ENH: Updating the FineRegistrationFilter to take into account a coarser grid...

ENH: Updating the FineRegistrationFilter to take into account a coarser grid step, an initial offset, and moving/fixed images with different size and resolutions

Code/DisparityMap/otbFineRegistrationImageFilter.h

@@ -37,7 +37,7 @@ namespace otb
  *
  * 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.
+ * minimum or maximum depending on the metric (default is On).
  *
  * 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.
@@ -46,11 +46,17 @@ namespace otb
  * the GetOutputDeformationField() method returns the corresponding offset.
  *
  * If the UseSpacingOn() flag is used, the output deformation field takes the input image spacing into account.
- * otherwise, the deformation field is expressed in pixels.
+ * otherwise, the deformation field is expressed in pixels (default ins On).
  *
- * 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.
+ * This filter accepts fixed and moving images with different sizes and spacing. Metric and search windows radius
+ * are expressed in terms of number of pixels in the fixed image.
+ *
+ * An initial offset can be used to reduce computation time in case of input and moving images with a significant
+ * offset. This offset is taken into account in the output deformation field.
+ *
+ * It is possible to generate an output metric map and deformation field at a coarser resolution by setting
+ * grid step to value higher than 1 (grid step is expressed in terms of number of fixed image pixels).
+ * Default value is 1.
  *
  * The FineRegistrationImageFilter allows to use the full range of itk::ImageToImageMetric provided by itk.
  *
@@ -82,6 +88,8 @@ public:
   typedef typename TInputImage::SizeType                          SizeType;
   typedef typename TInputImage::IndexType                         IndexType;
   typedef typename TInputImage::SpacingType                       SpacingType;
+  typedef typename TInputImage::PointType                         PointType;
+  typedef typename TInputImage::OffsetType                        OffsetType;
   typedef typename itk::InterpolateImageFunction
   <TInputImage, double>                                           InterpolatorType;
   typedef typename InterpolatorType::Pointer                      InterpolatorPointerType;
@@ -133,18 +141,32 @@ public:
   itkSetMacro(UseSpacing,bool);
   itkBooleanMacro(UseSpacing);
 
+  /** Set default offset between the two images */
+  itkSetMacro(InitialOffset,SpacingType);
+  itkGetConstReferenceMacro(InitialOffset,SpacingType);
+
+  /** Set the grid step */
+  itkSetMacro(GridStep,OffsetType);
+  itkGetConstReferenceMacro(GridStep,OffsetType);
+
   /** Set unsigned int radius */
   void SetRadius(unsigned int radius)
   {
     m_Radius.Fill(radius);
   }
   
-/** Set unsigned int radius */
+ /** Set unsigned int radius */
   void SetSearchRadius(unsigned int radius)
   {
     m_SearchRadius.Fill(radius);
   }
 
+  /** Set unsigned int grid step */
+  void SetGridStep(unsigned int step)
+  {
+    m_GridStep.Fill(step);
+  }
+
 protected:
   /** Constructor */
   FineRegistrationImageFilter();
@@ -157,6 +179,9 @@ protected:
   /** Generate the input requested regions  */
   virtual void GenerateInputRequestedRegion(void);
 
+  /** Generate output information */
+  virtual void GenerateOutputInformation(void);
+
 private:
   FineRegistrationImageFilter(const Self&); //purposely not implemented
   void operator=(const Self&); //purposely not implemented
@@ -185,6 +210,12 @@ private:
   /** The translation */
   TranslationPointerType        m_Translation;
 
+  /** Default offset */
+  SpacingType                   m_InitialOffset;
+
+  /** Grid step */
+  OffsetType                    m_GridStep;
+
 };
 
 } // end namespace otb

Code/DisparityMap/otbFineRegistrationImageFilter.txx

@@ -59,6 +59,12 @@ FineRegistrationImageFilter<TInputImage,T0utputCorrelation,TOutputDeformationFie
 
   // Translation
   m_Translation = TranslationType::New();
+
+  // Grid Step
+  m_GridStep.Fill(1);
+
+  // Default offset
+  m_InitialOffset.Fill(0);
  }
 
 template <class TInputImage, class T0utputCorrelation, class TOutputDeformationField>
@@ -115,6 +121,41 @@ FineRegistrationImageFilter<TInputImage,T0utputCorrelation,TOutputDeformationFie
   return static_cast<TOutputDeformationField *>(this->itk::ProcessObject::GetOutput(1));
  }
 
+template <class TInputImage, class TOutputCorrelation, class TOutputDeformationField>
+void
+FineRegistrationImageFilter<TInputImage,TOutputCorrelation,TOutputDeformationField>
+::GenerateOutputInformation()
+ {
+  // Call superclass implementation
+  Superclass::GenerateOutputInformation();
+
+  // Retrieve output pointers
+  TOutputCorrelation * outputPtr = this->GetOutput();
+  TOutputDeformationField *outputFieldPtr = this->GetOutputDeformationField();
+
+  // Update size and spacing according to grid step
+  InputImageRegionType largestRegion  = outputPtr->GetLargestPossibleRegion();
+  SizeType outputSize       = largestRegion.GetSize();
+  SpacingType outputSpacing = outputPtr->GetSpacing();
+
+  for(unsigned int dim = 0; dim < TOutputCorrelation::ImageDimension;++dim)
+    {
+    outputSize[dim] /= m_GridStep[dim];
+    outputSpacing[dim] *= m_GridStep[dim];
+    }
+
+  // Set spacing
+  outputPtr->SetSpacing(outputSpacing);
+  outputFieldPtr->SetSpacing(outputSpacing);
+
+  // Set largest region size
+  largestRegion.SetSize(outputSize);
+  outputPtr->SetLargestPossibleRegion(largestRegion);
+  outputFieldPtr->SetLargestPossibleRegion(largestRegion);
+
+  std::cout<<"Output region: "<<largestRegion<<std::endl;
+  std::cout<<"Output spacing:"<<outputSpacing<<std::endl;
+ }
 
 template <class TInputImage, class TOutputCorrelation, class TOutputDeformationField>
 void
@@ -137,17 +178,68 @@ FineRegistrationImageFilter<TInputImage,TOutputCorrelation,TOutputDeformationFie
 
   // get a copy of the fixed requested region (should equal the output
   // requested region)
-  typename TInputImage::RegionType fixedRequestedRegion;
+  InputImageRegionType fixedRequestedRegion, movingRequestedRegion;
   fixedRequestedRegion = fixedPtr->GetRequestedRegion();
+
+  // Apply grid step
+  SizeType fixedRequestedSize = fixedRequestedRegion.GetSize();
+  IndexType fixedRequestedIndex = fixedRequestedRegion.GetIndex();
+
+  for(unsigned int dim = 0; dim < TOutputCorrelation::ImageDimension;++dim)
+      {
+      fixedRequestedSize [dim] *= m_GridStep[dim];
+      fixedRequestedIndex[dim] *= m_GridStep[dim];
+      }
+
+  fixedRequestedRegion.SetSize(fixedRequestedSize);
+  fixedRequestedRegion.SetIndex(fixedRequestedIndex);
+
   // pad the input requested region by the operator radius
   fixedRequestedRegion.PadByRadius( m_Radius );
 
+  std::cout<<"Fixed requested region: "<<fixedRequestedRegion<<std::endl;
+
   // 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 );
+  InputImageRegionType searchFixedRequestedRegion = fixedRequestedRegion;
+  searchFixedRequestedRegion.PadByRadius(m_Radius);
+
+
+  // Find corners of the search window
+   IndexType ulIndex = searchFixedRequestedRegion.GetIndex();
+
+   IndexType lrIndex;
+   for(unsigned int dim = 0; dim < TInputImage::ImageDimension;++dim)
+     {
+     lrIndex[dim]= searchFixedRequestedRegion.GetIndex()[dim]
+                 + searchFixedRequestedRegion.GetSize()[dim]-1;
+     }
+
+   // Transform to physical space
+   PointType ulPoint, lrPoint;
+   fixedPtr->TransformIndexToPhysicalPoint(lrIndex,lrPoint);
+   fixedPtr->TransformIndexToPhysicalPoint(ulIndex,ulPoint);
+
+   // Apply default offset
+   lrPoint+=m_InitialOffset;
+   ulPoint+=m_InitialOffset;
+
+   // Transform back into moving region index space
+   IndexType movingIndex1, movingIndex2, movingIndex;
+   movingPtr->TransformPhysicalPointToIndex(ulPoint,movingIndex1);
+   movingPtr->TransformPhysicalPointToIndex(lrPoint,movingIndex2);
+
+   // Find requested region
+   SizeType movingSize;
+
+   for(unsigned int dim = 0; dim < TInputImage::ImageDimension;++dim)
+     {
+       movingIndex[dim] = std::min(movingIndex1[dim],movingIndex2[dim]);
+       movingSize[dim] = std::max(movingIndex1[dim],movingIndex2[dim]) - movingIndex[dim] + 1;
+     }
+
+   movingRequestedRegion.SetIndex(movingIndex);
+   movingRequestedRegion.SetSize(movingSize);
 
   // crop the fixed region at the fixed's largest possible region
   if ( fixedRequestedRegion.Crop(fixedPtr->GetLargestPossibleRegion()))
@@ -180,21 +272,16 @@ FineRegistrationImageFilter<TInputImage,TOutputCorrelation,TOutputDeformationFie
   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 );
+    // possible region). This case might happen so we do not throw any exception but
+    // request a null region instead
+    movingSize.Fill(0);
+    movingRequestedRegion.SetSize(movingSize);
+    movingIndex.Fill(0);
+    movingRequestedRegion.SetIndex(movingIndex);
 
-    // 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;
+    // store what we tried to request (prior to trying to crop)
+    movingPtr->SetRequestedRegion(movingRequestedRegion);
     }
-
   return;
  }
 
@@ -261,7 +348,14 @@ FineRegistrationImageFilter<TInputImage,TOutputCorrelation,TOutputDeformationFie
 
     // Build the region on which to compute the currentMetric
     InputImageRegionType currentMetricRegion;
-    currentMetricRegion.SetIndex(outputIt.GetIndex());
+
+    // Apply grid step
+    IndexType currentIndex = outputIt.GetIndex();
+    for(unsigned int dim = 0; dim < TInputImage::ImageDimension;++dim)
+      {
+      currentIndex[dim] *= m_GridStep[dim];
+      }
+    currentMetricRegion.SetIndex(currentIndex);
     SizeType size;
     size.Fill(1);
     currentMetricRegion.SetSize(size);
@@ -275,8 +369,8 @@ FineRegistrationImageFilter<TInputImage,TOutputCorrelation,TOutputDeformationFie
       {
       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]);
+        params[0]=m_InitialOffset[0] + static_cast<double>(i*fixedSpacing[0]);
+        params[1]=m_InitialOffset[0] + static_cast<double>(j*fixedSpacing[1]);
 
         try
         {
@@ -346,7 +440,7 @@ FineRegistrationImageFilter<TInputImage,TOutputCorrelation,TOutputDeformationFie
     else
       {
       deformationValue[0] = optParams[0]/fixedSpacing[0];
-      deformationValue[1] = optParams[1]/FixedSpacing[1];
+      deformationValue[1] = optParams[1]/fixedSpacing[1];
       }
     outputDfIt.Set(deformationValue);
     // Update iterators

Testing/Code/DisparityMap/CMakeLists.txt

@@ -333,6 +333,9 @@ ADD_TEST(feTvFineRegistrationImageFilterTestWithCorrelation ${DISPARITYMAP_TESTS
 	    2 # sradius
 	    0.1 # precision
 	    0 # Correlation
+	    1 # Grid step
+	    0 # Initial offset x
+	    0 # Initial offset y
 )
 
 ADD_TEST(feTvFineRegistrationImageFilterTestWithNormalizedCorrelation ${DISPARITYMAP_TESTS3}
@@ -349,7 +352,10 @@ ADD_TEST(feTvFineRegistrationImageFilterTestWithNormalizedCorrelation ${DISPARIT
         3 # radius
 	    2 # sradius
 	    0.1 # precision
-	    1 # Normalized Correlation
+	    1 # Normalized Correlation    
+	    1 # Grid step
+	    0 # Initial offset x
+	    0 # Initial offset y
 )
 
 ADD_TEST(feTvFineRegistrationImageFilterTestWithMeanSquare ${DISPARITYMAP_TESTS3}
@@ -367,6 +373,9 @@ ADD_TEST(feTvFineRegistrationImageFilterTestWithMeanSquare ${DISPARITYMAP_TESTS3
 	    2 # sradius
 	    0.1 # precision
 	    2 # Mean square
+	    1 # Grid step
+	    0 # Initial offset x
+	    0 # Initial offset y
 )
 
 ADD_TEST(feTvFineRegistrationImageFilterTestWithMeanReciprocalDifference ${DISPARITYMAP_TESTS3}
@@ -384,6 +393,9 @@ ADD_TEST(feTvFineRegistrationImageFilterTestWithMeanReciprocalDifference ${DISPA
 	    2 # sradius
 	    0.1 # precision
 	    3 # Mean reciprocal difference
+	    1 # Grid step
+	    0 # Initial offset x
+	    0 # Initial offset y
 )
 
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Testing/Code/DisparityMap/otbFineRegistrationImageFilterTest.cxx

@@ -33,10 +33,11 @@
 
 int otbFineRegistrationImageFilterTest( int argc, char * argv[] )
 {
-  if(argc!=9)
+  if(argc!=12)
     {
     std::cerr<<"Usage: "<<argv[0]<<" fixed_fname moving_fname output_correl output_field radius search_radius ";
-    std::cerr<<"subpixPrecision metric(0=CC,1=NCC,2=MeanSquare,3=Mean reciprocal square difference)"<<std::endl;
+    std::cerr<<"subpixPrecision metric(0=CC,1=NCC,2=MeanSquare,3=Mean reciprocal square difference) ";
+    std::cerr<<"gridStep offsetX offsetY"<<std::endl;
     return EXIT_FAILURE;
     }
   const char * fixedFileName  = argv[1];
@@ -47,6 +48,9 @@ int otbFineRegistrationImageFilterTest( int argc, char * argv[] )
   const unsigned int sradius  = atoi(argv[6]);
   const double      precision = atof(argv[7]);
   const unsigned int metric   = atoi(argv[8]);
+  const unsigned int gridStep = atoi(argv[9]);
+  const double       offsetx  = atof(argv[10]);
+  const double       offsety  = atof(argv[11]);
 
   typedef double      PixelType;
   const unsigned int  Dimension = 2;
@@ -73,6 +77,13 @@ int otbFineRegistrationImageFilterTest( int argc, char * argv[] )
   registration->SetRadius(radius);
   registration->SetSearchRadius(sradius);
   registration->SetSubPixelAccuracy(precision);
+  registration->SetGridStep(gridStep);
+
+  RegistrationFilterType::SpacingType offset;
+  offset[0] = offsetx;
+  offset[1] = offsety;
+
+  registration->SetInitialOffset(offset);
 
   switch(metric)
   {