itk/itkRayCastInterpolateImageFunctionWithOrigin.txx

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