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