diff --git a/Applications/Projections/otbRigidTransformResample.cxx b/Applications/Projections/otbRigidTransformResample.cxx
index 7b869a8ffaca6f606f4e78c245026146629b9015..9629c8790470862fd60e600c2234b1d5e4ac492e 100644
--- a/Applications/Projections/otbRigidTransformResample.cxx
+++ b/Applications/Projections/otbRigidTransformResample.cxx
@@ -18,11 +18,35 @@
#include "otbWrapperApplication.h"
#include "otbWrapperApplicationFactory.h"
+#include "itkLinearInterpolateImageFunction.h"
+#include "otbBCOInterpolateImageFunction.h"
+#include "itkNearestNeighborInterpolateImageFunction.h"
+
+//Transform
+#include "otbCompositeTransform.h"
+#include "itkScalableAffineTransform.h"
#include "itkTranslationTransform.h"
+#include "itkIdentityTransform.h"
+#include "itkScaleTransform.h"
+#include "itkCenteredRigid2DTransform.h"
+
#include "otbStreamingResampleImageFilter.h"
namespace otb
{
+enum
+{
+ Interpolator_NNeighbor,
+ Interpolator_Linear,
+ Interpolator_BCO
+};
+
+enum
+{
+ Transform_Translation,
+ Transform_Rotation,
+};
+
namespace Wrapper
{
@@ -41,6 +65,15 @@ public:
typedef itk::TranslationTransform<double, FloatVectorImageType::ImageDimension> TransformType;
typedef otb::StreamingResampleImageFilter<FloatVectorImageType, FloatVectorImageType, double> ResampleFilterType;
+ //typedef itk::IdentityTransform<double, FloatVectorImageType::ImageDimension> IdentityTransformType;
+
+ typedef itk::ScalableAffineTransform<double, FloatVectorImageType::ImageDimension> ScalableTransformType;
+ typedef ScalableTransformType::OutputVectorType OutputVectorType;
+
+ /** Rotation transform */
+ typedef itk::CenteredRigid2DTransform< double > RotationTransformType;
+ typedef RotationTransformType::ScalarType ScalarType;
+
itkTypeMacro(RigidTransformResample, otb::Application);
private:
@@ -55,7 +88,7 @@ private:
SetDocLimitations("None");
SetDocAuthors("OTB-Team");
SetDocSeeAlso("Translation");
-
+
AddDocTag("Conversion");
AddDocTag(Tags::Geometry);
@@ -63,18 +96,64 @@ private:
SetParameterDescription("in","The input image to translate.");
AddParameter(ParameterType_OutputImage, "out", "Output image");
SetParameterDescription("out","The translated output image.");
- AddParameter(ParameterType_Float, "tx", "The X translation (in physical units)");
- SetParameterDescription("tx","The translation value along X axis (in physical units).");
- AddParameter(ParameterType_Float, "ty", "The Y translation (in physical units)");
- SetParameterDescription("ty","The translation value along Y axis (in physical units)");
+
+ //Transform
+ AddParameter(ParameterType_Group,"transform","Transform parameters");
+ SetParameterDescription("transform","This group of parameters allows to set the transformation to apply.");
+
+ AddParameter(ParameterType_Choice, "transform.type", "Type of transformation");
+ SetParameterDescription("transform.type","Type of transformation. Available transformations are translation and rotation with scaling factor");
+
+ AddChoice("transform.type.translation", "translation");
+ SetParameterDescription("transform.type.translation","translation");
+
+ AddParameter(ParameterType_Float,"transform.type.translation.tx", "The X translation (in physical units)");
+ SetParameterDescription("transform.type.translation.tx","The translation value along X axis (in physical units).");
+ SetDefaultParameterFloat("transform.type.translation.tx",0.);
+ AddParameter(ParameterType_Float,"transform.type.translation.ty", "The Y translation (in physical units)");
+ SetParameterDescription("transform.type.translation.ty","The translation value along Y axis (in physical units)");
+ SetDefaultParameterFloat("transform.type.translation.tx",0.);
+
+ AddChoice("transform.type.rotation", "rotation");
+ SetParameterDescription("transform.type.rotation","rotation");
+
+ AddParameter(ParameterType_Float, "transform.type.rotation.angle", "Rotation angle");
+ SetParameterDescription("transform.type.rotation.angle","The rotation angle in degree (values between -180 and 180)");
+ SetDefaultParameterFloat("transform.type.rotation.angle",0.);
+
+ AddParameter(ParameterType_Float, "transform.type.rotation.scalex", "X factor");
+ SetParameterDescription("transform.type.rotation.scalex","X factor");
+ SetDefaultParameterFloat("transform.type.rotation.scalex",1.);
+
+ AddParameter(ParameterType_Float, "transform.type.rotation.scaley", "Y factor");
+ SetParameterDescription("transform.type.rotation.scaley","Y factor");
+ SetDefaultParameterFloat("transform.type.rotation.scaley",1.);
+
+ // Interpolators
+ AddParameter(ParameterType_Choice, "interpolator", "Interpolation");
+ SetParameterDescription("interpolator","This group of parameters allows to define how the input image will be interpolated during resampling.");
+ AddChoice("interpolator.nn", "Nearest Neighbor interpolation");
+ SetParameterDescription("interpolator.nn","Nearest neighbor interpolation leads to poor image quality, but it is very fast.");
+ AddChoice("interpolator.linear", "Linear interpolation");
+ SetParameterDescription("interpolator.linear","Linear interpolation leads to average image quality but is quite fast");
+ AddChoice("interpolator.bco", "Bicubic interpolation");
+ AddParameter(ParameterType_Radius, "interpolator.bco.radius", "Radius for bicubic interpolation");
+ SetParameterDescription("interpolator.bco.radius","This parameter allows to control the size of the bicubic interpolation filter. If the target pixel size is higher than the input pixel size, increasing this parameter will reduce aliasing artefacts.");
+ SetDefaultParameterInt("interpolator.bco.radius", 2);
+ SetParameterString("interpolator","bco");
+
+ // RAM available
+ AddRAMParameter("ram");
+ SetParameterDescription("ram","This allows to set the maximum amount of RAM available for processing. As the writing task is time consuming, it is better to write large pieces of data, which can be achieved by increasing this parameter (pay attention to your system capabilities)");
// Doc example parameter settings
SetDocExampleParameterValue("in", "qb_toulouse_sub.tif");
SetDocExampleParameterValue("out", "rigitTransformImage.tif");
- SetDocExampleParameterValue("tx", "5");
- SetDocExampleParameterValue("ty", "5");
+ SetDocExampleParameterValue("transform.type", "rotation");
+ SetDocExampleParameterValue("transform.type.scalex", "4.");
+ SetDocExampleParameterValue("transform.type.scaley", "4.");
}
-
+
void DoUpdateParameters()
{
// Nothing to do here : all parameters are independent
@@ -82,32 +161,201 @@ private:
void DoExecute()
{
-
FloatVectorImageType* inputImage = GetParameterImage("in");
-
- m_Transform = TransformType::New();
- TransformType::OutputVectorType offset;
- offset[0] = GetParameterFloat("tx");
- offset[1] = GetParameterFloat("ty");
- otbAppLogINFO( << "Offset : " << offset );
- m_Transform->SetOffset(offset);
m_Resampler = ResampleFilterType::New();
m_Resampler->SetInput(inputImage);
- m_Resampler->SetTransform(m_Transform);
- m_Resampler->SetOutputParametersFromImage(inputImage);
+
+ // Get Interpolator
+ switch ( GetParameterInt("interpolator") )
+ {
+ case Interpolator_Linear:
+ {
+ typedef itk::LinearInterpolateImageFunction<FloatVectorImageType,
+ double> LinearInterpolationType;
+ LinearInterpolationType::Pointer interpolator = LinearInterpolationType::New();
+ m_Resampler->SetInterpolator(interpolator);
+ }
+ break;
+ case Interpolator_NNeighbor:
+ {
+ typedef itk::NearestNeighborInterpolateImageFunction<FloatVectorImageType,
+ double> NearestNeighborInterpolationType;
+ NearestNeighborInterpolationType::Pointer interpolator = NearestNeighborInterpolationType::New();
+ m_Resampler->SetInterpolator(interpolator);
+ }
+ break;
+ case Interpolator_BCO:
+ {
+ typedef otb::BCOInterpolateImageFunction<FloatVectorImageType> BCOInterpolationType;
+ BCOInterpolationType::Pointer interpolator = BCOInterpolationType::New();
+ interpolator->SetRadius(GetParameterInt("interpolator.bco.radius"));
+ m_Resampler->SetInterpolator(interpolator);
+ }
+ break;
+ }
+
+ // Get Transform
+ switch ( GetParameterInt("transform.type") )
+ {
+ case Transform_Translation:
+ {
+ TransformType::Pointer transform = TransformType::New();
+ TransformType::OutputVectorType offset;
+ offset[0] = GetParameterFloat("transform.type.translation.tx");
+ offset[1] = GetParameterFloat("transform.type.translation.ty");
+ otbAppLogINFO( << "Offset : " << offset );
+ transform->SetOffset(offset);
+ m_Resampler->SetTransform(transform);
+ m_Resampler->SetOutputParametersFromImage(inputImage);
+ }
+ break;
+
+ case Transform_Rotation:
+ {
+ ScalableTransformType::Pointer transform = ScalableTransformType::New();
+
+ FloatVectorImageType::IndexType origin = inputImage->GetLargestPossibleRegion().GetIndex();
+ 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;
+
+
+ FloatVectorImageType::PointType centerPoint;
+ inputImage->TransformIndexToPhysicalPoint(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);
+
+ // Scale Transform
+ OutputVectorType scale;
+ scale[0] = 1.0 / GetParameterFloat("transform.type.rotation.scalex");
+ scale[1] = 1.0 / GetParameterFloat("transform.type.rotation.scaley");
+
+ //angle of rotation
+ ScalarType rot_angle = GetParameterFloat("transform.type.rotation.angle");
+ if (rot_angle < -180 || rot_angle > 180)
+ {
+ itkExceptionMacro(<<"The rotation angle must value be between -180 and 180.");
+ }
+
+ transform->SetIdentity();
+ if(spacing[0] > 0 && spacing[1] > 0) transform->Rotate2D( rot_angle * CONST_PI_180 );
+ else transform->Rotate2D( - rot_angle * CONST_PI_180 );
+
+ transform->SetCenter( centerPoint );
+ transform->Scale( scale );
+
+ //inverse transform
+ ScalableTransformType::Pointer inverseTransform = ScalableTransformType::New();
+ transform->GetInverse(inverseTransform);
+ m_Resampler->SetTransform(transform);
+
+
+ FloatVectorImageType::PointType ULpointTrans, URpointTrans, LRpointTrans, LLpointTrans, CenterPointTrans;
+
+ 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;
+ voutput.push_back(ULpointTrans);
+ voutput.push_back(URpointTrans);
+ voutput.push_back(LRpointTrans);
+ voutput.push_back(LLpointTrans);
+
+ double minX = voutput[0][0];
+ double maxX = voutput[0][0];
+ double minY = voutput[0][1];
+ double maxY = voutput[0][1];
+
+ for(unsigned int i = 0; i<voutput.size(); i++)
+ {
+ // Origins
+ if ( minX > voutput[i][0] )
+ {
+ minX = voutput[i][0];
+ }
+ if ( minY > voutput[i][1] )
+ {
+ minY = voutput[i][1];
+ }
+
+ // Sizes
+ if ( maxX < voutput[i][0] )
+ {
+ maxX = voutput[i][0];
+ }
+ if ( maxY < voutput[i][1] )
+ {
+ maxY = voutput[i][1];
+ }
+ }
+
+ if( spacing[0] > 0 ) orig[0] = minX;
+ else orig[0] = maxX;
+
+ if( spacing[1] > 0 ) orig[1] = minY;
+ else orig[1] = maxY;
+
+ m_Resampler->SetOutputOrigin(orig);
+
+ //size of output image
+ ResampleFilterType::SizeType size;
+ size[0]=vcl_abs(maxX-minX);
+ size[1]=vcl_abs(maxY-minY);
+
+ // Evaluate spacing
+ FloatVectorImageType::SpacingType OutputSpacing;
+ OutputSpacing=spacing;
+
+ //FIXME Find a way to update the output image spacing after resampling
+ //OutputSpacing[0] = scale[0] * spacing[0];
+ //OutputSpacing[1] = scale[1] * spacing[1];
+
+ m_Resampler->SetOutputSpacing(OutputSpacing);
+
+ // Evaluate size
+ ResampleFilterType::SizeType recomputedSize;
+ recomputedSize[0] = static_cast<unsigned int>(vcl_floor(vcl_abs(size[0]/OutputSpacing[0])));
+ recomputedSize[1] = static_cast<unsigned int>(vcl_floor(vcl_abs(size[1]/OutputSpacing[1])));
+ m_Resampler->SetOutputSize( recomputedSize );
+ otbAppLogINFO( << "Output image size : " << recomputedSize );
+ }
+ break;
+ }
+
FloatVectorImageType::PixelType defaultValue;
itk::PixelBuilder<FloatVectorImageType::PixelType>::Zero(defaultValue,
inputImage->GetNumberOfComponentsPerPixel());
m_Resampler->SetEdgePaddingValue(defaultValue);
+
m_Resampler->UpdateOutputInformation();
-
// Output Image
SetParameterOutputImage("out", m_Resampler->GetOutput());
}
ResampleFilterType::Pointer m_Resampler;
- TransformType::Pointer m_Transform;
}; //class
diff --git a/Testing/Applications/Projections/CMakeLists.txt b/Testing/Applications/Projections/CMakeLists.txt
index 467ffe590dab63d1cb7980aa8ea8a31cf68d3c28..b4b01e9eb892d9d1ffe66941ffbb3c28b3a874b2 100644
--- a/Testing/Applications/Projections/CMakeLists.txt
+++ b/Testing/Applications/Projections/CMakeLists.txt
@@ -19,14 +19,16 @@ OTB_TEST_APPLICATION(NAME apTvPrOrthorectification_UTM
${BASELINE}/owTvOrthorectifTest_UTM.tif
${TEMP}/apTvPrOrthorectifTest_UTM.tif)
-ENDIF(OTB_DATA_USE_LARGEINPUT)
+ENDIF(OTB_DATA_USE_LARGEINPUT)
OTB_TEST_APPLICATION(NAME apTvPrRigidTransformResample
APP RigidTransformResample
OPTIONS -in ${INPUTDATA}/poupees.tif
-out ${TEMP}/apTvPrRigidTransformResampleTest.tif
- -tx 5
- -ty -5
+ -transform.type translation
+ -transform.type.translation.tx 5
+ -transform.type.translation.ty -5
+ -interpolator nn
VALID --compare-image ${NOTOL}
${BASELINE}/owTvRigidTransformResampleTest.tif
${TEMP}/apTvPrRigidTransformResampleTest.tif)
@@ -79,48 +81,48 @@ OTB_TEST_APPLICATION(NAME apTvPrBundleToPerfectSensor
VALID --compare-image ${EPSILON_7}
${BASELINE}/apTvPrBundleToPerfectSensor.tif
${TEMP}/apTvPrBundleToPerfectSensor.tif)
-
+
OTB_TEST_APPLICATION(NAME apTvPrVectorDataReprojectionFromMapToImage
APP VectorDataReprojection
OPTIONS -in.vd ${INPUTDATA}/ToulousePoints-examples.shp
-out.vd ${TEMP}/apTvPrVectorDataReprojectionFromMapToImage.shp
- -out.proj image
+ -out.proj image
-out.proj.image.in ${INPUTDATA}/QB_Toulouse_Ortho_XS.tif
VALID --compare-ogr ${NOTOL}
${BASELINE_FILES}/prTvVectorDataProjectionFilterFromMapToImage.shp
${TEMP}/apTvPrVectorDataReprojectionFromMapToImage.shp)
-# With QGIS the two vectordata are equal
+# With QGIS the two vectordata are equal
OTB_TEST_APPLICATION(NAME apTuPrVectorDataReprojectionFromImageToMap
APP VectorDataReprojection
OPTIONS -in.vd ${INPUTDATA}/ToulousePointsUTM31N.shp
-in.kwl ${INPUTDATA}/QB_Toulouse_Ortho_XS.tif
-out.vd ${TEMP}/apTvPrVectorDataReprojectionFromImageToMap.shp
- -out.proj user
+ -out.proj user
-out.proj.user.map utm
-out.proj.user.map.utm.northhem true
-out.proj.user.map.utm.zone 31)
# VALID --compare-ogr ${NOTOL}
# ${INPUTDATA}/ToulousePoints-examples.shp
-# ${TEMP}/apTvPrVectorDataReprojectionFromImageToMap.shp)
-
-
+# ${TEMP}/apTvPrVectorDataReprojectionFromImageToMap.shp)
+
+
OTB_TEST_APPLICATION(NAME apTvPrVectorDataReprojectionFromMapToGeo
APP VectorDataReprojection
OPTIONS -in.vd ${INPUTDATA}/ToulousePoints-examples.shp
-out.vd ${TEMP}/apTvPrVectorDataReprojectionFromMapToGeo.kml
- -out.proj user
+ -out.proj user
-out.proj.user.map wgs
VALID --compare-ogr ${NOTOL}
${BASELINE_FILES}/prTvVectorDataProjectionFilterFromMapToGeo.kml
${TEMP}/apTvPrVectorDataReprojectionFromMapToGeo.kml)
-
+
OTB_TEST_APPLICATION(NAME apTvPrVectorDataReprojectionFromMapToMap
APP VectorDataReprojection
OPTIONS -in.vd ${INPUTDATA}/ToulousePoints-examples.shp
-out.vd ${TEMP}/apTvPrVectorDataReprojectionFromMapToMap.kml
- -out.proj user
+ -out.proj user
-out.proj.user.map lambert93
VALID --compare-ascii ${NOTOL}
${BASELINE_FILES}/prTvVectorDataProjectionFilterFromMapToMap.kml
@@ -130,20 +132,20 @@ OTB_TEST_APPLICATION(NAME apTvPrVectorDataReprojectionFromMapToMap2
APP VectorDataReprojection
OPTIONS -in.vd ${BASELINE_FILES}/prTvVectorDataProjectionFilterFromMapToMap.kml
-out.vd ${TEMP}/apTvPrVectorDataReprojectionFromMapToMap2.kml
- -out.proj user
+ -out.proj user
-out.proj.user.map utm
-out.proj.user.map.utm.northhem true
-out.proj.user.map.utm.zone 31
VALID --compare-ogr ${NOTOL}
${BASELINE_FILES}/apTvPrVectorDataReprojectionFromMapToMap2.kml
${TEMP}/apTvPrVectorDataReprojectionFromMapToMap2.kml)
-
+
IF(OTB_DATA_USE_LARGEINPUT)
OTB_TEST_APPLICATION(NAME apTvPrVectorDataReprojectionFromMapToSensor
APP VectorDataReprojection
OPTIONS -in.vd ${INPUTDATA}/ToulousePoints-examples.shp
-out.vd ${TEMP}/apTvPrVectorDataReprojectionFromMapToSensor.shp
- -out.proj image
+ -out.proj image
-out.proj.image.in ${LARGEINPUT}/QUICKBIRD/TOULOUSE/000000128955_01_P001_PAN/02APR01105228-P1BS-000000128955_01_P001.TIF
VALID --compare-ogr ${NOTOL}
${BASELINE_FILES}/prTvVectorDataProjectionFilterFromMapToSensor.shp
@@ -151,23 +153,23 @@ OTB_TEST_APPLICATION(NAME apTvPrVectorDataReprojectionFromMapToSensor
ENDIF(OTB_DATA_USE_LARGEINPUT)
-# With QGIS the two vectordata are equal
+# With QGIS the two vectordata are equal
OTB_TEST_APPLICATION(NAME apTuPrVectorDataReprojectionFromGeoToMap
APP VectorDataReprojection
OPTIONS -in.vd ${INPUTDATA}/ToulousePointsWGS.kml
-out.vd ${TEMP}/apTvPrVectorDataReprojectionFromGeoToMap.shp
- -out.proj image
+ -out.proj image
-out.proj.image.in ${INPUTDATA}/QB_Toulouse_Ortho_PAN.tif)
# VALID --compare-ogr ${NOTOL}
# ${INPUTDATA}/ToulousePoints-examples.shp
# ${TEMP}/apTvPrVectorDataReprojectionFromGeoToMap.shp)
-
-
-
-
\ No newline at end of file
+
+
+
+