Estimating the age of death of unknown human skeletal remains represents one of the major tasks of biological anthropologists. Traditionally, the estimation is performed visually by assessing degenerative changes of join surfaces (e.g. [8], [10]), among which the pubic symphysis of the pelvis is widely used (e.g. [4], [11]). Figure 1 illustrates a few examples of human symphyseal surfaces of individuals with the age-at-death of 25, 35, 45, 65 and 85 years, respectively. Visual observation, however, has its limitations, e.g. it is subjective, dependent on observer experience, and last but not least, its applicability suffers from low accuracy and reliability of estimates (e.g. [7], [9]). To achieve both accurate and reliable age estimates, it is recommended to use three broad intervals [1], [5]. Moreover, single-indicator methods do not work equally well throughout the adult period, for example, it has been reported that the pubic symphysis is no longer suitable for age estimation after the age of 40 years [2], [6]. Currently, the research has shifted to imaging technologies and sophisticated data mining methods (e.g. [3], [15]) that could offer a more objective and accurate perspective on age estimation in adults. The Algee-Hewitt -Slice -Stoyanova team ( [12], [13] proposed the most prominent approach [14]) with the estimation error (RMSE) ranging between 13.7 and 16.6 years (based on the dataset consisting of 93 samples) [13]. In this paper, we did not follow that approach. Instead, we developed a novel age-at-death estimation model based on CNN. Our model takes an image representing the 3D structure of the bone (see the following sections) and predicts an individual's age-at-death using the CNN for pattern recognition.
Table 1 shows the age-at-death distribution of our collection. The mean age is 53.7 years, and the standard deviation is 17.1 years. Table 2 shows the structure of the osteological collection. The input dataset consists of 483 skeletal samples from adult (18-92 years) males and females. All skeletal samples were digitised using the HP 3D Structured Light Scanner PRO S2 or S3 scanner and exported in STL format. The STL format is a file format describing an unstructured triangulated surface using a 3D Cartesian coordinate system. In this format, the surface geometry of a 3D object is represented as a number of small adjacent triangles. Figure 2 shows an example of the 3D scan.
The input data coming from the scanner needs to be preprocessed before applying the proposed age estimation method. The STL format is a very convenient input data format, however, it is relatively unsuitable for direct use due to the irregular distribution of all vertices in the 3D space. There are several options how the input data could be represented and regularised. For instance, one can use voxel representation, where each voxel encodes one bit of information -the presence of the bone. However, for a typical scan dimension 50×15×15 mm (see Figure 2) with the resolution of 0.1 mm, this results in over 11 million of voxels per one 3D scan, i.e. in about 11.3 Mbit of information. On the other hand, one can use only the top view of the symphysis surface (see Figure 3), since this area has the highest age prediction capabilities [4], [11]. The surface height could be encoded in colour or in grayscale. Using 8-bit grayscale gives even better resolution for the surface height, compared to the previously described voxel representation. This approach reduces the overall size to 0.6 Mbit, while keeping the same resolution for the other two dimensions. However, this "top view" representation ignores the side walls of the scan, and therefore, eliminates potentially additional age-related information. Therefore, we decided to transform the input data from the Cartesian coordinate system to a new coordinate system in such a way that the side walls of the scan could be examined in a similar way to the "top view" representation. First, the position, size, and orientation of all scans needs to be standardised. Then, every point (x, y, z) in the Cartesian coordinate system is transformed to a point (λ , µ, ν) in the new coordinate system using the following equations
where θ = atan(y/x), ϕ = acos(z/ρ), ρ = x 2 + y 2 + z 2 .
(
Here, we should note that θ = atan(y/x) is the four quadrant arctangent of the elements of x and y such that −π ≤ atan(y/x) ≤ π. The concept of this coordinate system transformation is illustrated in Figure 4. Figure 4 shows some point P on the surface S, where the surface S represents the surface of the symphysis. Initially, the point P is transformed into the spherical coordinates using the equation ( 2). Consecutively, this point is further transformed into the point P' using the equation ( 1). If all points obtained from the scan are transformed in such a way, then the surface S is transformed to the new surface S'. The above-described transformation has several interesting properties and offers multiple advantages. Figure 5 helps to understand these properties. Figure 5 (a) shows a cube placed at the centre of the Cartesian coordinate system. Figure 5 (b) shows the same cube but in the new coordinate system. Since the transformation preserves ν = z, all points have the same height above the x-y plane, or λµ plane, respectively. The cube is virtually stretched out from the bottom side of the cube in that way that the entire cube can be described as a function of the two variables λ and µ, i.e., for each point (λ , µ) in a portion of the λ -µ plane (the domain of the function) we can assign a unique number f (λ , µ). This is very advantageous, since the complicated 3D shape can be transformed to a 2D image practically without the loss of information.
Another advantage of the proposed transformation is a consequence of preserving ν = z. As already mentioned, all points have the same height above the x-y plane, or λµ plane, respectively. This allows us to detect and analyse the disturbances in the symphyseal surface profile quite This can be seen in Figure 5 (b), where the bottom edge of the "square" is much wider compared to the upper edge of that "square". Moreover, the bottom side of the original cube is completely deformed and rather resembles a ring, as seen in Figures 5 (b) in the dark blue areas of the image. However, this disadvantage is of little significance for our purposes, since all pubic scans have no bottom (3D scan captures only the surface of the bone, not internal parts of the bone), and the coordinate system is located in such a way that the most important areas of the scan are deformed only slightly. Figure 6 shows the distribution of a set of points in the λ -µ plane, which were originally uniformly placed on the surface of the cube. This helps to visualise the deformation of the cube shape. The top side of the cube is located around the origin of the λ -µ plane. The edges of the front side of the cube are highlighted in green. As can be seen, the bottom edge of the cube is more stretched compared to the top edge. The second disadvantage of the proposed coordinate system is more fundamental for complicated shapes, as 3D scans can be. Namely, not for all shapes, we can assign a unique number f (λ , µ) in the λµ plane. All points with the same value of θ and ϕ (for instance, points P, Q and R in Figure 4 (a)) are projected into the same (λ , µ) coordinates. In this case, we can select the maximum, minimum, median or the average of all points mapped to the same (λ , µ) point. In this situation, some information from the original shape is lost and cannot be fully recovered anymore. This creates unwanted artefacts in the transformed data. However, we have experimentally observed that it occurs only occasionally for our dataset and affects only small portions of the whole area. In our case, these artefacts are partially suppressed by scaling down all x-coordinates by a factor of 2.5 before applying the above described transformation of the coordinate system. Figure 7 shows the symphysis surface from Figure 2 in the new coordinate system. The ϕ variable is plotted with a resolution of 2 • for better visualisation. The actual resolution is set to 0.5 • . The surface from Figure 7 is projected onto a regular mesh in the λ -µ plane, where the ν coordinate is encoded in 8-bit (or 16-bit) value, effectively creating a grayscale image as shown in Figure 8. The range of angle ϕ can be arbitrarily chosen, e.g., if chosen such that ϕ ∈< 0 • , 90 • >, then only points above the x-y plane (with a positive z value) are used. The grayscale image can be directly used as input to the age estimation model. Moreover, to increase the variability of the input training dataset and the robustness of our model, we have generated 41 projections (grayscale images) for each 3D scan with a slightly rotated and translated origin of the Cartesian coordinate system. Our age estimation model consists of several identical age predictors. Each predictor is based on convolutional neural network [16]. Figure 9 shows the main idea of the age estimation flow. First, the 3D scan is transformed into several grayscale images (see Figure 8 as an example). These images are consecutively directly used as input for individuals age predictors. Second, an aggregation function is applied in order to combine the results from all predictors, and thus, to provide the final prediction. We have chosen mean and median as two possible aggregation functions. There is a single 3D scan for which we build multiple projections (41 in this case). By applying the predictor for each of the projections, we obtain multiple age predictions that are finally aggregated to gain the final predicted age.
The application of convolutional neural networks was an obvious choice. We have experimented with various architectures. In the final experiments, we were mainly inspired by the setup used for X-Ray images processing [17].
We also experimented with topologies based on DenseNet [18,19], which exploits a specific topology that shortens layer connections by connecting each layer to every other layer in a feed-forward fashion. The final predictor structure is shown in Figure 10. The model consists of a total of 20 layers. In the first part, the input image is reduced and transformed into features using convolutional layers in combination with pooling, activation (Ramp) and regularisation layers. The second part represents densely con-nected feedforward networks that transform 1024 features into a single real value. For better generalisation, we used a dropout layer in the second part.
The model was implemented in Wolfram Mathematica using an in-house neural networks package built on the MXNet framework.
To improve the robustness of the network training, we incorporated an ImageAugmentationLayer to the input layer. This layer takes the input image 160x160 pixels and randomly crops it to 128x128 pixels during the training phase -this allows us to efficiently expand the training dataset without having to implement a custom batch function. During the evaluation phase, the ImageAugmentationLayer crops the input image around the centre in a deterministic way, so it does not affect it during evaluation. We have chosen slightly larger kernels (7x7) in the first two layers to better handle larger structures in blurred images. We use the batch normalisation layers which are proposed as a technique to help coordinate the updating of multiple layers in the model [16]. Figures 11 and 12 are used as examples illustrating the extracted patterns for selected layers of the network for 20-and 72-years-old individuals, respectively. These figures show the input image and the output from layers #3, #6, #10 and #19. As it can be seen, the model identifies the vertical structures and the edge of the symphyseal surface reasonably well from the input image. These vertical structures combined with the shape of the symphyseal edge are also used by experts to identify the age of the individual.
As already mentioned, we have multiple images for a single individual (41 images per bone). For some individuals, we have at our disposal left and right bones, so there are 82 images for a given individual. All models are designed to process a single image on input, so the training set is an unstructured list of pairs {image, age}. Since we need proper testing, all images for a single individual must be in the same fold. Cross-validation uses the information about the ID of the individual to split the dataset properly. The folds are then flattened, and we perform standard training for mapping images to real values (age).
We use a standard Adam optimiser [20] with a batch size 16 (experimentally chosen) running on the GPU for training. Each training runs approximately 300 rounds (we have performed many experiments from 100 to many thousands), representing almost 10 million processed records.
To compare the models and study the behaviour, we have performed two ways of model evaluation. Obviously, the main goal is to predict the age of an (unknown) individual based on the 3D scan of the bone, so we evaluated the model for each individual (41 or 82 images) and computed the predicted age using aggregation. Figure 13 shows the actual age vs. the predicted age per individual. Based on the size of the dataset, we have chosen 5-fold validation for all of our experiments.
Stoyanova et al. presented in their study [13] five age-estimation models with similar age estimation performance and provided an open source software called forAGE (available at http://morphlab.sc.fsu.edu/). The best model (according to their study) is "SAH&VC" (SAH+Outline) and provides MAE of approximately 10.79 for their entire dataset (93 samples). To compare our model with their state-of-the-art model, we used their software and evaluated the MEA for our dataset as well.
The results are summarised in Table 3. It should be noted here that for our models, MAE and RMSE are computed using 5-fold cross-validation, whereas for the model from [13] MAE and RMSE are computed directly. To determine whether our model contains any systematic error or whether any particular age intervals introduce some anomalies in the prediction, we processed the predicted ages per one-year age intervals. Figure 14 shows the variation in age predictions for each age class. As can be seen, our model generally overestimates younger individuals (under the age of 55 years) and underestimates mature ones (above the age of 55 years). This is a result of the tendency to predict the age towards the mean age of the sample. Similarly, we analysed the SAH&VC model from [13]. Figure 15 shows the variation of age predictions for particular age class. As can be seen, the SAH&VC model generally underestimates all individuals above 40-50 years. More specifically, for individuals over 45 years old, the average of all estimations reaches only 37.7 years (for our dataset). We believe that this primarily results from the unbalanced age distribution of the dataset used in [13]. Moreover, we observed that our model can estimate the age-at-death of an individual over the entire age interval (in our case between 19-92 years) -see Figure 14. This contrasts with e.g. [2] and [21], where pubic symphysis is considered appropriate for individuals up to 40 years, or 60 years, respectively. When the maturation process of pubic symphysis is complete, the morphological changes are degenerative and highly variable between individuals [2], [22]. However, we believe that our model can capture even such changes.
We have developed a novel age-at-death estimation model based on convolution neural networks. Our model provides a mean absolute error of approximately 10.6 years and is suitable for adult and mature individuals. Our results indicate that the pubic symphysis reflects the age of an individual throughout their entire adult life. In other words, we have observed no limitations in terms of age prediction capabilities of pubic symphysis of adult individuals.
Acknowledgments This research was supported by a research grant awarded by the Technology Agency of the Czech Republic; project number TL03000646.