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   inline unsigned int
 333   MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 334   ::GetNumberOfParameters(void) const
 335   {
 336     unsigned int sum = 0;
 337     for (unsigned i = 0; i < m_nLabels; ++i)
 338       sum += m_trans[i]->GetNumberOfParameters();
 339     return sum;
 340   }
 341
 342   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 343   inline unsigned int
 344   MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 345   ::GetNumberOfParametersPerDimension(void) const
 346   {
 347     unsigned int sum = 0;
 348     for (unsigned i = 0; i < m_nLabels; ++i)
 349       sum += m_trans[i]->GetNumberOfParametersPerDimension();
 350     return sum;
 351   }
 352
 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
 357   {
 358     if (NULL == m_InputParametersPointer)
 359       {
 360         itkExceptionMacro( <<"Cannot GetParameters() because m_InputParametersPointer is NULL. Perhaps SetCoefficientImage() has been called causing the NULL pointer." );
 361       }
 362
 363     return (*m_InputParametersPointer);
 364   }
 365
 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
 370   {
 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;
 377   }
 378
 379   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 380   inline void
 381   MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 382   ::PrintSelf(std::ostream &os, itk::Indent indent) const
 383   {
 384     this->Superclass::PrintSelf(os, indent);
 385     os << indent << "nLabels" << m_nLabels << std::endl;
 386     itk::Indent ind = indent.GetNextIndent();
 387
 388     for (unsigned i = 0; i < m_nLabels; ++i)
 389     {
 390       os << indent << "Label " << i << std::endl;
 391       m_trans[i]->Print(os, ind);
 392     }
 393   }
 394
 395   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 396   inline bool
 397   MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 398   ::InsideValidRegion( const ContinuousIndexType& index ) const
 399   {
 400     bool res = false;
 401     for (unsigned i = 0; i < m_nLabels; ++i)
 402       res |= m_trans[i]->InsideValidRegion(index);
 403     return res;
 404   }
 405
 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
 410   {
 411     int lidx = 0;
 412     if (m_labels)
 413       lidx = m_labelInterpolator->Evaluate(inputPoint) - 1;
 414     if (lidx == -1)
 415       return inputPoint;
 416     return m_trans[lidx]->TransformPoint(inputPoint);
 417   }
 418
 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
 423   {
 424     int lidx = 0;
 425     if (m_labels)
 426       lidx = m_labelInterpolator->Evaluate(inputPoint) - 1;
 427     if (lidx == -1)
 428       return inputPoint;
 429     return m_trans[lidx]->DeformablyTransformPoint(inputPoint);
 430   }
 431
 432 #if ITK_VERSION_MAJOR >= 4
 433   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 434   inline void
 435   MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 436   ::ComputeJacobianWithRespectToParameters (const InputPointType &point, JacobianType &jacobian) const
 437   {
 438     if (m_LastJacobian != -1)
 439       m_trans[m_LastJacobian]->ResetJacobian();
 440
 441     int lidx = 0;
 442     if (m_labels)
 443       lidx = m_labelInterpolator->Evaluate(point) - 1;
 444     if (lidx == -1)
 445     {
 446       m_LastJacobian = lidx;
 447       jacobian = this->m_SharedDataBSplineJacobian;
 448       return;
 449     }
 450
 451     JacobianType unused;
 452     m_trans[lidx]->ComputeJacobianWithRespectToParameters(point, unused);
 453     m_LastJacobian = lidx;
 454
 455     jacobian = this->m_SharedDataBSplineJacobian;
 456   }
 457
 458 #else
 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
 463   {
 464     if (m_LastJacobian != -1)
 465       m_trans[m_LastJacobian]->ResetJacobian();
 466
 467     int lidx = 0;
 468     if (m_labels)
 469       lidx = m_labelInterpolator->Evaluate(point) - 1;
 470     if (lidx == -1)
 471     {
 472       m_LastJacobian = lidx;
 473       return this->m_Jacobian;
 474     }
 475
 476     m_trans[lidx]->GetJacobian(point);
 477     m_LastJacobian = lidx;
 478
 479     return this->m_Jacobian;
 480   }
 481 #endif
 482
 483   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 484   inline void
 485   MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions,NOutputDimensions>
 486   ::TransformPointToContinuousIndex( const InputPointType & point, ContinuousIndexType & index ) const
 487   {
 488     m_trans[0]->TransformPointToContinuousIndex(point, index);
 489   }
 490
 491   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 492   inline void
 493   MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions,NOutputDimensions>::InitJacobian()
 494   {
 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());
 499 #else
 500     this->m_Jacobian.set_size(OutputDimension, numberOfParameters);
 501     JacobianPixelType * jacobianDataPointer = reinterpret_cast<JacobianPixelType *>(this->m_Jacobian.data_block());
 502 #endif
 503     memset(jacobianDataPointer, 0,  numberOfParameters * sizeof (JacobianPixelType));
 504
 505     unsigned tot = 0;
 506     for (unsigned int j = 0; j < OutputDimension; j++)
 507     {
 508       for (unsigned i = 0; i < m_nLabels; ++i)
 509       {
 510         unsigned numberOfPixels = m_trans[i]->SetJacobianImageData(jacobianDataPointer, j);
 511         jacobianDataPointer += numberOfPixels;
 512         tot += numberOfPixels;
 513       }
 514     }
 515     assert(tot == numberOfParameters);
 516   }
 517
 518 }  // namespace clitk
 519
 520 #endif // __clitkMultipleBSplineDeformableTransform_txx