Python: filter FDK/FBP projections in cache-sized blocks, on several threads

[Dev-next-gen]

The filtering step of algs.fdk / algs.fbp runs on one core and works on whole padded projections. #331 already noted that it "could be improved by multithreading" (it was 50 s of a 91 s FDK there), and #312 / #423 tried to move it to the GPU without a clear win. This PR keeps it on the CPU, in numpy and scipy.fft as it is now, and makes it 5 times faster on 4 cores and 14 times faster on 16 for a 1024 × 1024 detector. The filtered projections are bit-identical to what master produces.

What changes

filtering() packs two projections in one complex array of nDetector[0] × filt_len and runs fft, the product with the filter and ifft over the whole array, each of them allocating another array of that size. For a 1024 × 1024 detector that is 32 MiB per array, so the data does not stay in cache, and only one core works.

Detector rows are filtered independently of each other. Each pair of projections is now cut in blocks of rows that take about 2 MiB once padded, a block is filtered in place in a buffer owned by the thread, and the blocks are shared among a ThreadPoolExecutor (numpy and scipy.fft release the GIL). The rest is as before: two projections packed in a complex array, the same filter, the same scipy.fft calls on each row. No new dependency, and only Python/tigre/utilities/filtering.py is touched.

Same output

The arithmetic done on a row does not change. The one thing that could change the last bit is how scipy's FFT groups rows in a SIMD vector, so a block is always a multiple of 16 rows, which keeps the grouping it has on a whole projection (with blocks of 5 rows I measured a relative difference of 3e-7).

I compared master's filtering() and this one with np.array_equal in 574 cases: called directly (9 detector shapes, odd and even numbers of angles, the 6 filter settings, 1, 2, 3, 5 and the default number of workers, the default and the smallest block size, float64 input), and through algs.fdk (centred and offset detector, anisotropic pixels) and algs.fbp, comparing the array they hand to Atb. All 574 are bit-identical, with NumPy 2.5.3 / SciPy 1.18.1 and with NumPy 1.26.4 / SciPy 1.11.4, on x86-64.

Timings

filtering(proj, geo, angles, parker=False) on random float32 projections, best of 3, 2 × 20-core Xeon E5-2698 v4, Python 3.12, NumPy 2.5.3, SciPy 1.18.1. "before" is master's function loaded in the same process.

projections	before	1 thread	2	4	8	16 (default)
360 × 512 × 512	2.37 s	2.19 s	1.15 s	0.67 s	0.38 s	0.25 s (×9.5)
360 × 1024 × 1024	14.3 s	9.2 s	5.2 s	2.7 s	1.58 s	1.00 s (×14)
100 × 2048 × 2048	21.0 s	11.3 s	6.1 s	3.3 s	1.93 s	1.22 s (×17)

The machine is shared and the "before" time of the 1024 × 1024 case moved between 14 and 30 s from one run to the next; the table is the run where it was lowest. The "after" times were stable (1.0 s on 16 threads in every run).

Small problems do not get slower: 100 × 64 × 64 goes from 12 to 9 ms, 100 × 256 × 256 from 179 to 27 ms, and a fan-beam sinogram of 720 × 1 × 1024 from 48 to 41 ms (single rows are too small to share out, so that one barely moves).

algs.fdk() up to the call to Atb, that is the weights plus the filtering: 15.9 s → 2.25 s for 360 × 1024 × 1024, 23.1 s → 2.8 s for 100 × 2048 × 2048.

Memory allocated by the step on top of the projections: 80 MiB → 35 MiB for a 1024 × 1024 detector and 320 MiB → 36 MiB for 2048 × 2048 with 16 threads (about 2 MiB per thread, 2.6 MiB with one thread).

End to end on a GPU

On an RTX 3070 Laptop (CUDA 12.4), I timed the full algs.fdk() — the filtering plus the real GPU backprojection — with master's filtering and with this one, in the same process. The reconstructed volume is bit-identical between the two (equal sha256), so the change leaves the reconstruction untouched and only saves time. 360 angles, 1024 detector: a 512³ volume goes 13.22 s → 6.47 s (×2.04), a 768³ volume 16.29 s → 9.70 s (×1.68), peak VRAM 1.2 and 3.3 GiB. The end-to-end speedup is smaller than the filtering speedup and shrinks as the volume grows, because the GPU backprojection is unchanged and takes a growing share of the total — which is the point: the filtering is no longer the bottleneck.

workers

filtering() gets an optional workers argument. By default it uses the number of cores the process is allowed to run on (os.sched_getaffinity where it exists, os.cpu_count() elsewhere), capped at 16: on the machine above, 24 threads and more were slower than 16 (1.5 s with 24, 1.9 s with 32, 2.6 s with 40, against 1.0 s with 16, for 360 × 1024 × 1024). FDK() and fbp() are not touched; passing workers through their kwargs is a small follow-up if you want it exposed there.

One known limit: a CPU quota (a container started with --cpus=4) does not show in the affinity mask, so the default still starts 16 threads there. In that setting the 1024 × 1024 case took 3.6 s with the default against 2.6 s with workers=4.

Tests

I have put the test file in the PR: Python/tests/test_filtering_workers.py, 37 cases — against a plain per-projection float64 filter, same bits for any number of workers and block size, in-place behaviour, an exception in a worker thread reaching the caller, and default_workers. I kept tests out of #781 and #783 as you preferred there; here I would rather hand them over, since they are what backs the "same output" claim, but say the word and I will drop the file if you would rather keep it separate. The tests of Python/tests that run without a GPU give the same result with and without the change.

Not touched

For a detector width that is a power of two, nextpow2(2 * n) returns one more than MATLAB's nextpow2, so filt_len is 4 times the width where filtering.m uses 2 times (2048 against 1024 for 512 columns). Using MATLAB's length would halve the FFT work again, but it changes the output slightly (2.5e-7 relative with ram_lak, hamming and hann, up to 5e-5 with shepp_logan and cosine, on random data), so I left it for a separate decision.

AI Usage

• Type of assistance: agentic code analysis, implementation and benchmarking, run by me.
• Scope: read the FDK/FBP path and the discussions in Implement FDK filtering on GPU #312, optimisation of fdk python setup #331 and FBP filtering on GPU  #423, wrote the change to Python/tigre/utilities/filtering.py, the harness that compares it with master, the benchmark and the test file (not part of the PR), and drafted this description.
• How it was found: by looking for the step that dominates the CPU time of the Python FDK, then keeping only a change whose output could be shown identical to master's.
• Level of modification: the patch, the harness and the test were written by the agent; the comparisons and the timings above were run by the pipeline in a network-less container, on the bare machine, and end to end on an RTX 3070 GPU, and are reported as run.

Found by a pipeline I build and run (Dev-next-gen), using Claude Code with Anthropic's Claude Fable 5.1.

Nothing this pipeline produces is submitted without human approval. The whole pull request — the change, the proof and this description — was read and approved by a human (lu et approuvé) before submission.

(@Dev-next-gen)