658 lines
25 KiB
C++
658 lines
25 KiB
C++
|
|
/*
|
|
|
|
Calculate Ditigally Reconstructed Topogram from a CT dataset using a modified
|
|
version of incremental ray-tracing algorithm
|
|
|
|
gfattori 08.11.2021
|
|
|
|
*/
|
|
|
|
|
|
|
|
/*=========================================================================
|
|
*
|
|
* Copyright NumFOCUS
|
|
*
|
|
* Licensed under the Apache License, Version 2.0 (the "License");
|
|
* you may not use this file except in compliance with the License.
|
|
* You may obtain a copy of the License at
|
|
*
|
|
* http://www.apache.org/licenses/LICENSE-2.0.txt
|
|
*
|
|
* Unless required by applicable law or agreed to in writing, software
|
|
* distributed under the License is distributed on an "AS IS" BASIS,
|
|
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
|
* See the License for the specific language governing permissions and
|
|
* limitations under the License.
|
|
*
|
|
*=========================================================================*/
|
|
/*=========================================================================
|
|
The algorithm was initially proposed by Robert Siddon and improved by
|
|
Filip Jacobs etc.
|
|
|
|
-------------------------------------------------------------------------
|
|
References:
|
|
|
|
R. L. Siddon, "Fast calculation of the exact radiological path for a
|
|
threedimensional CT array," Medical Physics 12, 252-55 (1985).
|
|
|
|
F. Jacobs, E. Sundermann, B. De Sutter, M. Christiaens, and I. Lemahieu,
|
|
"A fast algorithm to calculate the exact radiological path through a pixel
|
|
or voxel space," Journal of Computing and Information Technology ?
|
|
CIT 6, 89-94 (1998).
|
|
|
|
=========================================================================*/
|
|
|
|
#ifndef itkgSiddonJacobsRayCastInterpolateImageFunction_hxx
|
|
#define itkgSiddonJacobsRayCastInterpolateImageFunction_hxx
|
|
|
|
#include "itkgSiddonJacobsRayCastInterpolateImageFunction.h"
|
|
#include "itkContinuousIndex.h"
|
|
|
|
#include "itkMath.h"
|
|
#include <cstdlib>
|
|
#include <unistd.h>
|
|
|
|
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_SourcePoint[0] = 0.;
|
|
m_SourcePoint[1] = 0.;
|
|
m_SourcePoint[2] = 0.;
|
|
|
|
m_InverseTransform = TransformType::New();
|
|
m_InverseTransform->SetComputeZYX(true);
|
|
|
|
m_ComposedTransform = TransformType::New();
|
|
m_ComposedTransform->SetComputeZYX(true);
|
|
|
|
m_GantryRotTransform = TransformType::New();
|
|
m_GantryRotTransform->SetComputeZYX(true);
|
|
m_GantryRotTransform->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_Threshold = 0;
|
|
}
|
|
|
|
|
|
template <typename TInputImage, typename TCoordRep>
|
|
void
|
|
gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::PrintSelf(std::ostream & os, Indent indent) const
|
|
{
|
|
this->Superclass::PrintSelf(os, indent);
|
|
|
|
os << indent << "Threshold: " << m_Threshold << std::endl;
|
|
os << indent << "Transform: " << m_Transform.GetPointer() << std::endl;
|
|
}
|
|
|
|
|
|
template <typename TInputImage, typename TCoordRep>
|
|
typename gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::OutputType
|
|
gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(const PointType & point) const
|
|
{
|
|
float rayVector[3];
|
|
IndexType cIndex;
|
|
|
|
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,valuetril;
|
|
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();
|
|
|
|
// 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);
|
|
}
|
|
|
|
PointType PointReq = point;
|
|
//std::cout<<"PointReq: "<<point[0] <<" "<<point[1] <<" "<<point[2] <<" "<<std::endl;
|
|
PointReq[0] += m_PanelOffset;
|
|
|
|
drrPixelWorld = m_InverseTransform->TransformPoint(PointReq);
|
|
//std::cout<<"drrPixelWorld: "<<drrPixelWorld[0] <<" "<<drrPixelWorld[1] <<" "<<drrPixelWorld[2] <<" "<<std::endl;
|
|
|
|
|
|
// 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.;
|
|
|
|
PointType SourceWorld = m_InverseTransform->TransformPoint(SlidingSourcePoint);
|
|
//std::cout<<"SourceWorld: "<<SourceWorld[0] <<" "<<SourceWorld[1] <<" "<<SourceWorld[2] <<" "<<std::endl;
|
|
|
|
PointType O(3);
|
|
O[0] = -ctPixelSpacing[0]/2.;
|
|
O[1] = -ctPixelSpacing[1]/2.;
|
|
O[2] = -ctPixelSpacing[2]/2.;
|
|
|
|
// 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];
|
|
|
|
|
|
// 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];
|
|
|
|
|
|
/* 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 (fabs(rayVector[0]) > EPSILON/* != 0*/)
|
|
{
|
|
// 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
|
|
{
|
|
alphaXmin = -2;
|
|
alphaXmax = 2;
|
|
}
|
|
|
|
if (fabs(rayVector[1]) > EPSILON/* != 0*/)
|
|
{
|
|
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 (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])))
|
|
{
|
|
// Calculate entry and exit points using alphaCmin and alphaCminPrev
|
|
value = static_cast<float>(inputPtr->GetPixel(cIndex));
|
|
// Accumulate the interpolated intensity along the ray path
|
|
if (value > m_Threshold) /* Ignore voxels whose intensities are below the threshold. */
|
|
{
|
|
|
|
// Move along the ray by alphaCminPrev to find the entry point of this voxel
|
|
PointType entryPoint;
|
|
entryPoint[0] = SourceWorld[0] + alphaCminPrev * rayVector[0] ;
|
|
entryPoint[1] = SourceWorld[1] + alphaCminPrev * rayVector[1] ;
|
|
entryPoint[2] = SourceWorld[2] + alphaCminPrev * rayVector[2] ;
|
|
|
|
// Move along the ray by alphaCmin to find the exit point of this voxel
|
|
PointType exitPoint;
|
|
exitPoint[0] = SourceWorld[0] + alphaCmin * rayVector[0] ;
|
|
exitPoint[1] = SourceWorld[1] + alphaCmin * rayVector[1] ;
|
|
exitPoint[2] = SourceWorld[2] + alphaCmin * rayVector[2] ;
|
|
|
|
// Get the mid-point of the voxel / ray interception
|
|
PointType midpoint;
|
|
midpoint[0]= (entryPoint[0] + exitPoint[0]) * 0.5;
|
|
midpoint[1]= (entryPoint[1] + exitPoint[1]) * 0.5;
|
|
midpoint[2]= (entryPoint[2] + exitPoint[2]) * 0.5;
|
|
|
|
// Ray is computed in 'Siddon' geometry with zero origin,
|
|
// whereas the continuous index will be computed from the input
|
|
// image. The origin and shifts to the voxel edge are to be account for
|
|
midpoint[0] += ctOrigin[0] ;//+ O[0];
|
|
midpoint[1] += ctOrigin[1] ;//+ O[1];
|
|
midpoint[2] += ctOrigin[2] ;//+ O[2];
|
|
|
|
// Convert mid-point phyisical point into continuous index
|
|
// We need to use this position to find the neighbouring voxels
|
|
// for trilinear interpolation - not the cIndex!
|
|
itk::ContinuousIndex <double,3> continuousIndex;
|
|
inputPtr->TransformPhysicalPointToContinuousIndex(midpoint,continuousIndex);
|
|
|
|
// Get the baseIndex by flooring the continuous index
|
|
IndexType baseIndex;
|
|
baseIndex[0]=static_cast<IndexValueType>(std::floor(continuousIndex[0]));
|
|
baseIndex[1]=static_cast<IndexValueType>(std::floor(continuousIndex[1]));
|
|
baseIndex[2]=static_cast<IndexValueType>(std::floor(continuousIndex[2]));
|
|
|
|
// Calculate fractional parts for interpolation at the midpoint
|
|
double x_frac = continuousIndex[0] - baseIndex[0];
|
|
double y_frac = continuousIndex[1] - baseIndex[1];
|
|
double z_frac = continuousIndex[2] - baseIndex[2];
|
|
|
|
// Perform boundary checks for trilinear interpolation
|
|
bool within_bounds = (baseIndex[0] >= 0 && baseIndex[0] < sizeCT[0] - 1) &&
|
|
(baseIndex[1] >= 0 && baseIndex[1] < sizeCT[1] - 1) &&
|
|
(baseIndex[2] >= 0 && baseIndex[2] < sizeCT[2] - 1);
|
|
if (within_bounds)
|
|
{
|
|
// Fetch intensities from neighboring voxels for trilinear interpolation
|
|
float c000 = static_cast<float>(inputPtr->GetPixel(baseIndex));
|
|
baseIndex[0] += 1;
|
|
float c100 = static_cast<float>(inputPtr->GetPixel(baseIndex));
|
|
baseIndex[1] += 1;
|
|
float c110 = static_cast<float>(inputPtr->GetPixel(baseIndex));
|
|
baseIndex[0] -= 1;
|
|
float c010 = static_cast<float>(inputPtr->GetPixel(baseIndex));
|
|
baseIndex[2] += 1;
|
|
float c011 = static_cast<float>(inputPtr->GetPixel(baseIndex));
|
|
baseIndex[0] += 1;
|
|
float c111 = static_cast<float>(inputPtr->GetPixel(baseIndex));
|
|
baseIndex[1] -= 1;
|
|
float c101 = static_cast<float>(inputPtr->GetPixel(baseIndex));
|
|
baseIndex[0] -= 1;
|
|
float c001 = static_cast<float>(inputPtr->GetPixel(baseIndex));
|
|
|
|
// NOTE: do not use baseIndex anymore after this. It's messed up.
|
|
|
|
// Perform trilinear interpolation at the midpoint
|
|
float c00 = c000 * (1 - x_frac) + c100 * x_frac;
|
|
float c01 = c001 * (1 - x_frac) + c101 * x_frac;
|
|
float c10 = c010 * (1 - x_frac) + c110 * x_frac;
|
|
float c11 = c011 * (1 - x_frac) + c111 * x_frac;
|
|
|
|
float c0 = c00 * (1 - y_frac) + c10 * y_frac;
|
|
float c1 = c01 * (1 - y_frac) + c11 * y_frac;
|
|
|
|
valuetril = c0 * (1 - z_frac) + c1 * z_frac;
|
|
|
|
d12 += (alphaCmin - alphaCminPrev) * (valuetril - m_Threshold);
|
|
|
|
} else {
|
|
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
|
|
{
|
|
OutputPointType point;
|
|
this->m_Image->TransformContinuousIndexToPhysicalPoint(index, point);
|
|
|
|
return this->Evaluate(point);
|
|
}
|
|
|
|
|
|
template <typename TInputImage, typename TCoordRep>
|
|
void
|
|
gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::ComputeInverseTransform() const
|
|
{
|
|
m_ComposedTransform->SetIdentity();
|
|
m_ComposedTransform->Compose(m_Transform, 0);
|
|
|
|
typename TransformType::InputPointType isocenter;
|
|
isocenter = m_Transform->GetCenter();
|
|
|
|
// //This is problably not necessary since we correct the ray.... to be checked..
|
|
// InputImageConstPointer inputPtr = this->GetInputImage();
|
|
|
|
// 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);
|
|
|
|
// 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);
|
|
|
|
// 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).)
|
|
|
|
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();
|
|
}
|
|
|
|
|
|
template <typename TInputImage, typename TCoordRep>
|
|
void
|
|
gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Initialize()
|
|
{
|
|
this->ComputeInverseTransform();
|
|
m_SourceWorld = m_InverseTransform->TransformPoint(m_SourcePoint);
|
|
}
|
|
|
|
} // namespace itk
|
|
|
|
#endif
|