Statistical Shape Model: comparison between ICP and CPD algorithms on medical applications

The increasing availability of 3D anatomical models obtained from diagnostic images exploiting Reverse Engineering techniques allows the application of statistical analysis in the quantitative investigation of anatomical shapes variability. Statistical Shape Models are a well-established method for representing such variability, especially for complex forms like the anatomical ones. Not by chance, these models are widely used for medical applications, such as guiding segmentation of the diagnostic image and virtual reconstruction of incomplete anatomic region. The application of a statistical analysis on a set of shapes representing the same anatomical region essentially requires that shapes must be in correspondence, i.e. constituted by the same number of points in corresponding position. This work aims to compare two established algorithms, namely a modified version of the Iterative Closest Point and the non-rigid version of the Coherent Point Drift, to solve the correspondences’ problem in the construction of a Statistical Shape Model of the human cranium. The comparison is carried out on the models using the standard evaluation criteria: Generalization, Specificity and Compactness. The modified version of the Iterative Closest Point delivers a better Statistical Shape Model in terms of Generalization and Specificity, but not for Compactness, than the Coherent Point Drift-based model.


Introduction
Reverse Engineering (RE) techniques are typically related to manufacturing industry applications, but, in the last few years, they are proving their effectiveness also in nontraditional fields, as biology and anatomy.The increasing ability to interpret and process the information acquired with diagnostic imaging techniques, such as Computed Tomog-raphy (CT) and Magnetic Resonance Imaging (MRI), is transforming the surgical approach: they deliver comprehensive in vivo information about the anatomical Region of Interest (ROI), consisting of multiple two-dimensional images stacked along the third dimension.This information allows a three-dimensional representation of the acquired anatomies, enabling complex and even more refined analysis and processing for purposes beyond the diagnosis.By exploiting RE techniques, it is possible to digitalize and process the 3D anatomical ROI from the patient's diagnostic images obtained by CT or MRI.The typical approach relies on the data extrapolation from the acquired image, involving the segmentation of the ROI.Starting from the binary image provided by the segmentation process (using commercial software such as Mimics, 3D Slicer etc.), a quantitative and interactive model can be obtained by exporting the model in a STL format.This is a widespread and powerful mean to describe complex shapes (such as the anatomical structures), representing the model as a set of points and their related connectivity list.The STL format is supported by many software packages and is widely used in RE applications and in Additive Manufacturing.Computer-Aided Technologies (CAx) software packages provide advanced tools to properly handle the so obtained 3D model allowing more effective preoperative simulation, complex-surgery planning, quantitative evaluation of asymmetry or dysmorphism and the design of the patient-specific devices.To this aim, the increased availability of real and interactive 3D anatomical models, allows the application of statistical analysis in the quantitative investigation and modelling of anatomical shapes variability.A statistical analysis (SA) is the most suitable approach dealing with the anatomical structures, because of the difficulty tied to their geometrical complexity along with a very wide interpersonal variability.

Statistical Shape Analysis
In the plethora of SA methods, the Statistical Shape Analysis (SSA) [1] is the most established methodology able to gather and interpret information of the same anatomical structure from an adequate number of individuals (hereafter referred to the term "samples").This approach provides a parametric model able to generate new valid shapes (i.e.still belonging to the same family of shape described by each sample), usually referred to as Statistical Shape Models (SSMs) [2].The SSM, introduced in the '90 s by Cootes et al. [3], accounts for the shape variations learned from the collection of n s samples organized in a proper designed dataset called Training Set (TS).These variations are represented by the principal components φ i (main modes of variation) and their respective variance values λ i , defined by applying the Principal Component Analysis (PCA) on the TS.
In order to apply the PCA the shapes in the TS are required to be properly represented and aligned in a global reference system.The simplest, but also the most generic, method to represent complex 3D shapes is to use a set of points spread across its surface.A 3D shape, consisting of n p points, can be therefore represented by a vector x of spatial coordinates (see Eq. 1): Using the SSMs, a new shape x * can be obtained as a deformed version of the mean shape x through a linear combination of deformations along the first c main directions of variation √ λ i φ i (see Eq. 2): where c represents the number of significant φ i and it is defined in order to obtain a cumulative variance between 90% and 98% of the total variability described by the TS (see Eq. 3): Thanks to its ability to encode the information of an adequate number of ROI and to summarize the information of the physiological and pathological shapes, the SSMs have already been widely used, for example, to guide the segmentation process and the reconstruction of deficient anatomical regions [2,4].
For the construction of the SSM it must be strictly assumed that the points of all the involved shapes (provided by a STL file) are in correspondence.This means that the correspondent points of all training samples must be located at the same positions.In the past, for simple 2D shapes with prominent features, the user manually defined these correspondences.However, dealing with a TS defined by complex 3D shapes with a wide variability among individuals (i.e. the anatomical ones), the manual definition of the correspondences is time consuming and prone to error.Therefore, several methods have been developed to automatically establish such correspondences.Among them, the most promising approach involves pairwise registration algorithms [2].In particular, a modified non-rigid version of the Iterative Closest Point (ICP) [5] and the non-rigid version of the Coherent Point Drift (CPD) [6] represents, to date, the most established alternatives.In order to detect the correspondences within a given TS, the pairwise registration algorithms perform a non-rigid registration between a previously defined template and each sample of such a TS.Usually, this registration is carried out using a first rigid alignment followed by an elastic registration.In this way, each shape of the TS is represented as a deformed version of the template, therefore constituted by the same number of points in corresponding positions.In [5] an ICP based non-rigid registration algorithm is presented, where the elastic registration is achieved by means of a global deformation, modelled as a sum of Gaussian Radial Basis Function, and refined by a weighted locally rigid transformation.The CPD algorithm, introduced by Myronenko and Song [6], it is a probabilistic method for rigid and non-rigid registration.In such a method, the template is represented as the centroid of a Gaussian Mixture Model (GMM) which is fitted, rigidly or not, on a certain sample by maximizing the likelihood.In the rigid case the GMM centroids is constrained to move coherently, in the non-rigid one the displacement field is forced to be regular.
The aim of this work is to compare the performances of the two algorithms presented above used for solving the correspondences' problem in the case of human cranium, by applying the evaluation criteria proposed by Davies in [7].This comparison is necessary because the correspondences' detection is often most challenging part in the SSM construction and the algorithm used has a great impact on the final model's quality.

Evaluation criteria
The evaluation criteria proposed by Davies [7] are considered as a gold standard for the comparison of different SSMs obtained by the same TS.Over the years, the methods adopted for evaluating these criteria have been frequently criticized; however, the concepts underlying them have proven to be very robust.The evaluation is based on three parameters that identify the ideal properties of a deformable model: Generalization (G), Specificity (S) and Compactness (C).The three parameters are defined as a function of the number of main modes c employed for the representation of the shape.In detail, the Generalization is defined as the ability of the model to reproduce shapes not included in the TS; it is estimated performing a series of n s leave-one-out tests on the TS, in order to obtain a reduced model, and measuring the mean similarity between the excluded shape x i and the reconstruction with the reduced model x * (c) (see Eq. 4): D represents the metric used to compute the similarity between two shapes and it will be discussed below.G can also be interpreted as the mean of the average reconstruction error between the model and an unseen shape.
The Specificity S represents the ability to generate only valid shapes; it is determined assessing the mean similarity among a set of n r randomly generated shapes x k * , using the Eq. 2 with α i ≤ ±3, and the closest shapes x i in the TS (see Eq. 5): The last parameter is the Compactness C that shows the cumulative variance included in the model normalized with the total variance of the TS (see Eq. 6): A model is compact if the associated variance is minimal and requires few main modes, thus parameters, to obtain a correct representation of the shape.
By varying the metric D the significance of similarity between two shapes changes.Such parameter should be chosen according to the comparison to be made or the specific application of the SSM.Since the aim of this work is to compare two algorithms of correspondences' detection, the selected metric D is based on these correspondences.Therefore, the similarity between two shapes is evaluated as the Mean Absolute Distance (MAD) [8] among corresponding points (see Eq. 7): where ρ * j (c) and ρ i j represent the jth corresponding point of x * and x i respectively.Using the MAD, the G, S and C values are independent of the number of points n p constituting each shape and are expressed in mm.In this way, it is also possible to compare SSMs made up of a different number of points.

Results
After CT images segmentation, starting from TS containing n s 100 pathologically unaffected adult crania two SSMs were built following the methodology described by Marzola et al. in [4], based on the previous work of Di Angelo [9].The latter defines a reliable procedure for the alignment and the scaling before the application of the correspondences' detection algorithms to ensure proper results.The two obtained SSMs differ only in the pairwise registration algorithm used, so by comparing them, it is possible to determine which of these two algorithms is the most suitable for the correspondences' detection in this specific application.The comparison was carried out by using the parameters described in the previous section.
As stated by Davies [7] lower value of G, S and C for all c generally characterize a higher quality model, therefore a more appropriate algorithm for the correspondence's detection.So, as shown in Fig. 1, the SSM obtained using the modified version of the ICP resulted better than the one resulting from the CPD in terms of G and S (n r 10000 randomly generated shapes), but not for C. In order to obtain a representation of 98% of the total variance the ICP-based SSM requires 71 main modes of variation, while the CPD-based SSM needs 77.However, since the values of Compactness are very similar and the value of G and S are lower for the ICP-based SSM for all c, it's possible to conclude that, in the case of human neurocrania, the modified ICP version provide a better SSM if compared to CPD version.In particular, the ICP-based SSM is able to reproduce unseen shapes with an average reconstruction error equal to 1.83 mm with c 71 (r 98%).In addition, the modified version of ICP algorithm is preferred to the CPD as it required less time to register the template to each sample, then to detect the cor- respondences among all the samples of the TS and to build the SSM.

Conclusions
This work was meant to compare both a modified version of the Iterative Closest Point and the non-rigid version of the Coherent Point Drift to solve the correspondences' problem in the construction of a Statistical Shape Model of the human cranium.Starting from 100 pathologically unaffected adult crania two SSMs, different only for the correspondences' detection algorithm employed, were built and compared using the standard evaluation criteria of Generalization, Specificity and Compactness.The ICP-based SSM has proven to be preferable than the CPD-based in terms of Generalization and Specificity for all the significant modes of variation, but not for Compactness.In addition, the ICPbased SSM is able to reproduce an unseen cranium with an average reconstruction error equal to 1.83 mm using 71 main modes of variation.
The final aim was to understand which of the two literature techniques are the best suited to be adopted for the construction of an SSM of the human cranium to use as template in solving complex defective skulls reconstruction [10].Future works will be addressed to analyse in detail the meaning of SSM quality and the factors affecting it, depending on SSM's application.
Funding Open access funding provided by Università degli Studi di Firenze within the CRUI-CARE Agreement.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material.If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

Fig. 1
Fig. 1 Generalization Specificity S and Compactness C of ICP-based model (square) and the CPD-based model (circle)