Logarithmic intensity and speckle-based motion contrast methods for human retinal vasculature visualization using swept source optical coherence tomography

Reza Motaghiannezam; Scott Fraser

doi:10.1364/BOE.3.000503

1. Introduction

The importance of retinal and choroidal vasculature visualization is paramount in diagnosing various eye diseases, such as diabetic retinopathy [1], age-related macular degeneration (AMD) [2], glaucoma [3], and anterior ischemic optic neuropathy (AION) [4].

Color fundus photography (CF) and fluorescein angiography (FA) have served as valuable imaging methods for retinal vasculature network visualization [5]. Deeper choroidal vessel imaging has also been realized by Indocyanine green angiography (ICGA) [5]. However, fluorescent dye injection can cause undesired side-effects [6,7], and the 2-D nature of these imaging techniques limits their applications for providing depth information and/or deep choroidal blood vessel visualization.

To meet the need for 3D retinal and choroidal vasculature assessment without the use of fluorescent dye injection, optical coherence tomography (OCT) has emerged as an attractive non-invasive depth-resolved imaging technology [8]. OCT realizes 3D visualization of retinal structure [9] and vasculature [10] by capturing interferograms of the light reflected back from the retina. Doppler OCT (D-OCT) [11] and Phase Contrast (PC)-OCT [12] are motion sensitive methods for vessel visualization, since these methods collect more than one depth scan of the same retinal loci and analyze Doppler phase shifts between depth scans at different depths. For example, a modified D-OCT method used autocorrelation and averaging techniques to capture retinal and choroidal vessels from complex OCT signals obtained in adjacent depth scans [13,14]. Although D-OCT captures the regions of high-velocity blood flow, such as in major vessels, it is unable to capture slow flow in retinal capillaries [11] or deep flows such as the choroidal circulation due to limited phase sensitivity and small time separation between two successive A-scans. Smart scanning protocols enhance sensitivity to the smaller signals expected from the microvasculature by increasing the time separation between two OCT depth scans and relying on the acquired phase [12,15–17] of OCT signals for contrast. However, these methods [15–17] such as PC-OCT [12] are highly sensitive to the phase instability of the system and environment by relying on the phase information. The required bulk motion removal [18] and timing-induced phase error correction algorithms [19] as well as extra optical module [19] add to the complexity of the phase sensitive swept source (SS)-OCT software and hardware.

To avoid complications of phase instability in retinal microvasculature visualization, several elegant intensity contrast imaging methods have relied on the temporal speckle fluctuations caused by mobile particles, such as red blood cells. The first method captured microvasculature through Hilbert and Fourier analyses of the spectral domain (SD)-OCT signal and a moving reference arm [20,21]. However, the system complexity and high sensitivity to hyper-reflective static layers may limit its clinical application for distinguishing these areas from regions of motion. Another approach used speckle signals [22,23], which have been studied extensively in the literature [24–27]. Although this speckle variance imaging method was able to capture vessels no smaller than 25 µm diameter size in the mouse dorsal skinfold [22,23], the multiplicative nature of the speckle made its variance highly sensitive to static regions similar to the first method [20,21], as shown in our results.

Thus, there is a need for a simple OCT method which does not rely on the phase information and provides highly motion-sensitive contrast for distinguishing regions of motion from stationary areas. The latter is especially important for detecting leakage and abnormal vessels in patients with abnormal retinal and choroidal structure.

Here, we purpose two motion-sensitive logarithmic intensity-based contrasts for 3D microvasculature visualization in an in vivo human retina. The proposed logarithm operation converts the multiplicative amplitude or intensity fluctuations (speckle) into the additive variations and recovers the motion contrasts by removing the speckle free signals (static regions) through statistical analysis. We develop two methods, called logarithmic intensity variance (LOGIV) and differential logarithmic intensity variance (DLOGIV), which provide these surrogate markers for motion by acquiring more than one depth scan of the same retinal voxel and analyzing logarithmic intensity fluctuations at the same voxel. To study and compare the proposed logarithmic contrasts with speckle amplitude and intensity contrasts, we formulate a theory for each. While theoretical and experimental results demonstrate the dependence of the speckle contrasts on the sample reflectivity strength and their dispensability for differentiating vessels from static regions, logarithmic contrasts show less sensitivity to the sample reflectivity strength. Application of LOGIV, DLOGIV, and speckle contrasts for 3D microvasculature imaging in the in vivo human retina is validated by employing a high-speed SS-OCT at 1060 nm. LOGIV and DLOGIV retinal en face views show the enhanced motion contrasts in comparison with speckle contrasts for capturing microvasculature that lies between hyper-reflective regions. Compared to the differential phase contrast method [12], these logarithmic intensity-based motion contrast methods are simpler, have similar performance, and do not require extra software and hardware [19].

2. Theory

The coherent summation of light waves with random path lengths causes destructive and constructive interferences called speckle [24]. The speckle pattern is stationary over time when the sample and optical imaging beam are static. However, backscattered light from moving particles, such as red blood cells, forms the speckle pattern that changes rapidly over time. Speckle has a multiplicative nature and its amplitude and intensity (OCT signal) follow Rayleigh and exponential distributions, respectively [25–27]. This highlights the potential application of speckle statistical analysis to capture blood vessels containing moving red blood cells [22,23]. However, the OCT signal sometimes contains specular components from strong (fixed) local scatterers, which drastically change the amplitude and intensity PDFs of speckle (OCT signal). In this scenario, amplitude and intensity statistics (such as mean and variance) depend not on only motion but also on hyper-reflective static regions due to the presence of these specular components, as shown in our formulated theory. Thus, amplitude and intensity-based (speckle) contrasts are not ideal for capturing and differentiating mobile scatterers from static regions. To remove the direct dependence of the speckle on the sample reflectivity, we test two different approaches: (1) mean scaling of the speckle contrasts (speckle contrast ratios) and (2) logarithmic intensity along with variance analysis. In order to study and compare these intensity-based contrasts, we develop a theory. Our theoretical results demonstrate high sensitivity of speckle contrast ratios to hyper-reflective static regions. However, the natural logarithm function and variance analysis separate the static signals (specular components) from speckle and cancel specular components (fixed terms), respectively. The developed theory shows that logarithmic intensity contrasts outperform both speckle contrasts and speckle contrast ratios in terms of motion sensitivity.

Here, we first derive the generalized PDFs of the amplitude and intensity speckle in Fourier domain OCT. The Fourier transform of the received interferograms from a Fourier domain OCT system is intrinsically a complex stationary random process, and represents the amplitude (A) and phase (ϕ) at a certain location and time as well as a zero mean complex additive white Gaussian wide-sense stationary random process (N(0,σ_n²)).

U (z, t) = U_{Re a l} (z, t) + j U_{Im a g} (z, t) = A (z, t) e^{j φ (z, t)} + N (z, t)

where U_Real and U_Imag are the real and imaginary parts of the complex OCT signals, respectively.

At a given depth (z) in the sample, the first term in Eq. (1) is proportional to a coherent sum of all backscattered light from N scatterers in a coherence volume with a size given by the probe beam area and coherence length of the light source.

U_{Re a l} (z, t) = A_{r} (t) \sum_{i = 1}^{N} A_{s} (x_{i}, y_{i}, z, t) \cos (φ (x_{i}, y_{i}, z, t)) + N_{Re a l} (z, t) = V_{1} (z, t) + N_{Re a l} (z, t)

U_{Im a g} (z, t) = A_{r} (t) \sum_{i = 1}^{N} A_{s} (x_{i}, y_{i}, z, t) \sin (φ (x_{i}, y_{i}, z, t)) + N_{Im a g} (z, t) = V_{2} (z, t) + N_{Im a g} (z, t)

A_r, A_s, N_real(0,σ_nr² = σ_n²/2), and N_imag(0,σ_ni² = σ_n²/2) are, respectively, the amplitude of light returned from the reference arm, the amplitude of backscattered light from a given scatterer in the sample arm located at (x_i, y_i, z), and zero mean real additive white Gaussian wide-sense stationary random processes. σ_nr² and σ_ni² are variances of the real and imaginary parts of a zero mean complex additive white Gaussian wide-sense stationary random process (N). In Eqs. (2) and (3), V₁(µ_v1,σ_v1² = σ_v²/2) and V₂(µ_v2, σ_v2² = σ_v²/2) are the sum of many independent and identical random processes which approach two independent real Gaussian random processes with equal variances of σ_v1² and σ_v2², respectively. Thus, U(z,t) can be expressed as summation of a constant complex value (µ_v1(z) + jµ_v2(z)) and a zero mean complex Gaussian wide-sense stationary random process D(0, σ² = σ_n² + σ_v²).

U (z, t) = (μ_{v 1} (z) + j μ_{v 2} (z)) + D (z, t)

Linear intensity (I) is defined as

I (z, t) = {| U (z, t) |}^{2} = U_{Re a l}^{2} (z, t) + U_{Im a g}^{2} (z, t)

Thus, the intensity (I) and amplitude (M = √I) of OCT signal have a non-central chi-square distribution with 2 degrees of freedom and a Ricean distribution [28], respectively. PDFs of I and M are given by

p_{I} (i) = \frac{1}{2 (σ_{n r}^{2} + σ_{v 1}^{2})} e^{\frac{- (s_{I}^{2} + i)}{2 (σ_{n r}^{2} + σ_{v 1}^{2})}} I_{0} (\frac{s_{I} \sqrt{i}}{σ_{n r}^{2} + σ_{v 1}^{2}}) = \frac{1}{σ_{}^{2}} e^{\frac{- (s_{I}^{2} + i)}{σ_{}^{2}}} I_{0} (2 \frac{s_{I} \sqrt{i}}{σ_{}^{2}})

p_{M} (m) = \frac{m}{(σ_{n r}^{2} + σ_{v 1}^{2})} e^{\frac{- (s_{I}^{2} + m^{2})}{2 (σ_{n r}^{2} + σ_{v 1}^{2})}} I_{0} (\frac{s_{I} m}{σ_{n r}^{2} + σ_{v 1}^{2}}) = \frac{2 m}{σ_{}^{2}} e^{\frac{- (s_{I}^{2} + m^{2})}{σ_{}^{2}}} I_{0} (2 \frac{s_{I} m}{σ_{}^{2}})

where the non-centrality parameter of distributions and σ² are

s_{I}^{2} = μ_{v 1}^{2} + μ_{v 2}^{2}

σ^{2} = σ_{n r}^{2} + σ_{n i}^{2} + σ_{v 1}^{2} + σ_{v 2}^{2} = 2 (σ_{n r}^{2} + σ_{v 1}^{2}) = σ_{n}^{2} + σ_{v}^{2}

and I₀(.), the zeroth-order modified Bessel function of the first kind, is

I_{0} (c) = \sum_{k = 0}^{\infty} \frac{1}{k! Γ (k + 1)} {(\frac{c}{2})}^{2 k}

The non-centrality parameter shows the presence of specular components in OCT signal. Equations (6a) and (6b) approach exponential and Rayleigh distributions when the non-centrality parameter is zero (no specular component). While the thermal noise, shot noise, and relative intensity noise (RIN) of the laser source mainly determine σ_n², laser intensity fluctuation, transversal displacements between the scanning optical beam and sample, and mobile scatterers are contributing factors in σ_v²_.

To find speckle amplitude variance (SAV), speckle intensity variance (SIV), speckle amplitude contrast ratio (SACR), speckle intensity contrast ratio (SICR), LOGIV, and DLOGIV as a function of OCT signal-to-noise ratio (SNR = s_I²/σ²), the mean and variance of the amplitude (M), intensity (I), and natural logarithm of intensity (log(I)) are calculated.

2.1 Speckle amplitude and intensity variances

The mean and variance of the chi-square-distributed random variable with 2 degrees of freedom (I) and the Ricean random variable (M) are [28]

μ_{I} = 2 (σ_{n r}^{2} + σ_{v 1}^{2}) + s_{I}^{2} = σ_{}^{2} + s_{I}^{2} = σ_{}^{2} (1 + S N R)

S I V = σ_{I}^{2} = 4 {(σ_{n r}^{2} + σ_{v 1}^{2})}^{2} + 4 (σ_{n r}^{2} + σ_{v 1}^{2}) s_{I}^{2} = σ_{}^{4} + 2 σ_{}^{2} s_{I}^{2} = σ_{}^{4} (1 + 2 S N R)

μ_{M} = \sqrt{\frac{π}{2} (σ_{n r}^{2} + σ_{v 1}^{2})} L_{0.5} (\frac{- s_{I}^{2}}{2 (σ_{n r}^{2} + σ_{v 1}^{2})}) = \frac{σ \sqrt{π}}{2} L_{0.5} (\frac{- s_{I}^{2}}{σ^{2}}) = \frac{σ \sqrt{π}}{2} L_{0.5} (- S N R)

S A V = σ_{M}^{2} = 2 (σ_{n r}^{2} + σ_{v 1}^{2}) + s_{I}^{2} - \frac{π}{2} (σ_{n r}^{2} + σ_{v 1}^{2}) {[L_{0.5} (\frac{- s_{I}^{2}}{2 (σ_{n r}^{2} + σ_{v 1}^{2})})]}^{2} = σ_{}^{2} {1 + S N R - \frac{π}{4} {[L_{0.5} (- S N R)]}^{2}}

where L_0.5 denotes a Laguerre polynomial.

Equations (9b) and (10b) clearly show that SIV (σ²_I) and SAV (σ²_M) depend not on only the noise power (σ_n²) and motion (σ_v²) but also on the sample reflectivity strength (s_I²). Both SIV and SAV are also monotonically increasing functions of SNR. Therefore, speckle variance contrasts measure the noise power and motion of scatterers amplified with SNR value at the associated region. This highlights the inability of SIV and SAV to differentiate motion from static areas (σ_v≈0), especially hyper-reflective static regions with high SNR values.

2.2 Speckle amplitude and intensity contrast ratios

The speckle intensity and amplitude contrast ratios (SICR and SACR) are the ratios of the standard deviation to the mean of intensity and amplitude, respectively [29]. Using the derived Eqs. (9a) and (9b), and Eqs. (10a) and (10b), SICR and SACR can be expressed simply as following monotonically decreasing functions of SNR:

S I C R = \frac{σ_{I}}{μ_{I}} = \frac{\sqrt{1 + 2 S N R}}{1 + S N R}

S A C R = \frac{σ_{M}}{μ_{M}} = \sqrt{\frac{4 (1 + S N R)}{π {[L_{0.5} (- S N R)]}^{2}} - 1}

Although these ratios cancel the direct dependence of speckle contrasts to the sample reflectivity, the derived Eqs. (11a) and (11b) demonstrate the SNR dependence of speckle contrast ratios and their mechanisms for capturing motion. Lower SNR is associated with higher speckle contrast ratios and motion. Monotonically decreasing nature of SICR and SACR makes these two contrasts less sensitive to hyper-reflective regions. However, their narrow dynamic ranges (1 and √4π⁻¹-1) and nonzero values of SICR and SACR at hyper-reflective static regions limit their application for capturing motion and differentiating it from structural signal, as described later in details.

2.3 Logarithmic intensity and differential logarithmic intensity variances

As described in appendix, the mean of the defined logarithm of a non-central chi-square random variable is

m_{\log (I)} = {\begin{matrix} \log (s_{I}^{2}) + Γ (0, \frac{s_{I}^{2}}{σ_{}^{2}}) = \log (s_{I}^{2}) + \int_{\frac{s_{I}^{2}}{σ_{n}^{2}}}^{+ \infty} \frac{e^{- t}}{t} d t \\ \log (σ_{}^{2}) - γ \end{matrix}

Therefore, the mean of the natural logarithm of the normalized intensity to the non-centrality parameter (s_I²) is equal to Γ(0, s_I²/σ²).

The variance of log(I), LOGIV, can also be given by

L O G I V = σ_{\log (I)}^{2} = \sum_{n = 1}^{\infty} \frac{Q (2 n, S N R)}{n^{2}}

where Q is the normalized incomplete Gamma function [30]

Q (a, b) = \frac{Γ (a, b)}{Γ (a)}

It follows that the variance of differences of logarithm intensities (two identical and independent random processes with equal variances) is

D L O G I V = 2 σ_{\log (I)}^{2} = 2 \sum_{n = 1}^{\infty} \frac{Q (2 n, S N R)}{n^{2}}

LOGIV and DLOGIV approach π²/6 and π²/3, respectively, when SNR (s_I²/σ²) goes to zero (Eqs. (46) and (15)).

To validate our theoretically derived Eqs. (12) and (13) for the mean and variance of the natural logarithm of intensity, one hundred thousand zero mean complex Gaussian random variables with a defined variance (σ²) were generated and added to the desired complex signal with different intensity values (µ_v1² + µ_v2²), as expressed in Eq. (4). The red solid and blue dotted curves in Fig. 1(a) depict the derived Eq. (12) and the simulated result for the mean of the logarithm of the normalized intensity to the non-centrality parameter, Γ(0,(s_I/σ)²), as a function of different SNR values, respectively. The derived Eq. (12) clearly follows the simulated results. The mean of the logarithm of intensity approaches the logarithm of the non-centrality parameter for SNR>1. Figure 1(b) shows the derived Eq. (13) (blue dotted curve) and simulated result (red sold curve) for the variance of the logarithm of intensity. A close fit is found between the theoretical and simulated results in Fig. 1(b). The variance of the logarithm of intensity is π²/6 at a SNR equal to zero. When SNR is increased, the variance of the logarithm of intensity approaches zero. Thus, the motion sensing is based on the intensity fluctuation in logarithmic contrasts which are monotonically decreasing functions of SNR. A higher incomplete gamma function value detects lower SNR caused by the motion.

Fig. 1 The theoretical (blue dotted curves) and simulated (red solid curves) (a) mean of logarithm of normalized intensity to the non-centrality parameter and (b) variance of logarithm of intensity as a function of SNR. (c) The theoretical SACR, SICR, LOGIV, and DLOGIV as function of SNR. The close-up view shows the normalized SACR, SICR, LOGIV, and DLOGIV as function of SNR in Fig. 1(c).

Download Full Size | PDF

2.4 Comparison of speckle- and logarithmic intensity-based contrasts

To capture motion caused by the intensity fluctuation of OCT signal, the ideal intensity-based motion contrasts need to be: (1) independent of the sample reflectivity strength (s_I²) and (2) sensitive to motion over a wide SNR range. Theoretical results (Eqs. (9b) and (10b)) clearly show that speckle variances, SIV and SAV, depend on the sample reflectivity strength and are incapable of differentiating the motion signals from structural signals.

Although mean scaling and logarithm operation along with variance analysis are able to remove the direct dependence of the speckle on the sample reflectivity (Eqs. (11a), (11b), (13), and (15)), speckle contrast ratios (SICR and SACR) and logarithmic-based contrasts (LOGIV and DLOGIV) vary with SNR (11a-b, 13, and 15). To compare SICR, SACR, LOGIV, and DLOGIV variations as a function of SNR, the derived Eqs. (11a), (11b), (13), and (15) were plotted in Fig. 1(c). The tail of normalized graph to the dynamic range describes the contrast sensitivity to hyper reflective-regions (Fig. 1(c), close-up view). While all contrasts decrease and approach zero with increasing SNR (Fig. 1(c)), logarithmic intensity-based contrasts demonstrate faster decline. Logarithmic intensity-based contrasts also show wider dynamic ranges than speckle contrast ratios (Fig. 1(c)). Comparisons of four normalized graph tails (Fig. 1(c), close-up view) demonstrate that LOGIV and DLOGIV are less sensitive to hyper-reflective (high SNR value) regions.

Thus, less sensitivity to hyper-reflective regions make the logarithmic intensity-based motion contrasts more attractive than speckle contrasts for capturing mobile scatterers and differentiating motion from static regions in OCT applications such as ophthalmic OCT with SNR values < 35 dB.

3. Experimental methods

3.1 SS-OCT System

To validate the proposed intensity-based methods for providing motion contrasts and compare them with each other and DPV method [12,17], we used a prototype 50 kHz phase sensitive SS-OCT system, incorporating a polygon-based 1060 nm (1015-1103) swept laser source [31], with ~5.9 μm axial resolution in tissue and 102 dB sensitivity (1.2 mW incident power). The SS-OCT system was comprised of the polygon-based swept-laser source [31], an interferometer, and a data acquisition (DAQ) unit (Fig. 2 ). The swept source output was coupled to the interferometer through an isolator where a 90/10 coupler was used to split light into a sample arm: reference arm. The sample arm light was split equally between the calibration arm and a slit lamp biomicroscope (Carl Zeiss Meditec) as shown in Fig. 2. A 50/50 coupler combined and directed the reflected light from the sample to the one port of the interferometer output coupler. The reference arm light passed through a pair of collimators and was directed to the second port of the interferometer output coupler. The resulting interference fringes were detected on both output ports using a dual balanced photodetector. To generate a trigger signal at the beginning of the first interference fringe for data acquisition, a portion of the reference arm light was directed to a three port circulator and a fiber Bragg grating. A reflected optical pulse was detected using a photodetector and converted into an electronically tunable TTL signal as a trigger signal. The spectral signals were continuously digitized by triggering an AD conversion board. A D/A board was used to generate the driving signals of the two-axis galvanometers. A user interface and data acquisition was developed in LabView to coordinate instrument control and enable user interaction.

Fig. 2 Schematic of a polygon-based swept-laser source (gray box), interferometer (yellow box), and SS-OCT data acquisition unit used for imaging. SOA, semiconductor optical amplifier; diffraction grating (830 lines/mm); telescope (f₁ = 40 mm, f₂ = 35 mm); polygonal scanning mirrors (72 facets); COL, collimator; PC, polarization controller; ISO, isolator; CIRC, circulator; CUP, coupler; PD, photodetector; BD, dual balanced photodetector (InGaAs, 80 MHz); FBG, fiber Bragg grating (0.1 nm). Data acquisition unit comprising of an AD conversion board (14-bit, GaGe CompuScope 14,200) for digitizing interference fringes and a D/A board (16 bit, National Instruments) for driving scanning mirrors (Cambridge Technology).

Download Full Size | PDF

3.2 Scanning protocols

The prototype SS-OCT instrument was used to image the eye of a healthy volunteer (37-year-old Caucasian male) at the California Institute of Technology (Caltech). This study was approved by the Caltech committee for the protection of human subjects and adhered to the tenets of the Declaration of Helsinki. Written informed consent was obtained from subjects prior to OCT imaging. Total exposure time and incident exposure level were kept less than 5.5 seconds and 1.2 mW in each imaging session, consistent with the safe exposure determined by American National Standards Institute (ANSI) and International Commission on Non-Ionizing Radiation Protection (ICNIRP). Patient interfaces were based on a Stratus OCT-3 system (Carl Zeiss Meditec) adapted with optics to support the 1060 nm wavelength range. A 60-D lens was used at the exit of the fundus camera with 13 mm working distance providing a beam diameter of 1.5 mm on the cornea (~15 μm transverse resolution).

Three scanning protocols were implemented. A 1D protocol generated 1900 depth scans (A-scans) over the same transverse position in 0.038 seconds. A 2D protocol acquired four horizontal tomograms (B-scans) with 201 depth scans (A-scans) spanning the same transverse slice (2 mm) across the foveal centralis in 0.02 seconds. In the third protocol, a 3D OCT data set was collected by acquiring several neighboring B-scans over the parafovea. The system magnification, SS-OCT speed (50400 Hz), speed of the fast scan axis (200 Hz) with fly-back time (1 ms), and data acquisition time (4 seconds) gave an image size of 201 × 200 pixels over a 2 mm × 2 mm field of view (FOV); each B-scan was repeated four times. In the 3D scanning protocol, the fast scan axis was sagittal (superior-inferior) and the slow axis was horizontal (nasal-temporal).

3.3 Image processing and motion contrast imaging

The digitized signals were divided into individual spectral sweeps in the post-processing algorithm. Equal sample spacing in wave number (k) was achieved using a calibration trace at 1.5 mm interferometer delay and numerical correction of the nonlinearly swept waveforms [32]. Image background subtraction and numeric compensation for second order dispersion [33] were performed. The complex SS-OCT data sets were upsampled by a factor of 4 and Fourier transformed. Axial motion correction was achieved on the obtained 2D and 3D SS-OCT data sets by cross correlating the consecutive horizontal tomograms. The estimated mean and variance of the linear intensity, and logarithm of intensity were calculated for all voxels through acquired depth scans.

To perform motion contrast analysis and imaging, multiple A (or B)-scans were acquired over the same transverse position (or slice). Time separations were T_A = 20 µs and T_B = 5 ms between A- and B-scans for capturing the same position, respectively. Multiple linear intensity and phase measurements were recorded over the same transverse point separated in time. Six different intensity-based approaches were tested: SIV, SAV, SICR, SACR, LOGIV, and DLOGIV. In this paper, the number of measured samples (the linear intensity, logarithm of intensity, and phase) to estimate each of the different motion contrasts was 1900 in the 1D protocol for each voxel through depth. In the 2D and 3D protocols, four samples were measured and used to estimate these motion contrasts.

In the SICR, SACR, SIV, and SAV methods, the estimated amplitude and linear intensity means (µ), variances (σ²) as well as the ratios between their estimated standard deviations and means (σ/µ) were calculated for the same transverse point acquired in successive A (or B)-scans (Figs. 3(b) –3(e)). LOGIV was realized by calculating the estimated variance of multiple logarithmic intensity measurements (LOG(I(z,T))) of the same transverse point acquired in successive A (or B)-scans separated in time (Fig. 3(f)). DLOGIV and DPV captured the differences between multiple logarithmic intensity (LOG(I(z,T))) and phase measurements (φ(z,T)) of the same transverse points (separated in time) and calculated the estimated variance of these changes (Figs. 3(a) and 3(g)), respectively. To measure and remove timing-induced phase error [19] due to the random delay between the trigger signal and the subsequent A-to-D conversion (sample clock), a calibration signal was generated using a stationary mirror in the calibration arm (Fig. 3) as described in [19]. The calibration signal was located at a depth of 2 mm in the OCT intensity image. The corrected phase differences between adjacent A (or B)-scans for the same transverse point at a given depth were calculated by subtracting the phase difference of the calibration signal, linearly scaled with the sample signal depth, from the measured phase differences (Fig. 3(a)) [19]. Phase unwrapping was performed on all measurements. A weighted mean algorithm [18] estimated and removed the bulk axial motion phase change error (Fig. 3(a)).

Fig. 3 Flowchart representing the data processing procedures for generating different motion contrasts of OCT: (a) DPV, (b) SICR, (c) SIV, (d) SACR, (e) SAV, (f) LOGIV, and (g) DLOGIV.

Download Full Size | PDF

The same described procedures were repeated for the adjacent transverse points in the same and neighboring B-scans to capture the retinal vasculature in 2D and 3D data sets. To remove SNR-limited intensity and phase change errors in 2D and 3D data sets for vasculature visualization, an average intensity threshold (10 dB above the mean value of the noise floor) was applied; all contrasts with average intensity values < mean (log₁₀(I_noise)) + 10 dB were set to zero in the corresponding images (Fig. 3). Additional processing was applied to the 3D contrast data sets using a moving average filter of ~5 μm in the axial direction and a median filter of ~8 μm radius in the transverse directions.

To create the retinal en face views, the inner/outer photoreceptor segments (IS/OS) and vitreoretinal interface were detected using a segmentation algorithm [10]. Several depth integrated motion contrast en face images were generated by integrating the SICR, SACR, SIV, SAV, LOGIV, DLOGIV, and DPV between three different regions in the inner retina relative to IS/OS and vitreoretinal interface.

4. Experimental results

4.1 Experimental analysis of the linear intensity through retinal and choroidal layers

In order to statistically analyze the linear intensity of the retina and choroid, a stationary optical beam was used to capture 1900 depth scans of a given transverse location through an in vivo human retina and choroid. The 2D foveal centralis tomogram (Fig. 4(a) ) shows the transverse location and depths (red line) through which these A-scans were acquired in ~38 ms. Figure 4(b) presents the retinal and choroidal depth profile as a function of time (A-scan number) without alignment for the defined transverse position. This figure depicts the negligible axial bulk motion during collection of this data set. Dash-dotted blue and solid red lines show the estimated mean and standard deviation of the linear intensity through the retina and choroid in Fig. 4(c), respectively. This figure clearly demonstrates SIV dependence on the layer reflectivity as described by Eq. (9b). The similar result can be found for the SAV.

Fig. 4 (a) Averaged intensity tomogram across the fovea centralis (5 mm). Red line depicts the targeted transverse location and depths by a stationary optical beam. (b) The retinal and choroidal structure was captured in 1900 repeated A-scans through depth (red line in Fig. 4(a)). (c) The estimated mean (dash-dotted blue line) and standard deviation (solid red line) of the linear intensity through the retinal and choroidal layers. Histograms and theoretically derived PDFs (black dashed lines) of the captured linear intensity at the level of (d) RNFL, (e) IS, (f) RPE, (g) CC, (h) SL, and (i) HL. Histograms and theoretically derived PDFs (black dashed lines) of the captured amplitude at the level of (j) RNFL, (k) IS, (l) RPE, (m) CC, (n) SL, and (o) HL.

Download Full Size | PDF

To study PDFs of the linear intensity and amplitude of OCT signal at different retinal and choroidal layers, the estimated non-centrality parameter (s_I²) and standard deviation of the measured linear intensity and amplitude (σ_I and σ_M) were calculated throughout the captured depth profile. Figures 4(d)–4(o) show histograms of the measured linear intensity and amplitude at the level of six different retinal and choroidal layers, including the retina nerve fiber layer (RNFL), IS (inner segment), retina pigment epithelium (RPE), choriocapillaris (CC), Sattler’s layer (SL), and Haller’s layer (HL). The theoretically derived PDFs of the linear intensity (Eq. (6a)) are depicted as black dashed lines in Figs. 4(d)–4(i). These histograms follow chi-square distributions with 2 degrees of freedom (Eq. (6a)) for the estimated non-centrality parameter and standard deviation of the measured linear intensity. As shown in Eq. (6a), the PDF shape varies with the non-centrality parameter. For small non-centrality parameter values, the zeroth-order modified Bessel function of the first kind in Eq. (6a) is approximately equal to one.

In this scenario, the PDF of the linear intensity (Eq. (6a) approaches an exponential, as described previously [25–27]. Figures 4(d)–4(i) demonstrate a good agreement between experimental and theoretical results. Hyper- (Figs. 4(d)–4(f)) and hypo-reflective (Figs. 4(g)–4(i)) layers can be distinguished from each other using their linear intensity PDFs. Dashed lines in Figs. 4(j)–4(o) present the theoretically derived PDFs of the amplitude (Eq. (6b)). A close fit was found between histograms (Figs. 4(j)–4(o)) and a Ricean distribution (Eq. (6b)). The PDF of the amplitude (Eq. (6b)) matches with a Rayleigh distribution [25–27] when the non-centrality parameter approaches small values.

4.2 Quantitative analysis of speckle contrast ratios and logarithmic intensity-based contrasts though the retinal and choroidal layers

To investigate the quantitative variation of SICR, SACR, LOGIV, and DLOGIV through an in vivo human retina and choroid, a SS-OCT patient interface provided an optical beam for illuminating ~a 15 µm spot region on the retina of a healthy subject. The backscattered light was combined with the reference arm optical beam and detected in ~38 ms. 1900 depth scans of a given transverse location through the retina and choroid were captured. SICR, SACR, LOGIV, and DLOGIV were calculated throughout 750 µm depth by measuring the estimated mean (µ) and standard deviation (σ) of the linear intensity and amplitude as well as the estimated variance of the logarithm of intensity and differences of the logarithm of intensities over the same data set as previously described in the motion contrast imaging section (3.3).

Figures 5(a) –5(d) depict quantitative variations of SICR, SACR, LOGIV, and DLOGIV through the retina and choroid, respectively. The captured dynamic ranges of SICR, SACR, LOGIV, and DLOGIV are ~1, √(4π⁻¹-1), π²/6, and π²/3 as predicted by the theoretical and simulation results (Figs. 1(b) and 1(c)). Although LOGIV and DLOGIV values of hyper-reflective stationary layers such as RNFL, IS, and RPE are approximately zero (Figs. 5(c) and 5(d)), the measured SICR and SACR values are nonzero for the corresponding regions (Figs. 5(a) and 5(b)). Nonzero values in Figs. 5(c) and 5(d) are associated with low SNR depths and/or regions of motion. These quantitative measurements (Figs. 5(a)–5(d)) are consistent with the theoretical results (Fig. 1(c)). LOGIV and DLOGIV demonstrate their potential application for capturing motion due to their wide dynamic range and low sensitivity to the sample reflectivity strength. Figure 5(c) shows that LOGIV values are ~π²/6 in the choroid and sclera (>680 µm). However, only DLOGIV exhibits different values over a wide range for different choroidal regions and sclera (Fig. 5(d)). This highlights DLOGIV capability for differentiating CC and SL from the outer choroid and capturing the choroid-sclera junction.

Fig. 5 Quantitative variations of (a) SICR, (b) SACR, (c) LOGIV, and (d) DLOGIV through the retina and choroid.

Download Full Size | PDF

4.3 2D tomogram and en face view visualization of the retina using motion contrast imaging methods

To study different motion contrast methods, four B-scans were acquired across the foveal centralis (2 mm). The averaged intensity of four obtained B-scans is depicted in Fig. 6(a) . 2D SICR and SACR tomograms (Figs. 6(b) and 6(c)) show that these speckle contrast ratios capture not only regions of motion (between blue lines) in the inner choroid and small vessels (white arrows) in the inner retina but also highly reflective stationary regions in IS/OS, RPE (between red lines), and the inner retina. This finding is consistent with our theoretical result (Fig. 1(c)). While the SIV (Fig. 6(d)) is able to capture the inner retina vessels (white arrow), it highlights the static regions of IS/OS and RPE (between red lines) as motion. Motion in the inner choroid is barely detected in this tomogram. Figures 6(e) and 6(f) show the enhanced motion contrast in 2D LOGIV and DLOGIV tomograms. White static areas (between red lines) captured in 2D speckle tomograms (Figs. 6(b)–7(d) ) are invisible in 2D LOGIV and DLOGIV tomograms (Figs. 6(e) and 6(f)). Regions of motion in the inner choroid (white band between blue lines) and the small vessels in the inner retina (white arrows) are detectable in these 2D tomograms (Figs. 6(e) and 6(f)).

Fig. 6 Foveal (a) average intensity, (b) SICR, (c) SACR, (d) SAV, (e) LOGIV, (f) DLOGIV, (g) DPV before phase compensation, and (h) DPV after phase compensation tomograms (2 mm). White regions correspond to regions with higher either motion or/and reflectivity. White arrows indicate the small vessel in Figs. 6(b)–6(h). IS/OS and RPE are signified between red lines (static regions). Regions of motion in the inner choroid are indicated between blue lines.

Download Full Size | PDF

Fig. 7 Parafoveal depth-integrated en face views over 4 mm² FOV acquired in 4 seconds. Inverted (a) averaged intensity, (b) SICR, (c) SIV, (d) SACR, (e) SAV, (f) LOGIV, (g) DLOGIV, and (h) DPV en face images of the inner retina.

Download Full Size | PDF

To compare the intensity-based contrasts with DPV contrast, 2D DPV tomograms are shown in Figs. 6(g) and 6(h) before and after compensation, respectively. Figures 6(g) demonstrate DPV is unable to capture motion without use of compensation algorithms [18,19] and an extra hardware module [19]. In addition, the calibration mirror image limits imaging depth [19]. Thus, the simplicity and motion sensitivity of LOGIV and DLOGIV may make these two contrast methods more attractive than other proposed phase- and intensity-based methods [12,15–17,20–23] for capturing motion and microvasculature.

Figures 7(a)–7(h) illustrate the inverted intensity, SICR, SIV, SACR, SAV, LOGIV, DLOGIV, and DPV en face views generated by integrating their values between the region 30 μm posterior to the vitreoretinal interface and the region 130 μm anterior to IS/OS. Figure 7(a) shows that the meshwork of capillaries is barely visible in the intensity en face view. Although small vessels and capillaries are seen in the SICR, SIV, SACR en face images (Figs. 7(b)–7(d)), the narrow dynamic range and high sensitivity to hyper-reflective static regions degrade retinal microvasculature enface visualization through contrast integration in the depth. Gray areas highlight the hyper-reflective stationary regions captured around the fovea avascular zone (FAZ) and between the interconnected microvasculature networks (Figs. 7(b)–7(d)). While SAV captures not only motion but also static regions (Fig. 6(d), between red lines), it is able to visualize the meshwork of capillaries through contrast summation in depth (Fig. 7(e)). Motion contrast enhancement is depicted in Figs. 7(f) and 7(g) using LOGIV and DLOGIV methods. Blood vessels in the ganglion cell layer and capillary meshwork of the inner plexiform layer are visualized in the LOGIV and DLOGIV en face views (Figs. 7(f) and 7(g)). FAZ is resolvable by considering the capillary network around it as shown in the LOGIV and DLOGIV images in Figs. 7(f) and 7(g). To compare retinal visualization using the proposed intensity-based motion contrast methods with the phase contrast method, the DPV en face image (Fig. 7(h)) is generated by summing DPVs over the same regions in the inner retina. Although LOGIV, DLOGIV, and DPV en face images (Figs. 7(f)–7(h)) achieve the similar contrast for foveal vasculature visualization, DPV is a complicated method due to its need for the compensation algorithms [18,19] and an extra optical module [19].

To show the capillary meshwork of the inner retina through depth using logarithmic intensity-based motion contrast methods, the LOGIV and DLOGIV en face views are generated by integrating their values through different depths. Figures 8(a) and 8(b) show the capillary network of the inner retina between the regions 255 μm and 216 μm anterior to IS/OS in the inverted LOGIV, and DLOGIV en face views. The inverted DPV en face view (Fig. 8(c)) depicts the similar capillary meshwork of the inner retina in the same region. Similar retinal microvasculature network is also detected between the regions 216 μm and 169 μm anterior to IS/OS (Figs. 8(d)–8(f)) in the inverted LOGIV, DLOGIV, and DPV en face views. Figures 8(a)–8(f) clearly reveal depth-related variations of capillary meshwork morphology through the inner retina.

Fig. 8 Parafoveal depth-integrated en face views over 4 mm² FOV acquired in 4 seconds. Inverted (a) LOGIV, (b) DLOGIV, and (c) DPV en face images of the retina between the regions 255 μm and 216 μm anterior to IS/OS. Inverted (d) LOGIV, (e) DLOGIV, and (f) DPV en face images of the retina between the regions 216 μm and 169 μm anterior to IS/OS.

Download Full Size | PDF

5. Discussion

Our formulated theoretical description of the Fourier-domain OCT signal shows that the high sensitivity of amplitude and intensity speckle-based contrasts to static regions degrades microvasculature visualization. Inability of these contrasts to distinguish the motion signals from structural signals would also limit their application in detecting leakage and abnormal vessels in patients with irregular retinal and choroidal vessels and structure. To improve upon these methods, we described logarithmic intensity-based motion contrast methods and demonstrated their superiority for in vivo human parafoveal microvasculature visualization. Our results showed the potential application of depth-resolved LOGIV and DLOGIV methods for replacing invasive FA and ICGA for diagnosing eye diseases in the future. Microvasculature visualization in the retina with these logarithmic intensity-based motion contrasts also demonstrated their superior capability over previously reported phase [12,15–17] and intensity-based motion OCT [20–23] techniques without the need for compensation algorithms and extra optical modules. Further improvements to the logarithmic intensity-based motion contrast imaging methods might enhance vascular visualization, reduce the effects of inherent laser source intensity variation, and increase the scanning angle of the retina.

6. Appendix

To calculate the mean and variance of the natural logarithm of the intensity (log (I)), we define a new chi-square random variable W (=I/σ²) with a non-centrality parameter of s_W=s_I/σ. Thus, log(I)can be written as

\log (I) = \log (σ_{}^{2} W) = \log (σ_{}^{2}) + \log (W)

where the PDF of W is given by

p_{W} (w) = e^{- (s_{W}^{2} + w)} I_{0} (2 s_{W} \sqrt{w})

The mean of the defined logarithm of a non-central chi-square random variable is

m_{\log (I)} = \log (σ_{}^{2}) + \int_{0}^{+ \infty} \log (w) e^{- (s_{W}^{2} + w)} I_{0} (2 s_{W} \sqrt{w}) d w

To simplify Eq. (18), we use Eq. (8) and the n^th derivative of the Gamma function given by [30]

Γ^{(n)} (k) = \frac{d^{n}}{d k^{n}} Γ (k) = \int_{0}^{+ \infty} w^{k - 1} e^{- w} {[\log (w)]}^{n} d w

Therefore

m_{\log (I)} = \log (σ_{}^{2}) + e^{- s_{w}^{2}} \sum_{k = 0}^{\infty} \frac{{(s_{W}^{2})}^{k}}{k! Γ (k + 1)} Γ^{(1)} (k + 1)

The first derivative of Gamma function can be expressed as [30]

\frac{Γ^{(1)} (k + 1)}{Γ (k + 1)} = ψ_{0} (k + 1) = - γ + H_{k}^{(1)} = - γ + \sum_{j = 1}^{k} \frac{1}{j}

where ψ₀, γ (=0.577), and H_k (or H_k⁽¹⁾) are the digamma function, Euler's constant, and harmonic number, respectively.

By using the digamma function (21), the mean of log(I) in Eq. (20) can be rewritten as

m_{\log (I)} = \log (σ_{}^{2}) - γ + e^{- s_{w}^{2}} \sum_{k = 0}^{\infty} \frac{{(s_{W}^{2})}^{k} H_{k}}{k!}

The following identity Eq. (23) given by Gosper [34] simplifies Eq. (22) to Eq. (24)

e^{- f} \sum_{k = 0}^{\infty} \frac{{(f)}^{k} H_{k}}{k!} = γ + Γ (0, f) + \log (f)

m_{\log (I)} = \log (σ_{}^{2}) + \log (s_{W}^{2}) + Γ (0, s_{W}^{2})

where Γ is the incomplete gamma function.

Since s_w=s_I/σ, the mean of log(I) is given by

m_{\log (I)} = {\begin{matrix} \log (s_{I}^{2}) + Γ (0, \frac{s_{I}^{2}}{σ_{}^{2}}) = \log (s_{I}^{2}) + \int_{\frac{s_{I}^{2}}{σ_{n}^{2}}}^{+ \infty} \frac{e^{- t}}{t} d t \\ \log (σ_{}^{2}) - γ \end{matrix}

Therefore, the mean of the natural logarithm of the normalized intensity to the non-centrality parameter (s_I²) is equal to Γ(0, s_I²/σ²).

The variance of log(I) is

σ_{\log (I)}^{2} = E {{[\log (w)]}^{2}} - {[\log (s_{W}^{2}) + Γ (0, s_{W}^{2})]}^{2}

Using Eqs. (8), (17), (19) and the 2nd derivative of the Gamma function, the second moment of log (W) can be rewritten as

E {\log^{2} (w)} = e^{- s_{w}^{2}} \sum_{k = 0}^{\infty} \frac{{(s_{W}^{2})}^{k}}{k! Γ (k + 1)} Γ^{(2)} (k + 1)

The second derivative of the Gamma function can be expressed as [30]

\frac{Γ^{(2)} (k + 1)}{Γ (k + 1)} = ψ_{0}^{2} (k + 1) + \frac{d ψ_{0} (k)}{d k} = ψ_{0}^{2} (k + 1) + ψ_{1} (k + 1)

where ψ₁(k+1) is the trigamma function and is related to the Riemann zeta function ζ and the generalized harmonic number of order k of 2 (H_k⁽²⁾) by [30]

ψ_{1} (k + 1) = ξ (2) - H_{k}^{(2)} = \frac{π^{2}}{6} - H_{k}^{(2)} = \frac{π^{2}}{6} - \sum_{j = 1}^{k} \frac{1}{j^{2}}

Equations (21), (28), (29) and Taylor series expansion of exponential function simplify Eq. (27) to Eq. (30)

E {\log^{2} (w)} = \frac{π^{2}}{6} + γ^{2} - (2 γ) e^{- s_{W}^{2}} \sum_{k = 0}^{\infty} \frac{{(s_{W}^{2})}^{k} H_{k}}{k!} + e^{- s_{W}^{2}} \sum_{k = 0}^{\infty} \frac{{(s_{W}^{2})}^{k} [{(H_{k})}^{2} - H_{k}^{(2)}]}{k!}

By using the Gosper identity (23), the second moment of log(w) is rewritten as

E {\log^{2} (w)} = \frac{π^{2}}{6} + γ^{2} - 2 γ (γ + \log (s_{W}^{2}) + Γ (0, s_{W}^{2})) + F (s_{W}^{2})

where F is defined by

F (s_{W}^{2}) = e^{- s_{w}^{2}} \sum_{k = 0}^{\infty} \frac{{(s_{W}^{2})}^{k} [{(H_{k})}^{2} - H_{k}^{(2)}]}{k!}

To simplify Eq. (32), we define the first derivative of F as

\frac{\partial F (g)}{\partial g} = e^{- g} \sum_{k = 0}^{\infty} \frac{g^{k} [{(H_{k + 1})}^{2} - {(H_{k})}^{2} - H_{k + 1}^{(2)} + H_{k}^{(2)}]}{k!}

where g= s_W² and F(0)=0.

Using the following recursion equations (34) and (35), the first derivative of F can be given by Eq. (36)

H_{k + 1} - \frac{1}{k + 1} = H_{k}

H_{k + 1}^{(2)} - H_{k}^{(2)} = \frac{1}{{(k + 1)}^{2}}

\frac{\partial F (g)}{\partial g} = \frac{2 e^{- g}}{g} \sum_{k = 1}^{\infty} {(g)}^{k} (\frac{H_{k}}{k!} - \frac{1}{k k!})

The entire functions (E_in) [30] given in Eqs. (37), (38) simplify Eq. (36) to Eq. (39) for g>=0

E_{i n} (g) = \int_{0}^{g} \frac{(1 - e^{- t})}{t} d t = e^{- g} \sum_{k = 1}^{\infty} {(g)}^{k} (\frac{H_{k}}{k!}) = γ + Γ (0, g) + \log (g) = \sum_{k = 1}^{\infty} {(- 1)}^{k + 1} (\frac{g^{k}}{k k!})

E_{i n} (- g) = - \sum_{k = 1}^{\infty} \frac{g^{k}}{k k!}

\frac{\partial F (g)}{\partial g} = \frac{2 (1 - e^{- g})}{g} E_{i n} (g) + \frac{2 e^{- g}}{g} (E_{i n} (g) + E_{i n} (- g)) = 2 \frac{\partial E_{i n} (g)}{\partial g} E_{i n} (g) + \frac{2 e^{- g}}{g} (E_{i n} (g) + E_{i n} (- g))

Therefore

F (s_{W}^{2}) = {(γ + \log (s_{W}^{2}) + Γ (0, s_{W}^{2}))}^{2} + 2 \int_{0}^{s_{W}^{2}} \frac{e^{- g}}{g} (E_{i n} (g) + E_{i n} (- g)) d g

Using Eqs. (26), (31), and (40), the variance of log(I) can be expressed as

σ_{\log (I)}^{2} = \frac{π^{2}}{6} + 2 \int_{0}^{s_{W}^{2}} \frac{e^{- g}}{g} (E_{i n} (g) + E_{i n} (- g)) d g

Now, we can simplify Eq. (41) using power series expansions of E_in (Eqs. (37) and (38))

σ_{\log (I)}^{2} = \frac{π^{2}}{6} - 2 \sum_{k = 1}^{\infty} \frac{(1 + {(- 1)}^{k}) γ (k, s_{W}^{2})}{k k!}

where incomplete Gamma functions γ(a,z) and Γ(a,z) are given by [30]

γ (a, z) = \int_{0}^{z} t^{a - 1} e^{- t} d t

Γ (a, z) = \int_{z}^{\infty} t^{a - 1} e^{- t} d t

Using equalities given in Eq. (45) [30], the variance of log(I) can rewritten as

Γ (a) = (a - 1)! = γ (a, z) + Γ (a, z)

σ_{\log (I)}^{2} = \frac{π^{2}}{6} - 2 \sum_{k = 1}^{\infty} \frac{(1 + {(- 1)}^{k}) (Γ (k) - Γ (k, s_{W}^{2}))}{k k!} = \frac{π^{2}}{6} - 2 \sum_{k = 1}^{\infty} \frac{(1 + {(- 1)}^{k})}{k^{2}} + 4 \sum_{n = 1}^{\infty} \frac{Γ (2 n, s_{W}^{2})}{(2 n) (2 n)!}

The following power series can be expressed as a Riemann zeta function [30]

\sum_{k = 1}^{\infty} \frac{(1 + {(- 1)}^{k})}{k^{2}} = \frac{ξ (2)}{2} = \frac{π^{2}}{12}

and the variance of log(I) can be given by

σ_{\log (I)}^{2} = 4 \sum_{n = 1}^{\infty} \frac{Γ (2 n, \frac{s_{I}^{2}}{σ_{}^{2}})}{(2 n) (2 n)!} = 4 \sum_{n = 1}^{\infty} \frac{Γ (2 n, S N R)}{(2 n) (2 n)!}

The variance of log(I), LOGIV, can also be given by

L O G I V = σ_{\log (I)}^{2} = \sum_{n = 1}^{\infty} \frac{Q (2 n, S N R)}{n^{2}}

where Q is the normalized incomplete Gamma function [30]

Q (a, b) = \frac{Γ (a, b)}{Γ (a)}

It follows that the variance of differences of logarithm intensities (two independent random processes with equal variances) is

D L O G I V = 2 σ_{\log (I)}^{2} = 2 \sum_{n = 1}^{\infty} \frac{Q (2 n, S N R)}{n^{2}}

Acknowledgments

This work was supported in part by research from the California Institute for Regenerative Medicine (CIRM). The authors thank Carl Zeiss Meditec for providing the patient interface part of their Stratus OCT system.

References and links

1. V. Patel, S. Rassam, R. Newsom, J. Wiek, and E. Kohner, “Retinal blood flow in diabetic retinopathy,” BMJ 305(6855), 678–683 (1992). [CrossRef] [PubMed]

2. E. Friedman, “A hemodynamic model of the pathogenesis of age-related macular degeneration,” Am. J. Ophthalmol. 124(5), 677–682 (1997). [PubMed]

3. J. Flammer, S. Orgül, V. P. Costa, N. Orzalesi, G. K. Krieglstein, L. M. Serra, J.-P. Renard, and E. Stefánsson, “The impact of ocular blood flow in glaucoma,” Prog. Retin. Eye Res. 21(4), 359–393 (2002). [CrossRef] [PubMed]

4. S. S. Hayreh, Anterior Ischemic Optic Neuropathy, (Springer-Verlag, 1975).

5. J. D. Gass, Stereoscopic Atlas of Macular Diseases, 4th ed. (Mosby, 1997).

6. L. A. Yannuzzi, K. T. Rohrer, L. J. Tindel, R. S. Sobel, M. A. Costanza, W. Shields, and E. Zang, “Fluorescein angiography complication survey,” Ophthalmology 93(5), 611–617 (1986). [PubMed]

7. M. Hope-Ross, L. A. Yannuzzi, E. S. Gragoudas, D. R. Guyer, J. S. Slakter, J. A. Sorenson, S. Krupsky, D. A. Orlock, and C. A. Puliafito, “Adverse reactions due to indocyanine green,” Ophthalmology 101(3), 529–533 (1994). [PubMed]

8. W. Drexler and J. G. Fujimoto, Optical Coherence Tomography: Technology and Applications, Biological and Medical Physics, Biomedical Engineering (Springer, 2008).

9. M. Wojtkowski, V. Srinivasan, J. G. Fujimoto, T. Ko, J. S. Schuman, A. Kowalczyk, and J. S. Duker, “Three-dimensional retinal imaging with high-speed ultrahigh-resolution optical coherence tomography,” Ophthalmology 112(10), 1734–1746 (2005). [CrossRef] [PubMed]

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

11. T. Schmoll, C. Kolbitsch, and R. A. Leitgeb, “Ultra-high-speed volumetric tomography of human retinal blood flow,” Opt. Express 17(5), 4166–4176 (2009). [CrossRef] [PubMed]

12. J. Fingler, R. J. Zawadzki, J. S. Werner, D. Schwartz, and S. E. Fraser, “Volumetric microvascular imaging of human retina using optical coherence tomography with a novel motion contrast technique,” Opt. Express 17(24), 22190–22200 (2009). [CrossRef] [PubMed]

13. G. Liu, W. Qi, L. Yu, and Z. Chen, “Real-time bulk-motion-correction free Doppler variance optical coherence tomography for choroidal capillary vasculature imaging,” Opt. Express 19(4), 3657–3666 (2011). [CrossRef] [PubMed]

14. G. Liu, L. Chou, W. Jia, W. Qi, B. Choi, and Z. Chen, “Intensity-based modified Doppler variance algorithm: application to phase instable and phase stable optical coherence tomography systems,” Opt. Express 19(12), 11429–11440 (2011). [CrossRef] [PubMed]

15. S. Makita, J. Franck, M. Yamanari, M. Miura, and Y. Yasuno, “Comprehensive in vivo micro-vascular imaging of the human eye by dual-beam-scan Doppler optical coherence angiography,” Opt. Express 19(2), 1271–1283 (2011). [CrossRef] [PubMed]

16. A. Szkulmowska, M. Szkulmowski, D. Szlag, A. Kowalczyk, and M. Wojtkowski, “Three-dimensional quantitative imaging of retinal and choroidal blood flow velocity using joint Spectral and Time domain Optical Coherence Tomography,” Opt. Express 17(13), 10584–10598 (2009). [CrossRef] [PubMed]

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

18. M. C. Pierce, B. Hyle Park, B. Cense, and J. F. de Boer, “Simultaneous intensity, birefringence, and flow measurements with high-speed fiber-based optical coherence tomography,” Opt. Lett. 27(17), 1534–1536 (2002). [CrossRef] [PubMed]

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

20. R. K. Wang, L. An, P. Francis, and D. J. Wilson, “Depth-resolved imaging of capillary networks in retina and choroid using ultrahigh sensitive optical microangiography,” Opt. Lett. 35(9), 1467–1469 (2010). [CrossRef] [PubMed]

21. R. K. Wang and L. An, “Multifunctional imaging of human retina and choroid with 1050-nm spectral domain optical coherence tomography at 92-kHz line scan rate,” J. Biomed. Opt. 16(5), 050503 (2011). [CrossRef] [PubMed]

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

23. A. Mariampillai, M. K. K. Leung, M. Jarvi, B. A. Standish, K. Lee, B. C. Wilson, A. Vitkin, and V. X. D. Yang, “Optimized speckle variance OCT imaging of microvasculature,” Opt. Lett. 35(8), 1257–1259 (2010). [CrossRef] [PubMed]

24. J. W. Goodman, Statistical Optics, (Wiley, 1985).

25. J. M. Schmitt, S. H. Xiang, and K. M. Yung, “Speckle in optical coherence tomography,” J. Biomed. Opt. 4(1), 95 (1999). [CrossRef]

26. M. Pircher, E. Gotzinger, R. Leitgeb, A. F. Fercher, and C. K. Hitzenberger, “Speckle reduction in optical coherence tomography by frequency compounding,” J. Biomed. Opt. 8(3), 565–569 (2003). [CrossRef] [PubMed]

27. B. Karamata, K. Hassler, M. Laubscher, and T. Lasser, “Speckle statistics in optical coherence tomography,” J. Opt. Soc. Am. A 22(4), 593–596 (2005). [CrossRef] [PubMed]

28. A. Papoulis and S. U. Pillai, Probability, Random Variables and Stochastic Processes, (McGraw-Hill, 2002).

29. J. D. Briers, “Laser speckle contrast imaging for measuring blood flow,” in Proceeding of the Symposium on Photonics Technologies for 7th Framework Program (2006), pp. 328–332.

30. M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, (Dover, 1964).

31. E. C. Lee, J. F. de Boer, M. Mujat, H. Lim, and S. H. Yun, “In vivo optical frequency domain imaging of human retina and choroid,” Opt. Express 14(10), 4403–4411 (2006). [CrossRef] [PubMed]

32. Y. Yasuno, V. D. Madjarova, S. Makita, M. Akiba, A. Morosawa, C. Chong, T. Sakai, K.-P. Chan, M. Itoh, and T. Yatagai, “Three-dimensional and high-speed swept-source optical coherence tomography for in vivo investigation of human anterior eye segments,” Opt. Express 13(26), 10652–10664 (2005). [CrossRef] [PubMed]

33. B. Cense, N. Nassif, T. Chen, M. Pierce, S. H. Yun, B. Park, B. Bouma, G. Tearney, and J. de Boer, “Ultrahigh-resolution high-speed retinal imaging using spectral-domain optical coherence tomography,” Opt. Express 12(11), 2435–2447 (2004). [CrossRef] [PubMed]

34. R. W. Gosper, “harmonic Summation and exponential gfs,” http://mathworld.wolfram.com/HarmonicNumber.html.

Logarithmic intensity and speckle-based motion contrast methods for human retinal vasculature visualization using swept source optical coherence tomography

Abstract

1. Introduction

2. Theory

2.1 Speckle amplitude and intensity variances

2.2 Speckle amplitude and intensity contrast ratios

2.3 Logarithmic intensity and differential logarithmic intensity variances

2.4 Comparison of speckle- and logarithmic intensity-based contrasts

3. Experimental methods

3.1 SS-OCT System

3.2 Scanning protocols

3.3 Image processing and motion contrast imaging

4. Experimental results

4.1 Experimental analysis of the linear intensity through retinal and choroidal layers

4.2 Quantitative analysis of speckle contrast ratios and logarithmic intensity-based contrasts though the retinal and choroidal layers

4.3 2D tomogram and en face view visualization of the retina using motion contrast imaging methods

5. Discussion

6. Appendix

Acknowledgments

References and links

Cited By

Figures (8)

Equations (58)

Biomedical Optics Express