Statistical, Morphometric, Anatomical Shape Model (Atlas) of Calcaneus

The aim was to develop a morphometric and anatomically accurate atlas (statistical shape model) of calcaneus. The model is based on 18 left foot and 18 right foot computed tomography studies of 28 male individuals aged from 17 to 62 years, with no known foot pathology. A procedure for automatic atlas included extraction and identification of common features, averaging feature position, obtaining mean geometry, mathematical shape description and variability analysis. Expert manual assistance was included for the model to fulfil the accuracy sought by medical professionals. The proposed for the first time statistical shape model of the calcaneus could be of value in many orthopaedic applications including providing support in diagnosing pathological lesions, pre-operative planning, classification and treatment of calcaneus fractures as well as for the development of future implant procedures.


Introduction
Improving the diagnosis and supporting the therapy is a fundamental task of medical imaging [1,2]. Planning the surgery and choosing the method of surgical treatment often requires deep knowledge of the shape and morphological characteristics of an object at hand. In case of traumatology the detailed assessment of the morphological characteristics of bone is important for the choice of stabilization method and subsequent healing process. Approximately 10% of all fractures occur in foot bones [3]. The calcaneus (lat.os calcis) fracture is the most common fracture of tarsal bones, with an incidence of about 60% [4,5]. Calcaneus fractures are estimated at 2% of all fractures of human bones and can have long-term consequences for patients in terms of comfort and disability [6,7]. This largest tarsal bone plays stabilizing function for human body and pathologies of this bone influence human mobility and quality of life [8][9][10]. The evaluation of calcaneus fractures and determination of the effect of treatment essentially uses three types of imaging techniques: X-ray CT, and MRI [11][12][13][14][15]. Although X-ray study (radiography) remains the method of choice for the initial assessment of fracture [16,17], CT imaging has become the current state-of-art in diagnostics of calcaneus fractures [18][19][20]. An invaluable advantage of CT is the possibility of reconstructing the scanned object. The location and configuration at a joint can be provided having 3D morphological and architectural information about the individual foot bones [21]. Several studies considered characterizing the calcaneal bone shape. Early works of Stindel et at. [21,22] focused on classification of foot type in terms of pathological deformations using MRI. Stephan et al. [23] obtained information on the joint surfaces and 3D orientations of calcaneus using CT to quantify the integrity of calcaneus joint surfaces. Guterkunst et al. [24] aimed at assessing the measurement precision of landmark-based 3D bone-to-bone orientations of hind foot and lesser tarsal bones to compare an atlas-based automated method to the expert rating. Qiang et al. [7] described morphological characteristics of the calcaneus based on CT images. The survey in the later the paper has an epidemiological value. They described calcaneus by size and distances between characteristic points on the surface that were manually outlined by an experienced operator. Hence, there are possibilities of bone shape description from which a model of calcaneus could be constructed. This subsequently can be used for building statistical anatomical atlases (statistical shape models, SSMs) [25,26]. Many tasks in medical imaging involve automatic systems which use prior knowledge about an object. Model based methods have become popular because they have the potential to resolve possible confusion associated with structural complexity, to improve tolerance for noisy data and deliver more accurate results in case of missing data [27][28][29][30][31][32]. In medical imaging, SSMs are mainly used in segmentation and recognition tasks. An atlas guide/ model-based segmentation approach is recognized as one of the most successful methods for image analysis [33][34][35][36]. So far, applications of an atlas based model include brain structures [37][38][39][40][41], soft tissue structures in the abdominal and pelvic area [32,42], cardiac structures [29,43,44], and several bone structures [45][46][47], but not of a calcaneus. Shape representation is one of the most important issues in SSM based methods of medical image processing. Several shape representation forms are known [48]. Zhang and Liu [49] classify shape representation and description techniques into region-based methods and contour-based methods. Region-based methods consist of global methods such as geometric moment invariants, algebraic moment invariants, orthogonal moments, generic Fourier descriptors, grid based methods, and shape matrix as well as structural methods such as convex hull and medial axis. On the other hand, contour-based techniques exploit shape boundary including global methods such as simple shape descriptors (area, eccentricity, major axis orientation, blending energy), correspondencebased shape matching in space domain, shape signature (centroidal profile, complex coordinates, centroid distance, tangent angle, cumulative angle, curvature, area), boundary moments, elastic matching, stochastic method (autoregressive modeling), scale space method, and spectral transform as well as structural methods such as chain code representation, polygon decomposition, smooth curve decomposition, scale space method, syntactic analysis, and shape invariants. Shape representation forms are not limited to the classification proposed by Zhang and Liu [49]. Cottes et al. [25,27,28,50] classified them as "Hand Crafted" models (built from specific scalable objects like lines, ellipses and arcs) [30], articulated models [51], active contour models [52], Fourier-based shape models [53], statistical models of shape, and final elements models [54,55]. Also, for the 3D shape representation, Heiman and Meizer [56] considered medial models (pioneered by Pizer et al. [57]) that allow rendering 3D solids and spherical harmonics models [58][59][60].
The purpose of this work was to develop an anatomically accurate SSM of calcaneus which includes its morphological characteristics and provides mathematical representation of its shape.

Measurement setup
Retrospective data from regular hospital records have been used in this study. All patient records were anonymized and de-identified prior to processing according to the standard data release procedures of the two hospitals involved in the study. The Review Board of the Department of Radiology, Wroclaw Medical University, Wroclaw has approved the study. Poland Volume data have been acquired with three commercially available CT scanners ( [61] to extract physical dimensions of calcaneus. Volume data of 18 left and 18 right feet of 28 male subjects were collected. Subjects, aged from 17 to 62 years (mean 36.8 ± std 21.2 years), had no known foot pathology in the scanned foot. Subjects were aged matched, so they were no statistically significant differences between the age group averages of the left and the right foot subjects.

Data processing
A scheme for automatically building a morphometric and anatomical atlas was described by Subsol et al. [2]. It consists of the following stages: feature extraction, common feature identification, averaging feature position to obtain mean geometry and variability analysis. The process for building a morphometric and anatomical atlas of calcaneus is shown in Fig 1. In the following, we describe in detail the particular steps of the process. All procedures have been developed in Matlab (MathWorks, Inc., Natick, MA, USA) and are available from authors for free upon request.
1. Image pre-processing [ Fig 1A and 1B] The volume image I is decomposed in the sagittal plane into a series of 2D images I k , k = 1, 2, . . ., n. Each I k is first normalized in two steps: and then contrast is enhanced using sigmoid function [62]: where g is a gain, and c is cut-off value, I kmax is maximal and I kmin is minimal pixel intensity value. Both parameters were set empirically to 0.4 and 0.7, respectively.

Contour extraction [Fig 1C]
Extracting the contour of calcaneus in each image I k can be achieved in a variety of ways.
For the purpose of this study we use the region growing algorithm [63,64] with manually assisted assignment of a seed point by an expert operator.
3. 3D point cloud to surface [ Fig 1D] The result of the region growing algorithm is a set of coordinates W = [x, y, z], with x = [x1, x2, . . ., xn] 0 , y = [y1, y2, . . ., yn] 0 , and z = [z1, z2, . . ., zn] 0 representing the contour of the calcaneus gathered respectively from each I k , k = 1, 2, . . ., n. For such a point-cloud set the oriented normals are estimated and then using Poisson surface reconstruction method the surface mesh is generated [65]. Meshlab (Pisa, Italy) software was used to generate meshes [66]. Surface representation is needed to further land-marking of the characteristic points of the calcaneus. The objects of study are not scanned in the same position so their orientation and location need to be normalized for the statistical bone atlas to be built. For this, the point-cloud set W and the characteristic points C, D, E are used. First, the plane π 1 that includes points C, D, E is obtained This is followed by calculating the angle α between the plane π 1 and the plane π: z = 0. Then is applied and the angle β between x axis and vector C¼ ½x C ; y C ; z C is calculated. Finally, W is rotated about the x axis by the angle β using the rotation matrix: 6. Spherical harmonics decomposition [ Fig 1G] The mathematical shape description is achieved by using the spherical harmonics (SPHARM) [58,[67][68][69][70]. Using linear-in-parameters least square procedure, data (cloud of extracted points) are first fitted with a best sphere: and then translated to a new coordinate system centred at the estimated origin ðx 0 ;ŷ 0 ;ẑ 0 Þ: The center of the new origin is set as a center of the mass: Then, transformation from the Cartesian to spherical coordinates is made: , where θ d is the azimuth, ϕ d is the elevation angle and R d is the radius of sphere.
Then the normalization to a unit sphere is performed: The complete complex-valued spherical harmonics are solutions to the Laplace equation: are given by: where i ¼ ffiffiffiffiffiffi ffi À1 p ; N m l is the normalization factor defined as: P m l are the associated Legendre functions: The real-valued spherical harmonics can be expressed as: where l = dde -1 and m alternates between -l, l + 1, . . ., 0, l-1, l. The coefficient estimates in the spherical harmonic expansionĉ q , q = 1, 2, . . . Q can be easily evaluated using the method of linear least-squares by concatenating the data into a single column vector as described earlier. The estimate of the spherical harmonics based calcaneus surface is thenR Transforming this result back to Cartesian coordinate system: and then translating it back to the starting origin (see Fig 1G and 1H) creates an estimate of the calcaneus and allows to evaluate the quality of the approximation at the original points (x d , y d , z d ), d = 1, 2, . . ., D.
7. Model and model order selection [ Fig 1H] The optimal model order of the SPHARM expansion can be evaluated using one of the standard information criteria such as Akaike Information Criterion or the Rissanen Minimum Description Length (MDL) criterion [71,72]. For the purpose of this works the latter criterion was used. Given the optimal model order, further reduction of the number of coefficients in the SPHARM expansion is achieved by applying a multiple hypotheses procedure to test each of the coefficient values for zero.

Results
The MDL criterion estimated the optimal radial order of SPHARM expansion at 11 resulting in 144 characteristic coefficientsĉ m l . Fig 2 shows   shapes is obviously very high (r 2 = 0.99) and statistically significant (p ( 0.01). This correlation remains moderate when the three highest coefficients omitted amounting to r 2 = 0.59, p ( 0.01 (see Fig 5). Application of a multiple hypotheses procedure to test each of the coefficient values for zero reduced the set of coefficient to 41 and 47 representative coefficients for the right and left calcaneus, respectively.   No statistically significant differences were found between the parameters of the left and right foot (t-test, p > 0.05).

Discussion and Conclusions
Building an anatomically accurate statistical shape model of calcaneus requires close collaboration of image processing engineers, radiologists and orthopaedists. Following a general scheme for building such an atlas in an automatic function, as proposed by Subsol et al. [2], we were able to bring down the characterization of the mean geometry of calcaneus to several dozen representative coefficients of the spherical harmonic expansion. It should be noted that the scheme was not fully automatic and that the assistance of the experienced operator was necessary in two critical steps of building the model. That is, in the process of contour extraction where the selection of the seed point was manually assisted and in the difficult and tedious characterization of the representative anatomical features of the calcaneus. It is worth noting that the later was in a close agreement with the results of Qiang et al. [7] who also manually   Statistical Atlas of Calcaneus identified the anatomic landmarks of calcaneus. Shapes may be represented in many different ways and an argument can be made whether the choice of spherical harmonic representation was the most appropriate. Given the many advantages of spherical harmonics [69] such as their orthogonality and completeness, as well as the ability to estimate the coefficients using linear-in-parameters least-squares method provides a sufficient motivation for using this specific representation. It is important to acknowledge that other forms of 3D basis functions could also be used in place [73]. However, comparing different representations in terms of goodness-of-fit and the minimum number of representative coefficients, although interesting, was not within the scope of this work. One of the limitations of the study was the lack of complete CT data for both feet. It is expected that otherwise the symmetry of the model, as estimated via the two sets of spherical harmonic expansions, would be more evident. Another drawback was the acquisition of a single CT scanning per subjects (to limit the dose of radiation and cost), which prevented assessing the reproducibility of the measurement and modelling procedure. Nevertheless, other studies related to foot bone imaging with CT show relatively high test-retest reliability [24] that can be anticipated to hold for imaging calcaneus. Finally our study examined a set of male calcanei but the exact procedure can be readily applied for the female counterparts. To the best of our knowledge, the proposed statistical anatomical atlas of the calcaneus is presented for the first time. Tools of image processing, 3D shape modelling, and expert manual assistance were necessary to achieve the goal. The manual assistance in building a statistical shape model may appear to be a hindrance, but could be attempted, after verification, with robust shape registration methods [74,75]. However, at the same time it provides an assurance of anatomical accuracy sought by medical professionals. In this view, the proposed statistical model of calcaneus including the quantitative mathematical shape description could be of value to diagnosing pathological changes by matching a model which contains information about the expected shape [76,77]. It can be applied to monitoring growth process in childhood, bone and surrounding tissue age degeneration and osteoarthritis progression, as well as to detection of lesion (cysts, tumour-like changes). The model of calcaneus finds the application in pre-operative planning [78,79], classifications [6] and follow-up of treatment of calcaneus fractures [8,80], assessing the implant procedures [81,82], and 3D reconstruction of bones [83]. Consequently, the model could also be used in image processing and computer vision software supporting those medical imaging tasks like segmentation or recognition.
As the summary demonstrates, there is already a considerable amount of 3D SSMs employed in medical image analysis. Extending this area to models of calcaneus could further support medicine in better diagnosis and treatment.
Supporting Information S1 Data Set. The full data set for all 144 SPHARM coefficients (Fig 2). (XLS) S2 Data Set. Data used to reconstruction of 3D shape for the right calcaneus (Fig 3). (XLS)