OTB  6.7.0
Orfeo Toolbox
otbEstimateInnerProductPCAImageFilter.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 otbEstimateInnerProductPCAImageFilter_hxx
22 #define otbEstimateInnerProductPCAImageFilter_hxx
23 
25 
26 #include <vnl/algo/vnl_generalized_eigensystem.h>
27 #include <vnl/vnl_math.h>
28 
29 namespace otb
30 {
31 
35 template <class TInputImage, class TOutputImage>
37 ::EstimateInnerProductPCAImageFilter(): m_NumberOfPrincipalComponentsRequired(1), m_CenterData(true)
38 {
39 }
40 
44 template <class TInputImage, class TOutputImage>
45 void
47 ::PrintSelf(std::ostream& os, itk::Indent indent) const
48 {
49  this->Superclass::PrintSelf(os, indent);
50 }
51 
55 template<class TInputImage, class TOutputImage>
56 void
59 {
60  Superclass::GenerateOutputInformation();
61  this->GetOutput()->SetNumberOfComponentsPerPixel(m_NumberOfPrincipalComponentsRequired);
62 }
64 
68 template <class TInputImage, class TOutputImage>
69 void
72 {
73  // Instantiation object
74  StreamingInnerProductPointer streamingInnerProduct = StreamingInnerProductType::New();
75  streamingInnerProduct->SetInput(const_cast<InputImageType*>(this->GetInput()));
76  streamingInnerProduct->SetCenterData(m_CenterData);
77  streamingInnerProduct->Update();
78  m_InnerProduct = streamingInnerProduct->GetInnerProduct();
80 
81  MatrixType identityMatrix(m_InnerProduct.rows(), m_InnerProduct.columns());
82  identityMatrix.set_identity();
83  vnl_generalized_eigensystem eigenVectors_eigenValues(m_InnerProduct, identityMatrix);
84  m_EigenVectorsOfInnerProductMatrix = eigenVectors_eigenValues.V;
85 
86 }
87 
88 template <class TInputImage, class TOutputImage>
89 void
91 ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId)
92 {
93  typename InputImageType::ConstPointer inputPtr = this->GetInput();
94  typename OutputImageType::Pointer outputPtr = this->GetOutput();
95 
96  // Define the portion of the input to walk for this thread, using
97  // the CallCopyOutputRegionToInputRegion method allows for the input
98  // and output images to be different dimensions
99  unsigned int numberOfTrainingImages = inputPtr->GetNumberOfComponentsPerPixel();
100 
101  InputImageRegionType inputRegionForThread;
102  this->CallCopyOutputRegionToInputRegion(inputRegionForThread, outputRegionForThread);
103 
104  // Define the iterators
105  itk::ImageRegionConstIterator<TInputImage> inputIt(inputPtr, inputRegionForThread);
106  itk::ImageRegionIterator<TOutputImage> outputIt(outputPtr, outputRegionForThread);
107 
108  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
109 
110  inputIt.GoToBegin();
111  outputIt.GoToBegin();
112 
113  while (!inputIt.IsAtEnd())
114  {
115  InputPixelType inputPixel = inputIt.Get();
116  OutputPixelType outputPixel;
117  outputPixel.SetSize(m_NumberOfPrincipalComponentsRequired);
118  outputPixel.Fill(0);
119 
120  for (unsigned int img_number = 0; img_number < numberOfTrainingImages; ++img_number)
121  {
122  unsigned int indexNumberOfTrainingImages = numberOfTrainingImages - 1;
123  for (unsigned int vec_number = 0;
124  vec_number < m_NumberOfPrincipalComponentsRequired;
125  ++vec_number, indexNumberOfTrainingImages--)
126  {
127  outputPixel[vec_number] += static_cast<OutputInternalPixelType>(static_cast<double>(
128  inputPixel[img_number]) *
129  static_cast<double>(
130  m_EigenVectorsOfInnerProductMatrix[img_number
131  ][indexNumberOfTrainingImages]));
132  }
133  }
134  outputIt.Set(outputPixel);
135  ++inputIt;
136  ++outputIt;
137  progress.CompletedPixel(); // potential exception thrown here
138  }
139 
140 }
141 
142 }
143 
144 #endif
void Set(const PixelType &value) const
void PrintSelf(std::ostream &os, itk::Indent indent) const override
InputImageType::RegionType InputImageRegionType
unsigned int ThreadIdType
OutputImageType::RegionType OutputImageRegionType
bool IsAtEnd(void) const
void ThreadedGenerateData(const OutputImageRegionType &outputRegionForThread, itk::ThreadIdType threadId) override
PixelType Get(void) const