Blood flow velocity quantification using split-spectrum amplitude-decorrelation angiography with optical coherence tomography

The split-spectrum amplitude-decorrelation angiography (SSADA) algorithm was recently developed as a method for imaging blood flow in the human retina without the use of phase information. In order to enable absolute blood velocity quantification, in vitro phantom experiments are performed to correlate the SSADA signal at multiple time scales with various preset velocities. A linear model relating SSADA measurements to absolute flow velocities is derived using the phantom data. The operating range for the linear model is discussed along with its implication for velocity quantification with SSADA in a clinical setting. ©2013 Optical Society of America OCIS codes: (170.4500) Optical coherence tomography; (170.3880) Medical and biological imaging; (170.4470) Ophthalmology; (999.999) Blood flow. References and links 1. D. Huang, E. A. Swanson, C. P. Lin, J. S. Schuman, W. G. Stinson, W. Chang, M. R. Hee, T. Flotte, K. Gregory, C. A. Puliafito, and J. G. Fujimoto, “Optical coherence tomography,” Science 254(5035), 1178–1181 (1991). 2. M. R. Hee, J. A. Izatt, E. A. Swanson, D. Huang, J. S. Schuman, C. P. Lin, C. A. Puliafito, and J. G. Fujimoto, “Optical coherence tomography of the human retina,” Arch. Ophthalmol. 113(3), 325–332 (1995). 3. M. Wojtkowski, R. Leitgeb, A. Kowalczyk, T. Bajraszewski, and A. F. Fercher, “In vivo human retinal imaging by Fourier domain optical coherence tomography,” J. Biomed. Opt. 7(3), 457–463 (2002). 4. 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). 5. 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). 6. E. Friedman, “A hemodynamic model of the pathogenesis of age-related macular degeneration,” Am. J. Ophthalmol. 124(5), 677–682 (1997). 7. Y. Zhao, Z. Chen, C. Saxer, S. Xiang, J. F. de Boer, and J. S. Nelson, “Phase-resolved optical coherence tomography and optical Doppler tomography for imaging blood flow in human skin with fast scanning speed and high velocity sensitivity,” Opt. Lett. 25(2), 114–116 (2000). 8. V. X. D. Yang, M. L. Gordon, A. Mok, Y. Zhao, Z. Chen, R. S. C. Cobbold, B. C. Wilson, and I. Alex Vitkin, “Improved phase-resolved optical Doppler tomography using the Kasai velocity estimator and histogram segmentation,” Opt. Commun. 208(4-6), 209–214 (2002). 9. R. Leitgeb, L. Schmetterer, W. Drexler, A. Fercher, R. Zawadzki, and T. Bajraszewski, “Real-time assessment of retinal blood flow with ultrafast acquisition by color Doppler Fourier domain optical coherence tomography,” Opt. Express 11(23), 3116–3121 (2003). 10. M. A. Choma, A. K. Ellerbee, S. Yazdanfar, and J. A. Izatt, “Doppler flow imaging of cytoplasmic streaming using spectral domain phase microscopy,” J. Biomed. Opt. 11(2), 024014 (2006). 11. B. Baumann, B. Potsaid, M. F. Kraus, J. J. Liu, D. Huang, J. Hornegger, A. E. Cable, J. S. Duker, and J. G. Fujimoto, “Total retinal blood flow measurement with ultrahigh speed swept source/Fourier domain OCT,” Biomed. Opt. Express 2(6), 1539–1552 (2011). 12. R. Leitgeb, C. Hitzenberger, and A. Fercher, “Performance of fourier domain vs. time domain optical coherence tomography,” Opt. Express 11(8), 889–894 (2003). (C) 2013 OSA 1 October 2013 | Vol. 4, No. 10 | DOI:10.1364/BOE.4.001909 | BIOMEDICAL OPTICS EXPRESS 1909 #193860 $15.00 USD Received 17 Jul 2013; revised 22 Aug 2013; accepted 27 Aug 2013; published 3 Sep 2013 13. 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). 14. R. K. Wang and Z. Ma, “Real-time flow imaging by removing texture pattern artifacts in spectral-domain optical Doppler tomography,” Opt. Lett. 31(20), 3001–3003 (2006). 15. Y. J. Hong, S. Makita, F. Jaillon, M. J. Ju, E. J. Min, B. H. Lee, M. Itoh, M. Miura, and Y. Yasuno, “Highpenetration swept source Doppler optical coherence angiography by fully numerical phase stabilization,” Opt. Express 20(3), 2740–2760 (2012). 16. H. Ren and X. Li, “Clutter rejection filters for optical Doppler tomography,” Opt. Express 14(13), 6103–6112 (2006). 17. J. Tokayer and D. Huang, “Effect of blood vessel diameter on relative blood flow estimates in Doppler optical coherence tomography algorithms,” Proc. SPIE 7889, Optical Coherence Tomography and Coherence Domain Optical Methods in Biomedicine XV, 78892X, 78892X-8 (2011). 18. Y. Wang, B. A. Bower, J. A. Izatt, O. Tan, and D. Huang, “Retinal blood flow measurement by circumpapillary Fourier domain Doppler optical coherence tomography,” J. Biomed. Opt. 13(6), 064003 (2008). 19. M. Szkulmowski, A. Szkulmowska, T. Bajraszewski, A. Kowalczyk, and M. Wojtkowski, “Flow velocity estimation using joint Spectral and Time domain Optical Coherence Tomography,” Opt. Express 16(9), 6008– 6025 (2008). 20. Y. K. Tao, A. M. Davis, and J. A. Izatt, “Single-pass volumetric bidirectional blood flow imaging spectral domain optical coherence tomography using a modified Hilbert transform,” Opt. Express 16(16), 12350–12361 (2008). 21. I. Grulkowski, I. Gorczynska, M. Szkulmowski, D. Szlag, A. Szkulmowska, R. A. Leitgeb, A. Kowalczyk, and M. Wojtkowski, “Scanning protocols dedicated to smart velocity ranging in spectral OCT,” Opt. Express 17(26), 23736–23754 (2009). 22. S. Makita, Y. Hong, M. Yamanari, T. Yatagai, and Y. Yasuno, “Optical coherence angiography,” Opt. Express 14(17), 7821–7840 (2006). 23. L. An and R. K. Wang, “In vivo volumetric imaging of vascular perfusion within human retina and choroids with optical micro-angiography,” Opt. Express 16(15), 11438–11452 (2008). 24. 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). 25. L. Yu and Z. Chen, “Doppler variance imaging for three-dimensional retina and choroid angiography,” J. Biomed. Opt. 15(1), 016029 (2010). 26. 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). 27. J. Fingler, D. Schwartz, C. Yang, and S. E. Fraser, “Mobility and transverse flow visualization using phase variance contrast with spectral domain optical coherence tomography,” Opt. Express 15(20), 12636–12653 (2007). 28. 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). 29. Y. Yasuno, Y. Hong, S. Makita, M. Yamanari, M. Akiba, M. Miura, and T. Yatagai, “In vivo high-contrast imaging of deep posterior eye by 1um swept source optical coherence tomography and scattering optical coherence angiography,” Opt. Express 15(10), 6121–6139 (2007). 30. Y. Hong, S. Makita, M. Yamanari, M. Miura, S. Kim, T. Yatagai, and Y. Yasuno, “Three-dimensional visualization of choroidal vessels by using standard and ultra-high resolution scattering optical coherence angiography,” Opt. Express 15(12), 7538–7550 (2007). 31. 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). 32. H. C. Hendargo, R. Estrada, S. J. Chiu, C. Tomasi, S. Farsiu, and J. A. Izatt, “Automated non-rigid registration and mosaicing for robust imaging of distinct retinal capillary beds using speckle variance optical coherence tomography,” Biomed. Opt. Express 4(6), 803–821 (2013). 33. E. Jonathan, J. Enfield, and M. J. Leahy, “Correlation mapping method for generating microcirculation morphology from optical coherence tomography (OCT) intensity images,” J Biophotonics 4(9), 583–587 (2011). 34. 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). 35. Y. Jia, J. C. Morrison, J. Tokayer, O. Tan, L. Lombardi, B. Baumann, C. D. Lu, W. Choi, J. G. Fujimoto, and D. Huang, “Quantitative OCT angiography of optic nerve head blood flow,” Biomed. Opt. Express 3(12), 3127– 3137 (2012). 36. E. Wei, Y. Jia, O. Tan, B. Potsaid, J. J. Liu, W. Choi, J. G. Fujimoto, and D. Huang, “Parafoveal retinal vascular response to pattern visual stimulation assessed with OCT angiography,” PLoS ONE. submitted. 37. W. Li, C. T. Lancee, E. I. Cespedes, A. F. W. van der Steen, and N. Bom, “Decorrelation properties of intravascular echo signals: Potentials for blood velocity estimation,” J. Acoust. Soc. Am. 71(14), 70D–86D (1993). (C) 2013 OSA 1 October 2013 | Vol. 4, No. 10 | DOI:10.1364/BOE.4.001909 | BIOMEDICAL OPTICS EXPRESS 1910 #193860 $15.00 USD Received 17 Jul 2013; revised 22 Aug 2013; accepted 27 Aug 2013; published 3 Sep 2013 38. W. Li, A. F. W. van der Steen, C. T. Lancée, I. Céspedes, and N. Bom, “Blood flow imaging and volume flow quantitation with intravascular ultrasound,” Ultrasound Med. Biol. 24(2), 203–214 (1998). 39. 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). 40. G. Liu, W. Jia, V. Sun, B. Choi, and Z. Chen, “High-resolut


Introduction
Optical coherence tomography (OCT) is a high-resolution modality used for non-invasive depth-resolved imaging of biological tissue [1].OCT has been used extensively in ophthalmology for imaging both structure and function in the human retina [2][3][4].Of particular importance for studying various retinal pathologies is the measurement of retinal blood flow [5,6].Doppler OCT is able to calculate absolute blood flow velocity by evaluating the phase differences between adjacent A-scans [7][8][9][10][11].The speed and sensitivity advantages that Fourier-domain systems, which includes both spectral-domain and swept-source systems, demonstrate relative to time-domain systems make them particularly suited for this purpose [12].Additionally, swept-source OCT has a larger Doppler dynamic range than spectraldomain OCT due to a higher fringe washout velocity and is thus suitable for measuring very fast flow speeds [13].However, Doppler OCT requires good phase stability for quantifying slow flow velocities [13][14][15] and suffers from reduced sensitivity for small blood vessels [16,17].Furthermore, the measured Doppler frequency shift in a blood vessel varies inversely with the cosine of the Doppler angle, or the angle between the incident beam and blood vessel [18].Yet many blood vessels in the human retina are nearly perpendicular to the incident beam, so that the Doppler frequency shifts in these vessels are small and difficult to detect.Several techniques, including joint spectral and time domain OCT [19], use of a modified Hilbert transform [20] and smart scanning protocols [21], have been introduced to overcome some of the limitations of Doppler OCT.However, techniques that do not inherently depend on the Doppler angle may be particularly useful for visualizing blood vessels in the human retina.
Several angiographic techniques have been introduced for imaging retinal blood flow.Some of these techniques are phase-based, while others are amplitude-based.Some of the phase-based techniques include optical coherence angiography [22], optical microangiography [23,24], Doppler variance [25,26] and phase variance [27,28].Some of the amplitude-based techniques include scattering optical coherence angiography [29,30], speckle variance [31,32], correlation mapping [33] and split-spectrum amplitude-decorrelation angiography [34][35][36].These methods have been implemented using both spectral-domain and swept-source OCT imaging systems.Most of these techniques can detect both transverse and axial flow, and have been successful in visualizing retinal and choroidal microvascular networks.While good qualitative imaging results have been shown for all of these methods, quantitative results that map angiograms to flow velocities are lacking.Decorrelation has been used for quantitative flow measurements in ultrasound [37,38] and thus has the potential to be useful for measuring flow velocities in OCT.Recently, a decorrelation-based method termed intensity-based Doppler variance (IBDV) was introduced and its relationship with velocity was established [39][40][41].This method computes decorrelation using only the amplitude signal.The authors in [41] showed that the IBDV signal increased by 87.5% as the Doppler angle increased from 0° to 18°, indicating a significant Doppler angle dependence.One potential advantage of SSADA is that the algorithm first digitally creates an isotropic coherence volume, or resolution cell, prior to computing decorrelation.This could make the algorithm equally sensitive to axial and transverse motion so that SSADA may be used to quantify flow independent of Doppler angle.
The purpose of this paper is to determine a relationship between SSADA measurements and flow velocity.We hypothesize that this relationship is linear within a certain range of velocities and use in vitro blood flow phantom experiments to test this hypothesis.Whole blood was used in order to to closely mimic in vivo retinal imaging.The SSADA algorithm computes decorrelation between two OCT amplitudes that are taken at the same scanning position but are separated in time.We first examine the dependence of SSADA measurements on Doppler angle after splitting the source spectrum in wavenumber so that the axial resolution of the OCT imaging system is equal to its transverse resolution.The concept of multi-timescale SSADA is introduced, where the time separation between amplitude measurements is varied, and used to examine the time-dependence of decorrelation.We derive an equation that can be used to calculate the absolute flow velocity from a measured decorrelation value and particular time separation between OCT amplitude measurements.We then define the saturation velocity for which the decorrelation values saturate and the linear relationship is no longer valid.

System setup
A detailed description of the experimental Fourier domain OCT system can be found in [20].This system, originally designed for retinal imaging, was modified for capillary tube measurements by adding a focusing lens in the sample arm.Briefly, it contains a superluminescent diode source with a center wavelength of 840 μm and a bandwidth of 49 μm.The theoretical axial resolution of the system was 6.4 μm in air, while the measured axial resolution was 8.8 μm in air.The deviation from the theoretical value is likely due to the non-Gaussian spectrum profile as well as aberrations that reduce the spectral resolution of the spectrometer.A collimator was placed in the sample arm that produced a 1/e 2 intensity diameter of approximately 1.0 mm, and was focused down to 20 μm spot size using a 20 mm focal length lens (Thorlabs LB/450-B Bi-convex lens with antireflection coating in 650 -1050 nm range).The sample arm probe beam was focused onto a glass capillary tube (Wale Apparatus) with an outer diameter of 330 μm and an inner diameter of 200 μm.No antireflection material was used to coat the glass.The power incident on the capillary was 500 μW.The capillary was placed on top of a piece of paper and attached to a ball and socket mounting platform (Thorlabs).A syringe pump (Harvard Apparatus) was used to control the flow of human blood (coauthor J.T.) through the tube.The recombined field was spectrally dispersed by a spectrometer and detected by 1024 pixels on a line-scan camera.The wavelengths acquired by the spectrometer ranged from 806 nm to 876 nm in steps of 0.068 nm.The data from the camera was transferred to a computer for data processing.The time between two sequential A-scan acquisitions was 56 μs which corresponds to a 17 kHz repetition rate.A system schematic is show in Fig. 1.Fig. 1.Spectral domain OCT system schematic.SLD -superluminescent diode, FC -2x2 fiber coupler, PC -personal computer, H 2 0 -water cell, Galvos -scanning mirror galvanometers, BS beamsplitter.This system was originally designed for retinal imaging so that a water cell is in the reference arm for dispersion matching and a beamsplitter is in the sample arm to allow for slit lamp illumination of the retina.The retinal imaging system was modified by adding a focusing lens in the sample arm.

Scan protocol
Data was collected using both dual-plane [42] and M-mode protocols.The dual-plane protocol was implemented with two parallel B-scans that were separated by 100 μm and repeated eight times.The B-scans each consisted of 700 A-lines and covered 700 μm in the transverse dimension, resulting in a transverse step size of 1 μm.M-mode scans consisted of 2600 A-lines at the same transverse position inside the capillary.Measurements were taken for five different preset flow rates and five different Doppler angles.The Doppler angle was controlled using the ball and socket mount.

Data processing
The acquisition software was provided by Bioptigen Inc. (Durham, NC), and the custom processing and analysis software was coded in MATLAB.The custom processing software included methods for resampling to k-space via spline interpolation and dispersion compensation.Mirror reflections at two different sample depths were acquired and used for resampling to k-space as described in [43].A dispersion mismatch between the reference and sample arms leads to a nonlinear phase dependence of the spectral interferogram in wavenumber [44].To compensate dispersion, a linear function was fit to the unwrapped phase of one of the mirror reflections and the residual between the phase and the linear fit, representing the nolinear dependence, was used as a phase correction factor ( ) ( ) ( )  The dual-plane B-scans were used to compute the Doppler angle between the probe beam and the glass capillary, while the M-mode scans were used for velocity and decorrelation calculations.The velocities were calculated by first computing Doppler frequency shifts using the phase-resolved Doppler OCT algorithm [45] and then converting to absolute velocity using the measured Doppler angle.The capillary tube was automatically segmented in the Mmode scans using intensity thresholding and morphological operations.Figure 3 illustrates that there are significant projection artifacts in the paper underneath the glass capillary, which is likely due to multiple scattering by the moving blood cells.To avoid the influence of these artifacts on the measured Doppler frequency shifts, we chose a subregion in the upper half of the capillary to use for all subsequent calculations.The SSADA algorithm theory has been explained previously [34].It is based on the concept of creating an isotropic voxel by degrading the axial resolution until the coherence volume has equal size in the axial and transverse dimensions.This is accomplished by splitting the interferogram into different spectral bands with reduced bandwidth, which leads to a broadening of the axial resolution.The details of the algorithm, as well as its multiscale extension, are explained in the next section.

Theory
The OCT spectral interferogram can be written as [44] ( ) where 2 / k π λ = denotes the wavenumber, Re{} denotes the real part operator, ( , ) x z denotes the coordinate of a reference frame fixed to the sample, ( ) S k is the source power spectral density, ( , )  r x z is the backscattering coefficient of the sample and 0 ω is the full-width at half-maximum (FWHM) of the beam intensity profile.Note that Eq. ( 1) assumes that the zero path length difference in the interferometer is located at 0 z = and that the beam is centered at transverse position 0 x = .For a Gaussian spectral shape we can write is the center wavenumber of the source and k Δ is its FWHM spectral width in wavenumber.
After Fourier transforming the spectral interferogram we obtain the complex OCT signal as a function of depth Z , which can be written as where 0 ln 2 4 / z k δ = Δ denotes the FWHM axial resolution.
Each pixel in an OCT image is formed by a coherent sum of all of the light backscattered from a 3D coherence volume within a sample.The dimensions of the coherence volume are determined by 0 ω and 0 z δ , as can be seen by examining the rightmost exponential in Eq. ( 3).
The explicit dependence on the FWHM 0 ω of the intensity profile, rather than the electric field profile has been explained in [46].Briefly, the sample arm fiber's mode field modulates the illumination field, which the combines with the backscattering coefficient to produce the scattering field.The portion of the scattering field that couples back into the fiber is given by an overlap integral between the scattered field and the fiber mode field.Thus, the backscattering coefficients in the coherent sum are modulated by a function that is proportional to the square of the fiber mode field, which equals the intensity at the sample.
In our current OCT system the FWHM axial resolution 0 6.5 z δ ≈ µm is approximately two times higher than the FWHM transverse resolution 0 11.8 ω ≈ µm, which leads to higher decorrelation sensitivity for axial motion [47].In the fundus, blood flow is primarily in the transverse direction, whereas bulk motion noise sources, such as pulsation related to heartbeat, occur primarily in the axial direction [34].In order to decrease the sensitivity to this axial motion noise and improve the signal-to-noise ratio of flow detection, the SSADA algorithm digitally broadens the axial resolution prior to computing speckle decorrelation.
In the SSADA implementation, a desired axial resolution 0 z δ ′ is chosen that is at least as large as the nominal axial resolution of the OCT imaging system, which is 6.5 μm in blood.In order to create an isotropic coherence volume we set 0 0 11.8 z δ ω ′ = = μm.A set of Gaussian filters in wavenumber is then created that, when multiplied with the interferometric signal in Eq. ( 1), increases the axial resolution from 0 z δ to 0 z δ ′ .The j th filter in the filterbank can be expressed where 0 j k is the center wavenumber and G k Δ is its FWHM bandwidth.In order to determine the necessary G k Δ a few Fourier transform and Gaussian function properties are leveraged [48].The source spectrum is approximated with an equivalent Gaussian spectrum to simplify the subsequent mathematical analysis.We note that when multiplying Eq. ( 1) by Eq. ( 4) we are taking the product of two Gaussians in wavenumber space.This implies that, by the convolution property of Fourier transforms, we are convolving the original Gaussian response in Eq. ( 3) with a Gaussian that has FWHM equal to 4 ln 2 / G k Δ .The convolution results in another Gaussian in z-space with FWHM equal to . This is the quantity that we want to set, as it defines the modified FWHM axial resolution resulting from a spectrum split.Note the first term in the square root is equal to the original axial resolution 0 z δ .Thus, in order to increase the axial resolution from 0 . For our system this bandwidth is 285 radians/mm.In order to complete the filterbank specification, we next define the distance in wavenumber between adjacent filters as 0.5 G k × Δ and then determine the center wavenumber of each filter so that all of them fit within the range of acquired wavenumbers.The wavenumbers acquired by the spectrometer in our system ranged from 717 to 779 radians/mm and we are able to fit five Gaussian filters in this range with centers 719, 733, 747, 761 and 776 radians/mm, respectively.
The SSADA algorithm computes decorrelation at a single voxel between OCT magnitudes

( ) | ( ) |
A Z I Z = separated in time.In the original SSADA implementation in [34], consecutive B-scans at the same slow axis scanning position (M-B scans) and separated by 2.0 t Δ = ms are used to compute a decorrelation frame.In our current experiments, on the other hand, we use M-mode scans to compute the decorrelation of a single A-line.The main purpose of the switch to M-mode imaging is that it allows us to study the time course of decorrelation.We let ( ) A Z denote the m th split spectra from the n th A-line in our M-mode scan that is taken at time ( 1)  t n τ = − , where 56 τ = µs is the A-line rate of our system.
Generalizing the average decorrelation equation given in [34] to consider different time separations, the average decorrelation between 2 A-lines taken at the same transverse position can be written as where 2000 N = is the number of individual decorrelations that are being averaged, 5 M = is the number of split-spectra/Gaussian filters that was previously described and t Δ is an integer multiple of the line rate τ .Using our M-mode scans, we can compute a multi- timescale SSADA (MSSADA) image by computing the decorrelation at time separations ranging from τ through Nτ .A sample multi-timescale SSADA image is illustrated in Fig. 4.
(C) 2013 OSA For all of the results presented in this work, the Doppler and SSADA signals were binned into four distinct equally space depth bins within the region of.The data within each bin was averaged so that each scan produced four Doppler measurements and four SSADA measurements.

Doppler angle dependence
Since SSADA is an amplitude-based method and we create an isotropic voxel by splitting the source spectrum, we expect that the measured SSADA signal will exhibit similar sensitivities in the axial and transverse directions.In order to test this claim, we analyze the variation of the decorrelation signal with Doppler angle.To avoid both noise and saturation artifacts, we choose time separations t Δ for which decorrelation values lie in a central region around 0.145 at zero Doppler angle.Specifically, time separations of t Δ = 784, 280 and 168 μs for respective flow speeds of 0.37, 1.29 and 2.00 mm/s were used to plot SSADA measurements versus Doppler angle.Figure 5 illustrates that while the decorrelation signal remains relatively constant for Doppler angles less than 10°, it increases above 10° and then appears to plateau above 20°.This indicates that while SSADA measurements are equally sensitive to axial and transverse flow for Doppler angles less than 10°, the measurements are more sensitive to axial flow for Doppler angles above 10°.ρ is a nonparametric measure of how well the relationship between decorrelation and Doppler angle can be described using a monotonic function.In this case the function appears to be sigmoidal.
The change in decorrelation due to an angular variation can be used as an quantitative indicator of the sensitivity of the decorrelation signal to Doppler angle.Specifically, the decorrelation increases from approximated 0.144 at near perpendicular incidence to 0.159 at a Doppler angle of 30°.The change in decorrelation due to an angular deviation of 30°, is then computed as 100 x (0.159 -0.144)/((0.159+ 0.144) / 2) = 9.9%, which indicates a small but significant dependence of decorrelation measurements on Doppler angle.

Saturation
SSADA effectively enumerates the dissimilarity between a pixel's amplitude at two different time instances.If the interval between the two measurements is long enough, the respective amplitudes will be independent, and the decorrelation signal will be fully saturated.This defines a state of complete decorrelation, and any increase in the time interval will not alter the SSADA measurement.Thus, only decorrelation values that are below the saturation level are useful for distinguishing between varying flow speeds.By visually inspecting Fig. 6, we can see that complete decorrelation occurs in approximately 500 μs for a flow speed equal to 2 mm/s.At this speed and time separation the red blood cells (RBCs) are displaced by only 1.0 μm, less than one-tenth of the coherence volume size.This suggests that the decorrelation reaches full saturation well before the RBCs move completely through the coherence volume, which indicates a high sensitivity to speckle variation caused by the RBCs moving through the coherence volume.Note that the curves in Fig. 6 asymptotically approach the complete decorrelation value of 0.21.This motivates us to define a threshold over which the decorrelation rate slows down considerably and the curves in Fig. 6. begin to flatten.We set this decorrelation saturation threshold to 85% of the asymptotic decorrelation value.The resulting threshold is 0.18, and all decorrelation values above this threshold are referred to as saturated.Fig. 6.Multi-timescale decorrelation (Eq.( 5)) for various flow speeds.The asymptotic decorrelation value for our experiments is 0.21.

Relationship between decorrelation and velocity
The goal of this work is to be able to relate the velocity of blood flow to a measured decorrelation value.In order to do so, we examine the velocity dependence of the decorrelation signal at various time separations/intervals.We expect that, for a fixed time separation, the decorrelation will increase with increasing velocity.However, there should be some minimum velocity beneath which measured decorrelations are indistinguishable from noise, and that this velocity should be related to various system noise sources.Furthermore, there should also be a maximum velocity above which measured decorrelations are saturated.As discussed in the section on saturation, the rate of decorrelation begins to slow before decorrelation values are fully saturated.This means that the relationship between decorrelation and velocity changes as the decorrelation approaches full saturation.Therefore, in order to determine a relationship between decorrelation and velocity that does not change with time, we exclude all decorrelation values above the previously defined decorrelation saturation threshold of 0.18 from further analysis.With the saturated points removed, we fit a linear model to the decorrelation versus velocity relationship for a given time separation.The data and fits are shown in Fig. 7.We summarize the fitting parameters in Table 1.From this table we see that the slope changes with time separation.On the other hand, the intercept values remain relatively constant with changing time separations.We thus establish the linear relationship (1) (1) where t D Δ is the measured decorrelation value at a particular time separation t Δ , v is the flow velocity, (1)   t m Δ is the slope parameter that is a function of t Δ and (1) 0.08 b ≈ is the intercept parameter.The significance of this intercept parameter, which is equal to the decorrelation when the velocity is zero, is treated in the Discussion section.We next examined the dependence of the slopes (1)   t m Δ on time separation t Δ and found it to be linear, as decribed by the equation (Fig. 8).Since the intercept was very small, the equation could be simplified to (1)   (2) 0.24  Plugging this relationship into Eq.( 6) gives us the decorrelation as a function velocity and time separation (2) ( In practice, we measure the decorrelation at a particular time separation and wish to find the flow velocity.Thus, we can invert Eq. ( 7) to solve for velocity.Substituting m for (2) 1/ m and b for (1)  b we can write This model is only valid for a specific range of velocities and time separations, which define an operating range for our model.Using Eq. ( 8), we can compute the saturation velocity SAT v , which is defined as the velocity at which the decorrelation reaches the saturation cutoff value of 0.18, for various time separations t Δ .The linear model in Eq. ( 8) does not hold for velocities above SAT v .Some time separation-saturation velocity pairs are illustrated in Table 2.

Model parameters
In order to study the parameters of our model, we first rewrite Eq. ( 8) as where x Δ is the distance that the RBCs move between scans separated by time interval t Δ and all other terms have been previously defined.The two parameters in this model are the slope m and the decorrelation intercept b .Since decorrelation is dimensionless, the parameter b must be dimensionless as well.Furthermore, we see from Eq. ( 9) that when the RBC displacement equals zero, the decorrelation equals b .Thus the parameter b is equal to the decorrelation in non-flow pixels, or the bulk decorrelation.It equals the minimum measurable decorrelation value and can be defined as the decorrelation noise floor.We expect that this parameter will vary inversely with the system signal-to-noise ratio, similar to the way that the phase noise floor does for Doppler OCT [49].Further experiments are required to verify if and how this parameter relates to the signal-to-noise ratio of a particular OCT imaging system.
The slope parameter m must have spatial units (eg μm) in order to ensure consistency.We can understand its significance by examining decorrelation saturation.Specifically, after applying the decorrelation saturation cutoff value ( 0.18 ) to Eq. ( 9) and solving for m we have 4.17 / 0.1 where 0.417 μm indicates the distance the RBCs must move in order to saturate the decorrelation measurement.
is proportional to the saturation velocity SAT v for a constant t Δ and should depend on a particular OCT imaging system.We expect that SAT x Δ , and thus the parameter m , should be proportional to a particular spatial scale related to the physical imaging parameters.It has yet to be determined, however, what the relevant spatial scale might be.One hypothesis is that , and consequently SAT v , scales with wavelength, so that the saturation velocities given in Table 2 increase by a factor of 1.24 for a 1050 nm OCT imaging system.Further experiments using different wavelengths and beam widths are required to test these hypotheses.

Model limitations
There are a number of limitations of the linear model in Eq. ( 8).First, in order to establish the linear relationship between decorrelation and velocity we had to exclude saturated points from our analysis.In practice, this establishes an upper limit on a velocity-time separation pairing for quantitative SSADA using our linear model.Furthermore, as we can see in Fig. 6, for very slow flow speeds (e.g.0.37 mm/s) the curve looks more s-shaped than linear.We expect that for slower speeds the curve will look even more s-shaped.It seems then that a more accurate model for the decorrelation to velocity relationship might be sigmoidal.Both of these models would also naturally handle the saturated data points.
Another limitation of our model is the assumption that the parameter b , related to the decorrelation of stationary scatterers, is independent of the time separation t Δ .For the time separations examined in our analysis, this parameter was indeed constant.However, we expect that the bulk decorrelation should increase with time and reach a saturation point when the time separation is very large.Furthermore, we treated the intercept as a small constant and consequently left it out of the derivation of Eq. ( 7).This parameter may be related to the decorrelation caused by shot noise, system vibrations or Brownian motion.Accounting for this parameter would cause a slight increase in the calculated saturation velocities.Lastly, we expect that the decorrelation should depend on scatterer size and possibly blood cell concentration.We used blood in our experiments to mimic in vivo measurements as closely as possible.However, we might gain some additional insight by varying scatterer size and hematocrit levels.Additional experiments are needed to test these hypotheses.

Comparison with previous work on intensity-based Doppler variance angiography
Liu and colleagues also performed similar flow phantom experiments [40,41] and found similar results for both intensity-based Doppler variance (IBDV) and Doppler phase variance.Note that they used IBDV which is similar to SSADA but does not apply spectrum splitting and uses a difference averaging procedure.Many of the results of the two works are similar, including establishing decorrelation saturation values and a linear range relating the calculated signal to velocity that depends on the time separation between measurements.However, there are a couple of important differences between the results in that work and those shown here.First, there is a significant difference regarding the dependence on Doppler angle.In order to compare the variation with Doppler angle, we first normalize our SSADA measurements by subtracting the background decorrelation of 0.08 found at zero flow.Because the IBDV background values were negligibly small compared to flow signal in [41] background subtraction was not necessary for IBDV.We compare the variation in the measurements over a Doppler angular range of approximately 18° (the largest angle tested by Liu et al.).Specifically, after background subtraction, the SSADA signal increases from 0.065 at perpendicular incidence to approximately 0.080 at a Doppler angle of 18° for the data presented in Fig. 5. Thus, our results in this work indicate that the SSADA signal increases by approximately 23% over an angular range of 18°.On the other hand, the IBDV signal in [41] increases from 80 at perpendicular incidence to approximately 150 at 18°, an 87.5% increase over the same angular range.Thus, the Doppler angle dependence of IBDV [41] was significantly higher than the angle dependence of SSADA reported here.Another important difference between this work and that found in [40] is that in [40] the authors showed a saturation velocity over 100 mm/s for a time separation of 0.02 ms, whereas our model predicts a saturation velocity of approximately 20 mm/s for that time scale.
We hypothesize that these significant differences are likely caused by the choice of the algorithms and the difference in the flowing phantom.Specifically, by creating an isotropic voxel after splitting the source spectrum, SSADA aims to reduce flow directional sensitivity over IBDV.This could explain the reduced directional dependence of SSADA over IBDV, which did not split the OCT signal spectrum and therefore had much finer axial resolution compared to transverse spot size.Additionally, intralipid solution composed of spherical particles of 0.356 μm diameter was used as the flowing phantom in [40,41].In contrast our flow phantom used whole blood where the predominant scattering particles are red blood cells that have an average diameter of 7.2 μm and a disc-like shape [50].Because Doppler variance decreases with increasing particle size [51], we expect that the saturation velocity decreases as well.Another difference between this work and the work in [40] and [41] is that a swept-source OCT system was used for the experiments in [40] and [41] but a spectral domain OCT system was used here.For the same A-line rate, the integration time in the spectrometer OCT system is longer than that in a swept-source system.The difference in integration time may also affect the saturation velocity.Additional effects of blood such as tumbling and high viscosity could cause these observed differences as well.

Clinical SSADA
The flow phantom experiments presented in this work have a number of implications for clinical SSADA.We have shown that SSADA measurements have little dependence on Doppler angle for Doppler angles less than 10 degrees.So for retinal imaging where the OCT beam is nearly perpendicular to retinal vessels, clinical SSADA may be effectively angle independent.The clinical SSADA scans previously published by our research group [34][35][36] used an interframe time separation 2.0 t Δ = ms, which is on the long side of the time scale investigated in this article.Referring back to Eq. ( 8) we see that for this time scale the saturation velocity is 0.2 SAT v = mm/s (0.3 mm/sec if adjusted for the longer wavelength of 1050 nm).Since human retinal capillary flow speeds have been estimated to be in the range of 0.4-3.0mm/s [52][53][54], this suggests that SSADA is well suited for detailed angiography down to the capillary level.However, decorrelation signal should be saturated even at the capillary level according to our phantom calibration.This does not entirely agree with the clinical retinal angiograms that we have observed, where there is a gradation of flow signals at the smallest retinal vessels as visualized by a false color scale [34,35], and this graduated flow signal increased with visual stimulation [36].This difference could be caused in part by the fact that the work in [34][35][36] used a swept-source OCT system while this work utilized a spectral domain OCT system.As described previously, the difference in integration time may affect the saturation velocity and consequently the flow signal gradation.Another explanation for the graduated decorrelation signal that we see in clinical SSADA is that there is a long sigmoidal tail above what we set a the saturation point (top of the linear region) where the decorrelation still increased with velocity, albeit at a shallow slope.We further hypothesize that the gradation could also be due to the fact that real blood capillaries are smaller than the diameter of the OCT probe beam, therefore both flow and stationary tissue existed within the same interrogation volume for SSADA, which has an expanded axial resolution due to spectral splitting.These factors may account for the proportional SSADA signal to capillary flow that is a response to either capillary diameter or velocity.Our flow phantom used a capillary tube that is much larger than the OCT beam diameter.This is one aspect in which our phantom setup differs significantly from real human capillaries.In other aspects such as the beam diameter, the use of whole blood, and the SSADA algorithm, the phantom results should simulate the clinical parameters well.
For the purpose of measuring flow velocity, a faster OCT system with a shorter interframe time scale would be better suited.Specifically, in order to bring the saturation velocity above 3.0 mm/s and thus enable capillary velocity quantification within the linear range, a time separation less than 139 t Δ = μs is suitable, according to Eq. ( 8).If our M-B mode imaging protocol calls for 200 A-lines per B-scan, then an imaging speed of 1.4 million Alines per second (1.4 MHz) is needed.Thus, megahertz OCT systems [55,56] may be useful for blood flow velocity quantification within the linear range of SSADA.

Conclusion
We performed in vitro flow phantom experiments in order to derive a linear relationship between SSADA measurements and absolute flow velocity.We hypothesized that SSADA measurements are independent of Doppler angle after splitting the spectrum to create an isotropic voxel.Contrary to our hypothesis, an angle dependence was found, but the variation due to Doppler angle was small and much reduced compared to previous work using nonsplit-spectrum intensity-based angiography [40,41].The phantom experiments established that the decorrelation signal is linearly related to velocity over a limited range, and that this range is dependent on the time scale of SSADA measurement.The saturation velocities were also proportional to the SSADA time scale.Extrapolation to OCT systems used for clinical retinal imaging with SSADA [34][35][36] indicate that the systems should be detecting flow down to the lowest range capillary velocity.Additional experiments need to be performed to understand the dependence of SSADA on signal-to-noise ratio and beam width diameter.

Financial interests
Jason Tokayer, Yali Jia and David Huang have a significant financial interest in Optovue, a company that may have a commercial interest in the results of this research and technology.These potential conflicts of interest have been reviewed and managed by Oregon Health & Science University.Al-Hafeez Dhalla has no financial interest in the subject of this article.
correction factor was used to compensate each A-line for dispersion by multiplying each interferogram by exp( subtraction, resampling to k-space, dispersion compensation and Fourier transformation we obtain complex OCT data.Sample OCT reflectance images are shown in Fig.2.

Fig. 2 .
Fig. 2. Sample intensity images.Capillary was placed on top of a piece of paper.Left panel: Bscan; Right panel: M-mode scan.

Fig. 4 .
Fig. 4. Multi-timescale SSADA M-mode image of blood flow through capillary tube.Time separation indicates the time between the two A-lines used to compute the decorrelation.

Fig. 5 .
Fig. 5. Dependence of measured SSADA decorrelation signal on Doppler angle.Spearman's correlation coefficient between the decorrelation measurements and Doppler angle is ρ = 0.87.ρ is a nonparametric measure of how well the relationship between decorrelation and Doppler angle can be described using a monotonic function.In this case the function appears to be sigmoidal.

Fig. 7 .
Fig. 7. Decorrelation vs velocity for various time separations ( t Δ ) between A-lines.Left panel: black dotted line indicates the saturation cutoff value of 0.18.Right Panel: Linear fits of decorrelation vs velocity after saturated points have been removed.
effect of this approximation is discussed in the model limitations section.
. It may also be the case that SAT x Δ