A Robust and Automatic Method for the Best Symmetry Plane Detection of Craniofacial Skeletons

: The accurate location of the mid-sagittal plane is fundamental for the assessment of craniofacial dysmorphisms and for a proper corrective surgery planning. To date, these elaborations are carried out by skilled operators within speciﬁc software environments. Since the whole procedure is based on the manual selection of speciﬁc landmarks, it is time-consuming, and the results depend on the operators’ professional experience. This work aims to propose a new automatic and landmark-independent technique which is able to extract a reliable mid-sagittal plane from 3D CT images. The algorithm has been designed to perform a robust evaluation, also in the case of large defect areas. The presented method is an upgraded version of a mirroring-and registration technique for the automatic symmetry plane detection of 3D asymmetrically scanned human faces, previously published by the authors. With respect to the published algorithm, the improvements here introduced concern both the objective function formulation and the method used to minimize it. The automatic method here proposed has been veriﬁed in the analysis of real craniofacial skeletons also with large defects, and the results have been compared with other recent technologies.


Introduction
The face is one of the most important means of interaction with other people.A high degree of symmetry between its two halves is positively associated with health, immune competence, and attractiveness [1][2][3][4][5][6][7][8][9][10][11].Therefore, the correct evaluation of the craniofacial skeleton asymmetries is critical from a functional point of view (e.g., an effective occlusion [12,13]), as well as for the aesthetic and social aspects related to the asymmetries themselves [14].It is known that an incorrect restoration in the face symmetry, due to injury, disease or congenital malformation, could have psychological implications that cannot be underestimated.
In recent years, the development of the Medical Imaging, and in particular, Computerized Tomography (CT) and Magnetic Resonance (MRI), has greatly improved the ability to assess the craniofacial symmetry.Nowadays, computer tomography is the gold standard for imaging the skull and the bone structure of the face when planning surgery [15].Due to the complexities of skull shape, the most effective way to analyze CTs is to interact with 3D models.This involves the segmentation of the anatomical regions of interest (e.g., skull) from the cross-sectional 2D images provided by the CT, its 3D volumetric reconstruction by the voxelization and, by using a specific algorithm (as the Marching Cube algorithm), the generation of the surface model in STL format.This 3D discrete model representation is the simplest data format to manipulate complex geometries, but it provides low-level information limited to point coordinates and the normal of triangles enclosing it.In other words, there is no high-level information such as the position of one of the fundamental features of the facial skeleton: the symmetry plane.
To date, in clinical practice, the symmetry plane of the skull, as well as of the facial skeleton, is approximated with the mid-sagittal plane (MSP) [16][17][18][19][20][21][22] that is calculated based on the locations of manually identified landmarks.The corresponding methods, referred to as Cephalometric methods, like all manual ones, are time-consuming, not reproducible and repeatable; in addition, the best choice of the landmarks is still under debate.
In order to improve the accuracy and reliability of the MSP localization, several techniques have been proposed, thereby moving away from the traditional cephalometric approach toward more complex methodologies [23][24][25][26].These alternative approaches are principally based on statistical techniques such as Morphometrics-based methods (e.g., Ordinary Procrustes Analysis (OPA) [24]) or Principal Component Analysis (PCA) followed by Iterative Closest Point (ICP) methods [23,25].Both are mirroring-and-registration methods, in which the midsagittal plane is established after a refined superimposition between the original skull and its mirrored configuration.Morphometric methods reduce, but do not eliminate, the dependence of landmark selection on the final results; in addition, they continue to require the presence of a skilled operator.ICP based methods are landmark-independent but they fail the in presence of strong asymmetries [27].
Di Angelo and Di Stefano [28] have recently proposed a fully-automatic mirroring and registration method based on an original weighted function.Among all the analyzed methods, the results obtained in [28] for the symmetry plane detection for asymmetrically scanned human faces seemed to be the most promising for developing a method for estimation of the symmetry plane of skulls without interactions with the operator.As will be discussed in detail in Section 3, in this paper, firstly, problems in its direct application to the determination of the plane of symmetry of skulls are studied, and then, the substantial improvements to be introduced to overcome the resulting limitations are proposed.
The new method here proposed has been tested for the symmetry plane detection of the human skulls, i.e., both synthetic with artificial asymmetries, and real also with large defects.The results derived from these experiments and the comparison with other recent innovations are discussed critically thereafter.

Related Works
The methods proposed in the literature for symmetry plane detection of the craniofacial skeleton can be classified into the following three categories: -Cephalometric methods -Morphometric methods -ICP-based methods

Cephalometric Methods
In clinical practice, the evaluation of facial asymmetries is based on a cephalometric (2D or 3D) approach, which simplifies the analyses by using a selection of standardized landmarks.A very large number of possible combinations of landmarks have been explored in the literature [21,23,24,29]; some of the most relevant ones are reported in Table 1.Depending on the published methods, the MSP is defined as the plane: passing through a central landmark (e.g., Nasion or Crista Galli, see Table 1) and perpendicular to a horizontal line crossing bilateral skeletal landmarks [30]; -of best fit passing through a selection of midline landmarks [31].
An example of the first type of MSP definition is proposed by Kim et al. [32], whose MSP is defined by Crista Galli (CG) and the left and right Frontozygomatic Sutures (FZS) (Figure 1a,b).The accuracy of this method was demonstrated in [23].Green et al. in [31] proposed the evaluation of MSP as the plane passing through the three central landmarks: Nasion (N), Basion (Ba) and the Incisive Foramen (IF) (Figure 1a,c).According to the authors' results, this choice leads to better results than using the lateral ones.Landmark and MSP definitions according to the methods proposed in [30] and [31].(a) landmark definitions; (b) MSP definition according to [30]; (c) MSP definition according to [31].
Table 1.Landmark definition.As shown qualitatively in Figure 1, the main advantage of using landmarks is that the result is not affected by all skeletal defects, and the set of reference points considered can be changed and adapted as needed.On the other hand, it has been demonstrated that a different set of landmarks could lead to very different MSP [24]; there is still no consensus on the best set of landmarks to consider, and on which is the most accurate cephalometric plane.In any case, some authors highlight  As shown qualitatively in Figure 1, the main advantage of using landmarks is that the result is not affected by all skeletal defects, and the set of reference points considered can be changed and adapted as needed.On the other hand, it has been demonstrated that a different set of landmarks could lead to very different MSP [24]; there is still no consensus on the best set of landmarks to consider, and on which is the most accurate cephalometric plane.In any case, some authors highlight that a symmetry plane defined using manually-identified cephalometric landmarks results in an inaccurate estimation of symmetry, principally for the lack of accuracy and reproducibility of landmark identification and selection strategies [23,25,[32][33][34].

Basion
Although there have been some attempts to automate the selection of landmarks [33], to date, the selection is still completely manual and, as such, the operation is tedious and must be completed by a skilled person.

Morphometric Methods
The morphometric method for defining the plane of symmetry of a craniofacial skeleton has been proposed by Damstra et al. in [24]; results demonstrate that the technique can potentially produce very accurate and reliable MSP.The method exploits the Ordinary Procrustes Analysis (OPA) applied to landmarks manually defined by the user.In particular, the proposed method requires the identification of visible facial anatomical landmarks in the supraorbital and nasal bridge region for the superimposition (SOF, FNM, FZS, and FOM as defined in Table 1 and shown in Figure 2a), and three central landmarks (S, ANS and Pog, as defined in Table 1 and Figure 2a) for the MSP definition.The original surface model and centroid landmarks are mirrored around an arbitrary plane.A shape alignment by means of partial OPA is then applied to attain the superimposition.The centroid markers of the original and mirrored surface models were superimposed using a rigid translation.Then, a rigid rotation of the mirrored landmarks configuration around the geometric midpoint of the superimposed centroid markers is performed until the best fit between all homologous supraorbital and nasal bridge landmarks has been achieved by means of the least squared point distance.The MSP is defined as the plane passing through three points, in Figure 2b  inaccurate estimation of symmetry, principally for the lack of accuracy and reproducibility of landmark identification and selection strategies [23,25,[32][33][34].
Although there have been some attempts to automate the selection of landmarks [33], to date, the selection is still completely manual and, as such, the operation is tedious and must be completed by a skilled person.

Morphometric methods
The morphometric method for defining the plane of symmetry of a craniofacial skeleton has been proposed by Damstra et al. in [24]; results demonstrate that the technique can potentially produce very accurate and reliable MSP.The method exploits the Ordinary Procrustes Analysis (OPA) applied to landmarks manually defined by the user.In particular, the proposed method requires the identification of visible facial anatomical landmarks in the supraorbital and nasal bridge region for the superimposition (SOF, FNM, FZS, and FOM as defined in Table1 and shown in Figure 2a), and three central landmarks (S, ANS and Pog, as defined in Table 1 and Figure 2a) for the MSP definition.The original surface model and centroid landmarks are mirrored around an arbitrary plane.A shape alignment by means of partial OPA is then applied to attain the superimposition.The centroid markers of the original and mirrored surface models were superimposed using a rigid translation.Then, a rigid rotation of the mirrored landmarks configuration around the geometric midpoint of the superimposed centroid markers is performed until the best fit between all homologous supraorbital and nasal bridge landmarks has been achieved by means of the least squared point distance.The MSP is defined as the plane passing through three points, in Figure 2b

ICP-based methods
The ICP-based approach is a mirroring-and-registration method.The mirroring of the original data is carried out with respect to an initial plane of symmetry, which is usually estimated by the PCA [23,25,26,28,35].Then, the source point cloud and the mirrored data are registered by means of an ICP algorithm [36].After alignment, the MPS is calculated by approximating, in the least-squares sense, the middle of the segments joining homologous points in the original and mirrored point cloud.The superimposition could be done using: -complete models [23,35]; -some homologous surface areas selected by the user on the left and right facial skeleton [25] or by specifically designed operators [26].

ICP-Based Methods
The ICP-based approach is a mirroring-and-registration method.The mirroring of the original data is carried out with respect to an initial plane of symmetry, which is usually estimated by the PCA [23,25,26,28,35].Then, the source point cloud and the mirrored data are registered by means of an ICP algorithm [36].After alignment, the MPS is calculated by approximating, in the least-squares Symmetry 2019, 11, 245 5 of 13 sense, the middle of the segments joining homologous points in the original and mirrored point cloud.The superimposition could be done using: complete models [23,35]; -some homologous surface areas selected by the user on the left and right facial skeleton [25] or by specifically designed operators [26].
Although the ICP-based approach can be easily automated, being landmarks-independent, its use for MSP localization in clinical practice is mainly limited by 2 factors: the final solution is strictly affected by the initial one; -the function to be minimized likewise considers all the points of the model, including local defectiveness, such as holes or bony prominences, that should not be considered in the evaluation of symmetry plane of the skull.
Figure 3a,b show, respectively, MSP evaluation by using our implementation of [35] of a perfectly symmetric synthetic skull (Figure 3a) and a synthetic skull with a local asymmetry (Figure 3b).It is evident that these approaches are not able to achieve an acceptable result for the skull in presence of asymmetries, even if the asymmetry is localized.Although the ICP-based approach can be easily automated, being landmarks-independent, its use for MSP localization in clinical practice is mainly limited by 2 factors: -the final solution is strictly affected by the initial one; -the function to be minimized likewise considers all the points of the model, including local defectiveness, such as holes or bony prominences, that should not be considered in the evaluation of symmetry plane of the skull.Figures 3a and 3b show, respectively, MSP evaluation by using our implementation of [35] of a perfectly symmetric synthetic skull (Figure 3a) and a synthetic skull with a local asymmetry (Figure 3b).It is evident that these approaches are not able to achieve an acceptable result for the skull in presence of asymmetries, even if the asymmetry is localized.

The Mirroring and Weighted Approach
The proposed method is an upgraded version of the mirroring and registration technique presented in [28] (in the following referred as the MaWR method: mirroring and weighted registration method), which was proposed for bilateral symmetry estimations of human faces.It is based on an iterative registration algorithm which minimizes an objective function properly designed to filter out any kind of asymmetries.
In the following sections, the published method, its limitations, and the proposed improvements to extend its applicability to skulls are presented.

The published method
The published method [28], whose flow-chart is depicted in Figure 4, starts from a 3D discrete model (PC) and evaluates the symmetry plane by a preliminary first-attempt (Π0), carried out with a PCA algorithm, which is then refined iteratively until its final estimation.
The final estimation Πf of the MSP is obtained by the Levenberg-Marquardt algorithm [37] that iteratively minimizes the following objective function (OFi), whose expression at the i-th step is: where: n is the number of points of the source point cloud (PC);

The Mirroring and Weighted Approach
The proposed method is an upgraded version of the mirroring and registration technique presented in [28] (in the following referred as the MaWR method: mirroring and weighted registration method), which was proposed for bilateral symmetry estimations of human faces.It is based on an iterative registration algorithm which minimizes an objective function properly designed to filter out any kind of asymmetries.
In the following sections, the published method, its limitations, and the proposed improvements to extend its applicability to skulls are presented.

The Published Method
The published method [28], whose flow-chart is depicted in Figure 4, starts from a 3D discrete model (PC) and evaluates the symmetry plane by a preliminary first-attempt (Π 0 ), carried out with a PCA algorithm, which is then refined iteratively until its final estimation.
The final estimation Π f of the MSP is obtained by the Levenberg-Marquardt algorithm [37] that iteratively minimizes the following objective function (OF i ), whose expression at the i-th step is: where: n is the number of points of the source point cloud (PC); -p j is the j-th point belonging to PC; -TS(PC m,i ) is the tessellated surface of the mirrored data at i-th iteration (PC m,i ); -d Hauss (p j ,TS(PC m,i )) is the Hausdorff distance between p j and TS(PC m,i ); -w i,j is the weight associated to p j at i-th iteration.
The Hausdorff distance is calculated between each point p j on the source PC and the tessellated surface TS(PC m,i ), instead of point-point, in order to avoid the asymmetry in surface sampling and to make the distance as independent as possible on points density.The Hausdorff distance is calculated between each point pj on the source PC and the tessellated surface TS(PCm,i), instead of point-point, in order to avoid the asymmetry in surface sampling and to make the distance as independent as possible on points density.( The weights ws,i,j and wr,i,j are both expressed according to the Leclerc function [38], which has its maximum in correspondence with the symmetry plane: where: -di,j is the distance between pj and the symmetry plane Πi (at the i-th iteration).
-σs and σr define, respectively, the distance and the radius values for which the weight is the 36.79% of its maximum value [38].The weight ws,i,j works to reduce the effects of the asymmetries in the acquisition process which are mainly located far from the symmetry plane.The weight wr,i,j works as a filter which excludes from the registration process any local asymmetries, whether they are near or far from the symmetry The weights w i,j play an important role in the functionalities of the method being presented.It is expressed as the product of two specific weights: w i,j = w s,i,j •w r,i,j . ( The weights w s,i,j and w r,i,j are both expressed according to the Leclerc function [38], which has its maximum in correspondence with the symmetry plane: w r,i,j = 1 where: Symmetry 2019, 11, 245 7 of 13 -d i,j is the distance between p j and the symmetry plane Π i (at the i-th iteration).
σ s and σ r define, respectively, the distance and the radius values for which the weight is the 36.79% of its maximum value [38].The weight w s,i,j works to reduce the effects of the asymmetries in the acquisition process which are mainly located far from the symmetry plane.The weight w r,i,j works as a filter which excludes from the registration process any local asymmetries, whether they are near or far from the symmetry plane.
More details for the original algorithm are reported in [28].

Main Limitations and Proposed Improvements
Although the method presented in [28] represents a valid tool for the estimation of the plane of symmetry for human faces acquired by a 3D scanner, it leads to questionable results when applied to a 3D skull model presenting local asymmetry, as shown in Figure 5a.This is mainly due to two different limitations affecting the algorithm: the weight definition w i,j (Equation ( 2)) and the objective function's (Equation ( 1)) minimization strategy.
The weight definition is designed to also avoid the effects of the asymmetries due to the acquisition process.Since for a 3D facial-cranial skull model coming from CT images: the acquisition quality is not affected by the distance from the symmetry plane; -the most symmetrical areas may be those that are farthest from the symmetry plane; -the weight w s,i,j proves to be not only useless but also negatively affecting the results.
For this reason, in the new version of the algorithm, the contribution of w s,i,j has been neglected and the weight w i,j works only as a filter which excludes from the registration process any local asymmetries: Large values of σ r guarantee a robust registration, whereas small values afford an extremely accurate registration of the symmetric parts, excluding all asymmetries.In the experiments described hereafter, since we analyze real skulls with large asymmetries, the value of σ r is assumed to be the 50 percent the maximal width of the skull.
Dealing with the minimization strategy, the Levenberg-Marquardt algorithm used to minimize the objective function proved to be incapable of geting out of the local minimum resulting after the PCA.For this reason, the new algorithm uses the Particle Swarm Optimization [39] which belongs to the class of Meta-heuristic algorithms, whose adoption is widespread in practical applications for the lower computational burden and lack of strict usage hypothesis with respect to exact methods.Although the convergence to a global minimum cannot be guaranteed, we observed that in all the considered scenarios, by using an appropriate research space, a good solution (considering the application here analyzed) was found.
The discussed improvements have been implemented in a new algorithm properly designed for the symmetry plane detection of the craniofacial skeleton.As shown in Figure 5b, the novel algorithm can potentially lead to much better results compared with those obtainable with the original version.lower computational burden and lack of strict usage hypothesis with respect to exact methods.Although the convergence to a global minimum cannot be guaranteed, we observed that in all the considered scenarios, by using an appropriate research space, a good solution (considering the application here analyzed) was found.
The discussed improvements have been implemented in a new algorithm properly designed for the symmetry plane detection of the craniofacial skeleton.As shown in Figure 5b, the novel algorithm can potentially lead to much better results compared with those obtainable with the original version.

Experiments and results
The proposed method (MaWR-method) has been implemented using original software, coded in MATLAB.In order to evaluate the MaWR-method's accuracy in the identification of the plane of symmetry of craniofacial skeletons, selected test cases have been designed and the resulting planes are compared with the ones provided by two Cephalometric methods (Cplm1-method proposed by Kim et al. [30] and the Cplm2-method proposed by Green et al. in [31]) and the Morphometric method (Mpm-method) proposed in [24].These three methods, already presented in sections 2.1 and 2.2, have been selected since, according to [23], [24] and [31], they tend to deliver more reliable results amongst those published in the literature.

Experiments and Results
The proposed method (MaWR-method) has been implemented using original software, coded in MATLAB.In order to evaluate the MaWR-method's accuracy in the identification of the plane of symmetry of craniofacial skeletons, selected test cases have been designed and the resulting planes are compared with the ones provided by two Cephalometric methods (Cplm1-method proposed by Kim et al. [30] and the Cplm2-method proposed by Green et al. in [31]) and the Morphometric method (Mpm-method) proposed in [24].These three methods, already presented in Sections 2.1 and 2.2, have been selected since, according to [23,24,31], they tend to deliver more reliable results amongst those published in the literature.
The following paragraph discusses five of the test cases used for the experimental analysis.The presented test cases are skulls, segmented from anonymous CT, provided by the Meyer Children Hospital of Florence (Italy): -TC#1 and TC#2 are two healthy real skulls (Figure 6a,b), to demonstrate the reliability of each method in real cases.For TC#2 an incomplete skull model has been chosen because, commonly, the TC images acquisition addresses only the region of interest to reduce the irradiation risks.-TC#3, TC#4 and TC#5 are real skulls, each with a large defect (so, large asymmetry) (Figure 6c-e), to demonstrate the reliability and the robustness of the new method compared with the other approaches.These three case studies include one unilateral defect, one bilateral defect and one defect crossing the MSP.
For each skull, a proper segmentation of the hard tissue of the skull has provided by a skilled operator using Materialise Mimics software.Then, the segmented region is exported as a surface model (i.e., in STL format) and analyzed by using the four previously mentioned methods for the symmetry plane estimation.For the cephalometric and the morphometric methods, all the landmarks in Table 1 are manually-selected by a skilled operator using the volume renderings and the cross-sectional slices in all three planes of the CT images.
Since no real skull is perfectly symmetrical and the plane of symmetry does not actually exist, in this paper the comparison among the four methods here considered is carried out using the Asymmetry Value Value (AV): where: The AV is defined as the median of AV point that is the Euclidean distances between each point of the source point cloud (p j ) with respect to its closer tessellated surface triangle (TS(PC m,i )) of the mirrored configuration.The mirrored configuration is obtained reflecting the PC upon the estimated symmetry plane.The AV is calculated as the median instead of the mean in order to reduce the weight of the little strongly asymmetrical regions in the calculation.The distance point-triangle makes it possible to avoid the asymmetries due to the surface model sampling.Preliminarily, it has been verified that for the previously shown synthetic mesh (Figure 1a), the AV value is zero.For real skulls, based on the previously mentioned considerations, the AV value cannot be zero; however, the method for the smaller AV performs the best localization of the symmetry plane.
The MSPs of 5 skull surface models are calculated using the proposed algorithm (MWR-method), a morphometric-based (Mpm-method) and two cephalometric-based methods (Cplm1-method and Cplm2-method), as exposed below.The calculated AV is then used to perform a direct comparison among the different approaches.Table 2 shows the Asymmetry Value related to the different MSPs calculated for each of the considered test cases.For the TC#5, the MSP is not evaluated for Cplm2-method since, being the TC incomplete, the landmark IF is missing.Table 2 shows that the proposed method is able to find the MPS with a comparable (in fact, it is always lower for all the test cases) with respect to the other methods without any user interaction.The presented approach has therefore proven itself to be capable of obtaining a clinically acceptable result in a fully automatic manner.It is worth noting that even small differences in terms of AV correspond to different three-dimensional distributions of asymmetries of the skull.As an example, in Figure 7, the AV point maps for the TC#5 and the four methods here compared are shown.It is evident that the estimation of the symmetry plane carried out with the proposed method determines a better localization of the asymmetries of the skull.

Conclusions
In recent years, the reliable localization of the MSP of the skull has attracted a great deal of research effort, since it is crucial to helping with the clinical diagnosis and treatment planning of craniofacial asymmetry, as well as with the planning of craniofacial reconstructive surgery.Current clinical practice involves the use of cephalometric planes, which is based on the manual selection of the landmarks by an experienced operator.Experience using this method has shown it to be timeconsuming and the results are rarely reproducible and repeatable; in addition, the best choice of the landmarks is still under debate.Although several approaches have been proposed in order to improve the localization of the MSP, to date, there are still many important limitations.The morphometric approach is still tied to the manual selection of the landmarks, which is timeconsuming and quite inaccurate; methods based on the ICP algorithm, which are computer automated, are too sensitive to asymmetries.
To overcome these drawbacks, a new, fully automatic technique which is able to extract an accurate and reliable symmetry plane is proposed.It is, indeed, an ICP-based method with a specifically designed objective function and minimum searching function, which has proven itself to be robust to local minima.The performances of the proposed method are quantified and compared with other recent techniques by analyzing 20 skulls, segmented from anonymous CT and provided by the Meyer Children Hospital of Florence (Italy): two healthy real skulls and 18 real defective skulls with large defects.This pool of 18 cases that includes (i) large unilateral defect, (ii) large bilateral defect and (iii) large defect crossing the MSP, has been deemed by hospital doctors as representative of the clinical cases they have faced over the last five years.
In Table 3 the Asymmetry Values are reported for all these test cases and for four methods compared.Experimental results demonstrate that the method leads to the accurate localization of the MSP, even in cases where there are very large missing areas and strong asymmetries without any user interaction and any landmark selection.These results were carried out with a computational time, which is closely related to the resolution of the surface model (number of points of the point

Conclusions
In recent years, the reliable localization of the MSP of the skull has attracted a great deal of research effort, since it is crucial to helping with the clinical diagnosis and treatment planning of craniofacial asymmetry, as well as with the planning of craniofacial reconstructive surgery.Current clinical practice involves the use of cephalometric planes, which is based on the manual selection of the landmarks by an experienced operator.Experience using this method has shown it to be time-consuming and the results are rarely reproducible and repeatable; in addition, the best choice of the landmarks is still under debate.Although several approaches have been proposed in order to improve the localization of the MSP, to date, there are still many important limitations.The morphometric approach is still tied to the manual selection of the landmarks, which is time-consuming and quite inaccurate; methods based on the ICP algorithm, which are computer automated, are too sensitive to asymmetries.
To overcome these drawbacks, a new, fully automatic technique which is able to extract an accurate and reliable symmetry plane is proposed.It is, indeed, an ICP-based method with a specifically designed objective function and minimum searching function, which has proven itself to be robust to local minima.The performances of the proposed method are quantified and compared with other recent techniques by analyzing 20 skulls, segmented from anonymous CT and provided by the Meyer Children Hospital of Florence (Italy): two healthy real skulls and 18 real defective skulls with large defects.This pool of 18 cases that includes (i) large unilateral defect, (ii) large bilateral defect and (iii) large defect crossing the MSP, has been deemed by hospital doctors as representative of the clinical cases they have faced over the last five years.
In Table 3 the Asymmetry Values are reported for all these test cases and for four methods compared.Experimental results demonstrate that the method leads to the accurate localization of the MSP, even in cases where there are very large missing areas and strong asymmetries without any user interaction and any landmark selection.These results were carried out with a computational time, which is closely related to the resolution of the surface model (number of points of the point cloud) and the complexity of the model itself.For the analyzed test cases, i.e., 3D discrete models with point numbers from 80,000 to 120,000, the computational time was between 10 and 25 minutes.The proposed method, as it has been presented, fails in the case in which the symmetry plane of the skull does not coincide with the MSP.Among all the analyzed test cases, this occurred only for a rare case of craniofacial dysmorphism due to hypertelorism combined with a severe form of plagiocephaly (Figure 8a).In this case, the user has to segment the area (highlighted in green in Figure 8b) for which the symmetry plane approximates the MSP one.However, even in this case, the selection is much faster and easier to perform than the landmark identification.
Symmetry 2019, 11, x FOR PEER REVIEW 11 of 13 cloud) and the complexity of the model itself.For the analyzed test cases, i.e., 3D discrete models with point numbers from 80,000 to 120,000, the computational time was between 10 and 25 minutes.The proposed method, as it has been presented, fails in the case in which the symmetry plane of the skull does not coincide with the MSP.Among all the analyzed test cases, this occurred only for a rare case of craniofacial dysmorphism due to hypertelorism combined with a severe form of plagiocephaly (Figure 8a).In this case, the user has to segment the area (highlighted in green in Figure 8b) for which the symmetry plane approximates the MSP one.However, even in this case, the selection is much faster and easier to perform than the landmark identification.Finally, future efforts will be addressed both to test this methodology in a much larger case study and to reduce the calculation time, before implementing it in a CT segmentation software.
Ba The most anterior point on the margin of the foramen magnum in the midsagittal plane Nasion N Most anterior point of the frontonasal suture in the mid-sagittal plane Crista Galli CG Most superior point of the Crista Galli Incisive Foramen IF The midpoint of the Incisive Foramen Frontozygomatic Suture FZS The most medial and anterior point of left (FZSL) or right (FZSR) frontozygomatic suture at the level of the lateral orbital rim Supraorbital Foramen SOF The midpoint of supraorbital foramen Fontorbitomaxillare FOM Lateral point of the frontomaxillary suture on the medial margin of the orbit Frontonasomaxillare FNM The intersection of the nasomaxillary, frontomaxillary, and frontonasal sutures Sella S Center of the sella turcica Pogonion Pog Most anterior point of the bony chin in the median plane Anterior Nasal Spine ANS Most anterior point midpoint of the anterior nasal spine of the maxilla
called S*, ANS* and Pog*, each of which is the midpoint of the segment joining the original landmarks S, ANS and Pog with their mirrored counterparts.Symmetry 2019, 11, x FOR PEER REVIEW 4 of 13 called S*, ANS* and Pog*, each of which is the midpoint of the segment joining the original landmarks S, ANS and Pog with their mirrored counterparts.

Figure 3 .
Figure 3. MSP resulting by methods proposed in [35] for (a) complete and (b) defective skull.

Figure 3 .
Figure 3. MSP resulting by methods proposed in [35] for (a) complete and (b) defective skull.

Figure 4 .
Figure 4.The flow-chart of the MaWR method as presented in [28].

Figure 4 .
Figure 4.The flow-chart of the MaWR method as presented in [28].

Figure 5 .
Figure 5.The MSP resulting from: (a) the MaWR method as presented in [28]; (b) the MaWR improved method.

Figure 5 .
Figure 5.The MSP resulting from: (a) the MaWR method as presented in [28]; (b) the MaWR improved method.

13 Figure 6 .
Figure 6.The real test cases here analyzed.(a) and (b) are two healthy skulls; (c) skull with large unilateral defect; (d) skull with large bilateral defect; (e) skull with large defect crossing the MSP.

Figure 6 .
Figure 6.The real test cases here analyzed.(a) and (b) are two healthy skulls; (c) skull with large unilateral defect; (d) skull with large bilateral defect; (e) skull with large defect crossing the MSP.

Figure 8 .
Figure 8. Skull with craniofacial dysmorphism due to hypertelorism combined with a severe form of plagiocephaly.(a) MSP resulting from the proposed method; (b) MSP resulting from the proposed method applied to the user selected area highlighted in green.

Figure 8 .
Figure 8. Skull with craniofacial dysmorphism due to hypertelorism combined with a severe form of plagiocephaly.(a) MSP resulting from the proposed method; (b) MSP resulting from the proposed method applied to the user selected area highlighted in green.

Table 2 .
The Asymmetry Value (AV) for the five test cases and four methods here discussed.

Table 2 .
The Asymmetry Value (AV) for the five test cases and four methods here discussed.

Table 3 .
The Asymmetry Value for the twenty test cases and four methods analyzed.

Table 3 .
The Asymmetry Value for the twenty test cases and four methods analyzed