OTB  6.1.0
Orfeo Toolbox
otbVcaImageFilter.txx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 1999-2011 Insight Software Consortium
3  * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES)
4  *
5  * This file is part of Orfeo Toolbox
6  *
7  * https://www.orfeo-toolbox.org/
8  *
9  * Licensed under the Apache License, Version 2.0 (the "License");
10  * you may not use this file except in compliance with the License.
11  * You may obtain a copy of the License at
12  *
13  * http://www.apache.org/licenses/LICENSE-2.0
14  *
15  * Unless required by applicable law or agreed to in writing, software
16  * distributed under the License is distributed on an "AS IS" BASIS,
17  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
18  * See the License for the specific language governing permissions and
19  * limitations under the License.
20  */
21 
22 #ifndef otbVCAImageFilter_txx
23 #define otbVCAImageFilter_txx
24 
25 #include "otbVcaImageFilter.h"
27 
28 namespace otb {
29 
30 template<class TImage>
32  : m_NumberOfEndmembers(0)
33 {
34 }
35 
36 template<class TImage>
38 {
39 }
40 
41 template<class TImage>
43 {
44  // call the superclass' implementation of this method
45  Superclass::GenerateOutputInformation();
46 
47  // get pointers to the input and output
48  typename VectorImageType::ConstPointer inputPtr = this->GetInput();
49  typename VectorImageType::Pointer outputPtr = this->GetOutput();
50 
51  if ( !inputPtr || !outputPtr )
52  {
53  return;
54  }
55 
56  const unsigned int numberOfBands = inputPtr->GetNumberOfComponentsPerPixel();
57  RegionType outputRegion;
58  IndexType outputOrigin;
59  SizeType outputSize;
60 
61  outputOrigin.Fill(0);
62  outputSize[0] = m_NumberOfEndmembers;
63  outputSize[1] = 1;
64  outputRegion.SetIndex(outputOrigin);
65  outputRegion.SetSize(outputSize);
66 
67  outputPtr->SetLargestPossibleRegion( outputRegion );
68  outputPtr->SetNumberOfComponentsPerPixel( numberOfBands );
69 }
70 
71 template<class TImage>
73 {
74  typedef typename ForwardPCAImageFilterType::NormalizeFilterType NormalizeFilterType;
75 
76  VectorImageType* input = const_cast<VectorImageType*>(this->GetInput());
77  const unsigned int nbBands = this->GetInput()->GetNumberOfComponentsPerPixel();
78 
79 
80  otbMsgDevMacro( "Computing image stats" )
82  StreamingStatisticsVectorImageFilterType::New();
83 
84  statsInput->SetInput(input);
85  statsInput->SetEnableMinMax(false);
86  statsInput->Update();
87 
88  double SNR, SNRth;
89  vnl_matrix<PrecisionType> Ud;
90 
91  // SNR computation
92  {
93  vnl_matrix<PrecisionType> R = statsInput->GetCorrelation().GetVnlMatrix();
94  vnl_svd<PrecisionType> svd(R);
95  vnl_matrix<PrecisionType> U = svd.U();
96  Ud = U.get_n_columns(0, m_NumberOfEndmembers);
97  vnl_matrix<PrecisionType> Udt = Ud.transpose();
98 
99  // To remove the mean
100  typename NormalizeFilterType::Pointer normalize = NormalizeFilterType::New();
101  normalize->SetInput(input);
102  normalize->SetMean(statsInput->GetMean());
103  normalize->SetUseMean(true);
104  normalize->SetUseStdDev(false);
105 
106  typename MatrixImageFilterType::Pointer mulUd = MatrixImageFilterType::New();
107  mulUd->MatrixByVectorOn();
108  mulUd->SetInput(normalize->GetOutput());
109  mulUd->SetMatrix(Udt);
110  mulUd->UpdateOutputInformation();
111 
112  typename StreamingStatisticsVectorImageFilterType::Pointer transformedStats =
113  StreamingStatisticsVectorImageFilterType::New();
114 
115  transformedStats->SetInput(mulUd->GetOutput());
116  transformedStats->SetEnableMinMax(false);
117  transformedStats->Update();
118 
119  double P_R = nbBands * statsInput->GetComponentCorrelation();
120  double P_Rp = m_NumberOfEndmembers * transformedStats->GetComponentCorrelation()
121  + statsInput->GetMean().GetSquaredNorm();
122  //const double qL = static_cast<double>(m_NumberOfEndmembers) / nbBands;
123  SNR = vcl_abs(10*vcl_log10( (P_Rp - (m_NumberOfEndmembers/nbBands)*P_R) / (P_R - P_Rp) ));
124  }
125 
126  SNRth = 15.0 + 10.0 * vcl_log( static_cast<double>(m_NumberOfEndmembers) ) + 8.0;
127 
128  typename VectorImageType::Pointer Xd;
129  typename VectorImageType::Pointer Y;
130 
131  std::vector<itk::ProcessObject::Pointer> refHolder;
132  typename ForwardPCAImageFilterType::Pointer pca = ForwardPCAImageFilterType::New();
133  typename NormalizeFilterType::Pointer normalize = NormalizeFilterType::New();
134 
135  if (SNR > SNRth)
136  {
137  otbMsgDevMacro( "Using projective projection for dimensionnality reduction" )
138 
139  otbMsgDevMacro( "Computing SVD of correlation matrix" )
140  // Take the correlation matrix
141  vnl_matrix<PrecisionType> R = statsInput->GetCorrelation().GetVnlMatrix();
142 
143  // Apply SVD
144  vnl_svd<PrecisionType> svd(R);
145  vnl_matrix<PrecisionType> U = svd.U();
146  Ud = U.get_n_columns(0, m_NumberOfEndmembers);
147  vnl_matrix<PrecisionType> UdT = Ud.transpose();
148 
149  otbMsgDevMacro( "Apply dimensionality reduction" )
150  // Xd = Ud.'*M;
151  typename MatrixImageFilterType::Pointer mulUd = MatrixImageFilterType::New();
152  mulUd->MatrixByVectorOn();
153  mulUd->SetInput(this->GetInput());
154  mulUd->SetMatrix(UdT);
155  mulUd->UpdateOutputInformation();
156  refHolder.push_back(mulUd.GetPointer());
157 
158  Xd = mulUd->GetOutput();
159 
160  // Compute mean(Xd)
161  otbMsgDevMacro( "Compute mean(Xd)" )
163  statsXd = StreamingStatisticsVectorImageFilterType::New();
164  statsXd->SetEnableSecondOrderStats(false);
165  statsXd->SetInput(Xd);
166  // statsXd->GetStreamer()->SetNumberOfDivisionsStrippedStreaming(10);
167  statsXd->Update();
168  typename VectorImageType::PixelType Xdmean = statsXd->GetMean();
169  otbMsgDevMacro( "mean(Xd) = " << Xdmean)
170 
171  // Projective projection
172  // Xd ./ repmat( sum( Xd .* repmat(u, [1 N]) ) , [d 1]);
173  otbMsgDevMacro( "Compute projective projection" )
174  typename ProjectiveProjectionImageFilterType::Pointer proj = ProjectiveProjectionImageFilterType::New();
175  proj->SetInput(Xd);
176  proj->SetProjectionDirection(Xdmean);
177  refHolder.push_back(proj.GetPointer());
178 
179  Xd = proj->GetOutput();
180  Y = Xd;
181 
182  }
183  else
184  {
185  otbMsgDevMacro( "Using PCA for dimensionnality reduction" )
186 
187  normalize->SetInput(input);
188  normalize->SetMean(statsInput->GetMean());
189  normalize->SetUseMean(true);
190  normalize->SetUseStdDev(false);
191 
192  // Take the correlation matrix
193  vnl_matrix<PrecisionType> R = statsInput->GetCovariance().GetVnlMatrix();
194 
195  // Apply SVD
196  vnl_svd<PrecisionType> svd(R);
197  vnl_matrix<PrecisionType> U = svd.U();
198  Ud = U.get_n_columns(0, m_NumberOfEndmembers - 1);
199  vnl_matrix<PrecisionType> UdT = Ud.transpose();
200 
201  typename MatrixImageFilterType::Pointer mulUd = MatrixImageFilterType::New();
202  mulUd->MatrixByVectorOn();
203  mulUd->SetInput(normalize->GetOutput());
204  mulUd->SetMatrix(UdT);
205  mulUd->UpdateOutputInformation();
206  refHolder.push_back(mulUd.GetPointer());
207 
208  Xd = mulUd->GetOutput();
209 
210  typename VectorImageToAmplitudeImageFilterType::Pointer normComputer = VectorImageToAmplitudeImageFilterType::New();
211  normComputer->SetInput(Xd);
212 
213  typename StreamingMinMaxImageFilterType::Pointer maxNormComputer = StreamingMinMaxImageFilterType::New();
214  maxNormComputer->SetInput(normComputer->GetOutput());
215  maxNormComputer->Update();
216  typename ImageType::PixelType maxNorm = maxNormComputer->GetMaximum();
217  otbMsgDevMacro( "maxNorm : " << maxNorm)
218 
219  typename ConcatenateScalarValueImageFilterType::Pointer concat = ConcatenateScalarValueImageFilterType::New();
220  concat->SetInput(Xd);
221  concat->SetScalarValue(maxNorm);
222  refHolder.push_back(concat.GetPointer());
223 
224  Y = concat->GetOutput();
225  Y->UpdateOutputInformation();
226  }
227 
228  // E : result, will contain the endmembers
229  vnl_matrix<PrecisionType> E(nbBands, m_NumberOfEndmembers);
230 
231  // A = zeros(q, q)
232  // A(q, 1) = 1
233  vnl_matrix<PrecisionType> A(m_NumberOfEndmembers, m_NumberOfEndmembers);
234  A.fill(0);
235  A(m_NumberOfEndmembers - 1, 0) = 1;
236  typename RandomVariateGeneratorType::Pointer randomGen = RandomVariateGeneratorType::New();
237 
238  for (unsigned int i = 0; i < m_NumberOfEndmembers; ++i)
239  {
240  otbMsgDevMacro( "----------------------------------------" )
241  otbMsgDevMacro( "Iteration " << i )
242 
243  // w = rand(q, 1)
244  otbMsgDevMacro( "Random vector generation " )
245  vnl_vector<PrecisionType> w(m_NumberOfEndmembers);
246  for (unsigned int j = 0; j < w.size(); ++j)
247  {
248  w(j) = randomGen->GetVariateWithOpenRange();
249  }
250 
251  // f = ((I - A*pinv(A))*w) / (norm(I - A*pinv(A))*w))
252  otbMsgDevMacro( "f = ((I - A*pinv(A))*w) /(norm(I - A*pinv(A))*w))" )
253  vnl_matrix<PrecisionType> tmpMat(m_NumberOfEndmembers, m_NumberOfEndmembers);
254  tmpMat.set_identity();
255  otbMsgDevMacro( "A" << std::endl << A )
256  vnl_svd<PrecisionType> Asvd(A);
257  tmpMat -= A * Asvd.inverse();
258 
259  vnl_vector<PrecisionType> tmpNumerator = tmpMat * w;
260  vnl_vector<PrecisionType> f = tmpNumerator / tmpNumerator.two_norm();
261 
262  // v = f.'*Y
263  otbMsgDevMacro( "f = " << f );
264  otbMsgDevMacro( "v = f.'*Y" );
265  typename DotProductImageFilterType::Pointer dotfY = DotProductImageFilterType::New();
266  dotfY->SetInput(Y);
267 
268  typename VectorImageType::PixelType fV(f.data_block(), f.size());
269  dotfY->SetVector(typename VectorImageType::PixelType(fV));
270  typename ImageType::Pointer v = dotfY->GetOutput();
271 
272  // abs(v)
273  otbMsgDevMacro( "abs(v)" );
274  typename AbsImageFilterType::Pointer absVFilter = AbsImageFilterType::New();
275  absVFilter->SetInput(v);
276 
277  // max(abs(v))
278  otbMsgDevMacro( "max(abs(v))" );
279  typename StreamingMinMaxImageFilterType::Pointer maxAbs = StreamingMinMaxImageFilterType::New();
280  maxAbs->SetInput(absVFilter->GetOutput());
281  maxAbs->Update();
282 
283  // k = arg_max( max(abs(v)) )
284  otbMsgDevMacro( "k = arg_max( max(abs(v)) )" )
285  IndexType maxIdx = maxAbs->GetMaximumIndex();
286  otbMsgDevMacro( "max : " << maxAbs->GetMaximum() )
287  otbMsgDevMacro( "maxIdx : " << maxIdx )
288 
289  // extract Y(:, k)
290  otbMsgDevMacro( "Y(:, k)" )
291  RegionType region;
292  region.SetIndex( maxIdx );
293  SizeType size;
294  size.Fill(1);
295  region.SetSize(size);
296  Y->SetRequestedRegion(region);
297  Y->Update();
298 
299  // store new endmember in A
300  // A(:, i) = Y(:, k)
301  otbMsgDevMacro( "A(:, i) = Y(:, k)" )
302  typename VectorImageType::PixelType e = Y->GetPixel(maxIdx);
303  otbMsgDevMacro( "e = " << e )
304  A.set_column(i, e.GetDataPointer());
305 
306  otbMsgDevMacro( "A" << std::endl << A )
307 
308  // reproject new endmember in original space
309  vnl_vector<PrecisionType> u;
310  if (SNR > SNRth)
311  {
312  // u = Ud * Xd(:, k)
313  otbMsgDevMacro( "u = Ud * Xd(:, k)" )
314  Xd->SetRequestedRegion(region);
315  Xd->Update();
316  typename VectorImageType::PixelType xd = Xd->GetPixel(maxIdx);
317  vnl_vector<PrecisionType> xdV(xd.GetDataPointer(), xd.GetSize());
318  u = Ud * xdV;
319  }
320  else
321  {
322  // u = invPCA( Xd(:, k) )
323  otbMsgDevMacro( "u = Ud * Xd(:, k)" )
324  Xd->SetRequestedRegion(region);
325  Xd->Update();
326  typename VectorImageType::PixelType xd = Xd->GetPixel(maxIdx);
327  vnl_vector<PrecisionType> xdV(xd.GetDataPointer(), xd.GetSize());
328  //Ud = pca->GetTransformationMatrix().GetVnlMatrix().transpose();
329 
330  vnl_vector<PrecisionType> mean(statsInput->GetMean().GetDataPointer(), statsInput->GetMean().GetSize());
331  u = Ud * xdV + mean;
332  }
333 
334  // E(:, i) = u
335  otbMsgDevMacro( "E(:, i) = u" )
336  otbMsgDevMacro( "u = " << u )
337  E.set_column(i, u);
338  }
339 
340  typename VectorImageType::Pointer output = this->GetOutput();
341  output->SetRegions(output->GetLargestPossibleRegion());
342  //output->SetNumberOfComponentsPerPixel(input->GetNumberOfComponentsPerPixel());
343  output->Allocate();
344 
345  itk::ImageRegionIteratorWithIndex<VectorImageType> it(output, output->GetLargestPossibleRegion());
346  unsigned int i;
347  for (it.GoToBegin(), i = 0; !it.IsAtEnd(); ++it, ++i)
348  {
349  vnl_vector<PrecisionType> e = E.get_column(i);
350  typename VectorImageType::PixelType pixel(input->GetNumberOfComponentsPerPixel());
351  for (unsigned int j = 0; j < e.size(); ++j)
352  {
353  pixel[j] = E(j, i);
354  }
355  it.Set(pixel);
356  }
357 }
358 
362 template<class TImage>
363 void VCAImageFilter<TImage>::PrintSelf(std::ostream& os, itk::Indent indent) const
364 {
365  Superclass::PrintSelf( os, indent );
366 }
367 
368 } //end namespace otb
369 
370 #endif
~VCAImageFilter() ITK_OVERRIDE
ObjectType * GetPointer() const
VectorImageType::IndexType IndexType
void GenerateData() ITK_OVERRIDE
TVectorImage VectorImageType
VectorImageType::RegionType RegionType
void GenerateOutputInformation() ITK_OVERRIDE
void PrintSelf(std::ostream &os, itk::Indent indent) const ITK_OVERRIDE
Superclass::PixelType PixelType
Definition: otbImage.h:63
Normalize an VectorImage by setting its mean to zero and possibly variance to one (band by band)...
VectorImageType::SizeType SizeType
#define otbMsgDevMacro(x)
Definition: otbMacro.h:98