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 #include "clitkCommon.h"
  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   for(clitk::uint i=0; i<4; i++)
 502     for(clitk::uint j=0; j<3; j++)
 503       cornerVect[i][j] = 0.0; // to avoid warning
 504   int cross[4][3], noInterFlag[6];
 505   int nSidesCrossed, crossFlag, c[4];
 506   double ax, ay, az, bx, by, bz;
 507   double cubeInter[6][3];
 508   double denom;
 509
 510   int i,j, k;
 511   int NoSides = 6;  // =6 to allow truncation: =4 to remove truncated rays
 512
 513   // Calculate intercept of ray with planes
 514
 515   double interceptx[6], intercepty[6], interceptz[6];
 516   double d[6];
 517
 518   for( j=0; j<NoSides; j++) {
 519
 520     denom = (  m_BoundingPlane[j][0]*m_RayDirectionInMM[0]
 521                + m_BoundingPlane[j][1]*m_RayDirectionInMM[1]
 522                + m_BoundingPlane[j][2]*m_RayDirectionInMM[2]);
 523
 524     if( (long)(denom*100) != 0 ) {
 525       d[j] = -(   m_BoundingPlane[j][3]
 526                   + m_BoundingPlane[j][0]*m_CurrentRayPositionInMM[0]
 527                   + m_BoundingPlane[j][1]*m_CurrentRayPositionInMM[1]
 528                   + m_BoundingPlane[j][2]*m_CurrentRayPositionInMM[2] ) / denom;
 529
 530       interceptx[j] = m_CurrentRayPositionInMM[0] + d[j]*m_RayDirectionInMM[0];
 531       intercepty[j] = m_CurrentRayPositionInMM[1] + d[j]*m_RayDirectionInMM[1];
 532       interceptz[j] = m_CurrentRayPositionInMM[2] + d[j]*m_RayDirectionInMM[2];
 533
 534       noInterFlag[j] = 1;  //OK
 535     } else {
 536       noInterFlag[j] = 0;  //NOT OK
 537     }
 538   }
 539
 540
 541   nSidesCrossed = 0;
 542   for( j=0; j<NoSides; j++ ) {
 543
 544     // Work out which corners to use
 545
 546     if( j==0 ) {
 547       c[0] = 0;
 548       c[1] = 1;
 549       c[2] = 3;
 550       c[3] = 2;
 551     } else if( j==1 ) {
 552       c[0] = 4;
 553       c[1] = 5;
 554       c[2] = 7;
 555       c[3] = 6;
 556     } else if( j==2 ) {
 557       c[0] = 1;
 558       c[1] = 5;
 559       c[2] = 7;
 560       c[3] = 3;
 561     } else if( j==3 ) {
 562       c[0] = 0;
 563       c[1] = 2;
 564       c[2] = 6;
 565       c[3] = 4;
 566     } else if( j==4 ) {
 567       //TOP
 568       c[0] = 0;
 569       c[1] = 1;
 570       c[2] = 5;
 571       c[3] = 4;
 572     } else if( j==5 ) {
 573       //BOTTOM
 574       c[0] = 2;
 575       c[1] = 3;
 576       c[2] = 7;
 577       c[3] = 6;
 578     }
 579
 580     // Calculate vectors from corner of ct volume to intercept.
 581     for( i=0; i<4; i++ ) {
 582       if( noInterFlag[j]==1 ) {
 583         cornerVect[i][0] = m_BoundingCorner[c[i]][0] - interceptx[j];
 584         cornerVect[i][1] = m_BoundingCorner[c[i]][1] - intercepty[j];
 585         cornerVect[i][2] = m_BoundingCorner[c[i]][2] - interceptz[j];
 586       } else if( noInterFlag[j]==0 ) {
 587         cornerVect[i][0] = 0;
 588         cornerVect[i][1] = 0;
 589         cornerVect[i][2] = 0;
 590       }
 591
 592     }
 593
 594     // Do cross product with these vectors
 595     for( i=0; i<4; i++ ) {
 596       if( i==3 ) {
 597         k = 0;
 598       } else {
 599         k = i+1;
 600       }
 601       ax = cornerVect[i][0];
 602       ay = cornerVect[i][1];
 603       az = cornerVect[i][2];
 604       bx = cornerVect[k][0];
 605       by = cornerVect[k][1];
 606       bz = cornerVect[k][2];
 607
 608       // The int and divide by 100 are to avoid rounding errors.  If
 609       // these are not included then you get values fluctuating around
 610       // zero and so in the subsequent check, all the values are not
 611       // above or below zero.  NB. If you "INT" by too much here though
 612       // you can get problems in the corners of your volume when rays
 613       // are allowed to go through more than one plane.
 614       cross[i][0] = (int)((ay*bz - az*by)/100);
 615       cross[i][1] = (int)((az*bx - ax*bz)/100);
 616       cross[i][2] = (int)((ax*by - ay*bx)/100);
 617     }
 618
 619     // See if a sign change occured between all these cross products
 620     // if not, then the ray went through this plane
 621
 622     crossFlag=0;
 623     for( i=0; i<3; i++ ) {
 624       if( (   cross[0][i]<=0
 625               && cross[1][i]<=0
 626               && cross[2][i]<=0
 627               && cross[3][i]<=0)
 628
 629           || (   cross[0][i]>=0
 630                  && cross[1][i]>=0
 631                  && cross[2][i]>=0
 632                  && cross[3][i]>=0) ) {
 633         crossFlag++;
 634       }
 635     }
 636
 637
 638     if( crossFlag==3 && noInterFlag[j]==1 ) {
 639       cubeInter[nSidesCrossed][0] = interceptx[j];
 640       cubeInter[nSidesCrossed][1] = intercepty[j];
 641       cubeInter[nSidesCrossed][2] = interceptz[j];
 642       nSidesCrossed++;
 643     }
 644
 645   } // End of loop over all four planes
 646
 647   m_RayStartCoordInMM[0] = cubeInter[0][0];
 648   m_RayStartCoordInMM[1] = cubeInter[0][1];
 649   m_RayStartCoordInMM[2] = cubeInter[0][2];
 650
 651   m_RayEndCoordInMM[0] = cubeInter[1][0];
 652   m_RayEndCoordInMM[1] = cubeInter[1][1];
 653   m_RayEndCoordInMM[2] = cubeInter[1][2];
 654
 655   if( nSidesCrossed >= 5 ) {
 656     std::cerr << "WARNING: No. of sides crossed equals: " << nSidesCrossed << std::endl;
 657   }
 658
 659   // If 'nSidesCrossed' is larger than 2, this means that the ray goes through
 660   // a corner of the volume and due to rounding errors, the ray is
 661   // deemed to go through more than two planes.  To obtain the correct
 662   // start and end positions we choose the two intercept values which
 663   // are furthest from each other.
 664
 665
 666   if( nSidesCrossed >= 3 ) {
 667     maxInterDist = 0;
 668     for( j=0; j<nSidesCrossed-1; j++ ) {
 669       for( k=j+1; k<nSidesCrossed; k++ ) {
 670         interDist = 0;
 671         for( i=0; i<3; i++ ) {
 672           interDist += (cubeInter[j][i] - cubeInter[k][i])*
 673                        (cubeInter[j][i] - cubeInter[k][i]);
 674         }
 675         if( interDist > maxInterDist ) {
 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     return true;
 693   } else {
 694     return false;
 695   }
 696 }
 697
 698
 699 /* -----------------------------------------------------------------------
 700    SetRay() - Set the position and direction of the ray
 701    ----------------------------------------------------------------------- */
 702
 703 template<class TInputImage, class TCoordRep>
 704 bool
 705 RayCastHelper<TInputImage, TCoordRep>
 706 ::SetRay(OutputPointType RayPosn, DirectionType RayDirn)
 707 {
 708
 709   // Store the position and direction of the ray
 710   typename TInputImage::SpacingType spacing=this->m_Image->GetSpacing();
 711   SizeType dim=this->m_Image->GetLargestPossibleRegion().GetSize();
 712
 713   // we need to translate the _center_ of the volume to the origin
 714   m_NumberOfVoxelsInX = dim[0];
 715   m_NumberOfVoxelsInY = dim[1];
 716   m_NumberOfVoxelsInZ = dim[2];
 717
 718   m_VoxelDimensionInX = spacing[0];
 719   m_VoxelDimensionInY = spacing[1];
 720   m_VoxelDimensionInZ = spacing[2];
 721
 722   // (Aviv) Incorporating a fix by Jian Wu
 723   // http://public.kitware.com/pipermail/insight-users/2006-March/017265.html
 724   m_CurrentRayPositionInMM[0] = RayPosn[0];
 725   m_CurrentRayPositionInMM[1] = RayPosn[1];
 726   m_CurrentRayPositionInMM[2] = RayPosn[2];
 727
 728   m_RayDirectionInMM[0] = RayDirn[0];
 729   m_RayDirectionInMM[1] = RayDirn[1];
 730   m_RayDirectionInMM[2] = RayDirn[2];
 731
 732   // Compute the ray path for this coordinate in mm
 733
 734   m_ValidRay = this->CalcRayIntercepts();
 735
 736   if (! m_ValidRay) {
 737     Reset();
 738     return false;
 739   }
 740
 741   // Convert the start and end coordinates of the ray to voxels
 742
 743   this->EndPointsInVoxels();
 744
 745   /* Calculate the ray direction vector in voxels and hence the voxel
 746      increment required to traverse the ray, and the number of
 747      interpolation points on the ray.
 748
 749      This routine also shifts the coordinate frame by half a voxel for
 750      two of the directional components (i.e. those lying within the
 751      planes of voxels being traversed). */
 752
 753   this->CalcDirnVector();
 754
 755
 756   /* Reduce the length of the ray until both start and end
 757      coordinates lie inside the volume. */
 758
 759   m_ValidRay = this->AdjustRayLength();
 760
 761   // Reset the iterator to the start of the ray.
 762
 763   Reset();
 764
 765   return m_ValidRay;
 766 }
 767
 768
 769 /* -----------------------------------------------------------------------
 770    EndPointsInVoxels() - Convert the endpoints to voxels
 771    ----------------------------------------------------------------------- */
 772
 773 template<class TInputImage, class TCoordRep>
 774 void
 775 RayCastHelper<TInputImage, TCoordRep>
 776 ::EndPointsInVoxels(void)
 777 {
 778   m_RayVoxelStartPosition[0] = m_RayStartCoordInMM[0]/m_VoxelDimensionInX;
 779   m_RayVoxelStartPosition[1] = m_RayStartCoordInMM[1]/m_VoxelDimensionInY;
 780   m_RayVoxelStartPosition[2] = m_RayStartCoordInMM[2]/m_VoxelDimensionInZ;
 781
 782   m_RayVoxelEndPosition[0] = m_RayEndCoordInMM[0]/m_VoxelDimensionInX;
 783   m_RayVoxelEndPosition[1] = m_RayEndCoordInMM[1]/m_VoxelDimensionInY;
 784   m_RayVoxelEndPosition[2] = m_RayEndCoordInMM[2]/m_VoxelDimensionInZ;
 785
 786 }
 787
 788
 789 /* -----------------------------------------------------------------------
 790    CalcDirnVector() - Calculate the incremental direction vector in voxels.
 791    ----------------------------------------------------------------------- */
 792
 793 template<class TInputImage, class TCoordRep>
 794 void
 795 RayCastHelper<TInputImage, TCoordRep>
 796 ::CalcDirnVector(void)
 797 {
 798   double xNum, yNum, zNum;
 799
 800   // Calculate the number of voxels in each direction
 801
 802   xNum = vcl_fabs(m_RayVoxelStartPosition[0] - m_RayVoxelEndPosition[0]);
 803   yNum = vcl_fabs(m_RayVoxelStartPosition[1] - m_RayVoxelEndPosition[1]);
 804   zNum = vcl_fabs(m_RayVoxelStartPosition[2] - m_RayVoxelEndPosition[2]);
 805
 806   // The direction iterated in is that with the greatest number of voxels
 807   // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 808
 809   // Iterate in X direction
 810
 811   if( (xNum >= yNum) && (xNum >= zNum) ) {
 812     if( m_RayVoxelStartPosition[0] < m_RayVoxelEndPosition[0] ) {
 813       m_VoxelIncrement[0] = 1;
 814
 815       m_VoxelIncrement[1]
 816       = (m_RayVoxelStartPosition[1]
 817          - m_RayVoxelEndPosition[1])/(m_RayVoxelStartPosition[0]
 818                                       - m_RayVoxelEndPosition[0]);
 819
 820       m_VoxelIncrement[2]
 821       = (m_RayVoxelStartPosition[2]
 822          - m_RayVoxelEndPosition[2])/(m_RayVoxelStartPosition[0]
 823                                       - m_RayVoxelEndPosition[0]);
 824     } else {
 825       m_VoxelIncrement[0] = -1;
 826
 827       m_VoxelIncrement[1]
 828       = -(m_RayVoxelStartPosition[1]
 829           - m_RayVoxelEndPosition[1])/(m_RayVoxelStartPosition[0]
 830                                        - m_RayVoxelEndPosition[0]);
 831
 832       m_VoxelIncrement[2]
 833       = -(m_RayVoxelStartPosition[2]
 834           - m_RayVoxelEndPosition[2])/(m_RayVoxelStartPosition[0]
 835                                        - m_RayVoxelEndPosition[0]);
 836     }
 837
 838     // This section is to alter the start position in order to
 839     // place the center of the voxels in there correct positions,
 840     // rather than placing them at the corner of voxels which is
 841     // what happens if this is not carried out.  The reason why
 842     // x has no -0.5 is because this is the direction we are going
 843     // to iterate in and therefore we wish to go from center to
 844     // center rather than finding the surrounding voxels.
 845
 846     m_RayVoxelStartPosition[1] += ( (int)m_RayVoxelStartPosition[0]
 847                                     - m_RayVoxelStartPosition[0])*m_VoxelIncrement[1]*m_VoxelIncrement[0]
 848                                   + 0.5*m_VoxelIncrement[1] - 0.5;
 849
 850     m_RayVoxelStartPosition[2] += ( (int)m_RayVoxelStartPosition[0]
 851                                     - m_RayVoxelStartPosition[0])*m_VoxelIncrement[2]*m_VoxelIncrement[0]
 852                                   + 0.5*m_VoxelIncrement[2] - 0.5;
 853
 854     m_RayVoxelStartPosition[0] = (int)m_RayVoxelStartPosition[0] + 0.5*m_VoxelIncrement[0];
 855
 856     m_TotalRayVoxelPlanes = (int)xNum;
 857
 858     m_TraversalDirection = TRANSVERSE_IN_X;
 859   }
 860
 861   // Iterate in Y direction
 862
 863   else if( (yNum >= xNum) && (yNum >= zNum) ) {
 864
 865     if( m_RayVoxelStartPosition[1] < m_RayVoxelEndPosition[1] ) {
 866       m_VoxelIncrement[1] = 1;
 867
 868       m_VoxelIncrement[0]
 869       = (m_RayVoxelStartPosition[0]
 870          - m_RayVoxelEndPosition[0])/(m_RayVoxelStartPosition[1]
 871                                       - m_RayVoxelEndPosition[1]);
 872
 873       m_VoxelIncrement[2]
 874       = (m_RayVoxelStartPosition[2]
 875          - m_RayVoxelEndPosition[2])/(m_RayVoxelStartPosition[1]
 876                                       - m_RayVoxelEndPosition[1]);
 877     } else {
 878       m_VoxelIncrement[1] = -1;
 879
 880       m_VoxelIncrement[0]
 881       = -(m_RayVoxelStartPosition[0]
 882           - m_RayVoxelEndPosition[0])/(m_RayVoxelStartPosition[1]
 883                                        - m_RayVoxelEndPosition[1]);
 884
 885       m_VoxelIncrement[2]
 886       = -(m_RayVoxelStartPosition[2]
 887           - m_RayVoxelEndPosition[2])/(m_RayVoxelStartPosition[1]
 888                                        - m_RayVoxelEndPosition[1]);
 889     }
 890
 891     m_RayVoxelStartPosition[0] += ( (int)m_RayVoxelStartPosition[1]
 892                                     - m_RayVoxelStartPosition[1])*m_VoxelIncrement[0]*m_VoxelIncrement[1]
 893                                   + 0.5*m_VoxelIncrement[0] - 0.5;
 894
 895     m_RayVoxelStartPosition[2] += ( (int)m_RayVoxelStartPosition[1]
 896                                     - m_RayVoxelStartPosition[1])*m_VoxelIncrement[2]*m_VoxelIncrement[1]
 897                                   + 0.5*m_VoxelIncrement[2] - 0.5;
 898
 899     m_RayVoxelStartPosition[1] = (int)m_RayVoxelStartPosition[1] + 0.5*m_VoxelIncrement[1];
 900
 901     m_TotalRayVoxelPlanes = (int)yNum;
 902
 903     m_TraversalDirection = TRANSVERSE_IN_Y;
 904   }
 905
 906   // Iterate in Z direction
 907
 908   else {
 909
 910     if( m_RayVoxelStartPosition[2] < m_RayVoxelEndPosition[2] ) {
 911       m_VoxelIncrement[2] = 1;
 912
 913       m_VoxelIncrement[0]
 914       = (m_RayVoxelStartPosition[0]
 915          - m_RayVoxelEndPosition[0])/(m_RayVoxelStartPosition[2]
 916                                       - m_RayVoxelEndPosition[2]);
 917
 918       m_VoxelIncrement[1]
 919       = (m_RayVoxelStartPosition[1]
 920          - m_RayVoxelEndPosition[1])/(m_RayVoxelStartPosition[2]
 921                                       - m_RayVoxelEndPosition[2]);
 922     } else {
 923       m_VoxelIncrement[2] = -1;
 924
 925       m_VoxelIncrement[0]
 926       = -(m_RayVoxelStartPosition[0]
 927           - m_RayVoxelEndPosition[0])/(m_RayVoxelStartPosition[2]
 928                                        - m_RayVoxelEndPosition[2]);
 929
 930       m_VoxelIncrement[1]
 931       = -(m_RayVoxelStartPosition[1]
 932           - m_RayVoxelEndPosition[1])/(m_RayVoxelStartPosition[2]
 933                                        - m_RayVoxelEndPosition[2]);
 934     }
 935
 936     m_RayVoxelStartPosition[0] += ( (int)m_RayVoxelStartPosition[2]
 937                                     - m_RayVoxelStartPosition[2])*m_VoxelIncrement[0]*m_VoxelIncrement[2]
 938                                   + 0.5*m_VoxelIncrement[0] - 0.5;
 939
 940     m_RayVoxelStartPosition[1] += ( (int)m_RayVoxelStartPosition[2]
 941                                     - m_RayVoxelStartPosition[2])*m_VoxelIncrement[1]*m_VoxelIncrement[2]
 942                                   + 0.5*m_VoxelIncrement[1] - 0.5;
 943
 944     m_RayVoxelStartPosition[2] = (int)m_RayVoxelStartPosition[2] + 0.5*m_VoxelIncrement[2];
 945
 946     m_TotalRayVoxelPlanes = (int)zNum;
 947
 948     m_TraversalDirection = TRANSVERSE_IN_Z;
 949   }
 950 }
 951
 952
 953 /* -----------------------------------------------------------------------
 954    AdjustRayLength() - Ensure that the ray lies within the volume
 955    ----------------------------------------------------------------------- */
 956
 957 template<class TInputImage, class TCoordRep>
 958 bool
 959 RayCastHelper<TInputImage, TCoordRep>
 960 ::AdjustRayLength(void)
 961 {
 962   bool startOK, endOK;
 963
 964   int Istart[3];
 965   int Idirn[3];
 966
 967   if (m_TraversalDirection == TRANSVERSE_IN_X) {
 968     Idirn[0] = 0;
 969     Idirn[1] = 1;
 970     Idirn[2] = 1;
 971   } else if (m_TraversalDirection == TRANSVERSE_IN_Y) {
 972     Idirn[0] = 1;
 973     Idirn[1] = 0;
 974     Idirn[2] = 1;
 975   } else if (m_TraversalDirection == TRANSVERSE_IN_Z) {
 976     Idirn[0] = 1;
 977     Idirn[1] = 1;
 978     Idirn[2] = 0;
 979   } else {
 980     itk::ExceptionObject err(__FILE__, __LINE__);
 981     err.SetLocation( ITK_LOCATION );
 982     err.SetDescription( "The ray traversal direction is unset "
 983                         "- AdjustRayLength().");
 984     throw err;
 985     return false;
 986   }
 987
 988
 989   do {
 990
 991     startOK = false;
 992     endOK = false;
 993
 994     Istart[0] = (int) vcl_floor(m_RayVoxelStartPosition[0]);
 995     Istart[1] = (int) vcl_floor(m_RayVoxelStartPosition[1]);
 996     Istart[2] = (int) vcl_floor(m_RayVoxelStartPosition[2]);
 997
 998     if( (Istart[0] >= 0) && (Istart[0] + Idirn[0] < m_NumberOfVoxelsInX) &&
 999         (Istart[1] >= 0) && (Istart[1] + Idirn[1] < m_NumberOfVoxelsInY) &&
1000         (Istart[2] >= 0) && (Istart[2] + Idirn[2] < m_NumberOfVoxelsInZ) ) {
1001
1002       startOK = true;
1003     } else {
1004       m_RayVoxelStartPosition[0] += m_VoxelIncrement[0];
1005       m_RayVoxelStartPosition[1] += m_VoxelIncrement[1];
1006       m_RayVoxelStartPosition[2] += m_VoxelIncrement[2];
1007
1008       m_TotalRayVoxelPlanes--;
1009     }
1010
1011     Istart[0] = (int) vcl_floor(m_RayVoxelStartPosition[0]
1012                                 + m_TotalRayVoxelPlanes*m_VoxelIncrement[0]);
1013
1014     Istart[1] = (int) vcl_floor(m_RayVoxelStartPosition[1]
1015                                 + m_TotalRayVoxelPlanes*m_VoxelIncrement[1]);
1016
1017     Istart[2] = (int) vcl_floor(m_RayVoxelStartPosition[2]
1018                                 + m_TotalRayVoxelPlanes*m_VoxelIncrement[2]);
1019
1020     if( (Istart[0] >= 0) && (Istart[0] + Idirn[0] < m_NumberOfVoxelsInX) &&
1021         (Istart[1] >= 0) && (Istart[1] + Idirn[1] < m_NumberOfVoxelsInY) &&
1022         (Istart[2] >= 0) && (Istart[2] + Idirn[2] < m_NumberOfVoxelsInZ) ) {
1023
1024       endOK = true;
1025     } else {
1026       m_TotalRayVoxelPlanes--;
1027     }
1028
1029   } while ( (! (startOK && endOK)) && (m_TotalRayVoxelPlanes > 1) );
1030
1031
1032   return (startOK && endOK);
1033 }
1034
1035
1036 /* -----------------------------------------------------------------------
1037    Reset() - Reset the iterator to the start of the ray.
1038    ----------------------------------------------------------------------- */
1039
1040 template<class TInputImage, class TCoordRep>
1041 void
1042 RayCastHelper<TInputImage, TCoordRep>
1043 ::Reset(void)
1044 {
1045   int i;
1046
1047   m_NumVoxelPlanesTraversed = -1;
1048
1049   // If this is a valid ray...
1050
1051   if (m_ValidRay) {
1052     for (i=0; i<3; i++) {
1053       m_Position3Dvox[i] = m_RayVoxelStartPosition[i];
1054     }
1055     this->InitialiseVoxelPointers();
1056   }
1057
1058   // otherwise set parameters to zero
1059
1060   else {
1061     for (i=0; i<3; i++) {
1062       m_RayVoxelStartPosition[i] = 0.;
1063     }
1064     for (i=0; i<3; i++) {
1065       m_RayVoxelEndPosition[i] = 0.;
1066     }
1067     for (i=0; i<3; i++) {
1068       m_VoxelIncrement[i] = 0.;
1069     }
1070     m_TraversalDirection = UNDEFINED_DIRECTION;
1071
1072     m_TotalRayVoxelPlanes = 0;
1073
1074     for (i=0; i<4; i++) {
1075       m_RayIntersectionVoxels[i] = 0;
1076     }
1077     for (i=0; i<3; i++) {
1078       m_RayIntersectionVoxelIndex[i] = 0;
1079     }
1080   }
1081 }
1082
1083
1084 /* -----------------------------------------------------------------------
1085    InitialiseVoxelPointers() - Obtain pointers to the first four voxels
1086    ----------------------------------------------------------------------- */
1087
1088 template<class TInputImage, class TCoordRep>
1089 void
1090 RayCastHelper<TInputImage, TCoordRep>
1091 ::InitialiseVoxelPointers(void)
1092 {
1093   IndexType index;
1094
1095   int Ix, Iy, Iz;
1096
1097   Ix = (int)(m_RayVoxelStartPosition[0]);
1098   Iy = (int)(m_RayVoxelStartPosition[1]);
1099   Iz = (int)(m_RayVoxelStartPosition[2]);
1100
1101   m_RayIntersectionVoxelIndex[0] = Ix;
1102   m_RayIntersectionVoxelIndex[1] = Iy;
1103   m_RayIntersectionVoxelIndex[2] = Iz;
1104
1105   switch( m_TraversalDirection ) {
1106   case TRANSVERSE_IN_X: {
1107
1108     if( (Ix >= 0) && (Ix     < m_NumberOfVoxelsInX) &&
1109         (Iy >= 0) && (Iy + 1 < m_NumberOfVoxelsInY) &&
1110         (Iz >= 0) && (Iz + 1 < m_NumberOfVoxelsInZ)) {
1111
1112       index[0]=Ix;
1113       index[1]=Iy;
1114       index[2]=Iz;
1115       m_RayIntersectionVoxels[0]
1116       = this->m_Image->GetBufferPointer() + this->m_Image->ComputeOffset(index);
1117
1118       index[0]=Ix;
1119       index[1]=Iy+1;
1120       index[2]=Iz;
1121       m_RayIntersectionVoxels[1]
1122       = ( this->m_Image->GetBufferPointer() + this->m_Image->ComputeOffset(index) );
1123
1124       index[0]=Ix;
1125       index[1]=Iy;
1126       index[2]=Iz+1;
1127       m_RayIntersectionVoxels[2]
1128       = ( this->m_Image->GetBufferPointer() + this->m_Image->ComputeOffset(index) );
1129
1130       index[0]=Ix;
1131       index[1]=Iy+1;
1132       index[2]=Iz+1;
1133       m_RayIntersectionVoxels[3]
1134       = ( this->m_Image->GetBufferPointer() + this->m_Image->ComputeOffset(index) );
1135     } else {
1136       m_RayIntersectionVoxels[0] =
1137         m_RayIntersectionVoxels[1] =
1138           m_RayIntersectionVoxels[2] =
1139             m_RayIntersectionVoxels[3] = NULL;
1140     }
1141     break;
1142   }
1143
1144   case TRANSVERSE_IN_Y: {
1145
1146     if( (Ix >= 0) && (Ix + 1 < m_NumberOfVoxelsInX) &&
1147         (Iy >= 0) && (Iy     < m_NumberOfVoxelsInY) &&
1148         (Iz >= 0) && (Iz + 1 < m_NumberOfVoxelsInZ)) {
1149
1150       index[0]=Ix;
1151       index[1]=Iy;
1152       index[2]=Iz;
1153       m_RayIntersectionVoxels[0] = ( this->m_Image->GetBufferPointer()
1154                                      + this->m_Image->ComputeOffset(index) );
1155
1156       index[0]=Ix+1;
1157       index[1]=Iy;
1158       index[2]=Iz;
1159       m_RayIntersectionVoxels[1] = ( this->m_Image->GetBufferPointer()
1160                                      + this->m_Image->ComputeOffset(index) );
1161
1162       index[0]=Ix;
1163       index[1]=Iy;
1164       index[2]=Iz+1;
1165       m_RayIntersectionVoxels[2] = ( this->m_Image->GetBufferPointer()
1166                                      + this->m_Image->ComputeOffset(index) );
1167
1168       index[0]=Ix+1;
1169       index[1]=Iy;
1170       index[2]=Iz+1;
1171       m_RayIntersectionVoxels[3] = ( this->m_Image->GetBufferPointer()
1172                                      + this->m_Image->ComputeOffset(index) );
1173     } else {
1174       m_RayIntersectionVoxels[0]
1175       = m_RayIntersectionVoxels[1]
1176         = m_RayIntersectionVoxels[2]
1177           = m_RayIntersectionVoxels[3] = NULL;
1178     }
1179     break;
1180   }
1181
1182   case TRANSVERSE_IN_Z: {
1183
1184     if( (Ix >= 0) && (Ix + 1 < m_NumberOfVoxelsInX)   &&
1185         (Iy >= 0) && (Iy + 1 < m_NumberOfVoxelsInY) &&
1186         (Iz >= 0) && (Iz     < m_NumberOfVoxelsInZ)) {
1187
1188       index[0]=Ix;
1189       index[1]=Iy;
1190       index[2]=Iz;
1191       m_RayIntersectionVoxels[0] = ( this->m_Image->GetBufferPointer()
1192                                      + this->m_Image->ComputeOffset(index) );
1193
1194       index[0]=Ix+1;
1195       index[1]=Iy;
1196       index[2]=Iz;
1197       m_RayIntersectionVoxels[1] = ( this->m_Image->GetBufferPointer()
1198                                      + this->m_Image->ComputeOffset(index) );
1199
1200       index[0]=Ix;
1201       index[1]=Iy+1;
1202       index[2]=Iz;
1203       m_RayIntersectionVoxels[2] = ( this->m_Image->GetBufferPointer()
1204                                      + this->m_Image->ComputeOffset(index) );
1205
1206       index[0]=Ix+1;
1207       index[1]=Iy+1;
1208       index[2]=Iz;
1209       m_RayIntersectionVoxels[3] = ( this->m_Image->GetBufferPointer()
1210                                      + this->m_Image->ComputeOffset(index) );
1211
1212     } else {
1213       m_RayIntersectionVoxels[0]
1214       = m_RayIntersectionVoxels[1]
1215         = m_RayIntersectionVoxels[2]
1216           = m_RayIntersectionVoxels[3] = NULL;
1217     }
1218     break;
1219   }
1220
1221   default: {
1222     itk::ExceptionObject err(__FILE__, __LINE__);
1223     err.SetLocation( ITK_LOCATION );
1224     err.SetDescription( "The ray traversal direction is unset "
1225                         "- InitialiseVoxelPointers().");
1226     throw err;
1227     return;
1228   }
1229   }
1230 }
1231
1232 /* -----------------------------------------------------------------------
1233    IncrementVoxelPointers() - Increment the voxel pointers
1234    ----------------------------------------------------------------------- */
1235
1236 template<class TInputImage, class TCoordRep>
1237 void
1238 RayCastHelper<TInputImage, TCoordRep>
1239 ::IncrementVoxelPointers(void)
1240 {
1241   double xBefore = m_Position3Dvox[0];
1242   double yBefore = m_Position3Dvox[1];
1243   double zBefore = m_Position3Dvox[2];
1244
1245   m_Position3Dvox[0] += m_VoxelIncrement[0];
1246   m_Position3Dvox[1] += m_VoxelIncrement[1];
1247   m_Position3Dvox[2] += m_VoxelIncrement[2];
1248
1249   int dx = ((int) m_Position3Dvox[0]) - ((int) xBefore);
1250   int dy = ((int) m_Position3Dvox[1]) - ((int) yBefore);
1251   int dz = ((int) m_Position3Dvox[2]) - ((int) zBefore);
1252
1253   m_RayIntersectionVoxelIndex[0] += dx;
1254   m_RayIntersectionVoxelIndex[1] += dy;
1255   m_RayIntersectionVoxelIndex[2] += dz;
1256
1257   int totalRayVoxelPlanes
1258   = dx + dy*m_NumberOfVoxelsInX + dz*m_NumberOfVoxelsInX*m_NumberOfVoxelsInY;
1259
1260   m_RayIntersectionVoxels[0] += totalRayVoxelPlanes;
1261   m_RayIntersectionVoxels[1] += totalRayVoxelPlanes;
1262   m_RayIntersectionVoxels[2] += totalRayVoxelPlanes;
1263   m_RayIntersectionVoxels[3] += totalRayVoxelPlanes;
1264 }
1265
1266
1267 /* -----------------------------------------------------------------------
1268    GetCurrentIntensity() - Get the intensity of the current ray point.
1269    ----------------------------------------------------------------------- */
1270
1271 template<class TInputImage, class TCoordRep>
1272 double
1273 RayCastHelper<TInputImage, TCoordRep>
1274 ::GetCurrentIntensity(void) const
1275 {
1276   double a, b, c, d;
1277   double y, z;
1278
1279   if (! m_ValidRay) {
1280     return 0;
1281   }
1282   a = (double) (*m_RayIntersectionVoxels[0]);
1283   b = (double) (*m_RayIntersectionVoxels[1] - a);
1284   c = (double) (*m_RayIntersectionVoxels[2] - a);
1285   d = (double) (*m_RayIntersectionVoxels[3] - a - b - c);
1286
1287   switch( m_TraversalDirection ) {
1288   case TRANSVERSE_IN_X: {
1289     y = m_Position3Dvox[1] - vcl_floor(m_Position3Dvox[1]);
1290     z = m_Position3Dvox[2] - vcl_floor(m_Position3Dvox[2]);
1291     break;
1292   }
1293   case TRANSVERSE_IN_Y: {
1294     y = m_Position3Dvox[0] - vcl_floor(m_Position3Dvox[0]);
1295     z = m_Position3Dvox[2] - vcl_floor(m_Position3Dvox[2]);
1296     break;
1297   }
1298   case TRANSVERSE_IN_Z: {
1299     y = m_Position3Dvox[0] - vcl_floor(m_Position3Dvox[0]);
1300     z = m_Position3Dvox[1] - vcl_floor(m_Position3Dvox[1]);
1301     break;
1302   }
1303   default: {
1304     itk::ExceptionObject err(__FILE__, __LINE__);
1305     err.SetLocation( ITK_LOCATION );
1306     err.SetDescription( "The ray traversal direction is unset "
1307                         "- GetCurrentIntensity().");
1308     throw err;
1309     return 0;
1310   }
1311   }
1312
1313   return a + b*y + c*z + d*y*z;
1314 }
1315
1316 /* -----------------------------------------------------------------------
1317    IncrementIntensities() - Increment the intensities of the current ray point
1318    ----------------------------------------------------------------------- */
1319
1320 template<class TInputImage, class TCoordRep>
1321 void
1322 RayCastHelper<TInputImage, TCoordRep>
1323 ::IncrementIntensities(double increment)
1324 {
1325   short inc = (short) vcl_floor(increment + 0.5);
1326
1327   if (! m_ValidRay) {
1328     return;
1329   }
1330   *m_RayIntersectionVoxels[0] += inc;
1331   *m_RayIntersectionVoxels[1] += inc;
1332   *m_RayIntersectionVoxels[2] += inc;
1333   *m_RayIntersectionVoxels[3] += inc;
1334
1335   return;
1336 }
1337
1338
1339 /* -----------------------------------------------------------------------
1340    IntegrateAboveThreshold() - Integrate intensities above a threshold.
1341    ----------------------------------------------------------------------- */
1342
1343 template<class TInputImage, class TCoordRep>
1344 bool
1345 RayCastHelper<TInputImage, TCoordRep>
1346 ::IntegrateAboveThreshold(double &integral, double threshold)
1347 {
1348   double intensity;
1349 //  double posn3D_x, posn3D_y, posn3D_z;
1350
1351   integral = 0.;
1352
1353   // Check if this is a valid ray
1354
1355   if (! m_ValidRay) {
1356     return false;
1357   }
1358   /* Step along the ray as quickly as possible
1359      integrating the interpolated intensities. */
1360
1361   for (m_NumVoxelPlanesTraversed=0;
1362        m_NumVoxelPlanesTraversed<m_TotalRayVoxelPlanes;
1363        m_NumVoxelPlanesTraversed++) {
1364     intensity = this->GetCurrentIntensity();
1365
1366     if (intensity > threshold) {
1367       integral += intensity - threshold;
1368     }
1369     this->IncrementVoxelPointers();
1370   }
1371
1372   /* The ray passes through the volume one plane of voxels at a time,
1373      however, if its moving diagonally the ray points will be further
1374      apart so account for this by scaling by the distance moved. */
1375
1376   integral *= this->GetRayPointSpacing();
1377
1378   return true;
1379 }
1380
1381 /* -----------------------------------------------------------------------
1382    ZeroState() - Set the default (zero) state of the object
1383    ----------------------------------------------------------------------- */
1384
1385 template<class TInputImage, class TCoordRep>
1386 void
1387 RayCastHelper<TInputImage, TCoordRep>
1388 ::ZeroState()
1389 {
1390   int i;
1391
1392   m_ValidRay = false;
1393
1394   m_NumberOfVoxelsInX = 0;
1395   m_NumberOfVoxelsInY = 0;
1396   m_NumberOfVoxelsInZ = 0;
1397
1398   m_VoxelDimensionInX = 0;
1399   m_VoxelDimensionInY = 0;
1400   m_VoxelDimensionInZ = 0;
1401
1402   for (i=0; i<3; i++) {
1403     m_CurrentRayPositionInMM[i] = 0.;
1404   }
1405   for (i=0; i<3; i++) {
1406     m_RayDirectionInMM[i] = 0.;
1407   }
1408   for (i=0; i<3; i++) {
1409     m_RayVoxelStartPosition[i] = 0.;
1410   }
1411   for (i=0; i<3; i++) {
1412     m_RayVoxelEndPosition[i] = 0.;
1413   }
1414   for (i=0; i<3; i++) {
1415     m_VoxelIncrement[i] = 0.;
1416   }
1417   m_TraversalDirection = UNDEFINED_DIRECTION;
1418
1419   m_TotalRayVoxelPlanes = 0;
1420   m_NumVoxelPlanesTraversed = -1;
1421
1422   for (i=0; i<4; i++) {
1423     m_RayIntersectionVoxels[i] = 0;
1424   }
1425   for (i=0; i<3; i++) {
1426     m_RayIntersectionVoxelIndex[i] = 0;
1427   }
1428 }
1429 }; // end of anonymous namespace
1430
1431
1432 namespace itk
1433 {
1434
1435 /**************************************************************************
1436  *
1437  *
1438  * Rest of this code is the actual RayCastInterpolateImageFunctionWithOrigin
1439  * class
1440  *
1441  *
1442  **************************************************************************/
1443
1444 /* -----------------------------------------------------------------------
1445    Constructor
1446    ----------------------------------------------------------------------- */
1447
1448 template<class TInputImage, class TCoordRep>
1449 RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
1450 ::RayCastInterpolateImageFunctionWithOrigin()
1451 {
1452   m_Threshold = 0.;
1453
1454   m_FocalPoint[0] = 0.;
1455   m_FocalPoint[1] = 0.;
1456   m_FocalPoint[2] = 0.;
1457 }
1458
1459
1460 /* -----------------------------------------------------------------------
1461    PrintSelf
1462    ----------------------------------------------------------------------- */
1463
1464 template<class TInputImage, class TCoordRep>
1465 void
1466 RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
1467 ::PrintSelf(std::ostream& os, Indent indent) const
1468 {
1469   this->Superclass::PrintSelf(os,indent);
1470
1471   os << indent << "Threshold: " << m_Threshold << std::endl;
1472   os << indent << "FocalPoint: " << m_FocalPoint << std::endl;
1473   os << indent << "Transform: " << m_Transform.GetPointer() << std::endl;
1474   os << indent << "Interpolator: " << m_Interpolator.GetPointer() << std::endl;
1475
1476 }
1477
1478 /* -----------------------------------------------------------------------
1479    Evaluate at image index position
1480    ----------------------------------------------------------------------- */
1481
1482 template<class TInputImage, class TCoordRep>
1483 typename RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
1484 ::OutputType
1485 RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
1486 ::Evaluate( const PointType& point ) const
1487 {
1488   double integral = 0;
1489
1490   OutputPointType transformedFocalPoint
1491   = m_Transform->TransformPoint( m_FocalPoint );
1492
1493   DirectionType direction = transformedFocalPoint - point;
1494
1495   RayCastHelper<TInputImage, TCoordRep> ray;
1496   ray.SetImage( this->m_Image );
1497   ray.ZeroState();
1498   ray.Initialise();
1499
1500   // (Aviv) Added support for images with non-zero origin
1501   // ray.SetRay(point, direction);
1502   ray.SetRay(point - this->m_Image->GetOrigin().GetVectorFromOrigin(), direction);
1503   ray.IntegrateAboveThreshold(integral, m_Threshold);
1504
1505   return ( static_cast<OutputType>( integral ));
1506 }
1507
1508 template<class TInputImage, class TCoordRep>
1509 typename RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
1510 ::OutputType
1511 RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
1512 ::EvaluateAtContinuousIndex( const ContinuousIndexType& index ) const
1513 {
1514   OutputPointType point;
1515   this->m_Image->TransformContinuousIndexToPhysicalPoint(index, point);
1516
1517   return this->Evaluate( point );
1518 }
1519
1520 } // namespace itk
1521
1522
1523 #endif