Computerized Tomography

High-Dynamic-Range Computerized Tomography. Computerized Tomography (CT) is among the most revolutionary technologies in medical imaging, allowing the creation of 3D images from X-ray measurements. At the core of CT systems is the Radon transform, which links the image parameters to the X-ray measurements. Beyond CT, the Radon transform finds applications in geophysical explorations and remote sensing, as well. The existing approach HDR CT is based on the fusion of multiple low-dynamic-range images. This approach may suffer from ghosting artefacts and miss-calibration of the exposure and sensor responses. In addition, it is a quantitative approach that lacks mathematical guarantees. For example, how many images are required for sufficient HDR reconstruction?

Forward Model and the Modulo Radon Transform. The US-CT is a single-shot HDR CT, based on modulo sampling. At the core of this approach stands the modulo Radon transform (MRT), arising from the hardware implementation of the system, where instead of recording the Radon projections, the modulo of Radon projections are computed and sampled. The forward model of the US-CT can be written as,

\[\mathscr{R}^{\lambda}f(\theta,t)=\mathscr{M}_{\lambda}(\mathscr{R}f(\theta, t)).\]

where \(\mathscr{M}_{\lambda}\) is the modulo operator, \(\mathscr{R}\) is te Radon operator, and \(\theta, t\) denote angle and time.

CT Forward Model.

Fig. 31 CT Forward Model.

Despite the non-linearity the many-to-one mapping due to the modulo operator, the MRT can be injective, and identifiable from semidiscrete MRT samples. This motivates the derivation of recovery conditions and algorithm.

The Inverse Problem. To reconstruct an image, the MRT samples need to be inverted. The conditions and algorithms that can be used for sampling that allows reconstruction can be classified based on the domain of operation.

Spatial-Domain Inversion. For a semi-discrete MRT model, where the angles are available on the continuum, exact recovery from MRT samples can be guaranteed under the given conditions.

Unlimited sampling–filtered back projection

Any \(f \in L^{1}(\mathbb{R}^2)\cap\mathcal{B}_{\Omega}(\mathbb{R}^2)\) can be uniquely recovered from semidiscrete MRT samples if \(T < 1 / \Omega\) and the reconstruction filter \(F_{\Omega}\) is chosen to be the Ram–Lak filter with window \(W = \mathbf{1}_{[-1,1]}\).

The recovery is based on US-FBP [17], a variant of the US-ALG.

  1. Choose \(N = \left\lfloor \frac{\log(\lambda) - \log(\beta_f)}{\log(T S \varrho)} \right\rfloor\), \(J = \frac{6 \beta_f}{\lambda}\)

  2. Compute \(s_{(0)}[k] = \left(\mathscr{M}_\lambda \left( \Delta^N y_\theta^{\lambda} \right) - \Delta^N y_\theta^{\lambda} \right) [k]\)

  3. For \(n = 0, \dots, N - 2\) do:

    1. \(s_{(n+1)}[k] = 2 \lambda \left\lfloor \frac{s_{s(n)}[k] / \lambda}{2} \right\rfloor\)

    2. \(\kappa_n = \left\lfloor \frac{s_{(n+1)}[1] - s_{(n+1)}[J+1]}{12 \beta_f} + \frac{1}{2} \right\rfloor\)

    3. \(s_{(n+1)}[k] = s_{(n+1)}[k] + 2 \lambda \kappa_n\)

  4. End For

  5. \(y_\theta[k] = y_\theta^{\lambda}[k] + (S s_{(N-1)})[k]\)

  6. \(\mathcal{R}_\theta f = \sum_{k \in \mathbb{Z}} y_\theta[k] \, \operatorname{sinc} \left( \frac{\pi}{T} (\cdot - k T) \right)\)

Output: US-FBP reconstruction \(f_\Omega^{\lambda} = \frac{1}{2} \mathcal{R}^\# \left( F_\Omega * \mathcal{R}_\theta f \right)\)

In the case when \(f\) is not bandlimited but belongs to the more general Sobolev spaces of fractional order \(\alpha > 0\), given by \(H^{\alpha} (\mathbb{R}^2) = \bigl\{ f \in L^2(\mathbb{R}^2) \bigm| \| f \|_{\alpha} < \infty \bigr\}\), the reconstruction error using the US-FBP can be upper bounded.

If only finitely many time samples and angles are available, the image can be approximated based on a sequential application of modulo unfolding and a discrete FBP reconstruction formula.

US-CT Acquisition and Recovery Pipeline.

Fig. 32 US-CT Acquisition and Recovery Pipeline.

Hardware Experiment The proposed US-CT was experimentally tested using the M-ADC. In the experiment, and image from The Walnut Dataset was re-digitized using the M-ADC over \(20\) angles, and then algorithmically recovered. The technique achieves a reconstruction error at the order of \(10^{-3}\), with dynamic range enhancement of \(AAAA\). The sampled waveforms and the reconstruction from modulo samples for a fixed angle are presented below.

US-CT Hardware Experiment

Fig. 33 US-CT Hardware Experiment

Fourier-Domain Inversion Implementing spatial-domain inversion poses several challenges for real-time applications. These include the numerical instability when inverting higher differences, the high oversampling factor, and the need for prior knowledge of \(\lambda\). A Fourier-domain approach offers advantages by harnessing a wide range of tomographic recovery techniques to address these issues. Yet, working with non-linearities is particularly difficult in the transform domain.

OMP-DFR The first algorithm, OMP-DFR (Orthogonal Matching Pursuit - Direct Fourier Reconstruction), is a mixture of a Fourier domain and spatial domain reconstruction approaches: #. Reconstruction of the Fourier Radon samples based on a Fourier-domain separation principle (see US-FP) and a sparsity property. This is done using the OMP algorithm. #. Reconstruction of the target function using filtered back projection formula.

The reconstruction in first part can be guaranteed, while the error due to the FBP can be bounded.

US-CT OMP-based Approaches.

Fig. 34 US-CT OMP-based Approaches.