registration/clitkBSplineDeformableTransform.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 __clitkBSplineDeformableTransform_txx
  19 #define __clitkBSplineDeformableTransform_tx
  20 #include "clitkBSplineDeformableTransform.h"
  21
  22 //itk
  23 #include "itkContinuousIndex.h"
  24 #include "itkImageRegionIterator.h"
  25 #include "itkImageRegionConstIterator.h"
  26 #include "itkImageRegionConstIteratorWithIndex.h"
  27 #include "itkIdentityTransform.h"
  28
  29 namespace clitk
  30 {
  31
  32   // Constructor with default arguments
  33   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
  34   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
  35   ::BSplineDeformableTransform():Superclass(OutputDimension,0)
  36   {
  37     unsigned int i;
  38
  39     //JV default spline order
  40     for ( i=0;i<InputDimension; i++)
  41       m_SplineOrders[i]=3;
  42
  43     //JV default LUTSamplingfactor
  44     for ( i=0;i<InputDimension; i++)
  45       m_LUTSamplingFactors[i]=20;
  46
  47     for (i=0;i<InputDimension; i++)
  48       m_SupportSize[i] = m_SplineOrders[i]+1;
  49
  50     //Instantiate a Vector Bspline Interpolator
  51     m_VectorInterpolator =VectorInterpolatorType::New();
  52     m_VectorInterpolator->SetLUTSamplingFactors(m_LUTSamplingFactors);
  53     m_VectorInterpolator->SetSplineOrders(m_SplineOrders);
  54
  55     // Set Bulk transform to NULL (no bulk performed)
  56     m_BulkTransform = NULL;
  57
  58     // Mask
  59     m_Mask=NULL;
  60
  61     // Default grid size is zero
  62     typename RegionType::SizeType size;
  63     typename RegionType::IndexType index;
  64     size.Fill( 0 );
  65     index.Fill( 0 );
  66     m_GridRegion.SetSize( size );
  67     m_GridRegion.SetIndex( index );
  68
  69     m_GridOrigin.Fill( 0.0 );  // default origin is all zeros
  70     m_GridSpacing.Fill( 1.0 ); // default spacing is all ones
  71     m_GridDirection.SetIdentity(); // default spacing is all ones
  72
  73     m_InternalParametersBuffer = ParametersType(0);
  74     // Make sure the parameters pointer is not NULL after construction.
  75     m_InputParametersPointer = &m_InternalParametersBuffer;
  76
  77     // Initialize coeffient images
  78     m_WrappedImage = CoefficientImageType::New();
  79     m_WrappedImage->SetRegions( m_GridRegion );
  80     m_WrappedImage->SetOrigin( m_GridOrigin.GetDataPointer() );
  81     m_WrappedImage->SetSpacing( m_GridSpacing.GetDataPointer() );
  82     m_WrappedImage->SetDirection( m_GridDirection );
  83     m_CoefficientImage = NULL;
  84
  85     // Variables for computing interpolation
  86     for (i=0; i <InputDimension;i++)
  87       {
  88         m_Offset[i] = m_SplineOrders[i] / 2;
  89         m_SplineOrderOdd[i] = m_SplineOrders[i] % 2;
  90       }
  91     m_ValidRegion = m_GridRegion;
  92
  93     // Initialize jacobian images
  94     for (i=0; i <OutputDimension;i++)
  95       {
  96         m_JacobianImage[i] = JacobianImageType::New();
  97         m_JacobianImage[i]->SetRegions( m_GridRegion );
  98         m_JacobianImage[i]->SetOrigin( m_GridOrigin.GetDataPointer() );
  99         m_JacobianImage[i]->SetSpacing( m_GridSpacing.GetDataPointer() );
 100         m_JacobianImage[i]->SetDirection( m_GridDirection );
 101       }
 102
 103     /** Fixed Parameters store the following information:
 104      *     Grid Size
 105      *     Grid Origin
 106      *     Grid Spacing
 107      *     Grid Direction */
 108     //JV we add the splineOrders, LUTsamplingfactor, m_Mask  and m_BulkTransform
 109     /*
 110       Spline orders
 111       Sampling factors
 112       m_Mask
 113       m_BulkTransform
 114       The size of these is equal to the  NInputDimensions
 115     *********************************************************/
 116     this->m_FixedParameters.SetSize ( NInputDimensions * (NInputDimensions + 5)+2 );
 117     this->m_FixedParameters.Fill ( 0.0 );
 118     for ( i=0; i<NInputDimensions; i++)
 119       {
 120         this->m_FixedParameters[2*NInputDimensions+i] = m_GridSpacing[i];
 121       }
 122     for (unsigned int di=0; di<NInputDimensions; di++)
 123       {
 124         for (unsigned int dj=0; dj<NInputDimensions; dj++)
 125           {
 126             this->m_FixedParameters[3*NInputDimensions+(di*NInputDimensions+dj)] = m_GridDirection[di][dj];
 127           }
 128       }
 129
 130     //JV add splineOrders
 131     for ( i=0; i<NInputDimensions; i++)
 132       {
 133         this->m_FixedParameters[ ( (3+NInputDimensions)*NInputDimensions)+i] = (this->GetSplineOrders())[i];
 134       }
 135
 136     //JV add LUTsamplingFactors
 137     for ( i=0; i<NInputDimensions; i++)
 138       {
 139         this->m_FixedParameters[( (4+NInputDimensions)*NInputDimensions)+i ] = this->GetLUTSamplingFactors()[i];
 140       }
 141
 142     // JV add the mask pointer
 143     this->m_FixedParameters[( (5+NInputDimensions)*NInputDimensions)]=(double)((size_t)m_Mask);
 144
 145     // JV add the bulkTransform pointer
 146     this->m_FixedParameters[( (5+NInputDimensions)*NInputDimensions) +1]=(double)((size_t)m_BulkTransform);
 147
 148
 149     // Calculate the PointToIndex matrices
 150     DirectionType scale;
 151     for( unsigned int i=0; i<OutputDimension; i++)
 152       {
 153         scale[i][i] = m_GridSpacing[i];
 154       }
 155
 156     m_IndexToPoint = m_GridDirection * scale;
 157     m_PointToIndex = m_IndexToPoint.GetInverse();
 158
 159     // Jacobian
 160     m_LastJacobianIndex = m_ValidRegion.GetIndex();
 161
 162     //JV initialize some variables for jacobian calculation
 163     m_SupportRegion.SetSize(m_SupportSize);
 164     m_SupportIndex.Fill(0);
 165     m_SupportRegion.SetIndex(m_SupportIndex);
 166     for (  i = 0; i < OutputDimension; i++ )
 167       m_ZeroVector[i]=itk::NumericTraits<JacobianValueType>::Zero;
 168   }
 169
 170
 171   // Destructor
 172   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 173   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 174   ::~BSplineDeformableTransform()
 175   {
 176
 177   }
 178
 179
 180   // JV set Spline Order
 181   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 182   void
 183   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 184   ::SetSplineOrder(const unsigned int & splineOrder)
 185   {
 186     SizeType splineOrders;
 187     for (unsigned int i=0;i<InputDimension;i++)splineOrders[i]=splineOrder;
 188
 189     this->SetSplineOrders(splineOrders);
 190   }
 191
 192
 193   // JV set Spline Orders
 194   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 195   void
 196   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 197   ::SetSplineOrders(const SizeType & splineOrders)
 198   {
 199     if (m_SplineOrders!=splineOrders)
 200       {
 201         m_SplineOrders=splineOrders;
 202
 203         //update the interpolation function
 204         m_VectorInterpolator->SetSplineOrders(m_SplineOrders);
 205
 206         //update the variables for computing interpolation
 207         for (unsigned int i=0; i <InputDimension;i++)
 208         {
 209           m_SupportSize[i] = m_SplineOrders[i]+1;
 210           m_Offset[i] = m_SplineOrders[i] / 2;
 211           m_SplineOrderOdd[i] = m_SplineOrders[i] % 2;
 212         }
 213         this->Modified();
 214       }
 215   }
 216
 217
 218   // JV set sampling factor
 219   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 220   void
 221   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 222   ::SetLUTSamplingFactor( const int & samplingFactor)
 223   {
 224     SizeType samplingFactors;
 225     for (unsigned int i=0; i<NInputDimensions; i++)
 226       samplingFactors[i]=samplingFactor;
 227
 228     //update the interpolation function
 229     this->SetLUTSamplingFactors(samplingFactors);
 230   }
 231
 232
 233   // JV set sampling factors
 234   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 235   void
 236   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 237   ::SetLUTSamplingFactors( const SizeType & samplingFactors)
 238   {
 239     if(m_LUTSamplingFactors!=samplingFactors)
 240       {
 241         for (unsigned int i=0; i<NInputDimensions; i++)
 242           m_LUTSamplingFactors[i]=samplingFactors[i];
 243
 244         //update the interpolation function
 245         m_VectorInterpolator->SetLUTSamplingFactors(m_LUTSamplingFactors);
 246
 247         this->Modified();
 248       }
 249   }
 250
 251
 252   // Get the number of parameters
 253   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 254   unsigned int
 255   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 256   ::GetNumberOfParameters(void) const
 257   {
 258     // The number of parameters equal OutputDimension * number of
 259     // of pixels in the grid region.
 260     return ( static_cast<unsigned int>( OutputDimension ) *
 261              static_cast<unsigned int>( m_GridRegion.GetNumberOfPixels() ) );
 262
 263   }
 264
 265
 266   // Get the number of parameters per dimension
 267   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 268   unsigned int
 269   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 270   ::GetNumberOfParametersPerDimension(void) const
 271   {
 272     // The number of parameters per dimension equal number of
 273     // of pixels in the grid region.
 274     return ( static_cast<unsigned int>( m_GridRegion.GetNumberOfPixels() ) );
 275
 276   }
 277
 278
 279   // Set the grid region
 280   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 281   void
 282   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 283   ::SetGridRegion( const RegionType & region )
 284   {
 285     if ( m_GridRegion != region )
 286       {
 287         m_GridRegion = region;
 288
 289         // set regions for each coefficient and jacobian image
 290         m_WrappedImage->SetRegions( m_GridRegion );
 291         for (unsigned int j=0; j <OutputDimension;j++)
 292           m_JacobianImage[j]->SetRegions( m_GridRegion );
 293
 294         // Set the valid region
 295         // If the grid spans the interval [start, last].
 296         // The valid interval for evaluation is [start+offset, last-offset]
 297         // when spline order is even.
 298         // The valid interval for evaluation is [start+offset, last-offset)
 299         // when spline order is odd.
 300         // Where offset = vcl_floor(spline / 2 ).
 301         // Note that the last pixel is not included in the valid region
 302         // with odd spline orders.
 303         typename RegionType::SizeType size = m_GridRegion.GetSize();
 304         typename RegionType::IndexType index = m_GridRegion.GetIndex();
 305         for ( unsigned int j = 0; j < NInputDimensions; j++ )
 306           {
 307             index[j] +=
 308               static_cast< typename RegionType::IndexValueType >( m_Offset[j] );
 309             size[j] -=
 310               static_cast< typename RegionType::SizeValueType> ( 2 * m_Offset[j] );
 311             m_ValidRegionLast[j] = index[j] +
 312               static_cast< typename RegionType::IndexValueType >( size[j] ) - 1;
 313           }
 314         m_ValidRegion.SetSize( size );
 315         m_ValidRegion.SetIndex( index );
 316
 317         // If we are using the default parameters, update their size and set to identity.
 318         // Input parameters point to internal buffer => using default parameters.
 319         if (m_InputParametersPointer == &m_InternalParametersBuffer)
 320           {
 321             // Check if we need to resize the default parameter buffer.
 322             if ( m_InternalParametersBuffer.GetSize() != this->GetNumberOfParameters() )
 323               {
 324                 m_InternalParametersBuffer.SetSize( this->GetNumberOfParameters() );
 325                 // Fill with zeros for identity.
 326                 m_InternalParametersBuffer.Fill( 0 );
 327               }
 328           }
 329
 330         this->Modified();
 331       }
 332   }
 333
 334
 335   // Set the grid spacing
 336   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 337   void
 338   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 339   ::SetGridSpacing( const SpacingType & spacing )
 340   {
 341     if ( m_GridSpacing != spacing )
 342       {
 343         m_GridSpacing = spacing;
 344
 345         // Set spacing for each coefficient and jacobian image
 346         m_WrappedImage->SetSpacing( m_GridSpacing.GetDataPointer() );
 347         for (unsigned int j=0; j <OutputDimension; j++) m_JacobianImage[j]->SetSpacing( m_GridSpacing.GetDataPointer() );
 348
 349         // Set scale
 350         DirectionType scale;
 351         for( unsigned int i=0; i<OutputDimension; i++)
 352           {
 353             scale[i][i] = m_GridSpacing[i];
 354           }
 355
 356         m_IndexToPoint = m_GridDirection * scale;
 357         m_PointToIndex = m_IndexToPoint.GetInverse();
 358
 359         this->Modified();
 360       }
 361
 362   }
 363
 364   // Set the grid direction
 365   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 366   void
 367   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 368   ::SetGridDirection( const DirectionType & direction )
 369   {
 370     if ( m_GridDirection != direction )
 371       {
 372         m_GridDirection = direction;
 373
 374         // Set direction for each coefficient and jacobian image
 375         m_WrappedImage->SetDirection( m_GridDirection );
 376         for (unsigned int j=0; j <OutputDimension; j++) m_JacobianImage[j]->SetDirection( m_GridDirection );
 377
 378         // Set scale
 379         DirectionType scale;
 380         for( unsigned int i=0; i<OutputDimension; i++)
 381           {
 382             scale[i][i] = m_GridSpacing[i];
 383           }
 384
 385         m_IndexToPoint = m_GridDirection * scale;
 386         m_PointToIndex = m_IndexToPoint.GetInverse();
 387
 388         this->Modified();
 389       }
 390
 391   }
 392
 393
 394   // Set the grid origin
 395   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 396   void
 397   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 398   ::SetGridOrigin( const OriginType& origin )
 399   {
 400     if ( m_GridOrigin != origin )
 401       {
 402         m_GridOrigin = origin;
 403
 404         // Set origin for each coefficient and jacobianimage
 405         m_WrappedImage->SetOrigin( m_GridOrigin.GetDataPointer() );
 406         for (unsigned int j=0; j <OutputDimension; j++) m_JacobianImage[j]->SetOrigin( m_GridOrigin.GetDataPointer() );
 407
 408         this->Modified();
 409       }
 410
 411   }
 412
 413
 414   // Set the parameters
 415   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 416   void
 417   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 418   ::SetIdentity()
 419   {
 420     if( m_InputParametersPointer )
 421       {
 422         ParametersType * parameters =
 423           const_cast<ParametersType *>( m_InputParametersPointer );
 424         parameters->Fill( 0.0 );
 425         this->Modified();
 426       }
 427     else
 428       {
 429         itkExceptionMacro( << "Input parameters for the spline haven't been set ! "
 430                            << "Set them using the SetParameters or SetCoefficientImage method first." );
 431       }
 432   }
 433
 434   // Set the parameters
 435   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 436   void
 437   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 438   ::SetParameters( const ParametersType & parameters )
 439   {
 440
 441     // Check if the number of parameters match the
 442     // Expected number of parameters
 443     if ( parameters.Size() != this->GetNumberOfParameters() )
 444       {
 445         itkExceptionMacro(<<"Mismatched between parameters size "
 446                           << parameters.size()
 447                           << " and region size "
 448                           << m_GridRegion.GetNumberOfPixels() );
 449       }
 450
 451     // Clean up buffered parameters
 452     m_InternalParametersBuffer = ParametersType( 0 );
 453
 454     // Keep a reference to the input parameters
 455     m_InputParametersPointer = &parameters;
 456
 457     // Wrap flat array as images of coefficients
 458     this->WrapAsImages();
 459
 460     //Set input to vector interpolator
 461     m_VectorInterpolator->SetInputImage(this->GetCoefficientImage());
 462
 463     // Modified is always called since we just have a pointer to the
 464     // parameters and cannot know if the parameters have changed.
 465     this->Modified();
 466   }
 467
 468
 469   // Set the Fixed Parameters
 470   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 471   void
 472   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 473   ::SetFixedParameters( const ParametersType & parameters )
 474   {
 475
 476     // JV number should be exact, no defaults for spacing
 477     if ( parameters.Size() != NInputDimensions * (5 + NInputDimensions)+2 )
 478       {
 479         itkExceptionMacro(<< "Mismatched between parameters size "
 480                           << parameters.size()
 481                           << " and number of fixed parameters "
 482                           << NInputDimensions * (5 + NInputDimensions)+2 );
 483       }
 484     /*********************************************************
 485     Fixed Parameters store the following information:
 486         Grid Size
 487         Grid Origin
 488         Grid Spacing
 489         Grid Direction */
 490     // JV we add the splineOrders, LUTsamplingfactor, mask pointer and bulktransform pointer
 491     /*
 492       Spline orders
 493       Sampling factors
 494       m_Mask
 495       m_BulkTransform
 496     The size of these is equal to the  NInputDimensions
 497     *********************************************************/
 498
 499     /** Set the Grid Parameters */
 500     SizeType   gridSize;
 501     for (unsigned int i=0; i<NInputDimensions; i++)
 502       {
 503         gridSize[i] = static_cast<int> (parameters[i]);
 504       }
 505     RegionType bsplineRegion;
 506     bsplineRegion.SetSize( gridSize );
 507
 508     /** Set the Origin Parameters */
 509     OriginType origin;
 510     for (unsigned int i=0; i<NInputDimensions; i++)
 511       {
 512         origin[i] = parameters[NInputDimensions+i];
 513       }
 514
 515     /** Set the Spacing Parameters */
 516     SpacingType spacing;
 517     for (unsigned int i=0; i<NInputDimensions; i++)
 518       {
 519         spacing[i] = parameters[2*NInputDimensions+i];
 520       }
 521
 522     /** Set the Direction Parameters */
 523     DirectionType direction;
 524     for (unsigned int di=0; di<NInputDimensions; di++)
 525       {
 526         for (unsigned int dj=0; dj<NInputDimensions; dj++)
 527           {
 528             direction[di][dj] = parameters[3*NInputDimensions+(di*NInputDimensions+dj)];
 529           }
 530       }
 531
 532     //JV add the spline orders
 533     SizeType splineOrders;
 534     for (unsigned int i=0; i<NInputDimensions; i++)
 535       {
 536         splineOrders[i]=parameters[(3+NInputDimensions)*NInputDimensions+i];
 537       }
 538
 539     //JV add the LUT sampling factor
 540     SizeType  samplingFactors;
 541     for (unsigned int i=0; i<NInputDimensions; i++)
 542       {
 543         samplingFactors[i]=parameters[(4+NInputDimensions)*NInputDimensions+i];
 544       }
 545
 546     //JV add the MaskPointer
 547     m_Mask=((MaskPointer)( (size_t)parameters[(5+NInputDimensions)*NInputDimensions]));
 548
 549     //JV add the MaskPointer
 550     m_BulkTransform=((BulkTransformPointer)( (size_t)parameters[(5+NInputDimensions)*NInputDimensions+1]));
 551
 552     //Set the members
 553     this->SetGridSpacing( spacing );
 554     this->SetGridDirection( direction );
 555     this->SetGridOrigin( origin );
 556     this->SetGridRegion( bsplineRegion );
 557
 558     //JV add the LUT parameters
 559     this->SetSplineOrders( splineOrders );
 560     this->SetLUTSamplingFactors( samplingFactors );
 561
 562   }
 563
 564
 565   // Wrap flat parameters as images
 566   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 567   void
 568   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 569   ::WrapAsImages()
 570   {
 571     //JV Wrap parameter array in vectorial image, changed parameter order: A1x A1y A1z, A2x ....
 572     PixelType * dataPointer =reinterpret_cast<PixelType *>( const_cast<double *>(m_InputParametersPointer->data_block() )) ;
 573     unsigned int numberOfPixels = m_GridRegion.GetNumberOfPixels();
 574
 575     m_WrappedImage->GetPixelContainer()->SetImportPointer( dataPointer,numberOfPixels);//InputDimension
 576     m_CoefficientImage = m_WrappedImage;
 577
 578     //JV Wrap jacobian into OutputDimension X Vectorial images
 579     this->m_Jacobian.set_size( OutputDimension, this->GetNumberOfParameters() );
 580
 581     // Use memset to set the memory
 582     JacobianPixelType * jacobianDataPointer = reinterpret_cast<JacobianPixelType *>(this->m_Jacobian.data_block());
 583     memset(jacobianDataPointer, 0,  OutputDimension*numberOfPixels*sizeof(JacobianPixelType));
 584     m_LastJacobianIndex = m_ValidRegion.GetIndex();
 585     m_NeedResetJacobian = false;
 586
 587     for (unsigned int j=0; j<OutputDimension; j++)
 588       {
 589         m_JacobianImage[j]->GetPixelContainer()->
 590           SetImportPointer( jacobianDataPointer, numberOfPixels );
 591         jacobianDataPointer += numberOfPixels;
 592      }
 593   }
 594
 595
 596   // Set the parameters by value
 597   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 598   void
 599   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 600   ::SetParametersByValue( const ParametersType & parameters )
 601   {
 602
 603     // Check if the number of parameters match the
 604     // Expected number of parameters
 605     if ( parameters.Size() != this->GetNumberOfParameters() )
 606       {
 607         itkExceptionMacro(<<"Mismatched between parameters size "
 608                           << parameters.size()
 609                           << " and region size "
 610                           << m_GridRegion.GetNumberOfPixels() );
 611       }
 612
 613     // Copy it
 614     m_InternalParametersBuffer = parameters;
 615     m_InputParametersPointer = &m_InternalParametersBuffer;
 616
 617     // Wrap flat array as images of coefficients
 618     this->WrapAsImages();
 619
 620     // Modified is always called since we just have a pointer to the
 621     // Parameters and cannot know if the parameters have changed.
 622     this->Modified();
 623
 624   }
 625
 626   // Get the parameters
 627   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 628   const
 629   typename BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 630   ::ParametersType &
 631   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 632   ::GetParameters( void ) const
 633   {
 634     /** NOTE: For efficiency, this class does not keep a copy of the parameters -
 635      * it just keeps pointer to input parameters.
 636      */
 637     if (NULL == m_InputParametersPointer)
 638       {
 639         itkExceptionMacro( <<"Cannot GetParameters() because m_InputParametersPointer is NULL. Perhaps SetCoefficientImage() has been called causing the NULL pointer." );
 640       }
 641
 642     return (*m_InputParametersPointer);
 643   }
 644
 645
 646   // Get the parameters
 647   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 648   const
 649   typename BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 650   ::ParametersType &
 651   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 652   ::GetFixedParameters( void ) const
 653   {
 654     RegionType resRegion = this->GetGridRegion(  );
 655
 656     for (unsigned int i=0; i<NInputDimensions; i++)
 657       {
 658         this->m_FixedParameters[i] = (resRegion.GetSize())[i];
 659       }
 660     for (unsigned int i=0; i<NInputDimensions; i++)
 661       {
 662         this->m_FixedParameters[NInputDimensions+i] = (this->GetGridOrigin())[i];
 663       }
 664     for (unsigned int i=0; i<NInputDimensions; i++)
 665       {
 666         this->m_FixedParameters[2*NInputDimensions+i] =  (this->GetGridSpacing())[i];
 667       }
 668     for (unsigned int di=0; di<NInputDimensions; di++)
 669       {
 670         for (unsigned int dj=0; dj<NInputDimensions; dj++)
 671           {
 672             this->m_FixedParameters[3*NInputDimensions+(di*NInputDimensions+dj)] = (this->GetGridDirection())[di][dj];
 673           }
 674       }
 675
 676     //JV add splineOrders
 677     for (unsigned int i=0; i<NInputDimensions; i++)
 678       {
 679         this->m_FixedParameters[(3+NInputDimensions)*NInputDimensions+i] = (this->GetSplineOrders())[i];
 680       }
 681
 682     //JV add LUTsamplingFactor
 683     for (unsigned int i=0; i<NInputDimensions; i++)
 684       {
 685         this->m_FixedParameters[(4+NInputDimensions)*NInputDimensions+i] = (this->GetLUTSamplingFactors())[i];
 686       }
 687
 688     //JV add the mask
 689     this->m_FixedParameters[(5+NInputDimensions)*NInputDimensions]=(double)((size_t) m_Mask);
 690
 691     //JV add the bulktransform pointer
 692     this->m_FixedParameters[(5+NInputDimensions)*NInputDimensions+1]=(double)((size_t) m_BulkTransform);
 693
 694     return (this->m_FixedParameters);
 695   }
 696
 697
 698   // Set the B-Spline coefficients using input images
 699   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 700   void
 701   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 702   ::SetCoefficientImage( CoefficientImagePointer image )
 703   {
 704
 705     this->SetGridRegion( image->GetBufferedRegion() );
 706     this->SetGridSpacing( image->GetSpacing() );
 707     this->SetGridDirection( image->GetDirection() );
 708     this->SetGridOrigin( image->GetOrigin() );
 709     m_CoefficientImage = image;
 710
 711     // Update the interpolator
 712     m_VectorInterpolator->SetInputImage(this->GetCoefficientImage());
 713
 714     // Clean up buffered parameters
 715     m_InternalParametersBuffer = ParametersType( 0 );
 716     m_InputParametersPointer  = NULL;
 717
 718   }
 719
 720
 721   // Print self
 722   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 723   void
 724   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 725   ::PrintSelf(std::ostream &os, itk::Indent indent) const
 726   {
 727
 728     this->Superclass::PrintSelf(os, indent);
 729
 730     os << indent << "GridRegion: " << m_GridRegion << std::endl;
 731     os << indent << "GridOrigin: " << m_GridOrigin << std::endl;
 732     os << indent << "GridSpacing: " << m_GridSpacing << std::endl;
 733     os << indent << "GridDirection: " << m_GridDirection << std::endl;
 734     os << indent << "IndexToPoint: " << m_IndexToPoint << std::endl;
 735     os << indent << "PointToIndex: " << m_PointToIndex << std::endl;
 736
 737     os << indent << "CoefficientImage: [ ";
 738     os << m_CoefficientImage.GetPointer() << " ]" << std::endl;
 739
 740     os << indent << "WrappedImage: [ ";
 741     os << m_WrappedImage.GetPointer() << " ]" << std::endl;
 742
 743     os << indent << "InputParametersPointer: "
 744        << m_InputParametersPointer << std::endl;
 745     os << indent << "ValidRegion: " << m_ValidRegion << std::endl;
 746     os << indent << "LastJacobianIndex: " << m_LastJacobianIndex << std::endl;
 747     os << indent << "BulkTransform: ";
 748     os << m_BulkTransform << std::endl;
 749
 750     if ( m_BulkTransform )
 751       {
 752         os << indent << "BulkTransformType: "
 753            << m_BulkTransform->GetNameOfClass() << std::endl;
 754       }
 755     os << indent << "VectorBSplineInterpolator: ";
 756     os << m_VectorInterpolator.GetPointer() << std::endl;
 757     os << indent << "Mask: ";
 758     os << m_Mask<< std::endl;
 759   }
 760
 761
 762   // Verify
 763   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 764   bool
 765   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 766   ::InsideValidRegion( const ContinuousIndexType& index ) const
 767   {
 768     bool inside = true;
 769
 770     if ( !m_ValidRegion.IsInside( index ) )
 771       {
 772         inside = false;
 773       }
 774     //JV verify for each input dimension
 775     if ( inside)
 776       {
 777         typedef typename ContinuousIndexType::ValueType ValueType;
 778         for( unsigned int j = 0; j < InputDimension; j++ )
 779           {
 780             if (m_SplineOrderOdd[j])
 781               {
 782                 if ( index[j] >= static_cast<ValueType>( m_ValidRegionLast[j] ) )
 783                   {
 784                     inside = false;
 785                     break;
 786                   }
 787               }
 788           }
 789       }
 790
 791     return inside;
 792   }
 793
 794
 795   // Transform a point
 796   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 797   typename BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 798   ::OutputPointType
 799   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 800   ::TransformPoint(const InputPointType &inputPoint) const
 801   {
 802
 803    InputPointType transformedPoint;
 804    OutputPointType outputPoint;
 805
 806    // BulkTransform
 807    if ( m_BulkTransform )
 808      {
 809        transformedPoint = m_BulkTransform->TransformPoint( inputPoint );
 810      }
 811    else
 812       {
 813         transformedPoint = inputPoint;
 814       }
 815
 816     // Deformable transform
 817     if ( m_CoefficientImage )
 818       {
 819         // Check if inside mask
 820         if(m_Mask &&  !(m_Mask->IsInside(inputPoint) ) )
 821           {
 822             // Outside: no (deformable) displacement
 823             return transformedPoint;
 824           }
 825
 826         // Check if inside valid region
 827         bool inside = true;
 828         ContinuousIndexType index;
 829         this->TransformPointToContinuousIndex( inputPoint, index );
 830         inside = this->InsideValidRegion( index );
 831         if ( !inside )
 832           {
 833             // Outside: no (deformable) displacement
 834             return transformedPoint;
 835           }
 836
 837         // Call the vector interpolator
 838         itk::Vector<TCoordRep,OutputDimension> displacement=m_VectorInterpolator->EvaluateAtContinuousIndex(index);
 839
 840         // Return the results
 841         outputPoint = transformedPoint+displacement;
 842
 843       }
 844
 845     else
 846       {
 847         itkWarningMacro( << "B-spline coefficients have not been set" );
 848         outputPoint = transformedPoint;
 849       }
 850
 851     return outputPoint;
 852   }
 853
 854
 855
 856   //JV Deformably transform a point
 857   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 858   typename BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 859   ::OutputPointType
 860   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 861   ::DeformablyTransformPoint(const InputPointType &inputPoint) const
 862   {
 863     OutputPointType outputPoint;
 864     if ( m_CoefficientImage )
 865       {
 866
 867         // Check if inside mask
 868         if(m_Mask && !(m_Mask->IsInside(inputPoint) ) )
 869           {
 870             // Outside: no (deformable) displacement
 871             return inputPoint;
 872           }
 873
 874         // Check if inside
 875         bool inside = true;
 876         ContinuousIndexType index;
 877         this->TransformPointToContinuousIndex( inputPoint, index );
 878         inside = this->InsideValidRegion( index );
 879
 880         if ( !inside )
 881           {
 882             //outside: no (deformable) displacement
 883              return inputPoint;
 884              //return outputPoint;
 885           }
 886
 887         // Call the vector interpolator
 888         itk::Vector<TCoordRep,OutputDimension> displacement=m_VectorInterpolator->EvaluateAtContinuousIndex(index);
 889
 890         // Return the results
 891         outputPoint = inputPoint+displacement;
 892       }
 893
 894     // No coefficients available
 895     else
 896       {
 897         itkWarningMacro( << "B-spline coefficients have not been set" );
 898         outputPoint = inputPoint;
 899       }
 900
 901     return outputPoint;
 902   }
 903
 904
 905   // JV weights are identical as for transformpoint, could be done simultaneously in metric!!!!
 906   // Compute the Jacobian in one position
 907   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 908   const
 909   typename BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 910   ::JacobianType &
 911   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
 912   ::GetJacobian( const InputPointType & point ) const
 913   {
 914     // Can only compute Jacobian if parameters are set via
 915     // SetParameters or SetParametersByValue
 916     //     if( m_InputParametersPointer == NULL )
 917     //       {
 918     //  itkExceptionMacro( <<"Cannot compute Jacobian: parameters not set" );
 919     //       }
 920
 921     if (m_NeedResetJacobian)
 922       ResetJacobian();
 923
 924     //========================================================
 925     // For each dimension, copy the weight to the support region
 926     //========================================================
 927
 928     // Check if inside mask
 929     if(m_Mask &&  !(m_Mask->IsInside(point) ) )
 930       {
 931         // Outside: no (deformable) displacement
 932         return this->m_Jacobian;
 933       }
 934
 935     //Get index
 936     this->TransformPointToContinuousIndex( point, m_Index );
 937
 938     // NOTE: if the support region does not lie totally within the grid
 939     // we assume zero displacement and return the input point
 940     if ( !this->InsideValidRegion( m_Index ) )
 941       {
 942         return this->m_Jacobian;
 943       }
 944
 945     //Compute interpolation weights
 946     const WeightsDataType *weights=NULL;
 947     m_VectorInterpolator->EvaluateWeightsAtContinuousIndex( m_Index, &weights, m_LastJacobianIndex);
 948     m_SupportRegion.SetIndex( m_LastJacobianIndex );
 949
 950     //Reset the iterators
 951     unsigned int j = 0;
 952     for ( j = 0; j < OutputDimension; j++ )
 953       m_Iterator[j] = IteratorType( m_JacobianImage[j], m_SupportRegion);
 954
 955     // For each dimension, copy the weight to the support region
 956     while ( ! (m_Iterator[0]).IsAtEnd() )
 957       {
 958         //copy weight to jacobian image
 959         for ( j = 0; j < OutputDimension; j++ )
 960           {
 961             m_ZeroVector[j]=*weights;
 962             (m_Iterator[j]).Set( m_ZeroVector);
 963             m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
 964             ++(m_Iterator[j]);
 965           }
 966         // go to next coefficient in the support region
 967         weights++;
 968       }
 969     m_NeedResetJacobian = true;
 970
 971     // Return the result
 972     return this->m_Jacobian;
 973
 974   }
 975
 976
 977   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
 978   void
 979   BSplineDeformableTransform<TCoordRep, NInputDimensions,NOutputDimensions>
 980   ::TransformPointToContinuousIndex( const InputPointType & point, ContinuousIndexType & index ) const
 981   {
 982     unsigned int j;
 983
 984     itk::Vector<double, OutputDimension> tvector;
 985
 986     for ( j = 0; j < OutputDimension; j++ )
 987       {
 988         tvector[j] = point[j] - this->m_GridOrigin[j];
 989       }
 990
 991     itk::Vector<double, OutputDimension> cvector;
 992
 993     cvector = m_PointToIndex * tvector;
 994
 995     for ( j = 0; j < OutputDimension; j++ )
 996       {
 997         index[j] = static_cast< typename ContinuousIndexType::CoordRepType >( cvector[j] );
 998       }
 999   }
1000
1001   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
1002   unsigned
1003   BSplineDeformableTransform<TCoordRep, NInputDimensions,NOutputDimensions>
1004   ::SetJacobianImageData(JacobianPixelType * jacobianDataPointer, unsigned dim)
1005   {
1006     unsigned int numberOfPixels = m_GridRegion.GetNumberOfPixels();
1007     m_JacobianImage[dim]->GetPixelContainer()->SetImportPointer(jacobianDataPointer, numberOfPixels);
1008     return numberOfPixels;
1009   }
1010
1011   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
1012   void
1013   BSplineDeformableTransform<TCoordRep, NInputDimensions,NOutputDimensions>
1014   ::ResetJacobian() const
1015   {
1016     //========================================================
1017     // Zero all components of jacobian
1018     //========================================================
1019     // JV not thread safe (m_LastJacobianIndex), instantiate N transforms
1020     // NOTE: for efficiency, we only need to zero out the coefficients
1021     // that got fill last time this method was called.
1022
1023     unsigned int j = 0;
1024
1025     //Define the  region for each jacobian image
1026     m_SupportRegion.SetIndex(m_LastJacobianIndex);
1027
1028     //Initialize the iterators
1029     for (j = 0; j < OutputDimension; j++)
1030       m_Iterator[j] = IteratorType(m_JacobianImage[j], m_SupportRegion);
1031
1032     //Set the previously-set to zero
1033     while (!(m_Iterator[0]).IsAtEnd())
1034       {
1035         for (j = 0; j < OutputDimension; j++)
1036           {
1037             m_Iterator[j].Set(m_ZeroVector);
1038             ++(m_Iterator[j]);
1039           }
1040       }
1041     m_NeedResetJacobian = false;
1042   }
1043 } // namespace
1044
1045 #endif