diff --git a/Common/CUDA/Siddon_projection.cu b/Common/CUDA/Siddon_projection.cu index 2a025f8c..87e7a56d 100644 --- a/Common/CUDA/Siddon_projection.cu +++ b/Common/CUDA/Siddon_projection.cu @@ -176,9 +176,25 @@ __global__ void kernelPixelDetector( Geometry geo, float aM=fminf(fminf(axM,ayM),azM); // line intersects voxel space -> am=aM) + // The ray never enters the volume. RETURN HERE. Without it, execution + // falls into the index block below, whose imin..kmax are only meaningful + // for a ray that does enter. For a ray the eps guard has snapped to run + // exactly along one voxel axis, with the source OUTSIDE the volume on + // one of its zero axes, the two quotients for that axis share a sign, so + // fminf and fmaxf collapse (axm == axM == -inf) and aM inherits -inf. + // The other zero component then gives aM*ray = (-inf)*0 = NaN, Np is + // NaN, and (unsigned long)NaN is 2^63 on CUDA: that thread iterates + // ~9e18 times and the kernel never returns -- 100% SM, idle memory + // controller, no error, no result. Measured: 2 rays in 261 million on an + // odd x odd detector with an off-axis object; see + // Python/tests/test_siddon_missing_ray.py and + // Frontispiece/siddon_missed_ray_hang.png. A NEGATIVE Np does not do + // this -- CUDA saturates negatives and -inf to 0, a zero-trip loop. + 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 float imin,imax,jmin,jmax,kmin,kmax; diff --git a/Common/CUDA/Siddon_projection_parallel.cu b/Common/CUDA/Siddon_projection_parallel.cu index abd648b3..2f3ea6e2 100644 --- a/Common/CUDA/Siddon_projection_parallel.cu +++ b/Common/CUDA/Siddon_projection_parallel.cu @@ -180,9 +180,17 @@ __global__ void kernelPixelDetector_parallel( Geometry geo, float aM=(fminf(axM,ayM)); // line intersects voxel space -> am=aM) + // See Siddon_projection.cu: without the return a missed ray falls + // through to an Np built from indices that only mean something for a + // hit; for an axis-parallel ray with the source outside the volume that + // is (-inf)*0 = NaN, and (unsigned long)NaN = 2^63 spins forever. The + // hang was NOT reproduced in this kernel (26 configurations tried); the + // return here is by inspection. + 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 float imin,imax,jmin,jmax; diff --git a/Frontispiece/siddon_missed_ray_hang.png b/Frontispiece/siddon_missed_ray_hang.png new file mode 100644 index 00000000..e1acd50a Binary files /dev/null and b/Frontispiece/siddon_missed_ray_hang.png differ diff --git a/Frontispiece/siddon_missed_ray_hang.py b/Frontispiece/siddon_missed_ray_hang.py new file mode 100644 index 00000000..d5327b4b --- /dev/null +++ b/Frontispiece/siddon_missed_ray_hang.py @@ -0,0 +1,197 @@ +# -*- coding: utf-8 -*- +"""The geometry that makes Siddon_projection.cu never return. + +Top: the x-y plane in the kernel's voxel coordinates, drawn with y HORIZONTAL +so the ray (which runs exactly along -Y) lies across the page and the two +X-planes are horizontal lines. Real captured numbers for the ray that hangs +(proj 53, row 619, col 313 of the test geometry). Bottom: why an odd +detector dimension supplies the first zero direction component. +""" +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +from matplotlib.patches import Rectangle, FancyArrowPatch + +INK, INK2, SURF = "#0b0b0b", "#52514e", "#fcfcfb" +BLUE, BLUE_L = "#2a78d6", "#cde2fb" +CRIT, CRIT_L, GOOD, GREY = "#d03b3b", "#fbeaea", "#0ca30c", "#9a9891" + +# captured on the GPU, voxel units. nVoxel is (Z,Y,X)-ordered: X=64, Y=128. +NX, NY = 64.0, 128.0 +S_X, S_Y = 280.44769, 497.73584 # source +P_X, P_Y = 280.44794, -225.03015 # detector pixel + +fig = plt.figure(figsize=(13.5, 9.6), facecolor=SURF) +gs = fig.add_gridspec(2, 1, height_ratios=[2.35, 1.0], hspace=0.30, + left=0.055, right=0.985, top=0.885, bottom=0.04) + +# ================================================================ panel A +# horizontal axis = y (voxel), vertical axis = x (voxel) +ax = fig.add_subplot(gs[0]) +ax.set_facecolor(SURF) +ax.set_aspect("equal") +ax.set_xlim(-300, 600) +ax.set_ylim(-75, 345) + +# volume and the two X-planes Siddon measures distance to +ax.add_patch(Rectangle((0, 0), NY, NX, facecolor=BLUE_L, edgecolor=BLUE, + lw=1.4, zorder=2)) +ax.text(NY / 2, 49, "volume\nx in [0, 64], y in [0, 128]", ha="center", + va="center", fontsize=8.6, color=INK, zorder=3) +for xp, lab in ((0.0, "X-plane x = 0"), (NX, "X-plane x = nVoxelX = 64")): + ax.plot([-300, 600], [xp, xp], color=BLUE, lw=0.9, ls=(0, (4, 3)), + zorder=1, alpha=0.85) + ax.text(-296, xp + 4, lab, ha="left", va="bottom", fontsize=8.5, + color=BLUE) +ax.text(-292, -44, r"Siddon: $a_{xm}, a_{xM}$ = ray parameter where it " + "crosses each X-plane = (plane $-$ source.x) / ray.x", + ha="left", va="top", fontsize=8.5, color=BLUE) + +# --- the ordinary miss: same source, slight tilt -> finite, harmless +ax.plot([S_Y, P_Y], [S_X, 222.0], color=GREY, lw=1.3, ls=(0, (2, 3)), + zorder=3) +ax.text(-292, 214, "ordinary miss: ray.x $\\neq$ 0, quotients finite,\n" + "Np saturates to 0 -- harmless (257 million of these)", + ha="left", va="top", fontsize=8.5, color=INK2) + +# --- the ray that hangs: exactly along -Y at x = 280.45 +ax.plot([S_Y, P_Y], [S_X, P_X], color=CRIT, lw=2.4, zorder=5) +ax.plot(S_Y, S_X, marker="o", ms=9, mfc=CRIT, mec=SURF, mew=1.5, zorder=6) +ax.plot(P_Y, P_X, marker="s", ms=7.5, mfc=SURF, mec=CRIT, mew=1.7, zorder=6) +ax.text(S_Y + 10, S_X + 6, "source (x 280.45, y 497.74)", fontsize=9, + color=INK, va="bottom", ha="left") +ax.text(P_Y - 4, P_X - 9, "detector pixel\n(x 280.45, y $-$225.03)", + fontsize=9, color=INK, va="top", ha="left") +ax.text(200, 300, + "THE RAY THAT HANGS -- exactly along $-$Y\n" + "ray.x = 280.44794 $-$ 280.44769 = 0.00025 < eps $\\to$ 0 " + "ray.z = 63.999996 $-$ 64 = $-$4e$-$6 < eps $\\to$ 0\n" + "runs 216 voxels above the volume and never enters it", + fontsize=9, color=CRIT, va="bottom", ha="center", + bbox=dict(boxstyle="round,pad=0.4", fc=SURF, ec=CRIT, lw=0.9)) + +# both quotients point the SAME way: both planes are below the source +for yy, xp, lab in ((400, NX, r"$(64-280.4)\,/\,0 = -\infty$"), + (330, 0.0, r"$(0-280.4)\,/\,0 = -\infty$")): + ar = FancyArrowPatch((yy, S_X), (yy, xp), arrowstyle="-|>", + mutation_scale=13, color=CRIT, lw=1.4, zorder=4) + ax.add_patch(ar) + ax.text(yy + 7, (S_X + xp) / 2 + 20, lab, ha="left", va="center", + fontsize=9.2, color=CRIT, rotation=90) +ax.text(70, 150, + "source.x = 280.4 is OUTSIDE [0, 64]: both X-planes lie on the\n" + "SAME side, so fminf and fmaxf collapse onto one value,\n" + r"$a_{xm} = a_{xM} = -\infty$ $\Rightarrow$ $a_M = \min(a_{xM}, a_{yM}, a_{zM}) = -\infty$" "\n" + r"$\Rightarrow$ $a_M \cdot$ ray.z $= (-\infty)\cdot 0 =$ NaN $\Rightarrow$ Np = NaN" "\n" + r"$\Rightarrow$ (unsigned long) NaN $= 2^{63}$ loop iterations", + ha="center", va="center", fontsize=9.2, color=INK, + bbox=dict(boxstyle="round,pad=0.5", fc=CRIT_L, ec=CRIT, lw=0.9)) + +# --- the safe contrast: source INSIDE the volume's x-extent +G_X = 32.0 +ax.plot([S_Y, P_Y], [G_X, G_X], color=GOOD, lw=1.8, zorder=4) +ax.plot(S_Y, G_X, marker="o", ms=8, mfc=GOOD, mec=SURF, mew=1.5, zorder=6) +ax.text(S_Y + 10, G_X - 6, "same ray, but with\nsource.x = 32, inside [0, 64]", + fontsize=8.6, color=GOOD, va="top", ha="left") +for yy, xp, lab in ((470, NX, r"$+\infty$"), (470, 0.0, r"$-\infty$")): + ar = FancyArrowPatch((yy, G_X), (yy, xp), arrowstyle="-|>", + mutation_scale=11, color=GOOD, lw=1.3, zorder=5) + ax.add_patch(ar) + ax.text(yy + 6, xp + (5 if xp else -5), lab, ha="left", + va="bottom" if xp else "top", fontsize=9, color=GOOD) +ax.text(-292, -22, + "planes on OPPOSITE sides: $a_{xm} = -\\infty$, $a_{xM} = +\\infty$ " + "stay apart, $a_M$ comes from y or z and is finite. " + "(And this ray hits the volume anyway.)", + ha="left", va="top", fontsize=8.6, color=GOOD) + +for s in ax.spines.values(): + s.set_visible(False) +ax.tick_params(colors=INK2, labelsize=8, length=0) +ax.set_xlabel("y (voxel units) -- the ray runs along $-$Y, right to left", + color=INK2, fontsize=9) +ax.set_ylabel("x (voxel units)", color=INK2, fontsize=9) +ax.grid(True, color="#e6e5e1", lw=0.6) +ax.set_title("A. The x-y plane in kernel voxel coordinates -- real numbers " + "captured from the ray at proj 53, row 619, col 313", + loc="left", fontsize=10.5, color=INK, pad=8) + +# ================================================================ panel B +axb = fig.add_subplot(gs[1]) +axb.set_facecolor(SURF) +axb.set_aspect("equal") +axb.set_xlim(-0.6, 23.4) +axb.set_ylim(-1.2, 7.9) +axb.axis("off") +axb.set_title("B. Why ODD detector dimensions supply the first zero " + "component", loc="left", fontsize=10.5, color=INK, pad=8) + + +def grid(x0, y0, n, label, odd): + for i in range(n): + for j in range(n): + on_row = odd and (j == n // 2) + on_col = odd and (i == n // 2) + fc = "#ffffff" + if on_row or on_col: + fc = CRIT_L + if on_row and on_col: + fc = "#f3bdbd" + axb.add_patch(Rectangle((x0 + i, y0 + j), 1, 1, facecolor=fc, + edgecolor="#bdbcb6", lw=0.7, zorder=2)) + cx, cy = x0 + n / 2.0, y0 + n / 2.0 # optical axis through centre + axb.plot([x0 - 0.35, x0 + n + 0.35], [cy, cy], color=BLUE, lw=1.1, + ls=(0, (4, 3)), zorder=3) + axb.plot([cx, cx], [y0 - 0.35, y0 + n + 0.35], color=BLUE, lw=1.1, + ls=(0, (4, 3)), zorder=3) + axb.plot(cx, cy, marker="+", ms=11, mew=1.6, color=BLUE, zorder=4) + axb.text(x0 + n / 2.0, 7.55, label, ha="center", va="top", fontsize=10, + color=INK) + return cx, cy + + +# align the two grids on a common axis height +cx1, cy1 = grid(0.0, 0.0, 7, "ODD (7 x 7)", True) +cx2, cy2 = grid(9.0, 0.5, 6, "EVEN (6 x 6)", False) + +axb.text(8.0, 6.6, "optical axis\n(dashed)", ha="center", va="center", + fontsize=8.3, color=BLUE) + +# the half-pixel offset on the even grid, drawn outside the crosshair +axb.annotate("", xy=(cx2 + 3.55, cy2), xytext=(cx2 + 3.55, cy2 + 0.5), + arrowprops=dict(arrowstyle="<->", color=INK2, lw=1.0, + shrinkA=0, shrinkB=0)) +axb.text(cx2 + 3.7, cy2 + 0.25, "$\\frac{1}{2}$ px", fontsize=8.5, + color=INK2, va="center", ha="left") +axb.text(cx2, 0.15, "nearest pixel is $\\frac{1}{2}$ px off the axis\n" + "= 0.27 voxels = 270 x eps", ha="center", va="top", fontsize=8.3, + color=INK2) + +axb.text(17.0, 5.9, + "ODD: a pixel sits EXACTLY on the axis, so the whole centre\n" + "row has ray.z = 0. That is 1179 columns x 179 views = 211 000\n" + "rays on which a second coincidence (ray.x = 0) can land.\n" + "Both hanging rays were on row 619 = (1239 $-$ 1) / 2.", + ha="left", va="top", fontsize=8.7, color=CRIT) +axb.text(17.0, 2.5, + "EVEN: no pixel on the axis. The nearest is half a pixel\n" + "off -- 0.21 mm, 0.27 voxels, 270 x above eps = 0.001 --\n" + "so ray.z is never snapped to 0 and nothing collapses.\n" + "Adding ONE pixel to either dimension is enough to fix it.", + ha="left", va="top", fontsize=8.7, color=GOOD) + +fig.text(0.055, 0.955, + "Siddon forward projection never returns: a missed ray falls " + "through to a loop bound of (unsigned long) NaN = 2$^{63}$", + fontsize=13, color=INK, weight="bold") +fig.text(0.055, 0.918, + "All four at once: two ray components under eps (the ray lies along " + "a voxel axis) + source outside the volume along a zeroed axis " + "+ the ray misses + no early return. 2 rays in 261 million.", + fontsize=9.0, color=INK2) + +import os +out = os.path.join(os.path.dirname(os.path.abspath(__file__)), + "siddon_missed_ray_hang.png") +fig.savefig(out, dpi=170, facecolor=SURF) +print("wrote", out) diff --git a/Python/tests/test_siddon_missing_ray.py b/Python/tests/test_siddon_missing_ray.py new file mode 100644 index 00000000..39c0a803 --- /dev/null +++ b/Python/tests/test_siddon_missing_ray.py @@ -0,0 +1,225 @@ +# -*- coding: utf-8 -*- +"""Siddon forward projection must TERMINATE when rays miss the volume. + +Regression test for a hang in `Common/CUDA/Siddon_projection.cu`, which is +compiled into both the Python extension and the MATLAB mex (see +`MATLAB/Compile.m`), so this covers both bindings. + +THE BUG + + // line intersects voxel space -> am=aM) + detector[idx]=0; // <-- no return; execution falls through + +`am>=aM` means the ray never enters the volume. Without a `return`, +execution continues to + + unsigned long Np=(imax-imin+1)+(jmax-jmin+1)+(kmax-kmin+1); + for (unsigned long ii=0; ii=aM){ + detector[idx]=0; + return; + } + +Verified value-neutral: every configuration that returned before the fix +returns bit-identical sums after it. + +THE PARALLEL KERNEL IS NOT COVERED HERE, AND NOT FOR WANT OF TRYING + +`Siddon_projection_parallel.cu` carries the identical hazard -- + + unsigned long Np=(imax-imin+1)+(jmax-jmin+1); + +-- and was given the same `return`, but the hang could NOT be reproduced in +it. 26 configurations on a pre-fix build all returned in ~0.75 s: both +detector parities, offOrigin from +-5 mm to +-2000 mm, spanning the band +that hangs the cone kernel and well past it. Dispatch was confirmed genuine +(`_Ax.pyx` sends mode=="parallel" to siddon_ray_projection_parallel, and the +sums differ from cone mode), so the parallel path really was exercised. + +Two plausible reasons it resists, neither established: + * only TWO terms sum into Np instead of three, so there is less room to + go negative; + * that kernel already carries defensive code the cone one lacks -- an + epsilon on the divisions and explicit copysignf(1e11,...) handling for + ray.x/ray.y == 0 -- added, per its own comment, because "in paralel + beam, often ray.y or ray.x=0; this leads to infinities propagating and + breaking everything". It has been hardened once already. + +So the parallel fix is justified BY INSPECTION and is UNVERIFIED. Do not +read the cone tests below as covering it. If you find a parallel trigger, +add it here; the search above is recorded so it need not be repeated. + +Each projection runs in a CHILD PROCESS with a timeout, so a regression +FAILS this test rather than hanging the suite. +""" +import multiprocessing as mp + +import numpy as np + +TIMEOUT_S = 90 + +# MUST match the configuration the hang was verified on. The offsets and the +# per-view DSD/DSO are interpolated across this many views, so changing it +# resamples the trajectory and the grazing view that trips the bug may not +# be visited at all: at 60 views this test PASSES on a known-broken build. +# Do not lower it to make the test faster. +N_ANGLES = 179 + + +def _geometry(n_det, off_mm): + import tigre + geo = tigre.geometry(mode="cone", default=True, + nVoxel=np.array([128, 128, 64])) + geo.sVoxel = np.array([100.0, 125.0, 40.0]) + geo.dVoxel = geo.sVoxel / geo.nVoxel + geo.nDetector = np.asarray(n_det) + geo.dDetector = np.array([0.42, 0.42]) + geo.sDetector = geo.nDetector * geo.dDetector + geo.DSD = np.linspace(580.0, 920.0, N_ANGLES) + geo.DSO = np.linspace(407.0, 514.0, N_ANGLES) + geo.offOrigin = np.zeros((N_ANGLES, 3)) + geo.offOrigin[:, 2] = np.linspace(-off_mm, off_mm, N_ANGLES) + geo.offDetector = np.zeros((N_ANGLES, 2)) + geo.offDetector[:, 1] = np.linspace(-1.7 * off_mm, 1.7 * off_mm, N_ANGLES) + angles = np.linspace(0.9163, 2.2253, N_ANGLES).astype(np.float32) + return geo, angles + + +def _project(n_det, off_mm, q): + import tigre + geo, angles = _geometry(n_det, off_mm) + img = np.zeros(tuple(int(v) for v in geo.nVoxel), dtype=np.float32) + img[32:96, 32:96, 16:48] = 1.0 + proj = tigre.Ax(img, geo, angles) + q.put((float(np.asarray(proj).sum()), + bool(np.all(np.isfinite(np.asarray(proj)))))) + + +def _run(n_det, off_mm): + """Project in a child process. Returns (sum, all_finite) or None on + timeout -- None IS the failure this test exists to catch.""" + ctx = mp.get_context("spawn") + q = ctx.Queue() + p = ctx.Process(target=_project, args=(n_det, off_mm, q)) + p.start() + p.join(TIMEOUT_S) + if p.is_alive(): + p.terminate() + p.join() + return None + return None if q.empty() else q.get() + + +def test_missed_rays_terminate_with_odd_detector(): + """The reported case: both detector dimensions odd, object off-axis.""" + for off_mm in (50.0, 100.0, 200.0): + got = _run(np.array([1239, 1179]), off_mm) + assert got is not None, ( + f"Ax did not return within {TIMEOUT_S}s at offOrigin +-{off_mm} mm " + f"with a 1239x1179 detector -- the missed-ray fall-through in " + f"Siddon_projection.cu has regressed") + total, finite = got + assert finite, f"non-finite projection values at +-{off_mm} mm" + assert total > 0.0, f"empty projection at +-{off_mm} mm: {total}" + + +def test_detector_parity_does_not_change_the_result(): + """Padding the detector by one pixel must not change what is measured. + + This is the other half: the fix has to make the hang go away WITHOUT + moving any value. Odd and even detectors are different samplings, so + the totals differ slightly -- but only by the added row/column, not by + a factor. + """ + ref = _run(np.array([1240, 1180]), 100.0) + assert ref is not None, "even x even detector did not return" + for n_det in ([1239, 1179], [1239, 1180], [1240, 1179]): + got = _run(np.array(n_det), 100.0) + assert got is not None, f"{n_det} did not return within {TIMEOUT_S}s" + rel = abs(got[0] - ref[0]) / ref[0] + assert rel < 0.01, (f"{n_det} total {got[0]:.6g} differs from " + f"{ref[0]:.6g} by {rel:.3%} -- more than the one " + f"pixel of extra detector can account for") + + +def test_centred_object_is_unaffected(): + """A centred object never trips the bug and must be untouched by the + fix -- the control that catches a 'fix' which changes physics.""" + got = _run(np.array([1239, 1179]), 0.0) + assert got is not None, "centred case did not return" + assert got[1] and got[0] > 0.0 + + +if __name__ == "__main__": + for fn in (test_missed_rays_terminate_with_odd_detector, + test_detector_parity_does_not_change_the_result, + test_centred_object_is_unaffected): + print(f"{fn.__name__} ... ", end="", flush=True) + try: + fn() + print("PASS") + except AssertionError as e: + print(f"FAIL\n {e}")