From b451247f3242d88f18bed827ad538d5617d937bd Mon Sep 17 00:00:00 2001
From: Guillaume Pasero <guillaume.pasero@c-s.fr>
Date: Fri, 6 Jun 2014 18:20:39 +0200
Subject: [PATCH] BUG: Mantis-942: pixel convention consistency for
 Code/Projections filters

---
 .../otbImageToGenericRSOutputParameters.txx   | 42 +++++++++----------
 ...otbVectorDataIntoImageProjectionFilter.txx | 31 ++++++++------
 2 files changed, 39 insertions(+), 34 deletions(-)

diff --git a/Code/Projections/otbImageToGenericRSOutputParameters.txx b/Code/Projections/otbImageToGenericRSOutputParameters.txx
index 3f266c7197..54a8997afa 100644
--- a/Code/Projections/otbImageToGenericRSOutputParameters.txx
+++ b/Code/Projections/otbImageToGenericRSOutputParameters.txx
@@ -21,6 +21,7 @@
 
 #include "otbImageToGenericRSOutputParameters.h"
 #include "itkMacro.h"
+#include "itkContinuousIndex.h"
 
 namespace otb {
 
@@ -56,15 +57,15 @@ ImageToGenericRSOutputParameters<TImage>
   // Estimate the Output image Extent
   this->EstimateOutputImageExtent();
 
-  // Estimate the Output Origin
-  this->EstimateOutputOrigin();
-
   // Estimate the Output Spacing
   if(!m_ForceSpacing)
     this->EstimateOutputSpacing();
 
   // Finally Estimate the Output Size
   this->EstimateOutputSize();
+
+  // Estimate the Output Origin
+  this->EstimateOutputOrigin();
 }
 
 
@@ -94,25 +95,24 @@ ImageToGenericRSOutputParameters<TImage>
   m_Transform->GetInverse(invTransform);
 
   // Compute the 4 corners in the cartographic coordinate system
-  std::vector<IndexType>       vindex;
+  std::vector< itk::ContinuousIndex<double,2> > vindex;
   std::vector<PointType> voutput;
 
-  IndexType index1, index2, index3, index4;
-  SizeType  size;
+  itk::ContinuousIndex<double,2> index1(m_Input->GetLargestPossibleRegion().GetIndex());
+  index1[0] += -0.5;
+  index1[1] += -0.5;
+  itk::ContinuousIndex<double,2> index2(index1);
+  itk::ContinuousIndex<double,2> index3(index1);
+  itk::ContinuousIndex<double,2> index4(index1);
 
   // Image size
-  size = m_Input->GetLargestPossibleRegion().GetSize();
+  SizeType size = m_Input->GetLargestPossibleRegion().GetSize();
 
   // project the 4 corners
-  index1 = m_Input->GetLargestPossibleRegion().GetIndex();
-  index2 = m_Input->GetLargestPossibleRegion().GetIndex();
-  index3 = m_Input->GetLargestPossibleRegion().GetIndex();
-  index4 = m_Input->GetLargestPossibleRegion().GetIndex();
-
-  index2[0] += size[0] - 1;
-  index3[0] += size[0] - 1;
-  index3[1] += size[1] - 1;
-  index4[1] += size[1] - 1;
+  index2[0] += size[0];
+  index3[0] += size[0];
+  index3[1] += size[1];
+  index4[1] += size[1];
 
   vindex.push_back(index1);
   vindex.push_back(index2);
@@ -122,7 +122,7 @@ ImageToGenericRSOutputParameters<TImage>
   for (unsigned int i = 0; i < vindex.size(); ++i)
     {
     PointType physicalPoint;
-    m_Input->TransformIndexToPhysicalPoint(vindex[i], physicalPoint);
+    m_Input->TransformContinuousIndexToPhysicalPoint(vindex[i], physicalPoint);
     voutput.push_back(invTransform->TransformPoint(physicalPoint));
     }
 
@@ -168,8 +168,8 @@ ImageToGenericRSOutputParameters<TImage>
   // Set the output orgin in carto
   // projection
   PointType   origin;
-  origin[0] = m_OutputExtent.minX;
-  origin[1] = m_OutputExtent.maxY;
+  origin[0] = m_OutputExtent.minX + 0.5 * this->GetOutputSpacing()[0];
+  origin[1] = m_OutputExtent.maxY + 0.5 * this->GetOutputSpacing()[1];
   this->SetOutputOrigin(origin);
 }
 
@@ -187,8 +187,8 @@ ImageToGenericRSOutputParameters<TImage>
   double sizeCartoY = vcl_abs(m_OutputExtent.minY - m_OutputExtent.maxY);
 
   PointType o, oX, oY;
-  o[0] = this->GetOutputOrigin()[0];
-  o[1] = this->GetOutputOrigin()[1];
+  o[0] = m_OutputExtent.minX;
+  o[1] = m_OutputExtent.maxY;
 
   oX = o;
   oY = o;
diff --git a/Code/Projections/otbVectorDataIntoImageProjectionFilter.txx b/Code/Projections/otbVectorDataIntoImageProjectionFilter.txx
index 636e6e67cb..cc0c0eef0a 100644
--- a/Code/Projections/otbVectorDataIntoImageProjectionFilter.txx
+++ b/Code/Projections/otbVectorDataIntoImageProjectionFilter.txx
@@ -119,6 +119,7 @@ VectorDataIntoImageProjectionFilter<TInputVectorData, TInputImage>
 
   typedef typename ImageType::IndexType       IndexType;
   typedef typename ImageType::PointType       PointType;
+  typedef typename ImageType::SizeType       SizeType;
 
   if (m_InputImage.IsNull())
     {
@@ -132,24 +133,28 @@ VectorDataIntoImageProjectionFilter<TInputVectorData, TInputImage>
   std::cout << "ProjRef of the input vector data: "<< this->GetInput()->GetProjectionRef() << std::endl; */
 
   // Get the index of the corner of the image
-  IndexType ul, ur, ll, lr;
   PointType pul, pur, pll, plr;
-  ul = m_InputImage->GetLargestPossibleRegion().GetIndex();
-  ur = ul;
-  ll = ul;
-  lr = ul;
-  ur[0] += m_InputImage->GetLargestPossibleRegion().GetSize()[0];
-  lr[0] += m_InputImage->GetLargestPossibleRegion().GetSize()[0];
-  lr[1] += m_InputImage->GetLargestPossibleRegion().GetSize()[1];
-  ll[1] += m_InputImage->GetLargestPossibleRegion().GetSize()[1];
+  itk::ContinuousIndex<double,2> ul(m_InputImage->GetLargestPossibleRegion().GetIndex());
+  ul[0] += -0.5;
+  ul[1] += -0.5;
+
+  itk::ContinuousIndex<double,2> ur(ul);
+  itk::ContinuousIndex<double,2> ll(ul);
+  itk::ContinuousIndex<double,2> lr(ul);
+
+  SizeType size = m_InputImage->GetLargestPossibleRegion().GetSize();
+  ur[0] += size[0];
+  lr[0] += size[0];
+  lr[1] += size[1];
+  ll[1] += size[1];
 
   //std::cout << "bounding box of the input image (pixel): "<< ur << ", " << ul << ", " << lr << ", " << ll << std::endl;
 
   // Transform to physical point
-  m_InputImage->TransformIndexToPhysicalPoint(ul, pul);
-  m_InputImage->TransformIndexToPhysicalPoint(ur, pur);
-  m_InputImage->TransformIndexToPhysicalPoint(ll, pll);
-  m_InputImage->TransformIndexToPhysicalPoint(lr, plr);
+  m_InputImage->TransformContinuousIndexToPhysicalPoint(ul, pul);
+  m_InputImage->TransformContinuousIndexToPhysicalPoint(ur, pur);
+  m_InputImage->TransformContinuousIndexToPhysicalPoint(ll, pll);
+  m_InputImage->TransformContinuousIndexToPhysicalPoint(lr, plr);
   //std::cout << "bounding box of the input image (physical): "<< pur << ", " << pul << ", " << plr << ", " << pll << std::endl;
 
   // Build the cartographic region
-- 
GitLab