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 ParallelProjectorFunction, ParallelBackprojectorFunction
  6
  7
  8def shepp_logan_2d(Nx, Ny):
  9    Nx = int(Nx)
 10    Ny = int(Ny)
 11    phantom = np.zeros((Ny, Nx), dtype=np.float64)
 12    ellipses = [
 13        (0.0, 0.0, 0.69, 0.92, 0, 1.0),
 14        (0.0, -0.0184, 0.6624, 0.8740, 0, -0.8),
 15        (0.22, 0.0, 0.11, 0.31, -18.0, -0.8),
 16        (-0.22, 0.0, 0.16, 0.41, 18.0, -0.8),
 17        (0.0, 0.35, 0.21, 0.25, 0, 0.7),
 18    ]
 19    cx = (Nx - 1) / 2
 20    cy = (Ny - 1) / 2
 21    for ix in range(Nx):
 22        for iy in range(Ny):
 23            xnorm = (ix - cx) / (Nx / 2)
 24            ynorm = (iy - cy) / (Ny / 2)
 25            val = 0.0
 26            for (x0, y0, a, b, angdeg, ampl) in ellipses:
 27                th = np.deg2rad(angdeg)
 28                xprime = (xnorm - x0) * np.cos(th) + (ynorm - y0) * np.sin(th)
 29                yprime = -(xnorm - x0) * np.sin(th) + (ynorm - y0) * np.cos(th)
 30                if xprime * xprime / (a * a) + yprime * yprime / (b * b) <= 1.0:
 31                    val += ampl
 32            phantom[iy, ix] = val
 33    phantom = np.clip(phantom, 0.0, 1.0)
 34    return phantom
 35
 36def ramp_filter(sinogram_tensor):
 37    device = sinogram_tensor.device
 38    num_views, num_det = sinogram_tensor.shape
 39    freqs = torch.fft.fftfreq(num_det, device=device)
 40    omega = 2.0 * torch.pi * freqs
 41    ramp = torch.abs(omega)
 42    ramp_2d = ramp.reshape(1, num_det)
 43    sino_fft = torch.fft.fft(sinogram_tensor, dim=1)
 44    filtered_fft = sino_fft * ramp_2d
 45    filtered = torch.real(torch.fft.ifft(filtered_fft, dim=1))
 46    
 47    return filtered
 48
 49def main():
 50    Nx, Ny = 256, 256
 51    phantom = shepp_logan_2d(Nx, Ny)
 52    num_angles = 360
 53    angles_np = np.linspace(0, 2*np.pi, num_angles, endpoint=False).astype(np.float32)
 54
 55    num_detectors = 512
 56    detector_spacing = 1.0
 57    voxel_spacing = 1.0
 58
 59    device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
 60    image_torch = torch.tensor(phantom, device=device, dtype=torch.float32, requires_grad=True)
 61    angles_torch = torch.tensor(angles_np, device=device, dtype=torch.float32)
 62
 63    sinogram = ParallelProjectorFunction.apply(image_torch, angles_torch,
 64                                               num_detectors, detector_spacing, voxel_spacing)
 65    
 66    sinogram_filt = ramp_filter(sinogram)
 67
 68    reconstruction = F.relu(ParallelBackprojectorFunction.apply(sinogram_filt, angles_torch,
 69                                                         detector_spacing, Ny, Nx, voxel_spacing)) # ReLU to ensure non-negativity
 70    
 71    # --- FBP normalization ---
 72    # The backprojection is a sum over all angles. To approximate the integral,
 73    # we need to multiply by the angular step d_theta.
 74    # The FBP formula also includes a factor of 1/2 when integrating over [0, 2*pi].
 75    # d_theta = 2 * pi / num_angles
 76    # Normalization factor = (1/2) * d_theta = pi / num_angles
 77    reconstruction = reconstruction * (np.pi / num_angles)
 78
 79    loss = torch.mean((reconstruction - image_torch)**2)
 80    loss.backward()
 81
 82    print("Loss:", loss.item())
 83    print("Phantom gradient center pixel:", image_torch.grad[Ny//2, Nx//2].item())
 84    print("Reconstruction shape:", reconstruction.shape)
 85
 86    sinogram_cpu = sinogram.detach().cpu().numpy()
 87    reco_cpu = reconstruction.detach().cpu().numpy()
 88
 89    plt.figure(figsize=(12,4))
 90    plt.subplot(1,3,1)
 91    plt.imshow(phantom, cmap='gray')
 92    plt.title("Phantom")
 93    plt.axis('off')
 94    plt.subplot(1,3,2)
 95    plt.imshow(sinogram_cpu, aspect='auto', cmap='gray')
 96    plt.title("Differentiable Sinogram")
 97    plt.axis('off')
 98    plt.subplot(1,3,3)
 99    plt.imshow(reco_cpu, cmap='gray')
100    plt.title("Differentiable Recon")
101    plt.axis('off')
102    plt.tight_layout()
103    plt.show()
104
105    # print data range of the phantom and reco 
106    print("Phantom range:", phantom.min(), phantom.max())
107    print("Reco range:", reco_cpu.min(), reco_cpu.max())
108
109if __name__ == "__main__":
110    main()