Fan Beam Iterative Reconstruction

This example demonstrates 2D fan beam iterative reconstruction using the differentiable FanProjectorFunction and FanBackprojectorFunction from diffct.

Overview

Fan beam iterative reconstruction extends the optimization approach to the more realistic fan beam geometry. This example shows how to:

  • Formulate fan beam CT reconstruction as an optimization problem

  • Handle geometric complexities of divergent ray geometry

  • Apply gradient-based optimization with fan beam operators

  • Monitor convergence and reconstruction quality

Mathematical Background

Fan Beam Iterative Formulation

The fan beam reconstruction problem is formulated as:

\[\hat{f} = \arg\min_f \|A_{\text{fan}}(f) - p\|_2^2 + \lambda R(f)\]

where \(A_{\text{fan}}\) is the fan beam forward projection operator accounting for divergent ray geometry.

Fan Beam Forward Model

The fan beam projection operator maps 2D image \(f(x,y)\) to sinogram \(p(\beta, u)\):

\[p(\beta, u) = \int_{\text{ray}} f(x,y) \, dl\]

where integration follows the ray from point source to detector element \(u\) at source angle \(\beta\).

Gradient Computation

The gradient involves the fan beam backprojection operator (adjoint):

\[\frac{\partial L}{\partial f} = 2A_{\text{fan}}^T(A_{\text{fan}}(f) - p_{\text{measured}})\]

where \(A_{\text{fan}}^T\) is the fan beam backprojection operator.

Geometric Considerations

Fan beam geometry introduces complexities compared to parallel beam:

  • Ray Divergence: Non-parallel rays affect sampling density and conditioning

  • Magnification Effects: Variable magnification across the field of view

  • Non-uniform Resolution: Spatial resolution varies with distance from rotation center

  • Geometric Distortion: Requires careful handling of coordinate transformations

Implementation Steps

  1. Geometry Setup: Configure fan beam parameters (SID, SDD)

  2. Problem Formulation: Define parameterized image and fan beam forward model

  3. Loss Computation: Calculate L2 distance using FanProjectorFunction

  4. Gradient Computation: Use automatic differentiation through fan beam operators

  5. Optimization: Apply Adam optimizer with appropriate learning rate

  6. Convergence Monitoring: Track reconstruction quality and loss evolution

Model Architecture

The fan beam reconstruction model consists of:

  • Parameterized Image: Learnable 2D tensor representing the unknown image

  • Fan Beam Forward Model: FanProjectorFunction with geometric parameters

  • Loss Function: Mean squared error between predicted and measured sinograms

Convergence Characteristics

Fan beam reconstruction typically exhibits:

  1. Initial Convergence (0-100 iterations): Rapid loss decrease, basic structure

  2. Detail Refinement (100-500 iterations): Fine features develop, slower progress

  3. Final Convergence (500+ iterations): Minimal improvement, convergence plateau

Challenges and Solutions

  • Conditioning: Fan beam system matrix may have different conditioning properties

  • Geometric Artifacts: Proper weighting and filtering help reduce artifacts

  • Parameter Tuning: Learning rate may need adjustment for optimal convergence

  • Memory Usage: Similar to parallel beam but with additional geometric computations

2D Fan Beam Iterative Example
  1import math
  2import torch
  3import numpy as np
  4import matplotlib.pyplot as plt
  5import torch.nn as nn
  6import torch.optim as optim
  7from diffct.differentiable import FanProjectorFunction
  8
  9
 10def shepp_logan_2d(Nx, Ny):
 11    Nx = int(Nx)
 12    Ny = int(Ny)
 13    phantom = np.zeros((Ny, Nx), dtype=np.float32)
 14    ellipses = [
 15        (0.0, 0.0, 0.69, 0.92, 0, 1.0),
 16        (0.0, -0.0184, 0.6624, 0.8740, 0, -0.8),
 17        (0.22, 0.0, 0.11, 0.31, -18.0, -0.8),
 18        (-0.22, 0.0, 0.16, 0.41, 18.0, -0.8),
 19        (0.0, 0.35, 0.21, 0.25, 0, 0.7),
 20    ]
 21    cx = (Nx - 1)*0.5
 22    cy = (Ny - 1)*0.5
 23    for ix in range(Nx):
 24        for iy in range(Ny):
 25            xnorm = (ix - cx)/(Nx/2)
 26            ynorm = (iy - cy)/(Ny/2)
 27            val = 0.0
 28            for (x0, y0, a, b, angdeg, ampl) in ellipses:
 29                th = np.deg2rad(angdeg)
 30                xprime = (xnorm - x0)*np.cos(th) + (ynorm - y0)*np.sin(th)
 31                yprime = -(xnorm - x0)*np.sin(th) + (ynorm - y0)*np.cos(th)
 32                if xprime*xprime/(a*a) + yprime*yprime/(b*b) <= 1.0:
 33                    val += ampl
 34            phantom[iy, ix] = val
 35    phantom = np.clip(phantom, 0.0, 1.0)
 36    return phantom
 37
 38class IterativeRecoModel(nn.Module):
 39    def __init__(self, volume_shape, angles, 
 40                 num_detectors, detector_spacing, 
 41                 sdd, sid, voxel_spacing):
 42        
 43        super().__init__()
 44        self.reco = nn.Parameter(torch.zeros(volume_shape))
 45        self.angles = angles
 46        self.num_detectors = num_detectors
 47        self.detector_spacing = detector_spacing
 48        self.sdd = sdd
 49        self.sid = sid
 50        self.relu = nn.ReLU() # non negative constraint
 51        self.voxel_spacing = voxel_spacing
 52
 53    def forward(self, x):
 54        updated_reco = x + self.reco
 55        current_sino = FanProjectorFunction.apply(updated_reco, self.angles, 
 56                                                  self.num_detectors, self.detector_spacing, 
 57                                                  self.sdd, self.sid, self.voxel_spacing)
 58        return current_sino, self.relu(updated_reco)
 59
 60class Pipeline:
 61    def __init__(self, lr, volume_shape, angles, 
 62                 num_detectors, detector_spacing, 
 63                 sdd, sid, voxel_spacing,
 64                 device, epoches=1000):
 65        
 66        self.epoches = epoches
 67        self.model = IterativeRecoModel(volume_shape, angles,
 68                                        num_detectors, detector_spacing, 
 69                                        sdd, sid, voxel_spacing).to(device)
 70        
 71        self.optimizer = optim.AdamW(list(self.model.parameters()), lr=lr)
 72        self.loss = nn.MSELoss()
 73
 74    def train(self, input, label):
 75        loss_values = []
 76        for epoch in range(self.epoches):
 77            self.optimizer.zero_grad()
 78            predictions, current_reco = self.model(input)
 79            loss_value = self.loss(predictions, label)
 80            loss_value.backward()
 81            self.optimizer.step()
 82            loss_values.append(loss_value.item())
 83
 84            if epoch % 10 == 0:
 85                print(f"Epoch {epoch}, Loss: {loss_value.item()}")
 86                
 87        return loss_values, self.model
 88
 89def main():
 90    Nx, Ny = 128, 128
 91    phantom_cpu = shepp_logan_2d(Nx, Ny)
 92
 93    num_views = 360
 94    angles_np = np.linspace(0, 2 * math.pi, num_views, endpoint=False).astype(np.float32)
 95
 96    num_detectors = 256
 97    detector_spacing = 1.0
 98    voxel_spacing = 1.0
 99    sdd = 600.0
100    sid = 400.0
101
102    device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
103    phantom_torch = torch.tensor(phantom_cpu, device=device, dtype=torch.float32)
104    angles_torch = torch.tensor(angles_np, device=device, dtype=torch.float32)
105
106    # Generate the "real" sinogram
107    real_sinogram = FanProjectorFunction.apply(phantom_torch, angles_torch,
108                                               num_detectors, detector_spacing,
109                                               sdd, sid, voxel_spacing)
110
111    pipeline_instance = Pipeline(lr=1e-1,
112                                 volume_shape=(Ny,Nx),
113                                 angles=angles_torch,
114                                 num_detectors=num_detectors,
115                                 detector_spacing=detector_spacing,
116                                 sdd=sdd, voxel_spacing=voxel_spacing,
117                                 sid=sid,
118                                 device=device, epoches=1000)
119
120    ini_guess = torch.zeros_like(phantom_torch)
121
122    loss_values, trained_model = pipeline_instance.train(ini_guess, real_sinogram)
123
124    reco = trained_model(ini_guess)[1].squeeze().cpu().detach().numpy()
125
126    plt.figure()
127    plt.plot(loss_values)
128    plt.title("Loss Curve")
129    plt.xlabel("Epoch")
130    plt.ylabel("Loss")
131    plt.show()
132
133    plt.figure(figsize=(12, 6))
134    plt.subplot(1, 2, 1)
135    plt.imshow(phantom_cpu, cmap="gray")
136    plt.title("Original Phantom")
137    plt.axis("off")
138
139    plt.subplot(1, 2, 2)
140    plt.imshow(reco, cmap="gray")
141    plt.title("Reconstructed")
142    plt.axis("off")
143    plt.show()
144
145if __name__ == "__main__":
146    main()