Orfeo Toolbox  4.0
otbLocalRxDetectorFilter.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 __otbLocalRxDetectorFilter_txx
19 #define __otbLocalRxDetectorFilter_txx
20 
22 
23 namespace otb
24 {
25 
29 template <class TInputImage, class TOutputImage>
32  : m_InternalRadius(1), m_ExternalRadius(2)
33 {
34 }
35 
39 template <class TInputImage, class TOutputImage>
40 void
42 ::PrintSelf(std::ostream& os, itk::Indent indent) const
43 {
44  Superclass::PrintSelf(os, indent);
45 
46  os << indent << "Internal Radius: " << m_InternalRadius << std::endl;
47  os << indent << "External Radius: " << m_ExternalRadius << std::endl;
48 }
49 
53 template <class TInputImage, class TOutputImage>
54 void
57 {
58  // Get the input and output pointers
59  OutputPointerType outputPtr = this->GetOutput();
60 
61  // Fill the buffer with black pixels
62  outputPtr->FillBuffer(0);
63 }
64 
68 template <class TInputImage, class TOutputImage>
69 void
71 ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
72  itk::ThreadIdType threadId)
73 {
74  // Get the input and output pointers
75  InputConstPointerType inputPtr = this->GetInput();
76  OutputPointerType outputPtr = this->GetOutput();
77 
78  // Support progress methods/callbacks
79  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
80 
81  // Compute input region for thread
82  typename TInputImage::RegionType inputRegionForThread;
83  inputRegionForThread = outputRegionForThread;
84  inputRegionForThread.PadByRadius(m_ExternalRadius);
85  inputRegionForThread.Crop(inputPtr->GetLargestPossibleRegion());
86 
87  // Iterator on input region
89  radius.Fill(m_ExternalRadius);
90 
91  VectorFaceCalculatorType vectorFaceCalculator;
92  typename VectorFaceCalculatorType::FaceListType vectorFaceList;
93  typename VectorFaceCalculatorType::FaceListType::iterator vectorFit;
94 
95  vectorFaceList = vectorFaceCalculator(inputPtr, inputRegionForThread, radius);
96 
97  vectorFit = vectorFaceList.begin(); // Only the first face is used
98 
99  ConstShapedNeighborhoodIteratorType inputIt(radius, inputPtr, *vectorFit);
100 
101  // Neighborhood Configuration
103 
104  for (int y = -m_ExternalRadius; y <= m_ExternalRadius; y++)
105  {
106  off[1] = y;
107  for (int x = -m_ExternalRadius; x <= m_ExternalRadius; x++)
108  {
109  off[0] = x;
110  if ((abs(x) > m_InternalRadius) || (abs(y) > m_InternalRadius))
111  {
112  inputIt.ActivateOffset(off);
113  }
114  }
115  }
116 
117 
118  // iterator on output region
119  FaceCalculatorType faceCalculator;
120  typename FaceCalculatorType::FaceListType faceList;
121  typename FaceCalculatorType::FaceListType::iterator fit;
122 
123  faceList = faceCalculator(outputPtr, inputRegionForThread, radius);
124  fit = faceList.begin(); // Only the first face is used
125 
126  ImageRegionIteratorType outputIt(outputPtr, *fit);
127 
128  // Run Input Image
129  for (inputIt.GoToBegin(), outputIt.GoToBegin(); !inputIt.IsAtEnd(); ++inputIt, ++outputIt)
130  {
131  // Create ListSample
132  typename ListSampleType::Pointer listSample = ListSampleType::New();
133  listSample->SetMeasurementVectorSize(inputPtr->GetNumberOfComponentsPerPixel());
134 
135  // Run neighborhood
137  for (ci = inputIt.Begin(); !ci.IsAtEnd(); ++ci)
138  {
139  // Pushback element in listSample
140  listSample->PushBack(ci.Get());
141  }
142 
143  // Compute mean & covariance matrix
144  typename CovarianceCalculatorType::Pointer covarianceCalculator = CovarianceCalculatorType::New();
145  covarianceCalculator->SetInput(listSample);
146  covarianceCalculator->Update();
147 
148  MeasurementVectorRealType meanVector = covarianceCalculator->GetMean();
149  MatrixType covarianceMatrix = covarianceCalculator->GetCovarianceMatrix();
150 
151  // Compute RX value
152  typename MatrixType::InternalMatrixType invCovMat = covarianceMatrix.GetInverse();
153 
154  VectorMeasurementType testPixVec;
155  testPixVec = inputPtr->GetPixel(inputIt.GetIndex());
156 
157  VectorMeasurementType meanVec(meanVector.GetNumberOfElements());
158  for(unsigned int i = 0; i < meanVector.GetNumberOfElements(); ++i)
159  {
160  meanVec.SetElement(i, meanVector.GetElement(i));
161  }
162 
163  typename MatrixType::InternalMatrixType centeredTestPixMat(meanVector.GetNumberOfElements(), 1);
164  for (unsigned int i = 0; i < centeredTestPixMat.rows(); ++i)
165  {
166  centeredTestPixMat.put(i, 0, (testPixVec.GetElement(i) - meanVector.GetElement(i)));
167  }
168 
169  typename MatrixType::InternalMatrixType rxValue = centeredTestPixMat.transpose() * invCovMat * centeredTestPixMat;
170 
171  outputIt.Set(rxValue.get(0, 0));
172  }
173 }
174 
178 template <class TInputImage, class TOutputImage>
179 void
182 {
183  // call the superclass' implementation of this method
184  Superclass::GenerateInputRequestedRegion();
185 
186  // get pointers to the input and output
187  InputPointerType inputPtr =
188  const_cast< InputImageType * >( this->GetInput());
189  OutputPointerType outputPtr = this->GetOutput();
190 
191  if ( !inputPtr || !outputPtr )
192  {
193  return;
194  }
195 
196  // get a copy of the input requested region (should equal the output
197  // requested region)
198  typename TInputImage::RegionType inputRequestedRegion;
199  inputRequestedRegion = inputPtr->GetRequestedRegion();
200 
201  // pad the input requested region by the operator radius
202  inputRequestedRegion.PadByRadius( m_ExternalRadius );
203 
204  // crop the input requested region at the input's largest possible region
205  if ( inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()) )
206  {
207  inputPtr->SetRequestedRegion( inputRequestedRegion );
208  return;
209  }
210  else
211  {
212  // Couldn't crop the region (requested region is outside the largest
213  // possible region). Throw an exception.
214 
215  // store what we tried to request (prior to trying to crop)
216  inputPtr->SetRequestedRegion( inputRequestedRegion );
217 
218  // build an exception
219  itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
220  e.SetLocation(ITK_LOCATION);
221  e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
222  e.SetDataObject(inputPtr);
223  throw e;
224  }
225 }
226 
227 } // end namespace otb
228 
229 #endif
230 

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