1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
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
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.
13 It is distributed under dual licence
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
22 #include "clitkMultipleBSplineDeformableTransformInitializer.h"
25 #include <itkRelabelComponentImageFilter.h>
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)
35 ::MultipleBSplineDeformableTransform() : Superclass(OutputDimension, 0)
40 m_labelInterpolator = 0;
42 m_trans[0] = TransformType::New();
43 m_parameters.resize(1);
45 this->m_FixedParameters.SetSize(NInputDimensions * (NInputDimensions + 5) + 4);
46 this->m_FixedParameters.Fill(0.0);
48 m_InternalParametersBuffer = ParametersType(0);
49 // Make sure the parameters pointer is not NULL after construction.
50 m_InputParametersPointer = &m_InternalParametersBuffer;
55 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
56 MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
57 ::~MultipleBSplineDeformableTransform()
61 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
63 MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
64 ::SetLabels(ImageLabelPointer labels)
66 GetFixedParameters(); // Update m_FixedParameters
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);
86 m_parameters.resize(1);
87 m_trans[0] = TransformType::New();
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);
94 #define LOOP_ON_LABELS(FUNC, ARGS) \
95 for (unsigned i = 0; i < m_nLabels; ++i) \
96 m_trans[i]->FUNC(ARGS);
98 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
100 MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
101 ::SetSplineOrder(const unsigned int & splineOrder)
103 LOOP_ON_LABELS(SetSplineOrder, splineOrder);
107 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
109 MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
110 ::SetSplineOrders(const SizeType & splineOrders)
112 LOOP_ON_LABELS(SetSplineOrders, splineOrders);
116 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
118 MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
119 ::SetLUTSamplingFactor( const int & samplingFactor)
121 LOOP_ON_LABELS(SetLUTSamplingFactor, samplingFactor);
125 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
127 MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
128 ::SetLUTSamplingFactors( const SizeType & samplingFactors)
130 LOOP_ON_LABELS(SetLUTSamplingFactors, samplingFactors);
134 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
136 MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
137 ::SetGridRegion( const RegionType & region )
139 LOOP_ON_LABELS(SetGridRegion, region);
143 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
145 MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
146 ::SetGridSpacing( const SpacingType & spacing )
148 LOOP_ON_LABELS(SetGridSpacing, spacing);
152 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
154 MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
155 ::SetGridDirection( const DirectionType & direction )
157 LOOP_ON_LABELS(SetGridDirection, direction);
161 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
163 MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
164 ::SetGridOrigin( const OriginType& origin )
166 LOOP_ON_LABELS(SetGridOrigin, origin);
170 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
172 MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
175 LOOP_ON_LABELS(SetIdentity, );
176 if ( m_InputParametersPointer )
178 ParametersType * parameters =
179 const_cast<ParametersType *>( m_InputParametersPointer );
180 parameters->Fill( 0.0 );
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
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;
196 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
198 MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
199 ::SetCoefficientImages(std::vector<CoefficientImagePointer>& images)
201 if (images.size() != m_nLabels)
203 itkExceptionMacro( << "Mismatched between the number of labels "
205 << " and the number of coefficient images "
208 for (unsigned i = 0; i < m_nLabels; ++i)
209 m_trans[i]->SetCoefficientImage(images[i]);
212 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
214 MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
215 ::SetParameters( const ParametersType & parameters )
217 if ( parameters.Size() != this->GetNumberOfParameters() )
219 itkExceptionMacro(<<"Mismatched between parameters size "
221 << " and region size "
222 << this->GetNumberOfParameters() );
225 // Clean up buffered parameters
226 m_InternalParametersBuffer = ParametersType( 0 );
228 // Keep a reference to the input parameters
229 m_InputParametersPointer = ¶meters;
231 double * dataPointer = const_cast<double *>(m_InputParametersPointer->data_block());
233 for (unsigned i = 0; i < m_nLabels; ++i)
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;
245 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
247 MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
248 ::SetFixedParameters( const ParametersType & parameters )
250 // BSplineSeformableTransform parameters + m_labels pointer
251 if ( parameters.Size() != NInputDimensions * (5 + NInputDimensions) + 4)
253 itkExceptionMacro(<< "Mismatched between parameters size "
255 << " and number of fixed parameters "
256 << NInputDimensions * (5 + NInputDimensions) + 4);
258 ImageLabelPointer tmp2 = ((ImageLabelPointer)( (size_t)parameters[(5+NInputDimensions)*NInputDimensions+2]));
259 if (tmp2 != m_labels)
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);
277 m_parameters.resize(1);
278 m_trans[0] = TransformType::New();
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;
290 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
292 MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
293 ::SetParametersByValue( const ParametersType & parameters )
295 if ( parameters.Size() != this->GetNumberOfParameters() )
297 itkExceptionMacro(<<"Mismatched between parameters size "
299 << " and region size "
300 << this->GetNumberOfParameters() );
304 m_InternalParametersBuffer = parameters;
305 m_InputParametersPointer = &m_InternalParametersBuffer;
307 for (unsigned i = 0, tot = 0; i < m_nLabels; ++i)
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);
320 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
322 MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
323 ::SetBulkTransform(BulkTransformPointer b)
326 LOOP_ON_LABELS(SetBulkTransform, b);
329 #undef LOOP_ON_LABELS
331 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
333 MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
334 ::GetNumberOfParameters(void) const
336 unsigned int sum = 0;
337 for (unsigned i = 0; i < m_nLabels; ++i)
338 sum += m_trans[i]->GetNumberOfParameters();
342 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
344 MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
345 ::GetNumberOfParametersPerDimension(void) const
347 unsigned int sum = 0;
348 for (unsigned i = 0; i < m_nLabels; ++i)
349 sum += m_trans[i]->GetNumberOfParametersPerDimension();
353 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
354 inline const typename MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>::ParametersType &
355 MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
356 ::GetParameters() const
358 if (NULL == m_InputParametersPointer)
360 itkExceptionMacro( <<"Cannot GetParameters() because m_InputParametersPointer is NULL. Perhaps SetCoefficientImage() has been called causing the NULL pointer." );
363 return (*m_InputParametersPointer);
366 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
367 inline const typename MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>::ParametersType &
368 MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
369 ::GetFixedParameters() const
371 const ParametersType& trans0Para = m_trans[0]->GetFixedParameters();
372 for (unsigned i = 0; i < trans0Para.Size() ; ++i)
373 this->m_FixedParameters.SetElement(i, trans0Para.GetElement(i));
374 this->m_FixedParameters[(5+NInputDimensions)*NInputDimensions + 2] = (double)((size_t) m_labels);
375 this->m_FixedParameters[(5+NInputDimensions)*NInputDimensions + 3] = m_nLabels;
376 return this->m_FixedParameters;
379 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
381 MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
382 ::PrintSelf(std::ostream &os, itk::Indent indent) const
384 this->Superclass::PrintSelf(os, indent);
385 os << indent << "nLabels" << m_nLabels << std::endl;
386 itk::Indent ind = indent.GetNextIndent();
388 for (unsigned i = 0; i < m_nLabels; ++i)
390 os << indent << "Label " << i << std::endl;
391 m_trans[i]->Print(os, ind);
395 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
397 MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
398 ::InsideValidRegion( const ContinuousIndexType& index ) const
401 for (unsigned i = 0; i < m_nLabels; ++i)
402 res |= m_trans[i]->InsideValidRegion(index);
406 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
407 inline typename MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>::OutputPointType
408 MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
409 ::TransformPoint(const InputPointType &inputPoint) const
413 lidx = m_labelInterpolator->Evaluate(inputPoint) - 1;
416 return m_trans[lidx]->TransformPoint(inputPoint);
419 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
420 inline typename MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>::OutputPointType
421 MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
422 ::DeformablyTransformPoint(const InputPointType &inputPoint) const
426 lidx = m_labelInterpolator->Evaluate(inputPoint) - 1;
429 return m_trans[lidx]->DeformablyTransformPoint(inputPoint);
432 #if ITK_VERSION_MAJOR >= 4
433 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
435 MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
436 ::ComputeJacobianWithRespectToParameters (const InputPointType &point, JacobianType &jacobian) const
438 if (m_LastJacobian != -1)
439 m_trans[m_LastJacobian]->ResetJacobian();
443 lidx = m_labelInterpolator->Evaluate(point) - 1;
446 m_LastJacobian = lidx;
447 jacobian = this->m_SharedDataBSplineJacobian;
452 m_trans[lidx]->ComputeJacobianWithRespectToParameters(point, unused);
453 m_LastJacobian = lidx;
455 jacobian = this->m_SharedDataBSplineJacobian;
459 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
460 inline const typename MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>::JacobianType &
461 MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
462 ::GetJacobian( const InputPointType & point ) const
464 if (m_LastJacobian != -1)
465 m_trans[m_LastJacobian]->ResetJacobian();
469 lidx = m_labelInterpolator->Evaluate(point) - 1;
472 m_LastJacobian = lidx;
473 return this->m_Jacobian;
476 m_trans[lidx]->GetJacobian(point);
477 m_LastJacobian = lidx;
479 return this->m_Jacobian;
483 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
485 MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions,NOutputDimensions>
486 ::TransformPointToContinuousIndex( const InputPointType & point, ContinuousIndexType & index ) const
488 m_trans[0]->TransformPointToContinuousIndex(point, index);
491 template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
493 MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions,NOutputDimensions>::InitJacobian()
495 unsigned numberOfParameters = this->GetNumberOfParameters();
496 #if ITK_VERSION_MAJOR >= 4
497 this->m_SharedDataBSplineJacobian.set_size(OutputDimension, numberOfParameters);
498 JacobianPixelType * jacobianDataPointer = reinterpret_cast<JacobianPixelType *>(this->m_SharedDataBSplineJacobian.data_block());
500 this->m_Jacobian.set_size(OutputDimension, numberOfParameters);
501 JacobianPixelType * jacobianDataPointer = reinterpret_cast<JacobianPixelType *>(this->m_Jacobian.data_block());
503 memset(jacobianDataPointer, 0, numberOfParameters * sizeof (JacobianPixelType));
506 for (unsigned int j = 0; j < OutputDimension; j++)
508 for (unsigned i = 0; i < m_nLabels; ++i)
510 unsigned numberOfPixels = m_trans[i]->SetJacobianImageData(jacobianDataPointer, j);
511 jacobianDataPointer += numberOfPixels;
512 tot += numberOfPixels;
515 assert(tot == numberOfParameters);
520 #endif // __clitkMultipleBSplineDeformableTransform_txx