itk/clitkVectorBSplineInterpolateImageFunction.txx

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