# In vivo laser Doppler holography of the human retina

L. Puyo, M. Paques, M. Fink, José-Alain Sahel, Michael Atlan

## ► To cite this version:

L. Puyo, M. Paques, M. Fink, José-Alain Sahel, Michael Atlan. In vivo laser Doppler holography of the human retina. Biomedical optics express, 2018, 9 (9), pp.4113-4129. 10.1364/BOE.9.004113. hal-01875560

HAL Id: hal-01875560

https://hal.sorbonne-universite.fr/hal-01875560

Submitted on 17 Sep 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# In vivo laser Doppler holography of the human retina

L. PUYO,<fnref n="1" />,* M. PAQUES,<fnref n="2" />,<fnref n="3" /> M. FINK,<fnref n="1" /> J.-A. SAHEL,<fnref n="2" />,<fnref n="3" /> AND M. ATLAN<fnref n="1" />

**Abstract:** The eye offers a unique opportunity for the non-invasive exploration of cardiovascular diseases. Optical angiography in the retina requires sensitive measurements, which hinders conventional full-field laser Doppler imaging schemes. To overcome this limitation, we used digital holography to perform laser Doppler perfusion imaging of human retina with near-infrared light. Two imaging channels with a slow and a fast CMOS camera were used simultaneously for real-time narrowband measurements, and offline wideband measurements, respectively. The beat frequency spectrum of optical interferograms recorded with the fast (up to 75 kHz) CMOS camera was analyzed by short-time Fourier transformation. Power Doppler images drawn from the Doppler power spectrum density qualitatively revealed blood flow in retinal vessels over 512 × 512 pixels covering 2.4 × 2.4 mm<fnref n="2" /> on the retina with a temporal resolution down to 1.6 ms. The sensitivity to lateral motion as well as the requirements in terms of sampling frequency are discussed.

© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

OCIS codes: (090.0090) Holography; (170.4460) Ophthalmic optics and devices; (170.3340) Laser Doppler velocimetry.

## References and links

1. S. H. Hardarson and E. Stefánsson, "Retinal oxygen saturation is altered in diabetic retinopathy," Br. J. Ophthalmol. **96**, 560–563 (2012).

2. B. Pemp and L. Schmetterer, "Ocular blood flow in diabetes and age-related macular degeneration," Can. J. Ophthalmol. Can. d'Ophtalmologie **43**, 295–301 (2008).

3. A. P. Cherecheanu, G. Garhofer, D. Schmidl, R. Werkmeister, and L. Schmetterer, "Ocular perfusion pressure and ocular blood flow in glaucoma," Curr. Opin. Pharmacol. **13**, 36–42 (2013).

4. J. A. Briers and S. Webster, "Quasi real-time digital version of single-exposure speckle photography for full-field monitoring of velocity or flow fields," Opt. Comm. **116**, 36–42 (1995).

5. A. H. Kashani, C.-L. Chen, J. K. Gahm, F. Zheng, G. M. Richter, P. J. Rosenfeld, Y. Shi, and R. K. Wang, "Optical coherence tomography angiography: A comprehensive review of current methods and clinical applications," Prog. Retin. Eye Res. (2017).

6. R. K. Wang, L. An, P. Francis, and D. J. Wilson, "Depth-resolved imaging of capillary networks in retina and choroid using ultrahigh sensitive optical microangiography," Opt. Lett. **35**, 1467–1469 (2010).

7. Y. Jia, J. C. Morrison, J. Tokayer, O. Tan, L. Lombardi, B. Baumann, C. D. Lu, W. Choi, J. G. Fujimoto, and D. Huang, "Quantitative OCT angiography of optic nerve head blood flow," Biomed. Opt. Express **3**, 3127–3137 (2012).

8. W. Choi, K. J. Mohler, B. Potsaid, C. D. Lu, J. J. Liu, V. Jayaraman, A. E. Cable, J. S. Duker, R. Huber, and J. G. Fujimoto, "Choriocapillaris and choroidal microvasculature imaging with ultrahigh speed OCT angiography," PloS One **8**, e81499 (2013).

9. Y. Jia, S. T. Bailey, T. S. Hwang, S. M. McClintic, S. S. Gao, M. E. Pennesi, C. J. Flaxel, A. K. Lauer, D. J. Wilson, J. Hornegger, J. G. Fujimoto, and D. Huang, "Quantitative optical coherence tomography angiography of vascular abnormalities in the living human eye," Proc. Natl. Acad. Sci. **112**, E2395–E2402 (2015).

10. Z. Chen, T. E. Milner, S. Srinivas, X. Wang, A. Malekafzali, M. J. van Gemert, and J. S. Nelson, "Non-invasive imaging of in vivo blood flow velocity using optical Doppler tomography," Opt. Lett. **22**, 1119–1121 (1997).

11. J. A. Izatt, M. D. Kulkarni, S. Yazdanfar, J. K. Barton, and A. J. Welch, "In vivo bidirectional color Doppler flow imaging of picoliter blood volumes using optical coherence tomography," Opt. Lett. **22**, 1439–1441 (1997).

<footnote n="1">¹Institut Langevin, Centre National de la Recherche Scientifique (CNRS), Paris Sciences & Lettres (PSL Research University), École Supérieure de Physique et de Chimie Industrielles (ESPCI Paris) - 1 rue Jussieu, 75005 Paris, France</footnote>

<footnote n="2">²Institut de la Vision, INSERM UMR-S 968, CNRS UMR 7210, UPMC, 17 rue Moreau, 75012 Paris, France</footnote>

<footnote n="3">³Centre d'Investigation Clinique (CIC) Centre Hospitalier National d'Ophtalmologie des Quinze-Vingts, INSERM, 28 rue de Charenton, 75012 Paris, France</footnote>

<footnote>*gl.puyo@gmail.com</footnote>
12. 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).

13. R. A. Leitgeb, R. M. Werkmeister, C. Blatter, and L. Schmetterer, "Doppler Optical Coherence Tomography," Prog. Retin. Eye Res. 41, 26-43 (2014).

14. H. Fujii, "Visualisation of retinal blood flow by laser speckle flowgraphy," Med. Biol. Eng. Comput. 32, 302-304 (1994).

15. T. Sugiyama, M. Araie, C. E. Riva, L. Schmetterer, and S. Orgul, "Use of laser speckle flowgraphy in ocular blood flow research," Acta Ophthalmol. 88, 723-729 (2010).

16. C. Riva, S. Harino, B. Petrig, and R. Shonat, "Laser Doppler flowmetry in the optic nerve," Exp. Eye Res. 55, 499-506 (1992).

17. M. Leahy, F. De Mul, G. Nilsson, and R. Maniewski, "Principles and practice of the laser-Doppler perfusion technique," Technol. Heal. Care 7, 143-162 (1999).

18. C. E. Riva, M. Geiser, and B. L. Petrig, "Ocular blood flow assessment using continuous laser Doppler flowmetry," Acta Ophthalmol. 88, 622-629 (2009).

19. V. Rajan, B. Varghese, T. G. van Leeuwen, and W. Steenbergen, "Review of methodological developments in laser Doppler flowmetry," Lasers Med. Sci. 24, 269-283 (2009).

20. G. Michelson, B. Schmauss, M. Langhans, J. Harazny, and M. Groh, "Principle, validity, and reliability of scanning laser Doppler flowmetry." J. Glaucoma. 5, 99-105 (1996).

21. L. Schmetterer and G. Garhofer, "How can blood flow be measured?" Surv. Ophthalmol. 52, S134-S138 (2007).

22. A. Serov, W. Steenbergen, and F. de Mul, "Laser Doppler perfusion imaging with complementary metal oxide semiconductor image sensor," Opt. Lett. 27, 300-302 (2002).

23. H. Spahr, D. Hillmann, C. Hain, C. Pfäffle, H. Sudkamp, G. Franke, and G. Hüttmann, "Imaging pulse wave propagation in human retinal vessels using full-field swept-source optical coherence tomography," Opt. letters 40, 4771-4774 (2015).

24. H. Sudkamp, P. Koch, H. Spahr, D. Hillmann, G. Franke, M. Münst, F. Reinholz, R. Birngruber, and G. Hüttmann, "In-vivo retinal imaging with off-axis full-field time-domain optical coherence tomography," Opt. Lett. 41, 4987-4990 (2016).

25. D. J. Fechtig, B. Grajciar, T. Schmoll, C. Blatter, R. M. Werkmeister, W. Drexler, and R. A. Leitgeb, "Line-field parallel swept source MHz OCT for structural and functional retinal imaging," Biomed. Opt. Express 6, 716-735 (2015).

26. D. Hillmann, H. Spahr, C. Hain, H. Sudkamp, G. Franke, C. Pfäffle, C. Winter, and G. Hüttmann, "Aberration-free volumetric high-speed imaging of in vivo retina," Sci. Reports 6, 35209 (2016).

27. L. Ginner, T. Schmoll, A. Kumar, M. Salas, N. Pricoupenko, L. M. Wurster, and R. A. Leitgeb, "Holographic line field en-face OCT with digital adaptive optics in the retina in vivo," Biomed. Opt. Express 9, 472-485 (2018).

28. M. Pellizzari, M. Simonutti, J. Degardin, J.-A. Sahel, M. Fink, M. Paques, and M. Atlan, "High speed optical holography of retinal blood flow," Opt. letters 41, 3503-3506 (2016).

29. H. Albrecht, M. Borys, N. Damaschke, and C. Tropea, *Laser Doppler and Phase Doppler Measurement Techniques* (Springer. 2003).

30. C. E. Riva, J. E. Grunwald, S. H. Sinclair, and B. Petrig, "Blood velocity and volumetric flow rate in human retinal vessels." Investig. Ophthalmol. & Vis. Sci. 26, 1124-1132 (1985).

31. C. Magnain, A. Castel, T. Boucneau, M. Simonutti, I. Ferezou, A. Rancillac, T. Vitalis, J.-A. Sahel, M. Paques, and M. Atlan, "Holographic laser Doppler imaging of microvascular blood flow," JOSA A 31, 2723-2735 (2014).

32. J. W. Goodman, *Introduction to Fourier Optics* (Roberts and Company Publishers, 2005).

33. M. Guizar-Sicairos, S. T. Thurman, and J. R. Fienup, "Efficient subpixel image registration algorithms," Opt. Lett. 33, 156-158 (2008).

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 ultra-high-speed spectral domain optical Doppler tomography," Opt. Express 11, 3490-3497 (2003).

35. H. N. Sabbah and P. D. Stein, "Valve origin of the aortic incisura," Am. J. Cardiol. 41, 32-38 (1978).

36. M. Stern, D. Lappe, P. Bowen, J. Chimosky, G. Holloway, H. Keiser, and R. Bowman, "Continuous measurement of tissue blood flow by laser-Doppler spectroscopy." Am. J. Physiol. 232, H441-H448 (1977).

37. J. F. Polak, J. M. Alessi-Chinetti, A. R. Patel, and J. M. Estes, "Association of common carotid artery Doppler-determined dicrotic notch velocity with the left ventricular ejection fraction," J. Ultrasound Medicine 34, 461-467 (2015).

38. V. Doblhoff-Dier, L. Schmetterer, W. Vilser, G. Garhöfer, M. Gröschl, R. A. Leitgeb, and R. M. Werkmeister, "Measurement of the total retinal blood flow using dual beam Fourier-domain Doppler optical coherence tomography with orthogonal detection planes," Biomed. Opt. Express 5, 630-642 (2014).

39. Z. Zhong, B. L. Petrig, X. Qi, and S. A. Burns, "In vivo measurement of erythrocyte velocity and retinal blood flow using adaptive optics scanning laser ophthalmoscopy," Opt. Express 16, 12746-12756 (2008).

40. R. Bonner and R. Nossal, "Model for laser Doppler measurements of blood flow in tissue." Appl. Opt. 20, 2097-2107 (1981).

41. R. F. Bonner and R. Nossal, "Principles of laser-Doppler flowmetry," in "Laser-Doppler Blood Flowmetry," (Springer,
1990), pp. 17–45.

# 1. Introduction

The implication of vascular disorders has been questioned for major diseases affecting the retina such as diabetic retinopathy [1], age-related macular denervation [2] and glaucoma [3]. Monitoring retinal blood flow appears crucial to understand the pathophysiology of these diseases. Most state of the art non-invasive instruments used for ocular blood flow assessment rely on speckle contrast analysis or laser Doppler methods. Both of these approaches use the intensity fluctuations of coherent light caused by moving objects in order to measure blood flow [4]. Nowadays, Optical Coherence Tomography Angiography (OCT-A) has emerged as a powerful vascular imaging technique [5]. OCT-A makes use of the speckle variations in the OCT signal caused by moving scatterers to generate an angiographic contrast by calculating the local variance of speckle realizations over a few repeated measurements [6, 7]. OCT-A instruments can map the retinal micro-vascular network with a micrometer axial resolution [8] and can be used to measure metrics relevant to the development of diabetic retinopathy such as the size and distribution of capillaries and the extent of the foveal avascular zone [9]. However, as it is necessary to reconstruct a full volume to obtain en-face images, the temporal resolution for swept-source OCT-A is of the order of a second, which is too slow for blood flow monitoring during a single cardiac cycle. Doppler OCT is a technique based on OCT that measures the Doppler frequency [10, 11] or phase [12] shift caused by the flow. The phase shift can be calculated between two consecutive measurements to extract a dynamic signal and reveal blood flow. However large flows causing a modulation higher than $2\pi$ rad cannot be unambiguously calculated [13]. Another limitation is that it cannot detect flow in vessels perpendicular to the incident OCT beam and only the axial component of the flow can be easily retrieved. Laser speckle flowgraphy is another technique based on speckle contrast whose purpose is to image retinal blood flow quantitatively, although in arbitrary units [14, 15]. En-face images of the retina are formed on the sensor and several images are combined to extract an angiographic contrast; this allows for full-field blood flow measurements with a high temporal resolution. Indeed a commercialized ophthalmic instrument based on this technology (LSFG NAVI, Softcare Ltd, Fukuoka, Japan) is able to image pulsatile flow in real-time with a temporal resolution of 33 ms. In contrast with OCT-A, laser speckle flowgraphy has no depth sectioning ability so the influence of choroidal flow is unclear. Inter-eye comparisons are not possible but intraindividual measurements are highly reproducible [15], which makes it an appropriate instrument to investigate ocular blood flow changes upon physiological or pharmacological stimuli. Laser Doppler flowmetry is based on measuring the Doppler broadening of coherent light induced by moving red blood cells in vascularized tissues [16–19]. The frequency shift is measured by means of interference between Doppler shifted light and non shifted light. This method allows monitoring of blood flow only in a single area, but it has been implemented on scanning systems to monitor blood flow over a full-field [20]. Blood flow and mean velocity in relatives units can be derived from calculations based on the first moments of the power spectrum density. As for laser speckle flowgraphy, interindividual comparisons are difficult whereas intraindividual measurements are reproducible [21]. Laser Doppler flowmetry was used in a direct image detection scheme [22] but this method is too slow for retinal imaging due ocular safety limits.

Lately, holographic detection schemes such as full-field swept-source OCT [23], off-axis full-field time-domain OCT [24] or off-axis holographic line-field en-face OCT [25] have emerged. These techniques allow parallelized imaging and have been demonstrated for in vivo human

<footnote n="42">42. J. Schmitt, A. Gandjbakhche, and R. Bonner, “Use of polarized light to discriminate short-path photons in a multiply scattering medium,” Appl. Opt. 31, 6535 (1992).</footnote>

<footnote n="43">43. L. F. Rojas-Ochoa, D. Lacoste, R. Lenke, P. Schurtenberger, and F. Scheffold, “Depolarization of backscattered linearly polarized light,” JOSA A 21, 1799–1804 (2004).</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Optical setup. L1, L2 and L3 are converging lenses. PBS: Polarizing Beam Splitter. BS: Beam Splitter. The light source is a single wavelength laser diode (SWL-7513-H-P, Newport). The 90% output of the fiber coupler is used for the object arm. The Doppler broadened light backscattered by the retina is combined with the reference field and interferograms are recorded on the two imaging channels to study the Doppler beat frequency. The data from the fast CMOS camera is processed offline while the 80 Hz camera is used for real-time monitoring.</figcaption></figure>

retinal imaging; they have been made possible by technological progress made by the camera industry which now allows for higher acquisition speed, and they seem to have opened the way for new possibilities especially in the realm of offline numerical processing [26, 27]. The laser Doppler holography technique we present here is conceptually close to laser Doppler flowmetry, but like the aforementioned methods it is derived from holography which allows for full-field imaging. In laser Doppler holography, the optical field serving as non-Doppler shifted light is a separate reference beam beating against the Doppler shifted optical field backscattered by the retina. This allows to increase as much as needed the power of the reference field impinging on the sensor and thus be able to work with very low exposure time. Interferograms are recorded with a high throughput camera and wideband measurements of the beat frequency of digitally reconstructed holograms are performed. The angiographic contrast is drawn from the Doppler spectral broadening of light backscattered by the retina and blood flow changes during cardiac cycles are revealed using a short-time Fourier transform analysis. We report here on the implementation of laser Doppler holographic perfusion imaging in the human retina. In contrast with previous attempts on rodent eyes [28], the holographic configuration is on-axis in order to maximize the field of view.
<figure><img src="image_3.png" /><figcaption>Fig. 2. Short-Time Fourier Transform analysis. A 3D sliding window of consecutive holograms is moved along the full hologram stack; for each window, the Doppler Power Spectral Density (DPSD) is estimated from the squared magnitude of the Fourier Transform in Eq. (3). A single image with blood flow contrast is obtained by integrating the DPSD over [ω₁, ω₂] in Eq. (4), and the resulting image M₀⁺ is referred to as power Doppler image. Examples of non-averaged power Doppler images are shown in Fig. 5, and examples of averaged power Doppler images are shown in Fig. 3(c) and Fig. 4(c).</figcaption></figure>

## 2. Methods

### 2.1. Optical setup

The optical setup designed and used for this study, shown in Fig. 1, consists of a fiber Mach-Zehnder optical interferometer. The light source used for the experiments is a 45 mW, single-mode, fiber diode laser (Newport SWL-7513-H-P) at wavelength λ = 785 nm. The retina is illuminated with a 1.5 mW constant exposure over 2.4 × 2.4 mm² area. Informed consent was obtained for each subject and experimental procedures adhered to the tenets of the Declaration of Helsinki. The laser beam is focused in the front focal plane of the eye so that the light is collimated inside the eye and illuminates the retina on an extended area. This irradiation level is compliant with the exposure levels of the international standard for ophthalmic instruments ISO 15004-2:2007.

A Polarizing Beam Splitter (PBS) cube is used in the object arm to illuminate the eye and collect the light backscattered by the retina. The PBS separates linear polarizations: it fully reflects the illumination beam and only the cross-polarized component of the light backscattered by the retina is transmitted to the camera. The reference beam ELO is collimated and considered to be monochromatic and stable. The object and reference waves are combined using a non-polarizing beam splitter cube and they interfere on the sensor plane. The polarization of the reference wave is adjusted with a half-wave plate and a polarizer to optimize fringe contrast.

The interference pattern is recorded simultaneously on two imaging channels: a slow CMOS camera (Ximea - xiQ, frame rate 80 Hz, 2048 × 2048 pixels, 8 bit pixel depth, pixel size 5.5 µm)
for real-time monitoring and a fast CMOS camera (Ametek - Phantom V2511, frame rate, up to 75 kHz for a 512 × 512 format, quantum efficiency 40%, 12-bit pixel depth, pixel size 28 µm) for offline processing.

## 2.2. Laser Doppler measurements

In laser Doppler measurement techniques, the Doppler effect is treated as a frequency shift of the optical field itself [29]; here the optical field E backscattered by the retina is considered to be the sum of the Doppler frequency shifted contributions. The interference pattern on the sensor is $I = |E_{LO} + E|^2$. In the cross-beating terms, the optical field frequencies cancel out: the interferogram beat frequency spectrum carries the optical Doppler broadening. Doppler frequency shifts for biological tissues are small compared to the laser optical frequency $3.8 \times 10^{14}$ Hz as the typical Doppler shift frequency spectra width for retinal vessels has been reported to be around 15-20 kHz [30]. However, as the spectral width corresponds to the Nyquist frequency of the fast camera, a major part of the Doppler broadening can be measured within the camera bandwidth in the interferogram beat frequency.

Acousto-optics modulators (not shown on the schematic) are used in conjunction with the slow camera to introduce a 40 Hz frequency shift in order to perform stroboscopic and narrowband signal measurements in time-averaging conditions [31]. However due to eye motion the quality of the signal is limited and this channel is mostly used for alignment purposes. As both methods are based on holography, the slow camera gives indications on the quality of the holographic reconstruction that are also valid for the fast camera; moreover it is able to yield a contrast related to blood flow that allows visualizing the retinal area being imaged during the acquisition, which is not possible on the fast camera. The real-time hologram rendering is done with our software (Holovibes) that runs calculations on graphics processing units; for this purpose a personal computer equipped with a Nvidia GeForce GTX Titan graphics card is used.

The fast camera is used for wideband measurements of the Doppler spectral broadening. Interferograms from this camera are recorded on its on-board memory and the data are processed offline using Matlab. Interferograms are recorded on the 512 × 512 pixel sensor with an acquisition rate of $f_S = 39$ kHz (except in section 4 where experiments are made to explore the effect of the sampling frequency), and a frame exposure time of 25.6 µs. Image rendering of complex-valued holograms $H(x, y)$ is performed using the angular spectrum propagation of the recorded interferograms $I(x, y)$, appropriate for near-field reconstructions [32]:

$$H(x, y, z) \approx \iint \hat{I}(k_x, k_y) \exp(ik_z z) \exp(ik_x x + ik_y y) dk_x dk_y \quad (1)$$

where $\hat{I}(k_x, k_y)$ is the two-dimensional Fourier transform of the interferogram $I(x, y)$. The optical axis is $z$, and $k_x, k_y$, and $k_z = \sqrt{k_x^2 - k_x^2 - k_y^2}$ are the wave vectors projections along $x, y$, and $z$, respectively. The lateral field of view of the reconstructed holograms of 512 × 512 pixels extends over $2.4 \times 2.4$ mm² on the retina.

## 2.3. Data processing

Data processing consists of measuring the local optical fluctuations of the holograms recorded with the high-speed CMOS camera. A Short-Time Fourier Transform (STFT) method is used to analyze the changes in blood flow over time through the changes in the Doppler broadened spectrum. The variations over time of the Doppler Power Spectrum Density (DPSD) are resolved with a sliding short time window constituted of $j_{win}$ consecutive holograms $H$. Given a sampling
<figure><img src="image_4.png" /><figcaption>Fig. 3. Optic Nerve Head (ONH) region imaged with commercialized instruments and Laser Doppler holography. (a) Scanning laser ophthalmoscope (Spectralis, Heidelberg). (b) Adaptive optics flood illumination (rtx1, Imagine Eyes). (c) Power Doppler images $M_{0+}$ calculated from holograms recorded at $f_S = 39$ kHz; $S(\omega)$ is integrated over $[f_1, f_2] = 4-19.5$ kHz. Multiple power Doppler images $M_{0+}$ are averaged over a total time of 80 ms (see Visualization 1 for blood flow movie). The peri-papillary crescent is visible on the edge of the ONH. (d) Asymmetry of the DPSD $M_{0-}$ (averaged over the same period of time) illustrating the resultant flow direction with respect to the optical axis.</figcaption></figure>

frequency $f_S$, the duration of the short-time window is thus:

$$ \sigma_{win} = \frac{j_{win}}{f_S} \qquad (2) $$

For each short-time window, a single image with a blood flow contrast is formed using the process illustrated in Fig. 2. First, the short-time window is fast Fourier transformed along the temporal dimension, then the power spectral density is estimated from the squared magnitude of the Fourier transform:

$$ S(x, y, t_n, \omega) = \left| \int_{t_n}^{t_n + \sigma_{win}} H(x, y, \tau) \exp(-i\omega\tau) d\tau \right|^2 \qquad (3) $$

Finally the quantity $M_{0\pm}$ is calculated from the integration of the DPSD:

$$ M_{0\pm}(x, y, t_n) = \int_{\omega_1}^{\omega_2} S(x, y, t_n, \omega) \pm S(x, y, t_n, -\omega) d\omega \qquad (4) $$

$M_{0+}$ is the zeroth moment of the bandpass filtered DPSD; it is a single image corresponding to the area under the curve of the DPSD over the frequency interval $[f_1, f_2] = (2\pi)^{-1} \times [\omega_1, \omega_2]$. $M_{0-}$ images reveal the asymmetry of the DPSD and are discussed in detail in section 3.4. Temporal
frequencies in this article are referred to indifferently as *f* or ω = 2π*f*. M₀⁺ images are referred to as Power Doppler images (PDIs) and reveal the local blood flow. The frequency band [f₁, f₂] is typically a high-pass filter (i.e. f₂ = f₁/2). The lower boundary f₁ is set between 4 kHz and 7 kHz in order to filter off the contributions that are not due to pulsatile flow in the vessels. The effect of the frequency band [f₁, f₂] is explored in subsection 4.1. Once the M₀⁺ image has been calculated, the sliding window is moved along the temporal dimension of the hologram stack. In this way, a sequence of PDIs revealing the blood flow for each moment of the cardiac cycle is made.

### 3. Results

#### 3.1. Blood flow contrast

Output images of the STFT analysis M₀⁺ are displayed in Fig. 3(c) and Fig. 4(c) and associated visualizations linked in corresponding captions. These images reveal the tissue perfusion occurring during the duration of the short-time window σwin. Indeed the Doppler frequency shift is greater when the velocity of the scattering particle is greater. When integrating the whole high-pass filtered DPSD, the pixels where the Doppler broadening is the greatest have the brightest intensity on M₀⁺ images. Thus areas of the retina where the local flow is greater are directly revealed in the intensity of PDIs.

As a result, the most noticeable features of the PDIs are the retinal vessels, due to the considerable flow of red blood cells moving at high speed. Blood vessels are clearly visible on PDIs and give a stronger signal than background tissues because of the larger flow in the vessels. Other features such as the nerve fiber layer are visible on reflectance images obtained with an adaptive optics flood illumination instrument (rtx1, Imagine Eyes, Orsay, France) in Fig. 3(b) and in Fig. 4(b), but these features are not visible on M₀⁺ images.

#### 3.2. STFT parameter σwin

The chosen width of the short-time window σwin results from a trade-off between temporal resolution and Signal-to-Noise Ratio (SNR). Integrating the signal over a longer time period improves the SNR insofar as the amplitude of lateral displacements stays limited during the duration of the short-time window. For excessively long short-time windows, involuntary eye movements limit the improvement in SNR. More importantly, as two events happening during the same short-time window cannot be separated, the sliding window duration σwin defines the temporal resolution of the instrument. In Fig. 5 and Fig. 6, a comparison is made to evaluate the importance of σwin on the SNR. The STFT analysis was carried out two times on the same dataset by using different parameters with an approximate tenfold increase in temporal resolution; jwin = 64 (σwin = 1.6 ms) in the first case and jwin = 512 (σwin = 13.1 ms) in the second case. As can be observed on individual PDIs in Fig. 5, the SNR greatly improves when increasing σwin. Nonetheless the vasculature is already visible with σwin = 1.6 ms. Besides when looking at pulsatile flow, the shape of the power Doppler signal with jwin = 64 can show cardiac cycle signal with a lower SNR than for jwin = 64.

For the STFT analysis conducted with σwin = 13.1 ms, a smaller ROI is required to spatially average the results to have a better SNR which illustrates as the SNR in the case where jwin = 64 is lower. Finally, it is interesting to notice that when averaging power Doppler images made with a STFT analysis performed with jwin = 64 and jwin = 512, the averaged image in the end are virtually the same as long as the total time on which PDIs are averaged is the same. This results are illustrated in images visible in Fig. 6: in this case, 800 PDIs with jwin = 64 and 100 PDIs with jwin = 512 over the same period of time, and the two averaged PDIs show the same features with an equal SNR.
<figure><img src="image_5.png" /><figcaption>Fig. 4. An artery intertwined with a vein is imaged with commercialized instruments and Laser Doppler holography. (a) Scanning laser ophthalmoscope (Spectralis, Heidelberg). (b) Adaptive optics flood illumination (rtx1, Imagine Eyes). (c) Power Doppler images $M_0^+$ calculated from holograms recorded at $f_S = 39$ kHz; $S(\omega)$ is integrated over $[f_1, f_2] = 7-19.5$ kHz. Multiple power Doppler images $M_0^+$ are averaged over a total time of 0.66 s (see Visualization 2 for blood flow movie). (d) Asymmetry of the DPSD $M_0^-$ (averaged over the same period of time) illustrating the resultant flow direction with respect to the optical axis.</figcaption></figure>

### 3.3. Pulsatile flow

The temporal resolution of the instrument is short enough to observe the changes in blood flow within cardiac cycles. When a heartbeat occurs, the increased speed and concentration of particles in vessels causes momentarily a larger Doppler spectral broadening of the optical field scattered in these structures, which is revealed by the STFT analysis. In Fig. 6, the changes in blood flow throughout cardiac cycles can be observed in a retinal artery and vein. The power Doppler signal is spatially averaged in regions of interest corresponding to an artery, a vein, and background tissues; PDIs have been registered to compensate for lateral eye movements using [33]. During systole (when the cardiac muscles contract and blood pressure is maximum) and diastole (when the cardiac muscles relax and the blood pressure is the lowest) the intensity of both vessels vary with the flow pulsatility. The artery gives a stronger signal than the vein during systole, because of the higher speed of circulating red blood cells which is in agreement with results reported in the literature [34]. The artery also shows a larger modulation depth of the flow than the vein.

Interestingly, for each cardiac cycle, a first peak is reached at the time where the perfusion is maximum, and is followed by a second one of lower amplitude. This second peak corresponds to the dicrotic notch which is thought to due to the stretch and recoil of the aorta [35]. This waveform has previously been observed in various blood flow monitoring systems such as laser Doppler spectroscopy [36], Doppler sonography [37], Doppler OCT [38], or adaptive optics scanning laser ophthalmoscopy [39]. It should be noted that the background tissue also exhibits a
<figure><img src="image_6.png" /><figcaption>Fig. 5. Comparison in Signal to Noise Ratio (SNR) of power Doppler images $M_{0+}$ and $M_{0-}$ calculated with $j_{win} = 64$ and 512 holograms (i.e. $\sigma_{win} = 1.6$ ms and 13.1 ms, respectively). The four images are non-averaged power Doppler images calculated from raw holograms acquired at 39 kHz: (a) $M_{0+}$, $j_{win} = 64$. (b) $M_{0+}$, $j_{win} = 512$. (c) $M_{0-}$, $j_{win} = 64$. (d) $M_{0-}$, $j_{win} = 512$. The SNR dramatically improves with the duration of the short-time window.</figcaption></figure>

flow behavior that we assume to be due to unresolved retinal or choroidal capillaries.

### 3.4. Laser Doppler spectral asymmetry

In the full-field CMOS Doppler imaging scheme proposed by Serov et al. [22], the light backscattered by the sample was not combined with a reference beam but was self-interfering. With this method, the non Doppler shifted light played the role of reference field to perform heterodyne measurements. In this configuration, a negative or positive Doppler frequency shift generates the same intensity fluctuations on the sensor: the signal recorded is real so its spectrum is exactly symmetric and the sign of the frequency shift cannot be recovered. In our case, as the reconstructed holograms are complex-valued and the phase of the optical field is known, the spectrum is no longer symmetric and it becomes possible to distinguish positive and negative frequency shifts. It can be observed that the difference between the negative and positive part of the DPSD contains information about the resultant flow direction when calculating the following quantity:

$$M_{0-}(x, y, t_n) = \int_{\omega_1}^{\omega_2} S(x, y, t_n, \omega) - S(x, y, t_n, -\omega) \, d\omega \quad (5)$$
<figure><img src="image_7.png" /><figcaption>Fig. 6. Comparison of pulsatile flow measurements in a vein and an artery for two STFT analysis performed on the same dataset with (a) $j_{win} = 64$ and (b) $j_{win} = 512$ holograms. For each analysis, a time averaged power Doppler images $M_{0+}$ is displayed on the left hand side and shows a red, blue and green ROI which mark an artery (A), a vein (V), and the background tissues (B), respectively. These regions are used to spatially average the power Doppler signal for the plots on the right hand side. Although noisier, the pulsatility signal is already well visible when using $j_{win} = 64$.</figcaption></figure>

The asymmetry of the spectrum is in agreement with the vessel geometry (cf. Fig. 3(d) and Fig. 4(d) considering the expected sign of the Doppler shift. Indeed the Doppler frequency shift $\omega_D$ occurring when the field is scattered by a moving object is calculated using the following formula:

$$ \omega_D = (\vec{k}_s - \vec{k}_i) \cdot \vec{v} \qquad (6) $$

where $\vec{v}$ is the velocity vector of the scattering particle and $\vec{k}_i$ and $\vec{k}_s$ are the incoming and scattered optical wave vectors, respectively. According to the direction of the particle with respect to the optical axis, the frequency shift can be positive or negative as schematized in Fig. 7(a) and Fig. 7(b). $M_{0+}$ thus reveals the predominant direction of the flow with respect to the optical axis. If the vessels are not perpendicular to the optical axis, the projection of the velocity vector on the optical axis is non-null and the sign of the DPSD asymmetry varies accordingly with vessel directions. This is well visible in Fig. 7(d) and Fig. 7(e). The power spectrum density has been calculated according to Eq.(3) in two regions that differ from the direction of the local flow. In the first case the vessel is going slightly upward with respect to the optical axis and $S(\omega) > S(-\omega)$ over the frequency band $[f_1, f_2]$ which is used to integrate DPSD when calculating the PDIs. In the second case, the vessel is going downward and $S(\omega) < S(-\omega)$ over the frequency band of interest. This discrepancy between negative and positive parts of the DPSD is imputable to the directly backscattered light. It can also be noticed that for lower frequencies (for example the peak occurring at about $\approx 1$ kHz), this inversion cannot be observed in two different regions as
<figure><img src="image_8.png" /><figcaption>Fig. 7. Cases of spectral asymmetry between the positive and negative parts of the DPSD calculated according to Eq. (3). (a, b) Schematic illustrating the sign of the Doppler frequency shift for the direct backscatter (cf. Eq. (6)). (c) $M_0$- showing the two regions where the spectra are calculated. (d) Negative and positive part of the DPSD measured in the ROI 1. (e) Negative and positive part of the DPSD measured in the ROI 2. For the range of frequencies 4-19.5 kHz used to calculate PDIs, $S_1(\omega) > S_1(-\omega)$ and on the contrary $S_2(\omega) < S_2(-\omega)$ depending on the vessel geometry.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Fig. 8. Multiple laser Doppler holography measurements of a same region are made with different sampling frequencies. For each of these sampling frequencies, PDIs $M_{0+}$ are calculated for different frequency bands indicated below the images. (a) $f_S$ = 10 kHz. (b) $f_S$ = 20 kHz. (c) $f_S$ = 39 kHz. (d) $f_S$ = 75 kHz.</figcaption></figure>

these frequencies correspond to Doppler shifts induced by global movements.

## 4. Role of the frequency band [$f_1, f_2$] and sampling frequency $f_S$

In this section, results from several measurements performed at the same retinal location using different sampling frequencies (10, 20, 39, 75 kHz) are exploited to investigate on the effect of the frequency band [$f_1, f_2$] and sampling frequency $f_S$ on the resulting power Doppler images $M_{0+}$. All measurements were done with the Phantom V2511 camera; the exposure time of the camera was set to the maximum for each acquisition (i.e. 100, 50, 25.6, 13.3 µs, respectively). The retina of a young and healthy subject was illuminated with the same exposure for all acquisitions, only the power of the reference arm was adjusted using neutral densities so that the camera was filled with the same amount of light (near saturation of the brightest pixels). The region imaged in these acquisitions is centered on the optic disc. The field of view is close to 3 × 3 mm², slightly larger than for the results shown in previous Figs. because we used a different combination of lenses.
## 4.1. Choice for [$f_1, f_2$]

The frequency band [$f_1, f_2$] used to integrate the DPSD in Eq. (4) is chosen so as to filter off the DC and autocorrelation terms as well as the signals due to global eye movements while keeping the contribution of Doppler shifted light due to pulsatile flow in the vessels. The lower boundary is generally set between 4 and 7 kHz depending on the retinal region being imaged while the upper boundary is set to the camera Nyquist frequency $f_s/2$. In Fig. 8, the displayed PDIs are calculated for different frequency bands for each sampling frequency. Each row of the Fig. corresponds to a sampling frequency, $f_s = 10, 20, 39, \text{and } 75$ kHz, for the rows (a), (b), (c), and (d), respectively. For each column, $M_{0+}$ are calculated for several bandpass interval: [$f_1, f_2$] = 1-3 kHz, 3-5 kHz, 5-10 kHz, 10-19.5 kHz and 19.5-37.5 kHz from left to right. Except for the acquisition with the highest sampling frequency, some intervals are above the Nyquist frequency of some sampling frequencies and consequently PDIs cannot be calculated for these frequencies. The visualization linked in caption Fig. 8 shows PDIs for a more precise range of frequencies in the case of the acquisition at 75 kHz.

It can be observed in Fig. 8 that for the lowest frequency band 1-3 kHz, the background tissue DPSD is greater than the vessels DPSD, so the vessels appear darker than the tissue on the PDIs. For the frequency band 3-5 kHz, the edges of the large retinal vessel appear brighter than the lumen, as the latter remains dark. That is because the flow is greater at the center of the vessels and the Doppler broadening occurring in these structures is still undersampled. However, for this frequency band, the Doppler broadening caused by the flow in smaller retinal vessels seems adequately sampled as they appear brightly. For frequencies superior to 10 kHz, the DPSD in the lumen seems adequately sampled. Thus integrating lower frequencies allows us to get the signal from smaller retinal vessels, although it also increases the contribution of global axial and lateral movements. It can be noticed that PDIs calculated using a certain frequency band look alike independently from the sampling frequency. For example PDIs calculated with the frequency band 1-3 kHz are virtually the same with the sampling frequencies 10 kHz and 75 kHz. So in the former case, all the signal above 5 kHz (Nyquist frequency for $f_s = 10$ kHz) which is aliased has no visible effect on the resulting $M_{0+}$ image.

This leads to the conclusion that to make a simple image of the vasculature, the aliased part of the DPSD does not make significant change on PDIs. It is also possible to conclude that the DPSD corresponding to the contributions of pulsatile flow in large vessels lies in frequency bands that are above 5 – 10 kHz and up to 37.5 kHz depending on the size of the vessel.

## 4.2. Sampling frequency $f_s$

In Fig. 9, the requirements in terms of sampling frequency to perform laser Doppler holography in the human retina are investigated. For this purpose, both the vasculature images of the retina $M_{0+}$ and the changes in power Doppler signal during the cardiac cycles are considered. STFT analysis are performed with $j_{win} = 128, 256, 512, 1024$ holograms for the acquisitions at $f_s = 10, 20, 39, \text{and } 75$ kHz respectively. In this way the temporal resolutions of the STFT analyses are approximately the same for all acquisitions ($\sigma_{win} \approx 13$ ms). Power Doppler images $M_{0+}$ are calculated using the highest possible frequency bands according to the sampling frequency: [$f_1, f_2$] = 3-5 kHz, 5-10 kHz, 10-19.5 and 19.5-37.5 kHz, for $f_s = 10, 20, 39, \text{and } 75$ kHz, respectively. For each sampling frequency, an image calculated from PDIs averaged over time is displayed on the left, and the plots on the right correspond to the signal over time of $M_{0+}$ averaged in the region indicated on the associated image. This ROI is centered on the arterial bifurcation of the central retinal artery where the flow is expected to be the greatest in the retina.

As observed in the previous subsection 4.1, PDIs of the vasculature at lower frequencies show a weaker signal inside the vessel lumen because the Doppler broadening is undersampled. The curve for the lowest frequency band 3-5 kHz displayed on the right side on the Fig. does not show a pattern related to the cardiac cycle. However for the frequency band above (5-10 kHz), the
<figure><img src="image_10.png" /><figcaption>Fig. 9. Requirements in terms of sampling frequency for laser Doppler holography in the central retinal artery. For each sampling frequency, power Doppler images are spatially averaged over the depicted ROIs on the left hand side and the result is displayed in the associated plot. The STFT parameters used for each acquisition are displayed on the left and have been chosen to have $\sigma_{win} \approx 13$ ms for all acquisitions. (a) $f_S = 10$ kHz, (b) $f_S = 20$ kHz, (c) $f_S = 39$ kHz, (d) $f_S = 75$ kHz. Visualization 3 shows PDIs as a function of the frequency for the acquisition with $f_S = 75$ kHz.</figcaption></figure>
typical curve of flow over time in the vessels is recognizable. Then for the two highest frequency bands, the difference is minor. Indeed when looking closely it is possible to see that around the maximum of the systolic peak for the 10-19.5 kHz plot, the curve is not well rounded and shows an abrupt decrease (indicated by the arrows) that is not visible for 19.5-37.5 kHz plots. We assume that a minor undersampling briefly occurs at this moment for large systolic flows. In Fig. 9(d), an overall increase over time of the power Doppler signal can be noticed: it is due to an involuntary eye movement which progressively shifts the pupil position.

Hence Fig. 9 shows that the Doppler shifted light corresponding to the pulsatile flow in large vessels is located further along the power spectrum density. Considering these results, performing measurements with a sampling frequency above 30 to 40 kHz seems to be a requirement for large retinal vessels. Otherwise, the part of the DPSD corresponding to the pulsatile flow is undersampled, aliased, and mixed with low frequencies of greater amplitude and is therefore not retrievable. However, for vessels of smaller size it is possible to use a smaller sampling frequency and have a properly sampled Doppler signal but there is a risk that this signal might be too weak due to the low number of circulating red cells especially compared to the Doppler shifts induced by global eye movements.

## 5. Discussion and conclusion

It is possible to record up to 259,329 frames in a 512 × 512 format on the on-board memory of the Phantom V2511 camera. This corresponds to approximately 6.6 s of acquisition at $f_S = 39$ kHz or 3.5 s of acquisition at $f_S = 75$ kHz. We showed in the previous section that with $f_S = 39$ kHz the Doppler broadening caused by pulsatile flow seems somewhat undersampled only when imaging the central retinal artery during systole. As a result, when imaging other retinal regions where the flow is weaker, this sampling frequency yields negligible aliasing. Consequently, for this work the sampling frequency was usually set to 39 kHz to be able to record longer acquisitions and maximize our chances to measure consecutive cardiac cycles unaffected by micro-saccades. Indeed, as with most ophthalmic instruments, the ability of laser Doppler holography to perform retinal imaging during micro-saccades is limited. The Doppler broadening induced by global movements is greater than the one induced by local pulsatile flows (for example the strong peak occurring at 2.8 s in Fig. 6 is a micro-saccade). Consequently, during a micro-saccade, the intensity of the PDIs is predominated by global movements. However drifts and tremor do not prevent blood flow imaging, as can be observed in movies because slow movements generate low frequency Doppler shifts that are removed when high pass filtering the DPSD. It should be noted that this process also filters off the Doppler signal due to the flow in small retinal vessels.

The technique presented in this article makes use of the backscattered light that has been highly Doppler shifted. This part of the light is selected when high-pass filtering the DPSD in Eq.(3). One of the assumptions of the theory used for laser Doppler flowmetry is that photons scattered by red blood cells (i.e. highly Doppler shifted photons) have been scattered by static diffusers before they impinge on erythrocytes [40, 41]. The result of this is a randomization of the wave vectors involved in the scalar product in Eq. (6) yielding the Doppler frequency shift: the situation is comparable to red blood cells being illuminated from all directions. Thus, selecting highly Doppler shifted light implies selecting multiply scattered light. Yet multiply scattered light loses the original linear polarization of the incident beam because light is depolarized in the multiple scattering regime [42, 43]. Consequently, it is possible to carry out laser Doppler holography with cross-polarized light because we work with the backscattered light that is highly Doppler shifted, thus multiply scattered and depolarized. For this reason we used a polarizing beam splitter cube which allows us to filter off the specular reflection occurring on the cornea. Indeed this reflection conserves the original linear polarization of the incident beam and is reflected by the polarizing beam splitter cube instead of being transmitted to the camera.

An important aspect of working with multiply scattered light is that the instrument is sensitive to
lateral motion. For conventional blood flow monitoring methods based on Doppler approaches, the flow is assessed in a simple manner only in vessels parallel to the optical axis. The flow in vessels perpendicular to the optical axis cannot be assessed as the Doppler shift is minimum, and for vessels with intermediate inclinations, the geometry must be known to calculate the corresponding flow. In our case, because we select multiply scattered light, vessels with any inclinations with respect to the optical axis contribute to the Doppler broadening of the backscattered light. Consequently the instrument is sensitive to lateral motion and able to image blood vessels in en-face planes as can be seen in the PDIs shown in Figs. 3, 4, and 8.

However in laser Doppler flowmetry, the retinal area where the blood flow is monitored is not a large blood vessel but tissues that are mainly vascularized by capillaries. In this medium the randomization of the light is a safer assumption because small vessels are oriented in all directions and are surrounded by a dense matrix of static diffusers. However in our case it can be argued that the multiple scattering assumption is less valid as the vessel walls and other structures around large vessels may not diffuse the incident beam sufficiently, and the limits of this assumption can be seen in the asymmetry of the DPSD $M_0$. Indeed in a situation of purely multiple scattered light there should not be any difference between the positive and negative parts of the spectrum despite the vessels geometry since the asymmetry exists because of the signed Doppler shifts of the directly backscattered light. However, as explained in section 3.4, a mild asymmetry does exist in the DPSD depending on the vessel geometry (less than 10 to 15% in the example shown). This asymmetry may cause a bias of flow sensitivity according to the vessels orientation with respect to the optical axis. However the impact seems limited and is mitigated by the process of averaging the positive and negative parts of the spectrum. Moreover in the parts of vessels where the flow direction is perpendicular to the optical axis (i.e. the areas where $M_0$ = 0), the flow sensitivity does not seem to be affected by the local inclination of the vessel with respect to the optical axis.

In conclusion, wideband laser Doppler holography can reveal retinal blood flow pulsatility in humans over a full-field with a temporal resolution of only 1.6 ms while keeping the retinal exposure under 1.5 mW at 785 nm. The angiographic contrast which is sensitive to local dynamics, including lateral motion, is drawn from the Doppler power spectrum density using a short-time Fourier transform analysis. Thanks to its high temporal resolution, laser Doppler holography could provide new insight in the dynamics of the vascular system as well as potential new ways to diagnose and follow-up major retinal or cardio-vascular diseases.

## Funding

LABEX WIFI (Laboratory of Excellence ANR-10-LABX-24) within the French Program Investments for the Future under Reference ANR-10-IDEX-0001-02 PSL; and European Research Council (ERC Synergy HELMHOLTZ, grant agreement #610110).

## Acknowledgments

The authors would like to thank S. Meimon, K. Grieve and J-P. Huignard for helpful discussions. The Titan Xp used for this research was donated by the NVIDIA Corporation.

## Disclosures

The authors declare that there are no conflicts of interest related to this article.