OTB  9.0.0
Orfeo Toolbox
otbRadiometricMomentsFunctor.h
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 otbRadiometricMomentsFunctor_h
22 #define otbRadiometricMomentsFunctor_h
23 
24 #include "itkVariableLengthVector.h"
25 #include "otbMath.h"
26 
27 namespace otb
28 {
29 namespace Functor
30 {
38 template <class TNeighIter, class TPrecision = float>
40 {
41 public:
43  typedef TNeighIter IteratorType;
44  typedef typename IteratorType::PixelType PixelType;
45  typedef TPrecision ScalarRealType;
46  typedef itk::VariableLengthVector<ScalarRealType> OutputType;
47 
49  {
50  }
52  {
53  }
54 
55  inline OutputType operator()(TNeighIter& it) const
56  {
57  OutputType moments;
58  moments.SetSize(4);
59  moments.Fill(itk::NumericTraits<ScalarRealType>::Zero);
60 
61  ScalarRealType sum1, sum2, sum3, sum4;
62  sum1 = itk::NumericTraits<ScalarRealType>::Zero;
63  sum2 = itk::NumericTraits<ScalarRealType>::Zero;
64  sum3 = itk::NumericTraits<ScalarRealType>::Zero;
65  sum4 = itk::NumericTraits<ScalarRealType>::Zero;
66 
67  // it.GoToBegin();
68  const unsigned int size = it.Size();
69  for (unsigned int i = 0; i < size; ++i)
70  {
71  ScalarRealType value = static_cast<ScalarRealType>(it.GetPixel(i));
72  ScalarRealType value2 = value * value;
73 
74  // Update cumulants
75  sum1 += value;
76  sum2 += value2;
77  sum3 += value * value2;
78  sum4 += value2 * value2;
79  }
80 
81  // final computations
82  // Mean
83  moments[0] = sum1 / size;
84  // Variance
85  moments[1] = (sum2 - (sum1 * moments[0])) / (size - 1);
86 
87  ScalarRealType sigma = std::sqrt(moments[1]);
88  ScalarRealType mean2 = moments[0] * moments[0];
89 
90  const double epsilon = 1E-10;
91  if (std::abs(moments[1]) > epsilon)
92  {
93  // Skewness
94  moments[2] = ((sum3 - 3.0 * moments[0] * sum2) / size + 2.0 * moments[0] * mean2) / (moments[1] * sigma);
95  // Kurtosis
96  moments[3] = ((sum4 - 4.0 * moments[0] * sum3 + 6.0 * mean2 * sum2) / size - 3.0 * mean2 * mean2) / (moments[1] * moments[1]) - 3.0;
97  }
98 
99  return moments;
100  }
101 };
102 }
103 }
104 
105 #endif
otb::Functor::RadiometricMomentsFunctor::ScalarRealType
TPrecision ScalarRealType
Definition: otbRadiometricMomentsFunctor.h:45
otb::Functor::RadiometricMomentsFunctor::operator()
OutputType operator()(TNeighIter &it) const
Definition: otbRadiometricMomentsFunctor.h:55
otbMath.h
otb
The "otb" namespace contains all Orfeo Toolbox (OTB) classes.
Definition: otbJoinContainer.h:32
otb::Functor::RadiometricMomentsFunctor::OutputType
itk::VariableLengthVector< ScalarRealType > OutputType
Definition: otbRadiometricMomentsFunctor.h:46
otb::Functor::RadiometricMomentsFunctor::RadiometricMomentsFunctor
RadiometricMomentsFunctor()
Definition: otbRadiometricMomentsFunctor.h:48
otb::Functor::RadiometricMomentsFunctor::PixelType
IteratorType::PixelType PixelType
Definition: otbRadiometricMomentsFunctor.h:44
otb::Functor::RadiometricMomentsFunctor::~RadiometricMomentsFunctor
~RadiometricMomentsFunctor()
Definition: otbRadiometricMomentsFunctor.h:51
otb::Functor::RadiometricMomentsFunctor::IteratorType
TNeighIter IteratorType
Definition: otbRadiometricMomentsFunctor.h:43
otb::Functor::RadiometricMomentsFunctor::Self
RadiometricMomentsFunctor Self
Definition: otbRadiometricMomentsFunctor.h:42
otb::Functor::RadiometricMomentsFunctor
Definition: otbRadiometricMomentsFunctor.h:39