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