OTB  6.7.0
Orfeo Toolbox
otbRealMomentsImageFunction.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 otbRealMomentsImageFunction_hxx
22 #define otbRealMomentsImageFunction_hxx
23 
26 #include "itkNumericTraits.h"
27 
28 namespace otb
29 {
30 
34 template <class TInputImage, class TCoordRep>
37 {
38  m_NeighborhoodRadius = 1;
39  m_Pmax = 4;
40  m_Qmax = 4;
41 }
42 
43 template <class TInputImage, class TCoordRep>
44 void
46 ::PrintSelf(std::ostream& os, itk::Indent indent) const
47 {
48  this->Superclass::PrintSelf(os, indent);
49  os << indent << " p indice maximum value : " << m_Pmax << std::endl;
50  os << indent << " q indice maximum value : " << m_Qmax << std::endl;
51  os << indent << " Neighborhood radius value : " << m_NeighborhoodRadius << std::endl;
52 }
53 
54 template <class TInputImage, class TCoordRep>
57 ::EvaluateAtIndex(const IndexType& index) const
58 {
59  // Build moments vector
60  OutputType moments;
61  moments.resize(m_Pmax+1);
62 
63  std::vector<ScalarRealType> valXpY, valXqY;
64  valXpY.resize(m_Pmax+1);
65  valXqY.resize(m_Qmax+1);
66 
67  // Initialize moments
68  for (unsigned int p = 0; p <= m_Pmax; p++)
69  {
70  moments.at(p).resize(m_Qmax+1);
71  valXpY.at(p) = 1.0;
72  for (unsigned int q = 0; q <= m_Qmax; q++)
73  {
74  moments.at(p).at(q) = 0.0;
75  valXqY.at(q) = 1.0;
76  }
77  }
78 
79  // Check for input image
80  if( !this->GetInputImage() )
81  {
82  return moments;
83  }
84 
85  // Check for out of buffer
86  if ( !this->IsInsideBuffer( index ) )
87  {
88  return moments;
89  }
90 
91  // Create an N-d neighborhood kernel, using a zeroflux boundary condition
92  typename InputImageType::SizeType kernelSize;
93  kernelSize.Fill( m_NeighborhoodRadius );
94 
96  it(kernelSize, this->GetInputImage(), this->GetInputImage()->GetBufferedRegion());
97 
98  // Set the iterator at the desired location
99  it.SetLocation(index);
100 
101  // Walk the neighborhood
102  const unsigned int size = it.Size();
103  for (unsigned int i = 0; i < size; ++i)
104  {
105  // Retrieve value, and centered-reduced position
106  ScalarRealType value = static_cast<ScalarRealType>(it.GetPixel(i));
107  ScalarRealType x = static_cast<ScalarRealType>(it.GetOffset(i)[0])/(2*m_NeighborhoodRadius+1);
108  ScalarRealType y = static_cast<ScalarRealType>(it.GetOffset(i)[1])/(2*m_NeighborhoodRadius+1);
109 
110  unsigned int pTmp = 1;
111  unsigned int qTmp = 1;
112 
113  while (pTmp <= m_Pmax)
114  {
115  valXpY.at(pTmp) = valXpY.at(pTmp-1) * x;
116  pTmp ++;
117  }
118  while (qTmp <= m_Qmax)
119  {
120  valXqY.at(qTmp) = valXqY.at(qTmp-1) * y;
121  qTmp ++;
122  }
123 
124 
125  // Update cumulants
126  for (unsigned int p = 0; p <= m_Pmax; p++)
127  {
128  for (unsigned int q= 0; q <= m_Qmax; q++)
129  {
130  moments.at(p).at(q) += valXpY.at(p) * valXqY.at(q) * value;
131  }
132  }
133  }
134 
135  // Normalisation
136  for (int p = m_Pmax; p >= 0; p--)
137  {
138  for (int q= m_Qmax; q >= 0; q--)
139  {
140  moments.at(p).at(q) /= moments.at(0).at(0);
141  }
142  }
143 
144  // Return result
145  return moments;
146 }
147 
148 } // namespace otb
149 
150 #endif
itk::Size< Monteverdi_DIMENSION > SizeType
Definition: mvdTypes.h:137
OutputType EvaluateAtIndex(const IndexType &index) const override
void SetLocation(const IndexType &position)
void PrintSelf(std::ostream &os, itk::Indent indent) const override
virtual PixelType GetPixel(NeighborIndexType i, bool &IsInBounds) const
void Fill(SizeValueType value)