diff --git a/Code/DisparityMap/otbStreamingWarpImageFilter.txx b/Code/DisparityMap/otbStreamingWarpImageFilter.txx
index 29d4c9a13f17ecd046d6d8291640c0386db5ddcf..ed4426fc2a0b7e4a9ae84c86881f385c5f94ba52 100644
--- a/Code/DisparityMap/otbStreamingWarpImageFilter.txx
+++ b/Code/DisparityMap/otbStreamingWarpImageFilter.txx
@@ -23,6 +23,7 @@
#include "otbStreamingWarpImageFilter.h"
#include "otbStreamingTraits.h"
+#include "itkImageRegionIteratorWithIndex.h"
namespace otb
{
@@ -30,51 +31,152 @@ namespace otb
template<class TInputImage, class TOutputImage, class TDeformationField>
StreamingWarpImageFilter<TInputImage, TOutputImage, TDeformationField>
::StreamingWarpImageFilter()
-{
+ {
// Fill the default maximum deformation
m_MaximumDeformation.Fill(1);
-}
+ }
template<class TInputImage, class TOutputImage, class TDeformationField>
void
StreamingWarpImageFilter<TInputImage, TOutputImage, TDeformationField>
::GenerateInputRequestedRegion()
-{
- // First, ca the superclass implementation
- Superclass::GenerateInputRequestedRegion();
-
+ {
// Get the input and deformation field pointers
- InputImageType * inputPtr = const_cast<InputImageType *>(this->GetInput());
- DeformationFieldType * deformationPtr = const_cast<DeformationFieldType*>(this->GetDeformationField());
-
+ InputImageType * inputPtr
+ = const_cast<InputImageType *>(this->GetInput());
+ DeformationFieldType * deformationPtr
+ = const_cast<DeformationFieldType*>(this->GetDeformationField());
+ OutputImageType * outputPtr
+ = this->GetOutput();
// Check if the input and the deformation field exist
- if (!inputPtr || !deformationPtr)
+ if (!inputPtr || !deformationPtr || !outputPtr)
{
return;
}
- // Compute the security margin radius
- typename InputImageType::SizeType radius;
- typename InputImageType::SpacingType spacing = inputPtr->GetSpacing();
+ // Here we are breaking traditional pipeline steps because we need to access the deformation field data
+ // so as to compute the input image requested region
+
+ // 1) First, evaluate the deformation field requested region corresponding to the output requested region
+ // (Here we suppose that the deformation field and the output image are in the same geometry/map projection)
+ typename OutputImageType::RegionType outputRequestedRegion = outputPtr->GetRequestedRegion();
+ typename OutputImageType::IndexType outIndexStart = outputRequestedRegion.GetIndex();
+ typename OutputImageType::IndexType outIndexEnd;
+ for(unsigned int dim = 0; dim<OutputImageType::ImageDimension;++dim)
+ outIndexEnd[dim]= outIndexStart[dim] + outputRequestedRegion.GetSize()[dim]-1;
+ typename OutputImageType::PointType outPointStart,outPointEnd;
+ outputPtr->TransformIndexToPhysicalPoint(outIndexStart,outPointStart);
+ outputPtr->TransformIndexToPhysicalPoint(outIndexEnd,outPointEnd);
+ // Here min/max allows for inverted spacing in deformation field
+ typename DeformationFieldType::PointType defPointStart, defPointEnd;
+ for(unsigned int dim = 0; dim<OutputImageType::ImageDimension;++dim)
+ {
+ defPointStart[dim] = std::min(outPointStart[dim], outPointEnd[dim]);
+ defPointEnd[dim] = std::max(outPointStart[dim],outPointEnd[dim]);
+ }
+ typename DeformationFieldType::IndexType defIndexStart,defIndexEnd;
+ deformationPtr->TransformPhysicalPointToIndex(defPointStart,defIndexStart);
+ deformationPtr->TransformPhysicalPointToIndex(defPointEnd,defIndexEnd);
+ typename DeformationFieldType::SizeType defRequestedSize;
+ typename DeformationFieldType::IndexType defRequestedIndex;
+ for(unsigned int dim = 0; dim<OutputImageType::ImageDimension;++dim)
+ {
+ defRequestedIndex[dim] = std::min(defIndexStart[dim],defIndexEnd[dim]);
+ defRequestedSize[dim] = std::max(defIndexStart[dim],defIndexEnd[dim]) - defRequestedIndex[dim] + 1;
+ }
- // Add the padding due to the interpolator
- unsigned int interpolatorRadius =
- StreamingTraits<typename Superclass::InputImageType>::CalculateNeededRadiusForInterpolator(this->GetInterpolator());
+ // Finally, build the deformation field requested region
+ typename DeformationFieldType::RegionType deformationRequestedRegion;
+ deformationRequestedRegion.SetIndex(defRequestedIndex);
+ deformationRequestedRegion.SetSize(defRequestedSize);
- // Compute the margin due to the maximum deformation value and interpolator radius
- for (unsigned int i = 0; i < InputImageType::ImageDimension; ++i)
+ // crop the input requested region at the input's largest possible region
+ if (deformationRequestedRegion.Crop(deformationPtr->GetLargestPossibleRegion()))
{
- radius[i] = interpolatorRadius + static_cast<unsigned int>(vcl_ceil(m_MaximumDeformation[i] / vcl_abs(spacing[i])));
+ deformationPtr->SetRequestedRegion(deformationRequestedRegion);
}
+ else
+ {
+ // Couldn't crop the region (requested region is outside the largest
+ // possible region). Throw an exception.
- otbMsgDevMacro(
- << "WarpImageFilter: MaximumDeformation: " << m_MaximumDeformation << ", interpolator radius: " <<
- interpolatorRadius << ", total radius: " << radius);
+ // store what we tried to request (prior to trying to crop)
+ deformationPtr->SetRequestedRegion(deformationRequestedRegion);
+
+ // build an exception
+ itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
+ e.SetLocation(ITK_LOCATION);
+ e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
+ e.SetDataObject(inputPtr);
+ throw e;
+ }
- // get a copy of the input requested region (should equal the output
- // requested region)
- typename DeformationFieldType::RegionType inputRequestedRegion;
- inputRequestedRegion = deformationPtr->GetRequestedRegion();
+ // 2) If we are still there, we have a correct deformation field requested region.
+ // This next step breaks pipeline rule but we need to do it to compute the input requested region,
+ // since it depends on deformation value.
+
+ // Trigger pipeline update on the deformation field
+
+ deformationPtr->PropagateRequestedRegion();
+ deformationPtr->UpdateOutputData();
+
+ // Walk the loaded deformation field to derive maximum and minimum
+ itk::ImageRegionIteratorWithIndex<DeformationFieldType> defIt(deformationPtr,deformationRequestedRegion);
+ defIt.GoToBegin();
+
+ typename InputImageType::PointType currentPoint;
+ typename InputImageType::PointType inputStartPoint, inputEndPoint;
+
+ // Initialise start and end points
+ deformationPtr->TransformIndexToPhysicalPoint(defIt.GetIndex(),currentPoint);
+ for(unsigned int dim = 0; dim<DeformationFieldType::ImageDimension;++dim)
+ {
+ currentPoint[dim]+=defIt.Get()[dim];
+ }
+ inputStartPoint = currentPoint;
+ inputEndPoint = currentPoint;
+
+ ++defIt;
+
+ // 3) Now Walk the field and compute the physical bounding box
+ while(!defIt.IsAtEnd())
+ {
+ deformationPtr->TransformIndexToPhysicalPoint(defIt.GetIndex(),currentPoint);
+ for(unsigned int dim = 0; dim<DeformationFieldType::ImageDimension;++dim)
+ {
+ currentPoint[dim]+=defIt.Get()[dim];
+ if(inputStartPoint[dim] > currentPoint[dim])
+ inputStartPoint[dim] = currentPoint[dim];
+ if(inputEndPoint[dim] < currentPoint[dim])
+ inputEndPoint[dim] = currentPoint[dim];
+ }
+ ++defIt;
+ }
+
+ // Convert physical bouding box to requested region
+ typename InputImageType::IndexType inputStartIndex,inputEndIndex;
+ inputPtr->TransformPhysicalPointToIndex(inputStartPoint,inputStartIndex);
+ inputPtr->TransformPhysicalPointToIndex(inputEndPoint,inputEndIndex);
+
+ typename InputImageType::SizeType inputFinalSize;
+ typename InputImageType::IndexType inputFinalIndex;
+
+ for(unsigned int dim = 0; dim<DeformationFieldType::ImageDimension;++dim)
+ {
+ inputFinalIndex[dim] = std::min(inputStartIndex[dim],inputEndIndex[dim]);
+ inputFinalSize[dim] = std::max(inputStartIndex[dim],inputEndIndex[dim])-inputFinalIndex[dim]+1;
+ }
+
+ typename InputImageType::RegionType inputRequestedRegion;
+ inputRequestedRegion.SetIndex(inputFinalIndex);
+ inputRequestedRegion.SetSize(inputFinalSize);
+
+ // Compute the security margin radius
+ typename InputImageType::SizeType radius;
+
+ // Compute the padding due to the interpolator
+ unsigned int interpolatorRadius =
+ StreamingTraits<typename Superclass::InputImageType>::CalculateNeededRadiusForInterpolator(this->GetInterpolator());
// pad the input requested region by the operator radius
inputRequestedRegion.PadByRadius(radius);
@@ -83,7 +185,6 @@ StreamingWarpImageFilter<TInputImage, TOutputImage, TDeformationField>
if (inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()))
{
inputPtr->SetRequestedRegion(inputRequestedRegion);
- return;
}
else
{
@@ -100,16 +201,16 @@ StreamingWarpImageFilter<TInputImage, TOutputImage, TDeformationField>
e.SetDataObject(inputPtr);
throw e;
}
-}
+ }
template<class TInputImage, class TOutputImage, class TDeformationField>
void
StreamingWarpImageFilter<TInputImage, TOutputImage, TDeformationField>
::PrintSelf(std::ostream& os, itk::Indent indent) const
-{
+ {
Superclass::PrintSelf(os, indent);
os << indent << "Maximum deformation: " << m_MaximumDeformation << std::endl;
-}
+ }
} // end namespace otb
#endif