Parallel Beam Filtered Backprojection (FBP)

This example demonstrates 2D parallel beam filtered backprojection (FBP) reconstruction using the ParallelProjectorFunction and ParallelBackprojectorFunction from diffct.

Overview

Parallel beam FBP is the standard analytical reconstruction method for 2D parallel beam CT. This example shows how to:

  • Configure parallel beam geometry parameters

  • Generate synthetic projection data using the Shepp-Logan phantom

  • Apply ramp filtering in the frequency domain

  • Perform backprojection to reconstruct the image

  • Visualize and evaluate reconstruction results

Mathematical Background

Parallel Beam Geometry

In parallel beam CT, X-rays are collimated into parallel beams. The projection at angle \(\theta\) and detector position \(t\) is given by the Radon transform:

\[p(t, \theta) = \int_{-\infty}^{\infty} f(t\cos\theta - s\sin\theta, t\sin\theta + s\cos\theta) \, ds\]

where \(f(x,y)\) is the 2D attenuation coefficient distribution.

Filtered Backprojection Algorithm

The FBP reconstruction consists of three sequential steps:

  1. Forward Projection: Compute sinogram using the Radon transform

  2. Ramp Filtering: Apply frequency domain filter \(H(\omega) = |\omega|\)

  3. Backprojection: Reconstruct using filtered projections

The complete FBP formula is:

\[f(x,y) = \int_0^\pi p_f(x\cos\theta + y\sin\theta, \theta) \, d\theta\]

where \(p_f(t, \theta)\) is the filtered projection:

\[p_f(t, \theta) = \mathcal{F}^{-1}\{|\omega| \cdot \mathcal{F}\{p(t, \theta)\}\}\]

Implementation Steps

  1. Phantom Generation: Create Shepp-Logan phantom with 5 ellipses

  2. Forward Projection: Generate sinogram using ParallelProjectorFunction

  3. Ramp Filtering: Apply \(H(\omega) = |\omega|\) filter in frequency domain

  4. Backprojection: Reconstruct using ParallelBackprojectorFunction

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

Shepp-Logan Phantom

The phantom consists of 5 ellipses representing brain tissue structures:

  • Outer skull: Large ellipse with low attenuation

  • Brain tissue: Medium ellipse with baseline attenuation

  • Ventricles: Small ellipses with fluid-like attenuation

  • Lesions: High-contrast features for reconstruction assessment

Each ellipse is defined by center position, semi-axes, rotation angle, and attenuation coefficient.

2D Parallel Beam FBP Example
 1import numpy as np
 2import torch
 3import matplotlib.pyplot as plt
 4import torch.nn.functional as F
 5from diffct.differentiable import (
 6    ParallelProjectorFunction,
 7    ParallelBackprojectorFunction,
 8    angular_integration_weights,
 9    ramp_filter_1d,
10)
11
12
13def shepp_logan_2d(Nx, Ny):
14    Nx = int(Nx)
15    Ny = int(Ny)
16    phantom = np.zeros((Ny, Nx), dtype=np.float64)
17    ellipses = [
18        (0.0, 0.0, 0.69, 0.92, 0, 1.0),
19        (0.0, -0.0184, 0.6624, 0.8740, 0, -0.8),
20        (0.22, 0.0, 0.11, 0.31, -18.0, -0.8),
21        (-0.22, 0.0, 0.16, 0.41, 18.0, -0.8),
22        (0.0, 0.35, 0.21, 0.25, 0, 0.7),
23    ]
24    cx = (Nx - 1) / 2
25    cy = (Ny - 1) / 2
26    for ix in range(Nx):
27        for iy in range(Ny):
28            xnorm = (ix - cx) / (Nx / 2)
29            ynorm = (iy - cy) / (Ny / 2)
30            val = 0.0
31            for (x0, y0, a, b, angdeg, ampl) in ellipses:
32                th = np.deg2rad(angdeg)
33                xprime = (xnorm - x0) * np.cos(th) + (ynorm - y0) * np.sin(th)
34                yprime = -(xnorm - x0) * np.sin(th) + (ynorm - y0) * np.cos(th)
35                if xprime * xprime / (a * a) + yprime * yprime / (b * b) <= 1.0:
36                    val += ampl
37            phantom[iy, ix] = val
38    phantom = np.clip(phantom, 0.0, 1.0)
39    return phantom
40
41def main():
42    Nx, Ny = 256, 256
43    phantom = shepp_logan_2d(Nx, Ny)
44    num_angles = 360
45    angles_np = np.linspace(0, 2*np.pi, num_angles, endpoint=False).astype(np.float32)
46
47    num_detectors = 512
48    detector_spacing = 1.0
49    voxel_spacing = 1.0
50
51    device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
52    image_torch = torch.tensor(phantom, device=device, dtype=torch.float32, requires_grad=True)
53    angles_torch = torch.tensor(angles_np, device=device, dtype=torch.float32)
54
55    sinogram = ParallelProjectorFunction.apply(image_torch, angles_torch,
56                                               num_detectors, detector_spacing, voxel_spacing)
57    
58    sinogram_filt = ramp_filter_1d(sinogram, dim=1)
59    d_theta = angular_integration_weights(angles_torch, redundant_full_scan=True).view(-1, 1)
60    sinogram_filt = sinogram_filt * d_theta
61
62    reconstruction = F.relu(ParallelBackprojectorFunction.apply(sinogram_filt, angles_torch,
63                                                         detector_spacing, Ny, Nx, voxel_spacing)) # ReLU to ensure non-negativity
64
65    loss = torch.mean((reconstruction - image_torch)**2)
66    loss.backward()
67
68    print("Loss:", loss.item())
69    print("Phantom gradient center pixel:", image_torch.grad[Ny//2, Nx//2].item())
70    print("Reconstruction shape:", reconstruction.shape)
71
72    sinogram_cpu = sinogram.detach().cpu().numpy()
73    reco_cpu = reconstruction.detach().cpu().numpy()
74
75    plt.figure(figsize=(12,4))
76    plt.subplot(1,3,1)
77    plt.imshow(phantom, cmap='gray')
78    plt.title("Phantom")
79    plt.axis('off')
80    plt.subplot(1,3,2)
81    plt.imshow(sinogram_cpu, aspect='auto', cmap='gray')
82    plt.title("Differentiable Sinogram")
83    plt.axis('off')
84    plt.subplot(1,3,3)
85    plt.imshow(reco_cpu, cmap='gray')
86    plt.title("Differentiable Recon")
87    plt.axis('off')
88    plt.tight_layout()
89    plt.show()
90
91    # print data range of the phantom and reco 
92    print("Phantom range:", phantom.min(), phantom.max())
93    print("Reco range:", reco_cpu.min(), reco_cpu.max())
94
95if __name__ == "__main__":
96    main()