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