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