1.

Introduction

Photoplethysmography (PPG) is a noninvasive optical technique for cardiovascular monitoring.¹ In its simplest form, a PPG device comprises an illumination source (e.g., LED) and detector (e.g., photodiode). By monitoring the temporal illumination changes, PPG devices measure the pulsatile blood pulse waveform from transient local arterial volume fluctuations. This information can be used to monitor cardiovascular characteristics such as heart rate, heart rate variability, blood pressure, and cardiac output.¹ However, existing contact-based systems are limited to single-location monitoring (e.g., finger), can only be used by one individual per device, and motion artifacts inhibit its use in ambulatory scenarios.

Photoplethysmographic imaging (PPGI) systems are biophotonic systems that have recently gained interest for noncontact widefield cardiovascular monitoring of cardiac parameters such as heart rate, breathing rate, pulse oxygen saturation, and heart rate variability.^2,3 Extending upon PPG theory, PPGI systems decouple the illumination source and detector. A camera is commonly used among PPGI systems as the illumination detector, enabling new types of monitoring such as spatial perfusion analysis^4,5 and multi-individual monitoring.⁶ However, many systems use coarse spatial averaging to estimate cardiovascular perfusion, such as averaging over the facial bounding box^6–8 and hardcoded regions.^9–11 One study attained increased accuracy by incorporating spatial pulsatility priors;¹² however, the method relies on a real-time estimate of the true physiological state based on the aforementioned coarse averaging techniques to achieve accurate prediction.

In this paper, we developed a continuous probabilistic pulsatility model for importance-weighted blood pulse waveform extraction. The continuous model can be used by PPGI systems of any resolution through appropriate sampling to extract robust blood pulse waveforms. The model was developed using a data-driven approach over a 23 participant sample with highly varying characteristics (11/12 female/male, age 11 to 60 years, body fat 10.5% to 42.3%, muscle 31.0% to 52.7%, BMI 16.4 to $35.1\mathrm{kg}·{\mathrm{m}}^{-2}$). Using blood pulse waveform spatial correlation priors, an importance weighting scheme was developed which assigned locations with consistently strong pulsatility a higher weight. Samples were aligned and aggregated in a common Cartesian space, and the continuous probabilistic pulsatility model was computed using a kernel density estimation approach. This method was compared against a whole-area uniform spatial averaging approach used by existing studies.^6–8 Results showed that signals extracted using the pulsatility model were statistically significantly stronger in temporal (correlation, $p<0.01$) and spectral (SNR, $p<0.01$) characteristics than uniform averaging, and heart rates were in tight agreement with ground-truth measurements ($r^2=0.9619$, error $\mu =0.52\mathrm{bpm}$, $\sigma =1.7\mathrm{bpm}$). Model visualization elucidated important arterial pathways, including the neck, malar regions, glabella regions, lips, and nose. This work significantly extends a previous case study¹³ by developing a data-driven multiparticipant model rather than participant-specific weights, and presents more rigorous analysis across a larger participant sample. We discuss how the model can be used and trained for custom applications.

2.

Methods

The goal was to compute a continuous probabilistic spatial pulsatility model for use as a priori information in PPGI systems. By computing a continuous model, it can be used by datasets of any resolution through appropriate discrete spatial sampling. Figure 1 shows the processing pipeline to generate this pulsatility model. This study was approved by a Research Ethics Committee at the University of Waterloo and performed in accordance with the Declaration of Helsinki. Informed consent was obtained from all participants prior to data collection. Additionally, informed consent was obtained from those individuals whose photos were used in this paper.

Fig. 1

Processing pipeline to derive the continuous probabilistic pulsatility model. Videos of 23 participants were recorded using a PPGI system synchronously with the ground-truth waveform (1). Correlation priors were computed by comparing pixelwise temporal signals against the ground-truth signal using Pearson’s linear correlation coefficient ($\rho$) (2). The correlation prior maps were aligned relative to a template and projected into a 2-D Cartesian space, and spatially aggregated (3). A 2-D kernel density estimation method was used in this space (4) to generate the data-driven continuous probabilistic pulsatility model (5).

2.2.

Probabilistic Pulsatility Model

Different spatial locations exhibit varying amounts of observed pulsatility. For example, a skin location directly above a superficial artery will exhibit high pulsatility due to transient changes in local tissue optical properties from the arterial pulse,^1,16 whereas occluding hair will inhibit the observed pulsatility. Let $p(x,y)$ be the probabilistic pulsatility model such that $p(x,y)$ quantifies the probability of observing arterial pulsatility at location $(x,y)$. Using a physiologically motivated data-driven approach, $p(x,y)$ can be computed by determining the locations that consistently exhibited pulsatility across a diverse sample of participants (see Table 1). Such a model can be used as an a priori model for new data when extracting cardiovascular properties.

Table 1

Sample demographics.

Demographic	Sample representation
$n$ (male/female)	23 (12/11)
Age (years)	11 to 60 ($29.6\pm 11.9$)a
Mass (kg)	42.6 to 107.0 ($73.4\pm 17.8$)a
Height (cm)	145 to 193 ($168.5\pm 11.4$)a
Body fat (%)	10.5 to 42.3 ($21.3\pm 7.8$)a
Muscle (%)	31.0 to 52.7 ($39.7\pm 4.5$)a
BMI ($\mathrm{kg}·{\mathrm{m}}^{-2}$)	16.4 to 35.1 ($25.7\pm 5.2$)a

^aμ±σ.

2.2.1.

Absorbance mapping

Let $f(x,y,t)$ be a set of frames. The blood pulse waveform signal is typically represented by transient changes in absorbance rather than reflectance, thus $f$ was transformed to absorbance

Eq. (1)

$$a(x,y,t)=-\log [f(x,y,t)].$$

Each signal was then temporally detrended¹⁷ to normalize the illumination and eliminate respiratory-induced artifacts.

2.2.2.

Correlation priors

A transformation $T$ was sought to map the set of absorbance frames $a(x,y,t)$ to a pulsatility strength map $C(x,y)$ describing pixelwise pulsatile components

Eq. (2)

$$C(x,y)=T[a(x,y,t)].$$

Rather than estimating the pulsatility based on heuristic information, which is participant-independent and may introduce uncontrolled sources of noise, this model can be augmented by incorporating prior information directly into the transformation

Eq. (3)

$$C(x,y)=T[a(x,y,t)|z],$$

where $z$ is the ground-truth blood pulse waveform. This was measured synchronously with the frames, and Pearson’s linear correlation coefficient was computed between the ground-truth signal and each pixel’s temporal signal

Eq. (4)

$$T[a(x,y,t)|z]=\frac{{\sigma }_{az}}{{\sigma }_a{\sigma }_z}·{\mathbb{1}}_{\mathbb{R}+}(\frac{{\sigma }_{az}}{{\sigma }_a{\sigma }_z}),$$

where ${\sigma }_{az}$ is the covariance between the pixel and ground-truth signals, ${\sigma }_a$ and ${\sigma }_z$ are the standard deviations of the pixel and ground-truth signal, respectively, and ${\mathbb{1}}_{{\mathbb{R}}^{+}}$ is the indicator function of positive real numbers. This computation is scale- and offset-independent, suitable for capturing the ratio nature of the blood pulse waveform.

2.2.3.

Importance-weighted kernel density estimation

A physiologically derived importance-weighted scheme was developed to quantify the spatial pulsatility strength. An anatomical location that exhibited strong pulsatility contributed a larger weight than those that contained weak or no pulsing. This system allows for continuous kernel-based probabilistic pulsatility density estimation later. The importance map for participant $i$ was computed as

Eq. (5)

$$V_i(x,y)=\max [C_i(x,y),0].$$

To infer pulsatility patterns across the whole sample, the primary anatomical locations must be coaligned. The camera was systematically and consistently set up for all participants; however, differences in anatomy, minor rotation (relative to the camera), and translation (relative to the frame region) were observed. To correct for these relative distortions, the problem was posed as a coordinate mapping problem, where each participant’s weight data $V_i(x,y)$ were projected into the coaligned pulsatility space $V(x^{\prime },y^{\prime })$. Mathematically

Eq. (6)

$$\left[\begin{array}{c}{x}^{\prime }\\ y^{\prime }\end{array}\right]=H_i\left[\right(\begin{array}{c}x\\ y\end{array}\left)\right],$$

where $H_i$ is a coordinate mapping function that maps $(x,y)$ from participant space $V_i$ to $(x^{\prime },y^{\prime })$ in the coaligned space $V$. Note that this transformation projects points directly into a Cartesian space. There is no need for interpolation, which may have caused local inaccuracies. Implementation details of $H_i$ are discussed later (Sec. 2.3). This aggregate coaligned pulsatility space, $V(x^{\prime },y^{\prime })$, was populated with weighted points from each participant

Eq. (7)

$$V(x^{\prime },y^{\prime })=\text{median}[V_i(x_i,y_i)],$$

where $(x_i,y_i)$ are the coordinates in the participant space that project to $(x^{\prime },y^{\prime })$ in the aggregate space according to the mapping function $H_i$.

A resolution-agnostic model can be computed by estimating a continuous probability density function, and sampling this density function according to the given system’s resolution. A modified Parzen–Rosenblatt kernel density estimation method^18,19 was used to estimate the continuous pulsatility probability density function

Eq. (8)

$$p(x,y)=\frac{1}{|V|}\sum _{k=1}^n\frac{V(x_k,y_k)}{w^2}\mathrm{\Phi }\left(\frac{\mathbf{v}-{\mathbf{v}}_k}{w}\right),$$

where $n$ is the total number of points, $w$ is the spatial window width, $\mathbf{v}=(x,y)$, $\mathrm{\Phi }$ is the window kernel, and $|V|$ is a normalization term such that ${\int }_x{\int }_{y}p(x,y)\mathrm{d}y\mathrm{d}x=0$. The kernel was scaled according to the datum’s pulsatility weight $V(x_k,y_k)$. Using the two-dimensional (2-D) Gaussian kernel,

Eq. (9)

$$\mathrm{\Phi }(\mathbf{u})=\frac{1}{2\pi }\exp (-\frac{{\mathbf{u}}^T\mathbf{u}}{2}),$$

the final probabilistic pulsatility model formulation becomes

Eq. (10)

$$p(x,y)=\frac{1}{|V|}\sum _{k=1}^n\frac{V(x_k,y_k)}{{(w\sqrt{2\pi })}^2}\exp (-\frac{1}{2}\frac{\Vert \mathbf{v}-{\mathbf{v}}_k\Vert }{w^2}),$$

where $w$ is modeled as the spatial standard deviation. This pulsatility model can be used by systems of any resolution through appropriate discrete sampling

Eq. (11)

$$\mathrm{\Omega }(x,y)=\sum _{n=-\infty }^{\infty }\sum _{m=-\infty }^{\infty }\delta (x-n{\tau }_x,y-m{\tau }_y)p(x,y),$$

where $\delta$ is the 2-D Dirac delta function and ${\tau }_x,{\tau }_y$ are the resolution periods in the coordinate space of $p(x,y)$. A physiologically derived blood pulse waveform can be extracted using the discretely sampled pulsatility map

Eq. (12)

$$z(t)=\sum _x\sum _y{a}_i(x,y,t)\mathrm{\Omega }(x,y).$$

Figure 2 shows a graphical depiction of this process.

Fig. 2

Extracting a blood pulse waveform using the probabilistic pulsatility model. The continuous model was discretely sampled to match the resolution of the target frames, and transformed to align to match the anatomical characteristics of the participant. Pixelwise weighted averaging results in a robust blood pulse waveform.

3.

Results

3.2.

Data Analysis

Table 1 provides a summary of the sample demographics measured using bioelectrical impedance analysis. The sample contained a wide range of ages (11 to 60 years), body compositions (BMI 16.4 to $35.1\mathrm{kg}·{\mathrm{m}}^{-2}$), and fair gender representation (11/13 female/male).

Visualizing the final probabilistic pulsatility model elucidated common arterial pathways. Figure 3 shows the (sampled) density across all participants, as well as split by gender. The faces appeared structurally similar to a typical human face, indicating accurate sample alignment based on anatomical anchor points. The forehead, eyebrows, eyes, nose, nostrils, lips, chin, and neck were all visually distinguishable. Characteristic differences existed between the female and male probabilistic faces; however, the primary anatomy was consistent, leading to a cohesive combined pulsatility map. Facial hair in the male data was observed as reduced pulsatility probability in those areas due to occlusion.

Fig. 3

(a) Discretely sampled probabilistic pulsatility model across all participants, as well as (b, c) split by gender. The effect of facial hair on reduced observed pulsatility is apparent around the mouth and chin of the male distribution (c). Both genders exhibited strong pulsatility in anatomically relevant locations, such as the neck (carotid artery) and cheeks (facial arteries).

The strongest observed arterial pathway traveled across both sides of the neck. Its anatomical location was consistent with the common carotid artery, which is the primary arterial pathway to the brain and face. The common carotid artery bifurcates into the external (facial) and internal (cerebral) carotid arteries at the top of the neck. The continued blood flow traveling through the external carotid artery across the jawline could be observed in female pulsatility map [Fig. 3(b)]. This phenomenon was not present in the male data [Fig. 3(c)], due to facial hair occlusion in some participants. Areas around mouth and chin exhibited reduced pulsatility compared to the female data due to facial hair occlusion. Common areas of with high pulsing probability across all participants in both genders included the carotid artery pathways, the lips, malar regions, nose, and glabella regions. Pulsing was additionally observed in females along the jaw line.

Blood pulse waveform signals extracted using the probabilistic pulsatility model in a leave-one-out cross-validation scheme exhibited stronger temporal and spectral characteristics when compared to those extracted with the FaceMean method. Figure 4 shows the distributions of temporal correlation and SNR results across the participant sample. Combined, these measures provide signal fidelity information in both the time and frequency domains. Pairwise comparison of the signals extracted using the proposed and FaceMean methods showed that the proposed method yielded statistically significant correlation ($W=39,p<0.01$) and SNR ($W=31,p<0.01$) compared to FaceMean signals. Figure 5 shows the pairwise differences in correlation and SNR for each of the 23 participants. Temporal and spectral fidelity improvements were observed in a large number of participants. These results indicate that the probabilistic pulsatility model identified areas of strong pulsatility.

Fig. 4

Comparison of signal fidelity using the proposed probabilistic pulsatility model versus spatial averaging. Signals extracted using the pulsatility model exhibited statistically significantly higher pairwise correlation to the ground-truth waveform ($W=39,p<0.01$) and spectral SNR ($W=31,p<0.01$). (**Wilcoxon signed rank test, $p<0.01$)

Fig. 5

Pairwise comparison of correlation and SNR between signals extracted using the proposed pulsatility model (black) and FaceMean (gray). Improvements were observed across most participants when using the probabilistic pulsatility model for extracting blood pulse waveform signals.

Many PPGI studies investigate heart rate extraction algorithms for remote monitoring. To validate the proposed system’s ability to extract heart rate, correlation and agreement to the ground-truth heart rate were investigated. Figure 6 shows the Bland–Altman plot²³ for this heart rate comparison. The proposed system not only attained high correlation to the true heart rate ($r^2=0.9619$), but it also achieved strong and tight agreement (error $\mu =0.52\mathrm{bpm}$, $\sigma =1.7\mathrm{bpm}$). This is consistent with the strong signal fidelity results in Fig. 4, which enable strong inferential properties of summary statistics such as heart rate.

Fig. 6

Heart rate estimation using the pulsatility density function attained strong correlation to the true heart rate ($r^2=0.9619$) and high degree of agreement (error $\mu =0.52\mathrm{bpm}$ and $\sigma =1.69\mathrm{bpm}$). Marker size indicates the number of points at that coordinate.

4.

Discussion

The results demonstrated that using a spatial probabilistic pulsatility model can enhance blood pulse waveform extraction in PPGI data. The data-driven pulsatility model, which was constructed using data from a 23 participant sample with widely varying demographics, provides a priori anatomical guidance based on consistently observed pulsatile regions.

The data-driven model is built using training data containing spatial pulsatility information. Thus, the new “test” data must sufficiently match the training data’s spatial perspective for proper alignment. Though advanced warping models can adjust for some inconsistencies, large deviations in rotation and perspective relative to the training data may produce erroneous alignment. In particular, different imaging and illumination perspectives may change the light–tissue interaction geometry (e.g., shadows), leading to varied photon migration paths. To address this challenge, the model was designed such that its methodology is agnostic to the type of data with which it is trained. Thus, independent training data sets can be used [$f(x,y,t)$ in Eq. (1)] to build custom pulsatility models suitable for the study’s specific test environment. For example, some systems may find that training based on rotated viewpoints or a certain class of demographic (e.g., gender or age) may yield increased results for their specific test environment. Additionally, this model can be trained on anatomical locations other than the head, enabling whole-body cardiovascular monitoring. The model use can be extended to detect abnormal perfusion patterns that may be early markers for disorders such as peripheral vascular disease or arteriosclerosis.

The model was designed as a continuous model so that it can be sampled by systems of any resolution. To be used, all that is required is a coordinate mapping function $H$ [from Eq. (6)] for spatial alignment to the template model. In offline systems, this can be accomplished by manual or semiautomatic spatial alignment through various methods.²⁴ In real-time systems, automatic alignment is required, and can be accomplished using methods such as automatic face fitting.²⁵ In low-motion scenarios (e.g., sleeping studies and controlled experiments), a single alignment operation may be sufficient. In scenarios with increased motion, one strategy could be to calibrate the first frame to the template, and align all subsequent frames to the calibrated source frame, as it would have more similarities to frames within the same video than the model template.

5.

Conclusions

Here, we have developed a continuous probabilistic pulsatility model that describes anatomical locations that consistently exhibited arterial pulsing across a 23 participant sample. This model can be used as a priori information to enhance signal extraction in PPGI systems. Since the model is a continuous model, it can be used by systems of any resolution via appropriate spatial sampling. Results showed that signals extracted using the pulsatility model exhibited statistically significant higher correlation and SNR versus unguided whole-area uniform averaging. Discretely sampled maps identified areas of consistently strong pulsing across the training data in the head, specifically the neck, malar regions, glabella regions, lips, and nose. Discussions demonstrated how the model may be used in custom environments for enhanced blood pulse waveform extraction.