Complex differential variance angiography with noise-bias correction for optical coherence tomography of the retina

Boy Braaf; Sabine Donner; Ahhyun S. Nam; Brett E. Bouma; Benjamin J. Vakoc

doi:10.1364/BOE.9.000486

1. Introduction

Angiography has become one of the primary imaging modes of optical coherence tomography (OCT) [1,2]. The high axial (2-15 μm) and lateral (5-30 μm) resolution of OCT angiography (OCTA) allows for three-dimensional visualization of vascular networks down to the capillary level [3]. Clinically, OCTA provides a non-invasive and label-free alternative to the more commonly used dye-based angiography techniques. Accordingly, OCTA has been adopted in fields such as ophthalmology [4,5], dermatology [6,7], neurology [8,9] and oncology [10,11]. In ophthalmology, OCTA has rapidly grown in popularity and has been integrated in several commercial devices for clinical investigation of the retinal vasculature [12].

The angiography signal in OCTA is derived from temporal changes in OCT data obtained from the same tissue location. Flowing blood causes the amplitude and phase of the complex OCT signal to change over time, while other tissues generate more time-stable complex signals. There are multiple methods to derive an angiography signature from either the OCT signal amplitude (or intensity) [13–15], the OCT signal phase [16–18], or from both simultaneously [8,19,20]. Because the phase signal in OCT is sensitive to instrument synchronization and timing errors, small axial displacements of the sample, and mechanical motion or vibrations, phase-based angiography has required methods to phase stabilize the measurements. These include optical clocking and triggering of the wavelength-swept laser and signal digitization [21,22] as well as post-processing methods that rely on a reference (mirror) signal for numerical phase correction of digitization errors [23–25]. In addition, stable galvanometer scanning patterns have been developed to minimize vibrations and positioning errors [18,26], while sample bulk motion is often numerically compensated with phase changes measured from stationary tissues [4,16].

Recently a phase-sensitive OCTA method has been developed based on a complex differential variance (CDV) algorithm that does not require phase stability [19]. The CDV algorithm obtains an angiography signature from the complex difference between A-scans acquired at two time points. Over a limited depth range within the complex differential A-scan, the effect of phase instabilities can be approximated by a common phase offset that can be ignored in the angiographic calculation. As such, OCTA performed with the CDV algorithm is highly immune to phase instabilities and does not require sophisticated solutions for phase stabilization. It was further shown that the CDV algorithm provides improved angiographic contrast over a set of amplitude-only angiography methods.

Like all angiographic methods, CDV can mistake measurement noise for angiographic signal. The dominant noise sources in phase-sensitive OCTA methods are shot noise and galvanometer scanning errors [23,27]. The latter causes positioning errors in the revisitation of the same location that lead to phase noise. As a consequence, CDV angiographic signals are correlated to the measurement signal-to-noise ratio (SNR) at low and moderate SNRs, and are limited by positioning errors for higher SNRs. This SNR noise-bias conflates sample structure and angiographic signal, resulting in images that are sometimes difficult to interpret.

In this paper the CDV algorithm is extended to correct for noise-bias. This is done by analytically deriving the contribution of the shot noise to the CDV parameter and the corresponding SNR normalization of the CDV measurement range. To correctly handle high SNR measurements, asymptotic limits on the noise-bias corrected CDV are included. The algorithm was implemented on an ophthalmic optical frequency domain imaging (OFDI) instrument with dynamic eye tracking. The improvement in angiography is shown by comparing noise-bias corrected and conventional CDV angiography, and by comparing noise-bias corrected CDV with phase-variance angiography. Finally, the eye tracking was used to acquire and compound multiple motion corrected data sets to obtain high quality en face angiograms of the retinal and choroidal vasculature.

2. Theory and CDV analysis

In this section, we first describe the origin of the noise bias in conventional CDV, present its impact on an angiographic cross-sectional image of the retina, and describe the principle of the noise-bias correction pursued in this work (section 2.1). Next, we present a generalized CDV signal model including shot noise effects (section 2.2) that is used to derive the upper and lower normalization bounds needed to correct the noise-bias of the CDV measurement range (section 2.3). The quality of the CDV images is further improved by the suppression of low-SNR signals and by spatial, temporal, and polarization averaging (section 2.4). Finally, practical implementation guidelines are presented for convenience (section 2.5).

2.1 Conventional CDV angiography and its noise-bias

The CDV algorithm is based on the calculation of complex differential A-scans by complex conjugate multiplication of OCT signals obtained at different time points from the same sample location [19]. The conventional implementation of CDV is given as:

f_{C D V c o n v} (z, x, t_{A}, t_{B}) = 1 - \frac{C (z, x, t_{A}, t_{B})}{I (z, x, t_{A}, t_{B})}

where f_CDVconv(z,x,t_A,t_B) is the CDV angiography value calculated between time points t_A and t_B at depth position z for lateral A-scan location x. C(z,x,t_A,t_B) describes the coherent average of the complex differential A-scans along depth for which the amplitude is given as:

C (z, x, t_{A}, t_{B}) = | (\frac{1}{2 L + 1}) \sum_{k = - L}^{L} w (k) \cdot E (z - k, x, t_{A}) \cdot E^{*} (z - k, x, t_{B}) | .

In Eq. (2) E(z,x,t) is the complex OCT signal at time point t, E(z,x,t_A)·E*(z,x,t_B) represents the complex differential A-scan, and w(k) is an averaging depth kernel that includes 2L + 1 pixels along the z direction. The coherent average is normalized by division with the incoherent average of the original A-scan data as given by:

I (z, x, t_{A}, t_{B}) = (\frac{1}{2 L + 1}) \sum_{k = - L}^{L} \frac{w (k)}{2} \cdot [{| E (z - k, x, t_{A}) |}^{2} + {| E (z - k, x, t_{B}) |}^{2}] .

In the case of stationary tissue the signal phasors (pixels) of the complex differential A-scan will be aligned, and the coherent average will result in a phasor that has the same length as the incoherent average. However, in the case of blood flow, motion will change the phase angles and amplitudes differently per pixel in the complex differential A-scan, and the coherent average phasor will decrease in length. The CDV data is normally displayed over a range from 0 (stationary tissue) to 1 (strong blood flow) encoded in gray scale from black to white respectively. In comparison with Nam et al. [19] the square-root transformation of the overall CDV function is left out of Eq. (1) for simplicity.

Intensity and CDV angiography images of the retina are presented in Fig. 1. The OCT intensity shows the layered structure of the retina in inverted gray scale. Although several tissue layers are clearly delineated only a few large blood vessels can be identified from their strong reflectivity as denoted by the red arrows. Figure 1(b) shows the corresponding angiogram calculated by the conventional CDV algorithm in which the large blood vessels have high CDV values. These vessels and several others are now easily identified, in part by their vertical shadow streaks that are cast to deeper depths. However, in Fig. 1(b), stationary tissue layers with a low SNR also show increased CDV values, which obscure a clear visualization of the smaller retinal vasculature. Additionally, outside the tissue noise fluctuations cause a false indication of blood flow and prevent a clear visualization of the choroidal vasculature below the retina.

Fig. 1 (a) OCT intensity image showing the structures of the retina. Several large blood vessels can be identified from the hyper-reflective spots within the upper layers of the retina as denoted by red arrows. (b) The retinal angiogram obtained with conventional CDV in which blood flow and stationary tissues are indicated in white and black respectively. The large blood vessels and several others are visible by their increased CDV value and their shadows cast onto deeper tissues. The noise-bias is seen as the SNR dependence of the conventional CDV method, which creates a layered appearance of the angiogram similar to the structural image. In addition, deep areas with limited SNR give a false indication of blood flow; this makes it difficult to appreciate the choroidal vasculature. Image sizes: 1.5 mm × 4.4 mm (height × width).

Download Full Size | PDF

In the CDV method, signal changes caused by blood flow are entirely described by the coherent average in the numerator of the fraction of Eq. (1). Although this numerator is mathematically normalized by the incoherent average of the denominator, this does not take the SNR or other noise sources into account. This leads to the CDV noise-bias as described above. To address this deficiency, this paper presents a CDV range normalization by SNR-dependent parameters as:

f_{C D V c o r} (z, x, t_{A}, t_{B}) = 1 - \frac{C (z, x, t_{A}, t_{B}) - {\hat{C}}_{f l o w} (z, x, t_{A}, t_{B})}{{\hat{C}}_{s t a t i o n a r y} (z, x, t_{A}, t_{B}) - {\hat{C}}_{f l o w} (z, x, t_{A}, t_{B})}

where Ĉ_stationary(z,x,t_A,t_B) and Ĉ_flow(z,x,t_A,t_B) describe the upper and lower limits of C(z,x,t_A,t_B) that correspond to the expected CDV signals for stationary tissue and flowing blood for the locally observed SNR. In contrast to conventional CDV, this normalization adjusts the floor and ceiling of the CDV to the local SNR to reduce the noise-bias in the angiography image. In the following sub-sections, we derive a statistical model including noise for C(z,x,t_A,t_B) and use this model to find expressions for Ĉ_stationary(z,x,t_A,t_B) and Ĉ_flow(z,x,t_A,t_B).

2.2 CDV statistical model

In order to understand the influence of noise on the CDV calculations a statistical model for the expected value of C(z,x,t_A,t_B) is derived from Eq. (2) and the definition of the complex OCT signals. Assuming shot noise limited measurements, the complex OCT signals E(z,x,t) are defined as complex phasors that are subject to additive noise:

E (z, x, t) = s (z, x, t) \cdot p (z, x, t) + n (z, x, t)

where s(z,x,t) is the true OCT signal phasor corresponding to the sample measurement, p(z,x,t) is a phasor with unity amplitude that describes the A-scan specific phase offset due to acquisition trigger instability or small-scale axial motion of the sample, and n(z,x,t) is a phasor that describes the shot noise. Subsequently the complex differential A-scan signal D(z,x,t_A,t_B) is calculated by complex conjugate multiplication between two OCT measurements obtained at time points t_A and t_B:

\begin{array}{l} D (z, x, t_{A}, t_{B}) = E (z, x, t_{A}) \cdot E^{*} (z, x, t_{B}) = s (z, x, t_{A}) \cdot s^{*} (z, x, t_{B}) \cdot p (z, x, t_{A}) \cdot p^{*} (z, x, t_{B}) \\ + s (z, x, t_{A}) \cdot p (z, x, t_{A}) \cdot n^{*} (z, x, t_{B}) \\ + s^{*} (z, x, t_{B}) \cdot p^{*} (z, x, t_{B}) \cdot n (z, x, t_{A}) \\ + n (z, x, t_{A}) \cdot n^{*} (z, x, t_{B}) . \end{array}

The complex differential A-scan signal of Eq. (6) consists of a signal phasor described by the first right-hand side term and a random noise phasor given by the sum of the other three terms that all include a shot noise component. The second and third right hand side terms of Eq. (6) both include a shot noise factor n(z,x,t) and can therefore be considered distributed as zero-mean complex Gaussians. The fourth term of Eq. (6) is mathematically described by a complex double-Gaussian distribution, which is approximated here by a normal (single) complex Gaussian distribution for simplicity. This approximation is possible in practice for the application to CDV angiography since this noise term is small compared to the other noise terms and its statistical variation can be described with tolerable error by the standard deviation of a normal (single) complex Gaussian distribution. Statistically, D(z,x,t_A,t_B) is therefore described by a complex Gaussian distribution with non-zero mean for which the magnitude can be modeled according to a Rice distribution [28]. Correspondingly, the expected value for C(z,x,t_A,t_B) is given by the mean amplitude of a Rice distribution, which is known to be described as:

\hat{C} (z, x, t_{A}, t_{B}) = σ_{C} (z, x, t_{A}, t_{B}) \cdot \sqrt{\frac{π}{2}} \cdot L_{1 / 2} (\frac{- S^{2} (z, x, t_{A}, t_{B})}{2 σ_{C}^{2} (z, x, t_{A}, t_{B})})

where S(z,x,t_A,t_B) describes the signal vector amplitude (i.e., the Rician noncentrality parameter), σ_C(z,x,t_A,t_B) describes the standard deviation of the Gaussian noise cloud (i.e. the Rician scale parameter), and L_1/2() is the Laguerre polynomial with order ½. For high SNR Eq. (7) can be approximated as: Ĉ(z,x,t_A,t_B) ≈S(z,x,t_A,t_B). Including the coherent average of Eq. (2), the signal vector amplitude S(z,x,t_A,t_B) is given as:

S (z, x, t_{A}, t_{B}) = | (\frac{1}{2 L + 1}) \sum_{k = - L}^{L} w (k) \cdot s (z - k, x, t_{A}) \cdot s^{*} (z - k, x, t_{B}) \cdot p (z, x, t_{A}) \cdot p^{*} (z, x, t_{B}) | .

Over the limited range of the averaging kernel low order phase instabilities (e.g., phase offsets and ramps) can be approximated by a common phase offset [19], and the phase instability terms p(z,x,t) can therefore be neglected:

S (z, x, t_{A}, t_{B}) = | (\frac{1}{2 L + 1}) \sum_{k = - L}^{L} w (k) \cdot s (z - k, x, t_{A}) \cdot s^{*} (z - k, x, t_{B}) | .

The random noise of the Rice distribution that is additive to S(z,x,t_A,t_B) results from the last three terms of Eq. (6). All include the shot noise phasor n(z,x,t) which is statistically described by an uncorrelated Gaussian noise variable with zero mean and an equal variance σ_n² for its real and imaginary components. Consequently the Rician random noise phasor also has a zero mean, and its statistical distribution can therefore be described by its standard deviation. Assuming Gaussian distributions for the individual noise terms of Eq. (6), the standard deviation for the noise on D(z,x,t_A,t_B) (a single pixel) is estimated in accordance with the summation of three independent random variables as:

σ_{D} (z, x, t_{A}, t_{B}) = \sqrt{{(| s (z, x, t_{A}) | \cdot σ_{n})}^{2} + {(| s (z, x, t_{B}) | \cdot σ_{n})}^{2} + σ_{AB}^{2}}

where the terms in the square-root describe the noise variances of the second, third and fourth right hand side terms of Eq. (6) respectively, and σ_AB represents the standard deviation of the complex double Gaussian distribution of the fourth right hand side term. In addition, the uniform phase distribution of the shot noise allowed the removal of p(z,x,t) and the phase of s(z,x,t) from Eq. (10). The standard deviation for the Rician noise cloud is then calculated by a summation of independent random variables over the pixels in the averaging kernel as:

σ_{C} (z, x, t_{A}, t_{B}) = \sqrt{\sum_{k = - L}^{L} {(\frac{1}{2 L + 1}) \cdot w (k) \cdot σ_{D} (z - k, x, t_{A}, t_{B})}^{2}} .

In Fig. 2 the coherent averaging process is schematically shown for stationary tissue as well as for an area with blood flow. Illustrated in Fig. 2(a) are five consecutive pixels in depth with a pair of complex OCT signals denoted by blue and green phasor arrows while red phasor arrows represent the shot noise. Although each pixel has a random phase offset, phase instabilities create a common phase angle change for all pixels between the complex OCTsignals. For stationary tissue this leads to aligned phasors in the complex differential A-scan. Blood flow however changes the phase angles and amplitudes over time differently per pixel and correspondingly the phasors of the complex differential A-scan are not aligned. In the coherent averaging over all pixels as shown in Fig. 2(b) the phasors of the complex differential A-scans are added and result in a long straight phasor path for stationary tissue, while for blood flow the phasor path is curved. The resulting phasor is therefore shorter for blood flow in comparison with stationary tissue, while its length is described by S(z,x,t_A,t_B). In both cases the shot noise results in a random noise cloud around the end of the signal phasor sum with standard deviation σ_C(z,x,t_A,t_B). This introduces a bias in the measured length for C(z,x,t_A,t_B) as shown by the dashed lines. This bias is inversely related to the SNR of the measurement, since the deviation induced by the shot noise is proportionally larger for weaker signals.

Fig. 2 (a) Complex A-scan pairs E(z,x,t_A) and E(z,x,t_B) with their complex differential A-scan D(z,x,t_A,t_B) for stationary tissue (left) as well as tissue with blood flow (right). The complex A-scans pairs are denoted as blue and green vector arrows while the shot noise is visualized as red vector arrows. (b) In the coherent averaging step of Eq. (2) the phasors of the complex differential A-scans are summed and result in a long phasor path for stationary tissue due to the well aligned phasors, while in the case of blood flow the path is curved and the resulting phasor is significantly smaller in length. In both cases the shot noise results in a random ‘noise cloud’ around the phasor path end and introduces a bias in the measured length compared to the pure signal length as shown by the dashed lines. Im and Re indicate the imaginary and real axes of the phasor space.

Download Full Size | PDF

2.3 CDV limits

The amplitude of the coherent average C(z,x,t_A,t_B) must be normalized to be interpretable. In this work, we propose normalization by factors Ĉ_stationary(z,x,t_A,t_B) and Ĉ_flow(z,x,t_A,t_B) that are calculated for each location as a function of the local SNR. In this section, we derive these limiting values and additionally include practical effects such as noise induced by galvanometer positioning errors and by the use of a limited number of pixels in the averaging process.

2.3.1 Estimating Ĉ_stationary(z,x,t_A,t_B)

Stationary tissue results in OCT signals that are constant over time, i.e., s(z,x,t_A) = s(z,x,t_B). In this case, the phase components cancel out in the calculation of the complex differential A-scan and the amplitude of an ideal Rician signal phasor can be defined from Eq. (9) as:

S_{i d e a l} (z, x) = (\frac{1}{2 L + 1}) \sum_{k = - L}^{L} w (k) \cdot {| s (z - k, x, t) |}^{2} .

Positioning errors in the revisitation of the same sample location however induce additional (differential) noise and will add a random phasor component for each pixel [23,27]. In general this additional noise is limited and leads to only a small reduction of S_ideal(z,x). The change in S_ideal(z,x) can therefore be approximated with a scalar factor α that describes the reduction from the ideal case as:

S_{s t a t i o n a r y} (z, x, α) \approx α \cdot S_{i d e a l} (z, x) .

The upper limit of C(z,x,t_A,t_B) is then given as:

{\hat{C}}_{s t a t i o n a r y} (z, x, t_{A}, t_{B}) = σ_{C} (z, x, t_{A}, t_{B}) \cdot \sqrt{\frac{π}{2}} \cdot L_{1 / 2} (\frac{- S_{s t a t i o n a r y}^{2} (z, x, α)}{2 σ_{C}^{2} (z, x, t_{A}, t_{B})})

where σ_C(z,x,t_A,t_B) is given by Eq. (11). The practical implementation of the estimation of α from experimental data is described in section 2.5.4.

2.3.2 Estimating Ĉ_flow(z,x,t_A,t_B)

The lower bound of C(z,x,t_A,t_B) results from imaging (fast) blood flow which temporally decorrelates the complex OCT signals. In this process scatterers are randomly displaced and therefore randomize the amplitude and phase of the complex OCT signals. The coherent averaging step of Eq. (2) can be viewed in this case as a two-dimensional random walk for the phasors of the complex differential A-scan. The steps of this random walk can be further characterized by a mean step length (amplitude) that is equal to S_ideal(z,x) and random directional (phase) angles. The amplitude of the Rician signal phasor can therefore be approximated as a scalar reduction from the ideal stationary tissue case as:

S_{f l o w} (z, x, β) \approx β \cdot S_{i d e a l} (z, x) .

Here, β represents the expected traveled distance for a two-dimensional random walk of phasors with unit average length that are weighted by the CDV averaging kernel. Calculating the expected traveled distance is a well-known statistical problem for which the solution converges to zero for a large number of steps, i.e. a large number of pixels in the CDV calculations. In CDV the averaging kernel is however limited to a small number of pixels for which these solutions do not apply. In practice it is therefore easier to obtain β via numerical simulation (see section 2.5.5). The random noise phasor that acts on S_flow(z,x,β) is calculated in accordance with Eq. (11) and the lower limit of C(z,x,t_A,t_B) is therefore given as:

{\hat{C}}_{f l o w} (z, x, t_{A}, t_{B}) = σ_{C} (z, x, t_{A}, t_{B}) \cdot \sqrt{\frac{π}{2}} \cdot L_{1 / 2} (\frac{- S_{f l o w}^{2} (z, x, β)}{2 σ_{C ​}^{2} (z, x, t_{A}, t_{B})}) .

2.4 CDV noise reduction

In order to improve the performance of the CDV angiography, several methods were implemented for noise reduction. These included spatial, temporal and polarization-diverse averaging and the suppression of low SNR signals.

2.4.1 CDV spatial and temporal averaging

In practice, the calculation of C(z,x,t_A,t_B) with Eq. (2) can be improved by additional spatial averaging along the lateral direction as well as by temporal averaging over multiple time differential measurements:

C (z, x) = (\frac{1}{M}) \sum_{t = 1}^{M} [(\frac{1}{2 O + 1}) \sum_{k = - O}^{O} y (k) \cdot C (z, x - k, t_{A}, t_{B})] .

In Eq. (17) y(k) represents an averaging kernel along the lateral direction and includes 2O + 1 pixels, while the temporal averaging includes an M number of time differential measurements. These additional averaging steps provide an improved noise reduction, but are applied after the coherent averaging step and therefore do not reduce the CDV noise-bias. In order to include the lateral and temporal averaging in the CDV calculations, Eq. (17) is applied to the result of Eq. (2). Similarly the lateral and temporal averaging can also be applied to Ĉ_stationary(z,x,t_A,t_B) and Ĉ_flow(z,x,t_A,t_B).

2.4.2 Suppression of angiographic signals from low SNR

The normalization strategy described by Eq. (4) adapts measured signals to SNR-dependent expected ranges. This works well for SNR values from which one could reasonably expect to make a meaningful measurement. For extremely low SNRs, angiographic measurements are however unreliable. In order to optimally suppress these low SNR signals in the CDV angiograms, a look-up table was defined that applied a posthoc correction to the corrected CDV signal f_CDVcor(z,x,t_A,t_B). This look-up table used the f_CDVcor(z,x,t_A,t_B) output and the OCT intensity signal SNR for a correction to an ideal CDV response as:

f_{C D V i d e a l} (z, x, t_{A}, t_{B}) = f_{C D V c o r} (z, x, t_{A}, t_{B}) \cdot ρ (S N R (z, x), f_{C D V c o r} (z, x, t_{A}, t_{B}))

where ρ(SNR(z,x), f_CDVcor(z,x,t_A,t_B)) is the look-up table function. In the noise-bias corrected CDV method, the CDV value and the SNR are decoupled, i.e., a pixel with a certain SNR can have any CDV value depending on the blood flow, and therefore the look-up table is two dimensional. The look-up table was defined to transform f_CDVcor(z,x,t_A,t_B) to a sigmoid curve response to gradually suppress the CDV values for low SNR signals. The 50% CDV signal cut-off was set at an SNR of 5 dB to effectively suppress all signals with a lower SNR. A sigmoid curve was chosen here since it provides a smooth suppression of low SNR signals that is less sensitive to noise fluctuations than in the conventional case of a step function for an SNR-based threshold of the OCT angiography data.

2.4.3 Polarization-diverse averaging

In polarization-sensitive OCT systems, multiple polarization channels are available for CDV analysis. In this case a combined CDV signal can be obtained by signal amplitude weighted averaging of the CDV signals over all polarization channels as:

\bar{f_{C D V} (z, x, t_{A}, t_{B})} = \frac{\sum_{p} {f_{C D V i d e a l} (z, x, t_{A}, t_{B}) \cdot {| s (z, x, t) |}^{2}}_{p}}{\sum_{p} {{| s (z, x, t) |}^{2}}_{p}}

where the subscript p denotes the polarization channel enumeration. The use of the complex OCT signal amplitude as a weight ensures that polarization channels that only contain noise do not contribute. The noise from low SNR polarization channels is therefore suppressed.

2.5 Implementation

Sections 2.1-4 describe the theory of the noise-bias corrected CDV method. In this section annotations are given in a step-wise order for its practical implementation.

2.5.1 Step-wise execution

The complex OCT data required to run the noise-bias corrected CDV method consists of multiple repeated B-scans that are obtained from the same sample location. These B-scans are then paired in various combinations to form multiple time differential pairs for the CDV calculation. For each time differential B-scan pair the coherent average C(z,x,t_A,t_B) is calculated using Eq. (2), from which the result is spatially and temporally averaged using Eq. (17) to obtain C(z,x). Afterwards, Ĉ_stationary(z,x,t_A,t_B) and Ĉ_flow(z,x,t_A,t_B) are calculated according to Eqs. (14) and (16) respectively, and averaged using Eq. (17). Finally C(z,x), Ĉ_stationary(z,x) and Ĉ_flow(z,x) are used to calculate f_CDVcor(z,x) with Eq. (4) for which pixels with a low SNR are suppressed by application of Eq. (18). The final CDV data $\bar{f_{C D V} (z, x)}$ is then obtained after polarization averaging using Eq. (19). For these calculations the parameters S_ideal(z,x), σ_C(z,x,t_A,t_B), α and β were calculated as described below.

2.5.2 Calculating S_ideal(z,x)

The ideal stationary tissue phasor S_ideal(z,x) is calculated according to Eq. (12) and uses |s(z,x,t)| as its input. The parameter |s(z,x,t)| is estimated from the OCT intensity which is calculated as the squared magnitude of the complex OCT signal and is approximated as [29]:

{| E (z, x, t) |}^{2} \approx {| s (z, x, t) |}^{2} + 2 \cdot σ_{n}^{2}

where the (random noise) cross-terms between s(z,x,t) and n(z,x,t) are neglected for simplicity. The intensity noise floor from an empty part of the OCT image is then used to determine σ_n. Subsequently |s(z,x,t)| is obtained from the square root after subtracting 2·σ_n² from |E(z,x,t)|².

2.5.3 Calculating σ_C(z,x,t_A,t_B)

The noise standard deviation σ_C(z,x,t_A,t_B) is calculated according to Eqs. (10) and (11) using the parameters |s(z,x,t)|, σ_n, and σ_AB. The determination of the first two parameters was described in section 2.5.2. The parameter σ_AB describes the standard deviation of the complex double Gaussian distribution of the noise floor of the complex differential A-scan data. It can therefore be obtained as the complex standard deviation of the complex differential data from an empty part of the OCT image.

2.5.4 Estimating α

The scalar factor α describes the reduction of S_ideal(z,x) due to positioning errors in repeatedly measuring the same sample location. The parameter α is best measured repeatedly throughout a data set since due to varying galvanometer scanner performance or lateral sample motion the positioning error and thus α can fluctuate in time. The factor α can be estimated from experimental data by selecting pixels with high SNR (>10dB) for which it can be assumed that Ĉ_stationary(z,x) ≈S_ideal(z,x) in case of an ideal measurement (α = 1). The fraction between C(z,x) and S_ideal(z,x) will now provide α for data from stationary tissue as the ceiling value of the fraction over all data.

2.5.5 Estimating β

The scalar factor β describes the reduction of S_ideal(z,x) due to OCT signal decorrelation from blood flow. β can be obtained from the two-dimensional random walk statistics over a limited number of steps, but the analytic derivation is complicated. It is therefore more practical to use numerical simulation to estimate β. For this purpose E(z,x,t) was simulated for each pixel by the coherent sum over depth weighted by the axial point-spread function for a sub-resolution particle density of 100 / pixel. The electric field measured from each sub-resolution particle was modeled with a uniform random distributed phase and a Rayleigh random distributed amplitude [30]. Shot noise was added as a random complex phasor with a complex Gaussian distribution according to a 37 dB intensity SNR. The simulation included 100,000 pixels for both E(z,x,t_A) and E(z,x,t_B) from which C(z,x,t_A,t_B) and S_ideal(z,x) were calculated following Eq. (2) and Eq. (12) respectively. The factor β was then calculated by the fraction of C(z,x,t_A,t_B) and S_ideal(z,x) as a function of the window size of a Gaussian CDV averaging kernel w(k). In Fig. 3 the average β over all pixels is given and shows that β becomes smaller when w(k) includes more pixels in depth. In addition it can be seen that for the limited lengths of w(k) that are normally used, β will not reach zero and the value for Ĉ_flow(z,x) will be significant. This shows that by choosing a limited length for w(k) and thus maintaining axial resolution in the CDV image, the sensitivity for observing CDV signal is limited. If unaccounted for, as is the case with the conventional CDV method, this can lead to a weak appearance of flow signals in the CDV data. In practice the value of β is chosen in accordancewith the length for the axial Gaussian window w(k), while w(k) is chosen small enough to retain the axial resolution to distinguish stacked blood vessels along the axial direction.

Fig. 3 The scalar factor β as a function of the length of Gaussian window w(k). The blue data points are the estimated values for β obtained with the simulation. The red curve is a polynomial fit to these data points to show the trend in the decay of β with an increasing length of w(k). It can be seen that β becomes smaller when w(k) includes more pixels in depth, but does not reach zero. The lower CDV limit Ĉ_flow(z,x) will therefore have a significant value.

Download Full Size | PDF

3. Experimental setup and angiography scanning protocol

In this study a polarization-sensitive OFDI (PS-OFDI) system was developed based on a passive-component depth-multiplexed polarization-sensitive OCT design (section 3.1) [31–33]. The PS-OFDI setup was integrated with the ophthalmic interface of a commercial Heidelberg Engineering Spectralis OCT device which provided eye tracking and fixation capabilities (section 3.2) [34]. Angiographic imaging was performed with the combined setup using a bi-directional segmented (back-stitched) scan protocol that imaged the same lateral location on the retina multiple times (section 3.3) [10,18]. The optical layout of the setup and an overview of the scan protocol are provided in Fig. 4.

Fig. 4 Experimental setup and angiography imaging protocol. (a) PS-OFDI interferometer based on a passive-component depth-multiplexed PS-OCT design. (b) Ophthalmic interface for retinal imaging. The optical paths are indicated in red for OFDI, in orange for SLO and in blue for the fixation target. (c) Angiography imaging protocol. The bidirectional triangle waveform is plotted for two segments in blue on top of the corresponding OCT intensity image of retinal structures acquired near the optic nerve head. Red and green dashed boxes denote two waveform turning points that are magnified in the inset images. In the inset images the orange line sections indicate the position offset for the backward ramp sections of each triangle wave to match their measurement positions with the forward ramp sections. Component abbreviations: FBG: fiber Bragg grating, PC: polarization controller, PBC: fiber-based polarization beam combiner, BS: beam splitter, LP: linear polarizer, PBS: polarizing beam splitter, GP: glass plate, GS: galvanometer scanner, DM: dichroic mirror, FT: fixation target, LD: laser diode.

Download Full Size | PDF

3.1 PS-OFDI interferometer

The PS-OFDI setup was based on a 1 μm wavelength-swept source (Axsun Technologies, Inc., MA, USA) with an A-scan repetition rate of 100 kHz and an average output power of 21 mW. The light from the swept-source was coupled into a single-mode fiber-based interferometer as shown in Fig. 4(a). A fiber-coupler with a 99/1 ratio was used to sample onepercent of the light to create an optically derived A-scan trigger signal. For this purpose a fiber Bragg grating (FBG, O/E land Inc., QC, Canada) was used to generate a narrow-band reflection at 1050 nm that was sent via a 50/50 fiber coupler to a photodetector (Trigger, PDB420C, Thorlabs Inc., NJ, USA) that triggered the acquisition of every A-scan. The remaining 99 percent of the source light was split over the interferometer sample and reference arms in a 90/10 ratio. In the sample arm, two incident polarization states were created by splitting the light equally using a 50/50 coupler, sending it through air paths of different lengths, and recombining them in a fiber-based polarization beam combiner (PBC, OZ Optics, ON, Canada). Polarization controllers (PC) were used to orient the polarization state of each path to the crystal axis on the PBC inputs. This created two orthogonal polarization states that were sent simultaneously to the eye for which the air path length difference multiplexed these two signals in depth in the OFDI image. A 20/80 coupler sent 20% of this light to the ophthalmic interface for imaging and passed 80% of the returning back-reflected sample light to a polarization-sensitive detection unit. In the reference arm the path length and chromatic dispersion were matched to the sample arm by transmitting the light through an open-air path. The light from the reference arm was recombined with the sample arm in a bulk optics polarization-sensitive detection unit. First the reference arm light was sent through a linear polarizer (LP) with a 45° orientation to provide equal power to both polarization-diverse detection channels. The light from both arms was recombined in a broadband beam splitter (BS) after which polarizing beam splitter (PBS) cubes split the light from both outputs into its orthogonal polarization components. The horizontal and vertical polarization components were recorded by separate balanced detectors (PDB460C, Thorlabs Inc.), denoted in Fig. 4 as H-det. and V-det. respectively. As such both detection channels recorded two depth-multiplexed OFDI images corresponding to the two incident polarization states. Together these four OFDI images represented the four Jones-matrix components that are necessary for polarimetry analysis. In one arm of the H-det. a microscope cover glass slide (GS) was placed to create a calibration signal that was used for phase-stabilization of the data-acquisition.

After detection, the signals from H-det. and V-det. were recorded by a 14-bit resolution dual-channel data-acquisition (DAQ) board (PX14400D, Signatec, IL, USA). Every A-scan was digitized with 2048 samples for an optical bandwidth of 91 nm centered at 1040 nm. The recorded signals were mapped to a linear wavenumber space [35], phase-stabilized [23,24], numerically corrected for chromatic dispersion [36], apodized with a hamming window [37], and background signals were subtracted [24] following previously developed methods. The system sensitivity was measured to be 92 dB from a mirror by calculating the incoherent intensity sum of the four polarization channels. The setup had a −6 dB signal roll-off with depth over 5.0 mm with an axial resolution of 10 μm in air.

3.2 Ophthalmic interface with eye tracking

The ophthalmic interface was provided by a commercial Heidelberg Engineering Spectralis OCT device (Heidelberg, Germany) that was modified with custom optics in the 1 μm wavelength range as shown in Fig. 4(b). The PS-OFDI setup replaced the native Spectralis OCT hardware and its light (red paths) was directed via a set of galvanometer scanners (GS_OCT) and a telescope to the human eye. The OFDI beam incident on the cornea had an optical power of 1.6 mW and provided a diffraction-limited spot-size of 18 μm on the retina. Dichroic mirrors (DM) were used to combine the optical paths of the Spectralis’ fixation target (FT, blue paths) and scanning laser ophthalmoscope (SLO, orange paths) with the PS-OFDI optical path. The FT consisted of a blue light LED grid and provided a fixation stimulus to improve the subject’s gaze stability as well as coarse alignment functionality to a specific retinal location. SLO imaging was performed simultaneous to the PS-OFDI imaging with the Spectralis built-in Heidelberg Retinal Angiograph from which the 815 nm infrared imaging channel was used. The SLO used a laser diode (LD) from which the light was sent in transmission through a beam splitter (BS). A set of scanners (GS_SLO) provided high-speed raster scanning independent from the OFDI imaging. In the return pass the reflected SLO light from the retina was (partially) reflected by the BS to a photodetector (Det.) that included a confocal pinhole. The SLO had a maximum optical power of 80 μW on the cornea and provided a diffraction-limited spot-size of 15 μm on the retina. The combined operation of the PS-OFDI setup and the SLO instrument was safe in its use according to the ANSI 2007 laser safety standard [38]. In this study the SLO was operated at a 15.4 Hz frame rate for a 15° field-of-view (FOV). Eye motion was measured in real-time from the SLO images by the Spectralis software by analyzing motion-induced affine image transformations compared to a reference SLO image. The obtained eye motion was converted into a correction signal and imported into the PS-OFDI acquisition control software for on-the-fly correction of the OFDI galvanometer waveforms in order to follow eye motions during the acquisition.

3.3 Angiography imaging protocol and processing

Angiographic imaging of the retina was performed with a scan protocol as visualized in Fig. 4(c). This scan protocol repeatedly scans the same locations on the retina using a bidirectional segmented triangle waveform that enables inter-B-scan comparison for angiographic detection of blood flows down to the capillary level [10,18]. In the upper panel of Fig. 4(c) an example of two waveform segments is shown that were measured from the retina near the optic nerve head. The waveform is plotted in blue and shows how the lateral position on the retina changed as a function of time. The corresponding OCT intensity image is shown on the background to emphasize that the same retinal structures were measured multiple times. In this study each waveform segment was configured to scan the same lateral locations five times resulting in two-and-a-half triangle waves per segment. The final half triangle wave of each segment forms the connection to the next segment by concatenating the same waveform with a (lateral) position offset.

As demonstrated by Duma et al. [26] the galvanometer scanner motion can deviate from a triangle waveform due to inertia of the scanner, especially at the turning points in the waveform. This effect causes non-uniform lateral sampling of the retina and results in distortions or missing data in the final angiographic images. As a solution a small lateral position overlap in between adjacent segments is therefore included [18]. In Fig. 4(c) red and green dashed boxes denote the positions for two turning points of connecting waveform segments, and are shown magnified as the inset images in the lower panels. The position overlap between the adjacent segments is clearly seen in the inset images as the position (y-axis) overlap between the waveform curves. Any surplus A-scans acquired in the overlap regions were omitted from the final reconstructed B-scan images.

In angiographic OCT imaging it is important to precisely position repeated measurements at the same location to minimize phase noise [18]. In case of a bidirectional triangle waveform this means that the sampling positions on the forward and backward scanning ramps of each triangle need to be matched exactly. This was achieved by defining each ramp section of a triangle with the same integer number of A-scan sampling positions. However, a time lag between the commanded position by the waveform and the actual position of the galvanometer mirror resulted in a shift and therefore a mismatch in the sampling positions between the two triangle ramps. Rough matching was obtained by circular shifting the array of waveform sampling points by an integer number, while fine matching was achieved by implementing a small position offset for the backward triangle ramp sections. The position offset for the backward ramp section is indicated in the inset images of Fig. 4(c) by the orange line segments. The waveform circular shift and position offset were found to be constant values of 17 A-scans and 0.25 A-scans respectively, which indicated that the time lag was a fixed delay in the galvanometer mirror response.

The angiography scan protocol used in this study was configured with segments for which each (half triangle) ramp consisted of 273 A-scans. Due to the bidirectional design of the waveform the time sampling for the five repeated measurements is not the same for every position in a waveform segment. The middle position of a segment experiences a constant time difference between successive measurements, while near the waveform turning points the time difference alternates between short and long. This property makes the bidirectional segmented triangle waveform suitable for angiography as long as the average time difference between the measurements provides enough sensitivity to observe blood flow down to the capillary level. The current configuration provides an average time difference of 2.73 ms, which was found sufficient for capillary imaging in the retina during previous studies [18].

Two scan protocol configurations were used in this study. A wide-field imaging protocol imaged the retina over a field-of-view of 8.8 mm (29°) with 8 segments that included 5 repeated measurements per location. The wide-field bidirectional B-scan included 11264 A-scans (including fly-backs) and was in post-processing deinterleaved into 5 normal B-scans of 2000 A-scans. Volumetric data sets were obtained with 500 B-scan locations and the acquisition time for a single data set was 56 seconds during which the SLO eye tracking stabilized the OFDI imaging onto the moving retinal target. Alternatively a small-field imaging protocol was used with a field-of-view of 4.4 mm (15°) using 4 segments with 5 repeated measurements. In this configuration 1000 B-scan locations were defined per data set with a total acquisition time of 56 sec. A small-field bidirectional B-scan included 6144 A-scans and was in post-processing deinterleaved into 5 normal B-scans of 1000 A-scans.

Large axial bulk sample motion was minimized numerically as described by Braaf et al. [18]. The axial displacement was estimated with cross-correlation from the (sub-)pixel shift in between the intensity images of the compared B-scans. The found shift in pixels was corrected by multiplication of the OCT signals with a complex phase ramp. No additional bulk motion phase correction was applied.

CDV angiography images were calculated by combining the 5 repeated B-scans into 10 time differential pairs. The CDV algorithm was configured with Gaussian averaging kernels having an axial size of 7 pixels (1.2 pixels std.) and a lateral size of 9 pixels (1.6 pixels std.) and the parameter β was set accordingly at a value of 0.66. In addition, OCT intensity images were calculated from the incoherent intensity sum over all B-scans.

4. Results

The measurements of human retinas in vivo adhered to the tenets of the Declaration of Helsinki and were approved by the Massachusetts General Hospital Institutional Review Board. Informed consent was obtained from the imaged subject.

4.1 Cross-sectional CDV angiography: improvements by noise-bias correction

CDV angiography imaging was performed on a healthy volunteer in the macula using the wide-field imaging protocol. Horizontal cross-sectional images through the fovea are shown in Fig. 5. In Fig. 5(a) the OCT intensity shows the tissue layers of the retina and the choroid below it. Several large blood vessels can be identified by their strong reflectivity and are indicated by red arrows. In Fig. 5(b) an angiogram obtained with conventional CDV is shown. As described for Fig. 1(b), the conventional CDV angiogram has a layered appearance due to its SNR dependence. This makes it difficult to observe smaller vasculature in the retina besides the large vessels that are indicated by the red arrows. The erroneous indication of blood flow outside the tissue further hampers a clear visualization of the choroidal vasculature. In Fig. 5(c) the noise-bias corrected CDV angiogram is shown in which these problems are resolved. In this angiogram all the stationary tissues have equally low CDV values. The small vasculature in the retina can therefore be seen in cross-section as small white dots against a black background of stationary tissue. The choroidal vasculature below the retina can now clearly be delineated from the retina and the background. In Fig. 5(d) the suppression of noise from lateral positioning errors was switched off by setting a constant value of α = 1 to show its impact on the noise-bias corrected CDV image. Additional noise is observed in this case as an increase in the background noise throughout the CDV image, which makes it more difficult to identify individual capillaries. The severity of the noise is also not uniform throughout the image, which makes the application of a simple threshold difficult. This emphasizes the importance of properly addressing this noise in the calculation of the CDV image as was described in section 2.3. In Figs. 5(e) – (h) zoomed in sections of Figs. 5(a) – (d) are shown as denoted by the blue frames in the latter. These zoomed in figures clearly show that it is hard to identify retinal capillaries in either the OCT intensity or conventional CDV angiography images, while they are clearly visible in the noise-bias corrected CDV images. It can further be seen that the suppression of noise from positioning errors in Fig. 5(g) results in a much clearer observation of the retinal capillaries than without this suppression as shown in Fig. 5(h).

Fig. 5 Cross-sectional CDV angiography in the macula. (a) OCT intensity showing the structures of the retina. (b) Conventional CDV angiogram. This angiogram suffers from noise-bias artifacts that hide all but the larger vessels denoted by red arrows. (c) Noise-bias corrected CDV angiogram. Correction of the noise-bias suppressed the artifacts of the conventional CDV angiogram and enabled a clear observation of the smaller retinal and the choroidal vasculature. (d) Noise-bias corrected CDV angiogram without the suppression of positioning errors. This lead to increased background noise, which made it hard to observe individual retinal capillaries. Image sizes: 1.5 mm (height) × 8.8 mm (width). (e)-(h) Zoomed in sections of Figs. (a)-(d) as denoted by the blue frames. Zoomed image sizes: 0.65 mm (height) × 1.1 mm (width).

Download Full Size | PDF

4.2 Cross-sectional CDV angiography: artifact immunity

Nam et al. [19] demonstrated that conventional CDV is immune to phase-instability trigger distortions that otherwise creates vertical artifact lines in phase-sensitive angiography images. CDV is similarly insensitive for sample bulk axial motion phase artifacts as is shown in Fig. 6 for noise-bias corrected CDV in comparison to phase-variance angiography as described earlier by Braaf et al. [18]. Figure 6(a) shows the OCT intensity with the retinal structures superior from the fovea. In most phase-based angiography implementations it is necessary to apply a sample bulk motion correction of the phase to mitigate axial motion artifacts in the angiograms. In Fig. 6(b) the phase-variance angiogram is shown computed with the same data that is used as the input for CDV. Bulk axial motion (indicated by red arrows) creates oscillating artifacts with increased phase values that obscure the vasculature. After application of a bulk axial motion phase correction most of these artifacts are reduced as seen in Fig. 6(c). However, locations with lateral motion (indicated by green arrows) still show an increased background noise. In areas next to large arteries the pulsation of the cardiac cycle can compress the tissue and create a blood vessel mimicking artifact (indicated by the blue arrow). It is clear from the intensity image that no large blood vessel is present at this location. In Fig. 6(d) the noise-bias corrected CDV angiogram is shown which does not suffer from these artifacts. The suppression of positioning errors reduced the increased backgroundnoise from small lateral motions (indicated by the green arrows). Local expansion or compression of tissue mostly induces phase offsets in the complex OCT signals similar to axial bulk motion. CDV angiography is therefore also immune for this artifact and shows stationary tissue at this location.

Fig. 6 Noise-bias corrected CDV angiography compared to phase variance angiography. (a) OCT intensity. (b) Phase variance angiography without axial bulk motion phase correction. Several motion artifacts are seen the oscillatory phase increases indicated by red arrows. (c) Phase variance angiography after axial bulk motion phase correction. Artifacts can still be seen such as an increased background noise from lateral motion (green arrows) and a local displacement of the tissue due to arterial pulsation (blue arrow). (d) Noise-bias corrected CDV angiogram. The artifacts observed with phase variance angiography are suppressed in the CDV method. Image sizes: 1.5 mm (height) × 8.8 mm (width).

Download Full Size | PDF

4.3 En face CDV angiography

Wide-field en face angiographic imaging was performed on a healthy volunteer. The real-time compensation of eye motion during and between acquisitions allowed the compounding of multiple data sets. Ten data sets were obtained from the same retinal area and individually processed for CDV angiography. Afterwards the data sets were segmented to detect the retinal surface and the retinal pigment epithelium [39]. The segmentation results were used to generate separate en face angiography maps for the retina and the choroid by depth-integration of the CDV images. Individual en face angiograms included residual motion artifacts that occurred during saccadic eye motion, which created decorrelated segments. Figure 7(a) shows a single retinal en face angiogram including a time trace of the detected motion during the acquisition in a graph on the left. In this figure the decorrelated segments are seen as white horizontal line segments which distorted ~8% of the data set on average. The recording of the eye position during the acquisition with the SLO eye tracker allowed for the detection of residual motion artifacts in post-processing. A threshold of 15 μm was used on the instantaneous position change in the motion traces to discard motion distorted angiogram segments. Remaining artifacts were discarded using a threshold on the total integrated CDV signal. In Fig. 7(b) the retinal en face angiogram is shown after all the distorted segments were discarded. Due to the random occurrence of saccadic motion other data sets could be used to fill up the missing segments. However, in order to compound multiple data sets first residual eye motion shifts were detected and removed by cross-correlation of each B-scan (each horizontal en face line) onto a reference image. This step improved the location stability of smaller vasculature, which was necessary to avoid blurring of small features in the compounding process. The reference image was calculated as the average retinal en face map over all data sets, which provided a (relatively) motion free map of the larger vasculature. Angiograms with minimal motion artifacts of both the retina and the choroid were then obtained from the median pixel value over all motion-corrected data sets (shown in Fig. 8).

Fig. 7 Motion in retinal angiograms. (a) Motion distorted retinal en face angiogram. Saccadic motion decorrelated several line segments which appear completely white thereof. On the left a graph shows the horizontal (blue) and vertical (red) position of the retina as recorded by the SLO eye tracker. A threshold on the position change was used to reject the majority of the motion artifacts. (b) Retinal en face angiogram after the segments with motion artifacts were discarded. Image sizes: 8.8 mm × 8.8 mm.

Download Full Size | PDF

Fig. 8 Wide-field CDV angiography obtained with eye tracking for motion-correction and data set compounding (10 × ). (a) Retinal angiogram showing the vasculature from the large central arteries and veins down to the capillary level. (b) Choroidal angiogram showing a dense vascular network in which the larger vessels can be distinguished from their decrease in signal. Image sizes: 8.8 mm × 8.8 mm (29° × 29°).

Download Full Size | PDF

In Fig. 8 wide-field compounded angiograms of the retina and choroid are shown. These high-quality images with minimal motion and decorrelation artifacts show the benefit of performing CDV angiography together with SLO eye tracking. In Fig. 8(a) the compounded retinal angiogram shows clearly how the larger vessels are positioned around the macula along arcs starting from the optic nerve head on the right side. Their branches can be followed towards the fovea and despite the relatively sparse sampling along the vertical direction (18 μm steps) parts of the capillary beds can be seen. The compounded choroidal angiogram is shown in Fig. 8(b) and shows several large vessels in the periphery that are radially oriented towards the fovea. In the fovea a dense mesh of vessels is observed. Note that in the choroidal angiogram a part of the large vessels is observed from their reduced signal compared to surrounding vasculature. This is caused in part by the low SNR of the signals obtained from within the large choroidal blood vessels and the subsequent suppression of their CDV signals.

In Fig. 9 small-field compounded angiograms are shown of the retina centered at the fovea. The small-field angiography data sets were obtained similarly as described above forthe large-field data sets. In order to efficiently visualize the different retinal capillary networks, the data was segmented in a superficial and a deeper retinal layer by integrating the CDV angiography signals from the inner limiting membrane to the inner plexiform layer and from the inner nuclear layer to the outer plexiform layer respectively. Figure 9(a) shows the superficial retinal vasculature and clearly demonstrates the dense capillary network that interconnects the larger arterioles and venules. In addition in the center of the angiogram the foveal avascular zone is clearly visualized with several long single capillaries crossing on the side of this area. Figure 9(b) shows the deeper retinal vasculature with its capillaries organized in lobular patterns around the fovea. Although the superficial and deeper vascular layers are connected they clearly show the branching into different network patterns.

Fig. 9 Small-field CDV angiography of the retina with eye tracking for motion-correction and data set compounding (10 × ). (a) Superficial retinal vasculature with capillaries interconnecting the arterioles and venules. (b) Deeper retinal vasculature showing lobular patterns in the capillary network. Image sizes: 4.4 mm × 4.4 mm (15° × 15°).

Download Full Size | PDF

5. Discussion

The CDV angiography method presented in this paper offers several advantages. It is designed to be inherently insensitive to the triggering and clocking phase-instabilities that can be present in OFDI systems. In addition, when properly normalized, the CDV method is minimally affected by small-scale axial and lateral sample motion. These advantages make the CDV method cross-platform compatible between spectrometer and swept source-based OCT systems, applicable to samples within a wide range of bulk motion, while retaining its angiographic detection sensitivity.

The integration of an SLO eye tracking system with the OFDI system minimized eye motion artifacts in the angiography data sets and allowed the recording of multiple data sets to replace distorted data. This enabled the acquisition of large densely sampled data sets for angiographic evaluation of the retina over a wide field-of-view although with significant acquisition times. Multiple data sets from the same retinal area could further be compounded for image enhancement. This provided the angiographic imaging quality that is desired in the clinic for the investigation of ocular pathologies over large retinal areas. In the future it is expected that with the implementation of novel high-speed swept laser designs [40,41] and eye tracking technologies [39] these wide-field angiography images can be obtained with significantly shorter acquisition times.

Recent work has shown that combining OCTA with polarization-sensitive OCT can provide promising additional insights in retinal pathology [42]. PS-OFDI systems however often suffer from a significant loss in sensitivity (i.e. SNR) compared to conventional OFDI systems due to the distribution of the interferometric signal over multiple incident states and detection channels. In general it is therefore challenging with these systems to obtain high-quality OCTA data. In this study, this difficulty was overcome by implementation of several data averaging procedures into the CDV angiography method and the implementation of SLO eye tracking to support data set compounding. The current work is therefore promising in providing both OFDI-based angiography and polarimetry analysis in high quality.

The CDV method presented in this paper was designed and described with well-defined signals in mind that have a significant SNR, but also works for pixels that only contain shot noise. In the latter case the dominant shot noise causes C(z,x,t_A,t_B) and Ĉ_stationary(z,x,t_A,t_B) to become equivalent in value, and hence f_CDVcor(z,x,t_A,t_B) to become zero. This convenient side effect is further strengthened by the suppression of low SNR signals by Eq. (18), and results in an efficient suppression of the CDV background noise.

Speckle is inherently present in OCT data and is mostly noticeable from pixels wherein destructive interference results in poor SNR. In general, speckle will alter the SNR of an individual pixel in an OCT image by a change in the amplitude of the true OCT signal phasor s(z,x,t) of Eq. (5). In the CDV calculations, speckle therefore results in amplitude changes in the coherent averaging for the signal phasor S(z,x,t_A,t_B) of Eq. (9). The effect of speckle on the CDV method can be therefore be described as an amplitude weighting effect on the pixel level during the coherent averaging step of Eq. (2). In order for the CDV method to be minimally affected by pixels with speckle-degraded SNR, the length of the CDV depth kernel w(k) should be chosen longer than the axial size of a single speckle to increase the likelihood that pixels with sufficient SNR are included.

6. Conclusions

In this paper a novel CDV angiography method was demonstrated that corrected the noise-bias that otherwise could lead to false indications of blood flow in angiography images. This resulted in clearer observation of the retinal and choroidal vasculature down to the capillary level. The inherent immunity of the CDV method against phase instabilities was further shown to suppress sample motion artifacts. In addition, the application of eye tracking enabled the acquisition and compounding of angiography volumes for motion artifact free en face visualization. CDV combined with eye tracking can lead to more reliable angiography of the retina and choroid, which will be beneficial for the future clinical investigations of ocular pathologies.

Funding

Center for Biomedical OCT Research and Translation through Grant Number P41EB015903 awarded by the National Institute of Biomedical Imaging and Bioengineering and Grant Number R01CA163528 awarded by the National Cancer Institute of the National Institutes of Health. Additional support was provided by Heidelberg Engineering.

Acknowledgments

The authors like to thank Ellen Ziyi Zhang for support on the OFDI hardware development and Jacqueline Namati for her work on the study protocol and the management of the IRB regulatory process.

Disclosures

BB: (F), SD: (E), ASN: no conflict of interest, BEB: (F,P), BJV: (F,P)

References and links

1. 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]

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

3. R. A. Leitgeb, R. M. Werkmeister, C. Blatter, and L. Schmetterer, “Doppler optical coherence tomography,” Prog. Retin. Eye Res. 41, 26–43 (2014). [CrossRef] [PubMed]

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

5. K. V. Chalam and K. Sambhav, “Optical coherence tomography angiography in retinal diseases,” J. Ophthalmic Vis. Res. 11(1), 84–92 (2016). [CrossRef] [PubMed]

6. U. Baran, W. J. Choi, and R. K. Wang, “Potential use of OCT-based microangiography in clinical dermatology,” Skin Res. Technol. 22(2), 238–246 (2016). [CrossRef] [PubMed]

7. M. Ulrich, L. Themstrup, N. de Carvalho, M. Manfredi, C. Grana, S. Ciardo, R. Kästle, J. Holmes, R. Whitehead, G. B. Jemec, G. Pellacani, and J. Welzel, “Dynamic optical coherence tomography in dermatology,” Dermatology (Basel) 232(3), 298–311 (2016). [CrossRef] [PubMed]

8. R. K. Wang and S. Hurst, “Mapping of cerebro-vascular blood perfusion in mice with skin and skull intact by Optical Micro-AngioGraphy at 1.3 µm wavelength,” Opt. Express 15(18), 11402–11412 (2007). [CrossRef] [PubMed]

9. V. J. Srinivasan, A. C. Chan, and E. Y. Lam, “Doppler OCT and OCT angiography for in vivo imaging of vascular physiology,” in Selected Topics in Optical Coherence Tomography (InTech, 2012).

10. B. J. Vakoc, R. M. Lanning, J. A. Tyrrell, T. P. Padera, L. A. Bartlett, T. Stylianopoulos, L. L. Munn, G. J. Tearney, D. Fukumura, R. K. Jain, and B. E. Bouma, “Three-dimensional microscopy of the tumor microenvironment in vivo using optical frequency domain imaging,” Nat. Med. 15(10), 1219–1223 (2009). [CrossRef] [PubMed]

11. B. J. Vakoc, D. Fukumura, R. K. Jain, and B. E. Bouma, “Cancer imaging by optical coherence tomography: preclinical progress and clinical potential,” Nat. Rev. Cancer 12(5), 363–368 (2012). [CrossRef] [PubMed]

12. M. R. Munk, H. Giannakaki-Zimmermann, L. Berger, W. Huf, A. Ebneter, S. Wolf, and M. S. Zinkernagel, “OCT-angiography: A qualitative and quantitative comparison of 4 OCT-A devices,” PLoS One 12(5), e0177059 (2017). [CrossRef] [PubMed]

13. A. Mariampillai, B. A. Standish, E. H. Moriyama, M. Khurana, N. R. Munce, M. K. Leung, J. Jiang, A. Cable, B. C. Wilson, I. A. Vitkin, and V. X. Yang, “Speckle variance detection of microvasculature using swept-source optical coherence tomography,” Opt. Lett. 33(13), 1530–1532 (2008). [CrossRef] [PubMed]

14. Y. Jia, O. Tan, J. Tokayer, B. Potsaid, Y. Wang, J. J. Liu, M. F. Kraus, H. Subhash, J. G. Fujimoto, J. Hornegger, and D. Huang, “Split-spectrum amplitude-decorrelation angiography with optical coherence tomography,” Opt. Express 20(4), 4710–4725 (2012). [CrossRef] [PubMed]

15. C. Blatter, J. Weingast, A. Alex, B. Grajciar, W. Wieser, W. Drexler, R. Huber, and R. A. Leitgeb, “In situ structural and microangiographic assessment of human skin lesions with high-speed OCT,” Biomed. Opt. Express 3(10), 2636–2646 (2012). [CrossRef] [PubMed]

16. B. White, M. Pierce, N. Nassif, B. Cense, B. Park, G. Tearney, B. Bouma, T. Chen, and J. de Boer, “In vivo dynamic human retinal blood flow imaging using ultra-high-speed spectral domain optical coherence tomography,” Opt. Express 11(25), 3490–3497 (2003). [CrossRef] [PubMed]

17. D. Y. Kim, J. Fingler, J. S. Werner, D. M. Schwartz, S. E. Fraser, and R. J. Zawadzki, “In vivo volumetric imaging of human retinal circulation with phase-variance optical coherence tomography,” Biomed. Opt. Express 2(6), 1504–1513 (2011). [CrossRef] [PubMed]

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

19. 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]

20. S. Makita, K. Kurokawa, Y.-J. Hong, M. Miura, and Y. Yasuno, “Noise-immune complex correlation for optical coherence angiography based on standard and Jones matrix optical coherence tomography,” Biomed. Opt. Express 7(4), 1525–1548 (2016). [CrossRef] [PubMed]

21. J. Xi, L. Huo, J. Li, and X. Li, “Generic real-time uniform K-space sampling method for high-speed swept-Source optical coherence tomography,” Opt. Express 18(9), 9511–9517 (2010). [CrossRef] [PubMed]

22. H. C. Hendargo, R. P. McNabb, A.-H. Dhalla, N. Shepherd, and J. A. Izatt, “Doppler velocity detection limitations in spectrometer-based versus swept-source optical coherence tomography,” Biomed. Opt. Express 2(8), 2175–2188 (2011). [CrossRef] [PubMed]

23. B. Vakoc, S. Yun, J. de Boer, G. Tearney, and B. Bouma, “Phase-resolved optical frequency domain imaging,” Opt. Express 13(14), 5483–5493 (2005). [CrossRef] [PubMed]

24. B. Braaf, K. A. Vermeer, V. A. D. Sicam, E. van Zeeburg, J. C. van Meurs, and J. F. de Boer, “Phase-stabilized optical frequency domain imaging at 1-µm for the measurement of blood flow in the human choroid,” Opt. Express 19(21), 20886–20903 (2011). [CrossRef] [PubMed]

25. Y.-J. Hong, S. Makita, F. Jaillon, M. J. Ju, E. J. Min, B. H. Lee, M. Itoh, M. Miura, and Y. Yasuno, “High-penetration swept source Doppler optical coherence angiography by fully numerical phase stabilization,” Opt. Express 20(3), 2740–2760 (2012). [CrossRef] [PubMed]

26. V.-F. Duma, K. S. Lee, P. Meemon, and J. P. Rolland, “Experimental investigations of the scanning functions of galvanometer-based scanners with applications in OCT,” Appl. Opt. 50(29), 5735–5749 (2011). [CrossRef] [PubMed]

27. 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 µm,” Opt. Express 13(11), 3931–3944 (2005). [CrossRef] [PubMed]

28. S. O. Rice, “Mathematical analysis of random noise,” Bell Syst. Tech. J. 24(1), 46–156 (1945). [CrossRef]

29. M. Szkulmowski and M. Wojtkowski, “Averaging techniques for OCT imaging,” Opt. Express 21(8), 9757–9773 (2013). [CrossRef] [PubMed]

30. Y. Cheng, L. Guo, C. Pan, T. Lu, T. Hong, Z. Ding, and P. Li, “Statistical analysis of motion contrast in optical coherence tomography angiography,” J. Biomed. Opt. 20(11), 116004 (2015). [CrossRef] [PubMed]

31. Y. Lim, Y.-J. Hong, L. Duan, M. Yamanari, and Y. Yasuno, “Passive component based multifunctional Jones matrix swept source optical coherence tomography for Doppler and polarization imaging,” Opt. Lett. 37(11), 1958–1960 (2012). [CrossRef] [PubMed]

32. B. Baumann, W. Choi, B. Potsaid, D. Huang, J. S. Duker, and J. G. Fujimoto, “Swept source/Fourier domain polarization sensitive optical coherence tomography with a passive polarization delay unit,” Opt. Express 20(9), 10229–10241 (2012). [CrossRef] [PubMed]

33. 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]

34. A. Unterhuber, B. Považay, A. Müller, O. B. Jensen, M. Duelk, T. Le, P. M. Petersen, C. Velez, M. Esmaeelpour, P. E. Andersen, and W. Drexler, “Simultaneous dual wavelength eye-tracked ultrahigh resolution retinal and choroidal optical coherence tomography,” Opt. Lett. 38(21), 4312–4315 (2013). [CrossRef] [PubMed]

35. S. Yun, G. Tearney, J. de Boer, N. Iftimia, and B. Bouma, “High-speed optical frequency-domain imaging,” Opt. Express 11(22), 2953–2963 (2003). [CrossRef] [PubMed]

36. N. V. Iftimia, D. X. Hammer, C. E. Bigelow, D. I. Rosen, T. Ustun, A. A. Ferrante, D. Vu, and R. D. Ferguson, “Toward noninvasive measurement of blood hematocrit using spectral domain low coherence interferometry and retinal tracking,” Opt. Express 14(8), 3377–3388 (2006). [CrossRef] [PubMed]

37. R. Tripathi, N. Nassif, J. S. Nelson, B. H. Park, and J. F. de Boer, “Spectral shaping for non-Gaussian source spectra in optical coherence tomography,” Opt. Lett. 27(6), 406–408 (2002). [CrossRef] [PubMed]

38. I. J. White and H. D. Dederich, “American National Standard for Safe Use of Lasers, ANSI Z 136.1–2007,” Laser Inst. Am. Orlando (2007).

39. B. Braaf, K. V. Vienola, C. K. Sheehy, Q. Yang, K. A. Vermeer, P. Tiruveedhula, D. W. Arathorn, A. Roorda, and J. F. de Boer, “Real-time eye motion correction in phase-resolved OCT angiography with tracking SLO,” Biomed. Opt. Express 4(1), 51–65 (2013). [CrossRef] [PubMed]

40. T. Klein, W. Wieser, L. Reznicek, A. Neubauer, A. Kampik, and R. Huber, “Multi-MHz retinal OCT,” Biomed. Opt. Express 4(10), 1890–1908 (2013). [CrossRef] [PubMed]

41. R. Khazaeinezhad, M. Siddiqui, and B. J. Vakoc, “16 MHz wavelength-swept and wavelength-stepped laser architectures based on stretched-pulse active mode locking with a single continuously chirped fiber Bragg grating,” Opt. Lett. 42(10), 2046–2049 (2017). [CrossRef] [PubMed]

42. Y.-J. Hong, M. Miura, M. J. Ju, S. Makita, T. Iwasaki, and Y. Yasuno, “Simultaneous Investigation of Vascular and Retinal Pigment Epithelial Pathologies of Exudative Macular Diseases by Multifunctional Optical Coherence Tomography,” Invest. Ophthalmol. Vis. Sci. 55(8), 5016–5031 (2014). [CrossRef] [PubMed]

Complex differential variance angiography with noise-bias correction for optical coherence tomography of the retina

Abstract

1. Introduction