diff --git a/Code/Projections/otbRationalTransform.h b/Code/Projections/otbRationalTransform.h
new file mode 100644
index 0000000000000000000000000000000000000000..1fd1b150fa3bac46c57936c03b943395bb9b536d
--- /dev/null
+++ b/Code/Projections/otbRationalTransform.h
@@ -0,0 +1,180 @@
+/*=========================================================================
+
+ Program: ORFEO Toolbox
+ Language: C++
+ Date: $Date$
+ Version: $Revision$
+
+
+ Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
+ See OTBCopyright.txt for details.
+
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the above copyright notices for more information.
+
+=========================================================================*/
+#ifndef __otbRationalTransform_h
+#define __otbRationalTransform_h
+
+#include "itkTransform.h"
+#include "itkExceptionObject.h"
+#include "itkMacro.h"
+
+namespace otb
+{
+/** \class RationalTransform
+ * \brief This class implements a rational transfom
+ *
+ * A rational transform is a quotient of two polynomial functions.
+ *
+ * The degree of the numerator and denominator polynomial functions
+ * can be set using the appropriate setters.
+ *
+ * The number of parameters is then the number of dimensions times
+ * the numerator degree plus one times the denominator degree plus
+ * one.
+ *
+ * Parameters in the parameters vector are in the following order:
+ * dim0num0 dim0num1 ... dim0numN dim0denom0 dim0denom1
+ * ... dim0denomM ... dim1num0 ... dimDdenomM.
+ *
+ * \ingroup Transform
+ **/
+
+template <class TScalarType = double,
+ unsigned int Dimension = 2>
+class ITK_EXPORT RationalTransform : public itk::Transform<TScalarType,Dimension,Dimension>
+{
+public:
+ /** Standard class typedefs. */
+ typedef itk::Transform<TScalarType,Dimension,
+ Dimension> Superclass;
+ typedef RationalTransform Self;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
+
+ typedef typename Superclass::ScalarType ScalarType;
+ typedef itk::Point<ScalarType, Dimension> InputPointType;
+ typedef itk::Point<ScalarType, Dimension> OutputPointType;
+
+ typedef typename Superclass::InverseTransformBasePointer InverseTransformBasePointer;
+
+ /** Method for creation through the object factory. */
+ itkNewMacro(Self);
+
+ /** Run-time type information (and related methods). */
+ itkTypeMacro(RationalTransform, itk::Transform);
+
+ itkStaticConstMacro(SpaceDimension, unsigned int, Dimension);
+
+ /** Set the numerator degree */
+ itkSetMacro(NumeratorDegree,unsigned int);
+
+ /** Get the numerator degree */
+ itkGetMacro(NumeratorDegree,unsigned int);
+
+ /** Set the numerator degree */
+ itkSetMacro(DenominatorDegree,unsigned int);
+
+ /** Get the denominator degree */
+ itkGetMacro(DenominatorDegree,unsigned int);
+
+
+ /** The transform point method */
+ virtual OutputPointType TransformPoint(const InputPointType& point) const
+ {
+ // Check for consistency
+ if(this->GetNumberOfParameters() != this->m_Parameters.size())
+ {
+ {
+ itkExceptionMacro(<<"Wrong number of parameters: found "<<this->m_Parameters.Size()<<", expected "<<this->GetNumberOfParameters());
+ }
+ }
+
+ // Build output
+ OutputPointType outputPoint;
+
+ unsigned int dimensionStride = (m_DenominatorDegree+1)+(m_NumeratorDegree+1);
+
+ // Compute RPC transform
+ for(unsigned int dim = 0; dim < SpaceDimension; ++dim)
+ {
+ // Initialize numerator and denominator
+ TScalarType num = itk::NumericTraits<TScalarType>::Zero;
+ TScalarType denom = itk::NumericTraits<TScalarType>::Zero;
+ TScalarType currentPower = 1.;
+
+ // Compute numerator
+ for(unsigned int numDegree = 0; numDegree <= m_NumeratorDegree; ++numDegree)
+ {
+ num+=this->m_Parameters[dim*dimensionStride+numDegree]*currentPower;
+ currentPower*=point[dim];
+ }
+
+ // Compute denominator
+ currentPower = 1.;
+ for(unsigned int denomDegree = 0; denomDegree <= m_DenominatorDegree; ++denomDegree)
+ {
+ denom+=this->m_Parameters[dim*dimensionStride+m_NumeratorDegree+denomDegree+1]*currentPower;
+ currentPower*=point[dim];
+ }
+
+ // Finally, fill the output
+ outputPoint[dim]=num/denom;
+ }
+
+ // Return the output point
+ return outputPoint;
+ }
+
+ // Get the number of parameters
+ virtual unsigned int GetNumberOfParameters() const
+ {
+ return (m_NumeratorDegree +1 + m_DenominatorDegree+1)*SpaceDimension;
+ }
+
+ // Set parameter method
+ virtual void SetParameters(const typename Superclass::ParametersType & params)
+ {
+ // Check for the appropriate size
+ if(params.Size() != this->GetNumberOfParameters())
+ {
+ itkExceptionMacro(<<"Wrong number of parameters: found "<<params.Size()<<", expected "<<this->GetNumberOfParameters());
+ }
+
+ // Set parameters
+ this->m_Parameters = params;
+ }
+
+protected:
+ RationalTransform() : Superclass(SpaceDimension,16), m_NumeratorDegree(3), m_DenominatorDegree(3)
+ {
+ this->m_Parameters.SetSize(this->GetNumberOfParameters());
+ this->m_Parameters.Fill(0);
+ this->m_Parameters[1] = 1.;
+ this->m_Parameters[4] = 1.;
+ }
+
+ virtual ~RationalTransform() {}
+
+ void PrintSelf(std::ostream& os, itk::Indent indent) const
+ {
+ Superclass::PrintSelf(os,indent);
+ }
+
+private:
+ RationalTransform(const Self &); //purposely not implemented
+ void operator =(const Self&); //purposely not implemented
+
+ // Degree of numerator
+ unsigned int m_NumeratorDegree;
+
+ // Degree of denominator
+ unsigned int m_DenominatorDegree;
+};
+
+} // namespace otb
+
+#endif
diff --git a/Testing/Code/Projections/CMakeLists.txt b/Testing/Code/Projections/CMakeLists.txt
index e3999f4f8b8024e148ad9c4e215fab916888e1e0..764f56b11b30737e777cdec3bffba9deb9c5953f 100644
--- a/Testing/Code/Projections/CMakeLists.txt
+++ b/Testing/Code/Projections/CMakeLists.txt
@@ -1006,7 +1006,21 @@ ADD_TEST(prTvImageToGenericRSOutputParameters ${PROJECTIONS_TESTS4}
)
ENDIF(OTB_DATA_USE_LARGEINPUT)
+#----- otb::RationalTransform ------------------------
+ADD_TEST(prTuRationalTransformNew ${PROJECTIONS_TESTS4}
+otbRationalTransformNew)
+ADD_TEST(prTvRationalTransform ${PROJECTIONS_TESTS4}
+ --compare-ascii ${NOTOL}
+ ${BASELINE_FILES}/otbRationalTransformOutput.txt
+ ${TEMP}/otbRationalTransformOutput.txt
+otbRationalTransform
+${TEMP}/otbRationalTransformOutput.txt
+0 0
+1 1
+10 10
+-10 -10
+)
#=======================================================================================
@@ -1067,6 +1081,7 @@ otbGenericRSResampleImageFilter.cxx
otbElevDatabaseHeightAboveMSLFunction.cxx
otbImageToEnvelopeVectorDataFilter.cxx
otbImageToGenericRSOutputParameters.cxx
+otbRationalTransform.cxx
)
OTB_ADD_EXECUTABLE(otbProjectionsTests1 "${Projections_SRCS1}" "OTBProjections;OTBIO;OTBTesting")
diff --git a/Testing/Code/Projections/otbProjectionsTests4.cxx b/Testing/Code/Projections/otbProjectionsTests4.cxx
index 54d72c19f153b428f000d11687d3b477244dc2e8..de3ea5481bd1c572266fc24df1f7c15242b23531 100644
--- a/Testing/Code/Projections/otbProjectionsTests4.cxx
+++ b/Testing/Code/Projections/otbProjectionsTests4.cxx
@@ -35,4 +35,6 @@ void RegisterTests()
REGISTER_TEST(otbImageToEnvelopeVectorDataFilter);
REGISTER_TEST(otbImageToGenericRSOutputParametersNew);
REGISTER_TEST(otbImageToGenericRSOutputParameters);
+ REGISTER_TEST(otbRationalTransformNew);
+ REGISTER_TEST(otbRationalTransform);
}
diff --git a/Testing/Code/Projections/otbRationalTransform.cxx b/Testing/Code/Projections/otbRationalTransform.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..21d721dbc28935a64680da56168f37844b03373b
--- /dev/null
+++ b/Testing/Code/Projections/otbRationalTransform.cxx
@@ -0,0 +1,101 @@
+/*=========================================================================
+
+ Program: ORFEO Toolbox
+ Language: C++
+ Date: $Date$
+ Version: $Revision$
+
+
+ Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
+ See OTBCopyright.txt for details.
+
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the above copyright notices for more information.
+
+=========================================================================*/
+#if defined(_MSC_VER)
+#pragma warning ( disable : 4786 )
+#endif
+
+#include "otbRationalTransform.h"
+#include <fstream>
+
+int otbRationalTransformNew(int argc, char* argv[])
+{
+ typedef otb::RationalTransform<> RationalTransformType;
+
+ // Instantiation
+ RationalTransformType::Pointer rt = RationalTransformType::New();
+
+ return EXIT_SUCCESS;
+}
+
+int otbRationalTransform(int argc, char* argv[])
+{
+ typedef otb::RationalTransform<> RationalTransformType;
+
+ // Instantiation
+ RationalTransformType::Pointer rt = RationalTransformType::New();
+ rt->SetNumeratorDegree(4);
+ rt->SetDenominatorDegree(4);
+
+ RationalTransformType::ParametersType params(rt->GetNumberOfParameters());
+ params.Fill(1.);
+
+ // Rational is
+ // fx(x,y) = (1+2*x+3*x^2+4*x^3+5*x^4)/(6+7*x+8*x^2+9*x^3+10*x^4)
+ // fy(x,y) = (11+12*y+13*y^2+14*y^3+15*y^4)/(16+17*y+18*y^2+19*y^3+20*y^4)
+ params[0]=1;
+ params[1]=2;
+ params[2]=3;
+ params[3]=4;
+ params[4]=5;
+ params[5]=6;
+ params[6]=7;
+ params[7]=8;
+ params[8]=9;
+ params[9]=10;
+ params[10]=11;
+ params[11]=12;
+ params[12]=13;
+ params[13]=14;
+ params[14]=15;
+ params[15]=16;
+ params[16]=17;
+ params[17]=18;
+ params[18]=19;
+ params[19]=20;
+
+ rt->SetParameters(params);
+
+ RationalTransformType::InputPointType inputPoint;
+ RationalTransformType::OutputPointType outputPoint;
+
+ std::ofstream ofs;
+ ofs.open(argv[1]);
+
+ // Set floatfield to format writing properly
+ ofs.setf(std::ios::fixed, std::ios::floatfield);
+ ofs.precision(10);
+
+ unsigned int idx = 2;
+
+ ofs<<"Rational function is: "<<std::endl;
+ ofs<<"fx(x,y) = (1+2*x+3*x^2+4*x^3+5*x^4)/(6+7*x+8*x^2+9*x^3+10*x^4)"<<std::endl;
+ ofs<<"fy(x,y) = (11+12*y+13*y^2+14*y^3+15*y^4)/(16+17*y+18*y^2+19*y^3+20*y^4)"<<std::endl;
+ while(idx+1<(unsigned int)argc)
+ {
+ inputPoint[0] = atof(argv[idx]);
+ inputPoint[1] = atof(argv[idx+1]);
+ outputPoint = rt->TransformPoint(inputPoint);
+ ofs<<inputPoint<<" -> "<<outputPoint<<std::endl;
+ idx+=2;
+ }
+
+ ofs.close();
+
+ return EXIT_SUCCESS;
+}
+