itk/clitkRayCastInterpolateImageFunctionWithOrigin.txx

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