Conversation
|
Hmm can you explain in more detail why/when this is failing? I don't see this PR causing any negative consequences anyway, but I would like to understand exactly what is going on. |
`if (am>=aM) detector[idx]=0;` has no return, so a ray that never enters
the volume falls through to the index block below it and on into
unsigned long Np=(imax-imin+1)+(jmax-jmin+1)+(kmax-kmin+1);
for (unsigned long ii=0; ii<Np; ii++){ ... }
whose imin..kmax only mean something for a ray that does enter. Measured,
by instrumenting the kernel to dump its intermediates for the ray that
hangs (proj 53, row 619, col 313 of the test geometry):
1. the eps guard snaps ray.x = 0.00025 and ray.z = -4e-6 to exactly 0,
so the ray lies along -Y to ~1e-6 rad;
2. source.x = 280.4 is OUTSIDE nVoxelX = 64, so -source.x/0 and
(nVoxelX-source.x)/0 are BOTH -inf; fminf and fmaxf collapse and
axm == axM == -inf;
3. aM = fminf(axM, ayM, azM) = -inf while am = 0.51 -- the miss is
detected correctly;
4. without the return the Z block runs anyway:
kmin = ceilf(source.z + aM*ray.z), and (-inf)*0 = NaN;
5. Np = NaN, and on CUDA (unsigned long)NaN = 2^63. That thread iterates
~9e18 times and the kernel never returns.
No error, no warning, no result. On the GPU: 100% SM occupancy with an
IDLE memory controller and no PCIe traffic, because one or two threads are
re-fetching one clamped texture address while every other thread finished.
It is 2 rays in 261,479,799 in the test geometry. Two production runs sat
in it for 71 minutes against the 0.9 s the same call takes when it works.
A NEGATIVE Np does not do this: measured, CUDA saturates negatives and
-inf to 0, a zero-trip loop. It is specifically NaN.
WHEN. One ray needs all of: two direction components under eps (the ray
lies along a voxel axis); the source OUTSIDE the volume along a zeroed
axis (inside, the quotients are -inf and +inf, stay apart, and are
harmless); and the ray misses. Why both nDetector entries odd: an odd
count puts a pixel EXACTLY on the detector's centre line, so the whole
centre row has ray.z == 0 -- 211,000 candidate rays here, and both hanging
rays sat on row 619 = (1239-1)/2. An even count straddles the centre by
half a pixel (0.27 voxels, 270x eps) and never snaps. Adding one pixel to
either axis removes the hang; removing one re-creates it. The object
offset then decides whether the second coincidence lands: centred does not
hang, +-50/100/200 mm hangs, +-300/433 mm does not -- a band, not a
threshold.
Circular CT is not immune in principle, but off-axis slab geometry with a
large detector sweeps ray directions widely enough that the ~1e-6 rad
coincidence actually lands, which is why it surfaces in tomosynthesis and
laminography. A stock centred cone-beam geometry with an odd x odd
detector did NOT reproduce it.
Verified on A6000/sm_86 and RTX 6000 Ada/sm_89 over 21 measurements,
flipping parity in both directions. Value-neutral: every configuration
that returned before returns BIT-IDENTICAL sums after (1.68256e+08,
4.18034e+07, 2.89213e+07 on the test geometry), and a 184-test downstream
suite is unchanged.
Common/CUDA/Siddon_projection.cu is compiled into BOTH the Python
extension and the MATLAB mex (MATLAB/Compile.m), so this fixes both
bindings; the only difference in that file is an IS_FOR_MATLAB_TIGRE ifdef
transposing the detector index. Siddon_projection_parallel.cu carries the
identical pattern and gets the same return, but the hang could NOT be
reproduced there (26 configurations); that half is by inspection and the
test docstring says so. The ray_interpolated_* kernels step a fixed sample
count and are unaffected.
Python/tests/test_siddon_missing_ray.py runs each projection in a child
process under a timeout, so a regression FAILS rather than hanging the
suite. It fails on a pre-fix build and passes on a fixed one, with a
centred-object control that passes on both. N_ANGLES is pinned at 179 and
must not be lowered: the offsets and per-view DSD/DSO are interpolated
across the view count, and at 60 views the test passes on a known-broken
build.
Frontispiece/siddon_missed_ray_hang.png draws the ray that hangs in the kernel's voxel coordinates, using the real numbers captured from the GPU (proj 53, row 619, col 313 of the test geometry): both X-plane quotients pointing the same way because the source sits outside the volume on a zeroed axis, against the safe case with the source inside; and why an odd detector dimension puts a pixel exactly on the axis. The generator is alongside so the numbers can be regenerated. Kept as its own commit so the binary can be dropped without touching the fix.
65838ff to
a4f9600
Compare
|
Sure — and digging in, my code comment turned out to be wrong about the mechanism, so I've force-pushed a correction (7549d41) along with a figure. The chain, measured by instrumenting the kernel to dump its intermediates for the ray that hangs (proj 53, row 619, col 313 in the test geometry): The eps guard snaps ray.x = 0.00025 and ray.z = −4e-6 to exactly 0 — the ray lies along −Y to ~1e-6 rad. [Image] Top: the x–y plane in voxel coordinates with the real captured numbers. The two red arrows pointing the same way are the whole bug — both X-planes lie on one side of the source, so the two quotients share a sign and the min/max collapse. Green is the same ray with the source inside the volume's x-extent: opposite sides, infinities stay apart, harmless. Bottom: an odd detector count puts a pixel exactly on the axis, so a whole row has ray.z = 0. When. One ray needs all of: two direction components under eps; the source outside the volume along a zeroed axis; and a miss. It's 2 rays in 261 million here, which is why the symptom is so odd — nvidia-smi shows 100% SM with an idle memory controller and no PCIe, because two threads are re-reading one clamped texture address while everything else finished. Why odd × odd. Both NaN rays sit on row 619 = (1239−1)/2, the centre row. An odd count puts a pixel exactly on the centre line, so the whole row has ray.z = 0 — 211k candidates for the second coincidence. An even count straddles by half a pixel (0.27 voxels, 270× eps) and never snaps. Adding one pixel to either axis fixes it; removing one re-breaks it. The object offset then decides whether the second coincidence lands: centred doesn't hang, ±50/100/200 mm does, ±300/433 mm doesn't — a band, not a threshold. Scope. Circular CT isn't immune in principle, but off-axis slab geometry with a large detector sweeps ray directions widely enough that the ~1e-6 rad coincidence actually lands — a stock centred cone-beam with an odd×odd detector did not reproduce it. Siddon_projection_parallel.cu has the same pattern and gets the same return, but I couldn't reproduce a hang there (26 configurations); that half is by inspection and the docstring says so. Verified on A6000/sm_86 and RTX 6000 Ada/sm_89; bit-identical sums for every configuration that previously returned; 184-test downstream suite unchanged. The figure and its generator are a separate commit (a4f9600) — happy to drop the binary if you'd rather not carry it. |
2 participants
What is wrong. In
Common/CUDA/Siddon_projection.cu:am >= aMmeans the ray never enters the volume. Without areturnthe kernel continues towith
imin/imaxinverted for that ray. The sum is a negative float; cast tounsigned longit wraps to ~1.8e19, and the loop runs effectively forever. No error, no warning, no result. On the GPU it reads as 100 % SM occupancy with an idle memory controller and no PCIe traffic (every iteration re-fetches the same clamped texture address). Two of our production runs sat in it for 71 minutes before anyone compared against the 0.9 s the same call takes when it works.When it trips. Two conditions together, both of which ordinary circular CT avoids by accident:
nDetectorentries odd: the detector origin isP.y = dU*(-(nDetecU/2.0)+0.5),P.z = dV*((nDetecV/2.0)-0.5), so an odd count shifts that axis by half a pixel, and odd on both nudges grazing rays onto the degenerate side;Off-axis geometry is ordinary in tomosynthesis and laminography and rare in circular CT, which is where the blind spot comes from. A stock centred cone geometry with an odd x odd detector does NOT reproduce it. Verified on two machines and two GPU architectures (A6000 / sm_86, RTX 6000 Ada / sm_89) over 21 measurements, flipping parity in both directions: adding one pixel to either detector axis rescued a hang, removing one induced it.
Fix.
returnafter zeroing the pixel, inSiddon_projection.cuandSiddon_projection_parallel.cu. Value-neutral: every configuration that returned before returns bit-identical sums after.Scope, stated precisely. The cone kernel (
Siddon_projection.cu) is reproduced, fixed and covered by the test.Siddon_projection_parallel.cucarries the identical unsigned-wrap hazard and gets the samereturn, but the hang could not be reproduced in it: 26 configurations on a pre-fix build all returned, over both detector parities andoffOriginfrom 5 to 2000 mm. That half is justified by inspection only (it may be protected by having two terms inNp, or by thecopysignf(1e11, ...)guards it already carries). Theray_interpolated_*kernels step a fixed sample count and are unaffected.Siddon_projection.cuis compiled into both the Python extension and the MATLAB mex (Compile.m), so the fix covers both bindings; the file's only binding difference is theIS_FOR_MATLAB_TIGREdetector-index transpose.Test (
Python/tests/test_siddon_missing_ray.py). Each projection runs in a child process under a timeout, so a regression FAILS instead of hanging the suite. Three tests: the odd x odd detector with the object displaced across the band must return (fails on a pre-fix build by timeout, passes fixed); detector parity must not change the projection sum (even x even reference vs the three odd variants, < 1 %); a centred object is unaffected.N_ANGLESis pinned at 179: the offsets and per-view DSD/DSO are interpolated across the view count and at 60 views the test passes on a known-broken build.Before / after, on current master. The hang is build-dependent: the same kernel source compiled with nvcc on Linux trips it, the Windows (MSVC + CUDA 13.3) build of the same source does not - which is why a reproducer is hard to hand over and why the test runs under a timeout.
test_missed_rays_terminate...test_detector_parity...test_centred_object...[1239, 1179]did not return within 90 sThe Linux "without the fix" row is the fork's build of the identical
Siddon_projection.cu(byte-identical to current master), GPUs idle, no co-tenants.