OTB  6.1.0
Orfeo Toolbox
otbNormalizeInnerProductPCAImageFilter.txx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2017 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 otbNormalizeInnerProductPCAImageFilter_txx
22 #define otbNormalizeInnerProductPCAImageFilter_txx
23 
25 #include "itkImageRegionIterator.h"
26 #include "itkProgressReporter.h"
27 #include "itkNumericTraits.h"
28 
29 namespace otb
30 {
31 
35 template <class TInputImage, class TOutputImage>
38 {
39  this->SetNumberOfRequiredInputs(1);
40  this->SetNumberOfRequiredOutputs(1);
41  this->InPlaceOff();
43 
44  m_UseUnbiasedEstimator = true;
45 }
46 
50 template<class TInputImage, class TOutputImage>
51 void
54 {
55  Superclass::GenerateOutputInformation();
56 }
57 
61 template <class TInputImage, class TOutputImage>
62 void
65 {
66  StreamingStatisticsVectorImageFilterPointer stats = StreamingStatisticsVectorImageFilterType::New();
67  stats->SetInput(const_cast<InputImageType*>(this->GetInput()));
68  //set the normalization method to compute covariance to the StreamingStatisticsVectorImage filter
69  stats->SetUseUnbiasedEstimator(m_UseUnbiasedEstimator);
71 
72  stats->Update();
73 
74  RealPixelType means = stats->GetMean();
75  MatrixType cov = stats->GetCovariance();
76  double NbPixels = static_cast<double>(
77  this->GetInput()->GetLargestPossibleRegion().GetSize()[0] *
78  this->GetInput()->GetLargestPossibleRegion().GetSize()[1]);
79  if ((cov.Rows() != means.Size()) || (cov.Cols() != means.Size()))
80  {
81  itkExceptionMacro(<< "Covariance matrix with size (" << cov.Rows() << "," <<
82  cov.Cols() << ") is incompatible with mean vector with size" << means.Size());
83  }
84  m_CoefNorm.SetSize(means.Size());
85  for (unsigned int i = 0; i < m_CoefNorm.Size(); ++i)
86  {
87  m_CoefNorm[i] = (1. / vcl_sqrt(NbPixels * (cov[i][i] + means[i] * means[i])));
88  }
89 }
93 template <class TInputImage, class TOutputImage>
94 void
96 ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId)
97 {
98  typename Superclass::OutputImagePointer outputPtr = this->GetOutput();
99  typename Superclass::InputImageConstPointer inputPtr = this->GetInput();
101 
102  // Define the iterators
103  itk::ImageRegionConstIterator<InputImageType> inputIt(inputPtr, outputRegionForThread);
104  itk::ImageRegionIterator<OutputImageType> outputIt(outputPtr, outputRegionForThread);
105 
106  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
107 
108  inputIt.GoToBegin();
109  outputIt.GoToBegin();
110 
111  // Null pixel construction
112  InputPixelType nullPixel;
113  nullPixel.SetSize(inputPtr->GetNumberOfComponentsPerPixel());
115  while (!inputIt.IsAtEnd())
116  {
117  InputPixelType inPixel = inputIt.Get();
118  OutputPixelType outPixel;
119  outPixel.SetSize(inputPtr->GetNumberOfComponentsPerPixel());
121  //outPixel = m_Means * inPixel;
122  for (unsigned int j = 0; j < inputPtr->GetNumberOfComponentsPerPixel(); ++j)
123  {
124  outPixel[j] = static_cast<OutputInternalPixelType>(m_CoefNorm[j] * static_cast<double>(inPixel[j]));
125  }
126 
127  outputIt.Set(outPixel);
128  ++inputIt;
129  ++outputIt;
130  progress.CompletedPixel(); // potential exception thrown here
131  }
132 }
133 
134 template <class TInputImage, class TOutputImage>
135 void
137 ::PrintSelf(std::ostream& os, itk::Indent indent) const
138 {
139  this->Superclass::PrintSelf(os, indent);
140 }
141 
142 } // end namespace otb
143 
144 #endif
void Set(const PixelType &value) const
void ThreadedGenerateData(const OutputImageRegionType &outputRegionForThread, itk::ThreadIdType threadId) ITK_OVERRIDE
StreamingStatisticsVectorImageFilterType::RealPixelType RealPixelType
void PrintSelf(std::ostream &os, itk::Indent indent) const ITK_OVERRIDE
unsigned int ThreadIdType
StreamingStatisticsVectorImageFilterType::MatrixType MatrixType
OutputImageType::RegionType OutputImageRegionType
bool IsAtEnd(void) const
PixelType Get(void) const