Fan Beam Filtered Backprojection (FBP)#

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

Fan-beam FBP of the Shepp-Logan phantom and of a patient slice

Ground truth: Shepp-Logan, and a chest slice of TCIA case C002 (see Data sources), lung window. 512², 720 views over 360°, magnification 1.6, Shepp-Logan window. PSNR in the lower corner.#

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 full-turn circular fan-beam scan with 360 views, sid = 2.5 * size, sdd = 4 * size, 3 * size detector cells and a pitch of 0.8. Distances use the same units as the unit-spaced voxels.

  1. Project a phantom slice with Projector(beam="fan").

  2. Apply fan_cosine_weights for the detector pitch and source-detector distance.

  3. Apply ramp_filter_1d along dim=1, with sample_spacing=pitch and pad_factor=2.

  4. Multiply by angular_integration_weights with redundant_full_scan=True.

  5. Call fan_weighted_backproject for the distance-weighted gather.

This example uses a full scan. Short-scan reconstruction needs Parker weights and suitable angular coverage. Changing the view count does not enable it.

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()