registration/clitkBSplineDeformableTransform.h

   1 /*=========================================================================
   2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
   3
   4   Authors belong to:
   5   - University of LYON              http://www.universite-lyon.fr/
   6   - Léon Bérard cancer center       http://www.centreleonberard.fr
   7   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
   8
   9   This software is distributed WITHOUT ANY WARRANTY; without even
  10   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  11   PURPOSE.  See the copyright notices for more information.
  12
  13   It is distributed under dual licence
  14
  15   - BSD        See included LICENSE.txt file
  16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
  17 ===========================================================================**/
  18 #ifndef __clitkBSplineDeformableTransform_h
  19 #define __clitkBSplineDeformableTransform_h
  20 #include "clitkVectorBSplineResampleImageFunctionWithLUT.h"
  21
  22 //itk include
  23 #include "itkTransform.h"
  24 #include "itkImage.h"
  25 #include "itkImageRegion.h"
  26 #include "itkSpatialObject.h"
  27
  28 namespace clitk
  29 {
  30   // Forward declaration needed for friendship
  31   template <class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
  32   class ITK_EXPORT MultipleBSplineDeformableTransform;
  33
  34   template <
  35     class TCoordRep = double,               // Data type for scalars, coordinate representation,vectors
  36     unsigned int NInputDimensions = 3,        // Number of input dimensions
  37     unsigned int NOutputDimensions = 3 >      // Number of output dimensions
  38   class ITK_EXPORT BSplineDeformableTransform :
  39     public itk::Transform< TCoordRep, NInputDimensions, NOutputDimensions >
  40   {
  41   public:
  42
  43     //====================================================================
  44     // Typedefs
  45     //====================================================================
  46     typedef BSplineDeformableTransform                         Self;
  47     typedef itk::Transform< TCoordRep, NInputDimensions, NOutputDimensions > Superclass;
  48     typedef itk::SmartPointer<Self>                                 Pointer;
  49     typedef itk::SmartPointer<const Self>                           ConstPointer;
  50
  51     /** New macro for creation of through the object factory.*/
  52     itkNewMacro( Self );
  53
  54     /** Run-time type information (and related methods). */
  55     itkTypeMacro( BSplineDeformableTransform, Transform );
  56
  57     /** Dimension of the domain space. */
  58     itkStaticConstMacro(OutputDimension, unsigned int, NOutputDimensions);
  59
  60     /** Dimension of the input model. */
  61     itkStaticConstMacro(InputDimension, unsigned int, NInputDimensions);
  62
  63     /** Standard scalar type for this class. */
  64     typedef typename Superclass::ScalarType ScalarType;
  65
  66     /** Standard parameters container. */
  67     typedef typename Superclass::ParametersType ParametersType;
  68
  69     /** Standard Jacobian container. */
  70     typedef typename Superclass::JacobianType JacobianType;
  71
  72     /** Standard vector type for this class. */
  73     typedef itk::Vector<TCoordRep,
  74                         InputDimension> InputVectorType;
  75     typedef itk::Vector<TCoordRep,
  76                         OutputDimension> OutputVectorType;
  77
  78     /** Standard covariant vector type for this class. */
  79     typedef itk::CovariantVector<TCoordRep,
  80                                  InputDimension> InputCovariantVectorType;
  81     typedef itk::CovariantVector<TCoordRep,
  82                                  OutputDimension> OutputCovariantVectorType;
  83
  84     /** Standard vnl_vector type for this class. */
  85     typedef vnl_vector_fixed<TCoordRep,
  86                              InputDimension> InputVnlVectorType;
  87     typedef vnl_vector_fixed<TCoordRep,
  88                              OutputDimension> OutputVnlVectorType;
  89
  90     /** Standard coordinate point type for this class. */
  91     typedef itk::Point<TCoordRep,
  92                        InputDimension> InputPointType;
  93     typedef itk::Point<TCoordRep,
  94                        OutputDimension> OutputPointType;
  95
  96     //JV Parameters as images with OutputDimension number of components per Pixel
  97     typedef typename ParametersType::ValueType                          ParametersValueType;
  98     typedef typename itk::Vector<ParametersValueType, OutputDimension>  PixelType;
  99     typedef itk::Image<PixelType, InputDimension>                       CoefficientImageType;
 100     typedef typename CoefficientImageType::Pointer                      CoefficientImagePointer;
 101
 102
 103     /** Typedefs for specifying the extend to the grid. */
 104     typedef itk::ImageRegion<InputDimension>                RegionType;
 105     typedef typename RegionType::IndexType                  IndexType;
 106     typedef typename RegionType::SizeType                   SizeType;
 107     typedef typename CoefficientImageType::SpacingType      SpacingType;
 108     typedef typename CoefficientImageType::DirectionType    DirectionType;
 109     typedef typename CoefficientImageType::PointType        OriginType;
 110     typedef itk::ContinuousIndex<TCoordRep, InputDimension> ContinuousIndexType;
 111
 112     //JV added for the BLUT interpolator
 113     //typedef itk::Vector<TCoordRep, InputDimension> OutputSpacingType;
 114
 115     //JV m_VectorInterpolator
 116     typedef VectorBSplineResampleImageFunctionWithLUT
 117     <CoefficientImageType, TCoordRep> VectorInterpolatorType;
 118     typedef  typename  VectorInterpolatorType::CoefficientDataType CoefficientDataType;
 119     typedef  typename  VectorInterpolatorType::CoefficientDataType WeightsDataType;
 120
 121     /** Typedef of the bulk transform. */
 122     typedef itk::Transform<ScalarType, InputDimension, OutputDimension> BulkTransformType;
 123     typedef BulkTransformType*                  BulkTransformPointer;
 124
 125     /** Typedef of the mask */
 126     typedef itk::SpatialObject<  InputDimension >   MaskType;
 127     typedef MaskType*    MaskPointer;
 128
 129     //====================================================================
 130     // Set et Gets
 131     //====================================================================
 132     //JV  added for the BLUT interpolator
 133     void SetSplineOrder(const unsigned int &  splineOrder);
 134     void SetSplineOrders(const SizeType &  splineOrders);
 135     itkGetMacro( SplineOrders, SizeType );
 136     itkGetConstMacro( SplineOrders, SizeType );
 137     void SetLUTSamplingFactor(const int &  samplingFactor);
 138     void SetLUTSamplingFactors(const SizeType &  samplingFactors);
 139     itkGetMacro( LUTSamplingFactors, SizeType );
 140     itkGetConstMacro( LUTSamplingFactors,SizeType );
 141     //void SetOutputSpacing(const OutputSpacingType & outputSpacing);
 142     //itkGetMacro( OutputSpacing, OutputSpacingType );
 143     //itkGetConstMacro( OutputSpacing, OutputSpacingType );
 144
 145     void SetParameters(const ParametersType & parameters);
 146
 147     void SetFixedParameters(const ParametersType & parameters);
 148
 149     void SetParametersByValue(const ParametersType & parameters);
 150
 151     void SetIdentity();
 152
 153     /** Get the Transformation Parameters. */
 154     virtual const ParametersType& GetParameters(void) const;
 155
 156     /** Get the Transformation Fixed Parameters. */
 157     virtual const ParametersType& GetFixedParameters(void) const;
 158
 159     // The coefficientImage
 160     virtual CoefficientImagePointer GetCoefficientImage()
 161     { return m_CoefficientImage; }
 162     virtual const CoefficientImagePointer  GetCoefficientImage() const
 163     { return m_CoefficientImage; }
 164     virtual void SetCoefficientImage(CoefficientImagePointer image);
 165
 166     /** This method specifies the region over which the grid resides. */
 167     virtual void SetGridRegion( const RegionType& region );
 168     itkGetMacro( GridRegion, RegionType );
 169     itkGetConstMacro( GridRegion, RegionType );
 170
 171     /** This method specifies the grid spacing or resolution. */
 172     virtual void SetGridSpacing( const SpacingType& spacing );
 173     itkGetMacro( GridSpacing, SpacingType );
 174     itkGetConstMacro( GridSpacing, SpacingType );
 175
 176     /** This method specifies the grid directions . */
 177     virtual void SetGridDirection( const DirectionType & spacing );
 178     itkGetMacro( GridDirection, DirectionType );
 179     itkGetConstMacro( GridDirection, DirectionType );
 180
 181     /** This method specifies the grid origin. */
 182     virtual void SetGridOrigin( const OriginType& origin );
 183     itkGetMacro( GridOrigin, OriginType );
 184     itkGetConstMacro( GridOrigin, OriginType );
 185
 186     // Set the bulk transform, real pointer
 187     // itkSetConstObjectMacro( BulkTransform, BulkTransformType );
 188     // itkGetConstObjectMacro( BulkTransform, BulkTransformType );
 189     void SetBulkTransform(BulkTransformPointer b){m_BulkTransform=b;}
 190     BulkTransformPointer GetBulkTransform(void) {return m_BulkTransform;}
 191
 192     // Set mask, inside transform applies, outside zero, real pointer
 193     void SetMask(MaskPointer m){m_Mask=m;}
 194     MaskPointer GetMask(void){return m_Mask;}
 195     // itkSetConstObjectMacro( Mask, MaskType );
 196     // itkGetConstObjectMacro( Mask, MaskType );
 197
 198     /** Transform points by a BSpline deformable transformation. */
 199     OutputPointType  TransformPoint(const InputPointType  &point ) const;
 200
 201     // JV added for just the deformable part, without bulk
 202     OutputPointType  DeformablyTransformPoint(const InputPointType  &point ) const;
 203
 204     /** Parameter index array type. */
 205     typedef itk::Array<unsigned long> ParameterIndexArrayType;
 206
 207     /** Transform points by a BSpline deformable transformation.
 208      * On return, weights contains the interpolation weights used to compute the
 209      * deformation and indices of the x (zeroth) dimension coefficient parameters
 210      * in the support region used to compute the deformation.
 211      * Parameter indices for the i-th dimension can be obtained by adding
 212      * ( i * this->GetNumberOfParametersPerDimension() ) to the indices array.
 213      */
 214
 215     // JV not implemented
 216     //  virtual void TransformPoint( const InputPointType & inputPoint,
 217     //                           OutputPointType & outputPoint,
 218     //                           WeightsType & weights,
 219     //                           ParameterIndexArrayType & indices,
 220     //                           bool & inside ) const;
 221     //     virtual void DeformablyTransformPoint( const InputPointType & inputPoint,
 222     //                                     OutputPointType & outputPoint,
 223     //                                     WeightsType & weights,
 224     //                                     ParameterIndexArrayType & indices,
 225     //                                     bool & inside ) const;
 226     //     virtual void GetJacobian( const InputPointType & inputPoint,
 227     //                                WeightsType & weights,
 228     //                                ParameterIndexArrayType & indices
 229     //                                ) const;
 230
 231     /** Method to transform a vector -
 232      *  not applicable for this type of transform. */
 233     virtual OutputVectorType TransformVector(const InputVectorType &) const
 234     {
 235       itkExceptionMacro(<< "Method not applicable for deformable transform." );
 236       return OutputVectorType();
 237     }
 238
 239     /** Method to transform a vnl_vector -
 240      *  not applicable for this type of transform */
 241     virtual OutputVnlVectorType TransformVector(const InputVnlVectorType &) const
 242     {
 243       itkExceptionMacro(<< "Method not applicable for deformable transform. ");
 244       return OutputVnlVectorType();
 245     }
 246
 247     /** Method to transform a CovariantVector -
 248      *  not applicable for this type of transform */
 249     virtual OutputCovariantVectorType TransformCovariantVector(
 250                                                                const InputCovariantVectorType &) const
 251     {
 252       itkExceptionMacro(<< "Method not applicable for deformable transfrom. ");
 253       return OutputCovariantVectorType();
 254     }
 255
 256     /** Compute the Jacobian Matrix of the transformation at one point */
 257 #if ITK_VERSION_MAJOR >= 4
 258     virtual void ComputeJacobianWithRespectToParameters (const InputPointType &p, JacobianType &jacobian) const;
 259     virtual void ComputeJacobianWithRespectToPosition (const InputPointType &p, JacobianType &jacobian) const
 260     {
 261       itkExceptionMacro( "ComputeJacobianWithRespectToPosition not yet implemented for " << this->GetNameOfClass() );
 262     }
 263 #else
 264     virtual const JacobianType& GetJacobian(const InputPointType  &point ) const;
 265 #endif
 266
 267     /** Return the number of parameters that completely define the Transfom */
 268     virtual unsigned int GetNumberOfParameters(void) const;
 269
 270     /** Return the number of parameters per dimension */
 271     unsigned int GetNumberOfParametersPerDimension(void) const;
 272
 273     /** Return the region of the grid wholly within the support region */
 274     itkGetConstReferenceMacro( ValidRegion, RegionType );
 275
 276     /** Indicates that this transform is linear. That is, given two
 277      * points P and Q, and scalar coefficients a and b, then
 278      *
 279      *           T( a*P + b*Q ) = a * T(P) + b * T(Q)
 280      */
 281     virtual bool IsLinear() const { return false; }
 282
 283    //unsigned int GetNumberOfAffectedWeights() const;
 284
 285   protected:
 286     /** Print contents of an BSplineDeformableTransform. */
 287     void PrintSelf(std::ostream &os, itk::Indent indent) const;
 288
 289
 290     BSplineDeformableTransform();
 291     virtual ~BSplineDeformableTransform();
 292
 293     /** Wrap flat array into images of coefficients. */
 294     void WrapAsImages();
 295
 296     /** Convert an input point to a continuous index inside the BSpline grid */
 297     void TransformPointToContinuousIndex( const InputPointType & point, ContinuousIndexType & index ) const;
 298
 299   private:
 300     BSplineDeformableTransform(const Self&); //purposely not implemented
 301     void operator=(const Self&); //purposely not implemented
 302
 303     /** The bulk transform. */
 304     BulkTransformPointer  m_BulkTransform;
 305     MaskPointer m_Mask;
 306
 307     // JV added for BLUT interpolator
 308     SizeType m_SplineOrders;
 309     SizeType m_LUTSamplingFactors;
 310
 311     /** Variables defining the coefficient grid extend. */
 312     RegionType    m_GridRegion;
 313     SpacingType   m_GridSpacing;
 314     DirectionType m_GridDirection;
 315     OriginType    m_GridOrigin;
 316
 317     DirectionType m_PointToIndex;
 318     DirectionType m_IndexToPoint;
 319
 320     RegionType    m_ValidRegion;
 321
 322     /** Variables defining the interpolation support region. */
 323     SizeType m_Offset;
 324     itk::FixedArray<bool,InputDimension>    m_SplineOrderOdd;
 325     SizeType      m_SupportSize;
 326     IndexType     m_ValidRegionLast;
 327
 328     /** Array holding images wrapped from the flat parameters. */
 329     CoefficientImagePointer   m_WrappedImage;
 330
 331     /** Vector image representing the B-spline coefficients
 332      *  in each dimension. */
 333     CoefficientImagePointer   m_CoefficientImage;
 334
 335     /** Jacobian as OutputDimension number of images. */
 336     typedef typename JacobianType::ValueType JacobianValueType;
 337     typedef typename itk::Vector<JacobianValueType,OutputDimension> JacobianPixelType;
 338     typedef itk::Image<JacobianPixelType, OutputDimension> JacobianImageType;
 339     typename JacobianImageType::Pointer m_JacobianImage[OutputDimension];
 340     typedef itk::ImageRegionIterator<JacobianImageType> IteratorType;
 341
 342     /** Keep track of last support region used in computing the Jacobian
 343      * for fast resetting of Jacobian to zero.
 344      */
 345     //JV for J calculation
 346     mutable RegionType    m_SupportRegion;
 347     mutable IndexType     m_SupportIndex;
 348     mutable IndexType m_LastJacobianIndex;
 349     mutable IteratorType m_Iterator[OutputDimension];
 350     mutable JacobianPixelType m_ZeroVector;
 351     mutable ContinuousIndexType m_Index;
 352     mutable bool m_NeedResetJacobian;
 353
 354     /** Keep a pointer to the input parameters. */
 355     const ParametersType *  m_InputParametersPointer;
 356
 357     /** Internal parameters buffer. */
 358     ParametersType          m_InternalParametersBuffer;
 359
 360     //JV  the BLUT interpolator
 361     typename VectorInterpolatorType::Pointer m_VectorInterpolator;
 362
 363     /** Check if a continuous index is inside the valid region. */
 364     bool InsideValidRegion( const ContinuousIndexType& index ) const;
 365
 366     // VD Use external data container for JacobianImage
 367     unsigned SetJacobianImageData(JacobianPixelType * jacobianDataPointer, unsigned dim);
 368
 369     // VD Reset Jacobian
 370     void ResetJacobian() const;
 371
 372     // VD Add MultipleBSplineDeformableTransform as friend to facilitate wrapping
 373     friend class MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>;
 374 #if ITK_VERSION_MAJOR >= 4
 375     mutable JacobianType                            m_SharedDataBSplineJacobian;
 376 #endif
 377
 378   }; //class BSplineDeformableTransform
 379
 380
 381 }  // namespace itk
 382
 383 #if ITK_TEMPLATE_TXX
 384 # include "clitkBSplineDeformableTransform.txx"
 385 #endif
 386
 387
 388 #endif // __clitkBSplineDeformableTransform_h