itk/clitkVectorBSplineInterpolateImageFunction.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 _clitkVectorBSplineInterpolateImageFunction_txx
  19 #define _clitkVectorBSplineInterpolateImageFunction_txx
  20 #include "itkConfigure.h"
  21
  22 // Second, redirect to the optimized version if necessary
  23 // #ifdef ITK_USE_OPTIMIZED_REGISTRATION_METHODS
  24 // #include "itkOptVectorBSplineInterpolateImageFunction.txx"
  25 // #else
  26
  27 #include "clitkVectorBSplineInterpolateImageFunction.h"
  28
  29
  30 #include "itkImageLinearIteratorWithIndex.h"
  31 #include "itkImageRegionConstIteratorWithIndex.h"
  32 #include "itkImageRegionIterator.h"
  33
  34 #include "itkVector.h"
  35
  36 #include "itkMatrix.h"
  37
  38 namespace clitk
  39 {
  40
  41 /**
  42  * Constructor
  43  */
  44 template <class TImageType, class TCoordRep, class TCoefficientType>
  45 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
  46 ::VectorBSplineInterpolateImageFunction()
  47 {
  48   m_SplineOrder = 0;
  49   unsigned int SplineOrder = 3;
  50   m_CoefficientFilter = CoefficientFilter::New();
  51   // ***TODO: Should we store coefficients in a variable or retrieve from filter?
  52   m_Coefficients = CoefficientImageType::New();
  53   this->SetSplineOrder(SplineOrder);
  54   this->m_UseImageDirection = false;
  55 }
  56
  57 /**
  58  * Standard "PrintSelf" method
  59  */
  60 template <class TImageType, class TCoordRep, class TCoefficientType>
  61 void
  62 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
  63 ::PrintSelf(
  64   std::ostream& os,
  65   itk::Indent indent) const
  66 {
  67   Superclass::PrintSelf( os, indent );
  68   os << indent << "Spline Order: " << m_SplineOrder << std::endl;
  69   os << indent << "UseImageDirection = "
  70      << (this->m_UseImageDirection ? "On" : "Off") << std::endl;
  71
  72 }
  73
  74
  75 template <class TImageType, class TCoordRep, class TCoefficientType>
  76 void
  77 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
  78 ::SetInputImage(const TImageType * inputData)
  79 {
  80   if ( inputData ) {
  81
  82     DD("calling decomposition filter");
  83     m_CoefficientFilter->SetInput(inputData);
  84
  85     // the Coefficient Filter requires that the spline order and the input data be set.
  86     // TODO:  We need to ensure that this is only run once and only after both input and
  87     //        spline order have been set. Should we force an update after the
  88     //        splineOrder has been set also?
  89
  90     m_CoefficientFilter->Update();
  91     m_Coefficients = m_CoefficientFilter->GetOutput();
  92
  93
  94     // Call the Superclass implementation after, in case the filter
  95     // pulls in  more of the input image
  96     Superclass::SetInputImage(inputData);
  97
  98     m_DataLength = inputData->GetBufferedRegion().GetSize();
  99
 100   } else {
 101     m_Coefficients = NULL;
 102   }
 103 }
 104
 105
 106 template <class TImageType, class TCoordRep, class TCoefficientType>
 107 void
 108 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
 109 ::SetSplineOrder(unsigned int SplineOrder)
 110 {
 111   if (SplineOrder == m_SplineOrder) {
 112     return;
 113   }
 114   m_SplineOrder = SplineOrder;
 115   m_CoefficientFilter->SetSplineOrder( SplineOrder );
 116
 117   //this->SetPoles();/
 118   m_MaxNumberInterpolationPoints = 1;
 119   for (unsigned int n=0; n < ImageDimension; n++) {
 120     m_MaxNumberInterpolationPoints *= ( m_SplineOrder + 1);
 121   }
 122   this->GeneratePointsToIndex( );
 123 }
 124
 125
 126 template <class TImageType, class TCoordRep, class TCoefficientType>
 127 typename
 128 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
 129 ::OutputType
 130 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
 131 ::EvaluateAtContinuousIndex( const ContinuousIndexType & x ) const
 132 {
 133   vnl_matrix<long>        EvaluateIndex(ImageDimension, ( m_SplineOrder + 1 ));
 134
 135   // compute the interpolation indexes
 136   this->DetermineRegionOfSupport(EvaluateIndex, x, m_SplineOrder);
 137
 138   // Determine weights
 139   vnl_matrix<double>        weights(ImageDimension, ( m_SplineOrder + 1 ));
 140   SetInterpolationWeights( x, EvaluateIndex, weights, m_SplineOrder );
 141
 142   // Modify EvaluateIndex at the boundaries using mirror boundary conditions
 143   this->ApplyMirrorBoundaryConditions(EvaluateIndex, m_SplineOrder);
 144
 145   // perform interpolation
 146   //JV
 147   itk::Vector<double, VectorDimension> interpolated;
 148   for (unsigned int i=0; i< VectorDimension; i++) interpolated[i]=itk::NumericTraits<double>::Zero;
 149
 150   IndexType coefficientIndex;
 151   // Step through eachpoint in the N-dimensional interpolation cube.
 152   for (unsigned int p = 0; p < m_MaxNumberInterpolationPoints; p++) {
 153     // translate each step into the N-dimensional index.
 154     //      IndexType pointIndex = PointToIndex( p );
 155
 156     double w = 1.0;
 157     for (unsigned int n = 0; n < ImageDimension; n++ ) {
 158
 159       w *= weights[n][ m_PointsToIndex[p][n] ];
 160       coefficientIndex[n] = EvaluateIndex[n][m_PointsToIndex[p][n]];  // Build up ND index for coefficients.
 161     }
 162     // Convert our step p to the appropriate point in ND space in the
 163     // m_Coefficients cube.
 164     //JV shouldn't be necessary
 165     for (unsigned int i=0; i<VectorDimension; i++)
 166       interpolated[i] += w * m_Coefficients->GetPixel(coefficientIndex)[i];
 167   }
 168
 169   /*  double interpolated = 0.0;
 170     IndexType coefficientIndex;
 171     // Step through eachpoint in the N-dimensional interpolation cube.
 172     for (unsigned int sp = 0; sp <= m_SplineOrder; sp++)
 173       {
 174       for (unsigned int sp1=0; sp1 <= m_SplineOrder; sp1++)
 175         {
 176
 177         double w = 1.0;
 178         for (unsigned int n1 = 0; n1 < ImageDimension; n1++ )
 179           {
 180           w *= weights[n1][ sp1 ];
 181           coefficientIndex[n1] = EvaluateIndex[n1][sp];  // Build up ND index for coefficients.
 182           }
 183
 184           interpolated += w * m_Coefficients->GetPixel(coefficientIndex);
 185         }
 186       }
 187   */
 188   return(interpolated);
 189
 190 }
 191
 192
 193 template <class TImageType, class TCoordRep, class TCoefficientType>
 194 typename
 195 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
 196 :: CovariantVectorType
 197 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
 198 ::EvaluateDerivativeAtContinuousIndex( const ContinuousIndexType & x ) const
 199 {
 200   vnl_matrix<long>        EvaluateIndex(ImageDimension, ( m_SplineOrder + 1 ));
 201
 202   // compute the interpolation indexes
 203   // TODO: Do we need to revisit region of support for the derivatives?
 204   this->DetermineRegionOfSupport(EvaluateIndex, x, m_SplineOrder);
 205
 206   // Determine weights
 207   vnl_matrix<double>        weights(ImageDimension, ( m_SplineOrder + 1 ));
 208   SetInterpolationWeights( x, EvaluateIndex, weights, m_SplineOrder );
 209
 210   vnl_matrix<double>        weightsDerivative(ImageDimension, ( m_SplineOrder + 1));
 211   SetDerivativeWeights( x, EvaluateIndex, weightsDerivative, ( m_SplineOrder ) );
 212
 213   // Modify EvaluateIndex at the boundaries using mirror boundary conditions
 214   this->ApplyMirrorBoundaryConditions(EvaluateIndex, m_SplineOrder);
 215
 216   const InputImageType * inputImage = this->GetInputImage();
 217   const typename InputImageType::SpacingType & spacing = inputImage->GetSpacing();
 218
 219   // Calculate derivative
 220   CovariantVectorType derivativeValue;
 221   double tempValue;
 222   IndexType coefficientIndex;
 223   for (unsigned int n = 0; n < ImageDimension; n++) {
 224     derivativeValue[n] = 0.0;
 225     for (unsigned int p = 0; p < m_MaxNumberInterpolationPoints; p++) {
 226       tempValue = 1.0 ;
 227       for (unsigned int n1 = 0; n1 < ImageDimension; n1++) {
 228         //coefficientIndex[n1] = EvaluateIndex[n1][sp];
 229         coefficientIndex[n1] = EvaluateIndex[n1][m_PointsToIndex[p][n1]];
 230
 231         if (n1 == n) {
 232           //w *= weights[n][ m_PointsToIndex[p][n] ];
 233           tempValue *= weightsDerivative[n1][ m_PointsToIndex[p][n1] ];
 234         } else {
 235           tempValue *= weights[n1][ m_PointsToIndex[p][n1] ];
 236         }
 237       }
 238       derivativeValue[n] += m_Coefficients->GetPixel(coefficientIndex) * tempValue ;
 239     }
 240     derivativeValue[n] /= spacing[n];   // take spacing into account
 241   }
 242
 243 #ifdef ITK_USE_ORIENTED_IMAGE_DIRECTION
 244   if( this->m_UseImageDirection ) {
 245     CovariantVectorType orientedDerivative;
 246     inputImage->TransformLocalVectorToPhysicalVector( derivativeValue, orientedDerivative );
 247     return orientedDerivative;
 248   }
 249 #endif
 250
 251   return(derivativeValue);
 252 }
 253
 254
 255 template <class TImageType, class TCoordRep, class TCoefficientType>
 256 void
 257 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
 258 ::SetInterpolationWeights( const ContinuousIndexType & x, const vnl_matrix<long> & EvaluateIndex,
 259                            vnl_matrix<double> & weights, unsigned int splineOrder ) const
 260 {
 261   // For speed improvements we could make each case a separate function and use
 262   // function pointers to reference the correct weight order.
 263   // Left as is for now for readability.
 264   double w, w2, w4, t, t0, t1;
 265
 266   switch (splineOrder) {
 267   case 3:
 268     for (unsigned int n = 0; n < ImageDimension; n++) {
 269       w = x[n] - (double) EvaluateIndex[n][1];
 270       weights[n][3] = (1.0 / 6.0) * w * w * w;
 271       weights[n][0] = (1.0 / 6.0) + 0.5 * w * (w - 1.0) - weights[n][3];
 272       weights[n][2] = w + weights[n][0] - 2.0 * weights[n][3];
 273       weights[n][1] = 1.0 - weights[n][0] - weights[n][2] - weights[n][3];
 274     }
 275     break;
 276   case 0:
 277     for (unsigned int n = 0; n < ImageDimension; n++) {
 278       weights[n][0] = 1; // implements nearest neighbor
 279     }
 280     break;
 281   case 1:
 282     for (unsigned int n = 0; n < ImageDimension; n++) {
 283       w = x[n] - (double) EvaluateIndex[n][0];
 284       weights[n][1] = w;
 285       weights[n][0] = 1.0 - w;
 286     }
 287     break;
 288   case 2:
 289     for (unsigned int n = 0; n < ImageDimension; n++) {
 290       /* x */
 291       w = x[n] - (double)EvaluateIndex[n][1];
 292       weights[n][1] = 0.75 - w * w;
 293       weights[n][2] = 0.5 * (w - weights[n][1] + 1.0);
 294       weights[n][0] = 1.0 - weights[n][1] - weights[n][2];
 295     }
 296     break;
 297   case 4:
 298     for (unsigned int n = 0; n < ImageDimension; n++) {
 299       /* x */
 300       w = x[n] - (double)EvaluateIndex[n][2];
 301       w2 = w * w;
 302       t = (1.0 / 6.0) * w2;
 303       weights[n][0] = 0.5 - w;
 304       weights[n][0] *= weights[n][0];
 305       weights[n][0] *= (1.0 / 24.0) * weights[n][0];
 306       t0 = w * (t - 11.0 / 24.0);
 307       t1 = 19.0 / 96.0 + w2 * (0.25 - t);
 308       weights[n][1] = t1 + t0;
 309       weights[n][3] = t1 - t0;
 310       weights[n][4] = weights[n][0] + t0 + 0.5 * w;
 311       weights[n][2] = 1.0 - weights[n][0] - weights[n][1] - weights[n][3] - weights[n][4];
 312     }
 313     break;
 314   case 5:
 315     for (unsigned int n = 0; n < ImageDimension; n++) {
 316       /* x */
 317       w = x[n] - (double)EvaluateIndex[n][2];
 318       w2 = w * w;
 319       weights[n][5] = (1.0 / 120.0) * w * w2 * w2;
 320       w2 -= w;
 321       w4 = w2 * w2;
 322       w -= 0.5;
 323       t = w2 * (w2 - 3.0);
 324       weights[n][0] = (1.0 / 24.0) * (1.0 / 5.0 + w2 + w4) - weights[n][5];
 325       t0 = (1.0 / 24.0) * (w2 * (w2 - 5.0) + 46.0 / 5.0);
 326       t1 = (-1.0 / 12.0) * w * (t + 4.0);
 327       weights[n][2] = t0 + t1;
 328       weights[n][3] = t0 - t1;
 329       t0 = (1.0 / 16.0) * (9.0 / 5.0 - t);
 330       t1 = (1.0 / 24.0) * w * (w4 - w2 - 5.0);
 331       weights[n][1] = t0 + t1;
 332       weights[n][4] = t0 - t1;
 333     }
 334     break;
 335   default:
 336     // SplineOrder not implemented yet.
 337     itk::ExceptionObject err(__FILE__, __LINE__);
 338     err.SetLocation( ITK_LOCATION );
 339     err.SetDescription( "SplineOrder must be between 0 and 5. Requested spline order has not been implemented yet." );
 340     throw err;
 341     break;
 342   }
 343
 344 }
 345
 346 template <class TImageType, class TCoordRep, class TCoefficientType>
 347 void
 348 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
 349 ::SetDerivativeWeights( const ContinuousIndexType & x, const vnl_matrix<long> & EvaluateIndex,
 350                         vnl_matrix<double> & weights, unsigned int splineOrder ) const
 351 {
 352   // For speed improvements we could make each case a separate function and use
 353   // function pointers to reference the correct weight order.
 354   // Another possiblity would be to loop inside the case statement (reducing the number
 355   // of switch statement executions to one per routine call.
 356   // Left as is for now for readability.
 357   double w, w1, w2, w3, w4, w5, t, t0, t1, t2;
 358   int derivativeSplineOrder = (int) splineOrder -1;
 359
 360   switch (derivativeSplineOrder) {
 361
 362     // Calculates B(splineOrder) ( (x + 1/2) - xi) - B(splineOrder -1) ( (x - 1/2) - xi)
 363   case -1:
 364     // Why would we want to do this?
 365     for (unsigned int n = 0; n < ImageDimension; n++) {
 366       weights[n][0] = 0.0;
 367     }
 368     break;
 369   case 0:
 370     for (unsigned int n = 0; n < ImageDimension; n++) {
 371       weights[n][0] = -1.0;
 372       weights[n][1] =  1.0;
 373     }
 374     break;
 375   case 1:
 376     for (unsigned int n = 0; n < ImageDimension; n++) {
 377       w = x[n] + 0.5 - (double)EvaluateIndex[n][1];
 378       // w2 = w;
 379       w1 = 1.0 - w;
 380
 381       weights[n][0] = 0.0 - w1;
 382       weights[n][1] = w1 - w;
 383       weights[n][2] = w;
 384     }
 385     break;
 386   case 2:
 387
 388     for (unsigned int n = 0; n < ImageDimension; n++) {
 389       w = x[n] + .5 - (double)EvaluateIndex[n][2];
 390       w2 = 0.75 - w * w;
 391       w3 = 0.5 * (w - w2 + 1.0);
 392       w1 = 1.0 - w2 - w3;
 393
 394       weights[n][0] = 0.0 - w1;
 395       weights[n][1] = w1 - w2;
 396       weights[n][2] = w2 - w3;
 397       weights[n][3] = w3;
 398     }
 399     break;
 400   case 3:
 401
 402     for (unsigned int n = 0; n < ImageDimension; n++) {
 403       w = x[n] + 0.5 - (double)EvaluateIndex[n][2];
 404       w4 = (1.0 / 6.0) * w * w * w;
 405       w1 = (1.0 / 6.0) + 0.5 * w * (w - 1.0) - w4;
 406       w3 = w + w1 - 2.0 * w4;
 407       w2 = 1.0 - w1 - w3 - w4;
 408
 409       weights[n][0] = 0.0 - w1;
 410       weights[n][1] = w1 - w2;
 411       weights[n][2] = w2 - w3;
 412       weights[n][3] = w3 - w4;
 413       weights[n][4] = w4;
 414     }
 415     break;
 416   case 4:
 417     for (unsigned int n = 0; n < ImageDimension; n++) {
 418       w = x[n] + .5 - (double)EvaluateIndex[n][3];
 419       t2 = w * w;
 420       t = (1.0 / 6.0) * t2;
 421       w1 = 0.5 - w;
 422       w1 *= w1;
 423       w1 *= (1.0 / 24.0) * w1;
 424       t0 = w * (t - 11.0 / 24.0);
 425       t1 = 19.0 / 96.0 + t2 * (0.25 - t);
 426       w2 = t1 + t0;
 427       w4 = t1 - t0;
 428       w5 = w1 + t0 + 0.5 * w;
 429       w3 = 1.0 - w1 - w2 - w4 - w5;
 430
 431       weights[n][0] = 0.0 - w1;
 432       weights[n][1] = w1 - w2;
 433       weights[n][2] = w2 - w3;
 434       weights[n][3] = w3 - w4;
 435       weights[n][4] = w4 - w5;
 436       weights[n][5] = w5;
 437     }
 438     break;
 439
 440   default:
 441     // SplineOrder not implemented yet.
 442     itk::ExceptionObject err(__FILE__, __LINE__);
 443     err.SetLocation( ITK_LOCATION );
 444     err.SetDescription( "SplineOrder (for derivatives) must be between 1 and 5. Requested spline order has not been implemented yet." );
 445     throw err;
 446     break;
 447   }
 448
 449 }
 450
 451
 452 // Generates m_PointsToIndex;
 453 template <class TImageType, class TCoordRep, class TCoefficientType>
 454 void
 455 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
 456 ::GeneratePointsToIndex( )
 457 {
 458   // m_PointsToIndex is used to convert a sequential location to an N-dimension
 459   // index vector.  This is precomputed to save time during the interpolation routine.
 460   m_PointsToIndex.resize(m_MaxNumberInterpolationPoints);
 461   for (unsigned int p = 0; p < m_MaxNumberInterpolationPoints; p++) {
 462     int pp = p;
 463     unsigned long indexFactor[ImageDimension];
 464     indexFactor[0] = 1;
 465     for (int j=1; j< static_cast<int>(ImageDimension); j++) {
 466       indexFactor[j] = indexFactor[j-1] * ( m_SplineOrder + 1 );
 467     }
 468     for (int j = (static_cast<int>(ImageDimension) - 1); j >= 0; j--) {
 469       m_PointsToIndex[p][j] = pp / indexFactor[j];
 470       pp = pp % indexFactor[j];
 471     }
 472   }
 473 }
 474
 475 template <class TImageType, class TCoordRep, class TCoefficientType>
 476 void
 477 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
 478 ::DetermineRegionOfSupport( vnl_matrix<long> & evaluateIndex,
 479                             const ContinuousIndexType & x,
 480                             unsigned int splineOrder ) const
 481 {
 482   long indx;
 483
 484 // compute the interpolation indexes
 485   for (unsigned int n = 0; n< ImageDimension; n++) {
 486     if (splineOrder & 1) {   // Use this index calculation for odd splineOrder
 487       indx = (long)vcl_floor((float)x[n]) - splineOrder / 2;
 488       for (unsigned int k = 0; k <= splineOrder; k++) {
 489         evaluateIndex[n][k] = indx++;
 490       }
 491     } else {                   // Use this index calculation for even splineOrder
 492       indx = (long)vcl_floor((float)(x[n] + 0.5)) - splineOrder / 2;
 493       for (unsigned int k = 0; k <= splineOrder; k++) {
 494         evaluateIndex[n][k] = indx++;
 495       }
 496     }
 497   }
 498 }
 499
 500 template <class TImageType, class TCoordRep, class TCoefficientType>
 501 void
 502 VectorBSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
 503 ::ApplyMirrorBoundaryConditions(vnl_matrix<long> & evaluateIndex,
 504                                 unsigned int splineOrder) const
 505 {
 506   for (unsigned int n = 0; n < ImageDimension; n++) {
 507     long dataLength2 = 2 * m_DataLength[n] - 2;
 508
 509     // apply the mirror boundary conditions
 510     // TODO:  We could implement other boundary options beside mirror
 511     if (m_DataLength[n] == 1) {
 512       for (unsigned int k = 0; k <= splineOrder; k++) {
 513         evaluateIndex[n][k] = 0;
 514       }
 515     } else {
 516       for (unsigned int k = 0; k <= splineOrder; k++) {
 517         // btw - Think about this couldn't this be replaced with a more elagent modulus method?
 518         evaluateIndex[n][k] = (evaluateIndex[n][k] < 0L) ? (-evaluateIndex[n][k] - dataLength2 * ((-evaluateIndex[n][k]) / dataLength2))
 519                               : (evaluateIndex[n][k] - dataLength2 * (evaluateIndex[n][k] / dataLength2));
 520         if ((long) m_DataLength[n] <= evaluateIndex[n][k]) {
 521           evaluateIndex[n][k] = dataLength2 - evaluateIndex[n][k];
 522         }
 523       }
 524     }
 525   }
 526 }
 527
 528 } // namespace itk
 529
 530 #endif
 531
 532 //#endif