itk/clitkVectorBSplineInterpolateImageFunction.txx

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