diff --git a/Code/DisparityMap/otbDisparityMapToDEMFilter.txx b/Code/DisparityMap/otbDisparityMapToDEMFilter.txx
index b4a8b44f18b4754422f9f222d803170bb1eace74..79e1709b7f7d8643c901c80683469fe71b89fd2c 100644
--- a/Code/DisparityMap/otbDisparityMapToDEMFilter.txx
+++ b/Code/DisparityMap/otbDisparityMapToDEMFilter.txx
@@ -263,39 +263,54 @@ DisparityMapToDEMFilter<TDisparityImage,TInputImage,TOutputDEMImage,TEpipolarGri
// left image
typename SensorImageType::SizeType inputSize = leftImgPtr->GetLargestPossibleRegion().GetSize();
- typename SensorImageType::PointType tmpPoint;
- tmpPoint = leftImgPtr->GetOrigin();
- RSTransform2DType::OutputPointType left_ul = leftToGroundTransform->TransformPoint(tmpPoint);
-
- tmpPoint[0] = (leftImgPtr->GetOrigin())[0] + (leftImgPtr->GetSpacing())[0] * static_cast<double>(inputSize[0]);
- tmpPoint[1] = (leftImgPtr->GetOrigin())[1];
- RSTransform2DType::OutputPointType left_ur = leftToGroundTransform->TransformPoint(tmpPoint);
-
- tmpPoint[0] = (leftImgPtr->GetOrigin())[0] + (leftImgPtr->GetSpacing())[0] * static_cast<double>(inputSize[0]);
- tmpPoint[1] = (leftImgPtr->GetOrigin())[1] + (leftImgPtr->GetSpacing())[1] * static_cast<double>(inputSize[1]);
- RSTransform2DType::OutputPointType left_lr = leftToGroundTransform->TransformPoint(tmpPoint);
-
- tmpPoint[0] = (leftImgPtr->GetOrigin())[0];
- tmpPoint[1] = (leftImgPtr->GetOrigin())[1] + (leftImgPtr->GetSpacing())[1] * static_cast<double>(inputSize[1]);
- RSTransform2DType::OutputPointType left_ll = leftToGroundTransform->TransformPoint(tmpPoint);
+ typename SensorImageType::PointType ulp, urp, llp, lrp;
+ itk::ContinuousIndex<double,2> ul(leftImgPtr->GetLargestPossibleRegion().GetIndex());
+ ul[0] += -0.5;
+ ul[1] += -0.5;
+
+ itk::ContinuousIndex<double,2> ur(ul);
+ itk::ContinuousIndex<double,2> lr(ul);
+ itk::ContinuousIndex<double,2> ll(ul);
+ ur[0] += static_cast<double>(inputSize[0]);
+ lr[0] += static_cast<double>(inputSize[0]);
+ lr[1] += static_cast<double>(inputSize[1]);
+ ll[1] += static_cast<double>(inputSize[1]);
+
+ leftImgPtr->TransformContinuousIndexToPhysicalPoint(ul,ulp);
+ leftImgPtr->TransformContinuousIndexToPhysicalPoint(ur,urp);
+ leftImgPtr->TransformContinuousIndexToPhysicalPoint(ll,llp);
+ leftImgPtr->TransformContinuousIndexToPhysicalPoint(lr,lrp);
+
+ RSTransform2DType::OutputPointType left_ul, left_ur, left_ll, left_lr;
+ left_ul = leftToGroundTransform->TransformPoint(ulp);
+ left_ur = leftToGroundTransform->TransformPoint(urp);
+ left_ll = leftToGroundTransform->TransformPoint(llp);
+ left_lr = leftToGroundTransform->TransformPoint(lrp);
// right image
inputSize = rightImgPtr->GetLargestPossibleRegion().GetSize();
- tmpPoint = rightImgPtr->GetOrigin();
- RSTransform2DType::OutputPointType right_ul = rightToGroundTransform->TransformPoint(tmpPoint);
-
- tmpPoint[0] = (rightImgPtr->GetOrigin())[0] + (rightImgPtr->GetSpacing())[0] * static_cast<double>(inputSize[0]);
- tmpPoint[1] = (rightImgPtr->GetOrigin())[1];
- RSTransform2DType::OutputPointType right_ur = rightToGroundTransform->TransformPoint(tmpPoint);
-
- tmpPoint[0] = (rightImgPtr->GetOrigin())[0] + (rightImgPtr->GetSpacing())[0] * static_cast<double>(inputSize[0]);
- tmpPoint[1] = (rightImgPtr->GetOrigin())[1] + (rightImgPtr->GetSpacing())[1] * static_cast<double>(inputSize[1]);
- RSTransform2DType::OutputPointType right_lr = rightToGroundTransform->TransformPoint(tmpPoint);
-
- tmpPoint[0] = (rightImgPtr->GetOrigin())[0];
- tmpPoint[1] = (rightImgPtr->GetOrigin())[1] + (rightImgPtr->GetSpacing())[1] * static_cast<double>(inputSize[1]);
- RSTransform2DType::OutputPointType right_ll = rightToGroundTransform->TransformPoint(tmpPoint);
-
+ ul = rightImgPtr->GetLargestPossibleRegion().GetIndex();
+ ul[0] += -0.5;
+ ul[1] += -0.5;
+ ur = ul;
+ lr = ul;
+ ll = ul;
+ ur[0] += static_cast<double>(inputSize[0]);
+ lr[0] += static_cast<double>(inputSize[0]);
+ lr[1] += static_cast<double>(inputSize[1]);
+ ll[1] += static_cast<double>(inputSize[1]);
+
+ rightImgPtr->TransformContinuousIndexToPhysicalPoint(ul,ulp);
+ rightImgPtr->TransformContinuousIndexToPhysicalPoint(ur,urp);
+ rightImgPtr->TransformContinuousIndexToPhysicalPoint(ll,llp);
+ rightImgPtr->TransformContinuousIndexToPhysicalPoint(lr,lrp);
+
+ RSTransform2DType::OutputPointType right_ul, right_ur, right_lr, right_ll;
+ right_ul = rightToGroundTransform->TransformPoint(ulp);
+ right_ur = rightToGroundTransform->TransformPoint(urp);
+ right_ll = rightToGroundTransform->TransformPoint(llp);
+ right_lr = rightToGroundTransform->TransformPoint(lrp);
+
double left_xmin = std::min(std::min(std::min(left_ul[0],left_ur[0]),left_lr[0]),left_ll[0]);
double left_xmax = std::max(std::max(std::max(left_ul[0],left_ur[0]),left_lr[0]),left_ll[0]);
double left_ymin = std::min(std::min(std::min(left_ul[1],left_ur[1]),left_lr[1]),left_ll[1]);
@@ -316,12 +331,6 @@ DisparityMapToDEMFilter<TDisparityImage,TInputImage,TOutputDEMImage,TEpipolarGri
itkExceptionMacro(<<"Wrong reconstruction area, images don't overlap, check image corners");
}
- // Choose origin
- typename TOutputDEMImage::PointType outOrigin;
- outOrigin[0] = box_xmin;
- outOrigin[1] = box_ymax;
- outputPtr->SetOrigin(outOrigin);
-
// Compute step :
// TODO : use a clean RS transform instead
typename TOutputDEMImage::SpacingType outSpacing;
@@ -329,6 +338,12 @@ DisparityMapToDEMFilter<TDisparityImage,TInputImage,TOutputDEMImage,TEpipolarGri
outSpacing[1] = -57.295779513 * m_DEMGridStep / 6378137.0;
outputPtr->SetSpacing(outSpacing);
+ // Choose origin
+ typename TOutputDEMImage::PointType outOrigin;
+ outOrigin[0] = box_xmin + 0.5 * outSpacing[0];
+ outOrigin[1] = box_ymax + 0.5 * outSpacing[1];
+ outputPtr->SetOrigin(outOrigin);
+
// Compute output size
typename DEMImageType::RegionType outRegion;
outRegion.SetIndex(0,0);
@@ -391,32 +406,32 @@ DisparityMapToDEMFilter<TDisparityImage,TInputImage,TOutputDEMImage,TEpipolarGri
// up left at elevation min
TDPointType corners[8];
- corners[0][0]= outOrigin[0] + outSpacing[0] * outRegion.GetIndex(0);
- corners[0][1]= outOrigin[1] + outSpacing[1] * outRegion.GetIndex(1);
+ corners[0][0]= outOrigin[0] + outSpacing[0] * (-0.5 + static_cast<double>(outRegion.GetIndex(0)));
+ corners[0][1]= outOrigin[1] + outSpacing[1] * (-0.5 + static_cast<double>(outRegion.GetIndex(1)));
corners[0][2]= m_ElevationMin;
// up left at elevation max
corners[1][0]= corners[0][0];
corners[1][1]= corners[0][1];
corners[1][2]= m_ElevationMax;
// up right at elevation min
- corners[2][0]= outOrigin[0] + outSpacing[0] * (outRegion.GetIndex(0) + outRegion.GetSize(0));
- corners[2][1]= outOrigin[1] + outSpacing[1] * outRegion.GetIndex(1);
+ corners[2][0]= corners[0][0] + outSpacing[0] * static_cast<double>(outRegion.GetSize(0));
+ corners[2][1]= corners[0][1];
corners[2][2]= m_ElevationMin;
// up right at elevation max
corners[3][0]= corners[2][0];
corners[3][1]= corners[2][1];
corners[3][2]= m_ElevationMax;
// low right at elevation min
- corners[4][0]= outOrigin[0] + outSpacing[0] * (outRegion.GetIndex(0) + outRegion.GetSize(0));
- corners[4][1]= outOrigin[1] + outSpacing[1] * (outRegion.GetIndex(1) + outRegion.GetSize(1));
+ corners[4][0]= corners[0][0] + outSpacing[0] * static_cast<double>(outRegion.GetSize(0));
+ corners[4][1]= corners[0][1] + outSpacing[1] * static_cast<double>(outRegion.GetSize(1));
corners[4][2]= m_ElevationMin;
// low right at elevation max
corners[5][0]= corners[4][0];
corners[5][1]= corners[4][1];
corners[5][2]= m_ElevationMax;
// low left at elevation min
- corners[6][0]= outOrigin[0] + outSpacing[0] * outRegion.GetIndex(0);
- corners[6][1]= outOrigin[1] + outSpacing[1] * (outRegion.GetIndex(1) + outRegion.GetSize(1));
+ corners[6][0]= corners[0][0];
+ corners[6][1]= corners[0][1] + outSpacing[1] * static_cast<double>(outRegion.GetSize(1));
corners[6][2]= m_ElevationMin;
// low left at elevation max
corners[7][0]= corners[6][0];
@@ -540,8 +555,8 @@ DisparityMapToDEMFilter<TDisparityImage,TInputImage,TOutputDEMImage,TEpipolarGri
}
typename DisparityMapType::RegionType inputDisparityRegion;
- inputDisparityRegion.SetIndex(0, static_cast<int>(vcl_floor(epiIndexMin[0])));
- inputDisparityRegion.SetIndex(1, static_cast<int>(vcl_floor(epiIndexMin[1])));
+ inputDisparityRegion.SetIndex(0, static_cast<int>(vcl_floor(epiIndexMin[0] + 0.5)));
+ inputDisparityRegion.SetIndex(1, static_cast<int>(vcl_floor(epiIndexMin[1] + 0.5)));
inputDisparityRegion.SetSize(0, 1 + static_cast<unsigned int>(vcl_floor(epiIndexMax[0] - epiIndexMin[0] + 0.5)));
inputDisparityRegion.SetSize(1, 1 + static_cast<unsigned int>(vcl_floor(epiIndexMax[1] - epiIndexMin[1] + 0.5)));
diff --git a/Code/DisparityMap/otbDisparityTranslateFilter.h b/Code/DisparityMap/otbDisparityTranslateFilter.h
index 6754a9fbd880205fff29b231699640a09d3e25e0..71fcfe7892e9da4c77fb14b8b88c3e16b901d0fd 100644
--- a/Code/DisparityMap/otbDisparityTranslateFilter.h
+++ b/Code/DisparityMap/otbDisparityTranslateFilter.h
@@ -69,6 +69,7 @@ public:
typedef typename DispMapType::IndexType IndexType;
typedef typename DispMapType::SpacingType SpacingType;
typedef typename DispMapType::PointType PointType;
+ typedef typename DispMapType::IndexValueType IndexValueType;
typedef typename GridType::RegionType GridRegionType;
diff --git a/Code/DisparityMap/otbDisparityTranslateFilter.txx b/Code/DisparityMap/otbDisparityTranslateFilter.txx
index 5a75e759807101675baea1782ebcc747f64806a1..c2078afe965607ef041a477ffb49c2bdfaf9b3d5 100644
--- a/Code/DisparityMap/otbDisparityTranslateFilter.txx
+++ b/Code/DisparityMap/otbDisparityTranslateFilter.txx
@@ -245,12 +245,12 @@ DisparityTranslateFilter<TDisparityImage,TGridImage,TSensorImage,TMaskImage>
IndexType corners[4];
corners[0] = requested.GetIndex();
corners[1] = requested.GetIndex();
- corners[1][0] += requested.GetSize()[0];
+ corners[1][0] += static_cast<IndexValueType>(requested.GetSize()[0]) - 1;
corners[2] = requested.GetIndex();
- corners[2][1] += requested.GetSize()[1];
+ corners[2][1] += static_cast<IndexValueType>(requested.GetSize()[1]) - 1;
corners[3] = requested.GetIndex();
- corners[3][0] += requested.GetSize()[0];
- corners[3][1] += requested.GetSize()[1];
+ corners[3][0] += static_cast<IndexValueType>(requested.GetSize()[0]) - 1;
+ corners[3][1] += static_cast<IndexValueType>(requested.GetSize()[1]) - 1;
for (unsigned int k=0; k<4; ++k)
{
PointType pointSensor;
diff --git a/Code/DisparityMap/otbMulti3DMapToDEMFilter.txx b/Code/DisparityMap/otbMulti3DMapToDEMFilter.txx
index 282521de3a5cd1ba0a2fdd221410742d91cc043b..e42ab2358d0d32255186919f06ee53e968c57494 100644
--- a/Code/DisparityMap/otbMulti3DMapToDEMFilter.txx
+++ b/Code/DisparityMap/otbMulti3DMapToDEMFilter.txx
@@ -193,21 +193,29 @@ void Multi3DMapToDEMFilter<T3DImage, TMaskImage, TOutputDEMImage>::SetOutputPara
typename InputMapType::SizeType inputSize = imgPtr->GetLargestPossibleRegion().GetSize();
- typename InputMapType::PointType tmpPoint;
- tmpPoint = imgPtr->GetOrigin();
- RSTransform2DType::OutputPointType ul = mapToGroundTransform->TransformPoint(tmpPoint);
-
- tmpPoint[0] = (imgPtr->GetOrigin())[0] + (imgPtr->GetSpacing())[0] * static_cast<double> (inputSize[0]);
- tmpPoint[1] = (imgPtr->GetOrigin())[1];
- RSTransform2DType::OutputPointType ur = mapToGroundTransform->TransformPoint(tmpPoint);
-
- tmpPoint[0] = (imgPtr->GetOrigin())[0] + (imgPtr->GetSpacing())[0] * static_cast<double> (inputSize[0]);
- tmpPoint[1] = (imgPtr->GetOrigin())[1] + (imgPtr->GetSpacing())[1] * static_cast<double> (inputSize[1]);
- RSTransform2DType::OutputPointType lr = mapToGroundTransform->TransformPoint(tmpPoint);
-
- tmpPoint[0] = (imgPtr->GetOrigin())[0];
- tmpPoint[1] = (imgPtr->GetOrigin())[1] + (imgPtr->GetSpacing())[1] * static_cast<double> (inputSize[1]);
- RSTransform2DType::OutputPointType ll = mapToGroundTransform->TransformPoint(tmpPoint);
+ typename InputMapType::PointType ulPt, urPt, llPt, lrPt;
+ itk::ContinuousIndex<double,2> ulIdx(imgPtr->GetLargestPossibleRegion().GetIndex());
+ ulIdx[0] += -0.5;
+ ulIdx[1] += -0.5;
+
+ itk::ContinuousIndex<double,2> urIdx(ulIdx);
+ itk::ContinuousIndex<double,2> lrIdx(ulIdx);
+ itk::ContinuousIndex<double,2> llIdx(ulIdx);
+ urIdx[0] += static_cast<double>(inputSize[0]);
+ lrIdx[0] += static_cast<double>(inputSize[0]);
+ lrIdx[1] += static_cast<double>(inputSize[1]);
+ llIdx[1] += static_cast<double>(inputSize[1]);
+
+ imgPtr->TransformContinuousIndexToPhysicalPoint(ulIdx,ulPt);
+ imgPtr->TransformContinuousIndexToPhysicalPoint(urIdx,urPt);
+ imgPtr->TransformContinuousIndexToPhysicalPoint(llIdx,llPt);
+ imgPtr->TransformContinuousIndexToPhysicalPoint(lrIdx,lrPt);
+
+ RSTransform2DType::OutputPointType ul, ur, lr, ll;
+ ul = mapToGroundTransform->TransformPoint(ulPt);
+ ur = mapToGroundTransform->TransformPoint(urPt);
+ ll = mapToGroundTransform->TransformPoint(llPt);
+ lr = mapToGroundTransform->TransformPoint(lrPt);
double xmin = std::min(std::min(std::min(ul[0], ur[0]), lr[0]), ll[0]);
double xmax = std::max(std::max(std::max(ul[0], ur[0]), lr[0]), ll[0]);
@@ -227,12 +235,6 @@ void Multi3DMapToDEMFilter<T3DImage, TMaskImage, TOutputDEMImage>::SetOutputPara
}
- // Choose origin
- typename TOutputDEMImage::PointType outOrigin;
- outOrigin[0] = box_xmin;
- outOrigin[1] = box_ymax;
- outputPtr->SetOrigin(outOrigin);
-
// Compute step :
// TODO : use a clean RS transform instead
typename TOutputDEMImage::SpacingType outSpacing;
@@ -240,6 +242,13 @@ void Multi3DMapToDEMFilter<T3DImage, TMaskImage, TOutputDEMImage>::SetOutputPara
outSpacing[0] = 57.295779513 * m_DEMGridStep / (6378137.0 * vcl_cos((box_ymin + box_ymax) * 0.5 * 0.01745329251));
outSpacing[1] = -57.295779513 * m_DEMGridStep / 6378137.0;
outputPtr->SetSpacing(outSpacing);
+
+ // Choose origin
+ typename TOutputDEMImage::PointType outOrigin;
+ outOrigin[0] = box_xmin + 0.5 * outSpacing[0];
+ outOrigin[1] = box_ymax + 0.5 * outSpacing[1];
+ outputPtr->SetOrigin(outOrigin);
+
// Compute output size
typename TOutputDEMImage::RegionType outRegion;
outRegion.SetIndex(0, 0);
@@ -346,32 +355,32 @@ void Multi3DMapToDEMFilter<T3DImage, TMaskImage, TOutputDEMImage>::GenerateInput
// up left at elevation min
TDPointType corners[8];
- corners[0][0] = outOrigin[0] + outSpacing[0] * outRegion.GetIndex(0);
- corners[0][1] = outOrigin[1] + outSpacing[1] * outRegion.GetIndex(1);
+ corners[0][0] = outOrigin[0] + outSpacing[0] * (-0.5 + static_cast<double>(outRegion.GetIndex(0)));
+ corners[0][1] = outOrigin[1] + outSpacing[1] * (-0.5 + static_cast<double>(outRegion.GetIndex(1)));
corners[0][2] = m_ElevationMin;
// up left at elevation max
corners[1][0] = corners[0][0];
corners[1][1] = corners[0][1];
corners[1][2] = m_ElevationMax;
// up right at elevation min
- corners[2][0] = outOrigin[0] + outSpacing[0] * (outRegion.GetIndex(0) + outRegion.GetSize(0));
- corners[2][1] = outOrigin[1] + outSpacing[1] * outRegion.GetIndex(1);
+ corners[2][0] = corners[0][0] + outSpacing[0] * static_cast<double>(outRegion.GetSize(0));
+ corners[2][1] = corners[0][1];
corners[2][2] = m_ElevationMin;
// up right at elevation max
corners[3][0] = corners[2][0];
corners[3][1] = corners[2][1];
corners[3][2] = m_ElevationMax;
// low right at elevation min
- corners[4][0] = outOrigin[0] + outSpacing[0] * (outRegion.GetIndex(0) + outRegion.GetSize(0));
- corners[4][1] = outOrigin[1] + outSpacing[1] * (outRegion.GetIndex(1) + outRegion.GetSize(1));
+ corners[4][0] = corners[0][0] + outSpacing[0] * static_cast<double>(outRegion.GetSize(0));
+ corners[4][1] = corners[0][1] + outSpacing[1] * static_cast<double>(outRegion.GetSize(1));
corners[4][2] = m_ElevationMin;
// low right at elevation max
corners[5][0] = corners[4][0];
corners[5][1] = corners[4][1];
corners[5][2] = m_ElevationMax;
// low left at elevation min
- corners[6][0] = outOrigin[0] + outSpacing[0] * outRegion.GetIndex(0);
- corners[6][1] = outOrigin[1] + outSpacing[1] * (outRegion.GetIndex(1) + outRegion.GetSize(1));
+ corners[6][0] = corners[0][0];
+ corners[6][1] = corners[0][1] + outSpacing[1] * static_cast<double>(outRegion.GetSize(1));
corners[6][2] = m_ElevationMin;
// low left at elevation max
corners[7][0] = corners[6][0];
@@ -633,10 +642,11 @@ void Multi3DMapToDEMFilter<T3DImage, TMaskImage, TOutputDEMImage>::ThreadedGener
itk::ContinuousIndex<double, 2> continuousIndex;
//std::cout<<"point2D "<<point2D<<std::endl;
+ // The DEM cell at index 'n' contains continuous indexes from 'n-0.5' to 'n+0.5'
outputPtr->TransformPhysicalPointToContinuousIndex(point2D, continuousIndex);
typename OutputImageType::IndexType cellIndex;
- cellIndex[0] = static_cast<int> (vcl_floor(continuousIndex[0]));
- cellIndex[1] = static_cast<int> (vcl_floor(continuousIndex[1]));
+ cellIndex[0] = static_cast<int> (vcl_floor(continuousIndex[0] + 0.5));
+ cellIndex[1] = static_cast<int> (vcl_floor(continuousIndex[1] + 0.5));
//std::cout<<"cellindex "<<cellIndex<<std::endl;
//index from physical
/** -Wunused-variable