ct_reconstruction.tex

\documentclass[a4paper, 11pt]{report}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{amsmath, amssymb}
\usepackage{graphicx}
\usepackage{siunitx}
\usepackage{booktabs}
\usepackage{xcolor}
\usepackage[makestderr]{pythontex}

\definecolor{phantom}{RGB}{52, 152, 219}
\definecolor{sinogram}{RGB}{231, 76, 60}
\definecolor{recon}{RGB}{46, 204, 113}

\title{Computed Tomography:\\
Image Reconstruction and Artifact Analysis}
\author{Department of Medical Physics\\Technical Report MP-2024-001}
\date{\today}

\begin{document}
\maketitle

\begin{abstract}
This report presents a comprehensive analysis of CT image reconstruction algorithms. We implement the Radon transform, filtered back-projection with various filters, analyze reconstruction artifacts, compare iterative methods, and demonstrate noise reduction techniques. All simulations use PythonTeX for reproducibility.
\end{abstract}

\tableofcontents

\chapter{Introduction}

Computed Tomography reconstructs cross-sectional images from X-ray projections using the Radon transform:
\begin{equation}
p(s, \theta) = \int_{-\infty}^{\infty}\int_{-\infty}^{\infty} f(x, y) \delta(x\cos\theta + y\sin\theta - s) \, dx \, dy
\end{equation}

\section{Central Slice Theorem}
The Fourier transform of a projection equals a slice through the 2D Fourier transform:
\begin{equation}
P(\omega, \theta) = F(\omega\cos\theta, \omega\sin\theta)
\end{equation}

\begin{pycode}
import numpy as np
from scipy import ndimage
from scipy.fft import fft, ifft, fftfreq
import matplotlib.pyplot as plt
plt.rc('text', usetex=True)
plt.rc('font', family='serif')

np.random.seed(42)

def create_shepp_logan(size):
    """Create simplified Shepp-Logan phantom."""
    phantom = np.zeros((size, size))
    y, x = np.ogrid[-size//2:size//2, -size//2:size//2]

    # Outer ellipse (skull)
    mask = (x/0.69/size*2)**2 + (y/0.92/size*2)**2 <= 1
    phantom[mask] = 2.0

    # Brain
    mask = (x/0.6624/size*2)**2 + (y/0.874/size*2)**2 <= 1
    phantom[mask] = -0.98

    # Ventricles
    mask = ((x+0.22*size/2)/0.11/size*2)**2 + (y/0.31/size*2)**2 <= 1
    phantom[mask] = -0.02
    mask = ((x-0.22*size/2)/0.16/size*2)**2 + (y/0.41/size*2)**2 <= 1
    phantom[mask] = -0.02

    # Small features
    mask = ((x-0.35*size/2)/0.21/size*2)**2 + (y/0.25/size*2)**2 <= 1
    phantom[mask] = 0.01
    mask = (x/0.046/size*2)**2 + (y/0.046/size*2)**2 <= 1
    phantom[mask] = 0.01
    mask = ((x+0.08*size/2)/0.046/size*2)**2 + (y/0.023/size*2)**2 <= 1
    phantom[mask] = 0.01

    return phantom + 1

def radon_transform(image, angles):
    """Compute Radon transform (sinogram)."""
    n_angles = len(angles)
    size = image.shape[0]
    sinogram = np.zeros((size, n_angles))

    for i, angle in enumerate(angles):
        rotated = ndimage.rotate(image, -angle, reshape=False, order=3)
        sinogram[:, i] = np.sum(rotated, axis=0)

    return sinogram

def create_filter(size, filter_type='ram-lak'):
    """Create reconstruction filter in frequency domain."""
    freq = fftfreq(size)

    if filter_type == 'ram-lak':
        H = np.abs(freq)
    elif filter_type == 'shepp-logan':
        H = np.abs(freq) * np.sinc(freq)
    elif filter_type == 'cosine':
        H = np.abs(freq) * np.cos(np.pi * freq)
    elif filter_type == 'hamming':
        H = np.abs(freq) * (0.54 + 0.46 * np.cos(2 * np.pi * freq))
    else:
        H = np.ones(size)

    return H

def filtered_backprojection(sinogram, angles, filter_type='ram-lak'):
    """Reconstruct image using filtered back-projection."""
    size = sinogram.shape[0]
    n_angles = len(angles)
    reconstruction = np.zeros((size, size))

    # Create filter
    H = create_filter(size, filter_type)

    # Filter projections
    filtered_sinogram = np.zeros_like(sinogram)
    for i in range(n_angles):
        proj_ft = fft(sinogram[:, i])
        filtered_sinogram[:, i] = np.real(ifft(proj_ft * H))

    # Back-projection
    for i, angle in enumerate(angles):
        back_proj = np.tile(filtered_sinogram[:, i], (size, 1)).T
        back_proj = ndimage.rotate(back_proj, angle, reshape=False, order=1)
        reconstruction += back_proj

    return reconstruction * np.pi / n_angles
\end{pycode}

\chapter{Radon Transform and Sinogram}

\begin{pycode}
# Create phantom
size = 128
phantom = create_shepp_logan(size)

# Different angular sampling
angles_180 = np.linspace(0, 180, 180, endpoint=False)
angles_90 = np.linspace(0, 180, 90, endpoint=False)
angles_45 = np.linspace(0, 180, 45, endpoint=False)

sinogram_180 = radon_transform(phantom, angles_180)
sinogram_90 = radon_transform(phantom, angles_90)
sinogram_45 = radon_transform(phantom, angles_45)

fig, axes = plt.subplots(2, 3, figsize=(14, 8))

# Original phantom
ax = axes[0, 0]
im = ax.imshow(phantom, cmap='gray')
ax.set_title('Shepp-Logan Phantom')
ax.set_xlabel('x (pixels)')
ax.set_ylabel('y (pixels)')
plt.colorbar(im, ax=ax, fraction=0.046)

# Sinogram
ax = axes[0, 1]
im = ax.imshow(sinogram_180.T, cmap='gray', aspect='auto',
               extent=[0, size, 180, 0])
ax.set_title('Sinogram (180 projections)')
ax.set_xlabel('Detector Position')
ax.set_ylabel('Angle (degrees)')
plt.colorbar(im, ax=ax, fraction=0.046)

# Single projection
ax = axes[0, 2]
for angle_idx, label in [(0, '0'), (45, '45'), (90, '90')]:
    ax.plot(sinogram_180[:, angle_idx], label=f'{label}$^\\circ$')
ax.set_xlabel('Detector Position')
ax.set_ylabel('Projection Value')
ax.set_title('Sample Projections')
ax.legend()
ax.grid(True, alpha=0.3)

# Reconstructions with different angles
recon_180 = filtered_backprojection(sinogram_180, angles_180)
recon_90 = filtered_backprojection(sinogram_90, angles_90)
recon_45 = filtered_backprojection(sinogram_45, angles_45)

titles = ['45 projections', '90 projections', '180 projections']
recons = [recon_45, recon_90, recon_180]

for ax, recon, title in zip(axes[1], recons, titles):
    im = ax.imshow(recon, cmap='gray')
    ax.set_title(f'Reconstruction ({title})')
    ax.set_xlabel('x (pixels)')
    ax.set_ylabel('y (pixels)')
    plt.colorbar(im, ax=ax, fraction=0.046)

plt.tight_layout()
plt.savefig('ct_sinogram.pdf', dpi=150, bbox_inches='tight')
plt.close()

# Calculate metrics
mse_180 = np.mean((phantom - recon_180/recon_180.max()*phantom.max())**2)
mse_90 = np.mean((phantom - recon_90/recon_90.max()*phantom.max())**2)
mse_45 = np.mean((phantom - recon_45/recon_45.max()*phantom.max())**2)
\end{pycode}

\begin{figure}[htbp]
\centering
\includegraphics[width=0.95\textwidth]{ct_sinogram.pdf}
\caption{CT reconstruction: (a) phantom, (b) sinogram, (c) projections, (d-f) reconstructions with different angular sampling.}
\end{figure}

\chapter{Reconstruction Filters}

\section{Filter Comparison}
The ramp filter (Ram-Lak) is modified to reduce high-frequency noise:
\begin{equation}
H_{RL}(\omega) = |\omega|, \quad H_{SL}(\omega) = |\omega|\text{sinc}(\omega)
\end{equation}

\begin{pycode}
# Compare different filters
fig, axes = plt.subplots(2, 3, figsize=(14, 8))

# Filter frequency responses
ax = axes[0, 0]
size_filter = 256
freq = fftfreq(size_filter)
freq_plot = np.fft.fftshift(freq)

for filter_type in ['ram-lak', 'shepp-logan', 'cosine', 'hamming']:
    H = create_filter(size_filter, filter_type)
    H_plot = np.fft.fftshift(H)
    ax.plot(freq_plot, H_plot, label=filter_type.replace('-', ' ').title())

ax.set_xlabel('Frequency')
ax.set_ylabel('Filter Response')
ax.set_title('Reconstruction Filters')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_xlim(-0.5, 0.5)

# Reconstruct with different filters
filter_types = ['ram-lak', 'shepp-logan', 'cosine', 'hamming']
recon_filters = {}
mse_filters = {}

for filter_type in filter_types:
    recon = filtered_backprojection(sinogram_180, angles_180, filter_type)
    recon_filters[filter_type] = recon
    mse = np.mean((phantom - recon/recon.max()*phantom.max())**2)
    mse_filters[filter_type] = mse

# Display reconstructions
for ax, filter_type in zip([axes[0, 1], axes[0, 2], axes[1, 0], axes[1, 1]],
                           filter_types):
    im = ax.imshow(recon_filters[filter_type], cmap='gray')
    ax.set_title(f'{filter_type.replace("-", " ").title()} Filter')
    ax.set_xlabel('x (pixels)')
    ax.set_ylabel('y (pixels)')
    plt.colorbar(im, ax=ax, fraction=0.046)

# MSE comparison
ax = axes[1, 2]
filter_names = [f.replace('-', ' ').title() for f in filter_types]
mse_values = [mse_filters[f] for f in filter_types]
bars = ax.bar(filter_names, mse_values, color=['#3498db', '#e74c3c', '#2ecc71', '#f39c12'])
ax.set_xlabel('Filter Type')
ax.set_ylabel('Mean Squared Error')
ax.set_title('Reconstruction Quality')
ax.grid(True, alpha=0.3, axis='y')
plt.setp(ax.xaxis.get_majorticklabels(), rotation=45, ha='right')

plt.tight_layout()
plt.savefig('ct_filters.pdf', dpi=150, bbox_inches='tight')
plt.close()
\end{pycode}

\begin{figure}[htbp]
\centering
\includegraphics[width=0.95\textwidth]{ct_filters.pdf}
\caption{Filter comparison: (a) frequency responses, (b-e) reconstructions, (f) MSE comparison.}
\end{figure}

\chapter{Reconstruction Artifacts}

\begin{pycode}
# Demonstrate common artifacts
fig, axes = plt.subplots(2, 3, figsize=(14, 8))

# Metal artifact
phantom_metal = phantom.copy()
y, x = np.ogrid[-size//2:size//2, -size//2:size//2]
metal_mask = ((x-20)**2 + (y+15)**2 <= 25)
phantom_metal[metal_mask] = 10  # High attenuation

sinogram_metal = radon_transform(phantom_metal, angles_180)
recon_metal = filtered_backprojection(sinogram_metal, angles_180)

ax = axes[0, 0]
im = ax.imshow(recon_metal, cmap='gray')
ax.set_title('Metal Artifact')
ax.set_xlabel('x (pixels)')
ax.set_ylabel('y (pixels)')
plt.colorbar(im, ax=ax, fraction=0.046)

# Limited angle (missing wedge)
angles_limited = np.linspace(0, 120, 120, endpoint=False)
sinogram_limited = radon_transform(phantom, angles_limited)
recon_limited = filtered_backprojection(sinogram_limited, angles_limited)

ax = axes[0, 1]
im = ax.imshow(recon_limited, cmap='gray')
ax.set_title('Limited Angle (0-120$^\\circ$)')
ax.set_xlabel('x (pixels)')
ax.set_ylabel('y (pixels)')
plt.colorbar(im, ax=ax, fraction=0.046)

# Sparse angle (streak artifacts)
angles_sparse = np.linspace(0, 180, 18, endpoint=False)
sinogram_sparse = radon_transform(phantom, angles_sparse)
recon_sparse = filtered_backprojection(sinogram_sparse, angles_sparse)

ax = axes[0, 2]
im = ax.imshow(recon_sparse, cmap='gray')
ax.set_title('Sparse Angle (18 views)')
ax.set_xlabel('x (pixels)')
ax.set_ylabel('y (pixels)')
plt.colorbar(im, ax=ax, fraction=0.046)

# Noisy sinogram
noise_level = 0.1
sinogram_noisy = sinogram_180 + noise_level * np.random.randn(*sinogram_180.shape) * sinogram_180.max()
recon_noisy = filtered_backprojection(sinogram_noisy, angles_180)

ax = axes[1, 0]
im = ax.imshow(recon_noisy, cmap='gray')
ax.set_title(f'Noisy Data ({int(noise_level*100)}\\% noise)')
ax.set_xlabel('x (pixels)')
ax.set_ylabel('y (pixels)')
plt.colorbar(im, ax=ax, fraction=0.046)

# Beam hardening simulation
sinogram_hardened = sinogram_180 * (1 - 0.1 * sinogram_180/sinogram_180.max())
recon_hardened = filtered_backprojection(sinogram_hardened, angles_180)

ax = axes[1, 1]
im = ax.imshow(recon_hardened, cmap='gray')
ax.set_title('Beam Hardening')
ax.set_xlabel('x (pixels)')
ax.set_ylabel('y (pixels)')
plt.colorbar(im, ax=ax, fraction=0.046)

# Profile comparison
ax = axes[1, 2]
center = size // 2
ax.plot(phantom[center, :], 'k-', label='Original', linewidth=2)
ax.plot(recon_180[center, :]/recon_180.max()*phantom.max(), 'b--',
        label='Clean', alpha=0.7)
ax.plot(recon_noisy[center, :]/recon_noisy.max()*phantom.max(), 'r:',
        label='Noisy', alpha=0.7)
ax.set_xlabel('Position (pixels)')
ax.set_ylabel('Intensity')
ax.set_title('Central Profile')
ax.legend()
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('ct_artifacts.pdf', dpi=150, bbox_inches='tight')
plt.close()
\end{pycode}

\begin{figure}[htbp]
\centering
\includegraphics[width=0.95\textwidth]{ct_artifacts.pdf}
\caption{CT artifacts: (a) metal artifact, (b) limited angle, (c) sparse angle, (d) noisy data, (e) beam hardening, (f) profile comparison.}
\end{figure}

\chapter{Numerical Results}

\begin{pycode}
results_table = [
    ('Image size', f'{size} $\\times$ {size}', 'pixels'),
    ('MSE (180 projections)', f'{mse_180:.6f}', ''),
    ('MSE (90 projections)', f'{mse_90:.6f}', ''),
    ('MSE (45 projections)', f'{mse_45:.6f}', ''),
    ('Best filter (lowest MSE)', min(mse_filters, key=mse_filters.get).replace('-', ' ').title(), ''),
    ('Ram-Lak MSE', f'{mse_filters["ram-lak"]:.6f}', ''),
]
\end{pycode}

\begin{table}[htbp]
\centering
\caption{CT reconstruction results}
\begin{tabular}{@{}lcc@{}}
\toprule
Parameter & Value & Units \\
\midrule
\py{results_table[0][0]} & \py{results_table[0][1]} & \py{results_table[0][2]} \\
\py{results_table[1][0]} & \py{results_table[1][1]} & \py{results_table[1][2]} \\
\py{results_table[2][0]} & \py{results_table[2][1]} & \py{results_table[2][2]} \\
\py{results_table[3][0]} & \py{results_table[3][1]} & \py{results_table[3][2]} \\
\py{results_table[4][0]} & \py{results_table[4][1]} & \py{results_table[4][2]} \\
\py{results_table[5][0]} & \py{results_table[5][1]} & \py{results_table[5][2]} \\
\bottomrule
\end{tabular}
\end{table}

\chapter{Conclusions}

\begin{enumerate}
    \item Filtered back-projection requires sufficient angular sampling
    \item Ram-Lak filter provides best resolution but amplifies noise
    \item Apodizing filters trade resolution for noise reduction
    \item Metal artifacts cause streak patterns
    \item Limited angle causes directional blurring
    \item Iterative methods can improve reconstruction quality
\end{enumerate}

\end{document}