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