MOMBAT: Heart Rate Monitoring from Face Video using Pulse Modeling and Bayesian Tracking

Typically, any face videos based HR estimation method consists of the following three stages; preprocessing, HR estimation, and post-processing.

HR monitoring continuously performs HR estimations at different small time intervals and, usually, facial deformations affect a small number of frames. Due to this, some HR estimates in the monitoring can be spurious due to the inevitable facial deformations. Better estimation can be expected when a large number of frames are considered [47]. But it results in the loss of HR variations which is highly useful for medical purposes [5]. Furthermore, it restricts user-friendliness due to high wait time. Typically, the number of frames in a time interval is chosen such that the cardiovascular pulse wave can complete at least two cycles.

HR monitoring is performed by [39] using green channel variations and band-pass filtering. Likewise, methods [47] and [35] perform HR monitoring using matrix completion and convolution neural network (CNN) respectively. Erroneous HR estimates are rectified by [17] to improve the HR monitoring using image registration and global HR. The global HR is estimated from all the video frames and thus, it can be spurious when temporal signals contain noise.

Detecting the discriminatory facial features is an extensively studied research topic. These are referred to as facial landmark points [50]. Usually, they are located around face boundaries, eyes, eyebrows, mouth and nose. Constrained Local Model (CLM) is highly useful for landmark detection. It consists of: i) point distribution model (PDM) that uses rigid and non-rigid shape transformations for modelling the global location of discriminatory points; ii) patch experts which models the behaviour of a landmark by analysing the appearance around its local neighbourhood; and iii) joint optimization which aims to fit PDM and the experts in the best possible way [42]. The unknown shape parameter $\boldsymbol{p}$ is estimated by the joint optimization, which is given by

where $R$ is the regularization term which restricts the introduction of unlikely shapes; $D_{i}$ is the misalignment in the location of $i^{th}$ landmark in the image $I$; and $\mathrm{\boldsymbol{x}_{i}}$ is the $i^{th}$ landmark location in 3-D which is given by

where $\mathrm{\bar{\boldsymbol{x}_{i}}}$ denotes the mean value of $i^{th}$ feature given by PDM; $\phi_{i}$ is the component matrix; and vector $\boldsymbol{q}$ is used to control the non-rigid shape [2]. Remaining parameters scaling $\mathit{\boldsymbol{s}}$, translation $\boldsymbol{t}$ and rotation $\boldsymbol{R}$ controls the rigid shape. In essence, shape parameters are given by $\boldsymbol{p}=\left[\mathit{\boldsymbol{s}},\boldsymbol{t},\boldsymbol{R},\boldsymbol{q}\right]$. The performance of CLM heavily relies on PDM, patch expert and joint optimization. In CLNF [3], patch experts are given by local neural field for modeling spatial relationships between pixels, while non-uniform regularized landmark mean-shift is proposed for joint optimization by taking into account the reliability of patch experts.

Each temporal signal contains a pulse signal along with the noise. In case of multiple temporal signals, the pulse signal is extracted using blind source separation by estimating the individual source components [34]. Amplitudes of pulse signal and noise in the temporal signals depend on the facial structure, user characteristics (like skin color) and environmental settings (like illumination). Hence, z-score normalization [40] is applied to normalize the temporal signals. Moreover, the temporal signal, $F^{i}$ contains noise, $\boldsymbol{\eta}$ and actual pulse signal, $\boldsymbol{X}_{a}$ but modified by the facial structure. That is,

where $n$ and $A$ denote the time instant and matrix incorporating the effects of facial structure respectively. Further, the actual pulse signal is not known and it requires estimation from the temporal signal, that is,

where $\boldsymbol{X}_{e}$ and $B$ denote the estimated actual pulse and the transformation matrix respectively. It can be observed from Equations (3) and (4) that:

such that $T$ = $BA$ and $\hat{\boldsymbol{\eta}}$ = $B{\boldsymbol{\eta}}$.

For accurate HR monitoring, $\boldsymbol{X}_{e}$ should be similar to $\boldsymbol{X}_{a}$. That is, magnitude of $T$ should be 1 and appropriate shape constraints should be imposed on the estimated pulse spectrum. Such shape constraints are imposed using higher order cumulants [38]. The highest order of cumulant is restricted to 4 because higher-order cumulants are easily affected by the tail of the distribution which makes them sensitive to outliers and they are slightly independent in the middle of the distribution containing useful information [24]. It is proved in [32] that constraints on cumulant similarities till fourth order can be fulfilled by defining the objective function as:

where $K\left[\bullet\right]$ denotes the Kurtosis [12] while $\left|\bullet\right|$ and ^∗ represent the absolute value and conjugate operations respectively. This Kurtosis based maximization is solved using [32] to obtain the estimated pulse signal because it quickly provides the global convergence.

Quality of pulse signal can be estimated using the peak signal to noise ratio (PSNR) [54]. Typically, the amplitude of the pulse spectrum obtained after converting the pulse signal into the frequency domain, should contain a peak at the HR frequency and negligible values at other frequencies. Unfortunately, in the pulse spectrum, the noise increases the amplitude at other frequencies. Thus, PSNR can be defined such that the signal can be interpreted as the amplitudes corresponding to HR frequency while noise can be thought of the amplitudes at the remaining frequencies. Mathematically, the quality given by PSNR, $q$ is given by:

where $\boldsymbol{S}_{p}$ denotes the spectrum of the estimated pulse signal; $\bar{sum}$ performs the sum over all the frequencies; $n_{p}$ represent the neighbourhood size; and $maxLoc$ returns the position containing the maximum value (thus, the location of HR frequency is given by $maxLoc(\boldsymbol{S}_{p})$). In equation (7), signal (or numerator) is obtained by adding the amplitude of HR frequency and its few neighbourhoods while noise (or denominator) is obtained by adding the amplitude of the remaining frequencies.

HR monitoring requires the estimation of multiple HR at various time intervals and eventually concatenation of all HR estimates. Hence, just like the existing HR monitoring methods, we divide the face video into multiple overlapping windows [17].

In this section, we analyze each extracted window to estimate the corresponding cardiovascular pulse, HR and quality. Initially, we detect ROIs from the window and mitigate in-plane face movements. Then, we extract the cardiovascular pulse from the ROIs using Eulerian technique followed by FFT based analysis. Subsequently, we estimate the quality of the pulse according to their out-of-plane deformations and utilize it to rectify the pulse using pulse modeling.

The spectrum of noise-free pulse signal should contain maximum amplitude at the HR frequency, but it is violated when the pulse signal contains noise. The modeled pulse signal obtained after applying the proposed pulse modeling technique, contains noise and thus, it can provide spurious HR. Usually, small HR change is observed between subsequent HR estimates. It motivates us to introduce the Bayesian framework which rectifies the spurious HR estimates. The framework consolidates the likelihood information derived from the current window and the prior information related to previously analyzed windows [13].

To leverage the observation that there is a small HR change between subsequent HR estimates, we want to define our prior information such that the large confidence value is provided when the fluctuation between the current HR and previous HR is small while a small confidence value is provided when the fluctuation is large. Furthermore, we need to provide low prior information about the previous HR whenever they are spurious, otherwise the error can be propagated in the subsequent HR estimates. An important characteristic of spurious HR is that its corresponding pulse spectrum contains multiple peaks, thus the spectrum has low PSNR [54]. These conditions are met by defining the prior information for $\left(i+1\right)^{th}$ window using:

where $h_{i}$ and $\gamma_{i}$ denote the HR frequency and PSNR of the modeled pulse spectrum in the previous window (that is, $i^{th}$ window) respectively; $c$ is a predefined constant; $\boldsymbol{\theta}$ is the set of all probable HR frequencies; and $\mathcal{N}$ denote the normal distribution. It can be observed that the normal distribution is used in Equation (26) such that mean and variance are given by $h_{i}$ and $\frac{c}{\gamma_{i}}$ respectively. Thus, small PSNR results in high variance which in turn results in low prior knowledge. Also, it is suggested by [55] that the fluctuations between subsequent HR estimates usually lie within the range of -11bpm to +11bpm. Thus, we set $c$ equal to 4, so that 3 times of the variance covers most of our permissible HR estimates.

Our definition of prior information in Equation (26) prohibits large fluctuation from the previous HR. Hence, if previous HR is spurious with low PSNR value, then the current HR values should be restricted with large range. But one should consider all the plausible HR frequencies when previous HR is spurious. To incorporate this intuition, we add a constant value in all the plausible HR frequencies, which are lying between 0.7 to 4Hz. That is, we modify the prior information using:

where $P_{e}$ denotes the modified prior information; $P_{g}$ is the distribution described in Equation (26); and $P_{u}$ is given by:

In essence, $P_{u}$ is the uniform distribution, defined in the HR frequency ranges of 0.7 to 4 Hz such that any frequency is equally probable with the value of $\hat{c}$. We describe the parameter selection of $\hat{c}$ in Section 4. Furthermore, when the first window is analyzed $P_{g}$ is set to zero for all the possible HR frequency ranges, so that all the values are equally likely and hence no useful prior information is utilized.

The likelihood function is denoted by $P_{l}\left(\boldsymbol{S_{i+1}}|\boldsymbol{\theta}\right)$ where $\boldsymbol{S_{i+1}}$ denotes the spectrum of reconstructed pulse signal $\hat{\boldsymbol{X}_{e}}$ corresponding to the $\left(i+1\right)^{th}$ window. We estimate it using:

The posterior probability, $P_{p}\left(\boldsymbol{\theta}|\boldsymbol{S_{i+1}}\right)$ is evaluated by applying the Bayes rule [13], that is,

where $P\left(\boldsymbol{S_{i+1}}\right)$ is the evidence factor. Equations (27), (29) and (30) can be combined in the following manner:

where $Z_{2}$ is a normalization coefficient given by:

An illustration of the prior information, likelihood function and their corresponding posterior probability is shown in Figure 6. We apply maximum a posteriori estimation for HR frequency estimation which provides the minimum-error-rate classifier based on zero-one loss function [13]. For brevity, the expected loss incurred on selecting a particular frequency, $d$ is given by:

where $R$ is the incurred loss and $L$ represent the loss function given by:

Further, it is obvious that the sum of likelihood function at all the possible values (which lies between 0.7 to 4Hz in our case) will be equal to one, that is,

It can be seen by combining Equations (33), (34) and (35) that:

Hence, expected loss, $R$ is minimized when $\boldsymbol{\theta}$ is set to the value that maximizes the posterior probability $P_{p}$, that is, $\boldsymbol{\theta}$ is set to HR frequency. Hence, we obtain the HR frequency corresponding to $\left(i+1\right)^{th}$ window, $h_{i+1}$ using:

The corresponding HR is given by:

where $round$ operator rounds off the value to the nearest integer. Some examples depicting the usefulness of the proposed HR tracking are shown in Figure 7. The figure depicts the actual HR monitoring along with the predicted HR monitoring when the proposed HR tracking is avoided and utilized. It demonstrates that the HR monitoring can be improved significantly when the proposed HR tracking is used.

The performance of our method, $MOMBAT$ is evaluated on Intel i5-2400 CPU 3.10 GHz. Total 65 face videos have been collected from 65 different subjects (34 males and 31 females), out of which 15 videos are used for parameter selection (or training) and the remaining 50 videos are used for performance evaluation (or testing). The videos are acquired from Logitech webcam C270 camera which is mounted on a laptop and the subjects are free to perform natural facial movements and head pose variations. The resolution of these acquired videos is 640$\times$480 pixels. Furthermore, we avoid any compression mechanism and save the videos in AVI raw format. These are acquired for 1 minute at 30 frames per second. The ground truth is obtained by simultaneously acquiring the actual pulse from the right index fingertip using CMS 50D+ pulse oximeter. The percentage of distribution of ground truth HR estimation from the acquired database is shown in Figure 8.

The performance metrics used in our experiments are based on the predicted predicted HR error, $\left(\bar{P}\left(i,j\right)-\bar{A}\left(i,j\right)\right)$, where $\bar{P}\left(i,j\right)$ and $\bar{A}\left(i,j\right)$ denote the predicted and actual HR estimates respectively for $i^{th}$ subject in $j^{th}$ window. Accurate HR monitoring method requires that the prediction error is close to zero, alternatively, the mean $\mu$ and standard deviation $\sigma$ of the prediction error should be close to zero. Likewise, the percentage of samples with absolute error less than 5 bpm, $err_{5}$ should be close to 100% for correct HR monitoring. Another metric employed for the evaluation is mean average error, $MAE$ of all the subjects which is given by:

where $\left|\bullet\right|$ is the absolute operator; $n_{i}$ represents the number of windows for $i_{th}$ subject; and $z$ is the total number of subjects. Lower value of $MAE$ indicates that the predicted and estimated HR estimates are close to each other. Similarly, we also use total time, $t_{s}$ required for HR monitoring in seconds as a performance metric. Furthermore, we used the Pearson correlation coefficient, $\rho$ to evaluate the similarity between two variables in terms of linear relationship. It lies between -1 to 1 and is given by:

where $cov$ and $\sigma$ are the covariance and standard deviation operator respectively. Better HR estimation requires high similarity between the predicted and actual HR, that is, high $\rho$.

Our proposed method $MOMBAT$ requires proper selection of four parameters, which are: i) window size; ii) overlapping window size; iii) $\alpha$ representing the number of basis in pulse modeling; and iv) $\hat{c}$ which is required for defining prior in HR tracking (refer Equation (28)). Since minimum possible heart beats can be 42bpm and the proper window should contain at least two cycles of the cardiovascular pulse, we selected the window of 4 second video where slightly more than two cycles can be observed. Similarly, overlap between the successive windows is chosen in an application specific manner. We focus on frequently updating the previous HR from the new HR, typically, twice in a second. Hence, we set the overlap between successive windows at 3.5 seconds. The remaining parameters are set by the value, providing minimum $MAE$ on the training set. The $MAE$ for different parameter values are shown in Figure 9. We set $\alpha$ and $\hat{c}$ equal to 50 and 0.4 respectively, where minimum $MAE$ is attained on the training set.

For rigorous performance analysis, we create several other methods from our proposed method, $MOMBAT$ by avoiding or replacing its components. The following methods are considered for the performance analysis: a) $NorSysI$ which is obtained by avoiding image registration, pulse modeling and HR tracking in $MOMBAT$; b) $NorSysR$ which is same as $NorSysI$ except that it uses red light instead of green light for extracting the temporal signals; c) $NorSys$ which is obtained by avoiding pulse modeling and HR tracking in $MOMBAT$; d) $PulseModP$ and $PulseModL$ which avoid HR tracking, but utilize polynomial and Legendre basis [19] respectively, instead of Fourier basis for pulse modeling in $MOMBAT$; e) $PulseMod$ which is given by considering proposed (or Fourier basis based) pulse modeling, but avoiding HR tracking in $MOMBAT$; and f) $BayTrack$ which is given by avoiding pulse modeling, but considering HR tracking from $MOMBAT$. The description of these subversions of $MOMBAT$ is provided in Table 1. Furthermore, we compare our method with the following existing well known methods: [1]; CHROM [11]; POS [52]; [22]; [47]; [39]; [35]; and [17]. Pulse signal is extracted from the Lagrangian temporal signals using Principal Component Analysis in [1]. [22] registers the face and utilizes Eulerian temporal signals. Model based methods are utilized in [11] and [52] where temporal signals are extracted by fusing RGB color channels. Optical and physiological properties of skin reflection are used to perform such a fusion. Methods [1], [22], [11] and [52] provide one HR value. To conform these methods with $MOMBAT$, we extract the window and then analyze each window using these methods for HR monitoring. We are unable to conduct the comparative analysis with [27] for HR monitoring because it requires large window size as described in [47]. In [35], CNN trained on several windows is used for HR monitoring. The training and test sets contain different windows of the same subjects. For more rigorous analysis, we also perform the experimentation with another method $ModCNN$ where [35] is used, except that its training and testing sets do not contain windows of the same subjects.

Our experimental results on our dataset are presented in Table 2. It can be inferred from the table that [1] provides the most spurious HR monitoring because it requires the tracking of facial features, which is easily affected by expressions. Likewise, [22] exhibits lower performance than the other methods except [1] because it averages all the temporal signals for pulse extraction. This is error-prone because large noise in few temporal signals due to facial expressions, can tremendously affect the cardiovascular pulse after averaging. Both [1] and [22] employ highly time consuming feature tracking and BSS. In contrast, [39] extracts the pulse signal using only green channel intensity differences of full face and avoiding computationally expensive BSS step and feature tracking. It enables [39] to perform in the most computationally efficient manner, but such a method performs spuriously because it is easily affected by facial deformation, as shown in [17].

$NorSysR$ and $NorSysI$ are different only in the way that they utilize red and green light respectively for the temporal signal extraction. It can be observed from Table 2 that $NorSysI$ performs better than $NorSysR$, which indicates that green light is more effective in photo-plethysmographic imaging than red light. This observation is also mentioned in [48]. Furthermore, $NorSysI$ performs better than CHROM [11] and POS [52], which utilize optical and physiological properties of skin reflection to consolidate RGB color channels for temporal signal extraction. It is because CHROM and POS do not provide correct HR estimation when pulse signal and noise share similar amplitudes [52]. In addition, POS fails when face videos are illuminated by multiple light sources [52].

$NorSysI$ is the same as $NorSys$ except that $NorSysI$ avoids image registration and Table 2 points out that $NorSys$ performs better than $NorSysI$. It indicates that performance can be increased by utilizing image registration. Average computational time of $NorSysI$ and $NorSys$ are 6.8 second and 6.84 second, respectively, out of which, BSS is the most computationally expensive step requiring 5.13 seconds. Method $NorSys$ also performs better than [47] due to better ROI selection, image registration and proper BSS technique. [47] utilizes matrix completion to mitigate the noise, which increases the computation time significantly. Also, it can be observed from the table that $PulseModP$, $PulseModL$, $PulseMod$ and $BayTrack$ perform better HR monitoring than $NorSys$. $PulseModP$, $PulseModL$ and $PulseMod$ mitigate the problems of out-of-plane movements in $NorSys$ by modeling the pulse signal (refer Figure 5 for example) and in return, they incur an additional average time of 0.31 and 2.48 sec for the modeling and quality estimation respectively. $PulseMod$ performs better than $PulseModL$ and $PulseModP$ which demonstrate that Fourier basis is better suited for pulse signal modeling. Similarly, $BayTrack$ performs better monitoring than $NorSys$ because it rectifies the HR estimates by incorporating the prior knowledge of the HR estimates. Some of its examples are shown in Figure 7. It incurs an additional average time of 0.22 sec than $NorSys$ due to PSNR estimation.

Table 2 indicates that [35] exhibits good HR monitoring when the training and test sets contain different windows of the same subjects. But when the training and testing sets do not contain windows of the same subjects (that is, $ModCNN$) then there is significant performance degradation. It points out that CNN employed by [35], leverages the facial texture and skin color for HR monitoring. This is obviously a wrong way of performing HR monitoring. Likewise, [17] relying on only face reconstruction is incompetent to handle out-of-plane movements and hence, provide spurious HR estimates. But our method, $MOMBAT$ handles most of the spurious cases by utilizing the pulse modeling and Bayesian tracking. Furthermore, it provides the best HR monitoring amongst all the methods.But it incurs an additional average time of 2.94 sec when compared with $NorSys$ due to modeling and tracking. Such small time differences can be neglected to achieve significantly better HR monitoring for the face videos of 54 sec. A few cases where $MOMBAT$ has successfully performed the HR monitoring are shown in Figure 10. Just like other existing methods, $MOMBAT$ may perform spuriously when the face video contains noise that persists for long duration. Some such spurious monitoring cases by $MOMBAT$ are shown in Figure 11.

One major factor that hampers the progress in this realm of HR analysis using face video is the lack of appropriate datasets [10]. It is stated by [10] that several existing publicly available datasets that are extensively used in the literature, are not appropriate for HR estimation using face videos. One such example is MAHNOB-HCI dataset [44] which involves negligible illumination variations induced by the movie and thus, unable to cater complex real-world scenarios. COHFACE dataset is regarded as a more challenging dataset to cater more realistic conditions than MAHNOB-HCI by [21]. Thus, we have conducted our experiments on COHFACE dataset as well. But it lacks significant motion variations and thus, we create and conduct experiments on our dataset for better evaluation of our proposed method, $MOMBAT$.

The COHFACE dataset contains 160 face videos acquired from 40 subjects. Experiments are conducted on this dataset using the performance metrics, parameter selection and methods described in Section 4.2, Section 4.3 and Section 4.4 respectively. The corresponding results are shown in Table 3. It can be observed that these results are similar to the results on our dataset, that is, $MOMBAT$ performs best amongst the considered methods. Furthermore, it can be observed that the efficacy of $MOMBAT$ reduces slightly when COHFACE dataset is considered rather than our dataset. It is because the COHFACE dataset contains compressed videos which deteriorate the HR analysis [36].