Time-resolved quantitative inter-eye comparison of cardiac cycle-induced blood volume changes in the human retina

: We describe a low-cost, easy to use binocular instrument to acquire retinal video sequences of both eyes simultaneously. After image registration, cardiac cycle-induced pulsatile light attenuation changes can be measured quantitatively with high spatial and temporal resolution. Parameters such as amplitude, pulse form, and time shift between light attenuation changes can be calculated and compared between eye sides. Deviation from inter-eye symmetry can be not only an early sign of beginning eye diseases such as glaucoma but also a sign of pathological changes in the carotid arteries; hence, this method can improve the early detection of pathological changes. Important features compared to existing monocular instruments are a narrow band light source with the wavelength close to the peak of the blood extinction, and a proportional relationship of image intensity and light intensity, which are the main requirements for quantitative evaluation.


Introduction
The easily accessible retinal blood vessels, capillaries, and tissue offer possibilities for direct observation and a unique possibility to study vascular, neural, and ophthalmic disorders. Different parameters of the human eye change according to the cardiac cycle, e.g., the intraocular pressure (IOP), eye length, and retinal reflection. The pulsatile changing light reflection is due to the pulsatile changing amount of blood in the retinal tissue and vessels [1][2][3][4][5][6][7][8]. The pulse form and amplitude of the pulsatile changing reflection and the pulsatile changing part of the IOP (ophthalmic pulse amplitude) can be related to different factors, e.g., the IOP [4] and carotid stenosis [5,9].
In early studies describing pulsatile retinal reflection, non-imaging methods were used. Using state of the art CCD and CMOS camera technology, the acquisition of video sequences of the ocular fundus is also possible and adds high spatial resolution to the existing temporal resolution. Retinal video recording and image processing methods have been employed to assess different parameters of the eye, e.g., the reaction of the vessel caliber to flickering light stimulation [10,11], to assess spontaneous venous pulsation [8,[12][13][14][15], and to measure flickering light-induced reflectance changes of the human papilla and peripapillary region [16].
In the context of early detection of eye diseases, it is interesting that many parameters in glaucoma indicate inter-eye asymmetry. IOP asymmetry is a common finding in patients with glaucoma [17]. Inter-eye macular thickness asymmetry and total retinal nerve fiber layer (RNFL) thickness are also discussed as promising parameters for identifying early primary open-angle glaucoma [18]. De Leon et al. [19] reported differences of retinal arterial and venular caliber between both eyes in asymmetric glaucoma and Fansi et al. [20] submitted that the rim area to disc area asymmetry ratio (computed from values for the left and right eye) appeared to best distinguish between groups of healthy and three different groups of glaucoma progression. To compare these slowly changing or constant parameters between the eyes, they can be measured individually. The deviation from symmetry can be an early sign of beginning disease.
However, to compare cardiac cycle-induced fast changing parameters such as pulse form, amplitude, and especially time shift between sides, both eyes must be measured simultaneously. Whereas non-imaging binocular instruments were described in early experiments in photometric studies to measure arm-to-retina circulation times in both eyes simultaneously [1,4,21], simultaneous or binocular acquisition is not possible with the current retinal imaging technologies. Owing to their size, existing monocular imaging instruments cannot be combined to measure both eyes simultaneously.
The parallel video acquisition of both eyes simultaneously can enable the comparison of blood-volume-related parameters between both sides with high temporal and spatial resolution for the detection of deviation from symmetry. As the development of diseases such as glaucoma and carotid stenosis probably differs between sides at the beginning, the detection of any asymmetry could be an opportunity for early diagnosis.
In this paper, we describe an instrument for the acquisition of retinal video sequences of both eyes simultaneously with exact synchronization, a method to quantitatively calculate pulsatile attenuation changes from intensity values with high temporal and spatial resolution, and different examples demonstrating the characteristics of the presented method and possible clinical application.

Experimental system
There are certain indispensable requirements for such an instrument. First, it must be sufficiently small to allow both eyes to be measured simultaneously. Secondly, for quantitative measurements of the pulsatile attenuation changes, the illumination must be stable over time, and thirdly, there must be a proportional relationship between the light intensity I in the image plane and the resulting gray value GV of the image (I ∝ GV). Further, there should be high contrast for small changes in blood volume, which means the wavelength has to be selected so that the light absorption of blood is high. For stable fixation during the data acquisition, an effective fixation target must be present. Preferably, the instrument should be inexpensive and for a clinical application, it should be easy to use. The main components of the developed instrument and the principle of evaluation are described in the following sections.

Imaging
The developed instrument was based on the standard principle of fundus cameras [22]. Figure  1 displays the experimental setup. An ophthalmic lens (40D, Volk Optical Inc, USA) formed an intermediate aerial image of the retina in the image plane R2. This image was reimaged by a lens system (2 achromatic lenses 2 x 120 mm focal length, Qioptiq, Germany) to the sensor plane R3 of a CMOS camera (UI-3060 Rev 2, USB 3.0, IDS Imaging Development Systems GmbH, Germany) with 1936 × 1216 pixels. An area of interest of 1000 x 770 pixels was selected to match the field of view of the camera to formerly used CCD cameras with a smaller sensor [22], without changing the optics. With this setting, one pixel corresponded to 1.2 arcmin of the field of view. To ensure a linear transfer function, gamma was set to "1". The black level was automatically adjusted by the CMOS sensor.   3 , which is sufficiently small to place an ophthalmoscope in front of each eye. For alignment to the eyes, the ophthalmoscopes were attached to two independent slit lamp XYZ mounts.

Data acquisition
To ensure exact synchronization of the corresponding images of the video sequences, a hardware trigger for the cameras was employed [27]. Data acquisition software was developed to select the parameters of the cameras (frame rate, exposure time, gain), to move the position of the fixation targets on the OLED displays (to be able to direct the patients viewing direction), and to display live images of both cameras on the screen of a laptop. The data acquisition software executed on an HP ZBook 14 G2. During data acquisition, the time stamp was written to each frame (with millisecond resolution) indicating the moment of image capture by each camera to allow verification synchronization offline.

Image registration
Despite the use of fixation targets, the acquired video sequences indicated eye movements. To compensate for these eye movements, the video sequences were registered offline using a two-step process using the first frame of the sequence as a reference [28]. In the first stage, large movements were corrected via phase correlation, employing a Fourier shift theorem for the estimation of large inter-frame translation. The Lucas-Kanade [29] approach was used in the second stage of registration process to increase its precision. The tracking points were detected on a vascular tree, which was segmented using eigenvalues of the Hessian matrix of the pixel's neighborhood and Otsu thresholding. Small, segmented regions were consequently removed as we assumed that segmented blood vessels create longer structures. Finally, binary skeletonization was used to obtain the points on the vessel centerlines, which were used for tracking. The basic optical flow approach [29] was used for estimation of the global Euclidean transformation matrix for each frame n, employing translation and rotation. The final values of translations (Δx n , Δy n ) and rotations (Δφ n ) were computed for each frame by merging the data from both stages. These extracted signals were then used for the frame registration. Other applications are possible, e.g., eye fixation analysis [30].

Measurements
The study followed the tenets of the declaration of Helsinki for research involving human subjects and informed consent was obtained from all participants of the study. The clinical application of the method was performed in the framework of the Erlangen Glaucoma Registry (EGR), a clinical registry (ClinicalTrials.gov identifier: NCT00494923) for cross sectional and longitudinal observation of patients with open-angle glaucoma or glaucoma suspect, founded in 1991 [31]. Short video sequences of 10 s duration were captured from the ONH of subjects with dilated pupils. The video ophthalmoscope does function with undilated pupils with a minimum pupil diameter of approximately 2 mm. However, the subjects discussed herein were measured with dilated pupils owing to the order of the other clinical investigations. Furthermore, dilated pupils decrease the diameter fluctuations (and consequently intensity changes) during the acquisition and increase the range of tolerable eye movements without disturbing reflection and scattering from the pupil edge. First, the right instrument was aligned to the eye such that a clear image of the ONH was visible on the screen and reflections were removed or minimized by small XYZ movements of the instrument. Then, the subject was asked to maintain fixation, and the left eye was aligned in the same manner.
Finally, the v than 20 s. Aft

Calculati
For consisten expressed as f and an exposu

Pulsati
The temporal (AOI) within  resulting in an attenuation A(n) of the changing part of the blood layer (assuming T + A = 1) The calculated attenuation A(n) of the changing blood layer is independent of the constant parameters of the eye and the instrument (light intensity, camera gain, pupil size, fundus reflection) because changing intensity conditions influence both I(n) and I max . An essential condition for the validity of this equation for the resulting gray values of the images is a proportional relationship between the light intensity I in the image plane and the resulting gray value GV of the image (I ∝ GV).

Pulsatile attenuation amplitude (PAA)
The pulsatile attenuation amplitude PAA representing the maximum change in attenuation A max during one heartbeat can be calculated from Eq.
Equation (4) can also be applied to each pixel of a sequence of frames of an entire heartbeat such that the spatial distribution of PAA(x,y) can be calculated (see Fig. 6). The presented PAA spatial distribution was calculated for five pulses and the resulting images PAA(x,y) were averaged.

Trend correction
If an entire sequence (250 frames, approximately 10 heartbeats) is examined, the constant part of the light intensity could change slowly owing to small eye movements and loss of optimal alignment. Furthermore, the individual heartbeats are not identical and could have marginally different amplitude such that it is not useful to use the maximum value of a single heartbeat as reference, as in Eq.
The changing transmission T trend (n) compared to this trend signal is calculated with resulting in the attenuation

Presentation of the results
All results are displayed in the diagrams in the following manner. The reference value is either the value at the beginning of a heartbeat (I max ) or the average value (I avg ) according to Eq. (3) and Eq. (7), respectively. To avoid negative values and to display the distribution of PAA(x,y) as an image or a video sequence, the value "100" was added to all results. Hence, the reference value is "100" and a value of "105" indicates a change (increase) of attenuation of 5% and "95" a decrease of 5%. Visualization 3 ( Fig. 9) presents, as an example, a video sequence with calculated PAA(x,y,n) values for one pulse. In this case, the reference value is the beginning of the heartbeat at Frame 3 of the video sequence. This video sequence can thus be directly analyzed quantitatively at any area of interest using appropriate software (e.g., the live z-profile of ImageJ). To compare the pulse form between sides, the pulses can also be normalized to their minimum and maximum. If pulses are averaged, the standard deviation is displayed in the diagram. Patient data are covered by gray rectangles in the upper left corner of the images and video sequences.

Results
This section describes the evaluation of a typical measurement and includes specific cases that demonstrate the potential of the proposed instrument. However, the results displayed herein do not yet allow conclusions concerning the correlation of calculated parameters and diseases. Table 1 lists the clinical data of the subjects (peripapillary retinal nerve fiber layer (RNFL) thickness, IOP, and diagnosis). All subjects had clear media. For details on the clinical measurements see [32] and [33].   The PAA in the selected AOI is approximately 8% ( ± 4%); it changes to a small amount from pulse to pulse. The increasing edges of the pulses are considerably similar between both sides. The PAA(x,y) indicates a high value at the neuroretinal rim and also at and beside the vessels (see below).  es at the ues at the Fig. 6(C)) e ONH in Fig. 6(A)   ct with no eye ase, the single n normalized p d in Fig. 7

Compari
The ifferent positio eas with high be selected. ositions and bo eas on retinal v tion can be obs he spatial distr (A) and Fig. 9 however, the sm higher attenua ulses are displa areas for both the different p mpared to the O ian pulses and aque of the carotid yed in Fig. 6 (C)) attenuation A(n) in e pulses (displayed viation (bluish and ons and both PAA, and po Figure 9 disp oth eyes. Corre veins immedia served in this s ribution of the 9(B), respectiv mall AOIs on ation changes ayed in Fig. 9

Discussio
The develope changes of bo There are reason is the light. Recentl been demonst scattering and compression, Figure 6(A and glaucoma neuroretinal intensity chan magnitudes in number of aut RNFL thickn peripapillary e.g [35,36]. speckle flow described, e.g absorption, w The example attenuation th volume during However, It was shown 10 [39][40][41]. Using high-speed spectral domain optical coherence tomography, An et al. [42] found that the ONH experiences continuous oscillatory axial motion that is strongly correlated with simultaneously measured pulsatile blood flow in the central retinal artery. The amplitude of the ONH movement was approximately 4 µm. The increase of retinal thickness caused by expansion of the retinal vessels is in the range of 1 µm as measured by Spahr et al. using phase-sensitive full-field swept-source optical coherence tomography [43]. Schmetterer et al. [44] investigated the distribution of fundus pulsation amplitude in a region of −15° to + 15° around the macula using a laser interferometric method. The average fundus pulsation amplitude in the ONH was approximately 10 µm. These cardiac cycle induced sample movements in conjunction with the directionality of the tissue scattering could result in cardiac cycle induced intensity changes. Additionally, changes in scattering due to changes of orientation of blood constituents (e.g. red blood cells) over the cardiac cycle could also result in intensity changes.
Further work is necessary to clarify the contribution of the different effects described here on the origin of the cardiac cycle induced attenuation changes.
Changing attenuation due to the pulsation of large retinal blood vessels and even small vessels, which are visible owing to their changing blood volume or their movement, can also be observed. This is documented in Fig. 9. If the vessel diameter changes according to the cardiac cycle, there are areas with high PAA along the vessels. This is caused by the fact that the vessels are darker than the background and positions near the vessels have high reflection when the vessel caliber is small (at the beginning of the heartbeat) and become darker when the vessel caliber increases (at the peak of the heartbeat). This can be observed in the vein depicted as Vessel 1 in Fig. 9. In this case, the vein exhibits large diameter pulsation causing PAA changes around the vein.
Furthermore, vessels that bend and change their position, to a small degree, based on the cardiac cycle also change the image intensity and values of PAA(x,y) at corresponding locations. Reduced PAA(x,y) values result in areas where the vessels are located at the beginning of the heartbeat and high values where they are located at the peak of a heartbeat. This can be observed in Vessel 2 in Fig. 9(B) and the corresponding video sequence (Visualization 2). Whereas these results can help to locate moving vessels, it must be noted that precise image registration is necessary for this evaluation; otherwise, artifacts are created. Although in many cases the registration provides satisfactory results, other cases could fail for specific frames (due to eye movement, saccades, and blinking) [28]. The amount of vessel movement and pulsation has been observed as being highly variable among subjects. Two examples presented herein demonstrate a significant amount of vessel movement in Visualization 2, yet only minimal movement in Visualization 1.
Beyond the evaluation of PAA, the time shift of pulse parameters (i.e., the position of minima or maxima) between sides could be an important parameter to detect deviation from normal values due to beginning diseases. For the subject with macroangiopathy and plaque of the carotid communis displayed in Fig. 8, a time shift of pulse parameters between sides can be clearly observed. Carotid atherosclerosis could cause different ocular symptoms [45] and carotid stenosis has been proven to influence the arm to retina circulation time (ARCT), particularly ARCT delay, disparity of ARCT between the eyes, and delayed ciliary circulation at the optic disc [46]. As the level of carotid stenosis is not symmetric, it leads to an impaired blood flow in the ophthalmic artery and the posterior ciliary arteries causing differences in ocular pulse amplitude [3] and possibly differences in dynamic patterns in corresponding retina.
Another noticeable time shift was observed between A(n) of the entire ONH and the small AOIs placed on retinal veins (Fig. 9). We found these delays to be in the range of 2-3 frames, corresponding to 80 ms to 120 ms. Similar delays regarding the spontaneous vein pulsation phase following the arterial phase have been described in Morgan et al. [12] and Moret et al. [47].
The assumed time shift indicated in Fig. 8 is in the range of one frame. To better detect a time shift, the temporal resolution of the instrument could be improved by increasing the frame rate and reducing exposure time. In this case, the light intensity also must be increased. This can be increased up to a factor of "7" without leaving the safety range of retinal illumination for unlimited exposure.
Compared to existing (monocular) video ophthalmoscopes that are based on (modified) standard fundus cameras or slit lamps with digital SLR (single-lens reflex) cameras (e.g., Canon, Nikon) as the image detector and tungsten halogen lamps as the light source, e.g [12,13], the instrument proposed herein has several advantages. 1) The contrast for small blood volume changes is 60% higher because the narrow spectrum of the LED coincides with an absorption maximum of blood. Conversely, in commercial instruments used for video acquisition, a tungsten halogen lamp is used as the light source and the green channel of the SLR cameras as the detector. Hence, parts of the spectrum that are only marginally influenced by the light absorption of the blood reduce the contrast. 2) The intensity of illumination is significantly reduced because the selection of the wavelength is performed in the illumination beam pathway and not in the imaging beam pathway. In fundus camera-based systems, the fundus is illuminated with the entire spectrum and only a relatively small range of the spectrum contributes to the changing light absorption. 3) Digital SLR cameras have a nonlinear transfer function that is optimized for color photography, not for quantitative measurements. There is no proportional relationship between light input and the resulting image intensity. Even logarithmic correction does not lead to an exact proportional relation because the transfer function is typically not described by a power law expression, Output = Input γ , as it is designed for the "gamma correction" of computer monitors. Due to the nonproportional relationship, the calculated light attenuation depends on the light intensity. 4) The spectrum of tungsten lamps changes if the brightness is changed, e.g., for setting the image intensity. For quantitative evaluation (comparison of subjects, repeated measurements of the same subject), the tungsten lamp must be used at a fixed intensity level. 5) Over the data acquisition time of 10 s, an internal fixation ensures more stable fixation compared to external fixation with the contralateral eye.
It must be noted that the presented examples demonstrate the unique characteristics of the instrument and possible clinical applications; however, no conclusions should be drawn from the limited sample size. In the future, this instrument will be applied in a clinical study to assess normal values of different parameters (PAA, pulse form, and time shift) in a larger population to compare the results from patients with the results from normal subjects. It will also include other clinical parameters such as visual field, nerve fiber layer thickness, IOP, and the condition of the carotid arteries.

Conclusion
Because of the narrow band light source with a wavelength close to the peak of the blood extinction, proportional transfer function, internal fixation, and robust registration software, cardiac cycle-induced light attenuation changes caused by pulsatile changing blood volume can be quantitatively measured at any area of interest within the field of view of the video sequences. The parallel acquisition allows precise comparison of the parameters of pulsatile absorption changes and determination of the time shift between both eyes. The developed instrument is small, can function on battery power, and is low-cost and easy to use. It will be applied in a clinical study to determine the correlations of the new parameters, PAA, amplitude, and time shift with diseases such as glaucoma and stenosis of the carotid arteries.