OTB  6.7.0
Orfeo Toolbox
otbNAPCAImageFilter.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 otbNAPCAImageFilter_hxx
22 #define otbNAPCAImageFilter_hxx
23 #include "otbNAPCAImageFilter.h"
24 
25 #include "itkMacro.h"
26 
27 #include <vnl/vnl_matrix.h>
28 #include <vnl/algo/vnl_matrix_inverse.h>
29 #include <vnl/algo/vnl_symmetric_eigensystem.h>
30 #include <vnl/algo/vnl_generalized_eigensystem.h>
31 
32 namespace otb
33 {
34 
35 template <class TInputImage, class TOutputImage,
36  class TNoiseImageFilter,
37  Transform::TransformDirection TDirectionOfTransformation >
38 void
41 {
42  InternalMatrixType An = this->GetNoiseCovarianceMatrix().GetVnlMatrix();
44  vnl_vector<double> vectValPn;
45  vnl_symmetric_eigensystem_compute( An, Fn, vectValPn );
46 
47  /* We used normalized PCA */
48  InternalMatrixType valPn ( vectValPn.size(), vectValPn.size(), vnl_matrix_null );
49  for ( unsigned int i = 0; i < valPn.rows(); ++i )
50  {
51  if ( vectValPn[i] > 0. )
52  {
53  valPn(i, i) = 1. / std::sqrt( vectValPn[i] );
54  }
55  else if ( vectValPn[i] < 0. )
56  {
57  otbMsgDebugMacro( << "ValPn(" << i << ") neg : " << vectValPn[i] << " taking abs value" );
58  valPn(i, i) = 1. / std::sqrt( std::abs( vectValPn[i] ) );
59  }
60  else
61  {
62  throw itk::ExceptionObject( __FILE__, __LINE__,
63  "Null Eigen value !!", ITK_LOCATION );
64  }
65  }
66  Fn = Fn * valPn;
67 
69  = vnl_matrix_inverse< MatrixElementType > ( this->GetCovarianceMatrix().GetVnlMatrix() );
70  InternalMatrixType Aadj = Fn.transpose() * Ax * Fn;
71 
72  InternalMatrixType Fadj;
73  vnl_vector<double> vectValPadj;
74  vnl_symmetric_eigensystem_compute( Aadj, Fadj, vectValPadj );
75 
76  InternalMatrixType transf = Fn * Fadj;
77  transf.inplace_transpose();
78 
79  if ( this->GetNumberOfPrincipalComponentsRequired()
80  != this->GetInput()->GetNumberOfComponentsPerPixel() )
81  this->m_TransformationMatrix = transf.get_n_rows( 0, this->GetNumberOfPrincipalComponentsRequired() );
82  else
83  this->m_TransformationMatrix = transf;
84 
85  this->m_EigenValues.SetSize( this->GetNumberOfPrincipalComponentsRequired() );
86  for ( unsigned int i = 0; i < this->GetNumberOfPrincipalComponentsRequired(); ++i )
87  this->m_EigenValues[this->GetNumberOfPrincipalComponentsRequired()-1-i]
88  = static_cast< RealType >( vectValPadj[i] );
89 }
90 
91 
92 } // end of namespace otb
93 
94 #endif
95 
96 
void GenerateTransformationMatrix() override
CovarianceEstimatorFilterType::RealType RealType
#define otbMsgDebugMacro(x)
Definition: otbMacro.h:64
MatrixType::InternalMatrixType InternalMatrixType