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