1. Introduction

The performance of optical coherence tomography (OCT)-based angiographic imaging has improved dramatically over the past 5 years [1, 2]. However, confounding artifacts and high noise levels continue to complicate the interpretation and segmentation of OCT angiographic datasets. Therefore, it is necessary to improve the signal and image processing algorithms upon which OCT-based angiography depends. To this end, both numerical and in vivo testbeds are useful. In vivo testbeds provide ultimate confirmation of performance, while numerical testbeds allow a more controlled evaluation of algorithms and optimization of their parameters. Numerical models based on Monte Carlo (MC) approaches have been used in the development of anatomical and polarization-sensitive OCT [3, 4]. However, these MC OCT models are inadequate for simulating microvascular signals for two reasons. First, OCT-based angiography is based on the detection of the time-varying nature of OCT signals, and the general approach of MC does not recapitulate these temporal variations. Second, OCT models are based on layered geometries that are too simplistic to model complex three-dimensional vascular networks.

In this paper, we describe a MC-based stimulation of OCT-based angiographic imaging of arbitrary three-dimensional vascular networks. The model supports multiple tissue types organized spatially into arbitrary three-dimensional geometries. We have also implemented Fourier-domain coherence gating directly into this model. As such, the model recapitulates many of the critical features of OCT-angiographic images including point-spread function broadening and vessel shadows. Most significantly, we have implemented MC-based simulation of OCT signal temporal variations within the regime wherein the flow-based displacement of scatterers is larger than the OCT system resolution (larger of transverse or axial). This allows modeling of angiographic signals. In the first two sections of this report, the methods for performing MC transport within arbitrary three-dimensional geometries (section 2) and the method for simulating both transverse and axial signal localization (section 3) are overviewed. In section 4, we describe the simulation of time-varying signals from the MC model. Finally, in section 5, we validate the approach through comparison of simulated images with datasets acquired in vivo in dorsal skinfold chamber preparations in mice.

2. Monte Carlo modeling of light transport within regions of arbitrary geometry

In prior works, MC models imposed a homogeneous or layered geometry [5, 6]. To support modeling of microvascular networks, it is necessary to extend these methods to support more complex region geometries. Our approach uses the same validated MC engine described in [5] but modifies the methods used to detect the intersection of a photon packet and a region boundary. First, we summarize the MC engine. Photon packets are launched from one surface location orthogonally into the sample. Photon packets begin with an initial weight, W, equal to 1. Once a photon packet is launched, the step length is selected randomly from an exponential distribution, the scale of which is defined by the optical properties of the tissue. The photon packet is then advanced along this direction orthogonal to the first surface until it either hits a boundary between regions, or it reaches the step length. If the photon packet reaches the step length, its weight is decreased by the ratio of the medium absorption to total interaction coefficient. The remaining photon packet is then scattered with a deflection angle, θ, and an azimuthal angle, ψ, which are statistically sampled from the Henyey-Greenstein scattering phase function and between 0 and 2π, respectively. If the photon packet intersects a region boundary along its advancement before reaching the step length, the packet is moved to the boundary. Here it is internally reflected or transmitted according a random number generator with probabilities matched to the Fresnel’s equation (including regimes for total internal reflection). After a transmission or an internal reflection, the photon packet completes the remaining dimensionless step distance divided by the new region absorption and scattering coefficients (if transmitted). A photon packet is discarded if it crosses an outside boundary or by a Russian Roulette random process once weight values fall below a predefined threshold limit. In the method of [5], regions are defined by z-axis locations and it is trivial to identify boundary transections.

To support a more complex geometry, we subdivide the sample into cuboidal voxels. We note that because these cuboidal voxels have boundaries that are always aligned to x/y axes, this method would not accurately estimate the specular (i.e., Fresnel) reflection from a tilted surface. In this approach, each voxel is assigned a specific tissue type with corresponding optical properties. The resolution of the geometry is defined by the voxel size, and can be reduced arbitrarily (at a computational cost) to simulate samples with higher resolution. For microvascular networks, the smallest vessels of which are on the order of 4–9 μm, we adopt a voxel sizes of 3 μm (27 μm³) [7, 8].

A boundary intersection algorithm was implemented to first identify the voxel boundaries that are intersected by an advancing photon packet. If there are no cube boundary intersections or if there are intersections but each intersection occurs between voxels associated with the same region (and therefore identical optical properties), the photon packet is transported the full step length. If there are cube boundary intersections within the photon packet path and one of these intersections occurs between two voxels of different regions, the photon packet is transported to the first intersection between differing voxels where it is internally reflected or transmitted. Analogously to the multilayered approach, after a transmission or an internal reflection, the photon packet continues to travel as described above the step distance adjusted by the new region total interaction coefficient.

To validate that the developed MC model provides similar results as validated layered codes [5], we compared the absorption distribution of light in a three-dimensional medium. To record absorption in the cross-sectional plane x = 0, a k × l grid was defined where each element has a dimension of 5 × 5 μm. Photon packets were launched from one surface location at the origin of the medium and absorption was scored during propagation according to the appropriate grid element in a two-dimensional matrix A(y, z), where 1 ≤ y ≤ k and 1 ≤ z ≤ l. Figures 1(a) and 1(c) show the geometry of the models used to simulate the absorption distribution, and Figs. 1(b) and 1(d) show the absorption distribution for the multilayered MC code [5] and the newly developed code, respectively. Note that cylindrical regions were included within the new model to exercise the boundary condition code for arbitrary geometries, but that the optical properties of each cylindrical region are identical to its surrounding medium and as such the resulting absorption distribution should match that calculated by the multilayered MC code.

 

Fig. 1 Absorption distribution in a two-layered medium simulated using 10 million photon packets. (a) Model of the two layered medium. (b) Absorption distribution in the model simulated with an existing MC code for multilayer samples [5]. (c) Model of the two layered medium comprising two different cylindrical regions at a depth of 750 μm and 250 μm. Optical properties of the cylindrical regions are identical to surrounding tissue. (d) Absorption distribution in the model simulated with the new MC model. Scale bars correspond to 100 μm. Colorbar is in logarithmic scale (dB).

Download Full Size | PDF

3. Modeling OCT transverse and axial signal localization

The aforementioned methods allow MC-based calculation of light transport in arbitrary three-dimensional structures. To extend this model to simulate OCT imaging, it is necessary to include methods for proving both transverse and axial signal localization (e.g., the imaging component of OCT). Transverse resolution can be incorporated as has been done in prior models [3, 9]. Briefly, photon packets are launched with a location defined by a uniform circular distribution. The angular position of photon packets is sampled uniformly between 0 and 2π while the radial position is sampled between 0 and a 10-μm beam radius with a square root distribution. On collection, exiting photon packets are included only if they meet the spatial and angular constraints of a detection aperture. In our model, a spatial filter with a 10-μm radius (uniform intensity) and an angular filter limited to a +/− 5° detection angle [3] was adopted. The location of the launch and collection apertures is translated across the imaging surface to simulate beam scanning.

In our approach, we have modeled directly the generation of wavelength-dependent fringe signals, and we then use Fourier analysis to convert these fringes to depth-resolved A-lines. The simulation tool records the square root of the diffuse reflectance (weight) and the total path length of each photon packet that passes the detection aperture. At the conclusion of the MC simulation for a given beam location (e.g, the modeling of one A-line), we obtain a set of detected photons packets, each with a diffuse reflectance amplitude and a path length. We then construct the received sample signal (amplitude and phase) as a function of wavelength according to

The reference field is calculated as

We interfere the sample and reference arm fields in a balanced detection configuration to generate a fringe signal given by

We then follow conventional Fourier-domain processing by applying a window (Hamming) and Fourier transform to give the OCT signal as a function of depth. Figure 2(a) shows an OCT interference signal, I_D(k), as a function of wavelength obtained with a two-layered medium (air and cutaneous tissue), where the first layer starts at 0 μm and the second, at 50 μm. Figure 2(b) presents the corresponding i_D(z) signal as a function of depth following Fourier transform. Inspection of Fig. 2(b) reveals a typical OCT signal where the two regions are clearly visible.

 

Fig. 2 Simulated A-line signals. (a) Interference signal, I_D(k), as a function of wavelength and corresponding (b) OCT signal, i_D(z), as a function of depth following Fourier transform.

Download Full Size | PDF

To validate the OCT signal localization methods, we generated simulated OCT structural images of a simple phantom structure with a homogeneous background and two vessel inclusions. Figure 3(a) shows the geometry of the medium comprising the two blood vessels with a diameter of 150 μm and 100 μm. Figure 3(b) presents the corresponding structural OCT image showing the vessels with the typical shadowing effect and omnipresent speckle noise.

 

Fig. 3 Simulated blood vessels. (a) Two blood vessels modeled at a depth of 250 μm and 200 μm. (b) Structural OCT image of two cylindrical blood vessels (outlined in yellow) surrounded by a homogeneous sample. The parameter μ_s for blood is fixed to 650 cm⁻¹, μ_a, to 5 cm⁻¹, g, to 0.9888, and n, to 1.37. The parameter μ_s for the sample is set to 10 cm⁻¹, μ_a, to 1 cm⁻¹, g, to 0.7, and n, to 1.37. Scale bars correspond to 100 μm.

Download Full Size | PDF

4. Simulating blood flow

OCT angiography is based on the acquisition of repeated measurements from the same location and the analysis of these measurements to identify time-varying signals that result from flowing blood. To simulate angiographic OCT images, it is therefore necessary to derive time-series measurements that are constant or modulated according to their degree of interaction with intravascular scatterers. In this work, we make an assumption that the intravascular scatterers have translated by more than the larger of the axial or transverse resolutions. As a result there is no correlation between the OCT signals scattered from vessels across two time points. This accurately describes many of the angiographic OCT approaches that employ large time separations between measurements to achieve sensitivity for slower flowing vessels.

To model temporal decorrelation of intravascular scatterers, we used the following approach. First, during MC light transport, we flag each photon packet (using the vessel flag variable, F_n) that scatters within a voxel associated with the region defining the intravascular space (section 2). Then, during calculation of the OCT signal I_D1 (k), we include in each photon packet that is flagged a random phasor, ϕ_1n, between 0 and 2π within the expontential. To simulate a second measurement, I_D2 (k), at the same location, the same MC solution (R_n and l_n) is used but a new set of random phasors, ϕ_2n, are applied. In this way, the OCT signals derived from photon packets that have not interacted with intravascular scatterers are constant in time, but those signals that are derived, at least in part, from photon packets scattered from intravascular scatterers are proportionally modulated. In this work, the OCT angiography contrast is then calculated from the depth-directional variance of the differential signal between i_D1 (z) and i_D2 (z) [2]. However, it is also possible to model other angiographic algorithms that operate within the regime of full signal decorrelation between measurements. We note that the effect of vessel size on the OCT signal is accounted for by the inherent statistical nature of Monte Carlo methods and not, for example, by adjusting the magnitude of the phase (ϕ_n) in the random phasor according to vessel size. A smaller vessel will have a larger proportion of photon packets that traversed the vessel without a scattering (and are therefore not phase modulated) than would a larger vessel. We also note that the simulation of time-series data does not require that the MC simulation of light transport be recalculated, and thus adds minimal additional computational burden.

The MC model algorithm is summarized in Fig. 4. For clarity, the function of each step is represented by its box line style. Boxes outlined by a full line describe photon transport in current MC models [5], by a dashed line, photon propagation in a discretized medium, by a dotted line, OCT-based signal localization, and by a dash-dotted line, modeling of flowing blood. This algorithm was implemented in a C language environment and was developed to provide acceleration of simulation time by offering the option to run on multithreaded Central Processing Units (CPU) or CUDA Graphics Processing Units (GPU) [10] depending on available hardware. The use of parallel threads enables many photon packets to be simulated simultaneously. For instance, simulation time to generate Fig. 1(b) on an Intel Core i7 processor with the existing downloadable model was 152 seconds with a single thread, while simulation time to obtain Fig. 1(d) with the new model required 3 seconds with 12 threads. The improvements in algorithm speed are necessary to make complex three dimensional simulations practical.

 

Fig. 4 Flow chart of the MC algorithm.

Download Full Size | PDF

While simulation time is largely dependent upon hardware, it is also affected by the number of launched photon packets and cuboids implemented in the model. While the use of a large number of photons and cuboids in simulations improves signal precision and three-dimensional shape resolution, it may rapidly become computationally impractical. The user must therefore appropriately choose optimal parameters for a specific application to obtain desired data within a reasonable timeframe.

5. Validation

To validate the model, we compared the appearence and signal statistics of a simulated angiographic image with that of angiographic images acquired in vivo (Figs. 5(d) and 5(h)). The simulated phantom included vessels designed to match that of a reference in vivo image acquired with an existing OCT-based system on a mouse skin tumor through a dorsal skinfold window preparation. The in vivo data was acquired at an A-line rate of 50 kHz with a 1300 nm system. A smaller region of the OCT image containing two larger blood vessels (diameters of approximately 250 μm (left) and 70 μm) was extracted.

 

Fig. 5 Cross-sectional OCT-based images of two blood vessels simulated using different optical properties. Top row shows structural images while bottom row presents angiography images. Optical properties of the skin and blood were extracted from literature. For normal skin, values were set according to [11] for a wavelength of 1,200 nm. The scattering coefficient, μ_s, was fixed to 16 cm⁻¹, the absorption coefficient, μ_a, to 5 cm⁻¹, and the anisotropy factor, g, to 0.715. Values for blood were chosen according to prior studies [12, 13] where optical properties were measured for whole blood with physiological hematocrit percentages of 42% and 46%, respectively. (a, e) Simulated vessels obtained with μ_s for blood fixed to 650 cm⁻¹, μ_a, to 5 cm⁻¹, and g, to 0.97 [13]. (b, f) Simulated vessels generated with μ_s for blood set to 650 cm⁻¹, μ_a, to 5 cm⁻¹, and g, to 0.9888. (c, g) Simulated vessels obtained with μ_s for blood equal to 650 cm⁻¹, μ_a, to 5 cm⁻¹, and g, to 0.995 [12]. (d, h) Experimental images of skin vessels acquired with an existing OCT-based system on a murine model. Scale bars correspond to 100 μm.

Download Full Size | PDF

Figure 5 shows cross-sectional images of the two blood vessels simulated using different optical properties based on published values [12, 13]. While values for absorption and scattering coefficients are in good agreement for wavelengths above ∼ 1,050 nm [13], the anisotropy factor in [12] remains larger (g ≃ 0.995 [12] compared to g ≃ 0.97 [13]). Figures 5(a) and 5(e) were generated based on values published in [13] while Figs. 5(c) and 5(g) were obtained from data presented in [12]. Figures 5(b) and 5(f) were simulated using an intermediate anisotropy value. All results show the two blood vessels with the OCT shadowing effect below each vessel. Based on these simulations, the scattering coefficient for blood, μ_s, was set to 650 cm⁻¹, the absorption coefficient, μ_a, to 5 cm⁻¹, and the anisotropy factor, g, to 0.9888.

Figure 6 compares the angiographic signal statistics in simulated and in vivo datasets. Histograms were generated over a region of interest (ROI) in Figs. 6(a) and 6(b) with an axial depth of 100 μm and a width of 200 μm. Inspection reveals agreement in signal distributions (Figs. 6(c) and 6(d)).

 

Fig. 6 Comparison of angiography images showing two blood vessels. Histograms were generated for ROIs surrounded by a red box to quantitatively compare (a, c) simulated and (b, d) in vivo angiography datasets. Scale bars correspond to 100 μm.

Download Full Size | PDF

6. Discussion and conclusion

This paper describes a new MC-based simulation tool for computing OCT angiography images. The tool offers three key innovations relative to existing numerical models of OCT. First, a voxel-based region definition allows arbitrary three-dimensional structures including complex vascular networks. Second, MC methods are used to capture photon packets and their propagation distance, and Fourier-domain coherence gating methods are used to generate depth-resolved fringes. Third, a method allowing simulation of flow-induced temporal decorrelation was implemented which enables angiographic image generation.

While this tool provides simulation of angiographic signals, it does not accurately recapitulate the noise present in angiographic imaging of static scatterers. This will be a focus of future work, and would enable direct comparison of the contrast of angiographic signals relative to background for varying angiographic algorithms. The tool could be further optimized to include the Gaussian beam profile on the launched and collection apertures; the current system uses a circular aperture. While these changes to the aperture will allow a more accurate modeling of the system transverse resolution, they will not affect the statistics of the angiographic signals. This tool can be complemented by adding methods to mathematically generate realistic microvascular networks for simulation of three-dimensional OCT imaging. This would allow segmentation algorithms to be optimized on datasets in which the true vessel geometries are known a priori. Finally, this model may be extendable to allow the modeling of partial temporal decorrelation, and would allow the testing and development of quantitative flow measurement algorithms in addition to the morphological angiography described in this work.