1 /*=========================================================================
3 Program: Insight Segmentation & Registration Toolkit
6 Copyright (c) Insight Software Consortium. All rights reserved.
7 See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
9 This software is distributed WITHOUT ANY WARRANTY; without even
10 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
11 PURPOSE. See the above copyright notices for more information.
13 =========================================================================*/
15 #ifndef __itkRayCastInterpolateImageFunctionWithOrigin_txx
16 #define __itkRayCastInterpolateImageFunctionWithOrigin_txx
18 #include "itkRayCastInterpolateImageFunctionWithOrigin.h"
20 #include "vnl/vnl_math.h"
23 // Put the helper class in an anonymous namespace so that it is not
24 // exposed to the user
27 /** \class Helper class to maintain state when casting a ray.
28 * This helper class keeps the RayCastInterpolateImageFunctionWithOrigin thread safe.
30 template <class TInputImage, class TCoordRep = float>
34 /** Constants for the image dimensions */
35 itkStaticConstMacro(InputImageDimension, unsigned int,
36 TInputImage::ImageDimension);
39 * Type of the Transform Base class
40 * The fixed image should be a 3D image
42 typedef itk::Transform<TCoordRep,3,3> TransformType;
44 typedef typename TransformType::Pointer TransformPointer;
45 typedef typename TransformType::InputPointType InputPointType;
46 typedef typename TransformType::OutputPointType OutputPointType;
47 typedef typename TransformType::ParametersType TransformParametersType;
48 typedef typename TransformType::JacobianType TransformJacobianType;
50 typedef typename TInputImage::SizeType SizeType;
51 typedef itk::Vector<TCoordRep, 3> DirectionType;
52 typedef itk::Point<TCoordRep, 3> PointType;
54 typedef TInputImage InputImageType;
55 typedef typename InputImageType::PixelType PixelType;
56 typedef typename InputImageType::IndexType IndexType;
61 void SetImage(const InputImageType *input)
67 * Initialise the ray using the position and direction of a line.
69 * \param RayPosn The position of the ray in 3D (mm).
70 * \param RayDirn The direction of the ray in 3D (mm).
72 * \return True if this is a valid ray.
74 bool SetRay(OutputPointType RayPosn, DirectionType RayDirn);
78 * Integrate the interpolated intensities along the ray and
81 * This routine can be called after instantiating the ray and
82 * calling SetProjectionCoord2D() or Reset(). It may then be called
83 * as many times thereafter for different 2D projection
86 * \param integral The integrated intensities along the ray.
88 * \return True if a valid ray was specified.
90 bool Integrate(double &integral)
92 return IntegrateAboveThreshold(integral, 0);
97 * Integrate the interpolated intensities above a given threshold,
98 * along the ray and return the result.
100 * This routine can be called after instantiating the ray and
101 * calling SetProjectionCoord2D() or Reset(). It may then be called
102 * as many times thereafter for different 2D projection
105 * \param integral The integrated intensities along the ray.
106 * \param threshold The integration threshold [default value: 0]
108 * \return True if a valid ray was specified.
110 bool IntegrateAboveThreshold(double &integral, double threshold);
113 * Increment each of the intensities of the 4 planar voxels
114 * surrounding the current ray point.
116 * \parameter increment Intensity increment for each of the current 4 voxels
118 void IncrementIntensities(double increment=1);
120 /// Reset the iterator to the start of the ray.
123 /// Return the interpolated intensity of the current ray point.
124 double GetCurrentIntensity(void) const;
126 /// Return the ray point spacing in mm
127 double GetRayPointSpacing(void) const {
128 typename InputImageType::SpacingType spacing=this->m_Image->GetSpacing();
131 return vcl_sqrt(m_VoxelIncrement[0]*spacing[0]*m_VoxelIncrement[0]*spacing[0]
132 + m_VoxelIncrement[1]*spacing[1]*m_VoxelIncrement[1]*spacing[1]
133 + m_VoxelIncrement[2]*spacing[2]*m_VoxelIncrement[2]*spacing[2] );
138 /// Set the initial zero state of the object
141 /// Initialise the object
142 void Initialise(void);
145 /// Calculate the endpoint coordinats of the ray in voxels.
146 void EndPointsInVoxels(void);
149 * Calculate the incremental direction vector in voxels, 'dVoxel',
150 * required to traverse the ray.
152 void CalcDirnVector(void);
155 * Reduce the length of the ray until both start and end
156 * coordinates lie inside the volume.
158 * \return True if a valid ray has been, false otherwise.
160 bool AdjustRayLength(void);
163 * Obtain pointers to the four voxels surrounding the point where the ray
166 void InitialiseVoxelPointers(void);
168 /// Increment the voxel pointers surrounding the current point on the ray.
169 void IncrementVoxelPointers(void);
171 /// Record volume dimensions and resolution
172 void RecordVolumeDimensions(void);
174 /// Define the corners of the volume
175 void DefineCorners(void);
178 * Calculate the planes which define the volume.
180 * Member function to calculate the equations of the planes of 4 of
181 * the sides of the volume, calculate the positions of the 8 corners
182 * of the volume in mm in World, also calculate the values of the
183 * slopes of the lines which go to make up the volume( defined as
184 * lines in cube x,y,z dirn and then each of these lines has a slope
185 * in the world x,y,z dirn [3]) and finally also to return the length
186 * of the sides of the lines in mm.
188 void CalcPlanesAndCorners(void);
191 * Calculate the ray intercepts with the volume.
193 * See where the ray cuts the volume, check that truncation does not occur,
194 * if not, then start ray where it first intercepts the volume and set
195 * x_max to be where it leaves the volume.
197 * \return True if a valid ray has been specified, false otherwise.
199 bool CalcRayIntercepts(void);
202 * The ray is traversed by stepping in the axial direction
203 * that enables the greatest number of planes in the volume to be
207 UNDEFINED_DIRECTION=0, //!< Undefined
208 TRANSVERSE_IN_X, //!< x
209 TRANSVERSE_IN_Y, //!< y
210 TRANSVERSE_IN_Z, //!< z
212 } TraversalDirection;
214 // Cache the image in the structure. Skip the smart pointer for
215 // efficiency. This inner class will go in/out of scope with every
216 // call to Evaluate()
217 const InputImageType *m_Image;
219 /// Flag indicating whether the current ray is valid
223 * The start position of the ray in voxels.
225 * NB. Two of the components of this coordinate (i.e. those lying within
226 * the planes of voxels being traversed) will be shifted by half a
227 * voxel. This enables indices of the neighbouring voxels within the plane
228 * to be determined by simply casting to 'int' and optionally adding 1.
230 double m_RayVoxelStartPosition[3];
233 * The end coordinate of the ray in voxels.
235 * NB. Two of the components of this coordinate (i.e. those lying within
236 * the planes of voxels being traversed) will be shifted by half a
237 * voxel. This enables indices of the neighbouring voxels within the plane
238 * to be determined by simply casting to 'int' and optionally adding 1.
240 double m_RayVoxelEndPosition[3];
244 * The current coordinate on the ray in voxels.
246 * NB. Two of the components of this coordinate (i.e. those lying within
247 * the planes of voxels being traversed) will be shifted by half a
248 * voxel. This enables indices of the neighbouring voxels within the plane
249 * to be determined by simply casting to 'int' and optionally adding 1.
251 double m_Position3Dvox[3];
253 /** The incremental direction vector of the ray in voxels. */
254 double m_VoxelIncrement[3];
256 /// The direction in which the ray is incremented thorough the volume (x, y or z).
257 TraversalDirection m_TraversalDirection;
259 /// The total number of planes of voxels traversed by the ray.
260 int m_TotalRayVoxelPlanes;
262 /// The current number of planes of voxels traversed by the ray.
263 int m_NumVoxelPlanesTraversed;
265 /// Pointers to the current four voxels surrounding the ray's trajectory.
266 const PixelType *m_RayIntersectionVoxels[4];
269 * The voxel coordinate of the bottom-left voxel of the current
270 * four voxels surrounding the ray's trajectory.
272 int m_RayIntersectionVoxelIndex[3];
274 /// The dimension in voxels of the 3D volume in along the x axis
275 int m_NumberOfVoxelsInX;
276 /// The dimension in voxels of the 3D volume in along the y axis
277 int m_NumberOfVoxelsInY;
278 /// The dimension in voxels of the 3D volume in along the z axis
279 int m_NumberOfVoxelsInZ;
281 /// Voxel dimension in x
282 double m_VoxelDimensionInX;
283 /// Voxel dimension in y
284 double m_VoxelDimensionInY;
285 /// Voxel dimension in z
286 double m_VoxelDimensionInZ;
288 /// The coordinate of the point at which the ray enters the volume in mm.
289 double m_RayStartCoordInMM[3];
290 /// The coordinate of the point at which the ray exits the volume in mm.
291 double m_RayEndCoordInMM[3];
295 Planes which define the boundary of the volume in mm
296 (six planes and four parameters: Ax+By+Cz+D). */
297 double m_BoundingPlane[6][4];
298 /// The eight corners of the volume (x,y,z coordinates for each).
299 double m_BoundingCorner[8][3];
301 /// The position of the ray
302 double m_CurrentRayPositionInMM[3];
304 /// The direction of the ray
305 double m_RayDirectionInMM[3];
309 /* -----------------------------------------------------------------------
310 Initialise() - Initialise the object
311 ----------------------------------------------------------------------- */
313 template<class TInputImage, class TCoordRep>
315 RayCastHelper<TInputImage, TCoordRep>
318 // Save the dimensions of the volume and calculate the bounding box
319 this->RecordVolumeDimensions();
321 // Calculate the planes and corners which define the volume.
322 this->DefineCorners();
323 this->CalcPlanesAndCorners();
327 /* -----------------------------------------------------------------------
328 RecordVolumeDimensions() - Record volume dimensions and resolution
329 ----------------------------------------------------------------------- */
331 template<class TInputImage, class TCoordRep>
333 RayCastHelper<TInputImage, TCoordRep>
334 ::RecordVolumeDimensions(void)
336 typename InputImageType::SpacingType spacing=this->m_Image->GetSpacing();
337 SizeType dim=this->m_Image->GetLargestPossibleRegion().GetSize();
339 m_NumberOfVoxelsInX = dim[0];
340 m_NumberOfVoxelsInY = dim[1];
341 m_NumberOfVoxelsInZ = dim[2];
343 m_VoxelDimensionInX = spacing[0];
344 m_VoxelDimensionInY = spacing[1];
345 m_VoxelDimensionInZ = spacing[2];
349 /* -----------------------------------------------------------------------
350 DefineCorners() - Define the corners of the volume
351 ----------------------------------------------------------------------- */
353 template<class TInputImage, class TCoordRep>
355 RayCastHelper<TInputImage, TCoordRep>
356 ::DefineCorners(void)
358 // Define corner positions as if at the origin
360 m_BoundingCorner[0][0] =
361 m_BoundingCorner[1][0] =
362 m_BoundingCorner[2][0] =
363 m_BoundingCorner[3][0] = 0;
365 m_BoundingCorner[4][0] =
366 m_BoundingCorner[5][0] =
367 m_BoundingCorner[6][0] =
368 m_BoundingCorner[7][0] = m_VoxelDimensionInX*m_NumberOfVoxelsInX;
370 m_BoundingCorner[1][1] =
371 m_BoundingCorner[3][1] =
372 m_BoundingCorner[5][1] =
373 m_BoundingCorner[7][1] = m_VoxelDimensionInY*m_NumberOfVoxelsInY;
375 m_BoundingCorner[0][1] =
376 m_BoundingCorner[2][1] =
377 m_BoundingCorner[4][1] =
378 m_BoundingCorner[6][1] = 0;
380 m_BoundingCorner[0][2] =
381 m_BoundingCorner[1][2] =
382 m_BoundingCorner[4][2] =
383 m_BoundingCorner[5][2] =
384 m_VoxelDimensionInZ*m_NumberOfVoxelsInZ;
386 m_BoundingCorner[2][2] =
387 m_BoundingCorner[3][2] =
388 m_BoundingCorner[6][2] =
389 m_BoundingCorner[7][2] = 0;
393 /* -----------------------------------------------------------------------
394 CalcPlanesAndCorners() - Calculate the planes and corners of the volume.
395 ----------------------------------------------------------------------- */
397 template<class TInputImage, class TCoordRep>
399 RayCastHelper<TInputImage, TCoordRep>
400 ::CalcPlanesAndCorners(void)
405 // find the equations of the planes
407 int c1=0, c2=0, c3=0;
410 { // loop around for planes
412 { // which corners to take
434 double line1x, line1y, line1z;
435 double line2x, line2y, line2z;
437 // lines from one corner to another in x,y,z dirns
438 line1x = m_BoundingCorner[c1][0] - m_BoundingCorner[c2][0];
439 line2x = m_BoundingCorner[c1][0] - m_BoundingCorner[c3][0];
441 line1y = m_BoundingCorner[c1][1] - m_BoundingCorner[c2][1];
442 line2y = m_BoundingCorner[c1][1] - m_BoundingCorner[c3][1];
444 line1z = m_BoundingCorner[c1][2] - m_BoundingCorner[c2][2];
445 line2z = m_BoundingCorner[c1][2] - m_BoundingCorner[c3][2];
449 // take cross product
450 A = line1y*line2z - line2y*line1z;
451 B = line2x*line1z - line1x*line2z;
452 C = line1x*line2y - line2x*line1y;
455 D = -( A*m_BoundingCorner[c1][0]
456 + B*m_BoundingCorner[c1][1]
457 + C*m_BoundingCorner[c1][2] );
459 // initialise plane value and normalise
460 m_BoundingPlane[j][0] = A/vcl_sqrt(A*A + B*B + C*C);
461 m_BoundingPlane[j][1] = B/vcl_sqrt(A*A + B*B + C*C);
462 m_BoundingPlane[j][2] = C/vcl_sqrt(A*A + B*B + C*C);
463 m_BoundingPlane[j][3] = D/vcl_sqrt(A*A + B*B + C*C);
465 if ( (A*A + B*B + C*C) == 0 )
467 itk::ExceptionObject err(__FILE__, __LINE__);
468 err.SetLocation( ITK_LOCATION );
469 err.SetDescription( "Division by zero (planes) "
470 "- CalcPlanesAndCorners().");
478 /* -----------------------------------------------------------------------
479 CalcRayIntercepts() - Calculate the ray intercepts with the volume.
480 ----------------------------------------------------------------------- */
482 template<class TInputImage, class TCoordRep>
484 RayCastHelper<TInputImage, TCoordRep>
485 ::CalcRayIntercepts()
487 double maxInterDist, interDist;
488 double cornerVect[4][3];
489 int cross[4][3], noInterFlag[6];
490 int nSidesCrossed, crossFlag, c[4];
491 double ax, ay, az, bx, by, bz;
492 double cubeInter[6][3];
496 int NoSides = 6; // =6 to allow truncation: =4 to remove truncated rays
498 // Calculate intercept of ray with planes
500 double interceptx[6], intercepty[6], interceptz[6];
503 for( j=0; j<NoSides; j++)
506 denom = ( m_BoundingPlane[j][0]*m_RayDirectionInMM[0]
507 + m_BoundingPlane[j][1]*m_RayDirectionInMM[1]
508 + m_BoundingPlane[j][2]*m_RayDirectionInMM[2]);
510 if( (long)(denom*100) != 0 )
512 d[j] = -( m_BoundingPlane[j][3]
513 + m_BoundingPlane[j][0]*m_CurrentRayPositionInMM[0]
514 + m_BoundingPlane[j][1]*m_CurrentRayPositionInMM[1]
515 + m_BoundingPlane[j][2]*m_CurrentRayPositionInMM[2] ) / denom;
517 interceptx[j] = m_CurrentRayPositionInMM[0] + d[j]*m_RayDirectionInMM[0];
518 intercepty[j] = m_CurrentRayPositionInMM[1] + d[j]*m_RayDirectionInMM[1];
519 interceptz[j] = m_CurrentRayPositionInMM[2] + d[j]*m_RayDirectionInMM[2];
521 noInterFlag[j] = 1; //OK
525 noInterFlag[j] = 0; //NOT OK
531 for( j=0; j<NoSides; j++ )
534 // Work out which corners to use
538 c[0] = 0; c[1] = 1; c[2] = 3; c[3] = 2;
542 c[0] = 4; c[1] = 5; c[2] = 7; c[3] = 6;
546 c[0] = 1; c[1] = 5; c[2] = 7; c[3] = 3;
550 c[0] = 0; c[1] = 2; c[2] = 6; c[3] = 4;
554 c[0] = 0; c[1] = 1; c[2] = 5; c[3] = 4;
558 c[0] = 2; c[1] = 3; c[2] = 7; c[3] = 6;
561 // Calculate vectors from corner of ct volume to intercept.
564 if( noInterFlag[j]==1 )
566 cornerVect[i][0] = m_BoundingCorner[c[i]][0] - interceptx[j];
567 cornerVect[i][1] = m_BoundingCorner[c[i]][1] - intercepty[j];
568 cornerVect[i][2] = m_BoundingCorner[c[i]][2] - interceptz[j];
570 else if( noInterFlag[j]==0 )
572 cornerVect[i][0] = 0;
573 cornerVect[i][1] = 0;
574 cornerVect[i][2] = 0;
579 // Do cross product with these vectors
590 ax = cornerVect[i][0];
591 ay = cornerVect[i][1];
592 az = cornerVect[i][2];
593 bx = cornerVect[k][0];
594 by = cornerVect[k][1];
595 bz = cornerVect[k][2];
597 // The int and divide by 100 are to avoid rounding errors. If
598 // these are not included then you get values fluctuating around
599 // zero and so in the subsequent check, all the values are not
600 // above or below zero. NB. If you "INT" by too much here though
601 // you can get problems in the corners of your volume when rays
602 // are allowed to go through more than one plane.
603 cross[i][0] = (int)((ay*bz - az*by)/100);
604 cross[i][1] = (int)((az*bx - ax*bz)/100);
605 cross[i][2] = (int)((ax*by - ay*bx)/100);
608 // See if a sign change occured between all these cross products
609 // if not, then the ray went through this plane
629 if( crossFlag==3 && noInterFlag[j]==1 )
631 cubeInter[nSidesCrossed][0] = interceptx[j];
632 cubeInter[nSidesCrossed][1] = intercepty[j];
633 cubeInter[nSidesCrossed][2] = interceptz[j];
637 } // End of loop over all four planes
639 m_RayStartCoordInMM[0] = cubeInter[0][0];
640 m_RayStartCoordInMM[1] = cubeInter[0][1];
641 m_RayStartCoordInMM[2] = cubeInter[0][2];
643 m_RayEndCoordInMM[0] = cubeInter[1][0];
644 m_RayEndCoordInMM[1] = cubeInter[1][1];
645 m_RayEndCoordInMM[2] = cubeInter[1][2];
647 if( nSidesCrossed >= 5 )
649 std::cerr << "WARNING: No. of sides crossed equals: " << nSidesCrossed << std::endl;
652 // If 'nSidesCrossed' is larger than 2, this means that the ray goes through
653 // a corner of the volume and due to rounding errors, the ray is
654 // deemed to go through more than two planes. To obtain the correct
655 // start and end positions we choose the two intercept values which
656 // are furthest from each other.
659 if( nSidesCrossed >= 3 )
662 for( j=0; j<nSidesCrossed-1; j++ )
664 for( k=j+1; k<nSidesCrossed; k++ )
669 interDist += (cubeInter[j][i] - cubeInter[k][i])*
670 (cubeInter[j][i] - cubeInter[k][i]);
672 if( interDist > maxInterDist )
674 maxInterDist = interDist;
676 m_RayStartCoordInMM[0] = cubeInter[j][0];
677 m_RayStartCoordInMM[1] = cubeInter[j][1];
678 m_RayStartCoordInMM[2] = cubeInter[j][2];
680 m_RayEndCoordInMM[0] = cubeInter[k][0];
681 m_RayEndCoordInMM[1] = cubeInter[k][1];
682 m_RayEndCoordInMM[2] = cubeInter[k][2];
689 if (nSidesCrossed == 2 )
700 /* -----------------------------------------------------------------------
701 SetRay() - Set the position and direction of the ray
702 ----------------------------------------------------------------------- */
704 template<class TInputImage, class TCoordRep>
706 RayCastHelper<TInputImage, TCoordRep>
707 ::SetRay(OutputPointType RayPosn, DirectionType RayDirn)
710 // Store the position and direction of the ray
711 typename TInputImage::SpacingType spacing=this->m_Image->GetSpacing();
712 SizeType dim=this->m_Image->GetLargestPossibleRegion().GetSize();
714 // we need to translate the _center_ of the volume to the origin
715 m_NumberOfVoxelsInX = dim[0];
716 m_NumberOfVoxelsInY = dim[1];
717 m_NumberOfVoxelsInZ = dim[2];
719 m_VoxelDimensionInX = spacing[0];
720 m_VoxelDimensionInY = spacing[1];
721 m_VoxelDimensionInZ = spacing[2];
723 // (Aviv) Incorporating a fix by Jian Wu
724 // http://public.kitware.com/pipermail/insight-users/2006-March/017265.html
725 m_CurrentRayPositionInMM[0] = RayPosn[0];
726 m_CurrentRayPositionInMM[1] = RayPosn[1];
727 m_CurrentRayPositionInMM[2] = RayPosn[2];
729 m_RayDirectionInMM[0] = RayDirn[0];
730 m_RayDirectionInMM[1] = RayDirn[1];
731 m_RayDirectionInMM[2] = RayDirn[2];
733 // Compute the ray path for this coordinate in mm
735 m_ValidRay = this->CalcRayIntercepts();
743 // Convert the start and end coordinates of the ray to voxels
745 this->EndPointsInVoxels();
747 /* Calculate the ray direction vector in voxels and hence the voxel
748 increment required to traverse the ray, and the number of
749 interpolation points on the ray.
751 This routine also shifts the coordinate frame by half a voxel for
752 two of the directional components (i.e. those lying within the
753 planes of voxels being traversed). */
755 this->CalcDirnVector();
758 /* Reduce the length of the ray until both start and end
759 coordinates lie inside the volume. */
761 m_ValidRay = this->AdjustRayLength();
763 // Reset the iterator to the start of the ray.
771 /* -----------------------------------------------------------------------
772 EndPointsInVoxels() - Convert the endpoints to voxels
773 ----------------------------------------------------------------------- */
775 template<class TInputImage, class TCoordRep>
777 RayCastHelper<TInputImage, TCoordRep>
778 ::EndPointsInVoxels(void)
780 m_RayVoxelStartPosition[0] = m_RayStartCoordInMM[0]/m_VoxelDimensionInX;
781 m_RayVoxelStartPosition[1] = m_RayStartCoordInMM[1]/m_VoxelDimensionInY;
782 m_RayVoxelStartPosition[2] = m_RayStartCoordInMM[2]/m_VoxelDimensionInZ;
784 m_RayVoxelEndPosition[0] = m_RayEndCoordInMM[0]/m_VoxelDimensionInX;
785 m_RayVoxelEndPosition[1] = m_RayEndCoordInMM[1]/m_VoxelDimensionInY;
786 m_RayVoxelEndPosition[2] = m_RayEndCoordInMM[2]/m_VoxelDimensionInZ;
791 /* -----------------------------------------------------------------------
792 CalcDirnVector() - Calculate the incremental direction vector in voxels.
793 ----------------------------------------------------------------------- */
795 template<class TInputImage, class TCoordRep>
797 RayCastHelper<TInputImage, TCoordRep>
798 ::CalcDirnVector(void)
800 double xNum, yNum, zNum;
802 // Calculate the number of voxels in each direction
804 xNum = vcl_fabs(m_RayVoxelStartPosition[0] - m_RayVoxelEndPosition[0]);
805 yNum = vcl_fabs(m_RayVoxelStartPosition[1] - m_RayVoxelEndPosition[1]);
806 zNum = vcl_fabs(m_RayVoxelStartPosition[2] - m_RayVoxelEndPosition[2]);
808 // The direction iterated in is that with the greatest number of voxels
809 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
811 // Iterate in X direction
813 if( (xNum >= yNum) && (xNum >= zNum) )
815 if( m_RayVoxelStartPosition[0] < m_RayVoxelEndPosition[0] )
817 m_VoxelIncrement[0] = 1;
820 = (m_RayVoxelStartPosition[1]
821 - m_RayVoxelEndPosition[1])/(m_RayVoxelStartPosition[0]
822 - m_RayVoxelEndPosition[0]);
825 = (m_RayVoxelStartPosition[2]
826 - m_RayVoxelEndPosition[2])/(m_RayVoxelStartPosition[0]
827 - m_RayVoxelEndPosition[0]);
831 m_VoxelIncrement[0] = -1;
834 = -(m_RayVoxelStartPosition[1]
835 - m_RayVoxelEndPosition[1])/(m_RayVoxelStartPosition[0]
836 - m_RayVoxelEndPosition[0]);
839 = -(m_RayVoxelStartPosition[2]
840 - m_RayVoxelEndPosition[2])/(m_RayVoxelStartPosition[0]
841 - m_RayVoxelEndPosition[0]);
844 // This section is to alter the start position in order to
845 // place the center of the voxels in there correct positions,
846 // rather than placing them at the corner of voxels which is
847 // what happens if this is not carried out. The reason why
848 // x has no -0.5 is because this is the direction we are going
849 // to iterate in and therefore we wish to go from center to
850 // center rather than finding the surrounding voxels.
852 m_RayVoxelStartPosition[1] += ( (int)m_RayVoxelStartPosition[0]
853 - m_RayVoxelStartPosition[0])*m_VoxelIncrement[1]*m_VoxelIncrement[0]
854 + 0.5*m_VoxelIncrement[1] - 0.5;
856 m_RayVoxelStartPosition[2] += ( (int)m_RayVoxelStartPosition[0]
857 - m_RayVoxelStartPosition[0])*m_VoxelIncrement[2]*m_VoxelIncrement[0]
858 + 0.5*m_VoxelIncrement[2] - 0.5;
860 m_RayVoxelStartPosition[0] = (int)m_RayVoxelStartPosition[0] + 0.5*m_VoxelIncrement[0];
862 m_TotalRayVoxelPlanes = (int)xNum;
864 m_TraversalDirection = TRANSVERSE_IN_X;
867 // Iterate in Y direction
869 else if( (yNum >= xNum) && (yNum >= zNum) )
872 if( m_RayVoxelStartPosition[1] < m_RayVoxelEndPosition[1] )
874 m_VoxelIncrement[1] = 1;
877 = (m_RayVoxelStartPosition[0]
878 - m_RayVoxelEndPosition[0])/(m_RayVoxelStartPosition[1]
879 - m_RayVoxelEndPosition[1]);
882 = (m_RayVoxelStartPosition[2]
883 - m_RayVoxelEndPosition[2])/(m_RayVoxelStartPosition[1]
884 - m_RayVoxelEndPosition[1]);
888 m_VoxelIncrement[1] = -1;
891 = -(m_RayVoxelStartPosition[0]
892 - m_RayVoxelEndPosition[0])/(m_RayVoxelStartPosition[1]
893 - m_RayVoxelEndPosition[1]);
896 = -(m_RayVoxelStartPosition[2]
897 - m_RayVoxelEndPosition[2])/(m_RayVoxelStartPosition[1]
898 - m_RayVoxelEndPosition[1]);
901 m_RayVoxelStartPosition[0] += ( (int)m_RayVoxelStartPosition[1]
902 - m_RayVoxelStartPosition[1])*m_VoxelIncrement[0]*m_VoxelIncrement[1]
903 + 0.5*m_VoxelIncrement[0] - 0.5;
905 m_RayVoxelStartPosition[2] += ( (int)m_RayVoxelStartPosition[1]
906 - m_RayVoxelStartPosition[1])*m_VoxelIncrement[2]*m_VoxelIncrement[1]
907 + 0.5*m_VoxelIncrement[2] - 0.5;
909 m_RayVoxelStartPosition[1] = (int)m_RayVoxelStartPosition[1] + 0.5*m_VoxelIncrement[1];
911 m_TotalRayVoxelPlanes = (int)yNum;
913 m_TraversalDirection = TRANSVERSE_IN_Y;
916 // Iterate in Z direction
921 if( m_RayVoxelStartPosition[2] < m_RayVoxelEndPosition[2] )
923 m_VoxelIncrement[2] = 1;
926 = (m_RayVoxelStartPosition[0]
927 - m_RayVoxelEndPosition[0])/(m_RayVoxelStartPosition[2]
928 - m_RayVoxelEndPosition[2]);
931 = (m_RayVoxelStartPosition[1]
932 - m_RayVoxelEndPosition[1])/(m_RayVoxelStartPosition[2]
933 - m_RayVoxelEndPosition[2]);
937 m_VoxelIncrement[2] = -1;
940 = -(m_RayVoxelStartPosition[0]
941 - m_RayVoxelEndPosition[0])/(m_RayVoxelStartPosition[2]
942 - m_RayVoxelEndPosition[2]);
945 = -(m_RayVoxelStartPosition[1]
946 - m_RayVoxelEndPosition[1])/(m_RayVoxelStartPosition[2]
947 - m_RayVoxelEndPosition[2]);
950 m_RayVoxelStartPosition[0] += ( (int)m_RayVoxelStartPosition[2]
951 - m_RayVoxelStartPosition[2])*m_VoxelIncrement[0]*m_VoxelIncrement[2]
952 + 0.5*m_VoxelIncrement[0] - 0.5;
954 m_RayVoxelStartPosition[1] += ( (int)m_RayVoxelStartPosition[2]
955 - m_RayVoxelStartPosition[2])*m_VoxelIncrement[1]*m_VoxelIncrement[2]
956 + 0.5*m_VoxelIncrement[1] - 0.5;
958 m_RayVoxelStartPosition[2] = (int)m_RayVoxelStartPosition[2] + 0.5*m_VoxelIncrement[2];
960 m_TotalRayVoxelPlanes = (int)zNum;
962 m_TraversalDirection = TRANSVERSE_IN_Z;
967 /* -----------------------------------------------------------------------
968 AdjustRayLength() - Ensure that the ray lies within the volume
969 ----------------------------------------------------------------------- */
971 template<class TInputImage, class TCoordRep>
973 RayCastHelper<TInputImage, TCoordRep>
974 ::AdjustRayLength(void)
981 if (m_TraversalDirection == TRANSVERSE_IN_X)
987 else if (m_TraversalDirection == TRANSVERSE_IN_Y)
993 else if (m_TraversalDirection == TRANSVERSE_IN_Z)
1001 itk::ExceptionObject err(__FILE__, __LINE__);
1002 err.SetLocation( ITK_LOCATION );
1003 err.SetDescription( "The ray traversal direction is unset "
1004 "- AdjustRayLength().");
1016 Istart[0] = (int) vcl_floor(m_RayVoxelStartPosition[0]);
1017 Istart[1] = (int) vcl_floor(m_RayVoxelStartPosition[1]);
1018 Istart[2] = (int) vcl_floor(m_RayVoxelStartPosition[2]);
1020 if( (Istart[0] >= 0) && (Istart[0] + Idirn[0] < m_NumberOfVoxelsInX) &&
1021 (Istart[1] >= 0) && (Istart[1] + Idirn[1] < m_NumberOfVoxelsInY) &&
1022 (Istart[2] >= 0) && (Istart[2] + Idirn[2] < m_NumberOfVoxelsInZ) )
1029 m_RayVoxelStartPosition[0] += m_VoxelIncrement[0];
1030 m_RayVoxelStartPosition[1] += m_VoxelIncrement[1];
1031 m_RayVoxelStartPosition[2] += m_VoxelIncrement[2];
1033 m_TotalRayVoxelPlanes--;
1036 Istart[0] = (int) vcl_floor(m_RayVoxelStartPosition[0]
1037 + m_TotalRayVoxelPlanes*m_VoxelIncrement[0]);
1039 Istart[1] = (int) vcl_floor(m_RayVoxelStartPosition[1]
1040 + m_TotalRayVoxelPlanes*m_VoxelIncrement[1]);
1042 Istart[2] = (int) vcl_floor(m_RayVoxelStartPosition[2]
1043 + m_TotalRayVoxelPlanes*m_VoxelIncrement[2]);
1045 if( (Istart[0] >= 0) && (Istart[0] + Idirn[0] < m_NumberOfVoxelsInX) &&
1046 (Istart[1] >= 0) && (Istart[1] + Idirn[1] < m_NumberOfVoxelsInY) &&
1047 (Istart[2] >= 0) && (Istart[2] + Idirn[2] < m_NumberOfVoxelsInZ) )
1054 m_TotalRayVoxelPlanes--;
1057 } while ( (! (startOK && endOK)) && (m_TotalRayVoxelPlanes > 1) );
1060 return (startOK && endOK);
1064 /* -----------------------------------------------------------------------
1065 Reset() - Reset the iterator to the start of the ray.
1066 ----------------------------------------------------------------------- */
1068 template<class TInputImage, class TCoordRep>
1070 RayCastHelper<TInputImage, TCoordRep>
1075 m_NumVoxelPlanesTraversed = -1;
1077 // If this is a valid ray...
1083 m_Position3Dvox[i] = m_RayVoxelStartPosition[i];
1085 this->InitialiseVoxelPointers();
1088 // otherwise set parameters to zero
1094 m_RayVoxelStartPosition[i] = 0.;
1098 m_RayVoxelEndPosition[i] = 0.;
1102 m_VoxelIncrement[i] = 0.;
1104 m_TraversalDirection = UNDEFINED_DIRECTION;
1106 m_TotalRayVoxelPlanes = 0;
1110 m_RayIntersectionVoxels[i] = 0;
1114 m_RayIntersectionVoxelIndex[i] = 0;
1120 /* -----------------------------------------------------------------------
1121 InitialiseVoxelPointers() - Obtain pointers to the first four voxels
1122 ----------------------------------------------------------------------- */
1124 template<class TInputImage, class TCoordRep>
1126 RayCastHelper<TInputImage, TCoordRep>
1127 ::InitialiseVoxelPointers(void)
1133 Ix = (int)(m_RayVoxelStartPosition[0]);
1134 Iy = (int)(m_RayVoxelStartPosition[1]);
1135 Iz = (int)(m_RayVoxelStartPosition[2]);
1137 m_RayIntersectionVoxelIndex[0] = Ix;
1138 m_RayIntersectionVoxelIndex[1] = Iy;
1139 m_RayIntersectionVoxelIndex[2] = Iz;
1141 switch( m_TraversalDirection )
1143 case TRANSVERSE_IN_X:
1146 if( (Ix >= 0) && (Ix < m_NumberOfVoxelsInX) &&
1147 (Iy >= 0) && (Iy + 1 < m_NumberOfVoxelsInY) &&
1148 (Iz >= 0) && (Iz + 1 < m_NumberOfVoxelsInZ))
1151 index[0]=Ix; index[1]=Iy; index[2]=Iz;
1152 m_RayIntersectionVoxels[0]
1153 = this->m_Image->GetBufferPointer() + this->m_Image->ComputeOffset(index);
1155 index[0]=Ix; index[1]=Iy+1; index[2]=Iz;
1156 m_RayIntersectionVoxels[1]
1157 = ( this->m_Image->GetBufferPointer() + this->m_Image->ComputeOffset(index) );
1159 index[0]=Ix; index[1]=Iy; index[2]=Iz+1;
1160 m_RayIntersectionVoxels[2]
1161 = ( this->m_Image->GetBufferPointer() + this->m_Image->ComputeOffset(index) );
1163 index[0]=Ix; index[1]=Iy+1; index[2]=Iz+1;
1164 m_RayIntersectionVoxels[3]
1165 = ( this->m_Image->GetBufferPointer() + this->m_Image->ComputeOffset(index) );
1169 m_RayIntersectionVoxels[0] =
1170 m_RayIntersectionVoxels[1] =
1171 m_RayIntersectionVoxels[2] =
1172 m_RayIntersectionVoxels[3] = NULL;
1177 case TRANSVERSE_IN_Y:
1180 if( (Ix >= 0) && (Ix + 1 < m_NumberOfVoxelsInX) &&
1181 (Iy >= 0) && (Iy < m_NumberOfVoxelsInY) &&
1182 (Iz >= 0) && (Iz + 1 < m_NumberOfVoxelsInZ))
1185 index[0]=Ix; index[1]=Iy; index[2]=Iz;
1186 m_RayIntersectionVoxels[0] = ( this->m_Image->GetBufferPointer()
1187 + this->m_Image->ComputeOffset(index) );
1189 index[0]=Ix+1; index[1]=Iy; index[2]=Iz;
1190 m_RayIntersectionVoxels[1] = ( this->m_Image->GetBufferPointer()
1191 + this->m_Image->ComputeOffset(index) );
1193 index[0]=Ix; index[1]=Iy; index[2]=Iz+1;
1194 m_RayIntersectionVoxels[2] = ( this->m_Image->GetBufferPointer()
1195 + this->m_Image->ComputeOffset(index) );
1197 index[0]=Ix+1; index[1]=Iy; index[2]=Iz+1;
1198 m_RayIntersectionVoxels[3] = ( this->m_Image->GetBufferPointer()
1199 + this->m_Image->ComputeOffset(index) );
1203 m_RayIntersectionVoxels[0]
1204 = m_RayIntersectionVoxels[1]
1205 = m_RayIntersectionVoxels[2]
1206 = m_RayIntersectionVoxels[3] = NULL;
1211 case TRANSVERSE_IN_Z:
1214 if( (Ix >= 0) && (Ix + 1 < m_NumberOfVoxelsInX) &&
1215 (Iy >= 0) && (Iy + 1 < m_NumberOfVoxelsInY) &&
1216 (Iz >= 0) && (Iz < m_NumberOfVoxelsInZ))
1219 index[0]=Ix; index[1]=Iy; index[2]=Iz;
1220 m_RayIntersectionVoxels[0] = ( this->m_Image->GetBufferPointer()
1221 + this->m_Image->ComputeOffset(index) );
1223 index[0]=Ix+1; index[1]=Iy; index[2]=Iz;
1224 m_RayIntersectionVoxels[1] = ( this->m_Image->GetBufferPointer()
1225 + this->m_Image->ComputeOffset(index) );
1227 index[0]=Ix; index[1]=Iy+1; index[2]=Iz;
1228 m_RayIntersectionVoxels[2] = ( this->m_Image->GetBufferPointer()
1229 + this->m_Image->ComputeOffset(index) );
1231 index[0]=Ix+1; index[1]=Iy+1; index[2]=Iz;
1232 m_RayIntersectionVoxels[3] = ( this->m_Image->GetBufferPointer()
1233 + this->m_Image->ComputeOffset(index) );
1238 m_RayIntersectionVoxels[0]
1239 = m_RayIntersectionVoxels[1]
1240 = m_RayIntersectionVoxels[2]
1241 = m_RayIntersectionVoxels[3] = NULL;
1248 itk::ExceptionObject err(__FILE__, __LINE__);
1249 err.SetLocation( ITK_LOCATION );
1250 err.SetDescription( "The ray traversal direction is unset "
1251 "- InitialiseVoxelPointers().");
1258 /* -----------------------------------------------------------------------
1259 IncrementVoxelPointers() - Increment the voxel pointers
1260 ----------------------------------------------------------------------- */
1262 template<class TInputImage, class TCoordRep>
1264 RayCastHelper<TInputImage, TCoordRep>
1265 ::IncrementVoxelPointers(void)
1267 double xBefore = m_Position3Dvox[0];
1268 double yBefore = m_Position3Dvox[1];
1269 double zBefore = m_Position3Dvox[2];
1271 m_Position3Dvox[0] += m_VoxelIncrement[0];
1272 m_Position3Dvox[1] += m_VoxelIncrement[1];
1273 m_Position3Dvox[2] += m_VoxelIncrement[2];
1275 int dx = ((int) m_Position3Dvox[0]) - ((int) xBefore);
1276 int dy = ((int) m_Position3Dvox[1]) - ((int) yBefore);
1277 int dz = ((int) m_Position3Dvox[2]) - ((int) zBefore);
1279 m_RayIntersectionVoxelIndex[0] += dx;
1280 m_RayIntersectionVoxelIndex[1] += dy;
1281 m_RayIntersectionVoxelIndex[2] += dz;
1283 int totalRayVoxelPlanes
1284 = dx + dy*m_NumberOfVoxelsInX + dz*m_NumberOfVoxelsInX*m_NumberOfVoxelsInY;
1286 m_RayIntersectionVoxels[0] += totalRayVoxelPlanes;
1287 m_RayIntersectionVoxels[1] += totalRayVoxelPlanes;
1288 m_RayIntersectionVoxels[2] += totalRayVoxelPlanes;
1289 m_RayIntersectionVoxels[3] += totalRayVoxelPlanes;
1293 /* -----------------------------------------------------------------------
1294 GetCurrentIntensity() - Get the intensity of the current ray point.
1295 ----------------------------------------------------------------------- */
1297 template<class TInputImage, class TCoordRep>
1299 RayCastHelper<TInputImage, TCoordRep>
1300 ::GetCurrentIntensity(void) const
1309 a = (double) (*m_RayIntersectionVoxels[0]);
1310 b = (double) (*m_RayIntersectionVoxels[1] - a);
1311 c = (double) (*m_RayIntersectionVoxels[2] - a);
1312 d = (double) (*m_RayIntersectionVoxels[3] - a - b - c);
1314 switch( m_TraversalDirection )
1316 case TRANSVERSE_IN_X:
1318 y = m_Position3Dvox[1] - vcl_floor(m_Position3Dvox[1]);
1319 z = m_Position3Dvox[2] - vcl_floor(m_Position3Dvox[2]);
1322 case TRANSVERSE_IN_Y:
1324 y = m_Position3Dvox[0] - vcl_floor(m_Position3Dvox[0]);
1325 z = m_Position3Dvox[2] - vcl_floor(m_Position3Dvox[2]);
1328 case TRANSVERSE_IN_Z:
1330 y = m_Position3Dvox[0] - vcl_floor(m_Position3Dvox[0]);
1331 z = m_Position3Dvox[1] - vcl_floor(m_Position3Dvox[1]);
1336 itk::ExceptionObject err(__FILE__, __LINE__);
1337 err.SetLocation( ITK_LOCATION );
1338 err.SetDescription( "The ray traversal direction is unset "
1339 "- GetCurrentIntensity().");
1345 return a + b*y + c*z + d*y*z;
1348 /* -----------------------------------------------------------------------
1349 IncrementIntensities() - Increment the intensities of the current ray point
1350 ----------------------------------------------------------------------- */
1352 template<class TInputImage, class TCoordRep>
1354 RayCastHelper<TInputImage, TCoordRep>
1355 ::IncrementIntensities(double increment)
1357 short inc = (short) vcl_floor(increment + 0.5);
1363 *m_RayIntersectionVoxels[0] += inc;
1364 *m_RayIntersectionVoxels[1] += inc;
1365 *m_RayIntersectionVoxels[2] += inc;
1366 *m_RayIntersectionVoxels[3] += inc;
1372 /* -----------------------------------------------------------------------
1373 IntegrateAboveThreshold() - Integrate intensities above a threshold.
1374 ----------------------------------------------------------------------- */
1376 template<class TInputImage, class TCoordRep>
1378 RayCastHelper<TInputImage, TCoordRep>
1379 ::IntegrateAboveThreshold(double &integral, double threshold)
1382 // double posn3D_x, posn3D_y, posn3D_z;
1386 // Check if this is a valid ray
1392 /* Step along the ray as quickly as possible
1393 integrating the interpolated intensities. */
1395 for (m_NumVoxelPlanesTraversed=0;
1396 m_NumVoxelPlanesTraversed<m_TotalRayVoxelPlanes;
1397 m_NumVoxelPlanesTraversed++)
1399 intensity = this->GetCurrentIntensity();
1401 if (intensity > threshold)
1403 integral += intensity - threshold;
1405 this->IncrementVoxelPointers();
1408 /* The ray passes through the volume one plane of voxels at a time,
1409 however, if its moving diagonally the ray points will be further
1410 apart so account for this by scaling by the distance moved. */
1412 integral *= this->GetRayPointSpacing();
1417 /* -----------------------------------------------------------------------
1418 ZeroState() - Set the default (zero) state of the object
1419 ----------------------------------------------------------------------- */
1421 template<class TInputImage, class TCoordRep>
1423 RayCastHelper<TInputImage, TCoordRep>
1430 m_NumberOfVoxelsInX = 0;
1431 m_NumberOfVoxelsInY = 0;
1432 m_NumberOfVoxelsInZ = 0;
1434 m_VoxelDimensionInX = 0;
1435 m_VoxelDimensionInY = 0;
1436 m_VoxelDimensionInZ = 0;
1440 m_CurrentRayPositionInMM[i] = 0.;
1444 m_RayDirectionInMM[i] = 0.;
1448 m_RayVoxelStartPosition[i] = 0.;
1452 m_RayVoxelEndPosition[i] = 0.;
1456 m_VoxelIncrement[i] = 0.;
1458 m_TraversalDirection = UNDEFINED_DIRECTION;
1460 m_TotalRayVoxelPlanes = 0;
1461 m_NumVoxelPlanesTraversed = -1;
1465 m_RayIntersectionVoxels[i] = 0;
1469 m_RayIntersectionVoxelIndex[i] = 0;
1472 }; // end of anonymous namespace
1478 /**************************************************************************
1481 * Rest of this code is the actual RayCastInterpolateImageFunctionWithOrigin
1485 **************************************************************************/
1487 /* -----------------------------------------------------------------------
1489 ----------------------------------------------------------------------- */
1491 template<class TInputImage, class TCoordRep>
1492 RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
1493 ::RayCastInterpolateImageFunctionWithOrigin()
1497 m_FocalPoint[0] = 0.;
1498 m_FocalPoint[1] = 0.;
1499 m_FocalPoint[2] = 0.;
1503 /* -----------------------------------------------------------------------
1505 ----------------------------------------------------------------------- */
1507 template<class TInputImage, class TCoordRep>
1509 RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
1510 ::PrintSelf(std::ostream& os, Indent indent) const
1512 this->Superclass::PrintSelf(os,indent);
1514 os << indent << "Threshold: " << m_Threshold << std::endl;
1515 os << indent << "FocalPoint: " << m_FocalPoint << std::endl;
1516 os << indent << "Transform: " << m_Transform.GetPointer() << std::endl;
1517 os << indent << "Interpolator: " << m_Interpolator.GetPointer() << std::endl;
1521 /* -----------------------------------------------------------------------
1522 Evaluate at image index position
1523 ----------------------------------------------------------------------- */
1525 template<class TInputImage, class TCoordRep>
1526 typename RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
1528 RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
1529 ::Evaluate( const PointType& point ) const
1531 double integral = 0;
1533 OutputPointType transformedFocalPoint
1534 = m_Transform->TransformPoint( m_FocalPoint );
1536 DirectionType direction = transformedFocalPoint - point;
1538 RayCastHelper<TInputImage, TCoordRep> ray;
1539 ray.SetImage( this->m_Image );
1543 // (Aviv) Added support for images with non-zero origin
1544 // ray.SetRay(point, direction);
1545 ray.SetRay(point - this->m_Image->GetOrigin().GetVectorFromOrigin(), direction);
1546 ray.IntegrateAboveThreshold(integral, m_Threshold);
1548 return ( static_cast<OutputType>( integral ));
1551 template<class TInputImage, class TCoordRep>
1552 typename RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
1554 RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
1555 ::EvaluateAtContinuousIndex( const ContinuousIndexType& index ) const
1557 OutputPointType point;
1558 this->m_Image->TransformContinuousIndexToPhysicalPoint(index, point);
1560 return this->Evaluate( point );