[Insight-developers] inverse transformations
Paul Koshevoy
koshevoy at sci.utah.edu
Fri Oct 7 12:33:44 EDT 2005
Greetings,
Perhaps I wasn't clear enough in my original posting, I'll try to give
you some background information here. I have written an automated image
mosaicing application that requires some minor modification to the
itk::Transform and several of the derived classes. I've been
porting these changes from ITK 1.8.1 to 2.0.1 to 2.2.0 by hand.
When porting to ITK 2.2.0 I've noticed a comment in the
itkAffineTransform.h that the BackTransform API is deprecated. It seems
there has been an effort underway to standardize the inverse
transformation API which resulted in the GetInverse method.
My original posting has pointed out two of my problems with the
GetInverse method. You can compare my changes to itkTransform.h and
itkAffineTransform.h against ITK 2.2.0, my version of those files is in
the attachment.
I would like to know what are the chances of my changes being merged into
the main ITK tree.
Thank you,
Paul.
On Thu, 29 Sep 2005, Paul Koshevoy wrote:
> Date: Thu, 29 Sep 2005 18:08:13 -0600 (MDT)
> From: Paul Koshevoy <koshevoy at sci.utah.edu>
> To: insight-developers at itk.org
> Subject: [Insight-developers] inverse transformations
>
> Hi,
>
> I am working on an automatic image mosaicing problem with
> Ross Whitaker and Tolga Tasdizen.
>
> I am having some problems with the current implementation of the
> inverse transformation -- the GetInverse method.
>
> First, GetInverse is not a virtual function. Therefore, the
> transform type must be known in order for the right GetInverse
> to be called, which means it must be specified as a template
> parameter.
>
> Second, GetInverse is not meaningful for transforms where no general
> inverse exists. For example, the inverse of a radial distortion
> transform may require a different number of parameters than the
> forward transform.
>
> For my application I need a general method for inverse transformation
> of a point. The way I would go about this is by declaring
>
> virtual bool
> GetInverse(itk::Transform<TScalarType,
> NInputDimensions,
> NOutputDimensions>::Pointer & inverse) const;
>
> The requirement that the inverse transform and the forward transform
> must be of the same type is too restrictive. The transforms which do
> not have a general analytic inverse should be able to return an
> iterative numeric inverse wrapper class instead.
>
> An alternative solution is to provide BackTransformPoint,
> BackTransformVector and BackTransformCovariantVector as virtual member
> functions of the itk::Transform base class. These (currently deprecated)
> methods are already provided by a number of transforms. Transforms that do
> not have a general inverse may use an iterative numeric inverse method
> internally. This is the approached I have used so far (with minor
> alterations to ITK 1.8.1 and 2.0.1).
>
> Are these observations accurate and are my proposed solutions acceptable?
>
> Paul.
> _______________________________________________
> Insight-developers mailing list
> Insight-developers at itk.org
> http://www.itk.org/mailman/listinfo/insight-developers
>
-------------- next part --------------
/*
Program: Insight Segmentation & Registration Toolkit
Module: $RCSfile: itkAffineTransform.h,v $
Language: C++
Date: $Date: 2005/03/09 16:46:12 $
Version: $Revision: 1.60 $
Copyright (c) Insight Software Consortium. All rights reserved.
See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm 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 __itkAffineTransform_h
#define __itkAffineTransform_h
#include <iostream>
#include "itkMatrix.h"
#include "itkMatrixOffsetTransformBase.h"
#include "itkExceptionObject.h"
#include "itkMacro.h"
namespace itk
{
/**
* Affine transformation of a vector space (e.g. space coordinates)
*
* This class allows the definition and manipulation of affine
* transformations of an n-dimensional affine space (and its
* associated vector space) onto itself. One common use is to define
* and manipulate Euclidean coordinate transformations in two and
* three dimensions, but other uses are possible as well.
*
* An affine transformation is defined mathematically as a linear
* transformation plus a constant offset. If A is a constant n x n
* matrix and b is a constant n-vector, then y = Ax+b defines an
* affine transformation from the n-vector x to the n-vector y.
*
* The difference between two points is a vector and transforms
* linearly, using the matrix only. That is, (y1-y2) = A*(x1-x2).
*
* The AffineTransform class determines whether to transform an object
* as a point or a vector by examining its type. An object of type
* Point transforms as a point; an object of type Vector transforms as
* a vector.
*
* One common use of affine transformations is to define coordinate
* conversions in two- and three-dimensional space. In this
* application, x is a two- or three-dimensional vector containing the
* "source" coordinates of a point, y is a vector containing the
* "target" coordinates, the matrix A defines the scaling and rotation
* of the coordinate systems from the source to the target, and b
* defines the translation of the origin from the source to the
* target. More generally, A can also define anisotropic scaling and
* shearing transformations. Any good textbook on computer graphics
* will discuss coordinate transformations in more detail. Several of
* the methods in this class are designed for this purpose and use the
* language appropriate to coordinate conversions.
*
* Any two affine transformations may be composed and the result is
* another affine transformation. However, the order is important.
* Given two affine transformations T1 and T2, we will say that
* "precomposing T1 with T2" yields the transformation which applies
* T1 to the source, and then applies T2 to that result to obtain the
* target. Conversely, we will say that "postcomposing T1 with T2"
* yields the transformation which applies T2 to the source, and then
* applies T1 to that result to obtain the target. (Whether T1 or T2
* comes first lexicographically depends on whether you choose to
* write mappings from right-to-left or vice versa; we avoid the whole
* problem by referring to the order of application rather than the
* textual order.)
*
* There are two template parameters for this class:
*
* ScalarT The type to be used for scalar numeric values. Either
* float or double.
*
* NDimensions The number of dimensions of the vector space.
*
* This class provides several methods for setting the matrix and vector
* defining the transform. To support the registration framework, the
* transform parameters can also be set as an Array<double> of size
* (NDimension + 1) * NDimension using method SetParameters().
* The first (NDimension x NDimension) parameters defines the matrix in
* column-major order (where the column index) varies the fastest).
* The last NDimension parameters defines the translation
* in each dimensions.
*
* This class also supports the specification of a center of rotation (center)
* and a translation that is applied with respect to that centered rotation.
* By default the center of rotation is set to the origin.
*
* \ingroup Transforms
*
**/
template <
class TScalarType=double, // Data type for scalars
// (e.g. float or double)
unsigned int NDimensions=3> // Number of dimensions in the input space
class AffineTransform
: public MatrixOffsetTransformBase< TScalarType, NDimensions, NDimensions >
{
public:
/** Standard typedefs */
typedef AffineTransform Self;
typedef MatrixOffsetTransformBase< TScalarType,
NDimensions,
NDimensions > Superclass;
typedef SmartPointer<Self> Pointer;
typedef SmartPointer<const Self> ConstPointer;
/** Run-time type information (and related methods). */
itkTypeMacro( AffineTransform, MatrixOffsetTransformBase );
/** New macro for creation of through a Smart Pointer */
itkNewMacro( Self );
/** Dimension of the domain space. */
itkStaticConstMacro(InputSpaceDimension, unsigned int, NDimensions);
itkStaticConstMacro(OutputSpaceDimension, unsigned int, NDimensions);
itkStaticConstMacro(SpaceDimension, unsigned int, NDimensions);
itkStaticConstMacro(ParametersDimension, unsigned int,
NDimensions*(NDimensions+1));
/** Parameters Type */
typedef typename Superclass::ParametersType ParametersType;
typedef typename Superclass::JacobianType JacobianType;
typedef typename Superclass::ScalarType ScalarType;
typedef typename Superclass::InputPointType InputPointType;
typedef typename Superclass::OutputPointType OutputPointType;
typedef typename Superclass::InputVectorType InputVectorType;
typedef typename Superclass::OutputVectorType OutputVectorType;
typedef typename Superclass::InputVnlVectorType InputVnlVectorType;
typedef typename Superclass::OutputVnlVectorType OutputVnlVectorType;
typedef typename Superclass::InputCovariantVectorType
InputCovariantVectorType;
typedef typename Superclass::OutputCovariantVectorType
OutputCovariantVectorType;
typedef typename Superclass::MatrixType MatrixType;
typedef typename Superclass::InverseMatrixType InverseMatrixType;
typedef typename Superclass::CenterType CenterType;
typedef typename Superclass::OffsetType OffsetType;
typedef typename Superclass::TranslationType TranslationType;
/** Compose affine transformation with a translation
*
* This method modifies self to include a translation of the
* origin. The translation is precomposed with self if pre is
* true, and postcomposed otherwise.
* This updates Translation based on current center. */
void Translate(const OutputVectorType &offset, bool pre=0);
/** Compose affine transformation with a scaling
*
* This method modifies self to magnify the source by a given
* factor along each axis. If all factors are the same, or only a
* single factor is given, then the scaling is isotropic;
* otherwise it is anisotropic. If an odd number of factors are
* negative, then the parity of the image changes. If any of the
* factors is zero, then the transformation becomes a projection
* and is not invertible. The scaling is precomposed with self if
* pre is true, and postcomposed otherwise.
* Note that the scaling is applied centered at the origin. */
void Scale(const OutputVectorType &factor, bool pre=0);
void Scale(const TScalarType &factor, bool pre=0);
/** Compose affine transformation with an elementary rotation
*
* This method composes self with a rotation that affects two
* specified axes, replacing the current value of self. The
* rotation angle is in radians. The axis of rotation goes
* through the origin. The transformation is given by
*
* y[axis1] = cos(angle)*x[axis1] + sin(angle)*x[axis2]
* y[axis2] = -sin(angle)*x[axis1] + cos(angle)*x[axis2].
*
* All coordinates other than axis1 and axis2 are unchanged;
* a rotation of pi/2 radians will carry +axis1 into +axis2.
* The rotation is precomposed with self if pre is true, and
* postcomposed otherwise.
* Note that the rotation is applied centered at the origin. */
void Rotate(int axis1, int axis2, TScalarType angle, bool pre=0);
/** Compose 2D affine transformation with a rotation
*
* This method composes self, which must be a 2D affine
* transformation, with a clockwise rotation through a given angle
* in radians. The center of rotation is the origin. The
* rotation is precomposed with self if pre is true, and
* postcomposed otherwise.
* Note that the rotation is applied centered at the origin.
*
* \warning Only to be use in two dimensions
*
* \todo Find a way to generate a compile-time error
* is this is used with NDimensions != 2. */
void Rotate2D(TScalarType angle, bool pre=0);
/** Compose 3D affine transformation with a rotation
*
* This method composes self, which must be a 3D affine
* transformation, with a clockwise rotation around a specified
* axis. The rotation angle is in radians; the axis of rotation
* goes through the origin. The rotation is precomposed with self
* if pre is true, and postcomposed otherwise.
* Note that the rotation is applied centered at the origin.
*
* \warning Only to be used in dimension 3
*
* \todo Find a way to generate a compile-time error
* is this is used with NDimensions != 3. */
void Rotate3D(const OutputVectorType &axis, TScalarType angle, bool pre=0);
/** Compose affine transformation with a shear
*
* This method composes self with a shear transformation,
* replacing the original contents of self. The shear is
* precomposed with self if pre is true, and postcomposed
* otherwise. The transformation is given by
*
* y[axis1] = x[axis1] + coef*x[axis2]
* y[axis2] = x[axis2].
*
* Note that the shear is applied centered at the origin. */
void Shear(int axis1, int axis2, TScalarType coef, bool pre=0);
/** Back transform by an affine transformation
*
* This method finds the point or vector that maps to a given
* point or vector under the affine transformation defined by
* self. If no such point exists, an exception is thrown.
*
* \deprecated Please use GetInverseTransform and then call the
* forward transform function **/
inline InputPointType BackTransform(const OutputPointType &point ) const;
inline InputVectorType BackTransform(const OutputVectorType &vector) const;
inline InputVnlVectorType BackTransform(
const OutputVnlVectorType &vector) const;
inline InputCovariantVectorType BackTransform(
const OutputCovariantVectorType &vector) const;
/** virtual:
* Inverse transformations - If y = Transform(x), then x = BackTransform(y).
* If no mapping from y to x exists, then an exception is thrown.
*/
InputPointType
BackTransformPoint(const OutputPointType & y) const;
InputVectorType
BackTransformVector(const OutputVectorType & y) const
{ return BackTransform(y); }
InputVnlVectorType
BackTransformVector(const OutputVnlVectorType & y) const
{ return BackTransform(y); }
InputCovariantVectorType
BackTransformCovariantVector(const OutputCovariantVectorType & y) const
{ return BackTransform(y); }
/** Compute distance between two affine transformations
*
* This method computes a ``distance'' between two affine
* transformations. This distance is guaranteed to be a metric,
* but not any particular metric. (At the moment, the algorithm
* is to collect all the elements of the matrix and offset into a
* vector, and compute the euclidean (L2) norm of that vector.
* Some metric which could be used to estimate the distance between
* two points transformed by the affine transformation would be
* more useful, but I don't have time right now to work out the
* mathematical details.) */
ScalarType Metric(const Self * other) const;
/** This method computes the distance from self to the identity
* transformation, using the same metric as the one-argument form
* of the Metric() method. **/
ScalarType Metric(void) const;
protected:
/** Construct an AffineTransform object
*
* This method constructs a new AffineTransform object and
* initializes the matrix and offset parts of the transformation
* to values specified by the caller. If the arguments are
* omitted, then the AffineTransform is initialized to an identity
* transformation in the appropriate number of dimensions. **/
AffineTransform(const MatrixType &matrix,
const OutputVectorType &offset);
AffineTransform(unsigned int outputDims,
unsigned int paramDims);
AffineTransform();
/** Destroy an AffineTransform object **/
virtual ~AffineTransform();
/** Print contents of an AffineTransform */
void PrintSelf(std::ostream &s, Indent indent) const;
private:
AffineTransform(const Self & other);
const Self & operator=( const Self & );
}; //class AffineTransform
} // namespace itk
#ifndef ITK_MANUAL_INSTANTIATION
#include "itkAffineTransform.txx"
#endif
#endif /* __itkAffineTransform_h */
-------------- next part --------------
/*=========================================================================
Program: Insight Segmentation & Registration Toolkit
Module: $RCSfile: itkTransform.h,v $
Language: C++
Date: $Date: 2005/04/05 15:41:44 $
Version: $Revision: 1.47 $
Copyright (c) Insight Software Consortium. All rights reserved.
See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm 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 __itkTransform_h
#define __itkTransform_h
#include "itkTransformBase.h"
#include "itkPoint.h"
#include "itkVector.h"
#include "itkCovariantVector.h"
#include "vnl/vnl_vector_fixed.h"
#include "itkArray.h"
#include "itkArray2D.h"
#include "itkObjectFactory.h"
namespace itk
{
/** \class Transform
* \brief Transform points and vector from an input space to an output space.
*
* This abstract class define the generic interface for a geometrical
* transformation from one space to another. The class provides methods
* for mapping points, vectors and covariant vectors from the input space
* to the output space.
*
* Given that transformation are not necesarily invertible, this basic
* class does not provide the methods for back transfromation. Back transform
* methods are implemented in derived classes where appropriate.
*
* \par Registration Framework Support
* Typically a Transform class has several methods for setting its
* parameters. For use in the registration framework, the parameters must
* also be represented by an array of doubles to allow communication
* with generic optimizers. The Array of transformation parameters is set using
* the SetParameters() method.
*
* Another requirement of the registration framework is the computation
* of the transform Jacobian. In general, a ImageToImageMetric requires
* the knowledge of the Jacobian in order to compute the metric derivatives.
* The Jacobian is a matrix whose element are the partial derivatives
* of the output point with respect to the array of parameters that defines
* the transform.
*
* \ingroup Transforms
*
*/
template <class TScalarType,
unsigned int NInputDimensions=3,
unsigned int NOutputDimensions=3>
class ITK_EXPORT Transform : public TransformBase
{
public:
/** Standard class typedefs. */
typedef Transform Self;
typedef TransformBase Superclass;
typedef SmartPointer< Self > Pointer;
typedef SmartPointer< const Self > ConstPointer;
/** New method for creating an object using a factory. */
itkNewMacro(Self);
/** Run-time type information (and related methods). */
itkTypeMacro( Transform, TransformBase );
/** Dimension of the domain space. */
itkStaticConstMacro(InputSpaceDimension, unsigned int, NInputDimensions);
itkStaticConstMacro(OutputSpaceDimension, unsigned int, NOutputDimensions);
/** Get the size of the input space */
unsigned int GetInputSpaceDimension(void) const {return NInputDimensions;}
/** Get the size of the output space */
unsigned int GetOutputSpaceDimension(void) const {return NOutputDimensions;}
/** Type of the scalar representing coordinate and vector elements. */
typedef TScalarType ScalarType;
/** Type of the input parameters. */
typedef typename Superclass::ParametersType ParametersType;
/** Type of the Jacobian matrix. */
typedef Array2D< double > JacobianType;
/** Standard vector type for this class. */
typedef Vector<TScalarType, NInputDimensions> InputVectorType;
typedef Vector<TScalarType, NOutputDimensions> OutputVectorType;
/** Standard covariant vector type for this class */
typedef CovariantVector<TScalarType, NInputDimensions> InputCovariantVectorType;
typedef CovariantVector<TScalarType, NOutputDimensions> OutputCovariantVectorType;
/** Standard vnl_vector type for this class. */
typedef vnl_vector_fixed<TScalarType, NInputDimensions> InputVnlVectorType;
typedef vnl_vector_fixed<TScalarType, NOutputDimensions> OutputVnlVectorType;
/** Standard coordinate point type for this class */
typedef Point<TScalarType, NInputDimensions> InputPointType;
typedef Point<TScalarType, NOutputDimensions> OutputPointType;
/** Method to transform a point. */
virtual OutputPointType TransformPoint(const InputPointType & ) const
{ return OutputPointType(); }
/** Method to transform a vector. */
virtual OutputVectorType TransformVector(const InputVectorType &) const
{ return OutputVectorType(); }
/** Method to transform a vnl_vector. */
virtual OutputVnlVectorType TransformVector(const InputVnlVectorType &) const
{ return OutputVnlVectorType(); }
/** Method to transform a CovariantVector. */
virtual OutputCovariantVectorType TransformCovariantVector(
const InputCovariantVectorType &) const
{ return OutputCovariantVectorType(); }
/**
* Inverse transformations - If y = Transform(x), then x = BackTransform(y).
* If no mapping from y to x exists, then an exception is thrown.
*/
virtual InputPointType
BackTransformPoint(const OutputPointType & ) const
{
itkExceptionMacro(<< "BackTransformPoint not implemented. ");
return InputPointType();
}
/**
* Inverse transformations - If y = Transform(x), then x = BackTransform(y).
* If no mapping from y to x exists, then an exception is thrown.
*/
virtual InputVectorType
BackTransformVector(const OutputVectorType &) const
{
itkExceptionMacro(<< "BackTransformVector not implemented. ");
return InputVectorType();
}
/**
* Inverse transformations - If y = Transform(x), then x = BackTransform(y).
* If no mapping from y to x exists, then an exception is thrown.
*/
virtual InputVnlVectorType
BackTransformVector(const OutputVnlVectorType &) const
{
itkExceptionMacro(<< "BackTransformVector not implemented. ");
return InputVnlVectorType();
}
/**
* Inverse transformations - If y = Transform(x), then x = BackTransform(y).
* If no mapping from y to x exists, then an exception is thrown.
*/
virtual InputCovariantVectorType
BackTransformCovariantVector(const OutputCovariantVectorType &) const
{
itkExceptionMacro(<< "BackTransformCovariantVector not implemented. ");
return InputCovariantVectorType();
}
/** Set the transformation parameters and update internal transformation.
* SetParameters gives the transform the option to set it's
* parameters by keeping a reference to the parameters, or by
* copying. To force the transform to copy it's parameters call
* SetParametersByValue.
* \sa SetParametersByValue
*/
virtual void SetParameters( const ParametersType & )
{ itkExceptionMacro( << "Subclasses should override this method" ) };
/** Set the transformation parameters and update internal transformation.
* This method forces the transform to copy the parameters. The
* default implementation is to call SetParameters. This call must
* be overridden if the transform normally implements SetParameters
* by keeping a reference to the parameters.
* \sa SetParameters
*/
virtual void SetParametersByValue ( const ParametersType & p )
{ this->SetParameters ( p ); };
/** Get the Transformation Parameters. */
virtual const ParametersType& GetParameters(void) const
{ itkExceptionMacro( << "Subclasses should override this method" );
return m_Parameters; };
/** Set the fixed parameters and update internal transformation. */
virtual void SetFixedParameters( const ParametersType & )
{ itkExceptionMacro( << "Subclasses should override this method" ) };
/** Get the Fixed Parameters. */
virtual const ParametersType& GetFixedParameters(void) const
{ itkExceptionMacro( << "Subclasses should override this method" );
return m_Parameters; };
/** Compute the Jacobian of the transformation
*
* This method computes the Jacobian matrix of the transformation
* at a given input point. The rank of the Jacobian will also indicate
* if the transform is invertible at this point.
*
* The Jacobian is be expressed as a matrix of partial derivatives of the
* output point components with respect to the parameters that defined
* the transform:
*
* \f[
*
J=\left[ \begin{array}{cccc}
\frac{\partial x_{1}}{\partial p_{1}} &
\frac{\partial x_{1}}{\partial p_{2}} &
\cdots & \frac{\partial x_{1}}{\partial p_{m}}\\
\frac{\partial x_{2}}{\partial p_{1}} &
\frac{\partial x_{2}}{\partial p_{2}} &
\cdots & \frac{\partial x_{2}}{\partial p_{m}}\\
\vdots & \vdots & \ddots & \vdots \\
\frac{\partial x_{n}}{\partial p_{1}} &
\frac{\partial x_{n}}{\partial p_{2}} &
\cdots & \frac{\partial x_{n}}{\partial p_{m}}
\end{array}\right]
*
* \f]
* **/
virtual const JacobianType & GetJacobian(const InputPointType &) const
{ itkExceptionMacro( << "Subclass should override this method" );
return m_Jacobian; };
/** Return the number of parameters that completely define the Transfom */
virtual unsigned int GetNumberOfParameters(void) const
{ return m_Parameters.Size(); }
/** Return the inverse of the transform.
* The inverse is recomputed if it has been modified */
bool GetInverse(Self*) const {return false;}
/** Generate a platform independant name */
virtual std::string GetTransformTypeAsString() const;
protected:
Transform();
Transform(unsigned int Dimension, unsigned int NumberOfParameters);
virtual ~Transform() {};
mutable ParametersType m_Parameters;
mutable ParametersType m_FixedParameters;
mutable JacobianType m_Jacobian;
private:
Transform(const Self&); //purposely not implemented
void operator=(const Self&); //purposely not implemented
};
} // end namespace itk
#ifndef ITK_MANUAL_INSTANTIATION
#include "itkTransform.txx"
#endif
#endif
More information about the Insight-developers
mailing list