Novel data mining-based age-at-death estimation model using adult pubic symphysis 3D scans Zdeněk Buk1 , Michal Štepanovský1 , Anežka Kotěrová2 , Jana Velemínská2 , Jaroslav Brůžek2 , and Pavel Kordík1 1 Faculty of Information Technology, Czech Technical University in Prague, Thakurova 9, Prague 160 00, Czech Republic 2 Department of Anthropology and Human Genetics, Faculty of Science, Charles University, Vinicna 7, Prague 128 43, Czech Republic Abstract: The paper introduces a novel age-at-death es- representing the 3D structure of the bone (see the follow- timation model based on Convolutional Neural Network ing sections) and predicts an individual’s age-at-death us- (CNN). The model uses 3D scan of human pubic symph- ing the CNN for pattern recognition. ysis as an input and estimates the age-at-death of the in- dividual as an output. The Mean Absolute Error (MAE) of this model is about 10.6 years for individuals between 18 and 92 years of age-at-death. Moreover, the results of the study indicate that pubic symphysis can be used to es- timate the age of individuals across the entire age range. The study involved a sample of 483 bone scans collected from 374 individuals (from which 109 individuals pro- vided both left and right pubic symphysis). 1 Introduction Figure 1: From left to right, examples of symphyseal sur- face scans of individuals with the age-at-death of 25, 35, Estimating the age of death of unknown human skeletal 45, 65 and 85 years. remains represents one of the major tasks of biological an- thropologists. Traditionally, the estimation is performed visually by assessing degenerative changes of join surfaces 2 Input data (e.g. [8], [10]), among which the pubic symphysis of the pelvis is widely used (e.g. [4], [11]). Figure 1 illustrates Table 1 shows the age-at-death distribution of our collec- a few examples of human symphyseal surfaces of individ- tion. The mean age is 53.7 years, and the standard de- uals with the age-at-death of 25, 35, 45, 65 and 85 years, viation is 17.1 years. Table 2 shows the structure of the respectively. Visual observation, however, has its limita- osteological collection. tions, e.g. it is subjective, dependent on observer experi- ence, and last but not least, its applicability suffers from Table 1: Age distribution of the collection low accuracy and reliability of estimates (e.g. [7], [9]). Age-at-death: 18–29 30–39 40–49 50–59 60–69 70+ To achieve both accurate and reliable age estimates, it is Males 24 54 57 56 59 49 recommended to use three broad intervals [1], [5]. More- Females 10 32 33 39 21 49 over, single-indicator methods do not work equally well throughout the adult period, for example, it has been re- ported that the pubic symphysis is no longer suitable for Table 2: Structure of the osteological collection age estimation after the age of 40 years [2], [6]. Currently, Country: Portugal Switzerland Thailand Crete the research has shifted to imaging technologies and so- Males 129 45 114 10 phisticated data mining methods (e.g. [3], [15]) that could Females 91 21 68 5 offer a more objective and accurate perspective on age es- timation in adults. The Algee-Hewitt – Slice – Stoyanova The input dataset consists of 483 skeletal samples from team ([12], [13] proposed the most prominent approach adult (18–92 years) males and females. All skeletal sam- [14]) with the estimation error (RMSE) ranging between ples were digitised using the HP 3D Structured Light 13.7 and 16.6 years (based on the dataset consisting of 93 Scanner PRO S2 or S3 scanner and exported in STL for- samples) [13]. In this paper, we did not follow that ap- mat. The STL format is a file format describing an un- proach. Instead, we developed a novel age-at-death esti- structured triangulated surface using a 3D Cartesian co- mation model based on CNN. Our model takes an image ordinate system. In this format, the surface geometry of Copyright ©2021 for this paper by its authors. Use permitted under a 3D object is represented as a number of small adjacent Creative Commons License Attribution 4.0 International (CC BY 4.0). triangles. Figure 2 shows an example of the 3D scan. µ, ν) in the new coordinate system using the following equations λ = ϕcos(θ ), µ = ϕsin(θ ), (1) ν = z, where θ = atan(y/x), ϕ =pacos(z/ρ), (2) ρ = x2 + y2 + z2 . Figure 2: Example of the 3D scan of the symphysis pubica Here, we should note that θ = atan(y/x) is the four of a 25-year-old female quadrant arctangent of the elements of x and y such that −π ≤ atan(y/x) ≤ π. The concept of this coordinate sys- tem transformation is illustrated in Figure 4. Figure 4 3 Data preprocessing shows some point P on the surface S, where the surface S represents the surface of the symphysis. Initially, the The input data coming from the scanner needs to be pre- point P is transformed into the spherical coordinates us- processed before applying the proposed age estimation ing the equation (2). Consecutively, this point is further method. The STL format is a very convenient input data transformed into the point P’ using the equation (1). If all format, however, it is relatively unsuitable for direct use points obtained from the scan are transformed in such a due to the irregular distribution of all vertices in the 3D way, then the surface S is transformed to the new surface space. There are several options how the input data could S’. 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 infor- mation. 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. (a) First step (b) Second step Using 8-bit grayscale gives even better resolution for the surface height, compared to the previously described voxel Figure 4: Transformation from the Cartesian coordinate representation. This approach reduces the overall size to system 0.6 Mbit, while keeping the same resolution for the other two dimensions. However, this "top view" representation The above-described transformation has several inter- ignores the side walls of the scan, and therefore, eliminates esting properties and offers multiple advantages. Figure 5 potentially additional age-related information. helps to understand these properties. Figure 5 (a) shows a cube placed at the centre of the Cartesian coordinate sys- tem. Figure 5 (b) shows the same cube but in the new co- ordinate 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 en- tire cube can be described as a function of the two vari- Figure 3: Top view on the symphysis surface from Fig- ables λ and µ, i.e., for each point (λ , µ) in a portion of the ure 2. The high of the surface is encoded in colour. λ -µ plane (the domain of the function) we can assign a unique number f (λ , µ). This is very advantageous, since Therefore, we decided to transform the input data from the complicated 3D shape can be transformed to a 2D im- the Cartesian coordinate system to a new coordinate sys- age practically without the loss of information. tem in such a way that the side walls of the scan could be Another advantage of the proposed transformation is a examined in a similar way to the "top view" representa- consequence of preserving ν = z. As already mentioned, tion. First, the position, size, and orientation of all scans all points have the same height above the x-y plane, or λ - needs to be standardised. Then, every point (x, y, z) in the µ plane, respectively. This allows us to detect and analyse Cartesian coordinate system is transformed to a point (λ , the disturbances in the symphyseal surface profile quite some information from the original shape is lost and can- not be fully recovered anymore. This creates unwanted artefacts in the transformed data. However, we have ex- perimentally 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 coor- dinate system. (a) Cube in the Cartesian coordinate system (b) Cube in the proposed coordinate system Figure 6: Visualisation of a set of points originally uni- Figure 5: Example of the proposed coordinate system formly located on the cube surface in the Cartesian coor- transformation dinate system after the projection into the λ -µ plane. easily. This transformation, however, also has few draw- Figure 7 shows the symphysis surface from Figure 2 in backs. First, it does not preserve the global shape of the the new coordinate system. The ϕ variable is plotted with surface as seen from a perpendicular view to that surface. a resolution of 2◦ for better visualisation. The actual reso- For instance, perfectly square side-walls of a cube (Fig- lution is set to 0.5◦ . The surface from Figure 7 is projected ure 5 (a)) become increasingly stretched out as ϕ grows. onto a regular mesh in the λ -µ plane, where the ν coordi- This can be seen in Figure 5 (b), where the bottom edge nate is encoded in 8-bit (or 16-bit) value, effectively cre- of the “square” is much wider compared to the upper edge ating a grayscale image as shown in Figure 8. The range of that “square”. Moreover, the bottom side of the original of angle ϕ can be arbitrarily chosen, e.g., if chosen such cube is completely deformed and rather resembles a ring, that ϕ ∈< 0◦ , 90◦ >, then only points above the x-y plane as seen in Figures 5 (b) in the dark blue areas of the image. (with a positive z value) are used. The grayscale image However, this disadvantage is of little significance for our can be directly used as input to the age estimation model. purposes, since all pubic scans have no bottom (3D scan Moreover, to increase the variability of the input training captures only the surface of the bone, not internal parts of dataset and the robustness of our model, we have gener- the bone), and the coordinate system is located in such a ated 41 projections (grayscale images) for each 3D scan way that the most important areas of the scan are deformed with a slightly rotated and translated origin of the Carte- only slightly. Figure 6 shows the distribution of a set of sian coordinate system. 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 com- plicated 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 se- Figure 7: Symphyseal surface fom Figure 2 in the pro- lect the maximum, minimum, median or the average of all posed coordinate system with the resolution of 2◦ of ϕ points mapped to the same (λ , µ) point. In this situation, nected feedforward networks that transform 1024 features into a single real value. For better generalisation, we used a dropout layer in the second part. image NetChain  Input array ( size: 1 ×160 ×160 ) Augm ImageAugmentationLayer array ( size: 1 ×128×128) Model NetChain ( 20 nodes ) vector ( size: 1 ) Output scalar Model: NetChain Input array ( size: 1 ×128×128) 1 ConvolutionLayer array ( size: 8 ×122 ×122 ) 2 BatchNormalizationLayer array ( size: 8 ×122 ×122 ) Figure 8: Generated grayscale image 3 4 Ramp ConvolutionLayer array ( size: 8 ×122 ×122 ) array ( size: 16×116 ×116 ) 5 BatchNormalizationLayer array ( size: 16×116 ×116 ) 6 Ramp array ( size: 16×116 ×116 ) 7 PoolingLayer array ( size: 16×58 ×58 ) 4 Age estimation model 8 ConvolutionLayer array ( size: 16×56 ×56 ) 9 BatchNormalizationLayer array ( size: 16×56 ×56 ) 10 Ramp array ( size: 16×56 ×56 ) 11 PoolingLayer array ( size: 16×28 ×28 ) Our age estimation model consists of several identical age 12 ConvolutionLayer array ( size: 16×26 ×26 ) predictors. Each predictor is based on convolutional neu- 13 BatchNormalizationLayer array ( size: 16×26 ×26 ) 14 Ramp array ( size: 16×26 ×26 ) ral network [16]. Figure 9 shows the main idea of the age 15 PoolingLayer array ( size: 16×8 ×8 ) 16 FlattenLayer vector ( size: 1024 ) estimation flow. First, the 3D scan is transformed into sev- 17 DropoutLayer vector ( size: 1024 ) 18 LinearLayer vector ( size: 100 ) eral grayscale images (see Figure 8 as an example). These 19 Ramp vector ( size: 100 ) images are consecutively directly used as input for indi- 20 LinearLayer vector ( size: 1 ) Output vector ( size: 1 ) viduals age predictors. Second, an aggregation function is applied in order to combine the results from all predictors, Figure 10: Age predictor. The network is wrapped with and thus, to provide the final prediction. We have chosen an Image Augmentation layer, which implements random mean and median as two possible aggregation functions. transformations during training. 3D scan 41 projections Predictor age1 Aggregation function Predictor age2 4.2 Implementation details Predictor age3 Predicted Age Predictor age4 The model was implemented in Wolfram Mathemat- ica using an in-house neural networks package built Predictor age41 on the MXNet framework. To improve the ro- bustness of the network training, we incorporated an Figure 9: Age prediction for a single individual. There is ImageAugmentationLayer to the input layer. This layer a single 3D scan for which we build multiple projections takes the input image 160x160 pixels and randomly crops (41 in this case). By applying the predictor for each of it to 128x128 pixels during the training phase – this al- the projections, we obtain multiple age predictions that are lows us to efficiently expand the training dataset without finally aggregated to gain the final predicted age. having to implement a custom batch function. During the evaluation phase, the ImageAugmentationLayer crops the input image around the centre in a deterministic way, 4.1 Predictor structure so it does not affect it during evaluation. We have chosen slightly larger kernels (7x7) in the first two layers to bet- The application of convolutional neural networks was an ter handle larger structures in blurred images. We use the obvious choice. We have experimented with various ar- batch normalisation layers which are proposed as a tech- chitectures. In the final experiments, we were mainly in- nique to help coordinate the updating of multiple layers spired by the setup used for X-Ray images processing [17]. in the model [16]. Figures 11 and 12 are used as exam- We also experimented with topologies based on DenseNet ples illustrating the extracted patterns for selected layers of [18, 19], which exploits a specific topology that shortens the network for 20- and 72-years-old individuals, respec- layer connections by connecting each layer to every other tively. These figures show the input image and the output layer in a feed-forward fashion. The final predictor struc- from layers #3, #6, #10 and #19. As it can be seen, the ture is shown in Figure 10. The model consists of a total model identifies the vertical structures and the edge of the of 20 layers. In the first part, the input image is reduced symphyseal surface reasonably well from the input image. and transformed into features using convolutional layers These vertical structures combined with the shape of the in combination with pooling, activation (Ramp) and regu- symphyseal edge are also used by experts to identify the larisation layers. The second part represents densely con- age of the individual. Figure 11: An example of a 20-years-old individual (predicted age = 20.82) evaluation – the surface structure can be clearly recognised. Figure 12: An example of a 72-years-old individual (predicted age = 71.6) evaluation. Even though we can see almost no details in the input image, the model can identify vertical (in this orientation) structures, which seems to be a key for the age identification. 4.3 Training have performed many experiments from 100 to many thou- sands), representing almost 10 million processed records. As already mentioned, we have multiple images for a sin- gle 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 4.4 Evaluation to process a single image on input, so the training set is an unstructured list of pairs {image, age}. Since we need To compare the models and study the behaviour, we have proper testing, all images for a single individual must be in performed two ways of model evaluation. Obviously, the the same fold. Cross-validation uses the information about main goal is to predict the age of an (unknown) individual the ID of the individual to split the dataset properly. The based on the 3D scan of the bone, so we evaluated the folds are then flattened, and we perform standard training model for each individual (41 or 82 images) and computed for mapping images to real values (age). the predicted age using aggregation. Figure 13 shows the We use a standard Adam optimiser [20] with a batch actual age vs. the predicted age per individual. Based on size 16 (experimentally chosen) running on the GPU for the size of the dataset, we have chosen 5-fold validation training. Each training runs approximately 300 rounds (we for all of our experiments. Aggregation: Mean, RMSE = 12.92, 5-fold cross validation 90 be seen, our model generally overestimates younger in- dividuals (under the age of 55 years) and underestimates 80 mature ones (above the age of 55 years). This is a result of the tendency to predict the age towards the mean age of 70 the sample. Predicted age 60 90 90  50 80 80 40  70 70  Predicted age 30 60 60  20 50 50 20 30 40 50 60 70 80 90   Actual age 40  40 Figure 13: Age prediction with mean aggregation func- 30 30 tion. Each dot represents a single individual. 20 20 5 Results and discussion 20 30 40 50 60 70 80 90 Actual age Stoyanova et al. presented in their study [13] five Figure 14: Variation of age predictions for particular age age-estimation models with similar age estimation per- class (all individuals mixed) for our model. The image formance and provided an open source software called shows the minimum, maximum, q1 , q2 (median) and q3 forAGE (available at http://morphlab.sc.fsu.edu/). The quartiles. Black dots represent outliers defined by quartiles best model (according to their study) is "SAH&VC" and 1.5× interquartile range. The central line connects the (SAH+Outline) and provides MAE of approximately median values. Aggregation function: Mean 10.79 for their entire dataset (93 samples). To compare our model with their state-of-the-art model, we used their Similarly, we analysed the SAH&VC model from [13]. software and evaluated the MEA for our dataset as well. Figure 15 shows the variation of age predictions for partic- The results are summarised in Table 3. It should be noted ular age class. As can be seen, the SAH&VC model gen- here that for our models, MAE and RMSE are computed erally underestimates all individuals above 40-50 years. using 5-fold cross-validation, whereas for the model from More specifically, for individuals over 45 years old, the [13] MAE and RMSE are computed directly. average of all estimations reaches only 37.7 years (for our dataset). We believe that this primarily results from the Table 3: Age prediction results. The table compares MAE unbalanced age distribution of the dataset used in [13]. and RMSE of our models with the model designed by Stoyanova et al. [13] used on our dataset. 90 90 Age estimation model MAE RMSE 80 80 Our model (Median) 10.63 12.94 Our model (Mean) 10.60 12.92 70 70 SAH&VC 19.52 25.04  Predicted age 60 60  As the presented results indicate, our models outper- 50 50 form the model developed by Stoyanova et al. [13] in  terms of both MAE and RMSE. There is a significant dis- 40 40 crepancy between the MAE presented by Stoyanova et 30 30 al. and the MAE computed on our dataset (i.e., 10.79 vs.  19.52 years). This discrepancy is discussed below in the 20 20 text. To determine whether our model contains any system- 20 30 40 50 60 70 80 90 atic error or whether any particular age intervals introduce Actual age some anomalies in the prediction, we processed the pre- dicted ages per one-year age intervals. Figure 14 shows Figure 15: Variation of age predictions for particular age the variation in age predictions for each age class. As can class according to SAH&VC model by Stoyanova [13]. Moreover, we observed that our model can estimate the [9] Milner, G. R., Wood, J. W., & Boldsen, J. L. (2018). Paleode- age-at-death of an individual over the entire age interval mography: problems, progress, and potential. In: Biological (in our case between 19–92 years) – see Figure 14. This Anthropology of the Human Skeleton (Third, pp. 593–633). contrasts with e.g. [2] and [21], where pubic symphysis New York: John Wiley & Sons. is considered appropriate for individuals up to 40 years, [10] Nikita, E. (2017). Osteoarchaeology: A guide to the macro- or 60 years, respectively. When the maturation process of scopic study of human skeletal remains (First). London: pubic symphysis is complete, the morphological changes Academic Press. are degenerative and highly variable between individuals [11] Schmitt, A., Murail, P., Cunha, E., & Rougé, D. (2002). [2], [22]. However, we believe that our model can capture Variability of the pattern of aging on the human skeleton: Evidence from bone indicators and implication on age at even such changes. death estimation. J Forensic Sci, 47, 1203–1209. [12] Slice, D. E., & Algee-Hewitt, B. F. B. (2015). Modeling 6 Conclusion Bone Surface Morphology: A Fully Quantitative Method for Age-at-Death Estimation Using the Pubic Symphysis. Jour- We have developed a novel age-at-death estimation model nal of Forensic Sciences, 60, 835–843. based on convolution neural networks. Our model pro- [13] Stoyanova, Algee-Hewitt, B. F. B., Kim, J., & Slice, D. E. vides a mean absolute error of approximately 10.6 years (2017). A Computational Framework for Age-at-Death Esti- and is suitable for adult and mature individuals. Our re- mation from the Skeleton: Surface and Outline Analysis of sults indicate that the pubic symphysis reflects the age of 3D Laser Scans of the Adult Pubic Symphysis. J Forensic an individual throughout their entire adult life. In other Sci, 62, 1434–1444. words, we have observed no limitations in terms of age [14] Stoyanova, Algee-Hewitt, B. F. B., & Slice, D. E. (2015). prediction capabilities of pubic symphysis of adult indi- An enhanced computational method for age-at-death estima- viduals. tion based on the pubic symphysis using 3D laser scans and thin plate splines. American Journal of Physical Anthropol- ogy, 158, 431–440. Acknowledgments This research was supported by a re- search grant awarded by the Technology Agency of the Czech [15] Villa, C., Gaudio, D., Cattaneo, C., et al (2015). Surface Republic; project number TL03000646. Curvature of Pelvic Joints from Three Laser Scanners: Sep- arating Anatomy from Measurement Error. Journal of Foren- sic Sciences, 60, 374–381. References [16] Ian Goodfellow, Yoshua Bengio, and Aaron Courville. Deep Learning. MIT Press, 2016. [1] Adserias-Garriga, & Wilson-Taylor, R. (2019). Skeletal age [17] Christie Natashia Archie. Chest X-rays Pneumo- estimation in adults. In J. Adserias-Garriga (Ed.), Age Es- nia Detection using Convolutional Neural Network, timation: A Multidisciplinary Approach (First, pp. 55–73). 2020, https://towardsdatascience.com/chest-x-rays- London: Academic Press. pneumonia-detection-using-convolutional-neural-network- [2] Baccino, Sinfield, L., Colomb, S., Baum, T. P., & Martrille, 63d6ec2d1dee L. (2014). Technical note: The two step procedure (TSP) for [18] Huang, Gao, Zhuang Liu, and Kilian Q. Wein- the determination of age at death of adult human remains in berger. Densely Connected Convolutional Networks. CoRR forensic cases. Forensic Science Int., 244, 247–251. abs/1608.06993 (2016). http://arxiv.org/abs/1608.06993. [3] Biwasaka, H., Sato, K., Aoki, Y., et al. (2013). Three dimen- [19] Sultana, Farhana, Abu Sufian, and Paramartha Dutta. sional surface analyses of pubic symphyseal faces of con- Advancements in Image Classification Using Convolu- temporary Japanese reconstructed with 3D digitized scanner. tional Neural Network. CoRR abs/1905.03288 (2019). Legal Medicine, 15, 264–268. http://arxiv.org/abs/1905.03288. [4] Brooks, S., & Suchey, J. M. (1990). Skeletal age determi- [20] Diederik P. Kingma and Jimmy Ba. Adam: A Method for nation based on the os pubis: a comparison of the Acsádi- Stochastic Optimization. In: 3rd International Conference Nemeskéri and Suchey-Brooks methods. Human Evolution, for Learning Representations, San Diego, 2015. Revised 30 5, 227–238. Jan 2017. [5] Buk, Z., Kordík, P., Brůžek, J., Schmitt, A., & Šnorek, M. [21] Teixeira, F., Cunha, E. (2021). Aging the elderly: Does the (2012). The age at death assessment in a multi-ethnic sample skull tell us something about age at death? In B. F. B. Algee- of pelvic bones using nature-inspired data mining methods. Hewitt and J. Kim (Eds.), Remodeling Forensic Skeletal Age Forensic Science International, 220, 294.e1-e9. (pp. 75–97). Academic Press. [6] Dudzik, B., & Langley, N. R. (2015). Estimating age from [22] Meindl, R.S., Lovejoy, C.O., Mensforth, R.P., WalkerR.A., the pubic symphysis: A new component-based system. A revised method of age determination using the os pubis, Forensic Science International, 257, 98–105. with a review and tests of accuracy of other current methods [7] Hartnett, K. M. (2010). Analysis of age-at-death estimation of pubic symphyseal aging, Am. J. Phys. Anthropol. 68 (Sep using data from a new, modern autopsy sample - Part I: Pubic (1)) (1985) 29–45. bone. Journal of Forensic Sciences, 55, 1145–1151. [8] Langley, & Tersigni-Tarrant, M. A. (2017). Forensic Anthro- pology: A Comprehensive Introduction (2nd edn.). Boca Ra- ton: CRC Press.