1.1. Aim of Thesis

1.2. Rationale for Non-Contact Heart Rate Monitoring

1.3. Importance of the Heart Rate

1.4. Fundamental Theory of Photoplethysmography

Figure 1.1. AC component of PPG signal compared to ECG. Allen (2007). Reproduced by permission from IOP Science. 

Figure 1.1. AC component of PPG signal compared to ECG. Allen (2007). Reproduced by permission from IOP Science. 

Figure 1. 2. AC and DC components of PPG signal.

Figure 1. 2. AC and DC components of PPG signal.

1.5. Pulse Oximetry—An Example of Contact PPG

Figure 1. 3. Transmission (left) and reflection (right) mode contact pulse oximetry. Tamura et al. (2014).

Figure 1. 3. Transmission (left) and reflection (right) mode contact pulse oximetry. Tamura et al. (2014).

1.6. Non-Contact Photoplethysmography

Figure 1. 4. Representative setup for video-based rPPG. Poh et al. (2010).

Figure 1. 4. Representative setup for video-based rPPG. Poh et al. (2010).

1.7. Problems with Video-Based rPPG

1.8. Research Objectives

1.9. Thesis Outline

2.1. Previous Work

3.1. Research Materials

3.2. Experimental Setup

3.3. ROI Description and Classification for Associated rPPG Signals

Figure 3. 1. Various ROIs used for rPPG signal extraction. Clarification: ROI B represents the union of ROI A and the small rectangle above ROI A.

Figure 3. 1. Various ROIs used for rPPG signal extraction. Clarification: ROI B represents the union of ROI A and the small rectangle above ROI A.

Table 3. 1. ROI description and rPPG signal classification.

3.4. Mean Error, RMSE and SNR Definitions

4.1. Overview of Workflow

Figure 4. 1. Implementation workflow.

Figure 4. 1. Implementation workflow.

4.2. Video Input

4.3. ROI Selection and Tracking

4.3.1. Face Detection

Face detection is a type of object-class detection. In object-class detection, the objective is to identify every instance of an object in an image belonging to a particular class. Typical objects include pedestrians and cars. Face-detection algorithms are most often designed for the detection of frontal human faces. A widely used algorithm is the Viola-Jones face classifier (Viola and Jones 2001), which is the default option for the Cascade Object Detector in MATLAB’s Computer Vision Toolbox. The Viola-Jones algorithm is known for its robustness (very high detection rate), efficiency and utility in real-time applications. It is made up of four main components:

• Haar-like feature selection
• Integral image creation
• AdaBoost training
• Cascading classifiers

Figure 4. 2. Bounding box of face detected by Viola-Jones algorithm.

Figure 4. 2. Bounding box of face detected by Viola-Jones algorithm.

Haar-like features are digital image features used in object detection. They are so named because of their resemblance to Haar wavelets – a series of rescaled square-shaped functions which together comprise a wavelet family.

Figure 4. 3. A Haar wavelet. Wikimedia Commons contributors / CC-BY-SA-3.0.

Figure 4. 3. A Haar wavelet. Wikimedia Commons contributors / CC-BY-SA-3.0.

A Haar-like feature can be imagined as a crude approximation of the intensity profile of an image or image sub-window and attempts to capture the essential appearance attributes of an object class. For example, all human faces contain some commonalities which may be matched using Haar-like features. These include:

• Eye region darker than upper cheeks
• Nose bridge region brighter than eyes

The use of Haar-like features to evaluate these commonalities is illustrated in Figures 4.4 and 4.5

Figure 4. 4. Haar-like feature that looks like bridge of nose applied to face.

Figure 4. 4. Haar-like feature that looks like bridge of nose applied to face.

Figure 4. 5. Haar-like feature that looks like eye region applied to face.

Figure 4. 5. Haar-like feature that looks like eye region applied to face.

Haar-like feature evaluation is based on the principle of computing the difference between the sum of pixels in adjacent white and black rectangles in an image i.e., value = Σ (pixels in black area) – Σ (pixels in white area). In the original implementation of the Viola-Jones algorithm, feature evaluation was performed on 24 x 24-pixel sub-windows of an image (Viola and Jones 2001). Though the Haar-like features were limited to four types (shown in Figure 4.6), at all possible scales and positions, these features required more than 160,000 evaluations over the sub-window (far more than the number of pixels, which is 24 x 24 = 576). At 160,000 evaluations per sub-window, the summing over pixels required for each feature evaluation is computationally very expensive and of limited use in real-time applications. Viola and Jones proposed an novel solution to this problem by introducing the integral image, more generically known as the summed-area table (Crow 1984).

Figure 4. 6. Haar-like features used in the original implementation of the Viola-Jones algorithm. Prmorgan / Wikimedia Commons / Public Domain.

Figure 4. 6. Haar-like features used in the original implementation of the Viola-Jones algorithm. Prmorgan / Wikimedia Commons / Public Domain.

The integral image is a data structure which allows rapid computation of the sum of pixels in an image or image sub-region in constant time. The pixel intensity I(x,y) at any point (x,y) in an integral image is the sum of the intensity i(x,y) of all pixels above and to the left of (x, y), inclusive, in the original image. This may be mathematically expressed as:

$$\mathrm{I}\left(\mathrm{x},\mathrm{y}\right)=\sum _{x\text{'} \le x, y’ \le y}i\left(x^{\text{'}},y^{\text{'}}\right)$$

The integral image can be computed in a single pass over an image since the value I(x, y) in the integral at (x, y) is defined recursively as:

$$\mathrm{I}(\mathrm{x},\mathrm{y})=\mathrm{i}(\mathrm{x},\mathrm{y})+\mathrm{I}(\mathrm{x},\mathrm{y}-1)+\mathrm{I}(\mathrm{x}-1,\mathrm{y})-\mathrm{I}(\mathrm{x}-1,\mathrm{y}-1)$$

Once the integral image has been set up, computing the sum of pixels for any rectangular sub-region of an image requires only four references (independent of the size of the sub-region). For a rectangular sub-region defined by coordinates A, B, C, D, the sum of pixels is given by:

$$\sum _{x_o \le x\le x_{1,}y \le y\le y_{1,} }i(x,y) = I(A) + I(D) - I(B) - I(C)$$

where A=(x₀, y₀), B=(x₁, y₀), C=(x₀, y₁) and D=(x₁, y₁)

Figure 4. 7. Intensity values of an image and its integral counterpart. I(a) + I(d) – I(b) – I(c) = 20 + 74 – 39 – 33 = 22.

Figure 4. 7. Intensity values of an image and its integral counterpart. I(a) + I(d) – I(b) – I(c) = 20 + 74 – 39 – 33 = 22.

Even with the computational efficiency afforded by the integral image, evaluating 160,000+ Haar-like features for each 24 x 24 sub-window is still computationally prohibitive for real-time facial detection. Most Haar-like features are useless for a given object class. Therefore, a method for selecting the most useful Haar-like features could significantly improve the time complexity of the detector. Viola and Jones achieved this selectivity by using a variant of a pre-existing algorithm known as AdaBoost (Freund and Schapire 1997).

AdaBoost, short for Adaptive Boosting, is an example of a boosting algorithm – a class of machine-learning algorithms than can be used to improve the performance of other so-called weak learner algorithms. Within the Viola-Jones algorithm, Adaboost both identifies the weak learners and constructs a strong learner from a linear combination of weak learners. Each weak learner is based on one of a set of best-performing Haar-like features (i.e., those that yield consistently high detection rates when applied to a face). With a boosting algorithm, the individual weak learners need only perform slightly better than random guessing to converge to a strong learner. A weak learner h(x,f,p,θ) evaluated over a 24 x 24-pixel sub-window x consists of a feature (f), a threshold (θ, which decides whether x is classified as face or non-face) and a polarity (p, indicating the direction of the inequality). Both p and θ are determined during training.

$$\mathrm{h}(\mathrm{x},\mathrm{f},\mathrm{p},\mathsf{\theta })=\left\{\begin{array}{c}1 \\ 0\end{array}\right.\genfrac{}{}{0pt}{}{ \mathrm{i}\mathrm{f}pf(x)<p\theta }{\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}}$$

A strong learner H_T (a linear combination of weak learners) selected over T trials of testing may be mathematically defined as:

$${\mathrm{H}}_{\mathrm{T}}(\mathrm{x})=\sum _{i = 1}^T{h}_i(\mathrm{x})$$

where h_i represents the i^th weak learner in the T trials.

Apart from determining the best Haar-like features, the modified AdaBoost algorithm also sets the polarity and threshold for each weak classifier. There seems to be no satisfactory solution to this problem except brute force i.e., the determination of each new weak classifier in the T trials of testing involves evaluating each feature on all the training examples (at different thresholds and polarity). The best performing feature in a testing trial is the feature which minimizes the weighting error it produces. This weighting error is a function of the weights assigned to the training images. The weight of a correctly classified image is decreased, while that of a misclassified image is kept constant. As a result, the i^th feature added to the classifier is forced to “focus” harder on the test images misclassified by the (i – 1) previous features. The weighting system is, therefore, an integral aspect of AdaBoost.

Table 4. 1. Modified AdaBoost algorithm (Viola and Jones 2001).

Table 4. 1. Modified AdaBoost algorithm (Viola and Jones 2001).

With the integral image and modified AdaBoost, the Viola-Jones algorithm seems ready for facial detection. However, computations thus far have focused on a single 24 x 24-pixel sub-window. Most images have higher resolution, meaning that the Viola-Jones detector would have to evaluate every 24 x 24-pixel sub-window of that image. This is still computationally too expensive and also very inefficient as it is known in advance that most sub-windows contain no face. This calls for a method that can streamline evaluation by discarding sub-windows with no face with minimal computational effort and spending more time on potentially positive windows. Viola and Jones resolved this problem through yet another novel contribution known as the attentional cascade.

An attentional cascade or cascaded classifier is a multi-layered classifier architecture in which each level of the classifier consists of a strong classifier prepared by the modified AdaBoost. The objective at each layer is to determine whether a given sub-window in an image is not a face or possibly a face. A sub-window which fails at any layer is immediately discarded. A simple 2-feature classifier can be adjusted to achieve 100% detection rate with a 40% false positive rate, significantly reducing the number of times the cascade must be evaluated. This classifier can serve as the first layer (or attentional operator) of the cascade, followed by a more complex layer (made up of, say, 10 features) which can perform a stronger evaluation on potentially “harder” negative-windows which survived the first-layer evaluation, and so on. A classifier built on this logic achieves slightly lower detection rates with each successive layer but is compensated by a much faster reduction in false positive rate. This relationship can be evaluated through the detection rate and false positive rate equations below:

Given a trained cascade of K classifiers, the false positive rate F of the cascade is:

$$\mathrm{F}=\prod _{i = 1}^K{f}_i$$

where f_i is the false positive rate of the i^th classifier in the cascade on the examples that it is run on. The detection rate D is:

$$\mathrm{D}=\prod _{i = 1}^K{d}_i$$

where d_i is the detection rate of the i^th detector.

As an example, a detection rate of 0.9 and very low false positive rate of approximately 6 x 10⁻⁶ can be achieved by a 10-layer classifier if each layer can achieve a detection rate of 99% (since 0.9 ≈ 0.99¹⁰) and a false positive rate of about 30% (0.30¹⁰ ≈ 6 x 10⁻⁶).

Figure 4. 8. Schematic depiction of the attentional cascade.

Figure 4. 8. Schematic depiction of the attentional cascade.

Like the weak learners in the boosting stage, cascading classifiers are trained with hundreds of positive and negative samples of the object of interest. In the case of the face, the object of interest in this thesis, training was unnecessary as the MATLAB Cascade Object Detector is pre-trained on a large database of faces.

The process of optimal training of the classifiers in the cascade is not straightforward and requires the solution of an intractable optimization problem. Viola and Jones proposed, instead, a simple framework which nevertheless yields a highly efficient and effective classifier. Under this framework, the user selects the minimum acceptable detection rate and maximum acceptable false positive rate for each classifier. Each layer of the cascade is then trained by AdaBoost with the number of features increased until the target detection and false positive rates are met for that layer. A new layer is added to the cascade until the overall target false positive rate is achieved.

4.3.2. ROI Selection

A face detected by the Viola-Jones algorithm is bounded by a rectangular bounding box (bbox) which contains skin pixels (of the face) and non-skin pixels (in the space surrounding the face). The simplest method of ROI selection sets the modest objective of eliminating non-face pixels by cropping the image to a central fraction of the width and height of the bbox. One recommendation is 60% of the width and 80% of the height (Rahman et al. 2016).

Another method involves the use of a rectangular ROI centred on the forehead. The choice of the forehead as a ROI is based on the discovery by Verkruysse et al. that it provides the strongest rPPG signal (2008). It is also supported by the work of Lewandowska et al. which used thermal images from several subjects to establish the forehead as a particularly uniform region for rPPG signal extraction (2011). Both ROI selection methods were used in this thesis. The positional specifications of the ROIs implemented using the second method are based on the recommendation of Tofighi et al. (2014).

Figure 4. 9. Various ROIs for rPPG signal extraction.

Figure 4. 9. Various ROIs for rPPG signal extraction.

Smarter ROI selection methods such as those based on skin detection algorithms are also in use (Kumar, Veeraraghavan, and Sabharwal 2015). These algorithms attempt to take advantage of the higher signal strength associated with larger ROIs, while minimizing the noise contributed by non-skin pixels and the movement on non-rigid facial structures.

4.3.2. Motion Tracking

In the generation of the rPPG signal, it is imperative that the reflectance information be extracted from the same ROI. Inability to maintain this constancy leads to signal noise due to inhomogeneous spatial distribution of skin reflectance. The ideal way of attaining this constancy is by zero relative motion between the subject and camera. However, this is obviously unrealistic, even when attempted, as small variations in the location of the face associated with processes such as breathing and the ejection of blood from the heart (as well as small movements of the camera) invariably lead to error-inducing movements of the ROI. A natural technique for addressing this motional error is to track the movement of the face across frames so that relative motion between the ROI and face is minimized.

Motion tracking could easily be implemented by detecting the face in all frames. However, this is naïve not only because it is computationally inefficient but also because it is prone to significant errors when the face position in a frame is not properly detected.

More robust tracking systems may be designed by calculating the optical flow of image pixels from frame to frame. Optical flow is defined as the pattern of apparent motion of image objects between two consecutive frames caused by the movement of object or camera (Miura et al. 2017). Optical flow is best represented by a 2D displacement vector field in which every vector points from a point in a given frame to an equivalent point in the next frame. Optical flow algorithms may be classed into two categories – dense methods which compute optical flow by considering every pixel and sparse methods which make use of only select pixels. While providing greater tracking accuracy, dense methods are computationally expensive and therefore not viable for real-time applications. A very popular sparse method is the Lucas-Kanade algorithm which limits the optical flow computation to a small subset of image points known as feature points (Lucas and Kanade 1981).

The Lucas-Kanade algorithm is based on three assumptions:

• brightness constancy – the brightness of each pixel is constant between two consecutive frames.
• temporal persistence – movements of image objects are small.
• spatial coherence – neighbouring pixels belonging to a certain surface move in a similar way.

These assumptions are used to construct a system of linear equations describing the optical flow.

Deviation from any one of the three assumptions naturally affects the performance of the algorithm. Large motions, for example, cannot be captured by the optical flow. This limitation has been addressed by Bouguet in a pyramidal implementation of the algorithm (2001). With the pyramids, small motions are eliminated, and large motions are reduced to small motions. The pyramid is constructed by recursive scaling of the resolution of consecutive image frames. The number of frames over which this scaling is done is typically no more than four. If the resolution of the first image is 640 x 480, for example, then that of the next four images will be 320 x 240, 160 x 120, 80 x 60 and 40 x 30.

Figure 4. 10. Schematic depiction of pyramidal implementation of Lucas-Kanade algorithm. Tarasenko et al. (2016). Reproduced with permission from Research India Publications.

Figure 4. 10. Schematic depiction of pyramidal implementation of Lucas-Kanade algorithm. Tarasenko et al. (2016). Reproduced with permission from Research India Publications.

Even with the motion correction by Bouguet’s pyramid, the Lucas-Kanade algorithm is still far from ideal. This means that the optical flow motion estimation will not always match the actual motion of the subject. This manifests as drift in the ROI position which introduces noise in the associated rPPG signals. Garcia attempted to reduce this error by detected the face every second and readjusting the ROI (2013). However, he also reported that this technique may introduce spurious frequency information in the signal and compromise the accuracy of HR estimation.

Figure 4. 11. Tracking points (white) of the Lucas-Kanade algorithm.

Figure 4. 11. Tracking points (white) of the Lucas-Kanade algorithm.

4.4. Signal Extraction

Figure 4. 12. Raw red, green, blue and luminance (Y component of YCbCr) rPPG signals.

Figure 4. 12. Raw red, green, blue and luminance (Y component of YCbCr) rPPG signals.

4.5. Signal Pre-Processing

4.5.1. Signal Detrending

Detrending is used to remove undesirable trends in a signal that may affect its stationarity. A signal is described as stationary if its long-term statistics (such as mean and variance) are relatively invariant in time. Signals typically develop trends through drifting and noising of different kinds. In the case of rPPG signals, this may be caused by changes in physical parameters such as ambient light intensity and temperature. Spectral analysis algorithms typically assume that input signals are at least weakly stationary. The stationarity of an rPPG signal may be significantly compromised in the presence of motional noise. In this research, detrending is performed using the highly recommended smoothness priors approach (SPA), with smoothing parameter λ = 10 and cut-off frequency f_c = 0.059 Hz (Tarvainen, Ranta-aho, and Karjalainen 2002). SPA works by removing the low frequency aperiodic trend component from an input signal, leaving behind a nearly stationary signal of interest.

Figure 4. 13. Detrended red, green, blue and luminance rPPG signals.

Figure 4. 13. Detrended red, green, blue and luminance rPPG signals.

4.5.2. Signal Filtering

Signal filtering was performed with a fourth-order Butterworth bandpass filter with cut-off frequencies of 0.8 and 6 Hz, following the method of Tarassenko et al. (2014). Signal filtering amounts to a denoising process as it eliminates all frequency components outside the bandpass interval.

Figure 4. 14. Filtered red, green, blue and luminance rPPG signals.

Figure 4. 14. Filtered red, green, blue and luminance rPPG signals.

4.5.3. Signal Normalization

The last stage of pre-processing normalizes all signals of interest to zero mean and unit variance according to the method mentioned by Cochran (1988). The normalized signal X_i given by:

$$X_i(t) = \frac{Y_i(t) - {\mu }_i}{{\delta }_i}$$

for each i = R, G and B signals where μ_𝑖 is the mean and 𝛿_𝑖 is the standard deviation of signal 𝑌_𝑖.

Apart from establishing signal comparability, normalization improves algorithmic stability and is a necessary preparatory step for signal processing by procedures such as PCA and correlation analysis.

Figure 4. 15. Normalized red, green, blue and luminance rPPG signals.

Figure 4. 15. Normalized red, green, blue and luminance rPPG signals.

4.6. Spectral Analysis

4.6.1. Fast Fourier Transform (FFT)

The Fourier Transform is a popular method for spectral analysis i.e., the conversion of a signal from the time domain to a representation in the frequency domain. The Fourier transform may be continuous or discrete. Since the rPPG represents a discrete signal, Fourier analysis of the rPPG signal employs the discrete Fourier Transform (DFT).

FFT refers to a class of algorithmic implementations of the DFT. While a naïve (direct from definition) implementation of the DFT has time complexity O(n²), FFT implementation has complexity O(nlogn) and is the most efficient known method for implementing the DFT (Cochran et al. 1967). The DFT is defined by the formula:

$$X_k=\sum _{n = 0}^{N-1}{x}_n{e}^{-i2\pi kn/N}$$

where X_k represents a set of complex numbers of cardinality k, x_n represents the input signal of length N, and k = 0, …, N – 1.

In this thesis, the FFT was implemented using the directly available MATLAB function.

Figure 4. 16. FFT of red, green, blue and luminance rPPG signals.

Figure 4. 16. FFT of red, green, blue and luminance rPPG signals.

The FFT operates on the assumption that the input signal is periodic and is therefore sensitive to discontinuities at signal end points. In the frequency domain, discontinuities translate to the problem of spectral leakage in which the power of the actual frequency is leaked to surrounding frequencies. Leakage can be reduced by minimizing discontinuities (while preserving frequency information). This can be achieved by windowing in which a function is applied to the signal to produce a smoothed redistribution of amplitude with end points equal to zero.

The window function employed in this thesis is the Hann window (depicted in Figure 4.18). This window combines acceptable frequency resolution with moderate spectral leakage. Alternative window functions such as the Blackmann and Blackmann-Harris, though better at reducing spectral leakage, do so at the cost of frequency resolution. The Hann window is highly recommended for general-purpose applications where the nature of the signal is unknown. The mathematical definition of the Hann window (of length N) is given by:

$$w(n) = \frac{1}{2}\left[1 -cos\left(\frac{2\pi n}{N - 1}\right) \right]$$

where n = 0, …, N – 1.

Figure 4. 17. Hann window function (left) and its Fourier transform (right). Olli Niemitalo / Wikipedia Commons / Public Domain.

Figure 4. 17. Hann window function (left) and its Fourier transform (right). Olli Niemitalo / Wikipedia Commons / Public Domain.

4.6.2. Welch Periodogram

A second spectral analysis technique employed in this thesis is the Welch periodogram. It is an implementation of the Fourier transform optimized for unevenly sampled signals and periodic signals of different shapes. The Welch periodogram is an improvement over the standard periodogram as well as Bartlett’s periodogram. Bartlett’s periodogram – also known as the method of averaged periodograms – operates by calculating the average of a series of periodograms derived from non-overlapping segments of an input signal. Its advantage over the standard periodogram is that it reduces signal noise at the cost of frequency resolution. Welch’s method modifies Bartlett’s by averaging over a set of periodograms obtained from overlapping segments of in input signal, yielding greater noise reduction. The procedure for the Welch periodogram involves:

• splitting of signal into overlapping segments
• windowing of overlapped segments in the time domain
• computing of periodograms for individual segments by FFT, followed by squaring of the magnitude
• averaging of individual periodograms (resulting in the reduction of the noise/variance in the power of the frequency components)

Figure 4. 18. Welch periodogram of red, green, blue and luminance rPPG signals.

Figure 4. 18. Welch periodogram of red, green, blue and luminance rPPG signals.

4.7. Heart Rate Determination

4.8. Color Models and Spaces

5.1. Self-Captured Videos

Table 5. 1. Mean error and RMSE of HR estimates (max peak, ambient lighting).

Table 5. 2. Mean error and RMSE of HR estimates (signal peak, ambient lighting).

Table 5. 3. Mean error and RMSE of HR estimates (max peak, dedicated lighting).

Table 5. 4. Mean error and RMSE of HR estimates (signal peak, dedicated lighting).

Figure 5. 1. Bland-Altman plot for HR determined from rPPG signal 4.2 using max peak estimates (ambient lighting).

Figure 5. 1. Bland-Altman plot for HR determined from rPPG signal 4.2 using max peak estimates (ambient lighting).

Figure 5. 2. Bland-Altman plot for HR determined from rPPG signal 4.2 using signal peak estimates (ambient lighting).

Figure 5. 2. Bland-Altman plot for HR determined from rPPG signal 4.2 using signal peak estimates (ambient lighting).

Figure 5. 3. Bland-Altman plot for HR determined from rPPG signal 4.2 using max peak estimates (dedicated lighting).

Figure 5. 3. Bland-Altman plot for HR determined from rPPG signal 4.2 using max peak estimates (dedicated lighting).

Figure 5. 4. Bland-Altman plot for HR determined from rPPG signal 4.2 using signal peak estimates (dedicated lighting).

Figure 5. 4. Bland-Altman plot for HR determined from rPPG signal 4.2 using signal peak estimates (dedicated lighting).

Table 5. 5. Mean SNR of rPPG signals (ambient lighting).

Table 5. 6. Mean SNR of rPPG signals (dedicated lighting).

Table 5. 7. Correlation coefficients (r) and p-values for est. HR vs ref. HR (max peak, ambient lighting).

Table 5. 8. Correlation coefficients (r) and p-values for est. HR vs ref. HR (signal peak, ambient lighting).

Table 5. 9. Correlation coefficients (r) and p-values for est. HR vs ref. HR (max peak, dedicated lighting).

Table 5. 10. Correlation coefficients (r) and p values for est. HR vs ref. HR (signal peak, dedicated lighting).

Figure 5. 5. Scatter plot of est. HR versus ref. HR for rPPG signal 4.2 (max peak, ambient lighting).

Figure 5. 5. Scatter plot of est. HR versus ref. HR for rPPG signal 4.2 (max peak, ambient lighting).

Figure 5. 6. Scatter plot of est. HR versus ref. HR for rPPG signal 4.2 (signal, ambient lighting).

Figure 5. 6. Scatter plot of est. HR versus ref. HR for rPPG signal 4.2 (signal, ambient lighting).

Figure 5. 7. Scatter plot of est. HR versus ref. HR for rPPG signal 4.2 (max peak, dedicated lighting).

Figure 5. 7. Scatter plot of est. HR versus ref. HR for rPPG signal 4.2 (max peak, dedicated lighting).

Figure 5. 8. Scatter plot of est. HR versus ref. HR for rPPG signal 4.2 (signal peak, dedicated lighting).

Figure 5. 8. Scatter plot of est. HR versus ref. HR for rPPG signal 4.2 (signal peak, dedicated lighting).

Table 5. 11. Correlation coefficients (r) and p-values for error versus SNR (ambient lighting).

Table 5. 12. Correlation coefficients (r) and p-values for error versus SNR (dedicated lighting).

Figure 5. 9. Scatter plot of error versus SNR for rPPG signal 4.2 (ambient lighting).

Figure 5. 9. Scatter plot of error versus SNR for rPPG signal 4.2 (ambient lighting).

Figure 5. 10. Scatter plot of error versus SNR for rPPG signal 4.2 (dedicated lighting).

Figure 5. 10. Scatter plot of error versus SNR for rPPG signal 4.2 (dedicated lighting).

5.2. COHFACE Database Videos

Table 5. 13. Mean error and RMSE of HR estimates (ambient lighting).

Table 5. 14. Mean error and RMSE of HR estimates (dedicated lighting).

Figure 5. 11. Bland-Altman plot for HR determined from sig 2.1′ (ambient lighting).

Figure 5. 11. Bland-Altman plot for HR determined from sig 2.1′ (ambient lighting).

Figure 5. 12. Bland-Altman plot for HR determined from sig 3.1′ (ambient lighting).

Figure 5. 12. Bland-Altman plot for HR determined from sig 3.1′ (ambient lighting).

Figure 5. 13. Bland-Altman plot for HR determined from sig 2.1′ (dedicated lighting).

Figure 5. 13. Bland-Altman plot for HR determined from sig 2.1′ (dedicated lighting).

Figure 5. 14. Bland-Altman plot for HR determined from sig 3.1′ (dedicated lighting).

Figure 5. 14. Bland-Altman plot for HR determined from sig 3.1′ (dedicated lighting).

Table 5. 15. Mean SNR for rPPG signals (ambient lighting).

Table 5. 16. Mean SNR for rPPG signals (dedicated lighting).

Table 5. 17. Correlation coefficients (r) and p-values for est. HR vs ref. HR (ambient lighting).

Table 5. 18. Correlation coefficients (r) and p-values for est. HR vs ref. HR (dedicated lighting).

Table 5. 19. Correlation coefficients (r) and p-values for est. HR vs ref. HR (all videos).

Figure 5. 15. Scatter plot of est. HR versus ref. HR for sig 2.1′ (ambient lighting).

Figure 5. 15. Scatter plot of est. HR versus ref. HR for sig 2.1′ (ambient lighting).

Figure 5. 16. Scatter plot of est. HR versus ref. HR for sig 3.1′ (ambient lighting).

Figure 5. 16. Scatter plot of est. HR versus ref. HR for sig 3.1′ (ambient lighting).

Figure 5. 17. Scatter plot of est. HR versus ref. HR for sig 2.1′ (dedicated lighting).

Figure 5. 17. Scatter plot of est. HR versus ref. HR for sig 2.1′ (dedicated lighting).

Figure 5. 18. Scatter plot of est. HR versus ref. HR for sig 2.1′ (dedicated lighting).

Figure 5. 18. Scatter plot of est. HR versus ref. HR for sig 2.1′ (dedicated lighting).

6.1. Implement an Algorithmic Framework for HR Estimation by Video-Based Reflectance Photoplethysmography

Figure 6. 1. Representative example of signal peak in relation to maximum peak.

Figure 6. 1. Representative example of signal peak in relation to maximum peak.

6.2. Investigate the Effect of Presence/Absence of Dedicated Lighting on the SNR of rPPG Signals

6.3. Investigate the Effect of Choice/Size of Region of Interest (ROI) on the SNR of rPPG Signals

6.4. Compare the SNR of Green rPPG Signals with That of rPPG Signals Based on Luminance

6.5. Investigate the Viability of the Summary Autocorrelation Function (SACF) for Improving the SNR of rPPG Signals

6.6. Implement the HR Estimation Algorithm into a GUI

Figure 6. 2. GUI2.

Figure 6. 2. GUI2.

6.7. Prepare a Laboratory Practical Centred on HR Estimation by rPPG

7.1. Conclusions

7.2. Further Work

Appendix A. Raw Results

Table A1. Self-captured videos – rPPG HR estimates (max peak, ambient lighting).

Table A2. Self-captured videos – rPPG HR estimates (sig. peak, ambient lighting).

Table A3. Self-captured videos – rPPG HR estimates (max peak, dedicated lighting).

Table A4. Self-captured videos – rPPG HR estimates (sig. peak, dedicated lighting).

Table A5. COHFACE videos – rPPG HR estimates (ambient lighting).

Table A6. COHFACE videos – rPPG HR estimates (dedicated lighting).

Appendix B. HR Estimation Algorithm

Appendix C. Laboratory Practical

Figure 1. Characteristic PPG signal from pulse oximeter.

Figure 1. Characteristic PPG signal from pulse oximeter.

Figure 2. Transmission (left) and reflection (right) mode PPG.

Figure 2. Transmission (left) and reflection (right) mode PPG.

Figure 3. Various ROIs for rPPG signal extraction.

Figure 3. Various ROIs for rPPG signal extraction.

Table 1. ROI description and rPPG signal classification.

1. 

Open the MATLAB IDE by double clicking on the MATLAB shortcut icon in the desktop window

2. 

Wait for the completion of MATLAB initialization and the appearance of ready at the bottom left corner of the MATLAB window

3. 

Seat at approximately 0.5 m from the webcam (use a metre rule to measure the distance). Seat still facing the camera and ensure that field of view of the webcam is centered on your face.

4. 

Insert the tip of your index figure into the rubber case of the sensor of the pulse oximeter module (already connected to the computer). An IR LED illuminates the finger from the top while a photodetector detects transmitted light at the bottom. Place the hand on a level surface (such as table top or arm rest) and keep it steady. See Figure 4 for correct placement of finger into sensor.

5. 

Type GUI1 followed by the ENTER key into the MATLAB command window. A GUI will appear and will serve as the interface for rPPG signal extraction and the outputting of FFT results (which will be analysed for heart rate estimation).

6. 

Initiate video capture (as well as a timer and simultaneous capture of ground truth pulse oximeter signal) by pressing the START pushbutton at the top of the GUI window.

7. 

Terminate video capture after 1 min (by pressing the STOP pushbutton).

8. 

Press the SELECT VIDEO pushbutton at the top of the GUI window. This is a browse button which will allow you to select the video you just recorded and enter it into the GUI space for rPPG signal extraction.

9. 

Press the PROCESS VIDEO pushbutton adjacent to SELECT VIDEO. This will initiate the processing of the video for rPPG signal. This may take more than 5 min (depending on hardware). The immediate output of this processing will be the display of the FFT of rPPG signals 1 – 4 in four labelled panels.

10.  

The Welch periodograms of signals 1 – 4 can be accessed by pressing the WELCH pushbutton at the bottom of the GUI. Return to the FFT results is achieved by pressing the same pushbutton which is automatically reconfigured and re-labelled FFT.

11. 

Determine the heart rate (in bpm) for both the FFT and Welch periodogram results by clicking on the peak of maximum amplitude in the frequency spectrum.

12. 

An FFT and Welch periodogram of the pulse oximeter reference values is obtained by pressing the GROUND TRUTH pushbutton.

13. 

Record all results in a table like the one in the Recording and presenting your data section below (it is designed to accommodate Results from multiple subjects/videos).

Figure 4. Placement of index finger into pulse oximeter sensor.

Figure 4. Placement of index finger into pulse oximeter sensor.

• For each of the eight (8) heart rate estimates (one for each rPPG signal), you need a corresponding error and signal-to-noise ratio (SNR) measure. Error values are obtained by clicking the (static) drop-down list titled ERROR at the top of the GUI window (it allows display but no selection).
• SNR values are obtained by clicking the adjacent SNR drop-down list.
• Record Error and SNR values in a table like the one specified in the Recording and presenting your data section
• Close the GUI by pressing the EXIT pushbutton and repeat steps 3 – 13 of Exercise 1 and 1 – 4 of Exercise 2 with dedicated light source about 1 metre from the face.

1. 

In the MATLAB command window, type GUI2. A new GUI will appear. Using this GUI, heart rate is determined directly without the need for manual evaluation of the FFTs.

2. 

Adopt the sitting posture employed in Exercise 1 and keep the dedicated lighting.

3. 

Press the START button at the bottom of the GUI window. Video capture will be initiated after 5 s and will yield your HR every 15 s.

4. 

Let the video run for 5 min. By then you will have 20 estimated heart rate values and 20 corresponding pulse oximeter ground truth values. These values will appear sequentially in labelled static text fields in the GUI window and will also be automatically outputted to a file.

5. 

Record your estimated HR and reference HR values in a table like the one specified in the Recording and presenting your data section

• To provide students an introduction to the MATLAB Integrated Development Environment (IDE)
• To provide students experience with extracting a clinically relevant biosignals (rPPG) using a video-based method
• To provide students an introduction to signal processing and signal processing algorithms such as the FFT
• To provide students the opportunity to perform simple statistical analyses on data