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