From 5d22317a72ab90a8856dfd51c60595efcf512b5e Mon Sep 17 00:00:00 2001
From: Rashad Kanavath <mohammed.rashad-km@cnes.fr>
Date: Sun, 4 May 2014 21:21:22 -0400
Subject: [PATCH] ENH: Add dissimilarity texture in to advanced haralick
 textures

---
 .../otbHaralickTextureExtraction.cxx          |  3 +-
 .../otbScalarImageToAdvancedTexturesFilter.h  |  5 ++
 ...otbScalarImageToAdvancedTexturesFilter.txx | 81 +++++++++++++------
 3 files changed, 61 insertions(+), 28 deletions(-)

diff --git a/Applications/FeatureExtraction/otbHaralickTextureExtraction.cxx b/Applications/FeatureExtraction/otbHaralickTextureExtraction.cxx
index 21c4bd282e..d7664f249c 100644
--- a/Applications/FeatureExtraction/otbHaralickTextureExtraction.cxx
+++ b/Applications/FeatureExtraction/otbHaralickTextureExtraction.cxx
@@ -28,8 +28,6 @@
 #include "otbImageList.h"
 #include "otbImageListToVectorImageFilter.h"
 
-#include "itkTimeProbe.h"
-
 namespace otb
 {
 namespace Wrapper
@@ -239,6 +237,7 @@ void DoExecute()
     m_AdvTexFilter->SetNumberOfBinsPerAxis(GetParameterInt("parameters.nbbin"));
     m_AdvImageList->PushBack(m_AdvTexFilter->GetMeanOutput());
     m_AdvImageList->PushBack(m_AdvTexFilter->GetVarianceOutput());
+    m_AdvImageList->PushBack(m_AdvTexFilter->GetDissimilarityOutput());
     m_AdvImageList->PushBack(m_AdvTexFilter->GetSumAverageOutput());
     m_AdvImageList->PushBack(m_AdvTexFilter->GetSumVarianceOutput());
     m_AdvImageList->PushBack(m_AdvTexFilter->GetSumEntropyOutput());
diff --git a/Code/FeatureExtraction/otbScalarImageToAdvancedTexturesFilter.h b/Code/FeatureExtraction/otbScalarImageToAdvancedTexturesFilter.h
index 7b91d5d0e1..c7b6571e50 100644
--- a/Code/FeatureExtraction/otbScalarImageToAdvancedTexturesFilter.h
+++ b/Code/FeatureExtraction/otbScalarImageToAdvancedTexturesFilter.h
@@ -36,6 +36,8 @@ namespace otb
  *
  * "Sum of squares: Variance" \f$ = f_4 = \sum_{i, j}(i - \mu)^2 g(i, j) \f$
  *
+ * "Dissimilarity" \f$ = f_5 = \sum_{i, j}(i - j) g(i, j)^2 \f$
+ *
  * "Sum average" \f$ = f_6 = -\sum_{i}i g_{x+y}(i)
  *
  * "Sum Variance" \f$ = f_7 = \sum_{i}(i - f_8)^2 g_{x+y}(i) \f$
@@ -178,6 +180,9 @@ public:
   /** Get the variance output image */
   OutputImageType * GetVarianceOutput();
 
+  /** Get the dissimilarity output image */
+  OutputImageType * GetDissimilarityOutput();
+
   /** Get the sum average output image */
   OutputImageType * GetSumAverageOutput();
 
diff --git a/Code/FeatureExtraction/otbScalarImageToAdvancedTexturesFilter.txx b/Code/FeatureExtraction/otbScalarImageToAdvancedTexturesFilter.txx
index c2866ef261..0593b93754 100644
--- a/Code/FeatureExtraction/otbScalarImageToAdvancedTexturesFilter.txx
+++ b/Code/FeatureExtraction/otbScalarImageToAdvancedTexturesFilter.txx
@@ -37,10 +37,10 @@ ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
   m_InputImageMinimum(0),
   m_InputImageMaximum(255)
 {
-  // There are 9 outputs corresponding to the 8 textures indices
-  this->SetNumberOfRequiredOutputs(9);
+  // There are 10 outputs corresponding to the 9 textures indices
+  this->SetNumberOfRequiredOutputs(10);
 
-  // Create the 9 outputs
+  // Create the 10 outputs
   this->SetNthOutput(0, OutputImageType::New());
   this->SetNthOutput(1, OutputImageType::New());
   this->SetNthOutput(2, OutputImageType::New());
@@ -50,18 +50,18 @@ ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
   this->SetNthOutput(6, OutputImageType::New());
   this->SetNthOutput(7, OutputImageType::New());
   this->SetNthOutput(8, OutputImageType::New());
+  this->SetNthOutput(9, OutputImageType::New());
 }
 
 template <class TInputImage, class TOutputImage>
 ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
 ::~ScalarImageToAdvancedTexturesFilter()
 {}
-
 template <class TInputImage, class TOutputImage>
 typename ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
 ::OutputImageType *
 ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
-::GetVarianceOutput()
+::GetMeanOutput()
 {
   if (this->GetNumberOfOutputs() < 1)
     {
@@ -74,7 +74,7 @@ template <class TInputImage, class TOutputImage>
 typename ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
 ::OutputImageType *
 ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
-::GetSumAverageOutput()
+::GetVarianceOutput()
 {
   if (this->GetNumberOfOutputs() < 2)
     {
@@ -87,7 +87,7 @@ template <class TInputImage, class TOutputImage>
 typename ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
 ::OutputImageType *
 ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
-::GetSumVarianceOutput()
+::GetDissimilarityOutput()
 {
   if (this->GetNumberOfOutputs() < 3)
     {
@@ -100,7 +100,7 @@ template <class TInputImage, class TOutputImage>
 typename ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
 ::OutputImageType *
 ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
-::GetSumEntropyOutput()
+::GetSumAverageOutput()
 {
   if (this->GetNumberOfOutputs() < 4)
     {
@@ -113,7 +113,7 @@ template <class TInputImage, class TOutputImage>
 typename ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
 ::OutputImageType *
 ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
-::GetDifferenceEntropyOutput()
+::GetSumVarianceOutput()
 {
   if (this->GetNumberOfOutputs() < 5)
     {
@@ -126,7 +126,7 @@ template <class TInputImage, class TOutputImage>
 typename ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
 ::OutputImageType *
 ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
-::GetDifferenceVarianceOutput()
+::GetSumEntropyOutput()
 {
   if (this->GetNumberOfOutputs() < 6)
     {
@@ -139,7 +139,7 @@ template <class TInputImage, class TOutputImage>
 typename ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
 ::OutputImageType *
 ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
-::GetIC1Output()
+::GetDifferenceEntropyOutput()
 {
   if (this->GetNumberOfOutputs() < 7)
     {
@@ -152,7 +152,7 @@ template <class TInputImage, class TOutputImage>
 typename ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
 ::OutputImageType *
 ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
-::GetIC2Output()
+::GetDifferenceVarianceOutput()
 {
   if (this->GetNumberOfOutputs() < 8)
     {
@@ -165,15 +165,28 @@ template <class TInputImage, class TOutputImage>
 typename ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
 ::OutputImageType *
 ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
-::GetMeanOutput()
+::GetIC1Output()
 {
-  if (this->GetNumberOfOutputs() < 1)
+  if (this->GetNumberOfOutputs() < 9)
     {
     return 0;
     }
   return static_cast<OutputImageType *>(this->GetOutput(8));
 }
 
+template <class TInputImage, class TOutputImage>
+typename ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
+::OutputImageType *
+ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
+::GetIC2Output()
+{
+  if (this->GetNumberOfOutputs() < 10)
+    {
+    return 0;
+    }
+  return static_cast<OutputImageType *>(this->GetOutput(9));
+}
+
 template <class TInputImage, class TOutputImage>
 void
 ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
@@ -275,6 +288,7 @@ ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
   InputImagePointerType  inputPtr              =      const_cast<InputImageType *>(this->GetInput());
   OutputImagePointerType meanPtr               =      this->GetMeanOutput();
   OutputImagePointerType variancePtr           =      this->GetVarianceOutput();
+  OutputImagePointerType dissimilarityPtr      =      this->GetDissimilarityOutput();
   OutputImagePointerType sumAveragePtr         =      this->GetSumAverageOutput();
   OutputImagePointerType sumVariancePtr        =      this->GetSumVarianceOutput();
   OutputImagePointerType sumEntropytPtr        =      this->GetSumEntropyOutput();
@@ -286,6 +300,7 @@ ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
   // Build output iterators
   itk::ImageRegionIteratorWithIndex<OutputImageType> varianceIt(variancePtr, outputRegionForThread);
   itk::ImageRegionIterator<OutputImageType>          meanIt(meanPtr, outputRegionForThread);
+  itk::ImageRegionIterator<OutputImageType>          dissimilarityIt(dissimilarityPtr, outputRegionForThread);
   itk::ImageRegionIterator<OutputImageType>          sumAverageIt(sumAveragePtr, outputRegionForThread);
   itk::ImageRegionIterator<OutputImageType>          sumVarianceIt(sumVariancePtr, outputRegionForThread);
   itk::ImageRegionIterator<OutputImageType>          sumEntropytIt(sumEntropytPtr, outputRegionForThread);
@@ -297,6 +312,7 @@ ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
   // Go to begin
   varianceIt.GoToBegin();
   meanIt.GoToBegin();
+  dissimilarityIt.GoToBegin();
   sumAverageIt.GoToBegin();
   sumVarianceIt.GoToBegin();
   sumEntropytIt.GoToBegin();
@@ -313,6 +329,7 @@ ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
   // Iterate on outputs to compute textures
   while (!varianceIt.IsAtEnd()
          && !meanIt.IsAtEnd()
+         && !dissimilarityIt.IsAtEnd()
          && !sumAverageIt.IsAtEnd()
          && !sumVarianceIt.IsAtEnd()
          && !sumEntropytIt.IsAtEnd()
@@ -350,6 +367,8 @@ ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
     GreyLevelCooccurrenceListType glcList;
     TotalAbsoluteFrequencyType totalFrequency = 0;
 
+    const double minSizeHist = std::min (m_Histogram->GetSize()[0], m_Histogram->GetSize()[1]);
+
     neighborIt = NeighborhoodIteratorType(m_NeighborhoodRadius, inputPtr, inputRegion);
     for ( neighborIt.GoToBegin(); !neighborIt.IsAtEnd(); ++neighborIt )
     {
@@ -423,14 +442,15 @@ ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
       totalFrequency = totalFrequency + 2;
     }
 
-    MeasurementType m_Mean                                        = itk::NumericTraits< MeasurementType >::Zero;
-    MeasurementType m_Variance                                    = itk::NumericTraits< MeasurementType >::Zero;
-    MeasurementType m_SumAverage                                  = itk::NumericTraits< MeasurementType >::Zero;
-    MeasurementType m_SumEntropy                                  = itk::NumericTraits< MeasurementType >::Zero;
-    MeasurementType m_SumVariance                                 = itk::NumericTraits< MeasurementType >::Zero;
+    MeasurementType m_Mean            				= itk::NumericTraits< MeasurementType >::Zero;
+    MeasurementType m_Variance        				= itk::NumericTraits< MeasurementType >::Zero;
+    MeasurementType m_Dissimilarity    				= itk::NumericTraits< MeasurementType >::Zero;
+    MeasurementType m_SumAverage      				= itk::NumericTraits< MeasurementType >::Zero;
+    MeasurementType m_SumEntropy      				= itk::NumericTraits< MeasurementType >::Zero;
+    MeasurementType m_SumVariance     				= itk::NumericTraits< MeasurementType >::Zero;
     MeasurementType m_DifferenceEntropy       = itk::NumericTraits< MeasurementType >::Zero;
     MeasurementType m_DifferenceVariance      = itk::NumericTraits< MeasurementType >::Zero;
-    MeasurementType m_IC1                           = itk::NumericTraits< MeasurementType >::Zero;
+    MeasurementType m_IC1                     = itk::NumericTraits< MeasurementType >::Zero;
     MeasurementType m_IC2                     = itk::NumericTraits< MeasurementType >::Zero;
 
     double Entropy = 0;
@@ -488,8 +508,6 @@ ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
       ++vectorIt;
       }
 
-    double minSizeHist = std::min (m_Histogram->GetSize()[0], m_Histogram->GetSize()[1]);
-
     //second pass over normalized co-occurrence list to find variance and pipj.
     //pipj is needed to calculate f11
     vectorIt = glcList.begin();
@@ -539,14 +557,23 @@ ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
 
     /* pipj computed below is totally different from earlier one which was used
      * to compute hxy1. This need to force an iterator over entire histogram.
-       Processing time is propotional to the m_NumberOfBins */
+       Processing time is propotional to the histogram bin size */
     double hxy2 = 0;
-              for(long unsigned int i = 0; i < histSize; ++i)
+		for(int i = 0; i < histSize; ++i)
       {
-      for(long unsigned int j = 0; j < histSize; ++j)
+      for(int j = 0; j < histSize; ++j)
         {
         double pipj = hx[j] * hy[i];
         hxy2 -= (pipj > 0.0001) ? pipj * vcl_log(pipj) : 0.;
+        double frequency = 0;
+        int instanceId = i * m_Histogram->GetSize()[0] + j;
+        if (instanceId < lookupArray.size())
+          {
+          int findex = lookupArray[instanceId];
+          if(findex > -1 )
+            frequency = glcList[findex].frequency;
+          }
+        m_Dissimilarity+= ( j - i ) * (frequency * frequency);
         }
       }
 
@@ -556,8 +583,9 @@ ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
     m_IC2 = (m_IC2 >= 0) ? vcl_sqrt (m_IC2) : 0;
 
     // Fill outputs
-    varianceIt.Set(m_Variance);
     meanIt.Set(m_Mean);
+    varianceIt.Set(m_Variance);
+    dissimilarityIt.Set(m_Dissimilarity);
     sumAverageIt.Set(m_SumAverage);
     sumVarianceIt.Set(m_SumVariance);
     sumEntropytIt.Set(m_SumEntropy);
@@ -572,6 +600,7 @@ ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
     // Increment iterators
     ++varianceIt;
     ++meanIt;
+    ++dissimilarityIt;
     ++sumAverageIt;
     ++sumVarianceIt;
     ++sumEntropytIt;
-- 
GitLab