]> Creatis software - clitk.git/blob - registration/clitkMultipleBSplineDeformableTransform.txx
Merge branch 'VTK6_Qt5_Overlay4D' into VTK6_Qt5_Binarize
[clitk.git] / registration / clitkMultipleBSplineDeformableTransform.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 __clitkMultipleBSplineDeformableTransform_txx
19 #define __clitkMultipleBSplineDeformableTransform_txx
20
21 // clitk
22 #include "clitkMultipleBSplineDeformableTransformInitializer.h"
23
24 // itk
25 #include <itkRelabelComponentImageFilter.h>
26
27 namespace clitk
28 {
29   // Constructor with default arguments
30   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
31   MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
32 #if ITK_VERSION_MAJOR >= 4
33   ::MultipleBSplineDeformableTransform() : Superclass(0)
34 #else
35   ::MultipleBSplineDeformableTransform() : Superclass(OutputDimension, 0)
36 #endif
37   {
38     m_nLabels = 1;
39     m_labels = 0;
40     m_labelInterpolator = 0;
41     m_trans.resize(1);
42     m_trans[0] = TransformType::New();
43     m_parameters.resize(1);
44
45     this->m_FixedParameters.SetSize(NInputDimensions * (NInputDimensions + 5) + 4);
46     this->m_FixedParameters.Fill(0.0);
47
48     m_InternalParametersBuffer = ParametersType(0);
49     // Make sure the parameters pointer is not NULL after construction.
50     m_InputParametersPointer = &m_InternalParametersBuffer;
51     m_BulkTransform = 0;
52     m_LastJacobian = -1;
53   }
54
55   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
56   MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
57   ::~MultipleBSplineDeformableTransform()
58   {
59   }
60
61   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
62   void
63   MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
64   ::SetLabels(ImageLabelPointer labels)
65   {
66     GetFixedParameters(); // Update m_FixedParameters
67     if (labels)
68     {
69       typedef itk::RelabelComponentImageFilter<ImageLabelType, ImageLabelType> RelabelImageFilterType;
70       typename RelabelImageFilterType::Pointer relabelImageFilter = RelabelImageFilterType::New();
71       relabelImageFilter->SetInput(labels);
72       relabelImageFilter->Update();
73       m_labels = relabelImageFilter->GetOutput();
74       m_nLabels = relabelImageFilter->GetNumberOfObjects();
75       m_trans.resize(m_nLabels);
76       m_parameters.resize(m_nLabels);
77       for (unsigned i = 0; i < m_nLabels; ++i)
78         m_trans[i] = TransformType::New();
79       m_labelInterpolator = ImageLabelInterpolator::New();
80       m_labelInterpolator->SetInputImage(m_labels);
81     }
82     else
83     {
84       m_nLabels = 1;
85       m_trans.resize(1);
86       m_parameters.resize(1);
87       m_trans[0] = TransformType::New();
88     }
89     this->m_FixedParameters[(5+NInputDimensions)*NInputDimensions + 2] = (double)((size_t) m_labels);
90     this->m_FixedParameters[(5+NInputDimensions)*NInputDimensions + 3] = m_nLabels;
91     SetFixedParameters(this->m_FixedParameters);
92   }
93
94 #define LOOP_ON_LABELS(FUNC, ARGS)          \
95   for (unsigned i = 0; i < m_nLabels; ++i)  \
96     m_trans[i]->FUNC(ARGS);
97
98   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
99   inline void
100   MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
101   ::SetSplineOrder(const unsigned int & splineOrder)
102   {
103     LOOP_ON_LABELS(SetSplineOrder, splineOrder);
104     this->Modified();
105   }
106
107   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
108   inline void
109   MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
110   ::SetSplineOrders(const SizeType & splineOrders)
111   {
112     LOOP_ON_LABELS(SetSplineOrders, splineOrders);
113     this->Modified();
114   }
115
116   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
117   inline void
118   MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
119   ::SetLUTSamplingFactor( const int & samplingFactor)
120   {
121     LOOP_ON_LABELS(SetLUTSamplingFactor, samplingFactor);
122     this->Modified();
123   }
124
125   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
126   inline void
127   MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
128   ::SetLUTSamplingFactors( const SizeType & samplingFactors)
129   {
130     LOOP_ON_LABELS(SetLUTSamplingFactors, samplingFactors);
131     this->Modified();
132   }
133
134   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
135   inline void
136   MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
137   ::SetGridRegion( const RegionType & region )
138   {
139     LOOP_ON_LABELS(SetGridRegion, region);
140     this->Modified();
141   }
142
143   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
144   inline void
145   MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
146   ::SetGridSpacing( const SpacingType & spacing )
147   {
148     LOOP_ON_LABELS(SetGridSpacing, spacing);
149     this->Modified();
150   }
151
152   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
153   inline void
154   MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
155   ::SetGridDirection( const DirectionType & direction )
156   {
157     LOOP_ON_LABELS(SetGridDirection, direction);
158     this->Modified();
159   }
160
161   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
162   inline void
163   MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
164   ::SetGridOrigin( const OriginType& origin )
165   {
166     LOOP_ON_LABELS(SetGridOrigin, origin);
167     this->Modified();
168   }
169
170   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
171   inline void
172   MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
173   ::SetIdentity()
174   {
175     LOOP_ON_LABELS(SetIdentity, );
176     if ( m_InputParametersPointer )
177     {
178       ParametersType * parameters =
179         const_cast<ParametersType *>( m_InputParametersPointer );
180       parameters->Fill( 0.0 );
181     }
182     this->Modified();
183   }
184
185   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
186   inline const std::vector<typename MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>::CoefficientImagePointer>&
187   MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
188   ::GetCoefficientImages() const
189   {
190     m_CoefficientImages.resize(m_nLabels);
191     for (unsigned i = 0; i < m_nLabels; ++i)
192       m_CoefficientImages[i] = m_trans[i]->GetCoefficientImage();
193     return m_CoefficientImages;
194   }
195
196   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
197   inline void
198   MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
199   ::SetCoefficientImages(std::vector<CoefficientImagePointer>& images)
200   {
201     if (images.size() != m_nLabels)
202     {
203       itkExceptionMacro( << "Mismatched between the number of labels "
204           << m_nLabels
205           << " and the number of coefficient images "
206           << images.size());
207     }
208     for (unsigned i = 0; i < m_nLabels; ++i)
209       m_trans[i]->SetCoefficientImage(images[i]);
210   }
211
212   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
213   void
214   MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
215   ::SetParameters( const ParametersType & parameters )
216   {
217     if ( parameters.Size() != this->GetNumberOfParameters() )
218       {
219         itkExceptionMacro(<<"Mismatched between parameters size "
220                           << parameters.size()
221                           << " and region size "
222                           << this->GetNumberOfParameters() );
223       }
224
225     // Clean up buffered parameters
226     m_InternalParametersBuffer = ParametersType( 0 );
227
228     // Keep a reference to the input parameters
229     m_InputParametersPointer = &parameters;
230
231     double * dataPointer = const_cast<double *>(m_InputParametersPointer->data_block());
232     unsigned tot = 0;
233     for (unsigned i = 0; i < m_nLabels; ++i)
234     {
235       unsigned int numberOfParameters = m_trans[i]->GetNumberOfParameters();
236       m_parameters[i].SetData(dataPointer + tot, numberOfParameters);
237       m_trans[i]->SetParameters(m_parameters[i]);
238       tot += numberOfParameters;
239     }
240     InitJacobian();
241
242     this->Modified();
243   }
244
245   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
246   void
247   MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
248   ::SetFixedParameters( const ParametersType & parameters )
249   {
250     // BSplineSeformableTransform parameters + m_labels pointer
251     if ( parameters.Size() != NInputDimensions * (5 + NInputDimensions) + 4)
252       {
253         itkExceptionMacro(<< "Mismatched between parameters size "
254                           << parameters.size()
255                           << " and number of fixed parameters "
256                           << NInputDimensions * (5 + NInputDimensions) + 4);
257       }
258     ImageLabelPointer tmp2 = ((ImageLabelPointer)( (size_t)parameters[(5+NInputDimensions)*NInputDimensions+2]));
259     if (tmp2 != m_labels)
260     {
261       if (tmp2)
262       {
263         m_labels = tmp2;
264         m_nLabels = parameters[(5+NInputDimensions) * NInputDimensions + 3 ];
265         m_trans.resize(m_nLabels);
266         m_parameters.resize(m_nLabels);
267         for (unsigned i = 0; i < m_nLabels; ++i)
268           m_trans[i] = TransformType::New();
269         m_labelInterpolator = ImageLabelInterpolator::New();
270         m_labelInterpolator->SetInputImage(m_labels);
271       }
272       else
273       {
274         m_labels = tmp2;
275         m_nLabels = 1;
276         m_trans.resize(1);
277         m_parameters.resize(1);
278         m_trans[0] = TransformType::New();
279       }
280     }
281     unsigned int numberOfParameters = NInputDimensions * (5 + NInputDimensions) + 2;
282     ParametersType tmp(numberOfParameters);
283     for (unsigned j = 0; j < numberOfParameters; ++j)
284       tmp.SetElement(j, parameters.GetElement(j));
285     LOOP_ON_LABELS(SetFixedParameters, tmp);
286     this->m_FixedParameters = parameters;
287     this->Modified();
288   }
289
290   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
291   void
292   MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
293   ::SetParametersByValue( const ParametersType & parameters )
294   {
295     if ( parameters.Size() != this->GetNumberOfParameters() )
296       {
297         itkExceptionMacro(<<"Mismatched between parameters size "
298                           << parameters.size()
299                           << " and region size "
300                           << this->GetNumberOfParameters() );
301       }
302
303     // Copy it
304     m_InternalParametersBuffer = parameters;
305     m_InputParametersPointer = &m_InternalParametersBuffer;
306
307     for (unsigned i = 0, tot = 0; i < m_nLabels; ++i)
308     {
309       unsigned int numberOfParameters = m_trans[i]->GetNumberOfParameters();
310       ParametersType tmp(numberOfParameters);
311       for (unsigned j = 0; j < numberOfParameters; ++j, ++tot)
312         tmp.SetElement(j, parameters.GetElement(tot));
313       m_trans[i]->SetParametersByValue(tmp);
314     }
315     InitJacobian();
316     this->Modified();
317   }
318
319
320   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
321   inline void
322   MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
323   ::SetBulkTransform(BulkTransformPointer b)
324   {
325     m_BulkTransform = b;
326     LOOP_ON_LABELS(SetBulkTransform, b);
327   }
328
329 #undef LOOP_ON_LABELS
330
331   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
332 #if ITK_VERSION_MAJOR >= 4
333   inline typename MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>::NumberOfParametersType
334 #else
335   inline unsigned int
336 #endif
337   MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
338   ::GetNumberOfParameters(void) const
339   {
340     unsigned int sum = 0;
341     for (unsigned i = 0; i < m_nLabels; ++i)
342       sum += m_trans[i]->GetNumberOfParameters();
343     return sum;
344   }
345
346   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
347   inline unsigned int
348   MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
349   ::GetNumberOfParametersPerDimension(void) const
350   {
351     unsigned int sum = 0;
352     for (unsigned i = 0; i < m_nLabels; ++i)
353       sum += m_trans[i]->GetNumberOfParametersPerDimension();
354     return sum;
355   }
356
357   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
358   inline const typename MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>::ParametersType &
359   MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
360   ::GetParameters() const
361   {
362     if (NULL == m_InputParametersPointer)
363       {
364         itkExceptionMacro( <<"Cannot GetParameters() because m_InputParametersPointer is NULL. Perhaps SetCoefficientImage() has been called causing the NULL pointer." );
365       }
366
367     return (*m_InputParametersPointer);
368   }
369
370   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
371   inline const typename MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>::ParametersType &
372   MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
373   ::GetFixedParameters() const
374   {
375     const ParametersType& trans0Para = m_trans[0]->GetFixedParameters();
376     for (unsigned i = 0; i < trans0Para.Size() ; ++i)
377       this->m_FixedParameters.SetElement(i, trans0Para.GetElement(i));
378     this->m_FixedParameters[(5+NInputDimensions)*NInputDimensions + 2] = (double)((size_t) m_labels);
379     this->m_FixedParameters[(5+NInputDimensions)*NInputDimensions + 3] = m_nLabels;
380     return this->m_FixedParameters;
381   }
382
383   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
384   inline void
385   MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
386   ::PrintSelf(std::ostream &os, itk::Indent indent) const
387   {
388     this->Superclass::PrintSelf(os, indent);
389     os << indent << "nLabels" << m_nLabels << std::endl;
390     itk::Indent ind = indent.GetNextIndent();
391
392     for (unsigned i = 0; i < m_nLabels; ++i)
393     {
394       os << indent << "Label " << i << std::endl;
395       m_trans[i]->Print(os, ind);
396     }
397   }
398
399   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
400   inline bool
401   MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
402   ::InsideValidRegion( const ContinuousIndexType& index ) const
403   {
404     bool res = false;
405     for (unsigned i = 0; i < m_nLabels; ++i)
406       res |= m_trans[i]->InsideValidRegion(index);
407     return res;
408   }
409
410   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
411   inline typename MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>::OutputPointType
412   MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
413   ::TransformPoint(const InputPointType &inputPoint) const
414   {
415     int lidx = 0;
416     if (m_labels)
417       lidx = m_labelInterpolator->Evaluate(inputPoint) - 1;
418     if (lidx == -1)
419       return inputPoint;
420     return m_trans[lidx]->TransformPoint(inputPoint);
421   }
422
423   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
424   inline typename MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>::OutputPointType
425   MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
426   ::DeformablyTransformPoint(const InputPointType &inputPoint) const
427   {
428     int lidx = 0;
429     if (m_labels)
430       lidx = m_labelInterpolator->Evaluate(inputPoint) - 1;
431     if (lidx == -1)
432       return inputPoint;
433     return m_trans[lidx]->DeformablyTransformPoint(inputPoint);
434   }
435
436 #if ITK_VERSION_MAJOR >= 4
437   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
438   inline void
439   MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
440   ::ComputeJacobianWithRespectToParameters (const InputPointType &point, JacobianType &jacobian) const
441   {
442     if (m_LastJacobian != -1)
443       m_trans[m_LastJacobian]->ResetJacobian();
444
445     int lidx = 0;
446     if (m_labels)
447       lidx = m_labelInterpolator->Evaluate(point) - 1;
448     if (lidx == -1)
449     {
450       m_LastJacobian = lidx;
451       jacobian = this->m_SharedDataBSplineJacobian;
452       return;
453     }
454
455     JacobianType unused;
456     m_trans[lidx]->ComputeJacobianWithRespectToParameters(point, unused);
457     m_LastJacobian = lidx;
458
459     jacobian = this->m_SharedDataBSplineJacobian;
460   }
461
462 #else
463   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
464   inline const typename MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>::JacobianType &
465   MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
466   ::GetJacobian( const InputPointType & point ) const
467   {
468     if (m_LastJacobian != -1)
469       m_trans[m_LastJacobian]->ResetJacobian();
470
471     int lidx = 0;
472     if (m_labels)
473       lidx = m_labelInterpolator->Evaluate(point) - 1;
474     if (lidx == -1)
475     {
476       m_LastJacobian = lidx;
477       return this->m_Jacobian;
478     }
479
480     m_trans[lidx]->GetJacobian(point);
481     m_LastJacobian = lidx;
482
483     return this->m_Jacobian;
484   }
485 #endif
486
487   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
488   inline void
489   MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions,NOutputDimensions>
490   ::TransformPointToContinuousIndex( const InputPointType & point, ContinuousIndexType & index ) const
491   {
492     m_trans[0]->TransformPointToContinuousIndex(point, index);
493   }
494
495   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
496   inline void
497   MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions,NOutputDimensions>::InitJacobian()
498   {
499     unsigned numberOfParameters = this->GetNumberOfParameters();
500 #if ITK_VERSION_MAJOR >= 4
501     this->m_SharedDataBSplineJacobian.set_size(OutputDimension, numberOfParameters);
502     JacobianPixelType * jacobianDataPointer = reinterpret_cast<JacobianPixelType *>(this->m_SharedDataBSplineJacobian.data_block());
503 #else
504     this->m_Jacobian.set_size(OutputDimension, numberOfParameters);
505     JacobianPixelType * jacobianDataPointer = reinterpret_cast<JacobianPixelType *>(this->m_Jacobian.data_block());
506 #endif
507     memset(jacobianDataPointer, 0,  numberOfParameters * sizeof (JacobianPixelType));
508
509     unsigned tot = 0;
510     for (unsigned int j = 0; j < OutputDimension; j++)
511     {
512       for (unsigned i = 0; i < m_nLabels; ++i)
513       {
514         unsigned numberOfPixels = m_trans[i]->SetJacobianImageData(jacobianDataPointer, j);
515         jacobianDataPointer += numberOfPixels;
516         tot += numberOfPixels;
517       }
518     }
519     assert(tot == numberOfParameters);
520   }
521
522 }  // namespace clitk
523
524 #endif // __clitkMultipleBSplineDeformableTransform_txx