registration/clitkShapedBLUTSpatioTemporalDeformableTransform.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 __clitkShapedBLUTSpatioTemporalDeformableTransform_h
  19 #define __clitkShapedBLUTSpatioTemporalDeformableTransform_h
  20 #include "clitkVectorBSplineResampleImageFunctionWithLUT.h"
  21 #include "clitkExtractImageFilter.h"
  22 #include "clitkLinearCombinationImageFilter.h"
  23
  24 //itk include
  25 #include "itkTransform.h"
  26 #include "itkImage.h"
  27 #include "itkImageRegion.h"
  28 #include "itkSpatialObject.h"
  29 #include "itkPasteImageFilter.h"
  30
  31 namespace clitk
  32 {
  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 ShapedBLUTSpatioTemporalDeformableTransform :
  39     public itk::Transform< TCoordRep, NInputDimensions, NOutputDimensions >
  40   {
  41   public:
  42
  43     //====================================================================
  44     // Typedefs
  45     //====================================================================
  46     typedef ShapedBLUTSpatioTemporalDeformableTransform                         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( ShapedBLUTSpatioTemporalDeformableTransform, Transform );
  56
  57     /** Dimension of the output space. */
  58     itkStaticConstMacro(OutputDimension, unsigned int, NOutputDimensions);
  59
  60     /** Dimension of the input space. */
  61     itkStaticConstMacro(InputDimension, unsigned int, NInputDimensions);
  62
  63     //JV the number of spatial dimensions
  64     itkStaticConstMacro(SpaceDimension, unsigned int, NInputDimensions-1);
  65
  66     /** Standard scalar type for this class. */
  67     typedef typename Superclass::ScalarType ScalarType;
  68
  69     /** Standard parameters container. */
  70     typedef typename Superclass::ParametersType ParametersType;
  71 #if ITK_VERSION_MAJOR >= 4
  72     typedef typename Superclass::NumberOfParametersType NumberOfParametersType;
  73 #endif
  74
  75     /** Standard Jacobian container. */
  76     typedef typename Superclass::JacobianType JacobianType;
  77
  78     /** Standard vector type for this class. */
  79     typedef itk::Vector<TCoordRep,
  80                         InputDimension> InputVectorType;
  81     typedef itk::Vector<TCoordRep,
  82                         OutputDimension> OutputVectorType;
  83
  84     /** Standard covariant vector type for this class. */
  85     typedef itk::CovariantVector<TCoordRep,
  86                                  InputDimension> InputCovariantVectorType;
  87     typedef itk::CovariantVector<TCoordRep,
  88                                  OutputDimension> OutputCovariantVectorType;
  89
  90     /** Standard vnl_vector type for this class. */
  91     typedef vnl_vector_fixed<TCoordRep,
  92                              InputDimension> InputVnlVectorType;
  93     typedef vnl_vector_fixed<TCoordRep,
  94                              OutputDimension> OutputVnlVectorType;
  95
  96     /** Standard coordinate point type for this class. */
  97     typedef itk::Point<TCoordRep,
  98                        InputDimension> InputPointType;
  99     typedef itk::Point<TCoordRep,
 100                        OutputDimension> OutputPointType;
 101
 102     //JV Parameters as images with OutputDimension number of components per Pixel
 103     typedef typename ParametersType::ValueType                          ParametersValueType;
 104     typedef typename itk::Vector<ParametersValueType, SpaceDimension>   PixelType;
 105     typedef itk::Image<PixelType, InputDimension>                       CoefficientImageType;
 106     typedef typename CoefficientImageType::Pointer                      CoefficientImagePointer;
 107
 108
 109     /** Typedefs for specifying the extend to the grid. */
 110     typedef itk::ImageRegion<InputDimension>                RegionType;
 111     typedef typename RegionType::IndexType                  IndexType;
 112     typedef typename RegionType::SizeType                   SizeType;
 113     typedef typename CoefficientImageType::SpacingType      SpacingType;
 114     typedef typename CoefficientImageType::DirectionType    DirectionType;
 115     typedef typename CoefficientImageType::PointType        OriginType;
 116     typedef itk::ContinuousIndex<TCoordRep, InputDimension> ContinuousIndexType;
 117
 118     //JV added for the BLUT interpolator
 119     typedef itk::Vector<TCoordRep, InputDimension> OutputSpacingType;
 120
 121     //JV m_VectorInterpolator
 122     typedef VectorBSplineResampleImageFunctionWithLUT
 123     <CoefficientImageType, TCoordRep> VectorInterpolatorType;
 124     typedef  typename  VectorInterpolatorType::CoefficientDataType CoefficientDataType;
 125     typedef  typename  VectorInterpolatorType::CoefficientDataType WeightsDataType;
 126
 127     /** Typedef of the bulk transform. */
 128     typedef itk::Transform<ScalarType, InputDimension, OutputDimension> BulkTransformType;
 129     typedef BulkTransformType*                  BulkTransformPointer;
 130
 131     /** Typedef of the mask */
 132     typedef itk::SpatialObject<  InputDimension >   MaskType;
 133     typedef MaskType*    MaskPointer;
 134
 135     //====================================================================
 136     // Set et Gets
 137     //====================================================================
 138     //JV  added for the BLUT interpolator
 139     void SetSplineOrder(const unsigned int &  splineOrder);
 140     void SetSplineOrders(const SizeType &  splineOrders);
 141     itkGetMacro( SplineOrders, SizeType );
 142     itkGetConstMacro( SplineOrders, SizeType );
 143     void SetLUTSamplingFactor(const int &  samplingFactor);
 144     void SetLUTSamplingFactors(const SizeType &  samplingFactors);
 145     itkGetMacro( LUTSamplingFactors, SizeType );
 146     itkGetConstMacro( LUTSamplingFactors,SizeType );
 147
 148     void SetParameters(const ParametersType & parameters);
 149
 150     void SetFixedParameters(const ParametersType & parameters);
 151
 152     void SetParametersByValue(const ParametersType & parameters);
 153
 154     void SetIdentity();
 155
 156     /** Get the Transformation Parameters. */
 157     virtual const ParametersType& GetParameters(void) const;
 158
 159     /** Get the Transformation Fixed Parameters. */
 160     virtual const ParametersType& GetFixedParameters(void) const;
 161
 162     // The coefficientImage
 163     virtual CoefficientImagePointer GetCoefficientImage()
 164     { return m_CoefficientImage; }
 165     virtual const CoefficientImagePointer  GetCoefficientImage() const
 166     { return m_CoefficientImage; }
 167     virtual void SetCoefficientImage(CoefficientImagePointer image);
 168
 169     // The padded coefficient image
 170     virtual CoefficientImagePointer GetPaddedCoefficientImage()
 171     { return m_PaddedCoefficientImage; }
 172     virtual const CoefficientImagePointer  GetPaddedCoefficientImage() const
 173     { return m_PaddedCoefficientImage; }
 174     // virtual void SetPaddedCoefficientImage(CoefficientImagePointer image);
 175
 176     /** This method specifies the region over which the grid resides. */
 177     virtual void SetGridRegion( const RegionType& region );
 178     itkGetMacro( GridRegion, RegionType );
 179     itkGetConstMacro( GridRegion, RegionType );
 180
 181     /** This method specifies the grid spacing or resolution. */
 182     virtual void SetGridSpacing( const SpacingType& spacing );
 183     itkGetMacro( GridSpacing, SpacingType );
 184     itkGetConstMacro( GridSpacing, SpacingType );
 185
 186     /** This method specifies the grid directions . */
 187     virtual void SetGridDirection( const DirectionType & spacing );
 188     itkGetMacro( GridDirection, DirectionType );
 189     itkGetConstMacro( GridDirection, DirectionType );
 190
 191     /** This method specifies the grid origin. */
 192     virtual void SetGridOrigin( const OriginType& origin );
 193     itkGetMacro( GridOrigin, OriginType );
 194     itkGetConstMacro( GridOrigin, OriginType );
 195
 196     // Set the bulk transform, real pointer
 197     // itkSetConstObjectMacro( BulkTransform, BulkTransformType );
 198     // itkGetConstObjectMacro( BulkTransform, BulkTransformType );
 199     void SetBulkTransform(BulkTransformPointer b){m_BulkTransform=b;}
 200     BulkTransformPointer GetBulkTransform(void) {return m_BulkTransform;}
 201
 202     //Set mask, inside transform applies, outside zero, real pointer
 203     void SetMask(MaskPointer m){m_Mask=m;}
 204
 205     // JV the shape
 206     itkSetMacro( TransformShape , unsigned int  );
 207     itkGetMacro( TransformShape , unsigned int  );
 208     itkGetConstMacro( TransformShape, unsigned int  );
 209
 210     /** Transform points by a BSpline deformable transformation. */
 211     OutputPointType  TransformPoint(const InputPointType  &point ) const;
 212
 213     // JV added for just the deformable part, without bulk
 214     OutputPointType  DeformablyTransformPoint(const InputPointType  &point ) const;
 215
 216     /** Parameter index array type. */
 217     typedef itk::Array<unsigned long> ParameterIndexArrayType;
 218
 219     /** Transform points by a BSpline deformable transformation.
 220      * On return, weights contains the interpolation weights used to compute the
 221      * deformation and indices of the x (zeroth) dimension coefficient parameters
 222      * in the support region used to compute the deformation.
 223      * Parameter indices for the i-th dimension can be obtained by adding
 224      * ( i * this->GetNumberOfParametersPerDimension() ) to the indices array.
 225      */
 226
 227     // JV not implemented
 228     //  virtual void TransformPoint( const InputPointType & inputPoint,
 229     //                           OutputPointType & outputPoint,
 230     //                           WeightsType & weights,
 231     //                           ParameterIndexArrayType & indices,
 232     //                           bool & inside ) const;
 233     //     virtual void DeformablyTransformPoint( const InputPointType & inputPoint,
 234     //                                     OutputPointType & outputPoint,
 235     //                                     WeightsType & weights,
 236     //                                     ParameterIndexArrayType & indices,
 237     //                                     bool & inside ) const;
 238     //     virtual void GetJacobian( const InputPointType & inputPoint,
 239     //                                WeightsType & weights,
 240     //                                ParameterIndexArrayType & indices
 241     //                                ) const;
 242
 243     /** Method to transform a vector -
 244      *  not applicable for this type of transform. */
 245     virtual OutputVectorType TransformVector(const InputVectorType &) const
 246     {
 247       itkExceptionMacro(<< "Method not applicable for deformable transform." );
 248       return OutputVectorType();
 249     }
 250
 251     /** Method to transform a vnl_vector -
 252      *  not applicable for this type of transform */
 253     virtual OutputVnlVectorType TransformVector(const InputVnlVectorType &) const
 254     {
 255       itkExceptionMacro(<< "Method not applicable for deformable transform. ");
 256       return OutputVnlVectorType();
 257     }
 258
 259     /** Method to transform a CovariantVector -
 260      *  not applicable for this type of transform */
 261     virtual OutputCovariantVectorType TransformCovariantVector(
 262                                                                const InputCovariantVectorType &) const
 263     {
 264       itkExceptionMacro(<< "Method not applicable for deformable transfrom. ");
 265       return OutputCovariantVectorType();
 266     }
 267
 268     /** Compute the Jacobian Matrix of the transformation at one point */
 269 #if ITK_VERSION_MAJOR >= 4
 270     virtual void ComputeJacobianWithRespectToParameters (const InputPointType &p, JacobianType &jacobian) const;
 271     virtual void ComputeJacobianWithRespectToPosition (const InputPointType &p, JacobianType &jacobian) const
 272     {
 273       itkExceptionMacro( "ComputeJacobianWithRespectToPosition not yet implemented for " << this->GetNameOfClass() );
 274     }
 275 #else
 276     virtual const JacobianType& GetJacobian(const InputPointType  &point ) const;
 277 #endif
 278
 279     /** Return the number of parameters that completely define the Transfom */
 280 #if ITK_VERSION_MAJOR >= 4
 281     virtual NumberOfParametersType GetNumberOfParameters(void) const;
 282 #else
 283     virtual unsigned int GetNumberOfParameters(void) const;
 284 #endif
 285
 286     //JV Return the padded number of parameters
 287     virtual unsigned int GetPaddedNumberOfParameters(void) const;
 288
 289     /** Return the number of parameters per dimension */
 290     unsigned int GetNumberOfParametersPerDimension(void) const;
 291
 292     /** Return the region of the grid wholly within the support region */
 293     itkGetConstReferenceMacro( ValidRegion, RegionType );
 294
 295     /** Indicates that this transform is linear. That is, given two
 296      * points P and Q, and scalar coefficients a and b, then
 297      *
 298      *           T( a*P + b*Q ) = a * T(P) + b * T(Q)
 299      */
 300     virtual bool IsLinear() const { return false; }
 301
 302    //unsigned int GetNumberOfAffectedWeights() const;
 303
 304   protected:
 305     /** Print contents of an BSplineSpatioTemporalDeformableTransform. */
 306     void PrintSelf(std::ostream &os, itk::Indent indent) const;
 307
 308
 309     ShapedBLUTSpatioTemporalDeformableTransform();
 310     virtual ~ShapedBLUTSpatioTemporalDeformableTransform();
 311
 312     /** Wrap flat array into images of coefficients. */
 313     void WrapAsImages();
 314
 315     // JV Pad/Extract the coefficient image
 316     void PadCoefficientImage(void);
 317     typename  CoefficientImageType::Pointer ExtractTemporalRow(const typename CoefficientImageType::Pointer& coefficientImage, unsigned int temporalIndex);
 318     //void ExtractCoefficientImage(void);
 319
 320     /** Wrap flat array into images of coefficients. */
 321     inline void WrapRegion(const RegionType& support,
 322                            RegionType& first,
 323                            RegionType& second,
 324                            RegionType& third,
 325                            std::vector<RegionType>& bc,
 326                            std::vector<double>& bcValues,
 327                            std::vector<RegionType>& bc2,
 328                            std::vector<double>& bc2Values,
 329                            std::vector<RegionType>& bc3,
 330                            std::vector<double>& bc3Values,
 331                            unsigned int& m_InitialOffset ) const;
 332
 333     /** Convert an input point to a continuous index inside the BSpline grid */
 334     void TransformPointToContinuousIndex( const InputPointType & point, ContinuousIndexType & index ) const;
 335
 336   private:
 337     ShapedBLUTSpatioTemporalDeformableTransform(const Self&); //purposely not implemented
 338     void operator=(const Self&); //purposely not implemented
 339
 340     /** The bulk transform. */
 341     BulkTransformPointer  m_BulkTransform;
 342     MaskPointer m_Mask;
 343
 344     // JV added for BLUT interpolator
 345     SizeType m_SplineOrders;
 346     //OutputSpacingType m_OutputSpacing;
 347     SizeType m_LUTSamplingFactors;
 348
 349     /** Variables defining the coefficient grid extend. */
 350     RegionType    m_GridRegion;
 351     SpacingType   m_GridSpacing;
 352     DirectionType m_GridDirection;
 353     OriginType    m_GridOrigin;
 354
 355     // JV additional variables for the padded region
 356     RegionType    m_PaddedGridRegion;
 357     OriginType    m_PaddedGridOrigin;
 358
 359
 360     DirectionType m_PointToIndex;
 361     DirectionType m_IndexToPoint;
 362
 363     RegionType    m_ValidRegion;
 364
 365     /** Variables defining the interpolation support region. */
 366     SizeType m_Offset;
 367     itk::FixedArray<bool,InputDimension>    m_SplineOrderOdd;
 368     SizeType      m_SupportSize;
 369     IndexType     m_ValidRegionLast;
 370
 371     /** Array holding images wrapped from the flat parameters. */
 372     CoefficientImagePointer   m_WrappedImage;
 373
 374     /** Vector image representing the B-spline coefficients
 375      *  in each dimension. */
 376     CoefficientImagePointer   m_CoefficientImage;
 377     CoefficientImagePointer   m_PaddedCoefficientImage;
 378
 379     /** Jacobian as OutputDimension number of images. */
 380     typedef typename JacobianType::ValueType JacobianValueType;
 381     typedef typename itk::Vector<JacobianValueType,SpaceDimension> JacobianPixelType;
 382     typedef itk::Image<JacobianPixelType, OutputDimension> JacobianImageType;
 383     typename JacobianImageType::Pointer m_JacobianImage[OutputDimension];
 384     typename JacobianImageType::Pointer m_PaddedJacobianImage[OutputDimension];
 385     typedef itk::ImageRegionIterator<JacobianImageType> IteratorType;
 386
 387     //JV for J calculation
 388     IndexType m_NullIndex;
 389     SizeType m_NullSize;
 390
 391     mutable RegionType    m_SupportRegion;
 392     mutable IndexType    m_SupportIndex;
 393     mutable RegionType    m_FirstRegion;
 394     mutable RegionType    m_SecondRegion;
 395     mutable RegionType    m_ThirdRegion;
 396     mutable unsigned int m_ThirdSize;
 397     mutable unsigned int m_InitialOffset;
 398
 399     mutable std::vector<IteratorType>  m_BCIterators[SpaceDimension];
 400     mutable std::vector<double>   m_BCValues;
 401     mutable std::vector<RegionType>  m_BCRegions;
 402     mutable unsigned int m_BCSize;
 403
 404     mutable std::vector<IteratorType>  m_BC2Iterators[SpaceDimension];
 405     mutable std::vector<double>   m_BC2Values;
 406     mutable std::vector<RegionType>  m_BC2Regions;
 407     mutable unsigned int m_BC2Size;
 408
 409     mutable std::vector<IteratorType>  m_BC3Iterators[SpaceDimension];
 410     mutable std::vector<double>   m_BC3Values;
 411     mutable std::vector<RegionType>  m_BC3Regions;
 412     mutable unsigned int m_BC3Size;
 413
 414     mutable RegionType m_BCRegion;
 415     mutable IteratorType m_FirstIterator[SpaceDimension];
 416     mutable IteratorType m_SecondIterator[SpaceDimension];
 417     mutable IteratorType m_ThirdIterator[SpaceDimension];
 418     mutable ContinuousIndexType m_Index;
 419
 420     //JV add a padded Jacobian matrix
 421     mutable JacobianType m_PaddedJacobian;
 422     mutable JacobianPixelType m_ZeroVector;
 423
 424     /** Keep track of last support region used in computing the Jacobian
 425      * for fast resetting of Jacobian to zero.
 426      */
 427     mutable IndexType m_LastJacobianIndex;
 428
 429     /** Keep a pointer to the input parameters. */
 430     const ParametersType *  m_InputParametersPointer;
 431
 432     /** Internal parameters buffer. */
 433     ParametersType          m_InternalParametersBuffer;
 434
 435     //JV  the BLUT interpolator
 436     typename VectorInterpolatorType::Pointer m_VectorInterpolator;
 437
 438     // the coefficients to apply the BC
 439     std::vector<std::vector<double> > m_Weights;
 440     std::vector<std::vector<double> > m_WeightRatio;
 441
 442     /** Check if a continuous index is inside the valid region. */
 443     bool InsideValidRegion( const ContinuousIndexType& index ) const;
 444
 445     // JV Shape
 446     unsigned int m_TransformShape;
 447
 448 #if ITK_VERSION_MAJOR >= 4
 449     mutable JacobianType                            m_SharedDataBSplineJacobian;
 450 #endif
 451
 452   }; //class ShapedBLUTSpatioTemporalDeformableTransform
 453
 454
 455 }  // namespace itk
 456
 457 #if ITK_TEMPLATE_TXX
 458 # include "clitkShapedBLUTSpatioTemporalDeformableTransform.txx"
 459 #endif
 460
 461
 462 #endif // __clitkShapedBLUTSpatioTemporalDeformableTransform_h