itk/itkRayCastInterpolateImageFunctionWithOrigin.txx

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