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