OTB  6.7.0
Orfeo Toolbox
otbMaximumAutocorrelationFactorImageFilter.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2019 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 #ifndef otbMaximumAutocorrelationFactorImageFilter_hxx
22 #define otbMaximumAutocorrelationFactorImageFilter_hxx
23 
26 #include "otbMath.h"
27 #include "itkSubtractImageFilter.h"
28 
29 #include "vnl/algo/vnl_matrix_inverse.h"
30 #include "vnl/algo/vnl_generalized_eigensystem.h"
31 
33 #include "itkImageRegionIterator.h"
34 #include "itkProgressReporter.h"
35 
36 namespace otb
37 {
38 template <class TInputImage, class TOutputImage>
41 {
42  m_CovarianceEstimator = CovarianceEstimatorType::New();
43  m_CovarianceEstimatorH = CovarianceEstimatorType::New();
44  m_CovarianceEstimatorV = CovarianceEstimatorType::New();
45 }
46 
47 template <class TInputImage, class TOutputImage>
48 void
51 {
52  // Call superclass implementation
53  Superclass::GenerateOutputInformation();
54 
55  // Retrieve input images pointers
56  TInputImage * inputPtr = const_cast<TInputImage *>(this->GetInput());
57  //TOutputImage * outputPtr = this->GetOutput();
58 
59  // TODO: set the number of output components
60  unsigned int nbComp = inputPtr->GetNumberOfComponentsPerPixel();
61 
62  // Compute Dh and Dv
66 
67  InputImageRegionType largestInputRegion = inputPtr->GetLargestPossibleRegion();
68  InputImageRegionType referenceRegion;
69  InputImageSizeType size = largestInputRegion.GetSize();
70  size[0] -= 1;
71  size[1] -= 1;
72  referenceRegion.SetSize(size);
73  InputImageIndexType index = largestInputRegion.GetIndex();
74  referenceRegion.SetIndex(index);
75 
76  InputImageRegionType dhRegion;
77  InputImageRegionType dvRegion;
78 
79  index[0] += 1;
80 
81  dhRegion.SetSize(size);
82  dhRegion.SetIndex(index);
83 
84  index[0] -= 1;
85  index[1] += 1;
86 
87  dvRegion.SetSize(size);
88  dvRegion.SetIndex(index);
89 
90  typename ExtractFilterType::Pointer referenceExtract = ExtractFilterType::New();
91  referenceExtract->SetInput(inputPtr);
92  referenceExtract->SetExtractionRegion(referenceRegion);
93 
94  typename ExtractFilterType::Pointer dhExtract = ExtractFilterType::New();
95  dhExtract->SetInput(inputPtr);
96  dhExtract->SetExtractionRegion(dhRegion);
97 
99  dhExtractShift->SetInput(dhExtract->GetOutput());
100  dhExtractShift->SetReferenceImage(referenceExtract->GetOutput());
101  dhExtractShift->SetUseReferenceImage(true);
102  dhExtractShift->SetChangeOrigin(true);
103 
104  typename ExtractFilterType::Pointer dvExtract = ExtractFilterType::New();
105  dvExtract->SetInput(inputPtr);
106  dvExtract->SetExtractionRegion(dvRegion);
107 
109  dvExtractShift->SetInput(dvExtract->GetOutput());
110  dvExtractShift->SetReferenceImage(referenceExtract->GetOutput());
111  dvExtractShift->SetUseReferenceImage(true);
112  dvExtractShift->SetChangeOrigin(true);
113 
114 
115  typename DifferenceFilterType::Pointer diffh = DifferenceFilterType::New();
116  diffh->SetInput1(referenceExtract->GetOutput());
117  diffh->SetInput2(dhExtractShift->GetOutput());
118 
119  typename DifferenceFilterType::Pointer diffv = DifferenceFilterType::New();
120  diffv->SetInput1(referenceExtract->GetOutput());
121  diffv->SetInput2(dvExtractShift->GetOutput());
122 
123  //Compute pooled sigma (using sigmadh and sigmadv)
124  m_CovarianceEstimatorH->SetInput(diffh->GetOutput());
125  m_CovarianceEstimatorH->Update();
126  VnlMatrixType sigmadh = m_CovarianceEstimatorH->GetCovariance().GetVnlMatrix();
127 
128  m_CovarianceEstimatorV->SetInput(diffv->GetOutput());
129  m_CovarianceEstimatorV->Update();
130  VnlMatrixType sigmadv = m_CovarianceEstimatorV->GetCovariance().GetVnlMatrix();
131 
132  // Simple pool
133  VnlMatrixType sigmad = 0.5*(sigmadh+sigmadv);
134 
135  // Compute the original image covariance
136  referenceExtract->SetExtractionRegion(inputPtr->GetLargestPossibleRegion());
137  m_CovarianceEstimator->SetInput(referenceExtract->GetOutput());
138  m_CovarianceEstimator->Update();
139  VnlMatrixType sigma = m_CovarianceEstimator->GetCovariance().GetVnlMatrix();
140 
141  m_Mean = VnlVectorType(nbComp, 0);
142 
143  for(unsigned int i = 0; i<nbComp; ++i)
144  {
145  m_Mean[i] = m_CovarianceEstimator->GetMean()[i];
146  }
147 
148  vnl_generalized_eigensystem ges(sigmad, sigma);
149  VnlMatrixType d = ges.D;
150  m_V = ges.V;
151 
152  m_AutoCorrelation = VnlVectorType(nbComp, 1.);
153  m_AutoCorrelation -= 0.5 *d.get_diagonal();
154 
155  VnlMatrixType invstderr = VnlMatrixType(nbComp, nbComp, 0);
156  invstderr.set_diagonal(sigma.get_diagonal());
157  invstderr = invstderr.apply(&std::sqrt);
158  invstderr = invstderr.apply(&InverseValue);
159 
160  VnlMatrixType invstderrmaf = VnlMatrixType(nbComp, nbComp, 0);
161  invstderrmaf.set_diagonal((m_V.transpose() * sigma * m_V).get_diagonal());
162  invstderrmaf = invstderrmaf.apply(&std::sqrt);
163  invstderrmaf = invstderrmaf.apply(&InverseValue);
164 
165  VnlMatrixType aux1 = invstderr * sigma * m_V * invstderrmaf;
166 
167  VnlMatrixType sign = VnlMatrixType(nbComp, nbComp, 0);
168 
169  VnlVectorType aux2 = VnlVectorType(nbComp, 0);
170 
171  for(unsigned int i = 0; i < nbComp; ++i)
172  {
173  aux2=aux2 + aux1.get_row(i);
174  }
175 
176  sign.set_diagonal(aux2);
177  sign = sign.apply(&SignOfValue);
178 
179  // There is no need for scaling since vnl_generalized_eigensystem
180  // already gives unit variance
181  m_V = m_V * sign;
182 }
183 
184 template <class TInputImage, class TOutputImage>
185 void
187 ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId)
188 {
189  // Retrieve input images pointers
190  const TInputImage * inputPtr = this->GetInput();
191  TOutputImage * outputPtr = this->GetOutput();
192 
193 
194  typedef itk::ImageRegionConstIterator<InputImageType> ConstIteratorType;
195  typedef itk::ImageRegionIterator<OutputImageType> IteratorType;
196 
197  IteratorType outIt(outputPtr, outputRegionForThread);
198  ConstIteratorType inIt(inputPtr, outputRegionForThread);
199 
200  inIt.GoToBegin();
201  outIt.GoToBegin();
202 
203  // Get the number of components for each image
204  unsigned int outNbComp = outputPtr->GetNumberOfComponentsPerPixel();
205 
206 
207  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
208 
209  while(!inIt.IsAtEnd() && !outIt.IsAtEnd())
210  {
211  VnlVectorType x(outNbComp, 0);
212  VnlVectorType maf(outNbComp, 0);
213 
214  for(unsigned int i = 0; i < outNbComp; ++i)
215  {
216  x[i] = inIt.Get()[i];
217  }
218 
219  maf = (x-m_Mean)*m_V;
220 
221  typename OutputImageType::PixelType outPixel(outNbComp);
222 
223  for(unsigned int i = 0; i<outNbComp; ++i)
224  {
225  outPixel[i]=maf[i];
226  }
227 
228  outIt.Set(outPixel);
229 
230  ++inIt;
231  ++outIt;
232  progress.CompletedPixel();
233  }
234 }
235 }
236 
237 #endif
T SignOfValue(const T &value)
Definition: otbMath.h:99
Extract a spatial or spectral subset of a multi-channel image.
InputImageType::RegionType InputImageRegionType
T InverseValue(const T &value)
Definition: otbMath.h:94
unsigned int ThreadIdType
OutputImageType::RegionType OutputImageRegionType
void ThreadedGenerateData(const OutputImageRegionType &outputRegionForThread, itk::ThreadIdType threadId) override