fattori_g/reg23Topograms
0
0
Fork 0
Code Issues 1 Pull Requests Actions Packages Projects Releases Wiki Activity

Siddon

Browse Source 
revised code to correct ray calculation that uses boundaries instead of voxel centers as we were expecting
This commit is contained in:
Proton local user
2023-06-06 10:48:06 +02:00
parent 6f3a5d70b6
commit 90cfe3defd
2 changed files with 422 additions and 427 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
reg23Topograms/itkDTRrecon/itkgSiddonJacobsRayCastInterpolateImageFunction.h
View File

@ -131,6 +131,9 @@ public:
using ContinuousIndexType = typename Superclass::ContinuousIndexType;
/** Internal floating point comparison accuracy **/
static const double EPSILON;
/** \brief
* Interpolate the image at a point position.
*

846
reg23Topograms/itkDTRrecon/itkgSiddonJacobsRayCastInterpolateImageFunction.hxx
View File

@ -56,40 +56,45 @@ CIT 6, 89-94 (1998).
namespace itk
{
template <typename TInputImage, typename TCoordRep>
const double gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::EPSILON =
0.00001;
template <typename TInputImage, typename TCoordRep>
gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::gSiddonJacobsRayCastInterpolateImageFunction()
{
m_FocalPointToIsocenterDistance = 1000.; // Focal point to isocenter distance in mm.
m_ProjectionAngle = 0.; // Angle in radians betweeen projection central axis and reference axis
m_Threshold = 0.; // Intensity threshold, below which is ignored.
m_PanelOffset = 0.;
m_FocalPointToIsocenterDistance = 1000.; // Focal point to isocenter distance in mm.
m_ProjectionAngle = 0.; // Angle in radians betweeen projection central axis and reference axis
m_Threshold = 0.; // Intensity threshold, below which is ignored.
m_PanelOffset = 0.;
m_SourcePoint[0] = 0.;
m_SourcePoint[1] = 0.;
m_SourcePoint[2] = 0.;
m_SourcePoint[0] = 0.;
m_SourcePoint[1] = 0.;
m_SourcePoint[2] = 0.;
m_InverseTransform = TransformType::New();
m_InverseTransform->SetComputeZYX(true);
m_InverseTransform = TransformType::New();
m_InverseTransform->SetComputeZYX(true);
m_ComposedTransform = TransformType::New();
m_ComposedTransform->SetComputeZYX(true);
m_ComposedTransform = TransformType::New();
m_ComposedTransform->SetComputeZYX(true);
m_GantryRotTransform = TransformType::New();
m_GantryRotTransform->SetComputeZYX(true);
m_GantryRotTransform->SetIdentity();
m_GantryRotTransform = TransformType::New();
m_GantryRotTransform->SetComputeZYX(true);
m_GantryRotTransform->SetIdentity();
m_CamShiftTransform = TransformType::New();
m_CamShiftTransform->SetComputeZYX(true);
m_CamShiftTransform->SetIdentity();
m_CamShiftTransform = TransformType::New();
m_CamShiftTransform->SetComputeZYX(true);
m_CamShiftTransform->SetIdentity();
m_CamRotTransform = TransformType::New();
m_CamRotTransform->SetComputeZYX(true);
m_CamRotTransform->SetIdentity();
// constant for converting degrees into radians
const float dtr = (atan(1.0) * 4.0) / 180.0;
m_CamRotTransform->SetRotation(dtr * (-90.0), 0.0, 0.0);
m_CamRotTransform = TransformType::New();
m_CamRotTransform->SetComputeZYX(true);
m_CamRotTransform->SetIdentity();
// constant for converting degrees into radians
const float dtr = (atan(1.0) * 4.0) / 180.0;
m_CamRotTransform->SetRotation(dtr * (-90.0), 0.0, 0.0);
m_Threshold = 0;
m_Threshold = 0;
}
@ -97,10 +102,10 @@ template <typename TInputImage, typename TCoordRep>
void
gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::PrintSelf(std::ostream & os, Indent indent) const
{
this->Superclass::PrintSelf(os, indent);
this->Superclass::PrintSelf(os, indent);
os << indent << "Threshold: " << m_Threshold << std::endl;
os << indent << "Transform: " << m_Transform.GetPointer() << std::endl;
os << indent << "Threshold: " << m_Threshold << std::endl;
os << indent << "Transform: " << m_Transform.GetPointer() << std::endl;
}
@ -108,410 +113,388 @@ template <typename TInputImage, typename TCoordRep>
typename gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::OutputType
gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(const PointType & point) const
{
float rayVector[3];
IndexType cIndex;
float rayVector[3];
IndexType cIndex;
PointType drrPixelWorld; // Coordinate of a DRR pixel in the world coordinate system
OutputType pixval;
PointType drrPixelWorld; // Coordinate of a DRR pixel in the world coordinate system
OutputType pixval;
float firstIntersection[3];
float alphaX1, alphaXN, alphaXmin, alphaXmax;
float alphaY1, alphaYN, alphaYmin, alphaYmax;
float alphaZ1, alphaZN, alphaZmin, alphaZmax;
float alphaMin, alphaMax;
float alphaX, alphaY, alphaZ, alphaCmin, alphaCminPrev;
float alphaUx, alphaUy, alphaUz;
float alphaIntersectionUp[3], alphaIntersectionDown[3];
float d12, value;
float firstIntersectionIndex[3];
int firstIntersectionIndexUp[3], firstIntersectionIndexDown[3];
int iU, jU, kU;
float firstIntersection[3];
float alphaX1, alphaXN, alphaXmin, alphaXmax;
float alphaY1, alphaYN, alphaYmin, alphaYmax;
float alphaZ1, alphaZN, alphaZmin, alphaZmax;
float alphaMin, alphaMax;
float alphaX, alphaY, alphaZ, alphaCmin, alphaCminPrev;
float alphaUx, alphaUy, alphaUz;
float alphaIntersectionUp[3], alphaIntersectionDown[3];
float d12, value;
float firstIntersectionIndex[3];
int firstIntersectionIndexUp[3], firstIntersectionIndexDown[3];
int iU, jU, kU;
// Min/max values of the output pixel type AND these values
// represented as the output type of the interpolator
const OutputType minOutputValue = itk::NumericTraits<OutputType>::NonpositiveMin();
const OutputType maxOutputValue = itk::NumericTraits<OutputType>::max();
// Min/max values of the output pixel type AND these values
// represented as the output type of the interpolator
const OutputType minOutputValue = itk::NumericTraits<OutputType>::NonpositiveMin();
const OutputType maxOutputValue = itk::NumericTraits<OutputType>::max();
// If the volume was shifted, recalculate the overall inverse transform
unsigned long int interpMTime = this->GetMTime();
unsigned long int vTransformMTime = m_Transform->GetMTime();
// If the volume was shifted, recalculate the overall inverse transform
unsigned long int interpMTime = this->GetMTime();
unsigned long int vTransformMTime = m_Transform->GetMTime();
if (interpMTime < vTransformMTime)
{
this->ComputeInverseTransform();
// The m_SourceWorld should be computed here to avoid the repeatedly calculation
// for each projection ray. However, we are in a const function, which prohibits
// the modification of class member variables. So the world coordiate of the source
// point is calculated for each ray as below. Performance improvement may be made
// by using a static variable?
// m_SourceWorld = m_InverseTransform->TransformPoint(m_SourcePoint);
}
if (interpMTime < vTransformMTime)
{
this->ComputeInverseTransform();
// The m_SourceWorld should be computed here to avoid the repeatedly calculation
// for each projection ray. However, we are in a const function, which prohibits
// the modification of class member variables. So the world coordiate of the source
// point is calculated for each ray as below. Performance improvement may be made
// by using a static variable?
// m_SourceWorld = m_InverseTransform->TransformPoint(m_SourcePoint);
}
PointType PointReq = point;
PointReq[0] += m_PanelOffset;
PointType PointReq = point;
PointReq[0] += m_PanelOffset;
drrPixelWorld = m_InverseTransform->TransformPoint(PointReq);
drrPixelWorld = m_InverseTransform->TransformPoint(PointReq);
// Get ths input pointers
InputImageConstPointer inputPtr = this->GetInputImage();
typename InputImageType::SizeType sizeCT;
typename InputImageType::RegionType regionCT;
typename InputImageType::SpacingType ctPixelSpacing;
typename InputImageType::PointType ctOrigin;
ctPixelSpacing = inputPtr->GetSpacing();
ctOrigin = inputPtr->GetOrigin();
regionCT = inputPtr->GetLargestPossibleRegion();
sizeCT = regionCT.GetSize();
PointType SlidingSourcePoint = m_SourcePoint;
SlidingSourcePoint[0] = 0.;
// Simple sliding source model, the source follows the y of the pixel
// requested, which sets the Z position to look in the CT volume.
SlidingSourcePoint[1] = point[1];
SlidingSourcePoint[2] = 0.;
// Get ths input pointers
InputImageConstPointer inputPtr = this->GetInputImage();
PointType SourceWorld = m_InverseTransform->TransformPoint(SlidingSourcePoint);
typename InputImageType::SizeType sizeCT;
typename InputImageType::RegionType regionCT;
typename InputImageType::SpacingType ctPixelSpacing;
typename InputImageType::PointType ctOrigin;
PointType O(3);
O[0] = -ctPixelSpacing[0]/2.;
O[1] = -ctPixelSpacing[1]/2.;
O[2] = -ctPixelSpacing[2]/2.;
ctPixelSpacing = inputPtr->GetSpacing();
ctOrigin = inputPtr->GetOrigin();
regionCT = inputPtr->GetLargestPossibleRegion();
sizeCT = regionCT.GetSize();
// be sure that drr pixel position is inside the volume
float xSizeCT = O[0] + (sizeCT[0]) * ctPixelSpacing[0];
float ySizeCT = O[1] + (sizeCT[1]) * ctPixelSpacing[1];
float zSizeCT = O[2] + (sizeCT[2]) * ctPixelSpacing[2];
float xDrrPix = drrPixelWorld[0];
float yDrrPix = drrPixelWorld[1];
float zDrrPix = drrPixelWorld[2];
if(
xDrrPix<= O[0] || xDrrPix>=xSizeCT ||
yDrrPix<= O[1] || yDrrPix>=ySizeCT ||
zDrrPix<= O[2] || zDrrPix>=zSizeCT
){
pixval = static_cast<OutputType>(0.0);
return pixval;
}
// calculate the detector position for this pixel center by moving
// 2*m_FocalPointToIsocenterDistance from the source in the pixel
// directions
double p1[3], p2[3];
p1[0] = SourceWorld[0];
p1[1] = SourceWorld[1];
p1[2] = SourceWorld[2];
p2[0] = drrPixelWorld[0];
p2[1] = drrPixelWorld[1];
p2[2] = drrPixelWorld[2];
double P12D = sqrt(
(p2[0]-p1[0]) * (p2[0]-p1[0]) +
(p2[1]-p1[1]) * (p2[1]-p1[1]) +
(p2[2]-p1[2]) * (p2[2]-p1[2])
);
double p12V[3];
p12V[0] = p2[0] - p1[0];
p12V[1] = p2[1] - p1[1];
p12V[2] = p2[2] - p1[2];
double p12v [3];
p12v[0] = p12V[0]/P12D;
p12v[1] = p12V[1]/P12D;
p12v[2] = p12V[2]/P12D;
double pDest [3];
pDest[0] = p1[0] + m_FocalPointToIsocenterDistance * 2 * p12v[0];
pDest[1] = p1[1] + m_FocalPointToIsocenterDistance * 2 * p12v[1];
pDest[2] = p1[2] + m_FocalPointToIsocenterDistance * 2 * p12v[2];
PointType SlidingSourcePoint = m_SourcePoint;
SlidingSourcePoint[0] = 0.;
SlidingSourcePoint[1] = point[1] /*- ctPixelSpacing[2]/2.*/;
SlidingSourcePoint[2] = 0.;
PointType SourceWorld = m_InverseTransform->TransformPoint(SlidingSourcePoint);
//std::cout<<"SourceWorld: "<<SourceWorld<<std::endl;
// If Pixel position (in mm) is outside bounds of CT (zero-based)
// assign 0 at the pixel and move on. Ensures regular spacing of resulting
// DRR
float xSizeCT = -ctPixelSpacing[0]/2. + (sizeCT[0]) * ctPixelSpacing[0];
float ySizeCT = -ctPixelSpacing[1]/2. + (sizeCT[1]) * ctPixelSpacing[1];
float zSizeCT = -ctPixelSpacing[2]/2. + (sizeCT[2]) * ctPixelSpacing[2];
float xDrrPix = drrPixelWorld[0];
float yDrrPix = drrPixelWorld[1];
float zDrrPix = drrPixelWorld[2];
// if(zDrrPix < 0 /*|| yDrrPix < 0 || xDrrPix < 0*/)
// std::cout << drrPixelWorld[0]<<" " <<drrPixelWorld[1]<<" "<<drrPixelWorld[2]<<" / ";
if(xDrrPix>xSizeCT || yDrrPix>ySizeCT || zDrrPix>zSizeCT){
pixval = static_cast<OutputType>(0.0);
return pixval;
}
// std::cout<<point << " " <<PointReq <<" "<<drrPixelWorld<<std::endl;
// std::cout<< ctOrigin <<std::endl;
// std::cout<< m_Transform->GetCenter()<<std::endl;
// std::cout<<sizeCT << " " <<ctPixelSpacing <<" "<<drrPixelWorld<<std::endl;
// The following is the Siddon-Jacob fast ray-tracing algorithm
//rayVector[0] = drrPixelWorld[0] - SourceWorld[0];
//rayVector[1] = drrPixelWorld[1] - SourceWorld[1];
// correct ray to reach detector for topogram geometry
// Differs from conventional FPD geometry (above) that has
// the panel at the end of the ray.
rayVector[0] = pDest[0] - SourceWorld[0];
rayVector[1] = pDest[1] - SourceWorld[1];
rayVector[2] = pDest[2] - SourceWorld[2];
// exit(1);
// calculate the detector position for this pixel center by moving
// 2*m_FocalPointToIsocenterDistance from the source in the pixel
// directions
double p1[3], p2[3];
p1[0] = SourceWorld[0];
p1[1] = SourceWorld[1];
p1[2] = SourceWorld[2];
p2[0] = drrPixelWorld[0];
p2[1] = drrPixelWorld[1];
p2[2] = drrPixelWorld[2];
double P12D = sqrt(
(p2[0]-p1[0]) * (p2[0]-p1[0]) +
(p2[1]-p1[1]) * (p2[1]-p1[1]) +
(p2[2]-p1[2]) * (p2[2]-p1[2])
);
double p12V[3];
p12V[0] = p2[0] - p1[0];
p12V[1] = p2[1] - p1[1];
p12V[2] = p2[2] - p1[2];
double p12v [3];
p12v[0] = p12V[0]/P12D;
p12v[1] = p12V[1]/P12D;
p12v[2] = p12V[2]/P12D;
double pDest [3];
pDest[0] = p1[0] + m_FocalPointToIsocenterDistance * 2 * p12v[0];
pDest[1] = p1[1] + m_FocalPointToIsocenterDistance * 2 * p12v[1];
pDest[2] = p1[2] + m_FocalPointToIsocenterDistance * 2 * p12v[2];
// The following is the Siddon-Jacob fast ray-tracing algorithm
//rayVector[0] = drrPixelWorld[0] - SourceWorld[0];
//rayVector[1] = drrPixelWorld[1] - SourceWorld[1];
// correct ray to reach detector
rayVector[0] = pDest[0] - SourceWorld[0];
rayVector[1] = pDest[1] - SourceWorld[1];
rayVector[2] = pDest[2] - SourceWorld[2];
/* Calculate the parametric values of the first and the last
/* Calculate the parametric values of the first and the last
intersection points of the ray with the X, Y, and Z-planes that
define the CT volume. */
if (rayVector[0] != 0)
{
// alphaX1 = (0.0 - SourceWorld[0]) / rayVector[0];
// alphaXN = (sizeCT[0] * ctPixelSpacing[0] - SourceWorld[0]) / rayVector[0];
alphaX1 = ( (-ctPixelSpacing[0]/2.) - SourceWorld[0]) / rayVector[0];
alphaXN = ( (-ctPixelSpacing[0]/2.) + sizeCT[0] * ctPixelSpacing[0] - SourceWorld[0]) / rayVector[0];
alphaXmin = std::min(alphaX1, alphaXN);
alphaXmax = std::max(alphaX1, alphaXN);
}
else
{
alphaXmin = -2;
alphaXmax = 2;
}
if (rayVector[1] != 0)
{
alphaY1 = ( (-ctPixelSpacing[1]/2.) - SourceWorld[1]) / rayVector[1];
alphaYN = ( -ctPixelSpacing[1]/2. + sizeCT[1] * ctPixelSpacing[1] - SourceWorld[1]) / rayVector[1];
// alphaY1 = (0.0 - SourceWorld[1]) / rayVector[1];
// alphaYN = (sizeCT[1] * ctPixelSpacing[1] - SourceWorld[1]) / rayVector[1];
alphaYmin = std::min(alphaY1, alphaYN);
alphaYmax = std::max(alphaY1, alphaYN);
}
else
{
alphaYmin = -2;
alphaYmax = 2;
}
if (rayVector[2] != 0)
{
// alphaZ1 = (0.0 - SourceWorld[2]) / rayVector[2];
// alphaZN = (sizeCT[2] * ctPixelSpacing[2] - SourceWorld[2]) / rayVector[2];
alphaZ1 = ( (-ctPixelSpacing[2]/2.) - SourceWorld[2]) / rayVector[2];
alphaZN = ( -ctPixelSpacing[2]/2. + sizeCT[2] * ctPixelSpacing[2] - SourceWorld[2]) / rayVector[2];
alphaZmin = std::min(alphaZ1, alphaZN);
alphaZmax = std::max(alphaZ1, alphaZN);
}
else
{
alphaZmin = -2;
alphaZmax = 2;
}
/* Get the very first and the last alpha values when the ray
intersects with the CT volume. */
alphaMin = std::max(std::max(alphaXmin, alphaYmin), alphaZmin);
alphaMax = std::min(std::min(alphaXmax, alphaYmax), alphaZmax);
/* Calculate the parametric values of the first intersection point
of the ray with the X, Y, and Z-planes after the ray entered the
CT volume. */
firstIntersection[0] = SourceWorld[0] + alphaMin * rayVector[0] - ctPixelSpacing[0]/2.;
firstIntersection[1] = SourceWorld[1] + alphaMin * rayVector[1] - ctPixelSpacing[1]/2.;
firstIntersection[2] = SourceWorld[2] + alphaMin * rayVector[2] - ctPixelSpacing[2]/2.;
/* Transform world coordinate to the continuous index of the CT volume*/
firstIntersectionIndex[0] = firstIntersection[0] / ctPixelSpacing[0];
firstIntersectionIndex[1] = firstIntersection[1] / ctPixelSpacing[1];
firstIntersectionIndex[2] = firstIntersection[2] / ctPixelSpacing[2];
// now "correct" the indices (there are special cases):
if (alphaMin == alphaYmin) // j is special
{
if (alphaYmin == alphaY1)
firstIntersectionIndex [1] = 0;
else
// alpha_y_Ny
firstIntersectionIndex [1] = sizeCT[1] - 1;
}
else if (alphaMin == alphaXmin) // i is special
{
if (alphaXmin == alphaX1)
firstIntersectionIndex[0] = 0;
else
// alpha_x_Nx
firstIntersectionIndex[0] = sizeCT[0] - 1;
}
else if (alphaMin == alphaZmin) // k is special
{
if (alphaZmin == alphaZ1)
firstIntersectionIndex [2] = 0;
else
// alpha_z_Nz
firstIntersectionIndex [2] = sizeCT[2] - 1;
}
firstIntersectionIndexUp[0] = (int)ceil(firstIntersectionIndex[0]);
firstIntersectionIndexUp[1] = (int)ceil(firstIntersectionIndex[1]);
firstIntersectionIndexUp[2] = (int)ceil(firstIntersectionIndex[2]);
firstIntersectionIndexDown[0] = (int)floor(firstIntersectionIndex[0]);
firstIntersectionIndexDown[1] = (int)floor(firstIntersectionIndex[1]);
firstIntersectionIndexDown[2] = (int)floor(firstIntersectionIndex[2]);
if (rayVector[0] == 0)
{
alphaX = 2;
}
else
{
alphaIntersectionUp[0] = (firstIntersectionIndexUp[0] * ctPixelSpacing[0] -(ctPixelSpacing[0]/2.) - SourceWorld[0]) / rayVector[0];
alphaIntersectionDown[0] = (firstIntersectionIndexDown[0] * ctPixelSpacing[0] -(ctPixelSpacing[0]/2.) - SourceWorld[0]) / rayVector[0];
alphaX = std::max(alphaIntersectionUp[0], alphaIntersectionDown[0]);
}
if (rayVector[1] == 0)
{
alphaY = 2;
}
else
{
alphaIntersectionUp[1] = (firstIntersectionIndexUp[1] * ctPixelSpacing[1] -(ctPixelSpacing[1]/2.) - SourceWorld[1]) / rayVector[1];
alphaIntersectionDown[1] = (firstIntersectionIndexDown[1] * ctPixelSpacing[1] -(ctPixelSpacing[1]/2.) - SourceWorld[1]) / rayVector[1];
alphaY = std::max(alphaIntersectionUp[1], alphaIntersectionDown[1]);
}
if (rayVector[2] == 0)
{
alphaZ = 2;
}
else
{
alphaIntersectionUp[2] = (firstIntersectionIndexUp[2] * ctPixelSpacing[2] -(ctPixelSpacing[2]/2.) - SourceWorld[2]) / rayVector[2];
alphaIntersectionDown[2] = (firstIntersectionIndexDown[2] * ctPixelSpacing[2] -(ctPixelSpacing[2]/2.) - SourceWorld[2]) / rayVector[2];
alphaZ = std::max(alphaIntersectionUp[2], alphaIntersectionDown[2]);
}
/* Calculate alpha incremental values when the ray intercepts with x, y, and z-planes */
if (rayVector[0] != 0)
{
alphaUx = ctPixelSpacing[0] / fabs(rayVector[0]);
}
else
{
alphaUx = 999;
}
if (rayVector[1] != 0)
{
alphaUy = ctPixelSpacing[1] / fabs(rayVector[1]);
}
else
{
alphaUy = 999;
}
if (rayVector[2] != 0)
{
alphaUz = ctPixelSpacing[2] / fabs(rayVector[2]);
}
else
{
alphaUz = 999;
}
/* Calculate voxel index incremental values along the ray path. */
if (SourceWorld[0] < drrPixelWorld[0])
{
iU = 1;
}
else
{
iU = -1;
}
if (SourceWorld[1] < drrPixelWorld[1])
{
jU = 1;
}
else
{
jU = -1;
}
if (SourceWorld[2] < drrPixelWorld[2])
{
kU = 1;
}
else
{
kU = -1;
}
d12 = 0.0; /* Initialize the sum of the voxel intensities along the ray path to zero. */
/* Initialize the current ray position. */
alphaCmin = std::min(std::min(alphaX, alphaY), alphaZ);
/* Initialize the current voxel index. */
cIndex[0] = firstIntersectionIndexDown[0];
cIndex[1] = firstIntersectionIndexDown[1];
cIndex[2] = firstIntersectionIndexDown[2];
while (alphaCmin < alphaMax) /* Check if the ray is still in the CT volume */
{
/* Store the current ray position */
alphaCminPrev = alphaCmin;
if ((alphaX <= alphaY) && (alphaX <= alphaZ))
if (fabs(rayVector[0]) > EPSILON/* != 0*/)
{
/* Current ray front intercepts with x-plane. Update alphaX. */
alphaCmin = alphaX;
cIndex[0] = cIndex[0] + iU;
alphaX = alphaX + alphaUx;
}
else if ((alphaY <= alphaX) && (alphaY <= alphaZ))
{
/* Current ray front intercepts with y-plane. Update alphaY. */
alphaCmin = alphaY;
cIndex[1] = cIndex[1] + jU;
alphaY = alphaY + alphaUy;
// alphaX1 = (0.0 - SourceWorld[0]) / rayVector[0];
// alphaXN = (sizeCT[0] * ctPixelSpacing[0] - SourceWorld[0]) / rayVector[0];
alphaX1 = ( O[0] - SourceWorld[0]) / rayVector[0];
alphaXN = ( O[0] + sizeCT[0] * ctPixelSpacing[0] - SourceWorld[0]) / rayVector[0];
alphaXmin = std::min(alphaX1, alphaXN);
alphaXmax = std::max(alphaX1, alphaXN);
}
else
{
/* Current ray front intercepts with z-plane. Update alphaZ. */
alphaCmin = alphaZ;
cIndex[2] = cIndex[2] + kU;
alphaZ = alphaZ + alphaUz;
alphaXmin = -2;
alphaXmax = 2;
}
if ((cIndex[0] >= 0) && (cIndex[0] < static_cast<IndexValueType>(sizeCT[0])) && (cIndex[1] >= 0) &&
(cIndex[1] < static_cast<IndexValueType>(sizeCT[1])) && (cIndex[2] >= 0) &&
(cIndex[2] < static_cast<IndexValueType>(sizeCT[2])))
if (fabs(rayVector[1]) > EPSILON/* != 0*/)
{
/* If it is a valid index, get the voxel intensity. */
value = static_cast<float>(inputPtr->GetPixel(cIndex));
if (value > m_Threshold) /* Ignore voxels whose intensities are below the threshold. */
{
d12 += (alphaCmin - alphaCminPrev) * (value - m_Threshold);
}
alphaY1 = ( O[1] - SourceWorld[1]) / rayVector[1];
alphaYN = ( O[1] + sizeCT[1] * ctPixelSpacing[1] - SourceWorld[1]) / rayVector[1];
// alphaY1 = (0.0 - SourceWorld[1]) / rayVector[1];
// alphaYN = (sizeCT[1] * ctPixelSpacing[1] - SourceWorld[1]) / rayVector[1];
alphaYmin = std::min(alphaY1, alphaYN);
alphaYmax = std::max(alphaY1, alphaYN);
}
else
{
alphaYmin = -2;
alphaYmax = 2;
}
}
if (d12 < minOutputValue)
{
pixval = minOutputValue;
}
else if (d12 > maxOutputValue)
{
pixval = maxOutputValue;
}
else
{
pixval = static_cast<OutputType>(d12);
}
return (pixval);
if (fabs(rayVector[2]) > EPSILON/* != 0*/)
{
// alphaZ1 = (0.0 - SourceWorld[2]) / rayVector[2];
// alphaZN = (sizeCT[2] * ctPixelSpacing[2] - SourceWorld[2]) / rayVector[2];
alphaZ1 = ( O[2] - SourceWorld[2]) / rayVector[2];
alphaZN = ( O[2] + sizeCT[2] * ctPixelSpacing[2] - SourceWorld[2]) / rayVector[2];
alphaZmin = std::min(alphaZ1, alphaZN);
alphaZmax = std::max(alphaZ1, alphaZN);
}
else
{
alphaZmin = -2;
alphaZmax = 2;
}
/* Get the very first and the last alpha values when the ray
intersects with the CT volume. */
alphaMin = std::max(std::max(alphaXmin, alphaYmin), alphaZmin);
alphaMax = std::min(std::min(alphaXmax, alphaYmax), alphaZmax);
/* Calculate the parametric values of the first intersection point
of the ray with the X, Y, and Z-planes after the ray entered the
CT volume. */
firstIntersection[0] = SourceWorld[0] + alphaMin * rayVector[0] - O[0] ; //ctPixelSpacing[0]/2.;
firstIntersection[1] = SourceWorld[1] + alphaMin * rayVector[1] - O[1] ; //ctPixelSpacing[1]/2.;
firstIntersection[2] = SourceWorld[2] + alphaMin * rayVector[2] - O[2] ; //ctPixelSpacing[2]/2.;
/* Transform world coordinate to the continuous index of the CT volume*/
firstIntersectionIndex[0] = firstIntersection[0] / ctPixelSpacing[0];
firstIntersectionIndex[1] = firstIntersection[1] / ctPixelSpacing[1];
firstIntersectionIndex[2] = firstIntersection[2] / ctPixelSpacing[2];
// now "correct" the indices (there are special cases):
if (alphaMin == alphaYmin) // j is special
{
if (alphaYmin == alphaY1)
firstIntersectionIndex [1] = 0;
else
// alpha_y_Ny
firstIntersectionIndex [1] = sizeCT[1] - 1;
}
else if (alphaMin == alphaXmin) // i is special
{
if (alphaXmin == alphaX1)
firstIntersectionIndex[0] = 0;
else
// alpha_x_Nx
firstIntersectionIndex[0] = sizeCT[0] - 1;
}
else if (alphaMin == alphaZmin) // k is special
{
if (alphaZmin == alphaZ1)
firstIntersectionIndex [2] = 0;
else
// alpha_z_Nz
firstIntersectionIndex [2] = sizeCT[2] - 1;
}
firstIntersectionIndexUp[0] = (int)ceil(firstIntersectionIndex[0]);
firstIntersectionIndexUp[1] = (int)ceil(firstIntersectionIndex[1]);
firstIntersectionIndexUp[2] = (int)ceil(firstIntersectionIndex[2]);
firstIntersectionIndexDown[0] = (int)floor(firstIntersectionIndex[0]);
firstIntersectionIndexDown[1] = (int)floor(firstIntersectionIndex[1]);
firstIntersectionIndexDown[2] = (int)floor(firstIntersectionIndex[2]);
if (fabs(rayVector[0]) <= EPSILON/* == 0*/)
{
alphaX = 2;
}
else
{
alphaIntersectionUp[0] =
(O[0] + firstIntersectionIndexUp[0] * ctPixelSpacing[0] - SourceWorld[0]) / rayVector[0];
alphaIntersectionDown[0] =
(O[0] + firstIntersectionIndexDown[0] * ctPixelSpacing[0] - SourceWorld[0]) / rayVector[0];
alphaX = std::max(alphaIntersectionUp[0], alphaIntersectionDown[0]);
}
if (fabs(rayVector[1] ) <= EPSILON/* == 0*/)
{
alphaY = 2;
}
else
{
alphaIntersectionUp[1] =
(O[1] + firstIntersectionIndexUp[1] * ctPixelSpacing[1] - SourceWorld[1]) / rayVector[1];
alphaIntersectionDown[1] =
(O[1] + firstIntersectionIndexDown[1] * ctPixelSpacing[1] - SourceWorld[1]) / rayVector[1];
alphaY = std::max(alphaIntersectionUp[1], alphaIntersectionDown[1]);
}
if (fabs(rayVector[2]) <= EPSILON/* == 0*/)
{
alphaZ = 2;
}
else
{
alphaIntersectionUp[2] =
(O[2] + firstIntersectionIndexUp[2] * ctPixelSpacing[2] - SourceWorld[2]) / rayVector[2];
alphaIntersectionDown[2] =
(O[2] + firstIntersectionIndexDown[2] * ctPixelSpacing[2] - SourceWorld[2]) / rayVector[2];
alphaZ = std::max(alphaIntersectionUp[2], alphaIntersectionDown[2]);
}
/* Calculate alpha incremental values when the ray intercepts with x, y, and z-planes */
if (fabs(rayVector[0]) > EPSILON /*!= 0*/)
{
alphaUx = ctPixelSpacing[0] / fabs(rayVector[0]);
}
else
{
alphaUx = 999;
}
if (fabs(rayVector[1]) > EPSILON /*!= 0*/)
{
alphaUy = ctPixelSpacing[1] / fabs(rayVector[1]);
}
else
{
alphaUy = 999;
}
if (fabs(rayVector[2]) > EPSILON /*!= 0*/)
{
alphaUz = ctPixelSpacing[2] / fabs(rayVector[2]);
}
else
{
alphaUz = 999;
}
/* Calculate voxel index incremental values along the ray path. */
iU = (SourceWorld[0] < drrPixelWorld[0]) ? 1 : -1;
jU = (SourceWorld[1] < drrPixelWorld[1]) ? 1 : -1;
kU = (SourceWorld[2] < drrPixelWorld[2]) ? 1 : -1;
d12 = 0.0; /* Initialize the sum of the voxel intensities along the ray path to zero. */
/* Initialize the current ray position. */
alphaCmin = std::min(std::min(alphaX, alphaY), alphaZ);
/* Initialize the current voxel index. */
cIndex[0] = firstIntersectionIndexDown[0];
cIndex[1] = firstIntersectionIndexDown[1];
cIndex[2] = firstIntersectionIndexDown[2];
while (alphaCmin < alphaMax) /* Check if the ray is still in the CT volume */
{
/* Store the current ray position */
alphaCminPrev = alphaCmin;
if ((alphaX <= alphaY) && (alphaX <= alphaZ))
{
/* Current ray front intercepts with x-plane. Update alphaX. */
alphaCmin = alphaX;
cIndex[0] = cIndex[0] + iU;
alphaX = alphaX + alphaUx;
}
else if ((alphaY <= alphaX) && (alphaY <= alphaZ))
{
/* Current ray front intercepts with y-plane. Update alphaY. */
alphaCmin = alphaY;
cIndex[1] = cIndex[1] + jU;
alphaY = alphaY + alphaUy;
}
else
{
/* Current ray front intercepts with z-plane. Update alphaZ. */
alphaCmin = alphaZ;
cIndex[2] = cIndex[2] + kU;
alphaZ = alphaZ + alphaUz;
}
if ((cIndex[0] >= 0) && (cIndex[0] < static_cast<IndexValueType>(sizeCT[0])) && (cIndex[1] >= 0) &&
(cIndex[1] < static_cast<IndexValueType>(sizeCT[1])) && (cIndex[2] >= 0) &&
(cIndex[2] < static_cast<IndexValueType>(sizeCT[2])))
{
/* If it is a valid index, get the voxel intensity. */
value = static_cast<float>(inputPtr->GetPixel(cIndex));
if (value > m_Threshold) /* Ignore voxels whose intensities are below the threshold. */
{
d12 += (alphaCmin - alphaCminPrev) * (value - m_Threshold);
}
}
}
if (d12 < minOutputValue)
{
pixval = minOutputValue;
}
else if (d12 > maxOutputValue)
{
pixval = maxOutputValue;
}
else
{
pixval = static_cast<OutputType>(d12);
}
return (pixval);
}
template <typename TInputImage, typename TCoordRep>
typename gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::OutputType
gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::EvaluateAtContinuousIndex(
const ContinuousIndexType & index) const
const ContinuousIndexType & index) const
{
OutputPointType point;
this->m_Image->TransformContinuousIndexToPhysicalPoint(index, point);
OutputPointType point;
this->m_Image->TransformContinuousIndexToPhysicalPoint(index, point);
return this->Evaluate(point);
return this->Evaluate(point);
}
@ -519,43 +502,52 @@ template <typename TInputImage, typename TCoordRep>
void
gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::ComputeInverseTransform() const
{
m_ComposedTransform->SetIdentity();
m_ComposedTransform->Compose(m_Transform, 0);
m_ComposedTransform->SetIdentity();
m_ComposedTransform->Compose(m_Transform, 0);
typename TransformType::InputPointType isocenter;
isocenter = m_Transform->GetCenter();
typename TransformType::InputPointType isocenter;
isocenter = m_Transform->GetCenter();
// if(this->GetInputImage() != NULL)
//{
// isocenter[0] -= this->GetInputImage()->GetSpacing()[0];
// isocenter[1] -= this->GetInputImage()->GetSpacing()[1];
// isocenter[2] -= this->GetInputImage()->GetSpacing()[2];
//}
// An Euler 3D transform is used to rotate the volume to simulate the roation of the linac gantry.
// The rotation is about z-axis. After the transform, a AP projection geometry (projecting
// towards positive y direction) is established.
m_GantryRotTransform->SetRotation(0.0, 0.0, -m_ProjectionAngle);
m_GantryRotTransform->SetCenter(isocenter);
m_ComposedTransform->Compose(m_GantryRotTransform, 0);
// //This is problably not necessary since we correct the ray.... to be checked..
// InputImageConstPointer inputPtr = this->GetInputImage();
// An Euler 3D transfrom is used to shift the source to the origin.
typename TransformType::OutputVectorType focalpointtranslation;
focalpointtranslation[0] = -isocenter[0];
focalpointtranslation[1] = m_FocalPointToIsocenterDistance - isocenter[1];
focalpointtranslation[2] = -isocenter[2];
m_CamShiftTransform->SetTranslation(focalpointtranslation);
m_ComposedTransform->Compose(m_CamShiftTransform, 0);
// std::cout<<inputPtr<<std::endl;
// if(inputPtr != NULL)
// {
// isocenter[0] -= inputPtr->GetSpacing()[0]/2.;
// isocenter[1] -= inputPtr->GetSpacing()[1]/2.;
// isocenter[2] -= inputPtr->GetSpacing()[2]/2.;
// std::cout<<inputPtr->GetSpacing()<<std::endl;
// } else {
// isocenter[0] += 0.5;
// isocenter[1] += 0.5;
// isocenter[2] += 1;
// }
// An Euler 3D transform is used to rotate the volume to simulate the roation of the linac gantry.
// The rotation is about z-axis. After the transform, a AP projection geometry (projecting
// towards positive y direction) is established.
m_GantryRotTransform->SetRotation(0.0, 0.0, -m_ProjectionAngle);
m_GantryRotTransform->SetCenter(isocenter);
m_ComposedTransform->Compose(m_GantryRotTransform, 0);
// A Euler 3D transform is used to establish the standard negative z-axis projection geometry. (By
// default, the camera is situated at the origin, points down the negative z-axis, and has an up-
// vector of (0, 1, 0).)
// An Euler 3D transfrom is used to shift the source to the origin.
typename TransformType::OutputVectorType focalpointtranslation;
focalpointtranslation[0] = -isocenter[0];
focalpointtranslation[1] = m_FocalPointToIsocenterDistance - isocenter[1];
focalpointtranslation[2] = -isocenter[2];
m_CamShiftTransform->SetTranslation(focalpointtranslation);
m_ComposedTransform->Compose(m_CamShiftTransform, 0);
m_ComposedTransform->Compose(m_CamRotTransform, 0);
// A Euler 3D transform is used to establish the standard negative z-axis projection geometry. (By
// default, the camera is situated at the origin, points down the negative z-axis, and has an up-
// vector of (0, 1, 0).)
// The overall inverse transform is computed. The inverse transform will be used by the interpolation
// procedure.
m_ComposedTransform->GetInverse(m_InverseTransform);
this->Modified();
m_ComposedTransform->Compose(m_CamRotTransform, 0);
// The overall inverse transform is computed. The inverse transform will be used by the interpolation
// procedure.
m_ComposedTransform->GetInverse(m_InverseTransform);
this->Modified();
}
@ -563,8 +555,8 @@ template <typename TInputImage, typename TCoordRep>
void
gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Initialize()
{
this->ComputeInverseTransform();
m_SourceWorld = m_InverseTransform->TransformPoint(m_SourcePoint);
this->ComputeInverseTransform();
m_SourceWorld = m_InverseTransform->TransformPoint(m_SourcePoint);
}
} // namespace itk