OTB  9.0.0
Orfeo Toolbox
otbHuMomentsImageFunction.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2022 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 otbHuMomentsImageFunction_hxx
22 #define otbHuMomentsImageFunction_hxx
23 
25 #include "itkConstNeighborhoodIterator.h"
26 #include "itkNumericTraits.h"
27 #include <complex>
28 
29 namespace otb
30 {
31 
35 template <class TInputImage, class TCoordRep>
37 {
38  m_NeighborhoodRadius = 1;
39 }
40 
41 template <class TInputImage, class TCoordRep>
42 void HuMomentsImageFunction<TInputImage, TCoordRep>::PrintSelf(std::ostream& os, itk::Indent indent) const
43 {
44  this->Superclass::PrintSelf(os, indent);
45  os << indent << " Neighborhood radius value : " << m_NeighborhoodRadius << std::endl;
46 }
47 
48 template <class TInputImage, class TCoordRep>
51 {
52  // Build moments vector
53  OutputType moments;
54 
55  // Initialize moments
56  moments.Fill(itk::NumericTraits<ScalarRealType>::Zero);
57 
58  // Check for input image
59  if (!this->GetInputImage())
60  {
61  return moments;
62  }
63 
64  // Check for out of buffer
65  if (!this->IsInsideBuffer(index))
66  {
67  return moments;
68  }
69 
70  // Define complex type
71  typedef std::complex<ScalarRealType> ComplexType;
72 
73  // Define and initialize cumulants for complex moments
74  ComplexType c11, c20, c02, c30, c03, c21, c12;
75  c11 = itk::NumericTraits<ComplexType>::Zero;
76  c20 = itk::NumericTraits<ComplexType>::Zero;
77  c02 = itk::NumericTraits<ComplexType>::Zero;
78  c30 = itk::NumericTraits<ComplexType>::Zero;
79  c03 = itk::NumericTraits<ComplexType>::Zero;
80  c21 = itk::NumericTraits<ComplexType>::Zero;
81  c12 = itk::NumericTraits<ComplexType>::Zero;
82 
83  ScalarRealType c00 = itk::NumericTraits<ScalarRealType>::Zero;
84 
85  // Create an N-d neighborhood kernel, using a zeroflux boundary condition
86  typename InputImageType::SizeType kernelSize;
87  kernelSize.Fill(m_NeighborhoodRadius);
88 
89  itk::ConstNeighborhoodIterator<InputImageType> it(kernelSize, this->GetInputImage(), this->GetInputImage()->GetBufferedRegion());
90 
91  // Set the iterator at the desired location
92  it.SetLocation(index);
93 
94  // Walk the neighborhood
95  const unsigned int size = it.Size();
96  for (unsigned int i = 0; i < size; ++i)
97  {
98  // Retrieve value, and centered-reduced position
99  ScalarRealType value = static_cast<ScalarRealType>(it.GetPixel(i));
100  ScalarRealType x = static_cast<ScalarRealType>(it.GetOffset(i)[0]) / (2 * m_NeighborhoodRadius + 1);
101  ScalarRealType y = static_cast<ScalarRealType>(it.GetOffset(i)[1]) / (2 * m_NeighborhoodRadius + 1);
102  // Build complex value
103  ComplexType xpy(x, y), xqy(x, -y);
104 
105  // Update cumulants
106  c00 += value;
107  c11 += xpy * xqy * value;
108  c20 += xpy * xpy * value;
109  c02 += xqy * xqy * value;
110  c30 += xpy * xpy * xpy * value;
111  c03 += xqy * xqy * xqy * value;
112  c21 += xpy * xpy * xqy * value;
113  c12 += xpy * xqy * xqy * value;
114  }
115 
116  // Nomalisation
117  c11 /= c00;
118  c20 /= c00;
119  c02 /= c00;
120  c30 /= c00;
121  c03 /= c00;
122  c21 /= c00;
123  c12 /= c00;
124 
125  // Compute moments combinations
126  moments[0] = static_cast<ScalarRealType>(c11.real());
127  moments[1] = static_cast<ScalarRealType>((c20 * c02).real());
128  moments[2] = static_cast<ScalarRealType>((c30 * c03).real());
129  moments[3] = static_cast<ScalarRealType>((c21 * c12).real());
130  moments[4] = static_cast<ScalarRealType>((c30 * c12 * c12 * c12).real());
131  moments[5] = static_cast<ScalarRealType>((c20 * c12 * c12).real());
132  moments[6] = static_cast<ScalarRealType>((c30 * c12 * c12 * c12).imag());
133 
134  // Return result
135  return moments;
136 }
137 
138 } // namespace otb
139 
140 #endif
otb::HuMomentsImageFunction::OutputType
Superclass::OutputType OutputType
Definition: otbHuMomentsImageFunction.h:89
otb
The "otb" namespace contains all Orfeo Toolbox (OTB) classes.
Definition: otbJoinContainer.h:32
otb::HuMomentsImageFunction::PrintSelf
void PrintSelf(std::ostream &os, itk::Indent indent) const override
Definition: otbHuMomentsImageFunction.hxx:42
otb::HuMomentsImageFunction::EvaluateAtIndex
OutputType EvaluateAtIndex(const IndexType &index) const override
Definition: otbHuMomentsImageFunction.hxx:50
otb::HuMomentsImageFunction::HuMomentsImageFunction
HuMomentsImageFunction()
Definition: otbHuMomentsImageFunction.hxx:36
otb::HuMomentsImageFunction::ScalarRealType
OutputType::ValueType ScalarRealType
Definition: otbHuMomentsImageFunction.h:90
otb::HuMomentsImageFunction::IndexType
Superclass::IndexType IndexType
Definition: otbHuMomentsImageFunction.h:85
otbHuMomentsImageFunction.h