Skip to content
Snippets Groups Projects
Commit f5e36f4a authored by Rashad Kanavath's avatar Rashad Kanavath
Browse files

ENH: Improvement of ScalarImageToPanTexTextureFiler with otbGreyLevelCooccurrenceIndexedList

parent 0b331217
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Code/FeatureExtraction/otbScalarImageToPanTexTextureFilter.h
+ 32
− 9
View file @ f5e36f4a
......@@ -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
......
Code/FeatureExtraction/otbScalarImageToPanTexTextureFilter.txx
+ 90
− 60
View file @ f5e36f4a
......@@ -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)
......
Testing/Code/FeatureExtraction/otbScalarImageToPanTexTextureFilter.cxx
+ 1
− 3
View file @ f5e36f4a
......@@ -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;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment