Precision analysis and optimization in phase decorrelation OCT velocimetry

Maximilian G. O. Gräfe; Maude Gondre; Johannes F. de Boer

doi:10.1364/BOE.10.001297

1. Introduction

Common retinal diseases such as diabetic retinopathy, glaucoma and age-related macular degeneration have shown changes to the vasculature or blood flow during their progression [1–6]. Optical coherence tomography (OCT) has become a widely used technique for the studies of disease developments and clinical diagnosis while OCT angiography (OCTA) has been proven to show a particular strong contrast between static tissue and moving particles e.g. blood cells [7]. This makes OCTA already a powerful tool to detect morphological changes in the retinal vasculature and formation of abnormal vessels e.g. neovascularization [8–11]. Because early symptoms are often connected to reduction of healthy vessels [10] or reduced functionality of existing vessels [5] a quantification of flow is becoming increasingly important. OCT velocimetry is still in its infancy in clinical practice, but it is a prime candidate for further research towards the understanding and early diagnosis of ocular diseases.

The first established technique to determine the flow velocity of particles in OCT was Doppler OCT [12–14]. Doppler OCT has the disadvantage that only the velocity component parallel to the propagation direction of the light can be measured. Several attempts have been made in the past years to overcome this limitation, among these: Doppler OCT with multiple directions of illumination/detection [15,16], (complex) speckle decorrelation [17–19] and Doppler mean and bandwidth analysis [20,21]. Doppler OCT with multiple directions of illumination requires a significant increase in system complexity [15,16], which can be a limiting factor especially if other functionalities should be incorporated. (Complex) speckle decorrelation and Doppler mean and bandwidth analysis requires in general acquisition of large data ensembles for each voxel and velocity estimation [17–19]. Therefore, this approach was demonstrated in situations when samples can be fixated and bulk motion can be suppressed. The long acquisition times required to acquire the large data sets for these methods limits the applicability for ophthalmology.

An important constraint for the applicability of OCT velocimetry in vivo is the limited acquisition time and thus the number of samples at a tissue voxel that are available to achieve a certain flow estimation precision. An analysis investigating the minimally required number of samples and the optimal method or combination of methods to estimate the flow velocity precision has to our knowledge not been presented. Therefore we set out to provide a theoretical derivation of the probability density function (PDF) of changes in the OCT signal between repeated measurements at (nearly) the same location as a result of motion in individual voxels, and derive a prediction of the precision of flow velocity estimations based on the PDF.

Since scan patterns optimized for flow analysis may require different time delays between repeated measurements, we restricted the analysis to PDFs that only describe the signal change between two individual measurements. The motivation for that is the following: some methods [19–22] require all samples at (nearly) the same location belonging to one ensemble to be a continuous trace of equidistantly sampled data in time and are optimized for M-scan-like acquisition modes. In many methods for in vivo applications various scan patterns are used in which time intervals between scan repetitions are optimized for sensitivity in a particular velocity range or multiple time intervals are used within the same estimation [8,23]. Our analysis focuses on the latter data acquisition method using only pairs of data for a time interval. As a first demonstration of this analysis we choose the use of Doppler OCT in combination with phase decorrelation as previously introduced [20,21,24].

The manuscript is organized as follows: In section two the theoretical model is described, followed by a numerical simulation based on the theoretical model. An established method by Bouwens et al. [21] is compared to additionally proposed methods. The limit of the Cramér-Rao lower bound is predicted and compared to the achieved precision. In section three the experimental setup and the data processing is described. In section four the experimental measurements and theoretically predicted results are compared. The article will be concluded by a short summary and a discussion about potential future extensions of the theory.

2. Theory & simulations

2.1 Theory

In order to describe light incident on the sample and the collection of backscattered light in the sample arm of an OCT system, an experimental situation similar to the description by Bouwens et al. [21] is considered. An electrical field propagates from the principal plane of a lens or objective to the scattering volume of a sample and back. The incident field on the volume of scattering particles can be described as superposition of (monochromatic) plane waves with weight $m_{i n} (p_{i n})$ :

E_{i n} (r) = \int m_{i n} (p_{i n}) e^{i r p_{in}} d p_{⊥,}_{i n}

With

p

the wave vector of the electromagnetic field,

k = 2 π / λ = \sqrt{p_{x}^{2} + p_{y}^{2} + p_{z}^{2}}

,

λ

the wavelength,

r

defined with respect to the focal point of the lens as origin and

p_{⊥}

the wave vector components

p_{x}

and

p_{y}

perpendicular to the optic axis z.

E_{i n} (r)

is represented as a field incident from the positive

p_{z}

half space. This principle is also illustrated in Fig. 1. In our setup a Gaussian incident beam was used, in which case the plane wave weights

m_{i n} (p_{i n})

are given by

m_{i n} (p) = \exp [- w_{0}^{2} (p_{x}^{2} + p_{y}^{2}) / 4]

with

w_{0}

the radius of the beam waist. The backscattered field can be formulated as [25]:

E_{S} (r_{2}; p_{o u t}) = \sum_{n} \int G (r_{2}, r_{1}; p_{o u t}) F_{n} (r_{1}) E_{i n} (r_{1}) d r_{1}

with

G (r_{2}, r_{1}; p_{o u t})

the Green’s function for propagating from

r_{1}

to

r_{2}

with wave vector

p_{o u t}

and

F_{n} (r_{1}; k) = k^{2} (n^{2} (r_{1}) - 1) / 4 π

the scattering strength of the n-th particle in the first order Born approximation, similar to ref [21]. with

n

the refractive index. When the light is collected by a lens with focal length

f

, and the focal point of the lens is located at the origin, the Green’s function

G (r_{2}, r_{1}; p_{o u t})

for propagation to the principal plane of the lens can be approximated assuming that the focal length

f

is much larger than the focal volume from which the backscattered light is collected, i.e.

f ≫ | r_{1} |

,

G (r_{2}, r_{1}; p_{o u t}) = \frac{e^{i p_{o u t} (r_{1} - r)}}{| r_{1} - r_{2} |} \approx \frac{e^{i k f}}{f} e^{- i p_{o u t} r_{1}}

The detected field has to be weighted by the acceptance mode

m_{o u t} (p_{o u t})

with which light is collected, and

r_{2}

is replaced by the focal length

f

:

E (f, k) = \int m_{o u t} (p_{o u t}) E_{s} (f, p_{o u t}) d p_{⊥, o u t}

In case of a fiber based system the numerical aperture and possible propagation modes determine the accepted modes. For simplicity only point scatterers are considered,

F_{n} (r_{1}, k) \propto k^{2} δ (r_{n} - r_{1})

, where the scattering strength based on the refractive index contrast was omitted. When applying the above mentioned approximation, and combining Eqs. (1)-(4), the result for the detected field is given by:

\begin{matrix} E (r_{n}, k) = \frac{k^{2} e^{i k f}}{f} \sum_{n} \underset{E_{i n}}{\underset{︸}{\int m_{i n} (p_{i n}) \exp [- i r_{n} p_{i n}] d p_{⊥, i n}}} \\ \times \underset{E_{o u t}}{\underset{︸}{\int m_{o u t} (p_{o u t}) \exp [i r_{n} p_{o u t}] d p_{⊥, o u t}}} . \end{matrix}

Equation (5) expresses the field that is detected given a collection of point scatterers at locations

r_{n} .

In an OCT configuration, not only the lateral dimension of the beam plays a role in the intensity of the backscattered light, but also the coherence gate. Lateral dimension and coherence gate define a volume from which the OCT sample arm field is detected. The scatterers at locations

r_{n}

need to fall within this volume. To add the axial resolution a factor

\exp (- z^{2} / l_{c}^{2})

needs to be added to Eq. (5), with

l_{c} = \sqrt{2 \ln 2 / π} \cdot λ_{0}^{2} / Δ λ

the coherence length [26] with

λ_{0}

the central wavelength and

Δ λ

the full width half maximum bandwidth. Omitting the prefactor depending on the focal length, the ensemble averaged intensity

E (r, k)

that is backscattered from location r is now proportional to:

\begin{matrix} E (r, k) \propto e^{- z^{2} / l_{c}^{2}} \underset{E_{i n} (r)}{\underset{︸}{\int m_{i n} (p_{i n}) \exp [- i r p_{i n}] d p_{⊥, i n}}} \\ \times \underset{E_{o u t} (r)}{\underset{︸}{\int m_{o u t} (p_{o u t}) \exp [i r p_{o u t}] d p_{⊥, o u t}}} \end{matrix}

A very common realization of OCT which applies also in case of our setup (see 3. Experimental setup) consists of the output of a single mode fiber that is collimated and imaged on a sample (e.g. the retina) by a lens. The output mode of a single mode fiber is Gaussian. In this configuration

m_{i n} (p_{i n})

and

m_{o u t} (p_{o u t})

are determined by the optical system, and since the detected light is reflected back to the single mode fiber,

m_{i n} (p_{i n}) = m_{o u t} (p_{o u t})

. That makes Eq. (6) symmetric for the incident and detection field when we consider additionally that the integration over half spheres are opposite (incident and backscattered directions are opposite). A configuration in which the illumination and collection geometry is different can be accounted for by a different

m_{i n} (p_{i n})

and

m_{o u t} (p_{o u t})

, where

m_{i n} (p_{i n}) \neq m_{o u t} (p_{o u t})

.

Fig. 1 Schematic illustration of the incident and backpropagating field focused by a lens or objective onto a sample. The fields are considered as superpositions of plane waves with the wave vectors $p_{i n}$ and $p_{o u t}$ which are weighed by $m_{i n} (p_{i n})$ and $m_{o u t} (p_{o u t})$ , respectively.

Download Full Size | PDF

The integration in Eq. (6) has to be performed in the $(p_{x}, p_{y})$ -plane, where $p_{z}$ is determined by the relation $| p | = k$ . Under the approximation of small numerical apertures $(p_{z}^{2} ≫ p_{x}^{2}, p_{y}^{2})$ we find:

p_{z} = k \sqrt{1 - \frac{p_{x}^{2}}{k^{2}} - \frac{p_{y}^{2}}{k^{2}}} \approx k (1 - \frac{p_{x}^{2}}{2 k^{2}} - \frac{p_{y}^{2}}{2 k^{2}})

Now Eq. (7) can be used to solve Eq. (6). Due to the symmetry of Eq. (6) for

m_{i n}

and

m_{o u t}

(p_{i n, z} = - p_{o u t, z})

only the integral for the incident field is carried out explicitly:

E_{i n} (r) = e^{i z k} \underset{e^{^{- \frac{(x^{2} + y^{2})}{w {(z)}^{2}}}}}{\underset{︸}{e^{- \frac{k^{2} w_{0}^{2} (x^{2} + y^{2})}{k^{2} w_{0}^{4} + 2 z^{2}}}}} \underset{e^{^{- i \frac{k (x^{2} + y^{2})}{2 R (z)}}}}{\underset{︸}{e^{- i \frac{2 z k (x^{2} + y^{2})}{k^{2} w_{0}^{4} + 4 z^{2}}}}} \cdot \underset{\frac{w_{0}}{w (z)} e^{i ζ (z)}}{\underset{︸}{\frac{w_{0}^{2}}{w_{0}^{2} - i \frac{2 z}{k}}}}

with

w (z) = w_{0} \sqrt{1 + {(z / z_{R})}^{2}}

,

z_{R}

the Rayleigh length,

R (z)

the radius of curvature and

ζ (z)

the Gouy-phase. The field

E_{i n} (r)

in Eq. (8) contains all terms of a Gaussian beam [27].

With Eq. (5) and Eq. (6), the explicit solution for $E_{i n} (r)$ in Eq. (8) and the equivalence of $E_{i n} (r)$ and $E_{o u t} (r)$ , we have derived expressions for the detected complex field backscattered from a detection volume confined laterally by the beam diameter and in depth by the coherence length of the source containing scatterers (Eq. (5) with inclusion of the coherence gate) and the ensemble averaged field $E (r)$ (Eq. (6)) as a function of $r$ , respectively.

Next we will determine the change in the phase due to flow. In our analysis we ignore the effect of Brownian motion. Velocities and time differences between repeated scans are considered for which we assume Brownian motion has no significant influence on the signal. Furthermore, we assume bulk motion, a collective motion of all particles in the focal volume. Therefore, flow can be modelled as a translation of all scatterers in the focal volume by a vector $δ r$ . Following the approach of Vakoc et al. [24], the effect of flow is described by a displacement in a second measurement with a bulk transversal and/or axial translation of the detection volume. Two consecutive measurements of the sample result in phasor $Γ_{A}$ and phasor $Γ_{B}$ . The second phasor $Γ_{B}$ can be decomposed into a contribution of the first phasor $Γ_{A}$ and a random phasor $Γ_{C}$ representing the decorrelation. Phasor $Γ_{B}$ can then be expressed as:

Γ_{B} = μ_{A} Γ_{A} + μ_{C} Γ_{C},

where

μ_{A}

and

μ_{C}

are parameters reflecting the weight of

Γ_{A}

and

Γ_{C}

in the composition of

Γ_{B}

and should adhere to the normalization

μ_{A}^{2} + μ_{C}^{2} = 1

. The weight

μ_{A}

can be set equal to a parameter

α

(

μ_{A} = α

), where Vakoc demonstrates this

α

to be equal to the lateral overlap integral of the detected fields of consecutive measurements. Here we evaluate the volume overlap integral including the axial coherence gate. The field reflected from detection volume A is evaluated as,

E_{A} (r) = E_{i n} (r) E_{o u t} (r) \exp (- z^{2} / l_{c}^{2})

The field phasor

Γ_{B}

reflected from detection volume B after motion is defined by the following substitutions:

z \to z + δ z

for axial motion and

x \to x + δ x

for lateral motion. The overlap integral is defined as:

α = \frac{∭ E_{B} (x + δ x, y, z + δ z) E_{A}^{*} (x, y, z) d^{3} r}{\sqrt{∭ E_{B} (x + δ x, y, z + δ z) E_{B}^{*} (x + δ x, y, z + δ z) d^{3} r} \sqrt{∭ E_{A} (x, y, z) E_{A}^{*} (x, y, z) d^{3} r}}

with

d^{3} r = d x d y d z

. Equation (11) cannot be solved analytically for the most general case but we will demonstrate that a solution can be found under the condition that the Rayleigh length

z_{R}

is much larger than the coherence length

l_{c}

. In addition to the displacement by flow we also evaluate the effect when the sample is significantly out of focus. For this purpose an additional displacement in z is introduced: the substitution

z \to z - Δ z

is performed in

E_{A}

. Since the same location but with motion of the scattering particles is described the substitution

z \to z - Δ z + δ z

is performed in

E_{B}

. The numerator of Eq. (11) can then be written as:

\begin{matrix} ∭ E_{B} (r + δ r) E_{A}^{*} (r) d^{3} r = e^{i 2 k_{0} δ z} e^{- \frac{δ z^{2}}{2 l_{c}^{2}}} \\ \times \int \frac{w_{0}^{2}}{w {(z - Δ z + δ z)}^{2}} \cdot \frac{w_{0}^{2}}{w {(z - Δ z)}^{2}} e^{i (ζ (z - Δ z) - ζ (z - Δ z + δ z))} \\ \times \int \int e^{^{- [\frac{2 [{(x + δ x)}^{2} + y^{2}]}{w {(z - Δ z + δ z)}^{2}} + \frac{2 (x^{2} + y^{2})}{w {(z - Δ z)}^{2}}]}} e^{^{- i \frac{k [{(x + δ x)}^{2} + y^{2}]}{R (z - Δ z + δ z)} + i \frac{k (x^{2} + y^{2})}{R (z - Δ z)}}} d x d y \cdot e^{- \frac{2 {(z - Δ z + δ z / 2)}^{2}}{l_{c}^{2}}} d z \end{matrix}

Now the independence of the overlap integral on the defocus will be demonstrated, i.e., the independence on

Δ z

. First the effect of

δ z

has to be evaluated. Due to a rather rapid decorrelation caused by the axial resolution (first line of Eq. (12)) the effect of

δ z

on

w (z)

and

R (z)

for low numerical apertures will be negligible. The term of the Gouy-phase in the second line and the denominators in the exponents of the third line are then considered constant and independent of

δ z

. Next the substitution

z^{'} = z - Δ z

can be made in Eq. (12), resulting in:

\begin{matrix} ∭ E_{B} E_{A}^{*} d^{3} r = e^{i 2 k_{0} δ z} e^{- \frac{δ z^{2}}{2 l_{c}^{2}}} \int {(\frac{w_{0}^{2}}{w {(z^{'})}^{2}})}^{2} \cdot e^{- \frac{2 {(z^{'} + δ z / 2)}^{2}}{l_{c}^{2}}} \\ \times \int \int e^{^{- [\frac{2 [{(x + δ x)}^{2} + x^{2} + 2 y^{2}]}{w {(z^{'})}^{2}}]}} e^{^{- i \frac{k [{(x + δ x)}^{2} - x^{2}]}{R (z^{'})}}} d x d y d z^{'} . \end{matrix}

Now we will turn to the integration over

z^{'}

. In most common OCT applications the coherence length

l_{c}

is much smaller than the Rayleigh-length to be able to resolve many pixels along z with more or less uniform lateral resolution. The exponent

2 {(z^{'} + δ z / 2)}^{2} / l_{c}^{2}

determines the range over which the integral along

z^{'}

receives significant contributions. Within this range of order

l_{c}

,

w (z^{'})

and

R (z^{'})

do not change significantly and are therefore approximated by constants:

w (ξ)

and

R (ξ)

with

ξ

being the average value of

z^{'}

. This transforms Eq. (13) into:

\begin{matrix} ∭ E_{B} E_{A}^{*} d^{3} r = e^{i 2 k_{0} δ z} e^{- \frac{δ z^{2}}{2 l_{c}^{2}}} {(\frac{w_{0}^{2}}{w {(ξ)}^{2}})}^{2} \\ \times \int \int e^{^{- [\frac{2 [{(x + δ x)}^{2} + x^{2} + 2 y^{2}]}{w {(ξ)}^{2}}]}} e^{^{- i \frac{k [{(x + δ x)}^{2} - x^{2}]}{R (ξ)}}} d x d y \\ \times \int e^{- \frac{2 {(z^{'} + δ z)}^{2}}{l_{c}^{2}}} d z^{'} . \end{matrix}

Now only the integrant in the third line of Eq. (14) depends on

z^{'}

. The integrals are evaluated at the center of the coherence gate

(ξ = 0)

in analogy to the derivation previously described by Popov et al. [22]. In conclusion, the explicit solution of Eq. (11) is independent of the location of the coherence gate with respect to the beam focus and is given by

α = e^{i 2 k_{0} δ z} e^{- \frac{δ x^{2}}{w_{0}^{2}}} e^{- \frac{δ z^{2}}{2 l_{c}^{2}}}

Equation (15) describes the effect of flow velocity on the overlap integral if we replace the displacement

δ r

by

v \cdot τ

, with

v

the velocity and

τ

the time delay between consecutive measurements. Equation (15) is also used in methods for speckle decorrelation for which this factor describes the decay of the autocorrelation function [18,19,22]. It should be emphasized here that Eq. (15) shows that the overlap integral is independent of the distance of the scattering volume to the focal plane, as derived by Popov et al. [22] and previously demonstrated experimentally by Weiss et al. [28]. The remaining calculation of the expected distribution function of phase differences is analogous to Vakoc et al. [24] with one difference: In [24] only lateral motion was considered and investigated. To achieve the PDF of phasor

Γ_{B}

, the PDF of the sum of the phasors

Γ_{A}

and

Γ_{C}

has to be calculated (Eq. (9)), with the prefactor

α

determined by the overlap integral. The PDFs of the individual phasors are known to follow a Rayleigh and circular Gaussian distribution, respectively [24,29]:

P_{A} (x_{Γ}, y_{Γ} | μ_{A}) = \frac{2 x_{Γ}}{μ_{A}^{2}} \exp (- \frac{x_{Γ}^{2}}{μ_{A}^{2}})

P_{C} (x_{Γ}, y_{Γ} | μ_{C}) = \frac{1}{π μ_{C}^{2}} \exp (- \frac{x_{Γ}^{2} + y_{Γ}^{2}}{μ_{A}^{2}})

The PDF

P_{B}

for phasor B results from the convolution of

P_{A}

and

P_{C}

. It should be pointed out that due to a difference in definition of the fields in the overlap integral (Vakoc et al. defined only incident fields, we define incident fields and subsequent backpropagation) we have

μ_{A} = α

(Eq. (11)) and

μ_{C} = \sqrt{1 - α^{2}}

where Vakoc used

μ_{A} = α^{2}

. The complex part of Eq. (15) can be identified with the mean Doppler shift caused by axial motion. It represents a circular shift of the final PDF of phase differences by a phase

2 k_{0} δ z

. Since in this analysis we are interested in the shape of the PDF and not so much its center position, we neglect the phase

2 k_{0} δ z

in the further calculation. The convolution of

P_{A}

and

P_{C}

leads to [24]:

P_{B} = \frac{1}{π} e^{- ((x_{Γ}^{2} + y_{Γ}^{2}) / (1 - α^{2}))} + \frac{α^{2}}{\sqrt{π (1 - α^{2})}} x_{Γ} \cdot e^{- ((x_{Γ}^{2} + y_{Γ}^{2} - α^{2} x_{Γ}^{2}) / (1 - α^{2}))} e r f c (- x_{Γ} \frac{α^{2}}{\sqrt{1 - α^{2}}})

with

e r f c (...)

the complementary error function. To obtain the PDF to measure the phase difference

φ

given a value for the overlap integral of

α

, Eq. (18) is transformed to polar coordinates and integrated over the radial component (amplitude) [24]:

P (φ | α) = \frac{1 - α^{2}}{2 π (1 - α^{2} \cos^{2} φ)} [1 + \frac{α \cos φ}{\sqrt{1 - α^{2} \cos^{2} φ}} (\tan^{- 1} [\frac{α \cos φ}{\sqrt{1 - α^{2} \cos^{2} φ}}] + \frac{π}{2})]

The goal of any estimator is now to estimate the parameter $α$ given a set of measurements from which flow velocities can be calculated using the ratio of displacement and axial/ lateral resolution as given in Eq. (15).

2.2 Simulations

In this subsection we simulate multiple repeated phasor measurements for a given overlap integral value $α$ (and therefore a given flow velocity) to create an ensemble of measurements. Randomized ensembles for the same $α$ are created repeatedly to determine the theoretical estimation precision. In all following graphs the uncertainty plotted as standard deviation of the estimated velocities is used as an inverse measure for precision. Thus, high uncertainty indicates poor precision. Equation (19) is used as a starting point for any simulations within this subsection. Further assumptions are: 1. that similar to Srinivasan et al. [18] an isotropic voxel size is used making the broadening of the phase distribution equally sensitive to axial and lateral motion; 2. during prior processing steps the distribution has been shifted circularly to center the mean phase difference at 0 rad. This eliminates the axial Doppler shift $2 k_{0} δ z$ . Those two assumptions have the same effect as if only lateral motion is simulated.

One method to extract $α$ from phase sensitive measurements was adapted from Bouwens et al. [21]. The distribution of phase differences was fitted with a Gaussian function:

C = a \cdot \exp (- \frac{φ^{2}}{2 σ^{2}}) + h

where

σ

is the standard deviation of the Gaussian function,

a

is a factor to account for non-normalized data and

h

is an offset. The standard deviation is then interpreted as ratios

(σ = δ r / w_{0})

. If data sets are considered for which the axial Doppler shift

2 k_{0} δ z

has not been eliminated, it can be accounted for by replacing

φ

with

φ - 2 k δ z

. Bouwens et al. developed his method originally for Doppler frequencies. However, the goal of this paper is to apply his method to pairs of repeated scans without implying the condition that all samples of one ensemble belong to an equidistantly sampled time trace. Thus, the Doppler frequencies are not accessible directly. Figure 2 shows in blue graphs for the PDF of Eq. (19) for a) a small and b) a large ratio

δ r / w_{0}

. In red the fit of Eq. (20) is shown. As can be seen in Fig. 2(a) and 2(b) not just the Gaussian function broadens for larger

δ r / w_{0}

, but also the offset

h

increases. The ratio of the amplitude of the Gaussian

a

and the offset

h

increases with increasing velocity, as shown in Fig. 2(c). This information is not used in Bouwens method but can be used to gain additional information content about the velocity.

Fig. 2 In blue graphs of the PDF $P (φ | α)$ given in Eq. (19) for different ratios of displacement. a) $δ r / w_{0} = 0.05$ . b) $δ r / w_{0} = 1$ . In red the fits with a Gaussian (Eq. (20)) are shown. In b) $a$ and $h$ from Eq. (20) are indicated c) The ratio $h / a$ as a function of $δ r / w_{0}$ . The ratio increases with increasing displacement.

Download Full Size | PDF

In the following a simulation of measurements based on a finite number of phase difference samples and the effect on the flow precession for different estimation methods is described. For every velocity and therefore for every ratio $δ r / w_{0}$ which is subsequently referred to as normalized velocity a number of samples $N$ is drawn from the distribution. All computations are performed in Matlab 2014a (The Mathworks, Inc.). A set $A$ of $N$ randomized uniformly distributed numbers between 0 and 1 are generated (rand() function). The PDF $P (φ, α)$ is computed and the cumulative probability function $C (φ, α)$ is derived. The numbers of set $A$ are projected to a set $D (φ)$ of phase differences according to $C (φ, α)$ . Three estimation methods are used to estimate the velocity from the simulated phase difference measurements. 1. Gauss. fit: a histogram with 100 bins of phase differences is created from $D (φ)$ spanning from the minimum to the maximum phase difference, Eq. (20) is fit to it and $σ$ is related to the normalized velocity as discussed above. Of three determined parameters this method uses only one and therefore excludes information which could be used for a reduction of the estimation uncertainty. In order to use the additional information a metric must be found to combine data from the offset together with the data from the Gaussian curve. This requirement and the similarity to phase variance algorithms motivated the next estimation method. 2. second moment: the standard deviation as statistical second moment is computed. As the second moment also takes into account measurements which in the Gauss. fit contribute only to the offset, it is not linearly related to the velocity anymore. This non-linear relation was accounted for by a look-up-table (LUT) which was used to assign the second moment to a normalized velocity. The LUT was generated in advance by the evaluation of the second moment of $P (φ, α)$ for a defined range of velocities. Using the second moment is a method to characterize the PDF. To extract all remaining information, all properties of the PDF must be used which motivates the last estimation method. 3. MLE: the normalized velocity was determined by a maximum likelihood estimator (MLE) defined as:

{(\frac{δ r}{w_{0}})}_{M L E} = \underset{δ r / w_{0}}{\arg \max} \sum_{j} \log P (φ_{j} | δ r / w_{0})

where

j

is the index of the

j

-th phase difference. For the implementation used here, the likelihood for each phase difference was determined the following way: prior to the estimations, PDFs were created for a range of normalized velocities between 0.01 and 1.8 in 200 steps where each PDF had a resolution with 1024 datapoints in phase space from -π to π. For a phase difference

φ_{j}

the two nearest phase difference datapoints and the respective values of the PDF were found and the likelihood value was determined by linear interpolation between the values of the PDF. This was performed for every

φ_{j}

in an ensemble for all velocities and the one velocity was chosen as estimated velocity which maximized the joint likelihood function. The estimators were determined over many iterations where a new set

D (φ)

was created for each iteration. The accuracy and precision (as statistical fluctuation/ uncertainty) of each estimator was quantified by the mean and standard deviation of their results, respectively. Figure 3 shows the results of this simulation for N = 290, 20 different normalized velocities and 5000 iterations.

Fig. 3 Simulation of the estimated velocity and the velocity precision as a function of actual velocity determined by three estimation methods: Gauss. fit (blue), second moment (red) and MLE (black). Velocity is expressed as a ratio of $δ r / w_{0}$ . a) estimated velocity as a function of input velocity. Crosses show the average estimated velocity and the full error bar length indicates two times the standard deviation. b) zoom in from a) while omitting the error bars. The solid line represents a slope of 1 (indication when estimated is the same as set ratio). c) standard deviation of the fluctuation as a function of average velocity, showing the estimation precision. d) zoom in from c) for a better visibility of the differences between second moment and MLE. Additionally, the Cramér-Rao lower bound (CRLB, solid black line) was included. The MLE comes closest to the Cramér-Rao lower bound indicating a low estimation uncertainty. Each simulation was based on N = 290 samples, each simulation was repeated 5000 times for a given actual velocity.

Download Full Size | PDF

Figure 3(a) shows the mean of all estimated flow velocities as crosses and statistical fluctuation as error bar. All methods provide an accurate mean flow velocity over a large range of input velocities. Figure 3(b) shows the results from (a) enlarged on the upper range of set velocities. For accurate results the estimations should on average fall on a slope of 1 (solid line). More deviations are observed for the Gauss. fit, where the actual velocity is first underestimated, and for velocities approaching 1.5 the actual velocity is overestimated. The MLE provides an accurate estimate over the full velocity range. The results for the second moment method are not shown in Fig. 3(b) due to high overlap with the MLE. The statistical fluctuations are increasing the strongest for the Gaussian fit method. This is shown more quantitatively in Fig. 3(c) where the standard deviation of the velocity estimation is plotted as a function of the velocities. Here the Gaussian fit method shows a strong nonlinear increase with velocities. Especially for high velocities it is increasing much faster than for the other two estimation methods. To show the differences between the second moment and the MLE, the result of Fig. 3(b) is plotted on a smaller vertical axis in Fig. 3(c). In Fig. 3(c) the Cramér-Rao lower bound was included and is defined as [30]:

σ_{C R}^{2} = {- E [\frac{\partial^{2} \ln p (x | θ)}{\partial^{2} θ}]}^{- 1}

where

p (x | θ)

is the PDF for a one-dimensional variable and

E [...]

the expectation value of any variable inside the brackets. In Eq. (22) the PDF is

P (φ | α)

. Equation (22) describes the lowest standard deviation possible due to the sensitivity of

p (x | θ)

to a given parameter (here the normalized velocity). The simulated precision of the flow estimate of the MLE follows exactly the Cramér-Rao lower bound and deviates slightly from the Cramér-Rao lower bound only for the highest velocities plotted. The second moment method has a small offset in comparison to the MLE. The precision of the second moment method is worse than the Gauss. fit method for relative velocities up to 0.25. Hence, it does not perform better than the Gauss. fit method for small velocities (smaller than 0.25) but represents still a significant improvement for the largest part of the velocity range. The MLE clearly outperforms all other methods. Ideally, the statistical error in the determination of flow velocities should increase only linearly. This is the case for normalized velocities between 0 and 1, indicating the velocity range providing the best velocity estimate.

3. Experimental setup

3.1 OCT-system

The experimental system shown in Fig. 4 was described in [31] and was modified for this study to be used not only for polarization sensitive data acquisition but also for phase-sensitive and therefore Doppler measurements. The light from a 1-μm wavelength swept source (Axsun Technologies, Inc., MA, USA) was split into a lead for the sample arm and reference arm. After being sent through a polarization delay unit (PDU) it was split in the sample arm again by a 20/80 coupler. With respect to the earlier presented system [31] the calibration mirror in one lead of the 20/80 coupler was replaced by a combination of circulator, fiber-bragg-grating (FBG) and photodiode delivering a stable A-line trigger to the acquisition board (Alazar ATS 9360). 10% of the light in the reference arm was split off to a Mach-Zehnder interferometer for k-clocking. After recombination of sample and reference arm the signal was again split into two polarization channels and recorded via balanced detection.

Fig. 4 Interferometer for phase sensitive measurements as well as polarization sensitive acquisition. PDU: polarization delay unit, P: polarizer, PC: polarization controller, PBS: polarization beam splitter, FBG: fiber-bragg-grating, PD: photo diode

Download Full Size | PDF

3.2 Phantoms

After the signal was guided through the ophthalmic interface it was sent to a model eye shown in Fig. 5(a). The model eye included a plano-convex lens (15mm focal length, LA1540-B, Thorlabs Inc.) and a scattering medium (0.2% Titanium dioxide in cured silicone) with an embedded glass capillary (inner diameter 150 µm, Molex Inc.). The capillary was connected to a syringe pump (NE-4000, New Era Pump Systems Inc.) to create a constant flow rate of a 0.5% aqueous intralipid solution. The fluid was pumped with a flow rate of 0.371 ml/min. Figure 5(b) shows the structural image of a single B-scan of the phantom. The glass capillary is embedded at a depth of approximately 300 μm.

Fig. 5 a) Schematic drawing of the flow phantom setup in front of the ophthalmic interface (same as in [31]). f: lens of 15 mm focal length, Ph: phantom consisting of a scattering medium with embedded capillary with 150 µm inner diameter, R: reservoir for intralipid. b) Structural image of a single B-scan. The part of the static medium used for the analysis is delineated in yellow and the inner boundary circle of the capillary delineated in red.

Download Full Size | PDF

3.3 Data processing

The configuration of the system for polarization sensitivity creates four images in total by depth-multiplexing in the PDU and detection of two orthogonal polarization channels. A more detailed description can be found in [31]. First, a single phase image must be computed out of the individual images. After initial processing steps of fixed pattern removal and dispersion compensation the displacement between the depth-multiplexed images is identified and the images are overlapped. The four fields are averaged for each A-line according to [32]:

\bar{E_{o u t}} (z) = \frac{1}{4} [E_{11} (z) + e^{- i θ_{12}} E_{12} (z) + e^{- i θ_{21}} E_{21} (z) + e^{- i θ_{22}} E_{22} (z)]

In Eq. (23) a correction phase

θ_{m n} = \arg [\sum_{z} E_{11} (z) E_{m n}^{*} (z)]

is used to account for average phase differences between a reference image (here the indexes 1,1 were chosen) and the other three images. Bulk motion compensation has been applied to

\bar{E_{o u t}} (z)

similar to previously described in [33].

The calculation of normalized velocities and the standard deviation of their fluctuation has been processed in two ways. Intralipid: adjacent A-lines (0.77 μm apart) were used to calculate phase differences from which velocities are estimated. Each B-scan has been scanned repeatedly to create ensembles for each pixel. The mean Doppler shift was determined by circular averaging and afterwards the phase distribution was shifted circularly to bring the mean phase difference to zero. Under the condition of laminar flow a centrosymmetric flow profile is expected inside the capillary. Therefore pixels belonging to the same velocities are located on circles around the center or the capillary. The inner wall of the capillary was segmented manually and the isotachs are defined by a reduction of the diameter of the wall segmentation centered at the capillary. To avoid effects caused by multiple scattering such as those known as OCT angiography projection artifacts [34] only velocities on the upper half of each circle were used for the analysis of the precision. Static medium: adjacent A-lines (0.77 μm apart) were recorded with a high degree of beam overlap. By comparing A-lines with different displacements, flow can be mimicked because the static medium will look like a homogeneously moving object. For each displacement multiple measurements were collected in ensembles and the velocity was estimated for each ensemble. Afterwards the mean and standard deviation of the velocities were calculated for each displacement.

4. Results

Figure 6 shows the results for measurements of the static medium (a and c) and intralipid (b and d). According to theory the level of the precision of the velocity estimation depends on the number of samples per estimation. Thus, the measurements have been evaluated for a high number (N = 290, a and b) and a low number (N = 33, c and d) of phase differences. The uncertainty as standard deviation was plotted as a function of the experimentally determined normalized velocities. The results from the second moment are not plotted for the measurements of intralipid for better visibility of the other estimators. Overall, a good to excellent agreement was found between theory and experiment. In all cases, the MLE performs significantly better than the Gauss. fit or second moment method, and the second moment method performs better than the Gauss. fit except for small velocities. Figure 6 demonstrates that the MLE provides a velocity precision that even for a number of phase differences as low as 33 is in close agreement with the theoretical prediction. For a larger number of samples and flow in a capillary (Fig. 6(b)) the experimentally determined flow precision is poorer than predicted by theory, especially at lower velocities. We attribute this deviation to the precision with which the isotachs could be determined in the flow phantom. Close to the capillary wall where the flow velocity is low, the velocity gradient with radius is greatest, introducing the largest variation in velocity with voxel position.

Fig. 6 Analysis of the velocity estimation uncertainties for measurements (dots) compared to simulations (solid line) for the cases of 290 samples per estimation (a and b) and 33 samples per estimation (c and d). Flow is mimicked by a motion of the beam with respect to a static medium (a and c) and by pumping intralipid with a constant flow rate through a capillary (b and d). Three estimators were used (see 2. Theory & Simulations): Gauss. fit (blue), second moment (red) and MLE (black).

Download Full Size | PDF

Finally, we present the velocity precision as a function of required samples for both the Gauss. fit and the MLE. This allows selection of the number of samples required to reach a certain velocity precision. Figure 3(c), 3(d), Fig. 6(a)-6(d) all show that the velocity estimation uncertainty increases with relative velocity.

Figure 7 shows an additional simulation where we evaluated the velocity estimation uncertainty at a normalized velocity of 0.92 for the Gauss. fit and MLE as a function of the number of samples per estimation in a double logarithmic plot. The uncertainty was simulated from 8 samples to 400 samples per estimation in 40 steps. For uncorrelated data the velocity estimation uncertainty is expected to be inversely proportional to the square root of the number of samples per estimation. Linear functions with slopes of −0.5 were fit to the simulated uncertainties (dashed lines). The velocity estimation uncertainty follows closely the expected behavior, however for a low number of samples a small deviation can be seen. From the shift along the x-axis between the two fits a factor was derived that indicates by how much the number of samples can be reduced by using the MLE instead of the Gauss. fit achieving the same velocity estimation precision. The MLE achieves the same precision as the Gauss. fit with a factor of 8.9 less samples.

Fig. 7 Double logarithmic plot of the velocity estimation uncertainty for a set normalized velocity of 0.92 and varying samples per velocity estimation for the Gauss. fit method (blue) and the MLE (black). The dashed lines indicate a fit of a linear function with slope −0.5, indicating that the velocity estimation uncertainty is inversely proportional to the square root of the number of samples per estimation.

Download Full Size | PDF

5. Discussion

In the study the estimation precision of predominantly lateral flow velocities has been analyzed for an established method [21] and two proposed estimators. Both estimators have shown lower statistical fluctuations and therefore better precision for a large range of velocities in comparison to the established method. The MLE delivered the best results but already the use of the standard deviation as a statistical second moment represents a significant improvement. The latter is similar to approaches using phase variance for OCTA (e.g [35].) but a LUT as used in this work or calibration curve is necessary for accurate velocity estimations. A more practical demand on measurements, rather than the improvement of the precision, is that a specific precision is desired and ought to be reached with a minimal amount of measurements in order to keep the data acquisition to within clinical acceptable times. For this demand the more interesting aspect is by how much the number of samples can be reduced by choosing the optimal estimation method. This was tested for an exemplary case when the normalized velocity was set to a fixed value and a reduction in samples between the established method (Gauss. fit) and MLE by a factor of 8.9 was found. For in vivo measurements this can be especially helpful where short acquisition times are important to reduce the influence of motion artifacts.

The analysis presented here is based on phase sensitive velocimetry methods. Other methods have been shown, e.g. based on complex speckle decorrelation [19,22]. In the future we plan to extend the theoretical description and simulations to include the amplitude of the OCT signal and use it for a more general analysis that can be applied to methods using the full field of the OCT signal.

Velocities presented in this work are always reported as normalized displacements but the conversion to actual velocity is straight forward if $w_{0}$ and the time delay $τ$ over which phase differences are recorded are known. The determined beam overlap $δ r / w_{0}$ between neighboring pixels in a static medium was used to calculate the experimental beam radius as 9.59 μm, which was in good agreement to the theoretical diffraction limit of 8.98 μm. determined experimental flow velocity of the intralipid solution in the capillary was 864 mm/s. Assuming a parabolic velocity profile (laminar flow) the maximum velocity given by the set flow rate was 700 mm/s.

The feature of the MLE that the estimation uncertainty is approximately linearly proportional to the velocity itself means that the estimation uncertainty is a fixed percentage of the estimated velocity (up to a normalized velocity of 0.9-1.0). The displacement $δ r$ is a function of velocity $v$ and the time interval $τ$ $(δ r = v τ)$ which elapses between pairs of measurements from which phase differences are computed. A combination of time delays can be chosen to cover the full range of flow velocities in e.g., the retina, providing an estimation uncertainty which is a fixed percentage of the estimated velocity over the full range.

The optimization of time intervals between repeated scans for velocities expected in the human retina was demonstrated elsewhere [8,23] as well as an additional combination of different time intervals to extend the dynamic range of the velocity estimations. Alternatively, also the beam diameter could be changed.

The presented theory does not include any effects of particles of different sizes or the rotation of non-spherical objects. However, previous studies have shown that methods for velocimetry developed without including those factors can still remain translatable to measurements of blood particles [18,36]. This may still be important in future and require a higher level of numerical simulations but in this work we focused in the theoretical part on an analytical description of the distributions of the phase differences under influence of flow.

In the experiments presented here, the signal-to-noise ratio (SNR) was in the range of 30-35 dB. According to Park et al. [37] a SNR of 30 dB leads to a phase noise standard deviation of 0.022 rad and therefore approximately to a ratio $δ r / w_{0}$ of 0.022 which is significantly below all simulated and measured values. Therefore, we treated the effects to the velocity estimations in this work as negligible. However, for in vivo experiments this contribution will become more important as also areas with lower SNR will be present. Effects of low SNR to the determination of the Doppler frequency in phase sensitive OCT and joint spectral and time domain OCT have been investigated previously [38,39]. In such cases, noise from low SNR will change the PDF and therefore change the uncertainty of estimated velocities and the accuracy of the estimation because under such circumstances the measurements of phases contain partially information about decorrelation and partially about the SNR. Similar effects have earlier been treated in the design of an MLE for the Doppler frequency estimation [39].

6. Conclusion

The PDF was calculated describing the phase fluctuations expected for phasors reflected from a sample with axial and lateral flow. Using this PDF, the statistical limit for the estimation precision of phase based OCT velocimetry methods has been studied and used for the design of an estimator producing results approaching the Cramér-Rao lower bound. It has been compared to established techniques and an MLE based method using the derived PDF provided a significant improvement. The method is designed to be applicable to measurements of pairs of phases which gives it a larger versatility in comparison to methods e.g. restricted to M-scans. In future it is planned to extend this analysis to the full complex field for phasor based techniques.

Funding

The Dutch Technology Foundation STW (grant number 12822), which is part of the Netherlands Organisation for Scientific Research (NWO); the Ministry of Economic Affairs; the Health~Holland framework; the European Union's Horizon 2020 research and innovation program (654148); LaserLaB Europe and Heidelberg Engineering GmbH.

Acknowledgments

We thank Jelmer Weda for the preparation of phantoms and help with the implementation in the model eye.

Disclosures

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

References

1. M. H. M. Cuypers, J. S. Kasanardjo, and B. C. P. Polak, “Retinal blood flow changes in diabetic retinopathy measured with the Heidelberg scanning laser Doppler flowmeter,” Graefes Arch. Clin. Exp. Ophthalmol. 238(12), 935–941 (2000). [CrossRef] [PubMed]

2. J. E. Grunwald, T. I. Metelitsina, J. C. Dupont, G. S. Ying, and M. G. Maguire, “Reduced foveolar choroidal blood flow in eyes with increasing AMD severity,” Invest. Ophthalmol. Vis. Sci. 46(3), 1033–1038 (2005). [CrossRef] [PubMed]

3. R. Zeimer, “Nature is teaching us to be humble in our quest to measure structure and function in glaucoma,” Br. J. Ophthalmol. 91(1), 2–3 (2007). [CrossRef] [PubMed]

4. E. R. Muir, R. C. Rentería, and T. Q. Duong, “Reduced ocular blood flow as an early indicator of diabetic retinopathy in a mouse model of diabetes,” Invest. Ophthalmol. Vis. Sci. 53(10), 6488–6494 (2012). [CrossRef] [PubMed]

5. Q. Shao, F. M. Heussen, Y. Ouyang, and A. Hager, “Retinal vessel diameter changes in different severities of diabetic retinopathy by SD-OCT,” Eur. J. Ophthalmol. 26(4), 342–346 (2016). [CrossRef] [PubMed]

6. K. K. W. Chan, F. Tang, C. C. Y. Tham, A. L. Young, and C. Y. Cheung, “Retinal vasculature in glaucoma: a review,” BMJ Open Ophthalmol 1(1), e000032 (2017). [CrossRef] [PubMed]

7. A. H. Kashani, C. L. Chen, J. K. Gahm, F. Zheng, G. M. Richter, P. J. Rosenfeld, Y. Shi, and R. K. Wang, “Optical coherence tomography angiography: A comprehensive review of current methods and clinical applications,” Prog. Retin. Eye Res. 60, 66–100 (2017). [CrossRef] [PubMed]

8. 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] [PubMed]

9. A. S. Nam, I. Chico-Calero, and B. J. Vakoc, “Complex differential variance algorithm for optical coherence tomography angiography,” Biomed. Opt. Express 5(11), 3822–3832 (2014). [CrossRef] [PubMed]

10. T. E. de Carlo, A. Romano, N. K. Waheed, and J. S. Duker, “A review of optical coherence tomography angiography (OCTA),” Int. J. Retina Vitreous 1(1), 5 (2015). [CrossRef] [PubMed]

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

12. Z. 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] [PubMed]

13. Z. Chen, T. E. Milner, S. Srinivas, X. Wang, A. Malekafzali, M. J. C. van Gemert, and J. S. Nelson, “Noninvasive imaging of in vivo blood flow velocity using optical Doppler tomography,” Opt. Lett. 22(14), 1119–1121 (1997). [CrossRef] [PubMed]

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

15. 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] [PubMed]

16. 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] [PubMed]

17. 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] [PubMed]

18. 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] [PubMed]

19. 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 Stat. Nonlin. Soft Matter Phys. 88(4), 042312 (2013). [CrossRef] [PubMed]

20. H. Ren, K. M. Brecke, Z. Ding, Y. Zhao, J. S. Nelson, and Z. 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] [PubMed]

21. 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] [PubMed]

22. 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] [PubMed]

23. 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] [PubMed]

24. 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] [PubMed]

25. A. Lagendijk, M. P. van Albada, A. Lagendijk, MB van der Mark, and MP van Albada, “Light scattering in strongly scattering media: Multiple scattering and weak localization,” Phys. Rev. B Condens. Matter 37(7), 3575–3592 (1988). [CrossRef] [PubMed]

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

27. H. Kogelnik and T. Li, “Laser beams and resonators,” Appl. Opt. 5(10), 1550–1567 (1966). [CrossRef] [PubMed]

28. N. Weiss, T. G. van Leeuwen, and J. Kalkman, “Simultaneous and localized measurement of diffusion and flow using optical coherence tomography,” Opt. Express 23(3), 3448–3459 (2015). [CrossRef] [PubMed]

29. J. W. Goodman, Statistical optics, Wiley classics library ed., Wiley classics library (Wiley, New York, 2000), pp. xvii, 550 p.

30. 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.

31. 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] [PubMed]

32. M. J. Ju, Y. J. Hong, S. Makita, Y. Lim, K. Kurokawa, L. Duan, M. Miura, S. Tang, and Y. Yasuno, “Advanced multi-contrast Jones matrix optical coherence tomography for Doppler and polarization sensitive imaging,” Opt. Express 21(16), 19412–19436 (2013). [CrossRef] [PubMed]

33. S. Makita, Y. Hong, M. Yamanari, T. Yatagai, and Y. Yasuno, “Optical coherence angiography,” Opt. Express 14(17), 7821–7840 (2006). [CrossRef] [PubMed]

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

35. J. Fingler, D. Schwartz, C. Yang, and S. E. Fraser, “Mobility and transverse flow visualization using phase variance contrast with spectral domain optical coherence tomography,” Opt. Express 15(20), 12636–12653 (2007). [CrossRef] [PubMed]

36. 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] [PubMed]

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

38. 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] [PubMed]

39. 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] [PubMed]

Precision analysis and optimization in phase decorrelation OCT velocimetry

Abstract

1. Introduction

2. Theory & simulations

2.1 Theory

2.2 Simulations

3. Experimental setup

3.1 OCT-system

3.2 Phantoms

3.3 Data processing

4. Results

5. Discussion

6. Conclusion

Funding

Acknowledgments

Disclosures

References

Cited By

Figures (7)

Equations (23)

Biomedical Optics Express