From 364674a4ea7674190a83f8b72cab54543d9d265c Mon Sep 17 00:00:00 2001 From: Guillaume Pasero <guillaume.pasero@c-s.fr> Date: Tue, 2 Sep 2014 19:35:29 +0200 Subject: [PATCH] BUG: Mantis-942 : fix output physical space --- .../Projections/otbRigidTransformResample.cxx | 85 +++++++++++-------- 1 file changed, 48 insertions(+), 37 deletions(-) diff --git a/Applications/Projections/otbRigidTransformResample.cxx b/Applications/Projections/otbRigidTransformResample.cxx index 59fc3029d0..91532f1164 100644 --- a/Applications/Projections/otbRigidTransformResample.cxx +++ b/Applications/Projections/otbRigidTransformResample.cxx @@ -230,12 +230,18 @@ private: // Evaluate spacing FloatVectorImageType::SpacingType spacing = inputImage->GetSpacing(); FloatVectorImageType::SpacingType OutputSpacing; - OutputSpacing=spacing; - OutputSpacing[0] = spacing[0] * scale[0]; OutputSpacing[1] = spacing[1] * scale[1]; m_Resampler->SetOutputSpacing(OutputSpacing); + + FloatVectorImageType::PointType origin = inputImage->GetOrigin(); + FloatVectorImageType::PointType outputOrigin; + outputOrigin[0] = origin[0] + 0.5 * spacing[0] * (scale[0] - 1.0); + outputOrigin[1] = origin[1] + 0.5 * spacing[1] * (scale[1] - 1.0); + + m_Resampler->SetOutputOrigin(outputOrigin); + m_Resampler->SetTransform(transform); // Evaluate size @@ -267,13 +273,18 @@ private: // Evaluate spacing FloatVectorImageType::SpacingType spacing = inputImage->GetSpacing(); FloatVectorImageType::SpacingType OutputSpacing; - OutputSpacing=spacing; - OutputSpacing[0] = spacing[0] * scale[0]; OutputSpacing[1] = spacing[1] * scale[1]; m_Resampler->SetOutputSpacing(OutputSpacing); + FloatVectorImageType::PointType origin = inputImage->GetOrigin(); + FloatVectorImageType::PointType outputOrigin; + outputOrigin[0] = origin[0] + 0.5 * spacing[0] * (scale[0] - 1.0); + outputOrigin[1] = origin[1] + 0.5 * spacing[1] * (scale[1] - 1.0); + + m_Resampler->SetOutputOrigin(outputOrigin); + ResampleFilterType::SizeType recomputedSize; recomputedSize[0] = inputImage->GetLargestPossibleRegion().GetSize()[0] / scale[0]; recomputedSize[1] = inputImage->GetLargestPossibleRegion().GetSize()[1] / scale[1]; @@ -290,34 +301,33 @@ private: { ScalableTransformType::Pointer transform = ScalableTransformType::New(); - FloatVectorImageType::IndexType origin = inputImage->GetLargestPossibleRegion().GetIndex(); + FloatVectorImageType::SizeType inSize = inputImage->GetLargestPossibleRegion().GetSize(); FloatVectorImageType::SpacingType spacing = inputImage->GetSpacing(); - FloatVectorImageType::IndexType center; - center[0] = origin[0] + inputImage->GetLargestPossibleRegion().GetSize()[0] / 2.0; - center[1] = origin[1] + inputImage->GetLargestPossibleRegion().GetSize()[1] / 2.0; - + itk::ContinuousIndex<double,2> ULindex(inputImage->GetLargestPossibleRegion().GetIndex()); + ULindex[0] += -0.5; + ULindex[1] += -0.5; + itk::ContinuousIndex<double,2> center, URindex, LRindex, LLindex; + center[0] = ULindex[0] + static_cast<double>(inSize[0]) / 2.0; + center[1] = ULindex[1] + static_cast<double>(inSize[1]) / 2.0; + + URindex = ULindex; + LRindex = ULindex; + LLindex = ULindex; + URindex[0] += inSize[0]; + LRindex[0] += inSize[0]; + LRindex[1] += inSize[1]; + LLindex[1] += inSize[1]; FloatVectorImageType::PointType centerPoint; - inputImage->TransformIndexToPhysicalPoint(center, centerPoint); + inputImage->TransformContinuousIndexToPhysicalPoint(center, centerPoint); //image boundary - FloatVectorImageType::IndexType ULindex, URindex, LRindex, LLindex; - - ULindex[0]=origin[0]; - ULindex[1]=origin[1]; - URindex[0]=origin[0]+ inputImage->GetLargestPossibleRegion().GetSize()[0]; - URindex[1]=origin[1]; - LRindex[0]=origin[0]+ inputImage->GetLargestPossibleRegion().GetSize()[0]; - LRindex[1]=origin[1]+ inputImage->GetLargestPossibleRegion().GetSize()[1]; - LLindex[0]=origin[0]; - LLindex[1]=origin[1]+ inputImage->GetLargestPossibleRegion().GetSize()[1]; - - FloatVectorImageType::PointType orig, ULpoint, URpoint, LRpoint, LLpoint; - inputImage->TransformIndexToPhysicalPoint(ULindex, ULpoint); - inputImage->TransformIndexToPhysicalPoint(URindex, URpoint); - inputImage->TransformIndexToPhysicalPoint(LRindex, LRpoint); - inputImage->TransformIndexToPhysicalPoint(LLindex, LLpoint); + FloatVectorImageType::PointType ULpoint, URpoint, LRpoint, LLpoint; + inputImage->TransformContinuousIndexToPhysicalPoint(ULindex, ULpoint); + inputImage->TransformContinuousIndexToPhysicalPoint(URindex, URpoint); + inputImage->TransformContinuousIndexToPhysicalPoint(LRindex, LRpoint); + inputImage->TransformContinuousIndexToPhysicalPoint(LLindex, LLpoint); // Scale Transform OutputVectorType scale; @@ -344,13 +354,11 @@ private: m_Resampler->SetTransform(transform); - FloatVectorImageType::PointType ULpointTrans, URpointTrans, LRpointTrans, LLpointTrans, CenterPointTrans; - + FloatVectorImageType::PointType ULpointTrans, URpointTrans, LRpointTrans, LLpointTrans; ULpointTrans=inverseTransform->TransformPoint(ULpoint); URpointTrans=inverseTransform->TransformPoint(URpoint); LRpointTrans=inverseTransform->TransformPoint(LRpoint); LLpointTrans=inverseTransform->TransformPoint(LLpoint); - CenterPointTrans=inverseTransform->TransformPoint(centerPoint); //compute min and max std::vector<FloatVectorImageType::PointType> voutput; @@ -364,7 +372,7 @@ private: double minY = voutput[0][1]; double maxY = voutput[0][1]; - for(unsigned int i = 0; i<voutput.size(); i++) + for(unsigned int i = 1; i<voutput.size(); i++) { // Origins if ( minX > voutput[i][0] ) @@ -387,16 +395,15 @@ private: } } - if( spacing[0] > 0 ) orig[0] = minX; - else orig[0] = maxX; + FloatVectorImageType::PointType outputOrig; + if( spacing[0] > 0 ) outputOrig[0] = minX; + else outputOrig[0] = maxX; - if( spacing[1] > 0 ) orig[1] = minY; - else orig[1] = maxY; - - m_Resampler->SetOutputOrigin(orig); + if( spacing[1] > 0 ) outputOrig[1] = minY; + else outputOrig[1] = maxY; //size of output image - ResampleFilterType::SizeType size; + FloatVectorImageType::PointType size; size[0]=vcl_abs(maxX-minX); size[1]=vcl_abs(maxY-minY); @@ -407,6 +414,10 @@ private: m_Resampler->SetOutputSpacing(OutputSpacing); + outputOrig[0] += 0.5 * OutputSpacing[0]; + outputOrig[1] += 0.5 * OutputSpacing[1]; + m_Resampler->SetOutputOrigin(outputOrig); + // Evaluate size ResampleFilterType::SizeType recomputedSize; recomputedSize[0] = static_cast<unsigned int>(vcl_floor(vcl_abs(size[0]/OutputSpacing[0]))); -- GitLab