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://oncora1.lyon.fnclcc.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 ::BSplineDeformableTransform():Superclass(OutputDimension,0)
39 //JV default spline order
40 for ( i=0;i<InputDimension; i++)
43 //JV default LUTSamplingfactor
44 for ( i=0;i<InputDimension; i++)
45 m_LUTSamplingFactors[i]=20;
47 for (i=0;i<InputDimension; i++)
48 m_SupportSize[i] = m_SplineOrders[i]+1;
50 //Instantiate a Vector Bspline Interpolator
51 m_VectorInterpolator =VectorInterpolatorType::New();
52 m_VectorInterpolator->SetLUTSamplingFactors(m_LUTSamplingFactors);
53 m_VectorInterpolator->SetSplineOrders(m_SplineOrders);
55 // Set Bulk transform to NULL (no bulk performed)
56 m_BulkTransform = NULL;
61 // Default grid size is zero
62 typename RegionType::SizeType size;
63 typename RegionType::IndexType index;
66 m_GridRegion.SetSize( size );
67 m_GridRegion.SetIndex( index );
69 m_GridOrigin.Fill( 0.0 ); // default origin is all zeros
70 m_GridSpacing.Fill( 1.0 ); // default spacing is all ones
71 m_GridDirection.SetIdentity(); // default spacing is all ones
73 m_InternalParametersBuffer = ParametersType(0);
74 // Make sure the parameters pointer is not NULL after construction.
75 m_InputParametersPointer = &m_InternalParametersBuffer;
77 // Initialize coeffient images
78 m_WrappedImage = CoefficientImageType::New();
79 m_WrappedImage->SetRegions( m_GridRegion );
80 m_WrappedImage->SetOrigin( m_GridOrigin.GetDataPointer() );
81 m_WrappedImage->SetSpacing( m_GridSpacing.GetDataPointer() );
82 m_WrappedImage->SetDirection( m_GridDirection );
83 m_CoefficientImage = NULL;
85 // Variables for computing interpolation
86 for (i=0; i <InputDimension;i++)
88 m_Offset[i] = m_SplineOrders[i] / 2;
89 if ( m_SplineOrders[i] % 2 )
91 m_SplineOrderOdd[i] = true;
95 m_SplineOrderOdd[i] = false;
98 m_ValidRegion = m_GridRegion;
100 // Initialize jacobian images
101 for (i=0; i <OutputDimension;i++)
103 m_JacobianImage[i] = JacobianImageType::New();
104 m_JacobianImage[i]->SetRegions( m_GridRegion );
105 m_JacobianImage[i]->SetOrigin( m_GridOrigin.GetDataPointer() );
106 m_JacobianImage[i]->SetSpacing( m_GridSpacing.GetDataPointer() );
107 m_JacobianImage[i]->SetDirection( m_GridDirection );
110 /** Fixed Parameters store the following information:
115 //JV we add the splineOrders, LUTsamplingfactor, m_Mask and m_BulkTransform
121 The size of these is equal to the NInputDimensions
122 *********************************************************/
123 this->m_FixedParameters.SetSize ( NInputDimensions * (NInputDimensions + 5)+2 );
124 this->m_FixedParameters.Fill ( 0.0 );
125 for ( i=0; i<NInputDimensions; i++)
127 this->m_FixedParameters[2*NInputDimensions+i] = m_GridSpacing[i];
129 for (unsigned int di=0; di<NInputDimensions; di++)
131 for (unsigned int dj=0; dj<NInputDimensions; dj++)
133 this->m_FixedParameters[3*NInputDimensions+(di*NInputDimensions+dj)] = m_GridDirection[di][dj];
137 //JV add splineOrders
138 for ( i=0; i<NInputDimensions; i++)
140 this->m_FixedParameters[ ( (3+NInputDimensions)*NInputDimensions)+i] = (this->GetSplineOrders())[i];
143 //JV add LUTsamplingFactors
144 for ( i=0; i<NInputDimensions; i++)
146 this->m_FixedParameters[( (4+NInputDimensions)*NInputDimensions)+i ] = this->GetLUTSamplingFactors()[i];
149 // JV add the mask pointer
150 this->m_FixedParameters[( (5+NInputDimensions)*NInputDimensions)]=(double)((size_t)m_Mask);
152 // JV add the bulkTransform pointer
153 this->m_FixedParameters[( (5+NInputDimensions)*NInputDimensions) +1]=(double)((size_t)m_BulkTransform);
156 // Calculate the PointToIndex matrices
158 for( unsigned int i=0; i<OutputDimension; i++)
160 scale[i][i] = m_GridSpacing[i];
163 m_IndexToPoint = m_GridDirection * scale;
164 m_PointToIndex = m_IndexToPoint.GetInverse();
167 m_LastJacobianIndex = m_ValidRegion.GetIndex();
169 //JV initialize some variables for jacobian calculation
170 m_SupportRegion.SetSize(m_SupportSize);
171 m_SupportIndex.Fill(0);
172 m_SupportRegion.SetIndex(m_SupportIndex);
173 for ( i = 0; i < OutputDimension; i++ )
174 m_ZeroVector[i]=itk::NumericTraits<JacobianValueType>::Zero;
180 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
181 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
182 ::~BSplineDeformableTransform()
188 // JV set Spline Order
189 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
191 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
192 ::SetSplineOrder(const unsigned int & splineOrder)
194 SizeType splineOrders;
195 for (unsigned int i=0;i<InputDimension;i++)splineOrders[i]=splineOrder;
197 this->SetSplineOrders(splineOrders);
201 // JV set Spline Orders
202 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
204 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
205 ::SetSplineOrders(const SizeType & splineOrders)
207 if(m_SplineOrders!=splineOrders)
209 m_SplineOrders=splineOrders;
211 //update the interpolation function
212 m_VectorInterpolator->SetSplineOrders(m_SplineOrders);
214 //update the varaibles for computing interpolation
215 for (unsigned int i=0; i <InputDimension;i++)
217 m_SupportSize[i] = m_SplineOrders[i]+1;
218 m_Offset[i] = m_SplineOrders[i] / 2;
220 if ( m_SplineOrders[i] % 2 )
222 m_SplineOrderOdd[i] = true;
226 m_SplineOrderOdd[i] = false;
234 // JV set sampling factor
235 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
237 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
238 ::SetLUTSamplingFactor( const int & samplingFactor)
240 SizeType samplingFactors;
241 for (unsigned int i=0; i<NInputDimensions; i++)
242 samplingFactors[i]=samplingFactor;
244 //update the interpolation function
245 this->SetLUTSamplingFactors(samplingFactors);
249 // JV set sampling factors
250 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
252 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
253 ::SetLUTSamplingFactors( const SizeType & samplingFactors)
255 if(m_LUTSamplingFactors!=samplingFactors)
257 for (unsigned int i=0; i<NInputDimensions; i++)
258 m_LUTSamplingFactors[i]=samplingFactors[i];
260 //update the interpolation function
261 m_VectorInterpolator->SetLUTSamplingFactors(m_LUTSamplingFactors);
268 // Get the number of parameters
269 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
271 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
272 ::GetNumberOfParameters(void) const
275 // The number of parameters equal OutputDimension * number of
276 // of pixels in the grid region.
277 return ( static_cast<unsigned int>( OutputDimension ) *
278 static_cast<unsigned int>( m_GridRegion.GetNumberOfPixels() ) );
283 // Get the number of parameters per dimension
284 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
286 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
287 ::GetNumberOfParametersPerDimension(void) const
289 // The number of parameters per dimension equal number of
290 // of pixels in the grid region.
291 return ( static_cast<unsigned int>( m_GridRegion.GetNumberOfPixels() ) );
296 // Set the grid region
297 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
299 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
300 ::SetGridRegion( const RegionType & region )
302 if ( m_GridRegion != region )
304 m_GridRegion = region;
306 // set regions for each coefficient and jacobian image
307 m_WrappedImage->SetRegions( m_GridRegion );
308 for (unsigned int j=0; j <OutputDimension;j++)
309 m_JacobianImage[j]->SetRegions( m_GridRegion );
311 // Set the valid region
312 // If the grid spans the interval [start, last].
313 // The valid interval for evaluation is [start+offset, last-offset]
314 // when spline order is even.
315 // The valid interval for evaluation is [start+offset, last-offset)
316 // when spline order is odd.
317 // Where offset = vcl_floor(spline / 2 ).
318 // Note that the last pixel is not included in the valid region
319 // with odd spline orders.
320 typename RegionType::SizeType size = m_GridRegion.GetSize();
321 typename RegionType::IndexType index = m_GridRegion.GetIndex();
322 for ( unsigned int j = 0; j < NInputDimensions; j++ )
325 static_cast< typename RegionType::IndexValueType >( m_Offset[j] );
327 static_cast< typename RegionType::SizeValueType> ( 2 * m_Offset[j] );
328 m_ValidRegionLast[j] = index[j] +
329 static_cast< typename RegionType::IndexValueType >( size[j] ) - 1;
331 m_ValidRegion.SetSize( size );
332 m_ValidRegion.SetIndex( index );
334 // If we are using the default parameters, update their size and set to identity.
335 // Input parameters point to internal buffer => using default parameters.
336 if (m_InputParametersPointer == &m_InternalParametersBuffer)
338 // Check if we need to resize the default parameter buffer.
339 if ( m_InternalParametersBuffer.GetSize() != this->GetNumberOfParameters() )
341 m_InternalParametersBuffer.SetSize( this->GetNumberOfParameters() );
342 // Fill with zeros for identity.
343 m_InternalParametersBuffer.Fill( 0 );
352 // Set the grid spacing
353 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
355 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
356 ::SetGridSpacing( const SpacingType & spacing )
358 if ( m_GridSpacing != spacing )
360 m_GridSpacing = spacing;
362 // Set spacing for each coefficient and jacobian image
363 m_WrappedImage->SetSpacing( m_GridSpacing.GetDataPointer() );
364 for (unsigned int j=0; j <OutputDimension; j++) m_JacobianImage[j]->SetSpacing( m_GridSpacing.GetDataPointer() );
368 for( unsigned int i=0; i<OutputDimension; i++)
370 scale[i][i] = m_GridSpacing[i];
373 m_IndexToPoint = m_GridDirection * scale;
374 m_PointToIndex = m_IndexToPoint.GetInverse();
381 // Set the grid direction
382 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
384 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
385 ::SetGridDirection( const DirectionType & direction )
387 if ( m_GridDirection != direction )
389 m_GridDirection = direction;
391 // Set direction for each coefficient and jacobian image
392 m_WrappedImage->SetDirection( m_GridDirection );
393 for (unsigned int j=0; j <OutputDimension; j++) m_JacobianImage[j]->SetDirection( m_GridDirection );
397 for( unsigned int i=0; i<OutputDimension; i++)
399 scale[i][i] = m_GridSpacing[i];
402 m_IndexToPoint = m_GridDirection * scale;
403 m_PointToIndex = m_IndexToPoint.GetInverse();
411 // Set the grid origin
412 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
414 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
415 ::SetGridOrigin( const OriginType& origin )
417 if ( m_GridOrigin != origin )
419 m_GridOrigin = origin;
421 // Set origin for each coefficient and jacobianimage
422 m_WrappedImage->SetOrigin( m_GridOrigin.GetDataPointer() );
423 for (unsigned int j=0; j <OutputDimension; j++) m_JacobianImage[j]->SetOrigin( m_GridOrigin.GetDataPointer() );
431 // Set the parameters
432 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
434 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
437 if( m_InputParametersPointer )
439 ParametersType * parameters =
440 const_cast<ParametersType *>( m_InputParametersPointer );
441 parameters->Fill( 0.0 );
446 itkExceptionMacro( << "Input parameters for the spline haven't been set ! "
447 << "Set them using the SetParameters or SetCoefficientImage method first." );
452 // Set the parameters
453 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
455 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
456 ::SetParameters( const ParametersType & parameters )
459 // Check if the number of parameters match the
460 // Expected number of parameters
461 if ( parameters.Size() != this->GetNumberOfParameters() )
463 itkExceptionMacro(<<"Mismatched between parameters size "
465 << " and region size "
466 << m_GridRegion.GetNumberOfPixels() );
469 // Clean up buffered parameters
470 m_InternalParametersBuffer = ParametersType( 0 );
472 // Keep a reference to the input parameters
473 m_InputParametersPointer = ¶meters;
475 // Wrap flat array as images of coefficients
476 this->WrapAsImages();
478 //Set input to vector interpolator
479 m_VectorInterpolator->SetInputImage(this->GetCoefficientImage());
481 // Modified is always called since we just have a pointer to the
482 // parameters and cannot know if the parameters have changed.
487 // Set the Fixed Parameters
488 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
490 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
491 ::SetFixedParameters( const ParametersType & parameters )
494 // JV number should be exact, no defaults for spacing
495 if ( parameters.Size() != NInputDimensions * (5 + NInputDimensions)+2 )
497 itkExceptionMacro(<< "Mismatched between parameters size "
499 << " and number of fixed parameters "
500 << NInputDimensions * (5 + NInputDimensions)+2 );
502 /*********************************************************
503 Fixed Parameters store the following information:
508 // JV we add the splineOrders, LUTsamplingfactor, mask pointer and bulktransform pointer
514 The size of these is equal to the NInputDimensions
515 *********************************************************/
517 /** Set the Grid Parameters */
519 for (unsigned int i=0; i<NInputDimensions; i++)
521 gridSize[i] = static_cast<int> (parameters[i]);
523 RegionType bsplineRegion;
524 bsplineRegion.SetSize( gridSize );
526 /** Set the Origin Parameters */
528 for (unsigned int i=0; i<NInputDimensions; i++)
530 origin[i] = parameters[NInputDimensions+i];
533 /** Set the Spacing Parameters */
535 for (unsigned int i=0; i<NInputDimensions; i++)
537 spacing[i] = parameters[2*NInputDimensions+i];
540 /** Set the Direction Parameters */
541 DirectionType direction;
542 for (unsigned int di=0; di<NInputDimensions; di++)
544 for (unsigned int dj=0; dj<NInputDimensions; dj++)
546 direction[di][dj] = parameters[3*NInputDimensions+(di*NInputDimensions+dj)];
550 //JV add the spline orders
551 SizeType splineOrders;
552 for (unsigned int i=0; i<NInputDimensions; i++)
554 splineOrders[i]=parameters[(3+NInputDimensions)*NInputDimensions+i];
557 //JV add the LUT sampling factor
558 SizeType samplingFactors;
559 for (unsigned int i=0; i<NInputDimensions; i++)
561 samplingFactors[i]=parameters[(4+NInputDimensions)*NInputDimensions+i];
564 //JV add the MaskPointer
565 m_Mask=((MaskPointer)( (size_t)parameters[(5+NInputDimensions)*NInputDimensions]));
567 //JV add the MaskPointer
568 m_BulkTransform=((BulkTransformPointer)( (size_t)parameters[(5+NInputDimensions)*NInputDimensions+1]));
571 this->SetGridSpacing( spacing );
572 this->SetGridDirection( direction );
573 this->SetGridOrigin( origin );
574 this->SetGridRegion( bsplineRegion );
576 //JV add the LUT parameters
577 this->SetSplineOrders( splineOrders );
578 this->SetLUTSamplingFactors( samplingFactors );
583 // Wrap flat parameters as images
584 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
586 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
589 //JV Wrap parameter array in vectorial image, changed parameter order: A1x A1y A1z, A2x ....
590 PixelType * dataPointer =reinterpret_cast<PixelType *>( const_cast<double *>(m_InputParametersPointer->data_block() )) ;
591 unsigned int numberOfPixels = m_GridRegion.GetNumberOfPixels();
593 m_WrappedImage->GetPixelContainer()->SetImportPointer( dataPointer,numberOfPixels);//InputDimension
594 m_CoefficientImage = m_WrappedImage;
596 //JV Wrap jacobian into OutputDimension X Vectorial images
597 this->m_Jacobian.set_size( OutputDimension, this->GetNumberOfParameters() );
599 // Use memset to set the memory
600 JacobianPixelType * jacobianDataPointer = reinterpret_cast<JacobianPixelType *>(this->m_Jacobian.data_block());
601 memset(jacobianDataPointer, 0, OutputDimension*numberOfPixels*sizeof(JacobianPixelType));
602 m_LastJacobianIndex = m_ValidRegion.GetIndex();
604 for (unsigned int j=0; j<OutputDimension; j++)
606 m_JacobianImage[j]->GetPixelContainer()->
607 SetImportPointer( jacobianDataPointer, numberOfPixels );
608 jacobianDataPointer += numberOfPixels;
613 // Set the parameters by value
614 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
616 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
617 ::SetParametersByValue( const ParametersType & parameters )
620 // Check if the number of parameters match the
621 // Expected number of parameters
622 if ( parameters.Size() != this->GetNumberOfParameters() )
624 itkExceptionMacro(<<"Mismatched between parameters size "
626 << " and region size "
627 << m_GridRegion.GetNumberOfPixels() );
631 m_InternalParametersBuffer = parameters;
632 m_InputParametersPointer = &m_InternalParametersBuffer;
634 // Wrap flat array as images of coefficients
635 this->WrapAsImages();
637 // Modified is always called since we just have a pointer to the
638 // Parameters and cannot know if the parameters have changed.
643 // Get the parameters
644 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
646 typename BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
648 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
649 ::GetParameters( void ) const
651 /** NOTE: For efficiency, this class does not keep a copy of the parameters -
652 * it just keeps pointer to input parameters.
654 if (NULL == m_InputParametersPointer)
656 itkExceptionMacro( <<"Cannot GetParameters() because m_InputParametersPointer is NULL. Perhaps SetCoefficientImage() has been called causing the NULL pointer." );
659 return (*m_InputParametersPointer);
663 // Get the parameters
664 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
666 typename BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
668 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
669 ::GetFixedParameters( void ) const
671 RegionType resRegion = this->GetGridRegion( );
673 for (unsigned int i=0; i<NInputDimensions; i++)
675 this->m_FixedParameters[i] = (resRegion.GetSize())[i];
677 for (unsigned int i=0; i<NInputDimensions; i++)
679 this->m_FixedParameters[NInputDimensions+i] = (this->GetGridOrigin())[i];
681 for (unsigned int i=0; i<NInputDimensions; i++)
683 this->m_FixedParameters[2*NInputDimensions+i] = (this->GetGridSpacing())[i];
685 for (unsigned int di=0; di<NInputDimensions; di++)
687 for (unsigned int dj=0; dj<NInputDimensions; dj++)
689 this->m_FixedParameters[3*NInputDimensions+(di*NInputDimensions+dj)] = (this->GetGridDirection())[di][dj];
693 //JV add splineOrders
694 for (unsigned int i=0; i<NInputDimensions; i++)
696 this->m_FixedParameters[(3+NInputDimensions)*NInputDimensions+i] = (this->GetSplineOrders())[i];
699 //JV add LUTsamplingFactor
700 for (unsigned int i=0; i<NInputDimensions; i++)
702 this->m_FixedParameters[(4+NInputDimensions)*NInputDimensions+i] = (this->GetLUTSamplingFactors())[i];
706 this->m_FixedParameters[(5+NInputDimensions)*NInputDimensions]=(double)((size_t) m_Mask);
708 //JV add the bulktransform pointer
709 this->m_FixedParameters[(5+NInputDimensions)*NInputDimensions+1]=(double)((size_t) m_BulkTransform);
711 return (this->m_FixedParameters);
715 // Set the B-Spline coefficients using input images
716 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
718 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
719 ::SetCoefficientImage( CoefficientImagePointer image )
722 this->SetGridRegion( image->GetBufferedRegion() );
723 this->SetGridSpacing( image->GetSpacing() );
724 this->SetGridDirection( image->GetDirection() );
725 this->SetGridOrigin( image->GetOrigin() );
726 m_CoefficientImage = image;
728 // Update the interpolator
729 m_VectorInterpolator->SetInputImage(this->GetCoefficientImage());
731 // Clean up buffered parameters
732 m_InternalParametersBuffer = ParametersType( 0 );
733 m_InputParametersPointer = NULL;
739 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
741 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
742 ::PrintSelf(std::ostream &os, itk::Indent indent) const
745 this->Superclass::PrintSelf(os, indent);
747 os << indent << "GridRegion: " << m_GridRegion << std::endl;
748 os << indent << "GridOrigin: " << m_GridOrigin << std::endl;
749 os << indent << "GridSpacing: " << m_GridSpacing << std::endl;
750 os << indent << "GridDirection: " << m_GridDirection << std::endl;
751 os << indent << "IndexToPoint: " << m_IndexToPoint << std::endl;
752 os << indent << "PointToIndex: " << m_PointToIndex << std::endl;
754 os << indent << "CoefficientImage: [ ";
755 os << m_CoefficientImage.GetPointer() << " ]" << std::endl;
757 os << indent << "WrappedImage: [ ";
758 os << m_WrappedImage.GetPointer() << " ]" << std::endl;
760 os << indent << "InputParametersPointer: "
761 << m_InputParametersPointer << std::endl;
762 os << indent << "ValidRegion: " << m_ValidRegion << std::endl;
763 os << indent << "LastJacobianIndex: " << m_LastJacobianIndex << std::endl;
764 os << indent << "BulkTransform: ";
765 os << m_BulkTransform << std::endl;
767 if ( m_BulkTransform )
769 os << indent << "BulkTransformType: "
770 << m_BulkTransform->GetNameOfClass() << std::endl;
772 os << indent << "VectorBSplineInterpolator: ";
773 os << m_VectorInterpolator.GetPointer() << std::endl;
774 os << indent << "Mask: ";
775 os << m_Mask<< std::endl;
780 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
782 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
783 ::InsideValidRegion( const ContinuousIndexType& index ) const
787 if ( !m_ValidRegion.IsInside( index ) )
791 //JV verify for each input dimension
794 typedef typename ContinuousIndexType::ValueType ValueType;
795 for( unsigned int j = 0; j < InputDimension; j++ )
797 if (m_SplineOrderOdd[j])
799 if ( index[j] >= static_cast<ValueType>( m_ValidRegionLast[j] ) )
813 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
814 typename BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
816 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
817 ::TransformPoint(const InputPointType &inputPoint) const
820 InputPointType transformedPoint;
821 OutputPointType outputPoint;
824 if ( m_BulkTransform )
826 transformedPoint = m_BulkTransform->TransformPoint( inputPoint );
830 transformedPoint = inputPoint;
833 // Deformable transform
834 if ( m_CoefficientImage )
836 // Check if inside mask
837 if(m_Mask && !(m_Mask->IsInside(inputPoint) ) )
839 // Outside: no (deformable) displacement
840 return transformedPoint;
843 // Check if inside valid region
845 ContinuousIndexType index;
846 this->TransformPointToContinuousIndex( inputPoint, index );
847 inside = this->InsideValidRegion( index );
850 // Outside: no (deformable) displacement
851 return transformedPoint;
854 // Call the vector interpolator
855 itk::Vector<TCoordRep,OutputDimension> displacement=m_VectorInterpolator->EvaluateAtContinuousIndex(index);
857 // Return the results
858 outputPoint = transformedPoint+displacement;
864 itkWarningMacro( << "B-spline coefficients have not been set" );
865 outputPoint = transformedPoint;
873 //JV Deformably transform a point
874 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
875 typename BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
877 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
878 ::DeformablyTransformPoint(const InputPointType &inputPoint) const
880 OutputPointType outputPoint;
881 if ( m_CoefficientImage )
884 // Check if inside mask
885 if(m_Mask && !(m_Mask->IsInside(inputPoint) ) )
887 // Outside: no (deformable) displacement
893 ContinuousIndexType index;
894 this->TransformPointToContinuousIndex( inputPoint, index );
895 inside = this->InsideValidRegion( index );
899 //outside: no (deformable) displacement
901 //return outputPoint;
904 // Call the vector interpolator
905 itk::Vector<TCoordRep,OutputDimension> displacement=m_VectorInterpolator->EvaluateAtContinuousIndex(index);
907 // Return the results
908 outputPoint = inputPoint+displacement;
911 // No coefficients available
914 itkWarningMacro( << "B-spline coefficients have not been set" );
915 outputPoint = inputPoint;
922 // JV weights are identical as for transformpoint, could be done simultaneously in metric!!!!
923 // Compute the Jacobian in one position
924 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
926 typename BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
928 BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
929 ::GetJacobian( const InputPointType & point ) const
931 // Can only compute Jacobian if parameters are set via
932 // SetParameters or SetParametersByValue
933 // if( m_InputParametersPointer == NULL )
935 // itkExceptionMacro( <<"Cannot compute Jacobian: parameters not set" );
939 //========================================================
940 // Zero all components of jacobian
941 //========================================================
942 // JV not thread safe (m_LastJacobianIndex), instantiate N transforms
943 // NOTE: for efficiency, we only need to zero out the coefficients
944 // that got fill last time this method was called.
948 //Define the region for each jacobian image
949 m_SupportRegion.SetIndex( m_LastJacobianIndex );
951 //Initialize the iterators
952 for ( j = 0; j < OutputDimension; j++ )
953 m_Iterator[j] = IteratorType( m_JacobianImage[j], m_SupportRegion);
955 //Set the previously-set to zero
956 while ( ! (m_Iterator[0]).IsAtEnd() )
958 for ( j = 0; j < OutputDimension; j++ )
960 m_Iterator[j].Set( m_ZeroVector );
965 //========================================================
966 // For each dimension, copy the weight to the support region
967 //========================================================
969 // Check if inside mask
970 if(m_Mask && !(m_Mask->IsInside(point) ) )
972 // Outside: no (deformable) displacement
973 return this->m_Jacobian;
977 this->TransformPointToContinuousIndex( point, m_Index );
979 // NOTE: if the support region does not lie totally within the grid
980 // we assume zero displacement and return the input point
981 if ( !this->InsideValidRegion( m_Index ) )
983 return this->m_Jacobian;
986 //Compute interpolation weights
987 const WeightsDataType *weights=NULL;
988 m_VectorInterpolator->EvaluateWeightsAtContinuousIndex( m_Index, &weights, m_LastJacobianIndex);
989 m_SupportRegion.SetIndex( m_LastJacobianIndex );
991 //Reset the iterators
992 for ( j = 0; j < OutputDimension; j++ )
993 m_Iterator[j] = IteratorType( m_JacobianImage[j], m_SupportRegion);
995 // For each dimension, copy the weight to the support region
996 while ( ! (m_Iterator[0]).IsAtEnd() )
998 //copy weight to jacobian image
999 for ( j = 0; j < OutputDimension; j++ )
1001 m_ZeroVector[j]=*weights;
1002 (m_Iterator[j]).Set( m_ZeroVector);
1003 m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
1006 // go to next coefficient in the support region
1010 // Return the result
1011 return this->m_Jacobian;
1016 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
1018 BSplineDeformableTransform<TCoordRep, NInputDimensions,NOutputDimensions>
1019 ::TransformPointToContinuousIndex( const InputPointType & point, ContinuousIndexType & index ) const
1023 itk::Vector<double, OutputDimension> tvector;
1025 for ( j = 0; j < OutputDimension; j++ )
1027 tvector[j] = point[j] - this->m_GridOrigin[j];
1030 itk::Vector<double, OutputDimension> cvector;
1032 cvector = m_PointToIndex * tvector;
1034 for ( j = 0; j < OutputDimension; j++ )
1036 index[j] = static_cast< typename ContinuousIndexType::CoordRepType >( cvector[j] );