diff --git a/itkReg23DRT/itkgSiddonJacobsRayCastInterpolateImageFunction.h b/itkReg23DRT/itkgSiddonJacobsRayCastInterpolateImageFunction.h index 2df4f71..6093f4e 100644 --- a/itkReg23DRT/itkgSiddonJacobsRayCastInterpolateImageFunction.h +++ b/itkReg23DRT/itkgSiddonJacobsRayCastInterpolateImageFunction.h @@ -133,6 +133,7 @@ public: /** Internal floating point comparison accuracy **/ static const double EPSILON; + static const double TOL; /** \brief * Interpolate the image at a point position. diff --git a/itkReg23DRT/itkgSiddonJacobsRayCastInterpolateImageFunction.hxx b/itkReg23DRT/itkgSiddonJacobsRayCastInterpolateImageFunction.hxx index 4bf5335..896139e 100644 --- a/itkReg23DRT/itkgSiddonJacobsRayCastInterpolateImageFunction.hxx +++ b/itkReg23DRT/itkgSiddonJacobsRayCastInterpolateImageFunction.hxx @@ -56,6 +56,9 @@ CIT 6, 89-94 (1998). namespace itk { +template +const double gSiddonJacobsRayCastInterpolateImageFunction::TOL = + 1e-9; template const double gSiddonJacobsRayCastInterpolateImageFunction::EPSILON = @@ -113,46 +116,38 @@ template typename gSiddonJacobsRayCastInterpolateImageFunction::OutputType gSiddonJacobsRayCastInterpolateImageFunction::Evaluate(const PointType & point) const { - float rayVector[3]; + // --- Geometry / Siddon working variables use double --- + double 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; - float firstIntersectionIndex[3]; - int firstIntersectionIndexUp[3], firstIntersectionIndexDown[3]; - int iU, jU, kU; - + double firstIntersection[3]; + double alphaX1, alphaXN, alphaXmin, alphaXmax; + double alphaY1, alphaYN, alphaYmin, alphaYmax; + double alphaZ1, alphaZN, alphaZmin, alphaZmax; + double alphaMin, alphaMax; + double alphaX, alphaY, alphaZ; + double alphaUx, alphaUy, alphaUz; + double alphaIntersectionUp[3], alphaIntersectionDown[3]; + double d12, value; + double 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::NonpositiveMin(); const OutputType maxOutputValue = itk::NumericTraits::max(); - // If the volume was shifted, recalculate the overall inverse transform - unsigned long int interpMTime = this->GetMTime(); + 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; @@ -160,7 +155,7 @@ gSiddonJacobsRayCastInterpolateImageFunction::Evaluate(c drrPixelWorld = m_InverseTransform->TransformPoint(PointReq); - // Get ths input pointers + // Get the input pointers InputImageConstPointer inputPtr = this->GetInputImage(); typename InputImageType::SizeType sizeCT; @@ -169,13 +164,12 @@ gSiddonJacobsRayCastInterpolateImageFunction::Evaluate(c typename InputImageType::PointType ctOrigin; ctPixelSpacing = inputPtr->GetSpacing(); - ctOrigin = inputPtr->GetOrigin(); - regionCT = inputPtr->GetLargestPossibleRegion(); - sizeCT = regionCT.GetSize(); - + ctOrigin = inputPtr->GetOrigin(); + regionCT = inputPtr->GetLargestPossibleRegion(); + sizeCT = regionCT.GetSize(); PointType SlidingSourcePoint = m_SourcePoint; - SlidingSourcePoint[0] = 0.; + 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]; @@ -184,22 +178,22 @@ gSiddonJacobsRayCastInterpolateImageFunction::Evaluate(c PointType SourceWorld = m_InverseTransform->TransformPoint(SlidingSourcePoint); PointType O(3); - O[0] = -ctPixelSpacing[0]/2.; - O[1] = -ctPixelSpacing[1]/2.; - O[2] = -ctPixelSpacing[2]/2.; + 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 - ){ + double xSizeCT = O[0] + (sizeCT[0]) * ctPixelSpacing[0]; + double ySizeCT = O[1] + (sizeCT[1]) * ctPixelSpacing[1]; + double zSizeCT = O[2] + (sizeCT[2]) * ctPixelSpacing[2]; + double xDrrPix = drrPixelWorld[0]; + double yDrrPix = drrPixelWorld[1]; + double zDrrPix = drrPixelWorld[2]; + + if (xDrrPix <= O[0] || xDrrPix >= xSizeCT || + yDrrPix <= O[1] || yDrrPix >= ySizeCT || + zDrrPix <= O[2] || zDrrPix >= zSizeCT) + { pixval = static_cast(0.0); return pixval; } @@ -217,260 +211,351 @@ gSiddonJacobsRayCastInterpolateImageFunction::Evaluate(c 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 P12D = std::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]; + 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.0 * p12v[0]; + pDest[1] = p1[1] + m_FocalPointToIsocenterDistance * 2.0 * p12v[1]; + pDest[2] = p1[2] + m_FocalPointToIsocenterDistance * 2.0 * 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. + // Compute ray from source to "detector" point pDest rayVector[0] = pDest[0] - SourceWorld[0]; rayVector[1] = pDest[1] - SourceWorld[1]; rayVector[2] = pDest[2] - SourceWorld[2]; + const double rayLength = std::sqrt(rayVector[0]*rayVector[0] + + rayVector[1]*rayVector[1] + + rayVector[2]*rayVector[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*/) + // just in case (should never really be zero) + if (rayLength <= EPSILON) { - // 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]; + pixval = static_cast(0.0); + return pixval; + } + const double INF = std::numeric_limits::infinity(); + + // Convenience bounds for the CT volume in world coords + double xMin = O[0]; + double xMax = O[0] + sizeCT[0] * ctPixelSpacing[0]; + double yMin = O[1]; + double yMax = O[1] + sizeCT[1] * ctPixelSpacing[1]; + double zMin = O[2]; + double zMax = O[2] + sizeCT[2] * ctPixelSpacing[2]; + + /* Calculate the parametric values of the first and last + intersection points of the ray with the X, Y, and Z-planes + that define the CT volume. */ + + // ---- X slab ---- + if (std::fabs(rayVector[0]) > EPSILON) + { + alphaX1 = (xMin - SourceWorld[0]) / rayVector[0]; + alphaXN = (xMax - SourceWorld[0]) / rayVector[0]; alphaXmin = std::min(alphaX1, alphaXN); alphaXmax = std::max(alphaX1, alphaXN); } else { - alphaXmin = -2; - alphaXmax = 2; + // Ray is parallel to X planes: must lie within [xMin, xMax] or no intersection + if (SourceWorld[0] < xMin || SourceWorld[0] > xMax) + { + pixval = static_cast(0.0); + return pixval; + } + alphaXmin = -INF; + alphaXmax = INF; } - if (fabs(rayVector[1]) > EPSILON/* != 0*/) + // ---- Y slab ---- + if (std::fabs(rayVector[1]) > EPSILON) { - 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]; + alphaY1 = (yMin - SourceWorld[1]) / rayVector[1]; + alphaYN = (yMax - SourceWorld[1]) / rayVector[1]; alphaYmin = std::min(alphaY1, alphaYN); alphaYmax = std::max(alphaY1, alphaYN); } else { - alphaYmin = -2; - alphaYmax = 2; + if (SourceWorld[1] < yMin || SourceWorld[1] > yMax) + { + pixval = static_cast(0.0); + return pixval; + } + alphaYmin = -INF; + alphaYmax = INF; } - if (fabs(rayVector[2]) > EPSILON/* != 0*/) + // ---- Z slab ---- + if (std::fabs(rayVector[2]) > EPSILON) { - // 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]; + alphaZ1 = (zMin - SourceWorld[2]) / rayVector[2]; + alphaZN = (zMax - SourceWorld[2]) / rayVector[2]; alphaZmin = std::min(alphaZ1, alphaZN); alphaZmax = std::max(alphaZ1, alphaZN); } else { - alphaZmin = -2; - alphaZmax = 2; + if (SourceWorld[2] < zMin || SourceWorld[2] > zMax) + { + pixval = static_cast(0.0); + return pixval; + } + alphaZmin = -INF; + alphaZmax = INF; } /* Get the very first and the last alpha values when the ray - intersects with the CT volume. */ + intersects with the CT volume. */ alphaMin = std::max(std::max(alphaXmin, alphaYmin), alphaZmin); alphaMax = std::min(std::min(alphaXmax, alphaYmax), alphaZmax); + // No intersection at all + if (alphaMax <= alphaMin) + { + pixval = static_cast(0.0); + return pixval; + } + + // Restrict to the [0,1] segment (from SourceWorld to pDest) + alphaMax = std::min(alphaMax, 1.0); + alphaMin = std::max(alphaMin, 0.0); + + if (alphaMax <= alphaMin) + { + pixval = static_cast(0.0); + return pixval; + } + /* 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. */ + after the ray entered the CT volume. */ + firstIntersection[0] = SourceWorld[0] + alphaMin * rayVector[0] - O[0]; + firstIntersection[1] = SourceWorld[1] + alphaMin * rayVector[1] - O[1]; + firstIntersection[2] = SourceWorld[2] + alphaMin * rayVector[2] - O[2]; - 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*/ + /* Transform world coordinates to continuous indices 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 + // Helper for equality with tolerance (for on-grid cases) + auto nearlyEqual = [](double a, double b) { - if (alphaYmin == alphaY1) - firstIntersectionIndex [1] = 0; - else - // alpha_y_Ny - firstIntersectionIndex [1] = sizeCT[1] - 1; - } - else if (alphaMin == alphaXmin) // i is special + return std::fabs(a - b) < TOL; + }; + + // "Correct" the indices for special entry cases + if (nearlyEqual(alphaMin, alphaXmin)) // i is special { - if (alphaXmin == alphaX1) - firstIntersectionIndex[0] = 0; + if (nearlyEqual(alphaXmin, alphaX1)) + firstIntersectionIndex[0] = 0.0; // continuous index 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; + firstIntersectionIndex[0] = sizeCT[0] - 1.; // continuous index } - 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*/) + if (nearlyEqual(alphaMin, alphaYmin)) // j is special { - alphaX = 2; + if (nearlyEqual(alphaYmin, alphaY1)) + firstIntersectionIndex[1] = 0.0; + else + firstIntersectionIndex[1] = sizeCT[1] - 1.; + } + + if (nearlyEqual(alphaMin, alphaZmin)) // k is special + { + if (nearlyEqual(alphaZmin, alphaZ1)) + firstIntersectionIndex[2] = 0.0; + else + firstIntersectionIndex[2] = sizeCT[2] - 1.; + } + + firstIntersectionIndexUp[0] = static_cast(std::ceil(firstIntersectionIndex[0])); + firstIntersectionIndexUp[1] = static_cast(std::ceil(firstIntersectionIndex[1])); + firstIntersectionIndexUp[2] = static_cast(std::ceil(firstIntersectionIndex[2])); + firstIntersectionIndexDown[0] = static_cast(std::floor(firstIntersectionIndex[0])); + firstIntersectionIndexDown[1] = static_cast(std::floor(firstIntersectionIndex[1])); + firstIntersectionIndexDown[2] = static_cast(std::floor(firstIntersectionIndex[2])); + + // ---- Compute first plane hits and increments ---- + + // X + if (std::fabs(rayVector[0]) <= EPSILON) + { + alphaX = INF; + alphaUx = INF; } else { alphaIntersectionUp[0] = - (O[0] + firstIntersectionIndexUp[0] * ctPixelSpacing[0] - SourceWorld[0]) / rayVector[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]); + (O[0] + firstIntersectionIndexDown[0] * ctPixelSpacing[0] - SourceWorld[0]) / rayVector[0]; + alphaX = std::max(alphaIntersectionUp[0], alphaIntersectionDown[0]); + alphaUx = ctPixelSpacing[0] / std::fabs(rayVector[0]); } - if (fabs(rayVector[1] ) <= EPSILON/* == 0*/) + // Y + if (std::fabs(rayVector[1]) <= EPSILON) { - alphaY = 2; + alphaY = INF; + alphaUy = INF; } else { alphaIntersectionUp[1] = - (O[1] + firstIntersectionIndexUp[1] * ctPixelSpacing[1] - SourceWorld[1]) / rayVector[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]); + (O[1] + firstIntersectionIndexDown[1] * ctPixelSpacing[1] - SourceWorld[1]) / rayVector[1]; + alphaY = std::max(alphaIntersectionUp[1], alphaIntersectionDown[1]); + alphaUy = ctPixelSpacing[1] / std::fabs(rayVector[1]); } - if (fabs(rayVector[2]) <= EPSILON/* == 0*/) + // Z + if (std::fabs(rayVector[2]) <= EPSILON) { - alphaZ = 2; + alphaZ = INF; + alphaUz = INF; } else { alphaIntersectionUp[2] = - (O[2] + firstIntersectionIndexUp[2] * ctPixelSpacing[2] - SourceWorld[2]) / rayVector[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]); + (O[2] + firstIntersectionIndexDown[2] * ctPixelSpacing[2] - SourceWorld[2]) / rayVector[2]; + alphaZ = std::max(alphaIntersectionUp[2], alphaIntersectionDown[2]); + alphaUz = ctPixelSpacing[2] / std::fabs(rayVector[2]); } - /* Calculate alpha incremental values when the ray intercepts with x, y, and z-planes */ - if (fabs(rayVector[0]) > EPSILON /*!= 0*/) + // Make sure the first plane alphas are not before entry + + // Only advance if this axis is not parallel (alphaUx finite) + if (!std::isinf(alphaUx)) { - alphaUx = ctPixelSpacing[0] / fabs(rayVector[0]); + while (alphaX <= alphaMin + TOL) + alphaX += alphaUx; } - else + if (!std::isinf(alphaUy)) { - alphaUx = 999; + while (alphaY <= alphaMin + TOL) + alphaY += alphaUy; } - if (fabs(rayVector[1]) > EPSILON /*!= 0*/) + if (!std::isinf(alphaUz)) { - alphaUy = ctPixelSpacing[1] / fabs(rayVector[1]); - } - else - { - alphaUy = 999; - } - if (fabs(rayVector[2]) > EPSILON /*!= 0*/) - { - alphaUz = ctPixelSpacing[2] / fabs(rayVector[2]); - } - else - { - alphaUz = 999; + while (alphaZ <= alphaMin + TOL) + alphaZ += alphaUz; } + /* Calculate voxel index incremental values along the ray path. */ + // Use the direction of the ray + iU = (rayVector[0] > 0.0) ? 1 : -1; + jU = (rayVector[1] > 0.0) ? 1 : -1; + kU = (rayVector[2] > 0.0) ? 1 : -1; - iU = (SourceWorld[0] < drrPixelWorld[0]) ? 1 : -1; - jU = (SourceWorld[1] < drrPixelWorld[1]) ? 1 : -1; - kU = (SourceWorld[2] < drrPixelWorld[2]) ? 1 : -1; + // Initialize line integral + d12 = 0.0; - 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. */ + /* Initialize the current voxel index with the voxel entered at alphaMin. */ 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 */ + // Current parametric position and previous one + double alpha = alphaMin; + double alphaPrev = alphaMin; + + // Siddon traversal + while (alpha < alphaMax) { - /* Store the current ray position */ - alphaCminPrev = alphaCmin; + // Next plane intersection along the ray + double alphaNext = std::min(std::min(alphaX, alphaY), alphaZ); + if (alphaNext > alphaMax) + alphaNext = alphaMax; - if ((alphaX <= alphaY) && (alphaX <= alphaZ)) + double segmentLength = alphaNext - alphaPrev; + // don’t break immediately on a single zero-length segment. This would cause black pixels! + // if (segmentLength <= 0.0) + // break; // FP guard, especially for grid-aligned rays + if (segmentLength <= 0.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; - } - 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(sizeCT[0])) && (cIndex[1] >= 0) && - (cIndex[1] < static_cast(sizeCT[1])) && (cIndex[2] >= 0) && - (cIndex[2] < static_cast(sizeCT[2]))) - { - /* If it is a valid index, get the voxel intensity. */ - value = static_cast(inputPtr->GetPixel(cIndex)); - if (value > m_Threshold) /* Ignore voxels whose intensities are below the threshold. */ + // advance to the next plane and try again + if (nearlyEqual(alphaNext, alphaX)) { - d12 += (alphaCmin - alphaCminPrev) * (value - m_Threshold); + cIndex[0] += iU; + alphaX += alphaUx; } + else if (nearlyEqual(alphaNext, alphaY)) + { + cIndex[1] += jU; + alphaY += alphaUy; + } + else + { + cIndex[2] += kU; + alphaZ += alphaUz; + } + + alphaPrev = alphaNext; + alpha = alphaNext; + + if (alpha >= alphaMax) + break; + + continue; + } + + // Accumulate contribution if inside the volume + if (cIndex[0] >= 0 && cIndex[0] < static_cast(sizeCT[0]) && + cIndex[1] >= 0 && cIndex[1] < static_cast(sizeCT[1]) && + cIndex[2] >= 0 && cIndex[2] < static_cast(sizeCT[2])) + { + value = static_cast(inputPtr->GetPixel(cIndex)); + + if (value > m_Threshold) + { + // physical path length = segmentLength (in α) * |rayVector| + d12 += segmentLength * rayLength * (value - m_Threshold); + } + } + + if (alphaNext >= alphaMax) + break; + + alphaPrev = alphaNext; + alpha = alphaNext; + + // Advance along the axis whose plane we actually hit + if (nearlyEqual(alphaNext, alphaX)) + { + cIndex[0] += iU; + alphaX += alphaUx; + } + else if (nearlyEqual(alphaNext, alphaY)) + { + cIndex[1] += jU; + alphaY += alphaUy; + } + else // must be z + { + cIndex[2] += kU; + alphaZ += alphaUz; } } + // Clamp output value if (d12 < minOutputValue) { pixval = minOutputValue; @@ -483,7 +568,7 @@ gSiddonJacobsRayCastInterpolateImageFunction::Evaluate(c { pixval = static_cast(d12); } - return (pixval); + return pixval; }