I V ] 2 0 M ar 2 02 0 Spatio-temporal filtering in laser Doppler holography

Léo Puyo, 2, 3 Michel Paques, 3 and Michael Atlan Corresponding author: gl.puyo@gmail.com Centre Hospitalier National d’Ophtalmologie des Quinze-Vingts, INSERM-DHOS CIC 1423. 28 rue de Charenton, 75012 Paris France Institut de la Vision-Sorbonne Universités. 17 rue Moreau, 75012 Paris France Institut Langevin. Centre National de la Recherche Scientifique (CNRS). Paris Sciences & Lettres (PSL University). École Supérieure de Physique et de Chimie Industrielles (ESPCI Paris) 1 rue Jussieu. 75005 Paris France (Dated: March 24, 2020)


INTRODUCTION
Retinal motion caused by fixational eye movements and pulsatile swelling of the choroid prove challenging in ophthalmic imaging. Insufficient temporal resolution results in image blurring, or calls for complex registration schemes in the case of scanning methods. For interferometric imaging, eye motion additionally affect the backscattered light by the Doppler effect, which is especially problematic when the sought contrast is based on motion to measure blood flow or cellular activity. By now, hardware or software based methods have been developed to cope with the harmful effect of eye motion in many ophthalmic technologies. For example, in optical coherence tomography (OCT), it has been proposed to compensate eye motion using eye movements tracking with other modalities [1,2]. In OCT-angiography (OCT-A), the scanning pattern can be modified and motion artifacts can be detected and corrected for [3][4][5][6]. In Doppler-OCT, the bulk motion Doppler shift is also estimated and compensated for by various methods [7,8]. Axial motion can also be corrected for digitally in fullfield OCT by using a phase gradient autofocus algorithm [9]. Laser Doppler holography (LDH) is a digital holographic method where blood flow is measured from the interference of coherent light backscattered by the eye with a reference beam [10][11][12]. The coherent gain brought by the reference beam allows to use ultrahigh camera frame rates to perform full-field measurement of the local Doppler broadening. The power Doppler calculated pixel-wise as the integral of the high-pass filtered Doppler power spectrum density (DPSD) reveals blood flow thanks to the larger Doppler broadening of light scattered by red cells. By using a sliding short-time window, the variations over cardiac cycles of retinal blood flow measured in power Doppler units can be revealed with a few milliseconds of temporal resolution [13]. LDH is by nature very sensitive to eye motion as they affect all of the light backscattered by the eye which causes a strong contribution to the power Doppler. Global eye motion have for effect to degrade the blood flow contrast derived from the detection of light scattered by cells flowing in blood vessels. Intermittent motion also corrupt the measurement of blood flow variations. The pulsatile axial motion of the eye are especially challenging because they are partly synchronous to blood flow and occur in the direction maximizing the Doppler effect. In our previous work, we imaged blood flow in large vessels revealed by frequencies in the tens of kHz range, while removing eye motion occurring in the range of a few kHz [13]. By thresholding the power spectrum with a lower cutoff frequency of approximately 6-10 kHz, fast blood flow can be effectively isolated from eye motion. However when interested in slow blood flow or when in presence of strong eye motion, the Fourier space does not allow to separate these contributions. Yet, the detection of blood flow in smaller or pathological vessels is of crucial clinical interest. Interestingly, LDH shares the problem of slow flow detection with Doppler ultrasound imaging. The physical equations, biomarkers, and challenges to overcome are remarkably similar in optical and acoustic Doppler imaging, and one of these challenges is tissue motion. Ultrafast ultrasonic imaging, a technique based on the unfocused transmission of plane or diverging wave instead of the conventional focused emissions [14], has allowed the parallel measurement of blood flow, increasing the acquisition frame rate by more than an order of magnitude. This subsequently allowed to use spatio-temporal analysis of the Doppler fluctuations to filter tissue motion [15,16], as it has also been done in x-ray computed tomography (CT) [17], and magnetic resonance imaging (MRI) [18]. Eigen-based clutter filters make use of the differences in spatio-temporal characteristics of tissue and blood motion for a better discrimination. In Doppler ultrasound, it was found that a singular value decomposition (SVD) which is a eigendecomposition of a non-squared matrix, generates a basis where tissue and blood motion can be efficiently separated. A similar method has also been implemented for full-field interferometric measurement to remove motion artifacts [19]. We here report on the use of such type of spatio-temporal filter in the field of LDH.

METHOD
We use the LDH setup presented in [13]. Digital holograms resulting from the interference of light backscattered by the retina and an on-axis reference beam are recorded with a camera running at 60 to 75 kHz. Informed consent was obtained from the subjects, experimental procedures adhered to the tenets of the Declaration of Helsinki, and study authorization was sought from the appropriate local ethics review boards (clinical trial NCT04129021). The data processing applied to the reconstructed holograms consists of measuring the local optical temporal fluctuations of the holograms in order to reveal blood flow. In our previous work, this was done by calculating the temporal Fourier transform of the holograms over a short-time window of a few hundreds holograms. Then so-called power Doppler images were obtained by high-pass filtering of the DPSD. In order to implement a spatio-temporal filtering, an additional SVD filtering operation is here introduced prior to the calculation of the temporal Fourier transform. The subsequent data processing of the Doppler signal is left unchanged from our previous work [13,20,21]. For a given short-time window of n t consecutive holograms of spatial dimension n x × n y , the stack of holograms is reshaped into a 2D space-time matrix H of size (n x × n y , n t ), where each column of the matrix contains all the hologram pixels at a given time. A SVD allows to write H as a product of 3 matrices U , ∆, and V such that: Where t stands for the conjugate transpose. The process is represented in Fig. 1(a). ∆ is a rectangular and diagonal matrix of size (n x × n y , n t ) that contains n t singular values. Each singular value λ i is associated to both a spatial eigenvector U i (x, y) (i.e. a 2D image), and a temporal eigenvector V i (t) (i.e. a 1D pattern). The spatial and temporal eigenvectors are contained in the columns of the orthonormal matrices U and V of sizes (n x × n y , n x × n y ) and (n t , n t ), respectively. H can be written using these eigenvectors as: Where λ i are real-valued positive scalars ordered by decreasing value. Thanks to the SVD, H is written as the sum of n t independent terms, that are each determined by a specific spatial distribution of energy, and by a pattern of temporal variations, weighted by the corresponding singular value. The purpose of the SVD in this work is to constitute a basis where contributions modulated by the same pattern of temporal variations are regrouped in the same eigenvectors. It is assumed that retinal and corneal motion generate the same Doppler signal over the whole field of view, so their contributions is spatially very well correlated, and can thus be isolated in some singular values. That would lead to a better discrimination of the Doppler fluctuations originating of eye motion from those of blood flow than with the Fourier decomposition. The question is then how to identify the singular values associated to eye motion and blood flow. Because artifacts such as eye motion have a larger spatial coherence, it can be expected that their contribution lies in the singular values of highest energy. Conversely, blood flow induces Doppler shifts from the motion of individual red cells. This is a Doppler signal that comes from the uncorrelated motion of independent scatterers, which results in a low spatial coherence. Consequently, blood flow signal is distributed in more singular values of lower energy. Therefore, by truncating the first n c singular values from the matrix ∆ , it is possible to isolate blood flow from eye motion: Thus H f is the cube of holograms filtered from the singular contributions of highest energy. Afterwards, the beat frequency spectrum of the holograms filtered by SVD is analyzed pixel-wise by the conventional shorttime Fourier transform: (4) Where t win denotes the extent the short-time window, and f is the temporal frequency. Finally, the bandpass filtered DPSD is integrated over a given frequency range [f 1 , f 2 ] to obtain power Doppler images M 0 (x, y, t n ) that reveal blood flow: The process is then repeated for the next short-time window in order to obtain a time-resolved movie of blood flow fluctuations. In Fig. 1 is shown an example of SVD on LDH data. The LDH measurement was acquired at 75 kHz in n x n y = 512 × 512 pixels format, and n t = 1024 complex-valued holograms are used for the computation of the SVD. In Fig. 1(b) is plotted the energy of the n t ordered singular values in a logarithmic scale. The vertical red dashed line indicates the threshold n c = 100 that separates clutter from blood flow. In Fig. 1(c) is shown the power spectrum density of the ordered temporal eigenvectors weighted by singular values: each column of the image displays the magnitude of the DPSD. The horizontal dashed line indicates the filtering that would be obtained by Fourier thresholding of the DPSD, whereas the vertical dashed line shows the filtering that is obtained by truncating the singular value matrix. A few spatial eigenvectors are shown in Fig. 1(d): as postulated the first singular values originate from sources of clutter whereas the rest of the singular values reveal blood flow. It is interesting to observe that higher singular values tend to reveal faster blood flow: the spatial eigenvectors from 100 to 200 reveal smaller blood vessels than those from 200 to 500.
What is crucial to observe in Fig. 1(c) however, is that the SVD isolates in the first eigenvectors (cf. arrow) some contributions of very high frequencies that only reveal clutter, as it is visible from the corresponding spatial eigenvectors. These unwanted contributions are at very high frequency so they cannot be filtered by high-pass temporal filtering. This is why the SVD provides a better basis than the Fourier space to filter clutter from blood flow. Setting aside these singularities, the eigenvectors corresponding to singular values of high energy tend to have slower variations whereas singular values of lower energy are associated to faster variations.
The choice of the threshold to truncate the singular value matrix is delicate. In the field of Doppler ultrasound, some techniques have been proposed based on the characteristics of the spatial or temporal eigenvectors [22]. In this work we have used a fixed threshold based on the equivalent temporal frequency cutoff. For a window of n t holograms, whether with the SVD or fast Fourier transform (FFT), the cube of n t holograms is projected onto a set of n t vectors. Typically, for a measurement at f S = 60 kHz and a short time window of n t = 1024 holograms, a threshold equivalent to f 1 = 2 kHz is used: the SVD filtering is performed as described by setting n c to the integer nearest to n t .f 1 /f S . Then the power Doppler image is calculated with the f 1 lower cutoff frequency. This method amounts to rejecting the same number of elements with both decompositions. The choice of the 2 kHz threshold is the result of a compromise: a threshold too low comes with the risk of not rejecting all clutter, and a threshold too high prevents imaging slower blood flow.

Results in different situations
We first demonstrate the efficiency of the SVD filtering in Fig. 2 on LDH measurements strongly affected by motion artifacts when using the standard highpass Fourier filter. The power Doppler images calculated without and with the additional SVD filter are shown on the left and right, respectively. In Fig. 2(a), the 2-6 kHz power Doppler image obtained just with the Fourier highpass is largely dominated by strong retinal motion, which prevents blood flow from being revealed. However the SVD filtering applied to the same data allows to reject the retinal motion, so that blood flow can be revealed on the power Doppler image. In Fig. 2(b), the measurement was similarly affected by a strong corneal reflection (arrow) that could not be removed by Fourier filtering on the 2-6 kHz frequency range. Again, this spurious contribution can be efficiently removed with the SVD filtering. Finally, in the third example shown in Fig. 2(c) is demonstrated a case where the power Doppler image was strongly affected by camera jitter. Unlike previous examples, this time the parasitic signal affects the power Doppler image on the frequency range 6-37 kHz. This unwanted contribution takes the form of a diffraction pattern from the camera jittering pixel. Once again the SVD filtering allows to remove this artifactual contribution from the power Doppler image so that a clean image of blood flow can be produced. These examples demonstrate that a spatio-temporal filtering of the holograms is efficient to remove unwanted contributions characterized either by a strong spatial coherence in the case of eye motion, or by a specific temporal pattern in the case of camera jitter. Thanks to their higher spatial coherence, these clutter contributions are associated with singular values of higher energy, which are discarded by the SVD filter because it removes the first singular values. The SVD is able to provide a signal decomposition basis more adapted to the discrimination of clutter from blood flow than the Fourier basis. This way, the additional SVD filtering allows to access Doppler frequency ranges otherwise dominated by clutter. This digital filter is not only able to isolate the Doppler signal of blood flow from that of eye motion, but it is also able to compensate for flaws of the physical detection channel with a particular spatio-temporal pattern, such as the camera jitter.

SVD filtering of the eye motion
During the few seconds of a LDH measurement can occur fixational eye movements [23], and cardiac related axial retinal and corneal motion [24]. In Fig. 3, we show the effect of the SVD filtering in presence of eye motion. An acquisition performed at 66 kHz is considered, and the 2-33 kHz frequency band is used to reveal retinal blood flow. The short-time window is positioned at two moments corresponding to cases of mild and strong eye motion, where mild and strong motion are defined as generating Doppler signals below and above the 2 kHz threshold, respectively. All the plots displayed in the  Fig. 3(b).
We first investigate the effect of the SVD filtering on the DPSD. In Fig. 3(a) are represented the power spectra density of H and H f in presence of mild (left) and strong (right) eye motion. The vertical dashed line indicates the 2 kHz cutoff frequency used to calculate the power Doppler. In presence of mild motion, the DPSD with and without prior SVD filtering are essentially the same above the cutoff frequency: the SVD filtering only rejects signals of frequencies lower than the cutoff. As the power Doppler is calculated as the integral of the spectrum above the cutoff frequency, in presence of mild motion the power Doppler images with and without SVD filtering shown below in Fig. 3(b) are close to identical. In presence of strong eye motion however, the spectra now differ significantly above the cutoff frequency. Without SVD filtering, the DPSD shows a peak around 3 kHz corresponding to the eye movement. This contribution is integrated in the power Doppler, and the power Doppler image without SVD filtering is dominated by eye motion. For its part, the SVD filtered DPSD exhibits a dip around the same frequency where there should be the eye motion peak because the eye motion has been rejected. The corresponding power Doppler image is now able to The spectrograms of the holograms fluctuations with and without SVD filtering shown in Fig. 3(c) demonstrate how the SVD filter suppresses the power Doppler contribution due to eye motion occurring during cardiac cycles. The magnitude of the DPSDs for all short-time windows of the movie are represented in the column of the spectrogram in color tones. On the spectrogram on the left, showing the DPSD without prior SVD filtering, the eye motion occurring throughout cardiac cycles can be identified in red in the 0-6 frequency range. The frequency threshold required to filter eye motion is represented by the horizontal dashed line, at approximately 6 kHz: below is eye movement and above is fast blood flow. The spectrogram on the right shows the DPSD of the data filtered by SVD. It is possible to see that whenever they occur, global tissue motion are effectively rejected. This allows to lower the frequency threshold used for the calculation of the power Doppler down to 2 kHz (horizontal dashed line). It is also possible to observe that when there is no eye movement, the SVD only removes the signals below 2 kHz, because the SVD and the Fourier space provide very similar decompositions. This effect is also visible in Fig. 1(c) in the linear aspect of the Fourier decomposition of the ordered temporal eigenvectors.
Finally, the robustness with which the SVD filtering is able to reject eye motion as they occur throughout cardiac cycles is demonstrated in Fig. 3(d). Power Doppler variations are measured in the small artery, again with and without SVD filtering. These two plots show that without SVD filtering, the blood flow waveform cannot be revealed because there are too many eye movements. However when using the SVD filtering, the blood flow waveform can be revealed as eye motion contribution are mostly canceled. The two corresponding movies showing the variations of power Doppler without and with prior SVD filtering are shown on the left and right, respectively, in Visualization 1.
What can be concluded from Fig. 3 is that in presence of mild motion the SVD filtering does not modify the power Doppler image because the SVD and Fourier signal decomposition bases do not differ significantly. However in presence of strong motion, there is a contribution of high spatial coherence that leads to a difference between the SVD and Fourier basis. Then the SVD is able to reject those contributions that are of higher frequencies than the power Doppler cutoff frequency of the measurement band, and thus dramatically improves how power Doppler images reveal blood flow. The SVD filtering opens the possibility to exploit the low frequency ranges without being hindered by eye motion, i.e. LDH can measure lower flow velocities over cardiac cyces. This allows for example to combine low and high frequency ranges in composite images as detailed in [20] in blood flow movies, instead of doing it only on averaged images as before; this is shown in Visualization 2.

Retinal blood flow in pathology
Many retinal pathologies are associated with a decreased blood flow, which is challenging to image with LDH as it requires to work with lower frequencies. In Fig. 4, we show an example of blood flow imaging in a case of retinal vein occlusion. The superior hemivein is occluded, and preretinal hemorrhage has formed above it. As shown in Fig. 4(a), the 10-30 kHz power Doppler reveals all the large retinal vessels except the occluded vein because its blood flow is reduced. Without SVD filtering, the 1-6 kHz power Doppler image in Fig. 4(b) also fails to reveal the occluded vein because there are too many artifacts. However when applying the SVD filter in Fig. 4(c), the artifacts can be removed, and the 1-6 kHz power Doppler image is now able to reveal the flow in the occluded vein.
In Fig. 4(d) and Fig. 4(e) are displayed composite images obtained by stitching two LDH measurements. In Fig. 4(d), the 1-6 and 10-30 kHz frequency ranges are fused into a composite image that shows low flow in cyan and high flow in red. This image is very efficient to evidence the blood flow disruption in the superior hemivessels: healthy vessels would normally show, like the inferior hemivessels, a red lumen surrounded by an outer cyan sheath. However the occluded vein appears only in cyan, and the artery feeding the occluded vascular territory appears in red, with a weaker contrast. The systole and diastole timings of the high flow are fused in orange and blue in Fig. 4(e) and allow to easily identify arteries in orange and veins in blue. The image obtained in the same region by fluorescein angiography shown in Fig. 4(f) is not able to reveal the blood vessels through the preretinal hemorrhage.
This case of retinal vein occlusion shows that the SVD filter is necessary to detect an abnormal perfusion, because in a situation where blood flow is reduced the signal is revealed by the low frequencies subject to many parasitic Doppler contributions. It also demonstrates how LDH is relevant to study non-invasively retinal vascular diseases thanks to its ability to assess blood flow over an extended field of view. Furthermore, considering the specific dynamic fluctuations of blood flow that exist in cases of retinal vein occlusions [25], LDH could prove particularly useful for the study of this pathology.

Imaging the choroid
We previously showed that there are discrepancies of flow velocities between arteries and veins in the choroid which can be used to perform an arteriovenous differentiation [20]. This was previously demonstrated on a limited number of power Doppler images as it required manual filtering of images affected by motion artifacts. The SVD filtering is robust enough so that the low frequency range used to reveal the choroidal veins can now be accessed systematically, and manual filtering is no longer necessary. In the first row of Fig. 5, we show two LDH panoramas revealing the choroid obtained by stitching 5x5 individual images. Thanks to the SVD filter, the making of the 50 individual composite power Doppler images could be fully automated, and the low frequency range could be accessed in all measurements. These two LDH panorama reveal efficiently all the choroidal vessels, and the functional blood flow contrast offered by LDH allows to identify the vessels type. In the second and third rows of Fig. 5, we show the images obtained in the same eyes with optical coherence tomography (OCT) and indocyanin green angiography (ICG), which are the imaging modalities currently available in clinics to study the choroid en-face. Both OCT and ICG are not as efficient as LDH to reveal the choroidal arteries: the blood flow contrast brought by LDH allows to have a better understanding of the choroid anatomy. As shown earlier in this article, the SVD filtering is critical to access the low frequency range (here 2-6 kHz), which is the frequency range that reveals the choroidal veins and the smaller arteriolar vasculature. The LDH montages have been stitched from 5x5 images obtained by fusing the 2-6 and 10-30 kHz frequency ranges in cyan and red. The access to the low frequency range made routinely possible by SVD is critical to reveal choroidal veins.

CONCLUSION
The access to low frequencies in laser Doppler imaging is critical to improve the sensitivity to slower blood flow. Low frequency shifts are needed not only to reveal blood flow in smaller vessels, but also in large vessels where blood velocity is reduced because of vascular pathologies such as occlusions. It is also important to be able to use the low frequencies when imaging the choroid because veins are revealed by lower frequency shifts due to the lower velocities of flow in these vessels. Unfortunately, low Doppler frequency shifts from slow blood flow overlap in the Fourier space with the Doppler responses of unwanted contributions such as retinal and corneal motion. Adding the spatial dimensions into the analysis of the holograms fluctuations allows to take advantage of the greater spatio-temporal coherence of eye motion compared to the uncorrelated Doppler shifts of individual red cells. We have shown that rejecting the singularities of highest energy is a simple way to filter clutter from blood flow and greatly improve the quality of low frequency power Doppler images. Additionally, because the SVD generates a basis adapted to each short-time window, it is able to filter pulsatile or fixational eye movement during cardiac cycles. This dramatically improves the quality of the blood flow waveform measured at low frequency, and is very valuable to image blood flow in the eye of subjects with poor fixation ability. Finally, the SVD filter is also able to compensate for the camera jitter, thereby reducing the camera requirements. Overall, the SVD filtering allows to significantly lower the frequency threshold used to calculate the power Doppler and detect slower blood flow more robustly. This spatiotemporal analysis of the holograms fluctuations is effective because measurements are performed in full-field. Indeed it is necessary that all pixel receive the same signal so that it is correlated over the field of view and becomes a singularity in the SVD. Moreover, it is advantageous to dispose of a large number of spatial and temporal samples for the SVD to efficiently identify signals of high spatio-temporal coherence, thus the method seems particularly suitable for high-throughput imaging techniques. The spatio-temporal filtering demonstrated in this article confers LDH a decisive advantage compared to scanning laser Doppler techniques.
In conclusion, a spatio-temporal filtering can dramatically improve the LDH signal quality at low frequency in various situations. By doing so, it addresses the key issue of detecting slower blood flow, which is critical to reveal smaller vessels, vessels with a pathologically reduced blood flow, and choroidal veins. The ability to reveal blood flow in smaller retinal vessels will especially helpful for LDH to better reveal blood vessels in the optic nerve head, which could be of high clinical interest for the study of optic neuropathies. LDH appears as a promising method to study the ocular vasculature in healthy and diseased eyes.

Funding Information
This work was supported by the European Research Council (ERC Synergy HELMHOLTZ, grant agreement #610110). The Titan Xp graphics card used for this research was donated by the NVIDIA Corporation.