Perturbation Monte Carlo (pMC) has been previously proposed to rapidly recompute optical measurements when small perturbations of optical properties are considered, but it was largely restricted to changes associated with prior tissue segments or regions-of-interest. In this work, we expand pMC to compute spatially and temporally resolved sensitivity profiles, i.e. the Jacobians, for diffuse optical tomography (DOT) applications. By recording the pseudo random number generator (PRNG) seeds of each detected photon, we are able to “replay” all detected photons to directly create the 3D sensitivity profiles for both absorption and scattering coefficients. We validate the replay-based Jacobians against the traditional adjoint Monte Carlo (aMC) method, and demonstrate the feasibility of using this approach for efficient 3D image reconstructions using in vitro hyperspectral wide-field DOT measurements. The strengths and limitations of the replay approach regarding its computational efficiency and accuracy are discussed, in comparison with aMC, for point-detector systems as well as wide-field pattern-based and hyperspectral imaging systems. The replay approach has been implemented in both of our open-source MC simulators - MCX and MMC (http://mcx.space)
© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
Diffuse optical tomography (DOT) and fluorescence molecular tomography (FMT) with near-infrared light have been active research areas for both clinical and pre-clinical applications for over three decades. The appeal of these tomographic techniques resides in their high sensitivity, relatively low cost and unique utilization of intrinsic and exogenous functional biomarkers in charactering optically-thick tissue samples noninvasively . These methods aim to recover the 3D distributions of either fluorophores or intrinsic absorption and scattering contrasts inside turbid media by employing measurement data collected on the surface of tissues and match them with an accurate light propagation forward model through solving an inverse problem . Over the past years, DOT and FMT have been successfully applied to imaging human breast cancers , prostate cancers , muscle functions , joint conditions , peripheral vascular diseases , brain functions , small animals , and have also found applications in imaging ex vivo biological samples  as well as engineered tissues . Regardless of the application, building an accurate and efficient forward model remains crucial to the quantitative imaging performance of inverse model-based techniques.
It has been widely agreed that the radiative transfer equation (RTE) provides an accurate model for light propagation in complex bio-tissues but it is known to be difficult to solve. Therefore, the diffusion equation (DE), a low-order expansion to the RTE, is usually adopted as an alternative to generate approximated solutions . With the introduction of more sophisticated imaging strategies over the past decade, such as mesoscopic imaging , time-domain imaging with early photons  and whole-body imaging with structured light [15,16], DE becomes increasingly limited in both accuracy and flexibility. In these cases, the Monte Carlo (MC) method, a stochastic solver to the RTE, has gained increasing popularity.
Because MC solves for the RTE directly, in many studies, it was considered as the “gold standard” for modeling photon migration in biological tissues [17,18]. Nevertheless, MC was not deemed suitable for providing the forward solutions in DOT/FMT due primarily to two limitations. First, only limited types of domain configurations, such as layered and voxelated space, are supported by public MC packages, restricting their use in modeling complex tissues. This limitation was addressed by the recent development of mesh-based MC methods [19,20]. Second, MC has been associated with high computational expenses when simulating large number of photons especially in 3D applications. Thanks to the wide availability of modern accelerators such as general purpose graphics processing unit (GPGPU), the simulation time of MC has been reduced by orders of magnitude over the past years . Additionally, efficient MC formulations such as the perturbation Monte Carlo (pMC) have been proposed to speed up the forward computation of MC [22,23]. Hence, the combination of efficient MC formulations with increased computation power has led to a rejuvenated interest in using MC to compute the forward model for DOT/FMT [24,25].
A Jacobian matrix, also known as the sensitivity matrix, is a central element in model-based DOT/FMT image reconstruction. It represents a map associating the changes of localized optical properties with the changes of measurements from selected source-detector pairs. Traditionally, the adjoint Monte Carlo (aMC) is adopted as an efficient method to calculate the Jacobians, especially when optodes are point-sources and point-detectors [25, 26]. The adjoint method works by multiplying the light fields of a forward simulation, computed by propagating photons from the source position, and an adjoint simulation, computed by propagating photons from an imaginary source at the detector position . However, with the rise of more advanced DOT imaging architectures, such as those leveraging structured light illumination or detection , the aMC method becomes cumbersome with the rapid increase in data dimensionality via spatial, temporal and spectral multiplexing. The investigations to these expanded data dimensions have been driving the development of next generation DOT/FMT systems and resulted in significantly improved 3D imaging performances [29, 30]. In many of these studies, pMC was found to be well suited as it benefits greatly from rescaling strategies that do not require relaunching MC forward simulations. However, studies utilizing pMC have been primarily focused on analyzing domains with known segmentations; the perturbations being estimated are restricted to known tissue segments or regions-of-interest, and are not spatially distributed. Extending pMC to support spatially resolved image reconstructions can greatly benefit DOT in exploring additional measurement dimensions.
In this paper, we propose an approach to efficiently build the Jacobians via perturbation Monte Carlo with photon numerical “replay”, where the trajectory information of detected photons is first recorded from a “baseline” simulation and then utilized to build sensitivity profiles. By storing the pseudo random number generator (PRNG) seeds and recreating the photon trajectories in an MC simulation , we can avoid the prohibitive memory cost of storing photon “biographies” and efficiently generate spatially and temporally resolved Jacobians with only a small overhead. Additional optimization strategies are investigated for hyperspectral imaging and single-pixel detection to further accelerate the Jacobian calculation.
The remaining sections of the paper are organized as follows. In Section 2, we will briefly review perturbation Monte Carlo principles and derive spatial/temporal Jacobians for both absorption or scattering perturbations from pMC formulations. We will also provide a brief review on adjoint Monte Carlo. In Section 3, we will detail the “photon replay” algorithm for direct Jacobian constructions. We will particularly discuss two scenarios - camera-based wide-field detection and hyperspectral imaging. In Section 4, we validate the Jacobians from replay with those from the adjoint method - first provide a simple in silico DOT example as a demonstration, and then another example in solving an in vitro hyperspectral wide-field DOT problem. Detailed discussions and comparisons between the replay and adjoint methods are provided, followed by conclusions and future improvements in Section 5.
2. Perturbation Monte Carlo and adjoint Monte Carlo for DOT
2.1. Forward and inverse perturbation Monte Carlo
When a photon is launched from a source in a Monte Carlo simulation, it undergoes multiple scattering events before reaching a detector (or exits without being detected), and its weight decreases exponentially along the trajectory following the Beer-Lambert law . Sometimes, the optical properties, i.e. absorption coefficient (μa), scattering coefficient (μs), and anisotropy factor (g), can be perturbed, (μ̂a, μ̂s, ĝ), from their background values (μa, μs, g) in a small region, such as in a tumor. In such cases, instead of running a new MC simulation with perturbed optical properties, we can predict the measurements by processing the recorded photon “biographies” from the background simulation, known as the perturbation Monte Carlo method [22,23]. For simplicity, here we assume g is not perturbed. In such case, the detected photon weight can be expressed as:19,21], referred to as the “microscopic Beer-Lambert law” approach (mBLL) in . In this method, a photon’s trajectory is solely determined by μs and the photon weight loss is solely determined by μa. This is different from the “albedo-weight” approach  commonly used in pMC literature [23,24], but was shown to be equivalent to the latter . For a complex domain, we rewrite Eq. (1) into a spatially resolved form [24,34] as:
If we perturb the absorption or scattering coefficient one region (Ωj) at a time, we can obtain the corresponding change δW due to the perturbation as:4–5) and keep only the first order terms. Then we can compute the Jacobians in the form of limits as: 6–7) can be further simplified to
2.2. The Adjoint Monte CarloEq. (10) when using the Born normalized approach. For simplicity, we will focus on DOT hereinafter but interested readers are referred to .
In many realistic DOT systems, the quantity being measured is not the fluence, but the total diffuse reflectance or transmittance (R) at a detector immediately outside the domain , i.e.
If the measurement is the fluence inside the diffusive media, i.e. Φ(rd, rj) = G(rd, rj), it can then be simply replaced by the adjoint forward solution of the fluence G(rj, rd) as a result of reciprocity of the Green’s function, i.e. Φ(rd, rj) = G(rd, rj) = G(rj, rd). However, when the measured quantity is diffuse reflectance or normal flux on a refractive index mismatched boundary, reciprocity is not valid any more. In such cases, the relationship between R(rd, r), or Jn, and the adjoint fluence Green’s function G(rj, rd) is more complex, and is related to 1) the boundary condition, 2) the exiting photon angular distribution , 3) the asymmetry between the emitting profile of the adjoint source and the receiving profile of the detector, and 4) surface specular reflection. Although deriving an analytical relationship is possible in such cases for a planar boundary within the diffusion regime [12,39], the measurement can typically be expressed as a scaled version of the adjoint fluence solution, i.e. Φ(rd, rj) = αG(rj, rd). The scaling factor α can be numerically calibrated by two sets of MC simulations. Also note that ∇Φ (rd, rj) = −α∇G(rj, rd) because the two gradients have opposite directions. Thus, Eqs. (10–11) can be ultimately re-written as
2.3. Solving the DOT inverse problem
After computing the sensitivity matrices from either the adjoint or perturbation MC, the DOT inverse problem typically requires solving the below equation35].
Due to the highly scattering nature of near-infrared light in the tissue, Eq. (16) is known to be ill-posed . As a result, one must apply regularization to stabilize the solution. The standard Tikhonov regularization is often used. When there are more measurements than unknowns, i.e. over-determined, solving the below equation is more efficient :
3. Building Jacobian matrices via photon replay
3.1. Photon replay algorithm
Based on Eqs. (8–9), the main challenge for building Jacobians with pMC is to find an efficient way to calculate the weighted average scattering event counts p̄ and photon trajectory lengths L̄, both in time and space. A simple approach is to allocate a memory buffer with twice the numbers of voxels or elements, B, multiplied by the number of time gates, T, for each simulated photon, to store p̄(Ωj, T) and L̄(Ωj, T) during the propagation. However, the drawbacks of this method are apparent. First, the memory for storing 2B×T values needs to be pre-allocated for each detector. When the number of detectors is large, which is typical for DOT applications, or multiple photons are simulated simultaneously, such as GPU-based MC, the memory usage becomes unrealistic, especially for GPUs where memory utility is highly restrictive. Secondly, we do not know in advance if a photon will arrive at the detector. In many cases, only a small portion of the photons are captured, resulting in significant redundant calculations.
The “photon replay” method was inspired by the “reseeded” pMC method proposed by Sassaroli . By taking advantage of the deterministic pseudo-random number generator (RNG) sequence followed by a seed, one only needs to store the photon’s initial RNG state at launch-time when it is captured by a detector. By reinitializing the photon with the saved RNG state, the photon will repeat its trajectory as in the baseline simulation and is guaranteed to be detected again at the same position. We therefore propose a photon “replay” algorithm combining the reseeded MC with the Jacobian formulations in Eqs. (8–9), with a schematic shown in Fig. 1.
The photon replay algorithm is performed in two steps. First, a baseline Monte Carlo simulation is performed where a large number of photons with random seeds are launched. If a photon is captured by one of the specified detectors, its arrival time t, final weight w, RNG seed value s at launch, and the detector index are stored. After the number of photons reach the target, a second Monte Carlo simulation, i.e. the photon replay, is performed for each detector. This time, we initialize each new photon with the stored seed value sk, from the kth detected photon in the baseline simulation, so that they are guaranteed to be captured again. During its re-propagation, instead of depositing photon weight losses as in the baseline MC, we utilize the stored detected photon weight wk and deposit either weighted trajectory lengths wk Lj (“wl” mode) or weighted scattering numbers wk pj (“wp” mode) at any given discretized region (Ωj) and time gate (binned by photon’s total time-of-flight tk), as described in Eqs. (6–7). After all detected photons are replayed, the weighted average values L̄(rj, t) and p̄(rj, t) are generated and used to compute the temporal and spatial resolved Jacobians according to Eqs. (8–9).
3.2. Improved replay efficiency in wide-field and hyperspectral imaging
With the development of spatial light modulators (SLMs), structured light illumination such as spatial frequency domain imaging (SFDI) and camera based single-pixel detection strategies have gained increasing popularity . This novel system design has shown great potential to improve data acquisition speed as well as signal-to-noise ratio (SNR) for DOT/FMT, making full-body small-animal imaging more efficient. Because the SNR of the Jacobians from replay largely depends on the number of detected photons, a single-pixel detection over a large field-of-view (FOV) can potentially benefit more from replay compared to point detectors.
Typically, dozens of detection patterns belong to the same class are used, such as 36 patterns in the quantized low-frequency base and 64 patterns in the Hadamard base . Since they share the same detection area, we can run only one replay MC to generate Jacobians for all detection patterns with any given source configuration. To achieve this, the only additional data that need to be stored is a photon’s exiting bilinear position (xk, yk) in the detector aperture. Then, when depositing wk Lj and wk pj according to Eqs. (6–7), we simply loop over all detection patterns and modify wk to , where is the corresponding “detector weight”, defined as the ith detection pattern intensity (0 to 1) at (xk, yk). This approach is well suited for emerging methodologies in DOT/FMT such as compressive-sensing based or adaptive DOT where “optimal” patterns can be obtained either theoretically  or experimentally .
Another scenario where we can leverage the benefits of photon replay is multispectral or hyperspectral imaging, which provides enriched dataset at different wavelengths for unmixing various types of chromophores. The optical properties, either μa or μs, are usually slightly different as the wavelength changes, so we have to perform separate MC simulations at each wavelength if adjoint MC were to be used to build Jacobians. With the proposed replay method, we only need to run one baseline simulation and one replay for all wavelength channels, with the assumption that μs is a constant over the measured wavelength range. Both of our published MC simulators, MCX [21,44] and MMC [19,45], utilize the mBLL method. Therefore, as long as μs remains constant, we can exactly repeat the photon trajectory with the stored seeds and recalculate the detected photon weight for a new μa value using Eq. (1), without needing to rerun the baseline simulation . Applying the rescaled wk in the replay will obtain the Jacobians of different wavelengths as proposed in .
4.1. Validation of replay-derived Jacobians
We first compare the Jacobians calculated by the adjoint method (aMC) and photon replay for a single source-detector pair in a homogeneous domain. All simulations are performed with our GPU-accelerated MC simulator “MCX”  using an NVIDIA GTX TITAN V GPU. The domain has a volume of 60×60×60 mm3 and a voxel size of 1×1×1 mm3, with the optical properties μa= 0.005 mm−1, μs = 1 mm−1, g = 0 and n = 1.37. A pencil beam light source towards +z direction is placed at [29.5, 29.5, 0] mm and a detector with a radius 2.5 mm is located at [29.5, 39.5, 0] mm. To match the configurations of the two simulations, the light source for the adjoint MC simulation is set to a uniform disk source with 2.5 mm radius at the detector position, where the photons are launched uniformly within the disk towards −z direction. The constant α in Eqs. (14–15) is determined to be 8.47 through calibration.
For aMC, 108 photons were launched for both the forward and adjoint simulations, each with a runtime of ∼ 1.90 s. For the replay approach, the baseline simulation with 109 photons took 22.6 s, in which around 7.30×106 photons were detected. The replay steps to obtain the “weighted trajectory length” (wl) and “weighted scattering number” (wp) took 0.39 s and 0.30 s, respectively. All runtimes reported are the average value of 5 repeated tests. To reduce stochastic noise, we apply Anscombe transform on the raw outputs, i.e. G(rj, rs), G(rj, rd) in aMC and L̄(Ωj), p̄(Ωj) in replay, before filtering them with a 3-D Gaussian kernel with σ = 1, followed by an exact inverse Anscombe transform [46,47].
In Fig. 2, we show the absorption and scattering Jacobians from the adjoint approach in (a–c), those from the replay approach are in (d–f), along with contour plot comparisons between the two methods in (h–i). Distributions of weighted average pathlengths L̄(Ωj) (same as the replay adjoint Jacobian by definition) and weighted average scattering numbers p̄(Ωj)/μs (Ωj) in (d, g), respectively. All plots are in log-scale with base 10 extracted at plane x = 29.5 mm.
Several remarks can be made from Fig. 2. First, the distributions of L̄(rj) in Fig. 2(d) and p̄(Ωj) / μs (Ωj) in Fig. 2(g) are very similar to each other. This can be explained by the following relationship: the average distance traveled by a single photon between two scattering events can be calculated as l = L/p, and the photon “mean free path” happens to be the inverse of scattering coefficient (assuming g = 0) in an mBLL MC photon simulation. Therefore, the accumulated L̄(Ωj) and p̄(Ωj)/μs (Ωj) values for a collection of photons are also very close numerically at the same spatial location. Secondly, the overall magnitude of the scattering Jacobian is about 1–2 orders of magnitude lower than the absorption Jacobian. This is expected because the scattering Jacobian is calculated by subtracting images in (g) by (d), according to Eq. (9). For the same reason, the scattering Jacobian is much nosier than the absorption Jacobian, as the stochastic noise is relatively amplified with the smaller mean value.
Finally, in the scattering Jacobian, only a small region near the source-detector is negative while in majority of the space, it has a positive value. This is also expected. Mathematically, since ∇Φ(rj, rs) and ∇Φ (rj, rd) vectors are along the and directions, from Eq. (15) the sign of Jμs (Ωj) is determined by the inner product of these two vectors, so we can see the points with negative scattering sensitivity values are inside the circle with the diameter of |rs − rd|. Physically, for photons launched from the source, a higher scattering coefficient along or close to the source-detector path tends to raise the possibility of deflecting photons away from the detector, thus reducing the total detected photon weight; on the other hand, a higher scattering coefficient far from the source-detector path tends to reflect photons back to the detector, thus increasing the total detected photon weight. Note here that, although g = 0 is used in this example, aMC and replay Jacobians are also expected to match when g is non-zero.
In addition, we also show the time-resolved absorption Jacobians Jμa (Ωj, T), built with the replay method for a transmission measurement setup, in Fig. 3. The size of the phantom is 60×60×30 mm3 with the same optical properties as above. The source remains at [29.5, 29.5, 0] mm and the detector locates at [29.5, 39.5, 30] mm. We recorded over 5 ns with 0.1 ns time step after launching 109 photons in the MC simulation, where the baseline and replay MC with 7.46×105 took 29.4 s and 0.11 s, respectively. We see that Jμa in an early time-gate (Fig. 3(a)) has a significantly narrower span compared to a late time gate with the same accumulated photon weights W(T), shown in Fig. 3(b), due to the ballistic photons with fewer scatterings . Moreover, in Fig. 3(c), we compare the change of temporal point-spread functions (TPSFs) predicted by pMC and replay Jacobians, for a 20×20×10 mm3 inclusion with δµa = 10−4 mm−1 at the center of the phantom (20 ≤ x, y ≤ 40 mm, 10 ≤ z ≤ 20 mm). The TPSFs derived from pMC (as Eq. (3)) and the replay Jacobian, calculated as Jμa · Δμa, match almost exactly. This is not surprising as both were derived from the same set of detected photons.
4.2. Sample reconstructions using the replay Jacobians
In this section, we test the reconstruction performance of both absorption and scattering Jacobians obtained from photon replay with simulated data. The tested slab domain contains 40×40×20 isotropic 1×1×1 mm3 voxels. The optical properties for the background medium are μa = 0.01 mm−1, μs = 1 mm−1, g = 0 and n = 1.37. Although g is close to 0.9 in many biological tissues, in the diffusion regime, we used the equivalent settings, i.e. g = 0 and μs as the reduced scattering coefficient (μ′s), to speed up MC computation. Two 8 mm diameter tubes are placed side-by-side with their axes 10 mm below the surface and 16 mm apart. We varied the optical properties of the two cylindrical inclusions and ran multiple simulations. The light source is a collimated Gaussian beam with a waist radius of 2 mm, and 8×8 raster scanning positions are applied to the phantom bottom plane (z = 0 mm) within a rectangle region bounded by [6, 6, 0] mm and [34, 34, 0] mm. A total of 36 disk-shaped detectors of 3 mm radius are placed on the top plane (z = 20 mm) with centers located between 5 to 35 mm, with 6 mm spacing, in both x and y axes.
To generate the absorption and scattering Jacobians, we launched 109 photons from every source position in the baseline simulation and performed photon replays for each detector. On an NVIDIA TITAN V GPU, the average run-time for each baseline simulation was 9.63 s, where ∼ 4.1×107 photons on average were detected; the runtimes for replaying 36 detectors were 10.55 s and 7.20 s for the “wl” and “wp” steps, respectively, resulting in a total run-time of ∼ 29.2 minutes. The raw outputs were denoised as shown in Section 4.1 before constructing the full Jacobian of 2,304 (source-detector pairs) and 32,000 columns (voxel count). Meanwhile, we also generated 3 sets of simulated CW measurements (Φ0, Φ1 and Φ2), where Φ0 : the medium is homogeneous with background optical properties; Φ1 : the left tube has 3 times the background absorption μ̂a = 0.03 mm−1 and Φ2 : the right tube has twice the background scattering μ̂s = 2 mm−1. For each measurement set, 108 photons were used for simulations at each source position, taking an average of 1.05 s per simulation.
The inverse problems for absorption and scattering perturbations were solved based on Eq. (18) and optimal regularization parameters were chosen based on an L-curve approach. In Fig. 4, we show the reconstructed μa cross-sections with Born approximation in (a–c) and those with Rytov approximation in (d–f). The absorption and scattering reconstructions (positive component) on z = 10 mm and y = 20 mm planes are displayed in the first and second columns separately, and the comparison between the ground truth and reconstructed mean values along y dimension on z = 10 mm plane are plotted in (c,f). Note that only one of μa or μs was perturbed in each sample reconstruction; the case where both μa and μs are perturbed is more complicated, normally requiring time-resolved measurements, matrix rescaling or other techniques, and is thus beyond the scope of this paper.
From Fig. 4, with either approximation, we have successfully reconstructed the true locations of the inclusions. However, the results from the Rytov approach appears to be more accurate, as evident when comparing Fig. 4(f) with (c). We also notice that the reconstructed scattering perturbation is noisier than absorption, as shown both on the surface and interior. This is expected because the scattering Jacobian is noisy as shown in Section 4.1. To suppress artefacts and obtain more accurate images, one can apply regularizations, adopt more advanced inverse solvers or utilize a priori structural information, which again are not the focus of this paper. We have demonstrated in silico that building Jacobians and solving inverse problem via photon replay are feasible for a typical DOT problem with many source-detector pairs.
4.3. Replay for hyperspectral wide-field DOT
Recently, we reported a novel hyperspectral time-resolved wide-field imaging system, capable of collecting in vivo small-animal whole-body data within 14 minutes . One digital mirror device (DMD) is fiber-coupled with a pulsed laser to generate structured illumination patterns and a second DMD is coupled with a spectrophotometer-based photomultiplier tube detector (PMT) to enable single-pixel detection. In this section, we show that one can use the photon replay approach to efficiently generate Jacobians for this type of hyperspectral DOT system.
A liquid phantom of size 80×50×22.5 mm3 was prepared as the imaging target. The background medium contains 0.8% Intralipid, 0.008‰ Epolight 2735 (Ep), and 0.024‰ India Ink (Ink). Two tubular inclusions with 5 mm radius were suspended 10 mm below the surface, as shown with the black dash lines in Figs. 5(a–b). Both were filled with the same Intralipid concentration as the background, while the left one is 2.25× higher in Ep and the right one is 1.83× higher in Ink concentration. Therefore, both tubes have the same μ′s as the background but different μa values, and their ratio with respect to the unit concentrations is 1.23.
The experiment was performed under a transmission configuration, where the field-of-view (FOV) was set to 40×30 mm2 for both illumination (bottom) and detection (top). A total of 36 quantized low-frequency (QLF) illumination patterns and detection patterns were used in this experiment. These patterns was chosen for their high SNR at low photon counts .
CW data from 5 spectral channels, with the center wavelength ranging from 750 to 770 nm, were selected for unmixing two absorbers. Wavelength dependent optical properties of the background and two absorbing materials were calibrated for every channel. The μ′s value remained approximately constant at 0.84 mm−1, while μa of background medium increased from 0.015 mm−1 to 0.016 mm−1, μa of 0.008‰ Ep increased from 7.2×10−3 mm−1 to 8.6×10−3 mm−1, and μa of 0.024‰ Ink decreased from 9.5×10−3 mm−1 to 9.1×10−3 mm−1. We assume that all μa values change linearly with the wavelength, as shown in Fig. 5(c).
In this example, we used MMC  for both the baseline and replay calculations using 36 CPU nodes in parallel on the Blue Gene/Q cluster, each consisting of a 16-core 1.6 GHz A2 processor at the Center for Computational Innovations (CCI) at Rensselaer Polytechnic Institute (RPI). The phantom mesh consisted of 17,073 nodes and 100,241 elements . The baseline MC with 108 photons and μa = 0.015 mm−1 took 16.8 minutes on average on the CPU cluster, where 5.3×106 photons were collected by the detector. With the technique introduced in Section 3.2, a single replay MC session created Jacobians for all 36 detection patterns in 15.7 minutes. Furthermore, under a constant scattering coefficient, we ran five replay simulations for different wavelength channels. Therefore, the total simulation time is 16.8 + 15.7×5 = 95.3 min.
Rytov approximation was again used to compute the measurement vector, and an optode calibration method was adopted  to reduce artifacts caused by fluctuations of pattern intensities in the experiment. The inverse problem was solved with a MATLAB built-in solver CGS with the default tolerance . The reconstructed Cep and Cink at 50% isovolume are shown in Figs. 5(a–b). Despite distortions along the z-axis - an expected result of the fixed projection angle - the two absorbers were successfully differentiated with accurate locations and low surface artefacts. We calculate the reconstructed absorber concentrations as the mean of elements larger than 50% of the maximum value, and define the crosstalk of Ep as , where C̄ep and C̄ink are normalized to their own maximum and Ωep is the ground truth volume of Ep (vice versa for Ink). From Fig. 5(d), we found that with measurements from more wavelength channels, the crosstalk between Epolight 2735 and India Ink kept decreasing (from 100% to 50.4% and 18.0%), whereas the concentration ratio between two absorbers kept approaching the ground truth (123%). More details about further improving resolution/ratio and reducing crosstalk can be found at .
5. Discussions and conclusions
In this section, we discuss the pros and cons of the replay and adjoint approaches. First, we discuss the accuracy of the Jacobians from the two approaches. The replay Jacobians was derived based on pMC, built upon the direct connection between the individual detected photon data and the actual measurement changes. By contrast, the adjoint method relies on the reciprocity of light . Thus, when sufficient number of photons are detected, we consider the replay to be more accurate than the adjoint method. We recognize that several factors may contribute to the derivations between the two Jacobians. Firstly, the adjoint method requires assuming an imaginary source that is located at the detector position. Isotropic or pencil-beam sources have been widely used as the adjoint sources, but the light emission profiles from these source forms are obviously different from the light reception profile of common detectors, such as an optical-fiber. To accurately simulate an adjoint source, a spatially and angularly distributed light source must be used, which is not commonly available. Secondly, the formulation of the adjoint Jacobian depends on the actual type of the measurement, such as fluence, flux or diffuse reflectance. Choosing different measurement quantities can result in different scaling factors, α, as we discussed in Section 2.2. In reality, the scaling factor due to the chosen measurement is effectively removed as part of the data calibration .
Next, we want to compare the two methods in terms of computational efficiency. If we follow the procedures described above, if launching M photons at each of the Ns sources and Nd non-overlapping detectors, the run-time for the adjoint method satisfies tad ∝ [(Ns + Nd) M]; that for the replay method satisfies tre ∝ [(1 + f) NsM], where f is the fraction of the photons being detected in the baseline simulation. However, this comparison does not consider the differences of the two Jacobians in terms of SNRs due to their differences in the number of photons involved. Assuming a shot-noise model, we found that the SNR of the adjoint Jacobian is proportional to while that for the replay is proportional to . Therefore, to produce a Jacobian that matches the SNR of the adjoint computation, the replay method needs to launch times more photons in the baseline simulation. This leads to an equivalent replay time of tre ∝ [2NdNs (1 + f) M/f].
In the cases where the detected photons are shared among the detectors, such as in the single-pixel camera setup, if we assume that every photon is on average shared between 1 ≤ gd ≤ Nd detectors (we call gd as the detection reusability gain), the average SNR of the replay becomes , and the resulting replay time is with the assumption that the memory operation cost is neglectable compared to photon propagation. Similarly, in multi- /hyper-spectral imaging with only absorption perturbations, the detected photons can be replayed once for all Nλ spectral channels, leading to a new running time of . Lastly, in applications where iterative reconstruction is performed, the baseline simulation is only needed once at the first iteration if the replay method is used. The forward model outputs and Jacobians of the subsequent iterations can be computed rapidly using only the detected photons, yielding an overall for K iterations. In comparison, the adjoint method requires repeated simulations at all sources and detectors in all iterations, giving tad ∝ (Ns + Nd) KM.
In Fig. 6, we show 3 sample configurations to demonstrate the pros and cons of the replay algorithm in terms of computational expense. The contour plots of the runtime ratios between the replay and the adjoint methods for building Jacobians, i.e. tre/tad, are presented, with a range of source and detector numbers 1 ≤ Ns, Nd ≤ 20. In Fig. 6(a), we show a case with single wavelength (Nλ = 1), non-overlapping point detectors (gd = 1) and f = 0.1. In Fig. 6(b), we show the runtime ratios for a multispectral single-pixel camera detection system with gd = Nd, Nλ = 10 and f = 0.1. In Fig. 6(c), we show an example of iterative reconstruction (K = 10) for a hyperspectral (Nλ = 10) point-detector (gd = 1) system with f = 0.3.
From Fig. 6(a), it appears that if the Jacobian SNR is the major concern, the replay method performs poorly in terms of speed compared with the adjoint method in point-source/detector configurations: when Ns = Nd = 10, the replay takes 100–120 times longer than the adjoint to build a Jacobian with same SNR. In fact, the larger the source and detector numbers, the more efficient is the adjoint method. However, according to Fig. 6(b), when more detected photons are shared between detection channels, either in the form of spatial patterns or spectral channels, the replay method can be more efficient under the condition that . Even for point-detection systems, utilizing iterative reconstruction and multispectral data can also make the replay method more effective compared to the adjoint, as suggested by Fig. 6(c).
In summary, we have proposed a novel approach named “photon replay” to directly and conveniently build time-resolved Jacobians for both absorption and scattering perturbations for diffuse optical tomography. The replay Jacobians were validated with the traditional adjoint approach and show efficiency in several scenarios, such as hyperspectral and single-pixel detection systems. We also demonstrated that the absorption Jacobians from replay can be practically useful for solving real-world DOT reconstructions; however, the scattering Jacobian suffers from stochastic noise despite large photon counts and additional denoising. We are currently working on new strategies to improve the quality of the scattering Jacobian.
The “replay” feature has been implemented and tested in both our voxel- and mesh-based MC platforms. The updated software, MCX and MMC, are freely available at http://mcx.space/. We strongly believe that these tools could contribute to a wide spectrum of explorations in diffuse optical imaging.
National Institute of Health (NIH) (R01-GM114365, R01-CA204443 (QF), and R01-EB19443, R01 BRG CA207725 (XI)); National Science Foundation (NSF) Career Award CBET 1149407 (XI); NVIDIA GPU Grant Program.
The authors declare that there are no known conflicts of interest related to this article.
References and links
2. S. R. Arridge and J. C. Schotland, “Optical tomography: forward and inverse problems,” Inverse Probl. 25, 123010 (2009). [CrossRef]
3. D. Grosenick, H. Rinneberg, R. Cubeddu, and P. Taroni, “Review of optical breast imaging and spectroscopy,” J. of Biomed. Optics 21, 091311 (2016). [CrossRef]
4. D. Piao, K. E. Bartels, Z. Jiang, G. R. Holyoak, J. W. Ritchey, G. Xu, C. F. Bunting, and G. Slobodov, “Alternative transrectal prostate imaging: a diffuse optical tomography method,” IEEE J. Sel. Top. Quantum Electron. 16, 715–729 (2010). [CrossRef]
5. D. Contini, L. Zucchelli, L. Spinelli, M. Caffini, R. Re, A. Pifferi, R. Cubeddu, and A. Torricelli, “Brain and muscle near infrared spectroscopy/imaging techniques,” J. Near Infrared Spectrosc. 20, 15–27 (2012). [CrossRef]
6. A. H. Hielscher, H. K. Kim, L. D. Montejo, S. Blaschke, U. J. Netz, P. A. Zwaka, G. Illing, G. A. Muller, and J. Beuthan, “Frequency-domain optical tomographic imaging of arthritic finger joints,” IEEE Trans. Med. Imaging 30, 1725–1736 (2011). [CrossRef] [PubMed]
7. M. Khalil, H. Kim, J. Hoi, I. Kim, R. Dayal, G. Shrikhande, and A. Hielscher, “Detection of peripheral arterial disease within the foot using vascular optical tomographic imaging: a clinical pilot study,” Eur. J. Vasc. Endovasc. Surg. 49, 83–89 (2015). [CrossRef]
8. A. T. Eggebrecht, S. L. Ferradal, A. Robichaux-Viehoever, M. S. Hassanpour, H. Dehghani, A. Z. Snyder, T. Hershey, and J. P. Culver, “Mapping distributed brain function and networks with diffuse optical tomography,” Nat. Photonics 8, 448 (2014). [CrossRef] [PubMed]
9. A. H. Hielscher, “Optical tomographic imaging of small animals,” Current Opinion in Biotech. 16, 79–88 (2005). [CrossRef]
10. M. Pimpalkhare, J. Chen, V. Venugopal, and X. Intes, “Ex vivo fluorescence molecular tomography of the spine,” J. of Biomed. Imaging 2012, 942326 (2012).
11. M. S. Ozturk, V. K. Lee, L. Zhao, G. Dai, and X. Intes, “Mesoscopic fluorescence molecular tomography of reporter genes in bioprinted thick tissue,” J. of Biomed. Optics 18, 100501 (2013). [CrossRef]
12. R. C. Haskell, L. O. Svaasand, T.-T. Tsay, T.-C. Feng, M. S. McAdams, and B. J. Tromberg, “Boundary conditions for the diffusion equation in radiative transfer,” J. Opt. Soc. Am. A 11, 2727–2741 (1994). [CrossRef]
13. M. S. Ozturk, C.-W. Chen, R. Ji, L. Zhao, B.-N. B. Nguyen, J. P. Fisher, Y. Chen, and X. Intes, “Mesoscopic fluorescence molecular tomography for evaluating engineered tissues,” Ann. Biomed. Eng. 44, 667–679 (2016). [CrossRef]
14. M. J. Niedre, R. H. de Kleine, E. Aikawa, D. G. Kirsch, R. Weissleder, and V. Ntziachristos, “Early photon tomography allows fluorescence detection of lung carcinomas and disease progression in mice in vivo,” Proc Natl Acad Sci U S A 105,19126–19131 (2008). [CrossRef] [PubMed]
16. N. Ducros, C. D’andrea, G. Valentini, T. Rudge, S. Arridge, and A. Bassi, “Full-wavelet approach for fluorescence diffuse optical tomography with structured illumination,” Opt. Lett. 35, 3676–3678 (2010). [CrossRef] [PubMed]
18. C. Zhu and Q. Liu, “Review of Monte Carlo modeling of light transport in tissues,” J. of Biomed. Optics 18, 050902 (2013). [CrossRef]
20. H. Shen and G. Wang, “A tetrahedron-based inhomogeneous Monte Carlo optical simulator,” Phys. Med. & Biol. 55, 947 (2010). [CrossRef]
22. A. Sassaroli, C. Blumetti, F. Martelli, L. Alianelli, D. Contini, A. Ismaelli, and G. Zaccanti, “Monte Carlo procedure for investigating light propagation and imaging of highly scattering media,” Appl. Opt. 37, 7392–7400 (1998). [CrossRef]
23. C. K. Hayakawa, J. Spanier, F. Bevilacqua, A. K. Dunn, J. S. You, B. J. Tromberg, and V. Venugopalan, “Perturbation Monte Carlo methods to solve inverse photon migration problems in heterogeneous tissues,” Opt. Lett. 26, 1335–1337 (2001). [CrossRef]
26. A. R. Gardner, C. K. Hayakawa, and V. Venugopalan, “Coupled forward-adjoint Monte Carlo simulation of spatial-angular light fields to determine optical sensitivity in turbid media,” J. of Biomed. Optics 19, 065003 (2014). [CrossRef]
27. Q. Fang, P. M. Meaney, S. D. Geimer, A. V. Streltsov, and K. D. Paulsen, “Microwave image reconstruction from 3-D fields coupled to 2-D parameter estimation,” IEEE Trans. on Medical Imaging 23, 475–484 (2004). [CrossRef] [PubMed]
28. J. P. Angelo, S.-J. Chen, M. Ochoa, U. Sunar, S. Gioux, and X. Intes, “Review of Structured Light in Diffuse Optical Imaging,” J. of Biomed. Optics (in press).
29. Q. Pian, R. Yao, L. Zhao, and X. Intes, “Hyperspectral time-resolved wide-field fluorescence molecular tomography based on structured light and single-pixel detection,” Opt. Lett. 40, 431–434 (2015). [CrossRef] [PubMed]
32. E. Alerstam, S. Andersson-Engels, and T. Svensson, “White Monte Carlo for time-resolved photon migration,” J. of Biomed. Optics 13, 041304 (2008). [CrossRef]
33. A. Sassaroli and F. Martelli, “Equivalence of four Monte Carlo methods for photon migration in turbid media,” J. Opt. Soc. Am. A 29, 2110–2117 (2012). [CrossRef]
34. J. Chen, V. Venugopal, and X. Intes, “Monte Carlo based method for fluorescence tomographic imaging with lifetime multiplexing using time gates,” Biomed. Opt. Express 2, 871–886 (2011). [CrossRef] [PubMed]
35. M. A. O’Leary, “Imaging with diffuse photon density waves,” Ph.D. thesis, University of Pennsylvania, Philadelphia (1996).
36. L. V. Wang and H.-I. Wu, Biomedical Optics: Principles and Imaging (John Wiley & Sons, 2012).
37. D. A. Boas, J. Culver, J. Stott, and A. Dunn, “Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head,” Opt. Express 10, 159–170 (2002). [CrossRef] [PubMed]
38. F. Martelli, A. Sassaroli, G. Zaccanti, and Y. Yamada, “Properties of the light emerging from a diffusive medium: angular dependence and flux at the external boundary,” Phys. Med. & Biol. 44, 1257 (1999). [CrossRef]
39. J. Ripoll, M. Nieto-Vesperinas, S. R. Arridge, and H. Dehghani, “Boundary conditions for light propagation in diffusive media with nonscattering regions,” J. Opt. Soc. Am. A 17, 1671–1681 (2000). [CrossRef]
40. P. K. Yalavarthy, D. R. Lynch, B. W. Pogue, H. Dehghani, and K. D. Paulsen, “Implementation of a computationally efficient least-squares algorithm for highly under-determined three-dimensional diffuse optical tomography problems,” Med. Phys. 35, 1682–1697 (2008). [CrossRef] [PubMed]
41. M. I. Ochoa, Q. Pian, and X. Intes, “Comparison of Compressive Basis for Quantitative Single-Pixel Fluorescence Lifetime Imaging,” in “Optical Tomography and Spectroscopy,” (Optical Society of America, 2018), pp. OTu4D–4.
43. V. Venugopal and X. Intes, “Adaptive wide-field optical tomography,” J. of Biomed. Optics 18, 036006 (2013). [CrossRef]
44. L. Yu, F. Nina-Paravecino, D. Kaeli, and Q. Fang, “Scalable and massively parallel Monte Carlo photon transport simulations for heterogeneous computing platforms,” J. of Biomed. Optics 23, 010504 (2018). [CrossRef]
46. M. Makitalo and A. Foi, “Optimal inversion of the Anscombe transformation in low-count Poisson image denoising,” IEEE Trans. on Image Processing 20, 99–109 (2011). [CrossRef]
47. V. Venugopal, J. Chen, and X. Intes, “Robust imaging strategies in time-resolved optical tomography,” in “Optical Tomography and Spectroscopy of Tissue X,” (Int. Society for Optics and Photonics, 2013), p. 857827.
48. Q. Pian, R. Yao, and X. Intes, “Hyperspectral Single-Pixel Wide-Field Time Domain Diffuse Optical Tomography,” in “Bio-Optics: Design and Application,” (Optical Society of America, 2015), pp. BM2A–6.
49. Q. Fang and D. A. Boas, “Tetrahedral mesh generation from volumetric binary and grayscale images,” in “IEEE Int. Symp. on Biomed. Imaging: From Nano to Macro, 2009. ISBI’09,” (IEEE, 2009), pp. 1142–1145.
51. P. Sonneveld, “CGS, a fast Lanczos-type solver for nonsymmetric linear systems,” SIAM J. Sci. Statist. Comput. 10, 36–52 (1989). [CrossRef]
52. A. Corlu, R. Choe, T. Durduran, M. A. Rosen, M. Schweiger, S. R. Arridge, M. D. Schnall, and A. G. Yodh, “Three-dimensional in vivo fluorescence diffuse optical tomography of breast cancer in humans,” Opt. Express 15, 6696–6716 (2007). [CrossRef] [PubMed]