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://oncora1.lyon.fnclcc.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
29 /** \class Helper class to maintain state when casting a ray.
30 * This helper class keeps the RayCastInterpolateImageFunctionWithOrigin thread safe.
32 template <class TInputImage, class TCoordRep = float>
36 /** Constants for the image dimensions */
37 itkStaticConstMacro(InputImageDimension, unsigned int,
38 TInputImage::ImageDimension);
41 * Type of the Transform Base class
42 * The fixed image should be a 3D image
44 typedef itk::Transform<TCoordRep,3,3> TransformType;
46 typedef typename TransformType::Pointer TransformPointer;
47 typedef typename TransformType::InputPointType InputPointType;
48 typedef typename TransformType::OutputPointType OutputPointType;
49 typedef typename TransformType::ParametersType TransformParametersType;
50 typedef typename TransformType::JacobianType TransformJacobianType;
52 typedef typename TInputImage::SizeType SizeType;
53 typedef itk::Vector<TCoordRep, 3> DirectionType;
54 typedef itk::Point<TCoordRep, 3> PointType;
56 typedef TInputImage InputImageType;
57 typedef typename InputImageType::PixelType PixelType;
58 typedef typename InputImageType::IndexType IndexType;
63 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)
94 return IntegrateAboveThreshold(integral, 0);
99 * Integrate the interpolated intensities above a given threshold,
100 * along the ray and return the result.
102 * This routine can be called after instantiating the ray and
103 * calling SetProjectionCoord2D() or Reset(). It may then be called
104 * as many times thereafter for different 2D projection
107 * \param integral The integrated intensities along the ray.
108 * \param threshold The integration threshold [default value: 0]
110 * \return True if a valid ray was specified.
112 bool IntegrateAboveThreshold(double &integral, double threshold);
115 * Increment each of the intensities of the 4 planar voxels
116 * surrounding the current ray point.
118 * \parameter increment Intensity increment for each of the current 4 voxels
120 void IncrementIntensities(double increment=1);
122 /// Reset the iterator to the start of the ray.
125 /// Return the interpolated intensity of the current ray point.
126 double GetCurrentIntensity(void) const;
128 /// Return the ray point spacing in mm
129 double GetRayPointSpacing(void) const {
130 typename InputImageType::SpacingType spacing=this->m_Image->GetSpacing();
133 return vcl_sqrt(m_VoxelIncrement[0]*spacing[0]*m_VoxelIncrement[0]*spacing[0]
134 + m_VoxelIncrement[1]*spacing[1]*m_VoxelIncrement[1]*spacing[1]
135 + m_VoxelIncrement[2]*spacing[2]*m_VoxelIncrement[2]*spacing[2] );
140 /// Set the initial zero state of the object
143 /// Initialise the object
144 void Initialise(void);
147 /// Calculate the endpoint coordinats of the ray in voxels.
148 void EndPointsInVoxels(void);
151 * Calculate the incremental direction vector in voxels, 'dVoxel',
152 * required to traverse the ray.
154 void CalcDirnVector(void);
157 * Reduce the length of the ray until both start and end
158 * coordinates lie inside the volume.
160 * \return True if a valid ray has been, false otherwise.
162 bool AdjustRayLength(void);
165 * Obtain pointers to the four voxels surrounding the point where the ray
168 void InitialiseVoxelPointers(void);
170 /// Increment the voxel pointers surrounding the current point on the ray.
171 void IncrementVoxelPointers(void);
173 /// Record volume dimensions and resolution
174 void RecordVolumeDimensions(void);
176 /// Define the corners of the volume
177 void DefineCorners(void);
180 * Calculate the planes which define the volume.
182 * Member function to calculate the equations of the planes of 4 of
183 * the sides of the volume, calculate the positions of the 8 corners
184 * of the volume in mm in World, also calculate the values of the
185 * slopes of the lines which go to make up the volume( defined as
186 * lines in cube x,y,z dirn and then each of these lines has a slope
187 * in the world x,y,z dirn [3]) and finally also to return the length
188 * of the sides of the lines in mm.
190 void CalcPlanesAndCorners(void);
193 * Calculate the ray intercepts with the volume.
195 * See where the ray cuts the volume, check that truncation does not occur,
196 * if not, then start ray where it first intercepts the volume and set
197 * x_max to be where it leaves the volume.
199 * \return True if a valid ray has been specified, false otherwise.
201 bool CalcRayIntercepts(void);
204 * The ray is traversed by stepping in the axial direction
205 * that enables the greatest number of planes in the volume to be
209 UNDEFINED_DIRECTION=0, //!< Undefined
210 TRANSVERSE_IN_X, //!< x
211 TRANSVERSE_IN_Y, //!< y
212 TRANSVERSE_IN_Z, //!< z
214 } TraversalDirection;
216 // Cache the image in the structure. Skip the smart pointer for
217 // efficiency. This inner class will go in/out of scope with every
218 // call to Evaluate()
219 const InputImageType *m_Image;
221 /// Flag indicating whether the current ray is valid
225 * The start position of the ray in voxels.
227 * NB. Two of the components of this coordinate (i.e. those lying within
228 * the planes of voxels being traversed) will be shifted by half a
229 * voxel. This enables indices of the neighbouring voxels within the plane
230 * to be determined by simply casting to 'int' and optionally adding 1.
232 double m_RayVoxelStartPosition[3];
235 * The end coordinate of the ray in voxels.
237 * NB. Two of the components of this coordinate (i.e. those lying within
238 * the planes of voxels being traversed) will be shifted by half a
239 * voxel. This enables indices of the neighbouring voxels within the plane
240 * to be determined by simply casting to 'int' and optionally adding 1.
242 double m_RayVoxelEndPosition[3];
246 * The current coordinate on the ray in voxels.
248 * NB. Two of the components of this coordinate (i.e. those lying within
249 * the planes of voxels being traversed) will be shifted by half a
250 * voxel. This enables indices of the neighbouring voxels within the plane
251 * to be determined by simply casting to 'int' and optionally adding 1.
253 double m_Position3Dvox[3];
255 /** The incremental direction vector of the ray in voxels. */
256 double m_VoxelIncrement[3];
258 /// The direction in which the ray is incremented thorough the volume (x, y or z).
259 TraversalDirection m_TraversalDirection;
261 /// The total number of planes of voxels traversed by the ray.
262 int m_TotalRayVoxelPlanes;
264 /// The current number of planes of voxels traversed by the ray.
265 int m_NumVoxelPlanesTraversed;
267 /// Pointers to the current four voxels surrounding the ray's trajectory.
268 const PixelType *m_RayIntersectionVoxels[4];
271 * The voxel coordinate of the bottom-left voxel of the current
272 * four voxels surrounding the ray's trajectory.
274 int m_RayIntersectionVoxelIndex[3];
276 /// The dimension in voxels of the 3D volume in along the x axis
277 int m_NumberOfVoxelsInX;
278 /// The dimension in voxels of the 3D volume in along the y axis
279 int m_NumberOfVoxelsInY;
280 /// The dimension in voxels of the 3D volume in along the z axis
281 int m_NumberOfVoxelsInZ;
283 /// Voxel dimension in x
284 double m_VoxelDimensionInX;
285 /// Voxel dimension in y
286 double m_VoxelDimensionInY;
287 /// Voxel dimension in z
288 double m_VoxelDimensionInZ;
290 /// The coordinate of the point at which the ray enters the volume in mm.
291 double m_RayStartCoordInMM[3];
292 /// The coordinate of the point at which the ray exits the volume in mm.
293 double m_RayEndCoordInMM[3];
297 Planes which define the boundary of the volume in mm
298 (six planes and four parameters: Ax+By+Cz+D). */
299 double m_BoundingPlane[6][4];
300 /// The eight corners of the volume (x,y,z coordinates for each).
301 double m_BoundingCorner[8][3];
303 /// The position of the ray
304 double m_CurrentRayPositionInMM[3];
306 /// The direction of the ray
307 double m_RayDirectionInMM[3];
311 /* -----------------------------------------------------------------------
312 Initialise() - Initialise the object
313 ----------------------------------------------------------------------- */
315 template<class TInputImage, class TCoordRep>
317 RayCastHelper<TInputImage, TCoordRep>
320 // Save the dimensions of the volume and calculate the bounding box
321 this->RecordVolumeDimensions();
323 // Calculate the planes and corners which define the volume.
324 this->DefineCorners();
325 this->CalcPlanesAndCorners();
329 /* -----------------------------------------------------------------------
330 RecordVolumeDimensions() - Record volume dimensions and resolution
331 ----------------------------------------------------------------------- */
333 template<class TInputImage, class TCoordRep>
335 RayCastHelper<TInputImage, TCoordRep>
336 ::RecordVolumeDimensions(void)
338 typename InputImageType::SpacingType spacing=this->m_Image->GetSpacing();
339 SizeType dim=this->m_Image->GetLargestPossibleRegion().GetSize();
341 m_NumberOfVoxelsInX = dim[0];
342 m_NumberOfVoxelsInY = dim[1];
343 m_NumberOfVoxelsInZ = dim[2];
345 m_VoxelDimensionInX = spacing[0];
346 m_VoxelDimensionInY = spacing[1];
347 m_VoxelDimensionInZ = spacing[2];
351 /* -----------------------------------------------------------------------
352 DefineCorners() - Define the corners of the volume
353 ----------------------------------------------------------------------- */
355 template<class TInputImage, class TCoordRep>
357 RayCastHelper<TInputImage, TCoordRep>
358 ::DefineCorners(void)
360 // Define corner positions as if at the origin
362 m_BoundingCorner[0][0] =
363 m_BoundingCorner[1][0] =
364 m_BoundingCorner[2][0] =
365 m_BoundingCorner[3][0] = 0;
367 m_BoundingCorner[4][0] =
368 m_BoundingCorner[5][0] =
369 m_BoundingCorner[6][0] =
370 m_BoundingCorner[7][0] = m_VoxelDimensionInX*m_NumberOfVoxelsInX;
372 m_BoundingCorner[1][1] =
373 m_BoundingCorner[3][1] =
374 m_BoundingCorner[5][1] =
375 m_BoundingCorner[7][1] = m_VoxelDimensionInY*m_NumberOfVoxelsInY;
377 m_BoundingCorner[0][1] =
378 m_BoundingCorner[2][1] =
379 m_BoundingCorner[4][1] =
380 m_BoundingCorner[6][1] = 0;
382 m_BoundingCorner[0][2] =
383 m_BoundingCorner[1][2] =
384 m_BoundingCorner[4][2] =
385 m_BoundingCorner[5][2] =
386 m_VoxelDimensionInZ*m_NumberOfVoxelsInZ;
388 m_BoundingCorner[2][2] =
389 m_BoundingCorner[3][2] =
390 m_BoundingCorner[6][2] =
391 m_BoundingCorner[7][2] = 0;
395 /* -----------------------------------------------------------------------
396 CalcPlanesAndCorners() - Calculate the planes and corners of the volume.
397 ----------------------------------------------------------------------- */
399 template<class TInputImage, class TCoordRep>
401 RayCastHelper<TInputImage, TCoordRep>
402 ::CalcPlanesAndCorners(void)
407 // find the equations of the planes
409 int c1=0, c2=0, c3=0;
412 { // loop around for planes
414 { // which corners to take
436 double line1x, line1y, line1z;
437 double line2x, line2y, line2z;
439 // lines from one corner to another in x,y,z dirns
440 line1x = m_BoundingCorner[c1][0] - m_BoundingCorner[c2][0];
441 line2x = m_BoundingCorner[c1][0] - m_BoundingCorner[c3][0];
443 line1y = m_BoundingCorner[c1][1] - m_BoundingCorner[c2][1];
444 line2y = m_BoundingCorner[c1][1] - m_BoundingCorner[c3][1];
446 line1z = m_BoundingCorner[c1][2] - m_BoundingCorner[c2][2];
447 line2z = m_BoundingCorner[c1][2] - m_BoundingCorner[c3][2];
451 // take cross product
452 A = line1y*line2z - line2y*line1z;
453 B = line2x*line1z - line1x*line2z;
454 C = line1x*line2y - line2x*line1y;
457 D = -( A*m_BoundingCorner[c1][0]
458 + B*m_BoundingCorner[c1][1]
459 + C*m_BoundingCorner[c1][2] );
461 // initialise plane value and normalise
462 m_BoundingPlane[j][0] = A/vcl_sqrt(A*A + B*B + C*C);
463 m_BoundingPlane[j][1] = B/vcl_sqrt(A*A + B*B + C*C);
464 m_BoundingPlane[j][2] = C/vcl_sqrt(A*A + B*B + C*C);
465 m_BoundingPlane[j][3] = D/vcl_sqrt(A*A + B*B + C*C);
467 if ( (A*A + B*B + C*C) == 0 )
469 itk::ExceptionObject err(__FILE__, __LINE__);
470 err.SetLocation( ITK_LOCATION );
471 err.SetDescription( "Division by zero (planes) "
472 "- CalcPlanesAndCorners().");
480 /* -----------------------------------------------------------------------
481 CalcRayIntercepts() - Calculate the ray intercepts with the volume.
482 ----------------------------------------------------------------------- */
484 template<class TInputImage, class TCoordRep>
486 RayCastHelper<TInputImage, TCoordRep>
487 ::CalcRayIntercepts()
489 double maxInterDist, interDist;
490 double cornerVect[4][3];
491 int cross[4][3], noInterFlag[6];
492 int nSidesCrossed, crossFlag, c[4];
493 double ax, ay, az, bx, by, bz;
494 double cubeInter[6][3];
498 int NoSides = 6; // =6 to allow truncation: =4 to remove truncated rays
500 // Calculate intercept of ray with planes
502 double interceptx[6], intercepty[6], interceptz[6];
505 for( j=0; j<NoSides; j++)
508 denom = ( m_BoundingPlane[j][0]*m_RayDirectionInMM[0]
509 + m_BoundingPlane[j][1]*m_RayDirectionInMM[1]
510 + m_BoundingPlane[j][2]*m_RayDirectionInMM[2]);
512 if( (long)(denom*100) != 0 )
514 d[j] = -( m_BoundingPlane[j][3]
515 + m_BoundingPlane[j][0]*m_CurrentRayPositionInMM[0]
516 + m_BoundingPlane[j][1]*m_CurrentRayPositionInMM[1]
517 + m_BoundingPlane[j][2]*m_CurrentRayPositionInMM[2] ) / denom;
519 interceptx[j] = m_CurrentRayPositionInMM[0] + d[j]*m_RayDirectionInMM[0];
520 intercepty[j] = m_CurrentRayPositionInMM[1] + d[j]*m_RayDirectionInMM[1];
521 interceptz[j] = m_CurrentRayPositionInMM[2] + d[j]*m_RayDirectionInMM[2];
523 noInterFlag[j] = 1; //OK
527 noInterFlag[j] = 0; //NOT OK
533 for( j=0; j<NoSides; j++ )
536 // Work out which corners to use
540 c[0] = 0; c[1] = 1; c[2] = 3; c[3] = 2;
544 c[0] = 4; c[1] = 5; c[2] = 7; c[3] = 6;
548 c[0] = 1; c[1] = 5; c[2] = 7; c[3] = 3;
552 c[0] = 0; c[1] = 2; c[2] = 6; c[3] = 4;
556 c[0] = 0; c[1] = 1; c[2] = 5; c[3] = 4;
560 c[0] = 2; c[1] = 3; c[2] = 7; c[3] = 6;
563 // Calculate vectors from corner of ct volume to intercept.
566 if( noInterFlag[j]==1 )
568 cornerVect[i][0] = m_BoundingCorner[c[i]][0] - interceptx[j];
569 cornerVect[i][1] = m_BoundingCorner[c[i]][1] - intercepty[j];
570 cornerVect[i][2] = m_BoundingCorner[c[i]][2] - interceptz[j];
572 else if( noInterFlag[j]==0 )
574 cornerVect[i][0] = 0;
575 cornerVect[i][1] = 0;
576 cornerVect[i][2] = 0;
581 // Do cross product with these vectors
592 ax = cornerVect[i][0];
593 ay = cornerVect[i][1];
594 az = cornerVect[i][2];
595 bx = cornerVect[k][0];
596 by = cornerVect[k][1];
597 bz = cornerVect[k][2];
599 // The int and divide by 100 are to avoid rounding errors. If
600 // these are not included then you get values fluctuating around
601 // zero and so in the subsequent check, all the values are not
602 // above or below zero. NB. If you "INT" by too much here though
603 // you can get problems in the corners of your volume when rays
604 // are allowed to go through more than one plane.
605 cross[i][0] = (int)((ay*bz - az*by)/100);
606 cross[i][1] = (int)((az*bx - ax*bz)/100);
607 cross[i][2] = (int)((ax*by - ay*bx)/100);
610 // See if a sign change occured between all these cross products
611 // if not, then the ray went through this plane
631 if( crossFlag==3 && noInterFlag[j]==1 )
633 cubeInter[nSidesCrossed][0] = interceptx[j];
634 cubeInter[nSidesCrossed][1] = intercepty[j];
635 cubeInter[nSidesCrossed][2] = interceptz[j];
639 } // End of loop over all four planes
641 m_RayStartCoordInMM[0] = cubeInter[0][0];
642 m_RayStartCoordInMM[1] = cubeInter[0][1];
643 m_RayStartCoordInMM[2] = cubeInter[0][2];
645 m_RayEndCoordInMM[0] = cubeInter[1][0];
646 m_RayEndCoordInMM[1] = cubeInter[1][1];
647 m_RayEndCoordInMM[2] = cubeInter[1][2];
649 if( nSidesCrossed >= 5 )
651 std::cerr << "WARNING: No. of sides crossed equals: " << nSidesCrossed << std::endl;
654 // If 'nSidesCrossed' is larger than 2, this means that the ray goes through
655 // a corner of the volume and due to rounding errors, the ray is
656 // deemed to go through more than two planes. To obtain the correct
657 // start and end positions we choose the two intercept values which
658 // are furthest from each other.
661 if( nSidesCrossed >= 3 )
664 for( j=0; j<nSidesCrossed-1; j++ )
666 for( k=j+1; k<nSidesCrossed; k++ )
671 interDist += (cubeInter[j][i] - cubeInter[k][i])*
672 (cubeInter[j][i] - cubeInter[k][i]);
674 if( interDist > maxInterDist )
676 maxInterDist = interDist;
678 m_RayStartCoordInMM[0] = cubeInter[j][0];
679 m_RayStartCoordInMM[1] = cubeInter[j][1];
680 m_RayStartCoordInMM[2] = cubeInter[j][2];
682 m_RayEndCoordInMM[0] = cubeInter[k][0];
683 m_RayEndCoordInMM[1] = cubeInter[k][1];
684 m_RayEndCoordInMM[2] = cubeInter[k][2];
691 if (nSidesCrossed == 2 )
702 /* -----------------------------------------------------------------------
703 SetRay() - Set the position and direction of the ray
704 ----------------------------------------------------------------------- */
706 template<class TInputImage, class TCoordRep>
708 RayCastHelper<TInputImage, TCoordRep>
709 ::SetRay(OutputPointType RayPosn, DirectionType RayDirn)
712 // Store the position and direction of the ray
713 typename TInputImage::SpacingType spacing=this->m_Image->GetSpacing();
714 SizeType dim=this->m_Image->GetLargestPossibleRegion().GetSize();
716 // we need to translate the _center_ of the volume to the origin
717 m_NumberOfVoxelsInX = dim[0];
718 m_NumberOfVoxelsInY = dim[1];
719 m_NumberOfVoxelsInZ = dim[2];
721 m_VoxelDimensionInX = spacing[0];
722 m_VoxelDimensionInY = spacing[1];
723 m_VoxelDimensionInZ = spacing[2];
725 // (Aviv) Incorporating a fix by Jian Wu
726 // http://public.kitware.com/pipermail/insight-users/2006-March/017265.html
727 m_CurrentRayPositionInMM[0] = RayPosn[0];
728 m_CurrentRayPositionInMM[1] = RayPosn[1];
729 m_CurrentRayPositionInMM[2] = RayPosn[2];
731 m_RayDirectionInMM[0] = RayDirn[0];
732 m_RayDirectionInMM[1] = RayDirn[1];
733 m_RayDirectionInMM[2] = RayDirn[2];
735 // Compute the ray path for this coordinate in mm
737 m_ValidRay = this->CalcRayIntercepts();
745 // Convert the start and end coordinates of the ray to voxels
747 this->EndPointsInVoxels();
749 /* Calculate the ray direction vector in voxels and hence the voxel
750 increment required to traverse the ray, and the number of
751 interpolation points on the ray.
753 This routine also shifts the coordinate frame by half a voxel for
754 two of the directional components (i.e. those lying within the
755 planes of voxels being traversed). */
757 this->CalcDirnVector();
760 /* Reduce the length of the ray until both start and end
761 coordinates lie inside the volume. */
763 m_ValidRay = this->AdjustRayLength();
765 // Reset the iterator to the start of the ray.
773 /* -----------------------------------------------------------------------
774 EndPointsInVoxels() - Convert the endpoints to voxels
775 ----------------------------------------------------------------------- */
777 template<class TInputImage, class TCoordRep>
779 RayCastHelper<TInputImage, TCoordRep>
780 ::EndPointsInVoxels(void)
782 m_RayVoxelStartPosition[0] = m_RayStartCoordInMM[0]/m_VoxelDimensionInX;
783 m_RayVoxelStartPosition[1] = m_RayStartCoordInMM[1]/m_VoxelDimensionInY;
784 m_RayVoxelStartPosition[2] = m_RayStartCoordInMM[2]/m_VoxelDimensionInZ;
786 m_RayVoxelEndPosition[0] = m_RayEndCoordInMM[0]/m_VoxelDimensionInX;
787 m_RayVoxelEndPosition[1] = m_RayEndCoordInMM[1]/m_VoxelDimensionInY;
788 m_RayVoxelEndPosition[2] = m_RayEndCoordInMM[2]/m_VoxelDimensionInZ;
793 /* -----------------------------------------------------------------------
794 CalcDirnVector() - Calculate the incremental direction vector in voxels.
795 ----------------------------------------------------------------------- */
797 template<class TInputImage, class TCoordRep>
799 RayCastHelper<TInputImage, TCoordRep>
800 ::CalcDirnVector(void)
802 double xNum, yNum, zNum;
804 // Calculate the number of voxels in each direction
806 xNum = vcl_fabs(m_RayVoxelStartPosition[0] - m_RayVoxelEndPosition[0]);
807 yNum = vcl_fabs(m_RayVoxelStartPosition[1] - m_RayVoxelEndPosition[1]);
808 zNum = vcl_fabs(m_RayVoxelStartPosition[2] - m_RayVoxelEndPosition[2]);
810 // The direction iterated in is that with the greatest number of voxels
811 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
813 // Iterate in X direction
815 if( (xNum >= yNum) && (xNum >= zNum) )
817 if( m_RayVoxelStartPosition[0] < m_RayVoxelEndPosition[0] )
819 m_VoxelIncrement[0] = 1;
822 = (m_RayVoxelStartPosition[1]
823 - m_RayVoxelEndPosition[1])/(m_RayVoxelStartPosition[0]
824 - m_RayVoxelEndPosition[0]);
827 = (m_RayVoxelStartPosition[2]
828 - m_RayVoxelEndPosition[2])/(m_RayVoxelStartPosition[0]
829 - m_RayVoxelEndPosition[0]);
833 m_VoxelIncrement[0] = -1;
836 = -(m_RayVoxelStartPosition[1]
837 - m_RayVoxelEndPosition[1])/(m_RayVoxelStartPosition[0]
838 - m_RayVoxelEndPosition[0]);
841 = -(m_RayVoxelStartPosition[2]
842 - m_RayVoxelEndPosition[2])/(m_RayVoxelStartPosition[0]
843 - m_RayVoxelEndPosition[0]);
846 // This section is to alter the start position in order to
847 // place the center of the voxels in there correct positions,
848 // rather than placing them at the corner of voxels which is
849 // what happens if this is not carried out. The reason why
850 // x has no -0.5 is because this is the direction we are going
851 // to iterate in and therefore we wish to go from center to
852 // center rather than finding the surrounding voxels.
854 m_RayVoxelStartPosition[1] += ( (int)m_RayVoxelStartPosition[0]
855 - m_RayVoxelStartPosition[0])*m_VoxelIncrement[1]*m_VoxelIncrement[0]
856 + 0.5*m_VoxelIncrement[1] - 0.5;
858 m_RayVoxelStartPosition[2] += ( (int)m_RayVoxelStartPosition[0]
859 - m_RayVoxelStartPosition[0])*m_VoxelIncrement[2]*m_VoxelIncrement[0]
860 + 0.5*m_VoxelIncrement[2] - 0.5;
862 m_RayVoxelStartPosition[0] = (int)m_RayVoxelStartPosition[0] + 0.5*m_VoxelIncrement[0];
864 m_TotalRayVoxelPlanes = (int)xNum;
866 m_TraversalDirection = TRANSVERSE_IN_X;
869 // Iterate in Y direction
871 else if( (yNum >= xNum) && (yNum >= zNum) )
874 if( m_RayVoxelStartPosition[1] < m_RayVoxelEndPosition[1] )
876 m_VoxelIncrement[1] = 1;
879 = (m_RayVoxelStartPosition[0]
880 - m_RayVoxelEndPosition[0])/(m_RayVoxelStartPosition[1]
881 - m_RayVoxelEndPosition[1]);
884 = (m_RayVoxelStartPosition[2]
885 - m_RayVoxelEndPosition[2])/(m_RayVoxelStartPosition[1]
886 - m_RayVoxelEndPosition[1]);
890 m_VoxelIncrement[1] = -1;
893 = -(m_RayVoxelStartPosition[0]
894 - m_RayVoxelEndPosition[0])/(m_RayVoxelStartPosition[1]
895 - m_RayVoxelEndPosition[1]);
898 = -(m_RayVoxelStartPosition[2]
899 - m_RayVoxelEndPosition[2])/(m_RayVoxelStartPosition[1]
900 - m_RayVoxelEndPosition[1]);
903 m_RayVoxelStartPosition[0] += ( (int)m_RayVoxelStartPosition[1]
904 - m_RayVoxelStartPosition[1])*m_VoxelIncrement[0]*m_VoxelIncrement[1]
905 + 0.5*m_VoxelIncrement[0] - 0.5;
907 m_RayVoxelStartPosition[2] += ( (int)m_RayVoxelStartPosition[1]
908 - m_RayVoxelStartPosition[1])*m_VoxelIncrement[2]*m_VoxelIncrement[1]
909 + 0.5*m_VoxelIncrement[2] - 0.5;
911 m_RayVoxelStartPosition[1] = (int)m_RayVoxelStartPosition[1] + 0.5*m_VoxelIncrement[1];
913 m_TotalRayVoxelPlanes = (int)yNum;
915 m_TraversalDirection = TRANSVERSE_IN_Y;
918 // Iterate in Z direction
923 if( m_RayVoxelStartPosition[2] < m_RayVoxelEndPosition[2] )
925 m_VoxelIncrement[2] = 1;
928 = (m_RayVoxelStartPosition[0]
929 - m_RayVoxelEndPosition[0])/(m_RayVoxelStartPosition[2]
930 - m_RayVoxelEndPosition[2]);
933 = (m_RayVoxelStartPosition[1]
934 - m_RayVoxelEndPosition[1])/(m_RayVoxelStartPosition[2]
935 - m_RayVoxelEndPosition[2]);
939 m_VoxelIncrement[2] = -1;
942 = -(m_RayVoxelStartPosition[0]
943 - m_RayVoxelEndPosition[0])/(m_RayVoxelStartPosition[2]
944 - m_RayVoxelEndPosition[2]);
947 = -(m_RayVoxelStartPosition[1]
948 - m_RayVoxelEndPosition[1])/(m_RayVoxelStartPosition[2]
949 - m_RayVoxelEndPosition[2]);
952 m_RayVoxelStartPosition[0] += ( (int)m_RayVoxelStartPosition[2]
953 - m_RayVoxelStartPosition[2])*m_VoxelIncrement[0]*m_VoxelIncrement[2]
954 + 0.5*m_VoxelIncrement[0] - 0.5;
956 m_RayVoxelStartPosition[1] += ( (int)m_RayVoxelStartPosition[2]
957 - m_RayVoxelStartPosition[2])*m_VoxelIncrement[1]*m_VoxelIncrement[2]
958 + 0.5*m_VoxelIncrement[1] - 0.5;
960 m_RayVoxelStartPosition[2] = (int)m_RayVoxelStartPosition[2] + 0.5*m_VoxelIncrement[2];
962 m_TotalRayVoxelPlanes = (int)zNum;
964 m_TraversalDirection = TRANSVERSE_IN_Z;
969 /* -----------------------------------------------------------------------
970 AdjustRayLength() - Ensure that the ray lies within the volume
971 ----------------------------------------------------------------------- */
973 template<class TInputImage, class TCoordRep>
975 RayCastHelper<TInputImage, TCoordRep>
976 ::AdjustRayLength(void)
983 if (m_TraversalDirection == TRANSVERSE_IN_X)
989 else if (m_TraversalDirection == TRANSVERSE_IN_Y)
995 else if (m_TraversalDirection == TRANSVERSE_IN_Z)
1003 itk::ExceptionObject err(__FILE__, __LINE__);
1004 err.SetLocation( ITK_LOCATION );
1005 err.SetDescription( "The ray traversal direction is unset "
1006 "- AdjustRayLength().");
1018 Istart[0] = (int) vcl_floor(m_RayVoxelStartPosition[0]);
1019 Istart[1] = (int) vcl_floor(m_RayVoxelStartPosition[1]);
1020 Istart[2] = (int) vcl_floor(m_RayVoxelStartPosition[2]);
1022 if( (Istart[0] >= 0) && (Istart[0] + Idirn[0] < m_NumberOfVoxelsInX) &&
1023 (Istart[1] >= 0) && (Istart[1] + Idirn[1] < m_NumberOfVoxelsInY) &&
1024 (Istart[2] >= 0) && (Istart[2] + Idirn[2] < m_NumberOfVoxelsInZ) )
1031 m_RayVoxelStartPosition[0] += m_VoxelIncrement[0];
1032 m_RayVoxelStartPosition[1] += m_VoxelIncrement[1];
1033 m_RayVoxelStartPosition[2] += m_VoxelIncrement[2];
1035 m_TotalRayVoxelPlanes--;
1038 Istart[0] = (int) vcl_floor(m_RayVoxelStartPosition[0]
1039 + m_TotalRayVoxelPlanes*m_VoxelIncrement[0]);
1041 Istart[1] = (int) vcl_floor(m_RayVoxelStartPosition[1]
1042 + m_TotalRayVoxelPlanes*m_VoxelIncrement[1]);
1044 Istart[2] = (int) vcl_floor(m_RayVoxelStartPosition[2]
1045 + m_TotalRayVoxelPlanes*m_VoxelIncrement[2]);
1047 if( (Istart[0] >= 0) && (Istart[0] + Idirn[0] < m_NumberOfVoxelsInX) &&
1048 (Istart[1] >= 0) && (Istart[1] + Idirn[1] < m_NumberOfVoxelsInY) &&
1049 (Istart[2] >= 0) && (Istart[2] + Idirn[2] < m_NumberOfVoxelsInZ) )
1056 m_TotalRayVoxelPlanes--;
1059 } while ( (! (startOK && endOK)) && (m_TotalRayVoxelPlanes > 1) );
1062 return (startOK && endOK);
1066 /* -----------------------------------------------------------------------
1067 Reset() - Reset the iterator to the start of the ray.
1068 ----------------------------------------------------------------------- */
1070 template<class TInputImage, class TCoordRep>
1072 RayCastHelper<TInputImage, TCoordRep>
1077 m_NumVoxelPlanesTraversed = -1;
1079 // If this is a valid ray...
1085 m_Position3Dvox[i] = m_RayVoxelStartPosition[i];
1087 this->InitialiseVoxelPointers();
1090 // otherwise set parameters to zero
1096 m_RayVoxelStartPosition[i] = 0.;
1100 m_RayVoxelEndPosition[i] = 0.;
1104 m_VoxelIncrement[i] = 0.;
1106 m_TraversalDirection = UNDEFINED_DIRECTION;
1108 m_TotalRayVoxelPlanes = 0;
1112 m_RayIntersectionVoxels[i] = 0;
1116 m_RayIntersectionVoxelIndex[i] = 0;
1122 /* -----------------------------------------------------------------------
1123 InitialiseVoxelPointers() - Obtain pointers to the first four voxels
1124 ----------------------------------------------------------------------- */
1126 template<class TInputImage, class TCoordRep>
1128 RayCastHelper<TInputImage, TCoordRep>
1129 ::InitialiseVoxelPointers(void)
1135 Ix = (int)(m_RayVoxelStartPosition[0]);
1136 Iy = (int)(m_RayVoxelStartPosition[1]);
1137 Iz = (int)(m_RayVoxelStartPosition[2]);
1139 m_RayIntersectionVoxelIndex[0] = Ix;
1140 m_RayIntersectionVoxelIndex[1] = Iy;
1141 m_RayIntersectionVoxelIndex[2] = Iz;
1143 switch( m_TraversalDirection )
1145 case TRANSVERSE_IN_X:
1148 if( (Ix >= 0) && (Ix < m_NumberOfVoxelsInX) &&
1149 (Iy >= 0) && (Iy + 1 < m_NumberOfVoxelsInY) &&
1150 (Iz >= 0) && (Iz + 1 < m_NumberOfVoxelsInZ))
1153 index[0]=Ix; index[1]=Iy; index[2]=Iz;
1154 m_RayIntersectionVoxels[0]
1155 = this->m_Image->GetBufferPointer() + this->m_Image->ComputeOffset(index);
1157 index[0]=Ix; index[1]=Iy+1; index[2]=Iz;
1158 m_RayIntersectionVoxels[1]
1159 = ( this->m_Image->GetBufferPointer() + this->m_Image->ComputeOffset(index) );
1161 index[0]=Ix; index[1]=Iy; index[2]=Iz+1;
1162 m_RayIntersectionVoxels[2]
1163 = ( this->m_Image->GetBufferPointer() + this->m_Image->ComputeOffset(index) );
1165 index[0]=Ix; index[1]=Iy+1; index[2]=Iz+1;
1166 m_RayIntersectionVoxels[3]
1167 = ( this->m_Image->GetBufferPointer() + this->m_Image->ComputeOffset(index) );
1171 m_RayIntersectionVoxels[0] =
1172 m_RayIntersectionVoxels[1] =
1173 m_RayIntersectionVoxels[2] =
1174 m_RayIntersectionVoxels[3] = NULL;
1179 case TRANSVERSE_IN_Y:
1182 if( (Ix >= 0) && (Ix + 1 < m_NumberOfVoxelsInX) &&
1183 (Iy >= 0) && (Iy < m_NumberOfVoxelsInY) &&
1184 (Iz >= 0) && (Iz + 1 < m_NumberOfVoxelsInZ))
1187 index[0]=Ix; index[1]=Iy; index[2]=Iz;
1188 m_RayIntersectionVoxels[0] = ( this->m_Image->GetBufferPointer()
1189 + this->m_Image->ComputeOffset(index) );
1191 index[0]=Ix+1; index[1]=Iy; index[2]=Iz;
1192 m_RayIntersectionVoxels[1] = ( this->m_Image->GetBufferPointer()
1193 + this->m_Image->ComputeOffset(index) );
1195 index[0]=Ix; index[1]=Iy; index[2]=Iz+1;
1196 m_RayIntersectionVoxels[2] = ( this->m_Image->GetBufferPointer()
1197 + this->m_Image->ComputeOffset(index) );
1199 index[0]=Ix+1; index[1]=Iy; index[2]=Iz+1;
1200 m_RayIntersectionVoxels[3] = ( this->m_Image->GetBufferPointer()
1201 + this->m_Image->ComputeOffset(index) );
1205 m_RayIntersectionVoxels[0]
1206 = m_RayIntersectionVoxels[1]
1207 = m_RayIntersectionVoxels[2]
1208 = m_RayIntersectionVoxels[3] = NULL;
1213 case TRANSVERSE_IN_Z:
1216 if( (Ix >= 0) && (Ix + 1 < m_NumberOfVoxelsInX) &&
1217 (Iy >= 0) && (Iy + 1 < m_NumberOfVoxelsInY) &&
1218 (Iz >= 0) && (Iz < m_NumberOfVoxelsInZ))
1221 index[0]=Ix; index[1]=Iy; index[2]=Iz;
1222 m_RayIntersectionVoxels[0] = ( this->m_Image->GetBufferPointer()
1223 + this->m_Image->ComputeOffset(index) );
1225 index[0]=Ix+1; index[1]=Iy; index[2]=Iz;
1226 m_RayIntersectionVoxels[1] = ( this->m_Image->GetBufferPointer()
1227 + this->m_Image->ComputeOffset(index) );
1229 index[0]=Ix; index[1]=Iy+1; index[2]=Iz;
1230 m_RayIntersectionVoxels[2] = ( this->m_Image->GetBufferPointer()
1231 + this->m_Image->ComputeOffset(index) );
1233 index[0]=Ix+1; index[1]=Iy+1; index[2]=Iz;
1234 m_RayIntersectionVoxels[3] = ( this->m_Image->GetBufferPointer()
1235 + this->m_Image->ComputeOffset(index) );
1240 m_RayIntersectionVoxels[0]
1241 = m_RayIntersectionVoxels[1]
1242 = m_RayIntersectionVoxels[2]
1243 = m_RayIntersectionVoxels[3] = NULL;
1250 itk::ExceptionObject err(__FILE__, __LINE__);
1251 err.SetLocation( ITK_LOCATION );
1252 err.SetDescription( "The ray traversal direction is unset "
1253 "- InitialiseVoxelPointers().");
1260 /* -----------------------------------------------------------------------
1261 IncrementVoxelPointers() - Increment the voxel pointers
1262 ----------------------------------------------------------------------- */
1264 template<class TInputImage, class TCoordRep>
1266 RayCastHelper<TInputImage, TCoordRep>
1267 ::IncrementVoxelPointers(void)
1269 double xBefore = m_Position3Dvox[0];
1270 double yBefore = m_Position3Dvox[1];
1271 double zBefore = m_Position3Dvox[2];
1273 m_Position3Dvox[0] += m_VoxelIncrement[0];
1274 m_Position3Dvox[1] += m_VoxelIncrement[1];
1275 m_Position3Dvox[2] += m_VoxelIncrement[2];
1277 int dx = ((int) m_Position3Dvox[0]) - ((int) xBefore);
1278 int dy = ((int) m_Position3Dvox[1]) - ((int) yBefore);
1279 int dz = ((int) m_Position3Dvox[2]) - ((int) zBefore);
1281 m_RayIntersectionVoxelIndex[0] += dx;
1282 m_RayIntersectionVoxelIndex[1] += dy;
1283 m_RayIntersectionVoxelIndex[2] += dz;
1285 int totalRayVoxelPlanes
1286 = dx + dy*m_NumberOfVoxelsInX + dz*m_NumberOfVoxelsInX*m_NumberOfVoxelsInY;
1288 m_RayIntersectionVoxels[0] += totalRayVoxelPlanes;
1289 m_RayIntersectionVoxels[1] += totalRayVoxelPlanes;
1290 m_RayIntersectionVoxels[2] += totalRayVoxelPlanes;
1291 m_RayIntersectionVoxels[3] += totalRayVoxelPlanes;
1295 /* -----------------------------------------------------------------------
1296 GetCurrentIntensity() - Get the intensity of the current ray point.
1297 ----------------------------------------------------------------------- */
1299 template<class TInputImage, class TCoordRep>
1301 RayCastHelper<TInputImage, TCoordRep>
1302 ::GetCurrentIntensity(void) const
1311 a = (double) (*m_RayIntersectionVoxels[0]);
1312 b = (double) (*m_RayIntersectionVoxels[1] - a);
1313 c = (double) (*m_RayIntersectionVoxels[2] - a);
1314 d = (double) (*m_RayIntersectionVoxels[3] - a - b - c);
1316 switch( m_TraversalDirection )
1318 case TRANSVERSE_IN_X:
1320 y = m_Position3Dvox[1] - vcl_floor(m_Position3Dvox[1]);
1321 z = m_Position3Dvox[2] - vcl_floor(m_Position3Dvox[2]);
1324 case TRANSVERSE_IN_Y:
1326 y = m_Position3Dvox[0] - vcl_floor(m_Position3Dvox[0]);
1327 z = m_Position3Dvox[2] - vcl_floor(m_Position3Dvox[2]);
1330 case TRANSVERSE_IN_Z:
1332 y = m_Position3Dvox[0] - vcl_floor(m_Position3Dvox[0]);
1333 z = m_Position3Dvox[1] - vcl_floor(m_Position3Dvox[1]);
1338 itk::ExceptionObject err(__FILE__, __LINE__);
1339 err.SetLocation( ITK_LOCATION );
1340 err.SetDescription( "The ray traversal direction is unset "
1341 "- GetCurrentIntensity().");
1347 return a + b*y + c*z + d*y*z;
1350 /* -----------------------------------------------------------------------
1351 IncrementIntensities() - Increment the intensities of the current ray point
1352 ----------------------------------------------------------------------- */
1354 template<class TInputImage, class TCoordRep>
1356 RayCastHelper<TInputImage, TCoordRep>
1357 ::IncrementIntensities(double increment)
1359 short inc = (short) vcl_floor(increment + 0.5);
1365 *m_RayIntersectionVoxels[0] += inc;
1366 *m_RayIntersectionVoxels[1] += inc;
1367 *m_RayIntersectionVoxels[2] += inc;
1368 *m_RayIntersectionVoxels[3] += inc;
1374 /* -----------------------------------------------------------------------
1375 IntegrateAboveThreshold() - Integrate intensities above a threshold.
1376 ----------------------------------------------------------------------- */
1378 template<class TInputImage, class TCoordRep>
1380 RayCastHelper<TInputImage, TCoordRep>
1381 ::IntegrateAboveThreshold(double &integral, double threshold)
1384 // double posn3D_x, posn3D_y, posn3D_z;
1388 // Check if this is a valid ray
1394 /* Step along the ray as quickly as possible
1395 integrating the interpolated intensities. */
1397 for (m_NumVoxelPlanesTraversed=0;
1398 m_NumVoxelPlanesTraversed<m_TotalRayVoxelPlanes;
1399 m_NumVoxelPlanesTraversed++)
1401 intensity = this->GetCurrentIntensity();
1403 if (intensity > threshold)
1405 integral += intensity - threshold;
1407 this->IncrementVoxelPointers();
1410 /* The ray passes through the volume one plane of voxels at a time,
1411 however, if its moving diagonally the ray points will be further
1412 apart so account for this by scaling by the distance moved. */
1414 integral *= this->GetRayPointSpacing();
1419 /* -----------------------------------------------------------------------
1420 ZeroState() - Set the default (zero) state of the object
1421 ----------------------------------------------------------------------- */
1423 template<class TInputImage, class TCoordRep>
1425 RayCastHelper<TInputImage, TCoordRep>
1432 m_NumberOfVoxelsInX = 0;
1433 m_NumberOfVoxelsInY = 0;
1434 m_NumberOfVoxelsInZ = 0;
1436 m_VoxelDimensionInX = 0;
1437 m_VoxelDimensionInY = 0;
1438 m_VoxelDimensionInZ = 0;
1442 m_CurrentRayPositionInMM[i] = 0.;
1446 m_RayDirectionInMM[i] = 0.;
1450 m_RayVoxelStartPosition[i] = 0.;
1454 m_RayVoxelEndPosition[i] = 0.;
1458 m_VoxelIncrement[i] = 0.;
1460 m_TraversalDirection = UNDEFINED_DIRECTION;
1462 m_TotalRayVoxelPlanes = 0;
1463 m_NumVoxelPlanesTraversed = -1;
1467 m_RayIntersectionVoxels[i] = 0;
1471 m_RayIntersectionVoxelIndex[i] = 0;
1474 }; // end of anonymous namespace
1480 /**************************************************************************
1483 * Rest of this code is the actual RayCastInterpolateImageFunctionWithOrigin
1487 **************************************************************************/
1489 /* -----------------------------------------------------------------------
1491 ----------------------------------------------------------------------- */
1493 template<class TInputImage, class TCoordRep>
1494 RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
1495 ::RayCastInterpolateImageFunctionWithOrigin()
1499 m_FocalPoint[0] = 0.;
1500 m_FocalPoint[1] = 0.;
1501 m_FocalPoint[2] = 0.;
1505 /* -----------------------------------------------------------------------
1507 ----------------------------------------------------------------------- */
1509 template<class TInputImage, class TCoordRep>
1511 RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
1512 ::PrintSelf(std::ostream& os, Indent indent) const
1514 this->Superclass::PrintSelf(os,indent);
1516 os << indent << "Threshold: " << m_Threshold << std::endl;
1517 os << indent << "FocalPoint: " << m_FocalPoint << std::endl;
1518 os << indent << "Transform: " << m_Transform.GetPointer() << std::endl;
1519 os << indent << "Interpolator: " << m_Interpolator.GetPointer() << std::endl;
1523 /* -----------------------------------------------------------------------
1524 Evaluate at image index position
1525 ----------------------------------------------------------------------- */
1527 template<class TInputImage, class TCoordRep>
1528 typename RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
1530 RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
1531 ::Evaluate( const PointType& point ) const
1533 double integral = 0;
1535 OutputPointType transformedFocalPoint
1536 = m_Transform->TransformPoint( m_FocalPoint );
1538 DirectionType direction = transformedFocalPoint - point;
1540 RayCastHelper<TInputImage, TCoordRep> ray;
1541 ray.SetImage( this->m_Image );
1545 // (Aviv) Added support for images with non-zero origin
1546 // ray.SetRay(point, direction);
1547 ray.SetRay(point - this->m_Image->GetOrigin().GetVectorFromOrigin(), direction);
1548 ray.IntegrateAboveThreshold(integral, m_Threshold);
1550 return ( static_cast<OutputType>( integral ));
1553 template<class TInputImage, class TCoordRep>
1554 typename RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
1556 RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
1557 ::EvaluateAtContinuousIndex( const ContinuousIndexType& index ) const
1559 OutputPointType point;
1560 this->m_Image->TransformContinuousIndexToPhysicalPoint(index, point);
1562 return this->Evaluate( point );