From 619a1ae33081fd38e97aa1a541b1d54ce5874e0b Mon Sep 17 00:00:00 2001
From: Rashad Kanavath <mohammed.rashad-km@cnes.fr>
Date: Wed, 7 May 2014 14:56:26 +0200
Subject: [PATCH] ENH: Using GreyLevelCooccurrenceIndexedList in simple and
advanced haralick textures
---
.../otbGreyLevelCooccurrenceIndexedList.h | 3 +
.../otbGreyLevelCooccurrenceIndexedList.txx | 14 ++
.../otbScalarImageToAdvancedTexturesFilter.h | 41 ++--
.../otbScalarImageToTexturesFilter.h | 39 ++--
.../otbScalarImageToTexturesFilter.txx | 195 +++++++-----------
5 files changed, 115 insertions(+), 177 deletions(-)
diff --git a/Code/FeatureExtraction/otbGreyLevelCooccurrenceIndexedList.h b/Code/FeatureExtraction/otbGreyLevelCooccurrenceIndexedList.h
index bdffb96ff4..4028effe07 100644
--- a/Code/FeatureExtraction/otbGreyLevelCooccurrenceIndexedList.h
+++ b/Code/FeatureExtraction/otbGreyLevelCooccurrenceIndexedList.h
@@ -107,6 +107,9 @@ public:
/** Get std::vector containing non-zero co-occurrence pairs */
VectorType GetVector();
+ /** Normailze the GLCIL */
+ void Normalize();
+
/** Add the CooccurrencePairType to list. Next occurrence of same index is
* checked via m_LookupArray value and the correspoding frequency value is
* incremented. If m_Symmetry is true the co-occurrence pair is added again */
diff --git a/Code/FeatureExtraction/otbGreyLevelCooccurrenceIndexedList.txx b/Code/FeatureExtraction/otbGreyLevelCooccurrenceIndexedList.txx
index 1c4f11edc5..61df91f29a 100644
--- a/Code/FeatureExtraction/otbGreyLevelCooccurrenceIndexedList.txx
+++ b/Code/FeatureExtraction/otbGreyLevelCooccurrenceIndexedList.txx
@@ -43,6 +43,20 @@ GreyLevelCooccurrenceIndexedList<THistogram>
m_TotalFrequency = 0;
}
+template <class THistogram >
+void
+GreyLevelCooccurrenceIndexedList<THistogram>
+::Normalize()
+{
+ //Normalize the co-occurrence indexed list
+ typename VectorType::iterator it = m_Vector.begin();
+ while( it != m_Vector.end())
+ {
+ (*it).frequency = (*it).frequency / m_TotalFrequency;
+ ++it;
+ }
+}
+
template <class THistogram >
typename GreyLevelCooccurrenceIndexedList<THistogram>::RelativeFrequencyType
GreyLevelCooccurrenceIndexedList<THistogram>::
diff --git a/Code/FeatureExtraction/otbScalarImageToAdvancedTexturesFilter.h b/Code/FeatureExtraction/otbScalarImageToAdvancedTexturesFilter.h
index c7b6571e50..00cbdccbd4 100644
--- a/Code/FeatureExtraction/otbScalarImageToAdvancedTexturesFilter.h
+++ b/Code/FeatureExtraction/otbScalarImageToAdvancedTexturesFilter.h
@@ -18,10 +18,9 @@
#ifndef __otbScalarImageToAdvancedTexturesFilter_h
#define __otbScalarImageToAdvancedTexturesFilter_h
+#include "otbGreyLevelCooccurrenceIndexedList.h"
#include "itkImageToImageFilter.h"
#include "itkHistogram.h"
-#include "itkArray.h"
-#include <vector>
namespace otb
{
@@ -119,29 +118,16 @@ public:
typedef typename HistogramType::MeasurementVectorType MeasurementVectorType;
typedef typename HistogramType::RelativeFrequencyType RelativeFrequencyType;
typedef typename HistogramType::TotalAbsoluteFrequencyType TotalAbsoluteFrequencyType;
- typedef typename HistogramType::IndexType HistogramIndexType;
+ typedef typename HistogramType::SizeType HistogramSizeType;
+ typedef typename HistogramType::IndexType CooccurrenceIndexType;
- /** Lookup array used to store the index of the given pixel pair. This index
- * is calculated from the HistogramType */
- typedef itk::Array<int> LookupArrayType;
+ typedef GreyLevelCooccurrenceIndexedList< HistogramType > CooccurrenceIndexedListType;
+ typedef typename CooccurrenceIndexedListType::Pointer CooccurrenceIndexedListPointerType;
+ typedef typename CooccurrenceIndexedListType::ConstPointer CooccurrenceIndexedListConstPointerType;
- /** struct to hold the index of pixel pair and its frequency. The frequnecy is
- * incremented for the next occurrence of the same pair */
-
- typedef struct Cooccurrence {
- HistogramIndexType index;
- RelativeFrequencyType frequency;
- }CooccurrenceType;
-
- /** array to hold CooccurrenceType. Index of the array where the cooccurrence
- * is stored is saved in the LookupArrayType. */
- typedef std::vector<CooccurrenceType> GreyLevelCooccurrenceListType;
-
- /** std::vector iterator typedef for GreyLevelCooccurrenceListType */
- typedef typename std::vector<CooccurrenceType>::iterator GreyLevelCooccurrenceListIteratorType;
-
- /** std::vector const_iterator typedef for GreyLevelCooccurrenceListType */
- typedef typename std::vector<CooccurrenceType>::const_iterator GreyLevelCooccurrenceListConstIteratorType;
+ typedef typename CooccurrenceIndexedListType::VectorType VectorType;
+ typedef typename VectorType::iterator VectorIteratorType;
+ typedef typename VectorType::const_iterator VectorConstIteratorType;
itkStaticConstMacro(MeasurementVectorSize, int, 2 );
@@ -174,6 +160,8 @@ public:
/** Get the input image maximum */
itkGetMacro(InputImageMaximum, InputPixelType);
+ void SetBinsAndMinMax(unsigned int bins, InputPixelType min, InputPixelType max);
+
/** Get the mean output image */
OutputImageType * GetMeanOutput();
@@ -209,13 +197,8 @@ protected:
ScalarImageToAdvancedTexturesFilter();
/** Destructor */
~ScalarImageToAdvancedTexturesFilter();
-
- /** Before ThreadedGenerateData created a single histogram instance **/
- virtual void BeforeThreadedGenerateData( );
-
/** Generate the input requested region */
virtual void GenerateInputRequestedRegion();
-
/** Parallel textures extraction */
virtual void ThreadedGenerateData(const OutputRegionType& outputRegion, itk::ThreadIdType threadId);
@@ -247,6 +230,8 @@ private:
/** Histogram instance used to get index for the pixel pair */
HistogramPointer m_Histogram;
+ HistogramSizeType m_HistSize;
+
};
} // End namespace otb
diff --git a/Code/FeatureExtraction/otbScalarImageToTexturesFilter.h b/Code/FeatureExtraction/otbScalarImageToTexturesFilter.h
index fcc3de861a..6c3e445394 100644
--- a/Code/FeatureExtraction/otbScalarImageToTexturesFilter.h
+++ b/Code/FeatureExtraction/otbScalarImageToTexturesFilter.h
@@ -19,6 +19,7 @@ PURPOSE. See the above copyright notices for more information.
#ifndef __otbScalarImageToTexturesFilter_h
#define __otbScalarImageToTexturesFilter_h
+#include "otbGreyLevelCooccurrenceIndexedList.h"
#include "itkImageToImageFilter.h"
#include "itkHistogram.h"
@@ -130,30 +131,18 @@ public:
typedef typename HistogramType::MeasurementVectorType MeasurementVectorType;
typedef typename HistogramType::RelativeFrequencyType RelativeFrequencyType;
typedef typename HistogramType::TotalAbsoluteFrequencyType TotalAbsoluteFrequencyType;
- typedef typename HistogramType::IndexType HistogramIndexType;
+ typedef typename HistogramType::SizeType HistogramSizeType;
+ typedef typename HistogramType::IndexType CooccurrenceIndexType;
+ typedef GreyLevelCooccurrenceIndexedList< HistogramType > CooccurrenceIndexedListType;
+ typedef typename CooccurrenceIndexedListType::Pointer CooccurrenceIndexedListPointerType;
+ typedef typename CooccurrenceIndexedListType::ConstPointer CooccurrenceIndexedListConstPointerType;
- /** Lookup array used to store the index of the given pixel pair. This index
- * is calculated from the HistogramType */
- typedef itk::Array<int> LookupArrayType;
+ typedef typename CooccurrenceIndexedListType::VectorType VectorType;
+ typedef typename VectorType::iterator VectorIteratorType;
+ typedef typename VectorType::const_iterator VectorConstIteratorType;
- /** struct to hold the index of pixel pair and its frequency. The frequnecy is
- * incremented for the next occurrence of the same pair */
-
- typedef struct Cooccurrence {
- HistogramIndexType index;
- RelativeFrequencyType frequency;
- }CooccurrenceType;
-
- /** array to hold CooccurrenceType. Index of the array where the cooccurrence
- * is stored is saved in the LookupArrayType. */
- typedef std::vector<CooccurrenceType> GreyLevelCooccurrenceListType;
-
- /** std::vector iterator typedef for GreyLevelCooccurrenceListType */
- typedef typename std::vector<CooccurrenceType>::iterator GreyLevelCooccurrenceListIteratorType;
-
- /** std::vector const_iterator typedef for GreyLevelCooccurrenceListType */
- typedef typename std::vector<CooccurrenceType>::const_iterator GreyLevelCooccurrenceListConstIteratorType;
+ void SetBinsAndMinMax(unsigned int bins, InputPixelType min, InputPixelType max);
itkStaticConstMacro(MeasurementVectorSize, int, 2 );
@@ -217,10 +206,6 @@ protected:
~ScalarImageToTexturesFilter();
/** Generate the input requested region */
virtual void GenerateInputRequestedRegion();
-
- /** Before ThreadedGenerateData created a single histogram instance **/
- virtual void BeforeThreadedGenerateData( );
-
/** Parallel textures extraction */
virtual void ThreadedGenerateData(const OutputRegionType& outputRegion, itk::ThreadIdType threadId);
@@ -252,9 +237,11 @@ private:
/** Histogram instance used to get index for the pixel pair */
HistogramPointer m_Histogram;
+ HistogramSizeType m_HistSize;
+
/** This method is same as ComputeMeansAndVariance in
* itkHistogramToTexturesFilter. This has been modified for using new GLCIL */
- void ComputeMeansAndVariances(GreyLevelCooccurrenceListType& frequencyList, double & pixelMean,
+ void ComputeMeansAndVariances(VectorType& frequencyList, double & pixelMean,
double & marginalMean, double & marginalDevSquared, double & pixelVariance);
};
diff --git a/Code/FeatureExtraction/otbScalarImageToTexturesFilter.txx b/Code/FeatureExtraction/otbScalarImageToTexturesFilter.txx
index b1b21bcd74..499e605410 100644
--- a/Code/FeatureExtraction/otbScalarImageToTexturesFilter.txx
+++ b/Code/FeatureExtraction/otbScalarImageToTexturesFilter.txx
@@ -31,11 +31,13 @@ namespace otb
{
template <class TInputImage, class TOutputImage>
ScalarImageToTexturesFilter<TInputImage, TOutputImage>
-::ScalarImageToTexturesFilter() : m_Radius(),
- m_Offset(),
- m_NumberOfBinsPerAxis(8),
- m_InputImageMinimum(0),
- m_InputImageMaximum(255)
+::ScalarImageToTexturesFilter()
+: m_Radius()
+, m_Offset()
+, m_HistSize(2)
+, m_NumberOfBinsPerAxis(8)
+, m_InputImageMinimum(0)
+, m_InputImageMaximum(255)
{
// There are 8 outputs corresponding to the 8 textures indices
this->SetNumberOfRequiredOutputs(8);
@@ -160,6 +162,43 @@ ScalarImageToTexturesFilter<TInputImage, TOutputImage>
return static_cast<OutputImageType *>(this->GetOutput(7));
}
+
+template <class TInputImage, class TOutputImage>
+void
+ScalarImageToTexturesFilter<TInputImage, TOutputImage>
+::SetBinsAndMinMax(unsigned int numberOfBinsPerAxis,
+ InputPixelType inputImageMinimum,
+ InputPixelType inputImageMaximum)
+{
+ /** calulate minimum offset and set it as neigborhood radius **/
+ 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;
+ }
+ }
+ m_NeighborhoodRadius.Fill(minRadius);
+
+ /** Initalize m_Histogram with given min, max and number of bins **/
+ MeasurementVectorType lowerBound;
+ MeasurementVectorType upperBound;
+ m_InputImageMinimum = inputImageMinimum;
+ m_InputImageMaximum = inputImageMaximum;
+ m_NumberOfBinsPerAxis = numberOfBinsPerAxis;
+ m_Histogram = HistogramType::New();
+ m_Histogram->SetMeasurementVectorSize( MeasurementVectorSize );
+ lowerBound.SetSize( MeasurementVectorSize );
+ upperBound.SetSize( MeasurementVectorSize );
+ lowerBound.Fill(m_InputImageMinimum);
+ upperBound.Fill(m_InputImageMaximum+1);
+ m_HistSize.Fill(m_NumberOfBinsPerAxis);
+ m_Histogram->Initialize(m_HistSize, lowerBound, upperBound);
+}
+
+
template <class TInputImage, class TOutputImage>
void
ScalarImageToTexturesFilter<TInputImage, TOutputImage>
@@ -220,38 +259,6 @@ ScalarImageToTexturesFilter<TInputImage, TOutputImage>
}
}
-template <class TInputImage, class TOutputImage>
-void
-ScalarImageToTexturesFilter<TInputImage, TOutputImage>
-::BeforeThreadedGenerateData()
-{
- /** Initalize m_Histogram with given min, max and number of bins **/
- MeasurementVectorType m_LowerBound;
- MeasurementVectorType m_UpperBound;
- m_Histogram = HistogramType::New();
- m_Histogram->SetMeasurementVectorSize( MeasurementVectorSize );
- m_LowerBound.SetSize( MeasurementVectorSize );
- m_UpperBound.SetSize( MeasurementVectorSize );
- m_LowerBound.Fill(m_InputImageMinimum);
- m_UpperBound.Fill(m_InputImageMaximum+1);
- typename HistogramType::SizeType size( MeasurementVectorSize );
- size.Fill(m_NumberOfBinsPerAxis);
- m_Histogram->Initialize(size, m_LowerBound, m_UpperBound);
-
- /** calulate minimum offset and set it as neigborhood radius **/
- 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;
- }
- }
- m_NeighborhoodRadius.Fill(minRadius);
-
-}
-
template <class TInputImage, class TOutputImage>
void
ScalarImageToTexturesFilter<TInputImage, TOutputImage>
@@ -321,20 +328,12 @@ ScalarImageToTexturesFilter<TInputImage, TOutputImage>
inputRegion.SetSize(inputSize);
inputRegion.Crop(inputPtr->GetRequestedRegion());
+ CooccurrenceIndexedListPointerType m_GLCIList = CooccurrenceIndexedListType::New();
+ m_GLCIList->Initialize(m_HistSize);
+
typedef itk::ConstNeighborhoodIterator< InputImageType > NeighborhoodIteratorType;
NeighborhoodIteratorType neighborIt;
-
- //initalize lookup array.
- //unsigned int historamSize = m_Histogram->GetSize(0);
- LookupArrayType lookupArray(m_NumberOfBinsPerAxis * m_NumberOfBinsPerAxis);
-
- lookupArray.Fill(-1);
-
- GreyLevelCooccurrenceListType glcList;
- TotalAbsoluteFrequencyType totalFrequency = 0;
-
neighborIt = NeighborhoodIteratorType(m_NeighborhoodRadius, inputPtr, inputRegion);
-
for ( neighborIt.GoToBegin(); !neighborIt.IsAtEnd(); ++neighborIt )
{
@@ -363,68 +362,19 @@ ScalarImageToTexturesFilter<TInputImage, TOutputImage>
// is out-of-bounds.
}
- unsigned int instanceId = 0;
- HistogramIndexType instanceIndex;
+ CooccurrenceIndexType instanceIndex;
MeasurementVectorType measurement( MeasurementVectorSize );
measurement[0] = centerPixelIntensity;
measurement[1] = pixelIntensity;
-
//Get Index of the histogram for the given pixel pair;
m_Histogram->GetIndex(measurement, instanceIndex);
-
- //Find the 1D index of the historam index. This index is used to check the
- //entry in lookup array.
- instanceId = instanceIndex[1] * m_NumberOfBinsPerAxis + instanceIndex[0];
-
- if( lookupArray[instanceId] < 0)
- {
- lookupArray[instanceId] = glcList.size();
- CooccurrenceType cooccur;
- cooccur.index = instanceIndex;
- cooccur.frequency = 1;
- glcList.push_back(cooccur);
- }
- else
- {
- int vindex = lookupArray[instanceId];
- glcList[vindex].frequency++;
- }
-
- //For symmetry store the same pixel pair.
- measurement[1] = centerPixelIntensity;
- measurement[0] = pixelIntensity;
- m_Histogram->GetIndex(measurement, instanceIndex);
- instanceId = instanceIndex[1] * m_NumberOfBinsPerAxis + instanceIndex[0];
-
- if( lookupArray[instanceId] < 0)
- {
- lookupArray[instanceId] = glcList.size();
- CooccurrenceType cooccur;
- cooccur.index = instanceIndex;
- cooccur.frequency = 1;
- glcList.push_back(cooccur);
- }
- else
- {
- int vindex = lookupArray[instanceId];
- glcList[vindex].frequency = glcList[vindex].frequency + 1;
- }
-
- // Increment total frequency by two as we consider symmetry of
- // cooccurrence pairs
- totalFrequency = totalFrequency + 2;
+ m_GLCIList->AddPairToList(instanceIndex);
}
- GreyLevelCooccurrenceListIteratorType vectorIt;
- GreyLevelCooccurrenceListConstIteratorType covectorIt;
- //Normalize the GreyLevelCooccurrenceListType
- vectorIt = glcList.begin();
- while( vectorIt != glcList.end())
- {
- (*vectorIt).frequency = (*vectorIt).frequency / totalFrequency;
- ++vectorIt;
- }
+ m_GLCIList->Normalize();
+ VectorConstIteratorType constVectorIt;
+ VectorType glcList = m_GLCIList->GetVector();
double pixelMean;
double marginalMean;
@@ -456,11 +406,11 @@ ScalarImageToTexturesFilter<TInputImage, TOutputImage>
}
//Compute textures
- covectorIt = glcList.begin();
- while( covectorIt != glcList.end())
+ constVectorIt = glcList.begin();
+ while( constVectorIt != glcList.end())
{
- HistogramIndexType index = (*covectorIt).index;
- RelativeFrequencyType frequency = (*covectorIt).frequency;
+ CooccurrenceIndexType index = (*constVectorIt).index;
+ RelativeFrequencyType frequency = (*constVectorIt).frequency;
energy += frequency * frequency;
entropy -= ( frequency > 0.0001 ) ? frequency *vcl_log(frequency) / log2:0;
correlation += ( ( index[0] - pixelMean ) * ( index[1] - pixelMean ) * frequency ) / pixelVarianceSquared;
@@ -469,7 +419,7 @@ ScalarImageToTexturesFilter<TInputImage, TOutputImage>
clusterShade += vcl_pow( ( index[0] - pixelMean ) + ( index[1] - pixelMean ), 3 ) * frequency;
clusterProminence += vcl_pow( ( index[0] - pixelMean ) + ( index[1] - pixelMean ), 4 ) * frequency;
haralickCorrelation += index[0] * index[1] * frequency;
- ++covectorIt;
+ ++constVectorIt;
}
//todo avoid division by zero also here?
@@ -502,11 +452,10 @@ ScalarImageToTexturesFilter<TInputImage, TOutputImage>
template <class TInputImage, class TOutputImage>
void
-ScalarImageToTexturesFilter<TInputImage, TOutputImage>::ComputeMeansAndVariances(GreyLevelCooccurrenceListType& frequencyVector,
- double & pixelMean,
- double & marginalMean,
- double & marginalDevSquared,
- double & pixelVariance)
+ScalarImageToTexturesFilter<TInputImage, TOutputImage>::
+ComputeMeansAndVariances(VectorType& frequencyVector, double & pixelMean,
+ double & marginalMean, double & marginalDevSquared,
+ double & pixelVariance)
{
// 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
@@ -524,15 +473,15 @@ ScalarImageToTexturesFilter<TInputImage, TOutputImage>::ComputeMeansAndVariances
// Ok, now do the first pass through the histogram to get the marginal sums
// and compute the pixel mean
- GreyLevelCooccurrenceListConstIteratorType covectorIt;
- covectorIt = frequencyVector.begin();
- while( covectorIt != frequencyVector.end())
+ VectorConstIteratorType constVectorIt;
+ constVectorIt = frequencyVector.begin();
+ while( constVectorIt != frequencyVector.end())
{
- RelativeFrequencyType frequency = (*covectorIt).frequency;
- HistogramIndexType index = (*covectorIt).index;
+ RelativeFrequencyType frequency = (*constVectorIt).frequency;
+ CooccurrenceIndexType index = (*constVectorIt).index;
pixelMean += index[0] * frequency;
marginalSums[index[0]] += frequency;
- ++covectorIt;
+ ++constVectorIt;
}
/* Now get the mean and deviaton of the marginal sums.
@@ -564,13 +513,13 @@ ScalarImageToTexturesFilter<TInputImage, TOutputImage>::ComputeMeansAndVariances
marginalDevSquared = marginalDevSquared / m_NumberOfBinsPerAxis;
pixelVariance = 0;
- covectorIt = frequencyVector.begin();
- while( covectorIt != frequencyVector.end())
+ constVectorIt = frequencyVector.begin();
+ while( constVectorIt != frequencyVector.end())
{
- RelativeFrequencyType frequency = (*covectorIt).frequency;
- HistogramIndexType index = (*covectorIt).index;
+ RelativeFrequencyType frequency = (*constVectorIt).frequency;
+ CooccurrenceIndexType index = (*constVectorIt).index;
pixelVariance += ( index[0] - pixelMean ) * ( index[0] - pixelMean ) * frequency;
- ++covectorIt;
+ ++constVectorIt;
}
delete[] marginalSums;
--
GitLab