From d23622360f03ee1e7296534218c51d6e36e7a773 Mon Sep 17 00:00:00 2001 From: Koh Hashimoto Date: Mon, 5 Oct 2026 21:15:55 +0900 Subject: [PATCH] Fix Siddon projection kernel missing return on missed rays (cone and parallel) Add early return after zeroing detector on missed rays (am >= aM) to avoid invalid index ranges and unnecessary computation in subsequent steps. Add regression test test_siddon_missed_rays.py covering cone and parallel modes. --- Common/CUDA/Siddon_projection.cu | 4 ++- Common/CUDA/Siddon_projection_parallel.cu | 4 ++- Python/tests/test_siddon_missed_rays.py | 41 +++++++++++++++++++++++ 3 files changed, 47 insertions(+), 2 deletions(-) create mode 100644 Python/tests/test_siddon_missed_rays.py diff --git a/Common/CUDA/Siddon_projection.cu b/Common/CUDA/Siddon_projection.cu index 2a025f8c..b4cdfdbb 100644 --- a/Common/CUDA/Siddon_projection.cu +++ b/Common/CUDA/Siddon_projection.cu @@ -176,8 +176,10 @@ __global__ void kernelPixelDetector( Geometry geo, float aM=fminf(fminf(axM,ayM),azM); // line intersects voxel space -> am=aM) + if (am>=aM) { detector[idx]=0; + return; + } // Compute max/min image INDEX for intersection eq(11-19) // Discussion about ternary operator in CUDA: https://stackoverflow.com/questions/7104384/in-cuda-why-is-a-b010-more-efficient-than-an-if-else-version diff --git a/Common/CUDA/Siddon_projection_parallel.cu b/Common/CUDA/Siddon_projection_parallel.cu index abd648b3..af69ccb8 100644 --- a/Common/CUDA/Siddon_projection_parallel.cu +++ b/Common/CUDA/Siddon_projection_parallel.cu @@ -180,8 +180,10 @@ __global__ void kernelPixelDetector_parallel( Geometry geo, float aM=(fminf(axM,ayM)); // line intersects voxel space -> am=aM) + if (am>=aM) { detector[idx]=0.0f; + return; + } // Compute max/min image INDEX for intersection eq(11-19) // Discussion about ternary operator in CUDA: https://stackoverflow.com/questions/7104384/in-cuda-why-is-a-b010-more-efficient-than-an-if-else-version diff --git a/Python/tests/test_siddon_missed_rays.py b/Python/tests/test_siddon_missed_rays.py new file mode 100644 index 00000000..09aafc6e --- /dev/null +++ b/Python/tests/test_siddon_missed_rays.py @@ -0,0 +1,41 @@ +import numpy as np +import pytest +import tigre + + +@pytest.mark.parametrize("mode", ["cone", "parallel"]) +def test_siddon_missed_rays(mode): + """Test that Siddon projection handles missed rays properly for cone and parallel geometry.""" + geo = tigre.geometry_default(high_resolution=False) + geo.mode = mode + # Make detector large enough so that outer rays do not intersect the voxel volume + geo.nDetector = np.array([256, 256]) + geo.dDetector = np.array([3.0, 3.0]) + geo.sDetector = geo.nDetector * geo.dDetector + + img = np.ones(tuple(geo.nVoxel.astype(int)), dtype=np.float32) + angles = np.array([0.0, 0.5], dtype=np.float32) + + proj = tigre.Ax(img, geo, angles, "Siddon") + + # Missed rays must be exactly 0 + assert not np.isnan(proj).any(), f"NaNs detected in {mode} Siddon projection" + assert (proj == 0.0).sum() > 0, f"Expected missed rays in {mode} Siddon projection" + # Intersecting rays must be positive + assert (proj > 0.0).sum() > 0, f"Expected intersecting rays in {mode} Siddon projection" + + +@pytest.mark.parametrize("mode", ["cone", "parallel"]) +def test_siddon_completely_missed_volume(mode): + """Test that Siddon projection returns all zeros when the volume is completely outside the beam.""" + geo = tigre.geometry_default(high_resolution=False) + geo.mode = mode + # Offset the volume far away along Y (perpendicular to optical axis) + geo.offOrigin = np.array([0.0, 5000.0, 0.0]) + + img = np.ones(tuple(geo.nVoxel.astype(int)), dtype=np.float32) + angles = np.array([0.0], dtype=np.float32) + + proj = tigre.Ax(img, geo, angles, "Siddon") + + assert np.all(proj == 0.0), f"Expected all zeros for completely missed volume in {mode} Siddon"