A wavelet-based decomposition method for a robust extraction of pulse rate from video recordings

Related works

The feasibility of measuring PR by the means of digital camera was first shown on the time-lapse grayscale images (Takano & Ohta, 2007). In search of the most robust algorithm for the extraction of the pulse signal from video recordings, numerous methods have been proposed afterward: color-space-based, blind-source-separation-based (BSS-based), model-based and data-based. Color-space-based extracts the pulse signal from different channels of standard color spaces (Tsouri & Li, 2015; Verkruysse, Svaasand & Nelson, 2008). BSS-based methods can be divided into Independent-Component-Analysis-driven (ICA) (Poh, McDuff & Picard, 2010) and Principal-Component-Analysis-driven (PCA) (Lewandowska et al., 2011) approaches. Model-based methods include chrominance-based (CHROM) (De Haan & Jeanne, 2013), Blood-Volume-Pulse-vector-based (PBV) (De Haan & Van Leest, 2014) and Plane-Orthogonal-to-Skin (POS) (Wang et al., 2017a) algorithms. The last group, data-based methods, includes the Spatial-Subspace-Rotation (2SR) approach (Wang, Stuijk & De Haan, 2015a). Interested readers may refer to thorough reviews of rPPG methods (McDuff et al., 2015; Rouast et al., 2017; Sikdar, Behera & Dogra, 2016), to their algorithmic principles (Wang et al., 2017a) or to the comparison of different steps of PR extraction from the video recordings (Unakafov, 2017).

An important issue regarding the feasibility of the implementation of rPPG in real-world applications is its high sensitivity to different sources of noise. Most of the existing algorithms derive pulse signal from three-dimensional RGB signals and therefore only two distortions can be eliminated by linear combination of all three color channels which is usually not sufficient. This can either be solved by using a camera with multiple bands (McDuff, Gontarek & Picard, 2014) or by algorithmically increasing the dimensionality of the given RGB signals (Wang et al., 2017c). The latter is done by decomposing the raw RGB signals into multiple components by the means of Fourier Transform. The distortion signals are then suppressed in each of the components and the resulting signals are combined into a single one-dimensional rPPG signal. This algorithm is known as Sub-band rPPG (SB) (Wang et al., 2017c).

Outline of the paper

In this paper, a novel algorithm called Continuous-Wavelet-Transform-based Sub-Band rPPG (SB-CWT) is proposed. SB-CWT uses a different filter bank for the decomposition of the RGB signals than the SB algorithm, namely Continuous Wavelet Transform (CWT) using the analytic generalized Morse wavelet. In CWT, the target signal is correlated with the shifted (translated) and compressed/stretched (also termed dilated or scaled) versions of the analyzing wavelet. By varying the values of the scale and translation parameters the CWT coefficients are obtained. Multiplying these coefficients with their corresponding scaled and translated wavelets provides the constituent wavelets of the analyzed signal. Wavelet transform has the potential to increase the number of the degrees of freedom of the distortion elimination in comparison to the Short Time Fourier Transform (STFT) that was applied in SB. Additionally, wavelet transform provides more accurate information about the time and the frequency of an analyzed signal than STFT. This is achieved by varying the aspect ratio (i.e., ratio between the length of the time window and width of the frequency band) (Ganesan, Das & Venkataraman, 2004). The proposed algorithm is evaluated on the publicly available dataset MMSE-HR (Tulyakov et al., 2016) and compared with the state-of-the-art SB algorithm (Wang et al., 2017c).

The structure of the article is as follows. In the next section, Materials and Methods, the used MSSE-HR dataset is presented, the state-of-the-art SB algorithm used for evaluation of the proposed algorithm is described, the proposed SB-CWT is presented step-by-step, preprocessing steps of the reference signals are shown and the evaluation metrics are defined. In the Results and Discussion sections, the experimental results are presented and discussed. The conclusions are derived in the last section.

Dataset description

In our research the subset data from the Multimodal Spontaneous Emotion database (Zhang et al., 2016), which includes 2D frontal face video recordings and blood volume pulse (BVP) measurements—MMSE-HR (Tulyakov et al., 2016)—was used. This dataset consists of 102 RGB image sets and corresponding physiological signals of 40 subjects (23 females and 17 males of various ethnic/racial ancestries) performing six different tasks. Table 1 lists all the tasks together with the emotions arousing from them. Image sets (image size: 1,040 × 1,392 pixels) were recorded with RGB 2D color camera (frame rate of 25 fps) of Di3D dynamic imaging system (Dimensional Imaging Ltd, Glasgow, Scotland). BVP signals were measured using Biopac NIBP100D (Biopac System, Inc., Goleta, CA, USA) with the sampling rate of one kHz and collected by Biopac MP150 data acquisition system (Biopac System, Inc., Goleta, CA, USA). All the recorded data were synchronized.

Sub-Band rPPG (SB) method

The proposed SB-CWT method is benchmarked against the SB algorithm (Wang et al., 2017c). This method uses Fourier Transform for transforming the raw RGB signals into the frequency domain. The signals are then projected onto the plane termed as POS (Wang et al., 2017a), the axes of which most likely encapsulate the pulsatile region. The exact projection direction is then localized using the alpha-tuning (De Haan & Jeanne, 2013). Projection to POS together with alpha-tuning serves for the extraction of the rPPG (pulse) signal from the decomposed raw RGB signals. Next, the individual sub-band pulse signals are weighted using the ratio between the pulsatile amplitude and intensity variation amplitude as the weighting function. Finally, Inverse Fourier Transform is used for the transformation of the signal back into the time domain.

Choosing the analyzing wavelet

The general framework of the proposed SB-CWT algorithm is shown in Fig. 1 and the schematic representation of the algorithm is shown in Fig. 2. The framework consists of 13 steps all of which are explained throughout this section. The first question that emerges in any wavelet-based research is which wavelet is the most appropriate for a given application. In the case of using CWT there are several candidate analytic wavelets, such as Morlet, derivative of Gaussian, Shannon, log Gabor, etc. It was shown that there is a wavelet superfamily called generalized Morse wavelets in which all of the aforementioned wavelets can be subsumed (Lilly & Olhede, 2012). Furthermore, it is even possible to objectively define which wavelet within this family is the most appropriate one to be used (Lilly & Olhede, 2012). General Fourier-domain form of the generalized Morse wavelets is: (1) $${\widehat{\text{ψ}}}_{\text{β},\text{γ}}\left(\text{ω}\right)=U\left(\text{ω}\right)a_{\text{β}\text{γ}}{\text{ω}}^{\text{β}}{e}^{-\text{ωγ}}$$ where U(ω) is the unit step function, a_βγ is the normalizing constant, ω denotes frequency, β is the compactness or the order parameter and γ is the symmetry or the family parameter. Increasing β results in the increased number of oscillations in the time domain and consequent narrowing of the wavelet in the frequency domain. Symmetry parameter γ defines the shape of the wavelet in the frequency domain. The Eq. (1) holds for ω > 0 with ${\widehat{\text{ψ}}}_{\text{β},\text{γ}}\left(\text{ω}\right)=0$ elsewhere (due to the unit step function). Generalized Morse wavelets are controlled by two parameters, β and γ. The next step was therefore to adjust their values. Symmetry parameter γ was set to 3, because the family of wavelets with this γ value possesses a few desirable properties, which make it recommended for general use (Lilly & Olhede, 2012). The γ = 3 family is known as Airy family. Interested readers may refer to Lilly & Olhede (2012) for more details. The only parameter left that needed to be defined was β. This parameter can be defined as P²/γ, where P² is the time-bandwidth product. In our example, P² value was defined for each of the four window lengths l (32, 64, 128 and 256 frames) in such a way that the value of the minimum wavelet bandpass frequency corresponds to the minimum frequency within the human heart rate band (40, 240) BPM. The β value was determined accordingly. The relation between a scale and a pseudo-frequency ω_a (this notation is used because the mapping between a wavelet scale to a frequency can be only done in a general sense) is defined as: (2) $${\text{ω}}_a=\frac{1}{a\text{FF}}{f}_{\text{s}\_\text{cam}},$$ where $f_{\text{s}\_\text{cam}}$ is a camera frame rate, a denotes a scale and FF is the Fourier factor. The maximum scale (corresponds to the minimum wavelet bandpass frequency) can be calculated as N/2σ_t, where N is signal length and σ_t is time-domain standard deviation of generalized Morse wavelet. The latter can be approximated by $\sqrt{P^2/2}$. For the definition of σ_t refer to Eq. (21) in Lilly & Olhede (2009). Fourier factor is defined as $2\text{π}/{\left(\text{β}/\text{γ}\right)}^{1/\text{γ}}$. For example, given ω_a = 40 min⁻¹ and $f_{\text{s\_cam}}=25\text{\hspace{0.17em}fps}$, the maximum scale that corresponds to the ω_a equals 6.8224 is reached at P² = 11.

The next required step was to define how finely the scales were to be discretized. This step was important in our application, because discretization indirectly defines the number of sub-bands. In CWT the base scale is usually defined as 2^1/ν, where ν is the number of voices per octave parameter. The name of this parameter refers to the fact that ν intermediate steps (i.e., scales) are needed if the scale is increased by an octave (i.e., if it is doubled). By increasing the value of ν, the scales are discretized more finely, which results in more sub-bands and consequently in more degrees of freedom for noise suppression. In the present research ν was set to 10, which means that the base scale was set to 2^1/10. The reason for choosing ν = 10 (common values of ν are 10, 12, 14, 16 and 32) was that we did not want to increase the computational time. Increasing ν namely requires more computations, which results in longer processing time of the algorithm. The described procedure of the selection of the wavelet parameters is depicted by step 1 in Fig. 1.

Estimation of the pulse rate from video recordings

Viola-Jones face detector (Viola & Jones, 2001) was used to detect faces from the first image in each image set (step 2 in Fig. 1). The width of the box containing the detected face was resized to 60%, while the height was kept unchanged (Poh, McDuff & Picard, 2011) (step 3 in Fig. 1). This was done to decrease the number of non-facial pixels in the ROI. Position of the ROI was kept constant throughout the entire image set. Raw RGB signals were obtained by spatial averaging of the RGB pixel intensity values inside the ROI in each frame (De Haan & Jeanne, 2013; Verkruysse, Svaasand & Nelson, 2008) (step 4 in Fig. 1): (3) $$C_i=\frac{1}{mn}\sum _{j=1}^m\sum _{k=1}^n{I}_{i_{jk}},$$ where C_i denotes i-th spatial RGB mean color channel signal (also termed raw RGB signal), $I_{i_{jk}}$ denotes pixel intensity value I of a pixel located in j-th row and k-th column in i-th color channel, i = 1, 2, 3 and denotes R, G and B channels, respectively. Spatial averaging reduces the camera quantization noise. In the next step, the raw RGB signals were inverted (step 5 in Fig. 1) in order to enable a better comparison with the reference pulse signals derived from the BVP signals. Pulse signals were namely obtained in reflection rPPG mode and therefore the signal reaches its peaks in each end-diastole phase, when the intensity of the diffuse- or body-reflection component is maximal. On the other side, the BVP signals from the MMSE-HR database were acquired using the measurement device based on the Penaz’s principle (Penaz, 1973) or Vascular Unloading Technique. This is an approach where the digital blood flow is held constant by automatic inflation or deflation of the cuff attached to the finger. The applied pressure is directly related to the arterial pressure and is the highest during the ejection phase of the systole. Further processing was not carried out on the full-length raw RGB signals, but on the windowed signals using four different window lengths l = 32, l = 64, l = 128 and l = 256 (step 6 in Fig. 1). The windowed signals were temporally normalized (De Haan & Jeanne, 2013) (step 7 in Fig. 1) in order to remove the dependency of the raw RGB signals on the average skin reflection color. This normalization is defined as: (4) $$C_{n_{wi}}\left(t\right)=\frac{C_{wi}}{\text{μ}\left({\text{C}}_{wi}\right)}-1,$$ where $C_{n_{wi}}$ is the i-th temporally normalized windowed RGB color channel signal, C_wi is the i-th windowed spatial RGB mean color channel signal, μ(·) is the averaging operator, and i = 1, 2, 3 and denotes R, G and B channels, respectively.

Temporally normalized windowed signals were then decomposed using the CWT (step 8 in Fig. 1). Since all $C_{n_{wi}}$signals are discrete, and the generalized Morse wavelet is defined in Fourier domain, the discretized version of the CWT expressed as an inverse Fourier transform was used: (5) $$W_{a_i}\left(b\right)=\frac{1}{N}\sum _{k=0}^{N-1}{\widehat{C}}_{n_{wi}}\left(k\right){\widehat{\text{ψ}}}_a^{*}\left(k\right)e^{i2\text{π}kb/N}$$ where a and b (b = 0, 1, 2 … N−1) are dilation and translation parameters, respectively, ${\widehat{C}}_{n_{wi}}\left(k\right)$ is the discrete Fourier transform of the i-th temporally normalized windowed RGB signal, ${\widehat{\text{ψ}}}_a^{*}\left(k\right)$ is the complex conjugate of the Fourier-domain form of the wavelet ψ_a, N denotes signal length and k denotes index of frequency (k = 0, 1, 2 … N−1). To ensure unit energy at each scale, the following normalization was used: (6) $${\widehat{\text{ψ}}}_a\left(a{\omega }_k\right)=\sqrt{\frac{2\text{π}a}{\text{Δ}t}}\widehat{\text{ψ}}\left(a{\text{ω}}_k\right),$$ where $\text{Δ}t=1/f_{\text{s\_cam}}$ and ${\text{ω}}_k=2\text{π}k/N\text{Δ}t.$ Combining Eq. (5) with Eq. (6) the CWT can be expressed as: (7) $$W_{a_i}\left(b\right)=\frac{1}{N}\sqrt{\frac{2\text{π}a}{\text{Δ}t}}\sum _{k=0}^{N-1}{\widehat{C}}_{n_{wi}}\left(\frac{2\text{π}k}{N\text{Δ}t}\right){\widehat{\text{ψ}}}^{*}\left(a\frac{2\text{π}k}{N\text{Δ}t}\right)e^{i2\text{π}kb/N}.$$

The Eq. (7) was repeated for all scales a and all values b and the final result W_i, where i = 1, 2, 3, was a matrix of wavelet coefficients with the number of rows equal to the total number of scales and the number of columns equal to the signal length N.

In order to obtain the pulse signal all the wavelet coefficients were then projected on POS plane (Wang et al., 2017a) (step 9 in Fig. 1). POS is based on the physiological reasoning and uses a projection plane that is orthogonal to the temporally normalized skin-tone direction. The POS algorithm was chosen because it offers the best general performance in extracting PR from RGB signals (Wang et al., 2017a). The axes of the POS projection plane (P_p) are defined as: (8) $$P_p=\left(\begin{array}{c}\begin{array}{ccc}0& 1& -1\\ -2& 1& 1\end{array}\\ \end{array}\right).$$ Projection of the wavelet coefficients on P_p results in two projected components, S₁ and S₂: (9) $$\begin{array}{l}{S}_1=W_2-W_3\\ S_2=-2W_1+W_2+W_3.\end{array}$$ The components S₁ and S₂ have in-phase pulsatile and anti-phase specular components which can be algorithmically separated (Wang et al., 2017a). This can be done by the means of the so called alpha-tuning (De Haan & Jeanne, 2013) (step 10 in Fig. 1), which can, for the case of the applied POS algorithm, be expressed as (Wang et al., 2017a): (10) $$h=S_1+\alpha S_2,$$ where α = |S₁|/|S₂|, where |·| denotes complex magnitude operator. In the case when the pulsatile variation is stronger than the specular (S₁ and S₂ are in in-phase), the strength of the resulting signal h is enhanced (Wang et al., 2017a). In the case when specular variations prevail (S₁ and S₂ are in anti-phase), the α moves the specular variation strength of one signal closer to the other, which results in the specular variation being cancelled out (due to the addition of two signals that are in anti-phase) (Wang et al., 2017a). In the next step, the obtained sub-band pulse signals were weighted using the global (scale-dependent) energy distribution function (step 11 in Fig. 1) defined as the sum of the original wavelet coefficients (the ones obtained by the CWT) at each scale (Bousefsaf, Maaoui & Pruski, 2013). The i-th weighted sub-band (i.e., at i-th scale) is defined as: (11) $$P_{w_i}=\sum _{j=1}^n\left|h_{ij}\right|·h_{i*},$$ where $P_{w_i}$ denotes i-th weighted sub-band pulse signal, h_ij denotes the CWT coefficient at i-th scale and j-th translation parameter, h_i* denotes all CWT coefficients at i-th scale and n is the total number of translation parameters. After the weighting, the discretized Inverse Continuous Wavelet Transform (ICWT; step 12 in Fig. 1) together with scale-dependent thresholding was performed in order to transform the weighted sub-band pulse signals back into the time domain. Scale-dependent thresholding allows transformation over selected scales only. The inverses of the selected scales (i.e., frequencies) lie within the range of human heart rate band (i.e., 40–240 BPM). The transformation was performed by summing of the scaled weighted sub-band signals $P_{w\text{\_scaled}}$ over the selected scales only: (12) $$\tilde{P}=\sum _{i=i_{\text{min}}}^m{P}_{w{\text{\_scaled}}_i},$$ where i runs from i_min that is equal to the index of the scale that corresponds to the pseudo-frequency 240 BPM to the index of the last scale, which corresponds to the pseudo-frequency 40 BPM. The scaled weighted sub-band signals are defined as: (13) $$P_{w\_\text{scaled}}=2\log 2^{1/\nu }\frac{1}{C_{\psi }}\text{real}\left(P_w\right),$$ where real(·) denotes the operator that takes only the real part of the complex input and C_ψ is the admissibility constant defined as $a_{\text{β,γ}}/\text{Γ}\left(\text{β}/\text{γ}\right)$ (Γ denotes gamma function). This constant ensures that the wavelet of interest has zero mean. Final pulse signal was obtained by overlap adding zero-mean segments of the pulse signal $(\tilde{P})$ using a sliding window with a sliding step of a one time sample (De Haan & Jeanne, 2013; Wang et al., 2017a) (step 13 in Fig. 1).

Average PR estimation from the obtained rPPG pulse signals was based on the average length of interbeat intervals (IBIs) of the entire pulse signal. The IBIs are defined as the intervals between the consecutive systolic peaks. These were detected using the derivative-based method for finding local maxima. The minimum peak-to-peak distance was set to 0.25 s, which corresponds to the upper limit of the human heart rate band (i.e., 240 BPM). This was the only constraint applied in the rPPG signal peak detection algorithm. Extraction of the pulse signal using the state-of-the-art SB method closely followed its original algorithm (Wang et al., 2017a). The only difference is that before the average PR estimation, the pulse signal was inverted due to the same reasoning as explained above.

To summarize, in the proposed SB-CWT the raw RGB signals extracted from the facial images were decomposed applying the CWT. Wavelet transform coefficients were then projected to the POS plane, the axes of which most likely encapsulate the pulsatile region (Wang et al., 2017a). Next, alpha-tuning (De Haan & Jeanne, 2013) was used to localize the exact projection direction. After that the projected coefficients were weighted using the global (scale-dependent) energy distribution function and finally inverted using the ICWT. The comparison of the algorithmic steps of SB-CWT and SB is shown in Table 2. From this table it can be seen that the two methods differ in signal decomposition step, weighting of the sub-band signals and in the transformation of the sub-band signals back into the time domain.

Processing of the reference signals

Pre-processing of the reference signals included the following steps: bandpass filtering, clipping and squaring. Applied bandpass filter was zero phase first-order Butterworth filter with a lower cut-off frequency of 0.67 Hz and upper cut-off frequency of two Hz (these frequencies correspond to 40 and 120 BPM). Filtered signals were then clipped in order to keep only the parts of the signals that are above zero. By subsequent squaring of the signals, noise and diastolic waves were suppressed, while the systolic peaks were emphasized. All the pre-processing steps served for a robust detection of systolic peaks, needed for the estimation of the reference average PR. This was done in the same manner as described above for the rPPG signals, except for the minimum peak-to-peak distance between two consecutive peaks being empirically set to 0.3 s.

Evaluation metrics

In order to evaluate the performance of the proposed algorithm several metrics were used: Signal-to-Noise ratio (SNR), Mean Absolute Error (MAE), Mean Percentage Error (MPE), Root Mean Square Error (RMSE) and coefficient of determination r². SNR was defined as the ratio of the energy spectrum within five bins around the PR frequency (defined as the peak frequency in the spectrum of the reference BVP signal) and the remaining energy contained in the spectrum. Due to the fact that pulse signal is pseudo-periodic, SNR was calculated within a time interval defined by the 256-frames-long sliding window and sliding step of one frame. The chosen SNR metric is a modified version of the metric proposed by De Haan & Jeanne (2013). Mathematical expression of the chosen SNR metric is: (14) $$\text{SNR}=10{\log }_{10}\left(\frac{{\sum }_{f=40}^{240}{\left(w_t\left(f\right)S\left(f\right)\right)}^2}{{\sum }_{f=40}^{240}{\left(1-w_t\left(f\right)S\left(f\right)\right)}^2}\right)$$ where w_t (f) is a binary window, S (f) is the spectrum of the rPPG signal obtained by the Fourier transform and f is the frequency in BPM.

Mean absolute error measures the absolute average errors’ magnitude. Its formula is: (15) $$\text{MAE}=\frac{1}{n}\sum _{i=1}^n\left|{\text{PR}}_{{\text{rPPG}}_i}-{\text{PR}}_{{\text{ref}}_i}\right|,$$ where ${\text{PR}}_{{\text{rPPG}}_i}$ denotes the i-th average PR value obtained from rPPG signals, ${\text{PR}}_{{\text{ref}}_i}$ is the i-th reference average PR value and n denotes the total number of facial recordings.

Root mean square error is defined as the square root of the average of squared differences between the estimated $({\text{PR}}_{{\text{rPPG}}_i})$ and actual values ${\text{(PR}}_{{\text{ref}}_i})$: (16) $$\text{RMSE}=\sqrt{\frac{1}{n}\sum _{i=1}^n{\left({\text{PR}}_{{\text{rPPG}}_i}-{\text{PR}}_{{\text{ref}}_i}\right)}^2}.$$ Mean percentage error is the average of the percentage error by which the estimated values $({\text{PR}}_{{\text{rPPG}}_i})$ differ from the reference values ${\text{(PR}}_{{\text{ref}}_i})$: (17) $$\text{MPE}=\frac{100\%}{n}\sum _{i=1}^n\frac{{\text{PR}}_{{\text{ref}}_i}-{\text{PR}}_{{\text{rPPG}}_i}}{{\text{PR}}_{{\text{ref}}_i}}.$$

Coefficient of determination r² defines the proportion of the variance in the dependent variable x predictable from the independent variable y. In the case of a linear least square regression, r² is equal to the square of the Pearson correlation coefficient r, which is a measure of the linear dependence of two variables: (18) $$r=\frac{{\displaystyle {\sum }_{i=1}^n\left(P{R}_{rPP{G}_i}-\overline{P{R}_{rPPG}}\right)\left(P{R}_{re{f}_i}-\overline{P{R}_{ref}}\right)}}{\sqrt{{\displaystyle {\sum }_{i=1}^n{\left(P{R}_{rPP{G}_i}-\overline{P{R}_{rPPG}}\right)}^2}}\sqrt{{\displaystyle {\sum }_{i=1}^n{\left(P{R}_{re{f}_i}-\overline{PRref}\right)}^2}}},$$ where $\overline{{\text{PR}}_{\text{rPPG}}}$ and $\overline{{\text{PR}}_{\text{ref}}}$ are the mean values of the estimated and reference PR values, respectively.

The agreement between the reference pulse signal and rPPG signals was assessed using the Bland–Altman plots. In the original plot the differences between the results of the two methods are plotted against their averages (Bland & Altman, 1986). We however used the plots in which the differences between the results are plotted against the reference method (Krouwer, 2008).