]> Creatis software - clitk.git/blob - registration/clitkBSplineDeformableTransformInitializer.txx
remove debug DD
[clitk.git] / registration / clitkBSplineDeformableTransformInitializer.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 __clitkBSplineDeformableTransformInitializer_txx
19 #define __clitkBSplineDeformableTransformInitializer_txx
20 #include "clitkBSplineDeformableTransformInitializer.h"
21
22 namespace clitk
23 {
24
25
26   template < class TTransform, class TImage >
27   BSplineDeformableTransformInitializer<TTransform, TImage >
28   ::BSplineDeformableTransformInitializer()
29   {
30     this->m_NumberOfControlPointsInsideTheImage.Fill( 5 );
31     this->m_ControlPointSpacing.Fill(64.);
32     this->m_ChosenSpacing.Fill(64.);
33     m_NumberOfControlPointsIsGiven=false;
34     m_ControlPointSpacingIsGiven=false;
35     m_TransformRegionIsGiven=false;
36     m_SamplingFactorIsGiven=false;
37     m_SplineOrders.Fill(3);
38     m_Parameters=NULL;
39   }
40
41   template < class TTransform, class TImage >
42   void
43   BSplineDeformableTransformInitializer<TTransform, TImage >
44   ::InitializeTransform()
45   {
46     // Sanity check:
47     // The image is required for the region and spacing
48     if( ! this->m_Image )
49       {
50         itkExceptionMacro( "Reference Image has not been set" );
51         return;
52       }
53     
54     // A pointer to the transform is required
55     if( ! this->m_Transform )
56       {
57         itkExceptionMacro( "Transform has not been set" );
58         return;
59       }
60
61     // If the image come from a filter, then update that filter.
62     if( this->m_Image->GetSource() )
63       {
64         this->m_Image->GetSource()->Update();
65       }
66
67
68     //--------------------------------------
69     // Spacing & Size on image
70     //--------------------------------------
71     SpacingType fixedImageSpacing = m_Image->GetSpacing();
72     SizeType fixedImageSize=m_Image->GetLargestPossibleRegion().GetSize();
73     typename RegionType::SizeType   gridBorderSize;
74     typename RegionType::SizeType   totalGridSize;
75
76     // Only spacing given: adjust if necessary
77     if (m_ControlPointSpacingIsGiven  && !m_NumberOfControlPointsIsGiven)
78       {
79         for(unsigned int r=0; r<InputDimension; r++) 
80           {
81             // JV
82             m_ChosenSpacing[r]= m_ControlPointSpacing[r];
83             m_ControlPointSpacing[r]= ( round(m_ChosenSpacing[r]/fixedImageSpacing[r]) *fixedImageSpacing[r] ) ;
84             m_NumberOfControlPointsInsideTheImage[r] = ceil( (double)fixedImageSize[r]*fixedImageSpacing[r]/ m_ControlPointSpacing[r] );
85             if (  ( ceil( (double)fixedImageSize[r]*fixedImageSpacing[r]/ m_ControlPointSpacing[r] ) )
86                  == ( (double)fixedImageSize[r]*fixedImageSpacing[r]/ m_ControlPointSpacing[r] ) )
87               {
88                 m_NumberOfControlPointsInsideTheImage[r] +=1;
89                 DD("Adding control point");
90               }
91           }
92       }
93   
94     // Only number of CP given: adjust if necessary
95     // JV -1 ?
96     else if (!m_ControlPointSpacingIsGiven  && m_NumberOfControlPointsIsGiven)
97       {
98         for(unsigned int r=0; r<InputDimension; r++) 
99           {
100             m_ChosenSpacing[r]=fixedImageSpacing[r]*( (double)(fixedImageSize[r])  / 
101                                                       (double)(m_NumberOfControlPointsInsideTheImage[r]) );
102             m_ControlPointSpacing[r]= fixedImageSpacing[r]* ceil( (double)(fixedImageSize[r] - 1)  / 
103                                                                   (double)(m_NumberOfControlPointsInsideTheImage[r] - 1) );
104           }
105       }
106   
107     // Both or none of Spacing and number of CP given: no adjustment adjust, just warnings
108     else 
109       {
110         for(unsigned int r=0; r<InputDimension; r++) 
111           {
112             m_ChosenSpacing[r]= m_ControlPointSpacing[r];
113             if (m_NumberOfControlPointsInsideTheImage[r]*m_ControlPointSpacing[r]< fixedImageSize[r]*fixedImageSpacing[r]) 
114               {
115                 std::cout<<"WARNING: Specified control point region ("<<m_NumberOfControlPointsInsideTheImage[r]*m_ControlPointSpacing[r]
116                          <<"mm) does not cover the fixed image region ("<< fixedImageSize[r]*fixedImageSpacing[r]
117                          <<"mm) for dimension "<<r<<"!" <<std::endl
118                          <<"Specify only --spacing or --control for automatic adjustment..."<<std::endl;
119               }
120             if (  fmod(m_ControlPointSpacing[r], fixedImageSpacing[r]) ) 
121               {
122                 std::cout<<"WARNING: Specified control point spacing for dimension "<<r
123                          <<" does not allow exact representation of BLUT FFD!"<<std::endl
124                          <<"Spacing ratio is non-integer: "<<m_ControlPointSpacing[r]/ fixedImageSpacing[r]<<std::endl
125                          <<"Specify only --spacing or --control for automatic adjustment..."<<std::endl;  
126               }
127           }
128       }
129
130     if (m_Verbose) std::cout<<"The chosen control point spacing "<<m_ChosenSpacing<<"..."<<std::endl;
131     if (m_Verbose) std::cout<<"The control points spacing was adapted to "<<m_ControlPointSpacing<<"..."<<std::endl; 
132     if (m_Verbose) std::cout<<"The number of (internal) control points is "<<m_NumberOfControlPointsInsideTheImage<<"..."<<std::endl; 
133     
134
135     //--------------------------------------
136     // Border size 
137     //--------------------------------------
138     for(unsigned int r=0; r<InputDimension; r++) gridBorderSize[r]=m_SplineOrders[r]; // Border for spline order = 3 ( 1 lower, 2 upper )
139     totalGridSize = m_NumberOfControlPointsInsideTheImage + gridBorderSize;
140     if (m_Verbose) std::cout<<"The total control point grid size was set to "<<totalGridSize<<"..."<<std::endl;
141     RegionType gridRegion;
142     gridRegion.SetSize(totalGridSize);
143  
144   
145     //--------------------------------------
146     // Origin
147     //--------------------------------------
148     typedef typename TransformType::OriginType OriginType;
149     OriginType fixedImageOrigin, gridOrigin;
150
151     // In case of non-zero index
152     m_Image->TransformIndexToPhysicalPoint(m_Image->GetLargestPossibleRegion().GetIndex(),fixedImageOrigin);
153     typename ImageType::DirectionType gridDirection = m_Image->GetDirection();      
154     SizeType orderShift;
155
156     // Shift depends on order
157     for(unsigned int r=0; r<InputDimension; r++)
158       {
159         orderShift[r] = m_SplineOrders[r] / 2;
160         gridOrigin[r] = fixedImageOrigin[r]- (double) orderShift[r]* m_ControlPointSpacing[r];
161       } 
162     if (m_Verbose) std::cout<<"The control point grid origin was set to "<<gridOrigin<<"..."<<std::endl;
163
164
165     //--------------------------------------
166     // LUT sampling factors
167     //--------------------------------------
168     for (unsigned int i=0; i< InputDimension; i++)
169       {
170         if (m_Verbose) std::cout<<"For dimension "<<i<<", the ideal sampling factor (if integer) is a multitude of "
171                                 << m_ControlPointSpacing[i]/ fixedImageSpacing[i]<<"..."<<std::endl;
172         if ( !m_SamplingFactorIsGiven) m_SamplingFactors[i]= (int) ( m_ControlPointSpacing[i]/ fixedImageSpacing[i]);
173         if (m_Verbose) std::cout<<"Setting sampling factor "<<i<<" to "<<m_SamplingFactors[i]<<"..."<<std::endl;
174       }
175
176
177     //--------------------------------------
178     // Set
179     //--------------------------------------
180     this->m_Transform->SetSplineOrders(m_SplineOrders);
181     this->m_Transform->SetGridRegion( gridRegion );
182     this->m_Transform->SetGridOrigin( gridOrigin );
183     this->m_Transform->SetGridSpacing( m_ControlPointSpacing );
184     this->m_Transform->SetGridDirection( gridDirection );
185     this->m_Transform->SetLUTSamplingFactors(m_SamplingFactors);
186
187   }
188
189   template < class TTransform, class TImage >
190   void
191   BSplineDeformableTransformInitializer<TTransform, TImage >
192   ::SetInitialParameters( const std::string& s, ParametersType& params) 
193   {
194     //---------------------------------------
195     // Read Initial parameters
196     //---------------------------------------
197     typedef itk::ImageFileReader<CoefficientImageType> CoefficientReaderType;
198     typename CoefficientReaderType::Pointer coeffReader=CoefficientReaderType::New();
199     coeffReader->SetFileName(s);
200     if(m_Verbose) std::cout<<"Reading initial coefficients from file: "<<s<<"..."<<std::endl;
201     coeffReader->Update();
202     typename CoefficientImageType::Pointer coefficientImage=coeffReader->GetOutput();
203     this->SetInitialParameters(coefficientImage, params);
204   }
205   
206   template < class TTransform, class TImage >
207   void
208   BSplineDeformableTransformInitializer<TTransform, TImage >
209   ::SetInitialParameters(const typename CoefficientImageType::Pointer coefficientImage, ParametersType& params) 
210   {
211     // Keep a reference 
212     m_Parameters=&params;
213
214     // Resample
215     typedef ResampleBSplineDeformableTransformImageFilter<CoefficientImageType,CoefficientImageType> ResamplerType;
216     typename ResamplerType::Pointer resampler=ResamplerType::New();
217     resampler->SetInput(coefficientImage);      
218     resampler->SetOutputParametersFromImage(m_Transform->GetCoefficientImage()); 
219     if(m_Verbose) std::cout<<"Resampling initial coefficients..."<<std::endl;
220     resampler->Update();
221
222     // Copy parameters into the existing array, I know its crappy
223     typedef itk::ImageRegionConstIterator<CoefficientImageType> Iterator;
224     Iterator it (resampler->GetOutput(), resampler->GetOutput()->GetLargestPossibleRegion() );
225     it.GoToBegin();
226     unsigned int i=0;
227     while(! it.IsAtEnd())
228       {
229         for (unsigned int j=0; j<InputDimension;j++)
230           params[i+j]=it.Get()[j];
231
232         ++it;
233         i+=InputDimension;
234       }
235
236     // JV pass the reference to the array !!
237     this->m_Transform->SetParameters(params);
238   }
239
240
241   
242 }  // namespace itk
243
244 #endif