From 0bd70bd4beefabca6b0223b0ee787e19d330cce2 Mon Sep 17 00:00:00 2001
From: Aurelien Bricier <aurelien.bricier@c-s.fr>
Date: Tue, 12 Oct 2010 17:54:10 +0200
Subject: [PATCH] ENH: added one-pass computation of multiple complex moments
---
.../otbComplexMomentImageFunction.h | 58 ++++---
.../otbComplexMomentImageFunction.txx | 148 ++++++++++--------
.../ComplexMomentImageExample.cxx | 13 +-
Testing/Code/FeatureExtraction/CMakeLists.txt | 4 +-
.../otbComplexMomentImage.cxx | 22 +--
5 files changed, 141 insertions(+), 104 deletions(-)
diff --git a/Code/FeatureExtraction/otbComplexMomentImageFunction.h b/Code/FeatureExtraction/otbComplexMomentImageFunction.h
index 1b086db60a..a3223780e8 100644
--- a/Code/FeatureExtraction/otbComplexMomentImageFunction.h
+++ b/Code/FeatureExtraction/otbComplexMomentImageFunction.h
@@ -18,7 +18,7 @@
#ifndef __otbComplexMomentImageFunction_h
#define __otbComplexMomentImageFunction_h
-#include "otbGeometricMomentImageFunction.h"
+#include "itkImageFunction.h"
#include <complex>
@@ -45,43 +45,48 @@ namespace otb
*
* \ingroup ImageFunctions
*/
-template <class TInput,
- class TOutput = std::complex<double>,
- class TPrecision = double,
- class TCoordRep = float>
+
+template <class TInputImage, class TCoordRep = float>
class ITK_EXPORT ComplexMomentImageFunction :
- public GeometricMomentImageFunction<TInput, TOutput, TPrecision, TCoordRep>
+ public itk::ImageFunction <TInputImage,
+ std::vector< std::vector< std::complex<double> > >,
+ TCoordRep>
{
public:
/** Standard class typedefs. */
typedef ComplexMomentImageFunction Self;
- typedef GeometricMomentImageFunction<TInput, TOutput, TPrecision, TCoordRep> Superclass;
+ typedef itk::ImageFunction<TInputImage,
+ std::vector<
+ std::vector<
+ std::complex<double> > >,
+ TCoordRep> Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
/** Run-time type information (and related methods). */
- itkTypeMacro(ComplexMomentImageFunction, GeometricMomentImageFunction);
+ itkTypeMacro(ComplexMomentImageFunction, ImageFunction);
/** Method for creation through the object factory. */
itkNewMacro(Self);
/** InputImageType typedef support. */
- typedef typename Superclass::InputType InputType;
+ typedef TInputImage InputImageType;
typedef typename Superclass::IndexType IndexType;
typedef typename Superclass::ContinuousIndexType ContinuousIndexType;
typedef typename Superclass::PointType PointType;
- typedef typename Superclass::OutputType ComplexType;
-
- /** Type for calculation precision */
- typedef typename Superclass::PrecisionType PrecisionType;
+ typedef float ScalarRealType;
- /** ComplexType for calculation precision */
- typedef std::complex<PrecisionType> ComplexPrecisionType;
+ typedef typename Superclass::OutputType ComplexType;
+ typedef typename std::complex<ScalarRealType> ScalarComplexType;
+ /** Dimension of the underlying image. */
+ itkStaticConstMacro(ImageDimension, unsigned int,
+ InputImageType::ImageDimension);
+
/** Evalulate the function at specified index */
virtual ComplexType EvaluateAtIndex(const IndexType& index) const;
-
+
/** Evaluate the function at non-integer positions */
virtual ComplexType Evaluate(const PointType& point) const
{
@@ -96,11 +101,17 @@ public:
this->ConvertContinuousIndexToNearestIndex(cindex, index);
return this->EvaluateAtIndex(index);
}
-
- itkSetMacro(P, unsigned int);
- itkGetConstReferenceMacro(P, unsigned int);
- itkSetMacro(Q, unsigned int);
- itkGetConstReferenceMacro(Q, unsigned int);
+
+ /** Get/Set the radius of the neighborhood over which the
+ * statistics are evaluated
+ */
+ itkSetMacro( NeighborhoodRadius, unsigned int );
+ itkGetConstReferenceMacro( NeighborhoodRadius, unsigned int );
+
+ itkSetMacro(Pmax, unsigned int);
+ itkGetConstReferenceMacro(Pmax, unsigned int);
+ itkSetMacro(Qmax, unsigned int);
+ itkGetConstReferenceMacro(Qmax, unsigned int);
protected:
ComplexMomentImageFunction();
@@ -111,8 +122,9 @@ private:
ComplexMomentImageFunction(const Self &); //purposely not implemented
void operator =(const Self&); //purposely not implemented
- unsigned int m_P;
- unsigned int m_Q;
+ unsigned int m_Pmax;
+ unsigned int m_Qmax;
+ unsigned int m_NeighborhoodRadius;
};
diff --git a/Code/FeatureExtraction/otbComplexMomentImageFunction.txx b/Code/FeatureExtraction/otbComplexMomentImageFunction.txx
index d1eeebc2ec..96106d52c0 100644
--- a/Code/FeatureExtraction/otbComplexMomentImageFunction.txx
+++ b/Code/FeatureExtraction/otbComplexMomentImageFunction.txx
@@ -19,11 +19,9 @@
#define __otbComplexMomentImageFunction_txx
#include "otbComplexMomentImageFunction.h"
-#include "itkImageRegionIterator.h"
-#include "itkImage.h"
#include "itkConstNeighborhoodIterator.h"
#include "itkNumericTraits.h"
-#include "otbMacro.h"
+#include "itkMacro.h"
namespace otb
{
@@ -31,96 +29,124 @@ namespace otb
/**
* Constructor
*/
-template <class TInput, class TOutput, class TPrecision, class TCoordRep>
-ComplexMomentImageFunction<TInput, TOutput, TPrecision, TCoordRep>
+template <class TInputImage, class TCoordRep>
+ComplexMomentImageFunction<TInputImage, TCoordRep>
::ComplexMomentImageFunction()
{
- m_P = 0;
- m_Q = 0;
+ m_NeighborhoodRadius = 1;
+ m_Pmax = 4;
+ m_Qmax = 4;
}
/**
*
*/
-template <class TInput, class TOutput, class TPrecision, class TCoordRep>
+template <class TInputImage, class TCoordRep>
void
-ComplexMomentImageFunction<TInput, TOutput, TPrecision, TCoordRep>
+ComplexMomentImageFunction<TInputImage, TCoordRep>
::PrintSelf(std::ostream& os, itk::Indent indent) const
{
this->Superclass::PrintSelf(os, indent);
- os << indent << " p indice value : " << m_P << std::endl;
- os << indent << " q indice value : " << m_Q << std::endl;
+ os << indent << " p indice maximum value : " << m_Pmax << std::endl;
+ os << indent << " q indice maximum value : " << m_Qmax << std::endl;
+ os << indent << " Neighborhood radius value : " << m_NeighborhoodRadius << std::endl;
}
-template <class TInput, class TOutput, class TPrecision, class TCoordRep>
-typename ComplexMomentImageFunction<TInput, TOutput, TPrecision, TCoordRep>::ComplexType
-ComplexMomentImageFunction<TInput, TOutput, TPrecision, TCoordRep>
+template <class TInputImage, class TCoordRep>
+typename ComplexMomentImageFunction<TInputImage, TCoordRep>::ComplexType
+ComplexMomentImageFunction<TInputImage, TCoordRep>
::EvaluateAtIndex(const IndexType& index) const
{
- typename TInput::SizeType ImageSize;
- ComplexPrecisionType Sum;
- ComplexPrecisionType ValP;
- ComplexPrecisionType ValQ;
- PrecisionType Norm;
- IndexType IndexValue;
- IndexType indexPos = index;
- typename TInput::SizeType kernelSize;
-
- if (!this->GetInputImage())
+ // Build moments vector
+ ComplexType moments;//<--- pb VariableLengthVector : compute the length before
+ moments.resize(m_Pmax+1);
+
+ std::vector<ScalarComplexType> valXpY, valXqY;
+ valXpY.resize(m_Pmax+1);
+ valXqY.resize(m_Qmax+1);
+
+ // Initialize moments
+ for (unsigned int p = 0; p <= m_Pmax; p++)
{
- otbMsgDevMacro(<< "Pb with GetInputImage");
- return (ComplexType(itk::NumericTraits<PrecisionType>::Zero, itk::NumericTraits<PrecisionType>::Zero));
+ moments.at(p).resize(m_Qmax+1);
+ valXpY.at(p) = ScalarComplexType(1.0,0.0);
+ for (unsigned int q = 0; q <= m_Qmax; q++)
+ {
+ moments.at(p).at(q) = ScalarComplexType(0.0,0.0);
+ valXqY.at(q) = ScalarComplexType(1.0,0.0);
+ }
}
-
- if (this->GetNeighborhoodRadius() < 0)
+ std::cout << "moments initiated!" << std::endl;
+ // Check for input image
+ if( !this->GetInputImage() )
{
- ImageSize = this->GetInputImage()->GetBufferedRegion().GetSize();
-
- indexPos[0] = ImageSize[0] / 2;
- indexPos[1] = ImageSize[1] / 2;
-
- kernelSize[0] = indexPos[0];
- kernelSize[1] = indexPos[1];
+ return moments;
}
- else
+
+ // Check for out of buffer
+ if ( !this->IsInsideBuffer( index ) )
{
- kernelSize.Fill(this->GetNeighborhoodRadius());
+ return moments;
}
-
- itk::ConstNeighborhoodIterator<TInput>
- it(kernelSize, this->GetInputImage(), this->GetInputImage()->GetBufferedRegion());
-
+
+ // Create an N-d neighborhood kernel, using a zeroflux boundary condition
+ typename InputImageType::SizeType kernelSize;
+ kernelSize.Fill( m_NeighborhoodRadius );
+
+ itk::ConstNeighborhoodIterator<InputImageType>
+ it(kernelSize, this->GetInputImage(), this->GetInputImage()->GetBufferedRegion());
+
// Set the iterator at the desired location
- it.SetLocation(indexPos);
- Sum = ComplexPrecisionType(0.0, 0.0);
- Norm = 0;
-
+ it.SetLocation(index);
+
+ // Walk the neighborhood
const unsigned int size = it.Size();
for (unsigned int i = 0; i < size; ++i)
{
- IndexValue = it.GetIndex(i);
- ValP = ComplexPrecisionType(1.0, 0.0);
- ValQ = ComplexPrecisionType(1.0, 0.0);
- unsigned int p = m_P;
- while (p > 0)
+ // Retrieve value, and centered-reduced position
+ ScalarRealType value = static_cast<ScalarRealType>(it.GetPixel(i));
+ ScalarRealType x = static_cast<ScalarRealType>(it.GetOffset(i)[0]);
+ ScalarRealType y = static_cast<ScalarRealType>(it.GetOffset(i)[1]);
+
+ // Build complex value
+ ScalarComplexType xpy(x,y),xqy(x,-y);
+
+ unsigned int pTmp = 1;
+ unsigned int qTmp = 1;
+
+ while (pTmp <= m_Pmax)
{
- ValP *= ComplexPrecisionType(IndexValue[0] - indexPos[0], IndexValue[1] - indexPos[1]);
- --p;
+ valXpY.at(pTmp) = valXpY.at(pTmp-1) * xpy;
+ pTmp ++;
}
- unsigned int q = m_Q;
- while (q > 0)
+ while (qTmp <= m_Qmax)
{
- ValQ *= ComplexPrecisionType(IndexValue[0] - indexPos[0], -(IndexValue[1] - indexPos[1]));
- --q;
+ valXqY.at(qTmp) = valXqY.at(qTmp-1) * xqy;
+ qTmp ++;
}
+
- Sum += (ValP * ValQ * ComplexPrecisionType(static_cast<PrecisionType>(it.GetPixel(i)), 0.0));
- Norm += it.GetPixel(i);
+ // Update cumulants
+ for (unsigned int p = 0; p <= m_Pmax; p++)
+ {
+ for (unsigned int q= 0; q <= m_Qmax; q++)
+ {
+ moments.at(p).at(q) += valXpY.at(p) * valXqY.at(q) * value;
+ }
+ }
+ }
+
+ // Normalisation
+ for (unsigned int p = 0; p <= m_Pmax; p++)
+ {
+ for (unsigned int q= 0; q <= m_Qmax; q++)
+ {
+ moments.at(p).at(q) /= vcl_pow(moments.at(0).at(0), (p+q)/2);
+ }
}
- Norm = vcl_pow(Norm, ((PrecisionType) m_P + (PrecisionType) m_Q) / 2.);
-
- return (static_cast<ComplexType>(Sum / Norm));
+ // Return result
+ return moments;
}
} // namespace otb
diff --git a/Examples/FeatureExtraction/ComplexMomentImageExample.cxx b/Examples/FeatureExtraction/ComplexMomentImageExample.cxx
index aedbf9bf9e..1ba7425254 100644
--- a/Examples/FeatureExtraction/ComplexMomentImageExample.cxx
+++ b/Examples/FeatureExtraction/ComplexMomentImageExample.cxx
@@ -75,10 +75,9 @@ int main(int argc, char * argv[])
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
- typedef std::complex<float>
- ComplexType;
- typedef otb::ComplexMomentImageFunction<InputImageType, ComplexType> CMType;
-
+ typedef otb::ComplexMomentImageFunction<InputImageType> CMType;
+ typedef CMType::ComplexType ComplexType;
+
CMType::Pointer cmFunction = CMType::New();
// Software Guide : EndCodeSnippet
@@ -92,8 +91,8 @@ int main(int argc, char * argv[])
// Software Guide : BeginCodeSnippet
reader->Update();
cmFunction->SetInputImage(reader->GetOutput());
- cmFunction->SetQ(Q);
- cmFunction->SetP(P);
+ cmFunction->SetQmax(Q);
+ cmFunction->SetPmax(P);
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
@@ -126,7 +125,7 @@ int main(int argc, char * argv[])
ComplexType Result = cmFunction->EvaluateAtIndex(center);
std::cout << "The moment of order (" << P << "," << Q <<
- ") is equal to " << Result << std::endl;
+ ") is equal to " << Result.at(P).at(Q) << std::endl;
// Software Guide : EndCodeSnippet
return EXIT_SUCCESS;
diff --git a/Testing/Code/FeatureExtraction/CMakeLists.txt b/Testing/Code/FeatureExtraction/CMakeLists.txt
index 30506970fc..e95f12a706 100644
--- a/Testing/Code/FeatureExtraction/CMakeLists.txt
+++ b/Testing/Code/FeatureExtraction/CMakeLists.txt
@@ -69,8 +69,8 @@ ADD_TEST(feTvComplexMomentImage ${FEATUREEXTRACTION_TESTS1}
${TEMP}/feComplexMomentImage.txt
otbComplexMomentImage
${INPUTDATA}/TeteAToto.png
- 1
- 1
+ 4
+ 4
${TEMP}/feComplexMomentImage.txt)
ADD_TEST(feTvHuImage ${FEATUREEXTRACTION_TESTS1}
diff --git a/Testing/Code/FeatureExtraction/otbComplexMomentImage.cxx b/Testing/Code/FeatureExtraction/otbComplexMomentImage.cxx
index 2e37798187..3124331437 100644
--- a/Testing/Code/FeatureExtraction/otbComplexMomentImage.cxx
+++ b/Testing/Code/FeatureExtraction/otbComplexMomentImage.cxx
@@ -42,8 +42,8 @@ int otbComplexMomentImage(int argc, char * argv[])
typedef itk::Image<InputPixelType, Dimension> InputImageType;
typedef otb::ImageFileReader<InputImageType> ReaderType;
- typedef std::complex<float> ComplexType;
- typedef otb::ComplexMomentImageFunction<InputImageType, ComplexType> CMType;
+ typedef otb::ComplexMomentImageFunction<InputImageType> CMType;
+ typedef CMType::ComplexType ComplexType;
ReaderType::Pointer reader = ReaderType::New();
CMType::Pointer function = CMType::New();
@@ -53,8 +53,8 @@ int otbComplexMomentImage(int argc, char * argv[])
reader->Update();
function->SetInputImage(reader->GetOutput());
- function->SetQ(q);
- function->SetP(p);
+ function->SetQmax(q);
+ function->SetPmax(p);
InputImageType::IndexType index;
index[0] = 10;
@@ -64,15 +64,15 @@ int otbComplexMomentImage(int argc, char * argv[])
std::ofstream outputStream(outputFilename);
- Result = function->EvaluateAtIndex(index);
- outputStream << std::setprecision(10) << "function->EvaluateAtIndex( index ): " << Result << std::endl;
-
- outputStream << " With NeighborhoodRadius(3):" << std::endl;
-
function->SetNeighborhoodRadius(3);
Result = function->EvaluateAtIndex(index);
- outputStream << "function->EvaluateAtIndex( index ): " << Result << std::endl;
-
+ for (unsigned int k=0; k<=p; k++)
+ {
+ for (unsigned int l=0; l<=q; l++)
+ {
+ outputStream << "ComplexMoment c(" << k << l << ") : " << Result.at(k).at(l) << std::endl;
+ }
+ }
outputStream.close();
return EXIT_SUCCESS;
--
GitLab