Optical coherence angiography

Noninvasive angiography is demonstrated for the in vivo human eye. Three-dimensional flow imaging has been performe d with high-speed spectral-domain optical coherence tomography . Sample motion is compensated by two algorithms. Axial motion between adja cent A-lines within one OCT image is compensated by the Doppler shift due t o bulk sample motion. Axial displacements between neighboring im ages are compensated by a correlation-based algorithm. Three-dimensi onal vasculature of ocular vessels has been visualized. By integrating volum e sets of flow images, two-dimensional images of blood vessels are obtain ed. Retinal and choroidal blood vessel images are simultaneously obtained by separating the volume set into retinal part and choroidal parts. These a re corresponding to fluorescein angiogram and indocyanine angiogram. © 2006 Optical Society of America OCIS codes: (170.4470) Ophthalmology; (170.3890) Medical optics instr umentation; (110.6880) Three-dimensional image acquisition; (170.4500 ) Optical coherence tomography; (110.4500) Optical coherence tomography. References and links 1. J. Flammer, S. Orgul, V. P. Costa, N. Orzalesi, G. K. Krieglst ein, L. M. Serra, J.-P. Renard, and E. Stefansson, “The impact of ocular blood flow in glaucoma,” Prog. Retin. Eye R s.21, 359–393 (2002). 2. V. Patel, S. Rassam, R. Newsom, J. Wiek, and E. Kohner, “Retin al blood flow in diabetic retinopathy.” BMJ 305(6855), 678–683 (1992). 3. E. Friedman, “A hemodynamic model of the pathogenesis of age-r elated macular degeneration,” Am. J. Ophthalmol. 124, 677–682 (1997). 4. J. D. Gass, Stereoscopic atlas of macular diseases, 4th ed. (Mosby, 1997). 5. L. Yannuzzi, K. Rohrer, L. Tindel, R. Sobel, M. Costanza, W . Shields, and E. Zang, “Fluorescein angiography complication survey,” Ophthalmology 93, 611–617 (1986). 6. M. Hope-Ross, L. Yannuzzi, E. Gragoudas, D. Guyer, J. Slak ter, J. Sorenson, S. Krupsky, D. Orlock, and C. Puliafito, “Adverse reactions due to indocyanine green,” Opht halmology101, 529–533 (1994). 7. G. T. Feke, A. Yoshida, and C. L. Schepens, “Laser based ins truments for ocular blood flow assessment,” J. Biomed. Opt.3, 415–422 (1998). 8. G. T. Feke and C. E. Riva, “Laser Doppler measurements of bloo d velocity in human retinal vessels,” J. Opt. Soc. Am. 68, 526–531 (1978). 9. C. E. Riva, G. T. Feke, B. Eberli, and V. Benary, “Bidirecti onal LDV system for absolute measurement of blood speed in retinal vessels,” Appl. Opt. 18, 2301–2306 (1979). 10. C. E. Riva, S. Harino, B. Petrig, and R. Shonat, “Laser Dop pler flowmetry in the optic nerve,” Exp. Eye Res. 55, 499–506 (1992). 11. C. E. Riva, S. D. Cranstoun, J. E. Grunwald, and B. L. Petri g, “Choroidal blood flow in the foveal region of the human ocular fundus,” Invest. Ophthalmol. Vis. Sci. 35, 4273–4281 (1994). 12. R. Bonner and R. Nossal, “Model for laser Doppler measureme nts of blood flow in tissue,” Appl. Opt. 20, 2097– 2107 (1981). #71896 $15.00 USD Received 16 June 2006; revised 1 August 2006; accepted 5 August 2006 (C) 2006 OSA 21 August 2006 / Vol. 14, No. 17 / OPTICS EXPRESS 7821 13. G. Michelson, B. Schmauss, M. Langhans, J. Haraznv, and M. Groh, “Principle, validity, and reliability of scanning laser Doppler flowmetry,” J. Glaucoma 5, 99–105 (1996). 14. G. Michelson and B. Schmauss, “Two dimensional mapping of th e perfusion of the retina and optic nerve head.” Br. J. Ophthalmol. 79, 1126–1132 (1995). 15. J. Briers and A. Fercher, “Retinal blood-flow visualizat on by means of laser speckle photography,” Invest. Ophthalmol. Vis. Sci.22, 255–259 (1982). 16. Y. Tamaki, M. Araie, E. Kawamoto, S. Eguchi, and H. Fujii, “N oncontact, two-dimensional measurement of retinal microcirculation using laser speckle phenomenon,” I nvest. Ophthalmol. Vis. Sci. 35, 3825–3834 (1994). 17. X. Wang, T. Milner, and J. Nelson, “Characterization of fl uid flow velocity by optical Doppler tomography,” Opt. Lett. 20, 1337–1339 (1995). 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 h uman skin with fast scanning speed and high velocity sensitivity,” Opt. Lett.25, 114–116 (2000). 19. 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). 20. V. Westphal, S. Yazdanfar, A. M. Rollins, and J. A. Izatt, “Real-time, high velocity-resolution color Doppler optical coherence tomography,” Opt. Lett. 27, 34–36 (2002). 21. Z. Ding, Y. Zhao, H. Ren, J. S. Nelson, and Z. Chen, “Real-t ime phase-resolved optical coherence tomography and optical Doppler tomography,” Opt. Express 10, 236–245 (2002). 22. V. X. Yang, M. L. Gordon, E. Seng-Yue, S. Lo, B. Qi, J. Pekar , A. Mok, B. C. Wilson, and I. A. Vitkin, “High speed, wide velocity dynamic range Doppler optical coherenc e tomography(Part II): Imaging in vivo cardiac dynamics of Xenopus laevis,” Opt. Express 11), 1650–1658 (2003). 23. S. Yazdanfar, A. M. Rollins, and J. A. Izatt, “Imaging and v elocimetry of the human retinal circulation with color Doppler optical coherence tomography,” Opt. Lett. 25, 1448–1450 (2000). 24. S. Yazdanfar, A. M. Rollins, and J. A. Izatt, “In vivo imagi n of human retinal flow dynamics by color Doppler optical coherence tomography.” Arch. Ophthalmol. 121, 235–239 (2003). 25. D. Huang, E. A. Swanson, C. P. Lin, J. S. Schuman, W. G. Stins on, W. Chang, M. R. Hee, T. Flotte, K. Gregory, C. A. Puliafito, and J. G. Fujimoto, “Optical coherence tomogra phy,” Science254, 1178–1181 (1991). 26. J. K. Barton, J. A. Izatt, M. Kulkarni, S. Yazdanfar, and A . Welch, “Three-dimensional reconstruction of blood vessels from in vivo color Doppler optical coherence tomogra phy images,” Dermatology 198, 355–361 (1998). 27. T. Mitsui, “Dynamic range of optical reflectometry with spe ctral interferometry,” Jpn. J. Appl. Phys. 38, 6133– 6137 (1999). 28. R. Leitgeb, C. Hitzenberger, and A. Fercher, “Performanc e of fourier domain vs. time domain optical coherence tomography,” Opt. Express 11, 889–894 (2003). 29. J. F. de Boer, B. Cense, B. H. Park, M. C. Pierce, G. J. Tearn ey, and B. E. Bouma, “Improved signal-to-noise ratio in spectral-domain compared with time-domain optical coh erence tomography,” Opt. Lett. 28, 2067–2069 (2003). 30. N. Nassif, B. Cense, B. H. Park, S. H. Yun, T. C. Chen, B. E. B ouma, G. J. Tearney, and J. F. de Boer, “In vivo human retinal imaging by ultrahigh-speed spectral domain optical coherence tomography,” Opt. Lett. 29, 480–482 (2004). 31. M. Wojtkowski, V. Srinivasan, J. G. Fujimoto, T. Ko, J. S. S chuman, A. Kowalczyk, and J. S. Duker, “Threedimensional retinal imaging with high-speed ultrahigh-reso luti n optical coherence tomography,” Ophthalmology 112, 1734–1746 (2005). 32. U. Schmidt-Erfurth, R. A. Leitgeb, S. Michels, B. Povazay , S. Sacu, B. Hermann, C. Ahlers, H. Sattmann, C. Scholda, A. F. Fercher, and W. Drexler, “Three-dimensiona l ultrahigh-resolution optical coherence tomography of macular diseases,” Invest. Ophthalmol. Vis. Sci. 46, 3393–3402 (2005). 33. R. Leitgeb, L. Schmetterer, W. Drexler, A. Fercher, R. Zaw adzki, and T. Bajraszewski, “Real-time assessment of retinal blood flow with ultrafast acquisition by color Dop pler Fourier domain optical coherence tomography,” Opt. Express11, 3116–3121 (2003). 34. B. R. White, M. C. Pierce, N. Nassif, B. Cense, B. H. Park, G. J. Tearney, B. E. Bouma, T. C. Chen, and J. F. de Boer, “In vivo dynamic human retinal blood flow imaging using u ltra-high-speed spectral domain optical Doppler tomography,” Opt. Express 11, 3490–3497 (2003). 35. T. C. Chen, B. Cense, M. C. Pierce, N. Nassif, B. H. Park, S. H. Yun, B. R. White, B. E. Bouma, G. J. Tearney, and J. F. de Boer, “Spectral domain optical coherence tomography: ultra-high speed, ultra-high resolution ophthalmic imaging.” Arch. Ophthalmol. 123, 1715–1720 (2005). 36. American National Standards Institute, American National Standard for Safe Use of Lasers: ANSI Z136.1 (Laser Institute of America, Orlando, Florida, 2000). 37. 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-doma in optical coherence tomography of the human retina and optic nerve,” Opt. Express 12, 367–376 (2004). 38. Y. Yasuno, V. D. Madjarova, S. Makita, M. Akiba, A. Morosa w , C. Chong, T. Sakai, K.-P. Chan, M. Itoh, and T. Yatagai, “Three-dimensional and high-speed swept-sourc e optical coherence tomography for in vivo investi#71896 $15.00 USD Received 16 June 2006; revised 1 August 2006; accepted 5 August 2006 (C) 2006 OSA 21 August 2006 / Vol. 14, No. 17 / OPTICS EXPRESS 7822 gation of human anterior eye segments,” Opt. Express 13, 10,652–10,664 (2005). 39. B. H. Park, M. C. Pierce, B. Cense, S.-H. Yun, M. Mujat, G. J . Tearney, B. E. Bouma, and J. F. de Boer, “Realtime fiber-based multi-functional spectral-domain optical co herence tomography at 1.3 μm,” Opt. Express13, 3931–3944 (2005). 40. V. X. Yang, M. L. Gordon, A. Mok, Y. Zhao, Z. Chen, R. S. Cobb ld, B. C. Wilson, and I. A. Vitkin, “Improved phase-resolved optical Doppler tomography using Kasai velo city estimator and histogram segmentation,” Opt. Commun.208, 209–214 (2002). 41. V. X. Yang, M. L. Gordon, B. Qi, J. Pekar, S. Lo, E. Seng-Yue , A. Mok, B. C. Wilson, and I. A. Vitkin, “High speed, wide velocity dynamic range Doppler optical coherenc e tomography(Part I): System design, signal processing, and performance,” Opt. Express 11, 794–809 (2003). 42. D. Freedman and P. Diaconis, “On the histogram as a density est mator: L2 theory,” Z. Wahrscheinlichkeitstheorie verw. Gebiete57, 453–476 (1981). 43. D. W. Scott,Multivariate Density Estimation: Theory, Practice, and Visualization (John Wiley & Sons, Inc., 1992). 44. L. Ferman, H. Collewijn, T. C. Jansen, and A. V. V. den Berg, “


Introduction
Ocular circulation is very important not only for ophthalmic diagnosis but also for the study of eye diseases such as glaucoma [1], diabetic retinopathy [2], and age-related macular degeneration [3]. Fluorescein angiography (FA) and indocyanine green angiography (ICGA) [4] are the chief diagnostic methods for retinal diseases. Ocular vessels are contrasted by using the fluorescence of the dyes injected into a vein. Since the excitation wavelength of indocyanine green is relatively longer than that of sodium fluorescein, the choroidal vasculature is clearly visualized in ICGA rather than FA. However, these methods are invasive. The skin will be colored in yellow by the fluorescein dye. Intravenous administrations of sodium fluorescein and indocyanine green are considered as safe procedures but not no risk. There are some adverse effects due to dye injection. A frequencies of mild adverse effects, such as nausea, vomiting, and so forth, are between less than 1 % and 10 % for FA [5] and 0.15 % for ICGA [6]. Although rare, some severe adverse effects such as anaphylaxis also occur.
There are several noninvasive blood flow assessment techniques that use lasers, which have been applied in several investigations of ocular circulation [7]. Laser Doppler velocimetry (LDV) measures the blood flow velocity, and it has been applied to human retinal vessels [8,9]. Microcirculation in the optic nerve head (ONH) [10] and the choroidal blood flow [11] of the human eye have been investigated by laser Doppler flowmetry (LDF) [12], which acquires the volume, flux, and average velocity of the red blood cells in the small vessels surrounded by the scattering tissues. Two-dimensional perfusion mapping of blood flow can be obtained by using scanning laser Doppler flowmetry (SLDF) [13,14], laser speckle photography [15], and laser speckle flowgraphy [16]. Although these techniques are useful for investigating the ocular circulation, they demonstrate no or poor ability to resolve the depth structure. SLDF, which is based on scanning laser ophthalmoscopy, exhibit depth-resolving power inherently due to the confocal effect; however, its axial resolution is limited to ∼ 300 µm due to eye aberrations.
Optical Doppler tomography (ODT) [17] has been applied to the depth-resolved crosssectional flow imaging of in vivo tissues [18,19,20,21,22]. Further, ODT has been demonstrated for blood flow imaging of retinal vessels [23,24]. These ODTs are based on time- [25] -a cross-sectional imaging technique with high spatial resolution. Three-dimensional vessel structure extraction is possible using threedimensional scanning [26]. Three-dimensional retinal imaging is, however, restricted by a substantial and frequent head movement and exposure limit for the eye. Short imaging time is required to suppress the influence of the head movement, but the scanning speed of OCT is restricted to maintain sustainable sensitivity with low power of probing beam. In recent years, spectral-domain OCT (SD-OCT), which has a higher sensitivity than time-domain OCT [27, 28, 29], has emerged. SD-OCT has been performed at a 150 fold faster measurement speed than time-domain OCT [30]. This technique allows high-resolution and highly accurate threedimensional in vivo structural measurement of the human retina [31, 32]. High-speed retinal blood flow imaging by using SD-OCT has been demonstrated [33,34,35].
We present an in vivo noninvasive angiography of the human retina. Three-dimensional retinal imaging is performed by using high-speed SD-OCT. The bulk motion artifacts in the flow images are compensated by a histogram-based algorithm. The axial eye movements during three-dimensional imaging are compensated by using Doppler shifts of the bulk motion for axial displacement between adjacent A-lines and a simple and fast correlation-based algorithm for axial shifts between neighboring images. Three-dimensional retinal and choroidal vasculature are visualized, and two-dimensional vessel imaging that is comparable to FA and ICGA are demonstrated.

Method
For retinal imaging, we constructed a high-speed SD-OCT system with an 840-nm band light source. The schematic of the developed system is shown in Fig. 1. A superluminescent diode (Superlum, Russia, SLD-37) is used as the light source; this SLD has a center wavelength of 840 nm and FWHM Spectral Bandwidth Of 50 Nm. A Free-Space Optical Isolator Is Placed Immediately After The Light Source. The beam is split by a 20/80 fiber coupler, and 20% of the beam is guided to the sample arm. A 78D lens is used to deliver the beam to the posterior part of the human eye. The head of the volunteer is fixed on a handmade stage. There is a fixation target to stabilize the other eye that is not imaged by OCT. The optical power of the beam on the cornea is 700 µW, which is lower than the ANSI exposure limit [36]. The beam diameter on the cornea is set to ∼ 0.6 mm in FWHM in order to suppress the effect of eye aberrations and obtain a large confocal parameter. The spot size on the retina is calculated to be ∼ 30 µm, and the confocal parameter in the tissue is 2 mm (n = 1.38). In the detection arm, the backscattered signal light and reference light are combined and guided to the spectrometer that consists of a transmission grating (Wasatch Photonics, 1200 lp/mm), achromatic lens (f = 200 mm), and high-speed line-scan CCD camera (Basler, L103k-2k). The CCD camera detects the spectral interference signals.
First, the detected spectrum is processed to eliminate the periodic electronic noise with a period of 4 pixels. This noise appears in OCT images as a fixed-pattern noise. Since this noise is unstable in phase, it cannot be eliminated with fixed-pattern noise rejection by using an averaged spectrum, which is mentioned later. Alternatively, a Fourier harmonic of the spectrum at a frequency of 1/4 pixel −1 (period of 4 pixels) is obtained and subtracted from the spectrum. Other optical fixed-pattern noises are eliminated by using the averaged spectrum of all the spectra in one image [37]. The nonlinearity of the spectrometer is obtained by analyzing the spectral phase of the OCT signal from a single reflection plane. Then, the spectrum is linearly resampled in k-space [38]. Digital inverse Fourier transform of this spectrum yield a complex signal of an A-line profile.
This system performs OCT imaging with an A-line rate of 18.7 kHz. The exposure time of one A-line is 53.1 µs. Two-dimensional tomographic imaging with 1024 A-lines is achieved at a frame rate of 18 fps, and a three-dimensional volume comprising 1024 × 140 A-lines is obtained within 7.7 s. The axial resolution is estimated to be 6.2 µm in air, whereas the measured axial resolution is 8.8 µm in air. This correspond to 6.4 µm in tissues, where the refractive index is assumed to be 1.38. This discrepancy from the theoretical value may be due to the non-Gaussian shape of the source spectrum. The sensitivity of this system is measured to be 99.3 dB at a depth of 100 µm; this value reduces to 92.9 dB at a depth of 2 mm.

Flow imaging
In our method, the Doppler shift of the OCT signal is utilized to contrast the blood vessels. To achieve the flow images, first, the bulk motion artifacts are eliminated by a histogram-based algorithm. Then, two flow imaging implementations, bi-directional flow and power of Doppler shift, are applied. Finally, these two images are averaged using a moving average filter.
A bi-directional flow image f i is obtained by using the phase difference between adjacent A-lines ∆φ i [18] as follows: Here, T denotes the A-line period and i ∈ [0, N) denotes A-line number. By using the A-line period, the maximum detectable Doppler shift is estimated as f max = ± 1/2T . In our system, this value is calculated to be ± 9.38 kHz. Further, the minimum detectable Doppler shift is defined on the basis of the fluctuation of the phase difference. As mentioned by Park et.al, the minimum detectable Doppler shift depends on the signal-to-noise ratio and the ratio of the spacing of transversal scanning to the beam spot size [39]. Another implementation, which is power of Doppler shift imaging, is also applied [39]. The power of the Doppler shift image σ 2 i is obtained as and it exhibits similar properties to Doppler standard deviation imaging [19]. During retinal imaging, significant involuntary head movements cause large artifacts in the flow images. To compensate for these bulk motion-induced artifacts, we use a histogram-based bulk motion estimation method [40,41] and phase wrapping method [34]. The diagram of the algorithm is shown in Fig. 2. The phases of adjacent A-lines φ i (z), φ i+1 (z) are subtracted (φ i+1 (z) − φ i (z)) and wrapped into -π to π. From the differences, the phase difference ∆φ bulk corresponding to the bulk motion is obtained by the histogram-based bulk motion compensation method. ∆φ bulk i is subtracted from the raw phases of the A-line that follow φ i+1 (z); then, phases are wrapped . Then, the phase difference ∆φ i (z) is obtained by subtracting φ i (z) from the compensated phase φ ′ i+1 (z) and phase wrapping. φ ′ i+1 (z) is used for the compensation of next A-line.
When the vessel region is smaller than the other tissue region or the blood flow distribution exhibits a wide velocity range, the phase difference corresponding to the Doppler shift frequency due to the bulk motion will be the maximum count in the histogram of the phase differences between two adjacent A-lines. In other words, the Doppler shift of bulk motion is estimated by using the mode of the phase differences. Yang et al. determined the number of histogram bins from the phase stability [41]. However, the histogram shape also depends on the number of samples. If the number of samples is small, the histogram bin width will be relatively small, resulting in a noisy histogram. We apply the following data-based histogram bin width determination rule, which is known as the Freedman & Diaconis rule [42]: where h denotes the histogram bin width;, IQ, the interquartile range; and m, the number of samples. Prior to the formulation of histograms, the noise regions are eliminated by using an intensity-based thresholding. The threshold level is determined from the mean values of the noise floor in the OCT images. Although a smooth histogram can be obtained by using a large histogram bin width, another problem exists. Since the histogram bin width increases, the accuracy of the bulk motion calculation decreases. In order to increase the accuracy of measurement of the bulk motion Doppler shifts, we average several histograms with different origins of bins; the resulting histogram is referred to as an averaged shifted histogram [43].
The histograms of one depth profile of the phase differences in the region lacking vessels are shown in Fig. 3. The histogram in which the number of bins is determined by the phase-noisebased rule shows an almost zero-mean probability distribution ( Fig. 3(a)). Here, we estimate the phase noise by using the lateral spacing and beam spot size [39]. However, the maximum count is observed at approximately 1 radian. This error might be due to the extremely narrow histogram bin width. By using the Freedman & Diaconis rule ( Fig. 3(b)), such types of error peaks disappear. Instead, the histogram bin widths broaden. The accuracy of the bulk motion Doppler shift measurements decreases. By averaging the eight histograms with different origin of bins and identical histogram bin width as determined by the Freedman & Diaconis rule, a very smooth phase difference distribution is obtained, as shown in Fig. 3(c). For comparison,   Fig. 3(b). This histogram is extremely noisy, and the error in the maximum count is observed at approximately 1 radian. Figure 4 shows flow images of the retinal vessels near the ONH with bulk motion compensation by using different histogram bin width determination methods. Bi-directional flow images are obtained after bulk motion compensation, differentiation of phase, and 19 × 3 pixels (29 µm in the lateral direction × 10 µm in the axial direction) averaging filtering. In Fig. 4

(A)
(phase-noise-based method), bulk motion artifacts remain in the upper part of the blood vessels, which is indicated by the white arrow. These artifacts are suppressed in Fig. 4 the other for compensation for each cross section-are described here.

Axial motion compensation between adjacent A-lines by bulk motion Doppler shift
The Doppler shift in bulk motion causes artifacts in the flow images; they can be eliminated as mentioned in section 2.1. To eliminate these artifacts, the Doppler shifts of the bulk tissue motion at each A-line are obtained. From these signals, the displacement of the tissue can also be estimated. Since the Doppler shift is caused by the movement of a sample along the optical axis, only the axial motion will be compensated. The axial velocity of the sample is expressed using the Doppler shift of bulk motion as follows: where v zi is the averaged axial velocity of the sample between the (i − 1)-th and i-th A-line; λ 0 , the center wavelength; δ φ bulk i , the phase difference due to the axial bulk motion; and n, the refractive index of the tissue. Then, the relative axial displacement between the adjacent A-lines is calculated as When the origin is set to the 0-th A-line, the axial displacement for each A-line is expressed as where ∆φ bulk is the phase differences corresponding to the mode value of histogram (Sec. 2.1) after phase unwrapping. In eq. (6), the change in the shape of the entire eye is not considered, such as the volume change due to the pressure of perfusion and accommodations. The OCT images are restored by shifting each A-line along the depth by using ∆z i .

Motion compensation between neighboring frames
The lateral movement of retinal OCT images by head movements is quit small comparing to axial movement since volunteers are gazing a distant stationary target. The lateral displacement by the eye movement is estimated as ∼ 40 µm from about 0.13 degree (standard deviation) of the gaze stability with head free and stationary fixation [44] and 17 mm of the distance between back nodal point to retina of the human eye [45]. It is comparable to the beam spot size on the retina (30 µm), even though the large saccadic eye movement occurs sometimes which cause significant large displacement. Thus, we use only the axial correlation for motion compensation of three-dimensional images. The axial shifts due to motion are estimated by the following one-dimensional cross correlation between the i-th A-lines in the ( j + 1)-th and j-th frames: where χ is the cross-correlation function, Γ[]() is the cross-correlation operation, and I i,J (z) is the i-th A-line intensity profile in the j-th frame. The positions of the maximum peaks of each χ(z ′ ) denote the axial shifts. When a maximum correlation of χ k (z ′ ) is smaller than the average of maximum correlations, i.e. 1/N ∑ N−1 i=0 max[χ i (z ′ )], χ k (z ′ ) is discarded. The axial displacement between the j-th and ( j + 1)-th frames is determined as the median value of the maximum positions in each χ.
A displacement that includes the lateral shift between neighboring frames is simply determined by a two-dimensional cross correlation. The retina has a well-layered structure: therefore, the cross-correlation function between two retinal OCT images exhibits a fairly sharp peak along the axial direction. Motion compensation in the axial direction works. In contrast, the shape of the cross-correlation function along the transversal direction is broadened. Motion compensation in lateral direction depends on each case. When a characteristic landmark such as the ONH is present and imaging is performed over a large area, lateral motion compensation can be applied; otherwise, it will not be effective. Further improvement is required in the lateral motion compensation technique.

Performance of axial shift compensation
Cross-sectional imaging with a large number of A-lines can be distorted by the sample motion. Figure 5(A) shows a cross section of the human retina, which comprises 4000 A-lines. The acquisition time per image is 214 ms. Image distortion occurs due to sample motion at the region indicated by the black arrow. The axial displacement is compensated by using the motion compensation algorithm for adjacent A-lines described in Sec. 2.2.1, as shown in Fig. 5(B). The bi-directional flow and the power of Doppler shift images are shown in Figs. 5(C) and (D), respectively. Two small retinal vessels are evident. In the power of Doppler shift image, some choroidal vessels are observed. In the movie, the axial motion between each frame is compensated using the correlation-based method (section 2.2.2). Figure 6 shows a comparison of the motion compensation results for a three-dimensional OCT volume set. Figure 6 compensations, with the motion compensation between neighboring frames, and with motion compensations between adjacent A-lines and between neighboring frames, respectively. The cross-sectional image along the slow-scanning direction shows significant axial displacement without any motion compensation, as shown in Fig. 6(B). The imaging time of this image is extremely longer than that of images along to fast-scanning direction; roughly several hundred times. Involuntary axial head movement cause large distortions. Since the head is not bound to stage and volunteers are trying to gaze a target, low frequency distortions maybe caused by small vibrations of muscles due to tension. By axial motion compensation between neighboring frames described in Sec. 2.2.2, the axial displacement is reduced; however, a small fluctuation remains due to the distortion of each image ( Fig. 6(C)). Although cross-sectional imaging with the conventional density of an A-line does not exhibit an apparent distortion, some shearing does exist. This affects the performance of motion compensation between neighboring frames. This distortion is improved by axial motion compensation between adjacent A-lines, and a smoother structure is evident (Fig. 6(D)). The mean values of normalized cross-correlation between neighboring frames are calculated for a three-dimensional OCT volume set of the macular region, which comprises 1024 × 140 A-lines (5 mm × 5 mm). These are 0.25 without any motion compensation, 0.33 with the motion compensation between neighboring frames, and 0.38 with the motion compensations between adjacent A-lines and between neighboring frames. The increment of the normalized cross-correlation is + 0.08 ± 0.08 (mean ± standard deviation) through only the motion compensation between neighboring frames and increased more (+ 0.05 ± 0.04) with the motion compensation between adjacent A-lines.

Optical coherence angiography
Three-dimensional volume sets of the flow images are used to generate three-dimensional vasculature images and en-face vessel images. To eliminate the noise in the regions with no signal, only the areas of interest are extracted, i.e., from the vitreoretinal interface to the deepest penetration points, where signals are equal to noise floor. Two boundaries of the ares are segmented in the OCT images. Furthermore, a high reflectivity layer -corresponding to the inner/outer segment boundary (IS/OS), retinal pigment epithelium (RPE), and choriocapillaris (CC) -is segmented in order to separate the retinal and choroidal layers. The OCT signals exhibit low intensities in some retinal vessels since the high blood flow velocity in the vessel results in a fringe washout. The anterior boundary of the ares of interest is determined using an algorithm similar to that presented by Mujat et al. [46]. In order to suppress noise, the OCT images are blurred by using a Gaussian filter with a standard deviation of 3 × 3 pixels; gradient magnitude images are thus obtained. These images are indicated with the sign of the axial gradient to distinguish whether anterior edge or posterior edge.The thresholding of the gradient magnitude images yields binary edge images. The threshold level is determined from the statistical properties of noise to maintain the interfaces as a continuous curve. After eliminating the small particles in the binary edge images, the first exceeding points over the threshold level from the top of the images are detected. Then, positions that have a maximum gradient within ten pixels lower than the first unity points are obtained. In order to smooth the curve, an iterative smoothing algorithm is used. The large second derivative of axial positions of boundary are ignored, and then, a linear interpolation of the boundary positions is performed. These steps are applied until the curve is significantly smooth. Then, a median filter with a rank of 15 is applied.
The high reflectivity layer is segmented prior to posterior boundary segmentation to determine the retinal vessel regions. The axial positions of the maximum gradient lower than 25 pixels from the anterior boundary are detected. Some A-lines with no edges are ignored. The intensities of 15 pixels around the detected positions along the axial direction are averaged. Then, they are smoothed using a Savitzky-Golay filter. The points at which the intensity before smoothing significantly drops from the smoothed results are considered to form the vessel region. The signals under the blood vessels are weak due to the absorption of blood. Both the vessel regions and the no-edge regions are expanded by 10 pixels on both sides to completely eliminate the regions. The positions of the high reflectivity layer of both regions are then obtained by linear interpolation. The positions of this layer are smoothed in a similar manner as the segmentation of the anterior boundary.
For the posterior boundary segmentation, the boundaries between the signals and noises are determined. Binary images are created by thresholding for the Gaussian blurred intensity images. The threshold level is determined by using a histogram of logarithmic intensity. The probability density of noise exhibits a sharp peak; therefore, the threshold level is set to a value which exhibits a count of a fifth part of the maximum count in the histogram. Then, the small particles are eliminated. The threshold of the particle size is determined by the image size. The first unity points from the posterior edge of the images are detected. The posterior boundary curve is interpolated at no-edge and vessel regions. Then, the curve is smoothed using the same method to that is performed for the anterior boundary and high reflectivity layer.
All the segmentation parameters are automatically determined. By using a PC (Athlon 64 3500+ processor and a 2 GB RAM), the processing, which includes flow imaging, motion compensation, and segmentations, takes approximately 11.5 seconds per frame that comprises 1024 A-lines.
The three-dimensional vasculature is imaged by a volume rendering of the segmented volumes. The bi-directional flow images are squared to achieve a high dynamic range for visualizing the vasculature. For both the three-dimensional bi-directional and the power of Doppler shift volume sets, three-dimensional Gaussian filtering is performed to suppress the spurious noise and discontinuities in the vasculature.
Two-dimensional angiograms are produced by integrating the flow volume sets along the axial direction. Since the power of the Doppler shift images have a higher sensitivity than the bi-directional flow images, angiograms produced from the former exhibit a higher contrast. Furthermore, the retinal vessels can be separated from the choroidal vessels by segmenting the high reflectivity layer. Two angiograms of the retinal and choroidal vessels can be produced individually, which might correspond to FA and ICGA.

Results and discussion
In vivo three-dimensional blood flow imaging has been performed in a healthy human eye. The macular region of the retina is imaged with 1024 A-lines in one frame and 140 frames acquired over an area of 5 mm × 5 mm. The resulting images are shown in Fig. 7. Figure 7(A) shows the fundus projection image of the OCT volume set. Three boundaries segmented by using the algorithm that is described in Sec. 2.3 are indicated in the sequence of OCT images. In the bi-directional flow image ( Fig. 7(C)), the blood flow from some retinal vessels is observed. However, the flow in some vessels and the choroidal vessels is not visible. In contrast, the power of Doppler shift image (Fig. 7(D)) clearly shows the flow in the retinal vessels and some choroidal vessels. However, artifacts are observed in the low signal regions (i.e., vitreous, outer nuclear layer, and deeper region than penetration). Artifacts also appear under the blood vessels due to the absorption of blood and forward scattering by the red blood cells. Threedimensional ocular vasculature is presented using volume rendering images of both the squared bi-directional flow and the power of Doppler shift volume sets. These are volume rendered after three-dimensional Gaussian filtering with a standard deviation of 1 × 1.5 × 1 pixels (axial direction × B-scanning direction × C-scanning direction). The volume rendering image of the squared bi-directional flow volume (Fig. 7(E)) shows the vasculature of the retinal vessels; however, the choroidal vessels are not clearly observed. Artifact are observed at the surface that have high specular reflection. Although artifacts are observed under the retinal vessels due to the shadowing effect, and substantial noise occurs at the outer nuclear layer, the volume rendering image of the power of Doppler shift volume set ( Fig. 7(F)) shows the retinal vasculature and choroidal vessels.
The retinal and choroidal layers can be separated at the high reflectivity layer corresponding to IS/OS, RPE, and CC. A combination of the retinal part from Fig. 7 (E) and the choroidal part from Fig. 7(F) is shown in Fig. 7(G). The anterior part is encoded in cyan and the posterior part, yellow. Since the three-dimensional squared bi-directional flow images exhibit low noise and small artifacts and the three-dimensional power of Doppler shift images exhibit choroidal vessels, this combination results in a better visualization. Three-dimensional volume rendering images of the ONH region are shown in Fig. 8 and its movies.
Although the power of Doppler shift images exhibit some artifacts, they are highly sensitive. The projection of the power of Doppler shift volume set serves as two-dimensional angiography. In this case, the artifacts under the vessels enhance the contrast of the projection images. Twodimensional angiographies for both the macular and ONH regions are shown in Figs ). On the other hand, the projection images of the choroidal part reveal not only the choroidal vessels but also the retinal vessels due to shadowing effect (Figs. 9(E) and 10(E)). For practical applications, composite false-colored images are created from these projection image sets to distinguish the retinal vessels and the choroidal vessels. For the colored images, each of the three layers -red, green, and blue -are produced by the following manner. The images of the retina and the choroid are normalized; they correspond to the blue and green channels, respectively. The retinal images are inverted and multiplied with the images of the choroid in order to suppress the retinal vessel artifacts; this is the yellow channel. This process renders the retinal and choroidal vessels cyan and yellow, respectively (Figs. 9(F) and 10(F)).
Since there are no retinal and choroidal layers in the optic disk, the segmentation of the high reflectivity layer at the optic disk is inappropriate. However, the high reflectivity layer segmentation algorithm behaves to detect the posterior boundary of retinal vessels, so that the retinal vessels are included into the retinal part. Since the almost of all OCT signals in the retinal vessels at the optic disk are disappeared due to fringe washout, large rising edges at  the posterior boundary of the retinal vessels are detected. Nevertheless, the discontinuity of the retinal vessel inside the optic disk is observed in Fig. 10(B). Further improvement of the segmentation algorithm is required. Figure 11 shows a comparison between FA, ICGA, and two-dimensional OCA of the macular region. The angiograms of the ONH are also shown in Fig. 12. Intravenous injection has been administrated with sodium fluorescein (2.5 ml of 10 percent solution) and indocyanine green (2 ml of 1.25 percent solution). Each fluorescence angiograms are acquired at 2 minutes 43 seconds (Fig. 11(B)), 14 minutes 56 seconds (Fig. 12(B)), 27 seconds ( Fig. 11(C)), and 45 seconds (Fig. 12(C)) after dye injection. Both the two-dimensional OCA images of the macular region and the ONH shows the vasculature of the retinal vessels and choroidal vessels. Fine retinal vessels in the macular region are visible with FA, but not with OCA. It is expected that the density of the A-line is not sufficient for detecting these small vessels and/or their blood flow velocity is low. Some thick choroidal vessels that appear in ICGA do not appear in the OCA images. The reasons for this are assumed to be their low blood flow velocity and/or the low backscattered light from these vessels. However, these OCA images are in good agreement to FA and ICGA.
FA and ICGA are also capable of investigating circulation by dye-dilution method. Macroscopic blood flow dynamics, e.g. circulation time, will be investigated which is difficult for OCA. Instead of that, OCA are capable of measuring the local blood flow speed, not absolute value but relative value, that is hard for fluorescence angiographies.
There is an alternative en-face method to contrast only the retinal vessels from the OCT volume sets by utilizing the absorption of blood, which is referred to as the shadowgram [47]. It might be useful to create composite colored images of the macular region. The comparison between an OCA image of the retina and a shadowgram, composite images from each image and an OCA image of the choroid are shown in Fig. 13. A finer retinal vasculature can be observed in the shadowgram and the composite image than in the OCA image of the retina. On the other hand, a OCT intensity based choroidal vessel imaging requires significantly deep penetration to the choroid. The high-speed swept-source OCT in the 1 micron wavelength region has been performed the deeper penetration and the high contrast choroidal vessel imaging [48].

Conclusion
We present a new noninvasive angiography for the human eye. By using three-dimensional flow imaging with high-speed SD-OCT, the three-dimensional vasculature of the posterior part of the human eye is visualized. The retinal and choroidal parts are distinguished in the OCT images by segmentation. The two-dimensional vessel images of the retina and choroid represent the ocular vessels corresponding to FA and ICGA.