OTB  5.0.0
Orfeo Toolbox
otbNormalizeInnerProductPCAImageFilter.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 __otbNormalizeInnerProductPCAImageFilter_txx
19 #define __otbNormalizeInnerProductPCAImageFilter_txx
20 
22 #include "itkImageRegionIterator.h"
23 #include "itkProgressReporter.h"
24 #include "itkNumericTraits.h"
25 
26 namespace otb
27 {
28 
32 template <class TInputImage, class TOutputImage>
35 {
36  this->SetNumberOfRequiredInputs(1);
37  this->SetNumberOfRequiredOutputs(1);
38  this->InPlaceOff();
40 
41  m_UseUnbiasedEstimator = true;
42 }
43 
47 template<class TInputImage, class TOutputImage>
48 void
51 {
52  Superclass::GenerateOutputInformation();
53 }
54 
58 template <class TInputImage, class TOutputImage>
59 void
62 {
63  StreamingStatisticsVectorImageFilterPointer stats = StreamingStatisticsVectorImageFilterType::New();
64  stats->SetInput(const_cast<InputImageType*>(this->GetInput()));
65  //set the normalization method to compute covariance to the StreamingStatisticsVectorImage filter
66  stats->SetUseUnbiasedEstimator(m_UseUnbiasedEstimator);
68 
69  stats->Update();
70 
71  RealPixelType means = stats->GetMean();
72  MatrixType cov = stats->GetCovariance();
73  double NbPixels = static_cast<double>(
74  this->GetInput()->GetLargestPossibleRegion().GetSize()[0] *
75  this->GetInput()->GetLargestPossibleRegion().GetSize()[1]);
76  if ((cov.Rows() != means.Size()) || (cov.Cols() != means.Size()))
77  {
78  itkExceptionMacro(<< "Covariance matrix with size (" << cov.Rows() << "," <<
79  cov.Cols() << ") is incompatible with mean vector with size" << means.Size());
80  }
81  m_CoefNorm.SetSize(means.Size());
82  for (unsigned int i = 0; i < m_CoefNorm.Size(); ++i)
83  {
84  m_CoefNorm[i] = (1. / vcl_sqrt(NbPixels * (cov[i][i] + means[i] * means[i])));
85  }
86 }
90 template <class TInputImage, class TOutputImage>
91 void
93 ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId)
94 {
95  typename Superclass::OutputImagePointer outputPtr = this->GetOutput();
96  typename Superclass::InputImageConstPointer inputPtr = this->GetInput();
98 
99  // Define the iterators
100  itk::ImageRegionConstIterator<InputImageType> inputIt(inputPtr, outputRegionForThread);
101  itk::ImageRegionIterator<OutputImageType> outputIt(outputPtr, outputRegionForThread);
102 
103  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
104 
105  inputIt.GoToBegin();
106  outputIt.GoToBegin();
107 
108  // Null pixel construction
109  InputPixelType nullPixel;
110  nullPixel.SetSize(inputPtr->GetNumberOfComponentsPerPixel());
112  while (!inputIt.IsAtEnd())
113  {
114  InputPixelType inPixel = inputIt.Get();
115  OutputPixelType outPixel;
116  outPixel.SetSize(inputPtr->GetNumberOfComponentsPerPixel());
118  //outPixel = m_Means * inPixel;
119  for (unsigned int j = 0; j < inputPtr->GetNumberOfComponentsPerPixel(); ++j)
120  {
121  outPixel[j] = static_cast<OutputInternalPixelType>(m_CoefNorm[j] * static_cast<double>(inPixel[j]));
122  }
123 
124  outputIt.Set(outPixel);
125  ++inputIt;
126  ++outputIt;
127  progress.CompletedPixel(); // potential exception thrown here
128  }
129 }
130 
131 template <class TInputImage, class TOutputImage>
132 void
134 ::PrintSelf(std::ostream& os, itk::Indent indent) const
135 {
136  this->Superclass::PrintSelf(os, indent);
137 }
138 
139 } // end namespace otb
140 
141 #endif
void Set(const PixelType &value) const
void PrintSelf(std::ostream &os, itk::Indent indent) const
void ThreadedGenerateData(const OutputImageRegionType &outputRegionForThread, itk::ThreadIdType threadId)
StreamingStatisticsVectorImageFilterType::RealPixelType RealPixelType
StreamingStatisticsVectorImageFilterType::MatrixType MatrixType
OutputImageType::RegionType OutputImageRegionType
bool IsAtEnd(void) const
PixelType Get(void) const
unsigned int ThreadIdType