Fan Beam Filtered Backprojection (FBP)

This example demonstrates 2D fan-beam filtered backprojection on a Shepp-Logan phantom using the public call chain from diffct:

FanProjectorFunction.apply      -- differentiable fan forward projection
parker_weights (optional)       -- short-scan redundancy weighting
fan_cosine_weights              -- cos(gamma) pre-weighting
ramp_filter_1d                  -- ramp filter along the detector axis
angular_integration_weights     -- per-view integration weights
fan_weighted_backproject        -- voxel-driven FBP backprojection gather

The accompanying source file is examples/fbp_fan.py.

Overview

Fan-beam FBP is the analytical reconstruction method for 2D fan-beam CT on a circular source orbit. Compared with parallel beam it adds a divergent-ray geometry (source-to-isocenter SID and source-to-detector SDD distances), a cosine pre-weight, and a voxel-dependent distance weight inside the backprojection. This example shows how to:

  • Configure the 2D fan-beam geometry (source, detector, orbit).

  • Generate fan-beam projections of a 2D Shepp-Logan phantom.

  • Apply FBP cosine pre-weighting, ramp filtering (with optional zero-padding and frequency-domain windowing), and angular-integration weights.

  • Run the dedicated voxel-driven FBP gather kernel through fan_weighted_backproject to reconstruct the image, already amplitude-calibrated.

Mathematical Background

Fan-beam geometry

  • \(D\) (sid): Source-to-Isocenter Distance. Source position at angle \(\beta\) is \(\vec{r}_s(\beta) = (-D \sin\beta, D \cos\beta)\).

  • \(D_{sd}\) (sdd): Source-to-Detector Distance. The detector is at distance sdd from the source perpendicular to the central ray.

  • \(u\): detector in-plane coordinate; \(du\) is the detector cell pitch.

  • \(\gamma = \arctan(u / D_{sd})\): fan angle of a detector cell.

Forward projection

The fan-beam projection at source angle \(\beta\) and detector position \(u\) is

\[p(\beta, u) = \int_0^{\infty} f\!\left(\vec{r}_s(\beta) + t \cdot \vec{d}(\beta, u)\right) \, dt\]

where \(\vec{d}\) is the unit ray direction from source to the detector cell. FanProjectorFunction implements this integral via a Siddon ray-march with bilinear interpolation.

FBP algorithm

  1. Cosine pre-weighting of the raw projection:

    \[p_w(\beta, u) = \frac{D_{sd}}{\sqrt{D_{sd}^2 + u^2}}\, p(\beta, u) = p(\beta, u) \cos\gamma.\]

  2. Ramp filtering along the detector axis:

    \[p_f(\beta, u) = \mathcal{F}_u^{-1}\bigl\{\,|\omega_u|\, \mathcal{F}_u\{p_w(\beta, u)\}\bigr\}.\]

    ramp_filter_1d lets you choose between Ram-Lak, Hann, Hamming, Cosine and Shepp-Logan shapes.

  3. Voxel-driven backprojection with the classical FBP distance weight:

    \[f(x, y) = \frac{1}{2} \int_0^{2\pi} \frac{D^2}{U(\beta, x, y)^2}\, p_f\!\left(\beta, u_{xy}\right) d\beta,\]

    where \(U(\beta, x, y) = D + x \sin\beta - y \cos\beta\) is the distance from the source to the pixel along the central ray direction, and the pixel’s projected detector coordinate is

    \[u_{xy} = \frac{D_{sd}}{U}\,(x \cos\beta + y \sin\beta).\]

    fan_weighted_backproject dispatches to a dedicated voxel-driven gather kernel _fan_2d_fbp_backproject_kernel for this step. For each pixel it computes the projected \(u\), linearly samples the filtered sinogram, multiplies by \((D/U)^2\), and accumulates across views. The Siddon-based autograd fan backprojector is reserved for the differentiable path used by iterative reconstruction and is not touched by this analytical pipeline.

Ramp Filter Options

ramp_filter_1d accepts several high-precision parameters (all optional and fully backward compatible):

  • sample_spacing – physical detector pitch along dim. Set this to detector_spacing for fan FBP so the filter output is in physical units and the reconstructed amplitude is detector-pitch independent.

  • pad_factor – zero-pad the detector dimension to pad_factor * N before the FFT. 2 is recommended for FBP.

  • window – frequency-domain apodization. Options: None / "ram-lak" (unwindowed), "hann" (used in this example), "hamming", "cosine", "shepp-logan".

  • use_rfft – use torch.fft.rfft for real inputs (default True, ~2x faster than the complex path).

Analytical FBP scale

fan_weighted_backproject multiplies its output by the analytical constant sdd / (2 * pi * sid). The 1/(2*pi) comes from the Fourier convention used by ramp_filter_1d (the |omega| ramp in radian frequency), and the sdd/sid factor corrects for the fact that the ramp is applied on the physical detector plane rather than the virtual isocenter-plane detector. The analytical cone FDK helper uses the same constant for the same reason.

Implementation Steps

The main() function in examples/fbp_fan.py follows the same 8-step structure as the cone FDK example: volume geometry, source trajectory, detector geometry, CUDA transfer, forward projection, FBP pipeline (optional Parker / cosine pre-weight / ramp filter / angular weights / voxel-gather backprojection / optional ReLU), quantitative summary, visualization.

Short-scan support (Parker weighting)

For a short-scan trajectory (angular coverage of \(\pi + 2\gamma_{\max}\)), set apply_parker=True in the example. The parker_weights helper builds the standard smooth redundancy weight that tapers the two ends of the angular range. When Parker is enabled the redundancy factor of the full-scan FBP formula is dropped, so the example passes redundant_full_scan=False to angular_integration_weights.

Source

2D Fan Beam FBP Example
  1"""Circular-orbit fan-beam FBP reconstruction example.
  2
  3Pipeline (matches the cone-beam FDK example, one dimension lower):
  4
  5    FanProjectorFunction.apply      # forward projection (sinogram)
  6    parker_weights (optional)       # short-scan redundancy weighting
  7    fan_cosine_weights              # cos(gamma) pre-weighting
  8    ramp_filter_1d                  # ramp filter along detector axis
  9    angular_integration_weights     # per-view integration weights
 10    fan_weighted_backproject        # voxel-driven FBP gather
 11
 12Every geometry and reconstruction parameter below is annotated with
 13its meaning, units, and available options so the file can be used as
 14a reference for the other analytical 2D entry points.
 15"""
 16
 17import math
 18
 19import matplotlib.pyplot as plt
 20import numpy as np
 21import torch
 22import torch.nn.functional as F
 23
 24from diffct.differentiable import (
 25    FanProjectorFunction,
 26    angular_integration_weights,
 27    fan_cosine_weights,
 28    fan_weighted_backproject,
 29    parker_weights,
 30    ramp_filter_1d,
 31)
 32
 33
 34def shepp_logan_2d(Nx, Ny):
 35    """2D Shepp-Logan phantom clipped to ``[0, 1]``."""
 36    Nx = int(Nx)
 37    Ny = int(Ny)
 38    phantom = np.zeros((Ny, Nx), dtype=np.float32)
 39    # (x0, y0, a, b, angle_deg, amplitude)
 40    ellipses = [
 41        (0.0,    0.0,    0.69,   0.92,    0.0,  1.0),
 42        (0.0,   -0.0184, 0.6624, 0.8740,  0.0, -0.8),
 43        (0.22,   0.0,    0.11,   0.31,  -18.0, -0.8),
 44        (-0.22,  0.0,    0.16,   0.41,   18.0, -0.8),
 45        (0.0,    0.35,   0.21,   0.25,    0.0,  0.7),
 46    ]
 47    cx = (Nx - 1) * 0.5
 48    cy = (Ny - 1) * 0.5
 49    for ix in range(Nx):
 50        for iy in range(Ny):
 51            xnorm = (ix - cx) / (Nx / 2)
 52            ynorm = (iy - cy) / (Ny / 2)
 53            val = 0.0
 54            for (x0, y0, a, b, angdeg, ampl) in ellipses:
 55                th = np.deg2rad(angdeg)
 56                xprime = (xnorm - x0) * np.cos(th) + (ynorm - y0) * np.sin(th)
 57                yprime = -(xnorm - x0) * np.sin(th) + (ynorm - y0) * np.cos(th)
 58                if xprime * xprime / (a * a) + yprime * yprime / (b * b) <= 1.0:
 59                    val += ampl
 60            phantom[iy, ix] = val
 61    return np.clip(phantom, 0.0, 1.0)
 62
 63
 64def main():
 65    # ------------------------------------------------------------------
 66    # 1. Image geometry
 67    # ------------------------------------------------------------------
 68    # ``Nx`` / ``Ny`` are the reconstruction grid size in pixels. The
 69    # phantom tensor has shape ``(Ny, Nx)`` (rows, cols) which matches
 70    # the ``(H, W)`` layout every 2D routine in diffct expects.
 71    Nx, Ny = 256, 256
 72    phantom = shepp_logan_2d(Nx, Ny)
 73
 74    # ``voxel_spacing`` is the physical size of one pixel in the same
 75    # length unit used by ``detector_spacing``, ``sdd`` and ``sid``
 76    # (commonly millimeters). Internally, all physical quantities are
 77    # divided by ``voxel_spacing``, so only their *ratios* matter.
 78    voxel_spacing = 1.0
 79
 80    # ------------------------------------------------------------------
 81    # 2. Detector geometry
 82    # ------------------------------------------------------------------
 83    # (Listed before the trajectory so the short-scan coverage below can
 84    # use the detector fan angle to compute ``pi + 2*gamma_max``.)
 85    #
 86    # ``num_detectors`` is the number of detector cells along the
 87    # detector axis. ``detector_spacing`` is their physical pitch.
 88    # Make sure the detector is wide enough that no ray that intersects
 89    # the reconstructed field of view ever projects outside it - rays
 90    # that miss the detector are zero-filled and introduce truncation
 91    # artifacts at the image edges.
 92    num_detectors = 600
 93    detector_spacing = 1.0
 94
 95    # Principal-ray offset (shifts the whole detector sideways).
 96    # Non-zero values model a detector that is not perfectly centered
 97    # on the source-isocenter line.
 98    detector_offset = 0.0
 99
100    # ``sdd`` = source-to-detector distance, ``sid`` = source-to-
101    # isocenter distance, both in physical units. The magnification at
102    # the detector is ``sdd / sid`` (here = 1.6). Typical clinical fan
103    # beam systems have magnifications around 1.3 - 2.0.
104    sdd = 800.0
105    sid = 500.0
106
107    # ------------------------------------------------------------------
108    # 3. Source trajectory (circular orbit)
109    # ------------------------------------------------------------------
110    # ``apply_parker`` selects between two supported trajectories:
111    #
112    #   False -> full 2*pi scan, 1/2 FBP redundancy factor absorbed
113    #            inside ``redundant_full_scan=True``.
114    #
115    #   True  -> minimal short scan ``pi + 2*gamma_max`` with Parker
116    #            redundancy weighting. ``gamma_max`` is computed from
117    #            the detector half-width and ``sdd``.
118    #
119    # Both branches reuse the same ``num_angles`` so reconstruction
120    # runtime stays identical - only the angular *range* changes.
121    apply_parker = False
122
123    if apply_parker:
124        u_max = ((num_detectors - 1) * 0.5) * detector_spacing + abs(detector_offset)
125        gamma_max = math.atan(u_max / sdd)
126        scan_range = math.pi + 2.0 * gamma_max
127    else:
128        scan_range = 2.0 * math.pi
129
130    num_angles = 360
131    angles_np = np.linspace(
132        0.0, scan_range, num_angles, endpoint=False
133    ).astype(np.float32)
134
135    # ------------------------------------------------------------------
136    # 4. Move everything to CUDA
137    # ------------------------------------------------------------------
138    device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
139    image_torch = torch.tensor(phantom, device=device, dtype=torch.float32)
140    angles_torch = torch.tensor(angles_np, device=device, dtype=torch.float32)
141
142    # ------------------------------------------------------------------
143    # 4.5  Pick a forward projector backend
144    # ------------------------------------------------------------------
145    # ``FanProjectorFunction`` / ``FanBackprojectorFunction`` both accept a
146    # ``backend`` keyword that selects the underlying CUDA kernel family.
147    # The choice applies to both forward and adjoint, and each option ships
148    # with a matched scatter/gather pair so autograd and the standalone
149    # Backprojector Function work unchanged.
150    #
151    #   "siddon" (default) - Ray-driven Siddon traversal with bilinear
152    #                        interpolation on the image. One thread per
153    #                        (view, detector bin). Fastest forward. Pick
154    #                        this when you only need a forward projection
155    #                        and not a matched cell-integrated model.
156    #
157    #   "sf"               - Voxel-driven separable-footprint projector
158    #                        (SF-TR of Long-Fessler-Balter, IEEE TMI 2010).
159    #                        Each voxel's projection footprint is a
160    #                        trapezoid built from the four projected
161    #                        corners and closed-form integrated over each
162    #                        detector cell, so **mass is conserved per
163    #                        voxel**. About 3x slower forward than
164    #                        "siddon". On analytical FBP reconstructions
165    #                        with the matched SF gather backprojector
166    #                        (``fan_weighted_backproject(backend="sf")``),
167    #                        SF and VD produce visually identical edge
168    #                        profiles on Shepp-Logan at typical CBCT
169    #                        magnifications - the "SF is sharper"
170    #                        advantage that shows up in the SF / LEAP
171    #                        literature only manifests at extreme
172    #                        sub-nominal voxel sizes that are not hit
173    #                        in standard examples. The real reason to
174    #                        pick "sf" is the **forward** side: if you
175    #                        plan to use this projector inside an
176    #                        iterative reco, a learned prior, or any
177    #                        loss that compares sinograms directly,
178    #                        a cell-integrated mass-conserving forward
179    #                        is the right model.
180    #
181    # The default is kept at "sf" so the reader can see the SF path run
182    # end-to-end; switching it to "siddon" gives a visually equivalent
183    # reconstruction at this geometry.
184    projector_backend = "sf"
185
186    # ------------------------------------------------------------------
187    # 5. Forward projection: image -> fan sinogram
188    # ------------------------------------------------------------------
189    # ``FanProjectorFunction`` is the differentiable fan-beam forward
190    # projector. It returns a (num_angles, num_detectors) sinogram. The
191    # call is autograd-aware so the same function can be used inside an
192    # iterative reconstruction loop (see ``iterative_reco_fan.py``).
193    # ``backend`` selects the CUDA kernel family used for both the forward
194    # and its adjoint - see step 4.5 above for the trade-offs.
195    sinogram = FanProjectorFunction.apply(
196        image_torch,
197        angles_torch,
198        num_detectors,
199        detector_spacing,
200        sdd,
201        sid,
202        voxel_spacing,
203        detector_offset,
204        0.0,                # center_offset_x
205        0.0,                # center_offset_y
206        projector_backend,
207    )
208
209    # ==================================================================
210    # 6. FBP analytical reconstruction
211    # ==================================================================
212
213    # --- 6.1  Optional Parker redundancy weighting -------------------
214    # For short-scan trajectories, ``parker_weights`` returns a
215    # ``(num_angles, num_detectors)`` weight that smoothly tapers views
216    # near the two ends of the angular range so each ray contributes
217    # exactly once. For a full 2*pi scan this helper returns all-ones
218    # and is a no-op, so the ``if`` is really only for short scans.
219    if apply_parker:
220        parker = parker_weights(
221            angles_torch, num_detectors, detector_spacing, sdd, detector_offset
222        )
223        sinogram = sinogram * parker
224
225    # --- 6.2  Fan-beam cosine pre-weighting --------------------------
226    # Multiplies each detector cell by ``cos(gamma) = sdd / sqrt(sdd^2
227    # + u^2)``, i.e. the cosine of the fan angle. This compensates for
228    # the extra path length that off-center rays traverse relative to
229    # the principal ray.
230    weights = fan_cosine_weights(
231        num_detectors,
232        detector_spacing,
233        sdd,
234        detector_offset=detector_offset,
235        device=device,
236        dtype=image_torch.dtype,
237    ).unsqueeze(0)
238    sino_weighted = sinogram * weights
239
240    # --- 6.3  1D ramp filter along the detector axis -----------------
241    # See the fdk_cone.py example for the full list of options - the
242    # same ramp filter is used here. For fan FBP the recommended call
243    # is ``sample_spacing=detector_spacing``, ``pad_factor=2``,
244    # ``window="hann"``.
245    sinogram_filt = ramp_filter_1d(
246        sino_weighted,
247        dim=1,
248        sample_spacing=detector_spacing,
249        pad_factor=2,
250        window="hann",
251    ).contiguous()
252
253    # --- 6.4  Per-view angular integration weights -------------------
254    # For a full 2*pi scan uniformly sampled, each view contributes
255    # ``pi / num_angles`` to the FBP integral (the ``1/2`` redundancy
256    # factor of the full-scan formula is absorbed inside
257    # ``redundant_full_scan=True``). For a short scan with Parker
258    # weights we already handle redundancy there, so pass
259    # ``redundant_full_scan=False``.
260    d_beta = angular_integration_weights(
261        angles_torch, redundant_full_scan=(not apply_parker)
262    ).view(-1, 1)
263    sinogram_filt = sinogram_filt * d_beta
264
265    # --- 6.5  Voxel-driven FBP backprojection ------------------------
266    # ``fan_weighted_backproject`` dispatches to one of two fan FBP
267    # gather kernels based on ``backend``:
268    #
269    #   "siddon" (default) - bilinear voxel-driven gather: linearly
270    #                        sample the filtered sinogram at each
271    #                        pixel's projected u-coordinate, multiply
272    #                        by ``(sid/U)^2`` and accumulate.
273    #   "sf"               - LEAP-style chord-weighted separable-
274    #                        footprint gather: integrate the filtered
275    #                        sinogram over each pixel's trapezoidal
276    #                        footprint and weight by the in-plane
277    #                        chord through the voxel plus the fan
278    #                        ``sid/U`` first-power weight (matches the
279    #                        matched-adjoint form in LEAP's
280    #                        ``projectors_SF.cu``). Amplitude-calibrated
281    #                        against the Siddon path on Shepp-Logan;
282    #                        edge profiles at typical CBCT geometries
283    #                        are visually indistinguishable from the
284    #                        Siddon path. Pick this when you want a
285    #                        cell-integrated forward / backward model
286    #                        matched to the SF forward projector.
287    #
288    # Here we pass the same ``projector_backend`` we picked at step
289    # 4.5 so forward and backward stay consistent. Amplitude is
290    # calibrated by the wrapper so the returned image is ready to
291    # compare against ``image_torch`` directly.
292    reconstruction_raw = fan_weighted_backproject(
293        sinogram_filt,
294        angles_torch,
295        detector_spacing,
296        Ny,
297        Nx,
298        sdd,
299        sid,
300        voxel_spacing=voxel_spacing,
301        detector_offset=detector_offset,
302        backend=projector_backend,
303    )
304    reconstruction = F.relu(reconstruction_raw)
305
306    # ------------------------------------------------------------------
307    # 7. Quantitative summary
308    # ------------------------------------------------------------------
309    raw_loss = torch.mean((reconstruction_raw - image_torch) ** 2)
310    clamped_loss = torch.mean((reconstruction - image_torch) ** 2)
311
312    scan_label = "Parker short scan" if apply_parker else "full 2*pi scan"
313    print(f"Fan Beam FBP example ({scan_label}):")
314    print(f"  Raw MSE:              {raw_loss.item():.6f}")
315    print(f"  Clamped MSE:          {clamped_loss.item():.6f}")
316    print(f"  Reconstruction shape: {tuple(reconstruction.shape)}")
317    print(
318        "  Raw reco data range:  "
319        f"[{reconstruction_raw.min().item():.4f}, {reconstruction_raw.max().item():.4f}]"
320    )
321    print(
322        "  Clamped reco range:   "
323        f"[{reconstruction.min().item():.4f}, {reconstruction.max().item():.4f}]"
324    )
325    print(
326        "  Phantom data range:   "
327        f"[{float(phantom.min()):.4f}, {float(phantom.max()):.4f}]"
328    )
329
330    # ------------------------------------------------------------------
331    # 8. Visualization
332    # ------------------------------------------------------------------
333    sinogram_cpu = sinogram.detach().cpu().numpy()
334    reco_cpu = reconstruction.detach().cpu().numpy()
335
336    plt.figure(figsize=(12, 4))
337    plt.subplot(1, 3, 1)
338    plt.imshow(phantom, cmap="gray")
339    plt.title("Phantom")
340    plt.axis("off")
341    plt.subplot(1, 3, 2)
342    plt.imshow(sinogram_cpu, cmap="gray", aspect="auto")
343    plt.title("Fan Sinogram")
344    plt.axis("off")
345    plt.subplot(1, 3, 3)
346    plt.imshow(reco_cpu, cmap="gray")
347    plt.title("Fan Reconstruction")
348    plt.axis("off")
349    plt.tight_layout()
350    plt.show()
351
352
353if __name__ == "__main__":
354    main()