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 cell-constant Siddon ray-march (each ray accumulates each traversed pixel weighted by its exact intersection length, no sub-pixel 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    if not torch.cuda.is_available():
139        raise RuntimeError("This example requires CUDA.")
140    device = torch.device("cuda")
141    image_torch = torch.tensor(phantom, device=device, dtype=torch.float32)
142    angles_torch = torch.tensor(angles_np, device=device, dtype=torch.float32)
143
144    # ------------------------------------------------------------------
145    # 4.5  Pick a forward projector backend
146    # ------------------------------------------------------------------
147    # ``FanProjectorFunction`` / ``FanBackprojectorFunction`` both accept a
148    # ``backend`` keyword that selects the underlying CUDA kernel family.
149    # The choice applies to both forward and adjoint, and each option ships
150    # with a matched scatter/gather pair so autograd and the standalone
151    # Backprojector Function work unchanged.
152    #
153    #   "siddon" (default) - Ray-driven cell-constant Siddon traversal:
154    #                        each traversed pixel contributes its value
155    #                        weighted by the exact chord length, no sub-
156    #                        pixel interpolation. One thread per (view,
157    #                        detector bin). Fastest forward. Pick this
158    #                        when you only need a forward projection and
159    #                        not a matched cell-integrated model.
160    #
161    #   "sf"               - Voxel-driven separable-footprint projector
162    #                        (SF-TR of Long-Fessler-Balter, IEEE TMI 2010).
163    #                        Each voxel's projection footprint is a
164    #                        trapezoid built from the four projected
165    #                        corners and closed-form integrated over each
166    #                        detector cell, so **mass is conserved per
167    #                        voxel**. About 3x slower forward than
168    #                        "siddon". On analytical FBP reconstructions
169    #                        with the matched SF gather backprojector
170    #                        (``fan_weighted_backproject(backend="sf")``),
171    #                        SF and VD produce visually identical edge
172    #                        profiles on Shepp-Logan at typical CBCT
173    #                        magnifications - the "SF is sharper"
174    #                        advantage that shows up in the SF / LEAP
175    #                        literature only manifests at extreme
176    #                        sub-nominal voxel sizes that are not hit
177    #                        in standard examples. The real reason to
178    #                        pick "sf" is the **forward** side: if you
179    #                        plan to use this projector inside an
180    #                        iterative reco, a learned prior, or any
181    #                        loss that compares sinograms directly,
182    #                        a cell-integrated mass-conserving forward
183    #                        is the right model.
184    #
185    # The default is kept at "sf" so the reader can see the SF path run
186    # end-to-end; switching it to "siddon" gives a visually equivalent
187    # reconstruction at this geometry.
188    projector_backend = "sf"
189
190    # ------------------------------------------------------------------
191    # 5. Forward projection: image -> fan sinogram
192    # ------------------------------------------------------------------
193    # ``FanProjectorFunction`` is the differentiable fan-beam forward
194    # projector. It returns a (num_angles, num_detectors) sinogram. The
195    # call is autograd-aware so the same function can be used inside an
196    # iterative reconstruction loop (see ``iterative_reco_fan.py``).
197    # ``backend`` selects the CUDA kernel family used for both the forward
198    # and its adjoint - see step 4.5 above for the trade-offs.
199    sinogram = FanProjectorFunction.apply(
200        image_torch,
201        angles_torch,
202        num_detectors,
203        detector_spacing,
204        sdd,
205        sid,
206        voxel_spacing,
207        detector_offset,
208        0.0,                # center_offset_x
209        0.0,                # center_offset_y
210        projector_backend,
211    )
212
213    # ==================================================================
214    # 6. FBP analytical reconstruction
215    # ==================================================================
216
217    # --- 6.1  Optional Parker redundancy weighting -------------------
218    # For short-scan trajectories, ``parker_weights`` returns a
219    # ``(num_angles, num_detectors)`` weight that smoothly tapers views
220    # near the two ends of the angular range so each ray contributes
221    # exactly once. For a full 2*pi scan this helper returns all-ones
222    # and is a no-op, so the ``if`` is really only for short scans.
223    if apply_parker:
224        parker = parker_weights(
225            angles_torch, num_detectors, detector_spacing, sdd, detector_offset
226        )
227        sinogram = sinogram * parker
228
229    # --- 6.2  Fan-beam cosine pre-weighting --------------------------
230    # Multiplies each detector cell by ``cos(gamma) = sdd / sqrt(sdd^2
231    # + u^2)``, i.e. the cosine of the fan angle. This compensates for
232    # the extra path length that off-center rays traverse relative to
233    # the principal ray.
234    weights = fan_cosine_weights(
235        num_detectors,
236        detector_spacing,
237        sdd,
238        detector_offset=detector_offset,
239        device=device,
240        dtype=image_torch.dtype,
241    ).unsqueeze(0)
242    sino_weighted = sinogram * weights
243
244    # --- 6.3  1D ramp filter along the detector axis -----------------
245    # See the fdk_cone.py example for the full list of options - the
246    # same ramp filter is used here. For fan FBP the recommended call
247    # is ``sample_spacing=detector_spacing``, ``pad_factor=2``,
248    # ``window="hann"``.
249    sinogram_filt = ramp_filter_1d(
250        sino_weighted,
251        dim=1,
252        sample_spacing=detector_spacing,
253        pad_factor=2,
254        window="hann",
255    ).contiguous()
256
257    # --- 6.4  Per-view angular integration weights -------------------
258    # For a full 2*pi scan uniformly sampled, each view contributes
259    # ``pi / num_angles`` to the FBP integral (the ``1/2`` redundancy
260    # factor of the full-scan formula is absorbed inside
261    # ``redundant_full_scan=True``). For a short scan with Parker
262    # weights we already handle redundancy there, so pass
263    # ``redundant_full_scan=False``.
264    d_beta = angular_integration_weights(
265        angles_torch, redundant_full_scan=(not apply_parker)
266    ).view(-1, 1)
267    sinogram_filt = sinogram_filt * d_beta
268
269    # --- 6.5  Voxel-driven FBP backprojection ------------------------
270    # ``fan_weighted_backproject`` dispatches to one of two fan FBP
271    # gather kernels based on ``backend``:
272    #
273    #   "siddon" (default) - bilinear voxel-driven gather: linearly
274    #                        sample the filtered sinogram at each
275    #                        pixel's projected u-coordinate, multiply
276    #                        by ``(sid/U)^2`` and accumulate.
277    #   "sf"               - LEAP-style chord-weighted separable-
278    #                        footprint gather: integrate the filtered
279    #                        sinogram over each pixel's trapezoidal
280    #                        footprint and weight by the in-plane
281    #                        chord through the voxel plus the fan
282    #                        ``sid/U`` first-power weight (matches the
283    #                        matched-adjoint form in LEAP's
284    #                        ``projectors_SF.cu``). Amplitude-calibrated
285    #                        against the Siddon path on Shepp-Logan;
286    #                        edge profiles at typical CBCT geometries
287    #                        are visually indistinguishable from the
288    #                        Siddon path. Pick this when you want a
289    #                        cell-integrated forward / backward model
290    #                        matched to the SF forward projector.
291    #
292    # Here we pass the same ``projector_backend`` we picked at step
293    # 4.5 so forward and backward stay consistent. Amplitude is
294    # calibrated by the wrapper so the returned image is ready to
295    # compare against ``image_torch`` directly.
296    reconstruction_raw = fan_weighted_backproject(
297        sinogram_filt,
298        angles_torch,
299        detector_spacing,
300        Ny,
301        Nx,
302        sdd,
303        sid,
304        voxel_spacing=voxel_spacing,
305        detector_offset=detector_offset,
306        backend=projector_backend,
307    )
308    reconstruction = F.relu(reconstruction_raw)
309
310    # ------------------------------------------------------------------
311    # 7. Quantitative summary
312    # ------------------------------------------------------------------
313    raw_loss = torch.mean((reconstruction_raw - image_torch) ** 2)
314    clamped_loss = torch.mean((reconstruction - image_torch) ** 2)
315
316    scan_label = "Parker short scan" if apply_parker else "full 2*pi scan"
317    print(f"Fan Beam FBP example ({scan_label}):")
318    print(f"  Raw MSE:              {raw_loss.item():.6f}")
319    print(f"  Clamped MSE:          {clamped_loss.item():.6f}")
320    print(f"  Reconstruction shape: {tuple(reconstruction.shape)}")
321    print(
322        "  Raw reco data range:  "
323        f"[{reconstruction_raw.min().item():.4f}, {reconstruction_raw.max().item():.4f}]"
324    )
325    print(
326        "  Clamped reco range:   "
327        f"[{reconstruction.min().item():.4f}, {reconstruction.max().item():.4f}]"
328    )
329    print(
330        "  Phantom data range:   "
331        f"[{float(phantom.min()):.4f}, {float(phantom.max()):.4f}]"
332    )
333
334    # ------------------------------------------------------------------
335    # 8. Visualization
336    # ------------------------------------------------------------------
337    sinogram_cpu = sinogram.detach().cpu().numpy()
338    reco_cpu = reconstruction.detach().cpu().numpy()
339
340    plt.figure(figsize=(12, 4))
341    plt.subplot(1, 3, 1)
342    plt.imshow(phantom, cmap="gray")
343    plt.title("Phantom")
344    plt.axis("off")
345    plt.subplot(1, 3, 2)
346    plt.imshow(sinogram_cpu, cmap="gray", aspect="auto")
347    plt.title("Fan Sinogram")
348    plt.axis("off")
349    plt.subplot(1, 3, 3)
350    plt.imshow(reco_cpu, cmap="gray")
351    plt.title("Fan Reconstruction")
352    plt.axis("off")
353    plt.tight_layout()
354    plt.show()
355
356
357if __name__ == "__main__":
358    main()