Optical coherence tomography velocimetry based on decorrelation estimation of phasor pair ratios (DEPPAIR)

Maximilian G. O. Gräfe; Oleg Nadiarnykh; Johannes F. De Boer

doi:10.1364/BOE.10.005470

1. Introduction

In the course of many retinal diseases, such as age-related macular degeneration, diabetic retinopathy and glaucoma [1–3], the blood flow is altered and it is therefore considered to be an important indicator for abnormal developments in the eye. Optical coherence tomography (OCT) is an interferometric imaging technique to measure the coherent backscattering of photons from biological tissue. OCT angiography (OCTA) is a functional extension which has been very successful in recent years for the visualization [4] of blood flow and was even made available as commercial products [5]. However, OCTA only visualizes flow in binary maps. Because early symptoms are often related to small changes in the morphology and functionality of the vessels, quantification becomes increasingly important. Vessel density maps [6] can be used to monitor smaller vessels and capillaries but also the flow velocities and flow rates would be of great interest. Approaches such as SSADA [7] and VISTA [8] have been introduced to extract information about flow rate from modified OCTA measurements and processing.

Numerous approaches have been implemented to retrieve quantitative flow information (primarily velocity) from the OCT signal. The earliest implementation, Doppler OCT [9–11], measured the average frequency shift within A-scans. Significant work has been done to improve flow sensitivity by comparing adjacent A-scans [12] describe noise influences [13–15] and optimize estimation precision and accuracy [16–20]. However, Doppler-OCT reveals only information about the axial velocity component. In order to convert to absolute flow velocity the flow direction relative to the incident beam must be known, which can be challenging. Work has been reported to resolve the three-dimensional direction by e.g. multi-beam-illumination [21,22] or synthetic apertures in full field [23] and line field OCT [24] but at the cost of a significant increase of the complexity of the device or other limitations such as increased influence of multiply scattered photons [25]. A set of algorithms based on a statistical approach have been developed in which axial and transversal velocity are linked to the correlation between repeated A-line measurements. Examples are SSADA [7], VISTA [8], Intensity speckle decorrelation [26,27], complex decorrelation [28–30], Doppler frequency bandwidth analysis [31,32] and phase decorrelation [33]. Those algorithms can be applied using any point scanning OCT device but at the disadvantage that they often need a large number of measurements for each velocity estimation. Depending on the chosen method they can be applied to the amplitude, phase or the complex signal of OCT. The equivalence of the average estimation of some of those algorithms has been addressed in previous reports [15,17] but to our knowledge a quantitative comparison of the retrieved information content has been missing, so far. Additionally, some of those algorithms rely on the use of specific scan patterns, such as M-scans which require all measurements to be part of a trace with equidistant time intervals between individual measurements. This can be challenging when significant bulk motion should be avoided within those traces (rather than between pairs of measurements). It also has been shown that in order to make the dynamic range of the velocity estimations suitable for physiologically relevant velocities, the time interval between measurements and therefore the scan patterns require a certain variability [8,34,35].

In section 2 of this article we provide a derivation of the probability density function (PDF) of the changes of the complex OCT E-field under the influence of motion in the sample, combining the prediction of changes in amplitude, phase and the combined complex field (amplitude and phase) in one theory. This gives us the basis to compare the information for amplitude-based, phase-based and complex signal based algorithms. The Cramér-Rao lower bound (CRLB) predicts the best achievable uncertainty for these analysis bases. This provides us with a tool to evaluate the efficiency of each algorithm to extract information about flow from changes in the E-field. In order to provide best flexibility for the choice of scan patterns we restrict our analysis to the use of an ensemble of pairs of measurements where a pair is acquired sequentially with a time delay. Our analysis requires only a correlation between the measurements in one pair but not between the pairs. In section 3 the experimental system, measurements and data processing is described. Three maximum likelihood estimators (MLEs) are implemented to make the most efficient use of the available data. In section 4 the theory is compared to experimental data and the dynamic range of the estimation is investigated. In sections 5 and 6 the results are discussed and a conclusion is given. The results are an important step towards the understanding and optimization of flow estimations.

2. Theory

To model the PDF of E-field changes an approach was chosen which was previously introduced by Vakoc et al. [14]. OCT-measurements of the complex E-field are represented as vectors, which we will refer to as here phasors, in a two-dimensional plane where the x-axis represents the real and the y-axis the imaginary part. The phasors of two measurements separated by a time interval $\tau$ are defined as ${\Gamma _A}$ and ${\Gamma _B}.$ ${\Gamma _B}$ can be expressed as the sum of a fraction of ${\Gamma _A}$ and a random phasor ${\Gamma _C}$ [14,33]:

(1)$${\Gamma _B} = |\alpha |{\Gamma _A} + \sqrt {1 - \alpha {\alpha ^ \ast }} {\Gamma _C}$$

where $\alpha$ is the volume overlap integral between measurement A and B. In the case of continuous motion of the volume, $\Gamma $ will show a random walk restricted by the amplitude of the reflected field in the experiment and the average step size in this random walk is regulated by $\alpha$. A movie showing a simulation of the complex E-field backscattered from a collection of moving point scatterers is shown in Fig. 1. For the scope of this paper we chose to use the approximation of bulk flow i.e. all scatterers in the probed volume move with the same velocity. The volume overlap integral is then given by [33]:

(2)$$\alpha = {e^{i2{k_0}\delta z}}{e^{ - \frac{{\delta {x^2}}}{{{w_0}^2}}}}{e^{ - \frac{{\delta {z^2}}}{{2{l_c}^2}}}},$$

with $\delta x$ and $\delta z$ the lateral and axial displacements of the scatterers, respectively, ${w_0}$ the beam radius and ${l_c}$ the coherence length defined as ${l_c} = \sqrt {2\ln 2/\pi } \cdot {\lambda _0}^2/\Delta \lambda$ [36] with ${\lambda _0}$ the central wavelength and $\Delta \lambda$ the full width half maximum bandwidth. In Eq. (2) the phase of $\alpha$ represents the (average) Doppler shift and the absolute value represents the decorrelation by axial and transversal motion.

Fig. 1. Simulation of the complex E-field when a collection of point scatterers move (laterally) through the detection volume. Left: visualization of point scatterers. The sphere indicates where the collection efficiency of the OCT drops to 1/e². Units are displayed in w₀ (laterally) and $\sqrt 2$l_c (axially). Right: the E-field which is a coherent summation of the backscattering and collection of all point scatterers. Each step represents a motion of δx/w₀=0.5 (see Visualization 1)

Download Full Size | PDF

We address now the derivation of a joint PDF for ${\Gamma _A}$ and ${\Gamma _C}.$ For fully developed speckle patterns the two phasors in Eq. (1) follow individually a Gaussian distribution in the complex plane [14,37]. We express the real and imaginary part as coordinates in a two-dimensional plane, and the PDFs of $|\alpha |{\Gamma _A}$ and $\sqrt {1 - \alpha {\alpha ^ \ast }} {\Gamma _C}$ are given by,

(3)$${P_{|\alpha |{\Gamma _A}}}({x_1},{y_1}) = \frac{1}{{\pi \alpha {\alpha ^ \ast }}}\exp \left( { - \frac{{{x_1}^2 + {y_1}^2}}{{\alpha {\alpha^ \ast }}}} \right)$$

(4)$${P_{\sqrt {1 - \alpha {\alpha ^ \ast }} {\Gamma _C}}}(\Delta x,\Delta y) = \frac{1}{{\pi ({1 - \alpha {\alpha^ \ast }} )}}\exp \left( { - \frac{{\Delta {x^2} + \Delta {y^2}}}{{({1 - \alpha {\alpha^ \ast }} )}}} \right)$$

Here, only the decorrelation effect of the absolute value of $\alpha$ is considered. Later the joint PDF for ${\Gamma _A}$ and ${\Gamma _C}$ is transformed into the joint PDF for ${\Gamma _A}$ and ${\Gamma _B}$ where also the Doppler shift $2{k_0}\delta z$ of $\alpha$ is incorporated. Because the variables in Eq. (3) are independent of the variables in Eq. (4), the joint PDF can be expressed as the product of the individual PDFs:

(5)$$P({x_1},{y_1},\Delta x,\Delta y) = \frac{1}{{\pi \alpha {\alpha ^ \ast }}}\exp \left( { - \frac{{{x_1}^2 + {y_1}^2}}{{\alpha {\alpha^ \ast }}}} \right)\frac{1}{{\pi ({1 - \alpha {\alpha^ \ast }} )}}\exp \left( { - \frac{{\Delta {x^2} + \Delta {y^2}}}{{({1 - \alpha {\alpha^ \ast }} )}}} \right).$$

In order to turn the PDF for fractions of the phasors into a PDF for the phasors with full length, the substitutions ${x_A} = {x_1}/|\alpha |,$ ${y_A} = {y_1}/|\alpha |,$ ${x_C} = \Delta x/\sqrt {1 - \alpha {\alpha ^ \ast }}$ and ${y_C} = \Delta y/\sqrt {1 - \alpha {\alpha ^ \ast }}$ are made in Eq. (5):

(6)$$P({x_A},{y_A},{x_C},{y_C}) = \frac{1}{{{\pi ^2}}}{e^{ - {x_A}^2 - {y_A}^2}}{e^{ - {x_C}^2 - {y_C}^2}}.$$

In practice, ${\Gamma _C}$ cannot be measured directly but only the pair ${\Gamma _A}$ and ${\Gamma _B}$. Therefore, the joint PDF for the latter phasor pair is derived with the substitutions ${x_{B^{\prime}}} = {x_A}|\alpha |+ {x_C}\sqrt {1 - \alpha {\alpha ^ \ast }}$ and ${y_{B^{\prime}}} = {y_A}|\alpha |+ {y_C}\sqrt {1 - \alpha {\alpha ^ \ast }}$:

(7)$$P({x_A},{y_A},{x_{B^{\prime}}},{y_{B^{\prime}}}) = \frac{{{e^{ - {x_A}^2 - {y_A}^2}}}}{{{\pi ^2}({1 - \alpha {\alpha^ \ast }} )}}\exp \left[ { - {{\left( {\frac{{{x_{B^{\prime}}} - {x_A}|\alpha |}}{{({1 - \alpha {\alpha^ \ast }} )}}} \right)}^2} - {{\left( {\frac{{{y_{B^{\prime}}} - {y_A}|\alpha |}}{{({1 - \alpha {\alpha^ \ast }} )}}} \right)}^2}} \right].$$

Equation (7) depends only on the absolute value of the overlap integral $\alpha .$ Next, we incorporate the Doppler shift or phase term of Eq. (2). The effect of a factor ${e^{i\theta }}$ in the complex plane is a pure rotation around the origin by an angle $\theta$. To translate this effect to real valued vectors, which we use to describe the measurements ${\Gamma _A}$ and ${\Gamma _B},$ a rotation in a two-dimensional plane is generated by a rotation matrix:

(8)$$\left( {\begin{array}{{c}} {{x_{B^{\prime}}}}\\ {{y_{B^{\prime}}}} \end{array}} \right) = R \cdot \left( {\begin{array}{{c}} {{x_B}}\\ {{y_B}} \end{array}} \right) = \left[ {\begin{array}{{cc}} {\cos \theta }&{ - \sin \theta }\\ {\sin \theta }&{\cos \theta } \end{array}} \right] \cdot \left( {\begin{array}{{c}} {{x_B}}\\ {{y_B}} \end{array}} \right) = \left( {\begin{array}{{c}} {{x_B} \cdot \cos \theta - {y_B} \cdot \sin \theta }\\ {{x_B} \cdot \sin \theta + {y_B} \cdot \cos \theta } \end{array}} \right).$$

Note that the rotation is only applied to ${\Gamma _B}$ but not ${\Gamma _A}$ because ${\Gamma _A}$ represents the phasor before the phase shift. Replacing ${x_{B^{\prime}}}$ and ${y_{B^{\prime}}}$ in Eq. (8) and changing to polar coordinates yields:

(9)$$\begin{aligned} P({A_A},{\varphi _A},{A_B},{\varphi _B}) &= \frac{{{A_A}{A_B}{e^{ - {A_A}^2}}}}{{{\pi ^2}({1 - \alpha {\alpha^ \ast }} )}}\\ & \times \exp \left[ { - {{\left( {\frac{{{A_B}\cos ({\varphi_B} + \theta ) - |\alpha |{A_A}\cos {\varphi_A}}}{{({1 - \alpha {\alpha^ \ast }} )}}} \right)}^2}} \right]\\ & \times \exp \left[ { - {{\left( {\frac{{{A_B}\sin ({\varphi_B} + \theta ) - |\alpha |{A_A}\sin {\varphi_A}}}{{({1 - \alpha {\alpha^ \ast }} )}}} \right)}^2}} \right] \end{aligned}$$

with ${A_j}$ the radial and ${\varphi _j}$ the angular components of ${\Gamma _j}$ $(j = A,B)$ representing field amplitude and phase of the measured phasors, respectively.

Equation (9) now contains all information about the PDF for a measured pair of phasors. For practical applications not all information is of interest or even meaningful. Therefore, in the following paragraph parametrizations are introduced which help to derive the marginal PDFs that give information about directly accessible values of phasor changes. First, not the actual values of the phase but only the phase change between measurements can be used to determine motion. Thus, the variable for the phase difference $\varphi = {\varphi _B} - {\varphi _A}$ is introduced and Eq. (9) is integrated over ${\varphi _A}$:

(10)$$P({A_A},{A_B},\varphi ) = \frac{{2{A_A}{A_B}{e^{ - {A_A}^2}}}}{{\pi ({1 - \alpha {\alpha^ \ast }} )}}\exp \left[ { - \frac{{{A_B}^2 + {\alpha^2}{A_A}^2 - 2|\alpha |{A_A}{A_B}\cos (\varphi + \theta )}}{{({1 - \alpha {\alpha^ \ast }} )}}} \right].$$

To express changes of the amplitude, the coordinate transform $p = {A_A}{A_B},$ $q = {A_B}/{A_A}$ was chosen. This allows us to express phasor changes as a conjugate product ${\Gamma _B}{\Gamma _A}^ \ast{=} p{e^{i\varphi }}$ or a ratio ${\Gamma _B}/{\Gamma _A} = q{e^{i\varphi }}.$ Then Eq. (10) becomes:

(11)$$P(p,\;q,\;\varphi ) = \frac{{{e^{ - p/q}}}}{{\pi ({1 - \alpha {\alpha^ \ast }} )}}\frac{p}{q}\exp \left( { - \frac{{pq + |\alpha |\frac{p}{q} - 2|\alpha |p\cos ({\varphi + \theta } )}}{{({1 - \alpha {\alpha^ \ast }} )}}} \right).$$

From Eq. (11) the marginal PDFs for $P(p,\varphi )$ and $P(q,\varphi )$ can be derived by integrating over q and $p,$ respectively:

(12)$$P(p,\varphi ) = \frac{{2p}}{{\pi ({1 - \alpha {\alpha^ \ast }} )}}\exp \left( {\frac{{2p|\alpha |\cos ({\varphi + \theta } )}}{{({1 - \alpha {\alpha^ \ast }} )}}} \right){K_0}\left( {\frac{{2p}}{{({1 - \alpha {\alpha^ \ast }} )}}} \right)$$

(13)$$P(q,\varphi ) = \frac{{q({1 - \alpha {\alpha^ \ast }} )}}{{\pi {{({{q^2} + 1 - 2q|\alpha |\cos ({\varphi + \theta } )} )}^2}}}$$

with ${K_0}$ the modified Bessel function of the second kind.

We can verify that the here derived PDFs are in agreement with other theories on speckle decorrelation [27–30,38] which have in common that from the autocorrelation of a time trace of the electrical field E the overlap integral $\left\langle {{E_{l + 1}}{E_l}^ \ast } \right\rangle = \alpha$ can be calculated, with l being the index for the l-th measurement. The same expectation value of ${E_{l + 1}}{E_l}^ \ast $ can be calculated using Eq. (12), giving the same result:

(14)$$\left\langle {{E_{l + 1}}{E_l}^ \ast } \right\rangle = \int\!\!\!\int {p{e^{i\varphi }}P(p,\varphi )} dpd\varphi = \alpha .$$

For methods to estimate $\alpha$ the PDF $P(q,\varphi )$ is used since it has an easier practical implementation: ${A_A}$ and ${A_B}$ are currently modeled such that they show a Rayleigh distribution with a mean of $\sqrt \pi /2$. In an actual measurement the amplitude of the E-field will scale with the average reflectivity R of the scattering particles: $|E |= \sqrt R A.$ Therefore, the reflectivity will also influence the product p. Either R has to be estimated and the measurements have to be scaled by its inverse, or it must be included in the PDF $P(p,\varphi ).$ This issue does not exist for $P(q,\varphi )$ because q is insensitive to any scaling factors in the field amplitude. A second reason is that the numerical evaluation of $P(q,\varphi )$ is more stable than $P(p,\varphi )$ in further calculations due to the Bessel function.

The marginal PDFs for pure phase decorrelation and pure amplitude decorrelation can be derived from Eq. (13) by the integration over q and $\varphi ,$ respectively:

(15)$$P(\phi ) = \frac{{1 - \alpha {\alpha ^ \ast }}}{{2\pi (1 - \alpha {\alpha ^ \ast }{{\cos }^2}\phi )}}\left[ {1 + \frac{{|\alpha |\cos \phi }}{{\sqrt {1 - \alpha {\alpha^ \ast }{{\cos }^2}\phi } }}({\pi - {{\cos }^{ - 1}}[{|\alpha |\cos \phi } ]} )} \right],$$

(16)$$P(q) = \frac{{2q({1 + {q^2}} )({1 - \alpha {\alpha^ \ast }} )}}{{\sqrt {{{({1 + {q^2} + 2q|\alpha |} )}^3}{{({1 + {q^2} - 2q|\alpha |} )}^3}} }},$$

with $\phi = \varphi + \theta$. Equation (15) is identical to the expression derived earlier [14,33]. We can express $\alpha$ in the parameters $\beta ,$ representing lateral and axial motion and $\theta ,$ representing the Doppler shift:

(17)$$\alpha = \exp ({i\theta - {\beta^2}} ).$$

β can be understood as a term summarizing all influences which can lead to a decorrelation of the OCT-measurement. This can for instance be bulk motion, particle flow as depicted in the movie of Fig. 1, or even brownian motion. For bulk motion where all scatterers in the probed volume move with the same velocity along a specific axis (as assumed for Eq. (2)), $\beta$ can be identified with an axial $(\beta = \delta z/\sqrt 2 {l_c})$ or transversal $(\beta = \delta x/{w_0})$ velocity normalized on ${w_0}$ and ${l_c},$ respectively. The presented theory can also be used to describe more complex scenarios. It has been shown that diffusion [28–30,38], clutter [28,38] (the presence of moving and stationary particles in the same detection volume) and flow gradients [27,39] can be included in the parameter. Also different focusing and light collection geometries [32] can be included, affecting $\alpha$ but the PDFs for $P(q,\varphi ),\;P(q)$ and $P(\varphi )$ remain unchanged. We focus in our analysis on the estimation of $\alpha$ and the conversion to $\beta ,$ which can be identified as relative motion under the above mentioned conditions of bulk flow.

For an illustration of the behavior of phasor ratios, phase differences and amplitude ratios for bulk flow, the PDFs of Eq. (13), (15) and (16) are plotted in Fig. 2 for different values of relative motion $\beta .$ As can be seen in the graphs, the PDFs broaden for an increase in $\beta .$ The phase difference PDF in Fig. 2a is a circular symmetrical distribution around average Doppler shift $\theta$ (here the axial motion was zero, so Doppler shift $\theta = 0rad$) over the interval ${\pm} \pi$. The amplitude ratio PDF in 2b is an asymmetrical distribution between 0 and ∞ with half of the integral below $q = 1$ and the other half above $q = 1.$ Fig. 2c, 2d and 2e show the combined phase difference and amplitude ratio distributions for the same values of $\beta$ as in 2a and 2b. All numerical implementations were evaluated in Matlab 2014a (Mathworks Inc.).

Fig. 2. Plots of PDFs for different normalized motion β. (a) marginal PDF for phase differences; (b) marginal PDF for amplitude ratios; (c), (d) and (e) show the joint PDF for phase differences and amplitude ratios for the same values of β as in (a) and (b). The maximum of the color scale (“parula”, Mathworks Inc.) was set to the maximum of each PDF. When the motion becomes stronger, β becomes larger and therefore the phasors decorrelate more. This behavior is expressed as a broadening of the PDFs.

Download Full Size | PDF

With an exact expression for the PDFs describing the measurements, the calculation of the CRLB [40] is straight forward. The CRLB predicts the smallest achievable uncertainty (which we quantify here as standard deviation) with which a parameter can be estimated from a measurement. It depends on the sensitivity of the PDF with respect to a change of the parameter, in our case the normalized motion $\beta ,$ but is independent of the estimation method. Therefore we will use the terms phasor-based, phase-based and amplitude-based methods for the CRLBs which belong to the PDFs of Eq. (13), (15) and (16), respectively. The CRLBs for those methods are plotted in Fig. 3a as a function of $\beta .$ Because the phasor-based methods use the full information available by the OCT measurement, it is expected that the respective CRLB shows the lowest uncertainty. This is also confirmed by the black line in the graph. Very similar behavior is shown by the CRLB of the phase-based methods. In contrast to that, the amplitude-based analysis performs significantly worse and the uncertainty increases nonlinearly, especially for high velocities. The estimation uncertainty decreases with an increase in the number of measurements N (in our case pairs of measurements) which are collected per ensemble to perform each estimation increase. The variance is given as ${\sigma ^2} = {({NI} )^{ - 1}}$ with I the Fisher-Information [40,41].

(18)$$I = E\left[ {\frac{{\partial \log P(x|\alpha )}}{{\partial \alpha }}} \right],$$

with P(x) the PDF of the variable x and the function E […] in Eq. (18) the expectation value. This leads to the proportionality $\textrm{CRLB} \propto 1/\sqrt N$. This relation can also be reversed in order to calculate a factor describing how many measurements are necessary to reach the same uncertainty with the phase-based and amplitude-based analysis, relative to the phasor based analysis. The phasor-based analysis was chosen as a reference since it resulted in the lowest uncertainty. Fig. 3b shows the square of the ratio of the CRLBs of the phase-based and amplitude-based analysis with respect to the phasor-based analysis, which expresses how many more measurements are needed by these methods as a function of β to reach the same uncertainty as the phasor-based approach. The phase-based analysis performs nearly as good as the phasor-based analysis with ratios between two (small velocities) and one (high velocities). To reach the same uncertainty the amplitude-based analysis requires an order of magnitude more measurements than what is necessary with the phase- or phasor-based analysis.

Fig. 3. (a) CRLBs for phasor-based (black), purely phase-based (red) and purely amplitude-based (blue) decorrelation estimation methods for N = 1 measurements. The phasor-based methods have the lowest CRLB. The CRLBs decrease proportional to $\sqrt N$. b) The square of the ratio of the CRLB of the phase-based and amplitude-based methods with respect to the phasor-based method, which expresses how many more measurements are needed by these methods as a function of β to reach the same uncertainty as the phasor-based approach. The purely phase-based analysis performs nearly as well as the phasor-based but the purely amplitude-based needs an order of magnitude more measurements especially in the high velocities rang ($\beta \approx 1$) to reach the same uncertainty. The black horizontal line indicates a ratio of one.

Download Full Size | PDF

3. Experiments

3.1 Setup

A handheld OCT-system was used as described by Nadiarnykh et al. [42]. Briefly, a swept laser source (Axsun Inc.) with central wavelength at 1060nm and a repetition rate of 100 kHz was used with a fiber based system. The sample arm included a handheld piece which is designed to be held approximately 1cm above the cornea of a human eye with a beam diameter of 1.2mm, while the person is facing upwards. Instead of an actual eye, a model eye was used consisting of a lens (AC-127-019-C, Thorlabs Inc.) and a retina phantom as it was previously presented in [33]. A schematic drawing of the measurement configuration and a structural B-scan of the phantom is shown in Fig. 4. The phantom consisted of a bulk scattering material (TiO₂ particles in cured silicone) with an embedded glass vessel (inner diameter 150µm). A 0.5% aqueous intralipid solution was pumped by a syringe pump (NE-4000, New Era Pump Systems Inc.) through the vessel to create flow of a scattering fluid. Flow was measured in two different configurations. In the first configuration, the beam was translated over the static scattering medium, mimicking bulk flow. In the second configuration, flow was measured in the vessel. Data acquisition was performed in M-scan mode, for each lateral location 2025 A-lines were recorded, before laterally displacing the beam by 1.1 µm for the next M-scan. 2000 M-scans were acquired to create a B-scan. This generates a field of view of 2.2 mm and a acquisition time of 40.5 s per B-scan. After a phantom measurement, a noise measurement was performed for fixed pattern background removal by blocking the sample arm. The signal-to-noise ratio (SNR) in the phantom was between 20 and 30dB.

Fig. 4. a) schematic drawing of the setup. The scan module of an OCT system emitting a collimated beam was mounted above a lens and retina phantom. The phantom consist of a glass vessel embedded in a scattering medium. The capillary was connected to a syringe pump which pumped an aqueous intralipid with a constant flow rate. b) Structural B-scan of the phantom with indicated location of capillary and area of static scattering medium.

Download Full Size | PDF

3.2 Processing

The first A-line processing steps were background subtraction for fixed pattern removal, spectral shaping with a tapered cosine window and Fourier transform. Secondly, individual measurements in the M-scans were combined into phasor pairs. For the application of the MLEs to experimental data, the phasor pairs were first grouped by the $\beta$ they represent and then sorted in 2024/N ensembles with N phasor pairs per ensemble. Thus, the average and standard deviation of $\beta$ could be calculated for each group of ensembles. For the static medium, the phasor pairs were created from two A-lines at the same depth location separated by the same step size (a multiple of 1.1 µm) to represent a group of ensembles with the same $\beta .$ For measurements in intralipid each M-scan is the measurement of 2025 A-lines at the same location resulting in 2024 phasor pairs for each depth location which represents one group of ensembles.

Each phasor pair was corrected for bulk motion in the following way: the conjugate product was summed along each A-line in the bulk phantom: $E = \sum\nolimits_k {{E_k}{E_{k - 1}}^ \ast } .$ The angle of this sum was extracted and used as bulk phase shift for compensation. The glass vessel was embedded below a layer of bulk scattering medium. The bulk phase shift which was extracted in an A-line from the medium above the intralipid solution, was also applied to measurements of the intralipid. For the estimation of $\beta ,$ three MLEs were implemented according to:

(19)$${\beta _{MLE}} = \mathop {\arg \max }\limits_\beta \sum\limits_j {\log P({q_j},{\varphi _j}|\beta )} ,$$

where j={1,2,…, N} is the index of the j-th measurement. Equation (19) is the definition for the phasor-based MLE using the PDF of Eq. (13). For the phase-based and amplitude based MLEs it must be replaced by the PDFs in Eq. (15) and (16), respectively.

In the following paragraph we describe the implementation of the individual MLEs. The phase-based MLE was described in [33]. In brief: Before applying the MLE, the average axial Doppler shift $\theta$ in the phasor pair distribution was removed by circular averaging of the phase differences and afterwards the phase distribution was shifted circularly to bring the mean phase difference to zero. PDFs were created for $\beta$ ranging from zero to a maximum ${\beta _{\max }}$ of 2.5 in 200 steps. These PDFs were used to compute the joint likelihood function. In this implementation we added a spline interpolation of the joint likelihood function for a more accurate determination of the maximum likelihood. For the phasor-based MLE the estimation of the average axial Doppler shift is included in the MLE. $\beta$ and $\theta$ are estimated in two steps. In the first step, Eq. (19) is evaluated over a coarse grid of 200 samples for $\beta$ (0 to ${\beta _{\max }}$) and $\theta$ (-π to +π), respectively. It was implemented as a search for zero in the first derivative with respect to $\beta$ of the joint likelihood function (the function has only one extremum):

(20)$$- \frac{{2\textrm{N}{e^{ - {\beta _{MLE}}^2}}}}{{1 - {e^{ - 2{\beta _{MLE}}^2}}}} + \sum\limits_j {\frac{{4{q_j}\cos ({\varphi _j} + \theta )}}{{1 + {q_j}^2 - 2{q_j}{e^{ - {\beta _{MLE}}^2}}\cos ({\varphi _j} + \theta )}}} = 0.$$

The numerical evaluation of Eq. (20) on a grid provides pairs of $\theta$ and $\beta$ for which Eq. (20) is close to zeros. These pairs trace out a path in a two dimensional ($\theta$, $\beta$) space (see Fig. 5). Only for the smallest $\beta$ on this path a unique value for $\theta$ can be found. Finding this unique combination is equivalent to finding parameters which maximize the likelihood function and therefore yield ${\beta _{MLE}}$ and ${\theta _{MLE}}.$ In a second step the accuracy of the estimation is improved: the four combinations of $\theta$ and $\beta$ around the zero-crossing for ${\beta _{MLE}}$ and ${\theta _{MLE}}$ are used as boundaries to compute a second grid with 200 steps for each parameter.

Fig. 5. Illustration of the derivative of the joint likelihood which is used to find combinations of β and θ fulfilling Eq. (20) for one measured ensemble at β=1.2

Download Full Size | PDF

The implementation of the amplitude-based MLE is simpler in comparison to the phase-based and phasor-based MLEs because inherent to its definition the amplitude-based MLE is independent of $\varphi$ and $\theta .$ Similar to the phasor-based MLE it is also evaluated in two steps, first with a coarse search for ${\beta _{MLE}}$ and secondly with a finer grid (only two boundaries) for better accuracy.

4. Results

Before the MLEs were applied to measurement of intralipid flow in the glass vessel, the linear dependency of the estimation of $\beta$ on axial and lateral motion was verified and the dynamic range was investigated. This was done by controlled axial and lateral bulk motion of the static scattering medium. To mimic axial motion, a B-scan was created by selecting one A-scan out of every M-scan. Axial motion was generated digitally – a copy of the complex valued B-scan was Fourier transformed to spectral domain, a phase ramp $2\pi s\kappa /{\rm K}$ was applied with s the shift in pixels, $\kappa$ the wavenumber index and ${\rm K}$ the number of samples per A-line. The shifted image was back-transformed to depth domain and phasor pairs were assembled from corresponding pixels in the original and shifted copy of the B-scan. For lateral motion phasor pairs were generated from adjacent A-lines, mimicking a lateral translation of the sample. Due to a high degree of beam overlap, motion in fractions of the lateral resolution could be mimicked. The results for axial and lateral motion estimation are shown in Fig. 6 for the phasor-based MLE for different numbers (N) of phasor pairs per ensemble. The fit of a line shows the linear regime in both graphs (a and b). The upper limit of the linear regime, which defines the dynamic range depends on N. In Fig. 6(a) where the axial motion was mimicked, the linear regime extends on the lower end to zero. In Fig. 6(b) where the lateral motion was mimicked using adjacent A-lines, an offset can be noticed. Here, a line cannot be fitted which crosses the origin. Also close to the origin the measurements deviate from a linear slope. We attribute this deviation to the influence of shot noise and beam position jitter originating from the galvanometer scanners. The lowest measurable $\beta = 0.098$ represents the lower bound of the dynamic range and was obtained with a stationary beam (except for galvanometer scan jitter).

Fig. 6. Verification of linear estimation behavior and dynamic range of motion estimation with phasor-based MLE for (a) axial motion and (b) lateral motion. The vertical axis shows the estimated parameter β. Lines were fitted to the linear parts of the graphs. The upper limit of the dynamic range (linear regime) is dependent on the number of phasor pairs (N) per ensemble. The lower limit is independent on N and determined by noise sources.

Download Full Size | PDF

The slopes in Fig. 6 can be used to determine the coherence length and beam radius. According to Eq. (2) and (17) axial shifts (6a) are estimated as $\beta = \delta z/\sqrt 2 {l_c}.$ Assuming a Gaussian spectrum of the source the coherence length l_c was 14.6 µm. A fit to the slope in Fig. 6(a) resulted in a coherence length l_c of 13.3 µm. Lateral shifts (6b) are estimated as $\beta = \delta x/{w_0}$. The theoretically calculated diffraction limited beam radius in the sample is 11 µm, given by the beam radius at the last lens before the phantom and the focal length of the lens. The actual beam radius derived from the slope in 6b is 15 µm.

The estimation uncertainties of $\beta$ and the agreement with the CRLBs was verified by a measurement with intralipid in the glass vessel. The solution was pumped with a flow rate of 435µl/min. In every M-scan only measurements inside the vessel were used and every location along the depth inside the vessel represents a different $\beta .$ Thus, a group of ensembles was created from every depth location from every M-scan. The standard deviation of $\beta$ was then plotted as a function of the estimation of $\beta$ averaged over the ensembles. To reduce effects of multiple scattering only the upper half of the vessel was analyzed. The results are shown in Fig. 7(a) for an ensemble size of N = 200. The estimations of the measurements are displayed as dots while the CRLBs for the corresponding methods (phasor-based, phase-based, amplitude-based) are displayed as solid lines in the same color as the corresponding MLEs. Figure 7(a) shows that the measurements follow approximately the predicted CRLBs. To increase the number of ensembles from 10 to 90, binning was applied. For the binning the results were sorted according to the average estimations of $\beta$ and consecutive results were again averaged. This happened by a binning factor of 9 and is shown in Fig. 7(b).

Fig. 7. Uncertainty estimations as standard deviations for measurements (dots) of decorrelation for flow of intralipid through a glass vessel. The corresponding CRLBs are displayed as solid lines of the same color of the respective measurements. (a) all results for the uncertainty estimation with N = 200, (b) results from (a) binned by a binning factor of 9 for a more concise display, (c) effect on standard deviation when randomness of phasor pairs is increased (N = 50), (d) display of outliers in amplitude-based estimation when PDF is sampled insufficiently (N = 50).

Download Full Size | PDF

Figure 7(a) and 7(b) show that all three estimation methods achieve the theoretically predicted uncertainty for large ensembles (N = 200). For measurements in vivo, an ensemble of 200 measurements requires a significant amount of time. Therefore, we investigated the effect of using a smaller number of measurements in an ensemble. Estimations were also done for N = 50 in 7(c) and 7(d). Figure 7(c) shows the effect of selecting 50 phasor pairs created from 51 consecutive A-lines in an M-scan (Shuffle OFF), or selecting 50 phasor pairs out of consecutive A-lines randomly selected out of the 2025 A-lines in an M-scan (Shuffle ON). The effect of shuffling is shown in Fig. 7(c) as an example for the phasor-based estimations. Red dots represent estimations from phasor pairs as consecutive neighbors (no shuffling), black dots represent estimations from shuffled phasor pairs. Figure 7(c) demonstrates that the estimations from shuffled phasor pairs achieve lower uncertainty. This example demonstrates that phasor pairs constructed from 51 consecutive M-scans are not completely statistically independent. For all other sub-figures in Fig. 7 we used Shuffle ON.

While the phasor- and phase-based MLE achieve the uncertainties predicted by the CRLBs, large artifacts were found for the amplitude-based MLE in 7d (delineated in black). It was found that those artifacts have a higher chance of appearance for large $\beta$ and become reduced if N is increased. They exhibit a patch-like organization, indicated by the black ellipse. Each plotted point in the graph represents one group of estimations. In the first patch (bottom left in the ellipse) every group of ensembles contained one ensemble which had a constellation of measurements that resulted in a joint likelihood function for an unexpectedly high $\beta .$ In the second patch two outliers can be found in the group of ensembles, in the third patch three, etc..

If N is increased, more measurements are included in the ensemble, the expected PDF seems to be better resembled by the measurements and the joint likelihood function peaks again at a more reasonable $\beta .$ Those outliers of estimates for $\beta$ increase mean and standard deviation of the whole group of ensembles.

5. Discussion

In this study a new theory was developed to compare quantitatively the precision of flow estimation by analyzing the information content of the estimation of decorrelation and related motion by three categories of methods: (complex) phasor-based, phase-based and amplitude-based analysis. The PDF’s describing the decorrelation between two measurements were derived for these three methods. Based on the PDF’s, the CRLB was calculated, and experimental data was compared to the CRLB prediction. The most striking result is that the information on decorrelation is dominated by information from the phase over the amplitude of the OCT-signal. This is especially the case for the upper range of velocities / motion ($\beta$ close to 1). Because complex methods combine both information from phase and amplitude they have the best theoretically achievable precision. A maximum likelihood estimator was implemented for each of the above mentioned categories and was applied to measurements of a flowing intralipid solution. An excellent agreement of the experimental measurements and theory was found for the complex phasor-based and phase-based methods. For the amplitude-based MLE we found artifacts for the combination of strong decorrelation ($\beta$ close to 1) and small numbers of samples per ensemble. The dynamic range of the flow estimation was investigated and it was shown that at the upper end it depended on the number of samples per ensemble. For the numbers used here (N = 25 to N = 200) an upper range of $\beta$ of approximately 0.95 to 1.15 was found. The lower bound was independent of the number of samples and only determined by noise sources and was in the measurements here at approximately 0.1.

The dynamic range of the velocity calculation is an important aspect for applications to measure e.g. blood flow. As presented here, $\beta$ is for transversal motion a function of ${w_0}.$ the flow velocity v and the time delay $\tau$ between repeated measurements, $\beta = v \cdot \tau /{w_0}.$ The adjustment of the sensitivity of OCT to different flow velocities has been shown before through the adjustment of $\tau$ [8,35,43] This can be realized by the adjustment of scan patterns. A typical value for ${w_0}$ is 10µm. Under this assumption Table 1 shows three examples how $\tau$ can be designed to cover different flow regimes. This can be realized by the adjustment of scanning patterns. Alternatively, also scanning with multiple beam diameters might be possible as it has been reported before [44].

Table 1. Conversion of β under the assumption of a beam radius of 10µm

View Table

Figure 8 illustrates the combination of dynamic ranges as mentioned in Table 1 with the relative error (standard deviation divided by expectation value) for a number of N = 25 phasor pairs. One interesting aspect for the combination of multiple time intervals to cover a larger range of velocity estimations is that for the phase-based and phasor-based methods the statistical error is nearly linearly dependent on the expectation value of $\beta .$ Thus, for a combined range the error would be an approximately fixed fraction of $\beta$ over three decades of flow velocities. The amplitude-based analysis on the other hand results in strong discontinuities in the regions where two ranges are merged.

Fig. 8. Illustration of the combination of multiple dynamic ranges and their relative error (standard deviation divided by expectation value) for phasor-based (solid black), phase-based (red) and amplitude-based (blue) analysis for a number of N = 25 phasor pairs. The vertical dashed lines indicate point in the velocity regime where the individual dynamic ranges are merged. The phasor-based analysis delivers the most constant relative error making the standard deviation a fixed percentage of the expectation value

Download Full Size | PDF

The limiting noise sources for the lower end of the dynamic range are shot noise and motion artifacts. Shot noise can be related to the SNR of a measurement and has been shown to broaden the distributions of Doppler phase shifts [13] and Doppler frequencies [15] and the implementation of both shot noise and decorrelation noise for optimal average Doppler frequency estimation has been demonstrated [18]. In the future, an incorporation of the influence of shot noise into the model described here and estimation of decorrelation would be beneficial because it is expected that it can bias the estimation and influence the uncertainty. A modeling of bulk motion artifacts is more complicated because it can come from different sources e.g. positioning error (jitter in scanner position) and sample motion and will highly depend on the technical implementation and conditions of the measurement. OCT is often applied in ophthalmology where it is used for minimally invasive monitoring and diagnosis of diseases and where sample fixation is very limited. Eye motion as micro-saccades, drift and tremor [45,46] can significantly influence the repeatability of the measurement of the same location. Limiting the data analysis to ensembles of measurements where only pairs of measurements are correlated (as implemented here) rather than the requirement of traces of measurement might help to overcome this limitation.

A fundamental assumption for the theory presented here is that fully developed speckle are measured. A consequence is that measurements follow a random, two-dimensional Gaussian distribution in the complex plane and therefore the amplitude follows a Rayleigh distribution. It has been shown in the past [47] that this is not generally fulfilled in biological tissue. Then the amplitude distribution can be approximated by a K-distribution and different cell and tissue types can be identified by the fitting parameters of this distribution. This changes the speckle pattern and it can be expected that the distributions of phases and phase changes may be altered by this, as well. Since for this study only small scattering particles were used (TiO₂, intralipid) this effect did not play a role but in future an investigation of this effect to blood flow measurements might be important. The backscattering of blood cells is different from TiO₂ and intralipid and depends on the orientation of the cells [48,49]. Strong forward scattering can also influence the speckle pattern of an A-line profile. Light which is scattered in tissue below a blood vessel can carry amplitude and phase changes due to forward scattering in the blood vessel above.. In this theory the effects of multiple scattering was not included. Those can be important in the explanation of shadowing artifacts known in OCT angiography [50] and the consideration might be important for accurate velocity and flow rate quantification when multiple vessels are within one A-line profile. However, previous reports showed the applicability of similar conditions to blood flow measurements [7,28,39,51] which might indicate that estimations are still valid or only small changes are necessary.

6. Conclusion

A PDF was derived to describe the changes of the complex OCT E-field between individual measurements under a decorrelation influence such as axial and lateral flow. From this also the PDFs for pure phase and pure amplitude changes were derived. Using those PDFs the CRLBs for the individual categories of analysis methods were derived, predicting the uncertainty of the estimations of decorrelation under the condition of a given number of measurements. MLEs have been implemented for each methods category and were used to verify the theoretical predictions of the CRLBs. The most interesting result is that the estimation uncertainty is almost completely determined by the phase of the measurement while the respective MLE significantly outperforms the amplitude-based MLE. The method requires correlation between measurements in ensembles of pairs but does not require correlation between the individual pairs which provides high flexibility in the choice of scan patterns.

Funding

Stichting voor de Technische Wetenschappen (12822); Nederlandse Organisatie voor Wetenschappelijk Onderzoek; Ministerie van Economische Zaken; Health~Holland framework; Horizon 2020 Framework Programme (654148); Heidelberg Engineering.

Disclosures

The authors declare that there are no conflicts of interest related to this article.

References

1. R. Ehrlich, A. Harris, N. S. Kheradiya, D. M. Winston, T. A. Ciulla, and B. Wirostko, “Age-related macular degeneration and the aging eye,” Clin. Interventions Aging 3, 473–482 (2008). [CrossRef]

2. V. Patel, S. Rassam, R. Newsom, J. Wiek, and E. Kohner, “Retinal Blood-Flow in Diabetic-Retinopathy,” Br. Med. J. 305(6855), 678–683 (1992). [CrossRef]

3. M. Emre, S. Orgul, K. Gugleta, and J. Flammer, “Ocular blood flow alteration in glaucoma is related to systemic vascular dysregulation,” Br. J. Ophthalmol. 88(5), 662–666 (2004). [CrossRef]

4. C. L. Chen and R. K. Wang, “Optical coherence tomography based angiography [Invited],” Biomed. Opt. Express 8(2), 1056–1082 (2017). [CrossRef]

5. X. X. Li, W. Wu, H. Zhou, J. J. Deng, M. Y. Zhao, T. W. Qian, C. Yan, X. Xu, and S. Q. Yu, “A quantitative comparison of five optical coherence tomography angiography systems in clinical performance,” Int. J. Ophthalmol. (Engl. Ed.) 11(11), 1784–1795 (2018). [CrossRef]

6. A. Rabiolo, F. Gelormini, R. Sacconi, M. V. Cicinelli, G. Triolo, P. Bettin, K. Nouri-Mahdavi, F. Bandello, and G. Querques, “Comparison of methods to quantify macular and peripapillary vessel density in optical coherence tomography angiography,” PLoS One 13(10), e0205773 (2018). [CrossRef]

7. J. Tokayer, Y. Jia, A. H. Dhalla, and D. Huang, “Blood flow velocity quantification using split-spectrum amplitude-decorrelation angiography with optical coherence tomography,” Biomed. Opt. Express 4(10), 1909–1924 (2013). [CrossRef]

8. S. B. Ploner, E. M. Moult, W. Choi, N. K. Waheed, B. Lee, E. A. Novais, E. D. Cole, B. Potsaid, L. Husvogt, J. Schottenhamml, A. Maier, P. J. Rosenfeld, J. S. Duker, J. Hornegger, and J. G. Fujimoto, “Toward quantitative optical coherence tomography angiography: Visualizing Blood Flow Speeds in Ocular Pathology Using Variable Interscan Time Analysis,” Retina 36(Suppl 1), S118–S126 (2016). [CrossRef]

9. X. J. Wang, T. E. Milner, and J. S. Nelson, “Characterization of Fluid-Flow Velocity by Optical Doppler Tomography,” Opt. Lett. 20(11), 1337–1339 (1995). [CrossRef]

10. Z. P. Chen, T. E. Milner, D. Dave, and J. S. Nelson, “Optical Doppler tomographic imaging of fluid flow velocity in highly scattering media,” Opt. Lett. 22(1), 64–66 (1997). [CrossRef]

11. J. A. Izatt, M. D. Kulkami, S. Yazdanfar, J. K. Barton, and A. J. Welch, “In vivo bidirectional color Doppler flow imaging of picoliter blood volumes using optical coherence tomograghy,” Opt. Lett. 22(18), 1439–1441 (1997). [CrossRef]

12. Y. H. Zhao, Z. P. Chen, C. Saxer, S. H. Xiang, J. F. de Boer, and J. S. Nelson, “Phase-resolved optical coherence tomography and optical Doppler tomography for imaging blood flow in human skin with fast scanning speed and high velocity sensitivity,” Opt. Lett. 25(2), 114–116 (2000). [CrossRef]

13. B. H. Park, M. C. Pierce, B. Cense, S. H. Yun, M. Mujat, G. J. Tearney, B. E. Bouma, and J. F. de Boer, “Real-time fiber-based multi-functional spectral-domain optical coherence tomography at 1.3 mu m,” Opt. Express 13(11), 3931–3944 (2005). [CrossRef]

14. B. J. Vakoc, G. J. Tearney, and B. E. Bouma, “Statistical Properties of Phase-Decorrelation in Phase-Resolved Doppler Optical Coherence Tomography,” IEEE Trans. Med. Imaging 28(6), 814–821 (2009). [CrossRef]

15. J. Walther and E. Koch, “Relation of joint spectral and time domain optical coherence tomography (jSTdOCT) and phase-resolved Doppler OCT,” Opt. Express 22(19), 23129–23146 (2014). [CrossRef]

16. M. Szkulmowski, A. Szkulmowska, T. Bajraszewski, A. Kowalczyk, and M. Wojtkowski, “Flow velocity estimation using joint Spectral and Time domain Optical Coherence Tomography,” Opt. Express 16(9), 6008–6025 (2008). [CrossRef]

17. J. Walther and E. Koch, “Flow Measurement by Lateral Resonant Doppler Optical Coherence Tomography in the Spectral Domain,” Appl. Sci. 7(4), 382 (2017). [CrossRef]

18. A. C. Chan, V. J. Srinivasan, and E. Y. Lam, “Maximum Likelihood Doppler Frequency Estimation Under Decorrelation Noise for Quantifying Flow in Optical Coherence Tomography,” IEEE Trans. Med. Imaging 33(6), 1313–1323 (2014). [CrossRef]

19. Y. K. Tao, A. M. Davis, and J. A. Izatt, “Single-pass volumetric bidirectional blood flow imaging spectral domain optical coherence tomography using a modified Hilbert transform,” Opt. Express 16(16), 12350–12361 (2008). [CrossRef]

20. A. Akif, K. Walek, C. Polucha, and J. Lee, “Doppler OCT clutter rejection using variance minimization and offset extrapolation,” Biomed. Opt. Express 9(11), 5340–5352 (2018). [CrossRef]

21. R. Haindl, W. Trasischker, A. Wartak, B. Baumann, M. Pircher, and C. K. Hitzenberger, “Total retinal blood flow measurement by three beam Doppler optical coherence tomography,” Biomed. Opt. Express 7(2), 287–301 (2016). [CrossRef]

22. A. Wartak, R. Haindl, W. Trasischker, B. Baumann, M. Pircher, and C. K. Hitzenberger, “Active-passive path-length encoded (APPLE) Doppler OCT,” Biomed. Opt. Express 7(12), 5233–5251 (2016). [CrossRef]

23. H. Spahr, C. Pfaffle, P. Koch, H. Sudkamp, G. Huttmann, and D. Hillmann, “Interferometric detection of 3D motion using computational subapertures in optical coherence tomography,” Opt. Express 26(15), 18803–18816 (2018). [CrossRef]

24. L. Ginner, A. Wartak, M. Salas, M. Augustin, M. Niederleithner, L. M. Wurster, and R. A. Leitgeb, “Synthetic subaperture-based angle-independent Doppler flow measurements using single-beam line field optical coherence tomography in vivo,” Opt. Lett. 44(4), 967–970 (2019). [CrossRef]

25. B. Karamata, M. Leutenegger, M. Laubscher, S. Bourquin, T. Lasser, and P. Lambelet, “Multiple scattering in optical coherence tomography. II. Experimental and theoretical investigation of cross talk in wide-field optical coherence tomography,” J. Opt. Soc. Am. A 22(7), 1380–1388 (2005). [CrossRef]

26. N. Uribe-Patarroyo, M. Villiger, and B. E. Bouma, “Quantitative technique for robust and noise-tolerant speed measurements based on speckle decorrelation in optical coherence tomography,” Opt. Express 22(20), 24411–24429 (2014). [CrossRef]

27. N. Uribe-Patarroyo and B. E. Bouma, “Velocity gradients in spatially resolved laser Doppler flowmetry and dynamic light scattering with confocal and coherence gating,” Phys. Rev. E 94(2), 022604 (2016). [CrossRef]

28. V. J. Srinivasan, H. Radhakrishnan, E. H. Lo, E. T. Mandeville, J. Y. Jiang, S. Barry, and A. E. Cable, “OCT methods for capillary velocimetry,” Biomed. Opt. Express 3(3), 612–629 (2012). [CrossRef]

29. I. Popov, A. Weatherbee, and I. A. Vitkin, “Statistical properties of dynamic speckles from flowing Brownian scatterers in the vicinity of the image plane in optical coherence tomography,” Biomed. Opt. Express 8(4), 2004–2017 (2017). [CrossRef]

30. N. Weiss, T. G. van Leeuwen, and J. Kalkman, “Localized measurement of longitudinal and transverse flow velocities in colloidal suspensions using optical coherence tomography,” Phys. Rev. E 88(4), 042312 (2013). [CrossRef]

31. H. W. Ren, K. M. Brecke, Z. H. Ding, Y. H. Zhao, J. S. Nelson, and Z. P. Chen, “Imaging and quantifying transverse flow velocity with the Doppler bandwidth in a phase-resolved functional optical coherence tomography,” Opt. Lett. 27(6), 409–411 (2002). [CrossRef]

32. A. Bouwens, D. Szlag, M. Szkulmowski, T. Bolmont, M. Wojtkowski, and T. Lasser, “Quantitative lateral and axial flow imaging with optical coherence microscopy and tomography,” Opt. Express 21(15), 17711–17729 (2013). [CrossRef]

33. M. G. O. Gräfe, M. Gondre, and J. F. de Boer, “Precision analysis and optimization in phase decorrelation OCT velocimetry,” Biomed. Opt. Express 10(3), 1297–1314 (2019). [CrossRef]

34. B. Braaf, K. A. Vermeer, K. V. Vienola, and J. F. de Boer, “Angiography of the retina and the choroid with phase-resolved OCT using interval-optimized backstitched B-scans,” Opt. Express 20(18), 20516–20534 (2012). [CrossRef]

35. T. Park, S. J. Jang, M. Han, S. Ryu, and W. Y. Oh, “Wide dynamic range high-speed three-dimensional quantitative OCT angiography with a hybrid-beam scan,” Opt. Lett. 43(10), 2237–2240 (2018). [CrossRef]

36. C. Akcay, P. Parrein, and J. P. Rolland, “Estimation of longitudinal resolution in optical coherence imaging,” Appl. Opt. 41(25), 5256–5262 (2002). [CrossRef]

37. J. W. Goodman, Statistical Optics, Wiley classics library ed., Wiley classics library (Wiley, 2000), pp. xvii, 550 p.

38. J. Lee, W. Wu, J. Y. Jiang, B. Zhu, and D. A. Boas, “Dynamic light scattering optical coherence tomography,” Opt. Express 20(20), 22262–22277 (2012). [CrossRef]

39. J. Tang, S. E. Erdener, B. Li, B. Fu, S. Sakadzic, S. A. Carp, J. Lee, and D. A. Boas, “Shear-induced diffusion of red blood cells measured with dynamic light scattering-optical coherence tomography,” J. Biophotonics 11(2), e201700070 (2018). [CrossRef]

40. H. L. V. Trees, Detection, Estimation, and Modulation Theory: Radar-Sonar Signal Processing and Gaussian Signals in Noise (Krieger Publishing Co., Inc., 1992), p. 646.

41. R. A. Fisher, “Theory of statistical estimation,” Math. Proc. Cambridge Philos. Soc. 22(5), 700–725 (1925). [CrossRef]

42. O. Nadiarnykh, V. Davidoiu, M. G. O. Gräfe, M. Bosscha, A. C. Moll, and J. F. de Boer, “Phase-based OCT angiography in diagnostic imaging of pediatric retinoblastoma patients: abnormal blood vessels in post-treatment regression patterns,” Biomed. Opt. Express 10(5), 2213–2226 (2019). [CrossRef]

43. B. Braaf, K. A. Vermeer, M. de Groot, K. V. Vienola, and J. F. de Boer, “Fiber-based polarization-sensitive OCT of the human retina with correction of system polarization distortions,” Biomed. Opt. Express 5(8), 2736–2758 (2014). [CrossRef]

44. M. Salas, M. Augustin, F. Felberer, A. Wartak, M. Laslandes, L. Ginner, M. Niederleithner, J. Ensher, M. P. Minneman, R. A. Leitgeb, W. Drexler, X. Levecq, U. Schmidt-Erfurth, and M. Pircher, “Compact akinetic swept source optical coherence tomography angiography at 1060 nm supporting a wide field of view and adaptive optics imaging modes of the posterior eye,” Biomed. Opt. Express 9(4), 1871–1892 (2018). [CrossRef]

45. S. Martinez-Conde, S. L. Macknik, and D. H. Hubel, “The role of fixational eye movements in visual perception,” Nat. Rev. Neurosci. 5(3), 229–240 (2004). [CrossRef]

46. M. Rucci, P. V. McGraw, and R. J. Krauzlis, “Fixational eye movements and perception,” Vision Res. 118, 1–4 (2016). [CrossRef]

47. M. Sugita, R. A. Brown, I. Popov, and A. Vitkin, “K-distribution three-dimensional mapping of biological tissues in optical coherence tomography,” J. Biophotonics 11(3), e201700055 (2018). [CrossRef]

48. P. Cimalla, J. Walther, M. Mittasch, and E. Koch, “Shear flow-induced optical inhomogeneity of blood assessed in vivo and in vitro by spectral domain optical coherence tomography in the 1.3 mu m wavelength range,” J. Biomed. Opt. 16(11), 116020 (2011). [CrossRef]

49. M. T. Bernucci, C. W. Merkle, and V. J. Srinivasan, “Investigation of artifacts in retinal and choroidal OCT angiography with a contrast agent,” Biomed. Opt. Express 9(3), 1020–1040 (2018). [CrossRef]

50. R. F. Spaide, J. G. Fujimoto, and N. K. Waheed, “Image Artifacts in Optical Coherence Tomography Angiography,” Retina 35(11), 2163–2180 (2015). [CrossRef]

51. A. Bouwens, T. Bolmont, D. Szlag, C. Berclaz, and T. Lasser, “Quantitative cerebral blood flow imaging with extended-focus optical coherence microscopy,” Opt. Lett. 39(1), 37–40 (2014). [CrossRef]

$τ$	$v_{min} (β = 0.1)$	$v_{max} (β = 1.0)$
1ms	1mm/s	10mm/s
100µs	10mm/s	100mm/s
10µs	100mm/s	1000mm/s

Optical coherence tomography velocimetry based on decorrelation estimation of phasor pair ratios (DEPPAIR)

Abstract

1. Introduction

2. Theory

3. Experiments

3.1 Setup

3.2 Processing

4. Results

5. Discussion

6. Conclusion

Funding

Disclosures

References

Supplementary Material (1)

Cited By

Figures (8)

Tables (1)

Equations (20)

Biomedical Optics Express