MRG

Code/DisparityMap/otbFineRegistrationImageFilter.h

+/*=========================================================================
+
+  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 using a given metric.
+ *
+ * This filter tries to find at each location of the fixed image the corresponding best matching
+ * patch of radius set by SetRadius() method in the moving image within a search windows whose radius
+ * is defined by SetSearchRadius() method.
+ *
+ * To do so, it optimizes a metric set using the SetMetric() method, which accepts any itk metric
+ * deriving from the itk::ImageToImageMetric. The MinimizeOn()/MinimizeOff() flag allow to search for
+ * minimum or maximum depending on the metric.
+ *
+ * Once a coarse (pixel wise) offset has been found, this match is further refined using dichotomic search
+ * until sub-pixel accuracy given by the SetSubPixelAccuracy() is reached.
+ *
+ * The filter proposes two outputs: GetOutput() return the image of the metric optimum at each location, and
+ * the GetOutputDeformationField() method returns the corresponding offset.
+ *
+ * If the UseImageSapcingOn() flag is used, the output deformation field takes the input image spacing into account.
+ * otherwise, the deformation field is expressed in pixels.
+ *
+ * This filter provides similar functionality to the otb::FastCorrelationImageFilter.
+ * The otb::FastCorrelationImageFilter provides an optimized implementation of fine registration for correlation
+ * metric. It is faster but less flexible: images should have the same size, and the metric can not be changed.
+ *
+ * The FineRegistrationImageFilter allows to use the full range of itk::ImageToImageMetric provided by itk.
+ *
+ * \sa      FastCorrelationImageFilter
+ * \ingroup IntensityImageFilters, 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

Code/DisparityMap/otbFineRegistrationImageFilter.txx

+/*=========================================================================
+
+  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"
+#include "itkExceptionObject.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]);
+
+        try
+        {
+          // compute currentMetric
+          currentMetric = m_Metric->GetValue(params);
+
+          // Check for maximum
+          if((m_Minimize && (currentMetric < optMetric)) || (!m_Minimize && (currentMetric > optMetric)))
+            {
+            optMetric = currentMetric;
+            optParams = params;
+            }
+        }
+        catch(itk::ExceptionObject& err)
+        {
+          itkWarningMacro(<<err.GetDescription());
+        }
+        }
+      }
+
+    // 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]);
+
+          try
+          {
+            // compute currentMetric
+            currentMetric = m_Metric->GetValue(params);
+
+            // Check for maximum
+            if((m_Minimize && (currentMetric < optMetric)) || (!m_Minimize && (currentMetric > optMetric)))
+              {
+              optMetric = currentMetric;
+              optParams = params;
+              }
+          }
+          catch(itk::ExceptionObject& err)
+          {
+            itkWarningMacro(<<err.GetDescription());
+
+          }
+          }
+        }
+      }
+
+    // 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

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_10} 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_10} 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_10} 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_10} 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")

Testing/Code/DisparityMap/otbDisparityMapTests3.cxx

@@ -28,4 +28,6 @@ void RegisterTests()
 {
   REGISTER_TEST(otbStreamingWarpImageFilterNew);
   REGISTER_TEST(otbStreamingWarpImageFilter);
+  REGISTER_TEST(otbFineRegistrationImageFilterNew);
+  REGISTER_TEST(otbFineRegistrationImageFilterTest);
 }

Testing/Code/DisparityMap/otbFineRegistrationImageFilterNew.cxx

+/*=========================================================================
+
+  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;
+}

Testing/Code/DisparityMap/otbFineRegistrationImageFilterTest.cxx

+/*=========================================================================
+
+  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;
+}