Computationally efficient model of OCT scan formation by focused beams and its usage to demonstrate a novel principle of OCT-angiography

In this paper, we present a computationally highly efficient full-wave model of optical coherence tomography (OCT)-scan formation by focused beams in spectral-domain OCT. Similarly to some previous models, it is based on the summation of fields scattered by discrete sub-resolution scatterers, and enables one to account for the axial and lateral inhomogeneity of the illuminating Gaussian beam. The main feature ensuring the high computational efficiency of the described model is that, rather than numerical integration of the scattered signal over the receiving aperture, we apply an analytical description of both the illuminating-beam focusing and the collection of the scattered signals over the receiving aperture. Elimination of numerical integration over the receiving aperture increases the overall computation speed by a factor of ∼ 103 . This is of key importance in terms of the practical feasibility of simulations of 3D OCT data volumes for large amounts (∼ 105−107 ) of scatterers, corresponding to realistic densities of cells in biological tissues. We demonstrate the model’s possibilities by simulating the digital refocusing of strongly focused OCT beams in the presence moving scatterers, thereby presenting a novel principle of contrast-agent-free visualization of scatterer flows, with velocities typical of blood microcirculation.


Introduction
The development of analytical and numerical models of scan formation in optical coherence tomography (OCT) has been attracting a great deal of attention over the past two decades, 1 Author to whom any correspondence should be addressed since the publication of the first works [1,2] to describe the formation of OCT scans, including their peculiar speckle structure. The latter is sometimes considered as a negative factor, degrading the image quality, so that various methods of speckle suppression have been discussed in the literature (e.g. [3,4], and references therein). However, analysis of speckle structure can also be an informative source in the context of the development of multimodal OCT, with particular reference to elastographic [5][6][7] and angiographic applications [8,9]. A popular trend in the development of OCT-scan modeling is the Monte-Carlo approach [10][11][12]. Monte-Carlo methods are usually rather computationally demanding (since they intrinsically require the simulation of pseudo-random trajectories of huge numbers of photons in turbid media). Moreover, such an approach is not well-suited to the generation of numerous data volumes, corresponding to the evolution of OCT scans, due to the motion of scatterers, although the latter are indispensable for testing signal processing in OCT elastography and angiography, where the locations of scatterers are random, but their motions may correspond to either deterministic or random laws.
Other simulation methods have also been proposed, in which propagation/scattering is described as a regular process (although the medium is turbid). Such methods may consider fluctuations in the refractive index with a full-wave numerical simulation of beam propagation and scattering [13], which is also very demanding computationally. There are also models based on the utilization of discrete scatterers, either for the explicit representation of point-spread function for subresolution scatterers [14,15], or to obtain the point-spread function in the axial direction by summing the spectral components, which is similar to the process for spectral-domain OCT devices [16,17].
The approach [16] based on the summation of signals from discrete sub-resolution scatterers is highly computationally efficient, and is very convenient for generating a series of OCT scans, formed by arbitrarily moving scatterers; this process has been employed efficiently in simulations of elastographic processing in OCT (e.g. [7,[18][19][20],). However, only weakly focused OCT beams with depth-independent radii may be described by this model [16], leading to the proposal of an extension for describing the influence of pronounced beam focusing/defocusing in [21]. In [21], the influence of the divergence of the Gaussian illuminating beam is described, following rigorous analysis. The spherically diverged backscattered signal was convoluted with the receiving aperture, having the same phase-amplitude mask that initially formed the incident illuminating beam. This collection of the backscattered signal over the receiving aperture in [21] was performed using numerical integration, which, again, is rather computationally demanding.
Nevertheless, in comparison with Monte-Carlo methods, or finite-element full-wave approaches, the model demonstrated in [21] is significantly less demanding computationally, and makes it possible to simulate a series of 2D B-scans at intervals ranging from several minutes to several tens of minutes (depending on the number of scatterers) using Intel-I7 class CPUs, obviating the need for multi-core GPU-based calculations. However, even for such computationally low-demand models, such as that proposed in [21], the required computation times increase up to several hours, or even days, for the generation of a set of 3D OCT scans. These are indispensable in e.g. modeling the digital refocusing of OCT data acquired using strongly focused illuminating beams. Below, we present a computationally much faster modification of the semi-analytical full-wave model [21], in which the same basic principles are used to describe the illuminating beams and the signals scattered by localized sub-resolution scatterers. We emphasize that, similarly to the case for real tissue, scatterers may move along either regular or random trajectories, may demonstrate dependence of the scattering properties on the wavelength, etc. The key difference from [21] is that the integration of the backscattered spectral harmonics over the receiving aperture in the model described below is performed analytically, rather than numerically, and that only the summation of the complex-valued spectral amplitudes thus established should be performed numerically for all scatterers. Due to the elimination of numerical integration, the computational efficiently of the model described in this paper is greatly increased (by a factor of 10 2 . . . 10 3 depending on the number of scatterers) in comparison with [21].
In this paper, the new computationally efficient form of the semi-analytical full-wave model is applied to simulate 3D volumes of OCT data, consisting of hundreds of B-scans, in the presence of both static and moving scatterers, including flows imitating perfused blood vessels. The resulting acceleration of simulations paves the way for previously impractical operations, such as performing highly controllable and realistic massive numerical experiments on the formation of OCT scans in various conditions. In particular, we demonstrate the application of the developed model to studying refocusing procedures; these have attracted a great deal of attention in recent years, due to the attendant possibility of achieving increased lateral resolution in OCT over the entire imaged depth [22]. Based on simulations of the digital refocusing of strongly focused OCT beams in the presence of moving scatterers, we demonstrate a novel principle of contrast-agent-free visualization of scatterer flows, with velocities typical of blood microcirculation.

Optimized model of OCT-signal formation by a strongly focused beam
In model [21], used here as a basis, the focused illuminating field is described as a Gaussian beam incident onto the imaged tissue through an immersion layer, as shown in figure 1.
In spectral-domain OCT, an elementary step consists of the formation of a single one-dimensional A-scan in the direction of the z-axis, aligned with the probing-beam axis. The complex-valued amplitude A(z q )of the qth pixel in an A-scan consisting of N pixels, with the coordinate of the centers equal to z q , is given by the summation of all spectral harmonics with wavenumbers k n , received by the array of photodetectors subsequent to the backpropagation of the illuminating OCT beam from all scatterers with coordinates r s : Here, the visualized depth H is determined by the distance between the neighboring spectral components H = π/|k n+1 − k n |. The complex amplitudes B (⃗ r s , k n ) of the received spectral Schematic illustrating a focused illuminating Gaussian beam and a sub-resolution scatterer, with coordinates − → rs , producing a spherically-diverging scattered wave. Here, W 0 is the beam radius in the focal plane, and W(z) is its current radius; z d is the distance of the receiving aperture from the tissue surface (adapted from [21]. © 2019 Astro Ltd. All rights reserved). components are jointly determined by the forward propagation of the illuminating beam, the backward propagation of the spherically diverging fields scattered by the subresolution scatterers, and the collection of those backpropagated scattered waves over the receiving aperture (indicated in figure 1 by a dash-dotted line).
Equation (1) for the complex-valued pixel amplitudes can also be represented as where FFT −1 kn denotes the inverse Fourier transform, applied to the received complex-valued amplitudes .
Let us recall the integral expressions obtained in [21] for the illuminating beam amplitude u (⃗ r s , k n ), and the signal B (⃗ r s , k n ), collected at the receiving aperture. For a Gaussian beam incident at a localized (sub-wavelength) scatterer with coordinates − → r s = (x s , y s , z s ), the complex-valued amplitude, u − → r s , k ,of a wave with wavenumber k can be written in the form: Here, the focus position is − → r 0 = (x 0 , y 0 , z 0 ), where z 0 is the axial coordinate, and the lateral coordinates correspond to radius-vector − → r 0⊥ = (x 0 , y 0 ). Quantity W (z, z 0 ) determines the current amplitude at the axis of the illuminating beam, and has the following meaning of its current radius: where W 0 is the beam radius in the focal zone; λ = 2 · π/k represents the wavelength. Quantity Φ (z, z 0 ) = arctg λ · (z − z 0 ) /(π · W 2 0 ) is the phase difference between the Gaussian beam and an ideally plane wave front; the current radius of wave-front curvature in the beam is given by where R (z, z 0 ) tends to infinity in the focal region, corresponding to z = z 0 . Factor S(k) in equation (3) describes the spectral form of the probing beam. In principle, function S(k) may be considered to be jointly determined by the spectral properties of the source and the receiving system. The field u scat (⃗ r, ⇀ r s , k), scattered by a sub-resolution scatterer with coordinates⃗ r s and scattering strength K s (k), corresponds to a spherically diverging wave of the form The collection of the scattered field over the receiving aperture subsequent to its backward propagation and passage through the optical system can be expressed in terms of the following integral expression for the sought complex-valued amplitude B − → r s , k of the received signal for each spectral component k = k n , and each scatterer: where − → r d⊥ = (x 0 + ξ, y 0 + η) denotes the radius-vector of the integration point on the receiving aperture (represented by a dashed line in figure 1), where the integration coordinates ξ, η are counted from the coordinates (x 0 , y 0 ) of the illuminating-beam axis; vector − → r d has the following meaning, − → r d = (x 0 + ξ, y 0 + η, z d ), where z d is the axial coordinate of the receiving aperture. The last two exponential factors in the integral (7) describe transformations of the amplitude and phase of the scattered signal during its backward passage through the optical system that initially formed the illuminating Gaussian beam.
In study [21], integral expression (7) was evaluated numerically, this evaluation being the most time-consuming procedure of that model. To accelerate the evaluation of integral (7) without an appreciable loss of accuracy, the following observations can be made: taking into account that z d < 0 and z s > 0, we consider the scatterers satisfying the condition (z s − z d ) >> πW 2 0 /λ, so that it is then possible to approximately represent the argument in the first exponential function in equation (7): (8) Integral (7) can then be represented in a simplified form: For a sufficiently large receiving aperture, in comparison with the illuminating-beam cross section, simplified integral (9) can analytically be evaluated as: where the factor P is given by: The exponential decay of signals propagating though tissue is not yet accounted for in equations (3)-(11), and should, if necessary, be introduced in addition.
Overall, in comparison with the initial form proposed in [21], where integral (7) was numerically evaluated, the analytical representation (10) of this integral accelerates the computations required to simulate 3D volumes of OCT data by a factor ∼10 2 −10 3 (depending on the number of scatterers). This significant acceleration paves the way for previously impractical operations in terms of the generation of large amounts of 3D data, corresponding to the densities of scatterers typical of cell densities in real biological tissues. Indeed, for other models (including [21], which was faster in comparison with the alternative variants mentioned in the introduction), simulations similar to that presented below would require either practically unacceptable time intervals (several days or greater), or the necessity of utilizing very powerful and expensive computational means to reduce the computation time.

Examples of refocusing of simulated OCT scans including the presence of moving scatterers
Bearing in mind the level of interest in the digital refocusing of OCT scans obtained via a highly-focused illuminating beam, an important application of the developed model is its utilization for the development of refocusing procedures, and studying the effect of the motion of scatterers on the results of digital refocusing. By analogy with [21], in the section below we will use the following spectral transformations corresponding to the beam refocusing [22]: where A ∆z (x, y, z) is the optical field corresponding to the numerical shift of the focal plane by a distance of ∆z, (x, y) are the lateral coordinates, and z is the axial coordinate (depth). As in the above equations, FFT and FFT −1 denote the forward and inverse Fourier transforms for the respective pairs of variables indicated in equation (12), A 0 (x, y, z) denotes the complex-valued signal amplitude before the numerical shift of the focal plane, and (u, v) are the lateral coordinates in the Fourier domain for the 2D-Fourier-transformed function A 0 (x, y, z). For the purposes of the discussion below, it will be sufficient to use equation (12), based on subsequent forward and inverse Fourier transforms, although other refocusing algorithms can also be used (e.g. convolution of the signal with a specially constructed function [22]). Here, we consider discrete sub-resolution scatterers with a frequency-independent scattering coefficient, which is a reasonable approximation for many real situations. We assume that the illuminating beam has a central wavelength of 1310 nm, and a spectral width of ∼90 nm, corresponding to parameters widely used in real OCT applications. We also assume that the lens plane, over which integration is performed, is located at z d = 100µm above the tissue. The beam focus is located at a depth of z o = 128µm inside the tissue. The 1/e radius of the beam in the focal waist is W 0 = 3µm for a beam with pronounced focusing, such as that shown in figure 1.
In the following examples we will show relatively smallsize stacks, comprising 128 B-scans along the y-axis, covering 192 µm. Along the x-axis, each B-scan also covers 192 µm and comprises 128 pixels, so that the horizontal interpixel step in the (x,y) plane is 1.5 µm in both directions. Each A-scan has an in-depth length of 128 pixels, covering 512 µm in air, so that the vertical pixel size is close to 4 µm, and the imaged volume is therefore 192×192×512 microns in size.
To demonstrate that the numerical estimation of integral (7) used in [21] and its analytical estimate (given by equation (10)) coincide with a high degree of accuracy, we firstly consider images of seven well-separated localized scatterers, loc-ated along the vertical axis at depths z s = 32, 64, 128, 200, 300, 400 and 480 µm. The images of these scatterers prior to digital refocusing are shown in figure 2(a-1), where the left half of the images illustrates the direct numerical estimation of integral (7) and the right half depicts the analytical estimate (10) of this integral. The vertical profiles through the center of the scatterer images are shown in figure 2(a-2). these clearly demonstrate that the meaningful parts of the images coincide with a high degree of accuracy, and that only in the dark region between the scatterers (∼50-60 dB below the maximal intensities) is there some degree of discrepancy, but such small signal amplitude should be below the noise level in real images. Therefore, in the following examples, we will demonstrate only the results based on the utilization of the computationally much more efficient analytical estimate (10).
In figures 2(b-1) and (b-2), the result of refocusing via equation (12) over the entire depth is shown, from which it is clear that refocusing enables uniform lateral resolution similar to that in figure 2(a) for the scatterer located at the initial focal depth z s = z 0 = 128 µm. Figure 2(b) demonstrates, however, that the amplitudes of the refocused images of identical localized scatterers outside the focus have lower amplitudes than those in the focus.
To explain this fact, we must return to the structure of factor P that enters equation (10) for the received spectral amplitude B − → r s , k . It should be noted that in equation (11) for factor P, the first real-part term is related to purely geometrical variations in the amplitude of the illuminating beam. If this geometrical effect solely determined the received signal, then, performing refocusing via equation (12) by appropriately correcting the phase fronts described by the exponential factors with imaginary arguments in equation (10), it would be possible to collect the distance-dependent energy of the scattered signal and perform its in-phase summation over the receiving aperture. Consequently, for identical scatterers, the amplitudes of their refocused images would be identical. However, equation (11) for factor P also contains the other two imaginary terms. These terms are related to inconsistent structure of the phase in the received sphericalwaves from the scatterers (with the curvature radius |z s − z 0 |), and the curvature R (z d , z 0 ), corresponding to the phase transformation of a plane wave front passing through the optical system that collects the scattered signals over the receiving aperture. This optical system fully compensates the receivedwaveform curvature only for scatterers located in the focal plane, i.e. for z s = z 0 , so that the imaginary part in equation (11) disappears. For a scatterer with z s ̸ = z 0 , the spherical divergence of the scattered wave is not fully compensated by the lens, so that the summation is not perfectly in-phase over the receiving aperture, and the resultant amplitude of the refocused images is lower than for z s = z 0 . Therefore, in order to obtain identical amplitudes for the refocused images of identical scatterers with various depths z s , the following correction factor K corr should additionally be applied to equation (10): The additional application of the correcting factor (13) equalizes the amplitudes of the refocused images of the localized scatterers, as shown in figures 2(c-1) and (c-2).
Note that in figures 2(a)-(c) noise is absent. It is clear that application of the correcting factor (13) also increases the level of noise, together with the useful signal. However, refocusing collects the noise incoherently, whereas the useful scattered signal is collected coherently, so that the signalto-noise ratio after refocusing is increased; thus, even scatterers hidden by noise prior to refocusing may become visible after refocusing. This is demonstrated in figure 2(d-1), where the left part of the image is similar to the image in figure 2(a-1), but with added noise (simulating the noise at the receiving array) at an intensity of ∼ −50 dB with respect to the intensity of the scatterer image near the physical focus at z 0 = 128 µm. It is clear that in the lower part of the left half of the simulated scan in figure 2(d-1), the images of scatterers are almost hidden by the noise. In contrast, the right half of figure 2(d-1) shows that after refocusing and application of the correcting factor (13) the scatterers become clearly visible against the noise, and their amplitudes are nearly equal, as expected for identical scatterers. After such verification that the results of simulations of OCT scans and subsequent application of refocusing procedures correspond well to the results expected from the physical considerations for motionless scatterers, we next consider the influence of the motion of scatterers.
In the following examples, instead of using rarified scatterers, we increase their concentrations to make them comparable with the density of cells in real biological tissues, where we assume that the cells have characteristic diameters of ∼5-6 µm. This corresponds to effective densities of scatterers on the order of (4-6)x10 6 mm −3 . For the generation of 3D volumes of OCT data with such concentrations of scatterers, a numerical evaluation of integral (7) over the receiving aperture may require unacceptably long calculations (i.e. tens of hours or even days using the CPU of a desktop computer, or requiring specialized and expensive high-performance computational means to accelerate the calculations). However, the acceleration due to analytical estimate (10) of the signal collected over the receiving aperture opens up the possibility of performing such simulations using a 'typical' desktop computer (the results shown below were obtained using a 16-cores Intel I7 6950X CPU, and required several hours to simulate 3D volumes of 192×192×512 microns in size, containing 2 · 10 5 scatterers).
Based on the discussions above, it is anticipated that the refocusing of moving scatterers by collecting defocused images of scatterers via digital correction of the phase fronts (as represented in equation (13)) should be differently affected by the lateral and axial motions of scatterers. It  (7), as in model [21], and the right half of panel (a-1) is simulated using the analytical estimate (10) of equation (7). Panel (a-2) shows vertical profiles corresponding to the numerical (circles) and analytical (solid line) estimates of integral (7), corresponding to the left and right halves of panel (a-1). Panel (b-1) corresponds to refocused image (a-1), using equation (12) without correction of amplitudes of the refocused images of scatterers; the corresponding intensity profile through the scatterers is shown in (b-2). Panel (c-1) is the same refocused image as in (b-1), but with additional correction factor (12) the application of which results in the intensities of the identical scatterers having equal amplitudes, as shown in the corresponding vertical profile in (c-2). The left half of panel (d-1) is similar to the unfocused image (a-1), but with the amplitude-correction factor (13) and added noise − 20 dB with respect to the scatterer image at the focal depth, so that the out-of-focus scatterers are significantly masked. The right half of (d-1) is obtained by digital refocusing of the image in the left half of (d-1) with the addition of the amplitude-correction factor (13). The vertical profiles corresponding to image (d-1) before (solid line) and after refocusing (dotted line) are shown in panel (d-2). is clear that for purely geometrical reasons, lateral motions cause weaker phase distortions in comparison with axial motions of the same velocity. Note that by the same reasoning, the Doppler frequency shift of the OCT signal for axial motions of scatterers is much greater than for lateral motions.
In OCT-based angiography, the main focus of interest tends to be the visualization of microvasculature with diameters of vessels below ∼50 microns, whereas the smallest capillaries may have diameters comparable with the sizes of red blood cells (∼5-6 microns). The velocities of scatterers in such fairly thin vessels typically range from a few mm/s to fractions of mm/s. Since in real situations, the orientations of the vessels are not ideally orthogonal to the OCT-beam axis, there is usually an axial projection of the flow velocity of at least ∼several per cent of the total velocity. In view of this, we consider simulated examples below, in which a vertical flow of scatterers is introduced, so as to imitate a vessel with a 10 µm radius. Figure 3 shows an example of such a simulation for the same conditions as for individual scatterers in figures 2(a)-(c) (in the absence of additive noise), but for randomly distributed 2 · 10 5 scatterers, the density of which is close to the concentration of cells in biological tissue. In the simulation, the axially moved particles are displaced by 2.25 × 10 −4 µm between the consequent A-scans, corresponding to a velocity of 2.25 µm s −1 for the acquisition rate of the OCT scanner described in [9,23], and applied for angiographic mapping in [24][25][26][27] using the high-pass filtering principle. Figure 3 Figures 3(c) and (d) show the same B-scan and intensity profile, subsequent to the application of digital refocusing over the entire visualization depth, so that the presence of the vessel becomes clearly visible. Indeed, it is anticipated that, when collecting the scattered fields over the cross-section of the beam to perform refocusing, the phase distortions caused by the motion of scatterers should perturb the inter-scan phasing and reduce the amplitude of regions with moving scatterers in the resultant refocused image. This effect should be stronger in out-of-focus regions, where the beam is significantly wider than for in-focus regions. Consequently, the number of overlapped B-scans, from which the refocused image of a scatterer is formed, is greater than in the focal plane. Unlike in-focus regions, where only the nearest B-scans are partially overlapped, in out-of-focus regions, significantly larger stacks of overlapped B-scans are acquired during larger time intervals, so that the motion of scatterers causes stronger mutual de-phasing of the overlapped B-scans. During refocusing, this de-phasing impedes the in-phase summation of signals from moving scatterers, so that the flow region has a significantly lower intensity. However, closer to the focus, where the illuminating beam is insufficiently wide, the motion-induced dephasing is not yet able to visualize the flow in the refocused image.
With this observation in mind, in order to enable a higher degree of contrast of the so-visualized 'vessel' visualized over the desired depth range (0 < z < 500 µm), it makes  Examples of focusing/refocusing influence on OCT scans of 3D volumes, in which a cylindrical flow of vertically moving scatterers is introduced to imitate a blood vessel. In this example the vertical displacement of the particles is 2.5 × 10 −4 µm between A-scans. Here, the physical focus of the illuminating beam with radius W 0 = 3µm is intentionally placed below the depth of interest at z 0 = 628 µm. Panel (a-1) shows the B-scan passing through the 'vessel' axis before refocusing, and panel (b-1) depicts the refocused B-scan at the same position, where the vessel becomes clearly visible as a region of strongly reduced amplitude due to signal-phase perturbations caused by the motion of scatterers. Panel (c-1) shows a cross-sectional C-scan through the 'vessel' at depth z = 300 µm (represented by horizontal dashed lines in (a-1) and (b-1)) before refocusing, and (d-1) shows the same cross-section after refocusing. Panels (a-2) and (b-2) show the vertical profiles along the 'vessel' axis, indicated by dashed lines in B-scans (a-1) and (b-1). Panels (c-2) and (d-2) show the horizontal profiles through the 'vessel' along the dashed lines depicted in C-scans (c-1) and (d-1). sense to place the focus somewhat deeper. Figure 4 corresponds to such a case, where, similarly to figures 2 and 3, the depths z < 500 µm represent the main interest for imaging, but the location of the physical focus is intentionally selected to be deeper (at z 0 = 628 µm), so that within the imaged region, the beam is significantly wider than in the focal region. Such a configuration ensures sufficient sensitivity of refocusing to scatterer motions, and the illuminatingbeam convergence may also be helpful in compensating for signal attenuation due to absorption and scattering in real conditions. Figure 4(a-1) shows the initial B-scan, in which the 'vessel' is not visible against the background material, and figure 4(b-1) shows the refocused B-scan for the same position, where the 'vessel' can clearly be observed as a stripe with strongly reduced intensity, for the reasons given above. Figures 4(a-2) and (b-2) show the profiles of figures 4(a-1) and (b-1) along the vertical dashed lines, where the reduced intensity within the vessel's cross section after refocusing can also be clearly observed. Figures 4(c-1) and (d-1) depict the horizontal C-scans plotted through the dashed lines shown in figures 4(a-1) and (b-1). In the C-scan image, prior to refocusing, the vessel cannot be seen against the background, similarly to figure 4(a-1), whereas after refocusing, the vessel cross-section can clearly be observed in the C-scan in figure 4(d-1). The corresponding intensity profiles through the center of the vessel's cross sections in the initial and refocused C-scans are shown in figures 4(c-2) and (d-2), where, subsequent to refocusing, the vessel becomes clearly visible.
It has been verified that although in refocused crosssectional images (e.g. figure 4(d-1)) it is difficult to recognize an exceptionally narrow vessel with a diameter of 6 µm, i.e. close to the diameter of individual erythrocytes, in the axial sections (as in figure 4(b-1)) the trajectories of such narrow vessels were still clearly visible, which confirms that the refocusing principle opens up the possibility of resolving issues associated with extremely narrow capillaries. Subsequent to refocusing and singling out the vessels as 'dark channels', the so-constructed image can then be inverted to represent the microvasculature network in the conventional manner, i.e. as bright channels against a dark background.
With regard to the sensitivity of the refocusing-based method, it has been verified that, for the purpose of visually contrasting the flow in the refocused image, it is sufficient that, during the acquisition time of the entire B-scan-stack covering the cross-section of the defocused beam, the vertical positions of the scatterers should alter by ∼λ/4 . . . λ/2, so that the constructive interference of the summed B-scans during the refocusing procedure becomes significantly perturbed. For the 'typical' OCT system discussed above, this corresponds to an axial component of particle velocities of ∼several µm/s, which should generally be sufficient to visualize blood microcirculation. The additional advantage is that the lateral resolution due to focusing can be significantly improved as compared to systems with weakly focused beams with diameters several times larger.

Conclusions
This paper describes an optimized version of the full-wave model of OCT-scan formation [16], in which the medium is represented as an ensemble of discrete scatterers. Similarly to the earlier variant [16], our version of the model uses a Gaussian illuminating beam, allowing for pronounced focusing. The OCT scans are formed by collecting the field of every spectral component scattered by discrete sub-resolution particles over the receiving aperture. However, unlike [16], where the integral (7) over the receiving aperture was evaluated numerically, in this work we present an analytical expression (10) for the result of this integration, so that only the summation of the analytically found expressions for every spectral component and every scatter requires to be performed numerically. This elimination of the numerical integration allowed us to significantly increase the operation's computation speed, by ∼10 2 −10 3 times for the same computer configuration. Therefore, simulations of 3D OCT data sets for regions containing ∼10 5 −10 6 scatterers can now be performed at intervals on the order of several hours, using a CPU Intel I7 6950X. For the typical density of cells in real biological tissue, this corresponds to tissue volumes ∼10 8 −10 9 µm 3 , which is quite sufficient for a broad range of problems related to testing and the perfection of signalprocessing methods in OCT. This model is very convenient for modeling situations with arbitrarily moving scatterers, allowing for flexible control of all main parameters characterizing the OCT setup and scattering tissue. This facilitates the study of various factors affecting OCT-based imaging of biological tissues, at a level which has been impractical in the past. The realization of similar highly controllable physical experiments would be very challenging: expensive, laborious, or even impossible in practice.
In particular, in the above sections we demonstrated the results of such highly controllable numerical experiments in relation to the feasibility of the numerical refocusing of OCT scans formed by highly-focused beams in the presence of motion of scatterers. The results demonstrated interesting prospects for the visualization of particle flows with parameters similar to the motion of erythrocytes in blood microvasculature. Although the suggested method of vessel visualization, similarly to other contrast-agent-free approaches, is based on the simulated motion of the scattering particles, the principles of the visualization described above is significantly different. The proposed computationally efficient modeling of large series of realistic OCT scans can be used to study in detail the proposed new angiographic approach, as well as potential deployment in a broad range of applications where the realization of physical experiments is challenging or even impossible in sufficiently controllable conditions.