1. Introduction

OCT signal simulation opens unique possibilities for development/optimization of various OCT modalities in highly controllable, easily modifiable conditions, which cannot be enabled in physical experiments with real biological tissues and even phantoms. Numerous research groups make significant effort in developing various tools to simulate OCT signals for general and specific purposes, e.g. [1–16]. Main approaches to OCT-scan modeling roughly can be divided in such groups as Monte-Carlo methods [2,4,5], utilization of the point-spread function concept [12,13] and wave-based methods [6–9] (or “full-wave” according to the terminology used in [6]) describing in various approximations the forward propagation of the optical wave towards scatterers and then its backward propagation and collection over the receiving aperture

In many models the scatterers are approximated by discrete scattering centers, properties of which can be more complex than point scatterers (e.g., accounting for frequency dependence) to better imitate real tissues. The use of discrete scatterers concerns all three groups of simulation methods (e.g. [10,17,18,19]). This is quite a meaningful approximation, even if real scatterers are not literally point-like objects; nevertheless, they still can be treated as discrete sub-resolution particles in OCT scans and be reasonably corresponded to cells in biological tissues. Also it should be recalled that the very fact of successful realization of OCT imaging is based on utilization of ballistic back-scattering. Even within the single-scattering approximation one can study many important specific features of OCT images (such as evolution of speckle patterns in the presence of scatterer motions, statistical properties of OCT speckles, etc.)

Multiple scattering as a degrading factor in OCT (the ability to describe which is intrinsic to Monte-Carlo methods) often represents secondary interest, because the very principle of OCT is based on single-scattering and many important OCT-scan features can be well captured even in this approximation. It is noteworthy that, depending on the situation, speckles in OCT may be either considered as a degrading factor or play a useful role of high diagnostic value [1,20]. For instance, motion-induced speckle variability is the basis of OCT-angiography [21]. Therefore, it is reasonable that the majority of scan-formation models in OCT are based on ballistic approximation.

For the illuminating-beam description in such models, various additional simplifications are usually made, which may be more or less realistic. In simpler variants like [1,17], the illuminating beams may have uniform phase in the cross section (which still may be a reasonable approximation for OCT-systems with weakly focused beams without pronounced curvature of wave fronts, whereas the amplitude profile may be variable, e.g. Gaussian). Even using such simplified beam geometries as in [14,17], it is already possible to simulate important features of real OCT-scans formed by beams with weak focusing/divergence and obtain useful results for developing new OCT modalities such as elastography [22–25]. For weakly focused beams, the validity of such simplified models is also confirmed by similarity with physical experiments and results of more advanced models like [18,19] that belong to the full-wave approach and are based on rigorous consideration of Gaussian beams with allowance for pronounced focusing.

It should be pointed out that for simulations of realistic densities of scatterers (imitating the density of cells in biological tissues) and biologically relevant sizes of visualized regions (∼hundreds of microns and greater), sufficiently high computational efficacy of the simulations is of key importance, especially when large series of OCT-data volumes are required to study effects of moving scatterers. In this context utilization of Gaussian beams in modeling is rather attractive because of the possibility to analytically describe their forward propagation. Collection of the back-propagated waves over the receiving aperture can be made either using numerical integration (like in [18]), or even analytically for non-diaphragmed beams [19]. Certainly, passing from numerical integration to analytical description of the back-propagated wave collection over the receiving aperture like in [19] dramatically improves the computational efficiency. However, this improvement is enabled at the expense of a limitation: the Gaussian beam diameter should be significantly smaller than the illuminating/receiving aperture. For non-Gaussian beams (and even for a Gaussian beam with a sharp diaphragm or somewhat distorted by phase aberrations), analytical description of the illuminating-beam propagation becomes either impossible or requires a radical reformulation of all basic equations. Numerical realizations (concerning both the forward/backward propagation and collection of the backscattered fields over the receiving aperture) is very computationally demanding. Consequently, development of realistic and flexible models of OCT-scan formation requiring minimal simplifications/approximations and at the same time enabling high computationally efficiency remains challenging and highly demanded.

In what follows, we demonstrate that using the same basic approximations (i.e., the use of discrete scatterers and ballistic scattering) as in [10,17–19] and many other models it is possible to reformulate the OCT-scan formation model, such that the above-mentioned limitations can be essentially lifted. As a result, the modified model acquires the following previously unavailable features: (i) it may simulate OCT scans using arbitrary forms of illuminating beams (Gaussian, Bessel, etc.), including the possibility to introduce sharp diaphragms; (ii) arbitrary aberrations of the illuminating beam can be introduced (such that forward and backward propagation through the aberrating screen is rigorously accounted for); (iii) the widely accepted paraxial approximation is not required, so that wide-angle beams can be used; (iv) arbitrary motions of scatterers between subsequent scans can be readily introduced in the model; (v) although the widely used approximation of ballistic scattering at discrete scatterers is used, there are no other limitations on the density of scatterers, their distribution in space and scattering strengths including possible frequency-dependence.

The modified model described in the next section uses full-spectral representation of the illuminating/scattered fields (both in the axial direction, like is physically used in spectral-domain OCT, and in the lateral directions the beam is also represented in the spectral form instead of the spatial-domain representation). Due to such modifications, the new formulation of the model becomes very straightforward and strikingly compact. The possibility of using arbitrary beam profiles is not made at the expense of computational efficiency; on the contrary, the model computationally is very fast. To the best of our knowledge, there were no previously discussed models that could simultaneously enable such a range of possibilities.

We firstly demonstrate the possibilities of the proposed model by validating it via comparison with the reference model [19] for the same set of scatterers. Model [19] is based on exact analytical solutions for non-diaphragmed focused Gaussian beams. Then the universality of the new model is illustrated by comparing OCT images obtained for Gaussian and Bessel beams (again for identical sets of scatterers). Next, in view of high importance of accounting for aberrations for some applications of OCT (e.g., in ophthalmology [26]), the unique flexibility of the new model is demonstrated by generating images with various aberrations represented by Zernike functions that are often used for descriptions of eye-lens aberrations in ophthalmology. Finally, the model usage in configurations with moving scatterers (typical of angiographic and elastographic applications) is demonstrated.

2. Method of full-wave angular-spectrum description of OCT-scan formation

In the proposed description, each spectral component of the illuminating-beam with wavenumber k is characterized by a complex-amplitude distribution ${U_L}(x,y,k)$ at the output lens (with allowance for a sharp diaphragm). The lens is located in the (x,y) plane with the axial coordinate $z = {z_d}$ and z-axis is directed to the tissue depth along the axis of the illuminating beam as shown in Fig. 1. Coordinate ${z_d}$ may either directly correspond to the tissue boundary or the boundary of the immersion layer (in which scatterers are absent). The lateral coordinates $({x_0},{y_0})$ of the beam axis define the position of a single A-scan. The lateral scanning required for obtaining 3D volumes of OCT data corresponds to variation in $({x_0},{y_0})$. Each discrete scatterer is characterized by coordinates $({x_s},{y_s},{z_s})$. In models of spectral-domain OCT with the illuminating beam represented as a set of spectral components with wavenumbers k, the forward propagation of the axial spectral components and their back-propagation from scatterers is often described in the spatial domain [18,19]. However, for the same purpose, the angular-spectrum approach [27,28] can be used to describe the optical-wave evolution. For example, after obtaining 3D images in [18,19], the numerical refocusing was made by recalculating the wave fronts represented by their angular (lateral) spatial spectra. In what follows, we apply such a spectral-domain description (both in the axial and lateral directions) to describe the wave fields for all stages: first, for illuminating-beam propagation towards the scatterers and then for backward propagation of scattered fields and their collection over the illuminating/receiving aperture.

 

Fig. 1. Schematically shown geometry of the illuminating beam (not necessarily focused like in the figure) with the lateral coordinates $({x_0},{y_0})$ of the axis directed to the tissue depth along the z-axis of the coordinate system. Coordinates $({x_s},{y_s},{z_s})$ characterize positions of scatterers. The lens shown in the figure depicts not only the output window of the OCT probe, but also other internal optical elements of the OCT system coupled with the illuminating/receiving aperture.

Download Full Size | PDF

In such a representation, the complex-valued field distribution ${U_L}(x,y,k)$ over the input/output aperture is characterized by the following spatial (angular) spectrum $g({k_x},{k_y},k)$ (for simplicity, below we omit constant factors before the integrals):

After propagation of the illuminating beam between the input plane $z = {z_d}$ and plane $z = {z_s}$ corresponding to the depth of a scatterer with coordinate ${z_s}$, the illuminating-wave spectrum can be written as ${h_s}({k_x},{k_y},k)g({k_x},{k_y},k)$, where the propagator function ${h_s}({k_x},{k_y},k)$ has the following form (assuming that $\exp (ikz)$ corresponds to the in-depth propagation along z-axis):

Note that a similar method for describing the axial evolution of the transversal spectrum of the illuminating beam using the propagator (2) was recently applied in paper [16] to describe degradation of the point-spread function in the presence of refractive-index fluctuations. According to Eqs. (1) and (2), for the scatterer with coordinates $({x_s},{y_s},{z_s})$, the transversal spatial spectrum of the incident wave is ${h_s}({k_x},{k_y},k)g({k_x},{k_y},k)$. Then its Fourier transform over spectral coordinates $({k_x},{k_y})$ yields the following complex-valued amplitude ${U_s}({x_0},{y_0},k)$ of the field incident on the scatterer located at $({x_s},{y_s},{z_s})$ if lateral coordinates of the beam are $({x_0},{y_0})$:

Notice that in the left-hand side of Eq. (3) the subscript “s” means that ${U_s}({x_0},{y_0},k)$ depends not only on the position of the beam axis $({x_0},{y_0})$, but also on coordinates $({x_s},{y_s},{z_s})$ of the s-th scatterer.

The possibility to describe the field ${U_s}({x_0},{y_0},k)$ in the spectral form (3) is a very useful feature, because the transversal spectrum $g({k_x},{k_y},k)$ can readily be found numerically for an arbitrary beam profile ${U_L}(x,y,k)$. This should be made only once for all scatterers and does not require additional assumptions (such as a Gaussian profile, absence of a diaphragm, etc.) which are usually needed to simplify analytical considerations.

In can also be noted that by analogy with the spectral function $g({k_x},{k_y},k)$ related via Eq. (1) with the initial beam profile in spatial domain ${U_L}(x,y,k)$), for the spectral form of propagator ${h_s}({k_x},{k_y},k)$, its spatial-domain counterpart ${H_s}(x,y,k)$ can be introduced via the inverse Fourier transform over spectral coordinates $({k_x},{k_y})$:

Let us recall that according to Eq. (3), the transversal spectrum of ${U_s}({x_0},{y_0},k)$ is given by the product of the Fourier transforms $g({k_x},{k_y},k)$ and ${h_s}({k_x},{k_y},k)$ that correspond to the spatial-domain functions ${U_L}(x,y,k)$ and ${H_s}(x,y,k)$. Then according to the theorem about the product of spectra, function ${U_s}({x_0},{y_0},k)$ can also be expressed as the following convolution of functions ${U_L}(x,y,k)$ and ${H_s}(x,y,k)$:

In Eq. (5) the integral is found over the aperture $\Omega $ of the lens (see Fig. 1). In fact ${H_s}(x,y)$ is the Green’s function for a localized scatterer located at $({x_s},{y_s},{z_s})$ and Eq. (5) is a spatial-domain counterpart of the spectral representation of the incident wave given by Eq. (3).

The scattered wave for a discrete scatterer located at $({x_s},{y_s},{z_s})$ is proportional to the incident-wave amplitude ${U_s}({x_0},{y_0},k)$ and the scattering strength ${K_s}$ for this scatterer. This scattered field propagates back to the illuminating/receiving aperture and is collected by the lens. This stage is again described by the convolution of the Green’s function and receiving aperture similarly to Eq. (5). Then the complex-valued amplitude ${B_s}({x_0},{y_0},k)$ of the spectral component with wave number k collected at the receiving aperture is given by

Here, similarly to ${U_s}$, subscript “s” means that ${B_s}({x_0},{y_0},k)$ corresponds to the collected field from the scatterer with coordinates $({x_s},{y_s},{z_s})$. We point out that the pre-integral factor ${U_s}({x_0},{y_0},k)$ in Eq. (6) corresponds to the incident-field amplitude, whereas the integral factor describing the collection of the back-scattered wave is identical to integral in Eq. (5) that describes the forward propagation. In view of this, Eq. (6) can be rewritten in the very compact form as the incident field squared:

At this point it should be recalled that Eq. (6) describes the received complex-valued amplitude ${B_s}({x_0},{y_0},k)$ for a single spectral component with the wavenumber $k = {k_n}$, whereas in spectral-domain OCT, the optical-beam contains a large number of illuminating spectral components $\exp ({i{k_n}z} )$ with wavenumbers ${k_n}$. In previous papers [17,19] (where spectral amplitudes ${B_s}({x_0},{y_0},k)$ were found differently than above) we discussed in detail that superposition of $n = 1..N$ spectral components with wavenumbers ${k_n}$ forms a pixelated A-scan composed of $q = 1..N$ pixels. This superposition mathematically corresponds to Fourier transformation over the axial wavenumbers, which is performed in spectral-domain OCT scanners to form in-depth A-scans. Using the notations similar to those used in works [17–19] we can write the following expression for the complex-valued amplitude $A({x_0},{y_0},{z_q})$ of q-th pixel with the axial (depth) coordinate ${z_q}$ and coordinates $({x_0},{y_0})$ in the lateral plane:

Here, for the discussed two-quadrature signal that is enabled by dual-arm OCT setups, $D = \pi /\delta k$ is the maximal depth range that is unambiguously imaged by the used set of the illuminating spectral components ${k_n}$ separated by $\delta k = {k_{n + 1}} - {k_n}$. Notice, that for single-arm systems enabling only one-quadrature signals, the unambiguously imaged depth is twice smaller [11,29,30]. In Eq. (8) summation over $n = 1..N$ is made for the contributions of all spectral harmonics taking into account their amplitudes $S({k_n})$ in the illuminating-beam spectrum (the sign in the exponential function in (8) agrees with the chosen sign $\exp ( + ikz)$ for the in-depth propagation along z-axis). The second summation over index s in Eq. (8) is made over contributions of all scatterers located within the beam. In a more compact form Eq. (8) can be written as the Fourier transform over ${k_n}$

It should be emphasized that although Eq. (3) in the spectral form and Eq. (5) in the spatial-domain equivalently represent the same amplitude ${U_s}({x_0},{y_0},k)$, the possibility of its efficient evaluation in the spectral form (3) is the central point for computational efficiency of the described model. Indeed, spectrum $g({k_x},{k_y},k)$ readily allows one to represent an arbitrary phase-amplitude form ${U_L}(x,y,k)$ of the illuminating beam via Eq. (1). Therefore, the spectral representation confers the above-described model unprecedented computational speed for treating arbitrary initial beam profiles ${U_L}(x,y,k)$ without radical reformulation. In particular, not only Gaussian beams, but other beam forms can be used (e.g., Bessel or Laguerre-Gauss beams that attract high interest in recent years). Also aberrations of the illumination beam at the boundary of the visualized region can easily be introduced, which is of interest for OCT applications in ophthalmology, etc. For any beam profile ${U_L}(x,y,k)$, the transversal spectrum $g({k_x},{k_y},k)$ should be found only once and used for all A-scans and scatterers. Next, it should be pointed out that the spectral representation of the propagator ${h_s}({k_x},{k_y},k)$ can be used in the exact form (2) and does not require the use paraxial approximation, which is often necessary for obtaining analytical results in other models.

Finally, one can readily imitate motions of scatterers by varying coordinates of scatterers $({x_s},{y_s},{z_s})$ for subsequently obtained scans. Notice that for motionless scatterers, estimation of Fourier transform given by Eq. (3) can be made using rapid FFT-algorithms to obtain (via Eqs. (7) and (8)) spectral amplitudes ${B_s}({x_0},{y_0},k)$ for the entire set of A-scan positions $({x_0},{y_0})$ within the simulated 3D volume of OCT-data. For displaced scatterers, the data should be recalculated and the total computation time will depend on the total number of $({x_0},{y_0})$-coordinates of the illuminating beam within the 3D data set and the required number of recalculations for moving scatterers. Depending on their velocities, for certain cases, displacements may be sufficient to consider between entire 3D volumes. In other cases variations in $({x_s},{y_s},{z_s})$ may be significant even on inter-A-scan intervals requiring larger amount of recalculations. In the next section we present examples of utilization of the proposed model to illustrate its main new features/capabilities.

3. Simulation examples to demonstrate main model capabilities

In the below-presented examples we assume that the central wavelength of the OCT system is 1.3 μm and spectral width is 90 nm, which are “typical” parameters for OCT. Other parameters of the illuminating beam and scatterers are not the same and are specified for each of the examples. Here, we simulate the OCT data for rectangular volumes with sizes 128 × 128 × 128 pixels corresponding to biologically-relevant physical sizes of 512 μm in depth and 192 μm in the lateral directions. Thus, the lateral size of the pixel is equal to 1.5 μm, which corresponds to spatial step of the positions $x_0^{(i)}$ and $y_0^{(i)}$ of the simulated A-scans. This set of 128 × 128 positions in the lateral directions is used to perform 2D FFT operations for finding the angular spectrum $g({k_x},{k_y},k)$ defined by Eq. (1) for the chosen beam profile ${U_L}(x,y,k)$. Then, for a given position $({x_s},{y_s},{z_s})$ of a scatterer, the inverse 2D FFT operations represented by Eq. (3) yield the 2D matrix of the complex-valued amplitudes ${U_s}({x_0},{y_0},k)$ for the same set of $x_0^{(i)}$ and $y_0^{(i)}$. Next, utilizing Eq. (7) we obtain the spectral amplitudes ${B_s}$ for every lateral position $({x_0},{y_0})$ (in the examples presented below, in Eq. (7) the scattering strengths of the scatterers was supposed identical, ${K_s} = 1$). Finally, using summation represented by Eqs. (8) or (9) using conventional for spectral-domain OCT FFT operations, one obtains the entire 3D set of complex-valued amplitudes of the simulated OCT image of the chosen size 128 × 128 × 128 pixels. We also point out that for all simulated examples, the OCT probe is located above the simulated volume.

As the first example, for validation of the proposed model, it makes sense to compare the results of the new spectral formulation given by Eqs. (1), (3), (7) and (9) with the reference analytical model [19]. The latter is based on exact analytical solutions for the spectral amplitudes ${B_s}({x_0},{y_0},k)$ obtained for a non-diaphragmed Gaussian beam (i.e., assuming that the lens aperture is sufficiently big and can be considered infinite). For comparison of the present model with model [19], we accepted identical parameters of a focused Gaussian beam with a focus radius of $3$ µm and generated scans for a set of 7 scatterers with identical lateral coordinates and different depths, including the focal depth $z = 128$ µm. Figures 2(a) and 2(b) demonstrate that the described full spectral-domain formulation based on discrete Fourier transforms gives the spatial forms of point-spread functions very close to the results obtained using the exact analytical model [19]. The depth profiles shown in Fig. 2(c) confirm that the two methods demonstrate excellent coincidence of the signal amplitude, so that minor differences can be observed only in the zones sufficiently distant from the scatterers, where the signal has very low amplitudes (at a level about -60 dB below the maxima and lower). This coincidence with the exact analytical results for a pronouncedly focused Gaussian beam confirms the reliability of the proposed model, for which, however, the main new feature is the ability to describe OCT signals formed by arbitrary beams, not necessarily axially symmetric ones.

 

Fig. 2. Comparison of OCT images for the same set of scatterers and the same pronouncedly focused Gaussian beam simulated using: (a) reference model [19] (based on exact analytical solutions obtained in spatial domain for Gaussians beams) and (b) the present universal model formulated in spectral domain. (c) shows the depth profiles along the beam axis for the two models (asterisks are for model [19] and circles correspond to the present model).

Download Full Size | PDF

The following examples demonstrate that the ability of the model to describe arbitrary beams opens the possibility to efficiency simulate various aberrations. Such aberrations may characterize, for example, the eye cornea and lens, which is of interest for ophthalmic applications of OCT. In the following examples, the aberrations are modeled by introducing phase distortions of an initially Gaussian beam.

To clearer demonstrate how various aberration alter the PSF form, we again intentionally took a small number (ten) of scatterers located at a depth of the focal plane, which may be considered as imitation of focusing at the eye retina. Figure 3 demonstrates the simulated scans in (x,y)-plane (C-scans) without (Fig. 3(a)) and with aberrations of several characteristic types (Figs. 3(b), 3(c) and 3(d)). Such aberrations may arise in OCT-imaging of the eye retina surface through cornea and eye lens with imperfect refraction [26], phase distortions of which are often approximated by a set of Zernike polynomials [31]). It is clear from Fig. 3 that both worsening of the PSF-localization in the focal plane and aberration-induced translational shifts of the entire C-scan can readily be simulated. Such OCT images with a priori known aberrations can be used for development and comparison of aberration-correction procedures. The latter are important, e.g., for improving the quality of OCT-imaging for ophthalmic applications [26,32–34] (a special discussion of utilization of simulated OCT images for such applications will be published later).

 

Fig. 3. Examples of simulated OCT C-scans in $(x,y)$ plane corresponding to the focus depth $z = 128$ μm for a focused Gaussian beam similar to that, for which $(x,z)$ images in Fig. 2 were shown. Panel (a) with well-localized scatterer images is for the absence of aberrations; (b) is for the aberration corresponding to the Zernike polynomial $Z_3^1$ (coma); (c) is for the Zernike polynomial $Z_2^2$ (2^nd astigmatism); (d) is for Zernike polynomial $Z_6^0$ (2^nd spherical aberration).

Download Full Size | PDF

Next, the ability of the model to easily simulate non-Gaussian beams is demonstrated for a Bessel beam (such beam attract significant interest in OCT due to nearly non-diverging character of propagation [35,36]). In the modeling, we imitated formation of such a beam due to propagation of an initially planar beam with a Gaussian amplitude profile $\exp ( - r_ \bot ^2/{a^2})$ through a cylindrical axicone with the opening angle $2\gamma $ (see [37]). To this end we introduced in the input plane $z = {z_d}$ a radially-dependent factor $\exp [i{k_n}{r_ \bot }tg(\gamma )]$ describing the phase delay acquired by the spectral components of the beam after propagating through the axicone. In the B-scans shown in Fig. 4(a)–1 and 4(b)-1 we again intentionally used a small number (six) scatterers located at three depths. Figure 4 clearly demonstrates the difference between the propagation of the simulated Bessel beam (row (b) in Fig. 4) and a Gaussian beam with a similar radius of the focal area (row (a) in Fig. 4). The Gaussian beam demonstrates pronounced focusing and strong broadening of the PSF outside the focal depth. In contrast, the Bessel beam demonstrates good conservation of the beam-kernel radius, which is accompanied by slow evolution of the surrounding rings corresponding to the Bessel-function sidelobes (compare the C-scans in rows (a) and (b) in Fig. 4).

 

Fig. 4. Examples of OCT B-scans in $(x,z)$ plane (first column) and C-scans in $(x,y)$ plane (2^nd through 4^th columns) at three depths simulated for a Gaussian beam (panels (a-1) through (a-4)) and a Bessel beam (panels (b-1) through (b-4)). It is clear that, for the Gaussian beam, the minimal size of the PSF corresponds to the focal depth $z = 248$ μm, whereas outside the focal depth the PSF sizes are much greater, such that for $z = 452$ μm, the two scatterers are hardly resolved. In contrast, the Bessel beam at all three depth has nearly invariable size of the central kernel, whereas the outer ring-shape sidelobes described by the Bessel function gradually evolve remaining at least 19 dB lower than the amplitude at the axis. In this simulation, the axicone angle is $\gamma = \pi /40$ rad.

Download Full Size | PDF

The next examples in Fig. 5 demonstrate the intrinsic ease of the proposed model to simulate 3D data formed by illuminating beams with arbitrary degree of focusing, which can be used for development/verification of methods of numerical refocusing. In addition we demonstrate the ease of simulation of various noises (by adding a random complex-valued number to every pixel to imitate the noise at the receiving array). Figures 5(a) and 5(b) show examples of the initial OCT image for a highly-focused Gaussian beam before and after numerical refocusing (with the refocusing algorithm used in [38,18]). The focal plane is located closer to the bottom (the focus depth $z = 480$ μm and the focus radius is 3 μm). The density of scatterers is intentionally chosen fairly low to clearer demonstrate different forms of the point-spread function (PSF) outside the focus depth and near the focal plane at $z = 480$ μm. It is seen that outside the focal depth the weakly localized scatterer images are essentially masked by the measurement noise (the SNR is ∼12dB near the focal plane and gets lower beyond the focal depth). However, after performing numerical refocusing one obtains nearly identical highly localized PSF-forms independently of the depth, i.e., maximal possible resolution the same as at the focal depth is enabled everywhere. Consequently, the scatterers become well visible against the noise even outside the physical-focus depth ($z = 480$ μm).

 

Fig. 5. Examples of refocusing of simulated 3D scans for a focused Gaussian beam in the presence of noises and moving scatterers (vertically oriented flow of scatterers). (a) is a 3D image of a rarified set of scatterers; notice that above the focal depths, the scatterer images are strongly broadened and are nearly masked by the imitated Gaussian “measurement noise”. (b) is the same image after refocusing, due to which the initially defocused scatterer images became well localized and clearly seen everywhere. (c) is a 3D image for a density of scatterers closer to density of biological cells in tissues with an axial flow 10 μm in diameter imitating a blood capillary which is initially indistinguishable against the background; (d) the same refocused image, in which the capillary outside the focal depth becomes well visible with the maximal resolution and looks as a dark channel, which is due to flow-related breaking of the constructive interference in the refocused image.

Download Full Size | PDF

Besides simulations of 3D volumes of OCT data for motionless scatterers, the model allows one to readily re-simulate OCT images for arbitrarily displaced positions of scatterers. Such possibilities are important for testing algorithms of OCT-based elastography and angiography in highly controllable conditions with flexibly modifiable parameters. In physical experiments, even using phantoms, it is very challenging to realize such imaging in highly controllable conditions. In this context, the next examples demonstrate simulations of digital refocusing of 3D OCT-data in the presence of a scatterer flow to imitate the blood flow in a micro-vessel surrounded by the “solid” tissue. We demonstrate that even a fairly slow motion of scatterers may break conditions of constructive interference required to transform a broadened scatterer image out of focus into a well focused image. The refocusing is made by re-assembling signals from the same scatterer, the broadened image of which is present in a stack of overlapped neighboring scans. The refocusing algorithm numerically introduces appropriate phase corrections and performs reassembling of scattered signals. After such signal re-assembling for motionless scatterers, the signals collected over their broadened PSFs produce well-localized bright PSFs, like shown in Figs. 5(a) and 5(b). However, for moving scatterers, their motion introduces inter-frame phase shifts. These additional motion-related phase shifts perturb the conditions of constructive interference required for refocusing. Consequently, after refocusing these moving scatterers produce regions of very low intensity and become highly contrasted with bright motionless scatterers. This effect may be used for singling out moving scatterers against the background signal from motionless scatterers. As illustrated by Fig. 5(c), the moving and motionless scatterers are indistinguishable before refocusing, whereas Fig. 5(d) demonstrates that after refocusing the simulated vessel becomes well visible due to strongly reduced signal level. Notice that near the depths of the physical focus, where the PSF has the minimal size and scatterer images nearly do not overlap for neighboring A-scans, the refocusing only slightly changes the OCT image. Consequently, the images of moving scatterers do not change their intensity near the focal depth. By this reason, the flow after refocusing becomes visible as a dark channel only sufficiently far from the focal plane, where the initial size of PSF is much larger than in the focal area. Notice that after refocusing the PSF size becomes comparable with the focal-spot size everywhere, so that even quite narrow channels with flow can be visualized with a high resolution as illustrated in Fig. 5(d). In this example, particle displacement between the subsequent A-scans is $\Delta z = {10^{ - 4}}$ μm, which corresponds to the vertical velocity $\Delta z = 4 \cdot {10^{ - 3}}$ mm/s assuming 40 kHz A-scan rate. The total amount of scatterers is $2.5 \cdot {10^5}$ to imitate the density of real biological cells.

Finally, the next example relates to elastographic imaging of local strains in the simulated region, in which a prescribed strain-induced variation in the positions of scatterers is introduced. Strain mapping is the basic step in realization of compression optical coherence elastography (C-OCE) which in recent years opened a broad range of biomedical applications [25]. Such applications comprise elastographic delineation of normal versus tumorous tissue and even differentiation of molecular subtypes of tumors [39–41], demonstrations of detailed morphological segmentation of heterogeneous tumor tissues using differences in elasticity among various morphological components [42,43]; studying complex 3D strain fields of thermo-mechanical origin in corneal tissue subjected to laser treatment in the context of refraction correction [44]; detection of insufficient stability of laser-fabricated cartilaginous implants [45]; studying deformations of osmotic origin during tissue impregnation by non-isotonic liquids [46,25]. In view of these very different applications, numerical simulations of exactly controllable arbitrarily complex deformations with the possibility to flexibly vary OCT-system parameters, signal-to-noise ratios, etc., - open unprecedented prospects for detailed multi-parameter studies of elastographic processing in OCT.

For demonstration, in Fig. 6(a) we show a structural image of a simulated 3D volume 192 × 192 × 512 microns in size. Then to all scatterers axial displacements were prescribed to produce the spatially-inhomogeneous axial-strain distribution shown in Fig. 6(b), for which another OCT scan was simulated. The simulated inhomogeneous strain distribution resembles strains produced in experiments [44,49] in corneal samples subjected to heating by an infra-red laser beam. The simulated OCT images in Fig. 6 are obtained for the same highly-focused Gaussian illuminating beam as for the examples in Fig. 5 assuming the same focal depth $z = 480$ μm (near the bottom). The focal-waist radius is $3$ μm also like for Fig. 5. The total amount of scatterers is $2.5 \cdot {10^5}$. The inter-frame phase difference between the initial and deformed simulated scans was found and used for reconstructing strains by estimating axial phase-variation gradients as proposed in [47]. Figure 6(c) demonstrates the rather complex distribution of inter-frame phase variations corresponding to the prescribed strain shown in Fig. 6(b). These phase variations are found in the absence measurement noises other than the strain-induced “decorrelation noise”. However, a very noisy structure of the phase-variations is clearly seen above the focal depth in regions of pronouncedly laterally-inhomogeneous strains. In these regions, large groups of differently displaced scatterers within the OCT-beam cross section contribute to formation of individual A-scans. Consequently, superposition of waves back-propagated from these differently-displaced scatterers produce very noisy resultant phase variation. The next Fig. 6(d) shows the much cleaner phase-variation map obtained after preliminary refocusing, which suppresses this degrading interference of multiple differently-displaced scatterers.

 

Fig. 6. Examples of elastography-related simulations of OCT scans. (a) is simulated structural image in which the brighter lines separate regions subjected to different straining. (b) is the prescribed distribution of strains according to which displacements of scatterers were defined (the inset shows the cross section of the central part of the volume through the plane labeled by the black dashed line. (c) is the phase variation between the initial and deformed states. (d) is the phase-variation map with much better quality obtained with the preliminary performed refocusing of the compared OCT images. (e) reconstructed strain distribution based on the noisy phase-variation map (c) found without refocusing. (f) similar strain map corresponding to phase-variation map (d) found with preliminary refocusing. (g) schematically shown real experimental configuration in which thermomechanical strain distribution similar to the simulated example (b) was produced by laser heating (like in works [44,49]): 1- the OCT probe; 2- buffer silicone layers between which the studied samples of cartilaginous or corneal tissues (layer 3) were placed; 4 – infrared laser operating at 1.5 μm that produced thermo-mechanical deformations in the tissue. (h) a real example of interframe strains resembling the simulated strain distribution and visualized by processing phase-resolved OCT scans using the vector method proposed in [48].

Download Full Size | PDF

In phase-sensitive OCT, the sought axial strains are proportional to the axial gradients of interframe phase-variations, for estimation of which either conventional least-square fitting [47], or more advanced “vector” method [24,48] can be applied as in the present example. The so-found strain distribution corresponding to the straightforwardly found rather noisy interframe phase variations in Fig. 6(c) are shown in Fig. 6(e) that also demonstrates strong distortions/errors in strain estimation. These errors are especially strong outside the focus depth and in the regions, where the phase variations have strong lateral inhomogeneity and look very noisy in Fig. 6(c). In contrast, the strain map of much better quality shown in Fig. 6(f) is found using the phase-variation map shown in Fig. 6(d) that was obtained after preliminary refocusing. The reconstructed strain map in Fig. 6(f) is much closer to the exact (prescribed) strain distribution (Fig. 6(b)) with much better quality and uniform lateral resolution over the entire 3D scan.

To demonstrate that strain distributions similar to the simulated one do occur in real situations, Fig. 6 g schematically shows an experimental configuration, in which thermo-mechanical strains were produced by irradiating a sample of collagenous tissue (a piece of cartilage or eye cornea) by a beam of infra-red laser. Such studies relate to the development of new methods for non-surgical reshaping of cartilaginous or corneal tissues (see a more detailed discussion in [44,45]). The buffer silicone layers shown in Fig. 6(g) serve for performing quantitative compression elastography [25] to assess the irradiation-produced variations in the tissue elasticity and make conclusions about the microstructural alterations in the tissue [49]. Figure 6 h demonstrates a real approximately axially-symmetric experimental example of in-depth section of interframe thermo-mechanical alternating-sign strains that was laser-induced in a sandwich sample consisting of rabbit cornea and buffer silicone layers. The strain was reconstructed from sequentially acquired OCT scans using the vector method [48]. Such high-resolution OCT-based elastographic imaging based on approach [48] attracts much attention is the recent years [50–52], so that the possibility to realistically simulate OCT images of deformed tissues should be very helpful for further perfection of the OCT-based strain/elasticity imaging.

4. Discussion

The key new feature of the proposed compact, but rigorous angular-spectrum-based formulation of the wave-based model of OCT-scan formation, is that the model describes arbitrary-form beams, not requiring the paraxial approximation. This ability to utilize arbitrary phase-amplitude beam profiles opens very convenient possibilities for modeling arbitrary aberrations (see Fig. 3), which is important in the context of the development/comparison of various methods for their numerical correction. Also various degrees of focusing and qualitatively different types of illumination beams (like Gaussian or Bessel ones in Fig. 4) can efficiently be studied in a uniform manner using Eqs. (1),(3), (7) and (9), in which only the beam profile ${U_L}(x,y,k)$ should be changed.

Of no less importance is the ability of the proposed formulation to validate the applicability of previously developed approximate models of OCT-scan formation, in which special forms of the beam profile were assumed (such as the neglect of the beam-front curvature, assumption of non-diaphragmed Gaussian beams, paraxial approximation, etc.) It is important that such simplified models can enable incomparable computational efficiency for simulating large sequences of OCT data. Such large sequences are required, e.g., for studying OCT-scan evolution in the presence of flows or Brownian motion of scatterers, or for developing of numerical methods to compensate masking bulk motions of living tissues in OCT-angiography [53,54]. For such applications, more sophisticated models may require unacceptable in practice computation time (days and weeks) or need unreasonably expensive supercomputing facilities. By these reasons approximate models continue to attract high attention (see, e.g., recent work [55]). However, the drawback of approximate/simplified models is that, remaining within the framework of the used approximations, it is impossible to reliably verify the range of parameters, for which these simpler models are applicable. The above-described model readily allows one to validate the applicability of various simplified model formulations such as [17] or model [19] that have already been used for studies of OCT-scan evolution due to variation in the OCT-setup parameters and tissue deformation. For example, they were applied in comparative studies of the correlation-based and phase-resolved elastographic processing of OCT-scans [23] and perfection of strain and displacement reconstruction using phase-sensitive OCT [22,48]. Now the applicability of such models can be reliably validated in fairly arbitrary conditions. Such possibility as the intrinsic ease in imitating measurement noises (as illustrated in Fig. 5) is quite straightforward, so that we do not specially focus on this issue (procedures of introducing the measurement noises in more detail were discussed in [48]).

It can be noted that the model also allows one to incorporate such a degrading effect as roll-off by performing, similarly to work [14], convolution in the k-space of the spectral amplitudes ${B_s}({x_0},{y_0},{k_n})$ in Eq. (8) with a function describing the integration profile of the elements of the photodetector array

The model also can readily account for the material dispersion, because of which the equidistant (in dispersionless case) set $\{ {k_n}\} $ of wavenumbers becomes non-equidistant. For a dispersive material it is sufficient to use the set $\{ {k_n}\} $ corresponding to the dispersion law. The forms of the propagator function and the other expressions for finding the complex-valued amplitudes ${B_s}({x_0},{y_0},{k_n})$ collected at the receiving aperture remain unchanged, although they acquire dispersion-related phase distortions that should cause degradation of the image constructed via the same FFT procedures as in the absence of dispersion.

Next, the optical beam decay due to absorption and/or scattering can also be introduced by adding an imaginary part to the wavenumbers. Especially simple is the case of frequency-independent losses. In this situation it is sufficient to apply to the formed image a single exponential factor $\exp [ - 2({\mu _a} + {\mu _s})z]$ describing the beam decay during the forth-and-back propagation, where ${\mu _a}$ and ${\mu _s}$ are the decay coefficients due to absorption and scattering, respectively (see, e.g., [14,24]). Concerning properties of scatterers, arbitrary scattering strengths (including those dependent on ${k_n}$) can readily be introduced, so that we do not specially focus on these issue.

Less trivial is the situation with motion of scatterers, examples of which for a flow-type motion in Figs. 5 and displacements of scatterers due to strain in Figs. 6. For these examples, the images for subsequent scatterer configurations were obtained as if the signal acquisition is “instantaneous”, which implies that during the acquisition time the vertical displacements of scatterers are much smaller compared with the wavelength and lateral ones are much smaller than the beam diameter. For modern OCT setups with the acquisition rates ∼ several tens kHz and higher, these conditions are usually fulfilled for such applications as OCT-visualization of microcirculation [53] and visualization of strains in compression elastography [25]. However, for larger scatterer velocities, the finite acquisition time may lead to fringe washout (for fast axial motions) causing signal reduction and also (for lateral motions) apparent broadening of the beamwidth [56,57]. The described model also opens the possibility to simulate such effects by imitating the finite acquisition time using summation of a series of complex-valued scans. As a result of this superposition, the image parts with sufficiently rapidly moving scatterers will demonstrate the above-mentioned effects, such as apparent beam broadening or signal reduction, the physics of the latter being rather similar to the flow-induced signal reduction demonstrated in Fig. 5 for refocusing to regions with moving particles.

The main limitation of the described wave-based model is the use of single-scattering (ballistic) approximation. It should be pointed out, however, that the very principle of OCT-visualization is based on this approximation and the multiple-scattering influence in OCT is only one of various degrading factors. Multiple-scattering effect in many cases may be of secondary importance, although even degradation of PSF due to multiple forward scattering by refractive-index inhomogeneities can also be introduced in the described approach by analogy with study [16]. In fact for the majority of other model formulations, the neglect of multiple scattering is typical, except for Monte-Carlo-based models. The latter are intrinsically better suitable for considering “snake-like” photon trajectories including multiple scattering, but computationally the Monte-Carlo approach is known to be very demanding.

In this context, impressively high computational efficiently of the present model can be specially emphasized. The model does not use computationally demanding solvers to describe the wave propagation or numerical integration for collecting signals from scatterers. The main Eqs. (1),(3), (7), and (9) rely on Fourier transformations that can be performed using efficient algorithms of discrete fast Fourier transform (FFT). The discrete transversal spectrum $g({k_x},{k_y},k)$ of an arbitrary beam profile ${U_L}(x,y,k)$ should be found only once according to Eq. (1) for a chosen mesh on the $(x,y)$-plane to represent the positions of A-scans. Then for all of the chosen positions $({x_0},{y_0})$, the 2D Fourier transformation over wavenumbers $({k_x},{k_y})$ in Eq. (3) gives the amplitudes ${U_s}({x_0},{y_0},{k_n})$ of the waves incident onto the scatterers. Consequently, all spectral amplitudes ${B_s}({x_0},{y_0},{k_n})$ at the receiving array are readily found (see Eq. (7)). Then the computationally efficient FFT over the wave-numbers ${k_n}$ and straightforward summation over all scatterers are performed in Eq. (9), which yields the entire 3D OCT-image.

The required computation time linearly depends on the number of spectral harmonics in the source spectrum $S({k_n})$ (determining the number of pixels in the axial direction). Similarly, the linear proportionality holds for the total number of A-scan positions in $(x,y)$-plane (i.e., the number of pixels in lateral directions) and the total number of scatterers. To give an idea about the required time, consider a volume 128 × 128 × 128 pixels (corresponding to 192 × 192 microns laterally and 512 microns in depth) used in Figs. 2 through 6. For a small amount of motionless scatterers (say, <100) and one core of Intel I-7 6950X CPU, the 3D images require the computation time ∼several seconds. If the same 3D volume contains 250000 scatterers (which imitates the density of biological cells ∼4-6 microns in size), the same CPU requires ∼100 minutes (using 10 cores of the same CPU, although we did not specially optimized parallelization and the acceleration factor was only ∼3). Note that entire 3D images with densely located scatterers in many cases are not needed and simulation of a B-scan or C-scan formed by a much smaller amount of scatterers is sufficient. Then for the same density of scatterers as above, their number can be reduced by a factor of ∼100 and the simulation time drops down to ∼1 minute. These results are about twice faster than simulation for the same configuration using the reference model [19] based on the exact analytical expressions for ${B_s}({k_n})$ found in the special case of non-diaphragmed Gaussian beams. In that model the simulation of individual A-scans is reduced to FFT over wave-numbers ${k_n}$ and summation over contributions of all scatterers like in Eqs.(8) and (9). Nevertheless, for simulating 3D OCT-scans, the present model has demonstrated comparable and even higher computation efficiency than model [19], but without the requirement to mandatory usage of non-diaphragmed Gaussian beams. Certainly, utilization of multi-core GPUs may additionally dramatically reduce the above-indicated and already quite practically acceptable CPU-based computational time.

5. Conclusions

In this paper we presented a new form of computationally efficient full-wave model of OCT-scan formation. The model formulation is straightforward and compact enabling a rich and flexible set of capabilities/features: (i) the illuminating beam may have arbitrary phase-amplitude profile with allowance of a sharp diaphragm and presence of aberrations; (ii) the model does not use paraxial approximation that limits the degree of focusing/divergence of the illuminating beam; (iii) although the broadly used approximation of single (ballistic) scattering by discrete scatterers is assumed, there are no other limitations on the density of scatterers, their distribution in space and scattering strengths with possible frequency-dependence; (iv) besides rigorous accounting for focusing/divergence, factors describing the wave decay can readily be introduced by analogy with Monte-Carlo approaches; (v) arbitrary measurement noises also can easily be added to simulate required signal-to-noise ratios. Also the model readily allows one to re-simulate the data for arbitrary changed scatterer positions (e.g., random or flow/deformation-induced). Thus, within the framework of ballistic scattering by discrete scatterers, the above-listed intrinsic abilities of the model give reasons to characterize it as fairly comprehensive. This flexibility, rich set of abilities and high computational efficiency of the model open unique possibilities for studying OCT-scan properties and developing novel diagnostic OCT-based modalities for biomedical diagnostics.