Monte Carlo Methods for Simulation of Optical Coherence Tomography of Turbid Media

We describe two importance sampling techniques for a standard Monte Carlo (MC) method that could enable fast simulation of signals from optical coherence tomography (OCT) imaging systems. These OCT signals are generated due to diffusive reflections from either multilayered or arbitrary shaped, turbid media, for example, tissue. Such signals typically consist of ballistic and quasi-ballistic components, of scattered photons inside the medium, in addition to photons that undergo multiple scattering. We show that MC simulation of these OCT signals using importance sampling reduces its computation time on a serial processor by up to three orders of magnitude compared to its corresponding standard implementation. Therefore, these importance sampling techniques enable practical simulation of OCT B-scans of turbid media, for example, tissue, using commonly available workstations.


Introduction
Optical coherence tomography (OCT) is a non-invasive sub-surface imaging technique that has experienced significant growth in biomedical applications [1,2]. OCT systems could be implemented with a low-coherence light source and a mechanical scanning sub-system (time-domain OCT). More advanced systems use a low-coherence light source with a spectrometer or a wavelength-swept source (frequency-domain OCT). OCT has an imaging depth that could reach up to 3 mm, depending on the optical properties of the tissue, and it also has one to two orders of magnitude higher resolution than ultrasound imaging. OCT could also produce images inside the body when it is implemented using optical fiber probes. Unlike X-ray or gamma-ray imaging, OCT is safe for biological tissues because it utilizes non-ionizing radiation mainly in the infrared spectrum.

OCT signal simulation using a Monte Carlo method
The signal obtained by an OCT imaging system consists of ballistic and quasi-ballistic (Class I diffuse reflectance) photons, in addition to multiply scattered photons (Class II diffuse reflectance), that are reflected from tissue [3].
However, multiply scattered photons do not carry practically useful information about the imaged tissue; therefore, they result in a degradation of the OCT signal [4]. In addition, it has been shown that Class II diffuse reflectance represents a fundamental limit related to the imaging depth of OCT in tissue [5]. Therefore, it is important to account for both Class I and Class II diffusive reflectance in any practical simulation of OCT signals.
Since it is not practical to simulate light transport in turbid media, for example, tissue, using electromagnetic waves, especially due to diffusive scattering, a Monte Carlo (MC) method of simulating light transport in tissue has been typically used [6][7][8][9]. However, the computational cost of this MC-based simulation of OCT systems could be very high, as the probability of detecting diffusively reflected photons from tissue is very low [4,5].
To reduce the computational cost, thereby accelerate, this MC simulation, importance sampling could be used to speed up simulations by orders of magnitude. Importance sampling has been applied earlier to optical communications [10,11], confocal microscopy [12], atmospheric optics [13], and diffuse optical tomography [8].
To improve the computational efficiency of the MC-based simulation of OCT systems [6], Yao and Wang proposed the first importance sampling technique to simulate OCT signals from a multilayered turbid medium [3]. However, their method only enabled the simulation of OCT signals from a thin shallow layer in tissue, as the results obtained from deeper tissue regions were underestimated. In [14], we, the authors of this chapter, developed another more advanced importance sampling technique by implementing multiple biased scatterings per photon packet, and by developing a photon splitting procedure. Our advanced importance sampling resulted in more accurate and computationally efficient simulations of OCT signals due to ballistic and quasi-ballistic photons. However, it still underestimated OCT signals due to multiply scattered photons. To enable accurate simulation of OCT signals due to both Class I and Class II diffusive reflectances, we further developed our importance sampling technique to accurately and efficiently simulate diffusive reflectance due to photons that undergo multiple scattering [15]. In this method, additional biased scatterings were randomly applied, which enabled accurate simulation of both Class I and Class II diffusive reflectances, with a speed-up of three orders of magnitude compared to the standard MC method.
Our advanced importance sampling techniques above were implemented to simulate OCT of tissues with planar geometry [6]. To enable simulation of OCT of arbitrarily shaped turbid media, we used the mesh-based MC method of light transport in tissue proposed by Fang [16]. This method uses a Plücker coordinate system to efficiently calculate intersections between paths of light propagation with interfaces of the object regions that are modeled using tetrahedrons. We combined this mesh-based MC method with our importance sampling techniques to simulate OCT signals from tissue with arbitrarily shaped regions. However, since it was still computationally costly to simulate a full OCT B-scan using this method, we also developed a parallel implementation of this simulator that exploited the massively parallel computing capabilities of Graphics Processing Units (GPUs) to accelerate this simulator by two additional orders of magnitude [17,18]. This GPU-based implementation enabled simulation of OCT B-scans of arbitrarily shaped turbid media in a few minutes using commonly available workstations.
In Section 2 of this chapter, we present a standard MC method for simulating OCT signals. In Section 3, we present our first importance sampling implementation that enables the simulation of OCT signals from higher depths inside turbid media. In Section 4, we present our more advanced importance sampling implementation that accurately calculates both Class I and Class II diffusive reflectances, and is three orders of magnitude faster than the standard MC simulation.

Standard MC method for simulating OCT signals
Our implementation of the MC method to simulate OCT signals is based on Monte Carlo simulation of light transport in multilayered tissues (MCML) [6]. MCML simulates an ensemble of photon packets that are launched as a steady-state pencil beam, normal to the top surface of the medium. Within the tissue, each such photon packet undergoes a random walk whose step size is determined by an exponentially distributed random variable that is parameterized by the interaction coefficient of this tissue. This interaction coefficient is equal to the sum of the absorption μ a and the scattering μ s coefficients of this tissue. The scattering events that take place at the end of the random steps are characterized by two random angles that determine the next direction of the photon packet. To account for the photon packet scattering, given an anisotropy factor, g, of the tissue, MCML uses the Henyey-Greenstein probability density function that is given by where θ S is the angle between the propagation direction of the photon packetû before the current scattering andû 0 is the direction of the photon packet after the current scattering. The angle between the previous propagation direction and the new propagation direction is θ S . Therefore, cos(θ S ) =û Áû 0 . To ensure that the new propagation direction is statistically correct, the provisional scattering directionû 0 is rotated aroundû by an angle ϕ, which is randomly selected from a uniform probability density function with a range from 0 to 2π to generate the propagation directionû 0 after the current scattering. At the end of each scattering event, the photon packet weight W is reduced according to the step size and the local absorption coefficient μ a . The weight W, which is initially set at 1, is proportional to the number of photons in the photon packet. The photon packet is either removed with probability 1/m or is allowed to continue propagating with probability 1À1/m and a new weight equal to mÁW once the weight reaches W th ¼ 10 À4 . We use the value m = 10 in this work. This procedure, denoted Russian roulette, is an unbiased technique to end simulation of photon packets that have a negligibly low contribution to the Monte Carlo simulation, so that a new photon packet can be initiated and simulated.
The Class I diffuse reflectance at depth equal to z is obtained by calculating the mean value of an indicator function I 1 that represents a spatial and temporal filter of Class I diffuse reflectance for all simulated photon packets. The indicator function I 1 of such spatial and temporal filter for the ith photon packet is defined as where l c is the optical source's coherence length; r i is the distance to the origin of the ith reflected photon packet; d max and θ max are the maximum photon packet collecting diameter and angle, respectively; θ z,i is the angle with the z-axis (normal to the tissue interface); Δs i is the optical path; and z is the photon packet's maximum depth.
At any depth, the diffuse reflectance R 1 is the expected value of I 1 at this depth, and its standard deviation σ R,1 could be estimated by and where N is the number of photons packets used in MC-based simulations.

Importance sampling for simulation of Class I OCT signal
Our first importance sampling technique to simulate OCT signals aimed at increasing the number of photons collected at the detector. This algorithm uses the same method described in Section 2, where we also use the same square time gating given by [3].
Since most tissues are highly forward-scattering, their anisotropy factor is close to 1. Therefore, there is a very small probability that a simulated photon packet at any given depth in the tissue would undergo scattering in the backward direction toward the OCT probe. The probability of collecting Class I photons drops rapidly with depth in the tissue from which the photon is scattered in the backward direction. To allow faster simulation of Class I photons, we designed an importance sampling method that biases the directionû 0 of a scattered photon toward the tip of the light-collecting optical fiber,v, as the photon packet reaches the depth range of interest. By defining the origin of a Cartesian coordinate system at the center of the tip of the light-collecting optical fiber, the bias direction in which this fiber is located is defined asv ¼ ÀR= R j j, where R ¼ xx þ yŷ þ zẑ is the position vector of the scattering location in the tissue.
All photon packets propagating in a direction close tov will contribute to the simulated Class I diffuse reflectance with a higher probability. Therefore, this bias direction is more efficient than biasing only in the backward direction, which may not be consistent with the direction of the light-collecting optical fiber. This choice of the bias direction is particularly effective for photon packets propagating deep in the tissue, where such photon packets experience one or more scattering events before they are diffusively reflected.

Scattering angle due to first event of backscattering
As the photon packet reaches the depth range targeted, the propagation directionû 0 of the scattered photon packet is biased toward the bias directionv, as opposed to being most likely scattered close to the previous propagation directionû as in the practical case with anisotropy g close to 1 and different from the bias toward Àû, the opposite of the direction of propagation, as it is done in [3]. To randomly select the biased angle θ B between the new scattering directionû 0 and the biased directionv, we use the same probability density function in Eq. (1). However, the bias coefficient does not have to correspond to the anisotropy factor g. Therefore, the probability density function of the biased angle is given by where a is a bias coefficient. After randomly sampling a biased angle θ B away from the biased directionv, so that cos(θ B ) =v Áû 0 , the provisional scattering directionû 0 is rotated aroundv by an angle ϕ, which is randomly selected from a uniform probability density function with a range from 0 to 2π to generate the propagation directionû 0 after the current scattering. This procedure ensures a more accurate model of the light scattering in tissue. The difference in the rotation by ϕ between the model with importance sampling and the standard model is that the rotation in the model with importance sampling is done around the biased direction v, while the rotation in the standard model in is done around the directionû prior to the current scattering. After the first biased scattering, this procedure produces the new propagation directionû 0 of the photon packet. Afterward, the scattered photon packet is associated with a likelihood ratio as discussed in other applications of this method [10,11]. Using our biased angle's probability density function, the likelihood ratio of the photon packet, Eq. (5), is given by where cos(θ S ) =û Áû 0 is determined, after the biased scattering, from the randomly sampled values of θ B and ϕ. The ratio of the probability of the scattering angle appearing in the biased case with the standard case is the likelihood ratio that is shown in Eq. (6). In addition to depending on θ B , the likelihood ratio also depends also on θ S . Figure 1 shows a schematic drawing of these vectors and the angles used in this direction biasing procedure. Note that the choice of bias distribution only affects the speed of convergence of the simulation. Therefore, other biased probability function could also be used to randomly generate the biased scattering toward the bias directionv.

Scattering angles of additional events of backscattering
As a photon packet is biased toward the apparent position of the collecting optical fiber, at any given depth in the tissue, the photon packet becomes more likely to be collected at the tip of the fiber. However, the photon packet could be scattered several times after the first backscatter bias before reaching the optical collection system. These additional scatterings, according to Eq. (1), reduce the correlation between the biased direction and the event in which the photon packet is collected. We overcome this reduction in correlation by continuing to bias the scattering directionû 0 toward the directionv, pointing to the apparent position of the optical collection system, at every scattering point until the photon packet is removed. These additional biases still use both Eqs. (7) and (8). Since the random values drawn for the angle between the scattering direction and the biased direction are independent of each other and are also independent of the previous scattering events, the overall likelihood ratio of a collected photon packet results from the multiplication of all the likelihood ratios of all the biased scattering in that particular simulation.
Once a photon packet experiences the first biased scattering, that photon packet is biased at all additional scattering points until it is removed from the simulation, which can occur when the photon packet is removed by Russian roulette, as described in Section 2, or it leaves the tissue. After simulating N launched photon packets using importance sampling, the diffuse reflectance R 1 and its standard deviation σ R,1 could be calculated with and Eqs. (7) and (8) are similar to Eqs. (3) and (4), except that the indicator function is multiplied by its corresponding likelihood ratio. Using this method, a significantly larger number of photon packets are scattered from a specific depth range toward the collecting optical system than the number obtained using a standard MCML implementation. At the end of this biased simulation, each photon packet is weighted by its likelihood ratio, which adjusts the contribution of each packet to the estimation of the Class I diffuse reflectance. The estimated diffuse reflectance converges toward its true value faster, by several orders of magnitude, when compared to the standard Monte Carlo method.

Importance sampling effectiveness and depth of tissue
One drawback of previously existing bias methods, for example, [3,7,8] is an underestimation of the diffuse reflectance beyond the targeted depth range. The application of the first backward bias reduces the probability that this photon packet would propagate beyond that portion of the tissue. This would lead to a statistical bias to this importance sampling method similar to that in the angle biasing procedure used in [3] and the method used in [7,8], which limits the effectiveness of those methods to a thin target layer.
We make sure we obtain correct statistics by splitting the photon packet into two photon packets before the first biased scattering [14]. The first of these two photon packets is the one biased toward the collecting optical system. The second photon packet starts propagating from the location in which the biased backscattering occurred, where its initial direction calculated by applying the standard procedure to the previous directionû as shown in Section 3.1. To ensure that there is no statistical bias associated with the forward-propagating photon packet that was split, it will be assigned a likelihood ratio L 0 i ð Þ, which is a complement to the likelihood of the biased backward scattered photon packet L i ð Þ such that L 0 i ð Þ ¼ 1 À L i ð Þ to this second photon packet. This second photon packet, only generated if L i ð Þ < 1, also undergoes biased backscattering in the tissue at the end of the next step, which could result in another photon packet split, and successive additional biased scatterings toward the tip of the collecting optical fiber until the photon packet propagates beyond the simulation domain. In cases that we investigated, this procedure increased the computational time of each photon packet by five times when compared with a photon packet computed using the standard Monte Carlo method with the same number of launched photon packets N. The increase in the computational time of our importance sampling-based implementation, compared to the standard method with the same number of launched photon packets, depends on the average value of the mean free path, and on the width of the target depth range. We note that in our importance sampling implementation, we do not count split photon packets as additional photon packets when determining the value of N in Eqs. (7) and (8), as the use of their corresponding likelihood ratios will generate the correct result. As a photon packet propagates beyond the target region, the packet will propagate using the standard scattering procedure until it is terminated. Once this photon packet is terminated, a new photon packet will be simulated from the OCT probe, as it is the case in the standard MCML. Even though the splitting procedure implies that the cost of simulating a launched photon with this importance sampling method is higher than the computational cost of simulating launched photos using the standard MCML, the computation cost of the Class I diffuse reflectance in our Monte Carlo simulations with importance sampling required as little as one-thousandth of the computational cost required by the standard Monte Carlo method to achieve the same accuracy in the calculated diffuse reflectance.

Numerical results
We validate our importance sampling technique for simulation of OCT signals from multilayered tissue, with different refractive indices and scattering properties, by comparing its results with those obtained by the standard Monte Carlo method. As shown in Figure 2, light is emitted by an optical fiber probe that is reflected by a prism. The shown optical system has a focusing lens with a numerical aperture (NA) that allows collecting light at an angle of up to 4°and a diameter of 0.5 mm. Similar to the setups in [5,10,11], we assume a point source that emits in the vertical direction. Air is present between the center of the probe and the first layer of tissue, which is placed 2.12 mm from the center of the fiber. We simulate a three-layer turbid medium with refractive-index mismatch at its interfaces. The first layer, extending from 2.12 to 2.22 mm from the tip of the fiber, has absorption coefficient μ a = 1.5 cm À1 , scattering coefficient μ s = 60 cm À1 , and refractive index n = 1. The second layer, extending from 2.32 to 2.42 mm from the tip of the fiber, has the same absorption and scattering coefficients as the first layer, but its refractive index is n = 1.33. The third layer, extending from 2.42 to 2.62 mm from the tip of the fiber, has the following parameters: μ a = 1.5 cm À1 , μ s = 30 cm À1 , and n = 1. After the third layer, the medium was assumed to be air: μ a = 0 cm À1 , μ s = 0 cm À1 , and n = 1. The anisotropy factor was assumed g = 0.9 for the three diffusive layers.
From Figure 3, we note an excellent correspondence between results obtained with our new importance sampling method and results obtained using MCML, that is, standard Monte Carlo simulations. However, our results were obtained in one-thousandth of the time required by the standard method.

Importance sampling for simulation of Class I and Class II OCT signals
In this section, we further improve the importance sampling technique that was described in Section 3, so we can simulate Class II OCT signals more accurately and more efficiently [19].

Scattering angle of first backscattering event
In the MC simulation described in Section 3, we note that the bias function in (5) produces large values of the likelihood ratio (>>1) when photon packets are scattered in the then unlikely forward direction. These photon packets contributed to a slow decrease in the relative variation, which corresponds to relative error, with the increase in the number of photon packets launched for the calculation of the OCT signal. Referring to Figure 1, we could reduce this relative variation by choosing a distribution function for the scattering angle that limits it to the backward direction. This modified distribution is given by: Class I diffuse reflectance dependence on the distance from the center of the optical fiber for the simulation whose schematic is shown in Figure 4. The solid black line represents the result obtained with 2 Â 10 5 photon packets using the importance sampling technique presented in Section 3. The green long dashed line is the result obtained with 10 9 photon packets using standard MCML. The blue dots represent results obtained with 10 6 photon packets using MCML. The pink short dashed lines show the confidence interval of the importance sampling simulations with 2 Â 10 5 photon packets that were estimated using a much larger ensemble of 64 Â 10 5 simulations. Reprinted with permission from [14] © The Optical Society of America.
where a is the bias coefficient that can be selected between 0 and 1. Once a biased angle θ B is randomly selected, away from the direction of the center of the OCT probev, where cos θ B ð Þ ¼v Áû 0 , the provisional biased scattering directionû 0 is rotated aroundv by an angle ϕ randomly sampled from a uniform distribution in the range from 0 to 2π. These parameters are defined in the same manner as those used in the biased distribution presented in Section 3. The only difference is that the domain of cos θ B ð Þ is restricted to a maximum deviation from the biased angle by 90°. This ensures that there would not be packets with very large likelihood ratio that could reduce the efficiency of our importance sampling. The likelihood ratio of the photon packet that uses the biased probability density function in Eq. (9) is given by where cos θ S ð Þ ¼û Áû 0 . We note that cos θ B ð Þ is obtained using the probability density function in Eq. (9), where it is used to obtain the new propagation direction u 0 .
To sample angles according to the biased probability density function in (8), one could use any uniform pseudo-random number generator that would be typically available in scientific software libraries. For example, if u i is a random number distributed uniformly between 0 and 1, a random value for cos θ B ð Þ that satisfies Eq. (9) with bias coefficient a could be generated with the following conversion formula This conversion formula was derived using probability theory [20].

Scattering angles of additional biased backscatterings
A second enhancement that could be made to the importance sampling technique, described in Section 3, is to bias the additional scatterings toward the center of the OCT probev with probability 0 ≤ p ≤ 1. That contrasts with the technique presented in Section 3, in which p = 1 (all the additional scatterings were biased). An unbiased scattering is performed in case a bias scattering is not applied in a given point where scattering takes place. The likelihood ratio associated with this scattering is calculated according to the formula If the biased function f B cos θ B ð Þ ð Þis selected to sample a random value of cos θ S ð Þ, which is an event with probability p, cos θ S ð Þ ¼û Áû 0 is a function of cos θ B ð Þ that is statistically sampled from the probability density function in Eq. (9).
Otherwise, in the case of the complementary event with probability 1 À p, the unbiased function f HG cos θ S ð Þ ð Þis used to sample a random value of cos θ S ð Þ and cos θ B ð Þ ¼v Áû 0 depends on the value of cos θ S ð Þ. Since the two random angles associated to each scattering do not depend on the random angles selected in the previous scatterings, the likelihood ratio of each collected photon packet results from the multiplication of all the likelihood ratios of all the biased scatterings in that simulated photon packet.

Numerical results
We validate our importance sampling technique for simulation of OCT signals from multilayered tissue, with different refractive indices and scattering properties, by comparing its results with those obtained by the standard Monte Carlo method. We consider a tissue model that consists of multiple layers that could be imaged with an OCT system, as shown schematically in Figure 4. The modeled tissue extends from 0 to 1 mm, and consists primarily of a turbid layer with absorption coefficient μ a = 1.5 cm À1 and a scattering coefficient μ s = 60 cm À1 , and also contains five thin layers with absorption coefficient μ a = 3 cm À1 and a scattering coefficient μ s = 120 cm À1 . These five thin layers with higher scattering coefficient are located from 200 to 215 μm, from 365 to 395 μm, from 645 to 660 μm, from 760 to 775 μm, and from 900 to 915 μm. We assume that this tissue has the same refractive index n = 1 and an anisotropy factor g = 0.9. We note that our method is robust in the presence of refractive index mismatch along layer boundaries of the tissue [19]. We simulate an OCT system where the light is delivered/collected by the tip of an optical fiber having a radius of 10 μm and an acceptance angle of 5°. For simplicity, the light source is assumed to be a one-dimensional light beam propagating along the vertical direction as in [3,8], since the purpose of this example is to validate and demonstrate the effectiveness of our second importance sampling technique when it is applied to the standard MC simulation.
In Figures 5 and 6, we show results obtained with 10 8 Monte Carlo photon packets with importance sampling, which has a computational cost of simulating about 9 Â 10 8 photon packets using standard Monte Carlo. The computational cost of applying this importance sampling technique depends on the target depth range, and on the average photon mean free path in the given tissue. The target depths in the shown simulations were set from 0 to 1 mm. Therefore, every single photon scattering that occurs in the depth range from 0 to 1 mm would be biased. We used a bias coefficient a = 0.925, and an additional bias probability p = 0.5, to run the Monte Carlo simulations with importance sampling. The results shown in Figure 5. The Class I diffusive reflectance (thick solid black curve) and the Class II reflectance (thin solid red curve), as a function of the tissue depth using the importance sampling technique presented in Section 4 with 10 8 photon packets. The pink short dashed and the blue long dashed curves are results of simulating Class I reflectance and the Class II reflectance using standard Monte Carlo with 10 11 photon packets, respectively. Reprinted with permission from [15] © The Optical Society of America.  Figures 5 and 6 show that our new importance sampling technique reduces the computational cost for obtaining the Class I diffuse reflectance by approximately three orders of magnitude when compared to the standard Monte Carlo method. This algorithm is optimum when the additional bias probability is equal to p = 0.5. Since only half of the back-scatterings are biased, this choice contributes toward enabling an optimum number of Class II photons to be collected by the tip of the optical fiber.
We note that the results obtained with the MCML have confidence intervals that are noticeably larger than those of the results obtained with importance sampling shown in Figure 6, even though the standard Monte Carlo simulations have a computational cost 113 times larger than those obtained with importance sampling. In Figure 6, we also note that our importance sampling technique reduces the computational cost of calculating the Class II reflectance by more than two orders of magnitude.
In Figure 7 we show the relationship between the relative error in calculating Class I and the Class II reflectances at two different depths: 400 and 670 μm and the bias coefficient a for p = 0.5. The depths at 400 and 670 μm correspond to tissue regions near the second and third regions with high local reflectance due to the higher local scattering coefficient. The relative variation in the results is the ratio between the square root of the variance, shown in Eq. (4), and the reflectance in Eq. (3).
We note that Class I reflectance has a minimum relative error at 400 μm with a = 0.925, but the minimum error at 670 μm occurs at a = 0.95 μm at 670 μm. The deeper the tissue region, the stronger the required bias because of the increase in the number of scatterings with the depth. However, as the bias coefficient is increased toward 1, larger variations in the likelihood ratio lead to an increase in the relative error. We also note that Class II reflectance has its minimum relative error at 400 μm with a = 0.91, while its minimum relative error at 670 μm increased to only about a = 0.925 μm. The optimum amount of bias required by the Classs II diffusive reflectance in both wavelengths is lower than the optimum bias coefficient observed in the Class I reflectance because the number of ballistic and quasi-ballistic photons increases with the bias, which leads to a decrease in the number of collected photon packets that undergo multiple scatterings. Figure 7 also shows that there is a range for the bias parameter a between 0.9 and 0.95 that enables fast calculation of both Class I and Class II reflectances using our importance sampling-based implementation.

Conclusions
We described two importance sampling techniques for a standard Monte Carlo (MC) method that could enable fast simulation of signals from optical coherence tomography (OCT) imaging systems. These OCT signals are generated due to diffusive reflections from either multilayered or arbitrarily shaped, turbid media, for example, tissue. Such signals typically consist of ballistic and quasi-ballistic components, of scattered photons inside the medium, in addition to photons that undergo multiple scattering. We showed that MC simulation of these OCT signals using our importance sampling reduced its computation time on a serial processor by up to three orders of magnitude compared to its corresponding standard implementation. Therefore, our importance sampling techniques enable practical simulation of OCT B-scans of turbid media, for example, tissue, using commonly available workstations.