]> Creatis software - clitk.git/blob - registration/clitkShapedBLUTSpatioTemporalDeformableTransformInitializer.txx
Moved from repository clitk to clitk.private/tests_dav
[clitk.git] / registration / clitkShapedBLUTSpatioTemporalDeformableTransformInitializer.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 __clitkShapedBLUTSpatioTemporalDeformableTransformInitializer_txx
19 #define __clitkShapedBLUTSpatioTemporalDeformableTransformInitializer_txx
20 #include "clitkShapedBLUTSpatioTemporalDeformableTransformInitializer.h"
21
22 namespace clitk
23 {
24
25
26   template < class TTransform, class TImage >
27   ShapedBLUTSpatioTemporalDeformableTransformInitializer<TTransform, TImage >
28   ::ShapedBLUTSpatioTemporalDeformableTransformInitializer()
29   {
30     this->m_NumberOfControlPointsInsideTheImage.Fill( 5 );
31     this->m_ControlPointSpacing.Fill(64.);
32     m_ControlPointSpacing[InputDimension-1]=2;
33     this->m_ChosenSpacing.Fill(64.);
34     m_ChosenSpacing[InputDimension-1]=2;
35     m_ControlPointSpacing[InputDimension-1]=2;
36     m_NumberOfControlPointsIsGiven=false;
37     m_ControlPointSpacingIsGiven=false;
38     m_TransformRegionIsGiven=false;
39     m_SamplingFactorIsGiven=false;
40     m_SplineOrders.Fill(3);
41     m_BC1=true;
42     m_BC2=true;
43     m_TrajectoryShape=0;
44   }
45
46   template < class TTransform, class TImage >
47   void
48   ShapedBLUTSpatioTemporalDeformableTransformInitializer<TTransform, TImage >
49   ::InitializeTransform()
50   {
51     // Sanity check:
52     // The image is required for the region and spacing
53     if( ! this->m_Image )
54       {
55         itkExceptionMacro( "Reference Image has not been set" );
56         return;
57       }
58     
59     // A pointer to the transform is required
60     if( ! this->m_Transform )
61       {
62         itkExceptionMacro( "Transform has not been set" );
63         return;
64       }
65
66     // If the image come from a filter, then update that filter.
67     if( this->m_Image->GetSource() )
68       {
69         this->m_Image->GetSource()->Update();
70       }
71
72
73     //--------------------------------------
74     // Spacing & Size on image
75     //--------------------------------------
76     SpacingType fixedImageSpacing = m_Image->GetSpacing();
77     SizeType fixedImageSize=m_Image->GetLargestPossibleRegion().GetSize();
78     typename RegionType::SizeType   gridBorderSize;
79     typename RegionType::SizeType   totalGridSize;
80
81     // Only spacing given: adjust if necessary
82     if (m_ControlPointSpacingIsGiven  && !m_NumberOfControlPointsIsGiven)
83       {
84         for(unsigned int r=0; r<InputDimension; r++) 
85           {
86             // JV
87             m_ChosenSpacing[r]= m_ControlPointSpacing[r];
88             m_ControlPointSpacing[r]= ( itk::Math::Round<double>(m_ChosenSpacing[r]/fixedImageSpacing[r]) *fixedImageSpacing[r] ) ;
89             m_NumberOfControlPointsInsideTheImage[r] = ceil( (double)fixedImageSize[r]*fixedImageSpacing[r]/ m_ControlPointSpacing[r] );
90             if ( (  ( ceil( (double)fixedImageSize[r]*fixedImageSpacing[r]/ m_ControlPointSpacing[r] ) )
91                  == ( (double)fixedImageSize[r]*fixedImageSpacing[r]/ m_ControlPointSpacing[r] ) )
92                  &&  (r!= InputDimension-1)  )
93               {
94                 m_NumberOfControlPointsInsideTheImage[r] +=1;
95                 DD("Adding control point");
96               }
97           }
98       }
99   
100     // Only number of CP given: adjust if necessary
101     // JV -1 ?
102     else if (!m_ControlPointSpacingIsGiven  && m_NumberOfControlPointsIsGiven)
103       {
104         for(unsigned int r=0; r<InputDimension; r++) 
105           {
106             m_ChosenSpacing[r]=fixedImageSpacing[r]*( (double)(fixedImageSize[r])  / 
107                                                       (double)(m_NumberOfControlPointsInsideTheImage[r]) );
108             m_ControlPointSpacing[r]= fixedImageSpacing[r]* ceil( (double)(fixedImageSize[r] - 1)  / 
109                                                                   (double)(m_NumberOfControlPointsInsideTheImage[r] - 1) );
110           }
111       }
112   
113     // Both or none of Spacing and number of CP given: no adjustment adjust, just warnings
114     else 
115       {
116         for(unsigned int r=0; r<InputDimension; r++) 
117           {
118             m_ChosenSpacing[r]= m_ControlPointSpacing[r];
119             if (m_NumberOfControlPointsInsideTheImage[r]*m_ControlPointSpacing[r]< fixedImageSize[r]*fixedImageSpacing[r]) 
120               {
121                 std::cout<<"WARNING: Specified control point region ("<<m_NumberOfControlPointsInsideTheImage[r]*m_ControlPointSpacing[r]
122                          <<"mm) does not cover the fixed image region ("<< fixedImageSize[r]*fixedImageSpacing[r]
123                          <<"mm) for dimension "<<r<<"!" <<std::endl
124                          <<"Specify only --spacing or --control for automatic adjustment..."<<std::endl;
125               }
126             if (  fmod(m_ControlPointSpacing[r], fixedImageSpacing[r]) ) 
127               {
128                 std::cout<<"WARNING: Specified control point spacing for dimension "<<r
129                          <<" does not allow exact representation of BLUT FFD!"<<std::endl
130                          <<"Spacing ratio is non-integer: "<<m_ControlPointSpacing[r]/ fixedImageSpacing[r]<<std::endl
131                          <<"Specify only --spacing or --control for automatic adjustment..."<<std::endl;  
132               }
133           }
134       }
135
136     //--------------------------------------
137     // LUT sampling factors
138     //--------------------------------------
139     for (unsigned int i=0; i< InputDimension; i++)
140       {
141         if (m_Verbose) std::cout<<"For dimension "<<i<<", the ideal sampling factor (if integer) is a multitude of "
142                                 << m_ControlPointSpacing[i]/ fixedImageSpacing[i]<<"..."<<std::endl;
143         if ( !m_SamplingFactorIsGiven) m_SamplingFactors[i]= (int) ( m_ControlPointSpacing[i]/ fixedImageSpacing[i]);
144       }
145
146     //--------------------------------------
147     // Set the transform properties 
148     //--------------------------------------
149     RegionType gridRegion;
150     OriginType gridOrigin;
151     DirectionType gridDirection;
152
153     //--------------------------------------
154     // Border size 
155     //--------------------------------------
156     for(unsigned int r=0; r<InputDimension; r++) gridBorderSize[r]=m_SplineOrders[r]; // Border for spline order = 3 ( 1 lower, 2 upper )
157     totalGridSize = m_NumberOfControlPointsInsideTheImage + gridBorderSize;
158
159     //--------------------------------------
160     // Origin
161     //--------------------------------------
162     OriginType fixedImageOrigin;
163
164     // In case of non-zero index
165     m_Image->TransformIndexToPhysicalPoint(m_Image->GetLargestPossibleRegion().GetIndex(),fixedImageOrigin);
166     gridDirection = m_Image->GetDirection();      
167     itk::FixedArray<int, InputDimension> orderShift;
168     
169     // Spacing is 2.5 : manually modify the props
170     if( !(m_Transform->GetTransformShape()%2) )
171       {
172         m_ControlPointSpacing[InputDimension-1]=2.5;
173         m_SamplingFactors[InputDimension-1]=5;
174         totalGridSize[InputDimension-1]-=1;
175         gridRegion.SetSize(totalGridSize);
176       }
177
178
179     switch(m_Transform->GetTransformShape()){
180     
181       // The egg shape
182     case 0:
183     case 1:  
184       {
185         if (m_Verbose) std::cout<<"Using the egg shape..."<<std::endl;
186
187         //--------------------------------------
188         // Border size 
189         //--------------------------------------
190         // Boundary condition1: cyclic
191         if (m_BC1)
192           totalGridSize[InputDimension-1]-=3;
193         
194         // Boundary condition1: reference is on trajectory
195         if (m_BC2)
196           totalGridSize[InputDimension-1]-=1;
197         
198         gridRegion.SetSize(totalGridSize);
199                 
200         //--------------------------------------
201         // Origin
202         //--------------------------------------
203         // Shift depends on order
204         for(unsigned int r=0; r<InputDimension; r++)
205           {
206             orderShift[r] = m_SplineOrders[r] / 2;
207             if (r==InputDimension-1)
208               {
209                 // Boundary condition1: cyclic
210                 if (m_BC1)
211                   orderShift[r]-=1;
212                 
213                 // Boundary condition1: reference is on trajectory
214                 if (m_BC2)
215                   orderShift[r]-=1;
216               }
217             
218             gridOrigin[r] = fixedImageOrigin[r]- (double) orderShift[r]* m_ControlPointSpacing[r];
219           } 
220         break;
221       }
222
223       // The rabit: a diamond with fixed extreme values
224     case 2:
225     case 3:
226       {
227
228         if (m_Verbose) std::cout<<"Using the rabbit shape..."<<std::endl;
229         
230         //--------------------------------------
231         // Border size 
232         //--------------------------------------
233         // Boundary condition1: cyclic
234         if (m_BC1)
235           totalGridSize[InputDimension-1]-=1;
236         
237         // Boundary condition1: reference is on trajectory
238         if (m_BC2)
239           totalGridSize[InputDimension-1]-=1;
240
241         // Extreme values are constrained
242         totalGridSize[InputDimension-1]-=2;
243         
244         gridRegion.SetSize(totalGridSize);
245         
246         //--------------------------------------
247         // Origin
248         //--------------------------------------
249         // Shift depends on order
250         for(unsigned int r=0; r<InputDimension; r++)
251           {
252             orderShift[r] = m_SplineOrders[r] / 2;
253             if (r==InputDimension-1)
254               {
255                 // Boundary condition1: cyclic
256                 orderShift[r]=0;
257               }
258             gridOrigin[r] = fixedImageOrigin[r]- (double) orderShift[r]* m_ControlPointSpacing[r];
259           } 
260
261         break;
262       }
263
264       // The sputnik: a diamond with one indepent extreme value and 4DOF
265     case 4:
266     case 5:
267       {
268         if (m_Verbose) std::cout<<"Using the sputnik shape..."<<std::endl;
269
270         //--------------------------------------
271         // Border size 
272         //--------------------------------------
273         // Boundary condition1: cyclic
274         if (m_BC1)
275           totalGridSize[InputDimension-1]-=1;
276         
277         // Boundary condition1: reference is on trajectory
278         if (m_BC2)
279           totalGridSize[InputDimension-1]-=1;
280
281         // One extreme value is fixed
282         totalGridSize[InputDimension-1]-=1;
283         gridRegion.SetSize(totalGridSize);
284         
285         //--------------------------------------
286         // Origin
287         //--------------------------------------
288         // Shift depends on order
289         for(unsigned int r=0; r<InputDimension; r++)
290           {
291             orderShift[r] = m_SplineOrders[r] / 2;
292             if (r==InputDimension-1)
293               {
294                 // Boundary condition1: cyclic
295                 orderShift[r]=0;
296               }
297             gridOrigin[r] = fixedImageOrigin[r]- (double) orderShift[r]* m_ControlPointSpacing[r];
298           } 
299         break;
300       }
301
302       // The diamond
303     case 6:  
304     case 7:
305       {
306         if (m_Verbose) std::cout<<"Using the diamond shape..."<<std::endl;
307         
308         //--------------------------------------
309         // Border size 
310         //--------------------------------------
311         // Boundary condition1: cyclic
312         if (m_BC1)
313           totalGridSize[InputDimension-1]-=1;
314         
315         // Boundary condition1: reference is on trajectory
316         if (m_BC2)
317           totalGridSize[InputDimension-1]-=1;
318         
319         gridRegion.SetSize(totalGridSize);
320         
321         //--------------------------------------
322         // Origin
323         //--------------------------------------
324         // Shift depends on order
325         for(unsigned int r=0; r<InputDimension; r++)
326           {
327             orderShift[r] = m_SplineOrders[r] / 2;
328             gridOrigin[r] = fixedImageOrigin[r]- (double) orderShift[r]* m_ControlPointSpacing[r];
329           } 
330
331         break;
332
333       } //end case 7
334
335    // The sputnik: a diamond with one indepent extreme value and 4DOF
336     case 8:
337     case 9:
338       {
339         if (m_Verbose) std::cout<<"Using the sputnik shape..."<<std::endl;
340
341         //--------------------------------------
342         // Border size 
343         //--------------------------------------
344         // Boundary condition1: cyclic
345         if (m_BC1)
346           totalGridSize[InputDimension-1]-=1;
347         
348         // Boundary condition1: reference is on trajectory
349         if (m_BC2)
350           totalGridSize[InputDimension-1]-=1;
351
352         // One extreme value is fixed
353         totalGridSize[InputDimension-1]-=1;
354         gridRegion.SetSize(totalGridSize);
355         
356         //--------------------------------------
357         // Origin
358         //--------------------------------------
359         // Shift depends on order
360         for(unsigned int r=0; r<InputDimension; r++)
361           {
362             orderShift[r] = m_SplineOrders[r] / 2;
363             if (r==InputDimension-1)
364               {
365                 // Boundary condition1: cyclic
366                 orderShift[r]=0;
367               }
368             gridOrigin[r] = fixedImageOrigin[r]- (double) orderShift[r]* m_ControlPointSpacing[r];
369           } 
370         break;
371       }
372   
373     default:
374       {
375         DD("Shape not available");
376       }
377
378     }
379
380
381
382     if (m_Verbose) std::cout<<"The chosen control point spacing "<<m_ChosenSpacing<<"..."<<std::endl;
383     if (m_Verbose) std::cout<<"The control points spacing was adapted to "<<m_ControlPointSpacing<<"..."<<std::endl; 
384     if (m_Verbose) std::cout<<"The number of (internal) control points is "<<m_NumberOfControlPointsInsideTheImage<<"..."<<std::endl;
385     if (m_Verbose) std::cout<<"Setting sampling factors to "<<m_SamplingFactors<<"..."<<std::endl; 
386     if (m_Verbose) std::cout<<"The total control point grid size was set to "<<totalGridSize<<"..."<<std::endl;
387     if (m_Verbose) std::cout<<"The control point grid origin was set to "<<gridOrigin<<"..."<<std::endl;
388     
389     //--------------------------------------
390     // Set
391     //--------------------------------------
392     this->m_Transform->SetSplineOrders(m_SplineOrders);
393     this->m_Transform->SetGridSpacing( m_ControlPointSpacing );
394     this->m_Transform->SetGridOrigin( gridOrigin );
395     this->m_Transform->SetGridDirection( gridDirection );
396     this->m_Transform->SetGridRegion( gridRegion );
397     this->m_Transform->SetLUTSamplingFactors(m_SamplingFactors);
398
399   }
400
401   template < class TTransform, class TImage >
402   void
403   ShapedBLUTSpatioTemporalDeformableTransformInitializer<TTransform, TImage >
404   ::SetInitialParameters( const std::string& s, ParametersType& params) 
405   {
406     //---------------------------------------
407     // Read Initial parameters
408     //---------------------------------------
409     typedef itk::ImageFileReader<CoefficientImageType> CoefficientReaderType;
410     typename CoefficientReaderType::Pointer coeffReader=CoefficientReaderType::New();
411     coeffReader->SetFileName(s);
412     if(m_Verbose) std::cout<<"Reading initial coefficients from file: "<<s<<"..."<<std::endl;
413     coeffReader->Update();
414     typename CoefficientImageType::Pointer coefficientImage=coeffReader->GetOutput();
415     this->SetInitialParameters(coefficientImage, params);
416   }
417   
418   template < class TTransform, class TImage >
419   void
420   ShapedBLUTSpatioTemporalDeformableTransformInitializer<TTransform, TImage >
421   ::SetInitialParameters(const typename CoefficientImageType::Pointer coefficientImage, ParametersType& params) 
422   {
423     // Keep a reference 
424     m_Parameters=&params;
425
426     // Resample
427     typedef ResampleBSplineDeformableTransformImageFilter<CoefficientImageType,CoefficientImageType> ResamplerType;
428     typename ResamplerType::Pointer resampler=ResamplerType::New();
429     resampler->SetVerbose(m_Verbose);
430     resampler->SetInput(coefficientImage);      
431     resampler->SetOutputParametersFromImage(m_Transform->GetCoefficientImage()); 
432     if(m_Verbose) std::cout<<"Resampling initial coefficients..."<<std::endl;
433     resampler->Update();
434
435     // Copy parameters into the existing array, I know its crappy
436     typedef itk::ImageRegionConstIterator<CoefficientImageType> Iterator;
437     Iterator it (resampler->GetOutput(), resampler->GetOutput()->GetLargestPossibleRegion() );
438     it.GoToBegin();
439     unsigned int i=0;
440     while(! it.IsAtEnd())
441       {
442         for (unsigned int j=0; j<InputDimension-1;j++)
443           params[i+j]=it.Get()[j];
444
445         ++it;
446         i+=InputDimension-1;
447       }
448
449     // JV pass the reference to the array !!
450     // JV Update the padded version !!
451     this->m_Transform->SetParameters(params);
452   }
453
454
455   //   template < class TTransform, class TImage >
456   //   void
457   //   ShapedBLUTSpatioTemporalDeformableTransformInitializer<TTransform, TImage >
458   //   ::SetInitialPaddedParameters( const std::string& s, ParametersType& params) 
459   //   {
460   //     //---------------------------------------
461   //     // Read Initial parameters
462   //     //---------------------------------------
463   //     typedef itk::ImageFileReader<CoefficientImageType> CoefficientReaderType;
464   //     typename CoefficientReaderType::Pointer coeffReader=CoefficientReaderType::New();
465   //     coeffReader->SetFileName(s);
466   //     if(m_Verbose) std::cout<<"Reading initial padded coefficients from file: "<<s<<"..."<<std::endl;
467   //     coeffReader->Update();
468   //     typename CoefficientImageType::Pointer coefficientImage=coeffReader->GetOutput();
469   //     this->SetInitialPaddedParameters(coefficientImage, params);
470   //   }
471   
472   //   template < class TTransform, class TImage >
473   //   void
474   //   ShapedBLUTSpatioTemporalDeformableTransformInitializer<TTransform, TImage >
475   //   ::SetInitialPaddedParameters(const typename CoefficientImageType::Pointer coefficientImage, ParametersType& params) 
476   //   {
477   //     // Resample
478   //     typedef ResampleBSplineSpatioTemporalDeformableTransformImageFilter<CoefficientImageType,CoefficientImageType> ResamplerType;
479   //     typename ResamplerType::Pointer resampler=ResamplerType::New();
480   //     resampler->SetInput(coefficientImage); 
481   //     resampler->SetOutputParametersFromImage(m_Transform->GetCoefficientImage()); 
482   //     if(m_Verbose) std::cout<<"Resampling initial padded coefficients..."<<std::endl;
483   //     resampler->Update();
484     
485   //     // Copy parameters into the existing array, I know its crappy
486   //     typedef itk::ImageRegionConstIterator<CoefficientImageType> Iterator;
487   //     Iterator it (resampler->GetOutput(), resampler->GetOutput()->GetLargestPossibleRegion() );
488   //     it.GoToBegin();
489   //     unsigned int i=0;
490   //     while(! it.IsAtEnd())
491   //       {
492   //    for (unsigned int j=0; j<InputDimension;j++)
493   //      params[i+j]=it.Get()[j];
494
495   //    ++it;
496   //    i+=InputDimension;
497   //       }
498   //   }
499
500   
501 }  // namespace itk
502
503 #endif