Fan Beam Filtered Backprojection (FBP)

This example demonstrates 2D fan beam filtered backprojection (FBP) reconstruction using the FanProjectorFunction and FanBackprojectorFunction from diffct.

Overview

Fan beam FBP extends parallel beam reconstruction to the more realistic fan beam geometry used in clinical CT scanners. This example shows how to:

  • Configure fan beam geometry parameters (SID, SDD)

  • Generate fan beam projections with proper weighting

  • Apply ramp filtering with cosine weighting correction

  • Perform fan beam backprojection reconstruction

Mathematical Background

Fan Beam Geometry

Fan beam CT uses a point X-ray source creating a fan-shaped beam. Key geometric parameters:

  • SDD \(D_s\): Source-to-Detector Distance (distance from X-ray source to detector array)

  • SID \(D_{sid}\): Source-to-Isocenter Distance (distance from X-ray source to rotation center)

  • Fan angle \(\gamma\): Angle between central ray and detector element

The detector position \(u\) relates to fan angle \(\gamma\) by:

\[\gamma = \arctan\left(\frac{u}{D_s}\right)\]

Fan Beam 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 + t \cdot \vec{d}(\beta, u)\right) dt\]

where \(\vec{r}_s\) is the source position and \(\vec{d}(\beta, u)\) is the ray direction vector.

Fan Beam FBP Algorithm

Fan beam FBP reconstruction involves three sequential steps:

  1. Cosine Weighting: Compensate for ray divergence:

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

  2. Ramp Filtering: Apply frequency domain filter:

    \[p_f(\beta, u) = \mathcal{F}^{-1}\{|\omega| \cdot \mathcal{F}\{p_w(\beta, u)\}\}\]

  3. Fan Beam Backprojection: Reconstruct using weighted backprojection:

    \[f(x,y) = \int_0^{2\pi} \frac{D_s^2}{(D_s + x\cos\beta + y\sin\beta)^2} p_f(\beta, u_{xy}) d\beta\]

    where the detector coordinate \(u_{xy}\) for pixel \((x,y)\) is:

    \[u_{xy} = D_s \frac{-x\sin\beta + y\cos\beta}{D_s + x\cos\beta + y\sin\beta}\]

Implementation Steps

  1. Phantom Generation: Create Shepp-Logan phantom for testing

  2. Fan Beam Projection: Generate sinogram using FanProjectorFunction

  3. Cosine Weighting: Apply divergence correction weights

  4. Ramp Filtering: Filter each projection in frequency domain

  5. Fan Beam Backprojection: Reconstruct using FanBackprojectorFunction

  6. Normalization: Scale by \(\frac{\pi}{N_{\text{angles}}}\) factor

Advantages of Fan Beam Geometry

  • Clinical Relevance: Matches real CT scanner geometry

  • Higher Flux: Better X-ray utilization than parallel beam

  • Natural Magnification: Improved spatial resolution

  • Faster Acquisition: Wider coverage per projection angle

2D Fan Beam FBP Example
  1import math
  2import numpy as np
  3import torch
  4import matplotlib.pyplot as plt
  5import torch.nn.functional as F
  6from diffct.differentiable import (
  7    FanProjectorFunction,
  8    angular_integration_weights,
  9    fan_cosine_weights,
 10    fan_weighted_backproject,
 11    parker_weights,
 12    ramp_filter_1d,
 13)
 14
 15
 16def shepp_logan_2d(Nx, Ny):
 17    Nx = int(Nx)
 18    Ny = int(Ny)
 19    phantom = np.zeros((Ny, Nx), dtype=np.float32)
 20    ellipses = [
 21        (0.0, 0.0, 0.69, 0.92, 0, 1.0),
 22        (0.0, -0.0184, 0.6624, 0.8740, 0, -0.8),
 23        (0.22, 0.0, 0.11, 0.31, -18.0, -0.8),
 24        (-0.22, 0.0, 0.16, 0.41, 18.0, -0.8),
 25        (0.0, 0.35, 0.21, 0.25, 0, 0.7),
 26    ]
 27    cx = (Nx - 1)*0.5
 28    cy = (Ny - 1)*0.5
 29    for ix in range(Nx):
 30        for iy in range(Ny):
 31            xnorm = (ix - cx)/(Nx/2)
 32            ynorm = (iy - cy)/(Ny/2)
 33            val = 0.0
 34            for (x0, y0, a, b, angdeg, ampl) in ellipses:
 35                th = np.deg2rad(angdeg)
 36                xprime = (xnorm - x0)*np.cos(th) + (ynorm - y0)*np.sin(th)
 37                yprime = -(xnorm - x0)*np.sin(th) + (ynorm - y0)*np.cos(th)
 38                if xprime*xprime/(a*a) + yprime*yprime/(b*b) <= 1.0:
 39                    val += ampl
 40            phantom[iy, ix] = val
 41    phantom = np.clip(phantom, 0.0, 1.0)
 42    return phantom
 43
 44def main():
 45    Nx, Ny = 256, 256
 46    phantom = shepp_logan_2d(Nx, Ny)
 47    num_angles = 360
 48    angles_np = np.linspace(0, 2*math.pi, num_angles, endpoint=False).astype(np.float32)
 49
 50    num_detectors = 600
 51    detector_spacing = 1.0
 52    detector_offset = 0.0
 53    voxel_spacing = 1.0
 54    sdd = 800.0
 55    sid = 500.0
 56    apply_parker = False
 57
 58    device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
 59    image_torch = torch.tensor(phantom, device=device, dtype=torch.float32)
 60    angles_torch = torch.tensor(angles_np, device=device, dtype=torch.float32)
 61
 62    sinogram = FanProjectorFunction.apply(image_torch, angles_torch, num_detectors,
 63                                          detector_spacing, sdd, sid, voxel_spacing)
 64
 65    # --- FBP weighting and filtering ---
 66    # 1) Optional Parker redundancy weighting for short-scan trajectories
 67    if apply_parker:
 68        parker = parker_weights(angles_torch, num_detectors, detector_spacing, sdd, detector_offset)
 69        sinogram = sinogram * parker
 70
 71    # 2) Fan-beam cosine pre-weighting
 72    weights = fan_cosine_weights(
 73        num_detectors,
 74        detector_spacing,
 75        sdd,
 76        detector_offset=detector_offset,
 77        device=device,
 78        dtype=image_torch.dtype,
 79    ).unsqueeze(0)
 80    sino_weighted = sinogram * weights
 81
 82    # 3) Ramp filter along detector axis
 83    sinogram_filt = ramp_filter_1d(sino_weighted, dim=1)
 84
 85    # 4) Angle-integration weights
 86    d_beta = angular_integration_weights(angles_torch, redundant_full_scan=(not apply_parker)).view(-1, 1)
 87    sinogram_filt = sinogram_filt * d_beta
 88
 89    # 5) Weighted fan-beam backprojection
 90    reconstruction = F.relu(
 91        fan_weighted_backproject(
 92            sinogram_filt,
 93            angles_torch,
 94            detector_spacing,
 95            Ny,
 96            Nx,
 97            sdd,
 98            sid,
 99            voxel_spacing=voxel_spacing,
100            detector_offset=detector_offset,
101        )
102    )
103
104    loss = torch.mean((reconstruction - image_torch)**2)
105
106    print("Loss:", loss.item())
107
108    sinogram_cpu = sinogram.detach().cpu().numpy()
109    reco_cpu = reconstruction.detach().cpu().numpy()
110
111    plt.figure(figsize=(12,4))
112    plt.subplot(1,3,1)
113    plt.imshow(phantom, cmap='gray')
114    plt.title("Phantom")
115    plt.axis('off')
116    plt.subplot(1,3,2)
117    plt.imshow(sinogram_cpu, cmap='gray', aspect='auto')
118    plt.title("Fan Sinogram")
119    plt.axis('off')
120    plt.subplot(1,3,3)
121    plt.imshow(reco_cpu, cmap='gray')
122    plt.title("Fan Reconstruction")
123    plt.axis('off')
124    plt.tight_layout()
125    plt.show()
126
127    # print data range of the phantom and reco
128    print("Phantom range:", phantom.min(), phantom.max())
129    print("Reco range:", reco_cpu.min(), reco_cpu.max())
130
131if __name__ == "__main__":
132    main()