1 /*=========================================================================
3 Program: Insight Segmentation & Registration Toolkit
4 Module: $RCSfile: itkRayCastInterpolateImageFunctionWithOrigin.txx,v $
6 Date: $Date: 2010/01/06 13:32:01 $
7 Version: $Revision: 1.1 $
9 Copyright (c) Insight Software Consortium. All rights reserved.
10 See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
12 This software is distributed WITHOUT ANY WARRANTY; without even
13 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
14 PURPOSE. See the above copyright notices for more information.
16 =========================================================================*/
18 #ifndef __itkRayCastInterpolateImageFunctionWithOrigin_txx
19 #define __itkRayCastInterpolateImageFunctionWithOrigin_txx
21 #include "itkRayCastInterpolateImageFunctionWithOrigin.h"
23 #include "vnl/vnl_math.h"
26 // Put the helper class in an anonymous namespace so that it is not
27 // exposed to the user
30 /** \class Helper class to maintain state when casting a ray.
31 * This helper class keeps the RayCastInterpolateImageFunctionWithOrigin thread safe.
33 template <class TInputImage, class TCoordRep = float>
37 /** Constants for the image dimensions */
38 itkStaticConstMacro(InputImageDimension, unsigned int,
39 TInputImage::ImageDimension);
42 * Type of the Transform Base class
43 * The fixed image should be a 3D image
45 typedef itk::Transform<TCoordRep,3,3> TransformType;
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;
53 typedef typename TInputImage::SizeType SizeType;
54 typedef itk::Vector<TCoordRep, 3> DirectionType;
55 typedef itk::Point<TCoordRep, 3> PointType;
57 typedef TInputImage InputImageType;
58 typedef typename InputImageType::PixelType PixelType;
59 typedef typename InputImageType::IndexType IndexType;
64 void SetImage(const InputImageType *input)
70 * Initialise the ray using the position and direction of a line.
72 * \param RayPosn The position of the ray in 3D (mm).
73 * \param RayDirn The direction of the ray in 3D (mm).
75 * \return True if this is a valid ray.
77 bool SetRay(OutputPointType RayPosn, DirectionType RayDirn);
81 * Integrate the interpolated intensities along the ray and
84 * This routine can be called after instantiating the ray and
85 * calling SetProjectionCoord2D() or Reset(). It may then be called
86 * as many times thereafter for different 2D projection
89 * \param integral The integrated intensities along the ray.
91 * \return True if a valid ray was specified.
93 bool Integrate(double &integral)
95 return IntegrateAboveThreshold(integral, 0);
100 * Integrate the interpolated intensities above a given threshold,
101 * along the ray and return the result.
103 * This routine can be called after instantiating the ray and
104 * calling SetProjectionCoord2D() or Reset(). It may then be called
105 * as many times thereafter for different 2D projection
108 * \param integral The integrated intensities along the ray.
109 * \param threshold The integration threshold [default value: 0]
111 * \return True if a valid ray was specified.
113 bool IntegrateAboveThreshold(double &integral, double threshold);
116 * Increment each of the intensities of the 4 planar voxels
117 * surrounding the current ray point.
119 * \parameter increment Intensity increment for each of the current 4 voxels
121 void IncrementIntensities(double increment=1);
123 /// Reset the iterator to the start of the ray.
126 /// Return the interpolated intensity of the current ray point.
127 double GetCurrentIntensity(void) const;
129 /// Return the ray point spacing in mm
130 double GetRayPointSpacing(void) const {
131 typename InputImageType::SpacingType spacing=this->m_Image->GetSpacing();
134 return vcl_sqrt(m_VoxelIncrement[0]*spacing[0]*m_VoxelIncrement[0]*spacing[0]
135 + m_VoxelIncrement[1]*spacing[1]*m_VoxelIncrement[1]*spacing[1]
136 + m_VoxelIncrement[2]*spacing[2]*m_VoxelIncrement[2]*spacing[2] );
141 /// Set the initial zero state of the object
144 /// Initialise the object
145 void Initialise(void);
148 /// Calculate the endpoint coordinats of the ray in voxels.
149 void EndPointsInVoxels(void);
152 * Calculate the incremental direction vector in voxels, 'dVoxel',
153 * required to traverse the ray.
155 void CalcDirnVector(void);
158 * Reduce the length of the ray until both start and end
159 * coordinates lie inside the volume.
161 * \return True if a valid ray has been, false otherwise.
163 bool AdjustRayLength(void);
166 * Obtain pointers to the four voxels surrounding the point where the ray
169 void InitialiseVoxelPointers(void);
171 /// Increment the voxel pointers surrounding the current point on the ray.
172 void IncrementVoxelPointers(void);
174 /// Record volume dimensions and resolution
175 void RecordVolumeDimensions(void);
177 /// Define the corners of the volume
178 void DefineCorners(void);
181 * Calculate the planes which define the volume.
183 * Member function to calculate the equations of the planes of 4 of
184 * the sides of the volume, calculate the positions of the 8 corners
185 * of the volume in mm in World, also calculate the values of the
186 * slopes of the lines which go to make up the volume( defined as
187 * lines in cube x,y,z dirn and then each of these lines has a slope
188 * in the world x,y,z dirn [3]) and finally also to return the length
189 * of the sides of the lines in mm.
191 void CalcPlanesAndCorners(void);
194 * Calculate the ray intercepts with the volume.
196 * See where the ray cuts the volume, check that truncation does not occur,
197 * if not, then start ray where it first intercepts the volume and set
198 * x_max to be where it leaves the volume.
200 * \return True if a valid ray has been specified, false otherwise.
202 bool CalcRayIntercepts(void);
205 * The ray is traversed by stepping in the axial direction
206 * that enables the greatest number of planes in the volume to be
210 UNDEFINED_DIRECTION=0, //!< Undefined
211 TRANSVERSE_IN_X, //!< x
212 TRANSVERSE_IN_Y, //!< y
213 TRANSVERSE_IN_Z, //!< z
215 } TraversalDirection;
217 // Cache the image in the structure. Skip the smart pointer for
218 // efficiency. This inner class will go in/out of scope with every
219 // call to Evaluate()
220 const InputImageType *m_Image;
222 /// Flag indicating whether the current ray is valid
226 * The start position of the ray in voxels.
228 * NB. Two of the components of this coordinate (i.e. those lying within
229 * the planes of voxels being traversed) will be shifted by half a
230 * voxel. This enables indices of the neighbouring voxels within the plane
231 * to be determined by simply casting to 'int' and optionally adding 1.
233 double m_RayVoxelStartPosition[3];
236 * The end coordinate of the ray in voxels.
238 * NB. Two of the components of this coordinate (i.e. those lying within
239 * the planes of voxels being traversed) will be shifted by half a
240 * voxel. This enables indices of the neighbouring voxels within the plane
241 * to be determined by simply casting to 'int' and optionally adding 1.
243 double m_RayVoxelEndPosition[3];
247 * The current coordinate on the ray in voxels.
249 * NB. Two of the components of this coordinate (i.e. those lying within
250 * the planes of voxels being traversed) will be shifted by half a
251 * voxel. This enables indices of the neighbouring voxels within the plane
252 * to be determined by simply casting to 'int' and optionally adding 1.
254 double m_Position3Dvox[3];
256 /** The incremental direction vector of the ray in voxels. */
257 double m_VoxelIncrement[3];
259 /// The direction in which the ray is incremented thorough the volume (x, y or z).
260 TraversalDirection m_TraversalDirection;
262 /// The total number of planes of voxels traversed by the ray.
263 int m_TotalRayVoxelPlanes;
265 /// The current number of planes of voxels traversed by the ray.
266 int m_NumVoxelPlanesTraversed;
268 /// Pointers to the current four voxels surrounding the ray's trajectory.
269 const PixelType *m_RayIntersectionVoxels[4];
272 * The voxel coordinate of the bottom-left voxel of the current
273 * four voxels surrounding the ray's trajectory.
275 int m_RayIntersectionVoxelIndex[3];
277 /// The dimension in voxels of the 3D volume in along the x axis
278 int m_NumberOfVoxelsInX;
279 /// The dimension in voxels of the 3D volume in along the y axis
280 int m_NumberOfVoxelsInY;
281 /// The dimension in voxels of the 3D volume in along the z axis
282 int m_NumberOfVoxelsInZ;
284 /// Voxel dimension in x
285 double m_VoxelDimensionInX;
286 /// Voxel dimension in y
287 double m_VoxelDimensionInY;
288 /// Voxel dimension in z
289 double m_VoxelDimensionInZ;
291 /// The coordinate of the point at which the ray enters the volume in mm.
292 double m_RayStartCoordInMM[3];
293 /// The coordinate of the point at which the ray exits the volume in mm.
294 double m_RayEndCoordInMM[3];
298 Planes which define the boundary of the volume in mm
299 (six planes and four parameters: Ax+By+Cz+D). */
300 double m_BoundingPlane[6][4];
301 /// The eight corners of the volume (x,y,z coordinates for each).
302 double m_BoundingCorner[8][3];
304 /// The position of the ray
305 double m_CurrentRayPositionInMM[3];
307 /// The direction of the ray
308 double m_RayDirectionInMM[3];
312 /* -----------------------------------------------------------------------
313 Initialise() - Initialise the object
314 ----------------------------------------------------------------------- */
316 template<class TInputImage, class TCoordRep>
318 RayCastHelper<TInputImage, TCoordRep>
321 // Save the dimensions of the volume and calculate the bounding box
322 this->RecordVolumeDimensions();
324 // Calculate the planes and corners which define the volume.
325 this->DefineCorners();
326 this->CalcPlanesAndCorners();
330 /* -----------------------------------------------------------------------
331 RecordVolumeDimensions() - Record volume dimensions and resolution
332 ----------------------------------------------------------------------- */
334 template<class TInputImage, class TCoordRep>
336 RayCastHelper<TInputImage, TCoordRep>
337 ::RecordVolumeDimensions(void)
339 typename InputImageType::SpacingType spacing=this->m_Image->GetSpacing();
340 SizeType dim=this->m_Image->GetLargestPossibleRegion().GetSize();
342 m_NumberOfVoxelsInX = dim[0];
343 m_NumberOfVoxelsInY = dim[1];
344 m_NumberOfVoxelsInZ = dim[2];
346 m_VoxelDimensionInX = spacing[0];
347 m_VoxelDimensionInY = spacing[1];
348 m_VoxelDimensionInZ = spacing[2];
352 /* -----------------------------------------------------------------------
353 DefineCorners() - Define the corners of the volume
354 ----------------------------------------------------------------------- */
356 template<class TInputImage, class TCoordRep>
358 RayCastHelper<TInputImage, TCoordRep>
359 ::DefineCorners(void)
361 // Define corner positions as if at the origin
363 m_BoundingCorner[0][0] =
364 m_BoundingCorner[1][0] =
365 m_BoundingCorner[2][0] =
366 m_BoundingCorner[3][0] = 0;
368 m_BoundingCorner[4][0] =
369 m_BoundingCorner[5][0] =
370 m_BoundingCorner[6][0] =
371 m_BoundingCorner[7][0] = m_VoxelDimensionInX*m_NumberOfVoxelsInX;
373 m_BoundingCorner[1][1] =
374 m_BoundingCorner[3][1] =
375 m_BoundingCorner[5][1] =
376 m_BoundingCorner[7][1] = m_VoxelDimensionInY*m_NumberOfVoxelsInY;
378 m_BoundingCorner[0][1] =
379 m_BoundingCorner[2][1] =
380 m_BoundingCorner[4][1] =
381 m_BoundingCorner[6][1] = 0;
383 m_BoundingCorner[0][2] =
384 m_BoundingCorner[1][2] =
385 m_BoundingCorner[4][2] =
386 m_BoundingCorner[5][2] =
387 m_VoxelDimensionInZ*m_NumberOfVoxelsInZ;
389 m_BoundingCorner[2][2] =
390 m_BoundingCorner[3][2] =
391 m_BoundingCorner[6][2] =
392 m_BoundingCorner[7][2] = 0;
396 /* -----------------------------------------------------------------------
397 CalcPlanesAndCorners() - Calculate the planes and corners of the volume.
398 ----------------------------------------------------------------------- */
400 template<class TInputImage, class TCoordRep>
402 RayCastHelper<TInputImage, TCoordRep>
403 ::CalcPlanesAndCorners(void)
408 // find the equations of the planes
410 int c1=0, c2=0, c3=0;
413 { // loop around for planes
415 { // which corners to take
437 double line1x, line1y, line1z;
438 double line2x, line2y, line2z;
440 // lines from one corner to another in x,y,z dirns
441 line1x = m_BoundingCorner[c1][0] - m_BoundingCorner[c2][0];
442 line2x = m_BoundingCorner[c1][0] - m_BoundingCorner[c3][0];
444 line1y = m_BoundingCorner[c1][1] - m_BoundingCorner[c2][1];
445 line2y = m_BoundingCorner[c1][1] - m_BoundingCorner[c3][1];
447 line1z = m_BoundingCorner[c1][2] - m_BoundingCorner[c2][2];
448 line2z = m_BoundingCorner[c1][2] - m_BoundingCorner[c3][2];
452 // take cross product
453 A = line1y*line2z - line2y*line1z;
454 B = line2x*line1z - line1x*line2z;
455 C = line1x*line2y - line2x*line1y;
458 D = -( A*m_BoundingCorner[c1][0]
459 + B*m_BoundingCorner[c1][1]
460 + C*m_BoundingCorner[c1][2] );
462 // initialise plane value and normalise
463 m_BoundingPlane[j][0] = A/vcl_sqrt(A*A + B*B + C*C);
464 m_BoundingPlane[j][1] = B/vcl_sqrt(A*A + B*B + C*C);
465 m_BoundingPlane[j][2] = C/vcl_sqrt(A*A + B*B + C*C);
466 m_BoundingPlane[j][3] = D/vcl_sqrt(A*A + B*B + C*C);
468 if ( (A*A + B*B + C*C) == 0 )
470 itk::ExceptionObject err(__FILE__, __LINE__);
471 err.SetLocation( ITK_LOCATION );
472 err.SetDescription( "Division by zero (planes) "
473 "- CalcPlanesAndCorners().");
481 /* -----------------------------------------------------------------------
482 CalcRayIntercepts() - Calculate the ray intercepts with the volume.
483 ----------------------------------------------------------------------- */
485 template<class TInputImage, class TCoordRep>
487 RayCastHelper<TInputImage, TCoordRep>
488 ::CalcRayIntercepts()
490 double maxInterDist, interDist;
491 double cornerVect[4][3];
492 int cross[4][3], noInterFlag[6];
493 int nSidesCrossed, crossFlag, c[4];
494 double ax, ay, az, bx, by, bz;
495 double cubeInter[6][3];
499 int NoSides = 6; // =6 to allow truncation: =4 to remove truncated rays
501 // Calculate intercept of ray with planes
503 double interceptx[6], intercepty[6], interceptz[6];
506 for( j=0; j<NoSides; j++)
509 denom = ( m_BoundingPlane[j][0]*m_RayDirectionInMM[0]
510 + m_BoundingPlane[j][1]*m_RayDirectionInMM[1]
511 + m_BoundingPlane[j][2]*m_RayDirectionInMM[2]);
513 if( (long)(denom*100) != 0 )
515 d[j] = -( m_BoundingPlane[j][3]
516 + m_BoundingPlane[j][0]*m_CurrentRayPositionInMM[0]
517 + m_BoundingPlane[j][1]*m_CurrentRayPositionInMM[1]
518 + m_BoundingPlane[j][2]*m_CurrentRayPositionInMM[2] ) / denom;
520 interceptx[j] = m_CurrentRayPositionInMM[0] + d[j]*m_RayDirectionInMM[0];
521 intercepty[j] = m_CurrentRayPositionInMM[1] + d[j]*m_RayDirectionInMM[1];
522 interceptz[j] = m_CurrentRayPositionInMM[2] + d[j]*m_RayDirectionInMM[2];
524 noInterFlag[j] = 1; //OK
528 noInterFlag[j] = 0; //NOT OK
534 for( j=0; j<NoSides; j++ )
537 // Work out which corners to use
541 c[0] = 0; c[1] = 1; c[2] = 3; c[3] = 2;
545 c[0] = 4; c[1] = 5; c[2] = 7; c[3] = 6;
549 c[0] = 1; c[1] = 5; c[2] = 7; c[3] = 3;
553 c[0] = 0; c[1] = 2; c[2] = 6; c[3] = 4;
557 c[0] = 0; c[1] = 1; c[2] = 5; c[3] = 4;
561 c[0] = 2; c[1] = 3; c[2] = 7; c[3] = 6;
564 // Calculate vectors from corner of ct volume to intercept.
567 if( noInterFlag[j]==1 )
569 cornerVect[i][0] = m_BoundingCorner[c[i]][0] - interceptx[j];
570 cornerVect[i][1] = m_BoundingCorner[c[i]][1] - intercepty[j];
571 cornerVect[i][2] = m_BoundingCorner[c[i]][2] - interceptz[j];
573 else if( noInterFlag[j]==0 )
575 cornerVect[i][0] = 0;
576 cornerVect[i][1] = 0;
577 cornerVect[i][2] = 0;
582 // Do cross product with these vectors
593 ax = cornerVect[i][0];
594 ay = cornerVect[i][1];
595 az = cornerVect[i][2];
596 bx = cornerVect[k][0];
597 by = cornerVect[k][1];
598 bz = cornerVect[k][2];
600 // The int and divide by 100 are to avoid rounding errors. If
601 // these are not included then you get values fluctuating around
602 // zero and so in the subsequent check, all the values are not
603 // above or below zero. NB. If you "INT" by too much here though
604 // you can get problems in the corners of your volume when rays
605 // are allowed to go through more than one plane.
606 cross[i][0] = (int)((ay*bz - az*by)/100);
607 cross[i][1] = (int)((az*bx - ax*bz)/100);
608 cross[i][2] = (int)((ax*by - ay*bx)/100);
611 // See if a sign change occured between all these cross products
612 // if not, then the ray went through this plane
632 if( crossFlag==3 && noInterFlag[j]==1 )
634 cubeInter[nSidesCrossed][0] = interceptx[j];
635 cubeInter[nSidesCrossed][1] = intercepty[j];
636 cubeInter[nSidesCrossed][2] = interceptz[j];
640 } // End of loop over all four planes
642 m_RayStartCoordInMM[0] = cubeInter[0][0];
643 m_RayStartCoordInMM[1] = cubeInter[0][1];
644 m_RayStartCoordInMM[2] = cubeInter[0][2];
646 m_RayEndCoordInMM[0] = cubeInter[1][0];
647 m_RayEndCoordInMM[1] = cubeInter[1][1];
648 m_RayEndCoordInMM[2] = cubeInter[1][2];
650 if( nSidesCrossed >= 5 )
652 std::cerr << "WARNING: No. of sides crossed equals: " << nSidesCrossed << std::endl;
655 // If 'nSidesCrossed' is larger than 2, this means that the ray goes through
656 // a corner of the volume and due to rounding errors, the ray is
657 // deemed to go through more than two planes. To obtain the correct
658 // start and end positions we choose the two intercept values which
659 // are furthest from each other.
662 if( nSidesCrossed >= 3 )
665 for( j=0; j<nSidesCrossed-1; j++ )
667 for( k=j+1; k<nSidesCrossed; k++ )
672 interDist += (cubeInter[j][i] - cubeInter[k][i])*
673 (cubeInter[j][i] - cubeInter[k][i]);
675 if( interDist > maxInterDist )
677 maxInterDist = interDist;
679 m_RayStartCoordInMM[0] = cubeInter[j][0];
680 m_RayStartCoordInMM[1] = cubeInter[j][1];
681 m_RayStartCoordInMM[2] = cubeInter[j][2];
683 m_RayEndCoordInMM[0] = cubeInter[k][0];
684 m_RayEndCoordInMM[1] = cubeInter[k][1];
685 m_RayEndCoordInMM[2] = cubeInter[k][2];
692 if (nSidesCrossed == 2 )
703 /* -----------------------------------------------------------------------
704 SetRay() - Set the position and direction of the ray
705 ----------------------------------------------------------------------- */
707 template<class TInputImage, class TCoordRep>
709 RayCastHelper<TInputImage, TCoordRep>
710 ::SetRay(OutputPointType RayPosn, DirectionType RayDirn)
713 // Store the position and direction of the ray
714 typename TInputImage::SpacingType spacing=this->m_Image->GetSpacing();
715 SizeType dim=this->m_Image->GetLargestPossibleRegion().GetSize();
717 // we need to translate the _center_ of the volume to the origin
718 m_NumberOfVoxelsInX = dim[0];
719 m_NumberOfVoxelsInY = dim[1];
720 m_NumberOfVoxelsInZ = dim[2];
722 m_VoxelDimensionInX = spacing[0];
723 m_VoxelDimensionInY = spacing[1];
724 m_VoxelDimensionInZ = spacing[2];
726 // (Aviv) Incorporating a fix by Jian Wu
727 // http://public.kitware.com/pipermail/insight-users/2006-March/017265.html
728 m_CurrentRayPositionInMM[0] = RayPosn[0];
729 m_CurrentRayPositionInMM[1] = RayPosn[1];
730 m_CurrentRayPositionInMM[2] = RayPosn[2];
732 m_RayDirectionInMM[0] = RayDirn[0];
733 m_RayDirectionInMM[1] = RayDirn[1];
734 m_RayDirectionInMM[2] = RayDirn[2];
736 // Compute the ray path for this coordinate in mm
738 m_ValidRay = this->CalcRayIntercepts();
746 // Convert the start and end coordinates of the ray to voxels
748 this->EndPointsInVoxels();
750 /* Calculate the ray direction vector in voxels and hence the voxel
751 increment required to traverse the ray, and the number of
752 interpolation points on the ray.
754 This routine also shifts the coordinate frame by half a voxel for
755 two of the directional components (i.e. those lying within the
756 planes of voxels being traversed). */
758 this->CalcDirnVector();
761 /* Reduce the length of the ray until both start and end
762 coordinates lie inside the volume. */
764 m_ValidRay = this->AdjustRayLength();
766 // Reset the iterator to the start of the ray.
774 /* -----------------------------------------------------------------------
775 EndPointsInVoxels() - Convert the endpoints to voxels
776 ----------------------------------------------------------------------- */
778 template<class TInputImage, class TCoordRep>
780 RayCastHelper<TInputImage, TCoordRep>
781 ::EndPointsInVoxels(void)
783 m_RayVoxelStartPosition[0] = m_RayStartCoordInMM[0]/m_VoxelDimensionInX;
784 m_RayVoxelStartPosition[1] = m_RayStartCoordInMM[1]/m_VoxelDimensionInY;
785 m_RayVoxelStartPosition[2] = m_RayStartCoordInMM[2]/m_VoxelDimensionInZ;
787 m_RayVoxelEndPosition[0] = m_RayEndCoordInMM[0]/m_VoxelDimensionInX;
788 m_RayVoxelEndPosition[1] = m_RayEndCoordInMM[1]/m_VoxelDimensionInY;
789 m_RayVoxelEndPosition[2] = m_RayEndCoordInMM[2]/m_VoxelDimensionInZ;
794 /* -----------------------------------------------------------------------
795 CalcDirnVector() - Calculate the incremental direction vector in voxels.
796 ----------------------------------------------------------------------- */
798 template<class TInputImage, class TCoordRep>
800 RayCastHelper<TInputImage, TCoordRep>
801 ::CalcDirnVector(void)
803 double xNum, yNum, zNum;
805 // Calculate the number of voxels in each direction
807 xNum = vcl_fabs(m_RayVoxelStartPosition[0] - m_RayVoxelEndPosition[0]);
808 yNum = vcl_fabs(m_RayVoxelStartPosition[1] - m_RayVoxelEndPosition[1]);
809 zNum = vcl_fabs(m_RayVoxelStartPosition[2] - m_RayVoxelEndPosition[2]);
811 // The direction iterated in is that with the greatest number of voxels
812 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
814 // Iterate in X direction
816 if( (xNum >= yNum) && (xNum >= zNum) )
818 if( m_RayVoxelStartPosition[0] < m_RayVoxelEndPosition[0] )
820 m_VoxelIncrement[0] = 1;
823 = (m_RayVoxelStartPosition[1]
824 - m_RayVoxelEndPosition[1])/(m_RayVoxelStartPosition[0]
825 - m_RayVoxelEndPosition[0]);
828 = (m_RayVoxelStartPosition[2]
829 - m_RayVoxelEndPosition[2])/(m_RayVoxelStartPosition[0]
830 - m_RayVoxelEndPosition[0]);
834 m_VoxelIncrement[0] = -1;
837 = -(m_RayVoxelStartPosition[1]
838 - m_RayVoxelEndPosition[1])/(m_RayVoxelStartPosition[0]
839 - m_RayVoxelEndPosition[0]);
842 = -(m_RayVoxelStartPosition[2]
843 - m_RayVoxelEndPosition[2])/(m_RayVoxelStartPosition[0]
844 - m_RayVoxelEndPosition[0]);
847 // This section is to alter the start position in order to
848 // place the center of the voxels in there correct positions,
849 // rather than placing them at the corner of voxels which is
850 // what happens if this is not carried out. The reason why
851 // x has no -0.5 is because this is the direction we are going
852 // to iterate in and therefore we wish to go from center to
853 // center rather than finding the surrounding voxels.
855 m_RayVoxelStartPosition[1] += ( (int)m_RayVoxelStartPosition[0]
856 - m_RayVoxelStartPosition[0])*m_VoxelIncrement[1]*m_VoxelIncrement[0]
857 + 0.5*m_VoxelIncrement[1] - 0.5;
859 m_RayVoxelStartPosition[2] += ( (int)m_RayVoxelStartPosition[0]
860 - m_RayVoxelStartPosition[0])*m_VoxelIncrement[2]*m_VoxelIncrement[0]
861 + 0.5*m_VoxelIncrement[2] - 0.5;
863 m_RayVoxelStartPosition[0] = (int)m_RayVoxelStartPosition[0] + 0.5*m_VoxelIncrement[0];
865 m_TotalRayVoxelPlanes = (int)xNum;
867 m_TraversalDirection = TRANSVERSE_IN_X;
870 // Iterate in Y direction
872 else if( (yNum >= xNum) && (yNum >= zNum) )
875 if( m_RayVoxelStartPosition[1] < m_RayVoxelEndPosition[1] )
877 m_VoxelIncrement[1] = 1;
880 = (m_RayVoxelStartPosition[0]
881 - m_RayVoxelEndPosition[0])/(m_RayVoxelStartPosition[1]
882 - m_RayVoxelEndPosition[1]);
885 = (m_RayVoxelStartPosition[2]
886 - m_RayVoxelEndPosition[2])/(m_RayVoxelStartPosition[1]
887 - m_RayVoxelEndPosition[1]);
891 m_VoxelIncrement[1] = -1;
894 = -(m_RayVoxelStartPosition[0]
895 - m_RayVoxelEndPosition[0])/(m_RayVoxelStartPosition[1]
896 - m_RayVoxelEndPosition[1]);
899 = -(m_RayVoxelStartPosition[2]
900 - m_RayVoxelEndPosition[2])/(m_RayVoxelStartPosition[1]
901 - m_RayVoxelEndPosition[1]);
904 m_RayVoxelStartPosition[0] += ( (int)m_RayVoxelStartPosition[1]
905 - m_RayVoxelStartPosition[1])*m_VoxelIncrement[0]*m_VoxelIncrement[1]
906 + 0.5*m_VoxelIncrement[0] - 0.5;
908 m_RayVoxelStartPosition[2] += ( (int)m_RayVoxelStartPosition[1]
909 - m_RayVoxelStartPosition[1])*m_VoxelIncrement[2]*m_VoxelIncrement[1]
910 + 0.5*m_VoxelIncrement[2] - 0.5;
912 m_RayVoxelStartPosition[1] = (int)m_RayVoxelStartPosition[1] + 0.5*m_VoxelIncrement[1];
914 m_TotalRayVoxelPlanes = (int)yNum;
916 m_TraversalDirection = TRANSVERSE_IN_Y;
919 // Iterate in Z direction
924 if( m_RayVoxelStartPosition[2] < m_RayVoxelEndPosition[2] )
926 m_VoxelIncrement[2] = 1;
929 = (m_RayVoxelStartPosition[0]
930 - m_RayVoxelEndPosition[0])/(m_RayVoxelStartPosition[2]
931 - m_RayVoxelEndPosition[2]);
934 = (m_RayVoxelStartPosition[1]
935 - m_RayVoxelEndPosition[1])/(m_RayVoxelStartPosition[2]
936 - m_RayVoxelEndPosition[2]);
940 m_VoxelIncrement[2] = -1;
943 = -(m_RayVoxelStartPosition[0]
944 - m_RayVoxelEndPosition[0])/(m_RayVoxelStartPosition[2]
945 - m_RayVoxelEndPosition[2]);
948 = -(m_RayVoxelStartPosition[1]
949 - m_RayVoxelEndPosition[1])/(m_RayVoxelStartPosition[2]
950 - m_RayVoxelEndPosition[2]);
953 m_RayVoxelStartPosition[0] += ( (int)m_RayVoxelStartPosition[2]
954 - m_RayVoxelStartPosition[2])*m_VoxelIncrement[0]*m_VoxelIncrement[2]
955 + 0.5*m_VoxelIncrement[0] - 0.5;
957 m_RayVoxelStartPosition[1] += ( (int)m_RayVoxelStartPosition[2]
958 - m_RayVoxelStartPosition[2])*m_VoxelIncrement[1]*m_VoxelIncrement[2]
959 + 0.5*m_VoxelIncrement[1] - 0.5;
961 m_RayVoxelStartPosition[2] = (int)m_RayVoxelStartPosition[2] + 0.5*m_VoxelIncrement[2];
963 m_TotalRayVoxelPlanes = (int)zNum;
965 m_TraversalDirection = TRANSVERSE_IN_Z;
970 /* -----------------------------------------------------------------------
971 AdjustRayLength() - Ensure that the ray lies within the volume
972 ----------------------------------------------------------------------- */
974 template<class TInputImage, class TCoordRep>
976 RayCastHelper<TInputImage, TCoordRep>
977 ::AdjustRayLength(void)
984 if (m_TraversalDirection == TRANSVERSE_IN_X)
990 else if (m_TraversalDirection == TRANSVERSE_IN_Y)
996 else if (m_TraversalDirection == TRANSVERSE_IN_Z)
1004 itk::ExceptionObject err(__FILE__, __LINE__);
1005 err.SetLocation( ITK_LOCATION );
1006 err.SetDescription( "The ray traversal direction is unset "
1007 "- AdjustRayLength().");
1019 Istart[0] = (int) vcl_floor(m_RayVoxelStartPosition[0]);
1020 Istart[1] = (int) vcl_floor(m_RayVoxelStartPosition[1]);
1021 Istart[2] = (int) vcl_floor(m_RayVoxelStartPosition[2]);
1023 if( (Istart[0] >= 0) && (Istart[0] + Idirn[0] < m_NumberOfVoxelsInX) &&
1024 (Istart[1] >= 0) && (Istart[1] + Idirn[1] < m_NumberOfVoxelsInY) &&
1025 (Istart[2] >= 0) && (Istart[2] + Idirn[2] < m_NumberOfVoxelsInZ) )
1032 m_RayVoxelStartPosition[0] += m_VoxelIncrement[0];
1033 m_RayVoxelStartPosition[1] += m_VoxelIncrement[1];
1034 m_RayVoxelStartPosition[2] += m_VoxelIncrement[2];
1036 m_TotalRayVoxelPlanes--;
1039 Istart[0] = (int) vcl_floor(m_RayVoxelStartPosition[0]
1040 + m_TotalRayVoxelPlanes*m_VoxelIncrement[0]);
1042 Istart[1] = (int) vcl_floor(m_RayVoxelStartPosition[1]
1043 + m_TotalRayVoxelPlanes*m_VoxelIncrement[1]);
1045 Istart[2] = (int) vcl_floor(m_RayVoxelStartPosition[2]
1046 + m_TotalRayVoxelPlanes*m_VoxelIncrement[2]);
1048 if( (Istart[0] >= 0) && (Istart[0] + Idirn[0] < m_NumberOfVoxelsInX) &&
1049 (Istart[1] >= 0) && (Istart[1] + Idirn[1] < m_NumberOfVoxelsInY) &&
1050 (Istart[2] >= 0) && (Istart[2] + Idirn[2] < m_NumberOfVoxelsInZ) )
1057 m_TotalRayVoxelPlanes--;
1060 } while ( (! (startOK && endOK)) && (m_TotalRayVoxelPlanes > 1) );
1063 return (startOK && endOK);
1067 /* -----------------------------------------------------------------------
1068 Reset() - Reset the iterator to the start of the ray.
1069 ----------------------------------------------------------------------- */
1071 template<class TInputImage, class TCoordRep>
1073 RayCastHelper<TInputImage, TCoordRep>
1078 m_NumVoxelPlanesTraversed = -1;
1080 // If this is a valid ray...
1086 m_Position3Dvox[i] = m_RayVoxelStartPosition[i];
1088 this->InitialiseVoxelPointers();
1091 // otherwise set parameters to zero
1097 m_RayVoxelStartPosition[i] = 0.;
1101 m_RayVoxelEndPosition[i] = 0.;
1105 m_VoxelIncrement[i] = 0.;
1107 m_TraversalDirection = UNDEFINED_DIRECTION;
1109 m_TotalRayVoxelPlanes = 0;
1113 m_RayIntersectionVoxels[i] = 0;
1117 m_RayIntersectionVoxelIndex[i] = 0;
1123 /* -----------------------------------------------------------------------
1124 InitialiseVoxelPointers() - Obtain pointers to the first four voxels
1125 ----------------------------------------------------------------------- */
1127 template<class TInputImage, class TCoordRep>
1129 RayCastHelper<TInputImage, TCoordRep>
1130 ::InitialiseVoxelPointers(void)
1136 Ix = (int)(m_RayVoxelStartPosition[0]);
1137 Iy = (int)(m_RayVoxelStartPosition[1]);
1138 Iz = (int)(m_RayVoxelStartPosition[2]);
1140 m_RayIntersectionVoxelIndex[0] = Ix;
1141 m_RayIntersectionVoxelIndex[1] = Iy;
1142 m_RayIntersectionVoxelIndex[2] = Iz;
1144 switch( m_TraversalDirection )
1146 case TRANSVERSE_IN_X:
1149 if( (Ix >= 0) && (Ix < m_NumberOfVoxelsInX) &&
1150 (Iy >= 0) && (Iy + 1 < m_NumberOfVoxelsInY) &&
1151 (Iz >= 0) && (Iz + 1 < m_NumberOfVoxelsInZ))
1154 index[0]=Ix; index[1]=Iy; index[2]=Iz;
1155 m_RayIntersectionVoxels[0]
1156 = this->m_Image->GetBufferPointer() + this->m_Image->ComputeOffset(index);
1158 index[0]=Ix; index[1]=Iy+1; index[2]=Iz;
1159 m_RayIntersectionVoxels[1]
1160 = ( this->m_Image->GetBufferPointer() + this->m_Image->ComputeOffset(index) );
1162 index[0]=Ix; index[1]=Iy; index[2]=Iz+1;
1163 m_RayIntersectionVoxels[2]
1164 = ( this->m_Image->GetBufferPointer() + this->m_Image->ComputeOffset(index) );
1166 index[0]=Ix; index[1]=Iy+1; index[2]=Iz+1;
1167 m_RayIntersectionVoxels[3]
1168 = ( this->m_Image->GetBufferPointer() + this->m_Image->ComputeOffset(index) );
1172 m_RayIntersectionVoxels[0] =
1173 m_RayIntersectionVoxels[1] =
1174 m_RayIntersectionVoxels[2] =
1175 m_RayIntersectionVoxels[3] = NULL;
1180 case TRANSVERSE_IN_Y:
1183 if( (Ix >= 0) && (Ix + 1 < m_NumberOfVoxelsInX) &&
1184 (Iy >= 0) && (Iy < m_NumberOfVoxelsInY) &&
1185 (Iz >= 0) && (Iz + 1 < m_NumberOfVoxelsInZ))
1188 index[0]=Ix; index[1]=Iy; index[2]=Iz;
1189 m_RayIntersectionVoxels[0] = ( this->m_Image->GetBufferPointer()
1190 + this->m_Image->ComputeOffset(index) );
1192 index[0]=Ix+1; index[1]=Iy; index[2]=Iz;
1193 m_RayIntersectionVoxels[1] = ( this->m_Image->GetBufferPointer()
1194 + this->m_Image->ComputeOffset(index) );
1196 index[0]=Ix; index[1]=Iy; index[2]=Iz+1;
1197 m_RayIntersectionVoxels[2] = ( this->m_Image->GetBufferPointer()
1198 + this->m_Image->ComputeOffset(index) );
1200 index[0]=Ix+1; index[1]=Iy; index[2]=Iz+1;
1201 m_RayIntersectionVoxels[3] = ( this->m_Image->GetBufferPointer()
1202 + this->m_Image->ComputeOffset(index) );
1206 m_RayIntersectionVoxels[0]
1207 = m_RayIntersectionVoxels[1]
1208 = m_RayIntersectionVoxels[2]
1209 = m_RayIntersectionVoxels[3] = NULL;
1214 case TRANSVERSE_IN_Z:
1217 if( (Ix >= 0) && (Ix + 1 < m_NumberOfVoxelsInX) &&
1218 (Iy >= 0) && (Iy + 1 < m_NumberOfVoxelsInY) &&
1219 (Iz >= 0) && (Iz < m_NumberOfVoxelsInZ))
1222 index[0]=Ix; index[1]=Iy; index[2]=Iz;
1223 m_RayIntersectionVoxels[0] = ( this->m_Image->GetBufferPointer()
1224 + this->m_Image->ComputeOffset(index) );
1226 index[0]=Ix+1; index[1]=Iy; index[2]=Iz;
1227 m_RayIntersectionVoxels[1] = ( this->m_Image->GetBufferPointer()
1228 + this->m_Image->ComputeOffset(index) );
1230 index[0]=Ix; index[1]=Iy+1; index[2]=Iz;
1231 m_RayIntersectionVoxels[2] = ( this->m_Image->GetBufferPointer()
1232 + this->m_Image->ComputeOffset(index) );
1234 index[0]=Ix+1; index[1]=Iy+1; index[2]=Iz;
1235 m_RayIntersectionVoxels[3] = ( this->m_Image->GetBufferPointer()
1236 + this->m_Image->ComputeOffset(index) );
1241 m_RayIntersectionVoxels[0]
1242 = m_RayIntersectionVoxels[1]
1243 = m_RayIntersectionVoxels[2]
1244 = m_RayIntersectionVoxels[3] = NULL;
1251 itk::ExceptionObject err(__FILE__, __LINE__);
1252 err.SetLocation( ITK_LOCATION );
1253 err.SetDescription( "The ray traversal direction is unset "
1254 "- InitialiseVoxelPointers().");
1261 /* -----------------------------------------------------------------------
1262 IncrementVoxelPointers() - Increment the voxel pointers
1263 ----------------------------------------------------------------------- */
1265 template<class TInputImage, class TCoordRep>
1267 RayCastHelper<TInputImage, TCoordRep>
1268 ::IncrementVoxelPointers(void)
1270 double xBefore = m_Position3Dvox[0];
1271 double yBefore = m_Position3Dvox[1];
1272 double zBefore = m_Position3Dvox[2];
1274 m_Position3Dvox[0] += m_VoxelIncrement[0];
1275 m_Position3Dvox[1] += m_VoxelIncrement[1];
1276 m_Position3Dvox[2] += m_VoxelIncrement[2];
1278 int dx = ((int) m_Position3Dvox[0]) - ((int) xBefore);
1279 int dy = ((int) m_Position3Dvox[1]) - ((int) yBefore);
1280 int dz = ((int) m_Position3Dvox[2]) - ((int) zBefore);
1282 m_RayIntersectionVoxelIndex[0] += dx;
1283 m_RayIntersectionVoxelIndex[1] += dy;
1284 m_RayIntersectionVoxelIndex[2] += dz;
1286 int totalRayVoxelPlanes
1287 = dx + dy*m_NumberOfVoxelsInX + dz*m_NumberOfVoxelsInX*m_NumberOfVoxelsInY;
1289 m_RayIntersectionVoxels[0] += totalRayVoxelPlanes;
1290 m_RayIntersectionVoxels[1] += totalRayVoxelPlanes;
1291 m_RayIntersectionVoxels[2] += totalRayVoxelPlanes;
1292 m_RayIntersectionVoxels[3] += totalRayVoxelPlanes;
1296 /* -----------------------------------------------------------------------
1297 GetCurrentIntensity() - Get the intensity of the current ray point.
1298 ----------------------------------------------------------------------- */
1300 template<class TInputImage, class TCoordRep>
1302 RayCastHelper<TInputImage, TCoordRep>
1303 ::GetCurrentIntensity(void) const
1312 a = (double) (*m_RayIntersectionVoxels[0]);
1313 b = (double) (*m_RayIntersectionVoxels[1] - a);
1314 c = (double) (*m_RayIntersectionVoxels[2] - a);
1315 d = (double) (*m_RayIntersectionVoxels[3] - a - b - c);
1317 switch( m_TraversalDirection )
1319 case TRANSVERSE_IN_X:
1321 y = m_Position3Dvox[1] - vcl_floor(m_Position3Dvox[1]);
1322 z = m_Position3Dvox[2] - vcl_floor(m_Position3Dvox[2]);
1325 case TRANSVERSE_IN_Y:
1327 y = m_Position3Dvox[0] - vcl_floor(m_Position3Dvox[0]);
1328 z = m_Position3Dvox[2] - vcl_floor(m_Position3Dvox[2]);
1331 case TRANSVERSE_IN_Z:
1333 y = m_Position3Dvox[0] - vcl_floor(m_Position3Dvox[0]);
1334 z = m_Position3Dvox[1] - vcl_floor(m_Position3Dvox[1]);
1339 itk::ExceptionObject err(__FILE__, __LINE__);
1340 err.SetLocation( ITK_LOCATION );
1341 err.SetDescription( "The ray traversal direction is unset "
1342 "- GetCurrentIntensity().");
1348 return a + b*y + c*z + d*y*z;
1351 /* -----------------------------------------------------------------------
1352 IncrementIntensities() - Increment the intensities of the current ray point
1353 ----------------------------------------------------------------------- */
1355 template<class TInputImage, class TCoordRep>
1357 RayCastHelper<TInputImage, TCoordRep>
1358 ::IncrementIntensities(double increment)
1360 short inc = (short) vcl_floor(increment + 0.5);
1366 *m_RayIntersectionVoxels[0] += inc;
1367 *m_RayIntersectionVoxels[1] += inc;
1368 *m_RayIntersectionVoxels[2] += inc;
1369 *m_RayIntersectionVoxels[3] += inc;
1375 /* -----------------------------------------------------------------------
1376 IntegrateAboveThreshold() - Integrate intensities above a threshold.
1377 ----------------------------------------------------------------------- */
1379 template<class TInputImage, class TCoordRep>
1381 RayCastHelper<TInputImage, TCoordRep>
1382 ::IntegrateAboveThreshold(double &integral, double threshold)
1385 // double posn3D_x, posn3D_y, posn3D_z;
1389 // Check if this is a valid ray
1395 /* Step along the ray as quickly as possible
1396 integrating the interpolated intensities. */
1398 for (m_NumVoxelPlanesTraversed=0;
1399 m_NumVoxelPlanesTraversed<m_TotalRayVoxelPlanes;
1400 m_NumVoxelPlanesTraversed++)
1402 intensity = this->GetCurrentIntensity();
1404 if (intensity > threshold)
1406 integral += intensity - threshold;
1408 this->IncrementVoxelPointers();
1411 /* The ray passes through the volume one plane of voxels at a time,
1412 however, if its moving diagonally the ray points will be further
1413 apart so account for this by scaling by the distance moved. */
1415 integral *= this->GetRayPointSpacing();
1420 /* -----------------------------------------------------------------------
1421 ZeroState() - Set the default (zero) state of the object
1422 ----------------------------------------------------------------------- */
1424 template<class TInputImage, class TCoordRep>
1426 RayCastHelper<TInputImage, TCoordRep>
1433 m_NumberOfVoxelsInX = 0;
1434 m_NumberOfVoxelsInY = 0;
1435 m_NumberOfVoxelsInZ = 0;
1437 m_VoxelDimensionInX = 0;
1438 m_VoxelDimensionInY = 0;
1439 m_VoxelDimensionInZ = 0;
1443 m_CurrentRayPositionInMM[i] = 0.;
1447 m_RayDirectionInMM[i] = 0.;
1451 m_RayVoxelStartPosition[i] = 0.;
1455 m_RayVoxelEndPosition[i] = 0.;
1459 m_VoxelIncrement[i] = 0.;
1461 m_TraversalDirection = UNDEFINED_DIRECTION;
1463 m_TotalRayVoxelPlanes = 0;
1464 m_NumVoxelPlanesTraversed = -1;
1468 m_RayIntersectionVoxels[i] = 0;
1472 m_RayIntersectionVoxelIndex[i] = 0;
1475 }; // end of anonymous namespace
1481 /**************************************************************************
1484 * Rest of this code is the actual RayCastInterpolateImageFunctionWithOrigin
1488 **************************************************************************/
1490 /* -----------------------------------------------------------------------
1492 ----------------------------------------------------------------------- */
1494 template<class TInputImage, class TCoordRep>
1495 RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
1496 ::RayCastInterpolateImageFunctionWithOrigin()
1500 m_FocalPoint[0] = 0.;
1501 m_FocalPoint[1] = 0.;
1502 m_FocalPoint[2] = 0.;
1506 /* -----------------------------------------------------------------------
1508 ----------------------------------------------------------------------- */
1510 template<class TInputImage, class TCoordRep>
1512 RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
1513 ::PrintSelf(std::ostream& os, Indent indent) const
1515 this->Superclass::PrintSelf(os,indent);
1517 os << indent << "Threshold: " << m_Threshold << std::endl;
1518 os << indent << "FocalPoint: " << m_FocalPoint << std::endl;
1519 os << indent << "Transform: " << m_Transform.GetPointer() << std::endl;
1520 os << indent << "Interpolator: " << m_Interpolator.GetPointer() << std::endl;
1524 /* -----------------------------------------------------------------------
1525 Evaluate at image index position
1526 ----------------------------------------------------------------------- */
1528 template<class TInputImage, class TCoordRep>
1529 typename RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
1531 RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
1532 ::Evaluate( const PointType& point ) const
1534 double integral = 0;
1536 OutputPointType transformedFocalPoint
1537 = m_Transform->TransformPoint( m_FocalPoint );
1539 DirectionType direction = transformedFocalPoint - point;
1541 RayCastHelper<TInputImage, TCoordRep> ray;
1542 ray.SetImage( this->m_Image );
1546 // (Aviv) Added support for images with non-zero origin
1547 // ray.SetRay(point, direction);
1548 ray.SetRay(point - this->m_Image->GetOrigin().GetVectorFromOrigin(), direction);
1549 ray.IntegrateAboveThreshold(integral, m_Threshold);
1551 return ( static_cast<OutputType>( integral ));
1554 template<class TInputImage, class TCoordRep>
1555 typename RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
1557 RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
1558 ::EvaluateAtContinuousIndex( const ContinuousIndexType& index ) const
1560 OutputPointType point;
1561 this->m_Image->TransformContinuousIndexToPhysicalPoint(index, point);
1563 return this->Evaluate( point );