From d1c0b82b2b9c07fba0e8f4ee38ec6cd6018d7664 Mon Sep 17 00:00:00 2001
From: Aurelien Bricier <aurelien.bricier@c-s.fr>
Date: Mon, 12 Nov 2012 17:46:20 +0100
Subject: [PATCH] ENH: modified RadiometricMomentsImageFunction to use the
 functor

---
 .../otbRadiometricMomentsImageFunction.h      |  6 +++
 .../otbRadiometricMomentsImageFunction.txx    | 52 +++----------------
 2 files changed, 14 insertions(+), 44 deletions(-)

diff --git a/Code/FeatureExtraction/otbRadiometricMomentsImageFunction.h b/Code/FeatureExtraction/otbRadiometricMomentsImageFunction.h
index 8235798985..29416d75f5 100644
--- a/Code/FeatureExtraction/otbRadiometricMomentsImageFunction.h
+++ b/Code/FeatureExtraction/otbRadiometricMomentsImageFunction.h
@@ -19,6 +19,8 @@
 #define __otbRadiometricMomentsImageFunction_h
 
 #include "itkImageFunction.h"
+#include "otbRadiometricMomentsFunctor.h"
+#include "itkConstNeighborhoodIterator.h"
 #include "itkFixedArray.h"
 
 namespace otb
@@ -74,6 +76,9 @@ public:
 
   typedef TCoordRep                                CoordRepType;
 
+  typedef Functor::RadiometricMomentsFunctor< itk::ConstNeighborhoodIterator<InputImageType>, ScalarRealType>
+                                                   FunctorType;
+
   /** Dimension of the underlying image. */
   itkStaticConstMacro(ImageDimension, unsigned int,
                       InputImageType::ImageDimension);
@@ -112,6 +117,7 @@ private:
   void operator =(const Self&);  //purposely not implemented
 
   unsigned int m_NeighborhoodRadius;
+  FunctorType m_Functor;
 };
 
 } // namespace otb
diff --git a/Code/FeatureExtraction/otbRadiometricMomentsImageFunction.txx b/Code/FeatureExtraction/otbRadiometricMomentsImageFunction.txx
index b4bd675569..3062f05404 100644
--- a/Code/FeatureExtraction/otbRadiometricMomentsImageFunction.txx
+++ b/Code/FeatureExtraction/otbRadiometricMomentsImageFunction.txx
@@ -19,9 +19,9 @@
 #define __otbRadiometricMomentsImageFunction_txx
 
 #include "otbRadiometricMomentsImageFunction.h"
-#include "itkConstNeighborhoodIterator.h"
 #include "itkNumericTraits.h"
 #include "itkMacro.h"
+#include "itkVariableLengthVector.h"
 #include <complex>
 
 namespace otb
@@ -71,59 +71,23 @@ RadiometricMomentsImageFunction<TInputImage, TCoordRep>
     return moments;
     }
 
-  // Define and intialize cumulants
-  ScalarRealType sum1, sum2, sum3, sum4;
-  sum1 = itk::NumericTraits<ScalarRealType>::Zero;
-  sum2 = itk::NumericTraits<ScalarRealType>::Zero;
-  sum3 = itk::NumericTraits<ScalarRealType>::Zero;
-  sum4 = itk::NumericTraits<ScalarRealType>::Zero;
-
   // Create an N-d neighborhood kernel, using a zeroflux boundary condition
   typename InputImageType::SizeType kernelSize;
   kernelSize.Fill( m_NeighborhoodRadius );
-  
+
   itk::ConstNeighborhoodIterator<InputImageType>
-    it(kernelSize, this->GetInputImage(), this->GetInputImage()->GetBufferedRegion());
-  
+  it(kernelSize, this->GetInputImage(), this->GetInputImage()->GetBufferedRegion());
+
   // Set the iterator at the desired location
   it.SetLocation(index);
-  
-  // Walk the neighborhood
-  const unsigned int size = it.Size();
-  for (unsigned int i = 0; i < size; ++i)
-    {
-    ScalarRealType value = static_cast<ScalarRealType>(it.GetPixel(i));
-    ScalarRealType value2 = value * value;
 
-    // Update cumulants
-    sum1 += value;
-    sum2 += value2;
-    sum3 += value * value2;
-    sum4 += value2 * value2;
+  itk::VariableLengthVector<ScalarRealType> tmp(m_Functor(it));
 
-    }
-
-  // final computations
-
-  // Mean
-  moments[0] = sum1 / size;
-  // Variance
-  moments[1] = (sum2 - (sum1 * moments[0])) / (size - 1);
-  
-  
-  ScalarRealType sigma      = vcl_sqrt(moments[1]);
-  ScalarRealType mean2      = moments[0] * moments[0];
- 
-  const double epsilon = 1E-10;
-  if (vcl_abs(moments[1]) > epsilon)
+  for (unsigned int i=0; i<4;i++)
     {
-    // Skewness
-    moments[2] = ((sum3 - 3.0 * moments[0] * sum2) / size + 2.0 * moments[0] * mean2) / (moments[1] * sigma);
-    // Kurtosis
-    moments[3] = ((sum4 - 4.0 * moments[0] * sum3 + 6.0 * mean2 * sum2) / size - 3.0 * mean2 * mean2)
-      / (moments[1] * moments[1]) - 3.0;
+    moments[i] = tmp[i];
     }
-   
+
   // Return result
   return moments;
 }
-- 
GitLab