Computerized assembly of neurocranial fragments based on surface extrapolation

Abstract

Fossil crania are often fractured and fragmented due to compaction and diagenesis. To restore the antemortem appearance of a fossil cranium, it is necessary to correctly assemble the fragments into their original anatomical positions. In this study, we propose a concept for computerized reconstruction that employs surface extrapolation to aid the assembly of fossil neurocranial fragments. Specifically, we approximate the surface of each neurocranial fragment using a bicubic Bézier surface to extrapolate the surface and mathematically predict the shape of adjacent fragments. The positions and orientations of adjacent fragments were calculated by minimizing the fitting errors. To evaluate the usefulness of this concept, we virtually divided modern human and chimpanzee neurocrania into pieces and used the proposed method to reassemble the generated virtual fragments. The neurocranial fragments were smoothly and correctly assembled. Comparison of the results obtained using the proposed method and conventional manual assembly revealed that the proposed method delivered similar performance in terms of differences between the original and reassembled shapes. However, the accuracy of the reassembly was found to be worse in the chimpanzee case because the fragments were more curved than those for the human cranium. Although there are some methodological limitations, the proposed concept may be useful for development of digital reassembly of fossil neurocranial fragments.

Introduction

Fossil crania are often fractured and fragmented due to compaction and diagenesis. Consequently, they are usually recovered in multiple scattered pieces. To restore the ante-mortem appearance of a fossil cranium for morphological investigations, the fragments must be assembled like a jigsaw puzzle based on shape information. Conventionally, they are reconstructed manually based on the knowledge and experience of skilled anthropologists. However, manual reconstruction may not yield reproducible results since it is susceptible to subjective considerations. To realize more objective, precise, and reproducible reconstruction of cranial fossils, it is essential to develop digital reconstruction methods to aid reassembly of fossil fragments.

Recently, several cranial fossils have been reconstructed in a virtual environment. Cranial fossil fragments were first virtually reconstructed in 1995 when Kalvin et al. (1995) restored a mid-Pleistocene fossil cranium excavated in Morocco, and Zollikofer et al. (1995) reconstructed the cranium of a Neanderthal child from Devil’s Tower. Digital assembly of fossil crania has subsequently become widespread in the field of physical anthropology and has produced significant results and benefits (e.g. Ponce de León and Zollikofer, 1999). Recently discovered crania of the early hominids Sahelanthropus tchadensis and Ardipithecus ramidus have also been digitally restored (Zollikofer et al., 2005; Suwa et al., 2009). The importance of digital restoration is growing. However, current digital fossil restoration techniques essentially only facilitate manual assembly of fragments into anatomical positions by allowing each piece to be positioned and rotated using a computer mouse.

In the field of computational geometry, automated reassembly of fractured objects has been investigated, since this is essentially viewed as a pattern-matching problem. For example, methods of reassembling objects from fragments based on the geometric similarity of the fracture surfaces have been proposed (e.g. Papaioannou et al., 2002; Huang et al., 2006). Recently, Chowdhury et al. (2009) and Chowdhury and Bhandarkar (2011) attempted to reconstruct a highly fractured mandible based on the automatic registration of manually identified opposing fracture surfaces on the broken bone fragments. Yu et al. (2012) proposed an automatic reassembly method for cranial fragments that used a cranium template to guide the reassembly process. The heat kernel signature was used to perform the initial matching between the template and fragments, and the fragments were then repositioned based on fracture curve matching. However, as the fractured surfaces of fossil neurocranial fragments are typically thin and damaged by erosion, and some fragments of a fossil cranium are usually missing, such registration techniques based on the correspondence of fracture surfaces are not generally applicable. Willis and Cooper (2008) attempted to computationally reassemble fragments of ceramic pots by exploiting their axial symmetry. However, such simple geometrical constraints cannot usually be applied to human crania, and automatic reassembly of fossil skull fragments using a computer is thus a challenging problem.

In the present study, we propose a concept for computerized reconstruction based on surface extrapolation that may aid in the assembly of fossil neurocranial fragments. Specifically, we approximate the surface of each neurocranial fragment using a parametric Bézier surface to extrapolate the fragment and mathematically predict the surface shape of adjacent fragments. The positions and orientations of the fragments are then calculated by minimizing the fitting errors, making it possible to align the fragments continuously and smoothly. To evaluate the usefulness of this concept, we virtually divided modern human and chimpanzee crania into pieces and reassembled the generated virtual fragments using the proposed concept.

Methods

Representation and extrapolation of fragment surfaces

We used a bicubic Bézier surface to approximate the surface of a fragment. The Bézier surface can be described as a polynomial function of parameters u and v:

  

$$\begin{array}{l}\mathbf{p}(u,v)={[p_x\hspace{0.17em}{p}_y\hspace{0.17em}{p}_z]}^T=\sum _{i=0}^3\sum _{j=0}^3{\mathbf{c}}_{ij}\hspace{0.17em}{B}_i(u)B_j(v)\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }u,v\in [0,1]\\ B_i(u)=\frac{6}{i!(3-i)!}{u}^i{(1-u)}^{3-i},\end{array}$$(1)

where c_ij are the position vectors of the 16 control points that define the surface shape.

Let q_k, k = 1, …, K, be points on the surface and let both p and q be defined in the fragment-fixed coordinate system whose coordinate axes are the principal axes of the fragment and whose origin is the center of the axis-aligned minimum bounding box of the fragment. An approximation of the surface by a bicubic Bézier surface can be computed by solving the following least-squares minimization problem (Figure 1):

Figure 1

Representation and extrapolation of fragment surface. A bicubic Bézier surface is fitted to the surface of a fragment. The rectangular box indicates the axis-aligned minimum bounding box of the fragment.

  

$$E(\mathbf{c})=\sum _{k=1}^K{\{\mathbf{p}(u_k,v_k)-{\mathbf{q}}_k\}}^2+\lambda \iint \left\{{\left(\frac{{\partial }^2\mathbf{p}}{\partial u^2}\right)}^2+2{\left(\frac{{\partial }^2\mathbf{p}}{\partial u\partial v}\right)}^2+{\left(\frac{{\partial }^2\mathbf{p}}{\partial v^2}\right)}^2\right\}dudv\to \text{min},$$(2)

where E is the objective function, p(u_k, v_k) is the point on the Bézier surface closest to the kth point q_k, and λ is the weighting coefficient. The first term represents the fitting error between the points generated by the function and the corresponding points on the surface. The second term is the smoothness penalty function that represents the bending energy of the surface.

It is usually difficult to determine the correspondence between p and q in order to calculate the first term. However, cranial fragments are generally flat in the direction of the z axis, and hence the parametric coordinates (u_k, v_k) corresponding to q_k can be calculated by solving the following simultaneous equations:

  

$$\left[\begin{array}{c}{p}_x(u_k,v_k)\\ p_y(u_k,v_k)\end{array}\right]=\left[\begin{array}{c}{q}_{x,k}\\ q_{y,k}\end{array}\right]$$(3)

where q_k = [q_x,k, q_y,k, q_z,k]T. Therefore, the first term in equation (2) can be rewritten as:

  

$$E(c_z)=\sum _{k=1}^K{\left\{p_z(u_k,v_k)-q_{z,k}\right\}}^2+\lambda \iint \left\{{\left(\frac{{\partial }^2\mathbf{p}}{\partial u^2}\right)}^2+2{\left(\frac{{\partial }^2\mathbf{p}}{\partial u\partial v}\right)}^2+{\left(\frac{{\partial }^2\mathbf{p}}{\partial v^2}\right)}^2\right\}dudv\to \text{min},$$(4)

To predict the shape of the outer periphery of the fragment, the fragment must be approximated by the center of the domain of the two parameters u and v. For this reason, the x and y coordinates of the 16 control points (c_x,ij, c_y,ij) were predetermined as follows:

  

$$\left[\begin{array}{c}{c}_{x,ij}\\ c_{y,ij}\end{array}\right]=\left[\begin{array}{c}\mathrm{\Delta }x(-0.5+m_i)/0.4\\ \mathrm{\Delta }y(-0.5+m_j)/0.4\end{array}\right],$$(5)

where Δx and Δy are the side lengths of the axis-aligned minimum bounding box of the fragment and (m₁, m₂, m₃, m₄) = (0.0, 0.3, 0.7, 1.0). This equation implies that the x and y coordinates of the four interior control points, c₁₁, c₁₂, c₂₁, and c₂₂, are set to the vertices of the minimum boundary box so that the Bézier surface can be extrapolated to the fragment surface within the domain of the parameters. If the x and y coordinates of the 16 control points are predetermined, the parametric coordinates (u, v) corresponding to the control points c_ij can be calculated by solving the following simultaneous equations, which are analogous to equation (3):

  

$$\left[\begin{array}{c}{p}_x(u_{ij},v_{ij})\\ p_y(u_{ij},v_{ij})\end{array}\right]=\left[\begin{array}{c}{c}_{x,ij}\\ c_{y,ij}\end{array}\right]$$(6)

Hence, the parameters sought are the z coordinates of the 16 control points. An accurately and smoothly fitted surface can be obtained by finding the values of the 16 control points that minimize E. We used the quasi-Newton method to solve this minimization problem numerically. The nonlinear simultaneous equations (3) and (6) were also numerically solved by the Newton–Raphson method. As the prediction tends to get worse at either end of the domain of the parameters u and v, the Bézier surface is defined only when the parameters satisfy the following condition:

  

$$\sqrt{{(u-0.5)}^2+{(v-0.5)}^2}\le \mathrm{0.45.}$$(7)

This inequality implies that the fitted Bézier surface is elliptically rather than rectangularly distributed.

The approximation accuracy would be better if the smoothness penalty function (i.e. the second term in equation (4)) did not exist and a higher-order Bézier surface function was used. However, the approximation would then be more sensitive to noise and relatively unimportant details, resulting in overfitting, as has been previously suggested (Deane et al., 2005; Matsuura et al., 2010). Minimizing both the fitting error and the smoothness penalty is therefore essential for accurately capturing the global characteristics of each cranial fragment and for accurately extrapolating the fragment to predict the shape of adjacent fragments.

Pairwise assembly of fragments

Let A and B be two adjacent fragments. In the present study, the orientation and translation of B with respect to A is calculated by minimizing the following objective function (Figure 2):

Figure 2

Pairwise assembly of fragments. The fitting errors between the extrapolated Bézier surfaces of fragments A (^Ap) and the points on the surfaces of fragments B (^Bq) and vice versa are mutually minimized to calculate the position and orientation of fragment B with respect to fragment A. The distances between the manually defined corresponding points on the edges of A and B (^Ar and ^Br) are also simultaneously minimized for assembly.

  

$$\begin{array}{l}{E}_{AB}({\theta }_x,{\theta }_y,{\theta }_z,{}^A\mathbf{t})=\delta \sum _k{\{{}^A\mathbf{p}(u_k,v_k)-(\mathbf{R}·{{}^B\mathbf{q}}_k+{}^A\mathbf{t})\}}^2\\ +\varepsilon \sum _k{\{\mathbf{R}·{}^B\mathbf{p}(u_k,v_k)-({{}^A\mathbf{q}}_k-{}^A\mathbf{t})\}}^2\\ +\mu \sum _k{\{{{}^A\mathbf{r}}_k-(\mathbf{R}·{{}^B\mathbf{r}}_k+{}^A\mathbf{t})\}}^2+\kappa I_{AB},\end{array}$$(8)

where ^Ap and ^Bp are respectively the approximated Bézier surfaces of fragments A and B (the superscripts also refer to the coordinate system), R(θ_x, θ_y, θ_z) is a rotation matrix representing the orientation of B with respect to the A coordinate system, ^At is a translation vector connecting the origins of the A and B coordinate systems, ^Aq and ^Bq are respectively the points on the surfaces of fragments A and B, ^Ar and ^Br are respectively the manually defined corresponding points on the edges of A and B, I is the penetration penalty function, and δ, ɛ, μ, and κ are weighting coefficients. The first and second terms on the right-hand side respectively represent the mean square difference between ^Ap and ^Bq and between ^Bp and ^Aq, and the third term represents the mean square difference between the corresponding points. Therefore, B is spatially translated and rotated to align with A such as to minimize the fitting error between the extrapolated Bézier surfaces and the corresponding points on the surfaces. The weighting coefficients, δ and ɛ, are appropriately determined based on the number of data points in each calculation step. The penetration penalty function I is calculated by a collision detection library (SmartCollisionSDK, v. 2.3, I-NET Corporation, Japan), which can efficiently detect collisions and penetration depths of three-dimensional mesh models. The assembly results are affected by the values of the weighting coefficients in equation (8). Here, we systematically vary two coefficients, μ and κ, to evaluate how they affect the assembly accuracy; we determined them to be 0.3 and 500, respectively, so that the third term is balanced by the sum of the first and second terms, and the fourth becomes larger than the other three terms when the fragments penetrate each other.

Multi-piece assembly

The first three neighboring fragments, A, B, and C, were assembled by minimizing the following objective function:

  

$$E_{AB}+E_{BC}+E_{AC}\to \text{min}.$$(9)

Here, 12 parameters are being sought since six degrees of freedom of one of the three fragments are fixed by the global coordinates. From this calculation, the positions and orientations of the three fragments with respect to each other were fixed. The fourth fragment D, which is supposedly adjacent to B and C here, was assembled by minimizing the following function:

  

$$E_{BD}+E_{CD}\to \text{min}.$$(10)

In this case, six parameters are sought, namely the translation and orientation of D alone, as the positions and orientations of B and C have already been fixed. In other words, the relative position and orientation of each fragment were determined locally by the two neighboring fragments. By repeating this process for each fragment, multiple fragments were assembled.

Ideally, the position and orientation of all fragments should be simultaneously calculated. However, since this is computationally very expensive because of the large search space, the fragments are assembled one by one.

Evaluation of the proposed concept

To evaluate the usefulness of the proposed concept, virtual fragments were generated from human and chimpanzee cranial models, and then reassembled using the above method using software written by the authors. Because the original shapes of the virtually fragmented crania were known, the accuracy of the reassembly process could be quantitatively estimated. Three-dimensional human and chimpanzee cranial models were generated based on computed tomography of a modern Japanese human specimen housed at Kyoto University and a female chimpanzee specimen housed at the Japan Monkey Centre, respectively. The models were virtually divided into pieces, as shown in Figure 3. Specifically, the frontal and parietal portions of the crania were virtually clipped and each portion was manually divided into pieces. The edges of the fragments were then abraded to a depth of approximately 1 mm in order to mimic the fracture surfaces of actual fossil fragments. These virtual fragments were then used to test how well a Bézier surface could approximate the surface of neurocranial fragments and predict the surfaces of surrounding fragments. The four sets of virtual fragments were reassembled based on the proposed method, and the accuracy of the reconstruction was compared with that for manual reassembly by two postgraduate students majoring in mechanical engineering (biomechanics).

Figure 3

Virtual cranial fragments generated for evaluation of the reassembly method. (A) Frontal portion of a human cranium. (B) Parietal portion of a human cranium. (C) Frontal portion of a chimpanzee cranium. (B) Parietal portion of a chimpanzee cranium.

Results

Figure 4 shows approximated Bézier surfaces for the virtually generated fragments of the human cranium. The residual errors between the fragment surface and the Bézier surface were 0.03 ± 0.02 mm and 0.05 ± 0.04 mm (mean ± standard deviation) for the fragments in the frontal and parietal portions, respectively. For the chimpanzee fragments, the corresponding residual errors were 0.11 ± 0.09 mm and 0.12 ± 0.10 mm, respectively. Although the approximation errors were higher for the chimpanzee fragments, the proposed method can still approximate the surfaces of the cranial fragments with an accuracy of about 0.2 mm.

Figure 4

Approximated Bézier surfaces for the virtually generated neurocranial fragments. (A) Fragments of a human frontal portion. (B) Fragments of a human parietal portion. (C) Fragments of a chimpanzee frontal portion. (D) Fragments of a chimpanzee parietal portion.

We also calculated how well the Bézier surfaces could extrapolate the fragments to predict the surrounding surfaces. The average residual errors were 0.53 ± 0.28 mm and 0.81 ± 0.59 mm for the fragments in the frontal and parietal portions of the human cranium, respectively. This indicates that the surfaces around a fragment could be predicted with the accuracy of ~0.5–1 mm in the human case. In the chimpanzee case, however, the corresponding average residual errors were 1.52 ± 1.13 mm and 1.92 ± 1.38 mm, respectively, representing a prediction accuracy of ~2 mm.

Figure 5 shows the results of reassembling two human and two chimpanzee cranial portions using the proposed method. As shown in Figure 5, the fragments were assembled with smooth junctures. The results for manual assembly are also shown in the figure. The contours indicate the distributions of errors between the assembled fragments and the original shape of the cranium. For the proposed computerized method, the residual errors between the reassembled fragments and the original shape for the human neurocranium (Table 1) were 0.30 ± 0.30 mm and 0.30 ± 0.28 mm in the frontal and parietal portions, respectively, and those for the chimpanzee neurocranium were 0.60 ± 0.52 mm and 0.68 ± 0.62 mm, respectively. On the other hand, the corresponding errors for manual assembly by the first participant were 0.29 ± 0.27 mm and 0.50 ± 0.40 mm for the human neurocranium, and 0.17 ± 0.17 mm and 0.40 ± 0.39 mm for the chimpanzee neurocranium. For the second participant, the corresponding errors were 0.53 ± 0.43 mm and 0.65 ± 0.67 mm for the human neurocranium, and 0.87 ± 0.63 mm and 0.82 ± 0.68 mm for the chimpanzee neurocranium (Table 1). Thus, similar levels of accuracy were obtained using the automatic and manual reassembly methods.

Figure 5

Comparison of reassembly results. Contours indicate the distribution of errors between the assembled fragments and the original shape of the cranium. (A) Results of the computerized reassembly. (B) Results of manual reassembly by two postgraduate students.

Table 1 Comparison of residual errors (mean ± standard deviation) between proposed computerized reassembly and manual reassembly by two postgraduate students

		Computerized method [mm]	Manual assembly [mm]
			A	B
			Human	frontal	0.30 ± 0.30	0.29 ± 0.27	0.17 ± 0.17
	parietal	0.30 ± 0.28	0.50 ± 0.40	0.40 ± 0.39
Chimpanzee	frontal	0.60 ± 0.52	0.53 ± 0.43	0.87 ± 0.63
	parietal	0.68 ± 0.62	0.65 ± 0.67	0.82 ± 0.68

Discussion

The present study proposed a concept for computerized reassembly of fossil cranial fragments based on the tangent continuity (smoothness) of the joints between fragments. The surfaces of cranial fragments were approximated using Bézier surfaces to mathematically predict the anatomical positions and orientations of adjacent fragments, and the fragments were successfully assembled based on these predictions.

The results indicate that a bicubic Bézier surface can successfully approximate the overall surface morphology of human cranial fragments. The matching error for each point on the virtual human neurocranial fragments was less than 0.1 mm, demonstrating that the approximation is very accurate. However, the matching error for the virtual chimpanzee fragments was found to be larger, possibly due to the fact that the chimpanzee cranium is much smaller and the fragments are more curved. Therefore, although a bicubic Bezier surface may be suitable for approximating human neurocranial fragments, a higher-order parameterization may be necessary for fragments with higher curvature.

Reflecting this difference in the approximation accuracy for human and chimpanzee fragments, the prediction error for each point on adjacent fragments was also found to be much smaller (about ~0.5–1 mm) in the human case than in the chimpanzee case (~2 mm). For reassembling human cranial fragments, the proposed estimation method based on extrapolation is reasonably accurate and the calculated Bézier surface can be used as a template. However, these results again indicate the need for a higher-order surface function for the chimpanzee. The relationship between the prediction error at the outer periphery of fragments and the order of the surface function should be investigated in the future, so that an appropriate surface function can be chosen depending on the size and curvature of fragments.

If the prediction for the outer periphery is reasonably accurate, the present method can align the fragments even if their outer boundaries are eroded or broken off. However, the error is expected to increase if a fragment is relatively large or is highly elongated, because a larger outer periphery is being estimated using a Bézier surface based on a relatively small fragment surface. In such a case, the area of the extrapolated outer periphery should be reduced by modifying equation (5).

As shown in Figure 5, a similar predictive accuracy was obtained using computerized and manual assembly of the virtual fragments, which suggests that the proposed concept based on surface smoothness and continuity can be useful for automatic assembly of human cranial fragments. It is interesting to note that the residual errors are smaller than the prediction errors at the outer periphery of fragments. This is possibly due to the fact that the mean square difference between points on the fragments and the extrapolated Bézier surfaces are evaluated bidirectionally, as in equation (8). Such bidirectional evaluation thus seems to be an effective approach for assembling fossil cranial fragments based on smooth continuity.

Manual assembly of cranial fragments is a powerful method for performing virtual reconstruction based on the anatomical knowledge and experience of skilled anthropologists, and state-of-the-art virtual assemblies have been conducted for many well-known fossil crania (e.g. Zollikofer et al., 2005; Gunz et al., 2009; Suwa et al., 2009; Benazzi et al., 2011). However, it is not easy to determine the spatial position and triaxial orientation of any given fragment with respect to other fragments using a two-dimensional computer screen and a mouse. This procedure requires a lot of patience, particularly when there are many cranial fragments to be assembled. Furthermore, it is important to document how each cranial fragment was assembled in order to evaluate the reconstruction results and to ensure their reproducibility. However, this is generally difficult when manual reconstruction is performed. In contrast, computerized reconstruction can assemble cranial fragments while eliminating some of the subjectivity associated with conventional manual assembly. Moreover, the assembly process is reproducible if the assembly order of the fragments is recorded. Although reconstructing fossil crania usually involves dozens of steps depending on the state of preservation, and the proposed method can provide aid with only one of these steps, it may provide useful support for objective assembly of fossil cranial fragments.

The proposed concept does have some limitations. One is that the adjacency relations of the fragments should be provided in advance by defining rough correspondences among the fragments by manually digitizing points on their edges. Thus, this algorithm is not fully automatic. Finding compatible matches and consistent multi-piece adjacency relations by computer is very difficult because the number of pairwise combinations that must be checked becomes prohibitively large when there are many fragments; in addition, the fragments all have quite similar surface shapes. It is thus difficult to find anatomically consistent adjacency relations without any prior information. This task should be performed manually based on the anatomical knowledge and experience of skilled experts. Another limitation is that the assembly results are affected by the order in which the fragments are assembled. Although this situation is not ideal, a suitable order can be determined by performing the assembly several times and comparing the results obtained. Such inspection is necessary for anatomically consistent assembly. The proposed concept is also inapplicable to anatomical areas that are not as smooth as the frontal and parietal surfaces (e.g. cranial base) because the surfaces cannot be represented and extrapolated using a Bézier surface. A final limitation is that if the outer boundary of a fragment is highly eroded, correct placement of adjacent fragments is more difficult because the predicted shape of the outer periphery becomes less reliable as the distance from the boundary increases. Further studies of computerized reconstruction of fossil crania are therefore required in order to realize more precise, objective, and universal protocols for assembly of cranial fragments.

Acknowledgments

We wish to express our sincere gratitude to Takeru Akazawa (Kochi Institute of Technology), Hiromasa Suzuki (University of Tokyo), Osamu Kondo (University of Tokyo), Hajime Ishida (University of the Ryukyus), and Takashi Michikawa (University of Tokyo) for their continuous guidance and support throughout the course of the present study. This study was supported by a Grant-in-Aid for Scientific Research on Innovative Areas “Replacement of Neanderthals by Modern Humans: Testing Evolutionary Models of Learning” from the Japanese Ministry of Education, Culture, Sports, Science and Technology.

References

•  Benazzi  S.,  Bookstein  F.L.,  Strait  D.S., and  Weber  G.W. (2011) A new OH5 reconstruction with an assessment of its uncertainty. Journal of Human Evolution, 61: 75–88.
•  Chowdhury  A.S. and  Bhandarkar  S.M. (2011) Computer Vision-Guided Virtual Craniofacial Surgery: A Graph-Theoretic and Statistical Perspective. Springer, London.
•  Chowdhury  A.S.,  Bhandarkar  S.M.,  Robinson  R.W., and  Yu  J.C. (2009) Virtual multi-fracture craniofacial reconstruction using computer vision and graph matching. Computerized Medical Imaging and Graphics, 33: 333–342.
•  Deane  A.S.,  Kremer  E.P., and  Begun  D.R. (2005) New approach to quantifying anatomical curvatures using high-resolution polynomial curve fitting (HR-PCF). American Journal of Physical Anthropology, 128: 630–638.
•  Gunz  P.,  Mitteroecker  P.,  Neubauer  S.,  Weber  G.W., and  Bookstein  F.L. (2009) Principles for the virtual reconstruction of hominin crania. Journal of Human Evolution, 57: 48–62.
•  Huang  Q.X.,  Floery  S.,  Gelfand  N.,  Hofer  M., and  Pottmann  H. (2006) Reassembling fractured objects by geometric matching. ACM Transactions on Graphics, 25: 569–578.
•  Kalvin  A.D.,  Watson  T.J.,  Dean  D., and  Hublin  J.J. (1995) Reconstruction of human fossils. IEEE Computer Graphics & Applications, 15: 12–15.
•  Matsuura  Y.,  Ogihara  N., and  Nakatsukasa  M. (2010) A method for quantifying articular surface morphology of metacarpals using quadric surface approximation, International Journal of Primatology, 31: 263–274.
•  Papaioannou  G.,  Karabassi  E.A., and  Theoharis  T. (2002) Reconstruction of three-dimensional objects through matching of their parts. IEEE Transactions on Pattern Analysis and Machine Intelligence, 24: 114–124.
•  Ponce de León  M.S. and  Zollikofer  C.P.E. (1999) New evidence from Le Moustier 1: Computer-assisted reconstruction and morphometry of the skull. Anatomical Record, 254: 474–489.
•  Suwa  G.,  Asfaw  B.,  Kono  R.T.,  Kubo  D.,  Lovejoy  C.O., and  White  T.D. (2009) The Ardipithecus ramidus skull and its implications for hominid origins. Science, 326: 68e1–7.
•  Willis  A.R. and  Cooper  D.B. (2008) Computational reconstruction of ancient artifacts. IEEE Signal Processing Magazine, 25: 65–83.
•  Yu  W.,  Li  M., and  Li  X. (2012) Fragmented skull modeling using heat kernels. Graphical Models, 74: 140–151.
•  Zollikofer  C.P.E.,  Ponce de León  M.S.,  Martin  R.D., and  Stucki  P. (1995) Neanderthal computer skulls. Nature, 375: 283–285.
•  Zollikofer  C.P.E.,  Ponce de León  M.S.,  Lieberman  D.E.,  Guy  F.,  Pilbeam  D.,  Likius  A.,  Mackaye  H.T.,  Vignaud  P., and  Brunet  M. (2005) Virtual cranial reconstruction of Sahelanthropus tchadensis. Nature, 434: 755–759.