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 __clitkBSplineDeformableTransform_txx
19 #define __clitkBSplineDeformableTransform_tx
20 #include "clitkBSplineDeformableTransform.h"
23 #include "itkContinuousIndex.h"
24 #include "itkImageRegionIterator.h"
25 #include "itkImageRegionConstIterator.h"
26 #include "itkImageRegionConstIteratorWithIndex.h"
27 #include "itkIdentityTransform.h"
32 // Constructor with default arguments
33 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
34 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
35 #if ITK_VERSION_MAJOR >= 4
36 ::BSplineDeformableTransform():Superclass(0)
38 ::BSplineDeformableTransform():Superclass(OutputDimension,0)
43 //JV default spline order
44 for ( i=0;i<InputDimension; i++)
47 //JV default LUTSamplingfactor
48 for ( i=0;i<InputDimension; i++)
49 m_LUTSamplingFactors[i]=20;
51 for (i=0;i<InputDimension; i++)
52 m_SupportSize[i] = m_SplineOrders[i]+1;
54 //Instantiate a Vector Bspline Interpolator
55 m_VectorInterpolator =VectorInterpolatorType::New();
56 m_VectorInterpolator->SetLUTSamplingFactors(m_LUTSamplingFactors);
57 m_VectorInterpolator->SetSplineOrders(m_SplineOrders);
59 // Set Bulk transform to NULL (no bulk performed)
60 m_BulkTransform = NULL;
65 // Default grid size is zero
66 typename RegionType::SizeType size;
67 typename RegionType::IndexType index;
70 m_GridRegion.SetSize( size );
71 m_GridRegion.SetIndex( index );
73 m_GridOrigin.Fill( 0.0 ); // default origin is all zeros
74 m_GridSpacing.Fill( 1.0 ); // default spacing is all ones
75 m_GridDirection.SetIdentity(); // default spacing is all ones
77 m_InternalParametersBuffer = ParametersType(0);
78 // Make sure the parameters pointer is not NULL after construction.
79 m_InputParametersPointer = &m_InternalParametersBuffer;
81 // Initialize coeffient images
82 m_WrappedImage = CoefficientImageType::New();
83 m_WrappedImage->SetRegions( m_GridRegion );
84 m_WrappedImage->SetOrigin( m_GridOrigin.GetDataPointer() );
85 m_WrappedImage->SetSpacing( m_GridSpacing.GetDataPointer() );
86 m_WrappedImage->SetDirection( m_GridDirection );
87 m_CoefficientImage = NULL;
89 // Variables for computing interpolation
90 for (i=0; i <InputDimension;i++)
92 m_Offset[i] = m_SplineOrders[i] / 2;
93 m_SplineOrderOdd[i] = m_SplineOrders[i] % 2;
95 m_ValidRegion = m_GridRegion;
97 // Initialize jacobian images
98 for (i=0; i <OutputDimension;i++)
100 m_JacobianImage[i] = JacobianImageType::New();
101 m_JacobianImage[i]->SetRegions( m_GridRegion );
102 m_JacobianImage[i]->SetOrigin( m_GridOrigin.GetDataPointer() );
103 m_JacobianImage[i]->SetSpacing( m_GridSpacing.GetDataPointer() );
104 m_JacobianImage[i]->SetDirection( m_GridDirection );
107 /** Fixed Parameters store the following information:
112 //JV we add the splineOrders, LUTsamplingfactor, m_Mask and m_BulkTransform
118 The size of these is equal to the NInputDimensions
119 *********************************************************/
120 this->m_FixedParameters.SetSize ( NInputDimensions * (NInputDimensions + 5)+2 );
121 this->m_FixedParameters.Fill ( 0.0 );
122 for ( i=0; i<NInputDimensions; i++)
124 this->m_FixedParameters[2*NInputDimensions+i] = m_GridSpacing[i];
126 for (unsigned int di=0; di<NInputDimensions; di++)
128 for (unsigned int dj=0; dj<NInputDimensions; dj++)
130 this->m_FixedParameters[3*NInputDimensions+(di*NInputDimensions+dj)] = m_GridDirection[di][dj];
134 //JV add splineOrders
135 for ( i=0; i<NInputDimensions; i++)
137 this->m_FixedParameters[ ( (3+NInputDimensions)*NInputDimensions)+i] = (this->GetSplineOrders())[i];
140 //JV add LUTsamplingFactors
141 for ( i=0; i<NInputDimensions; i++)
143 this->m_FixedParameters[( (4+NInputDimensions)*NInputDimensions)+i ] = this->GetLUTSamplingFactors()[i];
146 // JV add the mask pointer
147 this->m_FixedParameters[( (5+NInputDimensions)*NInputDimensions)]=(double)((size_t)m_Mask);
149 // JV add the bulkTransform pointer
150 this->m_FixedParameters[( (5+NInputDimensions)*NInputDimensions) +1]=(double)((size_t)m_BulkTransform);
153 // Calculate the PointToIndex matrices
155 for( unsigned int i=0; i<OutputDimension; i++)
157 scale[i][i] = m_GridSpacing[i];
160 m_IndexToPoint = m_GridDirection * scale;
161 m_PointToIndex = m_IndexToPoint.GetInverse();
164 m_LastJacobianIndex = m_ValidRegion.GetIndex();
166 //JV initialize some variables for jacobian calculation
167 m_SupportRegion.SetSize(m_SupportSize);
168 m_SupportIndex.Fill(0);
169 m_SupportRegion.SetIndex(m_SupportIndex);
170 for ( i = 0; i < OutputDimension; i++ )
171 m_ZeroVector[i]=itk::NumericTraits<JacobianValueType>::Zero;
176 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
177 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
178 ::~BSplineDeformableTransform()
184 // JV set Spline Order
185 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
187 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
188 ::SetSplineOrder(const unsigned int & splineOrder)
190 SizeType splineOrders;
191 for (unsigned int i=0;i<InputDimension;i++)splineOrders[i]=splineOrder;
193 this->SetSplineOrders(splineOrders);
197 // JV set Spline Orders
198 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
200 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
201 ::SetSplineOrders(const SizeType & splineOrders)
203 if (m_SplineOrders!=splineOrders)
205 m_SplineOrders=splineOrders;
207 //update the interpolation function
208 m_VectorInterpolator->SetSplineOrders(m_SplineOrders);
210 //update the variables for computing interpolation
211 for (unsigned int i=0; i <InputDimension;i++)
213 m_SupportSize[i] = m_SplineOrders[i]+1;
214 m_Offset[i] = m_SplineOrders[i] / 2;
215 m_SplineOrderOdd[i] = m_SplineOrders[i] % 2;
222 // JV set sampling factor
223 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
225 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
226 ::SetLUTSamplingFactor( const int & samplingFactor)
228 SizeType samplingFactors;
229 for (unsigned int i=0; i<NInputDimensions; i++)
230 samplingFactors[i]=samplingFactor;
232 //update the interpolation function
233 this->SetLUTSamplingFactors(samplingFactors);
237 // JV set sampling factors
238 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
240 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
241 ::SetLUTSamplingFactors( const SizeType & samplingFactors)
243 if(m_LUTSamplingFactors!=samplingFactors)
245 for (unsigned int i=0; i<NInputDimensions; i++)
246 m_LUTSamplingFactors[i]=samplingFactors[i];
248 //update the interpolation function
249 m_VectorInterpolator->SetLUTSamplingFactors(m_LUTSamplingFactors);
256 // Get the number of parameters
257 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
259 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
260 ::GetNumberOfParameters(void) const
262 // The number of parameters equal OutputDimension * number of
263 // of pixels in the grid region.
264 return ( static_cast<unsigned int>( OutputDimension ) *
265 static_cast<unsigned int>( m_GridRegion.GetNumberOfPixels() ) );
270 // Get the number of parameters per dimension
271 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
273 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
274 ::GetNumberOfParametersPerDimension(void) const
276 // The number of parameters per dimension equal number of
277 // of pixels in the grid region.
278 return ( static_cast<unsigned int>( m_GridRegion.GetNumberOfPixels() ) );
283 // Set the grid region
284 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
286 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
287 ::SetGridRegion( const RegionType & region )
289 if ( m_GridRegion != region )
291 m_GridRegion = region;
293 // set regions for each coefficient and jacobian image
294 m_WrappedImage->SetRegions( m_GridRegion );
295 for (unsigned int j=0; j <OutputDimension;j++)
296 m_JacobianImage[j]->SetRegions( m_GridRegion );
298 // Set the valid region
299 // If the grid spans the interval [start, last].
300 // The valid interval for evaluation is [start+offset, last-offset]
301 // when spline order is even.
302 // The valid interval for evaluation is [start+offset, last-offset)
303 // when spline order is odd.
304 // Where offset = vcl_floor(spline / 2 ).
305 // Note that the last pixel is not included in the valid region
306 // with odd spline orders.
307 typename RegionType::SizeType size = m_GridRegion.GetSize();
308 typename RegionType::IndexType index = m_GridRegion.GetIndex();
309 for ( unsigned int j = 0; j < NInputDimensions; j++ )
312 static_cast< typename RegionType::IndexValueType >( m_Offset[j] );
314 static_cast< typename RegionType::SizeValueType> ( 2 * m_Offset[j] );
315 m_ValidRegionLast[j] = index[j] +
316 static_cast< typename RegionType::IndexValueType >( size[j] ) - 1;
318 m_ValidRegion.SetSize( size );
319 m_ValidRegion.SetIndex( index );
321 // If we are using the default parameters, update their size and set to identity.
322 // Input parameters point to internal buffer => using default parameters.
323 if (m_InputParametersPointer == &m_InternalParametersBuffer)
325 // Check if we need to resize the default parameter buffer.
326 if ( m_InternalParametersBuffer.GetSize() != this->GetNumberOfParameters() )
328 m_InternalParametersBuffer.SetSize( this->GetNumberOfParameters() );
329 // Fill with zeros for identity.
330 m_InternalParametersBuffer.Fill( 0 );
339 // Set the grid spacing
340 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
342 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
343 ::SetGridSpacing( const SpacingType & spacing )
345 if ( m_GridSpacing != spacing )
347 m_GridSpacing = spacing;
349 // Set spacing for each coefficient and jacobian image
350 m_WrappedImage->SetSpacing( m_GridSpacing.GetDataPointer() );
351 for (unsigned int j=0; j <OutputDimension; j++) m_JacobianImage[j]->SetSpacing( m_GridSpacing.GetDataPointer() );
355 for( unsigned int i=0; i<OutputDimension; i++)
357 scale[i][i] = m_GridSpacing[i];
360 m_IndexToPoint = m_GridDirection * scale;
361 m_PointToIndex = m_IndexToPoint.GetInverse();
368 // Set the grid direction
369 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
371 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
372 ::SetGridDirection( const DirectionType & direction )
374 if ( m_GridDirection != direction )
376 m_GridDirection = direction;
378 // Set direction for each coefficient and jacobian image
379 m_WrappedImage->SetDirection( m_GridDirection );
380 for (unsigned int j=0; j <OutputDimension; j++) m_JacobianImage[j]->SetDirection( m_GridDirection );
384 for( unsigned int i=0; i<OutputDimension; i++)
386 scale[i][i] = m_GridSpacing[i];
389 m_IndexToPoint = m_GridDirection * scale;
390 m_PointToIndex = m_IndexToPoint.GetInverse();
398 // Set the grid origin
399 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
401 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
402 ::SetGridOrigin( const OriginType& origin )
404 if ( m_GridOrigin != origin )
406 m_GridOrigin = origin;
408 // Set origin for each coefficient and jacobianimage
409 m_WrappedImage->SetOrigin( m_GridOrigin.GetDataPointer() );
410 for (unsigned int j=0; j <OutputDimension; j++) m_JacobianImage[j]->SetOrigin( m_GridOrigin.GetDataPointer() );
418 // Set the parameters
419 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
421 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
424 if( m_InputParametersPointer )
426 ParametersType * parameters =
427 const_cast<ParametersType *>( m_InputParametersPointer );
428 parameters->Fill( 0.0 );
433 itkExceptionMacro( << "Input parameters for the spline haven't been set ! "
434 << "Set them using the SetParameters or SetCoefficientImage method first." );
438 // Set the parameters
439 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
441 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
442 ::SetParameters( const ParametersType & parameters )
445 // Check if the number of parameters match the
446 // Expected number of parameters
447 if ( parameters.Size() != this->GetNumberOfParameters() )
449 itkExceptionMacro(<<"Mismatched between parameters size "
451 << " and region size "
452 << m_GridRegion.GetNumberOfPixels() );
455 // Clean up buffered parameters
456 m_InternalParametersBuffer = ParametersType( 0 );
458 // Keep a reference to the input parameters
459 m_InputParametersPointer = ¶meters;
461 // Wrap flat array as images of coefficients
462 this->WrapAsImages();
464 //Set input to vector interpolator
465 m_VectorInterpolator->SetInputImage(this->GetCoefficientImage());
467 // Modified is always called since we just have a pointer to the
468 // parameters and cannot know if the parameters have changed.
473 // Set the Fixed Parameters
474 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
476 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
477 ::SetFixedParameters( const ParametersType & parameters )
480 // JV number should be exact, no defaults for spacing
481 if ( parameters.Size() != NInputDimensions * (5 + NInputDimensions)+2 )
483 itkExceptionMacro(<< "Mismatched between parameters size "
485 << " and number of fixed parameters "
486 << NInputDimensions * (5 + NInputDimensions)+2 );
488 /*********************************************************
489 Fixed Parameters store the following information:
494 // JV we add the splineOrders, LUTsamplingfactor, mask pointer and bulktransform pointer
500 The size of these is equal to the NInputDimensions
501 *********************************************************/
503 /** Set the Grid Parameters */
505 for (unsigned int i=0; i<NInputDimensions; i++)
507 gridSize[i] = static_cast<int> (parameters[i]);
509 RegionType bsplineRegion;
510 bsplineRegion.SetSize( gridSize );
512 /** Set the Origin Parameters */
514 for (unsigned int i=0; i<NInputDimensions; i++)
516 origin[i] = parameters[NInputDimensions+i];
519 /** Set the Spacing Parameters */
521 for (unsigned int i=0; i<NInputDimensions; i++)
523 spacing[i] = parameters[2*NInputDimensions+i];
526 /** Set the Direction Parameters */
527 DirectionType direction;
528 for (unsigned int di=0; di<NInputDimensions; di++)
530 for (unsigned int dj=0; dj<NInputDimensions; dj++)
532 direction[di][dj] = parameters[3*NInputDimensions+(di*NInputDimensions+dj)];
536 //JV add the spline orders
537 SizeType splineOrders;
538 for (unsigned int i=0; i<NInputDimensions; i++)
540 splineOrders[i]=parameters[(3+NInputDimensions)*NInputDimensions+i];
543 //JV add the LUT sampling factor
544 SizeType samplingFactors;
545 for (unsigned int i=0; i<NInputDimensions; i++)
547 samplingFactors[i]=parameters[(4+NInputDimensions)*NInputDimensions+i];
550 //JV add the MaskPointer
551 m_Mask=((MaskPointer)( (size_t)parameters[(5+NInputDimensions)*NInputDimensions]));
553 //JV add the MaskPointer
554 m_BulkTransform=((BulkTransformPointer)( (size_t)parameters[(5+NInputDimensions)*NInputDimensions+1]));
557 this->SetGridSpacing( spacing );
558 this->SetGridDirection( direction );
559 this->SetGridOrigin( origin );
560 this->SetGridRegion( bsplineRegion );
562 //JV add the LUT parameters
563 this->SetSplineOrders( splineOrders );
564 this->SetLUTSamplingFactors( samplingFactors );
569 // Wrap flat parameters as images
570 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
572 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
575 //JV Wrap parameter array in vectorial image, changed parameter order: A1x A1y A1z, A2x ....
576 PixelType * dataPointer =reinterpret_cast<PixelType *>( const_cast<double *>(m_InputParametersPointer->data_block() )) ;
577 unsigned int numberOfPixels = m_GridRegion.GetNumberOfPixels();
579 m_WrappedImage->GetPixelContainer()->SetImportPointer( dataPointer,numberOfPixels);//InputDimension
580 m_CoefficientImage = m_WrappedImage;
582 //JV Wrap jacobian into OutputDimension X Vectorial images
583 #if ITK_VERSION_MAJOR >= 4
584 this->m_SharedDataBSplineJacobian.set_size( OutputDimension, this->GetNumberOfParameters() );
586 this->m_Jacobian.set_size( OutputDimension, this->GetNumberOfParameters() );
589 // Use memset to set the memory
590 #if ITK_VERSION_MAJOR >= 4
591 JacobianPixelType * jacobianDataPointer = reinterpret_cast<JacobianPixelType *>(this->m_SharedDataBSplineJacobian.data_block());
593 JacobianPixelType * jacobianDataPointer = reinterpret_cast<JacobianPixelType *>(this->m_Jacobian.data_block());
595 memset(jacobianDataPointer, 0, OutputDimension*numberOfPixels*sizeof(JacobianPixelType));
596 m_LastJacobianIndex = m_ValidRegion.GetIndex();
597 m_NeedResetJacobian = false;
599 for (unsigned int j=0; j<OutputDimension; j++)
601 m_JacobianImage[j]->GetPixelContainer()->
602 SetImportPointer( jacobianDataPointer, numberOfPixels );
603 jacobianDataPointer += numberOfPixels;
608 // Set the parameters by value
609 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
611 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
612 ::SetParametersByValue( const ParametersType & parameters )
615 // Check if the number of parameters match the
616 // Expected number of parameters
617 if ( parameters.Size() != this->GetNumberOfParameters() )
619 itkExceptionMacro(<<"Mismatched between parameters size "
621 << " and region size "
622 << m_GridRegion.GetNumberOfPixels() );
626 m_InternalParametersBuffer = parameters;
627 m_InputParametersPointer = &m_InternalParametersBuffer;
629 // Wrap flat array as images of coefficients
630 this->WrapAsImages();
632 // Modified is always called since we just have a pointer to the
633 // Parameters and cannot know if the parameters have changed.
638 // Get the parameters
639 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
641 typename BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
643 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
644 ::GetParameters( void ) const
646 /** NOTE: For efficiency, this class does not keep a copy of the parameters -
647 * it just keeps pointer to input parameters.
649 if (NULL == m_InputParametersPointer)
651 itkExceptionMacro( <<"Cannot GetParameters() because m_InputParametersPointer is NULL. Perhaps SetCoefficientImage() has been called causing the NULL pointer." );
654 return (*m_InputParametersPointer);
658 // Get the parameters
659 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
661 typename BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
663 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
664 ::GetFixedParameters( void ) const
666 RegionType resRegion = this->GetGridRegion( );
668 for (unsigned int i=0; i<NInputDimensions; i++)
670 this->m_FixedParameters[i] = (resRegion.GetSize())[i];
672 for (unsigned int i=0; i<NInputDimensions; i++)
674 this->m_FixedParameters[NInputDimensions+i] = (this->GetGridOrigin())[i];
676 for (unsigned int i=0; i<NInputDimensions; i++)
678 this->m_FixedParameters[2*NInputDimensions+i] = (this->GetGridSpacing())[i];
680 for (unsigned int di=0; di<NInputDimensions; di++)
682 for (unsigned int dj=0; dj<NInputDimensions; dj++)
684 this->m_FixedParameters[3*NInputDimensions+(di*NInputDimensions+dj)] = (this->GetGridDirection())[di][dj];
688 //JV add splineOrders
689 for (unsigned int i=0; i<NInputDimensions; i++)
691 this->m_FixedParameters[(3+NInputDimensions)*NInputDimensions+i] = (this->GetSplineOrders())[i];
694 //JV add LUTsamplingFactor
695 for (unsigned int i=0; i<NInputDimensions; i++)
697 this->m_FixedParameters[(4+NInputDimensions)*NInputDimensions+i] = (this->GetLUTSamplingFactors())[i];
701 this->m_FixedParameters[(5+NInputDimensions)*NInputDimensions]=(double)((size_t) m_Mask);
703 //JV add the bulktransform pointer
704 this->m_FixedParameters[(5+NInputDimensions)*NInputDimensions+1]=(double)((size_t) m_BulkTransform);
706 return (this->m_FixedParameters);
710 // Set the B-Spline coefficients using input images
711 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
713 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
714 ::SetCoefficientImage( CoefficientImagePointer image )
717 this->SetGridRegion( image->GetBufferedRegion() );
718 this->SetGridSpacing( image->GetSpacing() );
719 this->SetGridDirection( image->GetDirection() );
720 this->SetGridOrigin( image->GetOrigin() );
721 m_CoefficientImage = image;
723 // Update the interpolator
724 m_VectorInterpolator->SetInputImage(this->GetCoefficientImage());
726 // Clean up buffered parameters
727 m_InternalParametersBuffer = ParametersType( 0 );
728 m_InputParametersPointer = NULL;
734 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
736 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
737 ::PrintSelf(std::ostream &os, itk::Indent indent) const
740 this->Superclass::PrintSelf(os, indent);
742 os << indent << "GridRegion: " << m_GridRegion << std::endl;
743 os << indent << "GridOrigin: " << m_GridOrigin << std::endl;
744 os << indent << "GridSpacing: " << m_GridSpacing << std::endl;
745 os << indent << "GridDirection: " << m_GridDirection << std::endl;
746 os << indent << "IndexToPoint: " << m_IndexToPoint << std::endl;
747 os << indent << "PointToIndex: " << m_PointToIndex << std::endl;
749 os << indent << "CoefficientImage: [ ";
750 os << m_CoefficientImage.GetPointer() << " ]" << std::endl;
752 os << indent << "WrappedImage: [ ";
753 os << m_WrappedImage.GetPointer() << " ]" << std::endl;
755 os << indent << "InputParametersPointer: "
756 << m_InputParametersPointer << std::endl;
757 os << indent << "ValidRegion: " << m_ValidRegion << std::endl;
758 os << indent << "LastJacobianIndex: " << m_LastJacobianIndex << std::endl;
759 os << indent << "BulkTransform: ";
760 os << m_BulkTransform << std::endl;
762 if ( m_BulkTransform )
764 os << indent << "BulkTransformType: "
765 << m_BulkTransform->GetNameOfClass() << std::endl;
767 os << indent << "VectorBSplineInterpolator: ";
768 os << m_VectorInterpolator.GetPointer() << std::endl;
769 os << indent << "Mask: ";
770 os << m_Mask<< std::endl;
775 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
777 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
778 ::InsideValidRegion( const ContinuousIndexType& index ) const
782 if ( !m_ValidRegion.IsInside( index ) )
786 //JV verify for each input dimension
789 typedef typename ContinuousIndexType::ValueType ValueType;
790 for( unsigned int j = 0; j < InputDimension; j++ )
792 if (m_SplineOrderOdd[j])
794 if ( index[j] >= static_cast<ValueType>( m_ValidRegionLast[j] ) )
808 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
809 typename BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
811 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
812 ::TransformPoint(const InputPointType &inputPoint) const
815 InputPointType transformedPoint;
816 OutputPointType outputPoint;
819 if ( m_BulkTransform )
821 transformedPoint = m_BulkTransform->TransformPoint( inputPoint );
825 transformedPoint = inputPoint;
828 // Deformable transform
829 if ( m_CoefficientImage )
831 // Check if inside mask
832 if(m_Mask && !(m_Mask->IsInside(inputPoint) ) )
834 // Outside: no (deformable) displacement
835 return transformedPoint;
838 // Check if inside valid region
840 ContinuousIndexType index;
841 this->TransformPointToContinuousIndex( inputPoint, index );
842 inside = this->InsideValidRegion( index );
845 // Outside: no (deformable) displacement
846 return transformedPoint;
849 // Call the vector interpolator
850 itk::Vector<TCoordRep,OutputDimension> displacement=m_VectorInterpolator->EvaluateAtContinuousIndex(index);
852 // Return the results
853 outputPoint = transformedPoint+displacement;
859 itkWarningMacro( << "B-spline coefficients have not been set" );
860 outputPoint = transformedPoint;
868 //JV Deformably transform a point
869 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
870 typename BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
872 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
873 ::DeformablyTransformPoint(const InputPointType &inputPoint) const
875 OutputPointType outputPoint;
876 if ( m_CoefficientImage )
879 // Check if inside mask
880 if(m_Mask && !(m_Mask->IsInside(inputPoint) ) )
882 // Outside: no (deformable) displacement
888 ContinuousIndexType index;
889 this->TransformPointToContinuousIndex( inputPoint, index );
890 inside = this->InsideValidRegion( index );
894 //outside: no (deformable) displacement
896 //return outputPoint;
899 // Call the vector interpolator
900 itk::Vector<TCoordRep,OutputDimension> displacement=m_VectorInterpolator->EvaluateAtContinuousIndex(index);
902 // Return the results
903 outputPoint = inputPoint+displacement;
906 // No coefficients available
909 itkWarningMacro( << "B-spline coefficients have not been set" );
910 outputPoint = inputPoint;
917 // JV weights are identical as for transformpoint, could be done simultaneously in metric!!!!
918 // Compute the Jacobian in one position
919 #if ITK_VERSION_MAJOR >= 4
920 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
922 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
923 ::ComputeJacobianWithRespectToParameters(const InputPointType & point, JacobianType & jacobian) const
925 if (m_NeedResetJacobian)
928 //========================================================
929 // For each dimension, copy the weight to the support region
930 //========================================================
932 // Check if inside mask
933 if (m_Mask && !(m_Mask->IsInside(point) ) )
935 // Outside: no (deformable) displacement
936 jacobian = m_SharedDataBSplineJacobian;
941 this->TransformPointToContinuousIndex( point, m_Index );
943 // NOTE: if the support region does not lie totally within the grid
944 // we assume zero displacement and return the input point
945 if ( !this->InsideValidRegion( m_Index ) )
947 jacobian = m_SharedDataBSplineJacobian;
951 //Compute interpolation weights
952 const WeightsDataType *weights = NULL;
953 m_VectorInterpolator->EvaluateWeightsAtContinuousIndex( m_Index, &weights, m_LastJacobianIndex);
954 m_SupportRegion.SetIndex( m_LastJacobianIndex );
956 //Reset the iterators
958 for ( j = 0; j < OutputDimension; j++ )
959 m_Iterator[j] = IteratorType( m_JacobianImage[j], m_SupportRegion);
961 // For each dimension, copy the weight to the support region
962 while ( ! (m_Iterator[0]).IsAtEnd() )
964 //copy weight to jacobian image
965 for ( j = 0; j < OutputDimension; j++ )
967 m_ZeroVector[j]=*weights;
968 (m_Iterator[j]).Set( m_ZeroVector);
969 m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
972 // go to next coefficient in the support region
975 m_NeedResetJacobian = true;
978 jacobian = m_SharedDataBSplineJacobian;
982 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
984 typename BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
986 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
987 ::GetJacobian( const InputPointType & point ) const
989 // Can only compute Jacobian if parameters are set via
990 // SetParameters or SetParametersByValue
991 // if( m_InputParametersPointer == NULL )
993 // itkExceptionMacro( <<"Cannot compute Jacobian: parameters not set" );
996 if (m_NeedResetJacobian)
999 //========================================================
1000 // For each dimension, copy the weight to the support region
1001 //========================================================
1003 // Check if inside mask
1004 if(m_Mask && !(m_Mask->IsInside(point) ) )
1006 // Outside: no (deformable) displacement
1007 return this->m_Jacobian;
1011 this->TransformPointToContinuousIndex( point, m_Index );
1013 // NOTE: if the support region does not lie totally within the grid
1014 // we assume zero displacement and return the input point
1015 if ( !this->InsideValidRegion( m_Index ) )
1017 return this->m_Jacobian;
1020 //Compute interpolation weights
1021 const WeightsDataType *weights=NULL;
1022 m_VectorInterpolator->EvaluateWeightsAtContinuousIndex( m_Index, &weights, m_LastJacobianIndex);
1023 m_SupportRegion.SetIndex( m_LastJacobianIndex );
1025 //Reset the iterators
1027 for ( j = 0; j < OutputDimension; j++ )
1028 m_Iterator[j] = IteratorType( m_JacobianImage[j], m_SupportRegion);
1030 // For each dimension, copy the weight to the support region
1031 while ( ! (m_Iterator[0]).IsAtEnd() )
1033 //copy weight to jacobian image
1034 for ( j = 0; j < OutputDimension; j++ )
1036 m_ZeroVector[j]=*weights;
1037 (m_Iterator[j]).Set( m_ZeroVector);
1038 m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
1041 // go to next coefficient in the support region
1044 m_NeedResetJacobian = true;
1046 // Return the result
1047 return this->m_Jacobian;
1053 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
1055 BSplineDeformableTransform<TCoordRep, NInputDimensions,NOutputDimensions>
1056 ::TransformPointToContinuousIndex( const InputPointType & point, ContinuousIndexType & index ) const
1060 itk::Vector<double, OutputDimension> tvector;
1062 for ( j = 0; j < OutputDimension; j++ )
1064 tvector[j] = point[j] - this->m_GridOrigin[j];
1067 itk::Vector<double, OutputDimension> cvector;
1069 cvector = m_PointToIndex * tvector;
1071 for ( j = 0; j < OutputDimension; j++ )
1073 index[j] = static_cast< typename ContinuousIndexType::CoordRepType >( cvector[j] );
1077 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
1079 BSplineDeformableTransform<TCoordRep, NInputDimensions,NOutputDimensions>
1080 ::SetJacobianImageData(JacobianPixelType * jacobianDataPointer, unsigned dim)
1082 unsigned int numberOfPixels = m_GridRegion.GetNumberOfPixels();
1083 m_JacobianImage[dim]->GetPixelContainer()->SetImportPointer(jacobianDataPointer, numberOfPixels);
1084 return numberOfPixels;
1087 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
1089 BSplineDeformableTransform<TCoordRep, NInputDimensions,NOutputDimensions>
1090 ::ResetJacobian() const
1092 //========================================================
1093 // Zero all components of jacobian
1094 //========================================================
1095 // JV not thread safe (m_LastJacobianIndex), instantiate N transforms
1096 // NOTE: for efficiency, we only need to zero out the coefficients
1097 // that got fill last time this method was called.
1101 //Define the region for each jacobian image
1102 m_SupportRegion.SetIndex(m_LastJacobianIndex);
1104 //Initialize the iterators
1105 for (j = 0; j < OutputDimension; j++)
1106 m_Iterator[j] = IteratorType(m_JacobianImage[j], m_SupportRegion);
1108 //Set the previously-set to zero
1109 while (!(m_Iterator[0]).IsAtEnd())
1111 for (j = 0; j < OutputDimension; j++)
1113 m_Iterator[j].Set(m_ZeroVector);
1117 m_NeedResetJacobian = false;