diff --git a/Code/FeatureExtraction/otbScalarImageToAdvancedTexturesFilter.h b/Code/FeatureExtraction/otbScalarImageToAdvancedTexturesFilter.h
index ee767586604a81d5481321ab7285e963e4262aee..7b91d5d0e1eff19d3e4d3b6dd95f014b15d697da 100644
--- a/Code/FeatureExtraction/otbScalarImageToAdvancedTexturesFilter.h
+++ b/Code/FeatureExtraction/otbScalarImageToAdvancedTexturesFilter.h
@@ -19,14 +19,67 @@
#define __otbScalarImageToAdvancedTexturesFilter_h
#include "itkImageToImageFilter.h"
-#include "itkScalarImageToCooccurrenceMatrixFilter.h"
-#include "otbGreyLevelCooccurrenceMatrixAdvancedTextureCoefficientsCalculator.h"
+#include "itkHistogram.h"
+#include "itkArray.h"
+#include <vector>
namespace otb
{
/**
* \class ScalarImageToAdvancedTexturesFilter
- * \brief TODO
+ * \brief This class computes 9 haralick textures from the given input image.
+ *
+ * The textures are computed over a sliding window with user defined radius:
+ * (where \f$ g(i, j) \f$ is the element in cell i, j of a normalized GLCIL):
+ *
+ * "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)
+ *
+ * "Sum Variance" \f$ = f_7 = \sum_{i}(i - f_8)^2 g_{x+y}(i) \f$
+ *
+ * "Sum Entropy" \f$= f_8 = -\sum_{i}g_{x+y}(i) log (g_{x+y}(i)) \f$
+ *
+ * "Difference variance" \f$ = f_10 = variance of g_{x-y}(i)
+ *
+ * "Difference entropy" \f$ = f_11 = -\sum_{i}g_{x-y}(i) log (g_{x-y}(i)) \f$
+ *
+ * "Information Measures of Correlation IC1" \f$ = f_12 = \frac{f_9 - HXY1}{H} \f$
+ *
+ * "Information Measures of Correlation IC2" \f$ = f_13 = \sqrt{1 - \exp{-2}|HXY2 - f_9|} \f$
+ *
+ * Above, \f$ \mu = \f$ (weighted pixel average) \f$ = \sum_{i, j}i \cdot g(i, j) =
+ * \sum_{i, j}j \cdot g(i, j) \f$ (due to matrix summetry), and
+ *
+ * \f$ \g_{x+y}(k) = \sum_{i}\sum_{j}g(i)\f$ where \f$ i+j=k \f$ and \f$ k = 2, 3, .., 2N_[g} \f$ and
+ *
+ * \f$ \g_{x-y}(k) = \sum_{i}\sum_{j}g(i)\f$ where \f$ i-j=k \f$ and \f$ k = 0, 1, .., N_[g}-1 \f$
+ *
+ * Print references:
+ *
+ * Haralick, R.M., K. Shanmugam and I. Dinstein. 1973. Textural Features for
+ * Image Classification. IEEE Transactions on Systems, Man and Cybernetics.
+ * SMC-3(6):610-620.
+ *
+ * David A. Clausi and Yongping Zhao. 2002. Rapid extraction of image texture by
+ * co-occurrence using a hybrid data structure. Comput. Geosci. 28, 6 (July
+ * 2002), 763-774. DOI=10.1016/S0098-3004(01)00108-X
+ * http://dx.doi.org/10.1016/S0098-3004(01)00108-X
+ *
+ * de O.Bastos, L.; Liatsis, P.; Conci, A., Automatic texture segmentation based
+ * on k-means clustering and efficient calculation of co-occurrence
+ * features. Systems, Signals and Image Processing, 2008. IWSSIP 2008. 15th
+ * International Conference on , vol., no., pp.141,144, 25-28 June 2008
+ * doi: 10.1109/IWSSIP.2008.4604387
+ *
+ * Neighborhood size can be set using the SetRadius() method. Offset for co-occurence estimation
+ * is set using the SetOffset() method.
+ *
+ * \sa ScalarImageToCooccurrenceIndexedList
+ * \sa ScalarImageToTexturesFiler
+ *
*/
template<class TInpuImage, class TOutputImage>
class ScalarImageToAdvancedTexturesFilter : public itk::ImageToImageFilter
@@ -50,19 +103,45 @@ 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;
typedef TOutputImage OutputImageType;
typedef typename OutputImageType::Pointer OutputImagePointerType;
typedef typename OutputImageType::RegionType OutputRegionType;
- /** 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 GreyLevelCooccurrenceMatrixAdvancedTextureCoefficientsCalculator
- <HistogramType> TextureCoefficientsCalculatorType;
- typedef typename TextureCoefficientsCalculatorType::Pointer TextureCoefficientsCalculatorPointerType;
+ typedef typename itk::NumericTraits< InputPixelType >::RealType MeasurementType;
+ typedef itk::Statistics::DenseFrequencyContainer2 THistogramFrequencyContainer;
+ typedef itk::Statistics::Histogram< MeasurementType, THistogramFrequencyContainer > HistogramType;
+ typedef typename HistogramType::Pointer HistogramPointer;
+ typedef typename HistogramType::ConstPointer HistogramConstPointer;
+ typedef typename HistogramType::MeasurementVectorType MeasurementVectorType;
+ typedef typename HistogramType::RelativeFrequencyType RelativeFrequencyType;
+ typedef typename HistogramType::TotalAbsoluteFrequencyType TotalAbsoluteFrequencyType;
+ typedef typename HistogramType::IndexType HistogramIndexType;
+
+ /** Lookup array used to store the index of the given pixel pair. This index
+ * is calculated from the HistogramType */
+ typedef itk::Array<int> LookupArrayType;
+
+ /** 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;
+
+ itkStaticConstMacro(MeasurementVectorSize, int, 2 );
/** Set the radius of the window on which textures will be computed */
itkSetMacro(Radius, SizeType);
@@ -125,8 +204,13 @@ 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);
@@ -140,6 +224,9 @@ private:
/** Radius of the window on which to compute textures */
SizeType m_Radius;
+ /** Radius of the neighborhood iterator which is minumum of m_Radius */
+ SizeType m_NeighborhoodRadius;
+
/** Offset for co-occurence */
OffsetType m_Offset;
@@ -151,6 +238,10 @@ private:
/** Input image maximum */
InputPixelType m_InputImageMaximum;
+
+ /** Histogram instance used to get index for the pixel pair */
+ HistogramPointer m_Histogram;
+
};
} // End namespace otb
diff --git a/Code/FeatureExtraction/otbScalarImageToAdvancedTexturesFilter.txx b/Code/FeatureExtraction/otbScalarImageToAdvancedTexturesFilter.txx
index f59c3e4a54e51228095069a60574aae4c9e7725c..6d959997c6fd81f5973839dd2a45ceb029719e65 100644
--- a/Code/FeatureExtraction/otbScalarImageToAdvancedTexturesFilter.txx
+++ b/Code/FeatureExtraction/otbScalarImageToAdvancedTexturesFilter.txx
@@ -20,8 +20,12 @@
#include "otbScalarImageToAdvancedTexturesFilter.h"
#include "itkImageRegionIteratorWithIndex.h"
+#include "itkConstNeighborhoodIterator.h"
#include "itkImageRegionIterator.h"
#include "itkProgressReporter.h"
+#include "itkNumericTraits.h"
+#include "vnl/vnl_math.h"
+#include <algorithm>
namespace otb
{
@@ -31,7 +35,7 @@ ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
m_Offset(),
m_NumberOfBinsPerAxis(8),
m_InputImageMinimum(0),
- m_InputImageMaximum(256)
+ m_InputImageMaximum(255)
{
// There are 9 outputs corresponding to the 8 textures indices
this->SetNumberOfRequiredOutputs(9);
@@ -230,22 +234,54 @@ ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
}
}
+
+template <class TInputImage, class TOutputImage>
+void
+ScalarImageToAdvancedTexturesFilter<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
ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
::ThreadedGenerateData(const OutputRegionType& outputRegionForThread, itk::ThreadIdType threadId)
{
// 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();
- OutputImagePointerType sumEntropytPtr = this->GetSumEntropyOutput();
- OutputImagePointerType differenceEntropyPtr = this->GetDifferenceEntropyOutput();
+ InputImagePointerType inputPtr = const_cast<InputImageType *>(this->GetInput());
+ OutputImagePointerType meanPtr = this->GetMeanOutput();
+ OutputImagePointerType variancePtr = this->GetVarianceOutput();
+ OutputImagePointerType sumAveragePtr = this->GetSumAverageOutput();
+ OutputImagePointerType sumVariancePtr = this->GetSumVarianceOutput();
+ OutputImagePointerType sumEntropytPtr = this->GetSumEntropyOutput();
+ OutputImagePointerType differenceEntropyPtr = this->GetDifferenceEntropyOutput();
OutputImagePointerType differenceVariancePtr = this->GetDifferenceVarianceOutput();
- OutputImagePointerType ic1Ptr = this->GetIC1Output();
- OutputImagePointerType ic2Ptr = this->GetIC2Output();
+ OutputImagePointerType ic1Ptr = this->GetIC1Output();
+ OutputImagePointerType ic2Ptr = this->GetIC2Output();
// Build output iterators
itk::ImageRegionIteratorWithIndex<OutputImageType> varianceIt(variancePtr, outputRegionForThread);
@@ -269,15 +305,7 @@ ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
ic1It.GoToBegin();
ic2It.GoToBegin();
-
- // Build the co-occurence matrix generator
- CoocurrenceMatrixGeneratorPointerType coOccurenceMatrixGenerator = CoocurrenceMatrixGeneratorType::New();
- coOccurenceMatrixGenerator->SetOffset(m_Offset);
- coOccurenceMatrixGenerator->SetNumberOfBinsPerAxis(m_NumberOfBinsPerAxis);
- coOccurenceMatrixGenerator->SetPixelValueMinMax(m_InputImageMinimum, m_InputImageMaximum);
-
- // Build the texture calculator
- TextureCoefficientsCalculatorPointerType texturesCalculator = TextureCoefficientsCalculatorType::New();
+ const double log2 = vcl_log(2.0);
// Set-up progress reporting
itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
@@ -294,20 +322,15 @@ ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
&& !ic2It.IsAtEnd())
{
// 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;
- // First, apply offset
+ // First, create an window for neighborhood iterator based on m_Radius
+ // For example, if xradius and yradius is 2. window size is 5x5 (2 * radius + 1).
for (unsigned int dim = 0; dim < InputImageType::ImageDimension; ++dim)
{
- inputIndex[dim] = std::min(
- static_cast<int>(varianceIt.GetIndex()[dim] - m_Radius[dim]),
- static_cast<int>(varianceIt.GetIndex()[dim] - m_Radius[dim] + m_Offset[dim])
- );
- inputSize[dim] = 2 * m_Radius[dim] + 1 + std::abs(m_Offset[dim]);
-
- inputIndexWithTwiceOffset[dim] = static_cast<int>(varianceIt.GetIndex()[dim] - m_Radius[dim] - std::abs(m_Offset[dim]));
- inputSizeWithTwiceOffset[dim] = inputSize[dim] + std::abs(m_Offset[dim]);
+ inputIndex[dim] = varianceIt.GetIndex()[dim] - m_Radius[dim];
+ inputSize[dim] = 2 * m_Radius[dim] + 1;
}
// Build the input region
@@ -316,56 +339,232 @@ ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
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 *outputIt
- 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)
+ typedef itk::ConstNeighborhoodIterator< InputImageType > NeighborhoodIteratorType;
+ NeighborhoodIteratorType neighborIt;
+
+ //initalize lookup array.
+ 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 )
+ {
+ const InputPixelType centerPixelIntensity = neighborIt.GetCenterPixel();
+ if ( centerPixelIntensity < m_InputImageMinimum
+ || centerPixelIntensity > m_InputImageMaximum )
+ {
+ continue; // don't put a pixel in the histogram if the value
+ // is out-of-bounds.
+ }
+
+ 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.
+ }
+ if ( pixelIntensity < m_InputImageMinimum
+ || pixelIntensity > m_InputImageMaximum )
{
- itLocalInputImage.Set(itInputPtr.Get());
+ continue; // don't put a pixel in the histogram if the value
+ // is out-of-bounds.
}
- /*********************************************************************************/
- InputImagePointerType maskImage = InputImageType::New();
- maskImage->SetRegions(inputRegionWithTwiceOffset);
- maskImage->Allocate();
- maskImage->FillBuffer(0);
+ unsigned int instanceId = 0;
+ HistogramIndexType 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;
+ }
+
+ 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_DifferenceEntropy = itk::NumericTraits< MeasurementType >::Zero;
+ MeasurementType m_DifferenceVariance = itk::NumericTraits< MeasurementType >::Zero;
+ MeasurementType m_IC1 = itk::NumericTraits< MeasurementType >::Zero;
+ MeasurementType m_IC2 = itk::NumericTraits< MeasurementType >::Zero;
+
+ double Entropy = 0;
- ImageRegionIteratorType itMask(maskImage, inputRegion);
- for (itMask.GoToBegin(); !itMask.IsAtEnd(); ++itMask)
+ long unsigned int histSize = m_Histogram->GetSize()[0];
+ long unsigned int twiceHistSize = 2 * m_Histogram->GetSize()[0];
+
+ typedef itk::Array<double> DoubleArrayType;
+ DoubleArrayType hx(histSize);
+ DoubleArrayType hy(histSize);
+ DoubleArrayType pdxy(twiceHistSize);
+
+ for(long unsigned int i = 0; i < histSize; i++)
+ {
+ hx[i] = 0.0;
+ hy[i] = 0.0;
+ pdxy[i] = 0.0;
+ }
+ for(long unsigned int i = histSize; i < twiceHistSize; i++)
+ {
+ pdxy[i] = 0.0;
+ }
+
+ /* hx.Fill(0.0); hy.Fill(0.0); pdxy.Fill(0.0); */
+ double hxy1 = 0;
+
+ GreyLevelCooccurrenceListIteratorType vectorIt;
+ GreyLevelCooccurrenceListConstIteratorType covectorIt;
+ //Normalize the GreyLevelCooccurrenceListType
+ //Compute Mean, Entropy (f12), hx, hy, pdxy
+ vectorIt = glcList.begin();
+ while( vectorIt != glcList.end())
+ {
+ HistogramIndexType index = (*vectorIt).index;
+ double frequency = static_cast<double>((*vectorIt).frequency) / static_cast<double>(totalFrequency);
+ m_Mean += static_cast<double>(index[0]) * frequency;
+ Entropy -= (frequency > 0.0001) ? frequency * vcl_log(frequency) / log2 : 0.;
+ //update normalized frequency
+ (*vectorIt).frequency = frequency;
+
+ unsigned int i = index[1];
+ unsigned int j = index[0];
+ hx[j] += frequency;
+ hy[i] += frequency;
+
+ if( i+j > histSize-1)
+ {
+ pdxy[i+j] += frequency;
+ }
+ if( i <= j )
+ {
+ pdxy[j-i] += frequency;
+ }
+
+ ++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();
+ while( vectorIt != glcList.end())
{
- itMask.Set(1);
+ double frequency = (*vectorIt).frequency;
+ HistogramIndexType index = (*vectorIt).index;
+ unsigned int i = index[1];
+ unsigned int j = index[0];
+ double index0 = static_cast<double>(index[0]);
+ m_Variance += ((index0 - m_Mean) * (index0 - m_Mean)) * frequency;
+ double pipj = hx[j] * hy[i];
+ hxy1 -= (pipj > 0.0001) ? frequency * vcl_log(pipj) : 0.;
+ ++vectorIt;
}
- // Compute the co-occurence matrix
- coOccurenceMatrixGenerator->SetInput(localInputImage);
- coOccurenceMatrixGenerator->SetMaskImage(maskImage);
- coOccurenceMatrixGenerator->SetInsidePixelValue(1);
- coOccurenceMatrixGenerator->Update();
+ //iterate histSize to compute sumEntropy
+ double PSSquareCumul = 0;
+ for(long unsigned int k = histSize; k < twiceHistSize; k++)
+ {
+ m_SumAverage += k * pdxy[k];
+ m_SumEntropy -= (pdxy[k] > 0.0001) ? pdxy[k] * vcl_log(pdxy[k]) / log2 : 0;
+ PSSquareCumul += k * k * pdxy[k];
+ }
+ m_SumVariance = PSSquareCumul - m_SumAverage * m_SumAverage;
+
+ double PDSquareCumul = 0;
+ double PDCumul = 0;
+ double hxCumul = 0;
+ double hyCumul = 0;
+
+ for (long unsigned int i = 0; i < minSizeHist; ++i)
+ {
+ double pdTmp = pdxy[i];
+ PDCumul += i * pdTmp;
+ m_DifferenceEntropy -= (pdTmp > 0.0001) ? pdTmp * vcl_log(pdTmp) / log2 : 0;
+ PDSquareCumul += i * i * pdTmp;
+
+ //comput hxCumul and hyCumul
+ double marginalfreq = hx[i];
+ hxCumul += (marginalfreq > 0.0001) ? vcl_log (marginalfreq) * marginalfreq : 0;
+
+ marginalfreq = hy[i];
+ hyCumul += (marginalfreq > 0.0001) ? vcl_log (marginalfreq) * marginalfreq : 0;
+ }
+ m_DifferenceVariance = PDSquareCumul - PDCumul * PDCumul;
+
+ /* 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 */
+ double hxy2 = 0;
+ for(long unsigned int i = 0; i < histSize; ++i)
+ {
+ for(long unsigned int j = 0; j < histSize; ++j)
+ {
+ double pipj = hx[j] * hy[i];
+ hxy2 -= (pipj > 0.0001) ? pipj * vcl_log(pipj) : 0.;
+ }
+ }
- // Compute textures indices
- texturesCalculator->SetHistogram(const_cast<HistogramType*>(coOccurenceMatrixGenerator->GetOutput())); //FIXME const correctness
- texturesCalculator->Compute();
+ //Information measures of correlation 1 & 2
+ m_IC1 = (vcl_abs(std::max (hxCumul, hyCumul)) > 0.0001) ? (Entropy - hxy1) / (std::max (hxCumul, hyCumul)) : 0;
+ m_IC2 = 1 - vcl_exp (-2. * vcl_abs (hxy2 - Entropy));
+ m_IC2 = (m_IC2 >= 0) ? vcl_sqrt (m_IC2) : 0;
// Fill outputs
- varianceIt.Set(texturesCalculator->GetVariance());
- meanIt.Set(texturesCalculator->GetMean());
- sumAverageIt.Set(texturesCalculator->GetSumAverage());
- sumVarianceIt.Set(texturesCalculator->GetSumVariance());
- sumEntropytIt.Set(texturesCalculator->GetSumEntropy());
- differenceEntropyIt.Set(texturesCalculator->GetDifferenceEntropy());
- differenceVarianceIt.Set(texturesCalculator->GetDifferenceVariance());
- ic1It.Set(texturesCalculator->GetIC1());
- ic2It.Set(texturesCalculator->GetIC2());
+ varianceIt.Set(m_Variance);
+ meanIt.Set(m_Mean);
+ sumAverageIt.Set(m_SumAverage);
+ sumVarianceIt.Set(m_SumVariance);
+ sumEntropytIt.Set(m_SumEntropy);
+ differenceEntropyIt.Set(m_DifferenceEntropy);
+ differenceVarianceIt.Set(m_DifferenceVariance);
+ ic1It.Set(m_IC1);
+ ic2It.Set(m_IC2);
// Update progress
progress.CompletedPixel();
diff --git a/Testing/Code/FeatureExtraction/otbScalarImageToAdvancedTexturesFilter.cxx b/Testing/Code/FeatureExtraction/otbScalarImageToAdvancedTexturesFilter.cxx
index 41231a7b89b4b61e5645423dda5ae507da6d2619..9bdb81b9965e656727df7762e1567636fbd3493e 100644
--- a/Testing/Code/FeatureExtraction/otbScalarImageToAdvancedTexturesFilter.cxx
+++ b/Testing/Code/FeatureExtraction/otbScalarImageToAdvancedTexturesFilter.cxx
@@ -69,7 +69,7 @@ int otbScalarImageToAdvancedTexturesFilter(int argc, char * argv[])
otb::StandardFilterWatcher watcher(filter, "Textures filter");
filter->SetInputImageMinimum(0);
- filter->SetInputImageMaximum(256);
+ filter->SetInputImageMaximum(255);
// Write outputs
std::ostringstream oss;