ENH: updated texture computation based on new ScalarImageToTexturesFilter

Code/FeatureExtraction/otbHaralickTexturesImageFunction.h

@@ -18,10 +18,11 @@
 #ifndef __otbHaralickTexturesImageFunction_h
 #define __otbHaralickTexturesImageFunction_h
 
+#include "otbGreyLevelCooccurrenceIndexedList.h"
+#include "itkImageToImageFilter.h"
 #include "itkImageFunction.h"
 #include "itkFixedArray.h"
-#include "itkScalarImageToCooccurrenceMatrixFilter.h"
-#include "itkHistogramToTextureFeaturesFilter.h"
+#include "itkHistogram.h"
 
 namespace otb
 {
@@ -99,27 +100,23 @@ public:
   typedef typename InputImageType::Pointer         InputImagePointerType;
   typedef typename InputImageType::PixelType       InputPixelType;
   typedef typename InputImageType::RegionType      InputRegionType;
+  typedef typename InputImageType::OffsetType   OffsetType;
   typedef typename InputRegionType::SizeType       SizeType;
 
-  // Textures tools
-  /** Co-occurence matrix and textures calculator */
-    typedef itk::Statistics::ScalarImageToCooccurrenceMatrixFilter
-    <InputImageType>                               CoocurrenceMatrixGeneratorType;
-    typedef typename CoocurrenceMatrixGeneratorType
-        ::Pointer                                  CoocurrenceMatrixGeneratorPointerType;
-    typedef typename CoocurrenceMatrixGeneratorType
-        ::OffsetType                               OffsetType;
-    typedef typename CoocurrenceMatrixGeneratorType
-        ::HistogramType                            HistogramType;
-    typedef itk::Statistics::HistogramToTextureFeaturesFilter
-    <HistogramType>                                TextureCoefficientsCalculatorType;
-    typedef typename TextureCoefficientsCalculatorType
-        ::Pointer                                  TextureCoefficientsCalculatorPointerType;
-
   // Output typedef support
   typedef typename Superclass::OutputType          OutputType;
   typedef typename OutputType::ValueType           ScalarRealType;
-
+  // Textures tools
+  typedef GreyLevelCooccurrenceIndexedList< InputPixelType >   CooccurrenceIndexedListType;
+  typedef typename CooccurrenceIndexedListType::Pointer       CooccurrenceIndexedListPointerType;
+  typedef typename CooccurrenceIndexedListType::ConstPointer  CooccurrenceIndexedListConstPointerType;
+  typedef typename CooccurrenceIndexedListType::IndexType              CooccurrenceIndexType;
+  typedef typename CooccurrenceIndexedListType::PixelValueType         PixelValueType;
+  typedef typename CooccurrenceIndexedListType::RelativeFrequencyType  RelativeFrequencyType;
+  typedef typename CooccurrenceIndexedListType::VectorType             VectorType;
+
+  typedef typename VectorType::iterator                    VectorIteratorType;
+  typedef typename VectorType::const_iterator              VectorConstIteratorType;
   // Coord rep
   typedef TCoordRep                                CoordRepType;
 
@@ -154,7 +151,7 @@ public:
      * Get the radius of the neighborhood over which the
      * statistics are evaluated
      */
-    itkGetConstReferenceMacro( NeighborhoodRadius, unsigned int );
+  itkGetConstReferenceMacro( NeighborhoodRadius, unsigned int );
 
   /** Set the offset for co-occurence computation */
   itkSetMacro(Offset, OffsetType);
@@ -203,6 +200,9 @@ private:
 
   /** Input image maximum */
   InputPixelType m_InputImageMaximum;
+
+  static const double m_PixelValueTolerance = 0.0001;
+
 };
 
 } // namespace otb
@@ -212,4 +212,3 @@ private:
 #endif
 
 #endif
-

Code/FeatureExtraction/otbHaralickTexturesImageFunction.txx

@@ -38,7 +38,7 @@ HaralickTexturesImageFunction<TInputImage, TCoordRep>
  m_Offset(),
  m_NumberOfBinsPerAxis(8),
  m_InputImageMinimum(0),
- m_InputImageMaximum(256)
+ m_InputImageMaximum(255)
 {}
 
 template <class TInputImage, class TCoordRep>
@@ -77,24 +77,18 @@ HaralickTexturesImageFunction<TInputImage, TCoordRep>
     return textures;
     }
 
-
+  const double log2 = vcl_log(2.0);
   // Retrieve the input pointer
   InputImagePointerType inputPtr = const_cast<InputImageType *> (this->GetInputImage());
 
   // Compute the region on which co-occurence will be estimated
-  typename InputRegionType::IndexType inputIndex, inputIndexWithTwiceOffset;
-  typename InputRegionType::SizeType inputSize, inputSizeWithTwiceOffset;
+  typename InputRegionType::IndexType inputIndex;
+  typename InputRegionType::SizeType inputSize;
 
-  for(unsigned int dim = 0; dim<InputImageType::ImageDimension; ++dim)
+  for (unsigned int dim = 0; dim < InputImageType::ImageDimension; ++dim)
     {
-    inputIndex[dim] = std::min(
-                               (index[dim] - m_NeighborhoodRadius),
-                               (index[dim] - m_NeighborhoodRadius + m_Offset[dim])
-                              );
-    inputSize[dim] = 2 * m_NeighborhoodRadius + 1 + std::abs(m_Offset[dim]);
-
-    inputIndexWithTwiceOffset[dim] = index[dim] - m_NeighborhoodRadius - std::abs(m_Offset[dim]);
-    inputSizeWithTwiceOffset[dim] = inputSize[dim] + std::abs(m_Offset[dim]);
+    inputIndex[dim] = index[dim] - m_NeighborhoodRadius;
+    inputSize[dim] = 2 * m_NeighborhoodRadius + 1;
     }
 
   // Build the input  region
@@ -103,64 +97,141 @@ HaralickTexturesImageFunction<TInputImage, TCoordRep>
   inputRegion.SetSize(inputSize);
   inputRegion.Crop(inputPtr->GetRequestedRegion());
 
-  InputRegionType inputRegionWithTwiceOffset;
-  inputRegionWithTwiceOffset.SetIndex(inputIndexWithTwiceOffset);
-  inputRegionWithTwiceOffset.SetSize(inputSizeWithTwiceOffset);
-  inputRegionWithTwiceOffset.Crop(inputPtr->GetRequestedRegion());
-
-
-  /*********************************************************************************/
-  //Local copy of the input image around the processed pixel index
-  InputImagePointerType localInputImage = InputImageType::New();
-  localInputImage->SetRegions(inputRegionWithTwiceOffset);
-  localInputImage->Allocate();
-  typedef itk::ImageRegionIteratorWithIndex<InputImageType> ImageRegionIteratorType;
-  ImageRegionIteratorType itInputPtr(inputPtr, inputRegionWithTwiceOffset);
-  ImageRegionIteratorType itLocalInputImage(localInputImage, inputRegionWithTwiceOffset);
-  for (itInputPtr.GoToBegin(), itLocalInputImage.GoToBegin(); !itInputPtr.IsAtEnd(); ++itInputPtr, ++itLocalInputImage)
+  CooccurrenceIndexedListPointerType GLCIList = CooccurrenceIndexedListType::New();
+  GLCIList->Initialize(m_NumberOfBinsPerAxis, m_InputImageMinimum, m_InputImageMaximum);
+
+  // Next, find the minimum radius that encloses all the offsets.
+  unsigned int minRadius = 0;
+  for ( unsigned int i = 0; i < m_Offset.GetOffsetDimension(); i++ )
+    {
+    unsigned int distance = vcl_abs(m_Offset[i]);
+    if ( distance > minRadius )
+      {
+      minRadius = distance;
+      }
+    }
+  SizeType radius;
+  radius.Fill(minRadius);
+
+
+  typedef itk::ConstNeighborhoodIterator< InputImageType > NeighborhoodIteratorType;
+  NeighborhoodIteratorType neighborIt;
+  neighborIt = NeighborhoodIteratorType(radius, inputPtr, inputRegion);
+  for ( neighborIt.GoToBegin(); !neighborIt.IsAtEnd(); ++neighborIt )
     {
-    itLocalInputImage.Set(itInputPtr.Get());
+    const InputPixelType centerPixelIntensity = neighborIt.GetCenterPixel();
+    bool pixelInBounds;
+    const InputPixelType pixelIntensity =  neighborIt.GetPixel(m_Offset, pixelInBounds);
+    if ( !pixelInBounds )
+      {
+      continue; // don't put a pixel in the histogram if it's out-of-bounds.
+      }
+    GLCIList->AddPixelPair(centerPixelIntensity, pixelIntensity);
     }
-  /*********************************************************************************/
 
-  // Build the maskImage corresponding to inputRegion included in inputRegionWithTwiceOffset
-  InputImagePointerType maskImage = InputImageType::New();
-  maskImage->SetRegions(inputRegionWithTwiceOffset);
-  maskImage->Allocate();
-  maskImage->FillBuffer(0);
 
-  ImageRegionIteratorType itMask(maskImage, inputRegion);
-  for (itMask.GoToBegin(); !itMask.IsAtEnd(); ++itMask)
+  double pixelMean = 0.;
+  double marginalMean;
+  double marginalDevSquared = 0.;
+  double pixelVariance = 0.;
+
+  //Create and Intialize marginalSums
+  std::vector<double> marginalSums(m_NumberOfBinsPerAxis, 0);
+
+  //get co-occurrence vector and totalfrequency
+  VectorType glcVector = GLCIList->GetVector();
+  double totalFrequency = static_cast<double> (GLCIList->GetTotalFrequency());
+
+  //Normalize the co-occurrence indexed list and compute mean, marginalSum
+  typename VectorType::iterator it = glcVector.begin();
+  while( it != glcVector.end())
+    {
+    double frequency = (*it).second / totalFrequency;
+    CooccurrenceIndexType index = (*it).first;
+    pixelMean += index[0] * frequency;
+    marginalSums[index[0]] += frequency;
+    ++it;
+    }
+
+  /* Now get the mean and deviaton of the marginal sums.
+     Compute incremental mean and SD, a la Knuth, "The  Art of Computer
+     Programming, Volume 2: Seminumerical Algorithms",  section 4.2.2.
+     Compute mean and standard deviation using the recurrence relation:
+     M(1) = x(1), M(k) = M(k-1) + (x(k) - M(k-1) ) / k
+     S(1) = 0, S(k) = S(k-1) + (x(k) - M(k-1)) * (x(k) - M(k))
+     for 2 <= k <= n, then
+     sigma = vcl_sqrt(S(n) / n) (or divide by n-1 for sample SD instead of
+     population SD).
+  */
+
+  std::vector<double>::const_iterator msIt = marginalSums.begin();
+  marginalMean = *msIt;
+
+  //Increment iterator to start with index 1
+  ++msIt;
+  for(int k= 2; msIt != marginalSums.end(); ++k, ++msIt)
+    {
+    double M_k_minus_1 = marginalMean;
+    double S_k_minus_1 = marginalDevSquared;
+    double x_k = *msIt;
+    double M_k = M_k_minus_1 + ( x_k - M_k_minus_1 ) / k;
+    double S_k = S_k_minus_1 + ( x_k - M_k_minus_1 ) * ( x_k - M_k );
+    marginalMean = M_k;
+    marginalDevSquared = S_k;
+    }
+
+  marginalDevSquared = marginalDevSquared / m_NumberOfBinsPerAxis;
+
+  VectorConstIteratorType constVectorIt;
+  constVectorIt = glcVector.begin();
+  while( constVectorIt != glcVector.end())
+    {
+    RelativeFrequencyType frequency = (*constVectorIt).second / totalFrequency;
+    CooccurrenceIndexType        index = (*constVectorIt).first;
+    pixelVariance += ( index[0] - pixelMean ) * ( index[0] - pixelMean ) * frequency;
+    ++constVectorIt;
+    }
+
+  double pixelVarianceSquared = pixelVariance * pixelVariance;
+
+  // Variance is only used in correlation. If variance is 0, then (index[0] - pixelMean) * (index[1] - pixelMean)
+  // should be zero as well. In this case, set the variance to 1. in order to
+  // avoid NaN correlation.
+  if(pixelVarianceSquared < m_PixelValueTolerance)
+    {
+    pixelVarianceSquared = 1.;
+    }
+
+  //Initalize texture variables;
+  PixelValueType energy      = itk::NumericTraits< PixelValueType >::Zero;
+  PixelValueType entropy     = itk::NumericTraits< PixelValueType >::Zero;
+  PixelValueType correlation = itk::NumericTraits< PixelValueType >::Zero;
+  PixelValueType inverseDifferenceMoment      = itk::NumericTraits< PixelValueType >::Zero;
+  PixelValueType inertia             = itk::NumericTraits< PixelValueType >::Zero;
+  PixelValueType clusterShade        = itk::NumericTraits< PixelValueType >::Zero;
+  PixelValueType clusterProminence   = itk::NumericTraits< PixelValueType >::Zero;
+  PixelValueType haralickCorrelation = itk::NumericTraits< PixelValueType >::Zero;
+
+  //Compute textures
+  constVectorIt = glcVector.begin();
+  while( constVectorIt != glcVector.end())
     {
-    itMask.Set(1);
+    CooccurrenceIndexType index = (*constVectorIt).first;
+    RelativeFrequencyType frequency = (*constVectorIt).second / totalFrequency;
+    textures[0] += frequency * frequency;
+    textures[1] -= ( frequency > m_PixelValueTolerance ) ? frequency *vcl_log(frequency) / log2:0;
+    textures[2] += ( ( index[0] - pixelMean ) * ( index[1] - pixelMean ) * frequency ) / pixelVarianceSquared;
+    textures[3] += frequency / ( 1.0 + ( index[0] - index[1] ) * ( index[0] - index[1] ) );
+    textures[4] += ( index[0] - index[1] ) * ( index[0] - index[1] ) * frequency;
+    textures[5] += vcl_pow( ( index[0] - pixelMean ) + ( index[1] - pixelMean ), 3 ) * frequency;
+    textures[6] += vcl_pow( ( index[0] - pixelMean ) + ( index[1] - pixelMean ), 4 ) * frequency;
+    textures[7] += index[0] * index[1] * frequency;
+    ++constVectorIt;
     }
 
-  // Build the co-occurence matrix generator
-  CoocurrenceMatrixGeneratorPointerType coOccurenceMatrixGenerator = CoocurrenceMatrixGeneratorType::New();
-  coOccurenceMatrixGenerator->SetInput(localInputImage);
-  coOccurenceMatrixGenerator->SetOffset(m_Offset);
-  coOccurenceMatrixGenerator->SetNumberOfBinsPerAxis(m_NumberOfBinsPerAxis);
-  coOccurenceMatrixGenerator->SetPixelValueMinMax(m_InputImageMinimum, m_InputImageMaximum);
-  coOccurenceMatrixGenerator->SetMaskImage(maskImage);
-  coOccurenceMatrixGenerator->SetInsidePixelValue(1);
-  coOccurenceMatrixGenerator->Update();
-
-  // Build the texture calculator
-  TextureCoefficientsCalculatorPointerType texturesCalculator = TextureCoefficientsCalculatorType::New();
-
-  // Compute textures indices
-  texturesCalculator->SetInput(coOccurenceMatrixGenerator->GetOutput());
-  texturesCalculator->Update();
-
-  // Fill the output vector
-  textures[0]=texturesCalculator->GetEnergy();
-  textures[1]=texturesCalculator->GetEntropy();
-  textures[2]=texturesCalculator->GetCorrelation();
-  textures[3]=texturesCalculator->GetInverseDifferenceMoment();
-  textures[4]=texturesCalculator->GetInertia();
-  textures[5]=texturesCalculator->GetClusterShade();
-  textures[6]=texturesCalculator->GetClusterProminence();
-  textures[7]=texturesCalculator->GetHaralickCorrelation();
+  //Use epsilon check here?
+  textures[7] = (marginalDevSquared != 0) ? ( textures[7] - marginalMean * marginalMean )  / marginalDevSquared : 0;
+  // haralickCorrelation = ( haralickCorrelation - marginalMean * marginalMean )  / marginalDevSquared;
 
   // Return result
   return textures;

Code/ObjectDetection/otbHaralickTexturesIFFactory.txx

@@ -43,9 +43,9 @@ HaralickTexturesIFFactory<TImageType, TCoordRep, TPrecision>
 
   function->SetInputImage(image);
   function->GetInternalImageFunction()->SetNeighborhoodRadius(param[0]);
-  function->GetInternalImageFunction()->SetInputImageMinimum(param[1]);
-  function->GetInternalImageFunction()->SetInputImageMaximum(param[2]);
-  function->GetInternalImageFunction()->SetNumberOfBinsPerAxis(param[3]);
+  function->SetInputImageMinimum(param[1]);
+  function->SetInputImageMaximum(param[2]);
+  function->SetNumberOfBinsPerAxis(param[3]);
 
   OffsetType offset;
   offset.Fill(param[4]);

Testing/Code/FeatureExtraction/otbHaralickTexturesImageFunction.cxx

@@ -51,6 +51,9 @@ int otbHaralickTexturesImageFunction(int argc, char * argv[])
   HaralickTexturesImageFunctionType::Pointer haralick = HaralickTexturesImageFunctionType::New();
   haralick->SetInputImage(reader->GetOutput());
   haralick->SetNeighborhoodRadius(10);
+  haralick->SetNumberOfBinsPerAxis(8);
+  haralick->SetInputImageMinimum(0);
+  haralick->SetInputImageMaximum(255);
 
   HaralickTexturesImageFunctionType::OffsetType offset;
   offset.Fill(1);