diff --git a/Code/DisparityMap/otbFineRegistrationImageFilter.h b/Code/DisparityMap/otbFineRegistrationImageFilter.h new file mode 100644 index 0000000000000000000000000000000000000000..d117d94cb8cbc86e983345f64d09875e00ba69ea --- /dev/null +++ b/Code/DisparityMap/otbFineRegistrationImageFilter.h @@ -0,0 +1,197 @@ +/*========================================================================= + + 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. + +=========================================================================*/ +#ifndef __otbFineRegistrationImageFilter_h +#define __otbFineRegistrationImageFilter_h + +#include "itkImageToImageFilter.h" +#include "itkInterpolateImageFunction.h" +#include "itkContinuousIndex.h" + +#include "itkTranslationTransform.h" +#include "itkImageToImageMetric.txx" + +namespace otb +{ + +/** \class FineRegistrationImageFilter + * \brief Computes a deformation field between two images + * + * This class compute a correlation field and the associated + * deformation field between two images. + * + * The m_Radius parameter defines the size of the window to compute + * the normalized correlation. + * + * The m_SearchRadius parameter defines the size of the area where to + * search for a correlation peak in the moving image. + * + * Once the correlation peak has been found, it can be further refined + * by trying to fit a quadric surface in the premises of the + * maxima, thus obtaining sub-pixel precision. This option is + * activated by setting m_RefineLocation to true. + * + * TOutputCorrelation is supposed to be a scalar image with floating + * point precision. + * + * TOutputDeformationField is supposed to be an image whose pixel is a + * fixed size vector of size 2, with floating point precision. + * + * Images spacing is ignored during computation and in the output + * deformation field. + * + * All computation are done in double precision. + * + * \ingroup IntensityImageFilters Multithreaded, Streamed + * + */ +template <class TInputImage, class T0utputCorrelation, class TOutputDeformationField> +class ITK_EXPORT FineRegistrationImageFilter : public itk::ImageToImageFilter<TInputImage,T0utputCorrelation> +{ +public: + /** Standard class typedefs. */ + typedef FineRegistrationImageFilter Self; + typedef itk::ImageToImageFilter<TInputImage,T0utputCorrelation> Superclass; + typedef itk::SmartPointer<Self> Pointer; + typedef itk::SmartPointer<const Self> ConstPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro(FineRegistrationImageFilter, ImageToImageFilter); + + /** Some convenient typedefs. */ + typedef typename T0utputCorrelation::RegionType OutputImageRegionType; + typedef typename TOutputDeformationField::PixelType DeformationValueType; + typedef typename TInputImage::Pointer InputImagePointerType; + typedef typename TInputImage::RegionType InputImageRegionType; + typedef typename TInputImage::SizeType SizeType; + typedef typename TInputImage::IndexType IndexType; + typedef typename TInputImage::SpacingType SpacingType; + typedef typename itk::InterpolateImageFunction + <TInputImage, double> InterpolatorType; + typedef typename InterpolatorType::Pointer InterpolatorPointerType; + typedef typename itk::ContinuousIndex<double, 2> ContinuousIndexType; + typedef typename itk::ImageToImageMetric<TInputImage, + TInputImage> MetricType; + typedef typename MetricType::Pointer MetricPointerType; + typedef typename itk::TranslationTransform<double,2> TranslationType; + typedef typename TranslationType::Pointer TranslationPointerType; + + /** Set/Get the Metric used to compare images */ + itkSetObjectMacro(Metric,MetricType); + itkGetObjectMacro(Metric,MetricType); + + /** Set/Get the interpolator used to inpterolate moving image at non-grid positions */ + itkSetObjectMacro(Interpolator,InterpolatorType); + itkGetObjectMacro(Interpolator,InterpolatorType); + + /** Connect one of the operands for pixel-wise addition */ + void SetFixedInput( const TInputImage * image); + + /** Connect one of the operands for pixel-wise addition */ + void SetMovingInput( const TInputImage * image); + + /** Get the inputs */ + const TInputImage * GetFixedInput(); + const TInputImage * GetMovingInput(); + + /** Get the output deformation field */ + TOutputDeformationField * GetOutputDeformationField(); + + /** Set the radius of the area on which metric is evaluated */ + itkSetMacro(Radius, SizeType); + itkGetMacro(Radius, SizeType); + + /** Set the searh radius */ + itkSetMacro(SearchRadius,SizeType); + itkGetMacro(SearchRadius,SizeType); + + /** Set/Get subpixel accuracy */ + itkSetMacro(SubPixelAccuracy,double); + itkGetMacro(SubPixelAccuracy,double); + + /** True if metric should be minimized. False otherwise */ + itkSetMacro(Minimize,bool); + itkBooleanMacro(Minimize); + + /** True if deformation field takes spacing into account. False otherwise */ + itkSetMacro(UseSpacing,bool); + itkBooleanMacro(UseSpacing); + + /** Set unsigned int radius */ + void SetRadius(unsigned int radius) + { + m_Radius.Fill(radius); + } + +/** Set unsigned int radius */ + void SetSearchRadius(unsigned int radius) + { + m_SearchRadius.Fill(radius); + } + +protected: + /** Constructor */ + FineRegistrationImageFilter(); + /** Destructor */ + virtual ~FineRegistrationImageFilter() {}; + + /** Threaded generate data */ + virtual void GenerateData(); + + /** Generate the input requested regions */ + virtual void GenerateInputRequestedRegion(void); + +private: + FineRegistrationImageFilter(const Self&); //purposely not implemented + void operator=(const Self&); //purposely not implemented + + /** The radius for correlation */ + SizeType m_Radius; + + /** The search radius */ + SizeType m_SearchRadius; + + /** Minimize/maximize metric */ + bool m_Minimize; + + /** If true, deformation field uses spacing. Otherwise, uses pixel grid */ + bool m_UseSpacing; + + /** Search step */ + double m_SubPixelAccuracy; + + /** The interpolator */ + InterpolatorPointerType m_Interpolator; + + /** The metric */ + MetricPointerType m_Metric; + + /** The translation */ + TranslationPointerType m_Translation; + +}; + +} // end namespace otb + +#ifndef OTB_MANUAL_INSTANTIATION +#include "otbFineRegistrationImageFilter.txx" +#endif + +#endif diff --git a/Code/DisparityMap/otbFineRegistrationImageFilter.txx b/Code/DisparityMap/otbFineRegistrationImageFilter.txx new file mode 100644 index 0000000000000000000000000000000000000000..7d0d501a3d1a0a3dd2c0ef283dd95786ce51637d --- /dev/null +++ b/Code/DisparityMap/otbFineRegistrationImageFilter.txx @@ -0,0 +1,359 @@ +/*========================================================================= + + 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. + +=========================================================================*/ +#ifndef __otbFineRegistrationImageFilter_txx +#define __otbFineRegistrationImageFilter_txx + +#include "otbFineRegistrationImageFilter.h" + +#include "itkProgressReporter.h" +#include "itkLinearInterpolateImageFunction.h" +#include "itkImageRegionIteratorWithIndex.h" +#include "itkNormalizedCorrelationImageToImageMetric.h" + +namespace otb +{ +/** + * Constructor + */ +template <class TInputImage, class T0utputCorrelation, class TOutputDeformationField> +FineRegistrationImageFilter<TInputImage,T0utputCorrelation,TOutputDeformationField> +::FineRegistrationImageFilter() + { + this->SetNumberOfRequiredInputs( 2 ); + this->SetNumberOfOutputs(2); + this->SetNthOutput(1,TOutputDeformationField::New()); + + // Default radius + m_Radius.Fill(2); + m_SearchRadius.Fill(4); + + // Default sub-pixel precision + m_SubPixelAccuracy = 0.1; + + // Flags + m_UseSpacing = true; + m_Minimize = true; + + // Default currentMetric + m_Metric = itk::NormalizedCorrelationImageToImageMetric + <TInputImage,TInputImage>::New(); + + // Default interpolator + m_Interpolator = itk::LinearInterpolateImageFunction<TInputImage,double>::New(); + + // Translation + m_Translation = TranslationType::New(); + } + +template <class TInputImage, class T0utputCorrelation, class TOutputDeformationField> +void +FineRegistrationImageFilter<TInputImage,T0utputCorrelation,TOutputDeformationField> +::SetFixedInput( const TInputImage * image ) + { + // Process object is not const-correct so the const casting is required. + SetNthInput(0, const_cast<TInputImage *>( image )); + } + +template <class TInputImage, class T0utputCorrelation, class TOutputDeformationField> +void +FineRegistrationImageFilter<TInputImage,T0utputCorrelation,TOutputDeformationField> +::SetMovingInput( const TInputImage * image) + { + // Process object is not const-correct so the const casting is required. + SetNthInput(1, const_cast<TInputImage *>( image )); + } + +template <class TInputImage, class T0utputCorrelation, class TOutputDeformationField> +const TInputImage * +FineRegistrationImageFilter<TInputImage,T0utputCorrelation,TOutputDeformationField> +::GetFixedInput() + { + if (this->GetNumberOfInputs()<1) + { + return 0; + } + return static_cast<const TInputImage *>(this->itk::ProcessObject::GetInput(0)); + } + +template <class TInputImage, class T0utputCorrelation, class TOutputDeformationField> +const TInputImage * +FineRegistrationImageFilter<TInputImage,T0utputCorrelation,TOutputDeformationField> +::GetMovingInput() + { + if (this->GetNumberOfInputs()<2) + { + return 0; + } + return static_cast<const TInputImage *>(this->itk::ProcessObject::GetInput(1)); + } + +template <class TInputImage, class T0utputCorrelation, class TOutputDeformationField> +TOutputDeformationField * +FineRegistrationImageFilter<TInputImage,T0utputCorrelation,TOutputDeformationField> +::GetOutputDeformationField() + { + if (this->GetNumberOfOutputs()<2) + { + return 0; + } + return static_cast<TOutputDeformationField *>(this->itk::ProcessObject::GetOutput(1)); + } + + +template <class TInputImage, class TOutputCorrelation, class TOutputDeformationField> +void +FineRegistrationImageFilter<TInputImage,TOutputCorrelation,TOutputDeformationField> +::GenerateInputRequestedRegion() + { + // call the superclass' implementation of this method + Superclass::GenerateInputRequestedRegion(); + + // get pointers to the input and output + TInputImage * fixedPtr = const_cast< TInputImage * >( this->GetFixedInput()); + TInputImage * movingPtr = const_cast< TInputImage * >( this->GetMovingInput()); + + TOutputCorrelation * outputPtr = this->GetOutput(); + + if ( !fixedPtr || !movingPtr || !outputPtr ) + { + return; + } + + // get a copy of the fixed requested region (should equal the output + // requested region) + typename TInputImage::RegionType fixedRequestedRegion; + fixedRequestedRegion = fixedPtr->GetRequestedRegion(); + // pad the input requested region by the operator radius + fixedRequestedRegion.PadByRadius( m_Radius ); + + // get a copy of the moving requested region (should equal the output + // requested region) + typename TInputImage::RegionType movingRequestedRegion; + movingRequestedRegion = movingPtr->GetRequestedRegion(); + // pad the input requested region by the operator radius + movingRequestedRegion.PadByRadius( m_SearchRadius + m_Radius ); + + // crop the fixed region at the fixed's largest possible region + if ( fixedRequestedRegion.Crop(fixedPtr->GetLargestPossibleRegion())) + { + fixedPtr->SetRequestedRegion( fixedRequestedRegion ); + } + else + { + // Couldn't crop the region (requested region is outside the largest + // possible region). Throw an exception. + // store what we tried to request (prior to trying to crop) + fixedPtr->SetRequestedRegion( fixedRequestedRegion ); + + // build an exception + itk::InvalidRequestedRegionError e(__FILE__, __LINE__); + itk::OStringStream msg; + msg << this->GetNameOfClass() + << "::GenerateInputRequestedRegion()"; + e.SetLocation(msg.str().c_str()); + e.SetDescription("Requested region is (at least partially) outside the largest possible region of image 1."); + e.SetDataObject(fixedPtr); + throw e; + } + + // crop the moving region at the moving's largest possible region + if ( movingRequestedRegion.Crop(movingPtr->GetLargestPossibleRegion())) + { + movingPtr->SetRequestedRegion( movingRequestedRegion ); + } + else + { + // Couldn't crop the region (requested region is outside the largest + // possible region). Throw an exception. + // store what we tried to request (prior to trying to crop) + movingPtr->SetRequestedRegion( movingRequestedRegion ); + + // build an exception + itk::InvalidRequestedRegionError e(__FILE__, __LINE__); + itk::OStringStream msg; + msg << this->GetNameOfClass() + << "::GenerateInputRequestedRegion()"; + e.SetLocation(msg.str().c_str()); + e.SetDescription("Requested region is (at least partially) outside the largest possible region of image 1."); + e.SetDataObject(movingPtr); + throw e; + } + + return; + } + +template <class TInputImage, class TOutputCorrelation, class TOutputDeformationField> +void +FineRegistrationImageFilter<TInputImage,TOutputCorrelation,TOutputDeformationField> +::GenerateData() + { + // Allocate outputs + this->AllocateOutputs(); + + // Get the image pointers + const TInputImage * fixedPtr = this->GetFixedInput(); + const TInputImage * movingPtr = this->GetMovingInput(); + TOutputCorrelation * outputPtr = this->GetOutput(); + TOutputDeformationField * outputDfPtr = this->GetOutputDeformationField(); + + // Wire currentMetric + m_Interpolator->SetInputImage(this->GetMovingInput()); + m_Metric->SetTransform(m_Translation); + m_Metric->SetInterpolator(m_Interpolator); + m_Metric->SetFixedImage(this->GetFixedInput()); + m_Metric->SetMovingImage(this->GetMovingInput()); + m_Metric->SetComputeGradient(false); + + /** Output iterators */ + itk::ImageRegionIteratorWithIndex<TOutputCorrelation> outputIt(outputPtr,outputPtr->GetRequestedRegion()); + itk::ImageRegionIterator<TOutputDeformationField> outputDfIt(outputDfPtr,outputPtr->GetRequestedRegion()); + outputIt.GoToBegin(); + outputDfIt.GoToBegin(); + + // support progress methods/callbacks + itk::ProgressReporter progress(this, 0, outputPtr->GetRequestedRegion().GetNumberOfPixels()); + + // GoToBegin + outputIt.GoToBegin(); + outputDfIt.GoToBegin(); + + // Correl, max correl, maxPosition + double currentMetric, optMetric; + + // Optimal translation parameters + typename TranslationType::ParametersType params(2), optParams(2), tmpOptParams(2); + + // Final deformation value + DeformationValueType deformationValue; + + // Get fixed image spacing + SpacingType fixedSpacing = fixedPtr->GetSpacing(); + + // Walk the images + while (!outputIt.IsAtEnd() && !outputDfIt.IsAtEnd() ) + { + // Reset + if(m_Minimize) + { + optMetric = itk::NumericTraits<double>::max(); + } + else + { + optMetric = itk::NumericTraits<double>::NonpositiveMin(); + } + optParams.Fill(0); + + // Build the region on which to compute the currentMetric + InputImageRegionType currentMetricRegion; + currentMetricRegion.SetIndex(outputIt.GetIndex()); + SizeType size; + size.Fill(1); + currentMetricRegion.SetSize(size); + currentMetricRegion.PadByRadius(m_Radius); + currentMetricRegion.Crop(fixedPtr->GetLargestPossibleRegion()); + m_Metric->SetFixedImageRegion(currentMetricRegion); + m_Metric->Initialize(); + + // Compute the correlation at each location + for(int i =-(int)m_SearchRadius[0]; i<=(int)m_SearchRadius[0];++i) + { + for(int j=-(int)m_SearchRadius[1]; j<=(int)m_SearchRadius[1];++j) + { + params[0]=static_cast<double>(i*fixedSpacing[0]); + params[1]=static_cast<double>(j*fixedSpacing[1]); + + // compute currentMetric + currentMetric = m_Metric->GetValue(params); + + // Check for maximum + if((m_Minimize && (currentMetric < optMetric)) || (!m_Minimize && (currentMetric > optMetric))) + { + optMetric = currentMetric; + optParams = params; + } + } + } + +// // Exhaustive sub-pixel +// tmpOptParams = optParams; +// for(double dx = -vcl_abs(fixedSpacing[0]);dx<vcl_abs(fixedSpacing[0]);dx+=m_SubPixelAccuracy) +// { +// for(double dy = -vcl_abs(fixedSpacing[1]);dy<vcl_abs(fixedSpacing[1]);dy+=m_SubPixelAccuracy) +// { +// params[0] = tmpOptParams[0]+dx; +// params[1] = tmpOptParams[1]+dy; +// +// // compute currentMetric +// currentMetric = m_Metric->GetValue(params); +// +// // Check for maximum +// if((m_Minimize && (currentMetric < optMetric)) || (!m_Minimize && (currentMetric > optMetric))) +// { +// optMetric = currentMetric; +// optParams = params; +// } +// } +// } + + // Dichotomic sub-pixel + SpacingType subPixelSpacing = fixedSpacing; + while(subPixelSpacing[0]> m_SubPixelAccuracy || subPixelSpacing[1]> m_SubPixelAccuracy) + { + // Perform 1 step of dichotomic search + subPixelSpacing/=2.; + + // Store last opt params + tmpOptParams = optParams; + + for(int i =-1; i<=1;i+=2) + { + for(int j=-1; j<=1;j+=2) + { + params = tmpOptParams; + params[0]+=static_cast<double>(i*subPixelSpacing[0]); + params[1]+=static_cast<double>(j*subPixelSpacing[1]); + + // compute currentMetric + currentMetric = m_Metric->GetValue(params); + + // Check for maximum + if((m_Minimize && (currentMetric < optMetric)) || (!m_Minimize && (currentMetric > optMetric))) + { + optMetric = currentMetric; + optParams = params; + } + } + } + } + + // Store the offset and the correlation value + outputIt.Set(optMetric); + deformationValue[0] = optParams[0]; + deformationValue[1] = optParams[1]; + outputDfIt.Set(deformationValue); + // Update iterators + ++outputIt; + ++outputDfIt; + + // Update progress + progress.CompletedPixel(); + } + } +} // end namespace otb + +#endif diff --git a/Testing/Code/DisparityMap/CMakeLists.txt b/Testing/Code/DisparityMap/CMakeLists.txt index 93cf6e05dea5853f331341ed5aa7d368fef27f1d..7f71ba6af6e826f69e96e6dd9b03b7f1939fc7bf 100644 --- a/Testing/Code/DisparityMap/CMakeLists.txt +++ b/Testing/Code/DisparityMap/CMakeLists.txt @@ -5,6 +5,7 @@ SET(BASELINE ${OTB_DATA_ROOT}/Baseline/OTB/Images) SET(BASELINE_FILES ${OTB_DATA_ROOT}/Baseline/OTB/Files) SET(INPUTDATA ${OTB_DATA_ROOT}/Input) SET(TEMP ${OTBTesting_BINARY_DIR}/Temporary) +SET(EXAMPLEDATA ${OTB_DATA_ROOT}/Examples) # Tolerance for pixel difference SET(NOTOL 0.0) @@ -311,6 +312,79 @@ ADD_TEST(dmTvStreamingWarpImageFilter ${DISPARITYMAP_TESTS3} 5 ) +# ------- otb::FineRegistrationImageFilter ---------- + +ADD_TEST(feTuFineRegistrationImageFilterNew ${DISPARITYMAP_TESTS3} + otbFineRegistrationImageFilterNew +) + +ADD_TEST(feTvFineRegistrationImageFilterTestWithCorrelation ${DISPARITYMAP_TESTS3} + --compare-n-images ${EPSILON_8} 2 + ${BASELINE}/feTvFineRegistrationImageFilterTestWithCorrelationMetric.tif + ${TEMP}/feTvFineRegistrationImageFilterTestWithCorrelationMetric.tif + ${BASELINE}/feTvFineRegistrationImageFilterTestWithCorrelationField.tif + ${TEMP}/feTvFineRegistrationImageFilterTestWithCorrelationField.tif + otbFineRegistrationImageFilterTest + ${EXAMPLEDATA}/StereoFixed.png # fixedFileName + ${EXAMPLEDATA}/StereoMoving.png # movingFileName + ${TEMP}/feTvFineRegistrationImageFilterTestWithCorrelationMetric.tif # output correlFileName + ${TEMP}/feTvFineRegistrationImageFilterTestWithCorrelationField.tif # output fieldFileName + 3 # radius + 2 # sradius + 0.1 # precision + 0 # Correlation +) + +ADD_TEST(feTvFineRegistrationImageFilterTestWithNormalizedCorrelation ${DISPARITYMAP_TESTS3} + --compare-n-images ${EPSILON_8} 2 + ${BASELINE}/feTvFineRegistrationImageFilterTestWithNormalizedCorrelationMetric.tif + ${TEMP}/feTvFineRegistrationImageFilterTestWithNormalizedCorrelationMetric.tif + ${BASELINE}/feTvFineRegistrationImageFilterTestWithNormalizedCorrelationField.tif + ${TEMP}/feTvFineRegistrationImageFilterTestWithNormalizedCorrelationField.tif + otbFineRegistrationImageFilterTest + ${EXAMPLEDATA}/StereoFixed.png # fixedFileName + ${EXAMPLEDATA}/StereoMoving.png # movingFileName + ${TEMP}/feTvFineRegistrationImageFilterTestWithNormalizedCorrelationMetric.tif # output correlFileName + ${TEMP}/feTvFineRegistrationImageFilterTestWithNormalizedCorrelationField.tif # output fieldFileName + 3 # radius + 2 # sradius + 0.1 # precision + 1 # Normalized Correlation +) + +ADD_TEST(feTvFineRegistrationImageFilterTestWithMeanSquare ${DISPARITYMAP_TESTS3} + --compare-n-images ${EPSILON_8} 2 + ${BASELINE}/feTvFineRegistrationImageFilterTestWithMeanSquareMetric.tif + ${TEMP}/feTvFineRegistrationImageFilterTestWithMeanSquareMetric.tif + ${BASELINE}/feTvFineRegistrationImageFilterTestWithMeanSquareField.tif + ${TEMP}/feTvFineRegistrationImageFilterTestWithMeanSquareField.tif + otbFineRegistrationImageFilterTest + ${EXAMPLEDATA}/StereoFixed.png # fixedFileName + ${EXAMPLEDATA}/StereoMoving.png # movingFileName + ${TEMP}/feTvFineRegistrationImageFilterTestWithMeanSquareMetric.tif # output correlFileName + ${TEMP}/feTvFineRegistrationImageFilterTestWithMeanSquareField.tif # output fieldFileName + 3 # radius + 2 # sradius + 0.1 # precision + 2 # Mean square +) + +ADD_TEST(feTvFineRegistrationImageFilterTestWithMeanReciprocalDifference ${DISPARITYMAP_TESTS3} + --compare-n-images ${EPSILON_8} 2 + ${BASELINE}/feTvFineRegistrationImageFilterTestWithMeanReciprocalDifferenceMetric.tif + ${TEMP}/feTvFineRegistrationImageFilterTestWithMeanReciprocalDifferenceMetric.tif + ${BASELINE}/feTvFineRegistrationImageFilterTestWithMeanReciprocalDifferenceField.tif + ${TEMP}/feTvFineRegistrationImageFilterTestWithMeanReciprocalDifferenceField.tif + otbFineRegistrationImageFilterTest + ${EXAMPLEDATA}/StereoFixed.png # fixedFileName + ${EXAMPLEDATA}/StereoMoving.png # movingFileName + ${TEMP}/feTvFineRegistrationImageFilterTestWithMeanReciprocalDifferenceMetric.tif # output correlFileName + ${TEMP}/feTvFineRegistrationImageFilterTestWithMeanReciprocalDifferenceField.tif # output fieldFileName + 3 # radius + 2 # sradius + 0.1 # precision + 3 # Mean reciprocal difference +) # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ otbDisparityMapTests3 ~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -349,6 +423,8 @@ SET(BasicDisparityMap_SRCS3 otbDisparityMapTests3.cxx otbStreamingWarpImageFilterNew.cxx otbStreamingWarpImageFilter.cxx +otbFineRegistrationImageFilterNew.cxx +otbFineRegistrationImageFilterTest.cxx ) OTB_ADD_EXECUTABLE(otbDisparityMapTests1 "${BasicDisparityMap_SRCS1}" "OTBDisparityMap;OTBIO;OTBTesting") diff --git a/Testing/Code/DisparityMap/otbDisparityMapTests3.cxx b/Testing/Code/DisparityMap/otbDisparityMapTests3.cxx index b80af500feac7c755065cdfbb714ef2fe4f90dd4..c0f9baf9a4103fe2c33b701a3ce48ac20b008486 100644 --- a/Testing/Code/DisparityMap/otbDisparityMapTests3.cxx +++ b/Testing/Code/DisparityMap/otbDisparityMapTests3.cxx @@ -28,4 +28,6 @@ void RegisterTests() { REGISTER_TEST(otbStreamingWarpImageFilterNew); REGISTER_TEST(otbStreamingWarpImageFilter); + REGISTER_TEST(otbFineRegistrationImageFilterNew); + REGISTER_TEST(otbFineRegistrationImageFilterTest); } diff --git a/Testing/Code/DisparityMap/otbFineRegistrationImageFilterNew.cxx b/Testing/Code/DisparityMap/otbFineRegistrationImageFilterNew.cxx new file mode 100644 index 0000000000000000000000000000000000000000..e2f585f94f1f41b2a568d2e7ba5f5e15fa2f1567 --- /dev/null +++ b/Testing/Code/DisparityMap/otbFineRegistrationImageFilterNew.cxx @@ -0,0 +1,39 @@ +/*========================================================================= + + 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. + +=========================================================================*/ +#if defined(_MSC_VER) +#pragma warning ( disable : 4786 ) +#endif + +#include "otbImage.h" +#include "itkFixedArray.h" +#include "otbFineRegistrationImageFilter.h" + + +int otbFineRegistrationImageFilterNew( int argc, char * argv[] ) +{ + typedef double PixelType; + + typedef itk::FixedArray<PixelType> DeformationValueType; + typedef otb::Image< PixelType> ImageType; + typedef otb::Image<DeformationValueType,2> FieldImageType; + typedef otb::FineRegistrationImageFilter<ImageType,ImageType,FieldImageType> RegistrationFilterType; + + RegistrationFilterType::Pointer Registration = RegistrationFilterType::New(); + + return EXIT_SUCCESS; +} diff --git a/Testing/Code/DisparityMap/otbFineRegistrationImageFilterTest.cxx b/Testing/Code/DisparityMap/otbFineRegistrationImageFilterTest.cxx new file mode 100644 index 0000000000000000000000000000000000000000..5ee00b270ac2530f08eaf4d540529db21742d247 --- /dev/null +++ b/Testing/Code/DisparityMap/otbFineRegistrationImageFilterTest.cxx @@ -0,0 +1,141 @@ +/*========================================================================= + + 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. + +=========================================================================*/ +#if defined(_MSC_VER) +#pragma warning ( disable : 4786 ) +#endif + +#include "itkFixedArray.h" +#include "otbImage.h" +#include "otbImageFileReader.h" +#include "otbImageFileWriter.h" +#include "otbFineRegistrationImageFilter.h" +#include "otbStandardFilterWatcher.h" +#include "itkTimeProbe.h" + +#include "itkNormalizedCorrelationImageToImageMetric.h" +#include "itkMeanReciprocalSquareDifferenceImageToImageMetric.h" +#include "itkMeanSquaresImageToImageMetric.h" + +int otbFineRegistrationImageFilterTest( int argc, char * argv[] ) +{ + if(argc!=9) + { + std::cerr<<"Usage: "<<argv[0]<<" fixed_fname moving_fname output_correl output_field radius search_radius subpixPrecision metric(0=CC,1=NCC,2=MeanSquare,3=Mean reciprocal square difference)"<<std::endl; + return EXIT_FAILURE; + } + const char * fixedFileName = argv[1]; + const char * movingFileName = argv[2]; + const char * correlFileName = argv[3]; + const char * fieldFileName = argv[4]; + const unsigned int radius = atoi(argv[5]); + const unsigned int sradius = atoi(argv[6]); + const double precision = atof(argv[7]); + const unsigned int metric = atoi(argv[8]); + + typedef double PixelType; + const unsigned int Dimension = 2; + + typedef itk::FixedArray<PixelType,Dimension> DeformationValueType; + typedef otb::Image< PixelType, Dimension > ImageType; + typedef otb::Image<DeformationValueType,Dimension> FieldImageType; + typedef otb::ImageFileReader< ImageType > ReaderType; + typedef otb::ImageFileWriter< ImageType > CorrelWriterType; + typedef otb::ImageFileWriter< FieldImageType> FieldWriterType; + typedef otb::FineRegistrationImageFilter<ImageType,ImageType,FieldImageType> RegistrationFilterType; + + ReaderType::Pointer freader = ReaderType::New(); + freader->SetFileName(fixedFileName); + //freader->Update(); + + ReaderType::Pointer mreader = ReaderType::New(); + mreader->SetFileName(movingFileName); + //mreader->Update(); + + RegistrationFilterType::Pointer registration = RegistrationFilterType::New(); + registration->SetFixedInput(freader->GetOutput()); + registration->SetMovingInput(mreader->GetOutput()); + registration->SetRadius(radius); + registration->SetSearchRadius(sradius); + registration->SetSubPixelAccuracy(precision); + + switch(metric) + { + case 0: + { + std::cout<<"Metric: correlation"<<std::endl; + typedef itk::NormalizedCorrelationImageToImageMetric<ImageType,ImageType> NCCType; + NCCType::Pointer metricPtr = NCCType::New(); + metricPtr->SubtractMeanOff(); + registration->SetMetric(metricPtr); + registration->MinimizeOn(); + break; + } + case 1: + { + std::cout<<"Metric: normalized correlation"<<std::endl; + typedef itk::NormalizedCorrelationImageToImageMetric<ImageType,ImageType> NCCType; + NCCType::Pointer metricPtr = NCCType::New(); + metricPtr->SubtractMeanOn(); + registration->SetMetric(metricPtr); + registration->MinimizeOn(); + break; + } + case 2: + { + std::cout<<"Metric: mean squares"<<std::endl; + typedef itk::MeanSquaresImageToImageMetric<ImageType,ImageType> MeanSquareType; + MeanSquareType::Pointer metricPtr = MeanSquareType::New(); + registration->SetMetric(metricPtr); + registration->MinimizeOn(); + break; + } + case 3: + { + std::cout<<"Metric: mean reciprocal square difference"<<std::endl; + typedef itk::MeanReciprocalSquareDifferenceImageToImageMetric<ImageType,ImageType> MRSDType; + MRSDType::Pointer metricPtr = MRSDType::New(); + registration->SetMetric(metricPtr); + registration->MinimizeOff(); + break; + } + default: + { + std::cerr<<"Metric id should be between 0 and 3"<<std::endl; + return EXIT_FAILURE; + } + } + + otb::StandardFilterWatcher watcher(registration,"Registration"); + itk::TimeProbe chrono; + chrono.Start(); + registration->Update(); + chrono.Stop(); + std::cout<<chrono.GetMeanTime()<<"\t"; + + CorrelWriterType::Pointer correlWriter = CorrelWriterType::New(); + correlWriter->SetFileName(correlFileName); + correlWriter->SetInput(registration->GetOutput()); + correlWriter->Update(); + + FieldWriterType::Pointer fieldWriter = FieldWriterType::New(); + fieldWriter->SetFileName(fieldFileName); + fieldWriter->SetInput(registration->GetOutputDeformationField()); + fieldWriter->Update(); + + return EXIT_SUCCESS; +}