From: srit Date: Mon, 6 Dec 2010 15:40:02 +0000 (+0000) Subject: Removed unfinished work of Jef X-Git-Tag: v1.2.0~295 X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=commitdiff_plain;h=ce530f00de56df065ed5ce7155c75e9ba17cea44;p=clitk.git Removed unfinished work of Jef --- diff --git a/itk/clitkRayCastInterpolateImageFunctionWithOrigin.h b/itk/clitkRayCastInterpolateImageFunctionWithOrigin.h deleted file mode 100755 index 32f322e..0000000 --- a/itk/clitkRayCastInterpolateImageFunctionWithOrigin.h +++ /dev/null @@ -1,215 +0,0 @@ -/*========================================================================= - 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_h -#define __clitkRayCastInterpolateImageFunctionWithOrigin_h -#include "clitkBSplineDeformableTransform.h" - -#include "itkInterpolateImageFunction.h" -#include "itkTransform.h" -#include "itkVector.h" - -namespace clitk -{ - - -template -class ITK_EXPORT RayCastInterpolateImageFunctionWithOrigin : - public itk::InterpolateImageFunction -{ -public: - /** Standard class typedefs. */ - typedef RayCastInterpolateImageFunctionWithOrigin Self; - typedef itk::InterpolateImageFunction Superclass; - typedef itk::SmartPointer Pointer; - typedef itk::SmartPointer ConstPointer; - - /** 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 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 Superclass::InputPixelType PixelType; - - typedef typename TInputImage::SizeType SizeType; - - typedef itk::Vector DirectionType; - - /** Type of the Interpolator Base class */ - typedef itk::InterpolateImageFunction InterpolatorType; - - typedef typename InterpolatorType::Pointer InterpolatorPointer; - - //JV Type of the Interpolator for the transform - typedef itk::InterpolateImageFunction DeformableTransformInterpolatorType; - typedef typename DeformableTransformInterpolatorType::Pointer DeformableTransformInterpolatorPointer; - - - /** Run-time type information (and related methods). */ - itkTypeMacro(RayCastInterpolateImageFunctionWithOrigin, InterpolateImageFunction); - - /** Method for creation through the object factory. */ - itkNewMacro(Self); - - /** OutputType typedef support. */ - typedef typename Superclass::OutputType OutputType; - - /** InputImageType typedef support. */ - typedef typename Superclass::InputImageType InputImageType; - - /** RealType typedef support. */ - typedef typename Superclass::RealType RealType; - - /** Dimension underlying input image. */ - itkStaticConstMacro(ImageDimension, unsigned int,Superclass::ImageDimension); - - /** Point typedef support. */ - typedef typename Superclass::PointType PointType; - - /** Index typedef support. */ - typedef typename Superclass::IndexType IndexType; - - /** ContinuousIndex typedef support. */ - typedef typename Superclass::ContinuousIndexType ContinuousIndexType; - - /** \brief - * Interpolate the image at a point position. - * - * Returns the interpolated image intensity at a - * specified point position. No bounds checking is done. - * The point is assume to lie within the image buffer. - * - * ImageFunction::IsInsideBuffer() can be used to check bounds before - * calling the method. - */ - virtual OutputType Evaluate( const PointType& point ) const; - - /** Interpolate the image at a continuous index position - * - * Returns the interpolated image intensity at a - * specified index position. No bounds checking is done. - * The point is assume to lie within the image buffer. - * - * Subclasses must override this method. - * - * ImageFunction::IsInsideBuffer() can be used to check bounds before - * calling the method. - */ - virtual OutputType EvaluateAtContinuousIndex( - const ContinuousIndexType &index ) const; - - - /** Connect the Transform. */ - itkSetObjectMacro( Transform, TransformType ); - /** Get a pointer to the Transform. */ - itkGetObjectMacro( Transform, TransformType ); - - /** Connect the Interpolator. */ - itkSetObjectMacro( Interpolator, InterpolatorType ); - /** Get a pointer to the Interpolator. */ - itkGetObjectMacro( Interpolator, InterpolatorType ); - - /** Connect the focalPoint. */ - itkSetMacro( FocalPoint, InputPointType ); - /** Get a pointer to the focalpoint. */ - itkGetMacro( FocalPoint, InputPointType ); - - /** Connect the threshold. */ - itkSetMacro( Threshold, double ); - /** Get a pointer to the threshhold. */ - itkGetMacro( Threshold, double ); - - /** Check if a point is inside the image buffer. - * \warning For efficiency, no validity checking of - * the input image pointer is done. */ - inline bool IsInsideBuffer( const PointType & ) const - { - return true; - } - bool IsInsideBuffer( const ContinuousIndexType & ) const - { - return true; - } - bool IsInsideBuffer( const IndexType & ) const - { - return true; - } - - //JV transform is deformable - itkSetMacro( TransformIsDeformable, bool ); - itkGetMacro( TransformIsDeformable, bool ); - - //JV interpolator for deformable transform - itkSetObjectMacro( DeformableTransformInterpolator, DeformableTransformInterpolatorType ); - itkGetObjectMacro( DeformableTransformInterpolator, DeformableTransformInterpolatorType ); - -protected: - - /// Constructor - RayCastInterpolateImageFunctionWithOrigin(); - - /// Destructor - ~RayCastInterpolateImageFunctionWithOrigin(){}; - - /// Print the object - void PrintSelf(std::ostream& os, itk::Indent indent) const; - - /// Transformation used to calculate the new focal point position - TransformPointer m_Transform; - - /// The focal point or position of the ray source - InputPointType m_FocalPoint; - - /// The threshold above which voxels along the ray path are integrated. - double m_Threshold; - - /// Pointer to the interpolator - InterpolatorPointer m_Interpolator; - - //JV Pointer to the deformable transform interpolator - DeformableTransformInterpolatorPointer m_DeformableTransformInterpolator; - - - //JV Transform is deformable: fetch interpolated value at transformed intersect - bool m_TransformIsDeformable; - - -private: - RayCastInterpolateImageFunctionWithOrigin( const Self& ); //purposely not implemented - void operator=( const Self& ); //purposely not implemented - - -}; - -} // namespace clitk - -#ifndef ITK_MANUAL_INSTANTIATION -#include "clitkRayCastInterpolateImageFunctionWithOrigin.txx" -#endif - -#endif diff --git a/itk/clitkRayCastInterpolateImageFunctionWithOrigin.txx b/itk/clitkRayCastInterpolateImageFunctionWithOrigin.txx deleted file mode 100755 index 92cb8db..0000000 --- a/itk/clitkRayCastInterpolateImageFunctionWithOrigin.txx +++ /dev/null @@ -1,1706 +0,0 @@ -/*========================================================================= - 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 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 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 DirectionType; - typedef itk::Point PointType; - - typedef TInputImage InputImageType; - typedef typename InputImageType::PixelType PixelType; - typedef typename InputImageType::IndexType IndexType; - typedef itk::ContinuousIndex ContinuousIndexType; - - - //JV Add an interpolator for the deformable transform - typedef itk::InterpolateImageFunction 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 - void - RayCastHelper - ::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 - void - RayCastHelper - ::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 - void - RayCastHelper - ::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 - void - RayCastHelper - ::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 - bool - RayCastHelper - ::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=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 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 - bool - RayCastHelper - ::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 - void - RayCastHelper - ::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 - void - RayCastHelper - ::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 - bool - RayCastHelper - ::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 - void - RayCastHelper - ::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 - void - RayCastHelper - ::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 - void - RayCastHelper - ::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 - double - RayCastHelper - ::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 DeformableTransformType; - DeformableTransformType* bsplineTransform=dynamic_cast(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 - void - RayCastHelper - ::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 - bool - RayCastHelper - ::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_NumVoxelPlanesTraversedGetCurrentIntensity(); - - 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 - void - RayCastHelper - ::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 - 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 - 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 - 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 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( integral )); - } - - template - 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