BUG: #2046 Add an epsilon margin to compute the m_ReachableOutputRegion

The application Superimpose crashed with a segmentation fault during
resampling. Investigation revealed that the interpolated index was out
of image bounds. The application uses the PHR mode and resample the XS
image with a otb::GridResampleImageFilter. This filter doesn't check
explicitly that input indexes are in the input buffered region. It
uses a m_ReachableOutputRegion to crop the output region
processed. The problem appears because an output pixel is right on the
border of the input image extent, and when someone uses a nearest neighbor
interpolator, this one converts the continuous index to an
out-of-bound index.

Modules/Applications/AppProjection/app/otbSuperimpose.cxx

@@ -192,6 +192,7 @@ private:
       NNInterpolatorType::Pointer interpolator = NNInterpolatorType::New();
       m_Resampler->SetInterpolator(interpolator);
       m_BasicResampler->SetInterpolator(interpolator);
+      m_BasicResampler->SetInterpolationMargin(1e-9);
     }
     break;
     case Interpolator_BCO:

Modules/Applications/AppProjection/test/CMakeLists.txt

@@ -367,3 +367,13 @@ otb_test_application(NAME apTvPrSuperimpose_phr
                      VALID  --compare-image ${EPSILON_7}
                         ${BASELINE}/apTvPrSuperimposePHR.tif
                         ${TEMP}/apTvPrSuperimposePHR.tif)
+
+otb_test_application(NAME apTvPrSuperimpose_phr_nn
+                     APP Superimpose
+                     OPTIONS -inr  ${INPUTDATA}/phr_pan.tif
+                             -inm ${INPUTDATA}/phr_xs.tif
+                             -interpolator nn
+                             -out ${TEMP}/apTvPrSuperimposePHR_nn.tif int16
+                     VALID  --compare-image ${EPSILON_7}
+                        ${BASELINE}/apTvPrSuperimposePHR_nn.tif
+                        ${TEMP}/apTvPrSuperimposePHR_nn.tif)

Modules/Filtering/ImageManipulation/include/otbGridResampleImageFilter.h

@@ -124,6 +124,9 @@ public:
   itkGetMacro(CheckOutputBounds, bool);
   itkBooleanMacro(CheckOutputBounds);
 
+  itkSetMacro(InterpolationMargin, double);
+  itkGetMacro(InterpolationMargin, double);
+
   itkSetObjectMacro(Interpolator, InterpolatorType);
   itkGetObjectMacro(Interpolator, InterpolatorType);
 
@@ -190,6 +193,8 @@ private:
 
   OutputPixelType m_EdgePaddingValue; // Default pixel value
 
+  double m_InterpolationMargin;
+
   bool m_CheckOutputBounds; // Shall we check
                             // output bounds when
                             // casting?

Modules/Filtering/ImageManipulation/include/otbGridResampleImageFilter.hxx

@@ -42,6 +42,7 @@ GridResampleImageFilter<TInputImage, TOutputImage, TInterpolatorPrecision>::Grid
     m_OutputSpacing(),
     m_EdgePaddingValue(),
     m_CheckOutputBounds(true),
+    m_InterpolationMargin(0.0),
     m_Interpolator(),
     m_ReachableOutputRegion()
 {
@@ -227,7 +228,6 @@ void GridResampleImageFilter<TInputImage, TOutputImage, TInterpolatorPrecision>:
 
   // Compute the padding due to the interpolator
 
-
   IndexType inUL = this->GetInput()->GetBufferedRegion().GetIndex();
   IndexType inLR = this->GetInput()->GetBufferedRegion().GetIndex() + this->GetInput()->GetBufferedRegion().GetSize();
   inLR[0] -= 1;
@@ -246,9 +246,9 @@ void GridResampleImageFilter<TInputImage, TOutputImage, TInterpolatorPrecision>:
   this->GetInput()->TransformIndexToPhysicalPoint(inUL, inULp);
   this->GetInput()->TransformIndexToPhysicalPoint(inLR, inLRp);
 
-  inULp -= 0.5 * this->GetInput()->GetSignedSpacing();
-  inLRp += 0.5 * this->GetInput()->GetSignedSpacing();
-
+  inULp -= (0.5 - m_InterpolationMargin) * this->GetInput()->GetSignedSpacing();
+  inLRp += (0.5 - m_InterpolationMargin) * this->GetInput()->GetSignedSpacing();
+ 
   ContinuousInputIndexType outUL;
   ContinuousInputIndexType outLR;
   this->GetOutput()->TransformPhysicalPointToContinuousIndex(inULp, outUL);
@@ -265,6 +265,8 @@ void GridResampleImageFilter<TInputImage, TOutputImage, TInterpolatorPrecision>:
 
   m_ReachableOutputRegion.SetIndex(outputIndex);
   m_ReachableOutputRegion.SetSize(outputSize);
+
+  otbMsgDevMacro(<< "ReachableOutputRegion: " << m_ReachableOutputRegion);
 }
 
 template <typename TInputImage, typename TOutputImage, typename TInterpolatorPrecision>