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