]> Creatis software - clitk.git/blob - registration/clitkBSplineDeformableTransform.txx
small bugs corrections
[clitk.git] / registration / clitkBSplineDeformableTransform.txx
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 __clitkBSplineDeformableTransform_txx
19 #define __clitkBSplineDeformableTransform_tx
20 #include "clitkBSplineDeformableTransform.h"
21
22 //itk
23 #include "itkContinuousIndex.h"
24 #include "itkImageRegionIterator.h"
25 #include "itkImageRegionConstIterator.h"
26 #include "itkImageRegionConstIteratorWithIndex.h"
27 #include "itkIdentityTransform.h"
28
29 namespace clitk
30 {
31
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)
36   {
37     unsigned int i;
38     
39     //JV default spline order
40     for ( i=0;i<InputDimension; i++)
41       m_SplineOrders[i]=3;
42
43     //JV default LUTSamplingfactor
44     for ( i=0;i<InputDimension; i++)
45       m_LUTSamplingFactors[i]=20;
46
47     for (i=0;i<InputDimension; i++)
48       m_SupportSize[i] = m_SplineOrders[i]+1;
49
50     //Instantiate a Vector Bspline Interpolator
51     m_VectorInterpolator =VectorInterpolatorType::New();
52     m_VectorInterpolator->SetLUTSamplingFactors(m_LUTSamplingFactors);
53     m_VectorInterpolator->SetSplineOrders(m_SplineOrders);
54
55     // Set Bulk transform to NULL (no bulk performed)
56     m_BulkTransform = NULL;
57
58     // Mask
59     m_Mask=NULL;
60
61     // Default grid size is zero
62     typename RegionType::SizeType size;
63     typename RegionType::IndexType index;
64     size.Fill( 0 );
65     index.Fill( 0 );
66     m_GridRegion.SetSize( size );
67     m_GridRegion.SetIndex( index );
68
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
72
73     m_InternalParametersBuffer = ParametersType(0);
74     // Make sure the parameters pointer is not NULL after construction.
75     m_InputParametersPointer = &m_InternalParametersBuffer;
76
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;
84   
85     // Variables for computing interpolation
86     for (i=0; i <InputDimension;i++)
87       {
88         m_Offset[i] = m_SplineOrders[i] / 2;
89         m_SplineOrderOdd[i] = m_SplineOrders[i] % 2;
90       }
91     m_ValidRegion = m_GridRegion;
92         
93     // Initialize jacobian images 
94     for (i=0; i <OutputDimension;i++)
95       {
96         m_JacobianImage[i] = JacobianImageType::New();
97         m_JacobianImage[i]->SetRegions( m_GridRegion );
98         m_JacobianImage[i]->SetOrigin( m_GridOrigin.GetDataPointer() );
99         m_JacobianImage[i]->SetSpacing( m_GridSpacing.GetDataPointer() );
100         m_JacobianImage[i]->SetDirection( m_GridDirection );
101       }
102     
103     /** Fixed Parameters store the following information:
104      *     Grid Size
105      *     Grid Origin
106      *     Grid Spacing
107      *     Grid Direction */
108     //JV we add the splineOrders, LUTsamplingfactor, m_Mask  and m_BulkTransform
109     /*
110       Spline orders
111       Sampling factors
112       m_Mask
113       m_BulkTransform
114       The size of these is equal to the  NInputDimensions
115     *********************************************************/
116     this->m_FixedParameters.SetSize ( NInputDimensions * (NInputDimensions + 5)+2 );
117     this->m_FixedParameters.Fill ( 0.0 );
118     for ( i=0; i<NInputDimensions; i++)
119       {
120         this->m_FixedParameters[2*NInputDimensions+i] = m_GridSpacing[i];
121       }
122     for (unsigned int di=0; di<NInputDimensions; di++)
123       {
124         for (unsigned int dj=0; dj<NInputDimensions; dj++)
125           {
126             this->m_FixedParameters[3*NInputDimensions+(di*NInputDimensions+dj)] = m_GridDirection[di][dj];
127           }
128       }
129
130     //JV add splineOrders
131     for ( i=0; i<NInputDimensions; i++)
132       {
133         this->m_FixedParameters[ ( (3+NInputDimensions)*NInputDimensions)+i] = (this->GetSplineOrders())[i];
134       }
135     
136     //JV add LUTsamplingFactors
137     for ( i=0; i<NInputDimensions; i++)
138       {
139         this->m_FixedParameters[( (4+NInputDimensions)*NInputDimensions)+i ] = this->GetLUTSamplingFactors()[i];
140       }
141
142     // JV add the mask pointer
143     this->m_FixedParameters[( (5+NInputDimensions)*NInputDimensions)]=(double)((size_t)m_Mask);
144
145     // JV add the bulkTransform pointer
146     this->m_FixedParameters[( (5+NInputDimensions)*NInputDimensions) +1]=(double)((size_t)m_BulkTransform);
147
148
149     // Calculate the PointToIndex matrices
150     DirectionType scale;
151     for( unsigned int i=0; i<OutputDimension; i++)
152       {
153         scale[i][i] = m_GridSpacing[i];
154       }
155
156     m_IndexToPoint = m_GridDirection * scale;
157     m_PointToIndex = m_IndexToPoint.GetInverse();
158
159     // Jacobian
160     m_LastJacobianIndex = m_ValidRegion.GetIndex();
161
162     //JV initialize some variables for jacobian calculation
163     m_SupportRegion.SetSize(m_SupportSize);
164     m_SupportIndex.Fill(0);
165     m_SupportRegion.SetIndex(m_SupportIndex);
166     for (  i = 0; i < OutputDimension; i++ ) 
167       m_ZeroVector[i]=itk::NumericTraits<JacobianValueType>::Zero;
168   }
169     
170
171   // Destructor
172   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
173   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
174   ::~BSplineDeformableTransform()
175   {
176
177   }
178
179
180   // JV set Spline Order
181   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
182   void
183   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
184   ::SetSplineOrder(const unsigned int & splineOrder)
185   {
186     SizeType splineOrders;
187     for (unsigned int i=0;i<InputDimension;i++)splineOrders[i]=splineOrder;
188
189     this->SetSplineOrders(splineOrders);
190   }
191    
192
193   // JV set Spline Orders
194   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
195   void
196   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
197   ::SetSplineOrders(const SizeType & splineOrders)
198   {
199     if (m_SplineOrders!=splineOrders)
200       {
201         m_SplineOrders=splineOrders;
202         
203         //update the interpolation function
204         m_VectorInterpolator->SetSplineOrders(m_SplineOrders);
205         
206         //update the variables for computing interpolation
207         for (unsigned int i=0; i <InputDimension;i++)
208         {
209           m_SupportSize[i] = m_SplineOrders[i]+1;
210           m_Offset[i] = m_SplineOrders[i] / 2;
211           m_SplineOrderOdd[i] = m_SplineOrders[i] % 2;
212         }
213         this->Modified();
214       }
215   }
216
217
218   // JV set sampling factor
219   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
220   void
221   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
222   ::SetLUTSamplingFactor( const int & samplingFactor)
223   {
224     SizeType samplingFactors;
225     for (unsigned int i=0; i<NInputDimensions; i++)
226       samplingFactors[i]=samplingFactor;
227     
228     //update the interpolation function
229     this->SetLUTSamplingFactors(samplingFactors);
230   }
231
232
233   // JV set sampling factors
234   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
235   void
236   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
237   ::SetLUTSamplingFactors( const SizeType & samplingFactors)
238   {
239     if(m_LUTSamplingFactors!=samplingFactors)
240       {
241         for (unsigned int i=0; i<NInputDimensions; i++)
242           m_LUTSamplingFactors[i]=samplingFactors[i];
243         
244         //update the interpolation function
245         m_VectorInterpolator->SetLUTSamplingFactors(m_LUTSamplingFactors);
246         
247         this->Modified();
248       }
249   }
250
251
252   // Get the number of parameters
253   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
254   unsigned int
255   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
256   ::GetNumberOfParameters(void) const
257   {
258     // The number of parameters equal OutputDimension * number of
259     // of pixels in the grid region.
260     return ( static_cast<unsigned int>( OutputDimension ) *
261              static_cast<unsigned int>( m_GridRegion.GetNumberOfPixels() ) );
262
263   }
264
265
266   // Get the number of parameters per dimension
267   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
268   unsigned int
269   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
270   ::GetNumberOfParametersPerDimension(void) const
271   {
272     // The number of parameters per dimension equal number of
273     // of pixels in the grid region.
274     return ( static_cast<unsigned int>( m_GridRegion.GetNumberOfPixels() ) );
275
276   }
277
278
279   // Set the grid region
280   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
281   void
282   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
283   ::SetGridRegion( const RegionType & region )
284   {
285     if ( m_GridRegion != region )
286       {
287         m_GridRegion = region;
288
289         // set regions for each coefficient and jacobian image
290         m_WrappedImage->SetRegions( m_GridRegion );
291         for (unsigned int j=0; j <OutputDimension;j++)
292           m_JacobianImage[j]->SetRegions( m_GridRegion );
293
294         // Set the valid region
295         // If the grid spans the interval [start, last].
296         // The valid interval for evaluation is [start+offset, last-offset]
297         // when spline order is even.
298         // The valid interval for evaluation is [start+offset, last-offset)
299         // when spline order is odd.
300         // Where offset = vcl_floor(spline / 2 ).
301         // Note that the last pixel is not included in the valid region
302         // with odd spline orders.
303         typename RegionType::SizeType size = m_GridRegion.GetSize();
304         typename RegionType::IndexType index = m_GridRegion.GetIndex();
305         for ( unsigned int j = 0; j < NInputDimensions; j++ )
306           {
307             index[j] += 
308               static_cast< typename RegionType::IndexValueType >( m_Offset[j] );
309             size[j] -= 
310               static_cast< typename RegionType::SizeValueType> ( 2 * m_Offset[j] );
311             m_ValidRegionLast[j] = index[j] +
312               static_cast< typename RegionType::IndexValueType >( size[j] ) - 1;
313           }
314         m_ValidRegion.SetSize( size );
315         m_ValidRegion.SetIndex( index );
316
317         // If we are using the default parameters, update their size and set to identity.
318         // Input parameters point to internal buffer => using default parameters.
319         if (m_InputParametersPointer == &m_InternalParametersBuffer)
320           {
321             // Check if we need to resize the default parameter buffer.
322             if ( m_InternalParametersBuffer.GetSize() != this->GetNumberOfParameters() )
323               {
324                 m_InternalParametersBuffer.SetSize( this->GetNumberOfParameters() );
325                 // Fill with zeros for identity.
326                 m_InternalParametersBuffer.Fill( 0 );
327               }
328           }
329
330         this->Modified();
331       }
332   }
333
334
335   // Set the grid spacing
336   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
337   void
338   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
339   ::SetGridSpacing( const SpacingType & spacing )
340   {
341     if ( m_GridSpacing != spacing )
342       {
343         m_GridSpacing = spacing;
344
345         // Set spacing for each coefficient and jacobian image
346         m_WrappedImage->SetSpacing( m_GridSpacing.GetDataPointer() );
347         for (unsigned int j=0; j <OutputDimension; j++) m_JacobianImage[j]->SetSpacing( m_GridSpacing.GetDataPointer() );
348         
349         // Set scale
350         DirectionType scale;
351         for( unsigned int i=0; i<OutputDimension; i++)
352           {
353             scale[i][i] = m_GridSpacing[i];
354           }
355
356         m_IndexToPoint = m_GridDirection * scale;
357         m_PointToIndex = m_IndexToPoint.GetInverse();
358
359         this->Modified();
360       }
361
362   }
363
364   // Set the grid direction
365   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
366   void
367   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
368   ::SetGridDirection( const DirectionType & direction )
369   {
370     if ( m_GridDirection != direction )
371       {
372         m_GridDirection = direction;
373
374         // Set direction for each coefficient and jacobian image
375         m_WrappedImage->SetDirection( m_GridDirection );
376         for (unsigned int j=0; j <OutputDimension; j++) m_JacobianImage[j]->SetDirection( m_GridDirection );
377         
378         // Set scale
379         DirectionType scale;
380         for( unsigned int i=0; i<OutputDimension; i++)
381           {
382             scale[i][i] = m_GridSpacing[i];
383           }
384
385         m_IndexToPoint = m_GridDirection * scale;
386         m_PointToIndex = m_IndexToPoint.GetInverse();
387
388         this->Modified();
389       }
390
391   }
392
393
394   // Set the grid origin
395   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
396   void
397   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
398   ::SetGridOrigin( const OriginType& origin )
399   {
400     if ( m_GridOrigin != origin )
401       {
402         m_GridOrigin = origin;
403
404         // Set origin for each coefficient and jacobianimage
405         m_WrappedImage->SetOrigin( m_GridOrigin.GetDataPointer() );
406         for (unsigned int j=0; j <OutputDimension; j++) m_JacobianImage[j]->SetOrigin( m_GridOrigin.GetDataPointer() );
407         
408         this->Modified();
409       }
410
411   }
412
413
414   // Set the parameters
415   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
416   void
417   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
418   ::SetIdentity()
419   {
420     if( m_InputParametersPointer )
421       {
422         ParametersType * parameters =
423           const_cast<ParametersType *>( m_InputParametersPointer );
424         parameters->Fill( 0.0 );
425         this->Modified();
426       }
427     else 
428       {
429         itkExceptionMacro( << "Input parameters for the spline haven't been set ! "
430                            << "Set them using the SetParameters or SetCoefficientImage method first." );
431       }
432   }
433
434   // Set the parameters
435   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
436   void
437   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
438   ::SetParameters( const ParametersType & parameters )
439   {
440
441     // Check if the number of parameters match the
442     // Expected number of parameters
443     if ( parameters.Size() != this->GetNumberOfParameters() )
444       {
445         itkExceptionMacro(<<"Mismatched between parameters size "
446                           << parameters.size() 
447                           << " and region size "
448                           << m_GridRegion.GetNumberOfPixels() );
449       }
450
451     // Clean up buffered parameters
452     m_InternalParametersBuffer = ParametersType( 0 );
453
454     // Keep a reference to the input parameters
455     m_InputParametersPointer = &parameters;
456
457     // Wrap flat array as images of coefficients
458     this->WrapAsImages();
459
460     //Set input to vector interpolator
461     m_VectorInterpolator->SetInputImage(this->GetCoefficientImage());
462
463     // Modified is always called since we just have a pointer to the
464     // parameters and cannot know if the parameters have changed.
465     this->Modified();
466   }
467
468
469   // Set the Fixed Parameters
470   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
471   void
472   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
473   ::SetFixedParameters( const ParametersType & parameters )
474   {
475  
476     // JV number should be exact, no defaults for spacing
477     if ( parameters.Size() != NInputDimensions * (5 + NInputDimensions)+2 )
478       {
479         itkExceptionMacro(<< "Mismatched between parameters size "
480                           << parameters.size() 
481                           << " and number of fixed parameters "
482                           << NInputDimensions * (5 + NInputDimensions)+2 );
483       }
484     /********************************************************* 
485     Fixed Parameters store the following information:
486         Grid Size
487         Grid Origin
488         Grid Spacing
489         Grid Direction */
490     // JV we add the splineOrders, LUTsamplingfactor, mask pointer and bulktransform pointer
491     /*
492       Spline orders
493       Sampling factors
494       m_Mask
495       m_BulkTransform
496     The size of these is equal to the  NInputDimensions
497     *********************************************************/
498   
499     /** Set the Grid Parameters */
500     SizeType   gridSize;
501     for (unsigned int i=0; i<NInputDimensions; i++)
502       {
503         gridSize[i] = static_cast<int> (parameters[i]);
504       }
505     RegionType bsplineRegion;
506     bsplineRegion.SetSize( gridSize );
507   
508     /** Set the Origin Parameters */
509     OriginType origin;
510     for (unsigned int i=0; i<NInputDimensions; i++)
511       {
512         origin[i] = parameters[NInputDimensions+i];
513       }
514   
515     /** Set the Spacing Parameters */
516     SpacingType spacing;
517     for (unsigned int i=0; i<NInputDimensions; i++)
518       {
519         spacing[i] = parameters[2*NInputDimensions+i];
520       }
521
522     /** Set the Direction Parameters */
523     DirectionType direction;
524     for (unsigned int di=0; di<NInputDimensions; di++)
525       {
526         for (unsigned int dj=0; dj<NInputDimensions; dj++)
527           {
528             direction[di][dj] = parameters[3*NInputDimensions+(di*NInputDimensions+dj)];
529           }
530       }
531
532     //JV add the spline orders
533     SizeType splineOrders;
534     for (unsigned int i=0; i<NInputDimensions; i++)
535       {
536         splineOrders[i]=parameters[(3+NInputDimensions)*NInputDimensions+i];
537       }
538
539     //JV add the LUT sampling factor
540     SizeType  samplingFactors;
541     for (unsigned int i=0; i<NInputDimensions; i++)
542       {
543         samplingFactors[i]=parameters[(4+NInputDimensions)*NInputDimensions+i];
544       }
545
546     //JV add the MaskPointer
547     m_Mask=((MaskPointer)( (size_t)parameters[(5+NInputDimensions)*NInputDimensions]));
548
549     //JV add the MaskPointer
550     m_BulkTransform=((BulkTransformPointer)( (size_t)parameters[(5+NInputDimensions)*NInputDimensions+1]));
551
552     //Set the members
553     this->SetGridSpacing( spacing );
554     this->SetGridDirection( direction );
555     this->SetGridOrigin( origin );
556     this->SetGridRegion( bsplineRegion );
557
558     //JV add the LUT parameters
559     this->SetSplineOrders( splineOrders );
560     this->SetLUTSamplingFactors( samplingFactors );
561       
562   }
563
564
565   // Wrap flat parameters as images
566   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
567   void
568   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
569   ::WrapAsImages()
570   {
571     //JV Wrap parameter array in vectorial image, changed parameter order: A1x A1y A1z, A2x ....
572     PixelType * dataPointer =reinterpret_cast<PixelType *>( const_cast<double *>(m_InputParametersPointer->data_block() )) ;
573     unsigned int numberOfPixels = m_GridRegion.GetNumberOfPixels();
574     
575     m_WrappedImage->GetPixelContainer()->SetImportPointer( dataPointer,numberOfPixels);//InputDimension
576     m_CoefficientImage = m_WrappedImage;
577     
578     //JV Wrap jacobian into OutputDimension X Vectorial images
579     this->m_Jacobian.set_size( OutputDimension, this->GetNumberOfParameters() );
580
581     // Use memset to set the memory
582     JacobianPixelType * jacobianDataPointer = reinterpret_cast<JacobianPixelType *>(this->m_Jacobian.data_block());
583     memset(jacobianDataPointer, 0,  OutputDimension*numberOfPixels*sizeof(JacobianPixelType));
584     m_LastJacobianIndex = m_ValidRegion.GetIndex();
585     m_NeedResetJacobian = false;
586
587     for (unsigned int j=0; j<OutputDimension; j++)
588       {
589         m_JacobianImage[j]->GetPixelContainer()->
590           SetImportPointer( jacobianDataPointer, numberOfPixels );
591         jacobianDataPointer += numberOfPixels;
592      }
593   }
594
595
596   // Set the parameters by value
597   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
598   void
599   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
600   ::SetParametersByValue( const ParametersType & parameters )
601   {
602
603     // Check if the number of parameters match the
604     // Expected number of parameters
605     if ( parameters.Size() != this->GetNumberOfParameters() )
606       {
607         itkExceptionMacro(<<"Mismatched between parameters size "
608                           << parameters.size() 
609                           << " and region size "
610                           << m_GridRegion.GetNumberOfPixels() );
611       }
612
613     // Copy it
614     m_InternalParametersBuffer = parameters;
615     m_InputParametersPointer = &m_InternalParametersBuffer;
616
617     // Wrap flat array as images of coefficients
618     this->WrapAsImages();
619
620     // Modified is always called since we just have a pointer to the
621     // Parameters and cannot know if the parameters have changed.
622     this->Modified();
623
624   }
625
626   // Get the parameters
627   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
628   const 
629   typename BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
630   ::ParametersType &
631   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
632   ::GetParameters( void ) const
633   {
634     /** NOTE: For efficiency, this class does not keep a copy of the parameters - 
635      * it just keeps pointer to input parameters. 
636      */
637     if (NULL == m_InputParametersPointer)
638       {
639         itkExceptionMacro( <<"Cannot GetParameters() because m_InputParametersPointer is NULL. Perhaps SetCoefficientImage() has been called causing the NULL pointer." );
640       }
641
642     return (*m_InputParametersPointer);
643   }
644
645
646   // Get the parameters
647   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
648   const 
649   typename BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
650   ::ParametersType &
651   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
652   ::GetFixedParameters( void ) const
653   {
654     RegionType resRegion = this->GetGridRegion(  );
655   
656     for (unsigned int i=0; i<NInputDimensions; i++)
657       {
658         this->m_FixedParameters[i] = (resRegion.GetSize())[i];
659       }
660     for (unsigned int i=0; i<NInputDimensions; i++)
661       {
662         this->m_FixedParameters[NInputDimensions+i] = (this->GetGridOrigin())[i];
663       } 
664     for (unsigned int i=0; i<NInputDimensions; i++)
665       {
666         this->m_FixedParameters[2*NInputDimensions+i] =  (this->GetGridSpacing())[i];
667       }
668     for (unsigned int di=0; di<NInputDimensions; di++)
669       {
670         for (unsigned int dj=0; dj<NInputDimensions; dj++)
671           {
672             this->m_FixedParameters[3*NInputDimensions+(di*NInputDimensions+dj)] = (this->GetGridDirection())[di][dj];
673           }
674       }
675
676     //JV add splineOrders
677     for (unsigned int i=0; i<NInputDimensions; i++)
678       {
679         this->m_FixedParameters[(3+NInputDimensions)*NInputDimensions+i] = (this->GetSplineOrders())[i];
680       }
681     
682     //JV add LUTsamplingFactor
683     for (unsigned int i=0; i<NInputDimensions; i++)
684       {
685         this->m_FixedParameters[(4+NInputDimensions)*NInputDimensions+i] = (this->GetLUTSamplingFactors())[i];
686       }
687     
688     //JV add the mask
689     this->m_FixedParameters[(5+NInputDimensions)*NInputDimensions]=(double)((size_t) m_Mask);
690
691     //JV add the bulktransform pointer
692     this->m_FixedParameters[(5+NInputDimensions)*NInputDimensions+1]=(double)((size_t) m_BulkTransform);
693     
694     return (this->m_FixedParameters);
695   }
696
697   
698   // Set the B-Spline coefficients using input images
699   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
700   void 
701   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
702   ::SetCoefficientImage( CoefficientImagePointer image )
703   {
704     
705     this->SetGridRegion( image->GetBufferedRegion() );
706     this->SetGridSpacing( image->GetSpacing() );
707     this->SetGridDirection( image->GetDirection() );
708     this->SetGridOrigin( image->GetOrigin() );
709     m_CoefficientImage = image;
710
711     // Update the interpolator    
712     m_VectorInterpolator->SetInputImage(this->GetCoefficientImage());
713
714     // Clean up buffered parameters
715     m_InternalParametersBuffer = ParametersType( 0 );
716     m_InputParametersPointer  = NULL;
717       
718   }  
719
720
721   // Print self
722   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
723   void
724   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
725   ::PrintSelf(std::ostream &os, itk::Indent indent) const
726   {
727
728     this->Superclass::PrintSelf(os, indent);
729
730     os << indent << "GridRegion: " << m_GridRegion << std::endl;
731     os << indent << "GridOrigin: " << m_GridOrigin << std::endl;
732     os << indent << "GridSpacing: " << m_GridSpacing << std::endl;
733     os << indent << "GridDirection: " << m_GridDirection << std::endl;
734     os << indent << "IndexToPoint: " << m_IndexToPoint << std::endl;
735     os << indent << "PointToIndex: " << m_PointToIndex << std::endl;
736
737     os << indent << "CoefficientImage: [ ";
738     os << m_CoefficientImage.GetPointer() << " ]" << std::endl;
739
740     os << indent << "WrappedImage: [ ";
741     os << m_WrappedImage.GetPointer() << " ]" << std::endl;
742  
743     os << indent << "InputParametersPointer: " 
744        << m_InputParametersPointer << std::endl;
745     os << indent << "ValidRegion: " << m_ValidRegion << std::endl;
746     os << indent << "LastJacobianIndex: " << m_LastJacobianIndex << std::endl;
747     os << indent << "BulkTransform: ";
748     os << m_BulkTransform << std::endl;
749
750     if ( m_BulkTransform )
751       {
752         os << indent << "BulkTransformType: " 
753            << m_BulkTransform->GetNameOfClass() << std::endl;
754       }
755     os << indent << "VectorBSplineInterpolator: ";
756     os << m_VectorInterpolator.GetPointer() << std::endl;
757     os << indent << "Mask: ";
758     os << m_Mask<< std::endl;
759   }
760
761
762   // Verify 
763   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
764   bool 
765   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
766   ::InsideValidRegion( const ContinuousIndexType& index ) const
767   {
768     bool inside = true;
769
770     if ( !m_ValidRegion.IsInside( index ) )
771       {
772         inside = false;
773       }
774     //JV verify for each input dimension
775     if ( inside)
776       {
777         typedef typename ContinuousIndexType::ValueType ValueType;
778         for( unsigned int j = 0; j < InputDimension; j++ )
779           {
780             if (m_SplineOrderOdd[j])
781               {
782                 if ( index[j] >= static_cast<ValueType>( m_ValidRegionLast[j] ) )
783                   { 
784                     inside = false;
785                     break;
786                   }
787               }
788           }
789       }
790     
791     return inside;
792   }
793
794
795   // Transform a point
796   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
797   typename BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
798   ::OutputPointType
799   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
800   ::TransformPoint(const InputPointType &inputPoint) const 
801   {
802   
803    InputPointType transformedPoint;
804    OutputPointType outputPoint;
805
806    // BulkTransform
807    if ( m_BulkTransform )
808      {
809        transformedPoint = m_BulkTransform->TransformPoint( inputPoint );
810      }
811    else
812       {
813         transformedPoint = inputPoint;
814       }
815    
816     // Deformable transform
817     if ( m_CoefficientImage )
818       {
819         // Check if inside mask
820         if(m_Mask &&  !(m_Mask->IsInside(inputPoint) ) )
821           {
822             // Outside: no (deformable) displacement
823             return transformedPoint;
824           }             
825         
826         // Check if inside valid region
827         bool inside = true;
828         ContinuousIndexType index;
829         this->TransformPointToContinuousIndex( inputPoint, index );
830         inside = this->InsideValidRegion( index );
831         if ( !inside )
832           {
833             // Outside: no (deformable) displacement
834             return transformedPoint;
835           }
836         
837         // Call the vector interpolator
838         itk::Vector<TCoordRep,OutputDimension> displacement=m_VectorInterpolator->EvaluateAtContinuousIndex(index);
839
840         // Return the results
841         outputPoint = transformedPoint+displacement;
842
843       }
844
845     else
846       {
847         itkWarningMacro( << "B-spline coefficients have not been set" );
848         outputPoint = transformedPoint;
849       }
850     
851     return outputPoint;
852   }
853
854
855
856   //JV Deformably transform a point
857   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
858   typename BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
859   ::OutputPointType
860   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
861   ::DeformablyTransformPoint(const InputPointType &inputPoint) const 
862   {
863     OutputPointType outputPoint;
864     if ( m_CoefficientImage )
865       {
866
867         // Check if inside mask
868         if(m_Mask && !(m_Mask->IsInside(inputPoint) ) )
869           {
870             // Outside: no (deformable) displacement
871             return inputPoint;
872           }     
873         
874         // Check if inside
875         bool inside = true;
876         ContinuousIndexType index;
877         this->TransformPointToContinuousIndex( inputPoint, index );
878         inside = this->InsideValidRegion( index );
879
880         if ( !inside )
881           {
882             //outside: no (deformable) displacement
883              return inputPoint;
884              //return outputPoint;
885           }
886
887         // Call the vector interpolator
888         itk::Vector<TCoordRep,OutputDimension> displacement=m_VectorInterpolator->EvaluateAtContinuousIndex(index);
889         
890         // Return the results
891         outputPoint = inputPoint+displacement;
892       }
893
894     // No coefficients available
895     else
896       {
897         itkWarningMacro( << "B-spline coefficients have not been set" );
898         outputPoint = inputPoint;
899       }
900    
901     return outputPoint;
902   }
903   
904
905   // JV weights are identical as for transformpoint, could be done simultaneously in metric!!!!
906   // Compute the Jacobian in one position 
907   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
908   const 
909   typename BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
910   ::JacobianType & 
911   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
912   ::GetJacobian( const InputPointType & point ) const
913   {
914     // Can only compute Jacobian if parameters are set via
915     // SetParameters or SetParametersByValue
916     //     if( m_InputParametersPointer == NULL )
917     //       {
918     //  itkExceptionMacro( <<"Cannot compute Jacobian: parameters not set" );
919     //       }
920
921     if (m_NeedResetJacobian)
922       ResetJacobian();
923
924     //========================================================
925     // For each dimension, copy the weight to the support region
926     //========================================================
927
928     // Check if inside mask
929     if(m_Mask &&  !(m_Mask->IsInside(point) ) )
930       {
931         // Outside: no (deformable) displacement
932         return this->m_Jacobian;
933       } 
934
935     //Get index   
936     this->TransformPointToContinuousIndex( point, m_Index );
937
938     // NOTE: if the support region does not lie totally within the grid
939     // we assume zero displacement and return the input point
940     if ( !this->InsideValidRegion( m_Index ) )
941       {
942         return this->m_Jacobian;
943       }
944
945     //Compute interpolation weights
946     const WeightsDataType *weights=NULL;
947     m_VectorInterpolator->EvaluateWeightsAtContinuousIndex( m_Index, &weights, m_LastJacobianIndex);
948     m_SupportRegion.SetIndex( m_LastJacobianIndex );
949
950     //Reset the iterators
951     unsigned int j = 0;
952     for ( j = 0; j < OutputDimension; j++ ) 
953       m_Iterator[j] = IteratorType( m_JacobianImage[j], m_SupportRegion);
954
955     // For each dimension, copy the weight to the support region
956     while ( ! (m_Iterator[0]).IsAtEnd() )
957       {
958         //copy weight to jacobian image
959         for ( j = 0; j < OutputDimension; j++ )
960           {
961             m_ZeroVector[j]=*weights;
962             (m_Iterator[j]).Set( m_ZeroVector);
963             m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
964             ++(m_Iterator[j]);
965           }
966         // go to next coefficient in the support region
967         weights++;
968       }
969     m_NeedResetJacobian = true;
970
971     // Return the result
972     return this->m_Jacobian;
973
974   }
975
976
977   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
978   void
979   BSplineDeformableTransform<TCoordRep, NInputDimensions,NOutputDimensions>
980   ::TransformPointToContinuousIndex( const InputPointType & point, ContinuousIndexType & index ) const
981   {
982     unsigned int j;
983
984     itk::Vector<double, OutputDimension> tvector;
985
986     for ( j = 0; j < OutputDimension; j++ )
987       {
988         tvector[j] = point[j] - this->m_GridOrigin[j];
989       }
990
991     itk::Vector<double, OutputDimension> cvector;
992
993     cvector = m_PointToIndex * tvector;
994
995     for ( j = 0; j < OutputDimension; j++ )
996       {
997         index[j] = static_cast< typename ContinuousIndexType::CoordRepType >( cvector[j] );
998       }
999   }
1000
1001   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
1002   unsigned
1003   BSplineDeformableTransform<TCoordRep, NInputDimensions,NOutputDimensions>
1004   ::SetJacobianImageData(JacobianPixelType * jacobianDataPointer, unsigned dim)
1005   {
1006     unsigned int numberOfPixels = m_GridRegion.GetNumberOfPixels();
1007     m_JacobianImage[dim]->GetPixelContainer()->SetImportPointer(jacobianDataPointer, numberOfPixels);
1008     return numberOfPixels;
1009   }
1010
1011   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
1012   void
1013   BSplineDeformableTransform<TCoordRep, NInputDimensions,NOutputDimensions>
1014   ::ResetJacobian() const
1015   {
1016     //========================================================
1017     // Zero all components of jacobian
1018     //========================================================
1019     // JV not thread safe (m_LastJacobianIndex), instantiate N transforms
1020     // NOTE: for efficiency, we only need to zero out the coefficients
1021     // that got fill last time this method was called.
1022
1023     unsigned int j = 0;
1024
1025     //Define the  region for each jacobian image
1026     m_SupportRegion.SetIndex(m_LastJacobianIndex);
1027
1028     //Initialize the iterators
1029     for (j = 0; j < OutputDimension; j++)
1030       m_Iterator[j] = IteratorType(m_JacobianImage[j], m_SupportRegion);
1031
1032     //Set the previously-set to zero
1033     while (!(m_Iterator[0]).IsAtEnd())
1034       {
1035         for (j = 0; j < OutputDimension; j++)
1036           {
1037             m_Iterator[j].Set(m_ZeroVector);
1038             ++(m_Iterator[j]);
1039           }
1040       }
1041     m_NeedResetJacobian = false;
1042   }
1043 } // namespace
1044
1045 #endif