A computerized facial approximation method for archaic humans based on dense facial soft tissue thickness depths

Facial approximation (FA) is a common tool used to recreate the possible facial appearance of a deceased person based on the relationship between soft tissue and the skull. Although this technique has been primarily applied to modern humans in the realm of forensic science and archaeology, only a few studies have attempted to produce FAs for archaic humans. This study presented a computerized FA approach for archaic humans based on the assumption that the facial soft tissue thickness depths (FSTDs) of modern living humans are similar to those of archaic humans. Additionally, we employed geometric morphometrics (GM) to examine the geometric morphological variations between the approximated faces and modern human faces. Our method has been applied to the Jinniushan (JNS) 1 archaic human, which is one of the most important fossils of the Middle Pleistocene, dating back to approximately 260,000 BP. The overall shape of the approximated face has a relatively lower forehead and robust eyebrows; a protruding, wider, and elongated middle and upper face; and a broad and short nose. Results also indicate skull morphology and the distribution of FSTDs influence the approximated face. These experiments demonstrate that the proposed method can approximate a plausible and reproducible face of an archaic human.


Introduction
Facial approximation (FA) or craniofacial reconstruction aims at recreating a potential facial appearance from a dry skull. This technique is often the last hope in the realm of forensic science when no other clues and evidence support the investigation and identification (Wilkinson 2010). Based on the assumed relationship between soft tissue and the bony structure, FA has been applied in archaeology to reconstruct the portraits and facial appearances of people in the past (Kustar 2004;Benazzi et al. 2009;Hayes, et al. 2017;Shui and Wu 2018;Marić et al. 2020). It has sometimes been applied to named individuals, but more usually unnamed people from the past. Nonetheless, these applications always focused on modern humans, and they are less commonly applied to archaic humans, where differences in skull and facial morphologies make the approximation more challenging. In recent years, the approximated appearance of our fossil ancestor has become an area of study for anthropologists and has also captured the imagination of the general public, influencing perceptions of how "like us" and how "human" Neanderthals were. The visualization of the approximated face rather than an imaginary approximation provides an effective 3D presentation to help us perceive and understand the characteristic features of human fossils. In addition, FA offers a new insight to investigate the morphological shape variations between archaic humans and Homo sapiens.
3D manual facial approximation approaches have been widely used to recreate facial appearances (Hayes 2016). Anthropologists collaborated with artists to recreate a possible likeness by means of modeling clay or plasticine over the replica of the skull and adding the facial features, e.g. eyes, nose, and mouth. During this procedure, muscle structures and facial soft tissue thickness depths (FSTDs) at anatomical landmarks can be used to represent the craniofacial relationship between soft tissue and skull. The manual FA approaches can be divided into three main categories: the Russian anatomical approach, the American anthropometrical approach, and the combination Manchester approach (Verzé 2009). However, they are heavily dependent on the degree of anthropological interpretation and the practitioners' subjective experience. Under such circumstances, multiple approximated faces of the same skull can be produced. For instance, three portraits of Ferrante Gonzaga, an Italian nobleman of the Renaissance, have been recreated (Fatuzzo et al. 2016). Such various approximations with inconsistent facial features might probably lead to less public confidence when no convincing hypothesis, and evidence can be provided.
With the rapid progress in computer science and medical image acquisition, computerized FA technology has been gradually developed to increase the level of accuracy and reliability of the approximated face. The basic idea is to mimic the manual FA approach using the computer (Wilkinson 2005). Using the FSTDs at anatomical landmarks and knowledge of facial muscle, both 2D and 3D interactive graphics technologies which mimic the Manchester approach have been developed. In 2D interactive FA, the frontal and profile portraits were recreated using Adobe Photoshop software (Hayes et al. 2012). Likewise, 3D interactive FA was used to recreate a 3D probable likeness through a haptic feedback device and 3D software, e.g. Autodesk 3ds Max, ZBrush, and Blender (Wilkinson et al. 2006;Lee et al. 2014;Short et al. 2014;Miranda et al. 2018). In their work, the tissue depth pegs which represented the FSTDs at anatomical landmarks were attached to the correspondence vertices of the dry skull, and the facial muscles were revised and attached to the surface of the skull. Then, the facial features were added and sculpted to improve the reliability of the approximated face. However, all these technologies require both anatomical knowledge and expertise in modeling skills. Anthropologists have to invest great effort in manual modeling when they wish to produce a range of multiple candidate faces that use the FSTDs of different samples.
Pioneering work on 3D graphical computerized FA was first proposed by Vanezis (1989). The average FSTDs at a limited number of anatomical landmarks were used to produce a coarse mask, and then the generic face was deformed to recreate a facial appearance over the dry skull. It is acknowledged that the greater number of FSTDs is acquired, the greater reliability of the approximation is achieved. Another effective computerized FA employed the deformation-based approach based on the assumption the verified craniofacial relationship of the template model is similar to that of the dry skull, removing the skull morphology variations (Quatrehomme et al. 1997;Nelson and Michael 1998;Turner et al. 2005;Deng et al. 2011;De Buhan and Nardoni 2018). In this procedure, either a generic face or a specific face based on the properties of the dry skull, e.g. age, sex, and ethnic group, was chosen as the template model. Then, the template face was deformed following the same transformation that was calculated by deforming the template skull to the dry skull. This approach is simple and easy-to-use, because it does not require the FSTDs table at anatomical landmarks. In recent years, with the increasing availability of skull and face datasets of modern living humans, a regression-based method has been applied to study the craniofacial relationship based on principal components (PC) scores of every skull and face in the shape space (Paysan et al.2009;Berar et al, 2011;Deng et al. 2016). Then, this predicted craniofacial relationship can be used to recreate the facial appearance.
With regard to the similarities in the craniofacial relationship between relatively recent modern human remains and modern living humans, the majority of existing studies focused mainly on approximating the appearance of archaeological human fossils. In contrast, only a few publications concerned the investigation of archaic humans (Hayes 2016 (Hayes et al, 2013). In their work, the 2D profile outlines of the approximated face were created based on the FSTDs at landmarks. Then, muscle images were deformed and attached to the surface of the skull. Finally, the reliability of the reconstructed face was evaluated using geometric morphometrics (GM). Because 3D facial morphology might allow anthropologists to better elucidate the facial characteristics of archaic humans and investigate evolutionary changes in the face, a 3D computerized FA approach still needs to be further investigated.
The Jinniushan (JNS) 1 cranium, dating back to 260,000 years BP, was discovered in Yingkou County, Liaoning Province, in northeast China in 1984 (Wu 1988;Rosenberg et al. 2006). It is one of the most important fossils in East Asia and has been used to investigate morphological features and shape variations with other fossils (Hublin 2013; Athreya and Wu 2017). It appears that its supraorbital shape, superciliary arch thickness and shape, postorbital constriction, and paranasal inflation are somewhat closer to those of Dali and Maba individuals, who are considered to represent population immigration from outside of China, and to be the result of an admixture with archaic humans (Andrews, 1986;Rightmire 1998). Although a manual approach has been used to produce the facial appearance of JNS 1, considerable interest has been shown in investigating the approximated face of JNS 1 based on reasonable assumption and supporting data, rather than experience and imagination. This paper aims to provide a computerized FA method to approximate the plausible and reproducible face of the archaic human.

The archaic human fossil
The JNS 1 cranium retained most of the maxillary dentition although the bone has been broken into more than one hundred pieces (Wu 1988). In an attempt to perform FA successfully, the cranium required careful examination and restoration. It has been manually repaired by researchers from the Institute of Vertebrate Paleontology and Paleoanthropology (IVPP) in Beijing, China. The restoration procedures were as follows: firstly, every fragment of JNS 1 fossil was cleaned and strengthened. Secondly, the fractured fragments were carefully matched together based on the similarity of the boundary of every fragment following the experience of the researchers. Thirdly, super glue was used to adhere fragments to each other. Finally, plaster was used to fill in the missing region of the cranium guided by geometric constraints. Anthropologists predicted the sex and age of JNS 1 through the analysis of morphological features, sutures, and dental wear. In recent years, JNS 1 was suggested to be female because of two important features, the subpubic concavity and the medial aspect of the ischiopubic ramus (Rosenberg et al. 2006). Likewise, based on the comparison of dentition and the analysis of tooth wear, an early study suggested that JNS 1 was over 30 years old (Wu 1988), but more recently, it was suggested to have been approximately about 20-30 years old (Herrera and Garcia-Bertrand 2018).
Because only the JNS 1 cranium remained and the mandible was not preserved, a well-preserved late archaic human mandible was required to assemble the JNS 1 cranium. But it remains challenging to find a suitable mandible with similar age and features. We have to decide to use two archaic human mandibles whose ages covered the age of JNS 1, i.e. one mandible is more recent than JNS 1, and the other is older than JNS 1 to repair JNS 1.
The Tabun 2 mandible was found in stratigraphic layer C of the Tabun cave, one of the paleoanthropological sites in the Near East. It was reconstructed and virtually recovered in six fragments, but lacking the left condyle, part of basilar symphysis (Schwartz andTattersall 2000, Quam andSmith, 1998). Morphologically, Tabun 2 is relatively large and robust, and it indicates a strong development of anterior marginal tubercle and a triangular basal corpus profile at the symphysis and mandibular foramina. It exhibited a mixture of morphological features of Neanderthals and early modern humans (Harvati and Lopez 2017). In addition, the well-preserved and complete Mauer 1 mandible (Wagner et al. 2010), a holotype of Homo heidelbergensis, was found near Mauer, southeast of Heidelberg, Germany in 1907. It is the oldest hominin fossil reported to date from central and northern Europe. It is of note that Mauer 1 exhibits a mixture of both primitive and modern features (Mounier et al. 2009). In this study, Tabun 2 and Mauer 1 were selected to fit with the JNS 1 cranium. Then, these two reassembled skulls (called JNS 1 using Tabun 2 and JNS 1 using Mauer 1) were used to approximate the face of JNS 1. Figure 1a shows the JNS 1 cranium (peach color) and the Tabun 2 mandible (gray color). Figure 1b displays the JNS 1 cranium (peach color) and the Mauer 1 mandible (gray color).

Skull and face datasets of modern living humans
In order to obtain the craniofacial relationship between soft tissue and skull, a total of 60 modern Chinese living humans (30 females and 30 males aged 20-30 years old), who lived in Shaanxi province in northern China, were enrolled in this study. More details can be seen in our previous studies (Shui et al. 2017). Each individual had normal morphology and had never undergone any orthodontic treatment. Medical images were acquired by means of a clinical multi-slice CT scanner system (Siemens Sensation 16). The CT images of each individual were archived in standard DICOM 3.0 with a resolution of 512 × 512 . All participants were provided with full details of the study and written informed consent. This research was approved by the Ethics Review Committee of Department of Archaeology, University of York.
Our previous studies constructed dense corresponding vertices among skulls and faces, respectively. The procedure was as follows: firstly, image segmentation and the well-known marching cubes algorithm (Lorensen and Cline 1987) were used to convert a series of CT images to the digital skull (or face). Secondly, the external surface of every skull (or face) within our dataset was computed. Thirdly, anatomical landmarks of the skulls and faces were defined and placed. Fourthly, iterative closest point (ICP), thin-plate splines (TPS), and compact support radial basis function (CSRBF) algorithms were applied to register all skulls (or faces). Next, the closest points were considered as corresponding vertices, i.e. every skull (or face), had the same number of vertices and each vertex of every skull (or face) was located approximately in corresponding positions. To remove the effects of location, orientation, and scaling, generalized Procrustes analysis (GPA) and principal component analysis (PCA) were carried out to construct the skull and face statistical shape model. Every skull (or face) can be represented by the coordinates of the average skull (or face) and the linear combinations of PC scores and corresponding orthogonal PCs (Shui et al. 2020).

Anatomical landmark definitions
In order to estimate the overall shape of the facial appearance, an anatomical landmark and semi-landmark configuration was defined. A total of 91 anatomical landmarks were chosen, and their 3D coordinates were acquired using Landmark Editor software (Wiley et al. 2005), where 17 anatomical landmarks were located on the midline and 74 anatomical landmarks were bilateral, respectively ( Table 1). Most of these anatomical landmarks were defined according to Martin's definitions (Martin 1928). Then, 404 semilandmarks were placed on JNS 1, which were identified in 16 patches based on the given landmarks. Here the semilandmarks of each patch were equally spaced within a 3 × 3 frame as a patch, and each patch with less than 9 anatomical landmarks was replenished with the middle points of two adjacent anatomical landmarks (Table 2). Finally, these semi-landmarks were projected to the modern human skull, i.e. we established geometric correspondences between two skulls (Gunz and Mitteroecker 2013).  Ramus posterius, most concave point on the posterior margin of ramus* Bilateral 47 Ramus anterius, most concave point on the anterior margin of ramus* Bilateral 48 Vertical projection from lowest point of mandibular notch to lower margin of mandible along ramus Bilateral 49 Vertical projection from alveolare of lower m2 to lower margin of mandible* Bilateral 50 Temporale anterius, most anterior point of temporal squama* Bilateral Figure 2 summarizes the framework of the proposed FA method. Firstly, the JNS 1 cranium and the selected mandible were virtually reassembled, and the missing geometry was repaired. Secondly, a coarse-to-fine computerized FA approach was proposed to recreate the possible likeness of JNS 1 based on the assumption that the distribution of average FSTDs of the modern humans within the dataset is similar to that of JNS 1. This procedure comprised four steps: (a) a hybrid non-rigid registration approach was carried out to establish the dense geometric correspondences between the template skull and JNS 1, where 495 landmarks and semi-landmarks were used to guide the transformation mapping; (b) the dense FSTDs of the template were calculated and visualized in a graphical format; (c) the coarsely approximated face of JNS 1 was recreated by assigning dense FSTDs to the corresponding vertices of JNS 1. The TPS interpolation function was used to improve the approximation. Due to the absence of the actual mandible, multiple approximations can also be mathematically calculated by interpolating the surfaces of the approximated faces based upon Tabun 2 and Mauer 1. (d) Quantitative evaluation was used to validate the reliability of the approximation through comparison of the distributions of FSTDs. Finally, we employed GM to examine the morphological shape variations between the approximated faces and modern human faces. We examined the effects of skull morphology and FSTDs on the approximated faces. All these methods were programed * Anatomical landmarks are defined by the authors and the rest of anatomical landmarks without special notice are from Martin (1928)  Alveolare of upper P3* Bilateral 53

Methods
Alveolare of lower m1* Bilateral 54 Mental foramen Bilateral 55 Mental laterale, turning point from mental to mandibular body on the inferior margin* Bilateral

Fig. 2
The pipeline of the computerized facial approximation method using C + + and Matlab 2019. Taking the approximation of JNS 1 using Tabun 2 as an example, we introduce the proposed method.

The restoration of JNS 1
Because the temporal mandible joint (TMJ) that connected the JNS 1 cranium and Tabun 2 was insufficiently accurate, and the left condyle of Tabun 2 seemed incomplete, the first step was to predict the missing geometry of Tabun 2 and match the cranium and mandible closely. Before the virtual restoration, the external surface of JNS 1 was extracted. This comprised three steps: firstly, JNS 1 was transformed into the Frankfort coordinate system based on the left porion, right porion, left orbitale, and the glabella. Secondly, an external point cloud was generated based on the cylindrical sampling algorithm. In this procedure, a couple of cross-section planes were generated between the bottom and top of JNS 1. For every cross-section plane, the external points were obtained by calculating the intersection points between JNS 1 and a set of given rays, starting at centroid of every cross section along equally spaced angle vectors. Finally, the external point clouds were converted to a set of triangular meshes. Subsequently, we employed the mirror restoration method to repair the external surface of Tabun 2 (Gunz et al. 2009). Figure 3 shows the external surface of JNS 1 using Tabun 2 that comprised the anatomical landmark and semi-landmark configuration.

Computerized facial approximation
In anthropology, it is widely accepted that facial surface has a close relationship to the bony structure and that the overall shape of the face can be approximated based on skull morphology and the craniofacial relationship between soft tissue and skull (Wilkinson 2005). We proposed a coarseto-fine FA approach to produce a reproducible and objective approximation.
Geometric correspondences between the template skull and JNS 1 We employed two steps to establish high-quality geometric correspondences between JNS 1 and the template skull. We assumed that all the landmarks and semi-landmarks of the template skull and JNS 1 were represented by where l denoted the number of anatomical landmarks and semilandmarks. A popular non-rigid registration TPS function was first conducted to deform the template skull to JNS 1. During the deformation, it enabled the bending energy of the function f p i = q i minimized (Bookstein 1989). TPS can be represented by affine transformation parameters and nonaffine warping parameters as the following linear equation: where the radial basis kernel can be represented by was the regularization parameter that used to balance the smoothness. I denoted the l × l identity matrix; O denoted the 4 × 4 zero matrix; A denoted the 4 × 1 zero matrix; α = [a 0 a 1 a 2 a 3 ] T ; and w = [ϖ i ] T represented the affine and non-affine parameters, respectively.
Following the same transformation that was computed by warping the template skull to JNS 1, the template face was deformed to produce a possible likeness as a candidate face. Next, we employed non-rigid registration to allow The restoration of JNS 1 using Tabun 2 that comprised the anatomical landmark and semi-landmark configuration Archaeol Anthropol Sci (2021) 13: 186 the deformed template skull and JNS 1 to match closely by assigning an affine transformation to every vertex of the deformed template skull. Assumed affine transformations X = [X 1 X 2 X 3 …X n ] T of all the vertices, we defined the cost function, E(X), which consists of anatomical landmarks and semi-landmarks term E l (X), a local affine regularization E d (X), and a stiffness term E s (X). To evaluate the accuracy of the skull match, the geometric deviation between the deformed template skull and JNS 1 was quantitatively calculated and depicted in a graphical format.
The cost function was as follows: where α, β, and λ denoted the weights that guided the optimization process. E l (X) was used to initialize and guide the registration as follows: where m i was the i-th landmark and semi-landmark of JNS 1 and v i was the i-th corresponding landmark and semi-landmark of resulting deformation of the template skull via TPS. The local affine regularization term expressed the distance between a vertex of JNS 1 and the corresponding vertex of the resulting deformation of the template skull as follows: where dist() denoted the distances between the corresponding points of JNS 1 and the resulting deformation of the template skull, and i denoted the reliability of the correspondences between these two meshes. We assumed that the nearest points between two meshes were the correspondences, denoted by rq i and rp i . In this procedure, the angles between normal vectors of the corresponding points and the Euclidean distance of the corresponding points can be used to improve reliability and reject the outliers.
The stiffness term was applied to regularize the deformation as follows: where ‖ • ‖ F was the Frobenius norm. i and j were the transformations of neighboring vertices, which were connected by an edge that belonged to the resulting deformation of the template skull. G = diag(1,1,1,γ) denoted a weighting matrix.
Dense FSTDs of the template During the acquisition of FSTDs of the template, the normal vector that was almost perpendicular to the surface of the bony structure was considered to be the measurement direction. A ray that started at a vertex of the template skull along the normal vector often passed through the template face; thus, the intersection point can be calculated. The FSTDs were defined as the Euclidean distances between pairs of corresponding vertices. It is of note that the normal vector of a given point that was determined by the geometric coordinates and topologies of the neighboring vertices influenced the accuracy of the FSTDs measurement. When the surface contained noise and sharp features, e.g. boundary of the surface, normal estimation remained a challenge. We extracted stable regions with robust normal estimation from the whole skull and then used FSTDs of the vertices within these stable regions to accomplish FA. It comprised two steps: firstly, we calculated the FSTDs of all the vertices along the closest distance vectors (Huempfner-Hierl et al. 2015). For every vertex of the template skull, the nearest point on the template face was searched, and the FSTDs were defined as the Euclidean distances between every pair of corresponding vertices. Secondly, the discrepancies between FSTDs along the normal vectors and those along the closest distance vectors were calculated. Once the deviation was less than the threshold, the vertex was suggested to be a stable vertex. Figure 4a shows the stable regions (red points) and unstable regions (blue points). In addition, the boundary vertices of the skull often were not considered to belong to the stable region. They can be extracted from the triangle meshes based on the assumption that the onering adjacent points of every boundary vertex cannot form a closed loop (Shui et al. 2020). Figure 4b shows the boundary vertices of the template skull (green points). Neither the unstable regions nor the boundary vertices were used to generate the coarsely approximated face. A coarse-to-fine facial approximation The overall shape of the facial approximation can be coarsely produced based on the dense FSTDs of the template using the following equation: where f i and s i denoted the geometric coordinates of the i-th vertex of the approximated face and JNS 1, respectively. V i represented the normal vector of the i-th vertex of JNS 1, and d i denoted the soft tissue thickness of the corresponding vertex of the template.
The predictions of FSTDs and corresponding measurement directions were always inconsistent with the actual ones; thus, the approximated face would be unsmooth. In addition, there always existed some voids, such as the eyes, nose, and cheeks. We employed a TPS function to warp the candidate face to the coarsely approximated face to improve the approximation. Because the position of the control point located on the two faces will greatly influence the deformation, we calculated the corresponding intersection points as control points based on the known anatomical landmarks and semi-landmarks of JNS 1.
Additionally, we offered a tool to mathematically calculate multiple approximations that simulated the mandible morphology changes using the following equation: where represented the interpolated approximation.
Evaluation of the reliability We validated the reliability of the approximated face by means of examining whether or not the distribution of FSTDs of the template was consistent with that of the approximated face. The FSTDs defined along the closest distance vectors were used, because they were insensitive to measurement direction and data noise (Gietzen et al. 2019). The FSTDs deviation between the template and the approximation was calculated and visualized.

Morphological shape variations of facial approximation
GM was carried out to capture the main features of the approximated faces and examine the geometric morphological variations between the approximated faces and modern human faces. GPA was first used to register all the vertices of the approximated faces and modern human faces, removing translation, rotation, and scaling (O'Higgins and Jones 1998). Thus, all the faces can be represented in the non-linear Kendall's shape space. The centroid size (CS) which is defined as the square root of the summed squared distances between all corresponding vertices and their centroid was calculated. Then, PCA was conducted on the Procrustes aligned coordinates to construct a facial shape tangent space. In this shape space, every sample was represented by the average face and the linear combinations of PC scores and corresponding independent orthogonal PCs. Next, Student's t test was carried out to verify the significant level of PC of interest between the modern faces and the approximated faces. Finally, a visualization technique was used to investigate the extent to which PC greatly explained the main patterns of morphological variation. In this process, two new faces along the positive and negative PC of interest were generated as follows: where and denoted the average face and weighting coefficient (it was set to 1 or − 1), and i denoted the standard deviation of the i-th PC, and i represented the i-th PC.

The effects of FSTDs and skull morphology on the approximated face
It is noted that skull morphology and the distribution of FSTDs are the two fundamental components of FA. We examined how the choice of FSTDs affected the approximated face. The approximated faces of the same JNS 1 were produced based on the average FSTDs of the females and males within our dataset. Then, the FSTDs deviation between two different approximations was calculated and depicted in a graphical format.
As the real mandible of JNS 1 was not survived, we investigated the effect of different mandibles on the approximated faces. Different approximated faces of JNS 1 using Tabun 2 and JNS 1 using Mauer 1 were produced based on the same distributions of FSTDs, respectively. Then, the geometric deviation between the approximations was used to examine the shape difference.

Facial approximation of JNS 1
Since JNS 1 is suggested to be a female, the average skull and face of the female group (Fig. 5a) was chosen as the template to approximate the facial appearance. Based on 495 anatomical landmarks and semi-landmarks, we first used the TPS deformation approach to approximate the facial appearance of JNS 1 using Tabun 2. Figure 5b and c show the deformed template skull and the candidate face. Figure 5d shows the template skull (gray color) and JNS 1 (peach color). It can be seen that the deformed template skull (8) ( ) = + 3 i i (black points) does not match JNS 1 (peach color) closely, as shown in Fig. 5e. The geometric difference between the deformed template skull and JNS 1 is calculated and depicted in Fig. 5f, where the Euclidean distance of every corresponding vertex was calculated. The average distance value of the vertex was 0.98 mm. The smaller deviation regions were located around the top of the cranial vault, whereas the largest geometric deviation regions were found around superciliary arches, the bottom of frontal, zygomatic arch, lateral mastoid process, the bottom of the occipital bone, partial regions of maxilla and teeth, and the greater wing of the sphenoid. It is noted that the approximated face might be inaccurate around these regions, and the conventional deformation-based FA approach needs to be improved. We proposed a coarse-to-fine facial approximation approach by attaching the dense FSTDs of the template to the corresponding vertices of JNS 1. The non-rigid registration was used to warp the deformed template skull to JNS 1 to generate a better deformation result (Fig. 6a). Figure 6b shows JNS 1 (peach color) and the deformation result (black points). The geometric difference between the deformed template skull and JNS 1 is calculated and visualized in Fig. 6c, where the Euclidean distance of every corresponding vertex was calculated. Almost 96.9% of overall vertices on JNS 1 showed the deviation within a deformation error of 1.0 mm, and the average deviation value was 0.54 mm. The area (black points) where the deviation was greater than 2.0 mm can be mainly observed around the boundaries of the bony structure. The region where the geometric deviation was greater can be observed around the superciliary arches and nasal bone. Dense FSTDs of the template that were calculated along the normal vector are visualized in Fig. 7a. It appeared that the FSTDs were almost distributed symmetrically along the midsagittal plane. The thinner FSTDs were observed at vertices located around the forehead, superciliary arches, nasal bone, and the cranial vault, and the thickest FSTDs were observed around the cheek region, the greater wing of the sphenoid, and the bottom of the occipital bone. Figure 7b and c display the point clouds and the triangle meshes of the coarsely approximated face. For every landmark and semilandmark of JNS 1, we calculate the corresponding points on the candidate face (left figure) and the coarsely approximated face (right figure), as seen in Fig. 7d. Figure 7e and f show the improved resulting approximation that is created by warping the candidate face to the coarse approximation using TPS. The result shows that the approximated face has a lower forehead, and a protruding, wider, and elongated middle and upper face. Also, it has robust eyebrows, a broad and short nose, and a wide mouth.
We compared the approximated faces using two different methods based on the comparison of the distribution of FSTDs. Figure 8a-c illustrate the distributions of FSTDs of the candidate face, the coarsely approximated face, and the improved resulting approximation, respectively. Figure 8d shows the FSTDs deviation between the candidate face and the template (Table 3). Almost 62.9% of overall vertices on the candidate face showed geometric deviation within a discrepancy of 1.0 mm, and the average deviation value was 1.17 mm. The area where the deviation was greater than 2.5 mm can be observed around the coronoid process, ramus, greater wing of the sphenoid, and zygomatic. The deviation of the vertex within the superciliary arches and nasal bone was greater as well. Figure 8e shows the FSTDs deviation between the coarsely approximated face and the template. Almost 87.7% of overall vertices on the coarse approximation showed the deviation within a discrepancy of 1.0 mm, and the average deviation value was 0.39 mm. It can be seen that the distribution of average FSTDs of the template was very close to that of the coarse approximation, except for the side of the teeth and mandibular condyle. Figure 8f illustrates the FSTDs deviation between the resulting approximation and the template. Almost 76.4% of overall vertices on the improved approximation showed the deviation within a discrepancy of 1.0 mm, and the average deviation value was 0.85 mm. The area where the deviation was greater than 2.5 mm can be observed around the coronoid process, both sides of the teeth and the greater wing of the sphenoid.
Because the FSTDs are the fundamental basis of the facial approximation in our method, we employed the FSTDs that were calculated using different methods to predict the overall shape of the coarsely approximated faces. Figure 9 shows the point clouds of the coarse approximation based on the normal vectors of the overall vertices, rather than the stable vertices. It appeared that much noise and outliers occurred around the coarsely approximated face. Figure 10a displays the FSTDs that are computed using a cylindrical sampling method . Figure 10b shows the coarsely approximated faces that are recreated by attaching these FSTDs to JNS 1 along the cylindrical sampling  Figure 10c shows the FSTDs deviation between the template and the approximation. Almost 60.5% of overall vertices on the approximation showed the deviation within a discrepancy of 1.0 mm, and the average deviation value was 1.3 mm. The area where the deviation was greater than 2.5 mm can be observed around the side of teeth, greater   . 9 The coarsely approximated face using the normal vectors of overall the vertices wing of the sphenoid, the bottom of zygomatic region, mandibular condyle, and superciliary arches. These results indicated the proposed FA method can accurately assign the FSTDs to the corresponding vertices of JNS 1.

Facial approximation of modern humans
In order to evaluate the reliability of the proposed FA approach, we approximated faces of two modern female skulls (01 and 03 skulls) and two modern male skulls (02 and 04 skulls). Then, we compare the FSTDs deviation between the approximated and actual faces, and the geometric variations between the approximated and actual faces, as shown in Fig. 11. From left to right, the leftmost column displayed the skull, and the next two columns illustrated the approximated face and the distribution of FSTDs. The middle two columns displayed the actual face and the distribution of FSTDs. The sixth column depicted the FSTDs deviation for every vertex between the approximated and actual faces (Table 4). Additionally, the geometric deviation between the approximated and actual faces was calculated (Table 5). The rightmost column displayed the geometric discrepancy of every vertex between the approximated and actual faces. Figure 11a shows the approximation of the 01 skull based on the average FSTDs of the females. Almost 59.3% of overall vertices showed the FSTDs deviation within 2.5 mm, and the average deviation value was 2.12 mm. Figure 11b shows the approximation of the 02 skull based on the average FSTDs of the females. Almost 63.5% of overall vertices showed the FSTDs deviation within 2.5 mm, and the average deviation value was 2.19 mm. Figure 11c displays the approximation of the 03 skull based on the average FSTDs of the males. Almost 80.8% of overall vertices showed the FSTDs deviation within a discrepancy of 2.5 mm, and the average deviation value was 1.43 mm. Figure 11d shows the approximation of the 04 skull based on the average FSTDs of the males. Almost 94.4% of overall vertices showed the FSTDs deviation within a discrepancy of 2.5 mm, and the average deviation value was 0.95 mm. A potential reason is that the distribution of the average FSTDs of the males is consistent with that of the 04 skull. However, the approximated chin of the 04 skull bore little visual resemblance to the actual chin, because the mental protuberance and tubercle seemed to be incomplete.
These results indicate that the overall shapes of the approximated faces bore resemblances to the actual faces. The area where the deviation was thinner can be observed around the forehead, nasal bone, and maxillary bone. But the greater deviation can be found around the cheeks, nose tip, zygomatic region, parietal and temporal regions, and mental protuberance. The distribution of FSTDs deviation was relatively consistent with that of the geometric discrepancy, but the average value of the FSTDs deviation was less than the geometric deviation. Thus, the proposed FA approach can be applied to modern humans, and the FSTDs deviation is a promising evaluation indicator to access the reliability of the approximated face.

The effects of FSTDs and skull morphology on facial approximation
We can recreate multiple approximated faces based on the distributions of FSTDs of different templates. Figure 12a shows another approximation of JNS 1 based on the average FSTDs of the males within our dataset. Figure 12b displays the FSTDs of the approximated face. The FSTDs deviation between two approximations based on the average FSTDs of the females and males is calculated, as shown in Fig. 12c. Almost 69.9% of overall vertices showed the FSTDs deviation within a discrepancy of 1.0 mm, and the average FSTDs deviation value was 0.81 mm. The area (blue, red, and black points) where the FSTDs deviation was greater than 1.0 mm can be observed around the mouth, nasal base, zygomatic, parietal and occipital regions, and cheeks.
Using the above FA method, we recreated the approximated faces of JNS 1 using Mauer 1. Figure 13a shows the approximated face based on the average FSTDs of   the females. Figure 13b shows the geometric variations between this approximation and the approximation of JNS 1 using Tabun 2. We also use the distribution of average FSTDs of the males to recreate an approximation of JNS 1 using Mauer 1, as shown in Fig. 13c. Figure 13d shows the geometric deviation between this approximation and the approximation of JNS 1 using Tabun 2. The approximation of JNS 1 using Mauer 1 has a wider face and a more robust chin than the approximations of JNS 1 using Tabun 2. These results showed that the geometric shape of the skull greatly influences the overall shape of the approximated faces. Based on the distribution of average FSTDs of the females, we employed the linear interpolation method to mathematically approximate three faces. Figure 14 displays three approximations using η = 0.25, η = 0.50, and η = 0.75.

Shape analysis of the approximated faces
GM analysis was conducted to capture the geometric features of the approximated faces. A total of 27 PCs accounted for over 95% of the morphological variance in the shape space. The first PC (PC 1) accounted for 31.6% of the morphological variance, and the second PC (PC 2) accounted for 13.3% of the morphological variance. Figure 15 shows the plots of the first two PCs of four different approximations (Fig. S1), including the approximated face that used the average FSTDs of the females and Tabun 2 (green point), the approximated face that used the average FSTDs of the males and Tabun 2 (yellow point), the approximated face that used the average FSTDs of the females and Mauer 1 (black point), the approximated face that used the average FSTDs of the males and Mauer 1 (cyan point), 30 modern female faces (red points), and 30 It is of note that PC 1 (p < 0.005) and PC 2 (p < 0.05) have significant differences between modern human faces and the approximated faces.
To identify the main patterns of shape variance, four new faces were recreated along the positive and negative directions of PC 1 and PC 2. The positive PC 1 connected with the approximated face with a relatively lower forehead, and robust and wide eyebrows; a protruding, wider, and elongated middle and upper face; a broad and short nose; a wider mouth; and a robust chin. By contrast, the negative PC 1 represented a face with a prominent and protruding forehead, a narrower middle and upper face, and relatively narrower mouth, nose, and chin. All the approximated faces were located at the extreme positive end of PC 1. It indicated that the approximated face was greatly different from modern human faces and verified the characteristic features of the approximations. When the same FSTDs were performed, the approximated faces that used Mauer 1 were located on the right side along the positive PC 1. It indicates that these approximations have wider cheeks and robust chins.

Discussion
Many previous studies have employed manual and computerized FA approaches to modern humans in archaeology and anthropology (Claes et al. 2010;Wilkinson 2010). However, when FA is applied to archaic humans, the main challenges are the poor preservation of archaic human fossils and the absence of their craniofacial relationship and anatomical knowledge. In this study, we propose a computerized coarse-to-fine FA approach to attach the distribution of average FSTDs of modern humans to archaic humans. The resulting approximation is promising, objective, and repeatable, and the reliability can be evaluated through a quantitative comparison of the distributions of FSTDs. Furthermore, we investigate the effects of skull geometry and the distribution of FSTDs on the approximated face. The first stage of FA is to examine and restore the dry skull. When the skull had missing parts or distortion, the TPS function was always used to deform the reflection structure of the intact side to fill in the gaps based on a landmark and semi-landmark configuration. Other studies attempted to use computerized approaches to virtually reassemble the skull fractures together (Yu et al. 2012). In the worst case, the fossil specimen, as in the case of the mandible of JNS 1, had not survived. In these cases, there is a potential solution, using a mandible that has similar age and morphological features to match the cranium. But since the mandibles of archaic humans are rarely found, and the geometric shapes of the mandibles are always unique, it remains challenging to select an appropriate mandible to provide a good fit with the JNS 1 cranium. We attempted to use different mandibles to match the JNS 1 cranium and recreated different approximations based on the repaired skull. According to these approximated facial appearances, multiple approximations can further be mathematically recreated to provide some references to describe the overall shape. Although these approximations cannot be interpreted anatomically, they might provide a new perspective for researchers to illustrate the facial appearance.
During FA, the prediction of the craniofacial relationship between soft tissue and the dry skull and the assignment of the predicted craniofacial relationship are two fundamental questions. The average FSTDs at landmarks and muscle structures of modern humans are always suggested to be the craniofacial relationship of people in the past (Hamre et al. 2017). In the study previously mentioned, the FSTDs of chimpanzees can be used to depict a thinner mid-face of archaic humans because FSTDs around the cheek were almost half that of modern humans (Hayes et al. 2013). However, due to the lack of evidence, there is a particular challenge to decide which FSTDs are confidently reasonable for JNS 1. This study considers the distribution of average FSTDs of modern humans as that of the archaic human and strongly suggests that FSTDs of the reliable regions along the normal vectors are appropriate for approximating the facial appearances of archaic humans.
It is worth mentioning that the FSTDs make a great contribution to improving the accuracy and reliability of the approximated face (Claes et al. 2006;Starbuck and Ward 2007). It is acknowledged that many factors, such as sex, age, nutrition status, body mass index, and ethnic groups, will impact the FSTDs of landmarks. Previous studies always placed the landmarks on the skull manually and then acquired the FSTDs at these landmarks (Stephan 2017). We investigated the effects of the distribution of FSTDs on FA. Approximated faces of the same JNS 1 using different FSTDs shared a resemblance. The larger deviation areas of two approximated faces between using FSTDs of the males and females were almost consistent with the regions of the FSTDs discrepancy between males and females . Once a reliable and verified craniofacial relationship can be obtained, a more reliable and accurate approximation would be recreated.
Because of the greater shape variations between the modern and archaic human skulls, it remains challenging to assign the FSTDs of the template to those of JNS 1. Previous studies often employed the deformation-based FA approach (Nelson and Michael 1998;Turner et al. 2005;Deng et al. 2011). It is accepted that the closer the deformed template and the dry skulls matched, the more confidence there was in the reliability of the approximated face. However, the performance of the TPS deformation even incorporating a regularization is inadequate accuracy when two skulls have great geometric differences. To address this problem, this study employed a hybrid non-rigid registration algorithm to establish a high-quality set of geometric correspondences. Additionally, this approach can make use of FSTDs for every sample within the dataset to recreate a range of multiple approximated faces and then use PCA to construct a tailored approximation-space for JNS 1. In this context, the missing areas of the coarsely approximated face can be better predicted (Gietzen et al. 2019), and a range of possible approximation can be recreated by using appropriate coefficients of PCs of interests (Shui and Wu 2018).
As previous studies mentioned (Oxnard and O'Higgins 2009;Wärmländer et al. 2019), the anatomical landmarks and semi-landmarks should be very carefully designed with regard to the research question. The purpose of this study is to approximate the overall shape of the facial appearance. Previous studies employed different numbers of landmarks (Vandermeulen et al. 2006;Deng et al. 2011) and topographic features, e.g. crest lines (Turner et al. 2005), to guide the deformation. But there is no standard criterion for performing FA.
We recommend that landmarks and semi-landmarks need to cover the entire skull, particularly around the region where the template and dry skulls are quite different. Such definitions will improve the certainties of the establishment of dense geometric correspondences and so enhance the reliability of the approximated face. Additionally, we employed GM to capture the main features of the approximated faces. The dense corresponding vertices are used to provide better visualization and interpretation, rather than the use of landmarks and semilandmarks. But due to the complexity of biological structures, we need to carefully examine the effect of landmarks, semilandmarks, and high-density correspondences on the shape analysis and then decide which type of corresponding points can provide a reliable interpretation.
The geometric comparison of the actual and approximated faces has been used to validate the reliability of FA. In this process, the shell-to-shell deviation, surface-to-surface deviation, and the shortest distance between two point clouds have been conducted (Wilkinson et al. 2006;Short et al. 2014;Miranda et al. 2018). But, the resulting registration impacts the comparisons of discrepancy between the approximated and actual faces. Additionally, since no verified actual faces can be provided, the comparison of the actual and approximated faces cannot be used to archaic humans. Thus, the reliability of the approximated face is mainly evaluated based upon the experience and knowledge of the experts. Even though the experts marked the areas of the approximated faces that need to be improved, researchers still do not know how to revise them exactly. To tackle this problem, we can integrate the approximation of facial appearance and recognition of the less confidence in approximated regions iteratively. This method provides a promising tool to allow researchers to examine the extent to which the resulting approximation is unreliable and where it needs to be further revised to improve the reliability.
However, uncertainty still remains since there are no actual archaic human faces and we cannot know the actual relationship between soft tissue and skull of archaic human. Also, it remains challenging to understand the relationship between facial features, e.g. eyes, nose, mouth, and ears, and bony structures. Thus, the approximated features need to be further improved, e.g. the size and orientation of nasal aperture. In addition, the ages and features of Tabun 2 and Mauer 1 are different from JNS 1, and an appropriate mandible needs to be used to enhance the reliability of the approximated face.

Conclusion
The approximated face of JNS 1 is a typical case of interdisciplinary study that provides anthropologists and the general public with an improved visual interpretation of the facial morphology of archaic human. This study proposed a coarse-to-fine FA approach based on dense FSTDs of modern humans and presented an evaluation approach to validate the reliability of FA. We also investigated the effects of skull morphology and the distribution of FSTDs on the approximated faces. Since the mandibles of archaic humans are rarely found, it remains challenging to select an appropriate mandible for the JNS 1 cranium. In the future, we will attempt to collect the different mandibles of archaic humans to improve the reliability of the approximation of JNS 1.