From 174ca44afb935da5ba4aa1bc840e831dd70c636a Mon Sep 17 00:00:00 2001 From: Julien Michel <julien.michel@c-s.fr> Date: Tue, 17 Apr 2007 14:30:15 +0000 Subject: [PATCH] =?UTF-8?q?Correction=20erreurs=20dans=20l'estimation=20de?= =?UTF-8?q?=20la=20disparit=C3=A9.=20Ajout=20de=20tests=20de=20bout=20en?= =?UTF-8?q?=20bout=20pour=20l'estimation=20de=20la=20disparit=C3=A9.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- ...esInterpolateDeformationFieldGenerator.txx | 4 +- .../otbDisparityMapEstimationMethod.h | 60 +-- .../otbDisparityMapEstimationMethod.txx | 93 ++++- Testing/Code/DisparityMap/CMakeLists.txt | 133 +++++- ...enteredRigidDeformationFieldEstimation.cxx | 383 +++++++++++++++++ .../DisparityMap/otbDisparityMapTests.cxx | 3 + .../otbEuler2DDeformationFieldEstimation.cxx | 385 ++++++++++++++++++ ...bTranslationDeformationFieldEstimation.cxx | 371 +++++++++++++++++ 8 files changed, 1387 insertions(+), 45 deletions(-) create mode 100644 Testing/Code/DisparityMap/otbCenteredRigidDeformationFieldEstimation.cxx create mode 100644 Testing/Code/DisparityMap/otbEuler2DDeformationFieldEstimation.cxx create mode 100644 Testing/Code/DisparityMap/otbTranslationDeformationFieldEstimation.cxx diff --git a/Code/DisparityMap/otbBSplinesInterpolateDeformationFieldGenerator.txx b/Code/DisparityMap/otbBSplinesInterpolateDeformationFieldGenerator.txx index ff610948d2..4681f7ec6b 100644 --- a/Code/DisparityMap/otbBSplinesInterpolateDeformationFieldGenerator.txx +++ b/Code/DisparityMap/otbBSplinesInterpolateDeformationFieldGenerator.txx @@ -102,8 +102,8 @@ BSplinesInterpolateDeformationFieldGenerator<TPointSet, TDeformationField> inPixel = inIt.Get(); PixelType outPixel; outPixel.SetSize(2); - outPixel[0]=inPixel[0]; - outPixel[1]=inPixel[1]; + outPixel[0]=-inPixel[0]; + outPixel[1]=-inPixel[1]; outIt.Set(outPixel); } } diff --git a/Code/DisparityMap/otbDisparityMapEstimationMethod.h b/Code/DisparityMap/otbDisparityMapEstimationMethod.h index 45d6e3a107..19d344d4df 100644 --- a/Code/DisparityMap/otbDisparityMapEstimationMethod.h +++ b/Code/DisparityMap/otbDisparityMapEstimationMethod.h @@ -107,18 +107,6 @@ class ITK_EXPORT DisparityMapEstimationMethod /** Smart Pointer type to a DataObject. */ typedef typename itk::DataObject::Pointer DataObjectPointer; - /** Set/Get the Fixed image. */ - itkSetObjectMacro(FixedImage,FixedImageType); - itkGetObjectMacro(FixedImage,FixedImageType); - - /** Set/Get the Moving image. */ - itkSetObjectMacro(MovingImage,MovingImageType); - itkGetObjectMacro(MovingImage,MovingImageType); - - /** Set/Get the points set. */ - itkSetObjectMacro(PointSet,PointSetType); - itkGetObjectMacro(PointSet,PointSetType); - /** Set/Get the Optimizer. */ itkSetObjectMacro(Optimizer,OptimizerType); itkGetObjectMacro(Optimizer,OptimizerType); @@ -147,6 +135,42 @@ class ITK_EXPORT DisparityMapEstimationMethod itkSetMacro(InitialTransformParameters,ParametersType); itkGetConstReferenceMacro(InitialTransformParameters,ParametersType); + /** + * Set the source pointset. + * \param pointset The source pointset. + */ + void SetPointSet(const TPointSet * pointset); + /** + * Get the source pointset. + * \return The source pointset. + */ + const TPointSet * GetPointSet(void); + + /** + * Set the fixed image. + * \param image The fixed image. + **/ + void SetFixedImage(const TFixedImage * image); + + /** + * Get the fixed image. + * \return The fixed image. + **/ + const TFixedImage * GetFixedImage(void); + + /** + * Set the moving image. + * \param image The mobing image. + **/ + void SetMovingImage(const TMovingImage * image); + + /** + * Get the fixed image. + * \return The fixed image. + **/ + const TMovingImage * GetMovingImage(void); + + protected: /** * Constructor. @@ -184,18 +208,6 @@ private: * The interpolator used for local registration. */ InterpolatorPointerType m_Interpolator; - /** - * The input point set. - */ - PointSetPointerType m_PointSet; - /** - * The moving image. - */ - MovingImagePointerType m_MovingImage; - /** - * The fixed image. - */ - FixedImagePointerType m_FixedImage; /** * The initial transform parameters. */ diff --git a/Code/DisparityMap/otbDisparityMapEstimationMethod.txx b/Code/DisparityMap/otbDisparityMapEstimationMethod.txx index 727b5c05c3..3a0ff52a91 100644 --- a/Code/DisparityMap/otbDisparityMapEstimationMethod.txx +++ b/Code/DisparityMap/otbDisparityMapEstimationMethod.txx @@ -32,10 +32,10 @@ template < typename TFixedImage, typename TMovingImage, class TPointSet > DisparityMapEstimationMethod<TFixedImage,TMovingImage,TPointSet> ::DisparityMapEstimationMethod() { - this->SetNumberOfRequiredInputs(0); - this->SetNumberOfRequiredOutputs(1); - m_FixedImage = 0; // has to be provided by the user. - m_MovingImage = 0; // has to be provided by the user. + this->SetNumberOfRequiredInputs(3); + //this->SetNumberOfRequiredOutputs(1); + // this->SetReleaseDataFlag(false); + this->SetReleaseDataBeforeUpdateFlag(false); m_Transform = 0; // has to be provided by the user. m_Interpolator = 0; // has to be provided by the user. m_Metric = 0; // has to be provided by the user. @@ -48,34 +48,101 @@ DisparityMapEstimationMethod<TFixedImage,TMovingImage,TPointSet> /* * Destructor. */ -template < typename TFixedImage, typename TMovingImage, class TPointSet > +template < class TFixedImage, class TMovingImage, class TPointSet > DisparityMapEstimationMethod<TFixedImage,TMovingImage,TPointSet> ::~DisparityMapEstimationMethod() {} +/** + * Set the source pointset. + * \param pointset The source pointset. + */ +template < class TFixedImage, class TMovingImage, class TPointSet > +void +DisparityMapEstimationMethod<TFixedImage,TMovingImage,TPointSet> +::SetPointSet(const TPointSet * pointset) +{ + this->itk::ProcessObject::SetNthInput(0,const_cast<PointSetType *>(pointset)); +} +/** + * Get the source pointset. + * \return The source pointset. + */ +template < class TFixedImage, class TMovingImage, class TPointSet > +const TPointSet * +DisparityMapEstimationMethod<TFixedImage,TMovingImage,TPointSet> +::GetPointSet(void) +{ + return static_cast<const PointSetType *>(this->itk::ProcessObject::GetInput(0)); +} +/** + * Set the fixed image. + * \param image The fixed image. + **/ +template < class TFixedImage, class TMovingImage, class TPointSet > +void +DisparityMapEstimationMethod<TFixedImage,TMovingImage,TPointSet> +::SetFixedImage(const TFixedImage * image) +{ + this->itk::ProcessObject::SetNthInput(1,const_cast<FixedImageType *>(image)); +} +/** + * Get the fixed image. + * \return The fixed image. + **/ +template < class TFixedImage, class TMovingImage, class TPointSet > +const TFixedImage * +DisparityMapEstimationMethod<TFixedImage,TMovingImage,TPointSet> +::GetFixedImage(void) +{ + return static_cast<const FixedImageType *>(this->itk::ProcessObject::GetInput(1)); +} +/** + * Set the moving image. + * \param image The mobing image. + **/ +template < class TFixedImage, class TMovingImage, class TPointSet > +void +DisparityMapEstimationMethod<TFixedImage,TMovingImage,TPointSet> +::SetMovingImage(const TMovingImage * image) +{ + this->itk::ProcessObject::SetNthInput(2,const_cast<MovingImageType *>(image)); +} + /** + * Get the fixed image. + * \return The fixed image. + **/ +template < class TFixedImage, class TMovingImage, class TPointSet > +const TMovingImage * +DisparityMapEstimationMethod<TFixedImage,TMovingImage,TPointSet> +::GetMovingImage(void) +{ + return static_cast<const MovingImageType *>(this->itk::ProcessObject::GetInput(2)); +} + + /** * Main computation method. - * */ -template < typename TFixedImage, typename TMovingImage, class TPointSet > +template < class TFixedImage, class TMovingImage, class TPointSet > void DisparityMapEstimationMethod<TFixedImage,TMovingImage,TPointSet> ::GenerateData(void) { // inputs pointers - FixedImagePointerType fixed = this->GetFixedImage(); - MovingImagePointerType moving = this->GetMovingImage(); - PointSetPointerType pointSet = this->GetPointSet(); + const FixedImageType * fixed = this->GetFixedImage(); + const MovingImageType * moving = this->GetMovingImage(); + const PointSetType * pointSet = this->GetPointSet(); PointSetType* output = this->GetOutput(); // Typedefs typedef typename PointSetType::PointsContainer PointsContainer; - typedef typename PointsContainer::Iterator PointsIterator; + typedef typename PointsContainer::ConstIterator PointsIterator; typedef itk::ImageRegistrationMethod<FixedImageType,MovingImageType> RegistrationType; typedef otb::ExtractROI<FixedPixelType,FixedPixelType> FixedExtractType; typedef otb::ExtractROI<MovingPixelType,MovingPixelType> MovingExtractType; // points retrieving - typename PointsContainer::Pointer points = m_PointSet->GetPoints(); + typename PointsContainer::ConstPointer points = pointSet->GetPoints(); // Iterator set up PointsIterator pointIterator = points->Begin(); @@ -167,7 +234,7 @@ DisparityMapEstimationMethod<TFixedImage,TMovingImage,TPointSet> ++dataId; } } -template < typename TFixedImage, typename TMovingImage, class TPointSet > +template < class TFixedImage, class TMovingImage, class TPointSet > void DisparityMapEstimationMethod<TFixedImage,TMovingImage,TPointSet> ::PrintSelf(std::ostream& os, itk::Indent indent) const diff --git a/Testing/Code/DisparityMap/CMakeLists.txt b/Testing/Code/DisparityMap/CMakeLists.txt index 2c6353591b..6ee8e7e2d3 100644 --- a/Testing/Code/DisparityMap/CMakeLists.txt +++ b/Testing/Code/DisparityMap/CMakeLists.txt @@ -8,9 +8,11 @@ SET(INPUTDATA ${OTB_DATA_ROOT}/Input) SET(IMAGEDATA ${OTB_DATA_ROOT}/LargeInput ) SET(TEMP ${OTBTesting_BINARY_DIR}/Temporary) -#Tolerance sur diff pixel image +# Tolerance for pixel difference SET(TOL 0.0) +#EPSILON loose tolerance for multiplatform support. +SET(EPSILON 0.0000000001) SET(DISPARITYMAP_TESTS ${CXX_TEST_PATH}/otbDisparityMapTests) @@ -43,7 +45,7 @@ ADD_TEST(dmTuNearestPointDeformationFieldGeneratorNew ${DISPARITYMAP_TESTS} otbNearestPointDeformationFieldGeneratorNew) ADD_TEST(dmTvNearestPointDeformationFieldGenerator ${DISPARITYMAP_TESTS} - --compare-image ${TOL} + --compare-image ${EPSILON} ${BASELINE}/dmTvNearestPointDeformationField.hdr ${TEMP}/dmTvNearestPointDeformationField.hdr otbNearestPointDeformationFieldGenerator @@ -56,7 +58,7 @@ ADD_TEST(dmTuNNearestPointsLinearInterpolateDeformationFieldGeneratorNew ${DISPA otbNNearestPointsLinearInterpolateDeformationFieldGeneratorNew) ADD_TEST(dmTvNNearestPointsLinearInterpolateDeformationFieldGenerator ${DISPARITYMAP_TESTS} - --compare-image ${TOL} + --compare-image ${EPSILON} ${BASELINE}/dmTvNNearestPointsLinearInterpolateDeformationField.hdr ${TEMP}/dmTvNNearestPointsLinearInterpolateDeformationField.hdr otbNNearestPointsLinearInterpolateDeformationFieldGenerator @@ -70,7 +72,7 @@ ADD_TEST(dmTuBSplinesInterpolateDeformationFieldGeneratorNew ${DISPARITYMAP_TEST otbNNearestPointsLinearInterpolateDeformationFieldGeneratorNew) ADD_TEST(dmTvBSplinesInterpolateDeformationFieldGenerator ${DISPARITYMAP_TESTS} - --compare-image ${TOL} + --compare-image ${EPSILON} ${BASELINE}/dmTvBSplinesInterpolateDeformationField.hdr ${TEMP}/dmTvBSplinesInterpolateDeformationField.hdr otbBSplinesInterpolateDeformationFieldGenerator @@ -90,7 +92,7 @@ ADD_TEST(dmTuNearestTransformDeformationFieldGeneratorNew ${DISPARITYMAP_TESTS} otbNearestTransformDeformationFieldGeneratorNew) ADD_TEST(dmTvNearestTransformDeformationFieldGenerator ${DISPARITYMAP_TESTS} - --compare-image ${TOL} + --compare-image ${EPSILON} ${BASELINE}/dmTvNearestTransformDeformationField.hdr ${TEMP}/dmTvNearestTransformDeformationField.hdr otbNearestTransformDeformationFieldGenerator @@ -104,14 +106,130 @@ ADD_TEST(dmTuNNearestTransformsLinearInterpolateDeformationFieldGeneratorNew ${D otbNNearestTransformsLinearInterpolateDeformationFieldGeneratorNew) ADD_TEST(dmTvNNearestTransformsLinearInterpolateDeformationFieldGenerator ${DISPARITYMAP_TESTS} - --compare-image ${TOL} + --compare-image ${EPSILON} ${BASELINE}/dmTvNNearestTransformsLinearInterpolateDeformationField.hdr ${TEMP}/dmTvNNearestTransformsLinearInterpolateDeformationField.hdr otbNNearestTransformsLinearInterpolateDeformationFieldGenerator ${TEMP}/dmTvNNearestTransformsLinearInterpolateDeformationField.hdr ) +# ------- Additional tests for deformation fields estimation ---------- + +ADD_TEST(dmTvTranslationDeformationFieldEstimation ${DISPARITYMAP_TESTS} + --compare-n-images ${EPSILON} 10 + ${TEMP}/dmTvTranslationDeformationFieldEstimationOutput_nnt_oi.tif + ${BASELINE}/dmTvTranslationDeformationFieldEstimationOutput_nnt_oi.tif + ${TEMP}/dmTvTranslationDeformationFieldEstimationOutput_bs_df.hdr + ${BASELINE}/dmTvTranslationDeformationFieldEstimationOutput_bs_df.hdr + ${TEMP}/dmTvTranslationDeformationFieldEstimationOutput_bs_oi.tif + ${BASELINE}/dmTvTranslationDeformationFieldEstimationOutput_bs_oi.tif + ${TEMP}/dmTvTranslationDeformationFieldEstimationOutput_np_df.hdr + ${BASELINE}/dmTvTranslationDeformationFieldEstimationOutput_np_df.hdr + ${TEMP}/dmTvTranslationDeformationFieldEstimationOutput_np_oi.tif + ${BASELINE}/dmTvTranslationDeformationFieldEstimationOutput_np_oi.tif + ${TEMP}/dmTvTranslationDeformationFieldEstimationOutput_nnp_df.hdr + ${BASELINE}/dmTvTranslationDeformationFieldEstimationOutput_nnp_df.hdr + ${TEMP}/dmTvTranslationDeformationFieldEstimationOutput_nnp_oi.tif + ${BASELINE}/dmTvTranslationDeformationFieldEstimationOutput_nnp_oi.tif + ${TEMP}/dmTvTranslationDeformationFieldEstimationOutput_nt_df.hdr + ${BASELINE}/dmTvTranslationDeformationFieldEstimationOutput_nt_df.hdr + ${TEMP}/dmTvTranslationDeformationFieldEstimationOutput_nt_oi.tif + ${BASELINE}/dmTvTranslationDeformationFieldEstimationOutput_nt_oi.tif + ${TEMP}/dmTvTranslationDeformationFieldEstimationOutput_nnt_df.hdr + ${BASELINE}/dmTvTranslationDeformationFieldEstimationOutput_nnt_df.hdr + otbTranslationDeformationFieldEstimation + ${INPUTDATA}/ROI_IKO_PAN_LesHalles_sub.tif + ${INPUTDATA}/ROI_IKO_PAN_LesHalles_sub_warped_translation.tif + ${TEMP}/dmTvTranslationDeformationFieldEstimationOutput + 15 10 10 250 0.95 4 50 +) + +ADD_TEST(dmTvEuler2DDeformationFieldEstimation ${DISPARITYMAP_TESTS} + --compare-n-images ${EPSILON} 10 + ${TEMP}/dmTvEuler2DDeformationFieldEstimationOutput_nnt_oi.tif + ${BASELINE}/dmTvEuler2DDeformationFieldEstimationOutput_nnt_oi.tif + ${TEMP}/dmTvEuler2DDeformationFieldEstimationOutput_bs_df.hdr + ${BASELINE}/dmTvEuler2DDeformationFieldEstimationOutput_bs_df.hdr + ${TEMP}/dmTvEuler2DDeformationFieldEstimationOutput_bs_oi.tif + ${BASELINE}/dmTvEuler2DDeformationFieldEstimationOutput_bs_oi.tif + ${TEMP}/dmTvEuler2DDeformationFieldEstimationOutput_np_df.hdr + ${BASELINE}/dmTvEuler2DDeformationFieldEstimationOutput_np_df.hdr + ${TEMP}/dmTvEuler2DDeformationFieldEstimationOutput_np_oi.tif + ${BASELINE}/dmTvEuler2DDeformationFieldEstimationOutput_np_oi.tif + ${TEMP}/dmTvEuler2DDeformationFieldEstimationOutput_nnp_df.hdr + ${BASELINE}/dmTvEuler2DDeformationFieldEstimationOutput_nnp_df.hdr + ${TEMP}/dmTvEuler2DDeformationFieldEstimationOutput_nnp_oi.tif + ${BASELINE}/dmTvEuler2DDeformationFieldEstimationOutput_nnp_oi.tif + ${TEMP}/dmTvEuler2DDeformationFieldEstimationOutput_nt_df.hdr + ${BASELINE}/dmTvEuler2DDeformationFieldEstimationOutput_nt_df.hdr + ${TEMP}/dmTvEuler2DDeformationFieldEstimationOutput_nt_oi.tif + ${BASELINE}/dmTvEuler2DDeformationFieldEstimationOutput_nt_oi.tif + ${TEMP}/dmTvEuler2DDeformationFieldEstimationOutput_nnt_df.hdr + ${BASELINE}/dmTvEuler2DDeformationFieldEstimationOutput_nnt_df.hdr + otbEuler2DDeformationFieldEstimation + ${INPUTDATA}/ROI_IKO_PAN_LesHalles_sub.tif + ${INPUTDATA}/ROI_IKO_PAN_LesHalles_sub_warped_euler2D.tif + ${TEMP}/dmTvEuler2DDeformationFieldEstimationOutput + 15 10 0.01 250 0.95 4 50 128 128 +) +ADD_TEST(dmTvCenteredRigidDeformationFieldEstimation ${DISPARITYMAP_TESTS} + --compare-n-images ${EPSILON} 10 + ${TEMP}/dmTvCenteredRigidDeformationFieldEstimationOutput_nnt_oi.tif + ${BASELINE}/dmTvCenteredRigidDeformationFieldEstimationOutput_nnt_oi.tif + ${TEMP}/dmTvCenteredRigidDeformationFieldEstimationOutput_bs_df.hdr + ${BASELINE}/dmTvCenteredRigidDeformationFieldEstimationOutput_bs_df.hdr + ${TEMP}/dmTvCenteredRigidDeformationFieldEstimationOutput_bs_oi.tif + ${BASELINE}/dmTvCenteredRigidDeformationFieldEstimationOutput_bs_oi.tif + ${TEMP}/dmTvCenteredRigidDeformationFieldEstimationOutput_np_df.hdr + ${BASELINE}/dmTvCenteredRigidDeformationFieldEstimationOutput_np_df.hdr + ${TEMP}/dmTvCenteredRigidDeformationFieldEstimationOutput_np_oi.tif + ${BASELINE}/dmTvCenteredRigidDeformationFieldEstimationOutput_np_oi.tif + ${TEMP}/dmTvCenteredRigidDeformationFieldEstimationOutput_nnp_df.hdr + ${BASELINE}/dmTvCenteredRigidDeformationFieldEstimationOutput_nnp_df.hdr + ${TEMP}/dmTvCenteredRigidDeformationFieldEstimationOutput_nnp_oi.tif + ${BASELINE}/dmTvCenteredRigidDeformationFieldEstimationOutput_nnp_oi.tif + ${TEMP}/dmTvCenteredRigidDeformationFieldEstimationOutput_nt_df.hdr + ${BASELINE}/dmTvCenteredRigidDeformationFieldEstimationOutput_nt_df.hdr + ${TEMP}/dmTvCenteredRigidDeformationFieldEstimationOutput_nt_oi.tif + ${BASELINE}/dmTvCenteredRigidDeformationFieldEstimationOutput_nt_oi.tif + ${TEMP}/dmTvCenteredRigidDeformationFieldEstimationOutput_nnt_df.hdr + ${BASELINE}/dmTvCenteredRigidDeformationFieldEstimationOutput_nnt_df.hdr + otbCenteredRigidDeformationFieldEstimation + ${INPUTDATA}/ROI_IKO_PAN_LesHalles_sub.tif + ${INPUTDATA}/ROI_IKO_PAN_LesHalles_sub_warped_centered_rigid.tif + ${TEMP}/dmTvCenteredRigidDeformationFieldEstimationOutput + 15 5 0.01 250 0.95 4 50 127 65 +) + +ADD_TEST(dmTvSinusoidDeformationFieldEstimation ${DISPARITYMAP_TESTS} + --compare-n-images ${EPSILON} 10 + ${TEMP}/dmTvSinusoidDeformationFieldEstimationOutput_nnt_oi.tif + ${BASELINE}/dmTvSinusoidDeformationFieldEstimationOutput_nnt_oi.tif + ${TEMP}/dmTvSinusoidDeformationFieldEstimationOutput_bs_df.hdr + ${BASELINE}/dmTvSinusoidDeformationFieldEstimationOutput_bs_df.hdr + ${TEMP}/dmTvSinusoidDeformationFieldEstimationOutput_bs_oi.tif + ${BASELINE}/dmTvSinusoidDeformationFieldEstimationOutput_bs_oi.tif + ${TEMP}/dmTvSinusoidDeformationFieldEstimationOutput_np_df.hdr + ${BASELINE}/dmTvSinusoidDeformationFieldEstimationOutput_np_df.hdr + ${TEMP}/dmTvSinusoidDeformationFieldEstimationOutput_np_oi.tif + ${BASELINE}/dmTvSinusoidDeformationFieldEstimationOutput_np_oi.tif + ${TEMP}/dmTvSinusoidDeformationFieldEstimationOutput_nnp_df.hdr + ${BASELINE}/dmTvSinusoidDeformationFieldEstimationOutput_nnp_df.hdr + ${TEMP}/dmTvSinusoidDeformationFieldEstimationOutput_nnp_oi.tif + ${BASELINE}/dmTvSinusoidDeformationFieldEstimationOutput_nnp_oi.tif + ${TEMP}/dmTvSinusoidDeformationFieldEstimationOutput_nt_df.hdr + ${BASELINE}/dmTvSinusoidDeformationFieldEstimationOutput_nt_df.hdr + ${TEMP}/dmTvSinusoidDeformationFieldEstimationOutput_nt_oi.tif + ${BASELINE}/dmTvSinusoidDeformationFieldEstimationOutput_nt_oi.tif + ${TEMP}/dmTvSinusoidDeformationFieldEstimationOutput_nnt_df.hdr + ${BASELINE}/dmTvSinusoidDeformationFieldEstimationOutput_nnt_df.hdr + otbTranslationDeformationFieldEstimation + ${INPUTDATA}/ROI_IKO_PAN_LesHalles_sub.tif + ${INPUTDATA}/ROI_IKO_PAN_LesHalles_sub_warped_sinus.tif + ${TEMP}/dmTvSinusoidDeformationFieldEstimationOutput + 15 10 10 250 0.95 4 25 +) # ------- Fichiers sources CXX ----------------------------------- SET(BasicDisparityMap_SRCS otbDisparityMapEstimationMethodNew.cxx @@ -128,6 +246,9 @@ otbNearestTransformDeformationFieldGeneratorNew.cxx otbNearestTransformDeformationFieldGenerator.cxx otbNNearestTransformsLinearInterpolateDeformationFieldGeneratorNew.cxx otbNNearestTransformsLinearInterpolateDeformationFieldGenerator.cxx +otbTranslationDeformationFieldEstimation.cxx +otbEuler2DDeformationFieldEstimation.cxx +otbCenteredRigidDeformationFieldEstimation.cxx ) INCLUDE_DIRECTORIES("${OTBTesting_BINARY_DIR}") diff --git a/Testing/Code/DisparityMap/otbCenteredRigidDeformationFieldEstimation.cxx b/Testing/Code/DisparityMap/otbCenteredRigidDeformationFieldEstimation.cxx new file mode 100644 index 0000000000..df8b3d0ca0 --- /dev/null +++ b/Testing/Code/DisparityMap/otbCenteredRigidDeformationFieldEstimation.cxx @@ -0,0 +1,383 @@ +/*========================================================================= + + Program: ORFEO Toolbox + Language: C++ + Date: $Date$ + Version: $Revision$ + + + Copyright (c) Centre National d'Etudes Spatiales. All rights reserved. + See OTBCopyright.txt for details. + + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ +#include "otbDisparityMapEstimationMethod.h" +#include "otbImage.h" +#include "otbVectorImage.h" +#include "otbImageFileReader.h" +#include "otbImageFileWriter.h" +#include "itkCenteredRigid2DTransform.h" +#include "itkNormalizedCorrelationImageToImageMetric.h" +#include "itkWindowedSincInterpolateImageFunction.h" +#include "itkZeroFluxNeumannBoundaryCondition.h" +#include "itkGradientDescentOptimizer.h" +#include "itkRescaleIntensityImageFilter.h" + +#include "itkWarpImageFilter.h" +#include "otbNearestPointDeformationFieldGenerator.h" +#include "otbNNearestPointsLinearInterpolateDeformationFieldGenerator.h" +#include "otbBSplinesInterpolateDeformationFieldGenerator.h" +#include "otbNearestTransformDeformationFieldGenerator.h" +#include "otbNNearestTransformsLinearInterpolateDeformationFieldGenerator.h" + + +int otbCenteredRigidDeformationFieldEstimation(int argc, char* argv[]) +{ +try + { + if(argc!=13) + { + std::cout<<"usage: "<<argv[0]<<" fixedFileName movingFileName outputFileName explorationSize windowSize learningRate numberOfIterations metricThreshold nbPointToInterpolate pontsetStep centerx centery"<<std::endl; + return EXIT_SUCCESS; + } + + // Input parameters + const char* fixedFileName = argv[1]; + const char* movingFileName = argv[2]; + const char* outputFileNamePrefix = argv[3]; + const unsigned int exploSize = atoi(argv[4]); + const unsigned int winSize = atoi(argv[5]); + const double learningRate = atof(argv[6]); + const unsigned int niterations = atoi(argv[7]); + const double metricThreshold = atof(argv[8]); + const unsigned int nbPointsToInterpolate = atoi(argv[9]); + const unsigned int step = atoi(argv[10]); + const double centerx = atof(argv[11]); + const double centery = atof(argv[12]); + + // 0. DEFINTIONS + + // Images definition + const unsigned int Dimension=2; + typedef double PixelType; + typedef otb::Image<PixelType,Dimension> ImageType; + typedef otb::Image<unsigned char,Dimension> OutputImageType; + typedef ImageType::IndexType IndexType; + typedef ImageType::SizeType SizeType; + typedef otb::VectorImage<double,Dimension> DeformationFieldType; + + // Transform definition + typedef itk::CenteredRigid2DTransform<double> TransformType; + typedef TransformType::ParametersType ParametersType; + + // Pointset definition + typedef itk::PointSet<ParametersType,Dimension> PointSetType; + typedef PointSetType::PointType PointType; + typedef PointSetType::PointsContainer::Iterator PointSetIteratorType; + typedef PointSetType::PointsContainer PointsContainerType; + typedef PointSetType::PointDataContainer::Iterator PointDataIteratorType; + typedef PointSetType::PointDataContainer PointDataContainerType; + + // Disparity map estimator definition + typedef otb::DisparityMapEstimationMethod<ImageType,ImageType,PointSetType> DMEstimationType; + + // Metric definition + typedef itk::NormalizedCorrelationImageToImageMetric<ImageType,ImageType> MetricType; + + // Interpolator definition + typedef itk::Function::HammingWindowFunction<3> WindowFunctionType; + typedef itk::ZeroFluxNeumannBoundaryCondition<ImageType> ConditionType; + typedef itk::WindowedSincInterpolateImageFunction<ImageType,3,WindowFunctionType,ConditionType ,double> InterpolatorType; + + // Optimizer definition + typedef itk::GradientDescentOptimizer OptimizerType; + + // IO definition + typedef otb::ImageFileReader<ImageType> ReaderType; + typedef otb::ImageFileWriter<OutputImageType> WriterType; + typedef otb::ImageFileWriter<DeformationFieldType> DeformationFieldWriterType; + + // Additional filters + typedef itk::RescaleIntensityImageFilter<ImageType,OutputImageType> RescalerType; + + // Deformation fields generator + typedef otb::NearestPointDeformationFieldGenerator<PointSetType,DeformationFieldType> NearestPointGeneratorType; + typedef otb::NNearestPointsLinearInterpolateDeformationFieldGenerator<PointSetType,DeformationFieldType> NNearestPointGeneratorType; + typedef otb::BSplinesInterpolateDeformationFieldGenerator<PointSetType,DeformationFieldType> BSplinesGeneratorType; + typedef otb::NearestTransformDeformationFieldGenerator<PointSetType,DeformationFieldType> NearestTransformGeneratorType; + typedef otb::NNearestTransformsLinearInterpolateDeformationFieldGenerator<PointSetType,DeformationFieldType> NNearestTransformGeneratorType; + + // Warper + typedef itk::WarpImageFilter<ImageType,ImageType,DeformationFieldType> ImageWarperType; + + //Input images reading + ReaderType::Pointer fixedReader = ReaderType::New(); + ReaderType::Pointer movingReader = ReaderType::New(); + + fixedReader->SetFileName(fixedFileName); + movingReader->SetFileName(movingFileName); + fixedReader->Update(); + movingReader->Update(); + + // 1. POINTSET CREATION + SizeType fixedSize = fixedReader->GetOutput()->GetLargestPossibleRegion().GetSize(); + unsigned int NumberOfXNodes = (fixedSize[0]-2*exploSize-1)/step; + unsigned int NumberOfYNodes = (fixedSize[1]-2*exploSize-1)/step; + unsigned int NumberOfNodes = NumberOfXNodes*NumberOfYNodes; + + std::cout<<"PointSet size: "<<NumberOfNodes<<std::endl; + + IndexType firstNodeIndex; + firstNodeIndex[0] = exploSize; + firstNodeIndex[1] = exploSize; + + PointSetType::Pointer nodes = PointSetType::New(); + unsigned int nodeCounter = 0; + + std::cout << "Node coordinates : " << std::endl; + for(unsigned int x=0; x<NumberOfXNodes; x++) + for(unsigned int y=0; y<NumberOfYNodes; y++) + { + PointType p; + p[0] = firstNodeIndex[0]+x*step; // x coordinate + p[1] = firstNodeIndex[1]+y*step; // y coordinate + std::cout << "Id: " << nodeCounter << " -> " << p << std::endl; + nodes->SetPoint( nodeCounter++, p ); + } + + // Fix to avoid recomputing the disparity for each deformation field generation method. + nodes->SetBufferedRegion(0); + + // 2. DISPARITY MAP ESTIMATION + DMEstimationType::Pointer dmestimator = DMEstimationType::New(); + TransformType::Pointer transform = TransformType::New(); + OptimizerType::Pointer optimizer = OptimizerType::New(); + InterpolatorType::Pointer interpolator = InterpolatorType::New(); + MetricType::Pointer metric = MetricType::New(); + metric->SetSubtractMean(true); + optimizer->MinimizeOn(); + OptimizerType::ScalesType scales(transform->GetNumberOfParameters()); + scales[0] = 3.14; + scales[1] = 0.0001; + scales[2] = 0.0001; + scales[3] = 10; + scales[4] = 10; + optimizer->SetScales(scales); + + // Set up + dmestimator->SetTransform(transform); + dmestimator->SetOptimizer(optimizer); + dmestimator->SetInterpolator(interpolator); + dmestimator->SetMetric(metric); + + // For gradient descent + optimizer->SetLearningRate( learningRate ); + optimizer->SetNumberOfIterations( niterations ); + DMEstimationType::ParametersType initialParameters(transform->GetNumberOfParameters() ); + initialParameters[0] = 0.0; // Initial angle + initialParameters[1] = centerx; // x initial center + initialParameters[2] = centery; // y initial center + initialParameters[3] = 0.0; // Initial offset in mm along X + initialParameters[4] = 0.0; // Initial offset in mm along Y + //Initial parameter set up + dmestimator->SetInitialTransformParameters(initialParameters); + + // inputs wiring + ImageType::SizeType win,explo; + win.Fill(winSize); + explo.Fill(exploSize); + + dmestimator->SetFixedImage(fixedReader->GetOutput()); + dmestimator->SetMovingImage(movingReader->GetOutput()); + dmestimator->SetPointSet(nodes); + dmestimator->SetWinSize(win); + dmestimator->SetExploSize(explo); + dmestimator->SetInitialTransformParameters(initialParameters); + + // 3. DEFORMATION FIELDS COMPUTATION + WriterType::Pointer writer = WriterType::New(); + DeformationFieldWriterType::Pointer dfwriter = DeformationFieldWriterType::New(); + ImageWarperType::Pointer warper = ImageWarperType::New(); + RescalerType::Pointer rescaler = RescalerType::New(); + itk::OStringStream oss; + + //3.a Nearest point deformation field generator + NearestPointGeneratorType::Pointer generator1 = NearestPointGeneratorType::New(); + generator1->SetPointSet(dmestimator->GetOutput()); + generator1->SetOutputOrigin(fixedReader->GetOutput()->GetOrigin()); + generator1->SetOutputSpacing(fixedReader->GetOutput()->GetSpacing()); + generator1->SetOutputSize(fixedReader->GetOutput()->GetLargestPossibleRegion().GetSize()); + generator1->SetMetricThreshold(metricThreshold); + + warper->SetInput(fixedReader->GetOutput()); + warper->SetDeformationField(generator1->GetOutput()); + rescaler->SetInput(warper->GetOutput()); + rescaler->SetOutputMaximum(255); + rescaler->SetOutputMinimum(0); + rescaler->Update(); + + oss.str(""); + oss<<outputFileNamePrefix<<"_np_df.hdr"; + dfwriter->SetFileName(oss.str().c_str()); + dfwriter->SetInput(generator1->GetOutput()); + dfwriter->Update(); + + oss.str(""); + oss<<outputFileNamePrefix<<"_np_oi.tif"; + writer->SetFileName(oss.str().c_str()); + writer->SetInput(rescaler->GetOutput()); + writer->Update(); + + //3.b NNearest points deformation field linear interpolate generator + writer = WriterType::New(); + dfwriter = DeformationFieldWriterType::New(); + warper = ImageWarperType::New(); + rescaler = RescalerType::New(); + + NNearestPointGeneratorType::Pointer generator2 = NNearestPointGeneratorType::New(); + generator2->SetPointSet(dmestimator->GetOutput()); + generator2->SetOutputOrigin(fixedReader->GetOutput()->GetOrigin()); + generator2->SetOutputSpacing(fixedReader->GetOutput()->GetSpacing()); + generator2->SetOutputSize(fixedReader->GetOutput()->GetLargestPossibleRegion().GetSize()); + generator2->SetMetricThreshold(metricThreshold); + generator2->SetNumberOfPoints(nbPointsToInterpolate); + + + warper->SetInput(fixedReader->GetOutput()); + warper->SetDeformationField(generator2->GetOutput()); + rescaler->SetInput(warper->GetOutput()); + rescaler->SetOutputMaximum(255); + rescaler->SetOutputMinimum(0); + rescaler->Update(); + + oss.str(""); + oss<<outputFileNamePrefix<<"_nnp_df.hdr"; + dfwriter->SetFileName(oss.str().c_str()); + dfwriter->SetInput(generator2->GetOutput()); + dfwriter->Update(); + + oss.str(""); + oss<<outputFileNamePrefix<<"_nnp_oi.tif"; + writer->SetFileName(oss.str().c_str()); + writer->SetInput(rescaler->GetOutput()); + writer->Update(); + + //3.c Splines points deformation field linear interpolate generator + writer = WriterType::New(); + dfwriter = DeformationFieldWriterType::New(); + warper = ImageWarperType::New(); + rescaler = RescalerType::New(); + + BSplinesGeneratorType::Pointer generator3 = BSplinesGeneratorType::New(); + generator3->SetPointSet(dmestimator->GetOutput()); + generator3->SetOutputOrigin(fixedReader->GetOutput()->GetOrigin()); + generator3->SetOutputSpacing(fixedReader->GetOutput()->GetSpacing()); + generator3->SetOutputSize(fixedReader->GetOutput()->GetLargestPossibleRegion().GetSize()); + generator3->SetMetricThreshold(metricThreshold); + + + warper->SetInput(fixedReader->GetOutput()); + warper->SetDeformationField(generator3->GetOutput()); + rescaler->SetInput(warper->GetOutput()); + rescaler->SetOutputMaximum(255); + rescaler->SetOutputMinimum(0); + rescaler->Update(); + + oss.str(""); + oss<<outputFileNamePrefix<<"_bs_df.hdr"; + dfwriter->SetFileName(oss.str().c_str()); + dfwriter->SetInput(generator3->GetOutput()); + dfwriter->Update(); + + oss.str(""); + oss<<outputFileNamePrefix<<"_bs_oi.tif"; + writer->SetFileName(oss.str().c_str()); + writer->SetInput(rescaler->GetOutput()); + writer->Update(); + + //3.d Nearest transform deformation field generator + writer = WriterType::New(); + dfwriter = DeformationFieldWriterType::New(); + warper = ImageWarperType::New(); + rescaler = RescalerType::New(); + + NearestTransformGeneratorType::Pointer generator4 = NearestTransformGeneratorType::New(); + generator4->SetPointSet(dmestimator->GetOutput()); + generator4->SetOutputOrigin(fixedReader->GetOutput()->GetOrigin()); + generator4->SetOutputSpacing(fixedReader->GetOutput()->GetSpacing()); + generator4->SetOutputSize(fixedReader->GetOutput()->GetLargestPossibleRegion().GetSize()); + generator4->SetMetricThreshold(metricThreshold); + generator4->SetTransform(transform); + + warper->SetInput(fixedReader->GetOutput()); + warper->SetDeformationField(generator4->GetOutput()); + rescaler->SetInput(warper->GetOutput()); + rescaler->SetOutputMaximum(255); + rescaler->SetOutputMinimum(0); + rescaler->Update(); + + oss.str(""); + oss<<outputFileNamePrefix<<"_nt_df.hdr"; + dfwriter->SetFileName(oss.str().c_str()); + dfwriter->SetInput(generator4->GetOutput()); + dfwriter->Update(); + + oss.str(""); + oss<<outputFileNamePrefix<<"_nt_oi.tif"; + writer->SetFileName(oss.str().c_str()); + writer->SetInput(rescaler->GetOutput()); + writer->Update(); + + //3.e NNearest transforms deformation field linear interpolation generator + writer = WriterType::New(); + dfwriter = DeformationFieldWriterType::New(); + warper = ImageWarperType::New(); + rescaler = RescalerType::New(); + + NNearestTransformGeneratorType::Pointer generator5 = NNearestTransformGeneratorType::New(); + generator5->SetPointSet(dmestimator->GetOutput()); + generator5->SetOutputOrigin(fixedReader->GetOutput()->GetOrigin()); + generator5->SetOutputSpacing(fixedReader->GetOutput()->GetSpacing()); + generator5->SetOutputSize(fixedReader->GetOutput()->GetLargestPossibleRegion().GetSize()); + generator5->SetMetricThreshold(metricThreshold); + generator5->SetTransform(transform); + generator5->SetNumberOfPoints(nbPointsToInterpolate); + + warper->SetInput(fixedReader->GetOutput()); + warper->SetDeformationField(generator5->GetOutput()); + rescaler->SetInput(warper->GetOutput()); + rescaler->SetOutputMaximum(255); + rescaler->SetOutputMinimum(0); + rescaler->Update(); + + oss.str(""); + oss<<outputFileNamePrefix<<"_nnt_df.hdr"; + dfwriter->SetFileName(oss.str().c_str()); + dfwriter->SetInput(generator5->GetOutput()); + dfwriter->Update(); + + oss.str(""); + oss<<outputFileNamePrefix<<"_nnt_oi.tif"; + writer->SetFileName(oss.str().c_str()); + writer->SetInput(rescaler->GetOutput()); + writer->Update(); + +} +catch( itk::ExceptionObject & err ) + + { + std::cout << "Exception itk::ExceptionObject thrown !" << std::endl; + std::cout << err << std::endl; + return EXIT_FAILURE; + } +catch( ... ) + { + std::cout << "Unknown exception thrown !" << std::endl; + return EXIT_FAILURE; + } + return EXIT_SUCCESS; +} diff --git a/Testing/Code/DisparityMap/otbDisparityMapTests.cxx b/Testing/Code/DisparityMap/otbDisparityMapTests.cxx index aff01963a2..196aa8f004 100644 --- a/Testing/Code/DisparityMap/otbDisparityMapTests.cxx +++ b/Testing/Code/DisparityMap/otbDisparityMapTests.cxx @@ -40,4 +40,7 @@ REGISTER_TEST(otbNearestTransformDeformationFieldGeneratorNew); REGISTER_TEST(otbNearestTransformDeformationFieldGenerator); REGISTER_TEST(otbNNearestTransformsLinearInterpolateDeformationFieldGeneratorNew); REGISTER_TEST(otbNNearestTransformsLinearInterpolateDeformationFieldGenerator); +REGISTER_TEST(otbTranslationDeformationFieldEstimation); +REGISTER_TEST(otbEuler2DDeformationFieldEstimation); +REGISTER_TEST(otbCenteredRigidDeformationFieldEstimation); } diff --git a/Testing/Code/DisparityMap/otbEuler2DDeformationFieldEstimation.cxx b/Testing/Code/DisparityMap/otbEuler2DDeformationFieldEstimation.cxx new file mode 100644 index 0000000000..7213b45302 --- /dev/null +++ b/Testing/Code/DisparityMap/otbEuler2DDeformationFieldEstimation.cxx @@ -0,0 +1,385 @@ +/*========================================================================= + + Program: ORFEO Toolbox + Language: C++ + Date: $Date$ + Version: $Revision$ + + + Copyright (c) Centre National d'Etudes Spatiales. All rights reserved. + See OTBCopyright.txt for details. + + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ +#include "otbDisparityMapEstimationMethod.h" +#include "otbImage.h" +#include "otbVectorImage.h" +#include "otbImageFileReader.h" +#include "otbImageFileWriter.h" +#include "itkEuler2DTransform.h" +#include "itkNormalizedCorrelationImageToImageMetric.h" +#include "itkWindowedSincInterpolateImageFunction.h" +#include "itkZeroFluxNeumannBoundaryCondition.h" +#include "itkGradientDescentOptimizer.h" +#include "itkRescaleIntensityImageFilter.h" + +#include "itkWarpImageFilter.h" +#include "otbNearestPointDeformationFieldGenerator.h" +#include "otbNNearestPointsLinearInterpolateDeformationFieldGenerator.h" +#include "otbBSplinesInterpolateDeformationFieldGenerator.h" +#include "otbNearestTransformDeformationFieldGenerator.h" +#include "otbNNearestTransformsLinearInterpolateDeformationFieldGenerator.h" + + +int otbEuler2DDeformationFieldEstimation(int argc, char* argv[]) +{ +try + { + if(argc!=13) + { + std::cout<<"usage: "<<argv[0]<<" fixedFileName movingFileName outputFileName explorationSize windowSize learningRate numberOfIterations metricThreshold nbPointToInterpolate pontsetStep centerx centery"<<std::endl; + return EXIT_SUCCESS; + } + + // Input parameters + const char* fixedFileName = argv[1]; + const char* movingFileName = argv[2]; + const char* outputFileNamePrefix = argv[3]; + const unsigned int exploSize = atoi(argv[4]); + const unsigned int winSize = atoi(argv[5]); + const double learningRate = atof(argv[6]); + const unsigned int niterations = atoi(argv[7]); + const double metricThreshold = atof(argv[8]); + const unsigned int nbPointsToInterpolate = atoi(argv[9]); + const unsigned int step = atoi(argv[10]); + const double centerx = atof(argv[11]); + const double centery = atof(argv[12]); + + // 0. DEFINTIONS + + // Images definition + const unsigned int Dimension=2; + typedef double PixelType; + typedef otb::Image<PixelType,Dimension> ImageType; + typedef otb::Image<unsigned char,Dimension> OutputImageType; + typedef ImageType::IndexType IndexType; + typedef ImageType::SizeType SizeType; + typedef otb::VectorImage<double,Dimension> DeformationFieldType; + + // Transform definition + typedef itk::Euler2DTransform<double> TransformType; + typedef TransformType::ParametersType ParametersType; + + // Pointset definition + typedef itk::PointSet<ParametersType,Dimension> PointSetType; + typedef PointSetType::PointType PointType; + typedef PointSetType::PointsContainer::Iterator PointSetIteratorType; + typedef PointSetType::PointsContainer PointsContainerType; + typedef PointSetType::PointDataContainer::Iterator PointDataIteratorType; + typedef PointSetType::PointDataContainer PointDataContainerType; + + // Disparity map estimator definition + typedef otb::DisparityMapEstimationMethod<ImageType,ImageType,PointSetType> DMEstimationType; + + // Metric definition + typedef itk::NormalizedCorrelationImageToImageMetric<ImageType,ImageType> MetricType; + + // Interpolator definition + typedef itk::Function::HammingWindowFunction<3> WindowFunctionType; + typedef itk::ZeroFluxNeumannBoundaryCondition<ImageType> ConditionType; + typedef itk::WindowedSincInterpolateImageFunction<ImageType,3,WindowFunctionType,ConditionType ,double> InterpolatorType; + + // Optimizer definition + typedef itk::GradientDescentOptimizer OptimizerType; + + // IO definition + typedef otb::ImageFileReader<ImageType> ReaderType; + typedef otb::ImageFileWriter<OutputImageType> WriterType; + typedef otb::ImageFileWriter<DeformationFieldType> DeformationFieldWriterType; + + // Additional filters + typedef itk::RescaleIntensityImageFilter<ImageType,OutputImageType> RescalerType; + + // Deformation fields generator + typedef otb::NearestPointDeformationFieldGenerator<PointSetType,DeformationFieldType> NearestPointGeneratorType; + typedef otb::NNearestPointsLinearInterpolateDeformationFieldGenerator<PointSetType,DeformationFieldType> NNearestPointGeneratorType; + typedef otb::BSplinesInterpolateDeformationFieldGenerator<PointSetType,DeformationFieldType> BSplinesGeneratorType; + typedef otb::NearestTransformDeformationFieldGenerator<PointSetType,DeformationFieldType> NearestTransformGeneratorType; + typedef otb::NNearestTransformsLinearInterpolateDeformationFieldGenerator<PointSetType,DeformationFieldType> NNearestTransformGeneratorType; + + // Warper + typedef itk::WarpImageFilter<ImageType,ImageType,DeformationFieldType> ImageWarperType; + + //Input images reading + ReaderType::Pointer fixedReader = ReaderType::New(); + ReaderType::Pointer movingReader = ReaderType::New(); + + fixedReader->SetFileName(fixedFileName); + movingReader->SetFileName(movingFileName); + fixedReader->Update(); + movingReader->Update(); + + // 1. POINTSET CREATION + SizeType fixedSize = fixedReader->GetOutput()->GetLargestPossibleRegion().GetSize(); + unsigned int NumberOfXNodes = (fixedSize[0]-2*exploSize-1)/step; + unsigned int NumberOfYNodes = (fixedSize[1]-2*exploSize-1)/step; + unsigned int NumberOfNodes = NumberOfXNodes*NumberOfYNodes; + + std::cout<<"PointSet size: "<<NumberOfNodes<<std::endl; + + IndexType firstNodeIndex; + firstNodeIndex[0] = exploSize; + firstNodeIndex[1] = exploSize; + + PointSetType::Pointer nodes = PointSetType::New(); + unsigned int nodeCounter = 0; + + std::cout << "Node coordinates : " << std::endl; + for(unsigned int x=0; x<NumberOfXNodes; x++) + for(unsigned int y=0; y<NumberOfYNodes; y++) + { + PointType p; + p[0] = firstNodeIndex[0]+x*step; // x coordinate + p[1] = firstNodeIndex[1]+y*step; // y coordinate + std::cout << "Id: " << nodeCounter << " -> " << p << std::endl; + nodes->SetPoint( nodeCounter++, p ); + } + + // Fix to avoid recomputing the disparity for each deformation field generation method. + nodes->SetBufferedRegion(0); + + // 2. DISPARITY MAP ESTIMATION + DMEstimationType::Pointer dmestimator = DMEstimationType::New(); + TransformType::Pointer transform = TransformType::New(); + OptimizerType::Pointer optimizer = OptimizerType::New(); + InterpolatorType::Pointer interpolator = InterpolatorType::New(); + MetricType::Pointer metric = MetricType::New(); + metric->SetSubtractMean(true); + optimizer->MinimizeOn(); + OptimizerType::ScalesType scales(3); + scales[0] = 3.14; + scales[1] = 10; + scales[2] = 10; + optimizer->SetScales(scales); + + // Transform center + TransformType::CenterType center; + center[0] = centerx; + center[1] = centery; + transform->SetCenter(center); + + // Set up + dmestimator->SetTransform(transform); + dmestimator->SetOptimizer(optimizer); + dmestimator->SetInterpolator(interpolator); + dmestimator->SetMetric(metric); + + // For gradient descent + optimizer->SetLearningRate( learningRate ); + optimizer->SetNumberOfIterations( niterations ); + DMEstimationType::ParametersType initialParameters(transform->GetNumberOfParameters() ); + initialParameters[0] = 0.0; // Initial angle + initialParameters[1] = 0.0; // Initial offset in mm along X + initialParameters[2] = 0.0; // Initial offset in mm along Y + //Initial parameter set up + dmestimator->SetInitialTransformParameters(initialParameters); + + // inputs wiring + ImageType::SizeType win,explo; + win.Fill(winSize); + explo.Fill(exploSize); + + dmestimator->SetFixedImage(fixedReader->GetOutput()); + dmestimator->SetMovingImage(movingReader->GetOutput()); + dmestimator->SetPointSet(nodes); + dmestimator->SetWinSize(win); + dmestimator->SetExploSize(explo); + dmestimator->SetInitialTransformParameters(initialParameters); + + // 3. DEFORMATION FIELDS COMPUTATION + WriterType::Pointer writer = WriterType::New(); + DeformationFieldWriterType::Pointer dfwriter = DeformationFieldWriterType::New(); + ImageWarperType::Pointer warper = ImageWarperType::New(); + RescalerType::Pointer rescaler = RescalerType::New(); + itk::OStringStream oss; + + //3.a Nearest point deformation field generator + NearestPointGeneratorType::Pointer generator1 = NearestPointGeneratorType::New(); + generator1->SetPointSet(dmestimator->GetOutput()); + generator1->SetOutputOrigin(fixedReader->GetOutput()->GetOrigin()); + generator1->SetOutputSpacing(fixedReader->GetOutput()->GetSpacing()); + generator1->SetOutputSize(fixedReader->GetOutput()->GetLargestPossibleRegion().GetSize()); + generator1->SetMetricThreshold(metricThreshold); + + warper->SetInput(fixedReader->GetOutput()); + warper->SetDeformationField(generator1->GetOutput()); + rescaler->SetInput(warper->GetOutput()); + rescaler->SetOutputMaximum(255); + rescaler->SetOutputMinimum(0); + rescaler->Update(); + + oss.str(""); + oss<<outputFileNamePrefix<<"_np_df.hdr"; + dfwriter->SetFileName(oss.str().c_str()); + dfwriter->SetInput(generator1->GetOutput()); + dfwriter->Update(); + + oss.str(""); + oss<<outputFileNamePrefix<<"_np_oi.tif"; + writer->SetFileName(oss.str().c_str()); + writer->SetInput(rescaler->GetOutput()); + writer->Update(); + + //3.b NNearest points deformation field linear interpolate generator + writer = WriterType::New(); + dfwriter = DeformationFieldWriterType::New(); + warper = ImageWarperType::New(); + rescaler = RescalerType::New(); + + NNearestPointGeneratorType::Pointer generator2 = NNearestPointGeneratorType::New(); + generator2->SetPointSet(dmestimator->GetOutput()); + generator2->SetOutputOrigin(fixedReader->GetOutput()->GetOrigin()); + generator2->SetOutputSpacing(fixedReader->GetOutput()->GetSpacing()); + generator2->SetOutputSize(fixedReader->GetOutput()->GetLargestPossibleRegion().GetSize()); + generator2->SetMetricThreshold(metricThreshold); + generator2->SetNumberOfPoints(nbPointsToInterpolate); + + + warper->SetInput(fixedReader->GetOutput()); + warper->SetDeformationField(generator2->GetOutput()); + rescaler->SetInput(warper->GetOutput()); + rescaler->SetOutputMaximum(255); + rescaler->SetOutputMinimum(0); + rescaler->Update(); + + oss.str(""); + oss<<outputFileNamePrefix<<"_nnp_df.hdr"; + dfwriter->SetFileName(oss.str().c_str()); + dfwriter->SetInput(generator2->GetOutput()); + dfwriter->Update(); + + oss.str(""); + oss<<outputFileNamePrefix<<"_nnp_oi.tif"; + writer->SetFileName(oss.str().c_str()); + writer->SetInput(rescaler->GetOutput()); + writer->Update(); + + //3.c Splines points deformation field linear interpolate generator + writer = WriterType::New(); + dfwriter = DeformationFieldWriterType::New(); + warper = ImageWarperType::New(); + rescaler = RescalerType::New(); + + BSplinesGeneratorType::Pointer generator3 = BSplinesGeneratorType::New(); + generator3->SetPointSet(dmestimator->GetOutput()); + generator3->SetOutputOrigin(fixedReader->GetOutput()->GetOrigin()); + generator3->SetOutputSpacing(fixedReader->GetOutput()->GetSpacing()); + generator3->SetOutputSize(fixedReader->GetOutput()->GetLargestPossibleRegion().GetSize()); + generator3->SetMetricThreshold(metricThreshold); + + + warper->SetInput(fixedReader->GetOutput()); + warper->SetDeformationField(generator3->GetOutput()); + rescaler->SetInput(warper->GetOutput()); + rescaler->SetOutputMaximum(255); + rescaler->SetOutputMinimum(0); + rescaler->Update(); + + oss.str(""); + oss<<outputFileNamePrefix<<"_bs_df.hdr"; + dfwriter->SetFileName(oss.str().c_str()); + dfwriter->SetInput(generator3->GetOutput()); + dfwriter->Update(); + + oss.str(""); + oss<<outputFileNamePrefix<<"_bs_oi.tif"; + writer->SetFileName(oss.str().c_str()); + writer->SetInput(rescaler->GetOutput()); + writer->Update(); + + //3.d Nearest transform deformation field generator + writer = WriterType::New(); + dfwriter = DeformationFieldWriterType::New(); + warper = ImageWarperType::New(); + rescaler = RescalerType::New(); + + NearestTransformGeneratorType::Pointer generator4 = NearestTransformGeneratorType::New(); + generator4->SetPointSet(dmestimator->GetOutput()); + generator4->SetOutputOrigin(fixedReader->GetOutput()->GetOrigin()); + generator4->SetOutputSpacing(fixedReader->GetOutput()->GetSpacing()); + generator4->SetOutputSize(fixedReader->GetOutput()->GetLargestPossibleRegion().GetSize()); + generator4->SetMetricThreshold(metricThreshold); + generator4->SetTransform(transform); + + warper->SetInput(fixedReader->GetOutput()); + warper->SetDeformationField(generator4->GetOutput()); + rescaler->SetInput(warper->GetOutput()); + rescaler->SetOutputMaximum(255); + rescaler->SetOutputMinimum(0); + rescaler->Update(); + + oss.str(""); + oss<<outputFileNamePrefix<<"_nt_df.hdr"; + dfwriter->SetFileName(oss.str().c_str()); + dfwriter->SetInput(generator4->GetOutput()); + dfwriter->Update(); + + oss.str(""); + oss<<outputFileNamePrefix<<"_nt_oi.tif"; + writer->SetFileName(oss.str().c_str()); + writer->SetInput(rescaler->GetOutput()); + writer->Update(); + + //3.e NNearest transforms deformation field linear interpolation generator + writer = WriterType::New(); + dfwriter = DeformationFieldWriterType::New(); + warper = ImageWarperType::New(); + rescaler = RescalerType::New(); + + NNearestTransformGeneratorType::Pointer generator5 = NNearestTransformGeneratorType::New(); + generator5->SetPointSet(dmestimator->GetOutput()); + generator5->SetOutputOrigin(fixedReader->GetOutput()->GetOrigin()); + generator5->SetOutputSpacing(fixedReader->GetOutput()->GetSpacing()); + generator5->SetOutputSize(fixedReader->GetOutput()->GetLargestPossibleRegion().GetSize()); + generator5->SetMetricThreshold(metricThreshold); + generator5->SetTransform(transform); + generator5->SetNumberOfPoints(nbPointsToInterpolate); + + warper->SetInput(fixedReader->GetOutput()); + warper->SetDeformationField(generator5->GetOutput()); + rescaler->SetInput(warper->GetOutput()); + rescaler->SetOutputMaximum(255); + rescaler->SetOutputMinimum(0); + rescaler->Update(); + + oss.str(""); + oss<<outputFileNamePrefix<<"_nnt_df.hdr"; + dfwriter->SetFileName(oss.str().c_str()); + dfwriter->SetInput(generator5->GetOutput()); + dfwriter->Update(); + + oss.str(""); + oss<<outputFileNamePrefix<<"_nnt_oi.tif"; + writer->SetFileName(oss.str().c_str()); + writer->SetInput(rescaler->GetOutput()); + writer->Update(); + +} +catch( itk::ExceptionObject & err ) + + { + std::cout << "Exception itk::ExceptionObject thrown !" << std::endl; + std::cout << err << std::endl; + return EXIT_FAILURE; + } +catch( ... ) + { + std::cout << "Unknown exception thrown !" << std::endl; + return EXIT_FAILURE; + } + return EXIT_SUCCESS; +} diff --git a/Testing/Code/DisparityMap/otbTranslationDeformationFieldEstimation.cxx b/Testing/Code/DisparityMap/otbTranslationDeformationFieldEstimation.cxx new file mode 100644 index 0000000000..060e4adf23 --- /dev/null +++ b/Testing/Code/DisparityMap/otbTranslationDeformationFieldEstimation.cxx @@ -0,0 +1,371 @@ +/*========================================================================= + + Program: ORFEO Toolbox + Language: C++ + Date: $Date$ + Version: $Revision$ + + + Copyright (c) Centre National d'Etudes Spatiales. All rights reserved. + See OTBCopyright.txt for details. + + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ +#include "otbDisparityMapEstimationMethod.h" +#include "otbImage.h" +#include "otbVectorImage.h" +#include "otbImageFileReader.h" +#include "otbImageFileWriter.h" +#include "itkTranslationTransform.h" +#include "itkNormalizedCorrelationImageToImageMetric.h" +#include "itkWindowedSincInterpolateImageFunction.h" +#include "itkZeroFluxNeumannBoundaryCondition.h" +#include "itkGradientDescentOptimizer.h" +#include "itkRescaleIntensityImageFilter.h" + +#include "itkWarpImageFilter.h" +#include "otbNearestPointDeformationFieldGenerator.h" +#include "otbNNearestPointsLinearInterpolateDeformationFieldGenerator.h" +#include "otbBSplinesInterpolateDeformationFieldGenerator.h" +#include "otbNearestTransformDeformationFieldGenerator.h" +#include "otbNNearestTransformsLinearInterpolateDeformationFieldGenerator.h" + + +int otbTranslationDeformationFieldEstimation(int argc, char* argv[]) +{ +try + { + if(argc!=11) + { + std::cout<<"usage: "<<argv[0]<<" fixedFileName movingFileName outputFileName explorationSize windowSize learningRate numberOfIterations metricThreshold nbPointToInterpolate pontsetStep"<<std::endl; + return EXIT_SUCCESS; + } + + // Input parameters + const char* fixedFileName = argv[1]; + const char* movingFileName = argv[2]; + const char* outputFileNamePrefix = argv[3]; + const unsigned int exploSize = atoi(argv[4]); + const unsigned int winSize = atoi(argv[5]); + const double learningRate = atof(argv[6]); + const unsigned int niterations = atoi(argv[7]); + const double metricThreshold = atof(argv[8]); + const unsigned int nbPointsToInterpolate = atoi(argv[9]); + const unsigned int step = atoi(argv[10]); + + // 0. DEFINTIONS + + // Images definition + const unsigned int Dimension=2; + typedef double PixelType; + typedef otb::Image<PixelType,Dimension> ImageType; + typedef otb::Image<unsigned char,Dimension> OutputImageType; + typedef ImageType::IndexType IndexType; + typedef ImageType::SizeType SizeType; + typedef otb::VectorImage<double,Dimension> DeformationFieldType; + + // Transform definition + typedef itk::TranslationTransform<double,Dimension> TransformType; + typedef TransformType::ParametersType ParametersType; + + // Pointset definition + typedef itk::PointSet<ParametersType,Dimension> PointSetType; + typedef PointSetType::PointType PointType; + typedef PointSetType::PointsContainer::Iterator PointSetIteratorType; + typedef PointSetType::PointsContainer PointsContainerType; + typedef PointSetType::PointDataContainer::Iterator PointDataIteratorType; + typedef PointSetType::PointDataContainer PointDataContainerType; + + // Disparity map estimator definition + typedef otb::DisparityMapEstimationMethod<ImageType,ImageType,PointSetType> DMEstimationType; + + // Metric definition + typedef itk::NormalizedCorrelationImageToImageMetric<ImageType,ImageType> MetricType; + + // Interpolator definition + typedef itk::Function::HammingWindowFunction<3> WindowFunctionType; + typedef itk::ZeroFluxNeumannBoundaryCondition<ImageType> ConditionType; + typedef itk::WindowedSincInterpolateImageFunction<ImageType,3,WindowFunctionType,ConditionType ,double> InterpolatorType; + + // Optimizer definition + typedef itk::GradientDescentOptimizer OptimizerType; + + // IO definition + typedef otb::ImageFileReader<ImageType> ReaderType; + typedef otb::ImageFileWriter<OutputImageType> WriterType; + typedef otb::ImageFileWriter<DeformationFieldType> DeformationFieldWriterType; + + // Additional filters + typedef itk::RescaleIntensityImageFilter<ImageType,OutputImageType> RescalerType; + + // Deformation fields generator + typedef otb::NearestPointDeformationFieldGenerator<PointSetType,DeformationFieldType> NearestPointGeneratorType; + typedef otb::NNearestPointsLinearInterpolateDeformationFieldGenerator<PointSetType,DeformationFieldType> NNearestPointGeneratorType; + typedef otb::BSplinesInterpolateDeformationFieldGenerator<PointSetType,DeformationFieldType> BSplinesGeneratorType; + typedef otb::NearestTransformDeformationFieldGenerator<PointSetType,DeformationFieldType> NearestTransformGeneratorType; + typedef otb::NNearestTransformsLinearInterpolateDeformationFieldGenerator<PointSetType,DeformationFieldType> NNearestTransformGeneratorType; + + // Warper + typedef itk::WarpImageFilter<ImageType,ImageType,DeformationFieldType> ImageWarperType; + + //Input images reading + ReaderType::Pointer fixedReader = ReaderType::New(); + ReaderType::Pointer movingReader = ReaderType::New(); + + fixedReader->SetFileName(fixedFileName); + movingReader->SetFileName(movingFileName); + fixedReader->Update(); + movingReader->Update(); + + // 1. POINTSET CREATION + SizeType fixedSize = fixedReader->GetOutput()->GetLargestPossibleRegion().GetSize(); + unsigned int NumberOfXNodes = (fixedSize[0]-2*exploSize-1)/step; + unsigned int NumberOfYNodes = (fixedSize[1]-2*exploSize-1)/step; + unsigned int NumberOfNodes = NumberOfXNodes*NumberOfYNodes; + + std::cout<<"PointSet size: "<<NumberOfNodes<<std::endl; + + IndexType firstNodeIndex; + firstNodeIndex[0] = exploSize; + firstNodeIndex[1] = exploSize; + + PointSetType::Pointer nodes = PointSetType::New(); + unsigned int nodeCounter = 0; + + std::cout << "Node coordinates : " << std::endl; + for(unsigned int x=0; x<NumberOfXNodes; x++) + for(unsigned int y=0; y<NumberOfYNodes; y++) + { + PointType p; + p[0] = firstNodeIndex[0]+x*step; // x coordinate + p[1] = firstNodeIndex[1]+y*step; // y coordinate + std::cout << "Id: " << nodeCounter << " -> " << p << std::endl; + nodes->SetPoint( nodeCounter++, p ); + } + + // Fix to avoid recomputing the disparity for each deformation field generation method. + nodes->SetBufferedRegion(0); + + // 2. DISPARITY MAP ESTIMATION + DMEstimationType::Pointer dmestimator = DMEstimationType::New(); + TransformType::Pointer transform = TransformType::New(); + OptimizerType::Pointer optimizer = OptimizerType::New(); + InterpolatorType::Pointer interpolator = InterpolatorType::New(); + MetricType::Pointer metric = MetricType::New(); + metric->SetSubtractMean(true); + optimizer->MinimizeOn(); + + // Set up + dmestimator->SetTransform(transform); + dmestimator->SetOptimizer(optimizer); + dmestimator->SetInterpolator(interpolator); + dmestimator->SetMetric(metric); + + // For gradient descent + optimizer->SetLearningRate( learningRate ); + optimizer->SetNumberOfIterations( niterations ); + DMEstimationType::ParametersType initialParameters(transform->GetNumberOfParameters() ); + initialParameters[0] = 0.0; // Initial offset in mm along X + initialParameters[1] = 0.0; // Initial offset in mm along Y + //Initial parameter set up + dmestimator->SetInitialTransformParameters(initialParameters); + + // inputs wiring + ImageType::SizeType win,explo; + win.Fill(winSize); + explo.Fill(exploSize); + + dmestimator->SetFixedImage(fixedReader->GetOutput()); + dmestimator->SetMovingImage(movingReader->GetOutput()); + dmestimator->SetPointSet(nodes); + dmestimator->SetWinSize(win); + dmestimator->SetExploSize(explo); + dmestimator->SetInitialTransformParameters(initialParameters); + + // 3. DEFORMATION FIELDS COMPUTATION + WriterType::Pointer writer = WriterType::New(); + DeformationFieldWriterType::Pointer dfwriter = DeformationFieldWriterType::New(); + ImageWarperType::Pointer warper = ImageWarperType::New(); + RescalerType::Pointer rescaler = RescalerType::New(); + itk::OStringStream oss; + + //3.a Nearest point deformation field generator + NearestPointGeneratorType::Pointer generator1 = NearestPointGeneratorType::New(); + generator1->SetPointSet(dmestimator->GetOutput()); + generator1->SetOutputOrigin(fixedReader->GetOutput()->GetOrigin()); + generator1->SetOutputSpacing(fixedReader->GetOutput()->GetSpacing()); + generator1->SetOutputSize(fixedReader->GetOutput()->GetLargestPossibleRegion().GetSize()); + generator1->SetMetricThreshold(metricThreshold); + + warper->SetInput(fixedReader->GetOutput()); + warper->SetDeformationField(generator1->GetOutput()); + rescaler->SetInput(warper->GetOutput()); + rescaler->SetOutputMaximum(255); + rescaler->SetOutputMinimum(0); + rescaler->Update(); + + oss.str(""); + oss<<outputFileNamePrefix<<"_np_df.hdr"; + dfwriter->SetFileName(oss.str().c_str()); + dfwriter->SetInput(generator1->GetOutput()); + dfwriter->Update(); + + oss.str(""); + oss<<outputFileNamePrefix<<"_np_oi.tif"; + writer->SetFileName(oss.str().c_str()); + writer->SetInput(rescaler->GetOutput()); + writer->Update(); + + //3.b NNearest points deformation field linear interpolate generator + writer = WriterType::New(); + dfwriter = DeformationFieldWriterType::New(); + warper = ImageWarperType::New(); + rescaler = RescalerType::New(); + + NNearestPointGeneratorType::Pointer generator2 = NNearestPointGeneratorType::New(); + generator2->SetPointSet(dmestimator->GetOutput()); + generator2->SetOutputOrigin(fixedReader->GetOutput()->GetOrigin()); + generator2->SetOutputSpacing(fixedReader->GetOutput()->GetSpacing()); + generator2->SetOutputSize(fixedReader->GetOutput()->GetLargestPossibleRegion().GetSize()); + generator2->SetMetricThreshold(metricThreshold); + generator2->SetNumberOfPoints(nbPointsToInterpolate); + + + warper->SetInput(fixedReader->GetOutput()); + warper->SetDeformationField(generator2->GetOutput()); + rescaler->SetInput(warper->GetOutput()); + rescaler->SetOutputMaximum(255); + rescaler->SetOutputMinimum(0); + rescaler->Update(); + + oss.str(""); + oss<<outputFileNamePrefix<<"_nnp_df.hdr"; + dfwriter->SetFileName(oss.str().c_str()); + dfwriter->SetInput(generator2->GetOutput()); + dfwriter->Update(); + + oss.str(""); + oss<<outputFileNamePrefix<<"_nnp_oi.tif"; + writer->SetFileName(oss.str().c_str()); + writer->SetInput(rescaler->GetOutput()); + writer->Update(); + + //3.c Splines points deformation field linear interpolate generator + writer = WriterType::New(); + dfwriter = DeformationFieldWriterType::New(); + warper = ImageWarperType::New(); + rescaler = RescalerType::New(); + + BSplinesGeneratorType::Pointer generator3 = BSplinesGeneratorType::New(); + generator3->SetPointSet(dmestimator->GetOutput()); + generator3->SetOutputOrigin(fixedReader->GetOutput()->GetOrigin()); + generator3->SetOutputSpacing(fixedReader->GetOutput()->GetSpacing()); + generator3->SetOutputSize(fixedReader->GetOutput()->GetLargestPossibleRegion().GetSize()); + generator3->SetMetricThreshold(metricThreshold); + + + warper->SetInput(fixedReader->GetOutput()); + warper->SetDeformationField(generator3->GetOutput()); + rescaler->SetInput(warper->GetOutput()); + rescaler->SetOutputMaximum(255); + rescaler->SetOutputMinimum(0); + rescaler->Update(); + + oss.str(""); + oss<<outputFileNamePrefix<<"_bs_df.hdr"; + dfwriter->SetFileName(oss.str().c_str()); + dfwriter->SetInput(generator3->GetOutput()); + dfwriter->Update(); + + oss.str(""); + oss<<outputFileNamePrefix<<"_bs_oi.tif"; + writer->SetFileName(oss.str().c_str()); + writer->SetInput(rescaler->GetOutput()); + writer->Update(); + + //3.d Nearest transform deformation field generator + writer = WriterType::New(); + dfwriter = DeformationFieldWriterType::New(); + warper = ImageWarperType::New(); + rescaler = RescalerType::New(); + + NearestTransformGeneratorType::Pointer generator4 = NearestTransformGeneratorType::New(); + generator4->SetPointSet(dmestimator->GetOutput()); + generator4->SetOutputOrigin(fixedReader->GetOutput()->GetOrigin()); + generator4->SetOutputSpacing(fixedReader->GetOutput()->GetSpacing()); + generator4->SetOutputSize(fixedReader->GetOutput()->GetLargestPossibleRegion().GetSize()); + generator4->SetMetricThreshold(metricThreshold); + generator4->SetTransform(transform); + + warper->SetInput(fixedReader->GetOutput()); + warper->SetDeformationField(generator4->GetOutput()); + rescaler->SetInput(warper->GetOutput()); + rescaler->SetOutputMaximum(255); + rescaler->SetOutputMinimum(0); + rescaler->Update(); + + oss.str(""); + oss<<outputFileNamePrefix<<"_nt_df.hdr"; + dfwriter->SetFileName(oss.str().c_str()); + dfwriter->SetInput(generator4->GetOutput()); + dfwriter->Update(); + + oss.str(""); + oss<<outputFileNamePrefix<<"_nt_oi.tif"; + writer->SetFileName(oss.str().c_str()); + writer->SetInput(rescaler->GetOutput()); + writer->Update(); + + //3.e NNearest transforms deformation field linear interpolation generator + writer = WriterType::New(); + dfwriter = DeformationFieldWriterType::New(); + warper = ImageWarperType::New(); + rescaler = RescalerType::New(); + + NNearestTransformGeneratorType::Pointer generator5 = NNearestTransformGeneratorType::New(); + generator5->SetPointSet(dmestimator->GetOutput()); + generator5->SetOutputOrigin(fixedReader->GetOutput()->GetOrigin()); + generator5->SetOutputSpacing(fixedReader->GetOutput()->GetSpacing()); + generator5->SetOutputSize(fixedReader->GetOutput()->GetLargestPossibleRegion().GetSize()); + generator5->SetMetricThreshold(metricThreshold); + generator5->SetTransform(transform); + generator5->SetNumberOfPoints(nbPointsToInterpolate); + + warper->SetInput(fixedReader->GetOutput()); + warper->SetDeformationField(generator5->GetOutput()); + rescaler->SetInput(warper->GetOutput()); + rescaler->SetOutputMaximum(255); + rescaler->SetOutputMinimum(0); + rescaler->Update(); + + oss.str(""); + oss<<outputFileNamePrefix<<"_nnt_df.hdr"; + dfwriter->SetFileName(oss.str().c_str()); + dfwriter->SetInput(generator5->GetOutput()); + dfwriter->Update(); + + oss.str(""); + oss<<outputFileNamePrefix<<"_nnt_oi.tif"; + writer->SetFileName(oss.str().c_str()); + writer->SetInput(rescaler->GetOutput()); + writer->Update(); + +} +catch( itk::ExceptionObject & err ) + + { + std::cout << "Exception itk::ExceptionObject thrown !" << std::endl; + std::cout << err << std::endl; + return EXIT_FAILURE; + } +catch( ... ) + { + std::cout << "Unknown exception thrown !" << std::endl; + return EXIT_FAILURE; + } + return EXIT_SUCCESS; +} -- GitLab