Open Access
Add to BibSonomy
Abstract
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.
© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. 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.
2. 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.
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).
2.1 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/cm2. 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.
2.2 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].
Table 1. Main information about the presented subjects. Column Eye: R/L for right and left eye, respectively. Column Diagnosis: H - no eye disease, OHT - ocular hypertension, G - perimetric glaucoma. Column RNFL: the mean retinal nerve fiber layer thickness value along the circular scan. Column IOP: intraocular pressure value. Column MD: mean deviation value from perimetry.
2.3 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-to-frame 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)).
2.4 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 (Imax), 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:
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 Imin. 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 Iavg(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 (Iavg(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]:Fig. 3. A) An example of 10 plethysmographic signals of 10 neighbouring pixels directly extracted from a recorded video sequence from a specific location on the retina. B) corresponding averaged magnitude spectrum with visible 1st and 2nd harmonic components. 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 Iavg(n) (see Eq. (4)). Here, the light intensity is shown. For the blood volume, the signal has to be inverted (see Fig. 2).
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 Imin and maxima Imax 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.
2.5 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 (f1) is determined as the position of the 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 1st harmonic frequency (Amp1) is used for spectral-based estimation of the PAA parameter (see Eq. (7)). Furthermore, the position of the 2nd harmonic frequency (f2) is also estimated similarly, as a maximum in a small range defined by 3 spectral lines around 2 × f1, which corresponds to ±18bpm. This frequency analysis results in the following quantities: amplitudes and positions of the 1st and 2nd harmonic, denoted as Amp1, Amp2 and f1, f2, 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:This approach for PAAk 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 PAA1(x,y) and PAA2(x,y) of PAA spatial distribution for the first and the second harmonic frequency, respectively. In the following, these spatial distributions are called PAA1 or PAA2 maps, and the summation of these maps is called PAA1+2 map, i.e., PAA1+2(x,y) = PAA1(x,y) + PAA2(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 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.
Fig. 4. Results of a parallel measurement (right and left eye) for Subject 1 (without eye disease). A: color fundus image, B: average images of registered video sequences, C: calculated PAA1 map computed according to [Eq. (7)] for each pixel, D: calculated PAA1 map color coded. The yellow circle shows the position of the ONH for better orientation and comparison between different imaging modalities. The numbers in C correspond to the type of different regions, which can be typically found in all PAA maps (see text for more details).
Fig. 5. Results of a parallel measurement (right and left eye) for Subject 2 with OHT on the right eye and perimetric glaucoma on the left eye. A: color fundus image, B: average images of registered video sequences, C: calculated PAA1 map computed according to [Eq. (7)] for each pixel, D: calculated PAA1 map color coded. The yellow circle shows the position of the ONH for better orientation and comparison between different imaging modalities.
3. Results
3.1 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 PAA1 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 PAA1 maps of the subject with asymmetric glaucoma (Fig. 5) show clearly visible asymmetry between the left and right PAA1 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.
3.2 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., PAA2 map. The results of the same subjects as in Figs. 4 and 5 are shown in Fig. 6 together with the fundus image, PAA1 map, PAA2 map, and PAA1+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 2nd 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 1st harmonic component.
Fig. 6. Examples of parallel measurements for Subjects 1 and 2 (from top to bottom): fundus images; the 1st harmonic component PAA1(x,y); the 2nd harmonic component PAA2(x,y) (please note the halved magnitude of the colormap); sum of the 1st and 2nd with white inpainted blood vessels and ONH borders. The white and black arrows show the position of some main veins and arteries, respectively.
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.
Fig. 7. Examples of parallel measurements for two subjects 3 and 4 (from top to bottom): fundus images; the 1st harmonic component PAA1(x,y); the 2nd harmonic component PAA2(x,y) (please note the halved magnitude of the colormap); sum of the 1st and 2nd with white inpainted blood vessels and ONH borders. The white and black arrows show the position of some main veins and arteries, respectively. The blue arrows for the subject on the left indicate a moon-shaped area of peripapillary atrophy (zone beta) around ONH.
3.3 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 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).
Fig. 8. Examples of time-based and frequency-based methods. Five PAA spatial maps were computed by time-based method (second row) and frequency-based method (third row); the summation PAA1+PAA2 is presented. The first row shows the averaged image from the acquired sequence.
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 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.
Fig. 9. Bland-Altman plot for frequency-based and time-based methods for PAA estimation. A dependence of PAA differences on PAA average values is noticeable. The higher PAA difference occurs for higher PAA values. The total average PAA difference between the frequency-based and time-based approach is equal to 0.61%A. The units %A denotes percent attenuation.
4. Discussion
The described method allows to evaluate cardiac cycle-induced light intensity changes in the ONH and peripapillary region. The PAA1 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 PAA2 map correlates with the ONH and blood vessel morphology, which shows its relevance and ability to determine higher harmonic pulsation.
The PAA2 maps highly correlate with PAA1 maps, i.e., the regions with increased (decreased) PAA1 values correspond to regions with increased (decreased) PAA2 values. However, the PAA2 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 1st and 2nd 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 PAA1 and PAA2 maps show similar distributions of PAA values with noticeable noisier appearance of the PAA2 (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–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 PAA2 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 visualization of these two phenomena can be performed via fusion of fundus image and 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 PAA1+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.
Fig. 10. An example of efficient fusion of RGB fundus image with PAA map for Subject 1 on the left side and Subject 2 on the right. The red channel was replaced by the PAA1+2 map. 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.
5. 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 1st and 2nd harmonic components is enough to compute the PAA spatial maps (given as PAA1 + PAA2), 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.
Funding
Deutsche Forschungsgemeinschaft (TO 115/3-1); Grantová Agentura České Republiky (21-18578S).
Acknowledgments
We thank Robert Lämmer, Bettina Hohberger, Christian Mardin, and Folkert Horn for clinical support.
Disclosures
The authors declare no conflicts of interest.
Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.
References
1. S. K. Wagner, D. J. Fu, L. Faes, X. Liu, J. Huemer, H. Khalid, D. Ferraz, E. Korot, C. Kelly, K. Balaskas, A. K. Denniston, and P. A. Keane, “Insights into systemic disease through retinal imaging-based oculomics,” Transl. Vis. Sci. Technol.9(2), 6 (2020). [CrossRef]
2. A. Gruber, R. K. Wang, S. Hurst, S. R. Hanson, S. L. Jacques, and Z. Ma, “Three dimensional optical angiography,” Opt. Express 15(7), 4083–4097 (2007). [CrossRef]
3. A. M. Hagag, S. S. Gao, Y. Jia, and D. Huang, “Optical coherence tomography angiography: technical principles and clinical applications in ophthalmology,” Taiwan J. Ophthalmol. 7(3), 115–129 (2017). [CrossRef]
4. N. Konishi, Y. Tokimoto, K. Kohra, and H. Fujii, “New laser speckle flowgraphy system using CCD camera,” Opt. Rev. 9(4), 163–169 (2002). [CrossRef]
5. A. S. Mursch-Edlmayr, N. Luft, D. Podkowinski, M. Ring, L. Schmetterer, and M. Bolz, “Laser speckle flowgraphy derived characteristics of optic nerve head perfusion in normal tension glaucoma and healthy individuals: a Pilot study,” Sci. Rep. 8(1), 5343 (2018). [CrossRef]
6. R. A. Leitgeb, R. M. Werkmeister, C. Blatter, and L. Schmetterer, “Doppler optical coherence tomography,” Prog. Retin. Eye Res. 41, 26–43 (2014). [CrossRef]
7. A. Wartak, B. Baumann, C. K. Hitzenberger, M. Pircher, R. Haindl, and W. Trasischker, “Total retinal blood flow measurement by three beam Doppler optical coherence tomography,” Biomed. Opt. Express 7(12), 5233 (2016). [CrossRef]
8. J.-A. Sahel, L. Puyo, M. Atlan, M. Fink, and M. Paques, “In vivo laser Doppler holography of the human retina,” Biomed. Opt. Express 9(9), 4113–4129 (2018). [CrossRef]
9. L. Puyo, M. Atlan, M. Paques, M. Paques, and M. Atlan, “Spatio-temporal filtering in laser Doppler holography for retinal blood flow imaging,” Biomed. Opt. Express 11(6), 3274–3287 (2020). [CrossRef]
10. R.-P. Tornow, J. Odstrcilik, and R. Kolar, “Time-resolved quantitative inter-eye comparison of cardiac cycle-induced blood volume changes in the human retina,” Biomed. Opt. Express 9(12), 6237 (2018). [CrossRef]
11. R.-P. Tornow, R. Kolar, J. Odstrcilik, I. Labounkova, and F. Horn, “Imaging video plethysmography shows reduced signal amplitude in glaucoma patients in the area of the microvascular tissue of the optic nerve head,” Graefe’s Arch. Clin. Exp. Ophthalmol. 259(2), 483–494 (2021). [CrossRef]
12. J. Allen, “Photoplethysmography and its application in clinical physiological measurement,” Physiol. Meas. 28(3), R1–R39 (2007). [CrossRef]
13. D. Djeldjli, F. Bousefsaf, C. Maaoui, and F. Bereksi-Reguig, “Imaging photoplethysmography: signal waveform analysis,” in Proceedings of 10th IEEE International Conference on Intelligent Data Acquisition and Advanced Computing Systems: Technology and Applications (IEEE, 2019), pp. 830–834.
14. O. A. Lyubashina, O. V. Mamontov, M. A. Volynsky, V. V. Zaytsev, and A. A. Kamshilin, “Contactless assessment of cerebral autoregulation by photoplethysmographic imaging at green illumination,” Front. Neurosci. 13, 1235 (2019). [CrossRef]
15. H. Hassan, S. Jaidka, V. M. Dwyer, and S. Hu, “Assessing blood vessel perfusion and vital signs through retinal imaging photoplethysmography,” Biomed. Opt. Express 9(5), 2351 (2018). [CrossRef]
16. R. Kolar, R. P. Tornow, J. Odstrcilik, and I. Liberdova, “Registration of retinal sequences from new video-ophthalmoscopic camera,” Biomed. Eng. Online 15(1), 57 (2016). [CrossRef]
17. R. Kolar, I. Liberdova, J. Odstrcilik, M. Hracho, and R. P. Tornow, “Detection of distorted frames in retinal video-sequences via machine learning,” Proc. SPIE 10413, 104130A (2017). [CrossRef]
18. R. P. Tornow, A. Milczarek, J. Odstrcilik, and R. Kolar, “Binocular video ophthalmoscope for simultaneous recording of sequences of the human retina to compare dynamic parameters,” Proc. SPIE 10413, 1041309 (2017). [CrossRef]
19. B. Baumann, C. D. Lu, D. Huang, J. G. Fujimoto, J. Tokayer, J. C. Morrison, L. Lombardi, O. Tan, W. Choi, and Y. Jia, “Quantitative OCT angiography of optic nerve head blood flow,” Biomed. Opt. Express 3(12), 3127–3137 (2012). [CrossRef]
20. C.-L. L. Chen, K. D. Bojikian, D. Gupta, J. C. Wen, Q. Zhang, C. Xin, R. Kono, R. C. Mudumbai, M. A. Johnstone, P. P. Chen, and R. K. Wang, “Optic nerve head perfusion in normal eyes and eyes with glaucoma using optical coherence tomography-based microangiography,” Quant. Imaging Med. Surg. 6(2), 125–133 (2016). [CrossRef]
21. J. M. Bland and D. G. Altman, “Measuring agreement in method comparison studies:,” Stat. Methods Med. Res. 8(2), 135–160 (1999). [CrossRef]
22. T.-H. Nam, Y.-B. Park, Y.-J. Park, and S.-H. Shin, “Age-related changes of the finger photoplethysmogram in frequency domain analysis,” J. Soc. Korean Med. Diagnostics 12, 42–62 (2008).
23. H. Hsiu, C. L. Hsu, C. T. Chen, W. C. Hsu, H. F. Hu, and F. C. Chen, “Correlation of harmonic components between the blood pressure and photoplethysmography waveforms following local-heating stimulation,” Int. J. Biosci. Biochem. Bioinforma. 2(4), 248–253 (2012). [CrossRef]
24. A. A. Fedotov, “Techniques for the morphological analysis of the pulse wave,” Biomed. Eng. 53(4), 270–274 (2019). [CrossRef]
25. F. Kuroda, T. Iwase, K. Yamamoto, E. Ra, and H. Terasaki, “Correlation between blood flow on optic nerve head and structural and functional changes in eyes with glaucoma,” Sci. Rep. 10(1), 729 (2020). [CrossRef]
26. J.-A. Sahel, L. Puyo, M. Fink, M. Atlan, M. Paques, M. Fink, M. Fink, J.-A. Sahel, J.-A. Sahel, J.-A. Sahel, M. Atlan, and M. Atlan, “Waveform analysis of human retinal and choroidal blood flow with laser Doppler holography,” Biomed. Opt. Express 10(10), 4942–4963 (2019). [CrossRef]
27. H. Hou, S. Moghimi, L. M. Zangwill, T. Shoji, E. Ghahari, P. I. C. Manalastas, R. C. Penteado, and R. N. Weinreb, “Inter-eye asymmetry of optical coherence tomography angiography vessel density in bilateral glaucoma, glaucoma suspect, and healthy eyes,” Am. J. Ophthalmol. 190, 69–77 (2018). [CrossRef]
28. J. J. Kaluzny and M. Burduk, “Intereye asymmetry of optic nerve head parameters and retinal nerve fibre layer thickness in patients with open angle glaucoma detected by spectral domain optical coherence tomography,” Ophthalmol. J. 2(2), 35–41 (2017). [CrossRef]
29. A. L. Williams, S. Gatla, B. E. Leiby, I. Fahmy, A. Biswas, D. M. de Barros, R. Ramakrishnan, S. S. Bhardwaj, C. Wright, S. Dubey, J. F. Lynch, A. Bayer, R. Khandelwal, P. Ichhpujani, M. Gheith, G. Siam, R. M. Feldman, J. D. Henderer, and G. L. Spaeth, “The value of intraocular pressure asymmetry in diagnosing glaucoma,” J. Glaucoma 22(3), 215–218 (2013). [CrossRef]
30. F. Moret, C. M. Poloschek, W. A. Lagrèze, and M. Bach, “Visualization of fundus vessel pulsation using principal component analysis,” Invest. Ophthalmol. Vis. Sci.52(8), 5457 (2011). [CrossRef]
31. F. Moret, C. M. Reiff, W. A. Lagrèze, and M. Bach, “Quantitative analysis of fundus-image sequences reveals phase of spontaneous venous pulsations,” Transl. Vis. Sci. Technol.4(5), 3 (2015). [CrossRef]
32. A. Wartak, F. Beer, S. Desissaire, B. Baumann, M. Pircher, and C. K. Hitzenberger, “Investigating spontaneous retinal venous pulsation using Doppler optical coherence tomography,” Sci. Reports 9, 4237 (2019). [CrossRef]
33. K. Gugleta, A. Kochkorov, R. Katamay, C. Zawinka, J. Flammer, and S. Orgul, “On pulse-wave propagation in the ocular circulation,” Investig. Opthalmology Vis. Sci. 47(9), 4019 (2006). [CrossRef]
34. A. Kochkorov, K. Gugleta, D. Kavroulaki, R. Katamay, K. Weier, M. Mehling, L. Kappos, J. Flammer, and S. Orgül, “Rigidity of retinal vessels in patients with multiple sclerosis,” Klin. Monbl. Augenheilkd. 226(04), 276–279 (2009). [CrossRef]
35. K. Zuiderveld, “Contrast limited adaptive histogram equalization,” in Proceedings of Graphics Gems IV (Academic Professional, Inc., 1994), pp. 474–485.
