Cone Beam Feldkamp-Davis-Kress (FDK)#

This page describes the cone-beam part of examples/analytical_reconstruction.py. For filter options and helper signatures, see API Reference.

Cone-beam projection and FDK slices of a 3D phantom

3D Shepp-Logan phantom, 256³, 720 views on a circular orbit. The off-centre slice and the coronal section show the approximation error of FDK away from the orbit plane.#

Run it#

Run the commands from the root of a source checkout; pip install diffct does not install the examples/ directory. See the “From source” tab in Installation. --figure needs Matplotlib: pip install matplotlib.

python examples/analytical_reconstruction.py --size 64 --window shepp-logan
python examples/analytical_reconstruction.py --size 128 --figure analytical.png

The script runs parallel, fan and cone reconstruction on one CUDA GPU and prints PSNR against the synthetic phantom. --figure requires Matplotlib.

Pipeline#

The example uses a flat detector and a circular cone-beam scan with 360 views. It does not reconstruct native curved-detector data directly. To reuse this pipeline for a curved detector, first resample line integrals along matching rays onto a covered virtual flat grid. This adds interpolation error; FDK still retains its circular cone-beam approximation. For native curved projection and matched backprojection, see Parameterized detector surfaces and Examples.

The example projects a 3D phantom with Projector into a (views, U, V) sinogram. The helper examples/_common.py::fdk then calls diffct.reconstruction.fdk(), which performs these steps:

  1. Apply detector cosine weights computed from the actual per-view rays.

  2. Filter along the detector u axis (dim=1) with sample_spacing=du and pad_factor=2. The last axis is v, not u.

  3. Apply angular_integration_weights. A full turn includes the 1/2 redundancy factor.

  4. Call cone_weighted_backproject. The example clamps negative values to zero.

The same stages are public helpers. See Differentiable analytical reconstruction to replace or learn the weights and filter.

FDK is approximate for a circular cone-beam scan away from the central plane. Large cone angles increase artifacts. On non-circular trajectories, FDK is a heuristic baseline, not an exact reconstruction.

The analytical helpers are not the Projector.backproject() adjoint. Do not use the adjoint in this pipeline, and do not add extra scale factors.

Source#

Parallel FBP, fan FBP and cone FDK#
 1"""Analytical reconstruction: parallel-beam FBP, fan-beam FBP and cone-beam FDK.
 2
 3Each pipeline is the textbook one: optional cosine weighting, ramp filter with
 4a window, angular weights, voxel-driven weighted backprojection. The window
 5trades sharpness against noise: "ram-lak" keeps all frequencies, "shepp-logan"
 6and "cosine" damp the highest, "hamming" and "hann" damp the most.
 7
 8    python examples/analytical_reconstruction.py --window shepp-logan --figure fbp.png
 9"""
10
11import argparse
12import math
13from pathlib import Path
14
15import torch
16
17from diffct import (
18    Projector,
19    angular_integration_weights,
20    circular_trajectory_2d_fan,
21    circular_trajectory_2d_parallel,
22    fan_cosine_weights,
23    fan_weighted_backproject,
24    fbp,
25    parallel_weighted_backproject,
26    ramp_filter_1d,
27)
28from _common import Scan, fdk, psnr, save_slices, shepp_logan_3d
29
30
31def main():
32    parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
33    parser.add_argument("--size", type=int, default=128)
34    parser.add_argument("--window", default="shepp-logan",
35                        choices=["ram-lak", "shepp-logan", "cosine", "hamming", "hann"])
36    parser.add_argument("--figure", type=Path)
37    args = parser.parse_args()
38    if not torch.cuda.is_available():
39        raise RuntimeError("diffct needs a CUDA GPU")
40    n = args.size
41    phantom = torch.from_numpy(shepp_logan_3d(n)).cuda()
42    image = phantom[n // 2]
43    views = 360
44    full_turn = torch.linspace(0.0, 2 * math.pi, views + 1, device="cuda")[:-1]
45    half_turn = torch.linspace(0.0, math.pi, views + 1, device="cuda")[:-1]
46
47    # Parallel beam over 180 degrees; detector pitch 0.5 voxel.
48    n_det, pitch = 3 * n, 0.5
49    trajectory = circular_trajectory_2d_parallel(views, end_angle=math.pi, device="cpu")
50    sinogram = Projector(trajectory, (n, n), n_det, beam="parallel", detector_spacing=pitch).project(image)
51    filtered = ramp_filter_1d(sinogram, dim=1, sample_spacing=pitch, pad_factor=2, window=args.window)
52    filtered = filtered * angular_integration_weights(half_turn, redundant_full_scan=False).view(-1, 1)
53    parallel = parallel_weighted_backproject(filtered.contiguous(), *[t.cuda() for t in trajectory],
54                                             pitch, n, n).clamp_min(0)
55    parallel_api = fbp(sinogram, trajectory, (n, n), half_turn, detector_spacing=pitch,
56                       window=args.window).clamp_min(0)
57    parallel_difference = (parallel - parallel_api).abs().max().item()
58
59    # Fan beam over a full turn; magnification 1.6, detector pitch 0.8 (0.5 voxel at the isocentre).
60    sid, sdd, n_det, pitch = 2.5 * n, 4.0 * n, 3 * n, 0.8
61    trajectory = circular_trajectory_2d_fan(views, sid, sdd, device="cpu")
62    sinogram = Projector(trajectory, (n, n), n_det, beam="fan", detector_spacing=pitch).project(image)
63    weighted = sinogram * fan_cosine_weights(n_det, pitch, sdd, device="cuda").unsqueeze(0)
64    filtered = ramp_filter_1d(weighted, dim=1, sample_spacing=pitch, pad_factor=2, window=args.window)
65    filtered = filtered * angular_integration_weights(full_turn, redundant_full_scan=True).view(-1, 1)
66    fan = fan_weighted_backproject(filtered.contiguous(), *[t.cuda() for t in trajectory],
67                                   pitch, n, n).clamp_min(0)
68    fan_api = fbp(sinogram, trajectory, (n, n), full_turn, beam="fan", detector_spacing=pitch,
69                  window=args.window, isocenter=(0, 0)).clamp_min(0)
70    fan_difference = (fan - fan_api).abs().max().item()
71
72    # Cone beam FDK on a circular orbit.
73    scan = Scan(n, views)
74    trajectory = scan.trajectory("circular")
75    sinogram = Projector(trajectory, (n, n, n), scan.detector, detector_spacing=scan.pitch).project(phantom)
76    cone = fdk(scan, sinogram, trajectory, window=args.window)
77
78    print(f"window {args.window}: parallel FBP {psnr(parallel, image):.2f} dB, "
79          f"fan FBP {psnr(fan, image):.2f} dB, cone FDK {psnr(cone, phantom):.2f} dB")
80    print(f"manual/API maximum difference: parallel {parallel_difference:.3g}, fan {fan_difference:.3g}")
81    if args.figure:
82        save_slices(args.figure, [("phantom\nground truth", phantom),
83                                  (f"parallel FBP\n{psnr(parallel, image):.1f} dB", parallel),
84                                  (f"fan FBP\n{psnr(fan, image):.1f} dB", fan),
85                                  (f"cone FDK\n{psnr(cone, phantom):.1f} dB", cone)],
86                    f"{n}² / {n}³ Shepp-Logan, {args.window} window")
87
88
89if __name__ == "__main__":
90    main()