OTB  6.7.0
Orfeo Toolbox
otbPCAImageFilter.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 otbPCAImageFilter_hxx
22 #define otbPCAImageFilter_hxx
23 #include "otbPCAImageFilter.h"
24 
25 #include "itkMacro.h"
26 
27 #include <vnl/vnl_matrix.h>
28 #include <vnl/algo/vnl_symmetric_eigensystem.h>
29 #include <vnl/algo/vnl_generalized_eigensystem.h>
30 
31 namespace otb
32 {
33 
34 template < class TInputImage, class TOutputImage,
35  Transform::TransformDirection TDirectionOfTransformation >
38 {
39  this->SetNumberOfRequiredInputs(1);
40 
41  m_NumberOfPrincipalComponentsRequired = 0;
42 
43  m_UseNormalization = false;
44  m_UseVarianceForNormalization = false;
45  m_GivenMeanValues = false;
46  m_GivenStdDevValues = false;
47 
48  m_GivenCovarianceMatrix = false;
49  m_GivenTransformationMatrix = false;
50  m_IsTransformationMatrixForward = true;
51 
52  m_CovarianceEstimator = CovarianceEstimatorFilterType::New();
53  m_Transformer = TransformFilterType::New();
54  m_Transformer->MatrixByVectorOn();
55  m_Normalizer = NormalizeFilterType::New();
56 }
57 
58 template < class TInputImage, class TOutputImage,
59  Transform::TransformDirection TDirectionOfTransformation >
60 void
63 {
64  Superclass::GenerateOutputInformation();
65 
66  switch ( static_cast<int>(DirectionOfTransformation) )
67  {
68  case static_cast<int>(Transform::FORWARD):
69  {
70  if ( m_NumberOfPrincipalComponentsRequired == 0
71  || m_NumberOfPrincipalComponentsRequired
72  > this->GetInput()->GetNumberOfComponentsPerPixel() )
73  {
74  m_NumberOfPrincipalComponentsRequired =
75  this->GetInput()->GetNumberOfComponentsPerPixel();
76  }
77 
78  this->GetOutput()->SetNumberOfComponentsPerPixel(
79  m_NumberOfPrincipalComponentsRequired );
80  break;
81  }
82  case static_cast<int>(Transform::INVERSE):
83  {
84  unsigned int theOutputDimension = 0;
85  if ( m_GivenTransformationMatrix )
86  {
87  theOutputDimension = m_TransformationMatrix.Rows() >= m_TransformationMatrix.Cols() ?
88  m_TransformationMatrix.Rows() : m_TransformationMatrix.Cols();
89  }
90  else if ( m_GivenCovarianceMatrix )
91  {
92  theOutputDimension = m_CovarianceMatrix.Rows() >= m_CovarianceMatrix.Cols() ?
93  m_CovarianceMatrix.Rows() : m_CovarianceMatrix.Cols();
94  }
95  else
96  {
97  throw itk::ExceptionObject(__FILE__, __LINE__,
98  "Covariance or transformation matrix required to know the output size",
99  ITK_LOCATION);
100  }
101 
102  this->GetOutput()->SetNumberOfComponentsPerPixel( theOutputDimension );
103 
104  break;
105  }
106  default:
107  throw itk::ExceptionObject(__FILE__, __LINE__,
108  "Class should be templeted with FORWARD or INVERSE only...",
109  ITK_LOCATION );
110  }
111 
112 
113  switch ( static_cast<int>(DirectionOfTransformation) )
114  {
115  case static_cast<int>(Transform::FORWARD):
116  {
117  ForwardGenerateOutputInformation();
118  break;
119  }
120  case static_cast<int>(Transform::INVERSE):
121  {
122  ReverseGenerateOutputInformation();
123  break;
124  }
125  }
126 }
127 
128 template < class TInputImage, class TOutputImage,
129  Transform::TransformDirection TDirectionOfTransformation >
130 void
133 {
134  typename InputImageType::Pointer inputImgPtr
135  = const_cast<InputImageType*>( this->GetInput() );
136 
137  if ( !m_GivenTransformationMatrix )
138  {
139  if ( !m_GivenCovarianceMatrix )
140  {
141  if ( m_UseNormalization )
142  {
143  m_Normalizer->SetInput( inputImgPtr );
144  m_Normalizer->SetUseStdDev( m_UseVarianceForNormalization );
145 
146  if ( m_GivenMeanValues )
147  m_Normalizer->SetMean( m_MeanValues );
148 
149  if ( m_GivenStdDevValues )
150  m_Normalizer->SetStdDev( m_StdDevValues );
151 
152  m_Normalizer->GetOutput()->UpdateOutputInformation();
153 
154  if ( !m_GivenMeanValues )
155  {
156  m_MeanValues = m_Normalizer->GetCovarianceEstimator()->GetMean();
157 
158  if ( !m_GivenStdDevValues )
159  m_StdDevValues = m_Normalizer->GetFunctor().GetStdDev();
160 
161  if ( m_UseVarianceForNormalization )
162  m_CovarianceMatrix = m_Normalizer->GetCovarianceEstimator()->GetCorrelation();
163  else
164  m_CovarianceMatrix = m_Normalizer->GetCovarianceEstimator()->GetCovariance();
165  }
166  else
167  {
168  m_CovarianceEstimator->SetInput( m_Normalizer->GetOutput() );
169  m_CovarianceEstimator->UpdateOutputInformation();
170  m_CovarianceMatrix = m_CovarianceEstimator->GetCovariance();
171  }
172 
173  m_Transformer->SetInput( m_Normalizer->GetOutput() );
174  }
175  else
176  {
177  m_CovarianceEstimator->SetInput( inputImgPtr );
178  m_CovarianceEstimator->Update();
179 
180  m_CovarianceMatrix = m_CovarianceEstimator->GetCovariance();
181 
182  m_Transformer->SetInput( inputImgPtr );
183  }
184  }
185  else
186  {
187  m_Transformer->SetInput( inputImgPtr );
188  }
189 
190  GenerateTransformationMatrix();
191  }
192  else if ( !m_IsTransformationMatrixForward )
193  {
194  //m_TransformationMatrix = m_TransformationMatrix.GetTranspose();
195  m_TransformationMatrix = m_TransformationMatrix.GetInverse();
196 
197  m_Transformer->SetInput( inputImgPtr );
198  }
199 
200  if ( m_TransformationMatrix.GetVnlMatrix().empty() )
201  {
202  throw itk::ExceptionObject( __FILE__, __LINE__,
203  "Empty transformation matrix",
204  ITK_LOCATION);
205  }
206 }
207 
208 template < class TInputImage, class TOutputImage,
209  Transform::TransformDirection TDirectionOfTransformation >
210 void
213 {
214  if ( !m_GivenTransformationMatrix )
215  {
216  if ( !m_GivenCovarianceMatrix )
217  {
218  throw itk::ExceptionObject( __FILE__, __LINE__,
219  "No Covariance nor Transformation matrix given",
220  ITK_LOCATION );
221  }
222 
223  GenerateTransformationMatrix();
224  //m_TransformationMatrix = m_TransformationMatrix.GetTranspose();
225  m_TransformationMatrix = m_TransformationMatrix.GetInverse();
226  }
227  else if ( m_IsTransformationMatrixForward )
228  {
229  // prevents from multiple transpositions...
230  m_IsTransformationMatrixForward = false;
231  //m_TransformationMatrix = m_TransformationMatrix.GetTranspose();
232  m_TransformationMatrix = m_TransformationMatrix.GetInverse();
233  }
234 
235  if ( m_TransformationMatrix.GetVnlMatrix().empty() )
236  {
237  throw itk::ExceptionObject( __FILE__, __LINE__,
238  "Empty transformation matrix",
239  ITK_LOCATION);
240  }
241 
242  m_Transformer->SetInput( this->GetInput() );
243  m_Transformer->SetMatrix( m_TransformationMatrix.GetVnlMatrix() );
244  m_Normalizer->SetInput( m_Transformer->GetOutput() );
245 
246  if ( m_GivenStdDevValues || m_GivenMeanValues )
247  {
248  if ( m_GivenStdDevValues )
249  {
250  VectorType revStdDev ( m_StdDevValues.Size() );
251  for ( unsigned int i = 0; i < m_StdDevValues.Size(); ++i )
252  revStdDev[i] = 1. / m_StdDevValues[i];
253  m_Normalizer->SetStdDev( revStdDev );
254  }
255 
256  if ( m_GivenMeanValues )
257  {
258  VectorType revMean ( m_MeanValues.Size() );
259  if ( m_GivenStdDevValues )
260  {
261  for ( unsigned int i = 0; i < m_MeanValues.Size(); ++i )
262  revMean[i] = - m_MeanValues[i] / m_StdDevValues[i];
263  m_Normalizer->SetUseStdDev( true );
264  }
265  else
266  {
267  for ( unsigned int i = 0; i < m_MeanValues.Size(); ++i )
268  revMean[i] = -m_MeanValues[i];
269  m_Normalizer->SetUseStdDev( false );
270  }
271  m_Normalizer->SetMean( revMean );
272  }
273  }
274  else
275  {
276  m_Normalizer->SetUseMean(false);
277  m_Normalizer->SetUseStdDev(false);
278  }
279 }
280 
281 
282 template < class TInputImage, class TOutputImage,
283  Transform::TransformDirection TDirectionOfTransformation >
284 void
287 {
288  switch ( static_cast<int>(DirectionOfTransformation) )
289  {
290  case static_cast<int>(Transform::FORWARD):
291  return ForwardGenerateData();
292  case static_cast<int>(Transform::INVERSE):
293  return ReverseGenerateData();
294  default:
295  throw itk::ExceptionObject(__FILE__, __LINE__,
296  "Class should be templated with FORWARD or INVERSE only...",
297  ITK_LOCATION );
298  }
299 }
300 
301 template < class TInputImage, class TOutputImage,
302  Transform::TransformDirection TDirectionOfTransformation >
303 void
306 {
307  m_Transformer->SetMatrix( m_TransformationMatrix.GetVnlMatrix() );
308  m_Transformer->GraftOutput( this->GetOutput() );
309  m_Transformer->Update();
310  this->GraftOutput( m_Transformer->GetOutput() );
311 }
312 
313 template < class TInputImage, class TOutputImage,
314  Transform::TransformDirection TDirectionOfTransformation >
315 void
318 {
319 
320  if ( m_GivenStdDevValues || m_GivenMeanValues )
321  {
322  m_Normalizer->GraftOutput(this->GetOutput());
323  m_Normalizer->Update();
324  this->GraftOutput(m_Normalizer->GetOutput());
325  }
326  else
327  {
328  m_Transformer->GraftOutput( this->GetOutput() );
329  m_Transformer->Update();
330  this->GraftOutput( m_Transformer->GetOutput() );
331  }
332 }
333 
334 template < class TInputImage, class TOutputImage,
335  Transform::TransformDirection TDirectionOfTransformation >
336 void
339 {
340 #if 0
341  /*
342  * Old stuff but did work !
343  */
344 
345  MatrixType Id ( m_CovarianceMatrix );
346  Id.SetIdentity();
347 
348  typename MatrixType::InternalMatrixType A = m_CovarianceMatrix.GetVnlMatrix();
349  typename MatrixType::InternalMatrixType I = Id.GetVnlMatrix();
350 
351  vnl_generalized_eigensystem solver ( A, I );
352 
353  typename MatrixType::InternalMatrixType transf = solver.V;
354  transf.fliplr();
355  transf.inplace_transpose();
356 
357  vnl_vector< double > valP = solver.D.diagonal();
358  valP.flip();
359 
360  /*
361  * We used normalized PCA
362  */
363  for ( unsigned int i = 0; i < valP.size(); ++i )
364  {
365  if ( valP[i] != 0. )
366  valP[i] = 1. / std::sqrt( std::abs( valP[i] ) );
367  else
368  throw itk::ExceptionObject( __FILE__, __LINE__,
369  "Null Eigen value !!", ITK_LOCATION );
370  }
371  valP.post_multiply( transf );
372 
373  if ( m_NumberOfPrincipalComponentsRequired
374  != this->GetInput()->GetNumberOfComponentsPerPixel() )
375  m_TransformationMatrix = transf.get_n_rows( 0, m_NumberOfPrincipalComponentsRequired );
376  else
377  m_TransformationMatrix = transf;
378 
379 
380  m_EigenValues.SetSize( m_NumberOfPrincipalComponentsRequired );
381  for ( unsigned int i = 0; i < m_NumberOfPrincipalComponentsRequired; ++i )
382  m_EigenValues[i] = static_cast< RealType >( valP[i] );
383 #else
384  InternalMatrixType transf;
385  vnl_vector<double> vectValP;
386  vnl_symmetric_eigensystem_compute( m_CovarianceMatrix.GetVnlMatrix(), transf, vectValP );
387 
388  InternalMatrixType valP ( vectValP.size(), vectValP.size(), vnl_matrix_null );
389  for ( unsigned int i = 0; i < vectValP.size(); ++i )
390  valP(i, i) = vectValP[i];
391 
392  /* We used normalized PCA */
393  for ( unsigned int i = 0; i < valP.rows(); ++i )
394  {
395  if ( valP(i, i) > 0. )
396  {
397  valP(i, i) = 1. / std::sqrt( valP(i, i) );
398  }
399  else if ( valP(i, i) < 0. )
400  {
401  otbMsgDebugMacro( << "ValP(" << i << ") neg : " << valP(i, i) << " taking abs value" );
402  valP(i, i) = 1. / std::sqrt( std::abs( valP(i, i) ) );
403  }
404  else
405  {
406  throw itk::ExceptionObject( __FILE__, __LINE__,
407  "Null Eigen value !!", ITK_LOCATION );
408  }
409  }
410  transf = valP * transf.transpose();
411  transf.flipud();
412 
413  if ( m_NumberOfPrincipalComponentsRequired
414  != this->GetInput()->GetNumberOfComponentsPerPixel() )
415  m_TransformationMatrix = transf.get_n_rows( 0, m_NumberOfPrincipalComponentsRequired );
416  else
417  m_TransformationMatrix = transf;
418 
419  m_EigenValues.SetSize( m_NumberOfPrincipalComponentsRequired );
420  for ( unsigned int i = 0; i < m_NumberOfPrincipalComponentsRequired; ++i )
421  m_EigenValues[i] = static_cast< RealType >( valP(i, i) );
422 
423 #endif
424 }
425 
426 
427 template < class TInputImage, class TOutputImage,
428  Transform::TransformDirection TDirectionOfTransformation >
429 void
431 ::PrintSelf ( std::ostream& os, itk::Indent indent ) const
432 {
433  Superclass::PrintSelf( os, indent );
434 
435  os << indent << "m_UseNormalization = ";
436  if ( m_UseNormalization )
437  os << "true\n";
438  else
439  os << "false\n";
440 
441  if ( m_GivenMeanValues )
442  os << indent << "Given Mean : " << m_MeanValues << "\n";
443 
444  if ( m_GivenStdDevValues )
445  os << indent << "Given StdDev : " << m_StdDevValues << "\n";
446 
447  if ( !m_CovarianceMatrix.GetVnlMatrix().empty() )
448  {
449  os << indent << "Covariance matrix";
450  if ( m_GivenCovarianceMatrix )
451  os << " (given)";
452  os << "\n";
453 
454  m_CovarianceMatrix.GetVnlMatrix().print(os);
455 
456  if ( m_GivenCovarianceMatrix )
457  m_CovarianceEstimator->Print( os, indent.GetNextIndent() );
458  }
459 
460  if ( !m_TransformationMatrix.GetVnlMatrix().empty() )
461  {
462  os << indent;
463  if ( !m_IsTransformationMatrixForward )
464  os << "Invert ";
465  os << "Transformation matrix";
466  if ( m_GivenTransformationMatrix )
467  os << " (given)";
468  os << "\n";
469 
470  m_TransformationMatrix.GetVnlMatrix().print(os);
471  }
472 
473  if ( m_EigenValues.Size() > 0 )
474  {
475  os << indent << "Eigen value :";
476  for ( unsigned int i = 0; i < m_NumberOfPrincipalComponentsRequired; ++i )
477  os << " " << m_EigenValues[i];
478  os << "\n";
479  }
480 }
481 
482 } // end of namespace otb
483 
484 #endif // otbPCAImageFilter_hxx
485 
486 
virtual void ForwardGenerateOutputInformation()
void PrintSelf(std::ostream &os, itk::Indent indent) const override
MatrixObjectType::ComponentType MatrixType
void GenerateOutputInformation() override
virtual void ReverseGenerateOutputInformation()
#define otbMsgDebugMacro(x)
Definition: otbMacro.h:64
TInputImage InputImageType
virtual void ForwardGenerateData()
virtual void ReverseGenerateData()
CovarianceEstimatorFilterType::RealPixelType VectorType
MatrixType::InternalMatrixType InternalMatrixType
void GenerateData() override