Application of a Novel Automatic Method for Determining the Bilateral Symmetry Midline of the Facial Skeleton Based on Invariant Moments

: Assuming a symmetric pattern plays a fundamental role in the diagnosis and surgical treatment of facial asymmetry, for reconstructive craniofacial surgery, knowing the precise location of the facial midline is important since for most reconstructive procedures the intact side of the face serves as a template for the malformed side. However, the location of the midline is still a subjective procedure, despite its importance. This study aimed to automatically locate the bilateral symmetry midline of the facial skeleton based on an invariant moment technique using pseudo-Zernike moments. A total of 367 skull images were evaluated using the proposed technique. The technique was found to be reliable and provided good accuracy in the symmetry planes. This new technique will be utilized for subsequent studies to evaluate diverse craniofacial reconstruction techniques.


Introduction
Bilateral symmetry refers to a structure of interest having two sides, with one side being the mirror image of the other. A good example is the human skull, which has bilateral symmetry of its left and right sides. The goal of reconstructive surgery for clinical pathologies of the craniofacial skeleton, whether boney deformities, mandibular alterations, tumors, or trauma, is to restore bilateral symmetry in order to restore function. Similarly, clinical diagnosis and treatment for orthodontics [1], maxillofacial [2], cosmetic [3], and plastic reconstructive surgery [4] seek to restore craniofacial bilateral symmetry. In preceding studies, a method for finding the midline of the craniofacial skeleton was developed by the authors of [5]. This is now a generally accepted method for establishing the midline of the face. However, accurately locating the craniofacial midline is essential for correctly forming an intact template for correcting facial malformations or trauma and for the planning of procedures in the reconstruction process. To date, the most common method for locating the midline for a two-dimensional image, or the midsagittal plane (MSP) for a three-dimensional object, has been the method proposed by the authors of [6].
Midline symmetry planes have been calculated as a perpendicular midpoint for a stack of horizontal lines crossing bilaterally across the facial skeleton containing boney landmarks or by using boney landmarks to construct a midline that best represents facial symmetry [7]. A simple way to define the midline has been to manually select a number of cephalometric boney landmarks in the dataset, either directly on the plane or at equal distances on either side of the midline plane. However, this requires great attention and care during the selection process by an expert user. Manual selection of skeletal landmarks is time-consuming and unreliable, which results in less accurate estimations of symmetry. In addition, the results have been dependent on the user's ability to find appropriate landmarks, the landmark availability and the visibility of the anatomical landmarks [8].
Methods with some amount of automation have been attempted for finding the midline plane. The authors of [9] described a semi-automatic approach for calculating the symmetry plane of the facial skeleton using principal component analysis (PCA) and the iterative closest point (ICP) alignment method along with surface models reconstructed from computed tomography image data. The initial step was to determine the precise position of the mirror plane using PCA in order to approximately align the mirrored mesh and the original mesh. The ICP algorithm was then used to get a refined registration [10]. The disadvantage of this method was the need for a central point of the image for the calculation of the symmetric plane. This point was obtained using the average of the vertices of the facial mesh. If this point was not provided or if the point was in wrong location, in case of imperfect symmetry, this method could lead to an error in the construction of the symmetric plane. A second limitation lied in the lack of self-learning or the ability to learn. Once applied to an image, the subsequent image will not learn or improve its performance based on the previous data.
Most recently, an upgraded version of a mirroring and registration technique for automatic symmetry plane detection of 3D asymmetrically scanned human faces was discussed in [11]. This work described an ICP-based method that uses particle swarm optimization. This method starts from a 3D discrete model and evaluates the symmetry plane by a preliminary first-attempt, carried out with a PCA algorithm, which is then refined iteratively until its final estimation obtained by a Levenberg-Marquardt algorithm. This new proposed method was claimed to improve upon the limitations from the authors' previous algorithm [12]; however, the new model still has limitations and fails to incorporate self-learning to improve the model's outcome. For instance, in the case of craniofacial dysmorphism, this model requires the user to interactively segment the area to solve the asymmetry. However, if the algorithm is presented a similar image again, another intervention is required to select the area.
The aim of this study was to present a technique for the automatic calculation of the craniofacial symmetry midline using invariant moments. Figure 1 shows the steps of the proposed algorithm. First, based on the cephalometric landmarks, the image is rotated so that the midline passes through the center, i.e., 0 • . The image is then duplicated and a dataset of 30 images of 1 • resolution is created from −14 • to 15 • . This set of images is then passed through different feature extractors and then to the k-nearest neighbors classifier. In the classification phase, pseudo-Zernike moments (PZMs) and PCA were selected after going through the independent component analysis (ICA) feature extractor. A second set of images of 0.5 degrees resolution was taken from the centered images to test the PZMs and PCA. PZMs were selected for having achieved the best accuracy. After the classifier determines the rotation degree of these images, their midpoints are found. Finally, by joining the midpoints and grades described by the classifier, the midlines can be constructed.

Image Creation
To create the dataset, the IdentifyMe database [13] was used. This database contains 464 skull images composed of skull and digital face images (real-world examples of solved skull identification cases) and unlabeled supplementary skull images. The database also contains several unsuitable images that needed to be corrected before being added to the dataset. Additionally, this step intended to eliminate regions of non-interest through cropping and the correction of image rotation so that the processing efficiency of the system could be enhanced ( Figure 2). As the cropping process can be easily solved, the focus was to create a method to correct the rotation, since it will be necessary for training purpose. The first step was to manually identify six cephalometric landmarks (1-Crista Galli, 2-Frontozygomatic Suture, 3-Orbitale, 4-Anterior Nasal Spine, 5-Subspinale, and 6-Prosthion) ( Figure 3a). Then, using a grid as reference, the images were rotated using the nearest neighbor interpolation technique so that landmarks 1, 4, 5, and 6 were vertically aligned and landmarks 2 and 3 were horizontally aligned with their respective counterpart as in Figure 3b. These images were used as a gold standard to create the rotated images. After the images were vertically aligned, 367 images were selected as suitable for use. Based on each image, 30 images were created with inclination angles from −14 to 15 degrees with 1 degree of variation along the sagittal plane, totaling 11,010 images ( Figure 4). Those angles were also used as the image labels in the classification step. Finally, the images were resized to 128 × 128 using the nearest neighbor interpolation method. A resolution of 128 × 128 was selected so as to reduce the processing time and computational energy during the creation of PCA and ICA feature vectors. Due to the limitations of the i7-8850H CPU (2.60 GHz) with 16 GB RAM computer used, a resolution of 256 × 256 resulted in the computer crashing.

Feature Extractors
Three feature extraction methods (PZMs, ICA, and PCA) were compared to determine the method with the leading accuracy. The resultant method was then used for the algorithm to determine the midpoint to generate the final midlines.

Pseudo-Zernike Moments-PZMs
Zernike moments (ZMs) were first introduced in computer vision by the authors of [14] and are widely used to identify and highlight global features of an image in image processing and machine learning applications [15]. Zernike moments map an image into a set of complex Zernike polynomials. As these Zernike polynomials are orthogonal to each other, Zernike moments can represent the properties of an image without redundancy or overlap of information between moments. Pseudo-Zernike moments are a derivation of the Zernike moments that was shown to be more robust and less sensitive to image noise [16].
Pseudo-Zernike moments, i.e., pseudo-Zernike polynomials, are one of the best invariant image descriptors belonging to a family of circularly orthogonal moments due to their minimal information redundancy, robustness to image noise, provision of twice the number of moments, and having more low-order moments. In addition to being rotational invariants, they can be made into scale and translation invariants after certain geometric transformations. PZMs have been used in numerous machine learning and image analysis applications, such as cephalometric landmark detection [17], Alzheimer's disease detection ( [18][19][20]), medical image retrieval [21], detection of tumors and brain tumors ( [22][23][24]), facial expression [25] and facial age recognition [26], facial recognition ( [27,28]), and other industrial applications [29].
As described by the authors of [15], to calculate the PZMs, the image (or region of interest) is initially mapped into a unit disc, where the center of the image is the origin of the disc (Figure 5a). Pixels that are outside the disc are not used in the calculation. To include these pixels ( Figure 5b) the disc can be expanded so that the image function f(x,y) is completely enclosed inside the disc by performing the following: x, and y j = where x = y = 0, 1, ... (N − 1), x i , and y j are the image coordinates, and the image function f(x,y) is defined over the discrete square domain N × N. Obtaining the PZMs [15] of an image begins with the calculation of the pseudo-Zernike polynomials (PZP) V pq (x i ,y j ) of order p and repetition q: where p is a non-negative integer, 0 ≤ |q| ≤ p, θ = tan −1 (y/x), θ ∈ [0,2π], and r = x 2 i + y 2 j . The complex function V pq (x i ,y j ) has two separate parts: the radial polynomials R pq (r) and the angular functions e jqθ = (cosθ + j sinθ) q polynomials in cosθ and sinθ. The radial polynomials are expressed as: Since R pq (r) = R p,−q (r), we consider q ≥ 0 and rewrite Equation (3) as: where A p-recursive method focusing on reducing time complexity through the fast computing of the pseudo-Zernike radial polynomials using a hybrid method was presented by the authors of [30], where: R (q+1)q (r) = (2q + 3)r q+1 − 2(q + 1)r q , R pq (r) = (k 1 r + k 2 )R (p−1)q (r) + k 3 R (p−2)q (r), where The PZMs of order p and repetition q of an image function f(x i ,y j ) over a unit disc are represented in a discrete domain by:

Independent Component Analysis-ICA
ICA is a blind source separation or statistical signal processing technique where a given measurement is represented by a linear composition of statistically independent components (ICs) that aims to linearly decompose a random vector into components that are as independent as possible [31]. Given a set of random variables {x 1 (t), x 2 (t), ..., x n (t)}, where t is time or the sample index, it is assumed to be generated as a linear mixture of ICs {s 1 (t), s 2 (t), ..., s n (t)}: where A is an unknown mixture matrix A R n×n . The FastICA algorithm [31] was used to perform the task through the columns of its mixture matrix, which contain the main feature vectors.

Principal Component Analysis-PCA
PCA, first introduced by the authors of [32] and then independently developed by the authors of [33], is a technique that preserves the variation present in a large number of interrelated variables of a dataset while reducing the dimensionality. This is achieved by transforming the dataset into a new ordered set of uncorrelated variables or principal components (PCs), so that all the original variables can be represented by the first few variables that maintained the most important features [34].

Geometric Moments-Central Point
In current applications, the central point of the face is obtained manually. Thus, we propose the application of the geometric moments method for the automatic extraction of the central point. Geometric moments are a popular method of moments and have been used to identify the image centroid in a number of image processing tasks [35,36].
Two-dimensional (p + q)th order moments of a digitally sampled N × M image that has the gray function f (x, y), (x, y = 0, ...M − 1) [37] is given as: x p y q f (x, y), p, q = 0, 1, 2, 3, ... . (14) As described in [38], the mass and area of the zeroth order moment, M 00 , of a digital image is defined as: The center of mass of the image f(x,y) is represented by the two first moments: Thus, the centroid of an image can be calculated by: As best practice, the center of mass was chosen to represent the position of an image in the field of view. The centroid of the image f(x,y), given by Equation (18), can be used to describe the position of the image in space by using the point as a reference point.

Classification
PZMs feature vectors were generated through the first 20 PZM orders and repetitions, based on Equations (1)- (12), totalizing 121 features. After several tests, the 24 best features were selected to represent PZMs.
Before being presented to the ICA and PCA feature extractors, the images were vectorized. For ICA, the mixing matrix A, described in Equation (13), was used as its feature vectors. For PCA, the eigenvectors based on eigenvalue order on the covariance matrix or its principal components were used. In total, feature vectors of size 35 were used to represent ICA and PCA.
To make a reasonable comparison among the descriptors, the feature vectors were presented uniquely to the k-nearest neighbors classifier (k-NN) using the Euclidean distance and the eight nearest neighbors. In order to evaluate the performance of the feature extractors, the size of the training sets was varied from 10% to 80% of the available database and the rest was used for testing purposes. The k-NN outputs were the predicted images angles and Figure 6 presents the comparison of the accuracy rate among ICA, PCA, and PZMs. Accuracy was calculated using the following model: where TP-true positives, TN-true negatives, FP-false positives, and FN-false negatives. The extractor based on the ICA obtained a low accuracy rate with a maximum accuracy of 14.5% after training, making it unviable for this application. Thus, a further comparison could be carried out with PZMs versus PCA only to evaluate and finally select the best feature extractor.
Using the selected 367 images and the same inclination angles (−14 • to 15 • ), a new set of 59 images from the original images was created with a 0.5 degree resolution, totaling 21,653 images, and the results are presented in Figure 7. From these results it can be seen that: 1. PCA performed well, but it required more coefficients to achieve a performance similar to pseudo-Zernike moments (Figures 6 and 7); 2. In Figure 6, ICA estimation achieved a bad performance in the experiment. We attributed this to the fact that it is inherently affected by the rotation of the images, which has been already explored [39]; 3. In Figures 6 and 7, pseudo-Zernike moments outperformed ICA and PCA, as PZMs maintained almost all of the images' features in a few coefficients. An initial conclusion is that the rotation invariant feature descriptors in the image plane can be effectively developed using the pseudo-Zernike moments method, which has performed well in these scenarios. Thus, the superiority and choice of the feature extraction based on pseudo-Zernike moments in this application became obvious.

Midpoint Calculation
By using Equations (15)- (18), the center of the image was calculated. To validate the accuracy of the midpoint technique, visual correctness was used compared with the cephalometric landmarks. However, a calculation of the image center may be performed manually and compared with the midpoint. Figure 8 shows three images with rotations of 0 • , 3 • , 6 • , 9 • , 12 • , and 15 • with centers (midpoints) represented by red stars and the angle directions obtained in the process of classifying by the pseudo-Zernike moments represented by black dots. To calculate the PZMs, the 21,653 imagesdataset (obtained from a 0.5 degree resolution) was divided into 80% for training and 20% for testing, with a total of 4331 testing images. The PZMs method's accuracy was 98.64%, and 1.36% of the images received wrong labels (59 images). The images with wrong labels obtained an error of 0.5 • for 33 images, 1 • and 1.5 • for four images each, 2 • , 2.5 • , 10 • , 13.5 • , 14 • , and 15 • for one image each, 3.5 • for three images, and 3 • for two images.
Once the midpoints (red stars) and the angles predicted by k-NN based on PZMs (black dots) were calculated (Figure 8), the symmetrical line could be easily constructed. The proposed technique presented good results in obtaining the bilateral symmetry midline of the images. However, there are a few limitations to the proposed technique. A small deviation in obtaining the center line could be seen in 59 images. This error is likely due to the fact that some images suffer from a small rotation in the sagittal plane. Additionally, if the images suffer from deformation, incompleteness, or non-uniform brightness, the image center calculus becomes difficult to perform since the moment is a quantitative measure of the function or shape of an image. Furthermore, the algorithm was not tested on non-symmetrical skull images due to a lack of non-symmetrical skull datasets. Our lab will be collecting non-symmetrical skull images to test the algorithm further. Another limitation for the proposed method is the resolution size of the images (128 × 128). This resolution was selected as any resolution higher than this resulted in an error from the PC (not enough memory). A higher resolution may result in more accurate results. However, using the lower resolution of 128 × 128 created acceptable results.
From the viewpoint of the angles, PZMs performed with 98.64% accuracy in the angles and 1.36% (59 images) with wrong angles. However, 33 of the images (55.93%) had an error of only 0.5 • . The majority of this error could be found in −0.5 • , 0 • , and 0.5 • , where 14 images (42.42% of the 33 images related to 0.5 • ) had an error of 0.5 • .
It is possible to state that the required resolution of 0.5 • reduces the accuracy of the feature extractor, and this error can be related to the image size, which was resized to improve the processing of the image matrices. However, these errors can be disregarded as they would be imperceptible to the human eye.

Conclusions
This study proposed an automatic technique for determining the bilateral symmetry midline of the facial skeleton based on invariant moments. A total of 367 skull images were evaluated after preprocessing using the IdentifyMe database. A comparative study between pseudo-Zernike moments, independent component analysis, and principal component analysis as feature descriptors of images using k-nearest neighbors and Euclidean distance was performed and the study of the feature extraction step revealed that pseudo-Zernike moments for feature description had the best performance. PZMs offer an alternative to conventional landmark-based symmetry scores that depends on the general positions of cephalometric landmarks. PZMs are also an alternative to PCA-ICP techniques, which depend on the manual selection of the central point and cannot be improved. With the proposed technique, the central point could be found as the centroid of an image, and then the symmetrical midline could be constructed.
In this study, we have shown the proposed technique to be reliable and to provide the midline symmetry plane with great accuracy, which can be used to aid surgeons in reconstructive craniofacial surgeries.