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