Landmark-based homologous multi-point warping approach to 3D facial recognition using multiple datasets

Dataset & description

We used three datasets to validate the robustness of our method. The first was acquired from the Stirling/ESRC 3D face database, which was captured by a Di3D camera system (Stirling-ESRC, 2018). These images are in the format of wavefront obj files containing 101 subjects with 3D facial scans in a neutral position. The database was intended to facilitate research into face recognition, expression and perception, and we randomly selected 58 subjects for this study. The dataset was used as a test set for a competition involving 3D face reconstruction from 2D images, with the 3D scans acting as the ’ground truth’, at an IEEE conference. The second dataset was the Bosphorus database, which was intended for research on 3D and 2D human face processing tasks and contains 105 subjects. We randomly selected 57 subjects for this study. The dataset was acquired using structured-light-based 3D system, with the subjects being instructed to sit at a 1.5 m distance from the camera; the sensor resolution in the x, y and z directions was 0.3, 0.3, and 0.4 mm, respectively, and high-resolution color texture was used (Savran et al., 2008). The third dataset was the Face Recognition Grand Challenge Version 2.0 (FRGC v2) database, consisting of 466 facial images, of which we randomly selected 120 subjects for this study. Each subject image was captured under uniform illumination, with high resolution and fairly uncontrolled conditions (Phillips et al., 2005).

Template mesh

A template mesh of 92,995 vertices and 183,996 triangles was created by manually locating 16 anatomical points on the 3D face (Fig. 2) called fixed points, according to facial landmark standards (Caple & Stephan, 2016) (for more detail, see Table 1). The fixed landmarks were not subjected to sliding, but were used to establish the warping fields for minimizing the bending energy. Due to its ease of detection and pose correction (Creusot, Pears & Austin, 2010) and its invariance to facial expression (Colombo, Cusano & Schettini, 2006), the nose tip (pronasale) was selected as the most robust and prominent landmark point. The nose tip area can be approximated as a hemisphere on the human face, although any other facial anatomical point could be used. This is where the sliding points begin to spread across the facial surface. Using this fixed point (the pronasale), 484 semi-landmarks were automatically generated, with the overlapping on the pronasale shown in blue. These were first randomly placed on the facial mesh before being uniformly distributed on the selected facial surface (Fig. 3). This was done using the locational positions of the fixed anatomical points with 1.5 mm radius to accommodate all 500 points. To quantify the morphological data for the complex 3D traits of both reference and target shapes, we used geometric morphometric tools based on previously reported landmark-based methodologies (Halazonetis, 2014; Klingenberg & Zaklan, 2000; Kouli et al., 2018; Yong et al., 2018; Zelditch, Swiderski & Sheets, 2012), and the automatic point placement, semi-landmark sliding, and landmark acquisition were implemented in ViewBox 4.0 (Halazonetis, 2014).

Homologous multi-point warping

The geometry of curves and surfaces is easy in 2D or 3D, but it is less easy to define semi-landmarks for non-planar surfaces in 3D (Huanca Ghislanzoni et al., 2017) as they are not guaranteed to be homologous after the first placement. However, this could be achieved by subjecting the semi-landmarks to sliding in the direction that reduces shape variance, thus closely positioning the points at the same locations in the 3D space. The sliding step is important, as it places the landmarks in positions where they correspond better to each other between individuals (Mitteroecker et al., 2013). These semi-landmarks were allowed to slide on the curves and the surface mesh of each target using TPS warping of the template, which positioned the reference points on the target facial mesh by minimizing the bending energy.

According to Bookstein (1989), physical steel takes a bending form with a small displacement. This is because the function (x, y, z) is the configuration of lowest physical bending energy, which is consistent with the given constraints. In this 3D face deformation, the transformation of TPS is done mathematically via the interpolation of a smooth mapping of h from ℝ³ → ℝ³. This is a selected set of corresponding points P_Ri,P_Ti, i = 1, …, N on the faces of the reference object (template) and target (subject) that minimizes the bending energy function E(h)using the following interpolation conditions (Bookstein, 1989; Bookstein, 1997a; Corner, Lele & Richtsmeier, 1992): (1)${\displaystyle \mathrm{E}\left(h\right)={\iiint }_{{\mathbb{R}}^3}\left({\left(\frac{{\partial }^{2}h}{\partial x^2}\right)}^2+{\left(\frac{{\partial }^{2}h}{\partial y^2}\right)}^2+{\left(\frac{{\partial }^{2}h}{\partial z^2}\right)}^2+2{\left(\frac{{\partial }^{2}h}{\partial xy}\right)}^2+2{\left(\frac{{\partial }^{2}h}{\partial xz}\right)}^2+2{\left(\frac{{\partial }^{2}h}{\partial yz}\right)}^2\right)dxdydz,s.t.h\left({\mathrm{P}}_{Ti}\right)={\mathrm{P}}_{Ri},i=1,\dots ,M}$ where P_Ti is the target object, P_Ri is the reference object for the sets of corresponding points, and h is the bending energy function that minimizes the non-negative quantity of the interpolation of the integral bending norm or the integral quadratic variation E(h). The TPS method now decomposes each component into affine and non-affine components, such that (2)${\displaystyle h\left({\mathrm{P}}_h\right)=\Psi \left({\mathrm{P}}_h\right)\mathrm{K}+{\mathrm{P}}_h\Gamma }$ where P_h are the homogeneous coordinate points on the target 3D face, and Ψ(P_h) = Ψ₁(P_h), Ψ₂(P_h), …, Ψ_M(P_h) is a 1 × M kernel vector of TPS of the form: (3)${\displaystyle {\Psi }_w\left({\mathrm{P}}_h\right)=\parallel {\mathrm{P}}_h-{\mathrm{P}}_{Tw}\parallel }$ K is a M ×4 non-affine warping coefficient matrix, and Γ is a homogeneous affine transformation of a 4 × 4 matrix. The energy function is minimized to find the optimum solution to Eq. (4) if the interpolation condition in Eq. (1) is not met. (4)${\displaystyle \mathbb{E}\left(\beta ,K,\Psi \right)=\frac{1}{M}\sum _{J=1}^M\parallel h\left({\mathrm{P}}_{Tj}\right)-{\mathrm{P}}_{Rj}\parallel +\beta E\left(h\right).}$

The interpolation conditions in Eq. (1) are satisfied if the smoothing regularization term β is zero, where Γ and K are TPS parameters obtained by solving the linear equation: (5)${\displaystyle \left(\begin{array}{cc}\hfill \Psi \hfill & \hfill {\mathrm{P}}_R\hfill \\ \hfill {\mathrm{P}}_R^T\hfill & \hfill 0\hfill \end{array}\right)\left(\begin{array}{c}\hfill K\hfill \\ \hfill \Gamma \hfill \end{array}\right)=\left(\begin{array}{c}\hfill {\mathrm{P}}_T\hfill \\ \hfill 0\hfill \end{array}\right)}$

Ψ is a M ×M matrix with components Ψ_wl = ∥P_Tw − P_Tl ∥ and P_R is a M ×4 matrix in which each row is the homogeneous coordinate of the point P_Ri, i = 1 , …, M. Using Eq. (2), the target facial mesh P_Ti was deformed to the reference mesh P_Ri. The bending energy was applied, and the process was iterated for six cycles to achieve optimum sliding of the points on the facial surface which gives points relaxed. This changed the bending energy from the initial value E_i to the final value E_f after six complete iterations. This means that the semi-landmarks can be treated in the same way as homologous landmarks in downstream analyses. Since warping may result in points that do not lie directly on the facial surface on the target mesh, the transferred points were projected onto the closest point on the mesh surface using the ICP method (Creusot, Pears & Austin, 2010). The aim of using ICP is to iteratively minimize the mean square error between two point sets. If the distance between the two points is within an acceptable threshold, then the closest point is determined as the corresponding point (Mian, Bennamoun & Owens, 2008). The homologous landmark warping H_KΓ after six complete iterations is therefore: (6)${\displaystyle H_{K\Gamma }=E_{f-i}\left(\begin{array}{c}\hfill K\hfill \\ \hfill \Gamma \hfill \end{array}\right)}$

where (7)${\displaystyle \left(\begin{array}{c}\hfill K\hfill \\ \hfill \Gamma \hfill \end{array}\right)={\left(\begin{array}{cc}\hfill \Psi \hfill & \hfill {\mathrm{P}}_R\hfill \\ \hfill {\mathrm{P}}_R^T\hfill & \hfill 0\hfill \end{array}\right)}^{-1}\left(\begin{array}{c}\hfill {\mathrm{P}}_T\hfill \\ \hfill 0\hfill \end{array}\right)}$ is the linear TPS equation obtained during deformation of the surface of the target mesh to the reference mesh, before convergence was finally reached, and E_f−i = E_f − E_i after six complete iterations. The first iteration showed a partial distribution of sliding points on the target surface mesh (Fig. 4A). This was automatically repeated until the optimum homologous result was achieved, using an exponential decay sliding step of hundred to five percent. During relaxation of the spline, the semi-landmarks were slid along the surface and the curve tangent structures, rather than on the surfaces or the curves. This reduced the computational effort, as the minimization problem became linear since sliding along the tangents lets the semi-landmarks slip off the data (Gunz, Mitteroecker & Bookstein, 2005). The target surface mesh was then treated as a set of homologous points (Fig. 4B). Note that we did not construct a new deformable mathematical equation from scratch, but simply extended the standard deformable method established by Bookstein (1989). We added a minor extension to the computational process of projecting the surface semi-landmark from the template object to the target object and iteratively sliding the semi-landmark to a point relaxed.

After applying the step-by-step methods of digitization of the facial points and sliding of the semi-landmarks using ViewBox 4.0, we then applied Procrustes superimposition, error assessment with Procrustes ANOVA, canonical variate analysis, regression analysis, and shape visualization and variation with PCA, using MorphoJ 1.06d (Klingenberg, 2011). Boxplots for the size distribution and a MANOVA were performed in PAST 3.0 (Hammer, Harper & Ryan, 2001).

Sliding task comparison

In many studies dealing with sliding semi-landmarks, many iterations and tasks are required before optimum smoothness is reached, even with a small set of landmarks. Botton-Divet et al. (2015) compared Edgewarp (Bookstein & Green, 1994) and Morpho (Schlager, 2013) in terms of the time, workflow complexity, and computational efficiency required to slide semi-landmarks on long bones in different mustelids. The study used 27 manually placed 3D anatomical landmarks and 790 sliding semi-landmarks on curves and surfaces. Our study does not consider the time factor, due to the challenges arising from the different datasets and samples used. Furthermore, many studies that have implemented sliding semi-landmarks have not reported information on the sliding times, giving no scope for comparison. In our study, both the digitization of fixed landmarks and the sliding of semi-landmarks were performed in Viewbox. Unlike in Edgewarp, where digitization was first performed using a Breuckmann 3D white light fringe surface scanner (Botton-Divet et al., 2015), the remaining holes were filled and edges and spikes were removed using Geomagic. It should be noted that the task of surface pre-processing is not required in Viewbox, unlike in Edgewarp, where surface pre-processing, initial projection, and iterations against the Procrustes consensus must be accounted for (Botton-Divet et al., 2015).

Measurement error

The process of landmark coordinate extraction is always associated with some degree of measurement error, and this may be as a result of the non-coplanarity of landmarks, inconsistency of specimens relative to the plane of digitization, or difficulties in pinpointing the landmark locus (Webster & Sheets, 2010). Landmark digitization error can be minimized by careful landmark selection, but can never be totally eliminated. In assessing measurement error, three templates were designed for each population. The same individuals (five per dataset) were acquired three times each, using the three templates (Klingenberg, Barluenga & Meyer, 2002; O’Higgins & Jones, 1998; Robinson & Terhune, 2017). These were digitized for both manual and sliding semi-landmarks on the different reference objects, and this was followed by Procrustes superimposition on the landmark data using three partitions: fixed anatomical landmarks (FAL), sliding semi-landmarks (SSL), and combined landmarks (CL). Although several other error measurement methods were suggested by Fruciano (2016), the measurement error for this study was assessed using a Procrustes ANOVA. This technique (Klingenberg, Barluenga & Meyer, 2002; Klingenberg & McIntyre, 1998) was implemented in morphometrics to analyze measurement errors (Klingenberg et al., 2010; Leamy et al., 2015; Singh et al., 2012) using MorphoJ, which was achieved through the minimization of the squared sum of the distances of all objects and the consensus configuration (Fruciano, 2016). This approach uses a partition based on the sum of squares of the deviations from the average configuration of each coordinate in a two-factor ANOVA, which can then be summed across all the coordinates. Using a relevant number of degrees of freedom, computation of the mean squares is done by dividing the total sum of squares for an effect (Klingenberg & McIntyre, 1998). The Procrustes approach is most useful in computing mean shapes and in deciding whether two shapes [R_i] and [T_i] are random realizations of the same shape (Goodall, 1991).

The steps in this algorithm can be summarized as follows:

1. Anatomical fixed points (16) were digitized on the template facial mesh and a prominent point (the pronasale) was identified.

2. Semi-landmarks (484) were automatically generated and placed along the curves, located at a uniform distance along each curve for sliding in Step 5.

3. These semi-landmarks were first randomly placed and then uniformly distributed on the selected reference surface mesh, starting from the selected prominent point.

4. The reference facial model was warped to each target mesh configuration using a TPS transformation, and the surface semi-landmark was projected from the reference facial mesh to the target facial mesh.

5. The surface and curve semi-landmarks were then slid together in the direction that minimized the bending energy between each target configuration and the reference object. This was done iteratively in six complete cycles, in order to ensure convergence and optimum smoothness. This gave a homologous representation of the reference mesh.

6. A Procrustes superimposition of the landmark data was performed, and an error assessment was computed using a Procrustes ANOVA.

Shape and size variation

Since there may be an interaction between the size and shape in facial morphology due to changes in the shape associated with size differences (Klingenberg & McIntyre, 1998), we assessed the allometry by testing the statistical significant proportion of morphological variation in the shape components, using a multivariate regression of shape onto size. Canonical variate analysis (CVA) was also performed to test the group differences and an ordination plot was produced (plot not shown).

Differences in the effects and size were then examined by computing a non-parametric analysis of variance (MANOVA) in terms of Wilks’ lambda. Using the population as the group and the size as the covariate, the population by size interaction term was calculated. The MANOVA was recomputed after removing the interaction term (population by size) and the population effect tested for the difference in the regression intercept.

Sliding task

The tasks required for the initial projection and the relaxations against the Procrustes mean shape for the iterations are minimal in Viewbox, which requires a much simpler workflow compared to Edgewarp. Moreover, since a general Procrustes analysis was required in order to perform the sliding task, which is not implemented internally in Edgewarp, new input files were generated prior to performing the next iteration (Botton-Divet et al., 2015), creating greater complexity in the workflow.

PCA

After calculating the mean shape using Procrustes superimposition, we studied the variability in the shape using PCA. This was not only to reduce the number of dimensions but also to balance the fit of the model with the ease of analysis and potential loss of information (Cangelosi & Goriely, 2007). The PCA of the total sample yielded 168 principal components. When they were separately computed, the Stirling population yielded only 57 PCs, FRGC yielded 59 PCs, and Bosphorus 51 PCs, all with non-zero variability. In order to retain any component that accounted for a specific proportion or percentage, we used a broken stick approach to PCA selection (Cangelosi & Goriely, 2007; Klingenberg, 2013; Peres-Neto, Jackson & Somers, 2005). The first five PCs accounted for 98.05% of the total variation in the Stirling population, while the first four PCS accounted for 93.44% of the total variation in the FRGC population, the first three PCs accounted for 85.65% of the total variation in the Bosphorus population, and the first two PCs accounted for 97.17% of the total variation in the combined population (Cangelosi & Goriely, 2007). The first and second PCs (PC1 and PC2) of the Stirling population accounted for 65.40% and 20.03% of the variation, respectively; in the FRGC population, these accounted for 76.88% and 7.97%, respectively; for the Bosphorus population, these figures were 68.29% and 10.84%, respectively; and for the total population, they were 93.96% and 3.21%, respectively.

Shape variations for each population and the combined population are shown in Fig. 5. In the visualization of 3D dataset object, the number of landmarks is shown in red and the mean symmetry configuration is shown in light blue. A lollipop graph is shown of the first principal component of each population and the combined population, which indicates the difference in the face shape. Note that we only visualized the first PC of each population, since this accounted for the highest variation in the total shape.

Procrustes ANOVA and shape variation

A Procrustes superimposition of each set configuration produced a symmetric consensus configuration. When the variation was partitioned around this consensus using a Procrustes ANOVA (Table 2), the results showed that the variation in symmetric shape among individuals in each set accounted for the largest portion of the total variation, although this was not statistically significant for the manually placed landmarks. To assess the digitization errors of the manually placed landmarks, the sliding semi-landmarks and overall landmarks, the deviations for each were obtained by simply calculating the amount of displacement from the average position calculated from all digitization and the variation accounts for the smallest portion of the total variation. For manually placed (fixed anatomical) landmarks, the variation accounted for 4.31%, while for sliding semi-landmarks, this accounted for 3.71%, and for the overall landmarks, this was 3.67%.

The CVA indicated that each population studied was clearly distinct from the others. The Procrustes distances between populations (FRGC vs. Bosphorus = 0.040; Stirling vs. Bosphorus = 0.580; Stirling vs. FRGC = 0.061) were all statistically significant (p < 0.0001). All 10,000 pairwise permutation tests indicated that the mean shapes differed significantly in the population.

Figure 6A shows the first two principal components of the total shape in the three populations, accounting for the total variance. The scatterplots of the scores along the first two principal components for all the datasets (Fig. 6A) showed an apparent pattern of association between the manual anatomical landmarks and the sliding semi-landmarks for the same specimen. This pattern disappeared when the sliding semi-landmarks were removed, as the data for Bosphorus and FRGC were clustered more tightly.

Changes in the shape of the face as a function of size (allometry) were evaluated using multivariate regression of the effect of shape variables on the centroid size for all populations (Fig. 6B). This is usually of primary interest within a homogeneous population. Inter-population allometry explained only 13.28% of the shape differences according to size, and was significant (p < 0.0001). The box-plots in Fig. 6C showed that the Bosphorus dataset contained the largest sizes, followed by FRGC, while the Stirling dataset contained the smallest sizes.

The characteristics of the allometric trajectories of the population were then tested using a MANOVA (Table 3), which explained a significant proportion of the overall variation. The interaction term (test for slopes) was statistically significant. When the size effect was removed and the MANOVA was repeated, the result was still statistically significant, suggesting that the effect of size on shape for both the slope and intercept was strong, and that this was not the case in the population group.