From bc807bce5d20293faa6d593ce4686ee2a80e2096 Mon Sep 17 00:00:00 2001
From: Guillaume Pasero <guillaume.pasero@c-s.fr>
Date: Thu, 12 Jan 2017 11:25:13 +0100
Subject: [PATCH] ENH: implement sub-sampled computation in Advanced and
HigherOrder textures
---
.../otbScalarImageToAdvancedTexturesFilter.h | 19 ++++++
...otbScalarImageToAdvancedTexturesFilter.txx | 53 +++++++++++++++-
...tbScalarImageToHigherOrderTexturesFilter.h | 20 ++++++
...ScalarImageToHigherOrderTexturesFilter.txx | 61 +++++++++++++++++--
4 files changed, 148 insertions(+), 5 deletions(-)
diff --git a/Modules/Feature/Textures/include/otbScalarImageToAdvancedTexturesFilter.h b/Modules/Feature/Textures/include/otbScalarImageToAdvancedTexturesFilter.h
index 456cf4aca4..42c1a5f073 100644
--- a/Modules/Feature/Textures/include/otbScalarImageToAdvancedTexturesFilter.h
+++ b/Modules/Feature/Textures/include/otbScalarImageToAdvancedTexturesFilter.h
@@ -164,6 +164,18 @@ public:
/** Get the input image maximum */
itkGetMacro(InputImageMaximum, InputPixelType);
+ /** Set the sub-sampling factor */
+ itkSetMacro(SubsampleFactor, SizeType);
+
+ /** Get the sub-sampling factor */
+ itkGetMacro(SubsampleFactor, SizeType);
+
+ /** Set the sub-sampling offset */
+ itkSetMacro(SubsampleOffset, OffsetType);
+
+ /** Get the sub-sampling offset */
+ itkGetMacro(SubsampleOffset, OffsetType);
+
/** Get the mean output image */
OutputImageType * GetMeanOutput();
@@ -199,6 +211,8 @@ protected:
ScalarImageToAdvancedTexturesFilter();
/** Destructor */
~ScalarImageToAdvancedTexturesFilter() ITK_OVERRIDE;
+ /** Generate the output informations */
+ void GenerateOutputInformation() ITK_OVERRIDE;
/** Generate the input requested region */
void GenerateInputRequestedRegion() ITK_OVERRIDE;
/** Before Parallel textures extraction */
@@ -231,6 +245,11 @@ private:
/** Input image maximum */
InputPixelType m_InputImageMaximum;
+ /** Sub-sampling factor */
+ SizeType m_SubsampleFactor;
+
+ /** Sub-sampling offset */
+ OffsetType m_SubsampleOffset;
};
} // End namespace otb
diff --git a/Modules/Feature/Textures/include/otbScalarImageToAdvancedTexturesFilter.txx b/Modules/Feature/Textures/include/otbScalarImageToAdvancedTexturesFilter.txx
index e0463ab5c8..80cf17d041 100644
--- a/Modules/Feature/Textures/include/otbScalarImageToAdvancedTexturesFilter.txx
+++ b/Modules/Feature/Textures/include/otbScalarImageToAdvancedTexturesFilter.txx
@@ -37,6 +37,8 @@ ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
, m_NumberOfBinsPerAxis(8)
, m_InputImageMinimum(0)
, m_InputImageMaximum(255)
+, m_SubsampleFactor()
+, m_SubsampleOffset()
{
// There are 10 outputs corresponding to the 9 textures indices
this->SetNumberOfRequiredOutputs(10);
@@ -52,6 +54,9 @@ ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
this->SetNthOutput(7, OutputImageType::New());
this->SetNthOutput(8, OutputImageType::New());
this->SetNthOutput(9, OutputImageType::New());
+
+ this->m_SubsampleFactor.Fill(1);
+ this->m_SubsampleOffset.Fill(0);
}
template <class TInputImage, class TOutputImage>
@@ -188,6 +193,38 @@ ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
return static_cast<OutputImageType *>(this->GetOutput(9));
}
+template <class TInputImage, class TOutputImage>
+void
+ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
+::GenerateOutputInformation()
+{
+ // First, call superclass implementation
+ Superclass::GenerateOutputInformation();
+
+ // Compute output size, origin & spacing
+ InputRegionType inputRegion = this->GetInput()->GetLargestPossibleRegion();
+ OutputRegionType outputRegion;
+ outputRegion.SetIndex(0,0);
+ outputRegion.SetIndex(1,0);
+ outputRegion.SetSize(0, 1 + (inputRegion.GetSize(0) - 1 - m_SubsampleOffset[0]) / m_SubsampleFactor[0]);
+ outputRegion.SetSize(1, 1 + (inputRegion.GetSize(1) - 1 - m_SubsampleOffset[1]) / m_SubsampleFactor[1]);
+
+ typename OutputImageType::SpacingType outSpacing = this->GetInput()->GetSpacing();
+ outSpacing[0] *= m_SubsampleFactor[0];
+ outSpacing[1] *= m_SubsampleFactor[1];
+
+ typename OutputImageType::PointType outOrigin;
+ this->GetInput()->TransformIndexToPhysicalPoint(inputRegion.GetIndex()+m_SubsampleOffset,outOrigin);
+
+ for (unsigned int i=0 ; i<this->GetNumberOfOutputs() ; i++)
+ {
+ OutputImagePointerType outputPtr = this->GetOutput(i);
+ outputPtr->SetLargestPossibleRegion(outputRegion);
+ outputPtr->SetOrigin(outOrigin);
+ outputPtr->SetSpacing(outSpacing);
+ }
+}
+
template <class TInputImage, class TOutputImage>
void
ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
@@ -209,12 +246,19 @@ ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
// We use only the first output since requested regions for all outputs are enforced to be equal
// by the default GenerateOutputRequestedRegiont() implementation
OutputRegionType outputRequestedRegion = outputPtr->GetRequestedRegion();
+ InputRegionType inputLargest = inputPtr->GetLargestPossibleRegion();
typename OutputRegionType::IndexType outputIndex = outputRequestedRegion.GetIndex();
typename OutputRegionType::SizeType outputSize = outputRequestedRegion.GetSize();
typename InputRegionType::IndexType inputIndex;
typename InputRegionType::SizeType inputSize;
+ // Convert index and size to full grid
+ outputIndex[0] = outputIndex[0] * m_SubsampleFactor[0] + m_SubsampleOffset[0] + inputLargest.GetIndex(0);
+ outputIndex[1] = outputIndex[1] * m_SubsampleFactor[1] + m_SubsampleOffset[1] + inputLargest.GetIndex(1);
+ outputSize[0] = 1 + (outputSize[0] - 1) * m_SubsampleFactor[0];
+ outputSize[1] = 1 + (outputSize[1] - 1) * m_SubsampleFactor[1];
+
// First, apply offset
for (unsigned int dim = 0; dim < InputImageType::ImageDimension; ++dim)
{
@@ -311,6 +355,8 @@ ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
const unsigned int histSize = m_NumberOfBinsPerAxis;
const long unsigned int twiceHistSize = 2 * m_NumberOfBinsPerAxis;
+ InputRegionType inputLargest = inputPtr->GetLargestPossibleRegion();
+
// Set-up progress reporting
itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
@@ -330,11 +376,16 @@ ScalarImageToAdvancedTexturesFilter<TInputImage, TOutputImage>
typename InputRegionType::IndexType inputIndex;
typename InputRegionType::SizeType inputSize;
+ // Convert index to full grid
+ typename OutputImageType::IndexType outIndex;
+
// 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] = varianceIt.GetIndex()[dim] - m_Radius[dim];
+ outIndex[dim] = varianceIt.GetIndex()[dim] * m_SubsampleFactor[dim]
+ + m_SubsampleOffset[dim] + inputLargest.GetIndex(dim);
+ inputIndex[dim] = outIndex[dim] - m_Radius[dim];
inputSize[dim] = 2 * m_Radius[dim] + 1;
}
diff --git a/Modules/Feature/Textures/include/otbScalarImageToHigherOrderTexturesFilter.h b/Modules/Feature/Textures/include/otbScalarImageToHigherOrderTexturesFilter.h
index 4879c7565a..781d69f80b 100644
--- a/Modules/Feature/Textures/include/otbScalarImageToHigherOrderTexturesFilter.h
+++ b/Modules/Feature/Textures/include/otbScalarImageToHigherOrderTexturesFilter.h
@@ -148,6 +148,18 @@ public:
itkSetMacro(FastCalculations, bool);
itkBooleanMacro(FastCalculations);
+ /** Set the sub-sampling factor */
+ itkSetMacro(SubsampleFactor, SizeType);
+
+ /** Get the sub-sampling factor */
+ itkGetMacro(SubsampleFactor, SizeType);
+
+ /** Set the sub-sampling offset */
+ itkSetMacro(SubsampleOffset, OffsetType);
+
+ /** Get the sub-sampling offset */
+ itkGetMacro(SubsampleOffset, OffsetType);
+
/** Get the Short Run Emphasis output image */
OutputImageType * GetShortRunEmphasisOutput();
@@ -186,6 +198,8 @@ protected:
ScalarImageToHigherOrderTexturesFilter();
/** Destructor */
~ScalarImageToHigherOrderTexturesFilter() ITK_OVERRIDE;
+ /** Generate the output informations */
+ void GenerateOutputInformation() ITK_OVERRIDE;
/** Generate the input requested region */
void GenerateInputRequestedRegion() ITK_OVERRIDE;
/** Parallel textures extraction */
@@ -215,6 +229,12 @@ private:
/** Fast calculation */
bool m_FastCalculations;
+
+ /** Sub-sampling factor */
+ SizeType m_SubsampleFactor;
+
+ /** Sub-sampling offset */
+ OffsetType m_SubsampleOffset;
};
} // End namespace otb
diff --git a/Modules/Feature/Textures/include/otbScalarImageToHigherOrderTexturesFilter.txx b/Modules/Feature/Textures/include/otbScalarImageToHigherOrderTexturesFilter.txx
index 66d4a4b65a..fecb7ae6bb 100644
--- a/Modules/Feature/Textures/include/otbScalarImageToHigherOrderTexturesFilter.txx
+++ b/Modules/Feature/Textures/include/otbScalarImageToHigherOrderTexturesFilter.txx
@@ -32,7 +32,9 @@ ScalarImageToHigherOrderTexturesFilter<TInputImage, TOutputImage>
m_NumberOfBinsPerAxis(8),
m_InputImageMinimum(0),
m_InputImageMaximum(255),
- m_FastCalculations(false)
+ m_FastCalculations(false),
+ m_SubsampleFactor(),
+ m_SubsampleOffset()
{
// There are 11 outputs corresponding to the 8 textures indices
this->SetNumberOfRequiredOutputs(10);
@@ -71,6 +73,8 @@ ScalarImageToHigherOrderTexturesFilter<TInputImage, TOutputImage>
}
this->SetOffsets( offsets );
+ this->m_SubsampleFactor.Fill(1);
+ this->m_SubsampleOffset.Fill(0);
}
template <class TInputImage, class TOutputImage>
@@ -231,6 +235,38 @@ ScalarImageToHigherOrderTexturesFilter<TInputImage, TOutputImage>
this->SetOffsets( offsetVector );
}
+template <class TInputImage, class TOutputImage>
+void
+ScalarImageToHigherOrderTexturesFilter<TInputImage, TOutputImage>
+::GenerateOutputInformation()
+{
+ // First, call superclass implementation
+ Superclass::GenerateOutputInformation();
+
+ // Compute output size, origin & spacing
+ InputRegionType inputRegion = this->GetInput()->GetLargestPossibleRegion();
+ OutputRegionType outputRegion;
+ outputRegion.SetIndex(0,0);
+ outputRegion.SetIndex(1,0);
+ outputRegion.SetSize(0, 1 + (inputRegion.GetSize(0) - 1 - m_SubsampleOffset[0]) / m_SubsampleFactor[0]);
+ outputRegion.SetSize(1, 1 + (inputRegion.GetSize(1) - 1 - m_SubsampleOffset[1]) / m_SubsampleFactor[1]);
+
+ typename OutputImageType::SpacingType outSpacing = this->GetInput()->GetSpacing();
+ outSpacing[0] *= m_SubsampleFactor[0];
+ outSpacing[1] *= m_SubsampleFactor[1];
+
+ typename OutputImageType::PointType outOrigin;
+ this->GetInput()->TransformIndexToPhysicalPoint(inputRegion.GetIndex()+m_SubsampleOffset,outOrigin);
+
+ for (unsigned int i=0 ; i<this->GetNumberOfOutputs() ; i++)
+ {
+ OutputImagePointerType outputPtr = this->GetOutput(i);
+ outputPtr->SetLargestPossibleRegion(outputRegion);
+ outputPtr->SetOrigin(outOrigin);
+ outputPtr->SetSpacing(outSpacing);
+ }
+}
+
template <class TInputImage, class TOutputImage>
void
ScalarImageToHigherOrderTexturesFilter<TInputImage, TOutputImage>
@@ -252,7 +288,17 @@ ScalarImageToHigherOrderTexturesFilter<TInputImage, TOutputImage>
// We use only the first output since requested regions for all outputs are enforced to be equal
// by the default GenerateOutputRequestedRegiont() implementation
OutputRegionType outputRequestedRegion = outputPtr->GetRequestedRegion();
- InputRegionType inputRequestedRegion = outputRequestedRegion;
+ typename OutputRegionType::IndexType outputIndex = outputRequestedRegion.GetIndex();
+ typename OutputRegionType::SizeType outputSize = outputRequestedRegion.GetSize();
+ InputRegionType inputLargest = inputPtr->GetLargestPossibleRegion();
+
+ // Convert index and size to full grid
+ outputIndex[0] = outputIndex[0] * m_SubsampleFactor[0] + m_SubsampleOffset[0] + inputLargest.GetIndex(0);
+ outputIndex[1] = outputIndex[1] * m_SubsampleFactor[1] + m_SubsampleOffset[1] + inputLargest.GetIndex(1);
+ outputSize[0] = 1 + (outputSize[0] - 1) * m_SubsampleFactor[0];
+ outputSize[1] = 1 + (outputSize[1] - 1) * m_SubsampleFactor[1];
+
+ InputRegionType inputRequestedRegion(outputIndex,outputSize);
// Apply the radius
inputRequestedRegion.PadByRadius(m_Radius);
@@ -297,6 +343,8 @@ ScalarImageToHigherOrderTexturesFilter<TInputImage, TOutputImage>
inputPtr->TransformIndexToPhysicalPoint( outputImagesIterators[0].GetIndex() + m_Radius, bottomRightPoint );
double maxDistance = topLeftPoint.EuclideanDistanceTo(bottomRightPoint);
+ InputRegionType inputLargest = inputPtr->GetLargestPossibleRegion();
+
// Set-up progress reporting
itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
@@ -304,12 +352,17 @@ ScalarImageToHigherOrderTexturesFilter<TInputImage, TOutputImage>
while ( !outputImagesIterators[0].IsAtEnd() )
{
// Compute the region on which run-length matrix will be estimated
- typename InputRegionType::IndexType inputIndex = outputImagesIterators[0].GetIndex() - m_Radius;
+ typename InputRegionType::IndexType inputIndex;
typename InputRegionType::SizeType inputSize;
- // First, apply offset
+ // Convert index to full grid
+ typename OutputImageType::IndexType outIndex;
+
for (unsigned int dim = 0; dim < InputImageType::ImageDimension; ++dim)
{
+ outIndex[dim] = outputImagesIterators[0].GetIndex()[dim] * m_SubsampleFactor[dim]
+ + m_SubsampleOffset[dim] + inputLargest.GetIndex(dim);
+ inputIndex[dim] = outIndex[dim] - m_Radius[dim];
inputSize[dim] = 2 * m_Radius[dim] + 1;
}
--
GitLab