Exercises

ex01-fourier-slice

Easy

(a) Create a 128×128128 \times 128 Shepp-Logan phantom.

(b) Compute the Radon transform (projection) at angle θ=30∘\theta = 30^\circ.

(c) Take the 1D FFT of the projection and the 2D FFT of the phantom. Extract the radial slice of the 2D FFT at 30∘30^\circ.

(d) Plot both on the same axes. Verify they are identical (up to numerical precision).

(e) Repeat for θ=0∘,45∘,90∘,135∘\theta = 0^\circ, 45^\circ, 90^\circ, 135^\circ. At which angle is the numerical agreement worst? Why?

Show Hint

Use interpolation (e.g., scipy.ndimage.map_coordinates) to extract the radial slice from the 2D FFT grid.

Worst agreement typically occurs at 45 degrees where the interpolation path cuts diagonally across grid cells.

Solution
Setup

Use skimage.data.shepp_logan_phantom() or construct manually. Compute projections with skimage.transform.radon().

Verification

The 1D FFT of the projection at θ\theta and the radial slice of the 2D FFT at angle θ\theta should match to within numerical precision (<10−10< 10^{-10} relative error for 64-bit arithmetic). Diagonal angles have slightly worse agreement due to interpolation across the Cartesian FFT grid.

ex02-fbp-basic

Easy

(a) Implement FBP for X-ray CT using the Ram-Lak filter. Use the Shepp-Logan phantom (256×256256 \times 256).

(b) Reconstruct with Nv=180N_v = 180 projections. Compute the PSNR.

(c) Implement the Shepp-Logan and cosine filters. Reconstruct with Nv=36N_v = 36 and compare all three filters in PSNR and visual quality.

(d) Plot the three filter frequency responses h(ν)h(\nu). Explain why the cosine filter suppresses streak artifacts.

(e) What is the minimum NvN_v for PSNR >25> 25 dB with each filter?

Show Hint

The Ram-Lak filter is h(ν)=∣ν∣h(\nu) = |\nu|; discretize carefully to avoid the DC singularity.

The cosine filter rolls off high frequencies, reducing aliased energy that creates streaks.

Solution
Filter responses

Ram-Lak: h(ν)=∣ν∣h(\nu) = |\nu|. Shepp-Logan: h(ν)=∣ν∣ sinc(ν/(2νmax⁡))h(\nu) = |\nu|\,\text{sinc}(\nu / (2\nu_{\max})). Cosine: h(ν)=∣ν∣ cos⁡(πν/(2νmax⁡))h(\nu) = |\nu|\,\cos(\pi\nu/(2\nu_{\max})). The cosine filter suppresses high spatial frequencies where angular undersampling causes aliasing (streaks).

Minimum views for 25 dB

Typically: Ram-Lak needs Nv≈60N_v \approx 60; Shepp-Logan Nv≈40N_v \approx 40; cosine Nv≈30N_v \approx 30. The smoother filters tolerate sparser sampling.

ex03-ewald-geometry

Easy

(a) Write a function that computes the Ewald circle K(ϕs)=k0(k^s(ϕs)−k^i)\mathbf{K}(\phi_s) = k_0(\hat{\mathbf{k}}_s(\phi_s) - \hat{\mathbf{k}}_i) for a given incident direction k^i\hat{\mathbf{k}}_i and wavenumber k0k_0.

(b) Plot the Ewald circles for Nv=1,2,4,8N_v = 1, 2, 4, 8 views at f=10f = 10 GHz. Label each circle with its view angle.

(c) For each configuration, compute the coverage ratio η\eta (fraction of the 2k02k_0 disk that is covered).

(d) At what NvN_v does η\eta exceed 0.9 for a single frequency?

(e) Add Nf=5N_f = 5 frequencies (8--12 GHz). How does the coverage change?

Show Hint

Discretize each Ewald circle into approximately 360 points and compute coverage on a grid.

Solution
Ewald circle formula

For incident direction θi\theta_i: K(ϕs)=k0(cos⁡ϕs−cos⁡θi,  sin⁡ϕs−sin⁡θi)\mathbf{K}(\phi_s) = k_0(\cos\phi_s - \cos\theta_i,\; \sin\phi_s - \sin\theta_i). The circle has center (−k0cos⁡θi,−k0sin⁡θi)(-k_0\cos\theta_i, -k_0\sin\theta_i) and radius k0k_0.

Coverage analysis

With a single frequency, η>0.9\eta > 0.9 requires Nv≳30N_v \gtrsim 30 (the arcs are thin curves with limited area). Adding Nf=5N_f = 5 frequencies increases coverage substantially because each view now fills an annular band rather than a single arc.

ex04-diffraction-tomography

Medium

(a) Create a 2D object (64×6464 \times 64) with 3 circular inclusions of varying contrast (χ=0.1,0.3,0.5\chi = 0.1, 0.3, 0.5) in a homogeneous background.

(b) Simulate the scattered field under the Born approximation for Nv=18N_v = 18 views and Nf=10N_f = 10 frequencies (8--12 GHz). Add noise at SNR=30\text{SNR} = 30 dB.

(c) Reconstruct using: (i) gridding + IFFT with DCF, (ii) CGLS (100 iterations).

(d) Compare PSNR and visual quality. Which method best preserves the circular boundaries?

(e) Reduce to Nv=6N_v = 6 views. How do the results change? Which method degrades more gracefully?

Show Hint

For Born simulation, compute the forward integral via FFT convolution with the Green's function.

CGLS handles incomplete data better because it implicitly regularizes via early stopping.

Solution
Forward simulation

The Born model gives Es=k02F−1[G~(k) χ~(k−k0k^i)]E_s = k_0^2 \mathcal{F}^{-1}[ \tilde{G}(\mathbf{k})\,\tilde{\chi}(\mathbf{k} - k_0\hat{\mathbf{k}}_i)]. Implement via FFT-based convolution for each view and frequency.

Reconstruction comparison

At Nv=18N_v = 18, both methods produce reasonable reconstructions (PSNR ≈25\approx 25--2828 dB). CGLS slightly better preserves edges. At Nv=6N_v = 6, gridding shows severe streaks (PSNR ≈16\approx 16 dB) while CGLS remains acceptable (PSNR ≈20\approx 20 dB) due to implicit regularization from limited iterations.

ex05-coverage-optimization

Medium

(a) For a fixed budget of Nv=8N_v = 8 views, find the set of angles {θv}\{\theta_v\} that maximizes the coverage ratio η\eta (single frequency, f=10f = 10 GHz).

(b) Compare: (i) uniform spacing, (ii) random angles, (iii) greedy selection (add the angle that maximally increases η\eta).

(c) Plot η\eta vs. NvN_v for all three strategies.

(d) Repeat for wideband (8--12 GHz, Nf=10N_f = 10). Does the optimal angle selection change?

(e) For limited-angle imaging (θv∈[0∘,90∘]\theta_v \in [0^\circ, 90^\circ] only), which Fourier regions are always missing?

Show Hint

Greedy maximization of coverage ratio is a set cover problem; it gives a (1−1/e)(1 - 1/e)-optimal solution.

Solution
Greedy algorithm

Initialize K=∅\mathcal{K} = \emptyset. For v=1,…,Nvv = 1, \ldots, N_v: choose θv=arg⁡max⁡θ∣K∪E(θ)∣−∣K∣\theta_v = \arg\max_\theta |\mathcal{K} \cup \mathcal{E}(\theta)| - |\mathcal{K}| where E(θ)\mathcal{E}(\theta) is the Ewald arc for angle θ\theta.

Result analysis

Greedy and uniform give similar η\eta for single-frequency (both ≈0.35\approx 0.35 for Nv=8N_v = 8). For wideband, greedy slightly outperforms uniform because it can exploit the radial diversity to place views where angular gaps matter most. Limited-angle imaging always misses the Fourier quadrant opposite to the view range.

ex06-density-compensation

Medium

(a) For Nv=12N_v = 12 views and Nf=8N_f = 8 frequencies, compute the Ewald arc sample positions on a 128×128128 \times 128 Fourier grid.

(b) Implement the Voronoi-based DCF: compute the Voronoi cell area for each sample point.

(c) Implement the iterative (Pipe-Menon) DCF with 20 iterations.

(d) Reconstruct a test object with and without the DCF. Compare PSNR and visual quality.

(e) Plot the DCF as a function of ∣K∣|\mathbf{K}|. How does it compare to the theoretical ramp filter ∣K∣|\mathbf{K}|?

Show Hint

Use scipy.spatial.Voronoi; clip infinite cells at the 2kmax⁡2k_{\max} boundary.

Without DCF, the reconstruction is biased toward low spatial frequencies.

Solution
DCF effect

Without DCF, low spatial frequencies are overweighted (they appear in multiple overlapping arcs), producing a blurred reconstruction with PSNR ≈3\approx 3--55 dB lower. The Pipe-Menon DCF closely approximates the Voronoi DCF and produces reconstructions with uniform frequency weighting.

Comparison with ramp filter

The DCF approximately follows ∣K∣|\mathbf{K}| at low spatial frequencies (similar to the ramp filter in FBP) but deviates at high frequencies where the Ewald arc density is non-uniform.

ex07-nufft-reconstruction

Medium

(a) Implement gridding + FFT reconstruction using a Kaiser-Bessel interpolation kernel (W=6W = 6, oversampling α=2\alpha = 2).

(b) Implement NUFFT-based reconstruction using the finufft library (or equivalent).

(c) Compare both with direct summation on a 32×3232 \times 32 problem. Verify they produce the same result (within tolerance).

(d) Benchmark computation time for grid sizes N=32,64,128,256,512N = 32, 64, 128, 256, 512. At what NN does gridding become faster than direct summation?

(e) Vary the Kaiser-Bessel width W=2,4,6,8W = 2, 4, 6, 8 and measure reconstruction error vs. computation time.

Show Hint

The deapodization step divides by the Fourier transform of the Kaiser-Bessel kernel at each image point.

Gridding typically becomes faster than direct summation at N≈64N \approx 64.

Solution
Crossover point

Direct summation is O(MN)O(MN) where MM is the number of k-space samples. Gridding + FFT is O(MW2+Nlog⁡N)O(MW^2 + N\log N). For typical M∼104M \sim 10^4 and W=6W = 6, gridding wins above N≈64N \approx 64.

Kernel width tradeoff

W=2W = 2: fast but inaccurate (relative error ∼10−2\sim 10^{-2}). W=6W = 6: standard choice (error ∼10−6\sim 10^{-6}). W=8W = 8: marginally better accuracy but 80% more gridding operations. W=6W = 6 is the standard tradeoff.

ex08-born-validity

Hard

(a) Simulate scattering from a dielectric cylinder of radius aa and contrast χ\chi using the exact Mie series solution (2D).

(b) Compute the Born approximation prediction for the same cylinder.

(c) Plot the relative error ∥Esexact−EsBorn∥/∥Esexact∥\|E_s^{\text{exact}} - E_s^{\text{Born}}\| / \|E_s^{\text{exact}}\| as a function of the scattering parameter α=k0∣χ∣a\alpha = k_0 |\chi| a.

(d) At what α\alpha does the Born error exceed 10%? 20%?

(e) Reconstruct the object using FDT + gridding with exact and Born data. At what α\alpha does the reconstruction PSNR drop below 20 dB?

(f) Compare with the Rytov approximation. For which objects is Rytov more accurate than Born?

Show Hint

The 2D Mie series uses Bessel and Hankel functions of integer order.

Rytov is more accurate for large objects with small refractive index contrast.

Solution
Born validity threshold

The Born approximation typically breaks at α≈0.5\alpha \approx 0.5--1.01.0 (10% error) and fails badly at α>2\alpha > 2 (50%+ error). The Rytov approximation extends the validity to α≈2\alpha \approx 2--33 for smooth objects because it correctly handles phase accumulation even when the total phase shift is large.

Reconstruction impact

PSNR drops below 20 dB at α≈1.5\alpha \approx 1.5 for Born and α≈3\alpha \approx 3 for Rytov. Beyond these thresholds, iterative methods (DBIM) are needed.

ex09-near-field-tomography

Hard

(a) Simulate near-field scattering from 3 point scatterers using the exact Green's function (spherical waves) at f=10f = 10 GHz. Measurement array: 32 elements on a circle of radius 20 cm.

(b) Reconstruct using: (i) far-field FDT (plane-wave model), (ii) near-field model (spherical-wave back-propagation).

(c) Compute position error and PSNR for both methods.

(d) Vary the array radius from 10 cm to 2 m. Plot PSNR vs. distance for both methods. At what distance do they converge?

(e) Implement the NF-FF transformation for a planar measurement surface. Apply it to the near-field data, then reconstruct with the far-field FDT. Compare with direct near-field reconstruction.

Show Hint

The NF-FF transformation is a 2D spatial FFT, propagation phase correction, then IFFT.

Both methods should converge around r≈2D2/λr \approx 2D^2/\lambda.

Solution
Near-field advantage

At 20 cm range (r/dff≈0.3r/d_{\mathrm{ff}} \approx 0.3), the far-field model gives PSNR ≈14\approx 14 dB with 8 mm position errors. The near-field model gives PSNR ≈29\approx 29 dB with sub-millimeter position errors.

Convergence distance

Both methods converge (PSNR difference <1< 1 dB) at r≈dff=2D2/λr \approx d_{\mathrm{ff}} = 2D^2/\lambda. For D=20D = 20 cm at 10 GHz: dff≈2.7d_{\mathrm{ff}} \approx 2.7 m.

ex10-multi-freq-synthesis

Hard

(a) Implement bandwidth synthesis: combine Nf=20N_f = 20 narrowband measurements at frequencies uniformly spaced from 8 to 12 GHz.

(b) Verify that the range resolution matches δr=c/(2W)=3.75\delta_r = c/(2W) = 3.75 cm for W=4W = 4 GHz. Place two point scatterers at varying separations and find the minimum resolvable separation.

(c) What happens if the frequency spacing exceeds c/(2D)=500c/(2D) = 500 MHz (for a 30 cm object)? Demonstrate range ambiguity.

(d) Add 10% random phase errors across frequencies (simulating oscillator drift). How does this affect the range resolution?

(e) Implement phase-error calibration using a single known reference scatterer. Show that the phase errors can be estimated and corrected.

Show Hint

Range ambiguity appears as ghost targets at multiples of c/(2Δf)c/(2\Delta f).

Phase calibration uses the known scatterer position to estimate the phase offset at each frequency.

Solution
Bandwidth synthesis verification

With coherent combination of 20 frequencies over 4 GHz bandwidth, the range profile shows a sinc-like peak with mainlobe width matching c/(2W)=3.75c/(2W) = 3.75 cm. Two scatterers are resolvable at separations >3.75> 3.75 cm.

Range ambiguity

When Δf>c/(2D)=500\Delta f > c/(2D) = 500 MHz, ghost targets appear at range offsets of c/(2Δf)c/(2\Delta f). For Δf=1\Delta f = 1 GHz: ghosts at ±15\pm 15 cm -- within the 30 cm object.

Phase calibration

The reference scatterer at known position provides the expected phase at each frequency. The difference between expected and measured phase gives the frequency-dependent phase error, which can be subtracted from all measurements.

ex11-resolution-analysis

Hard

(a) For the configuration: f0=10f_0 = 10 GHz, W=4W = 4 GHz, Nf=20N_f = 20, Nv=36N_v = 36 (full 360∘360^\circ), compute and plot the 2D PSF.

(b) Measure the −3-3 dB mainlobe width in range and cross-range. Compare with the theoretical predictions δr=c/(2W)\delta_r = c/(2W) and δcr=λmin⁡/4\delta_{\text{cr}} = \lambda_{\min}/4.

(c) Measure the peak sidelobe level. Compare between Ram-Lak, Shepp-Logan, and cosine filters.

(d) Repeat for limited aperture: Δθ=30∘,60∘,120∘,180∘\Delta\theta = 30^\circ, 60^\circ, 120^\circ, 180^\circ. Plot resolution (range and cross-range) vs. aperture angle.

(e) Place two point scatterers at varying separations. Determine the minimum resolvable separation (Rayleigh criterion) for each configuration.

Show Hint

The PSF is the inverse FFT of the filtered coverage indicator function.

Compute the PSF on a fine grid (at least 4×4\times oversampled) for accurate mainlobe measurement.

Solution
PSF measurement

For W=4W = 4 GHz, full aperture: range mainlobe ≈3.8\approx 3.8 cm (close to theoretical 3.753.75 cm); cross-range mainlobe ≈0.64\approx 0.64 cm (close to λmin⁡/4=0.625\lambda_{\min}/4 = 0.625 cm at 12 GHz).

Filter comparison

Peak sidelobe levels: Ram-Lak ≈−13\approx -13 dB, Shepp-Logan ≈−18\approx -18 dB, cosine ≈−23\approx -23 dB. The tradeoff is sidelobe level vs. mainlobe width (cosine has ≈20%\approx 20\% wider mainlobe).

ex12-limited-angle-regularized

Hard

(a) Create a 64×6464 \times 64 scene with 5 objects (mix of point and extended targets). Simulate DT data with views limited to θ∈[0∘,90∘]\theta \in [0^\circ, 90^\circ], Nv=9N_v = 9, Nf=10N_f = 10, SNR=25\text{SNR} = 25 dB.

(b) Reconstruct with gridding + IFFT. Identify the artifact pattern and explain it in terms of missing Fourier data.

(c) Reconstruct with FISTA + ℓ1\ell_1 regularization.

(d) Reconstruct with ADMM + TV regularization.

(e) Compare all three methods in PSNR and visual quality. Which best fills the missing Fourier data?

(f) Increase the aperture to 180∘180^\circ, then 360∘360^\circ. At what aperture does gridding match the regularized methods?

Show Hint

Limited-angle tomography produces an elongated PSF along the missing directions.

TV regularization is particularly effective for extended targets with sharp edges.

Solution
Artifact analysis

Quarter-aperture (90∘90^\circ) leaves three-quarters of the angular Fourier data missing. The gridding reconstruction shows severe directional streaks and elongation of point targets along the missing directions.

Regularization comparison

FISTA (ℓ1\ell_1): PSNR ≈22\approx 22 dB; good for point targets, less effective for extended objects. ADMM (TV): PSNR ≈24\approx 24 dB; best for extended targets with edges; mild staircase artifacts. Gridding matches regularized methods at approximately 240∘240^\circ--300∘300^\circ aperture.

ex13-xlmimo-nearfield

Hard

(a) Model an XL-MIMO array with D=1D = 1 m aperture at 28 GHz (λ=1.07\lambda = 1.07 cm), 128 elements. Compute dffd_{\mathrm{ff}}.

(b) Place 5 point scatterers at range r=20r = 20 m. Simulate the received signal using the exact near-field Green's function.

(c) Reconstruct using: (i) far-field beamforming, (ii) near-field beamfocusing.

(d) Compute and compare the PSFs for both methods at ranges r=5,20,50,100,200r = 5, 20, 50, 100, 200 m.

(e) At what range does far-field beamforming become acceptable (PSNR within 1 dB of near-field)?

Show Hint

Near-field beamfocusing uses range-dependent phase compensation: the steering vector depends on both angle and range.

The transition occurs near r≈dff/4r \approx d_{\mathrm{ff}}/4 for practical purposes.

Solution
Far-field distance

dff=2(1)2/0.0107=187d_{\mathrm{ff}} = 2(1)^2/0.0107 = 187 m. At r=20r = 20 m, the target is at r/dff=0.11r/d_{\mathrm{ff}} = 0.11 -- deep in the Fresnel zone.

Transition analysis

Far-field beamforming becomes acceptable (PSNR within 1 dB) at approximately r≈50r \approx 50 m (r/dff≈0.27r/d_{\mathrm{ff}} \approx 0.27). This is earlier than the formal Fraunhofer criterion because imaging tolerates some PSF broadening.

ex14-full-pipeline

Challenge

Build an end-to-end diffraction tomography imaging system:

(a) Scene: A 20×2020 \times 20 cm region with 3 dielectric cylinders (radii 1, 2, 3 cm; contrasts χ=0.1,0.2,0.3\chi = 0.1, 0.2, 0.3) plus 5 point scatterers.

(b) System: 24 Tx and 24 Rx on a 30 cm circular array. Frequency: 6--14 GHz, 32 stepped frequencies.

(c) Forward model: Born approximation with near-field Green's functions (the array is in the Fresnel zone).

(d) Reconstruction pipeline: (i) Compute coverage map and PSF. (ii) Gridding + IFFT with DCF (fast preview). (iii) FISTA + ℓ1\ell_1 for sparse refinement. (iv) ADMM + TV for edge preservation of extended targets.

(e) Evaluation: PSNR, SSIM, and computation time for each stage.

(f) Sensitivity analysis: Vary SNR\text{SNR} from 10 to 40 dB. At what SNR\text{SNR} does gridding become acceptable (PSNR >20> 20 dB)? At what SNR\text{SNR} does regularization become essential?

(g) Near-field impact: Compare reconstruction with and without near-field correction. Quantify the PSNR loss from ignoring phase curvature.

Show Hint

Use the NUFFT for efficient Fourier-domain operations throughout the pipeline.

The near-field correction is most important for the innermost array elements.

Progressive refinement: gridding provides initialization for FISTA, which provides initialization for ADMM+TV.

Solution
Pipeline performance

Typical results at SNR=25\text{SNR} = 25 dB: Gridding: PSNR ≈22\approx 22 dB, 0.1 s. FISTA: PSNR ≈28\approx 28 dB, 5 s. ADMM+TV: PSNR ≈30\approx 30 dB, 15 s. Near-field correction adds ≈3\approx 3--55 dB at this configuration.

Sensitivity thresholds

Gridding acceptable (>20> 20 dB) at SNR≥20\text{SNR} \geq 20 dB. Regularization becomes essential (gridding PSNR <15< 15 dB) at SNR≤15\text{SNR} \leq 15 dB.

ex15-fbp-cnn

Challenge

(a) Train a simple CNN (3--5 convolutional layers) to post-process FBP reconstructions from limited-view DT.

(b) Training data: Generate 500 random phantoms (random circles and rectangles), simulate DT data (Nv=12N_v = 12, Nf=10N_f = 10, SNR=25\text{SNR} = 25 dB), reconstruct with gridding. Pairs: (gridding output, ground truth).

(c) Evaluate on 50 test phantoms. Compare: (i) gridding alone, (ii) gridding + CNN, (iii) FISTA (ℓ1\ell_1).

(d) How does the CNN perform on out-of-distribution phantoms (e.g., phantoms with triangular objects, not seen in training)?

(e) Discuss: is this approach (physics initialization + learned refinement) more or less robust than end-to-end learning? Connect to the FBPConvNet paradigm.

Show Hint

Use a U-Net architecture (3 levels) for the post-processing CNN.

Out-of-distribution generalization is a key limitation of purely data-driven methods.

Solution
Expected results

Gridding + CNN typically achieves PSNR ≈27\approx 27--3030 dB, comparable to or slightly better than FISTA, with much faster inference (∼10\sim 10 ms vs. ∼5\sim 5 s). The CNN effectively learns to suppress streak artifacts.

Generalization

On out-of-distribution phantoms, the CNN's PSNR drops by ≈3\approx 3--55 dB, while FISTA's performance is largely unchanged. This motivates hybrid approaches (Chapter 13) that combine physics-based algorithms with learned components.

Chapter SummaryReferences & Further Reading