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