diff --git a/Testing/Code/FeatureExtraction/CMakeLists.txt b/Testing/Code/FeatureExtraction/CMakeLists.txt
index 4b7c1899353fc7e322a3863d9dcea6bf70da44a5..18e402a7826a06eee1f124389c0dcd4063fcccc8 100644
--- a/Testing/Code/FeatureExtraction/CMakeLists.txt
+++ b/Testing/Code/FeatureExtraction/CMakeLists.txt
@@ -73,6 +73,16 @@ ADD_TEST(feTvComplexMomentImage ${FEATUREEXTRACTION_TESTS1}
4
${TEMP}/feComplexMomentImage.txt)
+ADD_TEST(feTvComplexMomentImageScaleInvariant ${FEATUREEXTRACTION_TESTS1}
+ --compare-ascii ${NOTOL} ${BASELINE_FILES}/feComplexMomentImageScaleInvariant.txt
+ ${TEMP}/feComplexMomentImageScaleInvariant.txt
+ otbComplexMomentImageScaleInvariant
+ ${INPUTDATA}/poupees.png
+ 4
+ 4
+ ${TEMP}/feComplexMomentImageScaleInvariant.txt
+ ${TEMP}/feComplexImageScaleInvariant.tiff)
+
ADD_TEST(feTvRealMomentImage ${FEATUREEXTRACTION_TESTS1}
--compare-ascii ${NOTOL} ${BASELINE_FILES}/feRealMomentImage.txt
${TEMP}/feRealMomentImage.txt
@@ -96,6 +106,21 @@ ADD_TEST(feTvHuImage ${FEATUREEXTRACTION_TESTS1}
${INPUTDATA}/poupees.png
${TEMP}/feHuImage.txt)
+ADD_TEST(feTuHuImageScaleInvariant ${FEATUREEXTRACTION_TESTS1}
+ otbHuImageScaleInvariant
+ ${INPUTDATA}/poupees.png
+ ${TEMP}/feHuImageScaleInvariant.txt
+ ${TEMP}/feHuImageScaleInvariant.tiff
+ )
+
+ADD_TEST(feTuHuImageRotationInvariant ${FEATUREEXTRACTION_TESTS1}
+ otbHuImageRotationInvariant
+ ${INPUTDATA}/poupees.png
+ ${TEMP}/feHuImageRotationInvariant.txt
+ ${TEMP}/feHuImageRotationInvariant.tiff
+ 90
+ )
+
ADD_TEST(feTuFlusserImage ${FEATUREEXTRACTION_TESTS1}
--compare-ascii ${EPSILON_8} ${BASELINE_FILES}/feFlusserImage.txt
${TEMP}/feFlusserImage.txt
@@ -103,6 +128,21 @@ ADD_TEST(feTuFlusserImage ${FEATUREEXTRACTION_TESTS1}
${INPUTDATA}/poupees.png
${TEMP}/feFlusserImage.txt)
+ADD_TEST(feTuFlusserImageScaleInvariant ${FEATUREEXTRACTION_TESTS1}
+ otbFlusserImageScaleInvariant
+ ${INPUTDATA}/poupees.png
+ ${TEMP}/feFlusserImageScaleInvariant.txt
+ ${TEMP}/feFlusserImageScaleInvariant.tiff
+ )
+
+ADD_TEST(feTuFlusserImageRotationInvariant ${FEATUREEXTRACTION_TESTS1}
+ otbFlusserImageRotationInvariant
+ ${INPUTDATA}/poupees.png
+ ${TEMP}/feFlusserImageRotationInvariant.txt
+ ${TEMP}/feFlusserImageRotationInvariant.tiff
+ 90
+ )
+
ADD_TEST(feTuComplexMomentPathNew ${FEATUREEXTRACTION_TESTS1}
otbComplexMomentPathNew)
diff --git a/Testing/Code/FeatureExtraction/otbComplexMomentImage.cxx b/Testing/Code/FeatureExtraction/otbComplexMomentImage.cxx
index 4b1c5cc5c645fdabe26896c225bbd676c243c979..02add720eb6ae24563cd4959af056621b8922a36 100644
--- a/Testing/Code/FeatureExtraction/otbComplexMomentImage.cxx
+++ b/Testing/Code/FeatureExtraction/otbComplexMomentImage.cxx
@@ -28,6 +28,11 @@
#include "otbImageFileReader.h"
#include "otbComplexMomentImageFunction.h"
+#include "otbBCOInterpolateImageFunction.h"
+#include "itkLinearInterpolateImageFunction.h"
+#include "otbStreamingResampleImageFilter.h"
+#include "otbStreamingImageFileWriter.h"
+#include "itkResampleImageFilter.h"
int otbComplexMomentImage(int argc, char * argv[])
{
@@ -78,3 +83,120 @@ int otbComplexMomentImage(int argc, char * argv[])
return EXIT_SUCCESS;
}
+
+
+int otbComplexMomentImageScaleInvariant(int argc, char * argv[])
+{
+ const char * inputFilename = argv[1];
+ unsigned int p((unsigned int) ::atoi(argv[2]));
+ unsigned int q((unsigned int) ::atoi(argv[3]));
+ const char * outputFilename = argv[4];
+ const char * outputImageName = argv[5];
+
+ typedef double InputPixelType;
+ const unsigned int Dimension = 2;
+ typedef otb::Image<InputPixelType, Dimension> InputImageType;
+ typedef otb::ImageFileReader<InputImageType> ReaderType;
+ typedef otb::StreamingResampleImageFilter<InputImageType,
+ InputImageType,
+ double> StreamingResampleImageFilterType;
+ typedef otb::BCOInterpolateImageFunction<InputImageType,
+ double> InterpolatorType;
+ //typedef itk::LinearInterpolateImageFunction<InputImageType,
+ // double> InterpolatorType;
+ typedef otb::ComplexMomentImageFunction<InputImageType> FunctionType;
+ typedef FunctionType::ComplexType ComplexType;
+ typedef otb::StreamingImageFileWriter<InputImageType> WriterType;
+
+ ReaderType::Pointer reader = ReaderType::New();
+ StreamingResampleImageFilterType::Pointer resampler = StreamingResampleImageFilterType::New();
+ InterpolatorType::Pointer interpolator = InterpolatorType::New();
+ FunctionType::Pointer function1 = FunctionType::New();
+ FunctionType::Pointer function2 = FunctionType::New();
+ WriterType::Pointer writer = WriterType::New();
+
+ reader->SetFileName(inputFilename);
+ reader->Update();
+
+ interpolator->SetInputImage(reader->GetOutput());
+ interpolator->SetRadius(2);
+ interpolator->SetAlpha(-0.5);
+
+ resampler->SetInput(reader->GetOutput());
+ resampler->SetInterpolator(interpolator);
+ StreamingResampleImageFilterType::SizeType size;
+ size[0] = 1024;
+ size[1] = 1024;
+ resampler->SetOutputSize(size);
+ resampler->SetOutputSpacing(0.5);
+ resampler->Update();
+ writer->SetInput(resampler->GetOutput());
+ writer->SetFileName(outputImageName);
+ writer->Update();
+
+ function1->SetInputImage(reader->GetOutput());
+ function1->SetQmax(q);
+ function1->SetPmax(p);
+ InputImageType::IndexType index1;
+ index1[0] = 100;
+ index1[1] = 100;
+ function1->SetNeighborhoodRadius(4);
+ ComplexType Result1 = function1->EvaluateAtIndex(index1);
+
+ function2->SetInputImage(resampler->GetOutput());
+ function2->SetQmax(q);
+ function2->SetPmax(p);
+ InputImageType::IndexType index2;
+ index2[0] = 200;
+ index2[1] = 200;
+ function2->SetNeighborhoodRadius(8);
+ ComplexType Result2 = function2->EvaluateAtIndex(index2);
+
+ std::ofstream outputStream(outputFilename);
+ outputStream << std::setprecision(10);
+
+ for (unsigned int k=0; k<=p; k++)
+ {
+ for (unsigned int l=0; l<=q; l++)
+ {
+ double error;
+ if ( (k+l) >= 2)
+ {
+ if (k == l)
+ {
+ error = 2*vcl_abs( Result1.at(k).at(l).real() - Result2.at(k).at(l).real())
+ / vcl_abs(Result1.at(k).at(l).real() + Result2.at(k).at(l).real());
+ }
+ else
+ {
+ error = (vcl_abs( Result1.at(k).at(l).real() - Result2.at(k).at(l).real())
+ / vcl_abs(Result1.at(k).at(l).real() + Result2.at(k).at(l).real())
+ + vcl_abs( Result1.at(k).at(l).imag() - Result2.at(k).at(l).imag())
+ / vcl_abs(Result1.at(k).at(l).imag() + Result2.at(k).at(l).imag()));
+ }
+
+ outputStream << "Error = " << error << std::endl
+ << "Original Image - Complex Moment c" << k << l
+ << " : " << Result1.at(k).at(l) << std::endl
+ << "Scaled Image - Complex Moment c" << k << l
+ << " : " << Result2.at(k).at(l) << std::endl
+ << std::endl;
+ /*
+ if (error > 1)
+ {
+ itkGenericExceptionMacro( <<std::endl
+ << "Error = " << error
+ << " > 1 -> TEST FAILLED" << std::endl
+ << "Original Image - Complex Moment c" << k << l
+ << " : " << Result1.at(k).at(l) << std::endl
+ << " Scaled Image - Complex Moment c" << k << l
+ << " : " << Result2.at(k).at(l) << std::endl );
+ }
+ */
+ }
+ }
+ }
+ outputStream.close();
+
+ return EXIT_SUCCESS;
+}
diff --git a/Testing/Code/FeatureExtraction/otbFeatureExtractionTests1.cxx b/Testing/Code/FeatureExtraction/otbFeatureExtractionTests1.cxx
index 1e5bb534e5280c3723dbcda504b73156178acd4e..9ad43cad066a42455c94c57f7288c3af2a95f9dd 100644
--- a/Testing/Code/FeatureExtraction/otbFeatureExtractionTests1.cxx
+++ b/Testing/Code/FeatureExtraction/otbFeatureExtractionTests1.cxx
@@ -30,10 +30,15 @@ void RegisterTests()
REGISTER_TEST(otbDrawPathDessinCarre);
REGISTER_TEST(otbDrawPathAlign);
REGISTER_TEST(otbComplexMomentImage);
+ REGISTER_TEST(otbComplexMomentImageScaleInvariant);
REGISTER_TEST(otbRealMomentImage);
REGISTER_TEST(otbRadiometricMomentsImage);
REGISTER_TEST(otbHuImage);
+ REGISTER_TEST(otbHuImageRotationInvariant);
+ REGISTER_TEST(otbHuImageScaleInvariant);
REGISTER_TEST(otbFlusserImage);
+ REGISTER_TEST(otbFlusserImageRotationInvariant);
+ REGISTER_TEST(otbFlusserImageScaleInvariant);
REGISTER_TEST(otbComplexMomentPathNew);
REGISTER_TEST(otbComplexMomentPath);
REGISTER_TEST(otbComplexMomentPathFloat);
diff --git a/Testing/Code/FeatureExtraction/otbFlusserImage.cxx b/Testing/Code/FeatureExtraction/otbFlusserImage.cxx
index 1b7e92c3de92133f7a4d39df4101a7873131a12a..4c8527a23b3223086592d32b3b64823f2a765b0c 100644
--- a/Testing/Code/FeatureExtraction/otbFlusserImage.cxx
+++ b/Testing/Code/FeatureExtraction/otbFlusserImage.cxx
@@ -24,16 +24,23 @@
#include <iomanip>
#include <fstream>
#include "itkExceptionObject.h"
-#include "itkImage.h"
+#include "otbImage.h"
#include "otbImageFileReader.h"
#include "otbFlusserImageFunction.h"
+#include "otbBCOInterpolateImageFunction.h"
+#include "itkLinearInterpolateImageFunction.h"
+#include "otbStreamingResampleImageFilter.h"
+#include "otbStreamingImageFileWriter.h"
+#include "itkResampleImageFilter.h"
+
int otbFlusserImage(int argc, char * argv[])
{
const char * inputFilename = argv[1];
const char * outputFilename = argv[2];
+
typedef unsigned char InputPixelType;
const unsigned int Dimension = 2;
@@ -69,3 +76,202 @@ int otbFlusserImage(int argc, char * argv[])
return EXIT_SUCCESS;
}
+
+int otbFlusserImageScaleInvariant(int argc, char * argv[])
+{
+ const char * inputFilename = argv[1];
+ const char * outputFilename = argv[2];
+ const char * outputImageName = argv[3];
+
+
+ typedef double InputPixelType;
+ const unsigned int Dimension = 2;
+ typedef otb::Image<InputPixelType, Dimension> InputImageType;
+ typedef otb::ImageFileReader<InputImageType> ReaderType;
+ typedef otb::StreamingResampleImageFilter<InputImageType,
+ InputImageType,
+ double> StreamingResampleImageFilterType;
+ typedef otb::BCOInterpolateImageFunction<InputImageType,
+ double> InterpolatorType;
+ typedef otb::FlusserImageFunction<InputImageType> FunctionType;
+ typedef FunctionType::RealType RealType;
+ typedef otb::StreamingImageFileWriter<InputImageType> WriterType;
+
+ ReaderType::Pointer reader = ReaderType::New();
+ StreamingResampleImageFilterType::Pointer resampler = StreamingResampleImageFilterType::New();
+ InterpolatorType::Pointer interpolator = InterpolatorType::New();
+ FunctionType::Pointer function1 = FunctionType::New();
+ FunctionType::Pointer function2 = FunctionType::New();
+ WriterType::Pointer writer = WriterType::New();
+
+ reader->SetFileName(inputFilename);
+ reader->Update();
+
+ interpolator->SetInputImage(reader->GetOutput());
+ interpolator->SetRadius(2);
+ interpolator->SetAlpha(-0.5);
+
+ resampler->SetInput(reader->GetOutput());
+ resampler->SetInterpolator(interpolator);
+ StreamingResampleImageFilterType::SizeType size;
+ size[0] = 1024;
+ size[1] = 1024;
+ resampler->SetOutputSize(size);
+ resampler->SetOutputSpacing(0.5);
+ resampler->Update();
+ writer->SetInput(resampler->GetOutput());
+ writer->SetFileName(outputImageName);
+ writer->Update();
+
+ function1->SetInputImage(reader->GetOutput());
+ InputImageType::IndexType index1;
+ index1[0] = 100;
+ index1[1] = 100;
+ function1->SetNeighborhoodRadius(2);
+ RealType Result1 = function1->EvaluateAtIndex(index1);
+
+ function2->SetInputImage(resampler->GetOutput());
+ InputImageType::IndexType index2;
+ index2[0] = 200;
+ index2[1] = 200;
+ function2->SetNeighborhoodRadius(4);
+ RealType Result2 = function2->EvaluateAtIndex(index2);
+
+ std::ofstream outputStream(outputFilename);
+ outputStream << std::setprecision(10);
+
+ for (unsigned int j = 1; j < 12; j++)
+ {
+ double error = 2*vcl_abs( Result1[j-1] - Result2[j-1]) / vcl_abs(Result1[j-1] + Result2[j-1]);
+
+ outputStream << "Error = " << error << std::endl
+ << "Original Image - Flusser Moment #" << j << " : " << Result1[j-1] << std::endl
+ << "Scaled Image - Flusser Moment #" << j << " : " << Result2[j-1] << std::endl
+ << std::endl;
+
+ if (error > 1)
+ {
+ itkGenericExceptionMacro( <<std::endl
+ << "Error = " << error
+ << " > 1 -> TEST FAILLED" << std::endl
+ << "Original Image - Flusser Moment #" << j << " : " << Result1[j-1] << std::endl
+ << "Scaled Image - Flusser Moment #" << j << " : " << Result2[j-1] << std::endl
+ << std::endl;)
+ }
+ }
+ outputStream.close();
+
+ return EXIT_SUCCESS;
+}
+
+int otbFlusserImageRotationInvariant(int argc, char * argv[])
+{
+ const char * inputFilename = argv[1];
+ const char * outputFilename = argv[2];
+ const char * outputImageName = argv[3];
+ const double angleInDegrees = atoi(argv[4]);
+
+
+ typedef double InputPixelType;
+ const unsigned int Dimension = 2;
+ typedef otb::Image<InputPixelType, Dimension> InputImageType;
+ typedef otb::ImageFileReader<InputImageType> ReaderType;
+ typedef itk::ResampleImageFilter<
+ InputImageType, InputImageType > FilterType;
+ typedef itk::LinearInterpolateImageFunction<InputImageType, double> InterpolatorType;
+ typedef otb::FlusserImageFunction<InputImageType> FunctionType;
+ typedef FunctionType::RealType RealType;
+ typedef otb::StreamingImageFileWriter<InputImageType> WriterType;
+ typedef itk::AffineTransform< double, Dimension > TransformType;
+
+
+ ReaderType::Pointer reader = ReaderType::New();
+ FilterType::Pointer filter = FilterType::New();
+ TransformType::Pointer transform = TransformType::New();
+ InterpolatorType::Pointer interpolator = InterpolatorType::New();
+ FunctionType::Pointer function1 = FunctionType::New();
+ FunctionType::Pointer function2 = FunctionType::New();
+ WriterType::Pointer writer = WriterType::New();
+
+ reader->SetFileName(inputFilename);
+ reader->Update();
+
+ filter->SetInterpolator(interpolator);
+ filter->SetDefaultPixelValue( 100 );
+
+ const InputImageType::SpacingType & spacing = reader->GetOutput()->GetSpacing();
+ const InputImageType::PointType & origin = reader->GetOutput()->GetOrigin();
+ InputImageType::SizeType size =
+ reader->GetOutput()->GetLargestPossibleRegion().GetSize();
+
+ filter->SetOutputOrigin( origin );
+ filter->SetOutputSpacing( spacing );
+ filter->SetOutputDirection( reader->GetOutput()->GetDirection() );
+ filter->SetSize( size );
+
+
+ filter->SetInput(reader->GetOutput());
+
+ TransformType::OutputVectorType translation1;
+ const double imageCenterX = origin[0] + spacing[0] * size[0] / 2.0;
+ const double imageCenterY = origin[1] + spacing[1] * size[1] / 2.0;
+ translation1[0] = -imageCenterX;
+ translation1[1] = -imageCenterY;
+ transform->Translate( translation1 );
+
+ const double degreesToRadians = vcl_atan(1.0) / 45.0;
+ const double angle = angleInDegrees * degreesToRadians;
+ transform->Rotate2D( -angle, false );
+
+ TransformType::OutputVectorType translation2;
+ translation2[0] = imageCenterX;
+ translation2[1] = imageCenterY;
+ transform->Translate( translation2, false );
+
+ filter->SetTransform( transform );
+ filter->Update();
+
+ writer->SetInput(filter->GetOutput());
+ writer->SetFileName(outputImageName);
+ writer->Update();
+
+ function1->SetInputImage(reader->GetOutput());
+ InputImageType::IndexType index1;
+ index1[0] = 256;
+ index1[1] = 256;
+ function1->SetNeighborhoodRadius(4);
+ RealType Result1 = function1->EvaluateAtIndex(index1);
+
+ function2->SetInputImage(filter->GetOutput());
+ InputImageType::IndexType index2;
+ index2[0] = 256;
+ index2[1] = 256;
+ function2->SetNeighborhoodRadius(4);
+ RealType Result2 = function2->EvaluateAtIndex(index2);
+
+ std::ofstream outputStream(outputFilename);
+ outputStream << std::setprecision(10);
+
+ for (unsigned int j = 1; j < 12; j++)
+ {
+ double error = vcl_abs( Result1[j-1] - Result2[j-1]);
+
+ outputStream << "Error = " << error << std::endl
+ << "Original Image - Flusser Moment #" << j << " : " << Result1[j-1] << std::endl
+ << "Rotated Image - Flusser Moment #" << j << " : " << Result2[j-1] << std::endl
+ << "Rotation angle : " << angleInDegrees << std::endl;
+
+ if (error > 1E-9)
+ {
+ itkGenericExceptionMacro( <<std::endl
+ << "Error = " << error
+ << " > 1E-9 -> TEST FAILLED" << std::endl
+ << "Original Image - Flusser Moment #" << j << " : " << Result1[j-1] << std::endl
+ << " Rotated Image - Flusser Moment #" << j << " : " << Result2[j-1] << std::endl
+ << "Rotation angle : " << angleInDegrees << std::endl );
+ }
+ }
+ outputStream.close();
+
+ return EXIT_SUCCESS;
+}
diff --git a/Testing/Code/FeatureExtraction/otbHuImage.cxx b/Testing/Code/FeatureExtraction/otbHuImage.cxx
index 2dbfaaad1fdd1ea29cf8b3623e32ec7415013623..acfd85ab5a4f9127ee60467202a65ea436cc72a3 100644
--- a/Testing/Code/FeatureExtraction/otbHuImage.cxx
+++ b/Testing/Code/FeatureExtraction/otbHuImage.cxx
@@ -24,10 +24,17 @@
#include <iomanip>
#include <fstream>
#include "itkExceptionObject.h"
-#include "itkImage.h"
+#include "otbImage.h"
#include "otbImageFileReader.h"
#include "otbHuImageFunction.h"
+#include "otbBCOInterpolateImageFunction.h"
+#include "itkLinearInterpolateImageFunction.h"
+#include "otbStreamingResampleImageFilter.h"
+#include "otbStreamingImageFileWriter.h"
+#include "itkResampleImageFilter.h"
+
+
int otbHuImage(int argc, char * argv[])
{
@@ -37,7 +44,7 @@ int otbHuImage(int argc, char * argv[])
typedef unsigned char InputPixelType;
const unsigned int Dimension = 2;
- typedef itk::Image<InputPixelType, Dimension> InputImageType;
+ typedef otb::Image<InputPixelType, Dimension> InputImageType;
typedef otb::ImageFileReader<InputImageType> ReaderType;
typedef otb::HuImageFunction<InputImageType> FunctionType;
typedef FunctionType::RealType RealType;
@@ -67,3 +74,201 @@ int otbHuImage(int argc, char * argv[])
return EXIT_SUCCESS;
}
+
+int otbHuImageScaleInvariant(int argc, char * argv[])
+{
+ const char * inputFilename = argv[1];
+ const char * outputFilename = argv[2];
+ const char * outputImageName = argv[3];
+
+
+ typedef double InputPixelType;
+ const unsigned int Dimension = 2;
+ typedef otb::Image<InputPixelType, Dimension> InputImageType;
+ typedef otb::ImageFileReader<InputImageType> ReaderType;
+ typedef otb::StreamingResampleImageFilter<InputImageType,
+ InputImageType,
+ double> StreamingResampleImageFilterType;
+ typedef otb::BCOInterpolateImageFunction<InputImageType,
+ double> InterpolatorType;
+ typedef otb::HuImageFunction<InputImageType> FunctionType;
+ typedef FunctionType::RealType RealType;
+ typedef otb::StreamingImageFileWriter<InputImageType> WriterType;
+
+ ReaderType::Pointer reader = ReaderType::New();
+ StreamingResampleImageFilterType::Pointer resampler = StreamingResampleImageFilterType::New();
+ InterpolatorType::Pointer interpolator = InterpolatorType::New();
+ FunctionType::Pointer function1 = FunctionType::New();
+ FunctionType::Pointer function2 = FunctionType::New();
+ WriterType::Pointer writer = WriterType::New();
+
+ reader->SetFileName(inputFilename);
+ reader->Update();
+
+ interpolator->SetInputImage(reader->GetOutput());
+ interpolator->SetRadius(2);
+ interpolator->SetAlpha(-0.5);
+
+ resampler->SetInput(reader->GetOutput());
+ resampler->SetInterpolator(interpolator);
+ StreamingResampleImageFilterType::SizeType size;
+ size[0] = 1024;
+ size[1] = 1024;
+ resampler->SetOutputSize(size);
+ resampler->SetOutputSpacing(0.5);
+ resampler->Update();
+ writer->SetInput(resampler->GetOutput());
+ writer->SetFileName(outputImageName);
+ writer->Update();
+
+ function1->SetInputImage(reader->GetOutput());
+ InputImageType::IndexType index1;
+ index1[0] = 100;
+ index1[1] = 100;
+ function1->SetNeighborhoodRadius(2);
+ RealType Result1 = function1->EvaluateAtIndex(index1);
+
+ function2->SetInputImage(resampler->GetOutput());
+ InputImageType::IndexType index2;
+ index2[0] = 200;
+ index2[1] = 200;
+ function2->SetNeighborhoodRadius(4);
+ RealType Result2 = function2->EvaluateAtIndex(index2);
+
+ std::ofstream outputStream(outputFilename);
+ outputStream << std::setprecision(10);
+
+ for (unsigned int j = 1; j < 8; j++)
+ {
+ double error = 2*vcl_abs( Result1[j-1] - Result2[j-1]) / vcl_abs(Result1[j-1] + Result2[j-1]);
+
+ outputStream << "Error = " << error << std::endl
+ << "Original Image - Hu Moment #" << j << " : " << Result1[j-1] << std::endl
+ << "Scaled Image - Hu Moment #" << j << " : " << Result2[j-1] << std::endl
+ << std::endl;
+
+ if (error > 1)
+ {
+ itkGenericExceptionMacro( <<std::endl
+ << "Error = " << error
+ << " > 1 -> TEST FAILLED" << std::endl
+ << "Original Image - Hu Moment #" << j << " : " << Result1[j-1] << std::endl
+ << "Scaled Image - Hu Moment #" << j << " : " << Result2[j-1] << std::endl
+ << std::endl);
+ }
+ }
+ outputStream.close();
+
+ return EXIT_SUCCESS;
+}
+
+int otbHuImageRotationInvariant(int argc, char * argv[])
+{
+ const char * inputFilename = argv[1];
+ const char * outputFilename = argv[2];
+ const char * outputImageName = argv[3];
+ const double angleInDegrees = atoi(argv[4]);
+
+
+ typedef double InputPixelType;
+ const unsigned int Dimension = 2;
+ typedef otb::Image<InputPixelType, Dimension> InputImageType;
+ typedef otb::ImageFileReader<InputImageType> ReaderType;
+ typedef itk::ResampleImageFilter<
+ InputImageType, InputImageType > FilterType;
+ typedef itk::LinearInterpolateImageFunction<InputImageType, double> InterpolatorType;
+ typedef otb::HuImageFunction<InputImageType> FunctionType;
+ typedef FunctionType::RealType RealType;
+ typedef otb::StreamingImageFileWriter<InputImageType> WriterType;
+ typedef itk::AffineTransform< double, Dimension > TransformType;
+
+
+ ReaderType::Pointer reader = ReaderType::New();
+ FilterType::Pointer filter = FilterType::New();
+ TransformType::Pointer transform = TransformType::New();
+ InterpolatorType::Pointer interpolator = InterpolatorType::New();
+ FunctionType::Pointer function1 = FunctionType::New();
+ FunctionType::Pointer function2 = FunctionType::New();
+ WriterType::Pointer writer = WriterType::New();
+
+ reader->SetFileName(inputFilename);
+ reader->Update();
+
+ filter->SetInterpolator(interpolator);
+ filter->SetDefaultPixelValue( 100 );
+
+ const InputImageType::SpacingType & spacing = reader->GetOutput()->GetSpacing();
+ const InputImageType::PointType & origin = reader->GetOutput()->GetOrigin();
+ InputImageType::SizeType size =
+ reader->GetOutput()->GetLargestPossibleRegion().GetSize();
+
+ filter->SetOutputOrigin( origin );
+ filter->SetOutputSpacing( spacing );
+ filter->SetOutputDirection( reader->GetOutput()->GetDirection() );
+ filter->SetSize( size );
+
+ filter->SetInput(reader->GetOutput());
+
+ TransformType::OutputVectorType translation1;
+ const double imageCenterX = origin[0] + spacing[0] * size[0] / 2.0;
+ const double imageCenterY = origin[1] + spacing[1] * size[1] / 2.0;
+ translation1[0] = -imageCenterX;
+ translation1[1] = -imageCenterY;
+ transform->Translate( translation1 );
+
+ const double degreesToRadians = vcl_atan(1.0) / 45.0;
+ const double angle = angleInDegrees * degreesToRadians;
+ transform->Rotate2D( -angle, false );
+
+ TransformType::OutputVectorType translation2;
+ translation2[0] = imageCenterX;
+ translation2[1] = imageCenterY;
+ transform->Translate( translation2, false );
+
+ filter->SetTransform( transform );
+ filter->Update();
+
+ writer->SetInput(filter->GetOutput());
+ writer->SetFileName(outputImageName);
+ writer->Update();
+
+ function1->SetInputImage(reader->GetOutput());
+ InputImageType::IndexType index1;
+ index1[0] = 256;
+ index1[1] = 256;
+ function1->SetNeighborhoodRadius(4);
+ RealType Result1 = function1->EvaluateAtIndex(index1);
+
+ function2->SetInputImage(filter->GetOutput());
+ InputImageType::IndexType index2;
+ index2[0] = 256;
+ index2[1] = 256;
+ function2->SetNeighborhoodRadius(4);
+ RealType Result2 = function2->EvaluateAtIndex(index2);
+
+ std::ofstream outputStream(outputFilename);
+ outputStream << std::setprecision(10);
+
+ for (unsigned int j = 1; j < 8; j++)
+ {
+ double error = vcl_abs( Result1[j-1] - Result2[j-1]);
+
+ outputStream << "Error = " << error << std::endl
+ << "Original Image - Hu Moment #" << j << " : " << Result1[j-1] << std::endl
+ << "Rotated Image - Hu Moment #" << j << " : " << Result2[j-1] << std::endl
+ << "Rotation angle : " << angleInDegrees << std::endl;
+
+ if (error > 1E-9)
+ {
+ itkGenericExceptionMacro( <<std::endl
+ << "Error = " << error
+ << " > 1E-9 -> TEST FAILLED" << std::endl
+ << "Original Image - Hu Moment #" << j << " : " << Result1[j-1] << std::endl
+ << "Rotated Image - Hu Moment #" << j << " : " << Result2[j-1] << std::endl
+ << "Rotation angle : " << angleInDegrees << std::endl );
+ }
+ }
+ outputStream.close();
+
+ return EXIT_SUCCESS;
+}