Three-dimensional Neumann-series approach to model light transport in nonuniform media

Image reconstruction in fluorescence molecular tomography with sparsity-initialized maximum-likelihood expectation maximization

We present a reconstruction method involving maximum-likelihood expectation maximization (MLEM) to model Poisson noise as applied to fluorescence molecular tomography (FMT). MLEM is initialized with the output from a sparse reconstruction-based approach, which performs truncated singular value decomposition-based preconditioning followed by fast iterative shrinkage-thresholding algorithm (FISTA) to enforce sparsity. The motivation for this approach is that sparsity information could be accounted for within the initialization, while MLEM would accurately model Poisson noise in the FMT system. Simulation experiments show the proposed method significantly improves images qualitatively and quantitatively. The method results in over 20 times faster convergence compared to uniformly initialized MLEM and improves robustness to noise compared to pure sparse reconstruction. We also theoretically justify the ability of the proposed approach to reduce noise in the background region compared to pure sparse reconstruction. Overall, these results provide strong evidence to model Poisson noise in FMT reconstruction and for application of the proposed reconstruction framework to FMT imaging.

Incorporating reflection boundary conditions in the Neumann series radiative transport equation: application to photon propagation and reconstruction in diffuse optical imaging

We propose a formalism to incorporate boundary conditions in a Neumann-series-based radiative transport equation. The formalism accurately models the reflection of photons at the tissue-external medium interface using Fresnel's equations. The formalism was used to develop a gradient descent-based image reconstruction technique. The proposed methods were implemented for 3D diffuse optical imaging. In computational studies, it was observed that the average root-mean-square error (RMSE) for the output images and the estimated absorption coefficients reduced by 38% and 84%, respectively, when the reflection boundary conditions were incorporated. These results demonstrate the importance of incorporating boundary conditions that model the reflection of photons at the tissue-external medium interface.

A radiative transfer equation-based image-reconstruction method incorporating boundary conditions for diffuse optical imaging

Developing reconstruction methods for diffuse optical imaging requires accurate modeling of photon propagation, including boundary conditions arising due to refractive index mismatch as photons propagate from the tissue to air. For this purpose, we developed an analytical Neumann-series radiative transport equation (RTE)-based approach. Each Neumann series term models different scattering, absorption, and boundary-reflection events. The reflection is modeled using the Fresnel equation. We use this approach to design a gradient-descent-based analytical reconstruction algorithm for a three-dimensional (3D) setup of a diffuse optical imaging (DOI) system. The algorithm was implemented for a three-dimensional DOI system consisting of a laser source, cuboidal scattering medium (refractive index > 1), and a pixelated detector at one cuboid face. In simulation experiments, the refractive index of the scattering medium was varied to test the robustness of the reconstruction algorithm over a wide range of refractive index mismatches. The experiments were repeated over multiple noise realizations. Results showed that by using the proposed algorithm, the photon propagation was modeled more accurately. These results demonstrated the importance of modeling boundary conditions in the photon-propagation model.

A three-step reconstruction method for fluorescence molecular tomography based on compressive sensing

Fluorescence molecular tomography (FMT) is a promising tool for real time in vivo quantification of neurotransmission (NT) as we pursue in our BRAIN initiative effort. However, the acquired image data are noisy and the reconstruction problem is ill-posed. Further, while spatial sparsity of the NT effects could be exploited, traditional compressive-sensing methods cannot be directly applied as the system matrix in FMT is highly coherent. To overcome these issues, we propose and assess a three-step reconstruction method. First, truncated singular value decomposition is applied on the data to reduce matrix coherence. The resultant image data are input to a homotopy-based reconstruction strategy that exploits sparsity via ℓ₁ regularization. The reconstructed image is then input to a maximum-likelihood expectation maximization (MLEM) algorithm that retains the sparseness of the input estimate and improves upon the quantitation by accurate Poisson noise modeling. The proposed reconstruction method was evaluated in a three-dimensional simulated setup with fluorescent sources in a cuboidal scattering medium with optical properties simulating human brain cortex (reduced scattering coefficient: 9.2 cm^-1, absorption coefficient: 0.1 cm^-1) and tomographic measurements made using pixelated detectors. In different experiments, fluorescent sources of varying size and intensity were simulated. The proposed reconstruction method provided accurate estimates of the fluorescent source intensity, with a 20% lower root mean square error on average compared to the pure-homotopy method for all considered source intensities and sizes. Further, compared with conventional ℓ₂ regularized algorithm, overall, the proposed method reconstructed substantially more accurate fluorescence distribution. The proposed method shows considerable promise and will be tested using more realistic simulations and experimental setups.

Quantitative photoacoustic tomography using forward and adjoint Monte Carlo models of radiance

Forward and adjoint Monte Carlo (MC) models of radiance are proposed for use in model-based quantitative photoacoustic tomography. A two-dimensional (2-D) radiance MC model using a harmonic angular basis is introduced and validated against analytic solutions for the radiance in heterogeneous media. A gradient-based optimization scheme is then used to recover 2-D absorption and scattering coefficients distributions from simulated photoacoustic measurements. It is shown that the functional gradients, which are a challenge to compute efficiently using MC models, can be calculated directly from the coefficients of the harmonic angular basis used in the forward and adjoint models. This work establishes a framework for transport-based quantitative photoacoustic tomography that can fully exploit emerging highly parallel computing architectures.

Radiance and photon noise: imaging in geometrical optics, physical optics, quantum optics and radiology

The statistics of detector outputs produced by an imaging system are derived from basic radiometric concepts and definitions. We show that a fundamental way of describing a photon-limited imaging system is in terms of a Poisson random process in spatial, angular, and wavelength variables. We begin the paper by recalling the concept of radiance in geometrical optics, radiology, physical optics, and quantum optics. The propagation and conservation laws for radiance in each of these domains are reviewed. Building upon these concepts, we distinguish four categories of imaging detectors that all respond in some way to the incident radiance, including the new category of photon-processing detectors (capable of measuring radiance on a photon-by-photon basis). This allows us to rigorously show how the concept of radiance is related to the statistical properties of detector outputs and to the information content of a single detected photon. A Monte-Carlo technique, which is derived from the Boltzmann transport equation, is presented as a way to estimate probability density functions to be used in reconstruction from photon-processing data.

Singular value decomposition for photon-processing nuclear imaging systems and applications for reconstruction and computing null functions

Recent advances in technology are enabling a new class of nuclear imaging systems consisting of detectors that use real-time maximum-likelihood (ML) methods to estimate the interaction position, deposited energy, and other attributes of each photon-interaction event and store these attributes in a list format. This class of systems, which we refer to as photon-processing (PP) nuclear imaging systems, can be described by a fundamentally different mathematical imaging operator that allows processing of the continuous-valued photon attributes on a per-photon basis. Unlike conventional photon-counting (PC) systems that bin the data into images, PP systems do not have any binning-related information loss. Mathematically, while PC systems have an infinite-dimensional null space due to dimensionality considerations, PP systems do not necessarily suffer from this issue. Therefore, PP systems have the potential to provide improved performance in comparison to PC systems. To study these advantages, we propose a framework to perform the singular-value decomposition (SVD) of the PP imaging operator. We use this framework to perform the SVD of operators that describe a general two-dimensional (2D) planar linear shift-invariant (LSIV) PP system and a hypothetical continuously rotating 2D single-photon emission computed tomography (SPECT) PP system. We then discuss two applications of the SVD framework. The first application is to decompose the object being imaged by the PP imaging system into measurement and null components. We compare these components to the measurement and null components obtained with PC systems. In the process, we also present a procedure to compute the null functions for a PC system. The second application is designing analytical reconstruction algorithms for PP systems. The proposed analytical approach exploits the fact that PP systems acquire data in a continuous domain to estimate a continuous object function. The approach is parallelizable and implemented for graphics processing units (GPUs). Further, this approach leverages another important advantage of PP systems, namely the possibility to perform photon-by-photon real-time reconstruction. We demonstrate the application of the approach to perform reconstruction in a simulated 2D SPECT system. The results help to validate and demonstrate the utility of the proposed method and show that PP systems can help overcome the aliasing artifacts that are otherwise intrinsically present in PC systems.

RADIANCE AND PHOTON NOISE: Imaging in geometrical optics, physical optics, quantum optics and radiology

A fundamental way of describing a photon-limited imaging system is in terms of a Poisson random process in spatial, angular and wavelength variables. The mean of this random process is the spectral radiance. The principle of conservation of radiance then allows a full characterization of the noise in the image (conditional on viewing a specified object). To elucidate these connections, we first review the definitions and basic properties of radiance as defined in terms of geometrical optics, radiology, physical optics and quantum optics. The propagation and conservation laws for radiance in each of these domains are reviewed. Then we distinguish four categories of imaging detectors that all respond in some way to the incident radiance, including the new category of photon-processing detectors. The relation between the radiance and the statistical properties of the detector output is discussed and related to task-based measures of image quality and the information content of a single detected photon.

Joint reconstruction of activity and attenuation map using LM SPECT emission data

Attenuation and scatter correction in single photon emission computed tomography (SPECT) imaging often requires a computed tomography (CT) scan to compute the attenuation map of the patient. This results in increased radiation dose for the patient, and also has other disadvantages such as increased costs and hardware complexity. Attenuation in SPECT is a direct consequence of Compton scattering, and therefore, if the scattered photon data can give information about the attenuation map, then the CT scan may not be required. In this paper, we investigate the possibility of joint reconstruction of the activity and attenuation map using list-mode (LM) SPECT emission data, including the scattered-photon data. We propose a path-based formalism to process scattered-photon data. Following this, we derive analytic expressions to compute the Cramér-Rao bound (CRB) of the activity and attenuation map estimates, using which, we can explore the fundamental limit of information-retrieval capacity from LM SPECT emission data. We then suggest a maximum-likelihood (ML) scheme that uses the LM emission data to jointly reconstruct the activity and attenuation map. We also propose an expectation-maximization (EM) algorithm to compute the ML solution.

Truncated Fourier-series approximation of the time-domain radiative transfer equation using finite elements

The radiative transfer equation (RTE) is widely accepted to accurately describe light transport in a medium with scattering particles, and it has been successfully applied as a light-transport model, for example, in diffuse optical tomography. Due to the computationally expensive nature of the RTE, most of these applications have been in the frequency domain. In this paper, an efficient solution method for the time-domain RTE is proposed. The method is based on solving the frequency-domain RTE at multiple modulation frequencies and using the Fourier-series representation of the radiance to obtain approximation of the time-domain solution. The approach is tested with simulations. The results show that the method can be used to obtain the solution of the time-domain RTE with good accuracy and with significantly fewer computational resources than are needed in the direct time-domain solution.

An ideal-observer framework to investigate signal detectability in diffuse optical imaging

With the emergence of diffuse optical tomography (DOT) as a non-invasive imaging modality, there is a requirement to evaluate the performance of the developed DOT systems on clinically relevant tasks. One such important task is the detection of high-absorption signals in the tissue. To investigate signal detectability in DOT systems for system optimization, an appropriate approach is to use the Bayesian ideal observer, but this observer is computationally very intensive. It has been shown that the Fisher information can be used as a surrogate figure of merit (SFoM) that approximates the ideal observer performance. In this paper, we present a theoretical framework to use the Fisher information for investigating signal detectability in DOT systems. The usage of Fisher information requires evaluating the gradient of the photon distribution function with respect to the absorption coefficients. We derive the expressions to compute the gradient of the photon distribution function with respect to the scattering and absorption coefficients. We find that computing these gradients simply requires executing the radiative transport equation with a different source term. We then demonstrate the application of the SFoM to investigate signal detectability in DOT by performing various simulation studies, which help to validate the proposed framework and also present some insights on signal detectability in DOT.

Simulating photon-transport in uniform media using the radiative transport equation: a study using the Neumann-series approach

We present the implementation, validation, and performance of a Neumann-series approach for simulating light propagation at optical wavelengths in uniform media using the radiative transport equation (RTE). The RTE is solved for an anisotropic-scattering medium in a spherical harmonic basis for a diffuse-optical-imaging setup. The main objectives of this paper are threefold: to present the theory behind the Neumann-series form for the RTE, to design and develop the mathematical methods and the software to implement the Neumann series for a diffuse-optical-imaging setup, and, finally, to perform an exhaustive study of the accuracy, practical limitations, and computational efficiency of the Neumann-series method. Through our results, we demonstrate that the Neumann-series approach can be used to model light propagation in uniform media with small geometries at optical wavelengths.