1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
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
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.
13 It is distributed under dual licence
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"
25 #include "itkTransform.h"
27 #include "itkImageRegion.h"
28 #include "itkSpatialObject.h"
29 #include "itkPasteImageFilter.h"
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 >
43 //====================================================================
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;
51 /** New macro for creation of through the object factory.*/
54 /** Run-time type information (and related methods). */
55 itkTypeMacro( ShapedBLUTSpatioTemporalDeformableTransform, Transform );
57 /** Dimension of the output space. */
58 itkStaticConstMacro(OutputDimension, unsigned int, NOutputDimensions);
60 /** Dimension of the input space. */
61 itkStaticConstMacro(InputDimension, unsigned int, NInputDimensions);
63 //JV the number of spatial dimensions
64 itkStaticConstMacro(SpaceDimension, unsigned int, NInputDimensions-1);
66 /** Standard scalar type for this class. */
67 typedef typename Superclass::ScalarType ScalarType;
69 /** Standard parameters container. */
70 typedef typename Superclass::ParametersType ParametersType;
72 /** Standard Jacobian container. */
73 typedef typename Superclass::JacobianType JacobianType;
75 /** Standard vector type for this class. */
76 typedef itk::Vector<TCoordRep,
77 InputDimension> InputVectorType;
78 typedef itk::Vector<TCoordRep,
79 OutputDimension> OutputVectorType;
81 /** Standard covariant vector type for this class. */
82 typedef itk::CovariantVector<TCoordRep,
83 InputDimension> InputCovariantVectorType;
84 typedef itk::CovariantVector<TCoordRep,
85 OutputDimension> OutputCovariantVectorType;
87 /** Standard vnl_vector type for this class. */
88 typedef vnl_vector_fixed<TCoordRep,
89 InputDimension> InputVnlVectorType;
90 typedef vnl_vector_fixed<TCoordRep,
91 OutputDimension> OutputVnlVectorType;
93 /** Standard coordinate point type for this class. */
94 typedef itk::Point<TCoordRep,
95 InputDimension> InputPointType;
96 typedef itk::Point<TCoordRep,
97 OutputDimension> OutputPointType;
99 //JV Parameters as images with OutputDimension number of components per Pixel
100 typedef typename ParametersType::ValueType ParametersValueType;
101 typedef typename itk::Vector<ParametersValueType, SpaceDimension> PixelType;
102 typedef itk::Image<PixelType, InputDimension> CoefficientImageType;
103 typedef typename CoefficientImageType::Pointer CoefficientImagePointer;
106 /** Typedefs for specifying the extend to the grid. */
107 typedef itk::ImageRegion<InputDimension> RegionType;
108 typedef typename RegionType::IndexType IndexType;
109 typedef typename RegionType::SizeType SizeType;
110 typedef typename CoefficientImageType::SpacingType SpacingType;
111 typedef typename CoefficientImageType::DirectionType DirectionType;
112 typedef typename CoefficientImageType::PointType OriginType;
113 typedef itk::ContinuousIndex<TCoordRep, InputDimension> ContinuousIndexType;
115 //JV added for the BLUT interpolator
116 typedef itk::Vector<TCoordRep, InputDimension> OutputSpacingType;
118 //JV m_VectorInterpolator
119 typedef VectorBSplineResampleImageFunctionWithLUT
120 <CoefficientImageType, TCoordRep> VectorInterpolatorType;
121 typedef typename VectorInterpolatorType::CoefficientDataType CoefficientDataType;
122 typedef typename VectorInterpolatorType::CoefficientDataType WeightsDataType;
124 /** Typedef of the bulk transform. */
125 typedef itk::Transform<ScalarType, InputDimension, OutputDimension> BulkTransformType;
126 typedef BulkTransformType* BulkTransformPointer;
128 /** Typedef of the mask */
129 typedef itk::SpatialObject< InputDimension > MaskType;
130 typedef MaskType* MaskPointer;
132 //====================================================================
134 //====================================================================
135 //JV added for the BLUT interpolator
136 void SetSplineOrder(const unsigned int & splineOrder);
137 void SetSplineOrders(const SizeType & splineOrders);
138 itkGetMacro( SplineOrders, SizeType );
139 itkGetConstMacro( SplineOrders, SizeType );
140 void SetLUTSamplingFactor(const int & samplingFactor);
141 void SetLUTSamplingFactors(const SizeType & samplingFactors);
142 itkGetMacro( LUTSamplingFactors, SizeType );
143 itkGetConstMacro( LUTSamplingFactors,SizeType );
145 void SetParameters(const ParametersType & parameters);
147 void SetFixedParameters(const ParametersType & parameters);
149 void SetParametersByValue(const ParametersType & parameters);
153 /** Get the Transformation Parameters. */
154 virtual const ParametersType& GetParameters(void) const;
156 /** Get the Transformation Fixed Parameters. */
157 virtual const ParametersType& GetFixedParameters(void) const;
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);
166 // The padded coefficient image
167 virtual CoefficientImagePointer GetPaddedCoefficientImage()
168 { return m_PaddedCoefficientImage; }
169 virtual const CoefficientImagePointer GetPaddedCoefficientImage() const
170 { return m_PaddedCoefficientImage; }
171 // virtual void SetPaddedCoefficientImage(CoefficientImagePointer image);
173 /** This method specifies the region over which the grid resides. */
174 virtual void SetGridRegion( const RegionType& region );
175 itkGetMacro( GridRegion, RegionType );
176 itkGetConstMacro( GridRegion, RegionType );
178 /** This method specifies the grid spacing or resolution. */
179 virtual void SetGridSpacing( const SpacingType& spacing );
180 itkGetMacro( GridSpacing, SpacingType );
181 itkGetConstMacro( GridSpacing, SpacingType );
183 /** This method specifies the grid directions . */
184 virtual void SetGridDirection( const DirectionType & spacing );
185 itkGetMacro( GridDirection, DirectionType );
186 itkGetConstMacro( GridDirection, DirectionType );
188 /** This method specifies the grid origin. */
189 virtual void SetGridOrigin( const OriginType& origin );
190 itkGetMacro( GridOrigin, OriginType );
191 itkGetConstMacro( GridOrigin, OriginType );
193 // Set the bulk transform, real pointer
194 // itkSetConstObjectMacro( BulkTransform, BulkTransformType );
195 // itkGetConstObjectMacro( BulkTransform, BulkTransformType );
196 void SetBulkTransform(BulkTransformPointer b){m_BulkTransform=b;}
197 BulkTransformPointer GetBulkTransform(void) {return m_BulkTransform;}
199 //Set mask, inside transform applies, outside zero, real pointer
200 void SetMask(MaskPointer m){m_Mask=m;}
203 itkSetMacro( TransformShape , unsigned int );
204 itkGetMacro( TransformShape , unsigned int );
205 itkGetConstMacro( TransformShape, unsigned int );
207 /** Transform points by a BSpline deformable transformation. */
208 OutputPointType TransformPoint(const InputPointType &point ) const;
210 // JV added for just the deformable part, without bulk
211 OutputPointType DeformablyTransformPoint(const InputPointType &point ) const;
213 /** Parameter index array type. */
214 typedef itk::Array<unsigned long> ParameterIndexArrayType;
216 /** Transform points by a BSpline deformable transformation.
217 * On return, weights contains the interpolation weights used to compute the
218 * deformation and indices of the x (zeroth) dimension coefficient parameters
219 * in the support region used to compute the deformation.
220 * Parameter indices for the i-th dimension can be obtained by adding
221 * ( i * this->GetNumberOfParametersPerDimension() ) to the indices array.
224 // JV not implemented
225 // virtual void TransformPoint( const InputPointType & inputPoint,
226 // OutputPointType & outputPoint,
227 // WeightsType & weights,
228 // ParameterIndexArrayType & indices,
229 // bool & inside ) const;
230 // virtual void DeformablyTransformPoint( const InputPointType & inputPoint,
231 // OutputPointType & outputPoint,
232 // WeightsType & weights,
233 // ParameterIndexArrayType & indices,
234 // bool & inside ) const;
235 // virtual void GetJacobian( const InputPointType & inputPoint,
236 // WeightsType & weights,
237 // ParameterIndexArrayType & indices
240 /** Method to transform a vector -
241 * not applicable for this type of transform. */
242 virtual OutputVectorType TransformVector(const InputVectorType &) const
244 itkExceptionMacro(<< "Method not applicable for deformable transform." );
245 return OutputVectorType();
248 /** Method to transform a vnl_vector -
249 * not applicable for this type of transform */
250 virtual OutputVnlVectorType TransformVector(const InputVnlVectorType &) const
252 itkExceptionMacro(<< "Method not applicable for deformable transform. ");
253 return OutputVnlVectorType();
256 /** Method to transform a CovariantVector -
257 * not applicable for this type of transform */
258 virtual OutputCovariantVectorType TransformCovariantVector(
259 const InputCovariantVectorType &) const
261 itkExceptionMacro(<< "Method not applicable for deformable transfrom. ");
262 return OutputCovariantVectorType();
265 /** Compute the Jacobian Matrix of the transformation at one point */
266 #if ITK_VERSION_MAJOR >= 4
267 virtual void ComputeJacobianWithRespectToParameters (const InputPointType &p, JacobianType &jacobian) const;
268 virtual void ComputeJacobianWithRespectToPosition (const InputPointType &p, JacobianType &jacobian) const
270 itkExceptionMacro( "ComputeJacobianWithRespectToPosition not yet implemented for " << this->GetNameOfClass() );
273 virtual const JacobianType& GetJacobian(const InputPointType &point ) const;
276 /** Return the number of parameters that completely define the Transfom */
277 virtual unsigned int GetNumberOfParameters(void) const;
279 //JV Return the padded number of parameters
280 virtual unsigned int GetPaddedNumberOfParameters(void) const;
282 /** Return the number of parameters per dimension */
283 unsigned int GetNumberOfParametersPerDimension(void) const;
285 /** Return the region of the grid wholly within the support region */
286 itkGetConstReferenceMacro( ValidRegion, RegionType );
288 /** Indicates that this transform is linear. That is, given two
289 * points P and Q, and scalar coefficients a and b, then
291 * T( a*P + b*Q ) = a * T(P) + b * T(Q)
293 virtual bool IsLinear() const { return false; }
295 //unsigned int GetNumberOfAffectedWeights() const;
298 /** Print contents of an BSplineSpatioTemporalDeformableTransform. */
299 void PrintSelf(std::ostream &os, itk::Indent indent) const;
302 ShapedBLUTSpatioTemporalDeformableTransform();
303 virtual ~ShapedBLUTSpatioTemporalDeformableTransform();
305 /** Wrap flat array into images of coefficients. */
308 // JV Pad/Extract the coefficient image
309 void PadCoefficientImage(void);
310 typename CoefficientImageType::Pointer ExtractTemporalRow(const typename CoefficientImageType::Pointer& coefficientImage, unsigned int temporalIndex);
311 //void ExtractCoefficientImage(void);
313 /** Wrap flat array into images of coefficients. */
314 inline void WrapRegion(const RegionType& support,
318 std::vector<RegionType>& bc,
319 std::vector<double>& bcValues,
320 std::vector<RegionType>& bc2,
321 std::vector<double>& bc2Values,
322 std::vector<RegionType>& bc3,
323 std::vector<double>& bc3Values,
324 unsigned int& m_InitialOffset ) const;
326 /** Convert an input point to a continuous index inside the BSpline grid */
327 void TransformPointToContinuousIndex( const InputPointType & point, ContinuousIndexType & index ) const;
330 ShapedBLUTSpatioTemporalDeformableTransform(const Self&); //purposely not implemented
331 void operator=(const Self&); //purposely not implemented
333 /** The bulk transform. */
334 BulkTransformPointer m_BulkTransform;
337 // JV added for BLUT interpolator
338 SizeType m_SplineOrders;
339 //OutputSpacingType m_OutputSpacing;
340 SizeType m_LUTSamplingFactors;
342 /** Variables defining the coefficient grid extend. */
343 RegionType m_GridRegion;
344 SpacingType m_GridSpacing;
345 DirectionType m_GridDirection;
346 OriginType m_GridOrigin;
348 // JV additional variables for the padded region
349 RegionType m_PaddedGridRegion;
350 OriginType m_PaddedGridOrigin;
353 DirectionType m_PointToIndex;
354 DirectionType m_IndexToPoint;
356 RegionType m_ValidRegion;
358 /** Variables defining the interpolation support region. */
360 itk::FixedArray<bool,InputDimension> m_SplineOrderOdd;
361 SizeType m_SupportSize;
362 IndexType m_ValidRegionLast;
364 /** Array holding images wrapped from the flat parameters. */
365 CoefficientImagePointer m_WrappedImage;
367 /** Vector image representing the B-spline coefficients
368 * in each dimension. */
369 CoefficientImagePointer m_CoefficientImage;
370 CoefficientImagePointer m_PaddedCoefficientImage;
372 /** Jacobian as OutputDimension number of images. */
373 typedef typename JacobianType::ValueType JacobianValueType;
374 typedef typename itk::Vector<JacobianValueType,SpaceDimension> JacobianPixelType;
375 typedef itk::Image<JacobianPixelType, OutputDimension> JacobianImageType;
376 typename JacobianImageType::Pointer m_JacobianImage[OutputDimension];
377 typename JacobianImageType::Pointer m_PaddedJacobianImage[OutputDimension];
378 typedef itk::ImageRegionIterator<JacobianImageType> IteratorType;
380 //JV for J calculation
381 IndexType m_NullIndex;
384 mutable RegionType m_SupportRegion;
385 mutable IndexType m_SupportIndex;
386 mutable RegionType m_FirstRegion;
387 mutable RegionType m_SecondRegion;
388 mutable RegionType m_ThirdRegion;
389 mutable unsigned int m_ThirdSize;
390 mutable unsigned int m_InitialOffset;
392 mutable std::vector<IteratorType> m_BCIterators[SpaceDimension];
393 mutable std::vector<double> m_BCValues;
394 mutable std::vector<RegionType> m_BCRegions;
395 mutable unsigned int m_BCSize;
397 mutable std::vector<IteratorType> m_BC2Iterators[SpaceDimension];
398 mutable std::vector<double> m_BC2Values;
399 mutable std::vector<RegionType> m_BC2Regions;
400 mutable unsigned int m_BC2Size;
402 mutable std::vector<IteratorType> m_BC3Iterators[SpaceDimension];
403 mutable std::vector<double> m_BC3Values;
404 mutable std::vector<RegionType> m_BC3Regions;
405 mutable unsigned int m_BC3Size;
407 mutable RegionType m_BCRegion;
408 mutable IteratorType m_FirstIterator[SpaceDimension];
409 mutable IteratorType m_SecondIterator[SpaceDimension];
410 mutable IteratorType m_ThirdIterator[SpaceDimension];
411 mutable ContinuousIndexType m_Index;
413 //JV add a padded Jacobian matrix
414 mutable JacobianType m_PaddedJacobian;
415 mutable JacobianPixelType m_ZeroVector;
417 /** Keep track of last support region used in computing the Jacobian
418 * for fast resetting of Jacobian to zero.
420 mutable IndexType m_LastJacobianIndex;
422 /** Keep a pointer to the input parameters. */
423 const ParametersType * m_InputParametersPointer;
425 /** Internal parameters buffer. */
426 ParametersType m_InternalParametersBuffer;
428 //JV the BLUT interpolator
429 typename VectorInterpolatorType::Pointer m_VectorInterpolator;
431 // the coefficients to apply the BC
432 std::vector<std::vector<double> > m_Weights;
433 std::vector<std::vector<double> > m_WeightRatio;
435 /** Check if a continuous index is inside the valid region. */
436 bool InsideValidRegion( const ContinuousIndexType& index ) const;
439 unsigned int m_TransformShape;
441 #if ITK_VERSION_MAJOR >= 4
442 mutable JacobianType m_SharedDataBSplineJacobian;
445 }; //class ShapedBLUTSpatioTemporalDeformableTransform
451 # include "clitkShapedBLUTSpatioTemporalDeformableTransform.txx"
455 #endif // __clitkShapedBLUTSpatioTemporalDeformableTransform_h