trilinear interpolation

Ray geometry and interpolation fixed.

reg23Topograms/itkDTRrecon/itkgSiddonJacobsRayCastInterpolateImageFunction.hxx

@@ -48,6 +48,7 @@ CIT 6, 89-94 (1998).
 #define itkgSiddonJacobsRayCastInterpolateImageFunction_hxx
 
 #include "itkgSiddonJacobsRayCastInterpolateImageFunction.h"
+#include "itkContinuousIndex.h"
 
 #include "itkMath.h"
 #include <cstdlib>
@@ -128,7 +129,7 @@ gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(c
     float alphaX, alphaY, alphaZ, alphaCmin, alphaCminPrev;
     float alphaUx, alphaUy, alphaUz;
     float alphaIntersectionUp[3], alphaIntersectionDown[3];
-    float d12, value;
+    float d12, value,valuetril;
     float firstIntersectionIndex[3];
     int   firstIntersectionIndexUp[3], firstIntersectionIndexDown[3];
     int   iU, jU, kU;
@@ -139,7 +140,6 @@ gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(c
     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();
@@ -156,9 +156,12 @@ gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(c
     }
 
     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();
@@ -182,6 +185,7 @@ gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(c
     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.;
@@ -431,6 +435,7 @@ gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(c
     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 */
@@ -458,17 +463,104 @@ gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(c
             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])))
+            (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. */
+            // 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. */
             {
-                d12 += (alphaCmin - alphaCminPrev) * (value - m_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)