Introduction

Due to the impact of the COVID-19 pandemic, face masks have been must-have accessories for the past year. Their effectiveness against the transmission of the virus has been a crucial research topic [1,2,3,4]. Although the complete picture of the transmission route of COVID-19 remains to be fully realized, airborne transmission with the spread of respiratory droplets and aerosols has been shown to play a significant role [5,6,7]. Face masks of all kinds, including respirator-type, surgical masks, and even homemade fabric masks, have been proven effective in filtering large droplets (defined by the WHO as particles larger than 5 μm in diameter) and aerosols going through the face mask [8, 9]. Various face mask materials have been tested against different filtration mechanisms involving the process of diffusion, interception, impaction, and electrostatic attraction [10,11,12]. Besides the materials, the face mask “fit” also affects the effectiveness of the face masks but has been less studied. Popular surgical or homemade face masks are often not tailored to the individual’s face, resulting in a loose or improper fit and allowing droplets and aerosols to leak out [13,14,15]. While the larger droplets have been shown to travel a limited distance, the smaller aerosols can be transported further, linger in the air for hours, and increase the risk of transmission to other people [16, 17]. The fitting of the face masks should be more thoroughly studied to learn how and where the leakage sites are present and to ensure higher protection efficacy of masks.

Research on the face mask fit problem has been conducted mainly through experiments using either human or 3D-printed faces as the subjects. Various metrics are used to quantify the leakage, such as if a wearer can smell lingering particles or direct measurement of particle counts inside and outside of a face mask [18,19,20]. Verma et al. [9] visualized the respiratory jets and reported that the mask fit and the mask material have an important impact on the mask’s effectiveness. With fluorescent tracers tracking particle trajectories, Oestenstad et al. [21, 22] concluded that the facial dimensions are significantly correlated to the leakage location. Tang et al. [23] studied the jet generated from coughing and suggested that a surgical mask, while effectively blocking the forward momentum of the jet, could allow air leakage through the top and side of the mask due to a loose fit. Koh et al. [24] showed that poorly fit N95 respirators can provide less adequate protection than a simple surgical mask, and Brooks et al. [25] from the Center for Diseases Control and Prevention (CDC) suggested that knotting the mask band and tucking in to achieve better fit can largely improve the protection efficiency. From these studies, we can conclude that mask fit is an imperative factor for mask efficiencies. Nevertheless, there has been a lack of knowledge about the relationship between masks perimeters leaks and facial features. The limited sample size in most experimental studies makes it difficult to derive statistically significant correlations, and the filtering efficiency reported often covers a widespread margin [26, 27]. A potential solution to this problem is through the use of 3D morphable face models. Recent research inspected the mask fit for seven participants and concluded that small facial differences could lead to significant fit discrepancies and suggested that 3D models could serve well in creating virtual samples for large-scale tests [28]. Another study by Kolewe et al. [13] used 3D-printed head models for testing the mask protection efficiency and emphasize the benefit of customizability and variability of using 3D models.

To bridge the knowledge gap between the face mask leakage and the facial features, we propose a framework that utilizes 3D-scanned face models to achieve large-scale experiments of fitting masks to faces with a wide variety of facial characteristics and mask designs. Our recent study [29] categorized the faces created from a database into different demographics based on weight, age, gender, and height and study the correlation between these demographic groups and the mask leakage. The results suggested that the CDC-recommended mask design is not an optimal choice for some demographics and the general public should have more options to receive better protection. It was also revealed that the mask material, specifically the mask’s elastic modulus, does not significantly affect the mask fit. In this study, we further utilize a systematic way of morphing the face based on the available populations to isolate the impact of different facial features on mask leakage. The correlation relation found is then tested with a virtual cohort of randomly generated faces. This work aims to build a framework to rapidly and systematically generate a series of test subjects for mask deployment and could be further extended to more complex mask designs. The goal is to enable personalized mask selection, similar to selecting personal clothes, based on some defining key metrics that ensure maximum protection for each member of society.

Materials and methods

Face database

The publicly available Basel Face Model (BFM) [30] is used to generate the face populations that will be used in the next step of modifying certain facial features. The BFM includes high-resolution face scans of more than 100 males and 100 females ranging from 8 to 62 years old, with weights ranging from 40 to 123 kg. The BFM included the principal component analysis (PCA) components of faces [31], which allows us to expand the samples to over 7000 faces with the random weighting of the PCA components. The benefits of creating a virtual cohort of faces are that these faces are created based on the actual face scan, and they share similar meshes, so it is simple to locate the characteristic points such as nose tip across all faces, which can later be used to calculate the characteristic lengths. The utilization of this database allows us to cover a wide range of samples easily and the PCA components can be used to extend the proposed framework to any face.

Geometrically weighted PCA

As previously stated, one of the major problems of face mask fit study is that even a seemingly insignificant facial detail could lead to a very different gap profile. However, it is nearly impossible to look at the effect of these facial features individually with human subjects. This is where the 3D models are unique as they can change only selected regions of the face. 3D morphing face models have long been a popular idea in facial recognition, animation and diverse industries since its introduction by Blanz and Vetter [31]. The most common approach of building a morphable face is through PCA [32], representing the face with combinations of shape vectors [33, 34]. Through different combinations of these vectors, the face can be morphed into different geometries. In the same pioneering paper, Blanz and Vetter further proposed the local model that segments the faces into different regions and employs PCA modification to different divisions to change only the specified region. Many studies later adopt this idea. For example, Paysan et al. [30] manually divided the face into four segments to enhance their facial recognition algorithm. Tena et al. [35] partitioned the face then modified them simultaneously by explicitly enforcing boundary constraints to create natural expressions, and Gilani et al. [36] separated the faces into Voronoi regions based on the geodesic distance for generating population-specific deformable models. This study adopts a novel geometrically weighted PCA to build the morphing faces as the test subjects to isolate and modify individual facial features.

Weighting the data matrix to emphasize particular metrics is common practice when applying modal analysis [37,38,39] and has been adopted to improve facial recognition [40, 41]. Here, we weigh the shape vectors based on the geodesic distance from the specific facial characteristic points to isolate the variation of specific facial features. Figure 1a shows the characteristic points chosen on the mean face of the exemplars based on the principle that each set of characteristic points should represent possible contact points with the mask. The six sets represent the eyes, nose, chin, ear, zygomatic arch, and cheeks region. The geodesic distance of all points on the face to the characteristic points is calculated by solving the Eikonal equation [42]. The geodesic distance from the chin is shown in Fig. 1b. The geometrical weighting is then calculated with the weighting functions weighted by the geodesic distance. The weighting function chosen here is a family of Sigmoidal function which is differentiable and smoothly transitions from 0 to 1:

$$f\left( {x;a_k,c_k} \right) \,=\, \frac{1}{{1 \,+\, e^{ - a_k(x \,-\, c_k)}}}.$$
Fig. 1: Facial features and geometrical weighting schematic.
figure 1

(a) The facial feature points selected for calculating the six facial regions (eyes—red circle; nose—green square; chin—orange triangle; ear—yellow star; zygomatic arch—magenta diamond; cheeks—cyan cross), (b) the geodesic distance from the chin and (c) the corresponding geodesic weighting, and (d) the facial characteristic length (A—face width; B—cheek width; C—nose breadth; D—nose height; E—nose–chin length).

Full size image

The weighting of each vertex on the face mesh is assigned by the relative geodesic distance from the landmark points, with the further point having zero weight and the points closest to the landmark points have unity weight. The width and the center of the weighting function transition area can be controlled with the parameter \(a_k\) and \(c_k\), respectively. These factors are chosen such that regions do not have overlap and the spatial weighting function drops to zero at the border of each region. Figure 1c shows the geodesic weighting based on the geodesic distance from the chin. Below we describe the whole procedure of the proposed “Geometrically weighted PCA”. The geometry of each face is represented by a shape vector \({{{{{{{\boldsymbol{s}}}}}}}}_{{{{{{{\boldsymbol{i}}}}}}}} \,=\, (X_{1,i},Y_{1,i},Z_{1,i},X_{2,i}, \ldots ,Y_{N,i},Z_{N,i}) \,\in\, R^{3N}\) which contains the Cartesian coordinates of its N vertices. For the current study, we adopted over 7000 (the number of faces will be denoted by M) faces randomly generated from the BFM database [30] as the exemplars. The exemplars are arranged into a data matrix A after removing the mean \({{{\bar{\boldsymbol s}}}}\) and introducing the geometric weighting matrix \({{{{{{{\boldsymbol{W}}}}}}}}_{{{{{{{\boldsymbol{i}}}}}}}} \,= (w_{1,i},w_{1,i},w_{1,i},w_{2,i}, \ldots ,w_{N,i},w_{N,i}) \,\in\, R^{3N}\). The covariance matrix C of the data matrix is constructed and its eigenvectors are computed using the singular value decomposition:

$$\begin{array}{*{20}{c}} {{{{\bar{\boldsymbol s}}}} \,=\, \frac{1}{M}\mathop {\sum }\limits_{i \,=\, 1}^M {{{{{{{\boldsymbol{s}}}}}}}}_i,{{{{{{{\boldsymbol{a}}}}}}}}_i \,=\, {{{{{{{\boldsymbol{W}}}}}}}}_i^{{{{{{{\mathbf{T}}}}}}}}\left( {{{{{{{{\boldsymbol{s}}}}}}}}_i - {{{\bar{\boldsymbol s}}}}} \right),{{{{{{{\boldsymbol{A}}}}}}}} \,=\, \left( {{{{{{{{\boldsymbol{a}}}}}}}}_1,{{{{{{{\boldsymbol{a}}}}}}}}_2, \ldots ,{{{{{{{\boldsymbol{a}}}}}}}}_M} \right) \,=\, {{{{{{{\boldsymbol{U}}}}}}}}\Lambda {{{{{{{\mathbf{V}}}}}}}}^{{{{{{{\mathrm{T}}}}}}}},} \\ {{{{{{{{\boldsymbol{C}}}}}}}} \,=\, \frac{1}{M}{{{{{{{\boldsymbol{AA}}}}}}}}^T \,=\, \frac{1}{M}{{{{{{{\boldsymbol{U}}}}}}}}\Lambda ^2{{{{{{{\boldsymbol{U}}}}}}}}^T} \end{array}$$

The columns of the orthogonal matrix U are the eigenvectors of the covariance matrix. The eigenvalues can be acquired with \(\sigma _i^2 \,=\, \frac{{\lambda _i^2}}{M}\), where the \(\lambda _i\)s are diagonal elements of the singular value matrix \({{{{{{{\mathbf{\Lambda }}}}}}}}\) arranged in decreasing order. The eigenvalues stand for the relative magnitude of each principal component. We can then represent the faces with the linear combination of the shape eigenvectors:

$${{{{{{{\boldsymbol{s}}}}}}}} \,=\, {{{\bar{\boldsymbol s}}}} \,+\, \mathop {\sum }\limits_{i \,=\, 1}^P \alpha _i\sigma _i{{{{{{{\boldsymbol{s}}}}}}}}_i,$$

where the number of components used P can be determined with the accumulated eigenvalues. The shape coefficients \(\alpha _i\) can be adjusted to determine the modification extent of a face. Characteristic lengths of the faces defined in Fig. 1d can be used to quantify the effect of changing the PCA weighting α. These characteristic lengths are defined based on the common metrics used for facial recognition but modified, accounting for the contact points of the mask.

Compared with the segmented face approach, the proposed geometrically weighted PCA has some advantages that make it more suitable for our intended use. First, there is no segmentation involved in the process; hence no merging sequence is required. The weighting is derived from the database and can be extended to a larger deformation extent naturally. With the Sigmoidal weighting function, we have significant control over the morphing region, and we can also assign multiple landmark points to achieve length-based morphing. Multiple facial regions can be modified simultaneously as:

$${{{{{{{\boldsymbol{s}}}}}}}} \,=\, {{{\bar{\boldsymbol s}}}} \,+\, \mathop {\sum }\limits_{i \,=\, 1}^{P_1} \alpha _{i,1}\sigma _{i,1}{{{{{{{\boldsymbol{s}}}}}}}}_{i,1} \,+\, \mathop {\sum }\limits_{i \,=\, 1}^{P_2} \alpha _{i,2}\sigma _{i,2}{{{{{{{\boldsymbol{s}}}}}}}}_{i,2} \,+\, \ldots ,$$

In conclusion, this novel method can systematically create morphable faces with natural weighting and is suitable for the next phase of this study: face mask deployment.

Face mask deployment model

Now that we have the virtual faces whose facial features can be modified individually, the next step is to deploy the mask on them to inspect the leakage produced by the change of the facial features. The deployment process starts with placing the mask model on top of the face with elastic bands surrounding the ears with no tension. The resting length of the band then gradually decreases, and the quasi-static equilibrium position is calculated each time until the mask rests in a steady state on the face.

A multi-scale approach that utilizes the minimum elastic energy concept is adopted to account for different mechanics of the mask-face interactions. The total elastic energy of the system is expressed as

$${\it{\epsilon }}_t\left( X \right) \,=\, {\it{\epsilon }}_{{{{{{{{\mathrm{cloth}}}}}}}}}^s \,+\, {\it{\epsilon }}_{{{{{{{{\mathrm{border}}}}}}}}}^s \,+\, {\it{\epsilon }}_{{{{{{{{\mathrm{border}}}}}}}}}^b \,+\, {\it{\epsilon }}_{{{{{{{{\mathrm{band}}}}}}}}}^s,$$

where the four terms in respective order are the extensional elastic energy stored in the cloth, the tension and compression energy of the border strip surrounding the cloth, the bending energy stored in the border strip, and the tension energy stored in the connecting bands. In addition to the energy minimization during the settling phase, the non-penetration contact force between the deformable mask and soft face tissue is also considered. The soft contact is modeled with nonconservative forces as

$$f_{{{{{{{{\mathrm{contact}}}}}}}}} \,=\, \left\{ {\begin{array}{*{20}{c}} {k_{con}\left\| {{{{{{{{\boldsymbol{X}}}}}}}} \,-\, {{{{{{{\boldsymbol{X}}}}}}}}_{{{{{{{{\mathrm{con}}}}}}}}}} \right\|{{{\hat{\boldsymbol n}}}}_{{{{{{{{\mathrm{con}}}}}}}}}} & {{{{{{{{\mathrm{for}}}}}}}}\left( {{{{{{{{\boldsymbol{X}}}}}}}} \,-\, {{{{{{{\boldsymbol{X}}}}}}}}_{{{{{{{{\mathrm{con}}}}}}}}}} \right) \cdot {{{\hat{\boldsymbol n}}}}_{{{{{{{{\mathrm{con}}}}}}}}} \, < \, 0} \\ 0 & {{{{{{{{\mathrm{otherwise}}}}}}}}} \end{array}} \right.,$$

where \({{{{{{{\boldsymbol{X}}}}}}}}_{{{{{{{{\mathrm{con}}}}}}}}}\) is the contact point (or the closest vortex) on the face to the point X on the mask, and \({{{\hat{\boldsymbol n}}}}_{{{{{{{{\mathrm{con}}}}}}}}}\) is the outward-facing normal at this point. The contact stiffness \(k_{{{{{{{{\mathrm{con}}}}}}}}}\) can be specified to take different materials into account. The equilibrium position at each deployment phase is calculated by relating the internal forces of the mask to the derivatives of the energy density function with a set of nonlinear equations for solving the vortices placements. The equilibrium equations are solved iteratively until the final position is reached. An example sequence of deploying a mask to a face is shown in Fig. 2a. Masks can be configured with different sizes and tuck-in ratios, defined as the ratio of the initial and final height of the mask deployment sequence, which is equivalent to the pleating ratio \(L/L_0\) of the side of the mask, shown in Fig. 2b. Readers are referred to our prior study [43] for mathematical details behind each energy term and the parameters adopted in the numerical simulation. This mask deployment procedure is highly robust at handling a wide range of mask materials, geometries, and placements, as well as different faces. Next, we will take a look at the proposed metric for quantifying mask leakage.

Fig. 2: Mask deployment process and leakage definition.
figure 2

(a) The mask deployment sequence from the simulation, (b) definition of the tuck-in ratio, (c) Selective calculated two-dimensional channels depicted with solid red lines, and (d) the schematic of such channel.

Full size image

Hydraulic resistance

For the mask fit, the gap between the mask and face is an obvious metric. However, in our prior research [43], neither the total leakage area nor maximum gap could show the full picture of the leakage. Hence, hydraulic resistance was proposed. The principal idea is to look at the gap between the mask and face as a combination of airflow channels. One end of the channel is at the mouth or nose, where there is a high-pressure due to breathing, coughing, or sneezing, and the other end is at the opening along the mask perimeter. In the current study, we consider the outward protection for talking or coughing respiratory events, and therefore, the channels are started at the mouth. The outflow from the mouth results in several magnitudes higher particle emission rate than breathing from the nose [44]. Figure 2c shows several of these calculated channels following the shortest distance from the high-pressure point under the mask. The Hagen–Poiseuille flow profile is used to approximate the velocity profile for permeable 2D channel flow [45] as

$$v\left( {x,y} \right) \,=\, \frac{{H^2}}{{2\mu }}\frac{{dp}}{{dx}}\left( {\frac{y}{H} \,-\, \frac{{y^2}}{{H^2}}} \right),$$

where x is the streamwise direction, y the cross-flow direction, μ the dynamic viscosity, H represents the channel height, and \(dp/dx\) is the pressure drop along the channel. Note that the velocity profile has been shown to not change significantly even with one side of the channel being porous [46], i.e., the mask. Figure 2d shows an example of such channels. Integrating the velocity profile across the height of the channel, we have the mass flow rate as

$$\dot m \,=\, \frac{{H^3}}{{12\mu }}\frac{{{{\Delta }}P}}{D},$$

where D is the channel length derived with the geodesic distance from the high-pressure point to the mask border. The hydraulic resistance R represents the relative resistance to the airflow that leaks out through the gap versus the cross-flow filtered out through the mask and can be calculated with

$$R \,\equiv\, \frac{{{{\Delta }}P}}{{\dot m}} \,=\, 12\mu \frac{D}{{H^3}}.$$

Since the geometrical parameters D and H are the only varying components, we will report the parameter \(\bar R \,\equiv\, \frac{D}{{H^3}}\). The larger the \(\bar R\) is, the less leakage mask has from its perimeter edges. The channels are assumed to be in a parallel configuration, and assuming the channel has finite width related to the number of the channel N and mask border perimeter P the total hydraulic resistance can be defined as

$$\bar R_{{{{{{{{\mathrm{total}}}}}}}}} \,=\, \left( {\mathop {\sum }\limits_i \frac{P}{N}\frac{1}{{\overline {R_i} }}} \right)^{ - 1}.$$

The mean hydraulic resistance can be defined as if the mask border is composed of N identical channels in a parallel configuration, hence \(\bar R_{{{{{{{{\mathrm{avg}}}}}}}}} \,\equiv\, N\bar R_{{{{{{{{\mathrm{total}}}}}}}}}\) . To better understand how the leakage relates to the facial features, we partition the mask perimeter into three regions: chin, cheeks, and nose and report the average hydraulic resistance within each region. Hydraulic resistance, besides the local gap size, also takes the airflow’s travel path into account and could reflect a more accurate filtration efficiency. Later, we will see a more distinct correlation between facial features and hydraulic resistance than with the simple gap measurement.

Results

In this section, we will show how the pipeline works to reveal the correlation between the facial features and the mask leakage. Figure 3 shows the faces constructed with the mean face and 200 leading modes with \(\alpha _i \,=\, \{ - 1,1\}\) to demonstrate how the geometrically weighted PCA modifies the faces within specified regions. The changes of PCA modes generally emphasize or dampen the corresponding face features, such as increasing the weighting of the cheek feature makes a face more “chubby” or the chin feature modification makes the chin longer or shorter. Note that for some features, the PCA modes modify multiple facial attributes simultaneously. For example, modifying the nose region also leads to changes in lip length. The faces are created from the PCA modes acquired from the BFM database [30], and with same weighting each feature is modified to different extent. For example, varying α from −1.5 to 1.5 for the chin modes modifies the nose–chin length by 12.4% while using the same weighting range for the nose region, the nose height changes by 34%. We first show the total hydraulic resistance with the medium-sized mask (CDC-recommended size 9 in. by 5.4 in.) and a tuck-in ratio of 0.5 deployed on each face across the PCA weighting range \(\alpha \,=\, [ - 1.5,1.5]\). For clarity, we plotted the inverse of the hydraulic resistance here, where larger values indicate more leakage at the edges. In Fig. 3, we overlay the gap height on the mask border, and we can see that the maximum gaps do not show much difference across the range of α for most of the features. However, distinct trends are observed in the hydraulic resistance with the changes in the eyes, nose, and zygomatic arch. The face and mask combinations show that while neither the maximum gap nor total opening area along the mask border has substantial changes, the opening area close to the nose increases with the weighting α of the eyes and nose features. This increase is amplified by the cube relation of hydraulic resistance and gap height which results in the increasing trend of the leakage with the increasing α of the eyes and nose features. In contrast, increasing the α weighting of the zygomatic arch area slightly decreases the total leakage. The reason is that with the increasing α the zygomatic arch flattens and the mask could fit better. Also, the total hydraulic resistance of the mask with nose-clip wire is shown in Fig. 3. The nose wire essentially eliminates the gap at the nose region, and we can see that the gap reduces significantly for eyes and nose features.

Fig. 3: The geodesic weighting of the six facial feature regions and the mean face modified with 200 modes using PCA weighting α  = −1 (left in each subplot) and α = 1 (right).
figure 3

The gap height is overlayed on top of the mask border. Under each face set is the corresponding total hydraulic resistance (blue circle—total hydraulic resistance; orange square—with nose-clip wire) with respect to the PCA weighting change. The mask size is 9 inches. By 5.4 in. and the tuck-in ratio is 0.5 for all figures.

Full size image

To gain more insight on how the leakage is affected by the facial feature change, we plot the inverse of the sectional hydraulic resistance of the eyes, chin, and zygomatic arch region in Fig. 4. We first discuss the relationship between hydraulic resistance and the changes in the eye region. With a medium mask size and tuck-in ratio of 0.5, the leakage increases with α, and the increase of the leakage in the nose section is more prominent than the other two sections. In contrast, changing the zygomatic arch region induces a larger change in the leakage of the chin section. With this mask configuration modifying the chin area does not show much impact on the hydraulic resistance. We can explain these observations by correlating the facial feature weighting α to the characteristic lengths described in Fig. 1d. The correlations are shown in Table 1, where it should be noted that all correlations provided have an R-squared score greater than 0.98. The correlation model is a simple linear regression model with \(\left\{ {\alpha _i} \right\} \,=\, b_0 \,+\, \mathop {\sum }\limits_{i \,=\, 1}^k b_kX_k\), where \(X_k\)s are the corresponding characteristic feature lengths. We can see that the eyes and the nose feature changes are highly correlated to the nose height, and both cases show a noticeable leakage increase when α increases. This suggests that the nose height determines the leakage around the nose. On the other hand, the zygomatic arch feature change is related to the increase of cheek breath, implying that the fuller cheek generates less leakage at the chin area. Interestingly, although changing the chin also greatly modifies the chin length, no noticeable leakage change is observed. To reach a thorough picture of how and where the leakage is formed, we will look at the impact of the mask dimension, decided by the mask size and tuck-in ratio.

Fig. 4: The sectional average hydraulic resistance with different mask sizes, tuck-in ratios, and modified facial regions.
figure 4

The (a) eyes, (b) chin, and (c) zygomatic arch region is modified, respectively.

Full size image
Table 1 Linear correlation coefficient between the PCA weighting of different facial features and the characteristic lengths.
Full size table

In Fig. 4, we show the sectional average hydraulic resistance for two different tuck-in ratios (0.3 and 0.5) for the large mask (10 in. by 6 in.). One immediate observation is that, with the larger mask, the leakage is overall more considerable, and increasing the tuck-in ratio also generally increases the leakage. The main reason for having a larger leakage is that a larger mask could overhang and thus fit less appropriately around the chin area, which is evident from the increase in the chin leakage. In addition, contrary to the nose leakage-dominant medium mask cases, the larger mask exhibits chin leakage that sometimes exceeds the leakage at the nose section by more than double. Another observation is that for the larger mask and larger tuck-in ratio, the chin feature change is highly correlated with the chin leakage profile, while the zygomatic arch feature change is related to the chin leakage for the smaller tuck-in ratio. This indicates that with a more tightly fit (smaller mask and/or smaller tuck-in ratio), the chin leakage is affected by the cheek breadth, while with a looser fit, the chin’s length controls the leakage of the chin. Also, no matter which size of the mask, modifying the nose height directly affects the nose leakage.

From these observation, we conclude that the leakage location and extent are dependent on a combination of the facial feature and the mask configurations. For the smaller, tighter mask, the nose leakage is dominant and could be modified by the nose height; for larger or looser masks, the gap at the chin area becomes increasingly important and correlated to multiple facial features. However, some of the correlations between the facial features and leakage are not intuitive, such as the zygomatic arch is more correlated with the chin leakage than the nose leakage. Two reasons might be behind such observation. First reason is that the facial region modifications involve multiple and even competing effects on the overall geometry of the face. We observe that for the zygomatic region chosen, increasing the weighting α does not result in a monotonic change of a single feature like what has been observed for the nose region. Since the PCA modes are acquired from real faces, they lead to nonlinear changes within the region. Specifically, increasing the weighting suppressed the arching zygomatic bone but protruding the supra-alar creases, which has a competing effect on the nose gap. The second reason is related to the nonlinear relation between the mask opening caused by the relatively large in-plane stiffness of the cloth. While the maks cloth can bend easily, it is almost inextensible and has a much larger stiffness in stretching directions. Therefore, slight changes in one edge can reorient the mask and its placement on the face and consequently modifies the leakage on the other edges. This nonlinear correlation effect is strongly linked to the shape and size of the face mask.

So far, the effects of each facial feature, as defined by α, on mask effectiveness have been studied independently. However, could this sectional approach work when multiple facial features change at the same time? As stated previously, a slight change in facial features could lead to a large leakage profile difference [28] and although a strong linear correlation is observed in our analysis, nonlinear effects could appear when multiple features change simultaneously. To further test whether nonlinear correlation terms can be incorporated to reach a more accurate prediction model of the total hydraulic resistance based on the six dominant PCA features, we employ a sparse regression model on mask deployment results of 1000 randomly generated faces with a medium mask and tuck-in ratio of 0.5. The 1000 faces were generated by randomly changing all 6 PCA features simultaneously. The LASSO (least absolute shrinkage and selection operator regression technique) [47] is used to reach a sparse relationship between the total hydraulic resistance and polynomials of α terms up to third order including all cross-terms. Bayesian information criterion (BIC) [48] is used to obtain the model’s sparsity-promoting regularization parameter and select the set of relevant terms in the regression model (Fig. 5). The model is trained with a random selection of 80% of the data and is tested with the remaining results. An adequate prediction of the hydraulic resistance can be achieved with 11 linear and nonlinear terms as listed in Fig. 5. The results show that the proposed model can be employed to reach a predictive model of the leakage based on facial features and the choice of the mask. From the results shown here and the previous study [29], we confirm the finding by Oestenstad et al. [21, 22] that facial features are a dominant factor in deciding the leakage site. With a large database of faces of different ethnicity, gender and age groups, we anticipate the predictions of the proposed model to be improved and it could provide a more precise assessment model of the leakage formation. Such a study would be otherwise costly and time-consuming to be done in the traditional experimental setups. Furthermore, the model can be employed to check mask designs and provide fitness metrics for different facial shapes.

Fig. 5: (Left) Standard error from the LASSO regression for the training and test data decreases with the inclusion of more polynomial terms.
figure 5

The inset in the left figure of the BIC shows that ~11 polynomial terms are sufficient for a good BIC. (Right) 11 polynomial terms and corresponding coefficients used for the regression model.

Full size image

This section demonstrates how geometrically weighted PCA can bridge the facial characteristic and the mask leakage. The proposed hydraulic resistance metric successfully reveals the connections between the leakage, facial features, and mask configurations. In the next section, we will talk about more aspects awaiting to be explored to build a complete pipeline for personalized mask selection.

Discussion

Now that we already know that the facial features can be directly related to mask leakage, the next step will be building a framework for rapid mask selection. The concept is based on measuring characteristic lengths of the individual’s face using standard digital face-scanning technology available in most mobile devices. Based on our discussion above, the information is then used to determine which mask size/tuck-in ratio fits the best. However, before making the mask choice as simple as picking a shirt from the rack, further studies are necessary to extend these results to different mask designs and demographical categories such as gender, weight, and age. The results are based on the BFM database [30]; hence the PCA modes are inherently restricted to faces in that database. With a more expanded database such as the BU-EEG multimodal facial action database [49], we can expect the leakage assessment model to be increasingly accurate or include the effects of expressions or mouth movements.

Furthermore, once we know how the leakage is formed, the question of what kind of leakage is “good” leakage arises. In other words, if leakage is inevitable, where should we allow leakage to occur such that the aerosol does not linger in the air for too long? Another question to be answered is related to the deformation of the mask due to breathing or sneezing. These respiratory episodes lead to higher pressure in the mask, which pushes the mask away from the face and forms larger gaps. We are working on answering these with an advanced fluid-structure interaction model that utilizes 3D morphing faces, deformable mask models and immersed boundary flow simulation technique.

To summarize, due to the interaction between the face, mask, and airflow, picking an appropriate mask for everyone is a more complex problem than picking outfits, and a lot more work is now getting done to build a streamlined process to provide better protection.

Conclusions

In this study, we aim to understand how facial feature changes lead to mask leakage. The geometrically weighted PCA algorithm is proposed to modify specific parts of the 3D face model based on real face scans. A minimum-energy based, quasi-static model is utilized to deploy masks with various configurations to the faces. Finally, a new mask leakage metric, hydraulic resistance, that considers the airflow traveling path from the mouth to the escaping perimeter on top of the opening is considered to capture the aerosol leakage better. From the results, we concluded that the leakage location and extent are jointly decided by the mask configurations and the facial features. We are now working on answering various questions regarding mask leakage, and we are excited to be able to bridge the knowledge gap for providing better protection to the public.