1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
5 - University of LYON http://www.universite-lyon.fr/
6 - Léon Bérard cancer center http://www.centreleonberard.fr
7 - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
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.
13 It is distributed under dual licence
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"
22 #include "vnl/vnl_math.h"
25 // Put the helper class in an anonymous namespace so that it is not
26 // 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) {
69 * Initialise the ray using the position and direction of a line.
71 * \param RayPosn The position of the ray in 3D (mm).
72 * \param RayDirn The direction of the ray in 3D (mm).
74 * \return True if this is a valid ray.
76 bool SetRay(OutputPointType RayPosn, DirectionType RayDirn);
80 * Integrate the interpolated intensities along the ray and
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
88 * \param integral The integrated intensities along the ray.
90 * \return True if a valid ray was specified.
92 bool Integrate(double &integral) {
93 return IntegrateAboveThreshold(integral, 0);
98 * Integrate the interpolated intensities above a given threshold,
99 * along the ray and return the result.
101 * This routine can be called after instantiating the ray and
102 * calling SetProjectionCoord2D() or Reset(). It may then be called
103 * as many times thereafter for different 2D projection
106 * \param integral The integrated intensities along the ray.
107 * \param threshold The integration threshold [default value: 0]
109 * \return True if a valid ray was specified.
111 bool IntegrateAboveThreshold(double &integral, double threshold);
114 * Increment each of the intensities of the 4 planar voxels
115 * surrounding the current ray point.
117 * \parameter increment Intensity increment for each of the current 4 voxels
119 void IncrementIntensities(double increment=1);
121 /// Reset the iterator to the start of the ray.
124 /// Return the interpolated intensity of the current ray point.
125 double GetCurrentIntensity(void) const;
127 /// Return the ray point spacing in mm
128 double GetRayPointSpacing(void) const {
129 typename InputImageType::SpacingType spacing=this->m_Image->GetSpacing();
132 return vcl_sqrt(m_VoxelIncrement[0]*spacing[0]*m_VoxelIncrement[0]*spacing[0]
133 + m_VoxelIncrement[1]*spacing[1]*m_VoxelIncrement[1]*spacing[1]
134 + m_VoxelIncrement[2]*spacing[2]*m_VoxelIncrement[2]*spacing[2] );
139 /// Set the initial zero state of the object
142 /// Initialise the object
143 void Initialise(void);
146 /// Calculate the endpoint coordinats of the ray in voxels.
147 void EndPointsInVoxels(void);
150 * Calculate the incremental direction vector in voxels, 'dVoxel',
151 * required to traverse the ray.
153 void CalcDirnVector(void);
156 * Reduce the length of the ray until both start and end
157 * coordinates lie inside the volume.
159 * \return True if a valid ray has been, false otherwise.
161 bool AdjustRayLength(void);
164 * Obtain pointers to the four voxels surrounding the point where the ray
167 void InitialiseVoxelPointers(void);
169 /// Increment the voxel pointers surrounding the current point on the ray.
170 void IncrementVoxelPointers(void);
172 /// Record volume dimensions and resolution
173 void RecordVolumeDimensions(void);
175 /// Define the corners of the volume
176 void DefineCorners(void);
179 * Calculate the planes which define the volume.
181 * Member function to calculate the equations of the planes of 4 of
182 * the sides of the volume, calculate the positions of the 8 corners
183 * of the volume in mm in World, also calculate the values of the
184 * slopes of the lines which go to make up the volume( defined as
185 * lines in cube x,y,z dirn and then each of these lines has a slope
186 * in the world x,y,z dirn [3]) and finally also to return the length
187 * of the sides of the lines in mm.
189 void CalcPlanesAndCorners(void);
192 * Calculate the ray intercepts with the volume.
194 * See where the ray cuts the volume, check that truncation does not occur,
195 * if not, then start ray where it first intercepts the volume and set
196 * x_max to be where it leaves the volume.
198 * \return True if a valid ray has been specified, false otherwise.
200 bool CalcRayIntercepts(void);
203 * The ray is traversed by stepping in the axial direction
204 * that enables the greatest number of planes in the volume to be
208 UNDEFINED_DIRECTION=0, //!< Undefined
209 TRANSVERSE_IN_X, //!< x
210 TRANSVERSE_IN_Y, //!< y
211 TRANSVERSE_IN_Z, //!< z
213 } TraversalDirection;
215 // Cache the image in the structure. Skip the smart pointer for
216 // efficiency. This inner class will go in/out of scope with every
217 // call to Evaluate()
218 const InputImageType *m_Image;
220 /// Flag indicating whether the current ray is valid
224 * The start position of the ray in voxels.
226 * NB. Two of the components of this coordinate (i.e. those lying within
227 * the planes of voxels being traversed) will be shifted by half a
228 * voxel. This enables indices of the neighbouring voxels within the plane
229 * to be determined by simply casting to 'int' and optionally adding 1.
231 double m_RayVoxelStartPosition[3];
234 * The end coordinate of the ray in voxels.
236 * NB. Two of the components of this coordinate (i.e. those lying within
237 * the planes of voxels being traversed) will be shifted by half a
238 * voxel. This enables indices of the neighbouring voxels within the plane
239 * to be determined by simply casting to 'int' and optionally adding 1.
241 double m_RayVoxelEndPosition[3];
245 * The current coordinate on the ray in voxels.
247 * NB. Two of the components of this coordinate (i.e. those lying within
248 * the planes of voxels being traversed) will be shifted by half a
249 * voxel. This enables indices of the neighbouring voxels within the plane
250 * to be determined by simply casting to 'int' and optionally adding 1.
252 double m_Position3Dvox[3];
254 /** The incremental direction vector of the ray in voxels. */
255 double m_VoxelIncrement[3];
257 /// The direction in which the ray is incremented thorough the volume (x, y or z).
258 TraversalDirection m_TraversalDirection;
260 /// The total number of planes of voxels traversed by the ray.
261 int m_TotalRayVoxelPlanes;
263 /// The current number of planes of voxels traversed by the ray.
264 int m_NumVoxelPlanesTraversed;
266 /// Pointers to the current four voxels surrounding the ray's trajectory.
267 const PixelType *m_RayIntersectionVoxels[4];
270 * The voxel coordinate of the bottom-left voxel of the current
271 * four voxels surrounding the ray's trajectory.
273 int m_RayIntersectionVoxelIndex[3];
275 /// The dimension in voxels of the 3D volume in along the x axis
276 int m_NumberOfVoxelsInX;
277 /// The dimension in voxels of the 3D volume in along the y axis
278 int m_NumberOfVoxelsInY;
279 /// The dimension in voxels of the 3D volume in along the z axis
280 int m_NumberOfVoxelsInZ;
282 /// Voxel dimension in x
283 double m_VoxelDimensionInX;
284 /// Voxel dimension in y
285 double m_VoxelDimensionInY;
286 /// Voxel dimension in z
287 double m_VoxelDimensionInZ;
289 /// The coordinate of the point at which the ray enters the volume in mm.
290 double m_RayStartCoordInMM[3];
291 /// The coordinate of the point at which the ray exits the volume in mm.
292 double m_RayEndCoordInMM[3];
296 Planes which define the boundary of the volume in mm
297 (six planes and four parameters: Ax+By+Cz+D). */
298 double m_BoundingPlane[6][4];
299 /// The eight corners of the volume (x,y,z coordinates for each).
300 double m_BoundingCorner[8][3];
302 /// The position of the ray
303 double m_CurrentRayPositionInMM[3];
305 /// The direction of the ray
306 double m_RayDirectionInMM[3];
310 /* -----------------------------------------------------------------------
311 Initialise() - Initialise the object
312 ----------------------------------------------------------------------- */
314 template<class TInputImage, class TCoordRep>
316 RayCastHelper<TInputImage, TCoordRep>
319 // Save the dimensions of the volume and calculate the bounding box
320 this->RecordVolumeDimensions();
322 // Calculate the planes and corners which define the volume.
323 this->DefineCorners();
324 this->CalcPlanesAndCorners();
328 /* -----------------------------------------------------------------------
329 RecordVolumeDimensions() - Record volume dimensions and resolution
330 ----------------------------------------------------------------------- */
332 template<class TInputImage, class TCoordRep>
334 RayCastHelper<TInputImage, TCoordRep>
335 ::RecordVolumeDimensions(void)
337 typename InputImageType::SpacingType spacing=this->m_Image->GetSpacing();
338 SizeType dim=this->m_Image->GetLargestPossibleRegion().GetSize();
340 m_NumberOfVoxelsInX = dim[0];
341 m_NumberOfVoxelsInY = dim[1];
342 m_NumberOfVoxelsInZ = dim[2];
344 m_VoxelDimensionInX = spacing[0];
345 m_VoxelDimensionInY = spacing[1];
346 m_VoxelDimensionInZ = spacing[2];
350 /* -----------------------------------------------------------------------
351 DefineCorners() - Define the corners of the volume
352 ----------------------------------------------------------------------- */
354 template<class TInputImage, class TCoordRep>
356 RayCastHelper<TInputImage, TCoordRep>
357 ::DefineCorners(void)
359 // Define corner positions as if at the origin
361 m_BoundingCorner[0][0] =
362 m_BoundingCorner[1][0] =
363 m_BoundingCorner[2][0] =
364 m_BoundingCorner[3][0] = 0;
366 m_BoundingCorner[4][0] =
367 m_BoundingCorner[5][0] =
368 m_BoundingCorner[6][0] =
369 m_BoundingCorner[7][0] = m_VoxelDimensionInX*m_NumberOfVoxelsInX;
371 m_BoundingCorner[1][1] =
372 m_BoundingCorner[3][1] =
373 m_BoundingCorner[5][1] =
374 m_BoundingCorner[7][1] = m_VoxelDimensionInY*m_NumberOfVoxelsInY;
376 m_BoundingCorner[0][1] =
377 m_BoundingCorner[2][1] =
378 m_BoundingCorner[4][1] =
379 m_BoundingCorner[6][1] = 0;
381 m_BoundingCorner[0][2] =
382 m_BoundingCorner[1][2] =
383 m_BoundingCorner[4][2] =
384 m_BoundingCorner[5][2] =
385 m_VoxelDimensionInZ*m_NumberOfVoxelsInZ;
387 m_BoundingCorner[2][2] =
388 m_BoundingCorner[3][2] =
389 m_BoundingCorner[6][2] =
390 m_BoundingCorner[7][2] = 0;
394 /* -----------------------------------------------------------------------
395 CalcPlanesAndCorners() - Calculate the planes and corners of the volume.
396 ----------------------------------------------------------------------- */
398 template<class TInputImage, class TCoordRep>
400 RayCastHelper<TInputImage, TCoordRep>
401 ::CalcPlanesAndCorners(void)
406 // find the equations of the planes
408 int c1=0, c2=0, c3=0;
410 for (j=0; j<6; j++) {
411 // loop around for planes
413 // which corners to take
447 double line1x, line1y, line1z;
448 double line2x, line2y, line2z;
450 // lines from one corner to another in x,y,z dirns
451 line1x = m_BoundingCorner[c1][0] - m_BoundingCorner[c2][0];
452 line2x = m_BoundingCorner[c1][0] - m_BoundingCorner[c3][0];
454 line1y = m_BoundingCorner[c1][1] - m_BoundingCorner[c2][1];
455 line2y = m_BoundingCorner[c1][1] - m_BoundingCorner[c3][1];
457 line1z = m_BoundingCorner[c1][2] - m_BoundingCorner[c2][2];
458 line2z = m_BoundingCorner[c1][2] - m_BoundingCorner[c3][2];
462 // take cross product
463 A = line1y*line2z - line2y*line1z;
464 B = line2x*line1z - line1x*line2z;
465 C = line1x*line2y - line2x*line1y;
468 D = -( A*m_BoundingCorner[c1][0]
469 + B*m_BoundingCorner[c1][1]
470 + C*m_BoundingCorner[c1][2] );
472 // initialise plane value and normalise
473 m_BoundingPlane[j][0] = A/vcl_sqrt(A*A + B*B + C*C);
474 m_BoundingPlane[j][1] = B/vcl_sqrt(A*A + B*B + C*C);
475 m_BoundingPlane[j][2] = C/vcl_sqrt(A*A + B*B + C*C);
476 m_BoundingPlane[j][3] = D/vcl_sqrt(A*A + B*B + C*C);
478 if ( (A*A + B*B + C*C) == 0 ) {
479 itk::ExceptionObject err(__FILE__, __LINE__);
480 err.SetLocation( ITK_LOCATION );
481 err.SetDescription( "Division by zero (planes) "
482 "- CalcPlanesAndCorners().");
490 /* -----------------------------------------------------------------------
491 CalcRayIntercepts() - Calculate the ray intercepts with the volume.
492 ----------------------------------------------------------------------- */
494 template<class TInputImage, class TCoordRep>
496 RayCastHelper<TInputImage, TCoordRep>
497 ::CalcRayIntercepts()
499 double maxInterDist, interDist;
500 double cornerVect[4][3];
501 int cross[4][3], noInterFlag[6];
502 int nSidesCrossed, crossFlag, c[4];
503 double ax, ay, az, bx, by, bz;
504 double cubeInter[6][3];
508 int NoSides = 6; // =6 to allow truncation: =4 to remove truncated rays
510 // Calculate intercept of ray with planes
512 double interceptx[6], intercepty[6], interceptz[6];
515 for( j=0; j<NoSides; j++) {
517 denom = ( m_BoundingPlane[j][0]*m_RayDirectionInMM[0]
518 + m_BoundingPlane[j][1]*m_RayDirectionInMM[1]
519 + m_BoundingPlane[j][2]*m_RayDirectionInMM[2]);
521 if( (long)(denom*100) != 0 ) {
522 d[j] = -( m_BoundingPlane[j][3]
523 + m_BoundingPlane[j][0]*m_CurrentRayPositionInMM[0]
524 + m_BoundingPlane[j][1]*m_CurrentRayPositionInMM[1]
525 + m_BoundingPlane[j][2]*m_CurrentRayPositionInMM[2] ) / denom;
527 interceptx[j] = m_CurrentRayPositionInMM[0] + d[j]*m_RayDirectionInMM[0];
528 intercepty[j] = m_CurrentRayPositionInMM[1] + d[j]*m_RayDirectionInMM[1];
529 interceptz[j] = m_CurrentRayPositionInMM[2] + d[j]*m_RayDirectionInMM[2];
531 noInterFlag[j] = 1; //OK
533 noInterFlag[j] = 0; //NOT OK
539 for( j=0; j<NoSides; j++ ) {
541 // Work out which corners to use
577 // Calculate vectors from corner of ct volume to intercept.
578 for( i=0; i<4; i++ ) {
579 if( noInterFlag[j]==1 ) {
580 cornerVect[i][0] = m_BoundingCorner[c[i]][0] - interceptx[j];
581 cornerVect[i][1] = m_BoundingCorner[c[i]][1] - intercepty[j];
582 cornerVect[i][2] = m_BoundingCorner[c[i]][2] - interceptz[j];
583 } else if( noInterFlag[j]==0 ) {
584 cornerVect[i][0] = 0;
585 cornerVect[i][1] = 0;
586 cornerVect[i][2] = 0;
591 // Do cross product with these vectors
592 for( i=0; i<4; i++ ) {
598 ax = cornerVect[i][0];
599 ay = cornerVect[i][1];
600 az = cornerVect[i][2];
601 bx = cornerVect[k][0];
602 by = cornerVect[k][1];
603 bz = cornerVect[k][2];
605 // The int and divide by 100 are to avoid rounding errors. If
606 // these are not included then you get values fluctuating around
607 // zero and so in the subsequent check, all the values are not
608 // above or below zero. NB. If you "INT" by too much here though
609 // you can get problems in the corners of your volume when rays
610 // are allowed to go through more than one plane.
611 cross[i][0] = (int)((ay*bz - az*by)/100);
612 cross[i][1] = (int)((az*bx - ax*bz)/100);
613 cross[i][2] = (int)((ax*by - ay*bx)/100);
616 // See if a sign change occured between all these cross products
617 // if not, then the ray went through this plane
620 for( i=0; i<3; i++ ) {
629 && cross[3][i]>=0) ) {
635 if( crossFlag==3 && noInterFlag[j]==1 ) {
636 cubeInter[nSidesCrossed][0] = interceptx[j];
637 cubeInter[nSidesCrossed][1] = intercepty[j];
638 cubeInter[nSidesCrossed][2] = interceptz[j];
642 } // End of loop over all four planes
644 m_RayStartCoordInMM[0] = cubeInter[0][0];
645 m_RayStartCoordInMM[1] = cubeInter[0][1];
646 m_RayStartCoordInMM[2] = cubeInter[0][2];
648 m_RayEndCoordInMM[0] = cubeInter[1][0];
649 m_RayEndCoordInMM[1] = cubeInter[1][1];
650 m_RayEndCoordInMM[2] = cubeInter[1][2];
652 if( nSidesCrossed >= 5 ) {
653 std::cerr << "WARNING: No. of sides crossed equals: " << nSidesCrossed << std::endl;
656 // If 'nSidesCrossed' is larger than 2, this means that the ray goes through
657 // a corner of the volume and due to rounding errors, the ray is
658 // deemed to go through more than two planes. To obtain the correct
659 // start and end positions we choose the two intercept values which
660 // are furthest from each other.
663 if( nSidesCrossed >= 3 ) {
665 for( j=0; j<nSidesCrossed-1; j++ ) {
666 for( k=j+1; k<nSidesCrossed; k++ ) {
668 for( i=0; i<3; i++ ) {
669 interDist += (cubeInter[j][i] - cubeInter[k][i])*
670 (cubeInter[j][i] - cubeInter[k][i]);
672 if( interDist > maxInterDist ) {
673 maxInterDist = interDist;
675 m_RayStartCoordInMM[0] = cubeInter[j][0];
676 m_RayStartCoordInMM[1] = cubeInter[j][1];
677 m_RayStartCoordInMM[2] = cubeInter[j][2];
679 m_RayEndCoordInMM[0] = cubeInter[k][0];
680 m_RayEndCoordInMM[1] = cubeInter[k][1];
681 m_RayEndCoordInMM[2] = cubeInter[k][2];
688 if (nSidesCrossed == 2 ) {
696 /* -----------------------------------------------------------------------
697 SetRay() - Set the position and direction of the ray
698 ----------------------------------------------------------------------- */
700 template<class TInputImage, class TCoordRep>
702 RayCastHelper<TInputImage, TCoordRep>
703 ::SetRay(OutputPointType RayPosn, DirectionType RayDirn)
706 // Store the position and direction of the ray
707 typename TInputImage::SpacingType spacing=this->m_Image->GetSpacing();
708 SizeType dim=this->m_Image->GetLargestPossibleRegion().GetSize();
710 // we need to translate the _center_ of the volume to the origin
711 m_NumberOfVoxelsInX = dim[0];
712 m_NumberOfVoxelsInY = dim[1];
713 m_NumberOfVoxelsInZ = dim[2];
715 m_VoxelDimensionInX = spacing[0];
716 m_VoxelDimensionInY = spacing[1];
717 m_VoxelDimensionInZ = spacing[2];
719 // (Aviv) Incorporating a fix by Jian Wu
720 // http://public.kitware.com/pipermail/insight-users/2006-March/017265.html
721 m_CurrentRayPositionInMM[0] = RayPosn[0];
722 m_CurrentRayPositionInMM[1] = RayPosn[1];
723 m_CurrentRayPositionInMM[2] = RayPosn[2];
725 m_RayDirectionInMM[0] = RayDirn[0];
726 m_RayDirectionInMM[1] = RayDirn[1];
727 m_RayDirectionInMM[2] = RayDirn[2];
729 // Compute the ray path for this coordinate in mm
731 m_ValidRay = this->CalcRayIntercepts();
738 // Convert the start and end coordinates of the ray to voxels
740 this->EndPointsInVoxels();
742 /* Calculate the ray direction vector in voxels and hence the voxel
743 increment required to traverse the ray, and the number of
744 interpolation points on the ray.
746 This routine also shifts the coordinate frame by half a voxel for
747 two of the directional components (i.e. those lying within the
748 planes of voxels being traversed). */
750 this->CalcDirnVector();
753 /* Reduce the length of the ray until both start and end
754 coordinates lie inside the volume. */
756 m_ValidRay = this->AdjustRayLength();
758 // Reset the iterator to the start of the ray.
766 /* -----------------------------------------------------------------------
767 EndPointsInVoxels() - Convert the endpoints to voxels
768 ----------------------------------------------------------------------- */
770 template<class TInputImage, class TCoordRep>
772 RayCastHelper<TInputImage, TCoordRep>
773 ::EndPointsInVoxels(void)
775 m_RayVoxelStartPosition[0] = m_RayStartCoordInMM[0]/m_VoxelDimensionInX;
776 m_RayVoxelStartPosition[1] = m_RayStartCoordInMM[1]/m_VoxelDimensionInY;
777 m_RayVoxelStartPosition[2] = m_RayStartCoordInMM[2]/m_VoxelDimensionInZ;
779 m_RayVoxelEndPosition[0] = m_RayEndCoordInMM[0]/m_VoxelDimensionInX;
780 m_RayVoxelEndPosition[1] = m_RayEndCoordInMM[1]/m_VoxelDimensionInY;
781 m_RayVoxelEndPosition[2] = m_RayEndCoordInMM[2]/m_VoxelDimensionInZ;
786 /* -----------------------------------------------------------------------
787 CalcDirnVector() - Calculate the incremental direction vector in voxels.
788 ----------------------------------------------------------------------- */
790 template<class TInputImage, class TCoordRep>
792 RayCastHelper<TInputImage, TCoordRep>
793 ::CalcDirnVector(void)
795 double xNum, yNum, zNum;
797 // Calculate the number of voxels in each direction
799 xNum = vcl_fabs(m_RayVoxelStartPosition[0] - m_RayVoxelEndPosition[0]);
800 yNum = vcl_fabs(m_RayVoxelStartPosition[1] - m_RayVoxelEndPosition[1]);
801 zNum = vcl_fabs(m_RayVoxelStartPosition[2] - m_RayVoxelEndPosition[2]);
803 // The direction iterated in is that with the greatest number of voxels
804 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
806 // Iterate in X direction
808 if( (xNum >= yNum) && (xNum >= zNum) ) {
809 if( m_RayVoxelStartPosition[0] < m_RayVoxelEndPosition[0] ) {
810 m_VoxelIncrement[0] = 1;
813 = (m_RayVoxelStartPosition[1]
814 - m_RayVoxelEndPosition[1])/(m_RayVoxelStartPosition[0]
815 - m_RayVoxelEndPosition[0]);
818 = (m_RayVoxelStartPosition[2]
819 - m_RayVoxelEndPosition[2])/(m_RayVoxelStartPosition[0]
820 - m_RayVoxelEndPosition[0]);
822 m_VoxelIncrement[0] = -1;
825 = -(m_RayVoxelStartPosition[1]
826 - m_RayVoxelEndPosition[1])/(m_RayVoxelStartPosition[0]
827 - m_RayVoxelEndPosition[0]);
830 = -(m_RayVoxelStartPosition[2]
831 - m_RayVoxelEndPosition[2])/(m_RayVoxelStartPosition[0]
832 - m_RayVoxelEndPosition[0]);
835 // This section is to alter the start position in order to
836 // place the center of the voxels in there correct positions,
837 // rather than placing them at the corner of voxels which is
838 // what happens if this is not carried out. The reason why
839 // x has no -0.5 is because this is the direction we are going
840 // to iterate in and therefore we wish to go from center to
841 // center rather than finding the surrounding voxels.
843 m_RayVoxelStartPosition[1] += ( (int)m_RayVoxelStartPosition[0]
844 - m_RayVoxelStartPosition[0])*m_VoxelIncrement[1]*m_VoxelIncrement[0]
845 + 0.5*m_VoxelIncrement[1] - 0.5;
847 m_RayVoxelStartPosition[2] += ( (int)m_RayVoxelStartPosition[0]
848 - m_RayVoxelStartPosition[0])*m_VoxelIncrement[2]*m_VoxelIncrement[0]
849 + 0.5*m_VoxelIncrement[2] - 0.5;
851 m_RayVoxelStartPosition[0] = (int)m_RayVoxelStartPosition[0] + 0.5*m_VoxelIncrement[0];
853 m_TotalRayVoxelPlanes = (int)xNum;
855 m_TraversalDirection = TRANSVERSE_IN_X;
858 // Iterate in Y direction
860 else if( (yNum >= xNum) && (yNum >= zNum) ) {
862 if( m_RayVoxelStartPosition[1] < m_RayVoxelEndPosition[1] ) {
863 m_VoxelIncrement[1] = 1;
866 = (m_RayVoxelStartPosition[0]
867 - m_RayVoxelEndPosition[0])/(m_RayVoxelStartPosition[1]
868 - m_RayVoxelEndPosition[1]);
871 = (m_RayVoxelStartPosition[2]
872 - m_RayVoxelEndPosition[2])/(m_RayVoxelStartPosition[1]
873 - m_RayVoxelEndPosition[1]);
875 m_VoxelIncrement[1] = -1;
878 = -(m_RayVoxelStartPosition[0]
879 - m_RayVoxelEndPosition[0])/(m_RayVoxelStartPosition[1]
880 - m_RayVoxelEndPosition[1]);
883 = -(m_RayVoxelStartPosition[2]
884 - m_RayVoxelEndPosition[2])/(m_RayVoxelStartPosition[1]
885 - m_RayVoxelEndPosition[1]);
888 m_RayVoxelStartPosition[0] += ( (int)m_RayVoxelStartPosition[1]
889 - m_RayVoxelStartPosition[1])*m_VoxelIncrement[0]*m_VoxelIncrement[1]
890 + 0.5*m_VoxelIncrement[0] - 0.5;
892 m_RayVoxelStartPosition[2] += ( (int)m_RayVoxelStartPosition[1]
893 - m_RayVoxelStartPosition[1])*m_VoxelIncrement[2]*m_VoxelIncrement[1]
894 + 0.5*m_VoxelIncrement[2] - 0.5;
896 m_RayVoxelStartPosition[1] = (int)m_RayVoxelStartPosition[1] + 0.5*m_VoxelIncrement[1];
898 m_TotalRayVoxelPlanes = (int)yNum;
900 m_TraversalDirection = TRANSVERSE_IN_Y;
903 // Iterate in Z direction
907 if( m_RayVoxelStartPosition[2] < m_RayVoxelEndPosition[2] ) {
908 m_VoxelIncrement[2] = 1;
911 = (m_RayVoxelStartPosition[0]
912 - m_RayVoxelEndPosition[0])/(m_RayVoxelStartPosition[2]
913 - m_RayVoxelEndPosition[2]);
916 = (m_RayVoxelStartPosition[1]
917 - m_RayVoxelEndPosition[1])/(m_RayVoxelStartPosition[2]
918 - m_RayVoxelEndPosition[2]);
920 m_VoxelIncrement[2] = -1;
923 = -(m_RayVoxelStartPosition[0]
924 - m_RayVoxelEndPosition[0])/(m_RayVoxelStartPosition[2]
925 - m_RayVoxelEndPosition[2]);
928 = -(m_RayVoxelStartPosition[1]
929 - m_RayVoxelEndPosition[1])/(m_RayVoxelStartPosition[2]
930 - m_RayVoxelEndPosition[2]);
933 m_RayVoxelStartPosition[0] += ( (int)m_RayVoxelStartPosition[2]
934 - m_RayVoxelStartPosition[2])*m_VoxelIncrement[0]*m_VoxelIncrement[2]
935 + 0.5*m_VoxelIncrement[0] - 0.5;
937 m_RayVoxelStartPosition[1] += ( (int)m_RayVoxelStartPosition[2]
938 - m_RayVoxelStartPosition[2])*m_VoxelIncrement[1]*m_VoxelIncrement[2]
939 + 0.5*m_VoxelIncrement[1] - 0.5;
941 m_RayVoxelStartPosition[2] = (int)m_RayVoxelStartPosition[2] + 0.5*m_VoxelIncrement[2];
943 m_TotalRayVoxelPlanes = (int)zNum;
945 m_TraversalDirection = TRANSVERSE_IN_Z;
950 /* -----------------------------------------------------------------------
951 AdjustRayLength() - Ensure that the ray lies within the volume
952 ----------------------------------------------------------------------- */
954 template<class TInputImage, class TCoordRep>
956 RayCastHelper<TInputImage, TCoordRep>
957 ::AdjustRayLength(void)
964 if (m_TraversalDirection == TRANSVERSE_IN_X) {
968 } else if (m_TraversalDirection == TRANSVERSE_IN_Y) {
972 } else if (m_TraversalDirection == TRANSVERSE_IN_Z) {
977 itk::ExceptionObject err(__FILE__, __LINE__);
978 err.SetLocation( ITK_LOCATION );
979 err.SetDescription( "The ray traversal direction is unset "
980 "- AdjustRayLength().");
991 Istart[0] = (int) vcl_floor(m_RayVoxelStartPosition[0]);
992 Istart[1] = (int) vcl_floor(m_RayVoxelStartPosition[1]);
993 Istart[2] = (int) vcl_floor(m_RayVoxelStartPosition[2]);
995 if( (Istart[0] >= 0) && (Istart[0] + Idirn[0] < m_NumberOfVoxelsInX) &&
996 (Istart[1] >= 0) && (Istart[1] + Idirn[1] < m_NumberOfVoxelsInY) &&
997 (Istart[2] >= 0) && (Istart[2] + Idirn[2] < m_NumberOfVoxelsInZ) ) {
1001 m_RayVoxelStartPosition[0] += m_VoxelIncrement[0];
1002 m_RayVoxelStartPosition[1] += m_VoxelIncrement[1];
1003 m_RayVoxelStartPosition[2] += m_VoxelIncrement[2];
1005 m_TotalRayVoxelPlanes--;
1008 Istart[0] = (int) vcl_floor(m_RayVoxelStartPosition[0]
1009 + m_TotalRayVoxelPlanes*m_VoxelIncrement[0]);
1011 Istart[1] = (int) vcl_floor(m_RayVoxelStartPosition[1]
1012 + m_TotalRayVoxelPlanes*m_VoxelIncrement[1]);
1014 Istart[2] = (int) vcl_floor(m_RayVoxelStartPosition[2]
1015 + m_TotalRayVoxelPlanes*m_VoxelIncrement[2]);
1017 if( (Istart[0] >= 0) && (Istart[0] + Idirn[0] < m_NumberOfVoxelsInX) &&
1018 (Istart[1] >= 0) && (Istart[1] + Idirn[1] < m_NumberOfVoxelsInY) &&
1019 (Istart[2] >= 0) && (Istart[2] + Idirn[2] < m_NumberOfVoxelsInZ) ) {
1023 m_TotalRayVoxelPlanes--;
1026 } while ( (! (startOK && endOK)) && (m_TotalRayVoxelPlanes > 1) );
1029 return (startOK && endOK);
1033 /* -----------------------------------------------------------------------
1034 Reset() - Reset the iterator to the start of the ray.
1035 ----------------------------------------------------------------------- */
1037 template<class TInputImage, class TCoordRep>
1039 RayCastHelper<TInputImage, TCoordRep>
1044 m_NumVoxelPlanesTraversed = -1;
1046 // If this is a valid ray...
1049 for (i=0; i<3; i++) {
1050 m_Position3Dvox[i] = m_RayVoxelStartPosition[i];
1052 this->InitialiseVoxelPointers();
1055 // otherwise set parameters to zero
1058 for (i=0; i<3; i++) {
1059 m_RayVoxelStartPosition[i] = 0.;
1061 for (i=0; i<3; i++) {
1062 m_RayVoxelEndPosition[i] = 0.;
1064 for (i=0; i<3; i++) {
1065 m_VoxelIncrement[i] = 0.;
1067 m_TraversalDirection = UNDEFINED_DIRECTION;
1069 m_TotalRayVoxelPlanes = 0;
1071 for (i=0; i<4; i++) {
1072 m_RayIntersectionVoxels[i] = 0;
1074 for (i=0; i<3; i++) {
1075 m_RayIntersectionVoxelIndex[i] = 0;
1081 /* -----------------------------------------------------------------------
1082 InitialiseVoxelPointers() - Obtain pointers to the first four voxels
1083 ----------------------------------------------------------------------- */
1085 template<class TInputImage, class TCoordRep>
1087 RayCastHelper<TInputImage, TCoordRep>
1088 ::InitialiseVoxelPointers(void)
1094 Ix = (int)(m_RayVoxelStartPosition[0]);
1095 Iy = (int)(m_RayVoxelStartPosition[1]);
1096 Iz = (int)(m_RayVoxelStartPosition[2]);
1098 m_RayIntersectionVoxelIndex[0] = Ix;
1099 m_RayIntersectionVoxelIndex[1] = Iy;
1100 m_RayIntersectionVoxelIndex[2] = Iz;
1102 switch( m_TraversalDirection ) {
1103 case TRANSVERSE_IN_X: {
1105 if( (Ix >= 0) && (Ix < m_NumberOfVoxelsInX) &&
1106 (Iy >= 0) && (Iy + 1 < m_NumberOfVoxelsInY) &&
1107 (Iz >= 0) && (Iz + 1 < m_NumberOfVoxelsInZ)) {
1112 m_RayIntersectionVoxels[0]
1113 = this->m_Image->GetBufferPointer() + this->m_Image->ComputeOffset(index);
1118 m_RayIntersectionVoxels[1]
1119 = ( this->m_Image->GetBufferPointer() + this->m_Image->ComputeOffset(index) );
1124 m_RayIntersectionVoxels[2]
1125 = ( this->m_Image->GetBufferPointer() + this->m_Image->ComputeOffset(index) );
1130 m_RayIntersectionVoxels[3]
1131 = ( this->m_Image->GetBufferPointer() + this->m_Image->ComputeOffset(index) );
1133 m_RayIntersectionVoxels[0] =
1134 m_RayIntersectionVoxels[1] =
1135 m_RayIntersectionVoxels[2] =
1136 m_RayIntersectionVoxels[3] = NULL;
1141 case TRANSVERSE_IN_Y: {
1143 if( (Ix >= 0) && (Ix + 1 < m_NumberOfVoxelsInX) &&
1144 (Iy >= 0) && (Iy < m_NumberOfVoxelsInY) &&
1145 (Iz >= 0) && (Iz + 1 < m_NumberOfVoxelsInZ)) {
1150 m_RayIntersectionVoxels[0] = ( this->m_Image->GetBufferPointer()
1151 + this->m_Image->ComputeOffset(index) );
1156 m_RayIntersectionVoxels[1] = ( this->m_Image->GetBufferPointer()
1157 + this->m_Image->ComputeOffset(index) );
1162 m_RayIntersectionVoxels[2] = ( this->m_Image->GetBufferPointer()
1163 + this->m_Image->ComputeOffset(index) );
1168 m_RayIntersectionVoxels[3] = ( this->m_Image->GetBufferPointer()
1169 + this->m_Image->ComputeOffset(index) );
1171 m_RayIntersectionVoxels[0]
1172 = m_RayIntersectionVoxels[1]
1173 = m_RayIntersectionVoxels[2]
1174 = m_RayIntersectionVoxels[3] = NULL;
1179 case TRANSVERSE_IN_Z: {
1181 if( (Ix >= 0) && (Ix + 1 < m_NumberOfVoxelsInX) &&
1182 (Iy >= 0) && (Iy + 1 < m_NumberOfVoxelsInY) &&
1183 (Iz >= 0) && (Iz < m_NumberOfVoxelsInZ)) {
1188 m_RayIntersectionVoxels[0] = ( this->m_Image->GetBufferPointer()
1189 + this->m_Image->ComputeOffset(index) );
1194 m_RayIntersectionVoxels[1] = ( this->m_Image->GetBufferPointer()
1195 + this->m_Image->ComputeOffset(index) );
1200 m_RayIntersectionVoxels[2] = ( this->m_Image->GetBufferPointer()
1201 + this->m_Image->ComputeOffset(index) );
1206 m_RayIntersectionVoxels[3] = ( this->m_Image->GetBufferPointer()
1207 + this->m_Image->ComputeOffset(index) );
1210 m_RayIntersectionVoxels[0]
1211 = m_RayIntersectionVoxels[1]
1212 = m_RayIntersectionVoxels[2]
1213 = m_RayIntersectionVoxels[3] = NULL;
1219 itk::ExceptionObject err(__FILE__, __LINE__);
1220 err.SetLocation( ITK_LOCATION );
1221 err.SetDescription( "The ray traversal direction is unset "
1222 "- InitialiseVoxelPointers().");
1229 /* -----------------------------------------------------------------------
1230 IncrementVoxelPointers() - Increment the voxel pointers
1231 ----------------------------------------------------------------------- */
1233 template<class TInputImage, class TCoordRep>
1235 RayCastHelper<TInputImage, TCoordRep>
1236 ::IncrementVoxelPointers(void)
1238 double xBefore = m_Position3Dvox[0];
1239 double yBefore = m_Position3Dvox[1];
1240 double zBefore = m_Position3Dvox[2];
1242 m_Position3Dvox[0] += m_VoxelIncrement[0];
1243 m_Position3Dvox[1] += m_VoxelIncrement[1];
1244 m_Position3Dvox[2] += m_VoxelIncrement[2];
1246 int dx = ((int) m_Position3Dvox[0]) - ((int) xBefore);
1247 int dy = ((int) m_Position3Dvox[1]) - ((int) yBefore);
1248 int dz = ((int) m_Position3Dvox[2]) - ((int) zBefore);
1250 m_RayIntersectionVoxelIndex[0] += dx;
1251 m_RayIntersectionVoxelIndex[1] += dy;
1252 m_RayIntersectionVoxelIndex[2] += dz;
1254 int totalRayVoxelPlanes
1255 = dx + dy*m_NumberOfVoxelsInX + dz*m_NumberOfVoxelsInX*m_NumberOfVoxelsInY;
1257 m_RayIntersectionVoxels[0] += totalRayVoxelPlanes;
1258 m_RayIntersectionVoxels[1] += totalRayVoxelPlanes;
1259 m_RayIntersectionVoxels[2] += totalRayVoxelPlanes;
1260 m_RayIntersectionVoxels[3] += totalRayVoxelPlanes;
1264 /* -----------------------------------------------------------------------
1265 GetCurrentIntensity() - Get the intensity of the current ray point.
1266 ----------------------------------------------------------------------- */
1268 template<class TInputImage, class TCoordRep>
1270 RayCastHelper<TInputImage, TCoordRep>
1271 ::GetCurrentIntensity(void) const
1279 a = (double) (*m_RayIntersectionVoxels[0]);
1280 b = (double) (*m_RayIntersectionVoxels[1] - a);
1281 c = (double) (*m_RayIntersectionVoxels[2] - a);
1282 d = (double) (*m_RayIntersectionVoxels[3] - a - b - c);
1284 switch( m_TraversalDirection ) {
1285 case TRANSVERSE_IN_X: {
1286 y = m_Position3Dvox[1] - vcl_floor(m_Position3Dvox[1]);
1287 z = m_Position3Dvox[2] - vcl_floor(m_Position3Dvox[2]);
1290 case TRANSVERSE_IN_Y: {
1291 y = m_Position3Dvox[0] - vcl_floor(m_Position3Dvox[0]);
1292 z = m_Position3Dvox[2] - vcl_floor(m_Position3Dvox[2]);
1295 case TRANSVERSE_IN_Z: {
1296 y = m_Position3Dvox[0] - vcl_floor(m_Position3Dvox[0]);
1297 z = m_Position3Dvox[1] - vcl_floor(m_Position3Dvox[1]);
1301 itk::ExceptionObject err(__FILE__, __LINE__);
1302 err.SetLocation( ITK_LOCATION );
1303 err.SetDescription( "The ray traversal direction is unset "
1304 "- GetCurrentIntensity().");
1310 return a + b*y + c*z + d*y*z;
1313 /* -----------------------------------------------------------------------
1314 IncrementIntensities() - Increment the intensities of the current ray point
1315 ----------------------------------------------------------------------- */
1317 template<class TInputImage, class TCoordRep>
1319 RayCastHelper<TInputImage, TCoordRep>
1320 ::IncrementIntensities(double increment)
1322 short inc = (short) vcl_floor(increment + 0.5);
1327 *m_RayIntersectionVoxels[0] += inc;
1328 *m_RayIntersectionVoxels[1] += inc;
1329 *m_RayIntersectionVoxels[2] += inc;
1330 *m_RayIntersectionVoxels[3] += inc;
1336 /* -----------------------------------------------------------------------
1337 IntegrateAboveThreshold() - Integrate intensities above a threshold.
1338 ----------------------------------------------------------------------- */
1340 template<class TInputImage, class TCoordRep>
1342 RayCastHelper<TInputImage, TCoordRep>
1343 ::IntegrateAboveThreshold(double &integral, double threshold)
1346 // double posn3D_x, posn3D_y, posn3D_z;
1350 // Check if this is a valid ray
1355 /* Step along the ray as quickly as possible
1356 integrating the interpolated intensities. */
1358 for (m_NumVoxelPlanesTraversed=0;
1359 m_NumVoxelPlanesTraversed<m_TotalRayVoxelPlanes;
1360 m_NumVoxelPlanesTraversed++) {
1361 intensity = this->GetCurrentIntensity();
1363 if (intensity > threshold) {
1364 integral += intensity - threshold;
1366 this->IncrementVoxelPointers();
1369 /* The ray passes through the volume one plane of voxels at a time,
1370 however, if its moving diagonally the ray points will be further
1371 apart so account for this by scaling by the distance moved. */
1373 integral *= this->GetRayPointSpacing();
1378 /* -----------------------------------------------------------------------
1379 ZeroState() - Set the default (zero) state of the object
1380 ----------------------------------------------------------------------- */
1382 template<class TInputImage, class TCoordRep>
1384 RayCastHelper<TInputImage, TCoordRep>
1391 m_NumberOfVoxelsInX = 0;
1392 m_NumberOfVoxelsInY = 0;
1393 m_NumberOfVoxelsInZ = 0;
1395 m_VoxelDimensionInX = 0;
1396 m_VoxelDimensionInY = 0;
1397 m_VoxelDimensionInZ = 0;
1399 for (i=0; i<3; i++) {
1400 m_CurrentRayPositionInMM[i] = 0.;
1402 for (i=0; i<3; i++) {
1403 m_RayDirectionInMM[i] = 0.;
1405 for (i=0; i<3; i++) {
1406 m_RayVoxelStartPosition[i] = 0.;
1408 for (i=0; i<3; i++) {
1409 m_RayVoxelEndPosition[i] = 0.;
1411 for (i=0; i<3; i++) {
1412 m_VoxelIncrement[i] = 0.;
1414 m_TraversalDirection = UNDEFINED_DIRECTION;
1416 m_TotalRayVoxelPlanes = 0;
1417 m_NumVoxelPlanesTraversed = -1;
1419 for (i=0; i<4; i++) {
1420 m_RayIntersectionVoxels[i] = 0;
1422 for (i=0; i<3; i++) {
1423 m_RayIntersectionVoxelIndex[i] = 0;
1426 }; // end of anonymous namespace
1432 /**************************************************************************
1435 * Rest of this code is the actual RayCastInterpolateImageFunctionWithOrigin
1439 **************************************************************************/
1441 /* -----------------------------------------------------------------------
1443 ----------------------------------------------------------------------- */
1445 template<class TInputImage, class TCoordRep>
1446 RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
1447 ::RayCastInterpolateImageFunctionWithOrigin()
1451 m_FocalPoint[0] = 0.;
1452 m_FocalPoint[1] = 0.;
1453 m_FocalPoint[2] = 0.;
1457 /* -----------------------------------------------------------------------
1459 ----------------------------------------------------------------------- */
1461 template<class TInputImage, class TCoordRep>
1463 RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
1464 ::PrintSelf(std::ostream& os, Indent indent) const
1466 this->Superclass::PrintSelf(os,indent);
1468 os << indent << "Threshold: " << m_Threshold << std::endl;
1469 os << indent << "FocalPoint: " << m_FocalPoint << std::endl;
1470 os << indent << "Transform: " << m_Transform.GetPointer() << std::endl;
1471 os << indent << "Interpolator: " << m_Interpolator.GetPointer() << std::endl;
1475 /* -----------------------------------------------------------------------
1476 Evaluate at image index position
1477 ----------------------------------------------------------------------- */
1479 template<class TInputImage, class TCoordRep>
1480 typename RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
1482 RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
1483 ::Evaluate( const PointType& point ) const
1485 double integral = 0;
1487 OutputPointType transformedFocalPoint
1488 = m_Transform->TransformPoint( m_FocalPoint );
1490 DirectionType direction = transformedFocalPoint - point;
1492 RayCastHelper<TInputImage, TCoordRep> ray;
1493 ray.SetImage( this->m_Image );
1497 // (Aviv) Added support for images with non-zero origin
1498 // ray.SetRay(point, direction);
1499 ray.SetRay(point - this->m_Image->GetOrigin().GetVectorFromOrigin(), direction);
1500 ray.IntegrateAboveThreshold(integral, m_Threshold);
1502 return ( static_cast<OutputType>( integral ));
1505 template<class TInputImage, class TCoordRep>
1506 typename RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
1508 RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
1509 ::EvaluateAtContinuousIndex( const ContinuousIndexType& index ) const
1511 OutputPointType point;
1512 this->m_Image->TransformContinuousIndexToPhysicalPoint(index, point);
1514 return this->Evaluate( point );