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