itk/itkRayCastInterpolateImageFunctionWithOrigin.txx

   1 #ifndef __itkRayCastInterpolateImageFunctionWithOrigin_txx
   2 #define __itkRayCastInterpolateImageFunctionWithOrigin_txx
   3 #include "itkRayCastInterpolateImageFunctionWithOrigin.h"
   4
   5 #include "vnl/vnl_math.h"
   6
   7
   8 // Put the helper class in an anonymous namespace so that it is not
   9 // exposed to the user
  10 namespace {
  11
  12 /** \class Helper class to maintain state when casting a ray.
  13  *  This helper class keeps the RayCastInterpolateImageFunctionWithOrigin thread safe.
  14  */
  15 template <class TInputImage, class TCoordRep = float>
  16 class RayCastHelper
  17 {
  18 public:
  19   /** Constants for the image dimensions */
  20   itkStaticConstMacro(InputImageDimension, unsigned int,
  21                       TInputImage::ImageDimension);
  22
  23   /**
  24    * Type of the Transform Base class
  25    * The fixed image should be a 3D image
  26    */
  27   typedef itk::Transform<TCoordRep,3,3> TransformType;
  28
  29   typedef typename TransformType::Pointer            TransformPointer;
  30   typedef typename TransformType::InputPointType     InputPointType;
  31   typedef typename TransformType::OutputPointType    OutputPointType;
  32   typedef typename TransformType::ParametersType     TransformParametersType;
  33   typedef typename TransformType::JacobianType       TransformJacobianType;
  34
  35   typedef typename TInputImage::SizeType             SizeType;
  36   typedef itk::Vector<TCoordRep, 3>                  DirectionType;
  37   typedef itk::Point<TCoordRep, 3>                   PointType;
  38
  39   typedef TInputImage                                InputImageType;
  40   typedef typename InputImageType::PixelType         PixelType;
  41   typedef typename InputImageType::IndexType         IndexType;
  42
  43   /**
  44    * Set the image class
  45    */
  46   void SetImage(const InputImageType *input)
  47     {
  48     m_Image = input;
  49     }
  50
  51   /**
  52    *  Initialise the ray using the position and direction of a line.
  53    *
  54    *  \param RayPosn       The position of the ray in 3D (mm).
  55    *  \param RayDirn       The direction of the ray in 3D (mm).
  56    *
  57    *  \return True if this is a valid ray.
  58    */
  59   bool SetRay(OutputPointType RayPosn, DirectionType RayDirn);
  60
  61
  62   /** \brief
  63    *  Integrate the interpolated intensities along the ray and
  64    *  return the result.
  65    *
  66    *  This routine can be called after instantiating the ray and
  67    *  calling SetProjectionCoord2D() or Reset(). It may then be called
  68    *  as many times thereafter for different 2D projection
  69    *  coordinates.
  70    *
  71    *  \param integral      The integrated intensities along the ray.
  72    *
  73    * \return True if a valid ray was specified.
  74    */
  75   bool Integrate(double &integral)
  76     {
  77     return IntegrateAboveThreshold(integral, 0);
  78     };
  79
  80
  81   /** \brief
  82    * Integrate the interpolated intensities above a given threshold,
  83    * along the ray and return the result.
  84    *
  85    * This routine can be called after instantiating the ray and
  86    * calling SetProjectionCoord2D() or Reset(). It may then be called
  87    * as many times thereafter for different 2D projection
  88    * coordinates.
  89    *
  90    * \param integral      The integrated intensities along the ray.
  91    * \param threshold     The integration threshold [default value: 0]
  92    *
  93    * \return True if a valid ray was specified.
  94    */
  95   bool IntegrateAboveThreshold(double &integral, double threshold);
  96
  97   /** \brief
  98    * Increment each of the intensities of the 4 planar voxels
  99    * surrounding the current ray point.
 100    *
 101    * \parameter increment      Intensity increment for each of the current 4 voxels
 102    */
 103   void IncrementIntensities(double increment=1);
 104
 105   /// Reset the iterator to the start of the ray.
 106   void Reset(void);
 107
 108   /// Return the interpolated intensity of the current ray point.
 109   double GetCurrentIntensity(void) const;
 110
 111   /// Return the ray point spacing in mm
 112   double GetRayPointSpacing(void) const {
 113     typename InputImageType::SpacingType spacing=this->m_Image->GetSpacing();
 114
 115     if (m_ValidRay)
 116       return vcl_sqrt(m_VoxelIncrement[0]*spacing[0]*m_VoxelIncrement[0]*spacing[0]
 117                     + m_VoxelIncrement[1]*spacing[1]*m_VoxelIncrement[1]*spacing[1]
 118                     + m_VoxelIncrement[2]*spacing[2]*m_VoxelIncrement[2]*spacing[2] );
 119     else
 120       return 0.;
 121   };
 122
 123   /// Set the initial zero state of the object
 124   void ZeroState();
 125
 126   /// Initialise the object
 127   void Initialise(void);
 128
 129 protected:
 130   /// Calculate the endpoint coordinats of the ray in voxels.
 131   void EndPointsInVoxels(void);
 132
 133   /**
 134    * Calculate the incremental direction vector in voxels, 'dVoxel',
 135    * required to traverse the ray.
 136    */
 137   void CalcDirnVector(void);
 138
 139   /**
 140    * Reduce the length of the ray until both start and end
 141    * coordinates lie inside the volume.
 142    *
 143    * \return True if a valid ray has been, false otherwise.
 144    */
 145   bool AdjustRayLength(void);
 146
 147   /**
 148    *   Obtain pointers to the four voxels surrounding the point where the ray
 149    *   enters the volume.
 150    */
 151   void InitialiseVoxelPointers(void);
 152
 153   /// Increment the voxel pointers surrounding the current point on the ray.
 154   void IncrementVoxelPointers(void);
 155
 156   /// Record volume dimensions and resolution
 157   void RecordVolumeDimensions(void);
 158
 159   /// Define the corners of the volume
 160   void DefineCorners(void);
 161
 162   /** \brief
 163    * Calculate the planes which define the volume.
 164    *
 165    * Member function to calculate the equations of the planes of 4 of
 166    * the sides of the volume, calculate the positions of the 8 corners
 167    * of the volume in mm in World, also calculate the values of the
 168    * slopes of the lines which go to make up the volume( defined as
 169    * lines in cube x,y,z dirn and then each of these lines has a slope
 170    * in the world x,y,z dirn [3]) and finally also to return the length
 171    * of the sides of the lines in mm.
 172    */
 173   void CalcPlanesAndCorners(void);
 174
 175   /** \brief
 176    *  Calculate the ray intercepts with the volume.
 177    *
 178    *  See where the ray cuts the volume, check that truncation does not occur,
 179    *  if not, then start ray where it first intercepts the volume and set
 180    *  x_max to be where it leaves the volume.
 181    *
 182    *  \return True if a valid ray has been specified, false otherwise.
 183    */
 184   bool CalcRayIntercepts(void);
 185
 186   /**
 187    *   The ray is traversed by stepping in the axial direction
 188    *   that enables the greatest number of planes in the volume to be
 189    *   intercepted.
 190    */
 191   typedef enum {
 192     UNDEFINED_DIRECTION=0,        //!< Undefined
 193     TRANSVERSE_IN_X,              //!< x
 194     TRANSVERSE_IN_Y,              //!< y
 195     TRANSVERSE_IN_Z,              //!< z
 196     LAST_DIRECTION
 197   } TraversalDirection;
 198
 199   // Cache the image in the structure. Skip the smart pointer for
 200   // efficiency. This inner class will go in/out of scope with every
 201   // call to Evaluate()
 202   const InputImageType *m_Image;
 203
 204   /// Flag indicating whether the current ray is valid
 205   bool m_ValidRay;
 206
 207   /** \brief
 208    * The start position of the ray in voxels.
 209    *
 210    * NB. Two of the components of this coordinate (i.e. those lying within
 211    * the planes of voxels being traversed) will be shifted by half a
 212    * voxel. This enables indices of the neighbouring voxels within the plane
 213    * to be determined by simply casting to 'int' and optionally adding 1.
 214    */
 215   double m_RayVoxelStartPosition[3];
 216
 217   /** \brief
 218    * The end coordinate of the ray in voxels.
 219    *
 220    * NB. Two of the components of this coordinate (i.e. those lying within
 221    * the planes of voxels being traversed) will be shifted by half a
 222    * voxel. This enables indices of the neighbouring voxels within the plane
 223    * to be determined by simply casting to 'int' and optionally adding 1.
 224    */
 225   double m_RayVoxelEndPosition[3];
 226
 227
 228   /** \brief
 229    * The current coordinate on the ray in voxels.
 230    *
 231    * NB. Two of the components of this coordinate (i.e. those lying within
 232    * the planes of voxels being traversed) will be shifted by half a
 233    * voxel. This enables indices of the neighbouring voxels within the plane
 234    * to be determined by simply casting to 'int' and optionally adding 1.
 235    */
 236   double m_Position3Dvox[3];
 237
 238   /** The incremental direction vector of the ray in voxels. */
 239   double m_VoxelIncrement[3];
 240
 241   /// The direction in which the ray is incremented thorough the volume (x, y or z).
 242   TraversalDirection m_TraversalDirection;
 243
 244   /// The total number of planes of voxels traversed by the ray.
 245   int m_TotalRayVoxelPlanes;
 246
 247   /// The current number of planes of voxels traversed by the ray.
 248   int m_NumVoxelPlanesTraversed;
 249
 250   /// Pointers to the current four voxels surrounding the ray's trajectory.
 251   const PixelType *m_RayIntersectionVoxels[4];
 252
 253   /**
 254    * The voxel coordinate of the bottom-left voxel of the current
 255    * four voxels surrounding the ray's trajectory.
 256    */
 257   int m_RayIntersectionVoxelIndex[3];
 258
 259   /// The dimension in voxels of the 3D volume in along the x axis
 260   int m_NumberOfVoxelsInX;
 261   /// The dimension in voxels of the 3D volume in along the y axis
 262   int m_NumberOfVoxelsInY;
 263   /// The dimension in voxels of the 3D volume in along the z axis
 264   int m_NumberOfVoxelsInZ;
 265
 266   /// Voxel dimension in x
 267   double m_VoxelDimensionInX;
 268   /// Voxel dimension in y
 269   double m_VoxelDimensionInY;
 270   /// Voxel dimension in z
 271   double m_VoxelDimensionInZ;
 272
 273   /// The coordinate of the point at which the ray enters the volume in mm.
 274   double m_RayStartCoordInMM[3];
 275   /// The coordinate of the point at which the ray exits the volume in mm.
 276   double m_RayEndCoordInMM[3];
 277
 278
 279   /** \brief
 280       Planes which define the boundary of the volume in mm
 281       (six planes and four parameters: Ax+By+Cz+D). */
 282   double m_BoundingPlane[6][4];
 283   /// The eight corners of the volume (x,y,z coordinates for each).
 284   double m_BoundingCorner[8][3];
 285
 286   /// The position of the ray
 287   double m_CurrentRayPositionInMM[3];
 288
 289   /// The direction of the ray
 290   double m_RayDirectionInMM[3];
 291
 292 };
 293
 294 /* -----------------------------------------------------------------------
 295    Initialise() - Initialise the object
 296    ----------------------------------------------------------------------- */
 297
 298 template<class TInputImage, class TCoordRep>
 299 void
 300 RayCastHelper<TInputImage, TCoordRep>
 301 ::Initialise(void)
 302 {
 303   // Save the dimensions of the volume and calculate the bounding box
 304   this->RecordVolumeDimensions();
 305
 306   // Calculate the planes and corners which define the volume.
 307   this->DefineCorners();
 308   this->CalcPlanesAndCorners();
 309 }
 310
 311
 312 /* -----------------------------------------------------------------------
 313    RecordVolumeDimensions() - Record volume dimensions and resolution
 314    ----------------------------------------------------------------------- */
 315
 316 template<class TInputImage, class TCoordRep>
 317 void
 318 RayCastHelper<TInputImage, TCoordRep>
 319 ::RecordVolumeDimensions(void)
 320 {
 321   typename InputImageType::SpacingType spacing=this->m_Image->GetSpacing();
 322   SizeType dim=this->m_Image->GetLargestPossibleRegion().GetSize();
 323
 324   m_NumberOfVoxelsInX = dim[0];
 325   m_NumberOfVoxelsInY = dim[1];
 326   m_NumberOfVoxelsInZ = dim[2];
 327
 328   m_VoxelDimensionInX = spacing[0];
 329   m_VoxelDimensionInY = spacing[1];
 330   m_VoxelDimensionInZ = spacing[2];
 331 }
 332
 333
 334 /* -----------------------------------------------------------------------
 335    DefineCorners() - Define the corners of the volume
 336    ----------------------------------------------------------------------- */
 337
 338 template<class TInputImage, class TCoordRep>
 339 void
 340 RayCastHelper<TInputImage, TCoordRep>
 341 ::DefineCorners(void)
 342 {
 343   // Define corner positions as if at the origin
 344
 345   m_BoundingCorner[0][0] =
 346     m_BoundingCorner[1][0] =
 347     m_BoundingCorner[2][0] =
 348     m_BoundingCorner[3][0] = 0;
 349
 350   m_BoundingCorner[4][0] =
 351     m_BoundingCorner[5][0] =
 352     m_BoundingCorner[6][0] =
 353     m_BoundingCorner[7][0] = m_VoxelDimensionInX*m_NumberOfVoxelsInX;
 354
 355   m_BoundingCorner[1][1] =
 356     m_BoundingCorner[3][1] =
 357     m_BoundingCorner[5][1] =
 358     m_BoundingCorner[7][1] = m_VoxelDimensionInY*m_NumberOfVoxelsInY;
 359
 360   m_BoundingCorner[0][1] =
 361     m_BoundingCorner[2][1] =
 362     m_BoundingCorner[4][1] =
 363     m_BoundingCorner[6][1] = 0;
 364
 365   m_BoundingCorner[0][2] =
 366     m_BoundingCorner[1][2] =
 367     m_BoundingCorner[4][2] =
 368     m_BoundingCorner[5][2] =
 369     m_VoxelDimensionInZ*m_NumberOfVoxelsInZ;
 370
 371   m_BoundingCorner[2][2] =
 372     m_BoundingCorner[3][2] =
 373     m_BoundingCorner[6][2] =
 374     m_BoundingCorner[7][2] = 0;
 375
 376 }
 377
 378 /* -----------------------------------------------------------------------
 379    CalcPlanesAndCorners() - Calculate the planes and corners of the volume.
 380    ----------------------------------------------------------------------- */
 381
 382 template<class TInputImage, class TCoordRep>
 383 void
 384 RayCastHelper<TInputImage, TCoordRep>
 385 ::CalcPlanesAndCorners(void)
 386 {
 387   int j;
 388
 389
 390   // find the equations of the planes
 391
 392   int c1=0, c2=0, c3=0;
 393
 394   for (j=0; j<6; j++)
 395     {                                // loop around for planes
 396     switch (j)
 397       {                // which corners to take
 398       case 0:
 399         c1=1; c2=2; c3=3;
 400         break;
 401       case 1:
 402         c1=4; c2=5; c3=6;
 403         break;
 404       case 2:
 405         c1=5; c2=3; c3=7;
 406         break;
 407       case 3:
 408         c1=2; c2=4; c3=6;
 409         break;
 410       case 4:
 411         c1=1; c2=5; c3=0;
 412         break;
 413       case 5:
 414         c1=3; c2=7; c3=2;
 415         break;
 416       }
 417
 418
 419     double line1x, line1y, line1z;
 420     double line2x, line2y, line2z;
 421
 422     // lines from one corner to another in x,y,z dirns
 423     line1x = m_BoundingCorner[c1][0] - m_BoundingCorner[c2][0];
 424     line2x = m_BoundingCorner[c1][0] - m_BoundingCorner[c3][0];
 425
 426     line1y = m_BoundingCorner[c1][1] - m_BoundingCorner[c2][1];
 427     line2y = m_BoundingCorner[c1][1] - m_BoundingCorner[c3][1];
 428
 429     line1z = m_BoundingCorner[c1][2] - m_BoundingCorner[c2][2];
 430     line2z = m_BoundingCorner[c1][2] - m_BoundingCorner[c3][2];
 431
 432     double A, B, C, D;
 433
 434     // take cross product
 435     A = line1y*line2z - line2y*line1z;
 436     B = line2x*line1z - line1x*line2z;
 437     C = line1x*line2y - line2x*line1y;
 438
 439     // find constant
 440     D = -(   A*m_BoundingCorner[c1][0]
 441              + B*m_BoundingCorner[c1][1]
 442              + C*m_BoundingCorner[c1][2] );
 443
 444     // initialise plane value and normalise
 445     m_BoundingPlane[j][0] = A/vcl_sqrt(A*A + B*B + C*C);
 446     m_BoundingPlane[j][1] = B/vcl_sqrt(A*A + B*B + C*C);
 447     m_BoundingPlane[j][2] = C/vcl_sqrt(A*A + B*B + C*C);
 448     m_BoundingPlane[j][3] = D/vcl_sqrt(A*A + B*B + C*C);
 449
 450     if ( (A*A + B*B + C*C) == 0 )
 451       {
 452       itk::ExceptionObject err(__FILE__, __LINE__);
 453       err.SetLocation( ITK_LOCATION );
 454       err.SetDescription( "Division by zero (planes) "
 455                           "- CalcPlanesAndCorners().");
 456       throw err;
 457       }
 458     }
 459
 460 }
 461
 462
 463 /* -----------------------------------------------------------------------
 464    CalcRayIntercepts() - Calculate the ray intercepts with the volume.
 465    ----------------------------------------------------------------------- */
 466
 467 template<class TInputImage, class TCoordRep>
 468 bool
 469 RayCastHelper<TInputImage, TCoordRep>
 470 ::CalcRayIntercepts()
 471 {
 472   double maxInterDist, interDist;
 473   double cornerVect[4][3];
 474   int cross[4][3], noInterFlag[6];
 475   int nSidesCrossed, crossFlag, c[4];
 476   double ax, ay, az, bx, by, bz;
 477   double cubeInter[6][3];
 478   double denom;
 479
 480   int i,j, k;
 481   int NoSides = 6;  // =6 to allow truncation: =4 to remove truncated rays
 482
 483   // Calculate intercept of ray with planes
 484
 485   double interceptx[6], intercepty[6], interceptz[6];
 486   double d[6];
 487
 488   for( j=0; j<NoSides; j++)
 489     {
 490
 491     denom = (  m_BoundingPlane[j][0]*m_RayDirectionInMM[0]
 492                + m_BoundingPlane[j][1]*m_RayDirectionInMM[1]
 493                + m_BoundingPlane[j][2]*m_RayDirectionInMM[2]);
 494
 495     if( (long)(denom*100) != 0 )
 496       {
 497       d[j] = -(   m_BoundingPlane[j][3]
 498                   + m_BoundingPlane[j][0]*m_CurrentRayPositionInMM[0]
 499                   + m_BoundingPlane[j][1]*m_CurrentRayPositionInMM[1]
 500                   + m_BoundingPlane[j][2]*m_CurrentRayPositionInMM[2] ) / denom;
 501
 502       interceptx[j] = m_CurrentRayPositionInMM[0] + d[j]*m_RayDirectionInMM[0];
 503       intercepty[j] = m_CurrentRayPositionInMM[1] + d[j]*m_RayDirectionInMM[1];
 504       interceptz[j] = m_CurrentRayPositionInMM[2] + d[j]*m_RayDirectionInMM[2];
 505
 506       noInterFlag[j] = 1;  //OK
 507       }
 508     else
 509       {
 510       noInterFlag[j] = 0;  //NOT OK
 511       }
 512     }
 513
 514
 515   nSidesCrossed = 0;
 516   for( j=0; j<NoSides; j++ )
 517     {
 518
 519     // Work out which corners to use
 520
 521     if( j==0 )
 522       {
 523       c[0] = 0; c[1] = 1; c[2] = 3; c[3] = 2;
 524       }
 525     else if( j==1 )
 526       {
 527       c[0] = 4; c[1] = 5; c[2] = 7; c[3] = 6;
 528       }
 529     else if( j==2 )
 530       {
 531       c[0] = 1; c[1] = 5; c[2] = 7; c[3] = 3;
 532       }
 533     else if( j==3 )
 534       {
 535       c[0] = 0; c[1] = 2; c[2] = 6; c[3] = 4;
 536       }
 537     else if( j==4 )
 538       { //TOP
 539       c[0] = 0; c[1] = 1; c[2] = 5; c[3] = 4;
 540       }
 541     else if( j==5 )
 542       { //BOTTOM
 543       c[0] = 2; c[1] = 3; c[2] = 7; c[3] = 6;
 544       }
 545
 546     // Calculate vectors from corner of ct volume to intercept.
 547     for( i=0; i<4; i++ )
 548       {
 549       if( noInterFlag[j]==1 )
 550         {
 551         cornerVect[i][0] = m_BoundingCorner[c[i]][0] - interceptx[j];
 552         cornerVect[i][1] = m_BoundingCorner[c[i]][1] - intercepty[j];
 553         cornerVect[i][2] = m_BoundingCorner[c[i]][2] - interceptz[j];
 554         }
 555       else if( noInterFlag[j]==0 )
 556         {
 557         cornerVect[i][0] = 0;
 558         cornerVect[i][1] = 0;
 559         cornerVect[i][2] = 0;
 560         }
 561
 562       }
 563
 564     // Do cross product with these vectors
 565     for( i=0; i<4; i++ )
 566       {
 567       if( i==3 )
 568         {
 569         k = 0;
 570         }
 571       else
 572         {
 573         k = i+1;
 574         }
 575       ax = cornerVect[i][0];
 576       ay = cornerVect[i][1];
 577       az = cornerVect[i][2];
 578       bx = cornerVect[k][0];
 579       by = cornerVect[k][1];
 580       bz = cornerVect[k][2];
 581
 582       // The int and divide by 100 are to avoid rounding errors.  If
 583       // these are not included then you get values fluctuating around
 584       // zero and so in the subsequent check, all the values are not
 585       // above or below zero.  NB. If you "INT" by too much here though
 586       // you can get problems in the corners of your volume when rays
 587       // are allowed to go through more than one plane.
 588       cross[i][0] = (int)((ay*bz - az*by)/100);
 589       cross[i][1] = (int)((az*bx - ax*bz)/100);
 590       cross[i][2] = (int)((ax*by - ay*bx)/100);
 591       }
 592
 593     // See if a sign change occured between all these cross products
 594     // if not, then the ray went through this plane
 595
 596     crossFlag=0;
 597     for( i=0; i<3; i++ )
 598       {
 599       if( (   cross[0][i]<=0
 600               && cross[1][i]<=0
 601               && cross[2][i]<=0
 602               && cross[3][i]<=0)
 603
 604           || (   cross[0][i]>=0
 605                  && cross[1][i]>=0
 606                  && cross[2][i]>=0
 607                  && cross[3][i]>=0) )
 608         {
 609         crossFlag++;
 610         }
 611       }
 612
 613
 614     if( crossFlag==3 && noInterFlag[j]==1 )
 615       {
 616       cubeInter[nSidesCrossed][0] = interceptx[j];
 617       cubeInter[nSidesCrossed][1] = intercepty[j];
 618       cubeInter[nSidesCrossed][2] = interceptz[j];
 619       nSidesCrossed++;
 620       }
 621
 622     } // End of loop over all four planes
 623
 624   m_RayStartCoordInMM[0] = cubeInter[0][0];
 625   m_RayStartCoordInMM[1] = cubeInter[0][1];
 626   m_RayStartCoordInMM[2] = cubeInter[0][2];
 627
 628   m_RayEndCoordInMM[0] = cubeInter[1][0];
 629   m_RayEndCoordInMM[1] = cubeInter[1][1];
 630   m_RayEndCoordInMM[2] = cubeInter[1][2];
 631
 632   if( nSidesCrossed >= 5 )
 633     {
 634     std::cerr << "WARNING: No. of sides crossed equals: " << nSidesCrossed << std::endl;
 635     }
 636
 637   // If 'nSidesCrossed' is larger than 2, this means that the ray goes through
 638   // a corner of the volume and due to rounding errors, the ray is
 639   // deemed to go through more than two planes.  To obtain the correct
 640   // start and end positions we choose the two intercept values which
 641   // are furthest from each other.
 642
 643
 644   if( nSidesCrossed >= 3 )
 645     {
 646     maxInterDist = 0;
 647     for( j=0; j<nSidesCrossed-1; j++ )
 648       {
 649       for( k=j+1; k<nSidesCrossed; k++ )
 650         {
 651         interDist = 0;
 652         for( i=0; i<3; i++ )
 653           {
 654           interDist += (cubeInter[j][i] - cubeInter[k][i])*
 655             (cubeInter[j][i] - cubeInter[k][i]);
 656           }
 657         if( interDist > maxInterDist )
 658           {
 659           maxInterDist = interDist;
 660
 661           m_RayStartCoordInMM[0] = cubeInter[j][0];
 662           m_RayStartCoordInMM[1] = cubeInter[j][1];
 663           m_RayStartCoordInMM[2] = cubeInter[j][2];
 664
 665           m_RayEndCoordInMM[0] = cubeInter[k][0];
 666           m_RayEndCoordInMM[1] = cubeInter[k][1];
 667           m_RayEndCoordInMM[2] = cubeInter[k][2];
 668           }
 669         }
 670       }
 671     nSidesCrossed = 2;
 672     }
 673
 674   if (nSidesCrossed == 2 )
 675     {
 676     return true;
 677     }
 678   else
 679     {
 680     return false;
 681     }
 682 }
 683
 684
 685 /* -----------------------------------------------------------------------
 686    SetRay() - Set the position and direction of the ray
 687    ----------------------------------------------------------------------- */
 688
 689 template<class TInputImage, class TCoordRep>
 690 bool
 691 RayCastHelper<TInputImage, TCoordRep>
 692 ::SetRay(OutputPointType RayPosn, DirectionType RayDirn)
 693 {
 694
 695   // Store the position and direction of the ray
 696   typename TInputImage::SpacingType spacing=this->m_Image->GetSpacing();
 697   SizeType dim=this->m_Image->GetLargestPossibleRegion().GetSize();
 698
 699   // we need to translate the _center_ of the volume to the origin
 700   m_NumberOfVoxelsInX = dim[0];
 701   m_NumberOfVoxelsInY = dim[1];
 702   m_NumberOfVoxelsInZ = dim[2];
 703
 704   m_VoxelDimensionInX = spacing[0];
 705   m_VoxelDimensionInY = spacing[1];
 706   m_VoxelDimensionInZ = spacing[2];
 707
 708   // (Aviv) Incorporating a fix by Jian Wu
 709   // http://public.kitware.com/pipermail/insight-users/2006-March/017265.html
 710   m_CurrentRayPositionInMM[0] = RayPosn[0];
 711   m_CurrentRayPositionInMM[1] = RayPosn[1];
 712   m_CurrentRayPositionInMM[2] = RayPosn[2];
 713
 714   m_RayDirectionInMM[0] = RayDirn[0];
 715   m_RayDirectionInMM[1] = RayDirn[1];
 716   m_RayDirectionInMM[2] = RayDirn[2];
 717
 718   // Compute the ray path for this coordinate in mm
 719
 720   m_ValidRay = this->CalcRayIntercepts();
 721
 722   if (! m_ValidRay)
 723     {
 724     Reset();
 725     return false;
 726     }
 727
 728   // Convert the start and end coordinates of the ray to voxels
 729
 730   this->EndPointsInVoxels();
 731
 732   /* Calculate the ray direction vector in voxels and hence the voxel
 733      increment required to traverse the ray, and the number of
 734      interpolation points on the ray.
 735
 736      This routine also shifts the coordinate frame by half a voxel for
 737      two of the directional components (i.e. those lying within the
 738      planes of voxels being traversed). */
 739
 740   this->CalcDirnVector();
 741
 742
 743   /* Reduce the length of the ray until both start and end
 744      coordinates lie inside the volume. */
 745
 746   m_ValidRay = this->AdjustRayLength();
 747
 748   // Reset the iterator to the start of the ray.
 749
 750   Reset();
 751
 752   return m_ValidRay;
 753 }
 754
 755
 756 /* -----------------------------------------------------------------------
 757    EndPointsInVoxels() - Convert the endpoints to voxels
 758    ----------------------------------------------------------------------- */
 759
 760 template<class TInputImage, class TCoordRep>
 761 void
 762 RayCastHelper<TInputImage, TCoordRep>
 763 ::EndPointsInVoxels(void)
 764 {
 765   m_RayVoxelStartPosition[0] = m_RayStartCoordInMM[0]/m_VoxelDimensionInX;
 766   m_RayVoxelStartPosition[1] = m_RayStartCoordInMM[1]/m_VoxelDimensionInY;
 767   m_RayVoxelStartPosition[2] = m_RayStartCoordInMM[2]/m_VoxelDimensionInZ;
 768
 769   m_RayVoxelEndPosition[0] = m_RayEndCoordInMM[0]/m_VoxelDimensionInX;
 770   m_RayVoxelEndPosition[1] = m_RayEndCoordInMM[1]/m_VoxelDimensionInY;
 771   m_RayVoxelEndPosition[2] = m_RayEndCoordInMM[2]/m_VoxelDimensionInZ;
 772
 773 }
 774
 775
 776 /* -----------------------------------------------------------------------
 777    CalcDirnVector() - Calculate the incremental direction vector in voxels.
 778    ----------------------------------------------------------------------- */
 779
 780 template<class TInputImage, class TCoordRep>
 781 void
 782 RayCastHelper<TInputImage, TCoordRep>
 783 ::CalcDirnVector(void)
 784 {
 785   double xNum, yNum, zNum;
 786
 787   // Calculate the number of voxels in each direction
 788
 789   xNum = vcl_fabs(m_RayVoxelStartPosition[0] - m_RayVoxelEndPosition[0]);
 790   yNum = vcl_fabs(m_RayVoxelStartPosition[1] - m_RayVoxelEndPosition[1]);
 791   zNum = vcl_fabs(m_RayVoxelStartPosition[2] - m_RayVoxelEndPosition[2]);
 792
 793   // The direction iterated in is that with the greatest number of voxels
 794   // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 795
 796   // Iterate in X direction
 797
 798   if( (xNum >= yNum) && (xNum >= zNum) )
 799     {
 800     if( m_RayVoxelStartPosition[0] < m_RayVoxelEndPosition[0] )
 801       {
 802       m_VoxelIncrement[0] = 1;
 803
 804       m_VoxelIncrement[1]
 805         = (m_RayVoxelStartPosition[1]
 806            - m_RayVoxelEndPosition[1])/(m_RayVoxelStartPosition[0]
 807                                         - m_RayVoxelEndPosition[0]);
 808
 809       m_VoxelIncrement[2]
 810         = (m_RayVoxelStartPosition[2]
 811            - m_RayVoxelEndPosition[2])/(m_RayVoxelStartPosition[0]
 812                                         - m_RayVoxelEndPosition[0]);
 813       }
 814     else
 815       {
 816       m_VoxelIncrement[0] = -1;
 817
 818       m_VoxelIncrement[1]
 819         = -(m_RayVoxelStartPosition[1]
 820             - m_RayVoxelEndPosition[1])/(m_RayVoxelStartPosition[0]
 821                                          - m_RayVoxelEndPosition[0]);
 822
 823       m_VoxelIncrement[2]
 824         = -(m_RayVoxelStartPosition[2]
 825             - m_RayVoxelEndPosition[2])/(m_RayVoxelStartPosition[0]
 826                                          - m_RayVoxelEndPosition[0]);
 827       }
 828
 829     // This section is to alter the start position in order to
 830     // place the center of the voxels in there correct positions,
 831     // rather than placing them at the corner of voxels which is
 832     // what happens if this is not carried out.  The reason why
 833     // x has no -0.5 is because this is the direction we are going
 834     // to iterate in and therefore we wish to go from center to
 835     // center rather than finding the surrounding voxels.
 836
 837     m_RayVoxelStartPosition[1] += ( (int)m_RayVoxelStartPosition[0]
 838         - m_RayVoxelStartPosition[0])*m_VoxelIncrement[1]*m_VoxelIncrement[0]
 839       + 0.5*m_VoxelIncrement[1] - 0.5;
 840
 841     m_RayVoxelStartPosition[2] += ( (int)m_RayVoxelStartPosition[0]
 842            - m_RayVoxelStartPosition[0])*m_VoxelIncrement[2]*m_VoxelIncrement[0]
 843       + 0.5*m_VoxelIncrement[2] - 0.5;
 844
 845     m_RayVoxelStartPosition[0] = (int)m_RayVoxelStartPosition[0] + 0.5*m_VoxelIncrement[0];
 846
 847     m_TotalRayVoxelPlanes = (int)xNum;
 848
 849     m_TraversalDirection = TRANSVERSE_IN_X;
 850     }
 851
 852   // Iterate in Y direction
 853
 854   else if( (yNum >= xNum) && (yNum >= zNum) )
 855     {
 856
 857     if( m_RayVoxelStartPosition[1] < m_RayVoxelEndPosition[1] )
 858       {
 859       m_VoxelIncrement[1] = 1;
 860
 861       m_VoxelIncrement[0]
 862         = (m_RayVoxelStartPosition[0]
 863            - m_RayVoxelEndPosition[0])/(m_RayVoxelStartPosition[1]
 864                                         - m_RayVoxelEndPosition[1]);
 865
 866       m_VoxelIncrement[2]
 867         = (m_RayVoxelStartPosition[2]
 868            - m_RayVoxelEndPosition[2])/(m_RayVoxelStartPosition[1]
 869                                         - m_RayVoxelEndPosition[1]);
 870       }
 871     else
 872       {
 873       m_VoxelIncrement[1] = -1;
 874
 875       m_VoxelIncrement[0]
 876         = -(m_RayVoxelStartPosition[0]
 877             - m_RayVoxelEndPosition[0])/(m_RayVoxelStartPosition[1]
 878                                          - m_RayVoxelEndPosition[1]);
 879
 880       m_VoxelIncrement[2]
 881         = -(m_RayVoxelStartPosition[2]
 882             - m_RayVoxelEndPosition[2])/(m_RayVoxelStartPosition[1]
 883                                          - m_RayVoxelEndPosition[1]);
 884       }
 885
 886     m_RayVoxelStartPosition[0] += ( (int)m_RayVoxelStartPosition[1]
 887                                     - m_RayVoxelStartPosition[1])*m_VoxelIncrement[0]*m_VoxelIncrement[1]
 888       + 0.5*m_VoxelIncrement[0] - 0.5;
 889
 890     m_RayVoxelStartPosition[2] += ( (int)m_RayVoxelStartPosition[1]
 891                                     - m_RayVoxelStartPosition[1])*m_VoxelIncrement[2]*m_VoxelIncrement[1]
 892       + 0.5*m_VoxelIncrement[2] - 0.5;
 893
 894     m_RayVoxelStartPosition[1] = (int)m_RayVoxelStartPosition[1] + 0.5*m_VoxelIncrement[1];
 895
 896     m_TotalRayVoxelPlanes = (int)yNum;
 897
 898     m_TraversalDirection = TRANSVERSE_IN_Y;
 899     }
 900
 901   // Iterate in Z direction
 902
 903   else
 904     {
 905
 906     if( m_RayVoxelStartPosition[2] < m_RayVoxelEndPosition[2] )
 907       {
 908       m_VoxelIncrement[2] = 1;
 909
 910       m_VoxelIncrement[0]
 911         = (m_RayVoxelStartPosition[0]
 912            - m_RayVoxelEndPosition[0])/(m_RayVoxelStartPosition[2]
 913                                         - m_RayVoxelEndPosition[2]);
 914
 915       m_VoxelIncrement[1]
 916         = (m_RayVoxelStartPosition[1]
 917            - m_RayVoxelEndPosition[1])/(m_RayVoxelStartPosition[2]
 918                                         - m_RayVoxelEndPosition[2]);
 919       }
 920     else
 921       {
 922       m_VoxelIncrement[2] = -1;
 923
 924       m_VoxelIncrement[0]
 925         = -(m_RayVoxelStartPosition[0]
 926             - m_RayVoxelEndPosition[0])/(m_RayVoxelStartPosition[2]
 927                                          - m_RayVoxelEndPosition[2]);
 928
 929       m_VoxelIncrement[1]
 930         = -(m_RayVoxelStartPosition[1]
 931             - m_RayVoxelEndPosition[1])/(m_RayVoxelStartPosition[2]
 932                                          - m_RayVoxelEndPosition[2]);
 933       }
 934
 935     m_RayVoxelStartPosition[0] += ( (int)m_RayVoxelStartPosition[2]
 936                                     - m_RayVoxelStartPosition[2])*m_VoxelIncrement[0]*m_VoxelIncrement[2]
 937       + 0.5*m_VoxelIncrement[0] - 0.5;
 938
 939     m_RayVoxelStartPosition[1] += ( (int)m_RayVoxelStartPosition[2]
 940                                     - m_RayVoxelStartPosition[2])*m_VoxelIncrement[1]*m_VoxelIncrement[2]
 941       + 0.5*m_VoxelIncrement[1] - 0.5;
 942
 943     m_RayVoxelStartPosition[2] = (int)m_RayVoxelStartPosition[2] + 0.5*m_VoxelIncrement[2];
 944
 945     m_TotalRayVoxelPlanes = (int)zNum;
 946
 947     m_TraversalDirection = TRANSVERSE_IN_Z;
 948     }
 949 }
 950
 951
 952 /* -----------------------------------------------------------------------
 953    AdjustRayLength() - Ensure that the ray lies within the volume
 954    ----------------------------------------------------------------------- */
 955
 956 template<class TInputImage, class TCoordRep>
 957 bool
 958 RayCastHelper<TInputImage, TCoordRep>
 959 ::AdjustRayLength(void)
 960 {
 961   bool startOK, endOK;
 962
 963   int Istart[3];
 964   int Idirn[3];
 965
 966   if (m_TraversalDirection == TRANSVERSE_IN_X)
 967     {
 968     Idirn[0] = 0;
 969     Idirn[1] = 1;
 970     Idirn[2] = 1;
 971     }
 972   else if (m_TraversalDirection == TRANSVERSE_IN_Y)
 973     {
 974     Idirn[0] = 1;
 975     Idirn[1] = 0;
 976     Idirn[2] = 1;
 977     }
 978   else if (m_TraversalDirection == TRANSVERSE_IN_Z)
 979     {
 980     Idirn[0] = 1;
 981     Idirn[1] = 1;
 982     Idirn[2] = 0;
 983     }
 984   else
 985     {
 986     itk::ExceptionObject err(__FILE__, __LINE__);
 987     err.SetLocation( ITK_LOCATION );
 988     err.SetDescription( "The ray traversal direction is unset "
 989                         "- AdjustRayLength().");
 990     throw err;
 991     return false;
 992     }
 993
 994
 995   do
 996     {
 997
 998     startOK = false;
 999     endOK = false;
1000
1001     Istart[0] = (int) vcl_floor(m_RayVoxelStartPosition[0]);
1002     Istart[1] = (int) vcl_floor(m_RayVoxelStartPosition[1]);
1003     Istart[2] = (int) vcl_floor(m_RayVoxelStartPosition[2]);
1004
1005     if( (Istart[0] >= 0) && (Istart[0] + Idirn[0] < m_NumberOfVoxelsInX) &&
1006         (Istart[1] >= 0) && (Istart[1] + Idirn[1] < m_NumberOfVoxelsInY) &&
1007         (Istart[2] >= 0) && (Istart[2] + Idirn[2] < m_NumberOfVoxelsInZ) )
1008       {
1009
1010       startOK = true;
1011       }
1012     else
1013       {
1014       m_RayVoxelStartPosition[0] += m_VoxelIncrement[0];
1015       m_RayVoxelStartPosition[1] += m_VoxelIncrement[1];
1016       m_RayVoxelStartPosition[2] += m_VoxelIncrement[2];
1017
1018       m_TotalRayVoxelPlanes--;
1019       }
1020
1021     Istart[0] = (int) vcl_floor(m_RayVoxelStartPosition[0]
1022                             + m_TotalRayVoxelPlanes*m_VoxelIncrement[0]);
1023
1024     Istart[1] = (int) vcl_floor(m_RayVoxelStartPosition[1]
1025                             + m_TotalRayVoxelPlanes*m_VoxelIncrement[1]);
1026
1027     Istart[2] = (int) vcl_floor(m_RayVoxelStartPosition[2]
1028                             + m_TotalRayVoxelPlanes*m_VoxelIncrement[2]);
1029
1030     if( (Istart[0] >= 0) && (Istart[0] + Idirn[0] < m_NumberOfVoxelsInX) &&
1031         (Istart[1] >= 0) && (Istart[1] + Idirn[1] < m_NumberOfVoxelsInY) &&
1032         (Istart[2] >= 0) && (Istart[2] + Idirn[2] < m_NumberOfVoxelsInZ) )
1033       {
1034
1035       endOK = true;
1036       }
1037     else
1038       {
1039       m_TotalRayVoxelPlanes--;
1040       }
1041
1042     } while ( (! (startOK && endOK)) && (m_TotalRayVoxelPlanes > 1) );
1043
1044
1045   return (startOK && endOK);
1046 }
1047
1048
1049 /* -----------------------------------------------------------------------
1050    Reset() - Reset the iterator to the start of the ray.
1051    ----------------------------------------------------------------------- */
1052
1053 template<class TInputImage, class TCoordRep>
1054 void
1055 RayCastHelper<TInputImage, TCoordRep>
1056 ::Reset(void)
1057 {
1058   int i;
1059
1060   m_NumVoxelPlanesTraversed = -1;
1061
1062   // If this is a valid ray...
1063
1064   if (m_ValidRay)
1065     {
1066     for (i=0; i<3; i++)
1067       {
1068       m_Position3Dvox[i] = m_RayVoxelStartPosition[i];
1069       }
1070     this->InitialiseVoxelPointers();
1071     }
1072
1073   // otherwise set parameters to zero
1074
1075   else
1076     {
1077     for (i=0; i<3; i++)
1078       {
1079       m_RayVoxelStartPosition[i] = 0.;
1080       }
1081     for (i=0; i<3; i++)
1082       {
1083       m_RayVoxelEndPosition[i] = 0.;
1084       }
1085     for (i=0; i<3; i++)
1086       {
1087       m_VoxelIncrement[i] = 0.;
1088       }
1089     m_TraversalDirection = UNDEFINED_DIRECTION;
1090
1091     m_TotalRayVoxelPlanes = 0;
1092
1093     for (i=0; i<4; i++)
1094       {
1095       m_RayIntersectionVoxels[i] = 0;
1096       }
1097     for (i=0; i<3; i++)
1098       {
1099       m_RayIntersectionVoxelIndex[i] = 0;
1100       }
1101     }
1102 }
1103
1104
1105 /* -----------------------------------------------------------------------
1106    InitialiseVoxelPointers() - Obtain pointers to the first four voxels
1107    ----------------------------------------------------------------------- */
1108
1109 template<class TInputImage, class TCoordRep>
1110 void
1111 RayCastHelper<TInputImage, TCoordRep>
1112 ::InitialiseVoxelPointers(void)
1113 {
1114   IndexType index;
1115
1116   int Ix, Iy, Iz;
1117
1118   Ix = (int)(m_RayVoxelStartPosition[0]);
1119   Iy = (int)(m_RayVoxelStartPosition[1]);
1120   Iz = (int)(m_RayVoxelStartPosition[2]);
1121
1122   m_RayIntersectionVoxelIndex[0] = Ix;
1123   m_RayIntersectionVoxelIndex[1] = Iy;
1124   m_RayIntersectionVoxelIndex[2] = Iz;
1125
1126   switch( m_TraversalDirection )
1127     {
1128     case TRANSVERSE_IN_X:
1129       {
1130
1131       if( (Ix >= 0) && (Ix     < m_NumberOfVoxelsInX) &&
1132           (Iy >= 0) && (Iy + 1 < m_NumberOfVoxelsInY) &&
1133           (Iz >= 0) && (Iz + 1 < m_NumberOfVoxelsInZ))
1134         {
1135
1136         index[0]=Ix; index[1]=Iy; index[2]=Iz;
1137         m_RayIntersectionVoxels[0]
1138           = this->m_Image->GetBufferPointer() + this->m_Image->ComputeOffset(index);
1139
1140         index[0]=Ix; index[1]=Iy+1; index[2]=Iz;
1141         m_RayIntersectionVoxels[1]
1142           = ( this->m_Image->GetBufferPointer() + this->m_Image->ComputeOffset(index) );
1143
1144         index[0]=Ix; index[1]=Iy; index[2]=Iz+1;
1145         m_RayIntersectionVoxels[2]
1146           = ( this->m_Image->GetBufferPointer() + this->m_Image->ComputeOffset(index) );
1147
1148         index[0]=Ix; index[1]=Iy+1; index[2]=Iz+1;
1149         m_RayIntersectionVoxels[3]
1150           = ( this->m_Image->GetBufferPointer() + this->m_Image->ComputeOffset(index) );
1151         }
1152       else
1153         {
1154         m_RayIntersectionVoxels[0] =
1155           m_RayIntersectionVoxels[1] =
1156           m_RayIntersectionVoxels[2] =
1157           m_RayIntersectionVoxels[3] = NULL;
1158         }
1159       break;
1160       }
1161
1162     case TRANSVERSE_IN_Y:
1163       {
1164
1165       if( (Ix >= 0) && (Ix + 1 < m_NumberOfVoxelsInX) &&
1166           (Iy >= 0) && (Iy     < m_NumberOfVoxelsInY) &&
1167           (Iz >= 0) && (Iz + 1 < m_NumberOfVoxelsInZ))
1168         {
1169
1170         index[0]=Ix; index[1]=Iy; index[2]=Iz;
1171         m_RayIntersectionVoxels[0] = ( this->m_Image->GetBufferPointer()
1172                                        + this->m_Image->ComputeOffset(index) );
1173
1174         index[0]=Ix+1; index[1]=Iy; index[2]=Iz;
1175         m_RayIntersectionVoxels[1] = ( this->m_Image->GetBufferPointer()
1176                                        + this->m_Image->ComputeOffset(index) );
1177
1178         index[0]=Ix; index[1]=Iy; index[2]=Iz+1;
1179         m_RayIntersectionVoxels[2] = ( this->m_Image->GetBufferPointer()
1180                                        + this->m_Image->ComputeOffset(index) );
1181
1182         index[0]=Ix+1; index[1]=Iy; index[2]=Iz+1;
1183         m_RayIntersectionVoxels[3] = ( this->m_Image->GetBufferPointer()
1184                                        + this->m_Image->ComputeOffset(index) );
1185         }
1186       else
1187         {
1188         m_RayIntersectionVoxels[0]
1189         = m_RayIntersectionVoxels[1]
1190         = m_RayIntersectionVoxels[2]
1191         = m_RayIntersectionVoxels[3] = NULL;
1192         }
1193       break;
1194       }
1195
1196     case TRANSVERSE_IN_Z:
1197       {
1198
1199       if( (Ix >= 0) && (Ix + 1 < m_NumberOfVoxelsInX)   &&
1200           (Iy >= 0) && (Iy + 1 < m_NumberOfVoxelsInY) &&
1201           (Iz >= 0) && (Iz     < m_NumberOfVoxelsInZ))
1202         {
1203
1204         index[0]=Ix; index[1]=Iy; index[2]=Iz;
1205         m_RayIntersectionVoxels[0] = ( this->m_Image->GetBufferPointer()
1206                                        + this->m_Image->ComputeOffset(index) );
1207
1208         index[0]=Ix+1; index[1]=Iy; index[2]=Iz;
1209         m_RayIntersectionVoxels[1] = ( this->m_Image->GetBufferPointer()
1210                                        + this->m_Image->ComputeOffset(index) );
1211
1212         index[0]=Ix; index[1]=Iy+1; index[2]=Iz;
1213         m_RayIntersectionVoxels[2] = ( this->m_Image->GetBufferPointer()
1214                                        + this->m_Image->ComputeOffset(index) );
1215
1216         index[0]=Ix+1; index[1]=Iy+1; index[2]=Iz;
1217         m_RayIntersectionVoxels[3] = ( this->m_Image->GetBufferPointer()
1218                                        + this->m_Image->ComputeOffset(index) );
1219
1220         }
1221       else
1222         {
1223         m_RayIntersectionVoxels[0]
1224         = m_RayIntersectionVoxels[1]
1225         = m_RayIntersectionVoxels[2]
1226         = m_RayIntersectionVoxels[3] = NULL;
1227         }
1228       break;
1229       }
1230
1231     default:
1232       {
1233       itk::ExceptionObject err(__FILE__, __LINE__);
1234       err.SetLocation( ITK_LOCATION );
1235       err.SetDescription( "The ray traversal direction is unset "
1236                           "- InitialiseVoxelPointers().");
1237       throw err;
1238       return;
1239       }
1240     }
1241 }
1242
1243 /* -----------------------------------------------------------------------
1244    IncrementVoxelPointers() - Increment the voxel pointers
1245    ----------------------------------------------------------------------- */
1246
1247 template<class TInputImage, class TCoordRep>
1248 void
1249 RayCastHelper<TInputImage, TCoordRep>
1250 ::IncrementVoxelPointers(void)
1251 {
1252   double xBefore = m_Position3Dvox[0];
1253   double yBefore = m_Position3Dvox[1];
1254   double zBefore = m_Position3Dvox[2];
1255
1256   m_Position3Dvox[0] += m_VoxelIncrement[0];
1257   m_Position3Dvox[1] += m_VoxelIncrement[1];
1258   m_Position3Dvox[2] += m_VoxelIncrement[2];
1259
1260   int dx = ((int) m_Position3Dvox[0]) - ((int) xBefore);
1261   int dy = ((int) m_Position3Dvox[1]) - ((int) yBefore);
1262   int dz = ((int) m_Position3Dvox[2]) - ((int) zBefore);
1263
1264   m_RayIntersectionVoxelIndex[0] += dx;
1265   m_RayIntersectionVoxelIndex[1] += dy;
1266   m_RayIntersectionVoxelIndex[2] += dz;
1267
1268   int totalRayVoxelPlanes
1269     = dx + dy*m_NumberOfVoxelsInX + dz*m_NumberOfVoxelsInX*m_NumberOfVoxelsInY;
1270
1271   m_RayIntersectionVoxels[0] += totalRayVoxelPlanes;
1272   m_RayIntersectionVoxels[1] += totalRayVoxelPlanes;
1273   m_RayIntersectionVoxels[2] += totalRayVoxelPlanes;
1274   m_RayIntersectionVoxels[3] += totalRayVoxelPlanes;
1275 }
1276
1277
1278 /* -----------------------------------------------------------------------
1279    GetCurrentIntensity() - Get the intensity of the current ray point.
1280    ----------------------------------------------------------------------- */
1281
1282 template<class TInputImage, class TCoordRep>
1283 double
1284 RayCastHelper<TInputImage, TCoordRep>
1285 ::GetCurrentIntensity(void) const
1286 {
1287   double a, b, c, d;
1288   double y, z;
1289
1290   if (! m_ValidRay)
1291     {
1292     return 0;
1293     }
1294   a = (double) (*m_RayIntersectionVoxels[0]);
1295   b = (double) (*m_RayIntersectionVoxels[1] - a);
1296   c = (double) (*m_RayIntersectionVoxels[2] - a);
1297   d = (double) (*m_RayIntersectionVoxels[3] - a - b - c);
1298
1299   switch( m_TraversalDirection )
1300     {
1301     case TRANSVERSE_IN_X:
1302       {
1303       y = m_Position3Dvox[1] - vcl_floor(m_Position3Dvox[1]);
1304       z = m_Position3Dvox[2] - vcl_floor(m_Position3Dvox[2]);
1305       break;
1306       }
1307     case TRANSVERSE_IN_Y:
1308       {
1309       y = m_Position3Dvox[0] - vcl_floor(m_Position3Dvox[0]);
1310       z = m_Position3Dvox[2] - vcl_floor(m_Position3Dvox[2]);
1311       break;
1312       }
1313     case TRANSVERSE_IN_Z:
1314       {
1315       y = m_Position3Dvox[0] - vcl_floor(m_Position3Dvox[0]);
1316       z = m_Position3Dvox[1] - vcl_floor(m_Position3Dvox[1]);
1317       break;
1318       }
1319     default:
1320       {
1321       itk::ExceptionObject err(__FILE__, __LINE__);
1322       err.SetLocation( ITK_LOCATION );
1323       err.SetDescription( "The ray traversal direction is unset "
1324                           "- GetCurrentIntensity().");
1325       throw err;
1326       return 0;
1327       }
1328     }
1329
1330   return a + b*y + c*z + d*y*z;
1331 }
1332
1333 /* -----------------------------------------------------------------------
1334    IncrementIntensities() - Increment the intensities of the current ray point
1335    ----------------------------------------------------------------------- */
1336
1337 template<class TInputImage, class TCoordRep>
1338 void
1339 RayCastHelper<TInputImage, TCoordRep>
1340 ::IncrementIntensities(double increment)
1341 {
1342   short inc = (short) vcl_floor(increment + 0.5);
1343
1344   if (! m_ValidRay)
1345     {
1346     return;
1347     }
1348   *m_RayIntersectionVoxels[0] += inc;
1349   *m_RayIntersectionVoxels[1] += inc;
1350   *m_RayIntersectionVoxels[2] += inc;
1351   *m_RayIntersectionVoxels[3] += inc;
1352
1353   return;
1354 }
1355
1356
1357 /* -----------------------------------------------------------------------
1358    IntegrateAboveThreshold() - Integrate intensities above a threshold.
1359    ----------------------------------------------------------------------- */
1360
1361 template<class TInputImage, class TCoordRep>
1362 bool
1363 RayCastHelper<TInputImage, TCoordRep>
1364 ::IntegrateAboveThreshold(double &integral, double threshold)
1365 {
1366   double intensity;
1367 //  double posn3D_x, posn3D_y, posn3D_z;
1368
1369   integral = 0.;
1370
1371   // Check if this is a valid ray
1372
1373   if (! m_ValidRay)
1374     {
1375     return false;
1376     }
1377   /* Step along the ray as quickly as possible
1378      integrating the interpolated intensities. */
1379
1380   for (m_NumVoxelPlanesTraversed=0;
1381        m_NumVoxelPlanesTraversed<m_TotalRayVoxelPlanes;
1382        m_NumVoxelPlanesTraversed++)
1383     {
1384     intensity = this->GetCurrentIntensity();
1385
1386     if (intensity > threshold)
1387       {
1388       integral += intensity - threshold;
1389       }
1390     this->IncrementVoxelPointers();
1391     }
1392
1393   /* The ray passes through the volume one plane of voxels at a time,
1394      however, if its moving diagonally the ray points will be further
1395      apart so account for this by scaling by the distance moved. */
1396
1397   integral *= this->GetRayPointSpacing();
1398
1399   return true;
1400 }
1401
1402 /* -----------------------------------------------------------------------
1403    ZeroState() - Set the default (zero) state of the object
1404    ----------------------------------------------------------------------- */
1405
1406 template<class TInputImage, class TCoordRep>
1407 void
1408 RayCastHelper<TInputImage, TCoordRep>
1409 ::ZeroState()
1410 {
1411   int i;
1412
1413   m_ValidRay = false;
1414
1415   m_NumberOfVoxelsInX = 0;
1416   m_NumberOfVoxelsInY = 0;
1417   m_NumberOfVoxelsInZ = 0;
1418
1419   m_VoxelDimensionInX = 0;
1420   m_VoxelDimensionInY = 0;
1421   m_VoxelDimensionInZ = 0;
1422
1423   for (i=0; i<3; i++)
1424     {
1425     m_CurrentRayPositionInMM[i] = 0.;
1426     }
1427   for (i=0; i<3; i++)
1428     {
1429     m_RayDirectionInMM[i] = 0.;
1430     }
1431   for (i=0; i<3; i++)
1432     {
1433     m_RayVoxelStartPosition[i] = 0.;
1434     }
1435   for (i=0; i<3; i++)
1436     {
1437     m_RayVoxelEndPosition[i] = 0.;
1438     }
1439   for (i=0; i<3; i++)
1440     {
1441     m_VoxelIncrement[i] = 0.;
1442     }
1443   m_TraversalDirection = UNDEFINED_DIRECTION;
1444
1445   m_TotalRayVoxelPlanes = 0;
1446   m_NumVoxelPlanesTraversed = -1;
1447
1448   for (i=0; i<4; i++)
1449     {
1450     m_RayIntersectionVoxels[i] = 0;
1451     }
1452   for (i=0; i<3; i++)
1453     {
1454     m_RayIntersectionVoxelIndex[i] = 0;
1455     }
1456 }
1457 }; // end of anonymous namespace
1458
1459
1460 namespace itk
1461 {
1462
1463 /**************************************************************************
1464  *
1465  *
1466  * Rest of this code is the actual RayCastInterpolateImageFunctionWithOrigin
1467  * class
1468  *
1469  *
1470  **************************************************************************/
1471
1472 /* -----------------------------------------------------------------------
1473    Constructor
1474    ----------------------------------------------------------------------- */
1475
1476 template<class TInputImage, class TCoordRep>
1477 RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
1478 ::RayCastInterpolateImageFunctionWithOrigin()
1479 {
1480   m_Threshold = 0.;
1481
1482   m_FocalPoint[0] = 0.;
1483   m_FocalPoint[1] = 0.;
1484   m_FocalPoint[2] = 0.;
1485 }
1486
1487
1488 /* -----------------------------------------------------------------------
1489    PrintSelf
1490    ----------------------------------------------------------------------- */
1491
1492 template<class TInputImage, class TCoordRep>
1493 void
1494 RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
1495 ::PrintSelf(std::ostream& os, Indent indent) const
1496 {
1497   this->Superclass::PrintSelf(os,indent);
1498
1499   os << indent << "Threshold: " << m_Threshold << std::endl;
1500   os << indent << "FocalPoint: " << m_FocalPoint << std::endl;
1501   os << indent << "Transform: " << m_Transform.GetPointer() << std::endl;
1502   os << indent << "Interpolator: " << m_Interpolator.GetPointer() << std::endl;
1503
1504 }
1505
1506 /* -----------------------------------------------------------------------
1507    Evaluate at image index position
1508    ----------------------------------------------------------------------- */
1509
1510 template<class TInputImage, class TCoordRep>
1511 typename RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
1512 ::OutputType
1513 RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
1514 ::Evaluate( const PointType& point ) const
1515 {
1516   double integral = 0;
1517
1518   OutputPointType transformedFocalPoint
1519     = m_Transform->TransformPoint( m_FocalPoint );
1520
1521   DirectionType direction = transformedFocalPoint - point;
1522
1523   RayCastHelper<TInputImage, TCoordRep> ray;
1524   ray.SetImage( this->m_Image );
1525   ray.ZeroState();
1526   ray.Initialise();
1527
1528   // (Aviv) Added support for images with non-zero origin
1529   // ray.SetRay(point, direction);
1530   ray.SetRay(point - this->m_Image->GetOrigin().GetVectorFromOrigin(), direction);
1531   ray.IntegrateAboveThreshold(integral, m_Threshold);
1532
1533   return ( static_cast<OutputType>( integral ));
1534 }
1535
1536 template<class TInputImage, class TCoordRep>
1537 typename RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
1538 ::OutputType
1539 RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
1540 ::EvaluateAtContinuousIndex( const ContinuousIndexType& index ) const
1541 {
1542   OutputPointType point;
1543   this->m_Image->TransformContinuousIndexToPhysicalPoint(index, point);
1544
1545   return this->Evaluate( point );
1546 }
1547
1548 } // namespace itk
1549
1550
1551 #endif