Orfeo Toolbox  4.0
otbLocalRxDetectorNonThreadFilter.txx
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Program: ORFEO Toolbox
4  Language: C++
5  Date: $Date$
6  Version: $Revision$
7 
8  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
9  See OTBCopyright.txt for details.
10 
11 
12  This software is distributed WITHOUT ANY WARRANTY; without even
13  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
14  PURPOSE. See the above copyright notices for more information.
15 
16 =========================================================================*/
17 
18 #ifndef __otbLocalRxDetectorNonThreadFilter_txx
19 #define __otbLocalRxDetectorNonThreadFilter_txx
20 
22 
23 namespace otb
24 {
25 
29 template <class TInputImage, class TOutputImage>
32 {
33  this->m_ExternalRadius = 0;
34  this->m_InternalRadius = 0;
35 }
36 
40 template <class TInputImage, class TOutputImage>
41 void
43 ::PrintSelf(std::ostream& os, itk::Indent indent) const
44 {
45  Superclass::PrintSelf(os, indent);
46 
47  os << indent << "Internal Radius: " << m_InternalRadius << std::endl;
48  os << indent << "External Radius: " << m_ExternalRadius << std::endl;
49 }
50 
54 template <class TInputImage, class TOutputImage>
55 void
58 {
59 
60  this->AllocateOutputs();
61 
62  // Get the input and output pointers
63  InputPointerType inputPtr =
64  const_cast< InputImageType * >( this->GetInput());
65  OutputPointerType outputPtr = this->GetOutput();
66 
67  inputPtr->Update();
68 
69  // Fill the output buffer with black pixels
70  outputPtr->FillBuffer(0);
71 
72  // Iterator on input region
74  radius.Fill(m_ExternalRadius);
75 
76  VectorFaceCalculatorType vectorFaceCalculator;
77  typename VectorFaceCalculatorType::FaceListType vectorFaceList;
78  typename VectorFaceCalculatorType::FaceListType::iterator vectorFit;
79 
80  vectorFaceList = vectorFaceCalculator(inputPtr, inputPtr->GetRequestedRegion(), radius);
81  vectorFit = vectorFaceList.begin(); // Only the first face is used
82  ConstShapedNeighborhoodIteratorType inputIt(radius, inputPtr, *vectorFit);
83 
84  // Neighborhood Configuration
86 
87  for (int y = -m_ExternalRadius; y <= m_ExternalRadius; y++)
88  {
89  off[1] = y;
90  for (int x = -m_ExternalRadius; x <= m_ExternalRadius; x++)
91  {
92  off[0] = x;
93  if ((abs(x) > m_InternalRadius) || (abs(y) > m_InternalRadius))
94  {
95  inputIt.ActivateOffset(off);
96  }
97  }
98  }
99 
100  // Iterator on output region
101  FaceCalculatorType faceCalculator;
102  typename FaceCalculatorType::FaceListType faceList;
103  typename FaceCalculatorType::FaceListType::iterator fit;
104 
105  faceList = faceCalculator(outputPtr, inputPtr->GetRequestedRegion(), radius);
106  fit = faceList.begin(); // Only the first face is used
107 
108 
109  ImageRegionIteratorType outputIt(outputPtr, *fit);
110 
111  //std::cout << "Region de la face list : " << *vectorFit << std::endl;
112  //std::cout << "Region dans le buffer : " << inputPtr->GetBufferedRegion() << std::endl;
113 
114  // Run Input Image
115  for (inputIt.GoToBegin(), outputIt.GoToBegin(); !inputIt.IsAtEnd(); ++inputIt, ++outputIt)
116  {
117  // Create ListSample
118  typename ListSampleType::Pointer listSample = ListSampleType::New();
119  listSample->SetMeasurementVectorSize(inputPtr->GetNumberOfComponentsPerPixel());
120 
121  // Run neighborhood
123  for (ci = inputIt.Begin(); !ci.IsAtEnd(); ++ci)
124  {
125  // Pushback element in listSample
126  //std::cout << "pixel of shaped iteror : " << ci.Get() << std::endl;
127  listSample->PushBack(ci.Get());
128  }
129 
130  // Compute mean & covariance matrix
131  typename CovarianceCalculatorType::Pointer covarianceCalculator = CovarianceCalculatorType::New();
132  covarianceCalculator->SetInput(listSample);
133  covarianceCalculator->Update();
134 
135  MeasurementVectorRealType meanVector = covarianceCalculator->GetMean();
136  MatrixType covarianceMatrix = covarianceCalculator->GetCovarianceMatrix();
137 
138  // Compute RX value
139  typename MatrixType::InternalMatrixType invCovMat = covarianceMatrix.GetInverse();
140 
141  VectorMeasurementType testPixVec;
142  testPixVec = inputPtr->GetPixel(inputIt.GetIndex());
143 
144  VectorMeasurementType meanVec(meanVector.GetNumberOfElements());
145  for(unsigned int i = 0; i < meanVector.GetNumberOfElements(); ++i)
146  {
147  meanVec.SetElement(i, meanVector.GetElement(i));
148  }
149 
150  typename MatrixType::InternalMatrixType centeredTestPixMat(meanVector.GetNumberOfElements(), 1);
151  for (unsigned int i = 0; i < centeredTestPixMat.rows(); ++i)
152  {
153  centeredTestPixMat.put(i, 0, (testPixVec.GetElement(i) - meanVector.GetElement(i)));
154  }
155 
156  typename MatrixType::InternalMatrixType rxValue = centeredTestPixMat.transpose() * invCovMat * centeredTestPixMat;
157 
158  outputIt.Set(rxValue.get(0, 0));
159  }
160 }
161 
165 template <class TInputImage, class TOutputImage>
166 void
169 {
170  // Call the superclass' implementation of this method
171  Superclass::GenerateInputRequestedRegion();
172 
173  // Get pointers to the input and output
174  InputPointerType inputPtr =
175  const_cast< InputImageType * >( this->GetInput());
176  OutputPointerType outputPtr = this->GetOutput();
177 
178  if ( !inputPtr || !outputPtr )
179  {
180  return;
181  }
182 
183  // Get a copy of the input requested region (should equal the output
184  // requested region)
185  typename TInputImage::RegionType inputRequestedRegion;
186  inputRequestedRegion = inputPtr->GetRequestedRegion();
187 
188  // Pad the input requested region by the operator radius
189  inputRequestedRegion.PadByRadius( m_ExternalRadius );
190 
191  // Crop the input requested region at the input's largest possible region
192  if ( inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()) )
193  {
194  inputPtr->SetRequestedRegion( inputRequestedRegion );
195  return;
196  }
197  else
198  {
199  // Couldn't crop the region (requested region is outside the largest
200  // possible region). Throw an exception.
201 
202  // Store what we tried to request (prior to trying to crop)
203  inputPtr->SetRequestedRegion( inputRequestedRegion );
204 
205  // Build an exception
206  itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
207  e.SetLocation(ITK_LOCATION);
208  e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
209  e.SetDataObject(inputPtr);
210  throw e;
211  }
212 }
213 
214 } // end namespace otb
215 
216 #endif

Generated at Sat Mar 8 2014 16:06:31 for Orfeo Toolbox with doxygen 1.8.3.1