Assessing blood vessel perfusion and vital signs through retinal imaging photoplethysmography

: One solution to the global challenge of increasing ocular disease is a cost-eﬀective technique for rapid screening and assessment. Current ophthalmic imaging techniques, e.g. scanning and ocular blood ﬂow systems, are expensive, complex to operate and utilize invasive contrast agents during assessment. The work presented here demonstrates a simple retinal imaging photoplethysmography (iPPG) system with the potential to provide screening, diagnosis, monitoring and assessment that is non-invasive, painless and radiationless. Time series of individual retinal blood vessel images, captured with an eye fundus camera, are processed using standard ﬁltering, amplitude demodulation and principle component analysis (PCA) methods to determine the values of the heart rate (HR) and respiration rate (RR), which are in compliance with simultaneously obtained measurements using commercial pulse oximetry. It also seems possible that some information on the dynamic changes in oxygen saturation levels ( SpO 2 ) in a retinal blood vessel may also be obtained. As a consequence, the retinal iPPG modality system demonstrates a potential avenue for rapid ophthalmic screening, and even early diagnosis, against ocular disease without the need for ﬂuorescent or contrast agents. compensation techniques.


Introduction
Ocular diseases (OD), including glaucoma, diabetes, cataracts and age-related macular degeneration (AMD) have experienced a rapid increase worldwide [1,2]. Studies drawing correlations with age, gender, genes, diet, lifestyle and exposure to sunlight, have neither been able to clearly indicate the causes of ODs, nor predict their development [3]. The ocular microcirculatory system consists of a complex network of arteries and veins driven by ocular perfusion pressure (OPP) [4]. Retinal vessel branches are connected to the optic nerve system via the lamina cribrosa structure (LCS), which is meshed with the central retina vein (CRV) and central retina artery (CRA), and so interact with the cardiovascular system. The intervals between heartbeats create fluctuations in the blood flow within the retina, attributed to continuous changes in the sympathetic-parasympathetic balance of the autonomous nervous system (ANS) [5]. We note from previous studies [6,7] that the quantification of such fluctuations is a vital factor in early detection of visual impairments resulting from cardiovascular disease, hypertension, diabetic autonomic dysfunction and some psychological disorders.
In recent decades digital ocular imaging techniques, associated with a number of screening modalities, have evolved in an attempt to satisfy the growing demand for a method of direct visualization of the vascular network which might allow for effective diagnosis and assessment in the prevention of OD development. These modalities tend to rely heavily on fundus photography (FP) and fluorescein angiography (FAG) as clinical gold standards for the direct assessment of the retinal vasculature [8,9]. FP is usually used to photograph the interior surface of eye using fluorescent dye on the retinal area, while FAG is used to determine the arteriovenous transit time of fluorescent dye within the retina. However, these two clinical bench techniques are still incapable of assessing changes in retinal vasculature. Moreover, the different manufacturing settings of these screening devices (FP and FAG) may result in inconsistent measurements, causing an unreliable assessment of such changes [10]. Imaging photoplethysmography (iPPG), an emerging technology that is being developed together with associated opto-physiological models [11], has the potential to provide an alternative, non-invasive and real-time solution for retinal vasculature assessment through the advances in photonics technology, including those of smartphone cameras. Recently, for example, optical imaging assessments have been applied in an ophthalmological clinical study to quantify blood volume changes in veins influenced by intraocular pressure (IOP) [12,13] with the aim of discovering phase correlations between intracranial pressure (ICP) and retinal vessel pulsation (RVP). Two further assessment studies [14,15] were carried out to investigate alterations to vascular pulsation properties due to glaucoma, ICP changes and retinal venous occlusion. Aligned 2D matrices from a particular region of interest (ROI) were generated by synchronizing an ophthalmodynamometer and a pulse oximeter. These ROIs, under green illumination, were segmented into grey scale images before extracting the pulsatile signals. The phase relations between retinal arteries and veins were quantified from peak (dilation) and trough (constriction) points, and found to be consistent with the cardiac cycle. Furthermore a USB-powered, lightweight ophthalmoscope connected to an optical coherent tomographic (OCT) system has also been developed to improve early screening of glaucoma and associated eye diseases [17]. The assessment was implemented with a 5mm diameter LED [16,17] to illuminate a dilated pupil which enabled the recording of a video at 25 frames per second. A two-stage procedure of frame-by-frame registration processing was adopted to overcome significant eye movements in the first stage, and a sub-pixel resolution in the second stage, allowing the pulse parameters and fixation pattern to be defined. Hence, the heart beat-induced reflection changes were easily measured and might be used to assist an ophthalmologist to reveal structural changes of the optic nerve head.
Incorporating the well documented challenges of motion artifacts [6,18], recent study has focused on how one might most effectively provide reliable vital physiological parameters derived from retinal blood pulsation changes [12,13,19] by the means of retinal iPPG. In the current work, the optical measurement settings for retinal vascular assessment and validation of the measured retinal iPPG have been designed using a protocol which correlates pulsatile changes aligned with cardiac cycles. Such changes were recorded simultaneously with a portable fundus camera, a contact photoplethysmographic (cPPG) sensor and a commercial pulse oximeter through the implementation of the designated protocol with five healthy participants. Here 2D matrix ROI tracking allowed sub-pixel motion compensation, using all raw RGB image channels captured from the fundus camera. The pulsatile waveforms from the series of video images are extracted through an algorithm which combines feature tracking, motion compensation and noise reduction. From these, the vital physiological parameters (heart rate (HR), respiration rate (RR) and oxygen saturation (SpO 2 )) in the eyes are derived. The study indicates that a practicable engineering solution, using existing and forthcoming photonics-based iPPG, may have the potential to provide an alternative modality for eye disease prognosis and rapid screening. Further details on the signal processing methods used to extract the imaging and cPPG signals and a tissue optics modeling framework, based around Radiative Transfer Equation (RTE), which allows further investigation of the suitability of this method as a diagnostic technique will be published separately.

Experimental setup
A schematic of the retinal assessment system is shown in Figs. 1 and 2. The system is composed of: 1) a Visuscount 100 (Carl Zeiss Meditec AG) fundus camera, mounted on a slit lamp with adapter assist, to record the fundus videos; 2) a commercial pulse oximeter sensor (CMS50DL, Contec Medical Systems Co. Ltd., China) connected via an e-health shield connection bridge to an Arduino UNO board (Arduino, Spain). Ethical approval of the experimental protocol was obtained (from Loughborough University Ethical Committee) to carry out retinal assessment on five healthy participants under the supervision of an ophthalmologist. Following pupil dilation, the video assessment was taken with the fundus camera at a rate of 10 frames per second (fps) and then transferred via a USB port to a laptop. Whilst recording the video, cPPG and oxygen saturation (SpO 2 ) signals were simultaneously recorded from, respectively, the PPG sensor amped (Pimoroni Ltd., UK) placed on the index finger of a participant's right hand or earlobe and from a pulse oximeter sensor placed on the index finger of participant's left hand. Data acquisition from both sensors was implemented on Arduino platforms (Arduino v1.0.5 and v1.6.6). The SpO 2 reading from the pulse oximeter and the PPG signals from the pulse sensor were amplified and then transferred to an open source serial port application running on a laptop for further signal processing. The image processing of the fundus video was performed in Matlab (2015a, MathWorks Inc., USA) using a number of built-in functions. The fundus camera was of CMOS design with 5 MP resolution, a ±20 Dioptre lens and 40°FOV. Level 2 to 3 brightness was used during the assessment to reduce eye motion and blinking. The images captured were presented on an 8-bit scale (0 to 255) RGB and exported in TIFF format.

The cPPG signal, c(t)
The contact photoplethysmography (cPPG) signal from the finger pulse device was used both as a validation of the final imaging signal and as a reference to assess the level of signal processing that is suitable for the current analysis. The formation of an iPPG signal depends on blood opacity and the cardiac cycle activity of the autonomous nerve system (ANS). The cardiac cycle and respiratory activity modulates the backscattered light through blood volume changes in the illuminated blood vessel, and this may be observed as a modulation of the blood opacity (intensity) measured. This picture is complicated somewhat by respiratory activity as, first, intra-thoracic pressure variations cause an exchange of blood between the pulmonary and the systemic circulation which results in a variation of perfusion baseline which adds an offset to the heart rate signal. In addition, a decrease in cardiac output due to reduced ventricular filling causes pulse amplitude variations that provide a respiratory induced amplitude modulation (AM) to the heart rate signal [20]. Here the former is removed by filtering and the latter by standard demodulation methods.
Following the recording of the contact (cPPG) signal, a Band Pass filter was used to remove out-of-band noise resulting from the lower frequency respiratory system and higher frequency noise. This filtered output is shown in Fig. 3(a). The resulting signal also suffers respiration induced amplitude variations as indicated above, and the upshot is the true cPPG signal c(t) is contaminated by a low frequency modulating signal m(t), to give This leads to in-band noise, which thus will have escaped the Band Pass filtering. The modulating signal m(t) can be removed by squaring x(t), low pass filtering the result and taking the square root. The idea is that, for a single component f m of m(t) and f c of c(t), the square has components at f = 0, 2 f m , 2 f c and 2( f c ± f m ). Low pass filtering, with a cut-off between 2 f m and 2( f c − f m ), leaves only the DC and 2 f m components. Taking the square root, the result is then c(t) to a good approximation. The modulating signal m(t) obtained in this way is shown in Fig.3(b). The signal can then be divided out as in Fig.3(c) and (d). A fuller explanation of this simple idea may easily be extended to include a sum of frequency components (for each of m(t) and c(t)) as a Fourier series.

The imaging PPG signal, s(t)
The video recording captured on the Fundus camera was imported into the MATLAB environment (in .avi format) where a computer vision cascade object detector classifier using the built-in Kanade-Lucas-Tomasi (KLT) algorithm was used to define the location of a Region of Interest (ROI). A group of features was selected in the first frame and the feature points tracker was used to track their positions through the video sequence. This allowed the video to be cropped to 32 × 32 pixels, positioned around a particular ROI in 8-bit RGB tiff format.
The red channel might be expected to suffer from saturation issues due to the way that MATLAB's colormaps work, and it may also contaminate the blue channel. Consequently, we largely concentrate on the green channel which was used to increase focus on the blood vessel highlighted in the ROI. The green component from two frames from one of the video sequences is shown in Fig.4 as both contour and surface plots. The dark (blue) feature in the surface plots corresponds to the blood vessel of interest, in this case a retinal vein. Despite the KLT tracker, there is still a certain amount of movement within the 32 × 32 block. As a result, we automated a search for blood vessel pixels by finding the 3 × 3 pixel mini-block which has the lowest mean pixel value, and this mean value was taken to define an iPPG signal g(t) and its location to define a pixel within the vein. The mean pixel value, of the same mini-block, in the red and blue channels then defines red and blue signals r(t) and b(t).
Following their generation, out-of-band and in-band noise was removed from the RGB signals r(t), g(t) and b(t) using the same methods as for the contact PPG signal, with the additional step that, to combine the three colored channels into a single signal and to remove additional in-band noise, we also implemented a principal component analysis (PCA). First the RGB components are z-scored by subtracting their mean value and normalizing by their standard deviation. The MATLAB pca function then returns principal components of the form The coefficient matrix for the PCA analysis is given in Table 1 and demonstrates as expected that the noisier red channel is reduced in importance in the first principal component (PC A 1 ) compared to the blue and green channels which are roughly equally weighted. The PCA analysis also has the additional advantage of removing some of the uncertainty in the camera sensitivity across the RGB spectrum. The second and third principal components were statistically tested for normality, as were the  differences of the z-scored red and green, and blue and green channels. In each case (using one and two sided Kolmogorov-Smirnov and F-tests) the null hypotheses that these distributions are drawn from the same normal distribution could safely be accepted at the level of 5% significance. This supports the assumption that the measured RGB signals carry the same PPG information together with some additional, contaminating Gaussian noise. The combined signal processing steps for the iPPG signal were: band pass filtering, principal component analysis and demodulation. The first principal component signal is shown in Fig. 5 where it is compared to the contact PPG signal from Fig.3(c).

Heart rate (HR)
Unfortunately, our simple system did not have the capacity to trigger the camera and the contact sensor to record from the same time reference; which leaves a certain amount of uncertainly in the line-up of the iPPG/PCA and the cPPG signals. However, there is clearly enough commonality in Fig.5 to observe these signals show a similar cardiac activity pattern; one located at a finger tip, and the other in a retinal vein. The sample frequency of the camera (10 frames per second (fps)) means that there are only around 8-9 sample points per cardiac cycle, compared to the 35-40 for the contact sensor. This makes frequency resolution, and the application of the filtering steps, more of a challenge; it also introduces limitations on the accuracy of information extracted from the signal. For example, a peak detector on each signal should reveal a dynamic heart rate, however the peak position is limited to an accuracy of ± half a same period, i.e. 0.05s for the iPPG signal. Similarly, it is not always clear whether a peak is a new cycle or some mid-cycle feature. A greater problem is that, in a very small number of frames, the motion compensation techniques are unable to track the vein under investigation correctly which results in the vein disappearing from the 32 × 32 block, in Fig. 4(a) or (c). The heart rate (HR) extraction from a peak detection of the two techniques (imaging and contact) are shown in Fig.6. The differences reflect the poorer frame rate, and it attendant problems, as well as occasional data point from the Fundus camera suffering from such motion artifacts. The data here is presented without any attempt to remove any such erroneous points.

Oxygen saturation (SpO 2 )
In principle it is possible to get an idea of the continuous variation in SpO 2 using the red and green signals, r(t) and g(t). The AC component, for each, can be obtained by using a peak-finder to record the position of the peaks and troughs in the plots, such as those shown in Fig.5. Resampling these curves at the sampling frequency of 10 fps and subtracting amplitudes, leads to a continuous sampling of the AC amplitudes for red and green, as in Fig.7. An assessment of oxygen saturation can be calculated as where A and B are fitting parameters. The derivation of Eq. 2 assumes that the camera signal provides a faithful (linear) representation of the light intensity. It is worth observing that many cameras employ a brightness resolution function (which also varies with wavelength,λ) so the camera signal V is a nonlinear function of intensity I, V = F(I, λ). However, as the DC component dominates the light intensity, the AC components may be regarded as small increments, thus the AC increments in camera signal and true intensity are related by δV = F (I, λ)δI, where F (I, λ) is the first derivative of F(I, λ) with respect to I. This only alters Eq. 2 by introducing a scale factor of F (I gr een , λ gr een )/F (I r ed , λ r ed ); a scaling which will be incorporated in the fitting parameter B. Essentially linearizing F(I, λ) around the green and red DC values simply modifies the fitting parameter. The upshot is that provided a temporal variation is discernible in the camera signal, any brightness resolution function will be accounted for in the fitting.
We compare the behavior of the venous SpO 2 extracted from the imaging signal to that (systemic value) generated by the contact sensor in Fig.8, although to facilitate comparison of the temporal variation in the two cases we have adjusted the value of A to fit the systemic response. The pulse oximeter observation is only able to provide SpO 2 levels with 1% accuracy, however the main features appear to be similar in the two cases: a dip after 1s, followed by a broad peak, a second dip after 11s, and a double peak in the interval 14-16s.

Respiratory rate (RR)
In order to obtain information of the respiration rate (RR) the passband of the filter, designed for HR extraction, was replaced by one with band edges at 0.1 and 0.4Hz (corresponding to 6 to 24 breaths/min), and a PCA signal was extracted as described in Section 2.3. As in Fig.5 a peaks/zeros/troughs detection scheme was employed to assess the respiration rate (RR). The result, shown in Fig.9, shows a slightly variable respiratory rate around a mean of roughly 20 breaths per minute. Fig.9 also compares this PCA respiration function to the two respiration controlled contributions of the cPPG signal, the low pass filtered signal RR c P PG (t) and the extracted envelope function (described as m(t) in section 2.1).

Blood vessel diameter d bv (t)
As something of a final sanity check, we use the surface/contour plots of the image sequence, such as those shown in Fig.4 to investigate the change in the diameter of the blood vessel. The current method of defining the PCA signal from the minimum average pixel value of a 3 × 3 mini-block tested across the 32 × 32 ROI regions, measures the time varying blood flow as the blood volume increase in the direction normal to the image. The same blood volume increase is    also identifiable in the in-plane direction perpendicular the line of blood vessel. Assuming that the position of the vessel is stationary relative to the camera, an edge detection of the sides of the blood vessel can be obtained from the contour plot. One would expect that the variation in the blood vessel diameter about some mean d 0 should, roughly, be proportional to the PCA signal s(t) as both are measuring the same increase in blood volume. Thus and consequently it is expected that the peaks in the variation of the blood vessels diameter should line up with those in the imaging PCA signal s(t), shown in Fig.5. A comparison of the two is shown in Fig.10.

Discussion
The aim of this study was to develop a protocol and equipment settings, a tracking ROI system and appropriate post signal processing in which vital signs (HR, RR and potentially SpO 2 ) are assessed (and which may be mapped) for healthy participants' eyes by adopting a non-invasive iPPG technique using a portable fundus camera at a 10 fps sampling rate.
In previous clinical studies, intra ocular pressure, (IOP, mmHg) and intra cranial pressure, (ICP, mmHg) variation have been observed and monitored for their relationship to the dynamic changes of blood vessel diameter (arterial or venous) during the cardiac cycle. Such previous studies have demonstrated various modalities for retinal screening and imaging, and the common ophthalmologic practice is to use fluorescence angiography. Other studies have observed the phase different between IOP and ICP signals using cardiac signals [12,13]. Some researchers have claimed that these signals were not in phase with cardiac cycle, while others have shown the opposite findings [19]. The relationship between retinal pulsation and IOP and ICP signals still remains uncertain.
Nevertheless it may be possible to identify ocular diseases such as hypertension, cardiovascular diseases, diabetic retinopathy, brain tumors, AMD and glaucoma using such measurements. Clearly early screening is an effective way to prevent severe damage to ocular tissue and to allow a suitable and specific treatment to be applied for these patients. The outcomes from the present study have demonstrated that vital signs may be obtained from iPPG based ophthalmologic measurement and are in agreement with cardiac cycle activity. The retinal pulsation extracted from these iPPG signals has been used to derive HR, RR and to investigate SpO 2 . This method may also be used to map variations of retinal blood perfusion in the retinal vascular network. A number of potential experimental issues remain: the pulsatile variation extracted from the iPPG signals is generally contaminated with motion artifacts; illumination is slightly uneven; electrical noise is generated by the fundus camera (which also has a low frame rate) and the photo sensor has a relatively low sensitivity.
The impact of motion contamination may be seen most clearly from the Bland Altman plot of the data from the five participants, Fig. 11. The plot shows the full data which has been subject only to the filtering, demodulation and PCA analysis steps, i.e. no data has been removed. Some eye motion artifacts still remain in the data (typically due to either blinking, rolling of the eyes or the blood vessel motion leaving the 32 pixel ×32 pixel window). An algorithm for excluding such events without impacting on data validity will be included in a future development. Motion such as this generally affects more than a single frame, leading to a small number of erroneous data points. This can cause the apparent introduction of two short cardiac cycles when there should only be a single longer one (causing HR doubling), or can cause the conflation of two cycles giving one extra-long cycle (HR halving). However, despite these limitations, which can be overcome, the framework of the retinal iPPG system presented in this study is potentially able to extract blood perfusion in retinal blood vessels and to interpret the vital signs (HR, RR and SpO 2 ) without the need to perform any angiographic procedures on the participant's eye. Fig. 11. Bland Altman representation of Heart Rate data obtained from imaging PPG compared to the value obtained from the contact sensor.

Conclusions
This paper attempts to address the need for a reliable and cost effective method of assessment and monitoring of retinal health status. Our aim has been to understand to what extent it is possible to access retinal blood vessel perfusion and physiological vital signs using only a simple fundus camera and inexpensive optical equipment. We have successfully developed appropriate signal processing methods and validated the final results using a cPPG signal obtained from a commercial pulse oximeter together with an PPG sensor. Our signal processing is simple and targeted to the removal of respiratory effects and contaminating noise. We use a band pass filter to remove out-of-band noise, amplitude demodulation to remove in-band respiratory effects (crucial for the cPPG signals to regularize the amplitude of the signals), and a principal component analysis s(t) to combine the signal content in the RGB outputs of the fundus camera. We justify the PCA analysis by performing a series of statistical tests to determine that the neglect of other principal components is equivalent to removal of in-band additive Gaussian noise. We are able to show a cardiac rhythm which is highly correlated with that from the pulse oximeter. However our automated assessment currently tends to underestimate Heart Rate (HR) given by the pulse oximeter. This is partly due to the low frame rate of the present fundus camera which contributes too few samples to each cycle for accurate timing information, and is partly due to the relatively simple motion compensation techniques which, on occasion, are unable to follow eye motion leading to erroneous signal points.
We further demonstrate how, in principal, it could be possible to measure SpO 2 levels, and to that extent compare our derived estimates with the discrete results (1% accuracy) from the pulse oximeter. The respiration rates extracted from iPPG signals are quite believable. In addition the dynamic variation of the width of the footprint that the blood vessel makes in the 32 × 32 pixel block agrees well with the PCA (iPPG) signal s(t). We are currently in the process of making the system more accurate through a using a camera with a higher sample rate (over 50 fps) camera and better sensitivity, and more sophisticated motion compensation techniques.