diff --git a/Code/FeatureExtraction/otbScalarImageToPanTexTextureFilter.h b/Code/FeatureExtraction/otbScalarImageToPanTexTextureFilter.h
index 01f2f7867865a452ee830c52266818512f035686..5ea92e62be6f08362fb976ed713ecf12293ea505 100644
--- a/Code/FeatureExtraction/otbScalarImageToPanTexTextureFilter.h
+++ b/Code/FeatureExtraction/otbScalarImageToPanTexTextureFilter.h
@@ -18,8 +18,9 @@
#ifndef __otbScalarImageToPanTexTextureFilter_h
#define __otbScalarImageToPanTexTextureFilter_h
+#include "otbGreyLevelCooccurrenceIndexedList.h"
#include "itkImageToImageFilter.h"
-#include "itkScalarImageToCooccurrenceMatrixFilter.h"
+#include "itkHistogram.h"
namespace otb
{
@@ -66,19 +67,35 @@ public:
typedef typename InputImageType::PixelType InputPixelType;
typedef typename InputImageType::RegionType InputRegionType;
typedef typename InputRegionType::SizeType SizeType;
+ typedef typename InputImageType::OffsetType OffsetType;
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 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::SizeType HistogramSizeType;
+ typedef typename HistogramType::IndexType CooccurrenceIndexType;
+
+ typedef GreyLevelCooccurrenceIndexedList< HistogramType > CooccurrenceIndexedListType;
+ typedef typename CooccurrenceIndexedListType::Pointer CooccurrenceIndexedListPointerType;
+ typedef typename CooccurrenceIndexedListType::ConstPointer CooccurrenceIndexedListConstPointerType;
+
+ typedef typename CooccurrenceIndexedListType::VectorType VectorType;
+ typedef typename VectorType::iterator VectorIteratorType;
+ typedef typename VectorType::const_iterator VectorConstIteratorType;
+
typedef typename std::vector<OffsetType> OffsetListType;
- typedef typename CoocurrenceMatrixGeneratorType::HistogramType HistogramType;
- typedef typename HistogramType::AbsoluteFrequencyType FrequencyType; //FIXME stat framework
- typedef typename HistogramType::MeasurementType MeasurementType;
- typedef typename HistogramType::Iterator HistogramIterator;
+
+ void SetBinsAndMinMax(unsigned int bins, InputPixelType min, InputPixelType max);
+
+ itkStaticConstMacro(MeasurementVectorSize, int, 2 );
/** Set the radius of the window on which textures will be computed */
itkSetMacro(Radius, SizeType);
@@ -134,6 +151,12 @@ private:
/** Input image maximum */
InputPixelType m_InputImageMaximum;
+
+ /** Histogram instance used to get index for the pixel pair */
+ HistogramPointer m_Histogram;
+
+ HistogramSizeType m_HistSize;
+
};
} // End namespace otb
diff --git a/Code/FeatureExtraction/otbScalarImageToPanTexTextureFilter.txx b/Code/FeatureExtraction/otbScalarImageToPanTexTextureFilter.txx
index da5e6856e9506a0732456300ad221b559ab4692d..b3515e841e2ce86e4e905335e7c236128038425d 100644
--- a/Code/FeatureExtraction/otbScalarImageToPanTexTextureFilter.txx
+++ b/Code/FeatureExtraction/otbScalarImageToPanTexTextureFilter.txx
@@ -20,8 +20,10 @@
#include "otbScalarImageToPanTexTextureFilter.h"
#include "itkImageRegionIteratorWithIndex.h"
+#include "itkConstNeighborhoodIterator.h"
#include "itkImageRegionIterator.h"
#include "itkProgressReporter.h"
+#include "itkNumericTraits.h"
namespace otb
{
@@ -29,8 +31,9 @@ template <class TInputImage, class TOutputImage>
ScalarImageToPanTexTextureFilter<TInputImage, TOutputImage>
::ScalarImageToPanTexTextureFilter() : m_Radius(),
m_NumberOfBinsPerAxis(8),
+ m_HistSize(2),
m_InputImageMinimum(0),
- m_InputImageMaximum(256)
+ m_InputImageMaximum(255)
{
// There are 1 output corresponding to the Pan Tex texture indice
this->SetNumberOfRequiredOutputs(1);
@@ -67,6 +70,29 @@ ScalarImageToPanTexTextureFilter<TInputImage, TOutputImage>
::~ScalarImageToPanTexTextureFilter()
{}
+template <class TInputImage, class TOutputImage>
+void
+ScalarImageToPanTexTextureFilter<TInputImage, TOutputImage>
+::SetBinsAndMinMax(unsigned int numberOfBinsPerAxis,
+ InputPixelType inputImageMinimum,
+ InputPixelType inputImageMaximum)
+{
+ /** 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
ScalarImageToPanTexTextureFilter<TInputImage, TOutputImage>
@@ -144,86 +170,90 @@ ScalarImageToPanTexTextureFilter<TInputImage, TOutputImage>
OffsetType currentOffset = *offIt;
// 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>(outputIt.GetIndex()[dim] - m_Radius[dim]),
- static_cast<int>(outputIt.GetIndex()[dim] - m_Radius[dim] + currentOffset[dim])
- );
- inputSize[dim] = 2 * m_Radius[dim] + 1 + std::abs(currentOffset[dim]);
-
- inputIndexWithTwiceOffset[dim] = static_cast<int>(outputIt.GetIndex()[dim] - m_Radius[dim] - std::abs(currentOffset[dim]));
- inputSizeWithTwiceOffset[dim] = inputSize[dim] + std::abs(currentOffset[dim]);
+ inputIndex[dim] = outputIt.GetIndex()[dim] - m_Radius[dim];
+ inputSize[dim] = 2 * m_Radius[dim] + 1;
}
-
// Build the input region
InputRegionType inputRegion;
inputRegion.SetIndex(inputIndex);
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)
+ CooccurrenceIndexedListPointerType m_GLCIList = CooccurrenceIndexedListType::New();
+ m_GLCIList->Initialize(m_HistSize);
+
+ SizeType neighborhoodRadius;
+ /** calulate minimum offset and set it as neigborhood radius **/
+ unsigned int minRadius = 0;
+ for ( unsigned int i = 0; i < currentOffset.GetOffsetDimension(); i++ )
{
- itLocalInputImage.Set(itInputPtr.Get());
+ unsigned int distance = vcl_abs(currentOffset[i]);
+ if ( distance > minRadius )
+ {
+ minRadius = distance;
+ }
}
- /*********************************************************************************/
-
- InputImagePointerType maskImage = InputImageType::New();
- maskImage->SetRegions(inputRegionWithTwiceOffset);
- maskImage->Allocate();
- maskImage->FillBuffer(0);
+ neighborhoodRadius.Fill(minRadius);
- ImageRegionIteratorType itMask(maskImage, inputRegion);
- for (itMask.GoToBegin(); !itMask.IsAtEnd(); ++itMask)
+ typedef itk::ConstNeighborhoodIterator< InputImageType > NeighborhoodIteratorType;
+ NeighborhoodIteratorType neighborIt;
+ neighborIt = NeighborhoodIteratorType(neighborhoodRadius, inputPtr, inputRegion);
+ for ( neighborIt.GoToBegin(); !neighborIt.IsAtEnd(); ++neighborIt )
{
- itMask.Set(1);
- }
+ 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(currentOffset, pixelInBounds);
- // Build the co-occurence matrix generator
- CoocurrenceMatrixGeneratorPointerType coOccurenceMatrixGenerator = CoocurrenceMatrixGeneratorType::New();
- coOccurenceMatrixGenerator->SetInput(localInputImage);
- coOccurenceMatrixGenerator->SetOffset(currentOffset);
- coOccurenceMatrixGenerator->SetNumberOfBinsPerAxis(m_NumberOfBinsPerAxis);
- coOccurenceMatrixGenerator->SetPixelValueMinMax(m_InputImageMinimum, m_InputImageMaximum);
+ if ( !pixelInBounds )
+ {
+ continue; // don't put a pixel in the histogram if it's out-of-bounds.
+ }
- // Compute the co-occurence matrix
- coOccurenceMatrixGenerator->SetMaskImage(maskImage);
- coOccurenceMatrixGenerator->SetInsidePixelValue(1);
- coOccurenceMatrixGenerator->Update();
+ if ( pixelIntensity < m_InputImageMinimum
+ || pixelIntensity > m_InputImageMaximum )
+ {
+ continue; // don't put a pixel in the histogram if the value
+ // is out-of-bounds.
+ }
+
+ 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);
+ m_GLCIList->AddPairToList(instanceIndex);
+ }
- typename HistogramType::ConstPointer histo = coOccurenceMatrixGenerator->GetOutput();
+ m_GLCIList->Normalize();
+ VectorConstIteratorType constVectorIt;
+ VectorType glcList = m_GLCIList->GetVector();
+ //Compute inertia aka contrast
double inertia = 0;
- typename HistogramType::TotalAbsoluteFrequencyType totalFrequency = histo->GetTotalFrequency();
- typename HistogramType::ConstIterator itr = histo->Begin();
- typename HistogramType::ConstIterator end = histo->End();
- for(; itr != end; ++itr )
+ constVectorIt = glcList.begin();
+ while( constVectorIt != glcList.end())
{
- MeasurementType frequency = itr.GetFrequency();
- if (frequency == 0)
- {
- continue; // no use doing these calculations if we're just multiplying by zero.
- }
- typename HistogramType::IndexType index = histo->GetIndex(itr.GetInstanceIdentifier());
- inertia += (index[0] - index[1]) * (index[0] - index[1]) * frequency / totalFrequency;
+ CooccurrenceIndexType index = (*constVectorIt).index;
+ RelativeFrequencyType frequency = (*constVectorIt).frequency;
+ inertia += ( index[0] - index[1] ) * ( index[0] - index[1] ) * frequency;
+ ++constVectorIt;
}
if (inertia < out)
diff --git a/Testing/Code/FeatureExtraction/otbScalarImageToPanTexTextureFilter.cxx b/Testing/Code/FeatureExtraction/otbScalarImageToPanTexTextureFilter.cxx
index 1535a2a2ac221b8f2d58ecb6a3f76de974ed48c6..2f5460b3a539d957fd5add9328050e9c58afa606 100644
--- a/Testing/Code/FeatureExtraction/otbScalarImageToPanTexTextureFilter.cxx
+++ b/Testing/Code/FeatureExtraction/otbScalarImageToPanTexTextureFilter.cxx
@@ -55,13 +55,11 @@ int otbScalarImageToPanTexTextureFilter(int argc, char * argv[])
sradius.Fill(radius);
filter->SetInput(reader->GetOutput());
- filter->SetNumberOfBinsPerAxis(nbBins);
filter->SetRadius(sradius);
otb::StandardFilterWatcher watcher(filter, "Textures filter");
- filter->SetInputImageMinimum(0);
- filter->SetInputImageMaximum(255);
+ filter->SetBinsAndMinMax(nbBins, 0, 255);
// Write outputs
std::ostringstream oss;