4D holographic microscopy of zebrafish larvae microcirculation

An original technique that combines digital holography, dual illumination of the sample and cleaning algorithm 3D reconstruction is proposed. It uses a standard transmission microscopy setup coupled with a digital holography detection. The technique is 4D, since it allows to determine, at each time step, the 3D locations (x,y,z) of many moving objects that scatter the dual illumination beam. The technique has been validated by imaging the microcirculation of blood in a fish larvae sample (the moving objects are thus red blood cells RBCs). Videos showing in 4D the moving RBCs superimposed with the perfused blood vessels are obtained. © 2016 Optical Society of America OCIS codes: (090.1995) Digital holography; (090.2880) Holographic interferometry; (180.3170) Interference microscopy; (290.5850) Scattering, particles. References and links 1. E. Friedman, S. Krupsky, A. Lane, S. Oak, E. Friedman, K. Egan, and E. Gragoudas, “Ocular blood flow velocity in age-related macular degeneration,” Ophthalmology 102, 640–646 (1995). 2. M. P. Pase, N. A. Grima, C. K. Stough, A. Scholey, and A. Pipingas, “Cardiovascular disease risk and cerebral blood flow velocity,” Stroke 43, 2803–2805 (2012). 3. J. Kur, E. A. Newman, and T. Chan-Ling, “Cellular and physiological mechanisms underlying blood flow regulation in the retina and choroid in health and disease,” Progress in retinal and eye research 31,377–406 (2012). 4. O. Sakurada, C. Kennedy, J. Jehle, J. Brown, G. L. Carbin, and L. Sokoloff, “Measurement of local cerebral blood flow with iodo [14c] antipyrine,” Am. J. Physiology-Heart and Circulatory Physiology 234, H59–H66 (1978). 5. I. Kanno, H. Iida, S. Miura, M. Murakami, K. Takahashi, H. Sasaki, A. Inugami, F. Shishido, and K. Uemura, “A system for cerebral blood flow measurement using an h215o autoradiographic method and positron emission tomography,” J. Cerebral Blood Flow & Metabolism 7, 143–153 (1987). 6. Y. Yeh and H. Cummins, “Localized fluid flow measurements with an he-ne laser spectrometer,” Appl. Phys. Lett. 4, 176–178 (1964). 7. J. D. Briers and S. Webster, “Laser speckle contrast analysis (lasca): a nonscanning, full-field technique for monitoring capillary blood flow,” J. Biomed. Opt. 1, 174–179 (1996). 8. A. K. Dunn, “Laser speckle contrast imaging of cerebral blood flow,” An. Biomedical Engineering 40, 367–377 (2012). 9. D. Briers, D. D. Duncan, E. Hirst, S. J. Kirkpatrick, M. Larsson, W. Steenbergen, T. Stromberg, and O. B. Thompson, “Laser speckle contrast imaging: theoretical and practical limitations,” J. Biomed. Opt. 18, 066018– 066018 (2013). 10. S. Yuan, A. Devor, D. A. Boas, and A. K. Dunn, “Determination of optimal exposure time for imaging of blood flow changes with laser speckle contrast imaging,” Appl. Opt. 44, 1823–1830 (2005). 11. Y. Zeng, M. Wang, G. Feng, X. Liang, and G. Yang, “Laser speckle imaging based on intensity fluctuation modulation,” Opt. Lett. 38, 1313–1315 (2013). 12. D. Gabor, “A new microscopic principle,” Nature 161, 777–778 (1948). 13. E. N. Leith and J. Upatnieks, “Reconstructed wavefronts and communication theory,” J. Opt. Soc. Am. A 52, 1123–1128 (1962). 14. U. Schnars and W. Jüptner, “Direct recording of holograms by a ccd target and numerical reconstruction,” Appl. Opt. 33, 179–181 (1994). 15. M. Atlan, M. Gross, P. Desbiolles, É. Absil, G. Tessier, and M. Coppey-Moisan, “Heterodyne holographic microscopy of gold particles,” Opt. lett. 33, 500–502 (2008). 16. S.-H. Lee, Y. Roichman, G.-R. Yi, S.-H. Kim, S.-M. Yang, A. van Blaaderen, P. van Oostrum, and D. G. Grier, “Characterizing and tracking single colloidal particles with video holographic microscopy,” Opt. Express 15, 18275–18282 (2007). 17. F. C. Cheong, B. J. Krishnatreya, and D. G. Grier, “Strategies for three-dimensional particle tracking with holographic video microscopy,” Opt. Express 18, 13563–13573 (2010). 18. J. Gao, J. A. Lyon, D. P. Szeto, and J. Chen, “In vivo imaging and quantitative analysis of zebrafish embryos by digital holographic microscopy,” Biomed. Opt. Express 3, 2623–2635 (2012). 19. N. Verrier, D. Alexandre, and M. Gross, “Laser doppler holographic microscopy in transmission: application to fish embryo imaging,” Opt. Express 22, 9368–9379 (2014). 20. F. Saglimbeni, S. Bianchi, A. Lepore, and R. Di Leonardo, ”Three-axis digital holographic microscopy for high speed volumetric imaging,” Opt. Express 22, 13710–13718 (2014). 21. J. Högbom, “Aperture synthesis with a non-regular distribution of interferometer baselines,” Astronomy and Astrophysics Supplement Series 15, 417 (1974). 22. F. Soulez, L. Denis, C. Fournier, É. Thiébaut, and C. Goepfert, “Inverse-problem approach for particle digital holography: accurate location based on local optimization,” J. Opt. Soc. Am. A 24, 1164–1171 (2007). 23. F. Soulez, L. Denis, E. Thiébaut, C. Fournier, and C. Goepfert, “Inverse problem approach in particle digital holography: out-of-field particle detection made possible,” J. Opt. Soc. Am. A 24, 3708–3716 (2007). 24. F. Le Clerc, L. Collot, and M. Gross, “Numerical heterodyne holography with two-dimensional photodetector arrays,” Opt. Lett. 25, 716–718 (2000). 25. M. Westerfield, “The Zebrafish Book: a Guide for the Laboratory use of Zebrafish (Danio rerio)” (Institute of Neuroscience. University of Oregon, 1995). 26. S. Isogai, M. Horiguchi, and B. M. Weinstein, “The vascular anatomy of the developing zebrafish: an atlas of embryonic and early larval development,” Developmental biology 230, 278–301 (2001). 27. N. Warnasooriya, F. Joud, P. Bun, G. Tessier, M. Coppey-Moisan, P. Desbiolles, M. Atlan, M. Abboud, and M. Gross, “Imaging gold nanoparticles in living cell environments using heterodyne digital holographic microscopy.” Opt. Express 18, 3264–73 (2010). 28. F. Verpillat, F. Joud, P. Desbiolles, and M. Gross, “Dark-field digital holographic microscopy for 3d-tracking of gold nanoparticles,” Opt. Express 19, 26044–26055 (2011). 29. N. Verrier, D. Alexandre, G. Tessier, and M. Gross, “Holographic microscopy reconstruction in both object and image half-spaces with an undistorted three-dimensional grid,” Appl. Opt. 54, 4672–4677 (2015). 30. E. Cuche, P. Marquet, and C. Depeursinge, “Spatial filtering for zero-order and twin-image elimination in digital off-axis holography,” Appl. Opt. 39, 4070–4075 (2000). 31. L. Yu and M. K. Kim, “Wavelength-scanning digital interference holography for tomographic three-dimensional imaging by use of the angular spectrum method,” Opt. lett. 30, 2092–2094 (2005).


Introduction
Blood flow i maging techniques are widely used in bio-medical studies, since they can assess physiological processes or can be used for early detection of disease [1][2][3]. However, many blood flow studies require, for imaging purposes, the use of a contrast agent, making blood flow characterization invasive [4,5]. Scanning Doppler imaging techniques can be considered to alleviate this issue, but due to the scanning step, acquisition of an image is a time consuming process [6]. By analyzing the spatial statistics of the dynamic speckle with the Laser Speckle Contrast Analysis/Imaging (LASCA/LSCI) technique [7] the blood flow can be imaged [8,9]. Improvement of the acquired contrast image has been achieved through exposure time optimization [10], or intensity fluctuation analysis [11] resulting in high quality perfusion images. However, these techniques are limited to perfusion monitoring in 2D and unable to reconstruct the structure of the blood flow in 3D.
Proposed by Gabor [12] and improved by Leith and Upatnieks [13] holography records the interference of the field scattered by an object with a known reference field. It is then possible to extract, from this interference pattern, the scattered field that reaches the holographic detector. This is particularly easy in digital holography that uses a camera detector. Holography is intrinsically a 3D technique since the Maxwell equations can be used to back propagate the field from the camera to any point of the 3D space [14]. This allows to localize and track a point-like object that scatters light, since its position coincides with the maximum of intensity of the reconstructed field [15]. The localization accuracy along the z axis depends strongly on the experimental conditions, in particular on the detection numerical aperture. Nanometric accuracy can be reached in holographic microscopy with a high numerical aperture (NA ≥ 1) objective [16,17].
Digital holography has been used to image and qualitatively analyze blood flows in 2D [18,19]. 3D holographic imaging of blood flow is challenging for a number of reasons. The biologically relevant size of the microcirculation makes it impossible to use a high NA objective, since the field of view would be too small. Blood cells are much larger than the wavelength, and their refractive index is close to the plasma one. The light is then scattered within a small angle in the forward direction. In transmission geometry, this corresponds to a small NA and to a low z accuracy. In reflection geometry, this leads to a low signal, which competes with light scattered by the surrounding living tissues, whose refractive index is not homogeneous in time and space.
In this paper, we have imaged in 3D moving red blood cells (RBCs) in a living fish larvae. We used an original technique that combines digital holography, illumination of the sample and reconstruction along several axis [20], and calculations that involve both standard holographic propagation and 3D reconstruction by a cleaning algorithm [21-23]. More specifically, the experiment was developed with an holographic microscopy setup working in transmission. To increase the angle diversity of the light scattered by the sample, the illumination was made by two beams that were angularly tilted with respect to the optical axis. At each time point, the difference of successive images was calculated yielding 2-frame holograms that contain information on the light that is scattered by the moving objects (i.e. the RBCs) at that time. Two holograms, that describe the field components of each illumination beam, were then extracted from the data. These two holograms were used to calculate the two scattered fields 3D maps corresponding to each illumination. Since these two maps describe the scattering made by the same RBCs, they must exhibit correlation maximums in locations (x, y, z) that correspond to RBCs. We used thus these correlation maximums to calculate the 3D map of the RBC locations. We got then 3D videos of the moving RBCs and, by averaging, 3D images of perfused blood vessels.

Experimental setup
Our holography set-up is similar to the one of the reference [19] (see Fig. 1). The set-up uses an upright microscope (Olympus® CX41) that has been modified to perform heterodyne holography [24] in transmission. The main laser (HL6545MG: 60 mW @ λ = 660 nm, optical frequency ω I ) is split into two arms (illumination and reference) by a beam splitter BS. Two acousto-optic modulators (AOM) at ω 1,2 /2π ≃ 80 MHz control the frequency ω LO of the reference (i.e. local oscillator) field E LO : The illumination arm (field E) is split into two branches so as to illuminate the sample in two directions, whose relative angle is 28.3 • . The sample is a zebrafish larvae (150 hours-old) between slide and cover slip. It is imaged with a microscope objective MO (NA= 0.30, G=10). To achieve off-axis holography, the beam splitter that recombines the signal (E) and reference (E LO ) fields is angularly tilted: (θ = 0 . Larvae were mounted on a standard glass slide under a 22 mm 1.5 round coverslip. A 0.5 mm thick caoutchouc ring served as a spacer and sealant. They were embedded in a drop of cooling 1.5% low melting point agarose at 37 • C and oriented before gelation. The chamber was filled with 100 mg/l tricaine in filtered tank water and imaging was performed at room temperature. Blood vessels nomenclature is according to [26].
3. Data filtering and reconstruction.

4 phase reconstruction of the MO pupil image and selection of the signal corresponding to each illumination direction
We first performed a test experience, in which the field scattered by the sample is detected by 4 phases heterodyne holography. The detection is then made at the illumination frequency, and the holographic signal corresponds to the field scattered by immobile or almost immobile scatterers of the sample (ground glass or zebrafish larvae). The resulting signal is strong and allows for precise optical adjustment of our holographic device. This test experiment is also useful to illustrate the reconstruction procedure. This procedure is similar to the one used in previous works [27, 28] (for details see [29]), but with small differences due to the double illumination, and to the need of selecting the signals corresponding to each illumination. In the test experiment, we tuned the LO frequency ω LO to perform 4 phase detection, and we considered 4 phase holograms H 4ϕ . We have thus: where I 0 , I 1 ...I 3 are consecutive camera frames, and j 2 = −1. To select the +1 holographic grating order and to compensate for the MO pupil phase curvature, we performed the holographic reconstruction of the MO pupil by the Schnars et al. method [14] yieldingH 4ϕ (k x , k y ) with: where FFT is the discrete Fourier transform operator, and k is the modulus of the wave-vector: k = 2π/λ . Here, e jk(x 2 +y 2 )/2d is the quadratic phase factor that determines the reconstruction plane. It depends on the curvature of the reference wave front, and on the distance from the camera to the reconstruction plane. Note that a phase quadratic in k x and k y that affects H 4ϕ is missing in Eq.
(2) (with respect to the Schnars et al. work [14]). This missing phase does not affect the intensity images |H 4ϕ (k x , k y )| 2 .

Figures 2(a) and 2(b)
show the reconstructed images of the pupil. The display is made in arbitrary log scale for intensity: |H 4ϕ | 2 . Images are obtained with a ground glass ( Fig. 2(a)) and with a living zebrafish sample ( Fig. 2(b)). Figure 2(a) exhibits a circular zone whose brightness is high and homogeneous. This zone corresponds to the image of the MO pupil that is back illuminated through the ground glass. The brightness is homogeneous because the ground glass scatters light over angles wider than the MO collection angle. Because of the off-axis configuration (angular tilt θ of BS), the pupil image is located in the upper right side of the calculation grid, but it is displayed in the center of Fig. 2(a). Note that diffusers do not move in the ground glass experiment. It results that the scattered field E does not vary with time, and that the 4-phase demodulation made by Eq. (2) selects the +1 grating order term, since we have: The -1 and 0 grating order images, that are zero, are not visible in Fig. 2(a). To get the images of Fig. 2, the reconstruction parameter d (used to calculate the phase factor e jk(x 2 +y 2 )/2d of Eq. (3) ) has been adjusted to perform the reconstruction in the MO pupil plane. As seen, the pupil image exhibits sharp edges (white dashed line circle). The useful holographic information was then selected by cropping the pupil zone and by filling the remaining of the calculation grid with zeros. This corresponds to the Cuche et al. spatial filtering of the +1 grating order term EE * LO [30]. To compensate for the off-axis angular tilt, the cropped zone was translated into the center of the calculation grid. Figure 2(b) was obtained by imaging a zebrafish sample. Because the sample was quite transparent, the two illumination beams yield two bright spots, which are well separated in the Fourier space (see arrows 1 and 2 in Fig. 2(b) ). Around the spots, we can see two blurred brighter regions corresponding to the light scattered by the sample. Since these brighter regions are well separated, we extracted, from the holographic dataH 4ϕ , two hologramsH purple and H green that describe the field scattered from each illumination direction. This was done by cropping, withinH 4ϕ , two circular zones of radius 240 pixels centered on each spot (dotted white circle in Fig. 2(c) ) yielding the two hologramsH purple andH green (see purple and green zones in Fig. 2(c) ). Note that all the scattered light has been caught by eitherH purple orH green . The selection of the illumination direction was thus done without loss of information.

Holographic reconstruction of the field scattered by the circulating RBCs
In order to select the signal from moving RBCs, we considered 2-phase holograms H 2ϕ recorded without frequency shift of the LO beam.
The signal resulting from an object that is not moving is not taken into account in H 2ϕ . The images that are reconstructed from H 2ϕ correspond then to the moving components of the sample, i.e. to RBCs.H 2ϕ is calculated similarly toH 4ϕ by: Because H 2ϕ (x, y) is real,H 2ϕ exhibits both the +1 and -1 order images. Nevertheless, because of the off-axis configuration, these two images do not overlap. The +1 image is thus selected by cropping the pupil zone (whose exact location has been determined previously by the ground glass experiment), and by translating this zone into the center of the calculation grid. By this way, the useful holographic information is selected and the off-axis angular tilt is compensated. To separate the signals of the two illuminations, we performed, like in section 3.1, two crops of radius 240 pixels centered on the two illumination spots of Fig. 2(c). We obtained the two Fourier space hologramsH purple (k x , k y ) andH green (k x , k y ), from which we have calculated the two holograms H purple (x, y) and H green (x, y): where FFT −1 is the reverse Fourier transform operator. The holograms H purple (x, y, z) and H green (x, y, z), reconstructed in planes z, are then calculated by the angular spectrum method [31]. This method involves two Fourier transforms from the real space holograms H purple (x, y) and H green (x, y), but only one Fourier transform from the Fourier space hologramsH purple (k x , k y ) andH green (k x , k y ). We have: H green (x, y, z) = FFT −1 e − j(k 2 x +k 2 y )z/2kH where FFT is the discrete Fourier transform operator. Here, x, y and k x , k y are discrete coordinates whose pitches are ∆x and ∆k, which obey to N ∆x ∆k = 2π where N = 1024 is the size of the calculation grid.
Up to now, we have considered that both acquisition and reconstruction were made in the image half-space (near the camera). The pitches are then ∆x = D pix and ∆k = 2π/(ND pix ) (where D pix = 14 µm is the size of the camera pixels). x, y and z are the coordinates of the image of the sample (that is conjugated with the sample by MO). This point of view, which is most often adopted in the literature, is not very convenient, since we want to reconstruct the scattered fields with coordinates related to the sample.
A better point of view [29] is to consider that both acquisition and reconstruction are made in the object half-space (i.e. near the sample). H 2ϕ (x, y) is then the hologram recorded in the plane of the object half-space that is conjugated with the camera by MO. All the mathematical transformations made to calculate H purple (x, y) and H green (x, y) yield holograms, whose phases are properly corrected (i.e. the phase of the hologram is the same as the phase of the field). Equation (7) performs then the reconstruction with orthogonal coordinates x, y and z that correspond to the real coordinates near the sample. The pitches are ∆x = D pix /G, ∆k = 2πG/(ND pix ), where G is the imaging gain from the plane z = 0 to the camera. Note that G is generally not equal to the objective nominal gain (×10 here), since it depends on the exact location of the camera with respect to MO. The best is to measure the gain G by a calibration procedure.
In analyzing our data, we have considered the second point of view. We have measured G by imaging a USAF target. This calibration yields G = 24.6, ∆x = 14 µm/24.6 = 0.569 µm and ∆k = 2π/(1024∆x). and H green,m (x, y, z) of the images were obtained by selecting the purple and the green zones of Fig. 2(c). Reconstruction was made with z = −26.7 (a), z = 0 (b) and z = +26.7 µm (c). The displayed images (a..c) correspond to image 77 of file Visualization 1 (a), Visualization 2 (b) and Visualization 3 (c). Dorsal side left, anterior to the top. Bar is 100 µm.

Dual illumination reconstructed images and videos
We imaged a zebrafish sample by recording, without frequency shift of the LO beam (ω LO = ω I ), a sequence of frames I m (with m = 1...129). With each couple of consecutive frames: I m and I m+1 , we calculated the 2-phase holograms H 2ϕ,m = I m − I m+1 . With each hologram H 2ϕ,m , we calculated H purple,m (x, y, z) and H green,m (x, y, z) by Eqs. (5)- (7). The video files Visualization 1, Visualization 2 and Visualization 3 show the intensity reconstructed images, whose components are |H purple,m (x, y, z)| 2 and |H green,m (x, y, z)| 2 , that were obtained, for m = 1 to 128, and for respectively z = 0, z = +26.7 and z = 53.5 µm.
The images 77 (time index m = 77) of these video files are displayed on Figs. 3(a)-3(c). In figure 3 (a), the reconstructed plane z = 0 corresponds to the locations of the caudal vein and caudal artery. Indeed, the purple and green components of the caudal vein and caudal artery images are superimposed yielding white patterns. In figure 3(b), with z = 26.7 µm, the large parallel longitudinal vessels (caudal vein and caudal artery) are out of the reconstruction plane, and the purple and green components are well separated. On the other hand, the dorsal capillary seen in the bottom left side of the image is within the reconstruction plane, and the purple and green images of the RBC pointed by the arrow coincides (see zoom in Fig. 3(b)). In figure 3 (b), with z = 53.5 µm, the purple and green components of both the caudal vein, caudal artery and dorsal capillaries are separated. To image the inside structure of the perfused blood vessels, corresponding to the entire trajectories of the RBCs along the observation time, we calculated the averaged intensity holograms |H purple | 2 and |H green | 2 : The video file Visualization 4 shows the averaged intensity reconstructed images, whose components are |H purple (x, y, z)| 2 and |H green (x, y, z)| 2 . These images are obtained for z = −53.5 to z = +53.5 µm in 101 steps of 1.07 µm. Figure 4 shows the averaged intensity reconstructed images obtained for z = 0 (a), z = +26.7 (b) and z = +53.5 µm (c).  5. Scheme of the 3D reconstruction process by the cleaning algorithm. In (a) the product between the two 3D grids is calculated. This allows to select, in (b), the point of highest correlation x 1 ; y 1 ; z 1 . Thus, in (c), the point is stored in S purple (x, y, z) and S green (x, y, z) and erased from |H purple (x, y, z 1 )| and |H green (x, y, z 1 )|. From this modified plane, in (d) the overall grid is recalculated. In (e) a new maximum correlation point x 2 ; y 2 ; z 2 is found. Again, in (f), new sources are stored in the associated space of sources and set to zero in the plane z = z 2 of the 3D grid. A new cycle will start and the operation will be repeated K times. the XY , Y Z and XZ cuts of the averaged reconstructed 3D maps obtained without and with the cleaning algorithm i.e. |H purple (x, y, z)| 2 and |H green (x, y, z)| 2 (without) and |S purple (x, y, z)| 2 and |S green (x, y, z)| 2 (with).
Without cleaning (see Fig. 6 (a) ), the z confinement is poor. In the XZ and Y Z cuts, the green and purple images of the blood vessels are angularly tilted, and the location of the vessels corresponds to the z location where the green and purple images coincide. With cleaning (see Fig. 6 (b) ), the z confinement is much better. Since the green and purple images coincide, the images are displayed in white. The 3D maps obtained without and with cleaning are also displayed in the video files Visualization 5 and Visualization 6 (x, y and z are swept in 128 steps of pitch 8∆x and ∆z). To image a quite large zone of the zebrafish blood vessel structure, the experiment was made with a NA=0.3 microscope objective. The purple versus green axis angle is thus lower than in [20] (see xz cut of Fig. 6 (a) ), and the z confinement is thus lower too.
Note that the RBCs are much larger than one voxel, but, to keep all useful information, the cleaning algorithm was made by considering one voxel scatterers. K = 25000 corresponds thus to the number of voxels that is considered in the cleaning algorithm. This figure is enough to handle most of the scattered energy, and thus to image all the RBCs, whose number is much lower than K.

Visualization of the data calculated by the cleaning algorithm
The visualization of the 4D maps S purple and S green calculated by the cleaning algorithm is challenging, because of the 4 coordinates. To perform visualization, we have projected the instantaneous intensities |S purple | 2 and |S green | 2 (that are 3D maps of real numbers) along the direction u θ = (cos θ , 1, sin θ ) yielding, for each time step m and each direction u θ , 2D images that can be displayed. The projections are defined by:   It should be stressed that the cleaning algorithm is nonlinear, in the sense that it emphasizes the moving scatterers, whose signal is highest, which happen to be here the RBCs. Lower contributions to the signal such as vasomotion have not been detected. , and θ = +20 • (c). The display is made in arbitrary logarithmic scale: |S purple,θ | 2 is displayed in purple, |S green,θ | 2 in green. Bar is 100 µm.
To visualize the entire trajectories of the RBCs, that correspond to the inside of the per-fused blood vessels, we have calculated the time averaged intensity holograms |H purple | 2 and |H green | 2 . We got: We have displayed |S purple,θ | 2 and |S green,θ | 2 in a video file of 162 images (Visualization 8), where the angle varies from θ = −40 • to +40 • in 81 steps, and from θ = +40 • to −40 • in 81 steps. Figure 7 shows the images 1 (a), 41 (b) and 61 (c) of file Visualization 8, which correspond to θ = −40 • (a), θ = 0 • (b), and θ = +20 • (c). The 3D character of the blood vessels is clearly seen in the video file since the images of the blood vessels depends on the projection angle θ . To better visualize the 3D character of the motion of the RBCs, we have superimposed in a video file (Visualization 9) the averaged projections |S purple,θ | 2 and |S green,θ | 2 that show the perfused blood vessels, with the instantaneous projections |S purple,θ | 2 and |S green,θ | 2 that show the motion of individuals RBCs. In the video file Visualization 9, the time index varies from m = 1 to 122 in 122 steps, while the angle varies from θ = −30 • to +30 • in 61 steps, and from θ = +30 • to −30 • in 61 steps. The projections are displayed in arbitrary logarithmic scale: the averaged projections |S purple,θ | 2 and |S green,θ | 2 are displayed in blue and red, while the instantaneous projections |S purple,θ | 2 and |S green,θ | 2 are both displayed in green. With this choice of colors, the blood vessels are seen in purple, while the individual RBCs are white. Figure 9 shows the images 1 (a) , 32 (b) and 46 (c) of file Visualization 9, which correspond to time index m = 1 and projection angle θ = −30 • (a), m = 32 and θ = 0 • (b), m = 46 and θ = +15 • (c). The 3D character of the blood vessels is clearly seen in the video file. The 3D character of the RBCs motion is seen visible. As expected, the RBCs follow the blood vessels paths.

Conclusion
In this paper, we proposed an original imaging technique that combines a digital holographic microscopy setup with dual illumination of the sample, and calculations involving both stan-dard holographic reconstruction and 4D calculations by a cleaning algorithm. The setup uses an upright microscope that has been transformed into an heterodyne holographic device [24] working in transmission. The holograms H purple and H green that correspond to the scattering of each illumination beams were extracted from the camera recorded data I m . By making the hypothesis that the two holograms H purple and H green were generated by the same scatterers (i.e. the moving red blood cells), the location of the scatterers (x k , y k and z k ) and their complex scattering amplitudes S purple and S green were calculated by a cleaning algorithm. The amplitudes S purple and S green are 4D data (x, y, z and t) that were used to generate video files showing in 4D the RBCs (Visualization 7), the blood vessels (Visualization 8) and the RBCs superimposed with the blood vessels (Visualization 9).
In comparison with other techniques dedicated to monitor in vivo blood flow, including laser Doppler and laser speckle, the proposed approach has the major advantage, specifically thanks to the recording of holograms, to be able to measure 3D aspects of the RBC movements. The proposed method can be still improved by using a more efficient cleaning algorithm for the 3D reconstruction and a setup allowing a larger angle of separation of the two illuminations beams. These improvements should yield faster calculations and higher z resolution. Future work should allow to calculate 3D maps of velocity vectors by analyzing H purple and H green at successive instants of time t m and t m+1 .