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