From c1b964b56f8d21c8e8ce0bd4ce53016de1a6137e Mon Sep 17 00:00:00 2001
From: Manuel Grizonnet <manuel.grizonnet@gmail.com>
Date: Fri, 25 Jun 2010 16:14:51 +0200
Subject: [PATCH] ENH:add MEAN texture in advanced texture coefficients
 computation

---
 ...rixAdvancedTextureCoefficientsCalculator.h |  7 +++--
 ...xAdvancedTextureCoefficientsCalculator.txx | 11 ++++----
 .../otbScalarImageToAdvancedTexturesFilter.h  |  3 +++
 ...otbScalarImageToAdvancedTexturesFilter.txx | 26 ++++++++++++++++---
 Testing/Code/FeatureExtraction/CMakeLists.txt |  4 ++-
 ...xAdvancedTextureCoefficientsCalculator.cxx | 12 ++++++++-
 ...otbScalarImageToAdvancedTexturesFilter.cxx |  6 +++++
 7 files changed, 56 insertions(+), 13 deletions(-)

diff --git a/Code/FeatureExtraction/otbGreyLevelCooccurrenceMatrixAdvancedTextureCoefficientsCalculator.h b/Code/FeatureExtraction/otbGreyLevelCooccurrenceMatrixAdvancedTextureCoefficientsCalculator.h
index 5b313456a2..b40d274d40 100644
--- a/Code/FeatureExtraction/otbGreyLevelCooccurrenceMatrixAdvancedTextureCoefficientsCalculator.h
+++ b/Code/FeatureExtraction/otbGreyLevelCooccurrenceMatrixAdvancedTextureCoefficientsCalculator.h
@@ -34,6 +34,8 @@ namespace otb {
  * The features calculated are as follows (where \f$ g(i, j) \f$ is the element in
  * cell i, j of a a normalized GLCM):
  *
+ * "Mean" \f$ = \sum_{i,j}i g(i,j) \f$
+ *
  * "Sum of squares: Variance" \f$ = f_4 = \sum_{i,j}(i -mu)^2 g(i,j) \f$
  *
  * "Sum average" \f$ = f_6 = -\sum_{i}i g_{x+y}(i)  
@@ -117,6 +119,7 @@ public:
   itkSetObjectMacro( Histogram, HistogramType );
   itkGetObjectMacro( Histogram, HistogramType );
   
+  itkGetMacro(Mean, double);
   itkGetMacro(Variance, double);
   itkGetMacro(SumAverage, double);
   itkGetMacro(SumVariance, double);
@@ -142,8 +145,8 @@ private:
   HistogramPointer m_Histogram;
   
   void NormalizeHistogram(void);
-  void ComputeMean( double &pixelMean );
-  double m_Variance, m_SumAverage, m_SumVariance,m_SumEntropy, m_DifferenceEntropy,
+  void ComputeMean();
+  double m_Mean, m_Variance, m_SumAverage, m_SumVariance,m_SumEntropy, m_DifferenceEntropy,
   m_DifferenceVariance, m_IC1, m_IC2;
 
   double ComputePS ( long unsigned int k );
diff --git a/Code/FeatureExtraction/otbGreyLevelCooccurrenceMatrixAdvancedTextureCoefficientsCalculator.txx b/Code/FeatureExtraction/otbGreyLevelCooccurrenceMatrixAdvancedTextureCoefficientsCalculator.txx
index 22bf0c22c4..c75af60b05 100644
--- a/Code/FeatureExtraction/otbGreyLevelCooccurrenceMatrixAdvancedTextureCoefficientsCalculator.txx
+++ b/Code/FeatureExtraction/otbGreyLevelCooccurrenceMatrixAdvancedTextureCoefficientsCalculator.txx
@@ -41,7 +41,7 @@ NormalizeHistogram( void )
 template< class THistogram >
 void
 GreyLevelCooccurrenceMatrixAdvancedTextureCoefficientsCalculator< THistogram >::
-ComputeMean( double &pixelMean )
+ComputeMean( )
 {
   // This function takes two passes through the histogram and two passes through
   // an array of the same length as a histogram axis. This could probably be
@@ -49,7 +49,7 @@ ComputeMean( double &pixelMean )
   typedef typename HistogramType::Iterator HistogramIterator;
       
   // Initialize everything
-  pixelMean = 0;
+  m_Mean = 0;
       
   // Ok, now do the first pass through the histogram to get the marginal sums
   // and compute the pixel mean
@@ -58,7 +58,7 @@ ComputeMean( double &pixelMean )
     {
     MeasurementType frequency = hit.GetFrequency();
     IndexType index = m_Histogram->GetIndex(hit.GetInstanceIdentifier());
-    pixelMean += index[0] * frequency;
+    m_Mean += index[0] * frequency;
     }
 } 
 
@@ -80,8 +80,7 @@ Compute ()
     }
       
   // Now get the pixel mean.
-  double pixelMean;
-  this->ComputeMean( pixelMean );
+  this->ComputeMean();
                                                                 
   m_SumAverage = 0;
   m_SumEntropy = 0;
@@ -151,7 +150,7 @@ Compute ()
     
     IndexType index = m_Histogram->GetIndex(hit.GetInstanceIdentifier());
     
-    m_Variance += ( (index[0] - pixelMean) * (index[0] - pixelMean) ) * frequency;
+    m_Variance += ( (index[0] - m_Mean) * (index[0] - m_Mean) ) * frequency;
     Entropy -= (frequency > 0.0001) ? frequency * vcl_log(frequency) / log2 : 0;
     
     double pipj = m_Histogram->GetFrequency (index[0] , 0) * m_Histogram->GetFrequency (index[1] , 1);
diff --git a/Code/FeatureExtraction/otbScalarImageToAdvancedTexturesFilter.h b/Code/FeatureExtraction/otbScalarImageToAdvancedTexturesFilter.h
index 46ad156c82..2f54b1bec8 100644
--- a/Code/FeatureExtraction/otbScalarImageToAdvancedTexturesFilter.h
+++ b/Code/FeatureExtraction/otbScalarImageToAdvancedTexturesFilter.h
@@ -91,6 +91,9 @@ public:
   /** Get the input image maximum */
   itkGetMacro(InputImageMaximum,InputPixelType);
 
+  /** Get the mean output image */
+  OutputImageType * GetMeanOutput();
+  
   /** Get the variance output image */
   OutputImageType * GetVarianceOutput();
 
diff --git a/Code/FeatureExtraction/otbScalarImageToAdvancedTexturesFilter.txx b/Code/FeatureExtraction/otbScalarImageToAdvancedTexturesFilter.txx
index 0849ce9b13..f58cbed03d 100644
--- a/Code/FeatureExtraction/otbScalarImageToAdvancedTexturesFilter.txx
+++ b/Code/FeatureExtraction/otbScalarImageToAdvancedTexturesFilter.txx
@@ -33,10 +33,10 @@ ScalarImageToAdvancedTexturesFilter<TInputImage,TOutputImage>
                                   m_InputImageMinimum(0),
                                   m_InputImageMaximum(256)
 {
-  // There are 8 outputs corresponding to the 8 textures indices
-  this->SetNumberOfOutputs(8);
+  // There are 9 outputs corresponding to the 8 textures indices
+  this->SetNumberOfOutputs(9);
 
-  // Create the 8 outputs
+  // Create the 9 outputs
   this->SetNthOutput(0,OutputImageType::New());
   this->SetNthOutput(1,OutputImageType::New());
   this->SetNthOutput(2,OutputImageType::New());
@@ -45,6 +45,7 @@ ScalarImageToAdvancedTexturesFilter<TInputImage,TOutputImage>
   this->SetNthOutput(5,OutputImageType::New());
   this->SetNthOutput(6,OutputImageType::New());
   this->SetNthOutput(7,OutputImageType::New());
+  this->SetNthOutput(8,OutputImageType::New());
 }
 
 template <class TInputImage,class TOutputImage>
@@ -156,6 +157,19 @@ ScalarImageToAdvancedTexturesFilter<TInputImage,TOutputImage>
   return static_cast<OutputImageType *>(this->GetOutput(7));
 }
 
+template <class TInputImage,class TOutputImage>
+typename ScalarImageToAdvancedTexturesFilter<TInputImage,TOutputImage>
+::OutputImageType *
+ScalarImageToAdvancedTexturesFilter<TInputImage,TOutputImage>
+::GetMeanOutput()
+{
+  if(this->GetNumberOfOutputs()<1)
+    {
+    return 0;
+    }
+  return static_cast<OutputImageType *>(this->GetOutput(8));
+}
+
 template <class TInputImage,class TOutputImage>
 void
 ScalarImageToAdvancedTexturesFilter<TInputImage,TOutputImage>
@@ -221,6 +235,7 @@ ScalarImageToAdvancedTexturesFilter<TInputImage,TOutputImage>
 {
   // Retrieve the input and output pointers
   InputImagePointerType  inputPtr             =      const_cast<InputImageType *>(this->GetInput());
+  OutputImagePointerType meanPtr            =      this->GetMeanOutput();
   OutputImagePointerType variancePtr            =      this->GetVarianceOutput();
   OutputImagePointerType sumAveragePtr           =      this->GetSumAverageOutput();
   OutputImagePointerType sumVariancePtr       =      this->GetSumVarianceOutput();
@@ -232,6 +247,7 @@ ScalarImageToAdvancedTexturesFilter<TInputImage,TOutputImage>
 
   // Build output iterators
   itk::ImageRegionIteratorWithIndex<OutputImageType> varianceIt(variancePtr,outputRegionForThread);
+  itk::ImageRegionIterator<OutputImageType>          meanIt(meanPtr,outputRegionForThread);
   itk::ImageRegionIterator<OutputImageType>          sumAverageIt(sumAveragePtr,outputRegionForThread);
   itk::ImageRegionIterator<OutputImageType>          sumVarianceIt(sumVariancePtr,outputRegionForThread);
   itk::ImageRegionIterator<OutputImageType>          sumEntropytIt(sumEntropytPtr,outputRegionForThread);
@@ -242,6 +258,7 @@ ScalarImageToAdvancedTexturesFilter<TInputImage,TOutputImage>
 
   // Go to begin
   varianceIt.GoToBegin();
+  meanIt.GoToBegin();
   sumAverageIt.GoToBegin();
   sumVarianceIt.GoToBegin();
   sumEntropytIt.GoToBegin();
@@ -265,6 +282,7 @@ ScalarImageToAdvancedTexturesFilter<TInputImage,TOutputImage>
 
   // Iterate on outputs to compute textures
   while(!varianceIt.IsAtEnd()
+      &&!meanIt.IsAtEnd()
       &&!sumAverageIt.IsAtEnd()
       &&!sumVarianceIt.IsAtEnd()
       &&!sumEntropytIt.IsAtEnd()
@@ -302,6 +320,7 @@ ScalarImageToAdvancedTexturesFilter<TInputImage,TOutputImage>
 
     // Fill outputs
     varianceIt.Set(texturesCalculator->GetVariance());
+    meanIt.Set(texturesCalculator->GetMean());
     sumAverageIt.Set(texturesCalculator->GetSumAverage());
     sumVarianceIt.Set(texturesCalculator->GetSumVariance());
     sumEntropytIt.Set(texturesCalculator->GetSumEntropy());
@@ -315,6 +334,7 @@ ScalarImageToAdvancedTexturesFilter<TInputImage,TOutputImage>
 
     // Increment iterators
     ++varianceIt;
+    ++meanIt;
     ++sumAverageIt;
     ++sumVarianceIt;
     ++sumEntropytIt;
diff --git a/Testing/Code/FeatureExtraction/CMakeLists.txt b/Testing/Code/FeatureExtraction/CMakeLists.txt
index fb6b1a2336..07e398aee7 100644
--- a/Testing/Code/FeatureExtraction/CMakeLists.txt
+++ b/Testing/Code/FeatureExtraction/CMakeLists.txt
@@ -1817,9 +1817,11 @@ ADD_TEST(feTuScalarImageToAdvancedTexturesFilterNew ${FEATUREEXTRACTION_TESTS15}
 )
 
 ADD_TEST(feTvScalarImageToAdvancedTexturesFilter ${FEATUREEXTRACTION_TESTS15}
-        --compare-n-images ${EPSILON_8} 8   
+        --compare-n-images ${EPSILON_8} 9   
             ${BASELINE}/feTvScalarImageToAdvancedTexturesFilterOutputVariance.tif
             ${TEMP}/feTvScalarImageToAdvancedTexturesFilterOutputVariance.tif
+            ${BASELINE}/feTvScalarImageToAdvancedTexturesFilterOutputMean.tif
+            ${TEMP}/feTvScalarImageToAdvancedTexturesFilterOutputMean.tif
             ${BASELINE}/feTvScalarImageToAdvancedTexturesFilterOutputSumAverage.tif
             ${TEMP}/feTvScalarImageToAdvancedTexturesFilterOutputSumAverage.tif
             ${BASELINE}/feTvScalarImageToAdvancedTexturesFilterOutputSumVariance.tif
diff --git a/Testing/Code/FeatureExtraction/otbGreyLevelCooccurrenceMatrixAdvancedTextureCoefficientsCalculator.cxx b/Testing/Code/FeatureExtraction/otbGreyLevelCooccurrenceMatrixAdvancedTextureCoefficientsCalculator.cxx
index 14e101c86c..604953efaf 100644
--- a/Testing/Code/FeatureExtraction/otbGreyLevelCooccurrenceMatrixAdvancedTextureCoefficientsCalculator.cxx
+++ b/Testing/Code/FeatureExtraction/otbGreyLevelCooccurrenceMatrixAdvancedTextureCoefficientsCalculator.cxx
@@ -79,6 +79,7 @@ int otbGreyLevelCooccurrenceMatrixAdvancedTextureCoefficientsCalculator(int argc
   glcmCalc->Compute();
   
   double trueVariance = 7.6475;
+  double trueMean = 3.95;
   double trueSumAverage = 0;
   double trueSumVariance = 0;
   double trueSumEntropy = 0;
@@ -88,6 +89,7 @@ int otbGreyLevelCooccurrenceMatrixAdvancedTextureCoefficientsCalculator(int argc
   double trueIC2 = 0.908663;
   
   double variance = glcmCalc->GetVariance();
+  double mean = glcmCalc->GetMean();
   double sumAverage = glcmCalc->GetSumAverage();
   double sumVariance = glcmCalc->GetSumVariance();
   double sumEntropy = glcmCalc->GetSumEntropy();
@@ -105,7 +107,15 @@ int otbGreyLevelCooccurrenceMatrixAdvancedTextureCoefficientsCalculator(int argc
       << variance << std::endl;
     passed = false;
     }
-   
+  
+  if( vnl_math_abs(mean - trueMean) > 0.001 )
+    {
+    std::cerr << "Error:" << std::endl;
+    std::cerr << "Mean calculated wrong. Expected: " << trueMean << ", got: " 
+      << mean << std::endl;
+    passed = false;
+    }
+    
   if( vnl_math_abs(sumAverage - trueSumAverage) > 0.001 )
     {
     std::cerr << "Error:" << std::endl;
diff --git a/Testing/Code/FeatureExtraction/otbScalarImageToAdvancedTexturesFilter.cxx b/Testing/Code/FeatureExtraction/otbScalarImageToAdvancedTexturesFilter.cxx
index 369b6a3ae7..a435e0f25d 100644
--- a/Testing/Code/FeatureExtraction/otbScalarImageToAdvancedTexturesFilter.cxx
+++ b/Testing/Code/FeatureExtraction/otbScalarImageToAdvancedTexturesFilter.cxx
@@ -88,6 +88,12 @@ int otbScalarImageToAdvancedTexturesFilter(int argc, char * argv[])
   oss<<outprefix<<"Variance.tif";
   writer->SetInput(filter->GetVarianceOutput());
   writer->SetFileName(oss.str());
+  writer->Update();
+  
+   oss.str("");
+  oss<<outprefix<<"Mean.tif";
+  writer->SetInput(filter->GetMeanOutput());
+  writer->SetFileName(oss.str());
   writer->Update();
   
   oss.str("");
-- 
GitLab