High-penetration swept source Doppler optical coherence angiography by fully numerical phase stabilization

A high-penetration swept-source optical coherence tomography (HP-SS-OCT) system based on a 1-μm short cavity laser is developed. Doppler OCT processing is applied, along with a custom-made numerical phase stabilization algorithm; this process does not require additional calibration hardware. Thus, our phase stabilization method is simple and can be employed in a variety of SS-OCT systems. The bidirectional blood flow and vasculature in the deep choroid was successfully imaged via two Doppler modes that use different time intervals for Doppler processing. En face projection image of squared power of Doppler shift is compared to ICGA, and the utility of our method is verified. © 2012 Optical Society of America OCIS codes: (170.4500) Optical coherence tomography; (170.4460) Ophthalmic optics and devices; (170.3340) Laser Doppler velocimetry (170.4470) Ophthalmology; (170.3880) Medical and biological imagin; (170.3890) Medical optics instrumentation. References and links 1. V. Patel, S. Rassam, R. Newsom, J. Wiek, and E. Kohner, “Retinal blood flow in diabetic retinopathy,” Br. Med. J. 305, 678–683 (1992). 2. E. Friedman, “A hemodynamic model of the pathogenesis of age-related macular degeneration,” Am. J. Ophthalmol. 124, 677–682 (1997). 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, 359–393 (2002). 4. 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, 529–533 (1994). 5. D. Huang, E. Swanson, C. Lin, J. Schuman, W. Stinson, W. Chang, M. Hee, T. Flotte, K. Gregory, C. Puliafito, and a. et, “Optical coherence tomography,” Science 254, 1178–1181 (1991). 6. X. J. Wang, T. E. Milner, and J. S. Nelson, “Characterization of fluid flow velocity by optical doppler tomography,” Opt. Lett. 20, 1337–1339 (1995), http://www.opticsinfobase.org/abstract.cfm?URI= ol-20-11-1337. 7. Z. Chen, T. E. Milner, D. Dave, and J. S. Nelson, “Optical doppler tomographic imaging of fluid flow velocity in highly scattering media,” Opt. Lett. 22, 64–66 (1997), http://www.opticsinfobase.org/abstract. cfm?URI=ol-22-1-64. #155068 $15.00 USD Received 21 Sep 2011; revised 9 Jan 2012; accepted 9 Jan 2012; published 23 Jan 2012 (C) 2012 OSA 30 January 2012 / Vol. 20, No. 3 / OPTICS EXPRESS 2740 8. A. F. Fercher, C. K. Hitzenberger, G. Kamp, and S. Y. El-Zaiat, “Measurement of intraocular distances by backscattering spectral interferometry,” Opt. Commun. 117, 43–48 (1995). 9. G. Häusler and M. W. Lindner, ““Coherence radar” and “spectral radar”—new tools for dermatological diagnosis,” J. Biomed. Opt. 3, 21–31 (1998). 10. 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, 457–463 (2002). 11. N. A. Nassif, B. Cense, B. H. Park, M. C. Pierce, S. H. Yun, B. E. Bouma, G. J. Tearney, T. C., Chen, and J. F. de Boer, “In vivo high-resolution video-rate spectral-domain optical coherence tomography of the human retina and optic nerve,” Opt. Express 12, 367–376 (2004), http://www.opticsinfobase.org/ oe/abstract.cfm?uri=oe-12-3-367. 12. R. A. Leitgeb, L. Schmetterer, W. Drexler, A. F. Fercher, R. J. 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, 3116–3121 (2003), http://www.opticsinfobase.org/abstract. cfm?uri=oe-11-23-3116. 13. 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, 3490–3497 (2003), http://www.opticsinfobase.org/oe/abstract.cfm?uri= oe-11-25-3490. 14. R. A. Leitgeb, L. Schmetterer, C. K. Hitzenberger, A. F. Fercher, F. Berisha, M. Wojtkowski, and T. Bajraszewski, “Real-time measurement of in vitro flow by fourier-domain color doppler optical coherence tomography,” Opt. Lett. 29, 171–173 (2004), http://www.opticsinfobase.org/abstract.cfm?URI= ol-29-2-171. 15. Y. Zhao, Z. Chen, C. Saxer, Q. Shen, S. Xiang, J. F. de Boer, and J. S. Nelson, “Doppler standard deviation imaging for clinical monitoring of in vivo human skin blood flow,” Opt. Lett. 25, 1358–1360 (2000), http: //www.opticsinfobase.org/abstract.cfm?URI=ol-25-18-1358. 16. D. Y. Kim, J. Fingler, J. S. Werner, D. M. Schwartz, S. E. Fraser, and R. J. Zawadzki, “In vivo volumetric imaging of human retinal circulation with phase-variance optical coherence tomography,” Biomed. Opt. Express 2, 1504–1513 (2011), http://www.opticsinfobase.org/abstract.cfm?uri=boe-2-6-1504. 17. 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, 6008–6025 (2008), http://www.opticsinfobase.org/abstract.cfm?URI=oe-16-9-6008. 18. 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, 114–116 (2000), http://www.opticsinfobase.org/abstract. cfm?URI=ol-25-2-114. 19. S. Makita, Y. Hong, M. Yamanari, T. Yatagai, and Y. Yasuno, “Optical coherence angiography,” Opt. Express 14, 7821–7840 (2006), http://www.opticsinfobase.org/abstract.cfm?URI=oe-14-17-7821. 20. Y. Ikuno and Y. Tano, “Retinal and choroidal biometry in highly myopic eyes with spectral-domain optical coherence tomography,” Invest. Ophthalmol. Vis. Sci. 50, 3876–3880 (2009). 21. Y. Imamura, T. Fujiwara, R. Margolis, and R. F. Spaide, “Enhanced depth imaging optical coherence tomography of the choroid in central serous chorioretinopathy,” Retina 29, 1469–1473 (2009). 22. U. Chakravarthy, J. Evans, and P. J. Rosenfeld, “Age related macular degeneration,” Br. Med. J. 340 (2010). 23. A. Unterhuber, B. Považay, B. Hermann, H. Sattmann, A. Chavez-Pirson, and W. Drexler, “In vivo retinal optical coherence tomography at 1040 nm enhanced penetration into the choroid,” Opt. Express 13, 3252–3258 (2005), http://www.opticsinfobase.org/oe/abstract.cfm?uri=oe-13-9-3252. 24. D. M. de Bruin, D. L. Burnes, J. Loewenstein, Y. Chen, S. Chang, T. C. Chen, D. D. Esmaili, and J. F. de Boer, “In vivo three-dimensional imaging of neovascular age-related macular degeneration using optical frequency domain imaging at 1050 nm,” Invest. Ophthalmol. Vis. Sci. 49, 4545–4552 (2008). 25. 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 1-um swept source optical coherence tomography andscattering optical coherence angiography,” Opt. Express 15, 6121–6139 (2007), http://www.opticsinfobase.org/oe/abstract. cfm?uri=oe-15-10-6121. 26. B. Povazay, B. Hermann, B. Hofer, V. Kají, E. Simpson, T. Bridgford, and W. Drexler, “Wide-field optical coherence tomography of the choroid in vivo,” Invest. Ophthalmol. Vis. Sci. 50, 1856–1863 (2009). 27. V. J. Srinivasan, D. C. Adler, Y. Chen, I. Gorczynska, R. Huber, J. S. Duker, J. S. Schuman, and J. G. Fujimoto, “Ultrahigh-speed optical coherence tomography for three-dimensional and en face imaging of the retina and optic nerve head,” Invest. Ophthalmol. Vis. Sci. 49, 5103–5110 (2008). 28. Y. Yasuno, M. Miura, K. Kawana, S. Makita, M. Sato, F. Okamoto, M. Yamanari, T. Iwasaki, T. Yatagai, and T. Oshika, “Visualization of sub-retinal pigment epithelium morphologies of exudative macular diseases by highpenetration optical coherence tomography,” Invest. Ophthalmol. Vis. Sci. 50, 405–413 (2009). 29. S. H. Yun, G. Tearney, J. de Boer, and B. Bouma, “Motion artifacts in optical coherence tomography with frequency-domain ranging,” Opt. Express 12, 2977–2998 (2004), http://www.opticsinfobase.org/ #155068 $15.00 USD Received 21 Sep 2011; revised 9 Jan 2012; accepted 9 Jan 2012; published 23 Jan 2012 (C) 2012 OSA 30 January 2012 / Vol. 20, No. 3 / OPTICS EXPRESS 2741 abstract.cfm?URI=oe-12-13-2977. 30. B. Potsaid, B. Baumann, D. Huang, S. Barry, A. E. Cable, J. S. Schuman, J. S. Duker, and J. G. Fujimoto, “Ultrahigh speed 1050nm swept source/fourier domain OCT retinal and anterior segment imaging at 100,000 to 400,000 axial scans per second,” Opt. Express 18, 20029–20048 (2010), http://www.opticsinfobase. org/abstract.cfm?URI=oe-18-19-20029. 31. W. Wieser, B. R. Biedermann, T. Klein, C. M. Eigenwillig, and R. Huber, “Multi-megahertz OCT: High quality 3D imaging at 20 million a-scans and 4.5 Gvoxels per second,” Opt. Express 18, 14685–14704 (2010), http: //www.opticsinfobase.org/abstract.cfm?URI=oe-18-14-14685. 32. T. Klein, W. Wieser, C. M. Eigenwillig, B. R. Biedermann, and R. Huber, “Megahertz OCT for ultrawide-field retinal imaging with a 1050nm fourier domain mode-locked laser,” Opt. Express 19, 3044–3062 (2011), http: //www.opticsinfobase.org/abstract.cfm?URI=oe-19-4-3044. 33. B. Vakoc, S. Yun, J. de Boer, G. Tearney, and B. Bouma, “Phase-resolved optical frequency domain imaging,” Opt. Express 13, 5483–5493 (2005), http://www.opticsinfobase.org/oe/abstract.cfm?uri= oe-13-14-5483. 34. J. Zhang and Z. Chen, “In vivo blood flow imaging by a swept laser source based fourier domain optical doppler tomography,” Opt. Express 13, 7449–7457 (2005), http://www.opticsinfobase.org/abstract. cfm?uri=oe-13-19-7449. 35. 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, 1539–1552 (2011), http://www.opticsinfobase.org/abstract.cfm? URI=boe-2-6-1539. 36. B. Braa


Introduction
Human tissues, including those in the eye, are fed and they excrete waste through blood circulation, and thus, the abnormality of the vasculature is an important indicator of malfunction in the eye [1][2][3]. Therefore, the examination of ocular vessels is important in ophthalmic diagnosis. Fluorescein angiography (FA) and indocyanine green angiography (ICGA) are current angiographic methods used in ophthalmology. However, they require dye injection to enhance the contrast of the vessels, they are uncomfortable for the patients, and are sometimes associated with severe adverse reactions [4].
OCT is a noninvasive imaging modality that produces cross-sectional images at a micrometer scale [5]. Doppler OCT is an extension of OCT that is capable of detecting the Doppler shift of a probe beam to provide blood flow contrast with a depth-resolved tissue structure [6,7]. On the other hand, Fourier-domain (FD) OCT, which is a second-generation OCT, enables rapid, threedimensional volumetric imaging [8][9][10][11]. In FD-OCT, the Doppler frequency shift is calculated from the phase difference between adjacent A-lines; for this phase difference calculation [12][13][14], high phase stability is required. Although there are several other methods of increasing the contrast of vessels using the Doppler effect, such as Doppler standard deviation imaging [15], phase-variance OCT [16], and joint spectral and time domain OCT [17], phase difference methods are the most common because of their simplicity of implementation [12,13,18].
To supplement current angiographic methods, optical coherence angiography (OCA) was developed; it is based on the 840-nm spectral domain (SD) Doppler OCT [19]. Recently, there were reports that several ocular diseases, myopia, central serous chorioretinopathy (CSC), and age-related macular degeneration (AMD), are related to choroidal properties [20][21][22]. Although the 800-nm band SD-OCT-based OCA was able to measure retinal vessels, its ability to examine choroid was limited. This is mainly because the 800-nm-band probing beam is strongly absorbed by both the melanin in retinal pigment epithelium (RPE) and the choroid, and thus, penetration into the deep choroid is low. To obtain high penetration into the choroid, in vivo retinal OCT at 1-μm was demonstrated [23-28]. The low absorption of melanin at a wavelength of around 1 μm enables high penetration into the choroid.
Recently, swept-source (SS) OCT became a common option to implement 1-μm OCT because of its advantages over SD-OCT, including robustness to motion [29], a long measurement range in depth due to short instantaneous line-width [30], k-linear sampling, compactness [30], and high imaging speed [31,32].
However, SS-OCT typically has jitter in synchronization between the wavelength sweep and data acquisition; this jitter causes random spectral shifts among interference spectra and results in low phase stability. This low phase stability is a critical problem for phase-sensitive OCT modalities including Doppler OCT. To resolve this issue, previously demonstrated SS-Doppler-OCT commonly utilized a stationary mirror or a glass plate in the sample arm [33][34][35]. In this scheme, the unstable phase is measured by monitoring the phase of the mirror or the glass plate and is then corrected numerically. Although this method successfully enhances the phase stability, it requires additional hardware in the sample arm, and thus, it makes the system more complicated and limits the depth measurement range. Recently, Braaf et al. demonstrated phase stabilized SS-OCT which stabilizes the phase by using a separate Mach-Zehnder interferometer from an OCT interferometer [36]. Although this system demonstrated excellent phase stability and did not sacrifice depth imaging range, it still requires additional hardware complexity.
In this study, we demonstrate a high-penetration Doppler OCT and OCA based on high-speed 1-μm SS-OCT with a custom-made, fully numerical OCT phase stabilization method without needing any additional calibration hardware. Bidirectional choroidal blood flow and vasculature imaging are then demonstrated. This method does not require any additional calibration hardware, and thus, can be used with any SS-OCT system with a simple software update.

Swept-source optical coherence tomography
High-penetration Doppler OCA (HP-DOCA) is based on an SS-OCT system with a fiber-based Michelson interferometer. A schematic of the SS-OCT is shown in Fig. 1. A commercially available short cavity laser (Axsun Technologies Inc., MA) with a center wavelength of 1060 nm, a spectral scanning width of 120 nm, a -3-dB bandwidth of 100 nm, and a scanning speed of 100,000 scans/s was utilized. The laser beam passes through a polarization controller and enters an polarization independent circulator (Agiltron, MA) with an operation wavelength of 1060 ± 30 nm. A filter coupler (Advanced Fiber Resources Ltd., Hong Kong) passes 50% of the light to the reference arm and reflects 50% into the sample arm. The reference beam power is adjusted using an iris diaphragm. A custom retinal scanner based on a commercial OCT device (3D-OCT-1000, TOPCON Corp., Tokyo, Japan) is connected to the sample arm. The galvano-  metric scanner steering the probe beam is controlled by a function generator board (PCI-6713, National Instruments, TX). The probing beam power was set as 1.85 mW, which is lower than the safe exposure limit of ANSI standard. The reflected beams from both arms are combined at the coupler and detected by a balanced photodetector with a bandwidth of DC to 350 MHz (PDB130C, Thorlabs Inc.). To avoid high-frequency noise, the detector output signal is passed through a low-pass filter with a cut-off frequency of 170 MHz (7th order type-I Chebyshev filter, R & K Corp., Japan). The interference signal is digitized by a data acquisition (DAQ) board with 12-bit quantum resolution and a maximum sampling rate of 500 MS/s (ATS9350, AlazarTech, QC, Canada).
Data sampling was performed using an optically generated clock signal that is equi-spacing in wavenumber, namely, a k-clock. With the light source utilized in this study, the k-clock becomes unstable during the backward sweep, although only the forward sweeps are utilized for imaging. This unstable k-clock electrically interferes with the triggering of data acquisition. To supply a stable clock signal to the DAQ board during the backward sweep and avoid this interference, a custom-made dummy clock generator is employed. The dummy clock generator switches the clock signal supplied to the DAQ board between the k-clock and a dummy clock signal with a constant frequency (240 MHz) during forward and backward sweeps, respectively. It is noteworthy that the recent versions of 1-um swept-source provided by AXSUN Technologies possess a similar dummy clock generation function, and thus, a custom-made dummy clock generator is not needed. Owing to this k-linear sampling, no rescaling process is required.
Synchronization between data acquisition and the galvanometric scanner is achieved using the sweep trigger signal of the light source, which is utilized for triggering the data acquisition of the DAQ for a single sweep and utilized as an update clock for the waveform generation for galvanometric scanning.
The short cavity laser possesses a very narrow instantaneous line width; this short line width provides a very low signal roll-off. The signal roll-off was measured to be -2.5 dB over 3.0 mm. The sensitivity decay was measured to be -7.4 dB over 3.0 mm in air. It should be noted that the sensitivity decay is higher than the signal roll-off. This is due to the noise floor increase in the deeper region. The sensitivity was measured to be 99.1 dB in which the recoupling loss at the fiber tip of the sample arm is compensated. The expected shot-noise-limit sensitivity was 106.1 dB. The optical loss of the system was measured to be 4 dB; this may partially account for the departure of the measured sensitivity from the shot-noise-limit sensitivity. The depth resolution defined by the -6-dB width was measured to be 11.0 μm in air, which corresponds to the resolution of 8.0 μm in tissue (n = 1.38).

Phase-1 numerical phase stabilization: rough spectral shift correction
SS-OCT is commonly associated with jitter in synchronization between the wavelength sweep and data acquisition. This jitter results in random shifts in spectral sampling, and thus, results in jitter in the phase of the OCT signal. As the Doppler frequency obtained by the phase sensitive Doppler OCT algorithm is defined as the phase difference between the two A-scans, this jitter causes a significant error in the Doppler OCT measurement. In this section, we describe a fully numerical method of cancelling this jitter and for stabilizing the phase of the OCT; we refer to this method as "phase-1 stabilization" in this paper.
The phase error under consideration is caused by a relative shift between two sampled spectra S 1 ( j) and S 2 ( j), where j is the index of the sampling point of the spectrum and is linear with wavenumber k. We assume that S 2 ( j) suffers from spectral shift by jitter and becomes where β is the spectral shift in the sampling points and * denotes convolution. According to the shift theorem of Fourier transform, the discrete Fourier transform of the shifted spectrum S 2 ( j) becomes the discrete Fourier transform of S 2 ( j) with an erroneous phase: where represents the discrete Fourier transform, N is the number of sampling points in a single spectrum, and ζ = 0, 1, · · ·N − 1 is the index of the discrete Fourier pair of j that also represents the depth position of the OCT signal. If the two spectra, and thus, the two A-lines yielded from these spectra were fully correlated except for the spectral shift, i.e., S 1 ( j) = S 2 ( j), the Doppler shift obtained from these two spectra would present this erroneous phase as ∠ S 2 (ζ ) S 1 * (ζ ) = −2πζ β /N, where the superscript * denotes the complex conjugate. Under this condition, the erroneous phase caused by the jitter could be determined and numerically completely canceled. However, in practice, S 1 ( j) and S 2 ( j) are not fully correlated because of the structural difference between the two A-scans and the Doppler shift caused by the flow and bulk motion of a sample. In our strategy, we obtain and correct this erroneous phase by excluding these disturbances by cropping the perfectly correlated component of the spectra, i.e., the spectrum of the reference beam.
In the practical OCT measurements, the spectral interferometric signal S 1 ( j) is expressed as where E r ( j) and E p1 ( j) are optical fields of the reference beam and probe beam of this A-scan. In a standard configuration, the reference arm length is adjusted so that the OCT signals, i.e., the third and the fourth terms in this equation, do not overlap with the zero-delay position (where the first and the second terms are located). Thus, we can filter out the third and fourth terms. In addition, in the standard OCT configuration, the reference power of the photodetector is significantly higher than the probe power; After performing this filtering and making this assumption, the spectrum becomes S 1 ( j) |E r ( j)| 2 ; its Fourier transform becomes Note that, in practical SS-OCT configuration, |E r ( j)| 2 to be measured is not a reference spectrum in a precise sense because of the common mode rejection capability of balanced photodetection. However, this effect is negligible as we discuss in Section 4.1.
Similarly, the shifted spectrum of S 2 ( j) in practical OCT measurement is expressed as where E p2 ( j) is the optical field of the probe beam of this A-scan. A similar filtering and assumption for S 1 ( j) provides the cleaned-up spectrum of It is evident that Eqs. (4) and (6) are identical except for the erroneous phase slope component caused by the spectral shift. This erroneous phase is easily obtained using the following operation and thus, is easily canceled numerically. where and is equivalent to the OCT signal intensity.
In practical calculations with measured spectra, the phase obtained from Eq. (7) does not always have a perfectly linear slope. Thus, we determine the phase slope by weighted leastsquare fitting which determines the value of β that minimizese the following error energy.
where ϕ(ζ ) is the differential phase defined as ∠ S 2 (ζ ) S * 1 (ζ ), and S 2 (ζ ) and S * 1 (ζ ) are derived from measured spectra S 1 ( j) and S 2 ( j). I(ζ ) is set to zero if it is smaller than the noise level.
A detailed overview of the numerical phase stabilization process used in our particular implementation is as follows. First, the DC offset for each spectrum is removed by subtracting the mean of each spectrum so as to set the balanced voltage to zero. Discrete Fourier transform is performed after applying a Gaussian window with a standard deviation of 0.2N where N = 1536. Offset cancelation is important, otherwise the interaction between the offset and the window function degrades the accuracy of phase-1 stabilization. Namely, the offset reformed by the window function numerically mimics a pseudo spectrum with the shape of the window function, which has no spectral shift and overlapped with the real spectrum, and hence it disturbers the phase-1 stabilization. It is noteworthy that numerical dispersion compensation is not applied because the signal of interest is formed only from the reference beam, and thus, it is not affected by dispersion. The ±30 pixels close to the zero-delay are utilized for the analysis. This ±30-pixel width corresponds 304 μm in air and roughly 28 times the measured depth resolution defined as the 6-dB width. The phase difference between an A-line under shift correction and a reference A-line (typically the first A-line in an OCT volume) is obtained by multiplying the complex conjugate of the reference A-line by the A-line under correction. Least square fitting, described above, is applied to obtain the phase slope, which is proportional to the spectral shift. Finally, cancelation of the spectral shift is performed by multiplying the complex conjugate of the estimated phase slope to the entire ζ -range of the Fourier-transformed spectrum and then taking the inverse Fourier transform. Standard Fourier domain OCT (FD-OCT) processing (including complex median subtraction to eliminate fixed pattern noise [37] and automatic numerical dispersion compensation [25]) is then applied to obtain a phase-stabilized OCT image.
It should be noted that, phase-1 stabilization uses a signal generated by a reference beam which occupies the region close to the zero delay. And hence OCT signal should not be overlapped with this region and hence this region cannot be utilized for OCT imaging despite of its high sensitivity. This might be a drawback of phase-1 stabilization. However, in the clinical practice, imaging a retina at a region close to zero delay is generally discouraged to prevent the image folding due to the eye motion. And hence the negative effect of this drawback would be minimal.

Phase-2 numerical phase stabilization: Correction of residual spectral shift and bulk motion
To obtain the Doppler OCT image, the phase difference between adjacent A-lines of the phasestabilized OCT image is calculated by multiplying an A-scan by the complex conjugate of the other A-scan. As described in Section 2.2, spectral shift cancelation (phase-1 stabilization) is performed using tens of data points, which is close to the zero delay. Because data points close to the zero delay are relatively insensitive to the phase error than the data points far from the zero delay, a certain level of residual fine spectral shift will still exist. In addition, the Doppler frequency measured from an in vivo subject, such as a human eye, is sensitive not only to the flow of interest but also to the bulk motion of the subject. Although bulk motion is a true signal in the broadest sense of term, it frequently disrupts flow measurement and should therefore be canceled. In this section, we describe the phase-2 phase stabilization algorithm, which eliminates the residual phase errors caused by residual spectral shift and bulk motion.
By accounting for these errors, the Doppler phase shift obtained from a sample is described as where τ is the time interval between two A-scans under Doppler calculation, λ c is the center wavelength, n is the refractive index of the sample, v z (ζ ) is an axial velocity of the flow of interest, v b is the axial velocity component of the sample bulk motion, and β is the residual spectral shift. It should be noted that v z (ζ ) is a function of depth (ζ ), whereas v b is a constant; therefore, the second term of this equation is a depth-independent phase offset. The third term is generated by the residual spectral shift and is a linear function of depth whose slope is in proportion to the residual spectral shift β . It is noteworthy that if β is very large, phase wrapping occurs and the third term can no longer be a linear function of depth. In our algorithm, β is warranted to be sufficiently small owing to phase-1 stabilization. In our phase-2 stabilization algorithm, weighted linear fitting is applied to each depthresolved Doppler signal Δϕ(ζ ). If the effect of the flow of interest is negligible (in other words, if we only consider the Doppler signal of a static part of the sample), the slope and intercept of the linear fitting line are respectively proportional to the residual spectral shift and the bulk motion as evident in Eq. (9). Therefore, these two errors can be canceled by subtracting this fitting line from the phase. However, in practical measurement, the flow of interest is not negligible. In the phase-2 stabilization algorithm, this disturbance from the flow of interest is eliminated by iterative fitting described below.
In this iteration, weighted linear fitting was applied to the measured Doppler signal, which determines the slope and intersept so as to minimize the following error energy: where Δϕ(ζ ) is the Doppler signal obtained from a measurement, a m and b m are the slope and the intersept determined by the i-th iteration, respectively, and W m (ζ ) is the weight of the fitting. For the first iteration (m = 0), W 0 (ζ ) is defined as where I(ζ ) is the OCT signal intensity of the corresponding A-line, and ε 2 is the noise level.
In this method, we iteratively exclude the part with flow from the fitting by relying on the assumption that the Doppler shift caused by the real flow is eccentric with respect to the phase slope defined by the second and third term of Eq. (9). To exclude the part that satisfies this assumption, the weight is updated as where m = 1, 2, · · ·. It should be noted that, according to this equation, the weight is not updated if the fitting error is less than or equal to a particular threshold (π/(2m)), otherwise, the weight is set to zero. The threshold becomes smaller as the iteration progresses. After a predefined number of iterations, the phase slope obtained is subtracted from the measured Doppler shift to cancel the effects of residual spectral shift and bulk motion. In our particular implementation, the maximum iteration was set to seven (m = 6).
To further suppress the residual phase error, a Kasai auto-correlation filter [38, 39] is applied with a three-dimensional kernel of 2 pixels (depth) × 4 A-lines × 2 B-scans for some cases.

Fast and slow Doppler modes
Because flow velocity is calculated from phase difference, the measurable minimum axial velocity is limited by phase noise σ ϕ [40]. Further, the measurable maximum axial velocity is limited by the 2π-ambiguity of the phase, i.e., the measured phase value is wrapped onto the range from to −π to π. Finally, the measurable velocity range is given by where τ is the time interval between the two A-scans utilized for the Doppler calculation. It should be noted that the velocity range is downward-shifted as τ increases. In other words, a larger value of τ provides slower Doppler velocity signal and vice versa.
In our device, we implemented two scanning modes to achieve two velocity ranges; so-called fast and slow Doppler modes. The fast Doppler mode, also referred to as the intra-B-scan mode, calculates the phase difference between adjacent A-lines in a B-scan as same as conventional Doppler OCT. In our implementation, τ = 10 μs. In this configuration, the theoretically estimated best velocity range is from 1.75 mm/s to 19.2 mm/s, where the phase noise is assumed to be limited by the signal-to-noise ratio (SNR) and structural decorrelation by the transverse scanning. The SNR was assumed to be 25 dB, and the structural decorrelation was estimated from the separation of A-lines at one-tenth the optical transverse resolution.
The slow Doppler mode, also referred to as the inter-B-scan mode, calculates the phase difference between adjacent B-scans at each A-line [41]. In our implementation, τ = 3.2 ms. The velocity range at this configuration is theoretically estimated from 5.2 μm/s (at 25-dB SNR) to 60 μm/s, where the phase noise is again assumed to be limited by the SNR and structural decorrelation. With this mode, almost all vessels are measured owing to the slow minimum measurable speed. However, the reduced maximum speed results in multiple phase-wrappings. Two Doppler modes are switched by changing the scanning protocol. The typical protocols for the two Doppler modes are 2048 A-lines (horizontal) × 256 A-lines (vertical) for the fast Doppler mode and 256 A-lines (horizontal) × 2048 A-lines (vertical) for the slow Doppler mode. Both scans of the two Doppler modes typically cover 6 mm × 6 mm area. To make the phase noise for both Doppler modes identical, the separation between A-lines utilized for the Doppler processing is set identically for both modes. The separation is about one tenth of the spot size of the probing beam: 2.9 μm for a spot size of 30 μm in our particular implementation.
For the visualization of Doppler images, two visualization modes are utilized. One is a bidirectional Doppler visualization mode in which a white-to-black colormap ranged from −π rad to π rad or a blue-to-red colormap ranged from −0.37π rad to 0.37π rad is assigned to Δϕ. The positive value of Δϕ represents a positive flow direction, which represents the flow direction from posterior to the anterior in a posterior eye. The other is a power Doppler visualization mode in which a black-to-white colormap or a gradual-red colormap (invisible for 0 rad 2 to red for π 2 rad 2 ) is assigned.

Phase stability
To quantitatively examine the enhancement of phase stability obtained from the phase-1 stabilization algorithm, a stationary mirror was used as a sample. In this measurement, lateral scanning was not performed and the mirror was measured at several distances from the zero delay. The standard deviation of the phase differences between adjacent A-lines (σ ϕ ) is measured at the mirror position. The standard deviations were obtained with and without phase-1 stabilization; 1024 A-lines were utilized for each calculation. As shown in Fig. 2, the phase stability is significantly enhanced, especially in the deep region.
As the phase-2 stabilization algorithm is designed to work with a practical sample that involves bulk motion and a flow of interest, it is not possible to quantitatively evaluate this algorithm using the mirror sample. To semi-quantitatively evaluate this algorithm, we measured an in vivo human retina and evaluated the histogram of a section of the retina that does not contain visible blood vessels. Figure 3 shows bidirectional Doppler phase shift images obtained from an in vivo human eye. Each pane represents the bidirectional Doppler signal without phase stabilization ( Fig. 3(a)), with phase-1 stabilization (Fig. 3(b)), with phase-1 and phase-2 stabilizations (Fig. 3(c)), and with phase-1 and phase-2 stabilizations and Kasai filter (Fig. 3(d)). The region indicated by the small yellow box is utilized for the phase stability evaluation. The insets indicated by the large yellow boxes are the magnified images of the small yellow box. Figure 4 shows histograms obtained from the regions indicated by yellow boxes in Fig. 3

. Figures 4(a)-4(d) correspond to Figs. 3(a)-3(d).
Abrupt spectral shift caused by erroneous synchronization between the starting trigger of the wavelength sweep and data acquisition resulted in several vertical line artifacts, as shown in Fig. 3(a) (arrows). This produced the side lobes in the corresponding histogram of Fig. 4(a)   a large standard deviation of Δϕ of 43.6 • . It should be noted that this standard deviation would be significantly larger than that measured using a static mirror possibly because of the lower SNR and structural decorrelation between A-lines, and thus, cannot be directly compared with the value presented in Fig. 2. By applying phase-1 stabilization, the vertical line artifacts were dramatically reduced, as shown in Fig. 3(b). This reduction corresponded with the disappearance of side lobes in Fig. 4(b). The standard deviation decreased to a value as small as 38.0 • . Despite the reduction of standard deviation, the constant phase offset induced by bulk motion still exists, as shown in the histogram. This phase offset was successfully eliminated by phase-2 stabilization, as shown in Fig. 4(c). It was also found that the histogram was further narrowed and that the standard deviation was decreased to 35.7 • . The Kasai filter further reduced the phase noise, as shown in Fig. 3(d), and further narrowed the histogram to the standard deviation of 9.8 • , as shown in Fig. 4(d).

Bidirectional blood flow imaging by fast Doppler mode
To demonstrate bidirectional in vivo flow imaging of a retina, we measured a healthy retina of a 31-year-old Asian male. 6 mm × 6 mm areas centered at a macula and an optic nerve head (ONH) were scanned with 2048 × 256 A-lines in the fast (intra-B-scan) Doppler mode. A single volumetric scan took 6.6 s. Figure 5 shows the cross-sectional OCT intensity and Doppler OCT images. In the macular region, several choroidal blood flows are visible; they are indicated by arrows in Fig. 5(b). In the ONH region, both the choroidal and the retinal vessels are visible; they are indicated by arrows in Fig. 5(d). The fly-through movies are available for  OCT intensity is rendered as green; the bidirectional blood flow is rendered using a red-blue colormap. Both in the macular and the ONH regions (Figs. 6(a) and 6(b)), several choroidal vessels are visualized. In the ONH region ( Fig. 6(b)), large retinal vessels are clearly seen, while the visibility of retinal vessels in macula (Figs. 6(a)) is relatively low. This may because the measurable minimum velocity of the fast Doppler mode is too high to visualize fine retinal vessels in the macular region. We will demonstrate the visualization of fine retinal vessels in the slow Doppler mode in Section 3.3.
To evaluate the performance of high-penetration bidirectional blood flow imaging of fast Doppler mode, we measured 12 eyes of 6 Asian subjects without marked posterior disorders. All subjects were male and their age was 29.6 ± 3.7 years (mean ± standard deviation); their ages ranged from 25 to 36 years. For this study, 3 mm × 3 mm areas of maculae and ONHs were scanned with 1024 × 256 A-lines. The scan of a single volume took 3.3 s. Some representative images are shown in Fig. 7.
In the ONH area, bidirectional blood flow was visible in the retina and choroid in all 12 eyes. Among them, blood flow beneath the lamina cribrosa was visible in 1 eye, as indicated by a yellow circle in Fig. 7(c) (case-1). In the macular area, blood flow was clearly observed in choroid in 11 out of the 12 eyes, as shown in Fig. 7(f) (yellow circles, case-2). Eight of the 12 eyes were myopic, which is defined as the spherical equivalent refractive error of more than -3 diopters; 6 of the 8 myopic eyes were appeared with myopic conus. As exemplified by Fig. 7(i) (yellow circle, case-3), 5 of the 6 eyes with myopic conus showed blood flow signals beneath the conus. The red circles in Figs. 7(b), 7(e), and 7(h) indicate the same locations as the yellow circles in Fig. 7(c), 7(f), and 7(i). It should be noted that here the vessels/flow were classified as "visible" as long as a portion of the vessel was observed in the bidirectional Doppler image.

Vasculature imaging in slow Doppler mode
The same subject as described as the first subject in Section 3.2 was examined in the slow Doppler mode. 6 mm × 6 mm areas of ONH and macula were scanned with 256 × 2048 Alines in 6.6 s. Figure 8 shows the cross-sectional OCT intensity (Fig. 8(a)) and power of the Doppler phase shift (Fig. 8(b)) images in the macula area, where the green arrows indicate the same location in these two figures. Instead of losing the bidirectionality of the Doppler phase shift because of the reduced measurable velocity range and its resulting multiple phasewrapping, many Doppler shift signals are observed. The retinal vessels are visible in both the OCT and the Doppler OCT images, as indicated by red arrows. Although it is not easy to    recognize individual vessel structures, densely visualized Doppler phase shift signals under the RPE are consistent with the anatomical property of the choroid. In Fig. 8(b), the power Doppler signal in the choroid shows a random appearance.
To provide a more comprehensive clinical understanding of vasculature, en face projections of the power Doppler signal, also referred to as OCA, were created, as shown in Fig. 9. For comparison, the power Doppler projections were created both in the fast Doppler mode ( Fig.  9(a)) and in the slow Doppler mode ( Fig. 9(b)). In the case of the fast Doppler mode ( Fig. 9(a)), the density of visualized vessels is small because of the relatively high minimum measurable velocity in the slow Doppler mode. In the case of the slow Doppler mode ( Fig. 9(a)    Owing to the depth resolution of OCT and automatic segmentation of RPE [19], we created two separate OCAs from a single slow Doppler volume; one is for retinal angiogram and the other is for choroidal angiogram, as shown in Fig. 10. It is noteworthy that this clear separation of depth layers can be performed neither by fluorescein angiography (FA) nor ICGA. The horizontal line artifacts indicated by red and yellow arrows are caused by imperfections in axial shift motion cancelation and saccadic motion, respectively. The details of the line artifacts are discussed in Section 4.4.

Clinical case: Polypoidal choroidal vasculopathy
To evaluate the clinical performance of this Doppler OCT system, one eye with polypoidal choroidal vasculopathy (PCV) was examined. The subject was a 71-year-old Japanese male. He was diagnosed with PCV after a set of standard clinical examination at the Ibaraki Medical Center of Tokyo Medical University and then transferred to the University of Tsukuba to be scanned by Doppler OCT. The protocol for examinations was approved by the institutional review board of the Tokyo Medical University and written informed consent was obtained before the examination. A 6-mm × 6-mm area of macula was scanned with 2048 × 256 A-lines for the fast Doppler mode and 256 × 2048 A-lines for the slow Doppler mode. The acquisition time is 6.6 s for each mode.
Several OCT and OCA images for this case are shown in Fig. 11; a comparison between ICGA and slow Doppler OCA images is shown in Fig. 12. A bright hyper-fluorescent spot is observed in the ICGA image; it is indicated by a yellow arrow. The corresponding region is  (Fig. 11(b)) and OCA image (Fig. 11(d)), respectively. In the OCT intensity image (11(a)), a moderately hyper scattering region is observed between the detached RPE and Bruch's membrane, as indicated by a yellow arrow. The appearance of this area is consistent with previous studies [28,42]; this hyper scattering tissue is hypothesized to be an abnormal vessel [42]. By comparing the intensity and Doppler B-scans, we found that the strong flow in the bidirectional Doppler Bscan corresponds to the hyper scattering region beneath the RPE. This appearance is consistent with a previous study [43], and may support the abovementioned hypothesis.
In the ICGA image ( Fig. 12(a)), abnormal choroidal vascular networks characterized by thick choroidal vessels with unclear boundaries were observed at the region around the yellow arrow. A similar appearance is observed in the slow Doppler OCA images (Figs. 11(f) and 12(b)). In the corresponding composite B-scan ( Fig. 11(g)), in which power Doppler signal is overlaid onto the OCT intensity image with a gradual-red colormap, the power Doppler signal is observed between the detached RPE and Bruch's membrane.

Effects of common mode rejection by dual-balanced photodetector
As described in Section 2.2, the phase-1 stabilization algorithm relies on the OCT signal that is generated only by the reference beam. However, this signal is a common mode for dualbalanced photo-detection, and thus, cannot be detected if the dual-balanced photo-detection is perfectly working. In this section, we show that a standard dual-balanced photodetector and standard power configuration of OCT provide sufficient signal intensity to allow the use of the phase-1 stabilization algorithm owing to the inevitable imperfection of the photodetector.
In a standard OCT configuration, the reference power at the detector is much higher than the probe power. In our configuration, for a probe power at the sample of 1.85 mW and a well known retinal reflectivity of -60 dB, the probe power at the detector is estimated to be 1.85 nW. In contrast, the reference power at the detector was configured to be approximately 230 μW. Although a huge portion of the reference power is rejected by the common-mode rejection nature of dual-balanced detection, an considerable portion remains. The typical common-mode rejection ratio of our balanced photodetector is -40 dB; therefore, the residual reference power should be 23 nW. This is still 12 times higher than that of the probe beam. Therefore, even with the dual-balanced photo-detection scheme, the signal intensity would be sufficiently high for the phase-1 stabilization.
The following analysis would provide more quantitative measure of the required signal intensity for phase-1 stabilization. The phase-1 stabilization relies on linear fitting of the phase around the zero-delay point. In our implementation, the fitting area is ±30 pixels around the zero delay. Since the phase is an odd function of a depth, this fitting is equivalent to the fitting of the area from 0th to 30th pixel with a constriction of null-intercept. Under this condition, a least square fitting provide a slope of the fitting line as ∑ i ζ i ϕ i / ∑ i ζ 2 i , where ζ i and ϕ i are the depth position and the phase value of the i-th data point to be fitted, respectively, and i = 0, 1, · · · , 30. Here we have assumed, for simplicity, the signal is uniformly distributed along the fitting area, and hence the weight of the fitting is constant.
By using this fitting line, the erroneous phase occurred by a spectral jitter at the depth of z is expected to be ϕ = ∑ i ζ i ϕ i / ∑ i ζ 2 i ζ . By partially differentiating this equation by ϕ i , we have the following equation which relates the tiny variations of ϕ and ϕ i as where dϕ and dϕ i are the tiny variations of ϕ and ϕ i . Here dϕ i is interpreted as an error in the phase measurement due to a stochastic phase noise of OCT and dϕ is a resulting error in the phase estimation. Although the stochastic phase error cannot be deterministically predicted, its statistical property can be predicted as , where σ i and SNR i are the standard deviation of the stochastic phase error and SNR at z i . Since the signal intensity is assumed to be uniformly distributed within the fitting region, σ 0 = σ 1 = · · · = σ 30 ≡ σ . In our specific configuration, σ is expected to be 0.025 rad.
By applying error propagation analysis to Eq. (14), it is found that the phase estimation at z might have a phase error with standard deviation of By substituting our experimental parameters, this equation gives a theoretical prediction of phase estimation error at each depth as shown in Fig. 2 (green line). This theoretical prediction is well agreed with the experimental results. According this prediction, the phase error at the deepest perimeter of the measurement range is predicted to be σ d = 0.19 rad. Since the condition which warrants a correct operation of phase-2 stabilization, i.e. non-existence of phase wrapping, is σ d < π, the phase-1 stabilization surely can provide a sufficient phase stabilization.

The role of phase-1 stabilization
The phase-2 stabilization aims at canceling both a Doppler phase slope occurred by spectral shift and Doppler phase offset occurred by bulk motion of the sample, while phase-1 stabilization is designed only to cancel the Doppler phase slope. Because of this inclusive property of phase-2 stabilization, phase-1 stabilization is not always required. However, phase-2 stabilization can work only if no phase wrapping exists in a static tissue which phase-2 stabilization uses to monitor the phase. This condition of no phase wrapping is not always warranted. In contrast, phase-1 stabilization is robust for this wrapping issue although its accuracy is less excellent than phase-2. This is because phase-1 stabilization relies on phase close to the zero delay where phase wrapping is unlikely happened. And hence, the purpose of phase-1 is to cancel the phase slope only with low accuracy but to surely eliminate phase wrapping in the region of the static tissue. This elimination assures successful application of the subsequent phase-2 stabilization.
An example showing this role of phase-1 stabilization is shown in Fig. 13. Figure 13(a) shows an example of a bidirectional Doppler image which is stabilized only by phase-2 stabilization. Evident artifacts, vertical lines, exist as indicated by yellow arrows. Figure 13(b) shows the same Doppler B-scan but stabilized both by phase-1 and phase-2 stabilizations. It is evident that the line artifacts are totally eliminated.

Weight for linear fit in numerical phase stabilization
In both the phase-1 and phase-2 stabilization algorithms, weighted least square fitting was utilized, as described in Sections 2.2 and 2.3. It should be noted that the OCT intensity is used as the weight in Eq. (8) of phase-1, whereas the OCT amplitude, i.e., the square root of the intensity, is utilized in Eq. (10) of phase-2.
It is known that the energy of the SNR-limited phase noise of OCT is inversely proportional to the OCT intensity as long as the OCT noise level is constant [40]. Therefore, using the OCT intensity as the weight in Eq. (8) is rational.
In contrast, the use of the square root of the intensity in Eq. (10) derives from the following compromise. The retinal blood vessels have strong scattering in OCT in general. Although the phase-2 stabilization algorithm is designed to avoid the effects of Doppler shift in the retinal vessel, the linear fit will be strongly affected by the blood flow if we utilize the OCT intensity as the weight. To reduce the effect of the blood flow and to give a larger weight to high-SNR regions, the square root of the intensity was utilized for the weight. The sufficient performance of the algorithm with this compromise was experimentally verified in Section 3.1.
In order to further demonstrate the adequacy of the square root of OCT intensity as a weight, an example is shown in Fig. 14. Figure 14(a) shows a bidirectional Doppler image which is stabilized by phase-1 stabilization followed by phase-2 stabilization but with the OCT intensity as a weight. Several vertical line artifacts are evident as indicated by arrows. On the other hand, standard phase-2 stabilization which uses the square root of the OCT intensity as the weight eliminates almost all of these artifacts as shown in Fig. 14(b).

Horizontal line artifacts in slow Doppler OCA images
Several horizontal line artifacts are observed in the slow Doppler OCA images, such as in Figs. 9(b) and 10. By observing these images, we can identify two types of line artifacts. The first type (a) (b) Fig. 14. Bidirectional Doppler images stabilized by phase-1 stabilization followed by nonstandard phase-2 stabilization which uses an OCT intensity as a weight (a) and by phase-1 stabilization followed by standard phase-2 stabilization. The arrows indicate phase artifacts occurred by the imperfection of the phase-2 stabilization. is a full-width artifact, which appears with the full horizontal transverse width, as indicated by yellow arrows in Fig. 10(a). The second type is a partial-width artifact, which appears partially in the horizontal transverse range, as indicated by red arrows in Fig. 10(a).
A full-width artifact could occur because of total decorrelation between the B-scans utilized for Doppler processing because of rapid saccadic eye motion. The disconnection of structures across the artifact line is the defining characteristic of this type of artifact.
The partial-width artifact results from the decorrelation of the A-scans utilized for Doppler calculation because of the imperfection of the axial shift cancelation. In our implementation, only axial motion between B-scans were corrected; neither transverse motion, rotation, nor skewing were corrected [19]. This property of the algorithm sometimes results in sufficient motion correction for only a part of the B-scan. In the region of insufficient correction, phase decorrelation is caused by the motion, and appears as a partial-width line artifact.

Conclusion
High-penetration Doppler optical coherence tomography and angiography were demonstrated based on a 1-μm SS-OCT. Our newly developed phase stabilization algorithm enables highsensitivity Doppler detection via SS-OCT, without the need for additional calibration hardware. As this method is purely a software method, it can be applied to most SS-OCT devices, without the need for hardware modification.
This system was utilized to examine normal eyes and a pathological case of PCV. The study with normal eyes demonstrated the effectiveness of high-penetration Doppler imaging in the choroid and ONH. The study of PCV demonstrated the clinical utility of this method.