Orfeo Toolbox  3.16
otbPCAImageFilter.txx
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Program: ORFEO Toolbox
4  Language: C++
5  Date: $Date$
6  Version: $Revision$
7 
8 
9  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
10  See OTBCopyright.txt for details.
11 
12 
13  This software is distributed WITHOUT ANY WARRANTY; without even
14  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
15  PURPOSE. See the above copyright notices for more information.
16 
17 =========================================================================*/
18 #ifndef __otbPCAImageFilter_txx
19 #define __otbPCAImageFilter_txx
20 #include "otbPCAImageFilter.h"
21 
22 #include "itkMacro.h"
23 
24 #include <vnl/vnl_matrix.h>
25 #include <vnl/algo/vnl_symmetric_eigensystem.h>
26 #include <vnl/algo/vnl_generalized_eigensystem.h>
27 
28 namespace otb
29 {
30 
31 template < class TInputImage, class TOutputImage,
32  Transform::TransformDirection TDirectionOfTransformation >
33 PCAImageFilter< TInputImage, TOutputImage, TDirectionOfTransformation >
34 ::PCAImageFilter ()
35 {
36  this->SetNumberOfRequiredInputs(1);
37 
38  m_NumberOfPrincipalComponentsRequired = 0;
39 
40  m_UseNormalization = false;
41  m_UseVarianceForNormalization = false;
42  m_GivenMeanValues = false;
43  m_GivenStdDevValues = false;
44 
45  m_GivenCovarianceMatrix = false;
46  m_GivenTransformationMatrix = false;
47  m_IsTransformationMatrixForward = true;
48 
49  m_CovarianceEstimator = CovarianceEstimatorFilterType::New();
50  m_Transformer = TransformFilterType::New();
51  m_Transformer->MatrixByVectorOn();
52  m_Normalizer = NormalizeFilterType::New();
53 }
54 
55 template < class TInputImage, class TOutputImage,
56  Transform::TransformDirection TDirectionOfTransformation >
57 void
58 PCAImageFilter< TInputImage, TOutputImage, TDirectionOfTransformation >
59 ::GenerateOutputInformation()
60 {
61  Superclass::GenerateOutputInformation();
62 
63  switch ( static_cast<int>(DirectionOfTransformation) )
64  {
65  case static_cast<int>(Transform::FORWARD):
66  {
67  if ( m_NumberOfPrincipalComponentsRequired == 0
68  || m_NumberOfPrincipalComponentsRequired
69  > this->GetInput()->GetNumberOfComponentsPerPixel() )
70  {
71  m_NumberOfPrincipalComponentsRequired =
72  this->GetInput()->GetNumberOfComponentsPerPixel();
73  }
74 
75  this->GetOutput()->SetNumberOfComponentsPerPixel(
76  m_NumberOfPrincipalComponentsRequired );
77  break;
78  }
79  case static_cast<int>(Transform::INVERSE):
80  {
81  unsigned int theOutputDimension = 0;
82  if ( m_GivenTransformationMatrix )
83  {
84  theOutputDimension = m_TransformationMatrix.Rows() >= m_TransformationMatrix.Cols() ?
85  m_TransformationMatrix.Rows() : m_TransformationMatrix.Cols();
86  }
87  else if ( m_GivenCovarianceMatrix )
88  {
89  theOutputDimension = m_CovarianceMatrix.Rows() >= m_CovarianceMatrix.Cols() ?
90  m_CovarianceMatrix.Rows() : m_CovarianceMatrix.Cols();
91  }
92  else
93  {
94  throw itk::ExceptionObject(__FILE__, __LINE__,
95  "Covariance or transformation matrix required to know the output size",
96  ITK_LOCATION);
97  }
98 
99  this->GetOutput()->SetNumberOfComponentsPerPixel( theOutputDimension );
100 
101  break;
102  }
103  default:
104  throw itk::ExceptionObject(__FILE__, __LINE__,
105  "Class should be templeted with FORWARD or INVERSE only...",
106  ITK_LOCATION );
107  }
108 }
109 
110 template < class TInputImage, class TOutputImage,
111  Transform::TransformDirection TDirectionOfTransformation >
112 void
113 PCAImageFilter< TInputImage, TOutputImage, TDirectionOfTransformation >
114 ::GenerateData ()
115 {
116  switch ( static_cast<int>(DirectionOfTransformation) )
117  {
118  case static_cast<int>(Transform::FORWARD):
119  return ForwardGenerateData();
120  case static_cast<int>(Transform::INVERSE):
121  return ReverseGenerateData();
122  default:
123  throw itk::ExceptionObject(__FILE__, __LINE__,
124  "Class should be templated with FORWARD or INVERSE only...",
125  ITK_LOCATION );
126  }
127 }
128 
129 template < class TInputImage, class TOutputImage,
130  Transform::TransformDirection TDirectionOfTransformation >
131 void
132 PCAImageFilter< TInputImage, TOutputImage, TDirectionOfTransformation >
133 ::ForwardGenerateData ()
134 {
135  typename InputImageType::Pointer inputImgPtr
136  = const_cast<InputImageType*>( this->GetInput() );
137 
138  if ( !m_GivenTransformationMatrix )
139  {
140  if ( !m_GivenCovarianceMatrix )
141  {
142  if ( m_UseNormalization )
143  {
144  m_Normalizer->SetInput( inputImgPtr );
145  m_Normalizer->SetUseStdDev( m_UseVarianceForNormalization );
146 
147  if ( m_GivenMeanValues )
148  m_Normalizer->SetMean( m_MeanValues );
149 
150  if ( m_GivenStdDevValues )
151  m_Normalizer->SetStdDev( m_StdDevValues );
152 
153  m_Normalizer->Update();
154 
155  if ( !m_GivenMeanValues )
156  {
157  m_MeanValues = m_Normalizer->GetCovarianceEstimator()->GetMean();
158 
159  if ( !m_GivenStdDevValues )
160  m_StdDevValues = m_Normalizer->GetFunctor().GetStdDev();
161 
162  if ( m_UseVarianceForNormalization )
163  m_CovarianceMatrix = m_Normalizer->GetCovarianceEstimator()->GetCorrelation();
164  else
165  m_CovarianceMatrix = m_Normalizer->GetCovarianceEstimator()->GetCovariance();
166  }
167  else
168  {
169  m_CovarianceEstimator->SetInput( m_Normalizer->GetOutput() );
170  m_CovarianceEstimator->Update();
171 
172  m_CovarianceMatrix = m_CovarianceEstimator->GetCovariance();
173  }
174 
175  m_Transformer->SetInput( m_Normalizer->GetOutput() );
176  }
177  else
178  {
179  m_CovarianceEstimator->SetInput( inputImgPtr );
180  m_CovarianceEstimator->Update();
181 
182  m_CovarianceMatrix = m_CovarianceEstimator->GetCovariance();
183 
184  m_Transformer->SetInput( inputImgPtr );
185  }
186  }
187  else
188  {
189  m_Transformer->SetInput( inputImgPtr );
190  }
191 
192  GenerateTransformationMatrix();
193  }
194  else if ( !m_IsTransformationMatrixForward )
195  {
196  //m_TransformationMatrix = m_TransformationMatrix.GetTranspose();
197  m_TransformationMatrix = m_TransformationMatrix.GetInverse();
198 
199  m_Transformer->SetInput( inputImgPtr );
200  }
201 
202  if ( m_TransformationMatrix.GetVnlMatrix().empty() )
203  {
204  throw itk::ExceptionObject( __FILE__, __LINE__,
205  "Empty transformation matrix",
206  ITK_LOCATION);
207  }
208 
209  m_Transformer->SetMatrix( m_TransformationMatrix.GetVnlMatrix() );
210  m_Transformer->GraftOutput( this->GetOutput() );
211  m_Transformer->Update();
212  this->GraftOutput( m_Transformer->GetOutput() );
213 }
214 
215 template < class TInputImage, class TOutputImage,
216  Transform::TransformDirection TDirectionOfTransformation >
217 void
218 PCAImageFilter< TInputImage, TOutputImage, TDirectionOfTransformation >
219 ::ReverseGenerateData ()
220 {
221  if ( !m_GivenTransformationMatrix )
222  {
223  if ( !m_GivenCovarianceMatrix )
224  {
225  throw itk::ExceptionObject( __FILE__, __LINE__,
226  "No Covariance nor Transformation matrix given",
227  ITK_LOCATION );
228  }
229 
230  GenerateTransformationMatrix();
231  //m_TransformationMatrix = m_TransformationMatrix.GetTranspose();
232  m_TransformationMatrix = m_TransformationMatrix.GetInverse();
233  }
234  else if ( m_IsTransformationMatrixForward )
235  {
236  // prevents from multiple transpositions...
237  m_IsTransformationMatrixForward = false;
238  //m_TransformationMatrix = m_TransformationMatrix.GetTranspose();
239  m_TransformationMatrix = m_TransformationMatrix.GetInverse();
240  }
241 
242  if ( m_TransformationMatrix.GetVnlMatrix().empty() )
243  {
244  throw itk::ExceptionObject( __FILE__, __LINE__,
245  "Empty transformation matrix",
246  ITK_LOCATION);
247  }
248 
249  m_Transformer->SetInput( this->GetInput() );
250  m_Transformer->SetMatrix( m_TransformationMatrix.GetVnlMatrix() );
251 
252  if ( m_GivenMeanValues || m_GivenStdDevValues )
253  {
254  m_Normalizer->SetInput( m_Transformer->GetOutput() );
255 
256  if ( m_GivenStdDevValues )
257  {
258  VectorType revStdDev ( m_StdDevValues.Size() );
259  for ( unsigned int i = 0; i < m_StdDevValues.Size(); ++i )
260  revStdDev[i] = 1. / m_StdDevValues[i];
261  m_Normalizer->SetStdDev( revStdDev );
262  }
263 
264  if ( m_GivenMeanValues )
265  {
266  VectorType revMean ( m_MeanValues.Size() );
267  if ( m_GivenStdDevValues )
268  {
269  for ( unsigned int i = 0; i < m_MeanValues.Size(); ++i )
270  revMean[i] = - m_MeanValues[i] / m_StdDevValues[i];
271  m_Normalizer->SetUseStdDev( true );
272  }
273  else
274  {
275  for ( unsigned int i = 0; i < m_MeanValues.Size(); ++i )
276  revMean[i] = -m_MeanValues[i];
277  m_Normalizer->SetUseStdDev( false );
278  }
279  m_Normalizer->SetMean( revMean );
280  }
281 
282  m_Normalizer->GraftOutput( this->GetOutput() );
283  m_Normalizer->Update();
284  this->GraftOutput( m_Normalizer->GetOutput() );
285  }
286  else
287  {
288  m_Transformer->GraftOutput( this->GetOutput() );
289  m_Transformer->Update();
290  this->GraftOutput( m_Transformer->GetOutput() );
291  }
292 }
293 
294 template < class TInputImage, class TOutputImage,
295  Transform::TransformDirection TDirectionOfTransformation >
296 void
297 PCAImageFilter< TInputImage, TOutputImage, TDirectionOfTransformation >
298 ::GenerateTransformationMatrix ()
299 {
300 #if 0
301  /*
302  * Old stuff but did work !
303  */
304 
305  MatrixType Id ( m_CovarianceMatrix );
306  Id.SetIdentity();
307 
308  typename MatrixType::InternalMatrixType A = m_CovarianceMatrix.GetVnlMatrix();
309  typename MatrixType::InternalMatrixType I = Id.GetVnlMatrix();
310 
311  vnl_generalized_eigensystem solver ( A, I );
312 
313  typename MatrixType::InternalMatrixType transf = solver.V;
314  transf.fliplr();
315  transf.inplace_transpose();
316 
317  vnl_vector< double > valP = solver.D.diagonal();
318  valP.flip();
319 
320  /*
321  * We used normalized PCA
322  */
323  for ( unsigned int i = 0; i < valP.size(); ++i )
324  {
325  if ( valP[i] != 0. )
326  valP[i] = 1. / vcl_sqrt( vcl_abs( valP[i] ) );
327  else
328  throw itk::ExceptionObject( __FILE__, __LINE__,
329  "Null Eigen value !!", ITK_LOCATION );
330  }
331  valP.post_multiply( transf );
332 
333  if ( m_NumberOfPrincipalComponentsRequired
334  != this->GetInput()->GetNumberOfComponentsPerPixel() )
335  m_TransformationMatrix = transf.get_n_rows( 0, m_NumberOfPrincipalComponentsRequired );
336  else
337  m_TransformationMatrix = transf;
338 
339 
340  m_EigenValues.SetSize( m_NumberOfPrincipalComponentsRequired );
341  for ( unsigned int i = 0; i < m_NumberOfPrincipalComponentsRequired; ++i )
342  m_EigenValues[i] = static_cast< RealType >( valP[i] );
343 #else
344  InternalMatrixType transf;
345  vnl_vector<double> vectValP;
346  vnl_symmetric_eigensystem_compute( m_CovarianceMatrix.GetVnlMatrix(), transf, vectValP );
347 
348  InternalMatrixType valP ( vectValP.size(), vectValP.size(), vnl_matrix_null );
349  for ( unsigned int i = 0; i < vectValP.size(); ++i )
350  valP(i, i) = vectValP[i];
351 
352  /* We used normalized PCA */
353  for ( unsigned int i = 0; i < valP.rows(); ++i )
354  {
355  if ( valP(i, i) > 0. )
356  {
357  valP(i, i) = 1. / vcl_sqrt( valP(i, i) );
358  }
359  else if ( valP(i, i) < 0. )
360  {
361  otbMsgDebugMacro( << "ValP(" << i << ") neg : " << valP(i, i) << " taking abs value" );
362  valP(i, i) = 1. / vcl_sqrt( vcl_abs( valP(i, i) ) );
363  }
364  else
365  {
366  throw itk::ExceptionObject( __FILE__, __LINE__,
367  "Null Eigen value !!", ITK_LOCATION );
368  }
369  }
370  transf = valP * transf.transpose();
371  transf.flipud();
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  m_EigenValues.SetSize( m_NumberOfPrincipalComponentsRequired );
380  for ( unsigned int i = 0; i < m_NumberOfPrincipalComponentsRequired; ++i )
381  m_EigenValues[i] = static_cast< RealType >( valP(i, i) );
382 
383 #endif
384 }
385 
386 
387 template < class TInputImage, class TOutputImage,
388  Transform::TransformDirection TDirectionOfTransformation >
389 void
390 PCAImageFilter< TInputImage, TOutputImage, TDirectionOfTransformation >
391 ::PrintSelf ( std::ostream& os, itk::Indent indent ) const
392 {
393  Superclass::PrintSelf( os, indent );
394 
395  os << indent << "m_UseNormalization = ";
396  if ( m_UseNormalization )
397  os << "true\n";
398  else
399  os << "false\n";
400 
401  if ( m_GivenMeanValues )
402  os << indent << "Given Mean : " << m_MeanValues << "\n";
403 
404  if ( m_GivenStdDevValues )
405  os << indent << "Given StdDev : " << m_StdDevValues << "\n";
406 
407  if ( !m_CovarianceMatrix.GetVnlMatrix().empty() )
408  {
409  os << indent << "Covariance matrix";
410  if ( m_GivenCovarianceMatrix )
411  os << " (given)";
412  os << "\n";
413 
414  m_CovarianceMatrix.GetVnlMatrix().print(os);
415 
416  if ( m_GivenCovarianceMatrix )
417  m_CovarianceEstimator->Print( os, indent.GetNextIndent() );
418  }
419 
420  if ( !m_TransformationMatrix.GetVnlMatrix().empty() )
421  {
422  os << indent;
423  if ( !m_IsTransformationMatrixForward )
424  os << "Invert ";
425  os << "Transformation matrix";
426  if ( m_GivenTransformationMatrix )
427  os << " (given)";
428  os << "\n";
429 
430  m_TransformationMatrix.GetVnlMatrix().print(os);
431  }
432 
433  if ( m_EigenValues.Size() > 0 )
434  {
435  os << indent << "Eigen value :";
436  for ( unsigned int i = 0; i < m_NumberOfPrincipalComponentsRequired; ++i )
437  os << " " << m_EigenValues[i];
438  os << "\n";
439  }
440 }
441 
442 } // end of namespace otb
443 
444 #endif // __otbPCAImageFilter_txx
445 
446 
Generated at Sun Feb 3 2013 00:41:22 for Orfeo Toolbox with doxygen 1.8.1.1