ENH : texture formalas correction + add sum entropy

Code/FeatureExtraction/otbContrastTextureFunctor.h

@@ -120,7 +120,7 @@ public:
 	    double rVal = (static_cast<double>(r)+0.5)*binsLength[1];
 	    for (unsigned s = 0; s<histo[r].size(); s++)
 	      { 
-		if( vcl_abs( rVal - (static_cast<double>(s)+0.5)*binsLength[0] - nCeil) < vcl_abs(binsLength[0]-binsLength[1]) )
+		if( vcl_abs((static_cast<double>(s)+0.5)*binsLength[0] - rVal - nCeil) < vcl_abs(binsLength[0]) )
 		  {
 		    double p =  static_cast<double>(histo[r][s])*areaInv;
 		    out += nCeilSquare * p;

Code/FeatureExtraction/otbDifferenceEntropyTextureFunctor.h

@@ -115,20 +115,19 @@ public:
       { 
 	double Px_y = 0.;
 	double nCeil = (static_cast<double>(sB)+0.5)*binsLength[0];
-	double nCeilSquare = vcl_pow( nCeil, 2);
 	for (unsigned r = 0; r<histo.size(); r++)
 	  {
 	    double rVal = (static_cast<double>(r)+0.5)*binsLength[1];
 	    for (unsigned s = 0; s<histo[r].size(); s++)
 	      { 
-		if( vcl_abs( rVal - (static_cast<double>(s)+0.5)*binsLength[0] - nCeil ) < vcl_abs( binsLength[0]+binsLength[1]) )
+		if( vcl_abs((static_cast<double>(s)+0.5)*binsLength[0] - rVal - nCeil) < vcl_abs(binsLength[0]) )
 		  {
 		    Px_y += static_cast<double>(histo[r][s])*areaInv;
 		  }
 	      }
 	  }
 	if(Px_y != 0.)
-	  out += Px_y * nCeil * vcl_log(Px_y);
+	  out += Px_y * vcl_log(Px_y);
       }
     
     if(out != 0)

Code/FeatureExtraction/otbSumAverageTextureFunctor.h

@@ -122,10 +122,11 @@ public:
 		double sVal = (static_cast<double>(s)+0.5)*binsLength[0];
 		// In theory don't have the abs but will deals with neighborhood and offset without the same histo
 		// thus loop over 2*Ng don't have sense
-		if( vcl_abs(rVal + sVal - nCeil) < vcl_abs(binsLength[0]+binsLength[1]) || vcl_abs(rVal + sVal - 2*nCeil) < vcl_abs(binsLength[0]+binsLength[1]) )
+		//if( vcl_abs(rVal + sVal - nCeil) < vcl_abs(binsLength[0]+binsLength[1]) || vcl_abs(rVal + sVal - 2*nCeil) < vcl_abs(binsLength[0]+binsLength[1]) )
+		if( vcl_abs(rVal + sVal - nCeil) < vcl_abs(binsLength[0]) || vcl_abs(rVal + sVal - 2*nCeil) < 2*vcl_abs(binsLength[0]) )
 		  {
 		    double p =  static_cast<double>(histo[r][s])*areaInv;
-		    out += sVal * p;
+		    out += nCeil * p;
 		  }
 	      }
 	  }

Code/FeatureExtraction/otbSumEntropyTextureFunctor.h

+/*=========================================================================
+
+  Program:   ORFEO Toolbox
+  Language:  C++
+  Date:      $Date$
+  Version:   $Revision$
+
+
+  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
+  See OTBCopyright.txt for details.
+
+
+     This software is distributed WITHOUT ANY WARRANTY; without even
+     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+     PURPOSE.  See the above copyright notices for more information.
+
+=========================================================================*/
+#ifndef __otbSumEntropyTextureFunctor_h
+#define __otbSumEntropyTextureFunctor_h
+
+#include "otbTextureFunctorBase.h"
+
+namespace otb
+{
+namespace Functor
+{
+/** \class SumEntropyTextureFunctor
+ *  \brief This functor calculates the inverse difference moment of an image
+ *
+ *   Computes joint histogram (neighborhood and offset neighborhood) 
+ *   which bins are computing using Scott formula.
+ *   Computes the probabiltiy p for each pair of pixel.
+ *   InverseSumMoment  is the sum 1/(1+(pi-poff)²)*p over the neighborhood.
+ *   TIterInput is an ietrator, TOutput is a PixelType.
+ *
+ *  \ingroup Functor
+ *  \ingroup 
+ *  \ingroup Statistics
+ */
+template <class TIterInput1, class TIterInput2, class TOutput>
+class ITK_EXPORT SumEntropyTextureFunctor : 
+public TextureFunctorBase<TIterInput1, TIterInput2, TOutput>
+{
+public:
+  SumEntropyTextureFunctor(){};
+  ~SumEntropyTextureFunctor(){};
+
+  typedef TIterInput1                           IterType1;
+  typedef TIterInput2                           IterType2;
+  typedef TOutput                               OutputType;
+  typedef typename IterType1::OffsetType        OffsetType;
+  typedef typename IterType1::RadiusType        RadiusType;
+  typedef typename IterType1::InternalPixelType InternalPixelType;
+  typedef typename IterType1::ImageType         ImageType;
+  typedef itk::Neighborhood<InternalPixelType,::itk::GetImageDimension<ImageType>::ImageDimension>    NeighborhoodType;
+  typedef std::vector<double>                   DoubleVectorType;
+  typedef std::vector<int>                      IntVectorType;
+  typedef std::vector<IntVectorType>            IntVectorVectorType;
+
+
+  virtual double ComputeOverSingleChannel(const NeighborhoodType &neigh, const NeighborhoodType &neighOff)
+  {
+    DoubleVectorType binsLength = this->StatComputation(neigh, neighOff);
+
+    RadiusType radius = neigh.GetRadius();
+    double area = static_cast<double>(neigh.GetSize()[0]*neigh.GetSize()[1]);
+    double areaInv = 1/area;
+    OffsetType offset;
+    offset.Fill(0);
+    OffsetType offsetOff;
+    OffsetType offsetOffInit;
+
+    offsetOffInit[0] = -radius[0]+this->GetOffset()[0]-1;
+    offsetOffInit[1] = -radius[1]+this->GetOffset()[1]-1;
+
+    int histoIdX = 0;
+    int histoIdY = 0;
+    double out = 0.;
+
+    IntVectorType histoTemp;
+    IntVectorVectorType histo;
+    if (binsLength[0] != 0)
+      histoTemp = IntVectorType( vcl_floor( static_cast<double>(this->GetMaxi()-this->GetMini())/binsLength[0])+1., 0);
+    else
+      histoTemp = IntVectorType( 1, 0 );
+
+    if (binsLength[1] != 0)
+        histo = IntVectorVectorType( vcl_floor(static_cast<double>(this->GetMaxiOff()-this->GetMiniOff())/binsLength[1])+1., histoTemp );
+    else
+      histo = IntVectorVectorType( 1, histoTemp );
+
+    offsetOff = offsetOffInit;
+    for ( int l = -static_cast<int>(radius[0]); l <= static_cast<int>(radius[0]); l++ )
+	{
+	  offsetOff[0]++;
+	  offsetOff[1] = offsetOffInit[1];
+	  offset[0] = l;
+	  for ( int k = -static_cast<int>(radius[1]); k <= static_cast<int>(radius[1]); k++)
+	    {
+	      offsetOff[1]++;
+	      offset[1] = k;
+	      histoIdX = 0;
+	      histoIdY = 0;
+	      if ( binsLength[1] != 0)
+		histoIdX = static_cast<int>(vcl_floor( (static_cast<double>(neighOff[offsetOff])-this->GetMiniOff()) / static_cast<double>(binsLength[1]) ));
+	      if ( binsLength[0] !=0 )
+		histoIdY = static_cast<int>(vcl_floor( (static_cast<double>(neigh[offset])-this->GetMini()) /static_cast<double>( binsLength[0]) ));
+	      
+	      histo[histoIdX][histoIdY]++;
+	      
+	    }
+	}
+    // loop over bin neighborhood values
+    for (unsigned sB = 0; sB<histo[0].size(); sB++)
+      { 
+	double Px_y = 0.;
+	double nCeil = (static_cast<double>(sB)+0.5)*binsLength[0];
+	for (unsigned r = 0; r<histo.size(); r++)
+	  {
+	    double rVal = (static_cast<double>(r)+0.5)*binsLength[1];
+	    for (unsigned s = 0; s<histo[r].size(); s++)
+	      { 
+		double sVal = (static_cast<double>(s)+0.5)*binsLength[0];
+		if( vcl_abs(rVal + sVal - nCeil) < vcl_abs(binsLength[0]) || vcl_abs(rVal + sVal - 2*nCeil) < vcl_abs(binsLength[0]) )
+		  {
+		    Px_y +=  static_cast<double>(histo[r][s])*areaInv;
+		  }
+	      }
+	  }
+	if(Px_y != 0.)
+	  out += Px_y * vcl_log(Px_y);
+      }
+    
+      
+    if(out != 0)
+      out = -out;
+
+  
+    return out;  
+  }
+  
+};
+ 
+ 
+} // namespace Functor
+} // namespace otb
+
+#endif
+

Testing/Code/FeatureExtraction/CMakeLists.txt

@@ -1187,6 +1187,21 @@ ADD_TEST(feTvDifferenceEntropyTextureFunctor ${FEATUREEXTRACTION_TESTS12}
     2    # offset[1]	
 )
 
+# -------   otb::SumEntropyTextureFunctor  -------------
+ADD_TEST(feTvSumEntropyTextureFunctor ${FEATUREEXTRACTION_TESTS12} 
+--compare-image ${EPS}
+                ${BASELINE}/feTvSumEntropyTextureFunctor.tif
+                ${TEMP}/feTvSumEntropyTextureFunctor.tif
+  otbTextureFunctor
+    SEN #difference entropy
+    ${INPUTDATA}/poupees_sub.png
+    ${TEMP}/feTvSumEntropyTextureFunctor.tif
+    5    # radius
+    -2   # offset[0]
+    2    # offset[1]	
+)
+
+
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ otbFeatureExtractionTests13 ~~~~~~~~~~~~~~~~~~~~~
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -1328,6 +1343,23 @@ ADD_TEST(feTvDifferenceEntropyTextureImageFunction ${FEATUREEXTRACTION_TESTS13}
     2    # offset[1]
  )
 
+
+# -------        otb::SumEntropyTextureImageFunction   -------------
+ADD_TEST(feTvSumEntropyTextureImageFunction ${FEATUREEXTRACTION_TESTS13} 
+--compare-image ${EPS}
+                ${BASELINE}/feTuSumEntropyTextureImageFunction.tif
+                ${TEMP}/feTvSumEntropyTextureImageFunction.tif
+   otbTextureImageFunction
+    SEN #difference entropy
+    ${INPUTDATA}/poupees_sub_c1.png
+    ${TEMP}/feTvSumEntropyTextureImageFunction.tif
+    5    # radius[0]
+    5    # radius[1]
+    -2   # offset[0]
+    2    # offset[1]
+ )
+
+
 # A enrichir
 SET(BasicFeatureExtraction_SRCS1
 otbAlignImageToPath.cxx

Testing/Code/FeatureExtraction/otbTextureFunctor.cxx

@@ -33,6 +33,9 @@
 #include "otbContrastTextureFunctor.h"
 #include "otbSumAverageTextureFunctor.h"
 #include "otbDifferenceEntropyTextureFunctor.h"
+#include "otbSumEntropyTextureFunctor.h"
+
+
 
 template<class TInputImage, class TOutputImage, class TFunctor>
 int generic_TextureFunctor(int argc, char * argv[])
@@ -125,6 +128,11 @@ int otbTextureFunctor(int argc, char * argv[])
       typedef otb::Functor::DifferenceEntropyTextureFunctor<IteratorType, IteratorType, PixelType> FunctorType;
       return( generic_TextureFunctor<ImageType, ImageType, FunctorType>(argc,argv) );
     }
+  else if ( strArgv == "SEN" )
+    {
+      typedef otb::Functor::SumEntropyTextureFunctor<IteratorType, IteratorType, PixelType> FunctorType;
+      return( generic_TextureFunctor<ImageType, ImageType, FunctorType>(argc,argv) );
+    }
   else
     {
       return EXIT_FAILURE;

Testing/Code/FeatureExtraction/otbTextureImageFunction.cxx

@@ -36,6 +36,7 @@
 #include "otbContrastTextureFunctor.h"
 #include "otbSumAverageTextureFunctor.h"
 #include "otbDifferenceEntropyTextureFunctor.h"
+#include "otbSumEntropyTextureFunctor.h"
 
 
 template<class TInputImage, class TOutputImage, class TFunctor>
@@ -98,7 +99,6 @@ int otbTextureImageFunction(int argc, char * argv[])
 
   if(strArgv == "ENJ")
     {
-      std::cout<<"ENEGRY"<<std::endl;
       typedef otb::Functor::EnergyTextureFunctor<IteratorType, IteratorType, VectorType> FunctorType;
       return( generic_TextureImageFunction<ImageType, ImageType, FunctorType>(argc,argv) );
     }
@@ -142,6 +142,11 @@ int otbTextureImageFunction(int argc, char * argv[])
       typedef otb::Functor::DifferenceEntropyTextureFunctor<IteratorType, IteratorType, VectorType> FunctorType;
       return( generic_TextureImageFunction<ImageType, ImageType, FunctorType>(argc,argv) );
     }
+  else if ( strArgv == "SEN" )
+    {
+      typedef otb::Functor::SumEntropyTextureFunctor<IteratorType, IteratorType, VectorType> FunctorType;
+      return( generic_TextureImageFunction<ImageType, ImageType, FunctorType>(argc,argv) );
+    }
   else
     {
       return EXIT_FAILURE;