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