2.1.

Implementation of Phase Variations

The ultrasonic modulation of light alters the phase, $\varphi$ , of the packet as it propagates through the turbid medium (it is recognized that alternative UL mechanisms may lead to modulation of the absorption coefficient of the medium and, hence, the weight of the packet, though such effects are not considered here). Leutz and Maret first described the phase variation ${\varphi }_{d,j}\left(t\right)$ caused by changes in the optical path length due to ultrasound-modulated particle displacement after the $\mathrm{j}\text{th}$ scattering event.¹¹ Wang incorporated this phase variation in his MC simulations and also introduced a new mechanism of phase variation, ${\varphi }_{n,j}\left(t\right)$ , which considers the ultrasound-induced change in refractive index over one free path⁷ between scattering events. These two mechanisms have been implemented in our MC simulations; the full definition of ${\varphi }_{d,j}\left(t\right)$ and ${\varphi }_{n,j}\left(t\right)$ can be found in the literature.⁷

2.2.

Implementation of the Autocorrelation Function

The autocorrelation function is often used in the analysis of optical phases where a coherent light source is used in a turbid medium (e.g., diffuse correlation spectroscopy^{21, 22} and, indeed, UL¹¹). The reason is that optical phases are often randomized in turbid media and the analysis of their statistical properties using second-order statistics, such as the autocorrelation function, is a convenient metric.

In UL, it has been suggested that the autocorrelation function for path length $s$ can be formed by averaging and summing the phase variations $\Delta {\varphi }_{n,j}(t+\tau )={\varphi }_{n,j}(t+\tau )-{\varphi }_{n,j}\left(t\right)$ over $N$ free paths, and $\Delta {\varphi }_{d,j}(t+\tau )={\varphi }_{d,j}(t+\tau )-{\varphi }_{d,j}\left(t\right)$ over $N-1$ scattering events along one whole photon path from the optical source to the detector,

Assuming a sinusoidal ultrasonic excitation, the autocorrelation function $G_1\left(\tau \right)$ is a cosine waveform from which the maximum variation is given by $1-G_1\left(0.5T_a\right)$ , where $T_a$ is the acoustic period. The maximum variation of $G_1\left(\tau \right)$ indicates the strength of the UOT signal,¹⁵ and the results in Sec. 3 will be presented in terms of this parameter.

Following the same approach as Wang,⁷ the MC simulations were performed on an infinite slab with light illumination on one side and a single-point detector on the other. To speed up the computation, the reciprocity of photon migration was exploited by simulating the injection of photons into the medium at a single point on one side, detecting the exiting photons at all locations on the other.

The parallel architecture of the GPU allows thousands of photon-packet propagations (threads) to be simulated at the same time. The maximum number of simultaneous threads is dependent on how many memory registers are required by the photon propagation algorithm. In the original CUDAMCML code, 26 registers are used in a $32\text{-bit}$ Windows XP compilation,¹⁸ which allows a total of 16,128 simultaneous photon-packet propagation threads. Because of the incorporation of ultrasound modulation, our GPU code is more complicated and requires 37 registers; this still allows 10,752 simultaneous photon packet propagation threads. Figure 1 depicts the basic structure of our GPU code. Although the core of the MC code is similar to that of CUDAMCML, our GPU code also involves the calculation and accumulation of phase variations $\Delta \varphi (t+\tau )=\sum \Delta {\varphi }_{d,j}(t+\tau )+\sum \Delta {\varphi }_{n,j}(t+\tau )$ as photon packets propagate. As a photon packet reaches the detection surface, the final photon weight $W$ and the real and imaginary parts of the phase variation (i.e., $\cos \left[\Delta \varphi (t+\tau )\right]$ and $\sin \left[\Delta \varphi (t+\tau )\right]$ ) are added to the global device memory, which is used to calculate the autocorrelation function, as described in Eq. 2, on completion of the photon propagation simulation.

Fig. 1

More than 10,000 photon packet propagations are simulated at the same time. The phase variations are calculated and accumulated as photon packets propagate. The final photon weight $W$ and the real and imaginary parts of the phase variation, i.e., $\cos \left[\Delta \varphi (t+\tau )\right]$ and $\sin \left[\Delta \varphi (t+\tau )\right]$ are added to the global device memory, which is used to calculate the autocorrelation function at the end of the simulation.

2.3.

Implementation of the Ultrasound Modulated Speckle Patterns

In practice, the UL signals collectively form a speckle pattern. By exporting the MC simulation results in an appropriate manner, an ultrasound-modulated speckle pattern can be generated in postprocessing. In this scenario, coherent light is applied at a single point on one side of the slab, and the ultrasound-modulated speckle pattern is modeled based on the light exiting on the other side. Despite the differing motivations, this is essentially the same process as discussed. The electric field of the $m\text{th}$ transmitted photon packet is

3.1.

Validation of the MC Simulations

Sakadzic and Wang provide an analytical solution for the autocorrelation function¹⁵ for the same scenario as simulated in this paper; this result has been used here to verify our MC simulations. Figure 2 shows the analytical results against our MC results for maximum variations of $G_1\left(\tau \right)$ over a range of scattering coefficients ${\mu }_{\mathrm{s}}$ , ultrasound induced displacement amplitudes $A$ and ultrasound frequencies $f_{\mathrm{a}}$ . The following parameters were used in the simulations: wavelength of light in vacuo ${\lambda }_0=500\mathrm{nm}$ , $n_0=1.33$ , ${\mu }_{\mathrm{a}}=0{\mathrm{cm}}^{-1}$ , anisotropy factor $g=0$ , $\eta =0.3211$ (a function of the adiabatic piezo-optical coefficient of the material, $\partial n∕\partial p$ , the density $\rho$ , and the acoustic velocity of $v_{\mathrm{a}}$ : $\eta =\partial n∕\partial p\cdot \rho \cdot v_{\mathrm{a}}^2$ ),⁷ and thickness of the slab $L=2\mathrm{cm}$ . These values were the same as those used by Sakadzic and Wang.¹⁵ Figure 2 shows agreement between the GPU-based MC results and those produced analytically, thus validating our GPU-based MC code. Similar validation was performed on our CPU-based MC code. It can be seen from Fig. 2 that there is certain discrepancy between the analytical and MC result, which is particularly apparent when $A$ is large (e.g., $A=1.7\mathrm{nm}$ ). Sakadzic and Wang, who derived the analytical solution, noted that an approximation has been made in the derivation of the autocorrelation function and the approximation is only valid when the accumulated phase variation is small $(⪡1)$ .¹⁵ As $A$ becomes larger, the phase variations grow and the approximation in the analytical solution becomes less accurate. The MC code presented here does not use any approximation, and hence, a discrepancy is observed.

Fig. 2

Dependence of the maximum variation, i.e., $1-G_1\left(0.5T_a\right)$ , on the scattering coefficient at different values of ultrasound frequency and displacement amplitude. This figure serves as a validation of our GPU MC results and presents the same information as in Fig. 3(d) of Sakadzic and Wang’s study.¹⁵ Their analytical model was used to calculate the solid lines and the symbols were results obtained by our GPU MC simulations. The following parameters were used in the simulations: ${\lambda }_0=500\mathrm{nm}$ , $n_0=1.33$ , ${\mu }_a=0{\mathrm{cm}}^{-1}$ , $g=0$ , $\eta =0.3211$ , and $L=2\mathrm{cm}$ .

3.2.

GPU Speedup

Two CPUs and two GPUs were each tested with the corresponding versions of the MC codes. The two CPUs were an Intel Core Quad, Q6700, $2.66\mathrm{GHz}$ (slower) and an Intel Core i7-920 (faster). The two GPUs used were an Nvidia GeForce 9800 GT with 112 cores (slower) and an Nvidia GTS250 with 128 cores (faster). The optical and acoustic parameters used were; ${\lambda }_0=500\mathrm{nm}$ , $n_0=1.33$ , ${\mu }_{\mathrm{a}}=0{\mathrm{cm}}^{-1}$ , $g=0$ , $\eta =0.3211$ , $L=2\mathrm{cm}$ , $f_{\mathrm{a}}=1\mathrm{MHz}$ , and $A=1.7\mathrm{nm}$ . All other optical and acoustic parameters were the same as those used in a previous study.¹⁵ The scattering coefficient ${\mu }_{\mathrm{s}}$ was varied over 2, 5, 10, 20, 50, 100, and $200{\mathrm{cm}}^{-1}$ . All computation times were taken as a ratio to that of the slower CPU (Core Quad). In comparison to the reference CPU (Core Quad), the faster CPU (Intel Core i7) offered a computational speedup of around 1.57 for all values of ${\mu }_{\mathrm{s}}$ . Figure 3 shows the ratio of the CPU (Core Quad) computation time to the GPU computation times for the same set of parameters. The speedup varied with the choice of ${\mu }_{\mathrm{s}}$ . The GPU version is faster than its CPU counterpart by up to a factor of 125 with ${\mu }_{\mathrm{s}}=20{\mathrm{cm}}^{-1}$ . As expected, the GTS 250 is, in general, faster than the GeForce 9800 especially for increasing ${\mu }_{\mathrm{s}}$ .

Fig. 3

GPU speedup, defined as the ratio of CPU computation time to GPU computation time, as a function of ${\mu }_s$ . As an example, the GPU computation time is faster than its CPU counterpart by a factor 125 for ${\mu }_s=20{\mathrm{cm}}^{-1}$ . The following parameters were used in the simulations: ${\lambda }_0=500\mathrm{nm}$ , $n_0=1.33$ , ${\mu }_a=0{\mathrm{cm}}^{-1}$ , $g=0$ , $\eta =0.3211$ , $L=2\mathrm{cm}$ , $f_a=1\mathrm{MHz}$ , and $A=1.7\mathrm{nm}$ .

3.3.

Additional Computation Time Due to Incorporation of Ultrasound

Because of the additional complexity introduced by the incorporation of ultrasound propagation the modeling of UL imposes an increased computational load in comparison to the modeling of conventional photon migration. Figure 4 shows the ratio of the GPU computation time with ultrasound to that without ultrasound as a function of ${\mu }_{\mathrm{s}}$ . All other optical and acoustic parameters used were the same as those used in Sec. 3.2. The MC code CUDAMCML developed by Alerstam ¹⁸ was used to obtain the computation times for conventional photon migration modeling (without ultrasound). The GPU used was the Nvidia GeForce 9800. It can be seen from Fig. 4 that the additional computation effort demanded by the incorporation of ultrasound increased the computation time by a factor of approximately 6 to 17, depending on the value of ${\mu }_{\mathrm{s}}$ .

Fig. 4

Additional computation ratio, defined as a ratio of the GPU computation time with ultrasound (US) to that without US, as a function of ${\mu }_{\mathrm{s}}$ . As an example, the incoporation of ultrasound in the MC simulations increased the computation time by a factor of $>6$ for ${\mu }_{\mathrm{s}}=10{\mathrm{cm}}^{-1}$ . The GPU used was the Nvidia GeForce 9800.

3.4.

Ultrasound-Modulated Speckle Patterns

Figure 5 is the video of the ultrasound-moudlated speckle pattern generated by postprocessing the MC simulation results, and it has a total of 12 frames played at $4\mathrm{fps}$ . Figure 6 depicts the ultrasound-modulated speckle pattern at six points in time over one acoustic cycle $\left(1\mu \mathrm{s}\right)$ . The two boxes inset in Figs. 6, 6, 6, 6, 6, 6 highlight two speckle grains that have different phases, i.e., the brightest moment of the speckle grain highlighted by the box on the left happened at (f) $0.83\mu \mathrm{s}$ whereas that highlighted by the box on the right happened at (c) $0.332\mu \mathrm{s}$ . The simulations were completed in $4\mathrm{s}$ using the GPU code compared to $290\mathrm{s}$ for the CPU version—a speedup of 72 times.

Fig. 5

Video 1 showing ultrasound-modulated speckles: speckle grains blink at the same frequency but with differing phases (QuickTime, $1.9\mathrm{MB}$ ). .

Fig. 6

Ultrasound-modulated speckle patterns varied over one acoustic cycle $\left(1\mu \mathrm{s}\right)$ computed by a GPU (Nvidia GeForce 9800). The two boxes in (a–f) highlight two speckle grains, respectively, that have different phases, e.g., the brightest moment of the speckle grain highlighted by the box on the left happened at (f) $0.83\mu \mathrm{s}$ , whereas that highlighted by the box on the right happened at (c) $0.332\mu \mathrm{s}$ . Computation $\text{time}=4\mathrm{s}$ with a GPU compared to $190\mathrm{s}$ with a CPU.

4.1.

Simulation Efficacy

The curves of Figs. 3 and 4 both demonstrate two distinct trends. Below a particular value of ${\mu }_{\mathrm{s}}$ , (e.g., ${\mu }_{\mathrm{s}}=20{\mathrm{cm}}^{-1}$ , in Fig. 3, and ${\mu }_{\mathrm{s}}=10{\mathrm{cm}}^{-1}$ , in Fig. 4) the efficacy of the GPU-based UL simulation increases sharply with ${\mu }_{\mathrm{s}}$ . Above this threshold, its efficacy gradually reduces. These trends are consequences of our GPU implementation rather than the underlying algorithm. Two aspects dominate the observed efficacy as discussed in the following. As ${\mu }_{\mathrm{s}}$ is increased, the first aspect of the implementation leads to an increase in efficacy while the second aspect a decrease. The overall change in efficacy depends on which of these two aspects dominate. The two apects are as follows:

4.2.

New Possibilities

The incorporation of ultrasound propagation in photon migration modeling increases the computational effort considerably (e.g., the computation time increased by a factor of at least 6, as shown in Sec. 3.3). The improvement in speed is therefore an important advancement in the modelling of UL. Using such techniques additional acousto-optic mechanisms can be simulated before the simulation speed becomes a limiting factor in such research. Furthermore, numerical approximations (which can be detrimental to both accuracy and the readability of the simulation code) may no longer be required to ensure tolerable computation time. Previously, the simulated geometry, medium, type of acoustic source, and locations of the optical source and detector were often simplified partly due to the issue of computation time. For example, one early scheme (as discussed in Sec. 2.2) required the optical source to be (theoretically) infinitely large and for the detector to be a single point on the other side of the slab. The reciprocity of photon migration was relied on such that the code was sufficiently efficient. With a general speedup in the code, a more realistic scenario can be modeled without incurring an unacceptable increase in computational time.

4.3.

Image Reconstruction

An alternative to MC modeling is the finite element method (FEM), which has been widely adopted in diffuse optical imaging²⁷ and has a relatively fast computational speed. The computational speed of forward modeling is especially important in image reconstruction because of the iterative nature of an inverse problem. Sakadzic and Wang have developed a correlation transfer equation for UL¹⁷ that can potentially be solved using FEM and used for image reconstruction. However, unlike solving the diffusion approximation, as in the case of diffuse optical imaging, solving the correlation transfer equation is not an easy task and may involve some difficult numerical issues. Moreover, the existing correlation transfer equation has been developed for specific acousto-optic mechanisms. This means that the incorporation of an additional mechanism will require a new form of the correlation transfer function to be derived. In comparison, the MC approach is more straightforward and flexible, allowing the easy incorporation of new mechanisms (e.g., the nonphase mechanism¹³). With the improvement of computational speed using a GPU, forward modeling using MC simulations is now a viable option for UOT image reconstruction. Although a GPU implementation of a carefully designed FEM scheme could potentially outperform its MC counterpart, the absolute time required for the MC implementation would nonetheless remain reasonable [e.g., $10\min$ with a GPU as compared to $21\mathrm{h}$ with a CPU (assuming a speedup of 125 times)].

4.4.

Ultrasound-Modulated Speckle Pattern

Figures 5 and 6 demonstrate an important characteristic of ultrasound-modulated speckles whereby individual speckle grains are modulated at the ultrasound frequency but with individual absolute phase (i.e., speckle grains blink at the ultrasound frequency but reach the brightest point at different times). In practice, the speckle grains will also be randomly modulated due to Brownian motion and this behavior has often been exploited to estimate particle size within the medium. Brownian motion was not modeled in our MC simulations such that the dynamics of the speckle grains are governed only by the applied ultrasound. We are currently investigating the incorporation of Brownian motion as another mechanism that may affect the UOT signal. One important use of the speckle pattern modeling described here is to help the design of an efficient detection scheme. Speckle patterns can be simulated using different combinations of acoustic and optical parameters, and the results used to test and compare the performance of different detection schemes in a standardized manner.