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