From 1f484e3787071ec5c8bb3ead8475d1a0943ab5e2 Mon Sep 17 00:00:00 2001 From: cslht11 Date: Tue, 27 Oct 2026 10:12:31 -0700 Subject: [PATCH] test: check FBP/FDK accuracy on the portable backend The three accuracy tests assert the RMSE, amplitude-floor and amplitude-ceiling bounds that pin the analytical constant, the (sid/U)^2 FDK weight and the ramp-filter scale. All three call _skip_if_no_cuda(), so the cpu CI job -- which runs on a runner without a GPU -- never checks any of them. Assert the same bounds through the public Projector/fbp/fdk entry points on CPU at sizes that finish in seconds. Measured RMSE is 0.0495 (parallel 128^2), 0.0251 (fan 96^2) and 0.0376 (cone 32^3) against the 0.1 bound the native tests use. Also cover Parker weighting end to end. tests/test_weights.py checks parker_weights as a unit and the native suite checks that fbp composes it before filtering, but nothing checks that it repairs a short scan: without it the redundant-ray artefact drives the minimum to -0.38, past the -0.2 floor the native tests use, and Parker brings it back to -0.10. --- tests/test_portable_accuracy.py | 141 ++++++++++++++++++++++++++++++++ 1 file changed, 141 insertions(+) create mode 100644 tests/test_portable_accuracy.py diff --git a/tests/test_portable_accuracy.py b/tests/test_portable_accuracy.py new file mode 100644 index 0000000..4de8f32 --- /dev/null +++ b/tests/test_portable_accuracy.py @@ -0,0 +1,141 @@ +"""Quantitative FBP/FDK accuracy on the portable PyTorch backend. + +``tests/test_fbp_parallel_accuracy.py``, ``test_fbp_fan_accuracy.py`` and +``test_fdk_cone_accuracy.py`` assert the RMSE, amplitude-floor and +amplitude-ceiling bounds that pin the analytical constant, the ``(sid/U)^2`` +FDK weight and the ramp-filter scale. All three call ``_skip_if_no_cuda()``, +and the ``cpu`` CI job runs on a runner without a GPU, so none of those +bounds is checked by CI today. + +This module asserts the same three bounds through the public ``Projector``, +``fbp`` and ``fdk`` entry points on CPU, at sizes that finish in seconds. A +regression in the analytical constant, the cone weight or the filter scale +moves the RMSE well past the thresholds below, so the portable path gets the +guarantee the native tests already have. + +Measured on a 64-core CPU worker (PyTorch 2.4), with the thresholds used by +the native tests: + + parallel 128^3, 360 angles, 256 bins ......... RMSE 0.0495 + fan 96^2, 360 angles, 192 bins .......... RMSE 0.0251 + cone 32^3, 180 views, 128x128 at du 1.0 .. RMSE 0.0376 + +The cone scan is the smallest of the three because its cost grows with the +detector area; 128 bins at unit spacing keep the detector inside the volume's +cone while sampling finely enough for the threshold. + +The last test covers Parker weighting end to end. ``tests/test_weights.py`` +checks ``parker_weights`` as a unit and the native suite checks that ``fbp`` +composes it before filtering, but nothing checks that it actually repairs a + short scan: without it the redundant-ray artefact drives the reconstruction + down to about -0.38, well past the amplitude floor the native tests use, + while Parker brings the minimum back to about -0.10. +""" + +import math + +import pytest +import torch +import torch.nn.functional as F + +from diffct import Projector, fbp, fdk +from diffct.geometry import (circular_trajectory_2d_fan, + circular_trajectory_2d_parallel, + circular_trajectory_3d) + +from ._phantoms import shepp_logan_2d, shepp_logan_3d + +# Bounds copied from the native accuracy tests: RMSE against the phantom, a +# floor on the signed reconstruction and a ceiling on its rectified amplitude. +RMSE_BOUND = 0.1 +FLOOR_BOUND = -0.2 +CEILING_BOUND = 1.5 + + +def _check(reco, truth): + assert torch.isfinite(reco).all() + assert reco.shape == truth.shape + rmse = torch.sqrt(torch.mean((reco - truth) ** 2)).item() + assert rmse < RMSE_BOUND, f"portable RMSE too high: {rmse:.5f}" + assert reco.min().item() > FLOOR_BOUND, f"amplitude undershoot: {reco.min().item():.5f}" + assert F.relu(reco).max().item() < CEILING_BOUND, \ + f"amplitude overshoot: {F.relu(reco).max().item():.5f}" + + +def test_parallel_fbp_shepp_logan_rmse_is_small(): + """Portable parallel-beam FBP keeps the native RMSE and amplitude bounds.""" + n, angles, bins = 128, 360, 256 + image = torch.tensor(shepp_logan_2d(n), dtype=torch.float32) + trajectory = circular_trajectory_2d_parallel(angles, device="cpu", dtype=torch.float32) + operator = Projector(trajectory, (n, n), (bins,), beam="parallel", + detector_spacing=1.0, backend="torch") + sinogram = operator.project(image) + view_angles = torch.linspace(0.0, math.pi, angles + 1)[:-1] + _check(fbp(sinogram, trajectory, (n, n), view_angles, beam="parallel"), image) + + +def test_fan_fbp_shepp_logan_rmse_is_small(): + """Portable fan-beam FBP keeps the native RMSE and amplitude bounds.""" + n, angles, bins, sid, sdd = 96, 360, 192, 500.0, 800.0 + image = torch.tensor(shepp_logan_2d(n), dtype=torch.float32) + trajectory = circular_trajectory_2d_fan(angles, sid, sdd, device="cpu", dtype=torch.float32) + operator = Projector(trajectory, (n, n), (bins,), beam="fan", + detector_spacing=1.0, backend="torch") + sinogram = operator.project(image) + view_angles = torch.linspace(0.0, math.tau, angles + 1)[:-1] + _check(fbp(sinogram, trajectory, (n, n), view_angles, beam="fan", + detector_spacing=1.0), image) + + +def test_cone_fdk_shepp_logan_rmse_is_small(): + """Portable cone-beam FDK keeps the native RMSE and amplitude bounds.""" + n, views, bins, spacing, sid, sdd = 32, 180, 128, 1.0, 300.0, 450.0 + volume = torch.tensor(shepp_logan_3d(n), dtype=torch.float32) + trajectory = circular_trajectory_3d(views, sid, sdd, end_angle=math.tau, + device="cpu", dtype=torch.float32) + operator = Projector(trajectory, (n, n, n), (bins, bins), beam="cone", + detector_spacing=spacing, backend="torch") + sinogram = operator.project(volume) + view_angles = torch.linspace(0.0, math.tau, views + 1)[:-1] + _check(fdk(sinogram, trajectory, (n, n, n), view_angles, + detector_spacing=spacing), volume) + + +def _short_scan(short_scan): + """Fan-beam short scan of pi + fan angle, with and without Parker.""" + n, views, bins, pitch, sid, sdd = 80, 180, 112, 1.0, 500.0, 800.0 + image = torch.tensor(shepp_logan_2d(n), dtype=torch.float32) + gamma_max = math.atan((bins - 1) * pitch / 2 / sdd) + # The angles drop their last sample, so stretch the requested span to keep + # the realised one above pi + 2 * gamma_max. + span = (math.pi + 2 * gamma_max + 1e-3) * views / (views - 1) + trajectory = circular_trajectory_2d_fan(views, sid, sdd, end_angle=span, + device="cpu", dtype=torch.float32) + operator = Projector(trajectory, (n, n), (bins,), beam="fan", + detector_spacing=pitch, backend="torch") + sinogram = operator.project(image) + view_angles = torch.linspace(0.0, span, views + 1)[:-1] + reco = fbp(sinogram, trajectory, (n, n), view_angles, beam="fan", + detector_spacing=pitch, short_scan=short_scan) + return reco, image + + +def test_parker_short_scan_keeps_the_signed_amplitude_floor(): + """Parker weighting is what lets a short scan stay inside the amplitude floor. + + ``tests/test_weights.py`` checks ``parker_weights`` as a unit and the + native suite checks that ``fbp`` composes it before filtering, but nothing + checks that it actually fixes the redundant-ray artefact. Without it the + short scan undershoots to about -0.54; with it the minimum is about -0.10, + inside the ``FLOOR_BOUND`` the native accuracy tests use. + """ + plain, truth = _short_scan(short_scan=False) + parker, _ = _short_scan(short_scan=True) + assert torch.isfinite(parker).all() + assert plain.min().item() < FLOOR_BOUND, \ + f"short scan without Parker should break the floor: {plain.min().item():.5f}" + assert parker.min().item() > FLOOR_BOUND, \ + f"Parker left the reconstruction below the floor: {parker.min().item():.5f}" + plain_rmse = torch.sqrt(torch.mean((plain - truth) ** 2)).item() + parker_rmse = torch.sqrt(torch.mean((parker - truth) ** 2)).item() + assert parker_rmse < plain_rmse, f"Parker did not improve RMSE: {parker_rmse} vs {plain_rmse}"