OTB  6.7.0
Orfeo Toolbox
otbFastICAInternalOptimizerVectorImageFilter.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 otbFastICAInternalOptimizerVectorImageFilter_hxx
22 #define otbFastICAInternalOptimizerVectorImageFilter_hxx
24 
25 #include <itkMacro.h>
26 #include <itkImageRegionIterator.h>
27 
28 #include <vnl/vnl_math.h>
29 #include <vnl/vnl_matrix.h>
30 
31 namespace otb
32 {
33 
34 template <class TInputImage, class TOutputImage>
37 {
38  this->SetNumberOfRequiredInputs(2);
39 
40  m_CurrentBandForLoop = 0;
41  m_Beta = 0.;
42  m_Den = 0.;
43 
44  m_NonLinearity = [](double x) {return std::tanh(x);};
45  m_NonLinearityDerivative = [](double x) {return 1-std::pow( std::tanh(x), 2. );};
46 
47  m_TransformFilter = TransformFilterType::New();
48 }
49 
50 template <class TInputImage, class TOutputImage>
51 void
54 {
55  Superclass::GenerateOutputInformation();
56 
57  this->GetOutput()->SetNumberOfComponentsPerPixel(
58  this->GetInput(0)->GetNumberOfComponentsPerPixel() );
59 }
60 
61 template <class TInputImage, class TOutputImage>
62 void
65 {
66  if ( m_W.empty() )
67  {
68  throw itk::ExceptionObject( __FILE__, __LINE__,
69  "Give the initial W matrix", ITK_LOCATION );
70  }
71 
72  m_BetaVector.resize( this->GetNumberOfThreads() );
73  m_DenVector.resize( this->GetNumberOfThreads() );
74  m_NbSamples.resize( this->GetNumberOfThreads() );
75 
76  std::fill(m_BetaVector.begin(), m_BetaVector.end(), 0.);
77  std::fill(m_DenVector.begin(), m_DenVector.end(), 0.);
78  std::fill(m_NbSamples.begin(), m_NbSamples.end(), 0.);
79 }
80 
81 template <class TInputImage, class TOutputImage>
82 void
84 ::ThreadedGenerateData ( const OutputRegionType & outputRegionForThread, itk::ThreadIdType threadId )
85 {
86  InputRegionType inputRegion;
87  this->CallCopyOutputRegionToInputRegion( inputRegion, outputRegionForThread );
88 
90  ( this->GetInput(0), inputRegion );
92  ( this->GetInput(1), inputRegion );
94  ( this->GetOutput(), outputRegionForThread );
95 
96  unsigned int nbBands = this->GetInput(0)->GetNumberOfComponentsPerPixel();
97  input0It.GoToBegin();
98  input1It.GoToBegin();
99  outputIt.GoToBegin();
100 
101  double beta = 0.;
102  double den = 0.;
103  double nbSample = 0.;
104 
105  while ( !input0It.IsAtEnd() && !input1It.IsAtEnd() && !outputIt.IsAtEnd() )
106  {
107  double x = static_cast<double>( input1It.Get()[GetCurrentBandForLoop()] );
108  double g_x = m_NonLinearity(x);
109 
110  double x_g_x = x * g_x;
111  beta += x_g_x;
112 
113  double gp = m_NonLinearityDerivative(x);
114  den += gp;
115 
116  nbSample += 1.;
117 
118  typename OutputImageType::PixelType z ( nbBands );
119  for ( unsigned int bd = 0; bd < nbBands; bd++ )
120  z[bd] = g_x * input0It.Get()[bd];
121  outputIt.Set(z);
122 
123  ++input0It;
124  ++input1It;
125  ++outputIt;
126  } // end while loop
127 
128  m_BetaVector[threadId] += beta;
129  m_DenVector[threadId] += den;
130  m_NbSamples[threadId] += nbSample;
131 }
132 
133 template <class TInputImage, class TOutputImage>
134 void
137 {
138  double beta = 0;
139  double den = 0.;
140  double nbSample = 0;
141 
142  for ( itk::ThreadIdType i = 0; i < this->GetNumberOfThreads(); ++i )
143  {
144  beta += m_BetaVector[i];
145  den += m_DenVector[i];
146  nbSample += m_NbSamples[i];
147  }
148 
149  m_Beta = beta / nbSample;
150  m_Den = den / nbSample - m_Beta;
151 }
152 
153 } // end of namespace otb
154 
155 #endif
156 
157 
void ThreadedGenerateData(const OutputRegionType &, itk::ThreadIdType) override
void Set(const PixelType &value) const
unsigned int ThreadIdType
bool IsAtEnd(void) const
PixelType Get(void) const