diff --git a/Code/DisparityMap/otbFineRegistrationImageFilter.h b/Code/DisparityMap/otbFineRegistrationImageFilter.h
index 929df1ef294dc4e7417711c2edc0ed7675f76824..fc15da4382cb6a783d7a708312e640a16d83e519 100644
--- a/Code/DisparityMap/otbFineRegistrationImageFilter.h
+++ b/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
diff --git a/Code/DisparityMap/otbFineRegistrationImageFilter.txx b/Code/DisparityMap/otbFineRegistrationImageFilter.txx
index 0ecb07fa55cc0ddeed45a6090af3a567dd8c1d28..e93a2e3a351ef0a67ec163191a4921e384f7bcea 100644
--- a/Code/DisparityMap/otbFineRegistrationImageFilter.txx
+++ b/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
diff --git a/Testing/Code/DisparityMap/CMakeLists.txt b/Testing/Code/DisparityMap/CMakeLists.txt
index 72f682a77c6779ea539a2f8eee766647f8662b41..9e5764970e0bda1aa5a4968ab740e1ce4b61bfc4 100644
--- a/Testing/Code/DisparityMap/CMakeLists.txt
+++ b/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
)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
diff --git a/Testing/Code/DisparityMap/otbFineRegistrationImageFilterTest.cxx b/Testing/Code/DisparityMap/otbFineRegistrationImageFilterTest.cxx
index b49b9b2a9eb741bb9a5a7f068cd3da17b225830f..a06714c6b6d73a7f9c2a1dc15fcc4f365ce613e0 100644
--- a/Testing/Code/DisparityMap/otbFineRegistrationImageFilterTest.cxx
+++ b/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)
{