Photoplethysmographic analysis of retinal videodata based on the Fourier domain approach.

Assessment of retinal blood flow inside the optic nerve head (ONH) and the peripapillary area is an important task in retinal imaging. For this purpose, an experimental binocular video ophthalmoscope that acquires precisely synchronized video sequences of the optic nerve head and peripapillary area from both eyes has been previously developed. It enables to compare specific characteristics of both eyes and efficiently detect the eye asymmetry. In this paper, we describe a novel methodology for the analysis of acquired video data using a photoplethysmographic approach. We describe and calculate the pulsatile attenuation amplitude (PAA) spatial map, which quantifies the maximum relative change of blood volume during a cardiac cycle using a frequency domain approach. We also describe in detail the origin of PAA maps from the fundamental (the first) and the second harmonic component of the pulsatile signal, and we compare the results obtained by time-based and frequency-based approaches. In several cases, we show the advantages and possibilities of this device and the appropriate image analysis approach - fast measurement and comparison of blood flow characteristics of both eyes at a glance, the robustness of this approach, and the possibility of easy detection of asymmetry.


Introduction
The retina is the only human tissue allowing direct noninvasive imaging of the microvascular circulation and central nervous system. Therefore, besides the diagnosis of various retinal diseases, retinal imaging systems provide a unique possibility for the characterization of systemic and neural diseases. The importance of retinal imaging is emphasized in a recent vision paper [1], where oculomics is introduced as a new important direction of retinal imaging, which enables the identification of ocular biomarkers of systemic disease on the macroscopic, microscopic, and molecular level. Retinal tissue perfusion and the assessment of blood flow inside the optic nerve head (ONH) and the peripapillary area is currently a hot topic in retinal imaging. New functional imaging approaches -optical coherence tomographic angiography (OCT-A) or laser speckle flowgraphy (LSFG) are being developed to find clinically relevant parameters. These technologies have become important in many research and clinical setups for retinal and systemic cardiovascular disease assessment as they enable the measurement of functional biomarkers (e.g., papers related to OCT-A [2,3] and LSFG [4,5]). Moreover, phase sensitive OCT [6] offers the possibility of high temporal and spatial resolution. This, so called Doppler OCT, can provide a volumetric information about the blood flow together with the arterial and vein anatomy. Furthermore, recent development in this field has provided solutions for precise blood velocity quantification using multiple beams to solve the problem with unknown Doppler angle [7]. Another new promising technology based on Doppler holography has been introduced by Puyo et al. [8,9]. This technique provides high spatial and temporal resolution to measure and visualize retinal blood flow.
Recently, we have contributed to this research field with the development of an experimental binocular video ophthalmoscope [10] and we showed its potential in glaucoma research [11]. In this paper, we describe a novel method for the processing of acquired retinal video sequences to estimate the spatial distribution of cardiac cycle-induced light attenuation in retinal tissue using the principle of photoplethysmography [12]. The underlying principle is simple -as the light passes through the tissue, it is attenuated by blood. In the eye, the light is reflected back to the detector (CMOS camera of the instrument). The reflected light intensity depends on the blood volume in the retinal blood vessels and microcapillaries, which changes during the cardiac cycle. Thus the fundamental frequency of the light intensity changes is driven by the heart rate. This approach has been applied in a form of imaging photoplethysmography for the spatial distribution of the light attenuation analysis in various applications [13,14] including retinal imaging [15]. Here, we review the previously introduced [10] pulsatile attenuation amplitude (PAA), which quantifies the maximum change of light reflected from the retina during each cardiac cycle, which corresponds to light attenuation in retinal tissue. This quantity can be estimated for each cardiac cycle and averaged to obtain more reliable and robust results. Nevertheless, the previously described time-based approach needs a precise detection of particular cardiac cycles, which can be a source of errors due to various artefacts occurring during video acquisition. Thus, here we introduce a novel estimation of PAA based on Fourier transform and we discuss its properties. This new frequency-based approach has several advantages -simple implementation, straightforward relation to PAA, and the possibility to estimate various harmonic components.
The paper is organized as follows. The Method section summarizes the experimental binocular video ophthalmoscope, dataset and sequence preprocessing. The description of the basic concept of pulsatile attenuation amplitude and its estimation via Fourier domain is also provided here. Section 3 presents our main findings and Section 4 discusses the results with respect to other imaging modalities and our previous results. The paper is concluded in Section 5.

Methods
In this section, we briefly describe the main principle and properties of the binocular experimental video ophthalmoscope and describe the methodology for estimation of the blood volume changes based on the pulsating light intensity signal using the Fourier spectra (see Fig. 1). In brief, the time stack is preprocessed to remove retinal movements and consequently the frequency stack is computed. The detailed information of particular steps is provided below.
Recently, we have contributed to this research field with the development of an 46 experimental binocular video ophthalmoscope [10] and we showed its potential in glaucoma research [11]. In this paper, we describe a novel method for the processing of acquired retinal 48 video sequences to estimate the spatial distribution of cardiac cycle-induced light attenuation 49 in retinal tissue using the principle of photoplethysmography [12]. The underlying principle is 50 simple -as the light passes through the tissue, it is attenuated by blood. In the eye, the light is 51 reflected back to the detector (CMOS camera of the instrument). The reflected light intensity 52 depends on the blood volume in the retinal blood vessels and microcapillaries, which changes 53 during the cardiac cycle. Thus the fundamental frequency of the light intensity changes is driven 54 by the heart rate. This approach has been applied in a form of imaging photoplethysmography 55 for the spatial distribution of the light attenuation analysis in various applications [13,14] 56 including retinal imaging [15]. Here, we review the previously introduced [10] pulsatile 57 attenuation amplitude (PAA), which quantifies the maximum change of light reflected from the 58 retina during each cardiac cycle, which corresponds to light attenuation in retinal tissue. This 59 quantity can be estimated for each cardiac cycle and averaged to obtain more reliable and robust 60 results. Nevertheless, the previously described time-based approach needs a precise detection 61 of particular cardiac cycles, which can be a source of errors due to various artefacts occurring 62 during video acquisition. Thus, here we introduce a novel estimation of PAA based on Fourier 63 transform and we discuss its properties. This new frequency-based method approach has several 64 advantages -simple implementation, straightforward relation to PAA, and the possibility to estimate various harmonic components.

66
The paper is organized as follows. The Method section summarizes the experimental 67 binocular video ophthalmoscope, dataset and sequence preprocessing. The description of the 68 basic concept of pulsatile attenuation amplitude and its estimation via Fourier domain is also 69 provided here. Section 3 presents our main findings and Section 4 discusses the results with 70 respect to other imaging modalities and our previous results. The paper is concluded in Section 71 5.

74
In this section, we briefly describe the main principle and properties of the binocular 75 experimental video ophthalmoscope and describe the methodology for estimation of the blood 76 volume changes based on the pulsating light intensity signal using the Fourier spectra (see Fig.   77 1). In brief, the time stack is preprocessed to remove retinal movements and consequently the 78 frequency stack is computed. The detailed information of particular steps is provided below. Fig. 1. The principle of the presented approach shows the main steps of the data processing -81 the acquired time stack is processed frame-by-frame [16] to align the sequence; distorted frames Fig. 1. The principle of the presented approach shows the main steps of the data processing -the acquired time stack is processed frame-by-frame [16] to align the sequence; distorted frames are detected [17] (providing indexes of frame numbers); the temporal signals for each x,y position are Fourier transformed to the spectral domain (missing values due to distorted frames are temporally interpolated and spatial filter is applied on each frame before the Fourier transform).

Video ophthalmoscope
All video sequences were acquired using a recently developed binocular video ophthalmoscope, which was described in detail elsewhere [10,18]. In brief, it is based on the principle of a fundus camera. A standard ophthalmical lens (40D, Volk Optical Inc, USA) creates an image of the retina (i.e. an aerial image), which is projected by an objective lens on the sensor of a CMOS camera. The illuminating LED (577 nm, 15 nm FWHM) is placed in the image of the pupil plane. The light intensity at the cornea is set to 12 µW, resulting in a retinal illumination of 30 µW/cm 2 . In the aerial image plane, a small OLED display allows to present a fixation target to reduce eye movements during the acquisition. The area of interest of 1000 × 770 pixels is selected to match the field of view 20°× 15°, defined by the field aperture placed in the aerial image plane. The frame rate of the CMOS camera is set to 25 fps resulting in an exposure time of 39 ms. The instruments for the left and right eyes can be aligned independently. Hardware trigger of the two cameras allows exact synchronization. Video sequences of 10 s (250 frames) are acquired for each subject.
Although the instrument was built binocular to record both sequences simultaneously, here we analyze these sequences of both eyes independently. Simultaneous acquisition offers a possibility to analyze and evaluate also the timing between both eyes (see [10] for more details), which will be evaluated in more detail in the future. Thus, the DFT based method described here can also be applied using only a monocular video ophthalmoscope.

Experimental dataset
Data and images used in this study are part of the Erlangen Glaucoma Registry cohort (http:// www.clinicaltrials.gov, NCT00494923) founded in 1991. All participants signed an informed consent form approved by the ethical committee at the Friedrich-Alexander University of Erlangen-Nürnberg, in concordance with the Declaration of Helsinki. We included in this methodical examination 8 eyes of 4 subjects: 4 eyes without eye disease (H), one eye with ocular hypertension (OHT), and 3 eyes with perimetric glaucoma (G). The video sequences were acquired with the described binocular video ophthalmoscope with dilated pupils. The basic clinical parameters of the examined subjects/eyes are summarized in Table 1. Details of the different glaucoma groups are explained in detail in [11].

Preprocessing of video sequences
Eye movements during the acquisition cause distortion of particular frames. Therefore, several preprocessing steps must be performed before spectral analysis. The first step is the frame-toframe alignment. We employed our already published approach [16] for the frame translation and rotation correction. This method corrects large and small eye movements using a combination of phase correlation and tracking of blood vessel centerlines. The registration error was evaluated [13] on 23 sequences and it was 0.78 ± 0.67 pixels inside the optic disc and 1.39 ± 0.63 pixels outside the optic disc, respectively. Due to the fast eye movements, which may occur also during the frame exposure time, some frames might be blurred. Furthermore, eye blinking also causes significant distortion. Such frames cannot be restored and have to be eliminated from further analysis. The detection of these distorted frames is based on a machine learning approach, described in our previously published paper [17]. The aligned video sequences with information about the indexes of distorted frames in the sequence are further processed by interpolation of the pulsation signal on these indexes. For each pixel, the temporal profile of intensity is extracted (i.e., the pulsation signal is obtained) and the distorted values are interpolated by 1D cubic interpolation (see Fig. 3(C)).

Calculation of the pulsatile blood volume changes from video sequences
In the registered video sequence, each pixel on the camera chip represents a single detector, which measures the reflected light intensity from a specific site on the retina. Thus, the basic principle of photoplethysmography can be used to analyze the video sequences. The pixel-wise light transmission of the changing blood volume can be obtained as a ratio of the reflected intensity detected on the specific pixel I(n) and the maximum intensity of I(n) during the cardiac cycle (I max ), corresponding to the minimum blood volume in the examined tissue (Fig. 2). Therefore, the transmission T for each pixel as a function of frame number n (corresponding to the sampled time) can be expressed as: (1)

165
To compensate for these changes, we compute the relative intensity I'(n) as [10,11]: where I avg (n) represents the low frequency component (i.e., trend signal) estimated by a Neglecting the scattering, the attenuation (A) can be written as The maximum intensity change during one single cardiac cycle can be easily expressed using minimum intensity value I min . This quantity can be expressed in percent, which defines the pulsatile attenuation amplitude PAA: The PAA value represents the maximum change in light absorption during one single heartbeat, because only the pulsating part of the intensity signal is considered. PAA can be calculated for a single pixel but also for larger areas. As the unit of PAA we use "percentage attenuation" [%A] to distinguish it from error specifications.
The reflected pulsatile intensity signal I(n) can be distorted mainly by small eye and head movements and pupil dilation, resulting in illumination intensity changes. This effect can be observed in pulsatile signals as a low-frequency component influencing the overall analysis. To compensate for these changes, we compute the relative intensity I'(n) as [10,11]: where I avg (n) represents the low frequency component (i.e., trend signal) estimated by a moving average finite impulse response filter with the length of impulse response N=25 (i.e., 1 second). An example of the trend signal (I avg (n)) is shown in Fig. 3(C). More information and examples of the trend correction is presented in [11]. This ratio can be considered as a normalization which results in an aligned signal with average value equal to 1. Using this relative intensity I'(n), the PAA can be expressed as [10,11]: This value can be estimated for each heartbeat and then averaged over several heartbeats to obtain a more robust PAA estimation. However, in the so far used time-domain approach, this would need reliable detection of minima I min and maxima I max for each pulse, which can be challenging due to the low signal-to-noise ratio (particularly if we estimate PAA for each pixel) and in the presence of eye movements and eye blinking artefacts. Thus, our new approach for PAA estimation is based on the frequency domain as described in the next section.

Frequency analysis of retinal plethysmographic signal
The pulsation signal can be considered as a periodic signal driven by heart rate frequency, which defines the first (or fundamental) harmonic component. Therefore, transforming the pulsating signal of relative intensity into the frequency domain enables one to estimate its spectral position and amplitude. The second harmonic frequency can be consequently also estimated. Thus, the spectrum of the normalized pulsating signal is given as: where w(n) represents the Hamming window function used for spectral leakage elimination and I'(n) represents the pulsation signal extracted from a specific region and from the entire sequence duration. Consequently, the heart rate frequency (f 1 ) is determined as the position of the limited our analysis only to the 1 st and 2 nd harmonics. Fig. 3 represents an averaged magnitude 203 spectrum of the 10 photoplethysmographic signals for ten neighbouring pixels. The first and 214 (see Eq. (4)). Here, the light intensity is shown. For the blood volume, the signal has to be This approach for PAA k estimation can be applied for a defined region of interest or 218 generally for each pixel. To decrease the noise and therefore the robustness and quality of the 219 estimation, we applied a spatial averaging filter to each frame (3x3 kernel) before taking the Fourier transformation. Finally, this pixel-by-pixel computation provides two images The strong low frequency components below 20 bpm represent the fluctuation of the mean value. C) An example of a pulsation signal extracted from a small region inside ONH is shown as yellow color (with interpolated missing values, before trend correction) and the original signal with artifacts is shown as a black curve. The red curve shows the low frequency trend signal I avg (n) (see Eq. (4)). Here, the light intensity is shown. For the blood volume, the signal has to be inverted (see Fig. 2). maximum amplitude in the spectral range, corresponding to the physiological values in the resting conditions: 50-100 beats per minute (bpm), i.e. approximately 0.83-1.67 Hz. The amplitude of this 1 st harmonic frequency (Amp 1 ) is used for spectral-based estimation of the PAA parameter (see Eq. (7)). Furthermore, the position of the 2 nd harmonic frequency (f 2 ) is also estimated similarly, as a maximum in a small range defined by 3 spectral lines around 2 × f 1 , which corresponds to ±18bpm. This frequency analysis results in the following quantities: amplitudes and positions of the 1 st and 2 nd harmonic, denoted as Amp 1 , Amp 2 and f 1 , f 2 , respectively. These amplitudes can be used for alternative PAA estimation, specifically for these two frequencies.
The equation for PAA can be generally written in the following form: where Amp k represents the amplitude of the k th harmonic. If we consider that the mean value of the relative intensity signal I'(n) is equal to 1 [Eq. (5)] then the term 1-Amp k represents the I' min and 1+Amp k represents the I' max in equation [Eq.
(3)], respectively. Generally, this equation can be used for higher harmonic frequencies. However, due to the limitation given by the noise level in our image data, which influences the level of higher harmonics, we have limited our analysis only to the 1 st and 2 nd harmonics. Figure 3 represents an averaged magnitude spectrum of the 10 photoplethysmographic signals for ten neighbouring pixels. The first and second harmonic components are easily recognizable while the third component cannot be clearly localized. This approach for PAA k estimation can be applied for a defined region of interest or generally for each pixel. To decrease the noise and therefore the robustness and quality of the estimation, we applied a spatial averaging filter to each frame (3 × 3 kernel) before taking the Fourier transformation. Finally, this pixel-by-pixel computation provides two images PAA 1 (x,y) and PAA 2 (x,y) of PAA spatial distribution for the first and the second harmonic frequency, respectively. In the following, these spatial distributions are called PAA 1 or PAA 2 maps, and the summation of these maps is called PAA 1+2 map, i.e., PAA 1+2 (x,y) = PAA 1 (x,y) + PAA 2 (x,y).
The images included in this paper are presented without any additional clipping or resizing. Black pixels at the borders are due to image registration (translation and rotation) and consequently

279
The strong pulsation of Subject 1 (Fig. 4) of retinal blood vessels is presented in 280 Visualization 2 (Subject 1, only left eye, 6 heart beats), which shows the entire field of view 281 together with two magnified areas (corresponding to the grey boxes in full field-of-view video).

282
The pulsating vessel diameter and the moving of the vessels can be clearly seen. A precise  reduced area for calculations. The arcuate structure on the border of some images (e.g., Fig. 4(B) and 5(B)) is due to a not properly aligned aperture in the instrument that has no influence on the results in areas outside this structure. The black boxes in the upper left corner of some images cover patient data.

Parallel measurement
The results of a parallel measurement of two subjects are shown in Fig. 4 and Fig. 5; the subject with no eye disease is shown in Fig. 4, while the glaucoma subject with asymmetric glaucoma diagnosis on the left and right eye is shown in Fig. 5. Color fundus images taken with a commercial fundus camera are included for comparison (A). The color fundus images are registered to the corresponding average images of the video sequence (B), so that a pixel-perfect comparison of all images of one eye is possible. Subfigure C shows the calculated PAA 1 map and D the same image color coded. The acquired video sequence for Subject 1 in Fig. 4 is shown in Visualization 1.
With the help of the PAA map, four different regions with higher PAA values can be identified: 1. Regions inside and near ONH without visible blood vessels -this pulsation is caused by blood volume changing inside the microcapillary network.
2. Regions inside the blood vessels -this is caused by changing blood volume inside particular blood vessels (arteries or veins).
3. Regions on both sides of the blood vessels -this is connected with diameter changes of particular blood vessels (arteries or veins). During the diameter expansion, the centerline of the vessel is fixed, only the diameter changes.
4. Regions on one side of the blood vessels -this is because of the bending of the blood vessels and it typically occurs in the arteries. This causes a small displacement of the blood vessels centerline, which we call bending. Additionally, the diameter can change.
Bending of a vessel together with diameter expansion could result in high PAA values only on one side of the blood vessel.
Examples of these different regions are marked in Fig. 4(C) left. Corresponding areas can be typically found in all PAA maps.
The PAA 1 maps of the subject with asymmetric glaucoma (Fig. 5) show clearly visible asymmetry between the left and right PAA 1 maps, respectively. The intensity is considerably reduced in the left eye (i.e., eye with perimetric glaucoma diagnosis) compared to the right eye, while the imaging modalities (fundus image and video sequence average image) look relatively similar. The lower ONH blood flow reduction was confirmed by several papers using different OCT-A and LSFG modality, e.g., [5,19,20] and we also showed lower ONH PAA values for glaucoma subjects in our previous study [11].
The pulsating areas and the differences between eye sides can also clearly be seen in the color coded images (Figs. 4(D) and 5(D)). The colormap has been designed such that it emphasizes differences in the range from 0%A to 8%A, which are the typical values of microcapillary pulsation. Strong pulsations connected with the main blood vessels are typically higher and reach up to 15%A.
The strong pulsation of Subject 1 (Fig. 4) of retinal blood vessels is presented in Visualization 2 (Subject 1, only left eye, 6 heartbeats), which shows the entire field of view together with two magnified areas (corresponding to the grey boxes in full field-of-view video). The pulsating vessel diameter and the moving of the vessels can be clearly seen. A precise observation even shows a small time shift between moving arteries and pulsating veins. For better observation, the Visualization 3 shows a zoomed video sequences of these two areas close together.

First and second harmonic
The Fourier analysis enables to visualize not only the spatial distribution of the first harmonic but also of the second harmonic amplitude, i.e., PAA 2 map. The results of the same subjects as in Figs. 4 and 5 are shown in Fig. 6 together with the fundus image, PAA 1 map, PAA 2 map, and PAA 1+2 map. The figures also contain black and white arrows, which depict some of the main arteries and veins, respectively. It can be observed that the pulsation phenomenon is also included in the 2 nd harmonic component (Fig. 6, third row). However, the magnitude of the second harmonic component is of course much lower. Therefore, to make it visible, the colormap is preserved, but the magnitude of the corresponding color is halved with respect to the 1 st harmonic component.  Other examples are shown in Fig. 7. Parallel measurement of the subject with perimetric glaucoma on both eyes is shown on the left side. The peripapillary atrophy (zone beta) is visible in fundus images on the temporal side of both eyes (marked as black arrows). The clear morphological correspondence of this pathology with PAA maps suggests a strong correlation between PAA and fundus images. Zone beta is typically adjacent to ONH and contains an atrophied retinal pigment epithelium, which changes the blood volume in this area. Another subject with no eye disease is shown on the right side (Fig. 7). A strong PAA magnitude inside ONH and a ring-shaped of increased PAA corresponding to the optic rim and close surroundings are clearly visible.

Comparison of the time-and frequency-based method
We compared the time-based and frequency-based PAA estimation methods using data from all eight eyes (Table 1). Figure 8 shows a few examples of PAA maps estimated by these approaches. The time-based PAA maps are shown in the second row, whereas frequency-based maps are in the third row. The time-based results have been obtained by careful manual selection of undistorted pulses for each recording to determine I' max and I' min for each of these pulses. Typically, from three to five pulses were selected from each sequence for this analysis. Finally, Eq. (5) was

347
Altman plot (Fig. 9). Because the PAA differences were not normally distributed (according to applied to estimate PAA in each pixel. Therefore, the maps generated by this time-based approach can be considered as a reference method, since we have avoided all distorted pulses. It can be seen that the spatial distribution of the PAA values for both approaches are very similar. However, the image contrast is higher for frequency-based approach (see especially the blood vessels in the first and the second column).
To compare the PAA values of these two approaches, we have automatically determined squared regions (15 × 15 pixels) on the approximately 30 × 30 grid in PAA maps for each eye (ignoring the border pixels), which provides almost 900 patches for each eye. In such a way, these regions uniformly covered ONH and peripapillary area, including areas with the blood vessels. Next, we determined an average value within these regions for the time-based and frequency-based PAA maps. The comparison of these averages was done with the Bland-Altman plot (Fig. 9). Because the PAA differences were not normally distributed (according to Shapiro-Wilk test), we estimated the limits of agreement (dotted lines in Fig. 9) with the help of regression as described in [21]. Finally, the Bland-Altman plot shows good agreement between these two approaches for the small PAA values up to approximately 8%A. For the higher values there are higher differences. However, these values typically occur for strong arterial and vein pulsation. In these regions quantitative analysis is not needed, because we are focusing on microcapillary perfusion, which has typically values below 8%A. The difference of PAA values shows only small

374
The described method allows to evaluate cardiac cycle-induced light intensity changes in the 375 ONH and peripapillary region. The PAA 1 maps visualize the maximum changes of reflected 376 light due to the heart rate frequency. It is obvious that blood volume changes occur at this 377 frequency and cause light absorption changes, which are consequently captured by the video 378 ophthalmoscope and appropriate image processing.

379
The second (and higher) harmonic components have shown to be promising indicators in 380 different clinical and experimental setups, including indicators of a vascular aging marker [22] 381 or relation to blood pressure [14]. Here, from the higher harmonic components, we focused on 382 the second harmonic component only and its spatial distribution in ONH and peripapillary area.

383
It is shown that the PAA 2 map correlates with the ONH and blood vessel morphology, which 384 shows its relevance and ability to determine higher harmonic pulsation.

385
The PAA 2 maps highly correlate with PAA 1 maps, i.e., the regions with increased 386 (decreased) PAA 1 values correspond to regions with increased (decreased) PAA 2 values. dependence on the average PAA as presented by the blue solid line in Fig. 9. The total average PAA difference between the frequency-based and time-based approach is equal to 0.61%A, with the frequency-based method providing higher PAA values.

Discussion
The described method allows to evaluate cardiac cycle-induced light intensity changes in the ONH and peripapillary region. The PAA 1 maps visualize the maximum changes of reflected light due to the heart rate frequency. It is obvious that blood volume changes occur at this frequency and cause light absorption changes, which are consequently captured by the video ophthalmoscope and appropriate image processing.
The second (and higher) harmonic components have shown to be promising indicators in different clinical and experimental setups, including indicators of a vascular aging marker [22] or relation to blood pressure [23]. Here, from the higher harmonic components, we focused on the second harmonic component only and its spatial distribution in ONH and peripapillary area. It is shown that the PAA 2 map correlates with the ONH and blood vessel morphology, which shows its relevance and ability to determine higher harmonic pulsation.
The PAA 2 maps highly correlate with PAA 1 maps, i.e., the regions with increased (decreased) PAA 1 values correspond to regions with increased (decreased) PAA 2 values. However, the PAA 2 values are noisier in the peripapillary, which is mainly due to the low signal-to-noise ratio at the second harmonic frequency. The summation of the 1 st and 2 nd components then represents the maximum changing blood absorption during the cardiac cycle due to these two most contributing frequency components. We observed from our data that the amplitude of these two harmonic components covers approximately 12.6 ± 2.1% of the magnitude spectrum in the range between 0.5 Hz and 12.5 Hz (half of the sampling rate, i.e. the frame rate). The low frequency components below 0.5 Hz were ignored due to the nonzero mean value. This value suggests a relatively high contribution of these two components.
As mentioned, the PAA 1 and PAA 2 maps show similar distributions of PAA values with noticeable noisier appearance of the PAA 2 ( Fig. 6 and Fig. 7). However, as the spectrum is influenced by the shape of the pulsation it also influences the magnitude of the harmonics. For example, it has been shown that the second harmonic of the plethysmographic signal corresponds to a dicrotic notch in the descending part of the particular pulses [24].
In our recent study [11], we have already shown that selected parameters derived from the plethysmographic signal extracted from ONH correlate with RNFL thickness and glaucoma progression. Here we observed that the highest PAA values are typically seen inside the ONH area and in the peripapillary area as a result of microcapillary perfusion. This is in line with some recent studies, which uses different modalities to assess the ONH perfusion. For example, it has been shown [19] that ONH perfusion assessed by OCT-A showed large flow index and vessel density reduction for three preperimetric glaucoma subjects. Another study based on OCT-A by Chen et al. [20] described significant ONH perfusion reduction in glaucomatous eyes of several perfusion characteristics (blood flux, vessel density area, normalized flux) from prelaminar tissues. The results obtained by LSFG show similar results for glaucoma subjects. Mursch-Edlmayr et al. [5] found a reduced blood flow measured by LSFG in the vessels supplying the retina assessed by mean blur rate (MBR) from ONH in NTG subjects. Recently, Kuroda et al. [25] showed that the MBR extracted from ONH tissue is a promising factor for assessing the blood flow perfusion in glaucomatous eye and its correlation with field loss. Another recently introduced technology in ophthalmology research is Doppler holography [26], which showed its potential to characterize ONH blood flow at microcapillary level using spatial coefficient of variation. For our cases, we showed that we can recognize PAA reduction for glaucomatous eyes using this novel spectral-based technique. Furthermore, a simple comparison of both eyes at a glance enables to effectively assess the asymmetry of PAA maps. As asymmetry is an important factor of glaucoma progression [27][28][29], the PAA-based asymmetry assessment will be the next step of our research.
It should be noted that the PAA values are influenced by the noise (due to low light level and high camera gain) and eye-movement artifacts. The noise will be more apparent in the regions with the lower PAA values and of course in PAA 2 maps. Here, the presented PAA maps are without any additional filtering and processing, to show the capability and real outputs of this approach. Additional processing, including noise level estimation, is left for further investigation. The only processing to decrease the influence of the noise and frequency sampling is the spectral averaging of three spectral lines as described in Section 2.5. Furthermore, the precision of the frame-to-frame registration can also influence the results. Although the registration precision is relatively high [16], small errors in the range of 1 to 2 pixels can create artificially high PAA values in the close vicinity of arteries and veins.
The PAA maps can be efficiently used for visualization of blood vessel pulsations. The presence of well-known, but not fully described, spontaneous vein pulsation (SVP) phenomena can be easily localized, particularly with the help of color fundus images, where the veins can be identified relatively easily. The SVPs have been visualized and analyzed by Moret et al. [30,31] using principal component analysis, which is also applicable for our modality. However, in Moret's papers the frame rate is only 9 frames/second, which provides only a few samples per cardiac cycle (depending on the heart rate) and might be insufficient for reliable analysis. Another recent paper [32] describes the simultaneous acquisition and analysis of blood flow and blood vessel diameter changes using Doppler OCT technique. This promising technique provides an advanced imaging tool with high temporal and spatial resolution to study retinal blood circulation. However, only twelve healthy subjects were included in the study and data from only two subjects were presented in the paper. Another phenomenon, which can be observed in PAA maps, is the bending of the pressured arteries. This is visible during each cardiac cycle, typically on one side of the artery. The bending changes can be used for the identification of the delay between pulsation of the retinal arterioles and venules, i.e., pulse wave analysis [33,34].
The  The red color inside and around the darker blood vessels represents the vein pulsation and the red color near the lighter blood vessels depicts the artery bending. The CLAHE method was applied on the green and blue components of the fundus image. The black arrows point to bending blood vessels (arteries in these examples) and the white arrows point to pulsating blood vessels (veins in these examples). Visualization 2 also shows bending and pulsating vessels with the same color-coded arrows.
PAA maps as a tool for morphological and functional image data merging. Two results of this fusion are shown in Fig. 10. The red channel of the registered RGB fundus image was replaced by the PAA 1+2 map. The values in this map were normalized into the range [0,255]. Finally, the green and blue channels of the original fundus image were equalized (using contrast-limited adaptive histogram equalization (CLAHE) [35] to emphasize the difference between arteries and veins. The strong pulsation can be identified inside and on both sides of the veins in and near the ONH. The bending of the arteries are slightly visible for Subject 2 near the upper artery.

Conclusion
In this methodology paper, we provided a description of an approach for local light attenuation assessment of retinal pulsatile changes via Fourier spectrum. This approach is an alternative to the time-based approach, because using only the 1 st and 2 nd harmonic components is enough to compute the PAA spatial maps (given as PAA 1 + PAA 2 ), which is comparable to the time-based approach. The introduced frequency-based approach has some advantages over the time-based approach. Firstly, this approach does not need any detection of particular pulses, because it works with the whole measured sequence. This makes it more robust, because the reliable pulse onset and end detection for the time analysis is a source of possible errors and depends on the evaluator. Therefore, the application in a clinical setting or population studies is easier. Furthermore, this approach enables us to estimate different harmonic frequencies.
Here we focus only on the fundamental and the second harmonic, but the analysis of higher harmonics and their combination is a promising direction to go. Another advantage is the straightforward relation of harmonic amplitudes to the PAA, which simplifies the interpretation with respect to plethysmographic principles. Finally, the implementation can be efficient with fast Fourier transform and thus the final PAA maps can be computed fast. Nevertheless, this Fourier-based approach is one possible way to evaluate the pulsatile signal. The shape of the waveform analyzed in the time domain is also a source of possible characteristics as we showed in [11], where the steepness of the ascending part of the pulsation waveform is correlated with the severity of glaucoma.
This method can be applied to assess the parameters of blood supply of capillaries of the ONH but also of the parameters of vessel movement and diameter change. An important advantage for clinical application is the very short measurement time (10 s) for both eyes together and the possibility to compare both sides.