OTB  6.5.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  // Run Input Image
116  for (inputIt.GoToBegin(), outputIt.GoToBegin(); !inputIt.IsAtEnd(); ++inputIt, ++outputIt)
117  {
118  // Create ListSample
119  typename ListSampleType::Pointer listSample = ListSampleType::New();
120  listSample->SetMeasurementVectorSize(inputPtr->GetNumberOfComponentsPerPixel());
121 
122  // Run neighborhood
124  for (ci = inputIt.Begin(); !ci.IsAtEnd(); ++ci)
125  {
126  // Pushback element in listSample
127  //std::cout << "pixel of shaped iteror : " << ci.Get() << std::endl;
128  listSample->PushBack(ci.Get());
129  }
130 
131  // Compute mean & covariance matrix
132  typename CovarianceCalculatorType::Pointer covarianceCalculator = CovarianceCalculatorType::New();
133  covarianceCalculator->SetInput(listSample);
134  covarianceCalculator->Update();
135 
136  MeasurementVectorRealType meanVector = covarianceCalculator->GetMean();
137  MatrixType covarianceMatrix = covarianceCalculator->GetCovarianceMatrix();
138 
139  // Compute RX value
140  typename MatrixType::InternalMatrixType invCovMat = covarianceMatrix.GetInverse();
141 
142  VectorMeasurementType testPixVec;
143  testPixVec = inputPtr->GetPixel(inputIt.GetIndex());
144 
145  VectorMeasurementType meanVec(meanVector.GetNumberOfElements());
146  for(unsigned int i = 0; i < meanVector.GetNumberOfElements(); ++i)
147  {
148  meanVec.SetElement(i, meanVector.GetElement(i));
149  }
150 
151  typename MatrixType::InternalMatrixType centeredTestPixMat(meanVector.GetNumberOfElements(), 1);
152  for (unsigned int i = 0; i < centeredTestPixMat.rows(); ++i)
153  {
154  centeredTestPixMat.put(i, 0, (testPixVec.GetElement(i) - meanVector.GetElement(i)));
155  }
156 
157  typename MatrixType::InternalMatrixType rxValue = centeredTestPixMat.transpose() * invCovMat * centeredTestPixMat;
158 
159  outputIt.Set(rxValue.get(0, 0));
160  }
161 }
162 
166 template <class TInputImage, class TOutputImage>
167 void
170 {
171  // Call the superclass' implementation of this method
172  Superclass::GenerateInputRequestedRegion();
173 
174  // Get pointers to the input and output
175  InputPointerType inputPtr =
176  const_cast< InputImageType * >( this->GetInput());
177  OutputPointerType outputPtr = this->GetOutput();
178 
179  if ( !inputPtr || !outputPtr )
180  {
181  return;
182  }
183 
184  // Get a copy of the input requested region (should equal the output
185  // requested region)
186  typename TInputImage::RegionType inputRequestedRegion;
187  inputRequestedRegion = inputPtr->GetRequestedRegion();
188 
189  // Pad the input requested region by the operator radius
190  inputRequestedRegion.PadByRadius( m_ExternalRadius );
191 
192  // Crop the input requested region at the input's largest possible region
193  if ( inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()) )
194  {
195  inputPtr->SetRequestedRegion( inputRequestedRegion );
196  return;
197  }
198  else
199  {
200  // Couldn't crop the region (requested region is outside the largest
201  // possible region). Throw an exception.
202 
203  // Store what we tried to request (prior to trying to crop)
204  inputPtr->SetRequestedRegion( inputRequestedRegion );
205 
206  // Build an exception
207  itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
208  e.SetLocation(ITK_LOCATION);
209  e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
210  e.SetDataObject(inputPtr);
211  throw e;
212  }
213 }
214 
215 } // end namespace otb
216 
217 #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