/*========================================================================= * * Copyright NumFOCUS * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. * You may obtain a copy of the License at * * https://www.apache.org/licenses/LICENSE-2.0.txt * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. * *=========================================================================*/ #ifndef itkTransform_h #define itkTransform_h #include // For std::enable_if #include "itkTransformBase.h" #include "itkVector.h" #include "itkSymmetricSecondRankTensor.h" #include "itkDiffusionTensor3D.h" #include "itkVariableLengthVector.h" #include "vnl/vnl_vector_fixed.h" #include "vnl/vnl_matrix_fixed.h" #include "itkMatrix.h" namespace itk { /** * \class Transform * \brief Transform points and vectors from an input space to an output space. * * This abstract class defines the generic interface for a geometric * 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 transformations are not necessarily invertible, this basic * class does not provide the methods for back transformation. 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, an 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. * * Subclasses must provide implementations for:
* virtual OutputPointType TransformPoint(const InputPointType &) const
* virtual OutputVectorType TransformVector(const InputVectorType &) const
* virtual OutputVnlVectorType TransformVector(const InputVnlVectorType &) const
* virtual OutputCovariantVectorType TransformCovariantVector(const InputCovariantVectorType &) const
* virtual void SetParameters(const ParametersType &)
* virtual void SetFixedParameters(const FixedParametersType &)
* virtual void ComputeJacobianWithRespectToParameters( * const InputPointType &, * JacobianType &) const
* virtual void ComputeJacobianWithRespectToPosition( * const InputPointType & x, * JacobianPositionType &jacobian ) const;
* * Since TransformVector and TransformCovariantVector have multiple * overloaded methods from the base class, subclasses must specify:
* using Superclass::TransformVector;
* using Superclass::TransformCovariantVector;
* * * \ingroup ITKTransform */ template class ITK_TEMPLATE_EXPORT Transform : public TransformBaseTemplate { public: ITK_DISALLOW_COPY_AND_MOVE(Transform); /** Standard class type aliases. */ using Self = Transform; using Superclass = TransformBaseTemplate; using Pointer = SmartPointer; using ConstPointer = SmartPointer; /** \see LightObject::GetNameOfClass() */ itkOverrideGetNameOfClassMacro(Transform); /** Dimension of the domain space. */ static constexpr unsigned int InputSpaceDimension = VInputDimension; static constexpr unsigned int OutputSpaceDimension = VOutputDimension; /** define the Clone method */ itkCloneMacro(Self); /** Get the size of the input space */ unsigned int GetInputSpaceDimension() const override { return VInputDimension; } /** Get the size of the output space */ unsigned int GetOutputSpaceDimension() const override { return VOutputDimension; } /* For storing the name of InputSpace */ itkSetMacro(InputSpaceName, std::string); itkGetConstReferenceMacro(InputSpaceName, std::string); /** For storing the name of InputSpace/OutputSpace. InputSpaceName, OutputSpaceName provide identifiers for the world spaces that the transform applied to and the direction of the spatial transformation. The direction of the transform goes from the input space to output space. Typical values include the names of an atlas or a dataset. */ itkSetMacro(OutputSpaceName, std::string); itkGetConstReferenceMacro(OutputSpaceName, std::string); /** Type of the input parameters. */ using typename Superclass::FixedParametersType; using typename Superclass::FixedParametersValueType; using typename Superclass::ParametersType; using typename Superclass::ParametersValueType; using DerivativeType = Array; /** Type of the scalar representing coordinate and vector elements. */ using ScalarType = ParametersValueType; /** Type of the Jacobian matrix. */ using JacobianType = Array2D; using JacobianPositionType = vnl_matrix_fixed; using InverseJacobianPositionType = vnl_matrix_fixed; /** Standard vector type for this class. */ using InputVectorType = Vector; using OutputVectorType = Vector; /** Standard variable length vector type for this class * this provides an interface for the VectorImage class */ using InputVectorPixelType = VariableLengthVector; using OutputVectorPixelType = VariableLengthVector; /* Standard symmetric second rank tensor type for this class */ using InputSymmetricSecondRankTensorType = SymmetricSecondRankTensor; using OutputSymmetricSecondRankTensorType = SymmetricSecondRankTensor; /* Standard tensor type for this class */ using InputDiffusionTensor3DType = DiffusionTensor3D; using OutputDiffusionTensor3DType = DiffusionTensor3D; /** Standard covariant vector type for this class */ using InputCovariantVectorType = CovariantVector; using OutputCovariantVectorType = CovariantVector; /** Standard vnl_vector type for this class. */ using InputVnlVectorType = vnl_vector_fixed; using OutputVnlVectorType = vnl_vector_fixed; /** Standard coordinate point type for this class */ using InputPointType = Point; using OutputPointType = Point; /** Base inverse transform type. This type should not be changed to the * concrete inverse transform type or inheritance would be lost. */ using InverseTransformBaseType = Transform; using InverseTransformBasePointer = typename InverseTransformBaseType::Pointer; using MatrixType = Matrix; using OutputDirectionMatrix = Matrix; using InputDirectionMatrix = Matrix; using DirectionChangeMatrix = Matrix; using typename Superclass::NumberOfParametersType; /** Method to transform a point. * \warning This method must be thread-safe. See, e.g., its use * in ResampleImageFilter. */ virtual OutputPointType TransformPoint(const InputPointType &) const = 0; /** Method to transform a vector. */ virtual OutputVectorType TransformVector(const InputVectorType &) const { itkExceptionMacro("TransformVector(const InputVectorType &)" "is unimplemented for " << this->GetNameOfClass()); } /** Method to transform a vector at a given location. * For global transforms, \c point is ignored and \c TransformVector( vector ) * is called. Local transforms (e.g. deformation * field transform) must override and provide required behavior. */ virtual OutputVectorType TransformVector(const InputVectorType & vector, const InputPointType & point) const; /** Method to transform a vnl_vector. */ virtual OutputVnlVectorType TransformVector(const InputVnlVectorType &) const { itkExceptionMacro("TransformVector( const InputVnlVectorType & ) is " "unimplemented for " << this->GetNameOfClass()); } /** Method to transform a vnl_vector, at a point. * For global transforms, \c point is ignored and \c TransformVector( vector ) * is called. Local transforms (e.g. deformation * field transform) must override and provide required behavior. */ virtual OutputVnlVectorType TransformVector(const InputVnlVectorType & vector, const InputPointType & point) const; /** Method to transform a vector stored in a VectorImage. */ virtual OutputVectorPixelType TransformVector(const InputVectorPixelType & itkNotUsed(vector)) const { itkExceptionMacro("TransformVector( const InputVectorPixelType & ) is " "unimplemented for " << this->GetNameOfClass()); } /** Method to transform a vector stored in a VectorImage, at a point. * For global transforms, \c point is ignored and \c TransformVector( vector ) * is called. Local transforms (e.g. deformation * field transform) must override and provide required behavior. */ virtual OutputVectorPixelType TransformVector(const InputVectorPixelType & vector, const InputPointType & point) const; /** Method to transform a CovariantVector. */ virtual OutputCovariantVectorType TransformCovariantVector(const InputCovariantVectorType &) const { itkExceptionMacro("TransformCovariantVector( const InputCovariantVectorType & ) is " "unimplemented for " << this->GetNameOfClass()); } /** Method to transform a CovariantVector, using a point. Global transforms * can ignore the \c point parameter. Local transforms (e.g. deformation * field transform) must override and provide required behavior. * By default, \c point is ignored and * \c TransformCovariantVector(vector) is called */ virtual OutputCovariantVectorType TransformCovariantVector(const InputCovariantVectorType & vector, const InputPointType & point) const; /** Method to transform a CovariantVector stored in a VectorImage. */ virtual OutputVectorPixelType TransformCovariantVector(const InputVectorPixelType & itkNotUsed(vector)) const { itkExceptionMacro("TransformCovariantVector(const InputVectorPixelType &)" "is unimplemented for " << this->GetNameOfClass()); } /** Method to transform a CovariantVector, using a point. Global transforms * can ignore the \c point parameter. Local transforms (e.g. deformation * field transform) must override and provide required behavior. * By default, \c point is ignored and \c TransformCovariantVector(vector) is * called */ virtual OutputVectorPixelType TransformCovariantVector(const InputVectorPixelType & vector, const InputPointType & point) const; /** Method to transform a diffusion tensor */ virtual OutputDiffusionTensor3DType TransformDiffusionTensor3D(const InputDiffusionTensor3DType & itkNotUsed(tensor)) const { itkExceptionMacro("TransformDiffusionTensor3D( const InputDiffusionTensor3DType & ) is " "unimplemented for " << this->GetNameOfClass()); } /** Method to transform a diffusion tensor at a point. Global transforms * can ignore the \c point parameter. Local transforms (e.g. deformation * field transform) must override and provide required behavior. * By default, \c point is ignored and \c TransformDiffusionTensor(tensor) is * called */ virtual OutputDiffusionTensor3DType TransformDiffusionTensor3D(const InputDiffusionTensor3DType & inputTensor, const InputPointType & point) const; /** Method to transform a diffusion tensor stored in a VectorImage */ virtual OutputVectorPixelType TransformDiffusionTensor3D(const InputVectorPixelType & itkNotUsed(tensor)) const { itkExceptionMacro("TransformDiffusionTensor( const InputVectorPixelType & ) is " "unimplemented for " << this->GetNameOfClass()); } virtual OutputVectorPixelType TransformDiffusionTensor3D(const InputVectorPixelType & inputTensor, const InputPointType & point) const; /** Method to transform a diffusion tensor at a point. Global transforms * can ignore the \c point parameter. Local transforms (e.g. deformation * field transform) must override and provide required behavior. * By default, \c point is ignored and \c TransformSymmetricSecondRankTensor(tensor) is * called */ virtual OutputSymmetricSecondRankTensorType TransformSymmetricSecondRankTensor(const InputSymmetricSecondRankTensorType & inputTensor, const InputPointType & point) const; /** Method to transform a ssr tensor stored in a VectorImage */ virtual OutputSymmetricSecondRankTensorType TransformSymmetricSecondRankTensor(const InputSymmetricSecondRankTensorType & itkNotUsed(tensor)) const { itkExceptionMacro("TransformSymmetricSecondRankTensor( const InputSymmetricSecondRankTensorType & ) is " "unimplemented for " << this->GetNameOfClass()); } /** Method to transform a ssr tensor stored in a VectorImage */ virtual OutputVectorPixelType TransformSymmetricSecondRankTensor(const InputVectorPixelType & itkNotUsed(tensor)) const { itkExceptionMacro("TransformSymmetricSecondRankTensor( const InputVectorPixelType & ) is " "unimplemented for " << this->GetNameOfClass()); } /** Method to transform a diffusion tensor stored in a VectorImage, at * a point. Global transforms * can ignore the \c point parameter. Local transforms (e.g. deformation * field transform) must override and provide required behavior. * By default, \c point is ignored and \c TransformDiffusionTensor(tensor) is * called */ virtual OutputVectorPixelType TransformSymmetricSecondRankTensor(const InputVectorPixelType & inputTensor, const InputPointType & point) const; /** 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 its parameters call * SetParametersByValue. * \sa SetParametersByValue */ void SetParameters(const ParametersType &) override = 0; /** 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 */ void SetParametersByValue(const ParametersType & p) override { this->SetParameters(p); } /** This function allow copying a range of values into the Parameters * The range of values must conform to std::copy(begin, end, m_Parameters) * requirements. */ void CopyInParameters(const ParametersValueType * const begin, const ParametersValueType * const end) override; /** This function allow copying a range of values into the FixedParameters * The range of values must conform to std::copy(begin, end, m_FixedParameters) * requirements. */ void CopyInFixedParameters(const FixedParametersValueType * const begin, const FixedParametersValueType * const end) override; /** Get the Transformation Parameters. */ const ParametersType & GetParameters() const override { return m_Parameters; } /** Set the fixed parameters and update internal transformation. */ void SetFixedParameters(const FixedParametersType &) override = 0; /** Get the Fixed Parameters. */ const FixedParametersType & GetFixedParameters() const override { return m_FixedParameters; } /** Update the transform's parameters by the values in \c update. * \param update must be of the same length as returned by * GetNumberOfParameters(). Throw an exception otherwise. * \param factor is a scalar multiplier for each value in \c update. * SetParameters is called at the end of this method, to allow the transform * to perform any required operations on the updated parameters - typically * a conversion to member variables for use in TransformPoint. */ virtual void UpdateTransformParameters(const DerivativeType & update, ParametersValueType factor = 1.0); /** Return the number of local parameters that completely defines the * Transform at an individual voxel. * For transforms with local support, this will enable downstream * computation of the jacobian wrt only the local support region. * For instance, in the case of a deformation field, this will be equal to * the number of image dimensions. If it is an affine transform, this will * be the same as the GetNumberOfParameters(). */ virtual NumberOfParametersType GetNumberOfLocalParameters() const { return this->GetNumberOfParameters(); } /** Return the number of parameters that completely define the Transform */ NumberOfParametersType GetNumberOfParameters() const override { return this->m_Parameters.Size(); } /** Return the number of parameters that define the constant elements of a Transform */ virtual NumberOfParametersType GetNumberOfFixedParameters() const { return this->m_FixedParameters.Size(); } /** Returns a boolean indicating whether it is possible or not to compute the * inverse of this current Transform. If it is possible, then the inverse of * the transform is returned in the inverseTransform variable passed by the * user. The inverse is recomputed if this current transform has been * modified. * This method is intended to be overridden as needed by derived classes. * */ bool GetInverse(Self * itkNotUsed(inverseTransform)) const { return false; } /** Return an inverse of this transform. If the inverse has not been * implemented, return nullptr. The type of the inverse transform * does not necessarily need to match the type of the forward * transform. This allows one to return a numeric inverse transform * instead. */ virtual InverseTransformBasePointer GetInverseTransform() const { return nullptr; } /** Generate a platform independent name */ std::string GetTransformTypeAsString() const override; using typename Superclass::TransformCategoryEnum; /** Indicates the category transform. * e.g. an affine transform, or a local one, e.g. a deformation field. */ TransformCategoryEnum GetTransformCategory() const override { return Superclass::TransformCategoryEnum::UnknownTransformCategory; } virtual bool IsLinear() const { return (this->GetTransformCategory() == Superclass::TransformCategoryEnum::Linear); } /** * 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] * * This is also used for efficient computation of a point-local jacobian * for dense transforms. * \c jacobian is assumed to be thread-local variable, otherwise memory corruption * will most likely occur during multi-threading. * To avoid repetitive memory allocation, pass in 'jacobian' with its size * already set. */ virtual void ComputeJacobianWithRespectToParameters(const InputPointType & itkNotUsed(p), JacobianType & itkNotUsed(jacobian)) const = 0; virtual void ComputeJacobianWithRespectToParametersCachedTemporaries(const InputPointType & p, JacobianType & jacobian, JacobianType & itkNotUsed(cachedJacobian)) const { // NOTE: default implementation is not optimized, and just falls back to original methods. this->ComputeJacobianWithRespectToParameters(p, jacobian); } /** This provides the ability to get a local jacobian value * in a dense/local transform, e.g. DisplacementFieldTransform. For such * transforms it would be unclear what parameters would refer to. * Generally, global transforms should return an identity jacobian * since there is no change with respect to position. */ virtual void ComputeJacobianWithRespectToPosition(const InputPointType & itkNotUsed(x), JacobianPositionType & itkNotUsed(jacobian)) const { itkExceptionMacro("ComputeJacobianWithRespectToPosition( InputPointType, JacobianType )" " is unimplemented for " << this->GetNameOfClass()); } itkLegacyMacro(virtual void ComputeJacobianWithRespectToPosition(const InputPointType & x, JacobianType & jacobian) const); /** This provides the ability to get a local jacobian value * in a dense/local transform, e.g. DisplacementFieldTransform. For such * transforms it would be unclear what parameters would refer to. * Generally, global transforms should return an identity jacobian * since there is no change with respect to position. */ virtual void ComputeInverseJacobianWithRespectToPosition(const InputPointType & pnt, InverseJacobianPositionType & jacobian) const; itkLegacyMacro(virtual void ComputeInverseJacobianWithRespectToPosition(const InputPointType & x, JacobianType & jacobian) const); /** Apply this transform to an image without resampling. * * Updates image metadata (origin, spacing, direction cosines matrix) in place. * * Only available when input and output space are of the same dimension. * Only works properly for linear transforms. * * The image parameter may be either a SmartPointer or a raw pointer. * */ template std::enable_if_t ApplyToImageMetadata(TImage * image) const; template std::enable_if_t ApplyToImageMetadata(SmartPointer image) const { this->ApplyToImageMetadata(image.GetPointer()); // Delegate to the raw pointer signature } protected: /** * Clone the current transform. * This does a complete copy of the transform * state to the new transform */ typename LightObject::Pointer InternalClone() const override; /** Default-constructor. Creates a transform, having empty `Parameters` and `FixedParameters`. */ Transform() = default; Transform(NumberOfParametersType numberOfParameters); #if defined(__GNUC__) // A bug in some versions of the GCC and Clang compilers // result in an ICE or linker error when "= default" is requested. // This was observed in at least gcc 4.8 and 5.4.0, and // AppleClang 7.0.2 and 8.0.0. Probably others too. // "= default" doesn't gain us much, so just don't use it here. ~Transform() override{}; #else ~Transform() override = default; #endif mutable ParametersType m_Parameters{}; mutable FixedParametersType m_FixedParameters{}; OutputDiffusionTensor3DType PreservationOfPrincipalDirectionDiffusionTensor3DReorientation(const InputDiffusionTensor3DType &, const InverseJacobianPositionType &) const; /** Returns the inverse of the specified transform. Returns null if it cannot invert the transform. Helper function * for the implementation of `GetInverseTransform()` in derived transform classes. */ template static InverseTransformBasePointer InvertTransform(const TTransform & transform) { const auto inverse = TTransform::New(); return transform.GetInverse(inverse) ? inverse.GetPointer() : nullptr; } private: std::string m_InputSpaceName{}; std::string m_OutputSpaceName{}; template static std::string GetTransformTypeAsString(TType *) { std::string rval("other"); return rval; } static std::string GetTransformTypeAsString(float *) { std::string rval("float"); return rval; } static std::string GetTransformTypeAsString(double *) { std::string rval("double"); return rval; } }; } // end namespace itk #ifndef ITK_MANUAL_INSTANTIATION # include "itkTransform.hxx" #endif #endif