+++ /dev/null
-/*=========================================================================
- Program: vv http://www.creatis.insa-lyon.fr/rio/vv
-
- Authors belong to:
- - University of LYON http://www.universite-lyon.fr/
- - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr
- - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
-
- This software is distributed WITHOUT ANY WARRANTY; without even
- the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
- PURPOSE. See the copyright notices for more information.
-
- It is distributed under dual licence
-
- - BSD See included LICENSE.txt file
- - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
-======================================================================-====*/
-#ifndef __clitkRayCastInterpolateImageFunctionWithOrigin_txx
-#define __clitkRayCastInterpolateImageFunctionWithOrigin_txx
-#include "clitkRayCastInterpolateImageFunctionWithOrigin.h"
-#include "itkContinuousIndex.h"
-#include "vnl/vnl_math.h"
-
-
-// Put the helper class in an anonymous namespace so that it is not
-// exposed to the user
-namespace {
-
- /** \class Helper class to maintain state when casting a ray.
- * This helper class keeps the RayCastInterpolateImageFunctionWithOrigin thread safe.
- */
- template <class TInputImage, class TCoordRep = float>
- class RayCastHelper
- {
- public:
- /** Constants for the image dimensions */
- itkStaticConstMacro(InputImageDimension, unsigned int,
- TInputImage::ImageDimension);
-
- /**
- * Type of the Transform Base class
- * The fixed image should be a 3D image
- */
- typedef itk::Transform<TCoordRep,3,3> TransformType;
-
- typedef typename TransformType::Pointer TransformPointer;
- typedef typename TransformType::InputPointType InputPointType;
- typedef typename TransformType::OutputPointType OutputPointType;
- typedef typename TransformType::ParametersType TransformParametersType;
- typedef typename TransformType::JacobianType TransformJacobianType;
-
- typedef typename TInputImage::SizeType SizeType;
- typedef itk::Vector<TCoordRep, 3> DirectionType;
- typedef itk::Point<TCoordRep, 3> PointType;
-
- typedef TInputImage InputImageType;
- typedef typename InputImageType::PixelType PixelType;
- typedef typename InputImageType::IndexType IndexType;
- typedef itk::ContinuousIndex<TCoordRep, 3> ContinuousIndexType;
-
-
- //JV Add an interpolator for the deformable transform
- typedef itk::InterpolateImageFunction<InputImageType, double> InterpolatorType;
- typedef typename InterpolatorType::Pointer InterpolatorPointer;
-
-
- /**
- * Set the image class
- */
- void SetImage(const InputImageType *input)
- {
- m_Image = input;
- }
-
- /**
- * Initialise the ray using the position and direction of a line.
- *
- * \param RayPosn The position of the ray in 3D (mm).
- * \param RayDirn The direction of the ray in 3D (mm).
- *
- * \return True if this is a valid ray.
- */
- bool SetRay(OutputPointType RayPosn, DirectionType RayDirn);
-
-
- /** \brief
- * Integrate the interpolated intensities along the ray and
- * return the result.
- *
- * This routine can be called after instantiating the ray and
- * calling SetProjectionCoord2D() or Reset(). It may then be called
- * as many times thereafter for different 2D projection
- * coordinates.
- *
- * \param integral The integrated intensities along the ray.
- *
- * \return True if a valid ray was specified.
- */
- bool Integrate(double &integral)
- {
- return IntegrateAboveThreshold(integral, 0);
- };
-
-
- /** \brief
- * Integrate the interpolated intensities above a given threshold,
- * along the ray and return the result.
- *
- * This routine can be called after instantiating the ray and
- * calling SetProjectionCoord2D() or Reset(). It may then be called
- * as many times thereafter for different 2D projection
- * coordinates.
- *
- * \param integral The integrated intensities along the ray.
- * \param threshold The integration threshold [default value: 0]
- *
- * \return True if a valid ray was specified.
- */
- bool IntegrateAboveThreshold(double &integral, double threshold);
-
- /** \brief
- * Increment each of the intensities of the 4 planar voxels
- * surrounding the current ray point.
- *
- * \parameter increment Intensity increment for each of the current 4 voxels
- */
- void IncrementIntensities(double increment=1);
-
- /// Reset the iterator to the start of the ray.
- void Reset(void);
-
- /// Return the interpolated intensity of the current ray point.
- double GetCurrentIntensity(void) const;
-
- /// Return the ray point spacing in mm
- double GetRayPointSpacing(void) const {
- typename InputImageType::SpacingType spacing=this->m_Image->GetSpacing();
-
- if (m_ValidRay)
- return vcl_sqrt(m_VoxelIncrement[0]*spacing[0]*m_VoxelIncrement[0]*spacing[0]
- + m_VoxelIncrement[1]*spacing[1]*m_VoxelIncrement[1]*spacing[1]
- + m_VoxelIncrement[2]*spacing[2]*m_VoxelIncrement[2]*spacing[2] );
- else
- return 0.;
- };
-
- /// Set the initial zero state of the object
- void ZeroState();
-
- /// Initialise the object
- void Initialise(void);
-
- //JV set transform and interpolator
- void SetTransformIsDeformable(bool m){m_TransformIsDeformable=m;}
- void SetTransform (TransformType* m){m_Transform=m;}
- void SetDeformableTransformInterpolator(InterpolatorType* m){m_DeformableTransformInterpolator=m;}
-
- protected:
- /// Calculate the endpoint coordinats of the ray in voxels.
- void EndPointsInVoxels(void);
-
- /**
- * Calculate the incremental direction vector in voxels, 'dVoxel',
- * required to traverse the ray.
- */
- void CalcDirnVector(void);
-
- /**
- * Reduce the length of the ray until both start and end
- * coordinates lie inside the volume.
- *
- * \return True if a valid ray has been, false otherwise.
- */
- bool AdjustRayLength(void);
-
- /**
- * Obtain pointers to the four voxels surrounding the point where the ray
- * enters the volume.
- */
- void InitialiseVoxelPointers(void);
-
- /// Increment the voxel pointers surrounding the current point on the ray.
- void IncrementVoxelPointers(void);
-
- /// Record volume dimensions and resolution
- void RecordVolumeDimensions(void);
-
- /// Define the corners of the volume
- void DefineCorners(void);
-
- /** \brief
- * Calculate the planes which define the volume.
- *
- * Member function to calculate the equations of the planes of 4 of
- * the sides of the volume, calculate the positions of the 8 corners
- * of the volume in mm in World, also calculate the values of the
- * slopes of the lines which go to make up the volume( defined as
- * lines in cube x,y,z dirn and then each of these lines has a slope
- * in the world x,y,z dirn [3]) and finally also to return the length
- * of the sides of the lines in mm.
- */
- void CalcPlanesAndCorners(void);
-
- /** \brief
- * Calculate the ray intercepts with the volume.
- *
- * See where the ray cuts the volume, check that truncation does not occur,
- * if not, then start ray where it first intercepts the volume and set
- * x_max to be where it leaves the volume.
- *
- * \return True if a valid ray has been specified, false otherwise.
- */
- bool CalcRayIntercepts(void);
-
- /**
- * The ray is traversed by stepping in the axial direction
- * that enables the greatest number of planes in the volume to be
- * intercepted.
- */
- typedef enum {
- UNDEFINED_DIRECTION=0, //!< Undefined
- TRANSVERSE_IN_X, //!< x
- TRANSVERSE_IN_Y, //!< y
- TRANSVERSE_IN_Z, //!< z
- LAST_DIRECTION
- } TraversalDirection;
-
- // Cache the image in the structure. Skip the smart pointer for
- // efficiency. This inner class will go in/out of scope with every
- // call to Evaluate()
- const InputImageType *m_Image;
-
- /// Flag indicating whether the current ray is valid
- bool m_ValidRay;
-
- /** \brief
- * The start position of the ray in voxels.
- *
- * NB. Two of the components of this coordinate (i.e. those lying within
- * the planes of voxels being traversed) will be shifted by half a
- * voxel. This enables indices of the neighbouring voxels within the plane
- * to be determined by simply casting to 'int' and optionally adding 1.
- */
- double m_RayVoxelStartPosition[3];
-
- /** \brief
- * The end coordinate of the ray in voxels.
- *
- * NB. Two of the components of this coordinate (i.e. those lying within
- * the planes of voxels being traversed) will be shifted by half a
- * voxel. This enables indices of the neighbouring voxels within the plane
- * to be determined by simply casting to 'int' and optionally adding 1.
- */
- double m_RayVoxelEndPosition[3];
-
-
- /** \brief
- * The current coordinate on the ray in voxels.
- *
- * NB. Two of the components of this coordinate (i.e. those lying within
- * the planes of voxels being traversed) will be shifted by half a
- * voxel. This enables indices of the neighbouring voxels within the plane
- * to be determined by simply casting to 'int' and optionally adding 1.
- */
- double m_Position3Dvox[3];
-
- /** The incremental direction vector of the ray in voxels. */
- double m_VoxelIncrement[3];
-
- /// The direction in which the ray is incremented thorough the volume (x, y or z).
- TraversalDirection m_TraversalDirection;
-
- /// The total number of planes of voxels traversed by the ray.
- int m_TotalRayVoxelPlanes;
-
- /// The current number of planes of voxels traversed by the ray.
- int m_NumVoxelPlanesTraversed;
-
- /// Pointers to the current four voxels surrounding the ray's trajectory.
- const PixelType *m_RayIntersectionVoxels[4];
-
- /**
- * The voxel coordinate of the bottom-left voxel of the current
- * four voxels surrounding the ray's trajectory.
- */
- int m_RayIntersectionVoxelIndex[3];
-
- /// The dimension in voxels of the 3D volume in along the x axis
- int m_NumberOfVoxelsInX;
- /// The dimension in voxels of the 3D volume in along the y axis
- int m_NumberOfVoxelsInY;
- /// The dimension in voxels of the 3D volume in along the z axis
- int m_NumberOfVoxelsInZ;
-
- /// Voxel dimension in x
- double m_VoxelDimensionInX;
- /// Voxel dimension in y
- double m_VoxelDimensionInY;
- /// Voxel dimension in z
- double m_VoxelDimensionInZ;
-
- /// The coordinate of the point at which the ray enters the volume in mm.
- double m_RayStartCoordInMM[3];
- /// The coordinate of the point at which the ray exits the volume in mm.
- double m_RayEndCoordInMM[3];
-
-
- /** \brief
- Planes which define the boundary of the volume in mm
- (six planes and four parameters: Ax+By+Cz+D). */
- double m_BoundingPlane[6][4];
- /// The eight corners of the volume (x,y,z coordinates for each).
- double m_BoundingCorner[8][3];
-
- /// The position of the ray
- double m_CurrentRayPositionInMM[3];
-
- /// The direction of the ray
- double m_RayDirectionInMM[3];
-
- //JV transform is deformable
- bool m_TransformIsDeformable;
- TransformType* m_Transform;
- InterpolatorType* m_DeformableTransformInterpolator;
-
- };
-
- /* -----------------------------------------------------------------------
- Initialise() - Initialise the object
- ----------------------------------------------------------------------- */
-
- template<class TInputImage, class TCoordRep>
- void
- RayCastHelper<TInputImage, TCoordRep>
- ::Initialise(void)
- {
- // Save the dimensions of the volume and calculate the bounding box
- this->RecordVolumeDimensions();
-
- // Calculate the planes and corners which define the volume.
- this->DefineCorners();
- this->CalcPlanesAndCorners();
- }
-
-
- /* -----------------------------------------------------------------------
- RecordVolumeDimensions() - Record volume dimensions and resolution
- ----------------------------------------------------------------------- */
-
- template<class TInputImage, class TCoordRep>
- void
- RayCastHelper<TInputImage, TCoordRep>
- ::RecordVolumeDimensions(void)
- {
- typename InputImageType::SpacingType spacing=this->m_Image->GetSpacing();
- SizeType dim=this->m_Image->GetLargestPossibleRegion().GetSize();
-
- m_NumberOfVoxelsInX = dim[0];
- m_NumberOfVoxelsInY = dim[1];
- m_NumberOfVoxelsInZ = dim[2];
-
- m_VoxelDimensionInX = spacing[0];
- m_VoxelDimensionInY = spacing[1];
- m_VoxelDimensionInZ = spacing[2];
- }
-
-
- /* -----------------------------------------------------------------------
- DefineCorners() - Define the corners of the volume
- ----------------------------------------------------------------------- */
-
- template<class TInputImage, class TCoordRep>
- void
- RayCastHelper<TInputImage, TCoordRep>
- ::DefineCorners(void)
- {
- // Define corner positions as if at the origin
-
- m_BoundingCorner[0][0] =
- m_BoundingCorner[1][0] =
- m_BoundingCorner[2][0] =
- m_BoundingCorner[3][0] = 0;
-
- m_BoundingCorner[4][0] =
- m_BoundingCorner[5][0] =
- m_BoundingCorner[6][0] =
- m_BoundingCorner[7][0] = m_VoxelDimensionInX*m_NumberOfVoxelsInX;
-
- m_BoundingCorner[1][1] =
- m_BoundingCorner[3][1] =
- m_BoundingCorner[5][1] =
- m_BoundingCorner[7][1] = m_VoxelDimensionInY*m_NumberOfVoxelsInY;
-
- m_BoundingCorner[0][1] =
- m_BoundingCorner[2][1] =
- m_BoundingCorner[4][1] =
- m_BoundingCorner[6][1] = 0;
-
- m_BoundingCorner[0][2] =
- m_BoundingCorner[1][2] =
- m_BoundingCorner[4][2] =
- m_BoundingCorner[5][2] =
- m_VoxelDimensionInZ*m_NumberOfVoxelsInZ;
-
- m_BoundingCorner[2][2] =
- m_BoundingCorner[3][2] =
- m_BoundingCorner[6][2] =
- m_BoundingCorner[7][2] = 0;
-
- }
-
- /* -----------------------------------------------------------------------
- CalcPlanesAndCorners() - Calculate the planes and corners of the volume.
- ----------------------------------------------------------------------- */
-
- template<class TInputImage, class TCoordRep>
- void
- RayCastHelper<TInputImage, TCoordRep>
- ::CalcPlanesAndCorners(void)
- {
- int j;
-
-
- // find the equations of the planes
-
- int c1=0, c2=0, c3=0;
-
- for (j=0; j<6; j++)
- { // loop around for planes
- switch (j)
- { // which corners to take
- case 0:
- c1=1; c2=2; c3=3;
- break;
- case 1:
- c1=4; c2=5; c3=6;
- break;
- case 2:
- c1=5; c2=3; c3=7;
- break;
- case 3:
- c1=2; c2=4; c3=6;
- break;
- case 4:
- c1=1; c2=5; c3=0;
- break;
- case 5:
- c1=3; c2=7; c3=2;
- break;
- }
-
-
- double line1x, line1y, line1z;
- double line2x, line2y, line2z;
-
- // lines from one corner to another in x,y,z dirns
- line1x = m_BoundingCorner[c1][0] - m_BoundingCorner[c2][0];
- line2x = m_BoundingCorner[c1][0] - m_BoundingCorner[c3][0];
-
- line1y = m_BoundingCorner[c1][1] - m_BoundingCorner[c2][1];
- line2y = m_BoundingCorner[c1][1] - m_BoundingCorner[c3][1];
-
- line1z = m_BoundingCorner[c1][2] - m_BoundingCorner[c2][2];
- line2z = m_BoundingCorner[c1][2] - m_BoundingCorner[c3][2];
-
- double A, B, C, D;
-
- // take cross product
- A = line1y*line2z - line2y*line1z;
- B = line2x*line1z - line1x*line2z;
- C = line1x*line2y - line2x*line1y;
-
- // find constant
- D = -( A*m_BoundingCorner[c1][0]
- + B*m_BoundingCorner[c1][1]
- + C*m_BoundingCorner[c1][2] );
-
- // initialise plane value and normalise
- m_BoundingPlane[j][0] = A/vcl_sqrt(A*A + B*B + C*C);
- m_BoundingPlane[j][1] = B/vcl_sqrt(A*A + B*B + C*C);
- m_BoundingPlane[j][2] = C/vcl_sqrt(A*A + B*B + C*C);
- m_BoundingPlane[j][3] = D/vcl_sqrt(A*A + B*B + C*C);
-
- if ( (A*A + B*B + C*C) == 0 )
- {
- itk::ExceptionObject err(__FILE__, __LINE__);
- err.SetLocation( ITK_LOCATION );
- err.SetDescription( "Division by zero (planes) "
- "- CalcPlanesAndCorners().");
- throw err;
- }
- }
-
- }
-
-
- /* -----------------------------------------------------------------------
- CalcRayIntercepts() - Calculate the ray intercepts with the volume.
- ----------------------------------------------------------------------- */
-
- template<class TInputImage, class TCoordRep>
- bool
- RayCastHelper<TInputImage, TCoordRep>
- ::CalcRayIntercepts()
- {
- double maxInterDist, interDist;
- double cornerVect[4][3];
- int cross[4][3], noInterFlag[6];
- int nSidesCrossed, crossFlag, c[4];
- double ax, ay, az, bx, by, bz;
- double cubeInter[6][3];
- double denom;
-
- int i,j, k;
- int NoSides = 6; // =6 to allow truncation: =4 to remove truncated rays
-
- // Calculate intercept of ray with planes
-
- double interceptx[6], intercepty[6], interceptz[6];
- double d[6];
-
- for( j=0; j<NoSides; j++)
- {
-
- denom = ( m_BoundingPlane[j][0]*m_RayDirectionInMM[0]
- + m_BoundingPlane[j][1]*m_RayDirectionInMM[1]
- + m_BoundingPlane[j][2]*m_RayDirectionInMM[2]);
-
- if( (long)(denom*100) != 0 )
- {
- d[j] = -( m_BoundingPlane[j][3]
- + m_BoundingPlane[j][0]*m_CurrentRayPositionInMM[0]
- + m_BoundingPlane[j][1]*m_CurrentRayPositionInMM[1]
- + m_BoundingPlane[j][2]*m_CurrentRayPositionInMM[2] ) / denom;
-
- interceptx[j] = m_CurrentRayPositionInMM[0] + d[j]*m_RayDirectionInMM[0];
- intercepty[j] = m_CurrentRayPositionInMM[1] + d[j]*m_RayDirectionInMM[1];
- interceptz[j] = m_CurrentRayPositionInMM[2] + d[j]*m_RayDirectionInMM[2];
-
- noInterFlag[j] = 1; //OK
- }
- else
- {
- noInterFlag[j] = 0; //NOT OK
- }
- }
-
-
- nSidesCrossed = 0;
- for( j=0; j<NoSides; j++ )
- {
-
- // Work out which corners to use
-
- if( j==0 )
- {
- c[0] = 0; c[1] = 1; c[2] = 3; c[3] = 2;
- }
- else if( j==1 )
- {
- c[0] = 4; c[1] = 5; c[2] = 7; c[3] = 6;
- }
- else if( j==2 )
- {
- c[0] = 1; c[1] = 5; c[2] = 7; c[3] = 3;
- }
- else if( j==3 )
- {
- c[0] = 0; c[1] = 2; c[2] = 6; c[3] = 4;
- }
- else if( j==4 )
- { //TOP
- c[0] = 0; c[1] = 1; c[2] = 5; c[3] = 4;
- }
- else if( j==5 )
- { //BOTTOM
- c[0] = 2; c[1] = 3; c[2] = 7; c[3] = 6;
- }
-
- // Calculate vectors from corner of ct volume to intercept.
- for( i=0; i<4; i++ )
- {
- if( noInterFlag[j]==1 )
- {
- cornerVect[i][0] = m_BoundingCorner[c[i]][0] - interceptx[j];
- cornerVect[i][1] = m_BoundingCorner[c[i]][1] - intercepty[j];
- cornerVect[i][2] = m_BoundingCorner[c[i]][2] - interceptz[j];
- }
- else if( noInterFlag[j]==0 )
- {
- cornerVect[i][0] = 0;
- cornerVect[i][1] = 0;
- cornerVect[i][2] = 0;
- }
-
- }
-
- // Do cross product with these vectors
- for( i=0; i<4; i++ )
- {
- if( i==3 )
- {
- k = 0;
- }
- else
- {
- k = i+1;
- }
- ax = cornerVect[i][0];
- ay = cornerVect[i][1];
- az = cornerVect[i][2];
- bx = cornerVect[k][0];
- by = cornerVect[k][1];
- bz = cornerVect[k][2];
-
- // The int and divide by 100 are to avoid rounding errors. If
- // these are not included then you get values fluctuating around
- // zero and so in the subsequent check, all the values are not
- // above or below zero. NB. If you "INT" by too much here though
- // you can get problems in the corners of your volume when rays
- // are allowed to go through more than one plane.
- cross[i][0] = (int)((ay*bz - az*by)/100);
- cross[i][1] = (int)((az*bx - ax*bz)/100);
- cross[i][2] = (int)((ax*by - ay*bx)/100);
- }
-
- // See if a sign change occured between all these cross products
- // if not, then the ray went through this plane
-
- crossFlag=0;
- for( i=0; i<3; i++ )
- {
- if( ( cross[0][i]<=0
- && cross[1][i]<=0
- && cross[2][i]<=0
- && cross[3][i]<=0)
-
- || ( cross[0][i]>=0
- && cross[1][i]>=0
- && cross[2][i]>=0
- && cross[3][i]>=0) )
- {
- crossFlag++;
- }
- }
-
-
- if( crossFlag==3 && noInterFlag[j]==1 )
- {
- cubeInter[nSidesCrossed][0] = interceptx[j];
- cubeInter[nSidesCrossed][1] = intercepty[j];
- cubeInter[nSidesCrossed][2] = interceptz[j];
- nSidesCrossed++;
- }
-
- } // End of loop over all four planes
-
- m_RayStartCoordInMM[0] = cubeInter[0][0];
- m_RayStartCoordInMM[1] = cubeInter[0][1];
- m_RayStartCoordInMM[2] = cubeInter[0][2];
-
- m_RayEndCoordInMM[0] = cubeInter[1][0];
- m_RayEndCoordInMM[1] = cubeInter[1][1];
- m_RayEndCoordInMM[2] = cubeInter[1][2];
-
- if( nSidesCrossed >= 5 )
- {
- std::cerr << "WARNING: No. of sides crossed equals: " << nSidesCrossed << std::endl;
- }
-
- // If 'nSidesCrossed' is larger than 2, this means that the ray goes through
- // a corner of the volume and due to rounding errors, the ray is
- // deemed to go through more than two planes. To obtain the correct
- // start and end positions we choose the two intercept values which
- // are furthest from each other.
-
-
- if( nSidesCrossed >= 3 )
- {
- maxInterDist = 0;
- for( j=0; j<nSidesCrossed-1; j++ )
- {
- for( k=j+1; k<nSidesCrossed; k++ )
- {
- interDist = 0;
- for( i=0; i<3; i++ )
- {
- interDist += (cubeInter[j][i] - cubeInter[k][i])*
- (cubeInter[j][i] - cubeInter[k][i]);
- }
- if( interDist > maxInterDist )
- {
- maxInterDist = interDist;
-
- m_RayStartCoordInMM[0] = cubeInter[j][0];
- m_RayStartCoordInMM[1] = cubeInter[j][1];
- m_RayStartCoordInMM[2] = cubeInter[j][2];
-
- m_RayEndCoordInMM[0] = cubeInter[k][0];
- m_RayEndCoordInMM[1] = cubeInter[k][1];
- m_RayEndCoordInMM[2] = cubeInter[k][2];
- }
- }
- }
- nSidesCrossed = 2;
- }
-
- if (nSidesCrossed == 2 )
- {
- return true;
- }
- else
- {
- return false;
- }
- }
-
-
- /* -----------------------------------------------------------------------
- SetRay() - Set the position and direction of the ray
- ----------------------------------------------------------------------- */
-
- template<class TInputImage, class TCoordRep>
- bool
- RayCastHelper<TInputImage, TCoordRep>
- ::SetRay(OutputPointType RayPosn, DirectionType RayDirn)
- {
-
- // Store the position and direction of the ray
- typename TInputImage::SpacingType spacing=this->m_Image->GetSpacing();
- SizeType dim=this->m_Image->GetLargestPossibleRegion().GetSize();
-
- // we need to translate the _center_ of the volume to the origin
- m_NumberOfVoxelsInX = dim[0];
- m_NumberOfVoxelsInY = dim[1];
- m_NumberOfVoxelsInZ = dim[2];
-
- m_VoxelDimensionInX = spacing[0];
- m_VoxelDimensionInY = spacing[1];
- m_VoxelDimensionInZ = spacing[2];
-
- // (Aviv) Incorporating a fix by Jian Wu
- // http://public.kitware.com/pipermail/insight-users/2006-March/017265.html
- m_CurrentRayPositionInMM[0] = RayPosn[0];
- m_CurrentRayPositionInMM[1] = RayPosn[1];
- m_CurrentRayPositionInMM[2] = RayPosn[2];
-
- m_RayDirectionInMM[0] = RayDirn[0];
- m_RayDirectionInMM[1] = RayDirn[1];
- m_RayDirectionInMM[2] = RayDirn[2];
-
- // Compute the ray path for this coordinate in mm
-
- m_ValidRay = this->CalcRayIntercepts();
-
- if (! m_ValidRay)
- {
- Reset();
- return false;
- }
-
- // Convert the start and end coordinates of the ray to voxels
-
- this->EndPointsInVoxels();
-
- /* Calculate the ray direction vector in voxels and hence the voxel
- increment required to traverse the ray, and the number of
- interpolation points on the ray.
-
- This routine also shifts the coordinate frame by half a voxel for
- two of the directional components (i.e. those lying within the
- planes of voxels being traversed). */
-
- this->CalcDirnVector();
-
-
- /* Reduce the length of the ray until both start and end
- coordinates lie inside the volume. */
-
- m_ValidRay = this->AdjustRayLength();
-
- // Reset the iterator to the start of the ray.
-
- Reset();
-
- return m_ValidRay;
- }
-
-
- /* -----------------------------------------------------------------------
- EndPointsInVoxels() - Convert the endpoints to voxels
- ----------------------------------------------------------------------- */
-
- template<class TInputImage, class TCoordRep>
- void
- RayCastHelper<TInputImage, TCoordRep>
- ::EndPointsInVoxels(void)
- {
- m_RayVoxelStartPosition[0] = m_RayStartCoordInMM[0]/m_VoxelDimensionInX;
- m_RayVoxelStartPosition[1] = m_RayStartCoordInMM[1]/m_VoxelDimensionInY;
- m_RayVoxelStartPosition[2] = m_RayStartCoordInMM[2]/m_VoxelDimensionInZ;
-
- m_RayVoxelEndPosition[0] = m_RayEndCoordInMM[0]/m_VoxelDimensionInX;
- m_RayVoxelEndPosition[1] = m_RayEndCoordInMM[1]/m_VoxelDimensionInY;
- m_RayVoxelEndPosition[2] = m_RayEndCoordInMM[2]/m_VoxelDimensionInZ;
-
- }
-
-
- /* -----------------------------------------------------------------------
- CalcDirnVector() - Calculate the incremental direction vector in voxels.
- ----------------------------------------------------------------------- */
-
- template<class TInputImage, class TCoordRep>
- void
- RayCastHelper<TInputImage, TCoordRep>
- ::CalcDirnVector(void)
- {
- double xNum, yNum, zNum;
-
- // Calculate the number of voxels in each direction
-
- xNum = vcl_fabs(m_RayVoxelStartPosition[0] - m_RayVoxelEndPosition[0]);
- yNum = vcl_fabs(m_RayVoxelStartPosition[1] - m_RayVoxelEndPosition[1]);
- zNum = vcl_fabs(m_RayVoxelStartPosition[2] - m_RayVoxelEndPosition[2]);
-
- // The direction iterated in is that with the greatest number of voxels
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
- // Iterate in X direction
-
- if( (xNum >= yNum) && (xNum >= zNum) )
- {
- if( m_RayVoxelStartPosition[0] < m_RayVoxelEndPosition[0] )
- {
- m_VoxelIncrement[0] = 1;
-
- m_VoxelIncrement[1]
- = (m_RayVoxelStartPosition[1]
- - m_RayVoxelEndPosition[1])/(m_RayVoxelStartPosition[0]
- - m_RayVoxelEndPosition[0]);
-
- m_VoxelIncrement[2]
- = (m_RayVoxelStartPosition[2]
- - m_RayVoxelEndPosition[2])/(m_RayVoxelStartPosition[0]
- - m_RayVoxelEndPosition[0]);
- }
- else
- {
- m_VoxelIncrement[0] = -1;
-
- m_VoxelIncrement[1]
- = -(m_RayVoxelStartPosition[1]
- - m_RayVoxelEndPosition[1])/(m_RayVoxelStartPosition[0]
- - m_RayVoxelEndPosition[0]);
-
- m_VoxelIncrement[2]
- = -(m_RayVoxelStartPosition[2]
- - m_RayVoxelEndPosition[2])/(m_RayVoxelStartPosition[0]
- - m_RayVoxelEndPosition[0]);
- }
-
- // This section is to alter the start position in order to
- // place the center of the voxels in there correct positions,
- // rather than placing them at the corner of voxels which is
- // what happens if this is not carried out. The reason why
- // x has no -0.5 is because this is the direction we are going
- // to iterate in and therefore we wish to go from center to
- // center rather than finding the surrounding voxels.
-
- m_RayVoxelStartPosition[1] += ( (int)m_RayVoxelStartPosition[0]
- - m_RayVoxelStartPosition[0])*m_VoxelIncrement[1]*m_VoxelIncrement[0]
- + 0.5*m_VoxelIncrement[1] - 0.5;
-
- m_RayVoxelStartPosition[2] += ( (int)m_RayVoxelStartPosition[0]
- - m_RayVoxelStartPosition[0])*m_VoxelIncrement[2]*m_VoxelIncrement[0]
- + 0.5*m_VoxelIncrement[2] - 0.5;
-
- m_RayVoxelStartPosition[0] = (int)m_RayVoxelStartPosition[0] + 0.5*m_VoxelIncrement[0];
-
- m_TotalRayVoxelPlanes = (int)xNum;
-
- m_TraversalDirection = TRANSVERSE_IN_X;
- }
-
- // Iterate in Y direction
-
- else if( (yNum >= xNum) && (yNum >= zNum) )
- {
-
- if( m_RayVoxelStartPosition[1] < m_RayVoxelEndPosition[1] )
- {
- m_VoxelIncrement[1] = 1;
-
- m_VoxelIncrement[0]
- = (m_RayVoxelStartPosition[0]
- - m_RayVoxelEndPosition[0])/(m_RayVoxelStartPosition[1]
- - m_RayVoxelEndPosition[1]);
-
- m_VoxelIncrement[2]
- = (m_RayVoxelStartPosition[2]
- - m_RayVoxelEndPosition[2])/(m_RayVoxelStartPosition[1]
- - m_RayVoxelEndPosition[1]);
- }
- else
- {
- m_VoxelIncrement[1] = -1;
-
- m_VoxelIncrement[0]
- = -(m_RayVoxelStartPosition[0]
- - m_RayVoxelEndPosition[0])/(m_RayVoxelStartPosition[1]
- - m_RayVoxelEndPosition[1]);
-
- m_VoxelIncrement[2]
- = -(m_RayVoxelStartPosition[2]
- - m_RayVoxelEndPosition[2])/(m_RayVoxelStartPosition[1]
- - m_RayVoxelEndPosition[1]);
- }
-
-
- m_RayVoxelStartPosition[0] += ( (int)m_RayVoxelStartPosition[1]
- - m_RayVoxelStartPosition[1])*m_VoxelIncrement[0]*m_VoxelIncrement[1]
- + 0.5*m_VoxelIncrement[0] - 0.5;
-
- m_RayVoxelStartPosition[2] += ( (int)m_RayVoxelStartPosition[1]
- - m_RayVoxelStartPosition[1])*m_VoxelIncrement[2]*m_VoxelIncrement[1]
- + 0.5*m_VoxelIncrement[2] - 0.5;
-
- m_RayVoxelStartPosition[1] = (int)m_RayVoxelStartPosition[1] + 0.5*m_VoxelIncrement[1];
-
- m_TotalRayVoxelPlanes = (int)yNum;
-
- m_TraversalDirection = TRANSVERSE_IN_Y;
- }
-
- // Iterate in Z direction
-
- else
- {
-
- if( m_RayVoxelStartPosition[2] < m_RayVoxelEndPosition[2] )
- {
- m_VoxelIncrement[2] = 1;
-
- m_VoxelIncrement[0]
- = (m_RayVoxelStartPosition[0]
- - m_RayVoxelEndPosition[0])/(m_RayVoxelStartPosition[2]
- - m_RayVoxelEndPosition[2]);
-
- m_VoxelIncrement[1]
- = (m_RayVoxelStartPosition[1]
- - m_RayVoxelEndPosition[1])/(m_RayVoxelStartPosition[2]
- - m_RayVoxelEndPosition[2]);
- }
- else
- {
- m_VoxelIncrement[2] = -1;
-
- m_VoxelIncrement[0]
- = -(m_RayVoxelStartPosition[0]
- - m_RayVoxelEndPosition[0])/(m_RayVoxelStartPosition[2]
- - m_RayVoxelEndPosition[2]);
-
- m_VoxelIncrement[1]
- = -(m_RayVoxelStartPosition[1]
- - m_RayVoxelEndPosition[1])/(m_RayVoxelStartPosition[2]
- - m_RayVoxelEndPosition[2]);
- }
-
- m_RayVoxelStartPosition[0] += ( (int)m_RayVoxelStartPosition[2]
- - m_RayVoxelStartPosition[2])*m_VoxelIncrement[0]*m_VoxelIncrement[2]
- + 0.5*m_VoxelIncrement[0] - 0.5;
-
- m_RayVoxelStartPosition[1] += ( (int)m_RayVoxelStartPosition[2]
- - m_RayVoxelStartPosition[2])*m_VoxelIncrement[1]*m_VoxelIncrement[2]
- + 0.5*m_VoxelIncrement[1] - 0.5;
-
- m_RayVoxelStartPosition[2] = (int)m_RayVoxelStartPosition[2] + 0.5*m_VoxelIncrement[2];
-
- m_TotalRayVoxelPlanes = (int)zNum;
-
- m_TraversalDirection = TRANSVERSE_IN_Z;
- }
- }
-
-
- /* -----------------------------------------------------------------------
- AdjustRayLength() - Ensure that the ray lies within the volume
- ----------------------------------------------------------------------- */
-
- template<class TInputImage, class TCoordRep>
- bool
- RayCastHelper<TInputImage, TCoordRep>
- ::AdjustRayLength(void)
- {
- bool startOK, endOK;
-
- int Istart[3];
- int Idirn[3];
-
- if (m_TraversalDirection == TRANSVERSE_IN_X)
- {
- Idirn[0] = 0;
- Idirn[1] = 1;
- Idirn[2] = 1;
- }
- else if (m_TraversalDirection == TRANSVERSE_IN_Y)
- {
- Idirn[0] = 1;
- Idirn[1] = 0;
- Idirn[2] = 1;
- }
- else if (m_TraversalDirection == TRANSVERSE_IN_Z)
- {
- Idirn[0] = 1;
- Idirn[1] = 1;
- Idirn[2] = 0;
- }
- else
- {
- itk::ExceptionObject err(__FILE__, __LINE__);
- err.SetLocation( ITK_LOCATION );
- err.SetDescription( "The ray traversal direction is unset "
- "- AdjustRayLength().");
- throw err;
- return false;
- }
-
-
- do
- {
-
- startOK = false;
- endOK = false;
-
- Istart[0] = (int) vcl_floor(m_RayVoxelStartPosition[0]);
- Istart[1] = (int) vcl_floor(m_RayVoxelStartPosition[1]);
- Istart[2] = (int) vcl_floor(m_RayVoxelStartPosition[2]);
-
- if( (Istart[0] >= 0) && (Istart[0] + Idirn[0] < m_NumberOfVoxelsInX) &&
- (Istart[1] >= 0) && (Istart[1] + Idirn[1] < m_NumberOfVoxelsInY) &&
- (Istart[2] >= 0) && (Istart[2] + Idirn[2] < m_NumberOfVoxelsInZ) )
- {
-
- startOK = true;
- }
- else
- {
- m_RayVoxelStartPosition[0] += m_VoxelIncrement[0];
- m_RayVoxelStartPosition[1] += m_VoxelIncrement[1];
- m_RayVoxelStartPosition[2] += m_VoxelIncrement[2];
-
- m_TotalRayVoxelPlanes--;
- }
-
- Istart[0] = (int) vcl_floor(m_RayVoxelStartPosition[0]
- + m_TotalRayVoxelPlanes*m_VoxelIncrement[0]);
-
- Istart[1] = (int) vcl_floor(m_RayVoxelStartPosition[1]
- + m_TotalRayVoxelPlanes*m_VoxelIncrement[1]);
-
- Istart[2] = (int) vcl_floor(m_RayVoxelStartPosition[2]
- + m_TotalRayVoxelPlanes*m_VoxelIncrement[2]);
-
- if( (Istart[0] >= 0) && (Istart[0] + Idirn[0] < m_NumberOfVoxelsInX) &&
- (Istart[1] >= 0) && (Istart[1] + Idirn[1] < m_NumberOfVoxelsInY) &&
- (Istart[2] >= 0) && (Istart[2] + Idirn[2] < m_NumberOfVoxelsInZ) )
- {
-
- endOK = true;
- }
- else
- {
- m_TotalRayVoxelPlanes--;
- }
-
- } while ( (! (startOK && endOK)) && (m_TotalRayVoxelPlanes > 1) );
-
-
- return (startOK && endOK);
- }
-
-
- /* -----------------------------------------------------------------------
- Reset() - Reset the iterator to the start of the ray.
- ----------------------------------------------------------------------- */
-
- template<class TInputImage, class TCoordRep>
- void
- RayCastHelper<TInputImage, TCoordRep>
- ::Reset(void)
- {
- int i;
-
- m_NumVoxelPlanesTraversed = -1;
-
- // If this is a valid ray...
-
- if (m_ValidRay)
- {
- for (i=0; i<3; i++)
- {
- m_Position3Dvox[i] = m_RayVoxelStartPosition[i];
- }
- this->InitialiseVoxelPointers();
- }
-
- // otherwise set parameters to zero
-
- else
- {
- for (i=0; i<3; i++)
- {
- m_RayVoxelStartPosition[i] = 0.;
- }
- for (i=0; i<3; i++)
- {
- m_RayVoxelEndPosition[i] = 0.;
- }
- for (i=0; i<3; i++)
- {
- m_VoxelIncrement[i] = 0.;
- }
- m_TraversalDirection = UNDEFINED_DIRECTION;
-
- m_TotalRayVoxelPlanes = 0;
-
- for (i=0; i<4; i++)
- {
- m_RayIntersectionVoxels[i] = 0;
- }
- for (i=0; i<3; i++)
- {
- m_RayIntersectionVoxelIndex[i] = 0;
- }
- }
- }
-
-
- /* -----------------------------------------------------------------------
- InitialiseVoxelPointers() - Obtain pointers to the first four voxels
- ----------------------------------------------------------------------- */
-
- template<class TInputImage, class TCoordRep>
- void
- RayCastHelper<TInputImage, TCoordRep>
- ::InitialiseVoxelPointers(void)
- {
-
- //JV don't do anything if deformable
- if(!m_TransformIsDeformable)
- {
- IndexType index;
-
- int Ix, Iy, Iz;
-
- Ix = (int)(m_RayVoxelStartPosition[0]);
- Iy = (int)(m_RayVoxelStartPosition[1]);
- Iz = (int)(m_RayVoxelStartPosition[2]);
-
- m_RayIntersectionVoxelIndex[0] = Ix;
- m_RayIntersectionVoxelIndex[1] = Iy;
- m_RayIntersectionVoxelIndex[2] = Iz;
-
- switch( m_TraversalDirection )
- {
- case TRANSVERSE_IN_X:
- {
-
- if( (Ix >= 0) && (Ix < m_NumberOfVoxelsInX) &&
- (Iy >= 0) && (Iy + 1 < m_NumberOfVoxelsInY) &&
- (Iz >= 0) && (Iz + 1 < m_NumberOfVoxelsInZ))
- {
- index[0]=Ix; index[1]=Iy; index[2]=Iz;
- m_RayIntersectionVoxels[0]
- = this->m_Image->GetBufferPointer() + this->m_Image->ComputeOffset(index);
-
- index[0]=Ix; index[1]=Iy+1; index[2]=Iz;
- m_RayIntersectionVoxels[1]
- = ( this->m_Image->GetBufferPointer() + this->m_Image->ComputeOffset(index) );
-
- index[0]=Ix; index[1]=Iy; index[2]=Iz+1;
- m_RayIntersectionVoxels[2]
- = ( this->m_Image->GetBufferPointer() + this->m_Image->ComputeOffset(index) );
-
- index[0]=Ix; index[1]=Iy+1; index[2]=Iz+1;
- m_RayIntersectionVoxels[3]
- = ( this->m_Image->GetBufferPointer() + this->m_Image->ComputeOffset(index) );
- }
- else
- {
- m_RayIntersectionVoxels[0] =
- m_RayIntersectionVoxels[1] =
- m_RayIntersectionVoxels[2] =
- m_RayIntersectionVoxels[3] = NULL;
- }
- break;
- }
-
- case TRANSVERSE_IN_Y:
- {
-
- if( (Ix >= 0) && (Ix + 1 < m_NumberOfVoxelsInX) &&
- (Iy >= 0) && (Iy < m_NumberOfVoxelsInY) &&
- (Iz >= 0) && (Iz + 1 < m_NumberOfVoxelsInZ))
- {
-
- index[0]=Ix; index[1]=Iy; index[2]=Iz;
- m_RayIntersectionVoxels[0] = ( this->m_Image->GetBufferPointer()
- + this->m_Image->ComputeOffset(index) );
-
- index[0]=Ix+1; index[1]=Iy; index[2]=Iz;
- m_RayIntersectionVoxels[1] = ( this->m_Image->GetBufferPointer()
- + this->m_Image->ComputeOffset(index) );
-
- index[0]=Ix; index[1]=Iy; index[2]=Iz+1;
- m_RayIntersectionVoxels[2] = ( this->m_Image->GetBufferPointer()
- + this->m_Image->ComputeOffset(index) );
-
- index[0]=Ix+1; index[1]=Iy; index[2]=Iz+1;
- m_RayIntersectionVoxels[3] = ( this->m_Image->GetBufferPointer()
- + this->m_Image->ComputeOffset(index) );
- }
- else
- {
- m_RayIntersectionVoxels[0]
- = m_RayIntersectionVoxels[1]
- = m_RayIntersectionVoxels[2]
- = m_RayIntersectionVoxels[3] = NULL;
- }
- break;
- }
-
- case TRANSVERSE_IN_Z:
- {
-
- if( (Ix >= 0) && (Ix + 1 < m_NumberOfVoxelsInX) &&
- (Iy >= 0) && (Iy + 1 < m_NumberOfVoxelsInY) &&
- (Iz >= 0) && (Iz < m_NumberOfVoxelsInZ))
- {
-
- index[0]=Ix; index[1]=Iy; index[2]=Iz;
- m_RayIntersectionVoxels[0] = ( this->m_Image->GetBufferPointer()
- + this->m_Image->ComputeOffset(index) );
-
- index[0]=Ix+1; index[1]=Iy; index[2]=Iz;
- m_RayIntersectionVoxels[1] = ( this->m_Image->GetBufferPointer()
- + this->m_Image->ComputeOffset(index) );
-
- index[0]=Ix; index[1]=Iy+1; index[2]=Iz;
- m_RayIntersectionVoxels[2] = ( this->m_Image->GetBufferPointer()
- + this->m_Image->ComputeOffset(index) );
-
- index[0]=Ix+1; index[1]=Iy+1; index[2]=Iz;
- m_RayIntersectionVoxels[3] = ( this->m_Image->GetBufferPointer()
- + this->m_Image->ComputeOffset(index) );
-
- }
- else
- {
- m_RayIntersectionVoxels[0]
- = m_RayIntersectionVoxels[1]
- = m_RayIntersectionVoxels[2]
- = m_RayIntersectionVoxels[3] = NULL;
- }
- break;
- }
-
- default:
- {
- itk::ExceptionObject err(__FILE__, __LINE__);
- err.SetLocation( ITK_LOCATION );
- err.SetDescription( "The ray traversal direction is unset "
- "- InitialiseVoxelPointers().");
- throw err;
- return;
- }
- }
-
- }//JV end if deformable
- }
-
- /* -----------------------------------------------------------------------
- IncrementVoxelPointers() - Increment the voxel pointers
- ----------------------------------------------------------------------- */
-
- template<class TInputImage, class TCoordRep>
- void
- RayCastHelper<TInputImage, TCoordRep>
- ::IncrementVoxelPointers(void)
- {
-
-
- double xBefore = m_Position3Dvox[0];
- double yBefore = m_Position3Dvox[1];
- double zBefore = m_Position3Dvox[2];
-
- m_Position3Dvox[0] += m_VoxelIncrement[0];
- m_Position3Dvox[1] += m_VoxelIncrement[1];
- m_Position3Dvox[2] += m_VoxelIncrement[2];
-
-
- //JV don't do anything to pointers if deformable
- if(!m_TransformIsDeformable)
- {
- int dx = ((int) m_Position3Dvox[0]) - ((int) xBefore);
- int dy = ((int) m_Position3Dvox[1]) - ((int) yBefore);
- int dz = ((int) m_Position3Dvox[2]) - ((int) zBefore);
-
- m_RayIntersectionVoxelIndex[0] += dx;
- m_RayIntersectionVoxelIndex[1] += dy;
- m_RayIntersectionVoxelIndex[2] += dz;
-
- int totalRayVoxelPlanes
- = dx + dy*m_NumberOfVoxelsInX + dz*m_NumberOfVoxelsInX*m_NumberOfVoxelsInY;
-
- m_RayIntersectionVoxels[0] += totalRayVoxelPlanes;
- m_RayIntersectionVoxels[1] += totalRayVoxelPlanes;
- m_RayIntersectionVoxels[2] += totalRayVoxelPlanes;
- m_RayIntersectionVoxels[3] += totalRayVoxelPlanes;
- }//JV end if deformable
- }
-
-
- /* -----------------------------------------------------------------------
- GetCurrentIntensity() - Get the intensity of the current ray point.
- ----------------------------------------------------------------------- */
-
- template<class TInputImage, class TCoordRep>
- double
- RayCastHelper<TInputImage, TCoordRep>
- ::GetCurrentIntensity(void) const
- {
-
- //DD("In get current intensity");
- //JV fetch the interpolated value by applying the deformable transform
- if(m_TransformIsDeformable)
- {
- //3DPositionVox is in voxels AND shifted half a voxel in plane
- typename InputImageType::PointType inputPoint;
- // switch( m_TraversalDirection )
- // {
- // case TRANSVERSE_IN_X:
- // {
- // double x = (int)m_Position3Dvox[0] - 0.5*m_VoxelIncrement[0];
- // inputPoint[0]=x*m_VoxelDimensionInX;
-
- // inputPoint[1] -= ( (int(x)-x)*m_VoxelIncrement[1]*m_VoxelIncrement[0]
- // + 0.5*m_VoxelIncrement[1] - 0.5)*m_VoxelDimensionInY;
-
- // inputPoint[2] -= ( (int(x)-x)*m_VoxelIncrement[2]*m_VoxelIncrement[0]
- // + 0.5*m_VoxelIncrement[2] - 0.5)*m_VoxelDimensionInZ;
-
- // // inputPoint[0]=m_Position3Dvox[0]*m_VoxelDimensionInX;
- // // inputPoint[1]=(m_Position3Dvox[1]+0.5)*m_VoxelDimensionInY;
- // // inputPoint[2]=(m_Position3Dvox[2]+0.5)*m_VoxelDimensionInZ;
- // // DD("X");
- // break;
- // }
- // case TRANSVERSE_IN_Y:
- // {
- // double y = (int)m_Position3Dvox[1] - 0.5*m_VoxelIncrement[1];
- // double y=m_Position3Dvox[1];
- // inputPoint[1]=y*m_VoxelDimensionInY;
-
- // inputPoint[0] = (m_Position3Dvox[0]- ((int)y-y)*m_VoxelIncrement[0]*m_VoxelIncrement[1]
- // + 0.5*m_VoxelIncrement[0] - 0.5) * m_VoxelDimensionInX;
-
- // inputPoint[2] = (m_Position3Dvox[2]- ((int)y-y)*m_VoxelIncrement[2]*m_VoxelIncrement[1]
- // + 0.5*m_VoxelIncrement[2] - 0.5) * m_VoxelDimensionInZ;
-
- // // inputPoint[0]=(m_Position3Dvox[0]+0)*m_VoxelDimensionInX;
- // // inputPoint[1]=m_Position3Dvox[1]*m_VoxelDimensionInY;
- // // inputPoint[2]=(m_Position3Dvox[2]+0.5)*m_VoxelDimensionInZ;
-
- // break;
- // }
- // case TRANSVERSE_IN_Z:
- // {
- // double z = (int)m_Position3Dvox[2] - 0.5*m_VoxelIncrement[2];
- // inputPoint[2]=z*m_VoxelDimensionInZ;
-
- // inputPoint[0] += -( z*m_VoxelIncrement[1]*m_VoxelIncrement[0]
- // - 0.5*m_VoxelIncrement[2] + 0.5)*m_VoxelDimensionInZ;
-
- // inputPoint[1] += -( z*m_VoxelIncrement[1]*m_VoxelIncrement[0]
- // - 0.5*m_VoxelIncrement[2] + 0.5)*m_VoxelDimensionInZ;
-
- // // inputPoint[0]=(m_Position3Dvox[0]+0.5)*m_VoxelDimensionInX;
- // // inputPoint[1]=(m_Position3Dvox[1]+0.5)*m_VoxelDimensionInY;
- // // inputPoint[2]=m_Position3Dvox[2]*m_VoxelDimensionInZ;
- // // DD("Z");
- // break;
- // }
-
- // default:
- // {
- // itk::ExceptionObject err(__FILE__, __LINE__);
- // err.SetLocation( ITK_LOCATION );
- // err.SetDescription( "The ray traversal direction is unset "
- // "- GetCurrentIntensity().");
- // throw err;
- // return 0;
- // }
- // }
-
-
- //JV there seems to remain an small offset
- inputPoint[0]=(m_Position3Dvox[0])*m_VoxelDimensionInX;
- inputPoint[1]=(m_Position3Dvox[1])*m_VoxelDimensionInY;
- inputPoint[2]=m_Position3Dvox[2]*m_VoxelDimensionInZ;
- //JV work around, call Deformably transform point (no bulk)
- typedef clitk::BSplineDeformableTransform<double,3,3> DeformableTransformType;
- DeformableTransformType* bsplineTransform=dynamic_cast<DeformableTransformType*>(m_Transform);
- typename DeformableTransformType::OutputPointType transformedPoint=bsplineTransform->DeformablyTransformPoint(inputPoint);
-
- //check wether inside
- ContinuousIndexType contIndex;
- if (m_Image->TransformPhysicalPointToContinuousIndex(transformedPoint, contIndex))
- //interpolate this position in the input image
- return this->m_DeformableTransformInterpolator->EvaluateAtContinuousIndex(contIndex);
- else return 0;//m_EdgePaddingValue;
-
- }//JV end if deformable
-
- else
- {
- double a, b, c, d;
- double y, z;
-
- if (! m_ValidRay)
- {
- return 0;
- }
- a = (double) (*m_RayIntersectionVoxels[0]);
- b = (double) (*m_RayIntersectionVoxels[1] - a);
- c = (double) (*m_RayIntersectionVoxels[2] - a);
- d = (double) (*m_RayIntersectionVoxels[3] - a - b - c);
-
- switch( m_TraversalDirection )
- {
- case TRANSVERSE_IN_X:
- {
- y = m_Position3Dvox[1] - vcl_floor(m_Position3Dvox[1]);
- z = m_Position3Dvox[2] - vcl_floor(m_Position3Dvox[2]);
- break;
- }
- case TRANSVERSE_IN_Y:
- {
- y = m_Position3Dvox[0] - vcl_floor(m_Position3Dvox[0]);
- z = m_Position3Dvox[2] - vcl_floor(m_Position3Dvox[2]);
- break;
- }
- case TRANSVERSE_IN_Z:
- {
- y = m_Position3Dvox[0] - vcl_floor(m_Position3Dvox[0]);
- z = m_Position3Dvox[1] - vcl_floor(m_Position3Dvox[1]);
- break;
- }
- default:
- {
- itk::ExceptionObject err(__FILE__, __LINE__);
- err.SetLocation( ITK_LOCATION );
- err.SetDescription( "The ray traversal direction is unset "
- "- GetCurrentIntensity().");
- throw err;
- return 0;
- }
- }
-
- return a + b*y + c*z + d*y*z;
- }
- }
-
- /* -----------------------------------------------------------------------
- IncrementIntensities() - Increment the intensities of the current ray point
- ----------------------------------------------------------------------- */
-
- template<class TInputImage, class TCoordRep>
- void
- RayCastHelper<TInputImage, TCoordRep>
- ::IncrementIntensities(double increment)
- {
- short inc = (short) vcl_floor(increment + 0.5);
-
- if (! m_ValidRay)
- {
- return;
- }
- *m_RayIntersectionVoxels[0] += inc;
- *m_RayIntersectionVoxels[1] += inc;
- *m_RayIntersectionVoxels[2] += inc;
- *m_RayIntersectionVoxels[3] += inc;
-
- return;
- }
-
-
- /* -----------------------------------------------------------------------
- IntegrateAboveThreshold() - Integrate intensities above a threshold.
- ----------------------------------------------------------------------- */
-
- template<class TInputImage, class TCoordRep>
- bool
- RayCastHelper<TInputImage, TCoordRep>
- ::IntegrateAboveThreshold(double &integral, double threshold)
- {
- double intensity;
- // double posn3D_x, posn3D_y, posn3D_z;
-
- integral = 0.;
-
- // Check if this is a valid ray
-
- if (! m_ValidRay)
- {
- return false;
- }
- /* Step along the ray as quickly as possible
- integrating the interpolated intensities. */
-
- for (m_NumVoxelPlanesTraversed=0;
- m_NumVoxelPlanesTraversed<m_TotalRayVoxelPlanes;
- m_NumVoxelPlanesTraversed++)
- {
- intensity = this->GetCurrentIntensity();
-
- if (intensity > threshold)
- {
- integral += intensity - threshold;
- }
- this->IncrementVoxelPointers();
- }
-
- /* The ray passes through the volume one plane of voxels at a time,
- however, if its moving diagonally the ray points will be further
- apart so account for this by scaling by the distance moved. */
-
- integral *= this->GetRayPointSpacing();
-
- return true;
- }
-
- /* -----------------------------------------------------------------------
- ZeroState() - Set the default (zero) state of the object
- ----------------------------------------------------------------------- */
-
- template<class TInputImage, class TCoordRep>
- void
- RayCastHelper<TInputImage, TCoordRep>
- ::ZeroState()
- {
- int i;
-
- m_ValidRay = false;
-
- m_NumberOfVoxelsInX = 0;
- m_NumberOfVoxelsInY = 0;
- m_NumberOfVoxelsInZ = 0;
-
- m_VoxelDimensionInX = 0;
- m_VoxelDimensionInY = 0;
- m_VoxelDimensionInZ = 0;
-
- for (i=0; i<3; i++)
- {
- m_CurrentRayPositionInMM[i] = 0.;
- }
- for (i=0; i<3; i++)
- {
- m_RayDirectionInMM[i] = 0.;
- }
- for (i=0; i<3; i++)
- {
- m_RayVoxelStartPosition[i] = 0.;
- }
- for (i=0; i<3; i++)
- {
- m_RayVoxelEndPosition[i] = 0.;
- }
- for (i=0; i<3; i++)
- {
- m_VoxelIncrement[i] = 0.;
- }
- m_TraversalDirection = UNDEFINED_DIRECTION;
-
- m_TotalRayVoxelPlanes = 0;
- m_NumVoxelPlanesTraversed = -1;
-
- for (i=0; i<4; i++)
- {
- m_RayIntersectionVoxels[i] = 0;
- }
- for (i=0; i<3; i++)
- {
- m_RayIntersectionVoxelIndex[i] = 0;
- }
- }
-}; // end of anonymous namespace
-
-
-namespace clitk
-{
-
- /**************************************************************************
- *
- *
- * Rest of this code is the actual RayCastInterpolateImageFunctionWithOrigin
- * class
- *
- *
- **************************************************************************/
-
- /* -----------------------------------------------------------------------
- Constructor
- ----------------------------------------------------------------------- */
-
- template<class TInputImage, class TCoordRep>
- RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
- ::RayCastInterpolateImageFunctionWithOrigin()
- {
- m_Threshold = 0.;
-
- m_FocalPoint[0] = 0.;
- m_FocalPoint[1] = 0.;
- m_FocalPoint[2] = 0.;
-
- //JV
- m_TransformIsDeformable=false;
- }
-
-
- /* -----------------------------------------------------------------------
- PrintSelf
- ----------------------------------------------------------------------- */
-
- template<class TInputImage, class TCoordRep>
- void
- RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
- ::PrintSelf(std::ostream& os, itk::Indent indent) const
- {
- this->Superclass::PrintSelf(os,indent);
-
- os << indent << "Threshold: " << m_Threshold << std::endl;
- os << indent << "FocalPoint: " << m_FocalPoint << std::endl;
- os << indent << "Transform: " << m_Transform.GetPointer() << std::endl;
- os << indent << "Interpolator: " << m_Interpolator.GetPointer() << std::endl;
-
- }
-
- /* -----------------------------------------------------------------------
- Evaluate at image index position
- ----------------------------------------------------------------------- */
-
- template<class TInputImage, class TCoordRep>
- typename RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
- ::OutputType
- RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
- ::Evaluate( const PointType& point ) const
- {
- double integral = 0;
-
- OutputPointType transformedFocalPoint
- = m_Transform->TransformPoint( m_FocalPoint );
-
- DirectionType direction = transformedFocalPoint - point;
-
- RayCastHelper<TInputImage, TCoordRep> ray;
-
- //JV
- ray.SetTransformIsDeformable(m_TransformIsDeformable);
- if (m_TransformIsDeformable)
- {
- ray.SetTransform(m_Transform);
- ray.SetDeformableTransformInterpolator(m_DeformableTransformInterpolator);
- }
- //JV
-
- ray.SetImage( this->m_Image );
- ray.ZeroState();
-
- ray.Initialise();
- // (Aviv) Added support for images with non-zero origin
- // ray.SetRay(point, direction);
- ray.SetRay(point - this->m_Image->GetOrigin().GetVectorFromOrigin(), direction);
- ray.IntegrateAboveThreshold(integral, m_Threshold);
-
- return ( static_cast<OutputType>( integral ));
- }
-
- template<class TInputImage, class TCoordRep>
- typename RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
- ::OutputType
- RayCastInterpolateImageFunctionWithOrigin< TInputImage, TCoordRep >
- ::EvaluateAtContinuousIndex( const ContinuousIndexType& index ) const
- {
- OutputPointType point;
- this->m_Image->TransformContinuousIndexToPhysicalPoint(index, point);
-
- return this->Evaluate( point );
- }
-
-} // namespace clitk
-
-
-#endif
git://git.creatis.insa-lyon.fr / clitk.git / commitdiff