OTB  6.1.0
Orfeo Toolbox
otbLocalRxDetectorNonThreadFilter.txx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES)
3  *
4  * This file is part of Orfeo Toolbox
5  *
6  * https://www.orfeo-toolbox.org/
7  *
8  * Licensed under the Apache License, Version 2.0 (the "License");
9  * you may not use this file except in compliance with the License.
10  * You may obtain a copy of the License at
11  *
12  * http://www.apache.org/licenses/LICENSE-2.0
13  *
14  * Unless required by applicable law or agreed to in writing, software
15  * distributed under the License is distributed on an "AS IS" BASIS,
16  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17  * See the License for the specific language governing permissions and
18  * limitations under the License.
19  */
20 
21 
22 #ifndef otbLocalRxDetectorNonThreadFilter_txx
23 #define otbLocalRxDetectorNonThreadFilter_txx
24 
26 
27 namespace otb
28 {
29 
33 template <class TInputImage, class TOutputImage>
36 {
37  this->m_ExternalRadius = 0;
38  this->m_InternalRadius = 0;
39 }
40 
44 template <class TInputImage, class TOutputImage>
45 void
47 ::PrintSelf(std::ostream& os, itk::Indent indent) const
48 {
49  Superclass::PrintSelf(os, indent);
50 
51  os << indent << "Internal Radius: " << m_InternalRadius << std::endl;
52  os << indent << "External Radius: " << m_ExternalRadius << std::endl;
53 }
54 
58 template <class TInputImage, class TOutputImage>
59 void
62 {
63 
64  this->AllocateOutputs();
65 
66  // Get the input and output pointers
67  InputPointerType inputPtr =
68  const_cast< InputImageType * >( this->GetInput());
69  OutputPointerType outputPtr = this->GetOutput();
70 
71  inputPtr->Update();
72 
73  // Fill the output buffer with black pixels
74  outputPtr->FillBuffer(0);
75 
76  // Iterator on input region
78  radius.Fill(m_ExternalRadius);
79 
80  VectorFaceCalculatorType vectorFaceCalculator;
81  typename VectorFaceCalculatorType::FaceListType vectorFaceList;
82  typename VectorFaceCalculatorType::FaceListType::iterator vectorFit;
83 
84  vectorFaceList = vectorFaceCalculator(inputPtr, inputPtr->GetRequestedRegion(), radius);
85  vectorFit = vectorFaceList.begin(); // Only the first face is used
86  ConstShapedNeighborhoodIteratorType inputIt(radius, inputPtr, *vectorFit);
87 
88  // Neighborhood Configuration
90 
91  for (int y = -m_ExternalRadius; y <= m_ExternalRadius; y++)
92  {
93  off[1] = y;
94  for (int x = -m_ExternalRadius; x <= m_ExternalRadius; x++)
95  {
96  off[0] = x;
97  if ((abs(x) > m_InternalRadius) || (abs(y) > m_InternalRadius))
98  {
99  inputIt.ActivateOffset(off);
100  }
101  }
102  }
103 
104  // Iterator on output region
105  FaceCalculatorType faceCalculator;
106  typename FaceCalculatorType::FaceListType faceList;
107  typename FaceCalculatorType::FaceListType::iterator fit;
108 
109  faceList = faceCalculator(outputPtr, inputPtr->GetRequestedRegion(), radius);
110  fit = faceList.begin(); // Only the first face is used
111 
112 
113  ImageRegionIteratorType outputIt(outputPtr, *fit);
114 
115  //std::cout << "Region de la face list : " << *vectorFit << std::endl;
116  //std::cout << "Region dans le buffer : " << inputPtr->GetBufferedRegion() << std::endl;
117 
118  // Run Input Image
119  for (inputIt.GoToBegin(), outputIt.GoToBegin(); !inputIt.IsAtEnd(); ++inputIt, ++outputIt)
120  {
121  // Create ListSample
122  typename ListSampleType::Pointer listSample = ListSampleType::New();
123  listSample->SetMeasurementVectorSize(inputPtr->GetNumberOfComponentsPerPixel());
124 
125  // Run neighborhood
127  for (ci = inputIt.Begin(); !ci.IsAtEnd(); ++ci)
128  {
129  // Pushback element in listSample
130  //std::cout << "pixel of shaped iteror : " << ci.Get() << std::endl;
131  listSample->PushBack(ci.Get());
132  }
133 
134  // Compute mean & covariance matrix
135  typename CovarianceCalculatorType::Pointer covarianceCalculator = CovarianceCalculatorType::New();
136  covarianceCalculator->SetInput(listSample);
137  covarianceCalculator->Update();
138 
139  MeasurementVectorRealType meanVector = covarianceCalculator->GetMean();
140  MatrixType covarianceMatrix = covarianceCalculator->GetCovarianceMatrix();
141 
142  // Compute RX value
143  typename MatrixType::InternalMatrixType invCovMat = covarianceMatrix.GetInverse();
144 
145  VectorMeasurementType testPixVec;
146  testPixVec = inputPtr->GetPixel(inputIt.GetIndex());
147 
148  VectorMeasurementType meanVec(meanVector.GetNumberOfElements());
149  for(unsigned int i = 0; i < meanVector.GetNumberOfElements(); ++i)
150  {
151  meanVec.SetElement(i, meanVector.GetElement(i));
152  }
153 
154  typename MatrixType::InternalMatrixType centeredTestPixMat(meanVector.GetNumberOfElements(), 1);
155  for (unsigned int i = 0; i < centeredTestPixMat.rows(); ++i)
156  {
157  centeredTestPixMat.put(i, 0, (testPixVec.GetElement(i) - meanVector.GetElement(i)));
158  }
159 
160  typename MatrixType::InternalMatrixType rxValue = centeredTestPixMat.transpose() * invCovMat * centeredTestPixMat;
161 
162  outputIt.Set(rxValue.get(0, 0));
163  }
164 }
165 
169 template <class TInputImage, class TOutputImage>
170 void
173 {
174  // Call the superclass' implementation of this method
175  Superclass::GenerateInputRequestedRegion();
176 
177  // Get pointers to the input and output
178  InputPointerType inputPtr =
179  const_cast< InputImageType * >( this->GetInput());
180  OutputPointerType outputPtr = this->GetOutput();
181 
182  if ( !inputPtr || !outputPtr )
183  {
184  return;
185  }
186 
187  // Get a copy of the input requested region (should equal the output
188  // requested region)
189  typename TInputImage::RegionType inputRequestedRegion;
190  inputRequestedRegion = inputPtr->GetRequestedRegion();
191 
192  // Pad the input requested region by the operator radius
193  inputRequestedRegion.PadByRadius( m_ExternalRadius );
194 
195  // Crop the input requested region at the input's largest possible region
196  if ( inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()) )
197  {
198  inputPtr->SetRequestedRegion( inputRequestedRegion );
199  return;
200  }
201  else
202  {
203  // Couldn't crop the region (requested region is outside the largest
204  // possible region). Throw an exception.
205 
206  // Store what we tried to request (prior to trying to crop)
207  inputPtr->SetRequestedRegion( inputRequestedRegion );
208 
209  // Build an exception
210  itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
211  e.SetLocation(ITK_LOCATION);
212  e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
213  e.SetDataObject(inputPtr);
214  throw e;
215  }
216 }
217 
218 } // end namespace otb
219 
220 #endif
virtual IndexType GetIndex(void) const
virtual bool IsAtEnd() const
void SetDataObject(DataObject *dobj)
virtual void SetDescription(const std::string &s)
virtual void SetLocation(const std::string &s)
Superclass::RadiusType RadiusType
void Set(const PixelType &value) const
virtual void ActivateOffset(const OffsetType &off)
CovarianceCalculatorType::MeasurementVectorRealType MeasurementVectorRealType
TInputImage InputImageType
vnl_matrix< T > GetInverse(void) const
void PrintSelf(std::ostream &os, itk::Indent indent) const ITK_OVERRIDE
bool abs(const bool x)
const ConstIterator & Begin() const
vnl_matrix< T > InternalMatrixType