OTB  6.7.0
Orfeo Toolbox
otbRadiometricMomentsFunctor.h
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2017 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 
25 #include "otbMath.h"
26 
27 namespace otb
28 {
29 namespace Functor
30 {
38 template <class TNeighIter, class TPrecision = float>
40 {
41 public:
42 
44  typedef TNeighIter IteratorType;
46  typedef TPrecision ScalarRealType;
48 
51 
52  inline OutputType operator ()(TNeighIter& it) const
53  {
54  OutputType moments;
55  moments.SetSize(4);
57 
58  ScalarRealType sum1, sum2, sum3, sum4;
63 
64  //it.GoToBegin();
65  const unsigned int size = it.Size();
66  for (unsigned int i = 0; i < size; ++i)
67  {
68  ScalarRealType value = static_cast<ScalarRealType>(it.GetPixel(i));
69  ScalarRealType value2 = value * value;
70 
71  // Update cumulants
72  sum1 += value;
73  sum2 += value2;
74  sum3 += value * value2;
75  sum4 += value2 * value2;
76  }
77 
78  // final computations
79  // Mean
80  moments[0] = sum1 / size;
81  // Variance
82  moments[1] = (sum2 - (sum1 * moments[0])) / (size - 1);
83 
84  ScalarRealType sigma = std::sqrt(moments[1]);
85  ScalarRealType mean2 = moments[0] * moments[0];
86 
87  const double epsilon = 1E-10;
88  if (std::abs(moments[1]) > epsilon)
89  {
90  // Skewness
91  moments[2] = ((sum3 - 3.0 * moments[0] * sum2) / size + 2.0 * moments[0] * mean2) / (moments[1] * sigma);
92  // Kurtosis
93  moments[3] = ((sum4 - 4.0 * moments[0] * sum3 + 6.0 * mean2 * sum2) / size - 3.0 * mean2 * mean2)
94  / (moments[1] * moments[1]) - 3.0;
95  }
96 
97  return moments;
98  }
99 
100 };
101 
102 }
103 }
104 
105 #endif
OutputType operator()(TNeighIter &it) const
itk::VariableLengthVector< ScalarRealType > OutputType
void Fill(TValue const &v) noexcept
void SetSize(unsigned int sz, TReallocatePolicy reallocatePolicy, TKeepValuesPolicy keepValues)