Gérard Subsol

Jean-Philippe Thirion

Nicholas Ayache

A scheme for automatically building three-dimensional morphometric anatomical atlases: application to a skull atlas

Keywords: automatic diagnosis, computer-aided surgery, crest lines, digital atlas, morphometry, non-rigid registration, skull

We present a general scheme for automatically building a morphometric anatomical atlas. We detail each stage of the method, including the non-rigid registration algorithm, three-dimensional line averaging and statistical processes. We apply the method to obtain a quantitative atlas of skull crest lines. Finally, we use the resulting atlas to study a craniofacial disease; we show how we can obtain qualitative and quantitative results by contrasting a skull affected by a mandible deformation with the atlas.

INTRODUCTION

To improve diagnosis, treatment planning, delivery and follow-up, a physician needs to compare three-dimensional (3-D) medical images from various modalities [computed tomography (CT), magnetic resonance imagery (MRI) or nuclear medicine (NM)] [START_REF] Ayache | Medical computer vision, virtual reality and robotics[END_REF]. We can distinguish three kinds of comparisons: comparison of a single patient's images to study the evolution of a disease; comparison of different patients' images to contrast a healthy and a sick person and registration of images with an anatomical atlas to facilitate the anatomical interpretation.

This last type of comparison is necessary to identify and locate precisely the various anatomical structures of the patient. It also allows one to study potential variations from 'standard' anatomy.

Limitations of conventional anatomical atlases

For years, medical doctors have used books such as Pernkopf (1983), nevertheless, such atlases are quite difficult to use, especially by a non-skilled person, due to:

• Two-dimensional (2-D) representation. Images are usually in 2-D and are taken from a single point of view.

The 3-D shape of the anatomical structure has to be reconstructed mentally, which is a task that requires a lot of experience.

• Ambiguous landmark definition. In general, such atlases are based on different kinds of features: points (e.g. apex), lines (e.g. sutura, crista) or areas (e.g. pars, foramen). Moreover, the definition of these features can vary according to the observer.

• Lack of a quantitative description. Most of these atlases provide only a qualitative description; they do not give much quantitative information about feature position. In fact, localization of the anatomical structures is based only on the relationship between features requiring expert anatomy. One exception is the stereotactic brain atlas developed by [START_REF] Talairach | Co-planar stereotaxic atlas of the human brain[END_REF]. Moreover, conventional atlases usually present the description of only one patient. So, it is impossible to estimate the statistical distribution of the position, the size or the topology of anatomical structures.

These limits are particularly emphasized when physicians want to use atlases with volumetric medical images. They are then obliged to compare either 2-D slices that are not taken in exactly the same position or a 2-D plate with a 3-D rendering of the medical image.

Current digital anatomical atlases

To overcome the limitations of conventional atlases, some digital anatomical atlases based on 3-D medical images have been developed over the last 10 years. These can be separated into three classes of application:

• 3-D database. These atlases, e.g. Voxel-Man [START_REF] Höhne | A volume-based anatomical atlas[END_REF]Schiemann et al., 1996b) which is already marketed, include powerful volume visualization techniques to display anatomical data from any viewpoint with functionalities such as cutting or transparency.

They also integrate a sophisticated structure labelling thanks to a database engine. These atlases can be used to consult reference cases and, above all, for teaching anatomy. Nevertheless, they remain a digital version of conventional atlases [or of several atlases as in [START_REF] Nowinski | Talairach-Tournoux/Schaltenbrand-Wahren based electronic brain atlas system[END_REF]] with a qualitative description of only one patient. • Normalized registration. With non-rigid 3-D registration tools, it becomes possible to register the patient images with the atlas in order to locate anatomical structures (Christensen et al., 1996a). In [START_REF] Marrett | A volume of interest (VOI) atlas for the analysis of neurophysiological image data[END_REF], a manual method is proposed to register an atlas with MRI data: first a global affine matching is performed manually then the user can choose one (or a set of) volumes of interest in order to apply an affine transformation locally. [START_REF] Greitz | A computerized brain atlas: construction, anatomical content and some applications[END_REF] introduced quadratic transformations (called pear, skew, asymmetry or scoliosis) in addition to linear ones (translation, rotation, scalings) in their 'computerized brain atlas'. The first automatic atlas registration method was introduced by [START_REF] Bajcsy | Multiresolution elastic matching[END_REF] to find the cortical and brain ventricle structures. Since then, other methods have been proposed (as in Thirion, 1995;[START_REF] Bro-Nielsen | Fast fluid registration of medical images[END_REF]Christensen et al., 1996b;[START_REF] Feldmar | Rigid, affine and locally affine registration of free-form surfaces[END_REF][START_REF] Szeliski | Matching 3-D anatomical surfaces with non-rigid deformations using octree-splines[END_REF]. • Morphometric study. Morphometry allows one to study covariances of biological shapes [START_REF] Bookstein | Morphometric Tools for Landmark Data[END_REF][START_REF] Bookstein | Landmark methods for forms without landmarks: morphometrics of group differences in outline shape[END_REF]. After the registration between the atlas and the patient data, a statistical analysis of landmark positions, or shape parameters of anatomical structures, is performed. The structures with coordinates, or parameters, which are outside of 'normal' statistical bounds are considered as 'abnormal'. A pathology diagnosis could then be inferred, for example, by using an expert system [START_REF] Suzuki | A supporting system for getting tomograms and screening with a computerized 3D brain atlas and a knowledge database[END_REF]. Such morphometric tools can also be used to characterize the evolution of an anatomical structure over time, or throughout history [START_REF] Dean | The Middle Pleistocene Homo Erectus/Homo Sapiens Transition: New Evidence from Space Curve Statistics[END_REF].

In fact, the study of morphometry opens up new atlas applications, in particular, in computer-aided diagnosis, e.g. characterization of Crouzon's disease (Cutting et al., 1995a), brain abnormality detection [START_REF] Thompson | High-resolution random mesh algorithms for creating a probabilistic 3D surface atlas of the human brain[END_REF][START_REF] Thompson | Detection, visualization and animation of abnormal anatomic structure with a deformable probabilistic brain atlas based on random vector field transformations[END_REF], sex differences in the morphology of the corpus callosum [START_REF] Davatzikos | A computerized approach for morphological analysis of the corpus callosum[END_REF] or computer-aided surgery, e.g. craniofacial operations (Cutting et al., 1995b).

Complexity of building a morphometric atlas

To build a computerized morphometric anatomical atlas, we have to address the following two major problems:

• Defining shape description parameters. Morphology results depend on the statistical study of shape. So we have to determine a set of parameters that characterize the shape of an anatomical structure that is quite complex. Thus, if we take an isosurface of a simple structure, e.g. the cerebral ventricles, extracted by the classical 'marching cubes' algorithm in a high-resolution MR image (voxel size of 1.0×1.0×1.5 mm 3 ), we obtain several tens of thousands of 3-D points to define the geometry. What we want is to compute an extremely condensed set of 'meaningful' shape parameters, say around ∼20-30, to obtain a more compact and easy to understand representation. Moreover, once these parameters are chosen, we can calculate their average and covariance which will be used in statistical tests. This can be done for points or frames [START_REF] Pennec | Validation of 3-D registration methods based on points and frames[END_REF] or for more abstract parameters such as vibration modes.

• Handling large amounts of data. To obtain meaningful statistical morphometric results, we have to process a large database of at least several dozen 3-D images, each one requiring several Mbytes. Currently, however, digital anatomical atlases are built by delineating manually anatomical structures in one medical 3-D image (Schiemann et al., 1996a): an anatomist uses a semiautomatic interactive segmentation tool to identify the voxels, very often slice by slice. This task takes too long and is too difficult to be generalized for a huge database. Moreover, it is not always possible to identify manually landmarks with a precision of one voxel in a huge 3-D medical image.

We conclude that only automatic tools could lead to the construction of morphometric anatomical atlases that both take into account the accuracy of new medical images and integrate the computation of quantitative parameters. These tools must integrate automatic segmentation of anatomical features, automatic non-rigid registration of these features between patients, automatic identification and statistical comparison of shape parameters and must be applicable to large databases of very high-resolution images.

Some related work

In [START_REF] Cutting | A spline-based approach for averaging three-dimensional curves and surfaces[END_REF]Cutting et al. ( , 1995a)), the authors propose building an average template of the skull composed of points, lines and surface patches based on a database of nine normal skulls segmented from CT scans. They use this average model to study the shape of skulls affected by Crouzon's disease. This work involves anatomists, morphometricians and surgeons. In particular, anatomy specialists are needed to define an a priori template of anatomically meaningful feature points and lines, and also to adapt this template to the patient data. Such manual intervention could limit or prevent the generalization of the technique to large databases, required to obtain a good statistical accuracy in shape variability studies. Moreover, the template is based on manually extracted points and lines that could appear too sparse (only 50 points and 100 lines) for a very precise shape parametrization. Also, the template is static and cannot be improved by adding points or lines which could evolve with the study of new data. Morphometric tools presented in this work were also applied in [START_REF] Boes | Generating a normalized geometric liver model using warping[END_REF] to build an average model of the liver based on six landmark points extracted in 15 CT images.

Our work can also be related to some research developed for brain study. In [START_REF] Royackkers | Feature extraction for cortical sulci identification[END_REF], a sulci statistical model is built from several MR images. The aim is to identify efficiently and robustly the superficial part of six major sulci in patient data. Here also, the structure of the atlas is fixed once by the user. In [START_REF] Mangin | From 3D magnetic resonance images to structural representations of the cortex topography using topology preserving deformations[END_REF], a high-level representation of the cortical topography is inferred from a brain MR image. This representation is very complete and integrates the topology and quantitative parameters such as length or depth. This model is intended to estimate statistically the inter-subject variability and to detect and recognize automatically the main cortical sulci.

Similar to these works, we would like to emphasize that building an atlas consists not only in computing the average parameters and standard deviations of features, but also in detecting and identifying in the images which features are common to all the subjects and will be used in the atlas.

Content

In this paper we describe a scheme for automatically building morphometric anatomical atlases. After an overview of the method in Section 2, we detail each stage in Sections 3-7. The scheme is fully automatic and has been tested on a database of six different skulls extracted from high-resolution CT scans. At the end of the process, we obtain an average skull model composed of common line features in their mean positions and their standard deviations. In Section 8, we present a sample morphometric study of a skull affected by a severe maxillary deformation. This study, in spite of its simplicity from a medical point of view, shows how an automatically built morphometric atlas could be useful for medical applications. Section 9 describes future work, in particular, the application of the scheme to other anatomical structures.

This paper must be considered as a long-term computer science research presentation to test a new concept and not as a direct medical application. In this spirit, we have given a greater importance to the presentation of a globally consistent scheme with working prototypes, being aware that each stage of this scheme could be improved.

THE ATLAS BUILDING SCHEME

Our scheme summarizes the method used by anatomists to draw up atlases: the study of different patients' structures allows one to identify which features appear visually common to all the data and in a stable position. This would correspond to the notion of 'biological homology'.

First, we need to collect a database of high-resolution 3-D medical images of different patients. In our skull example, we use six high-resolution CT scans a of dry skulls, without any artefacts and with a voxel size of 1.0 × 1.0 × 1.5 mm 3 .

We then apply a preprocessing stage to segment the anatomical structure we want to study. Various methods can be applied and one can refer to [START_REF] Ayache | Segmentation of complex three-dimensional medical objects: a challenge and a requirement for computerassisted surgery planning and performance[END_REF] for an overview of the methods, [START_REF] Mcinerney | Deformable models in medical image analysis: a survey[END_REF] about using deformable models or [START_REF] Kapur | Segmentation of brain tissue from magnetic resonance images[END_REF] for brain segmentation. From the segmented binary image, we extract the anatomical structure surface by using the 'marching cubes' algorithm [START_REF] Lorensen | Marching cubes: a high resolution 3D surface construction algorithm[END_REF].

In the case of the skull, a simple intensity thresholding gives a good segmentation of the bone as it appears very bright in CT scans. In Figure 1, we present the surfaces of the six skulls (A to F) which constitute our database. We notice a very important diversity in the skull orientation, size and shape.

The building scheme itself is composed of four stages (see Figure 2, left):

• Stage 1: Feature extraction. We extract some geometrical features automatically. We have to choose a feature type which combines a mathematical definition and an anatomical relevance. • Stage 2: Common feature identification. We find correspondences between the sets of features in different images by using a non-rigid registration algorithm. We then identify which features are common to all the data sets. These common feature subsets will form the structure of the atlas. These skulls were segmented from high-resolution CT scan images (acquired from two different devices, the voxel size is about 1 × 1 × 1.5 mm 3 ) by classic mathematical morphology and thresholding tools. We notice an important diversity in the skull size, orientation and shape.

• Stage 3: Average position. We average the common feature positions and then obtain the atlas mean geometry.

• Stage 4: Variability analysis. We analyse the variability of common feature positions with respect to the mean position and we compute some shape parameters. They describe and quantify, concisely and precisely, the shape of the common features and thus those of the anatomical structure.

When we want to study a patient (see Figure 2, right), we extract features from the 3-D image. We use the nonrigid registration algorithm to find correspondences between atlas and patient features (automatic labelling and normalized superimposition). We can then compare the shape parameters statistically. This leads to the detection of 'abnormal' shapes of the anatomical structure (shape analysis).

FEATURE EXTRACTION

Choice of feature

Once anatomical surfaces are extracted, one decomposes the surface into characteristic features. It can be surface patches, e.g. defined by their local shape [START_REF] Brady | Describing Surfaces[END_REF], line features, e.g. based on differential geometry [START_REF] Hosaka | Modeling of curves and surfaces in CAD/CAM[END_REF], or point features, e.g. 'extremal' points (Thirion, 1996a).

Surface features are difficult to handle because they integrate all the points of the structure. At the other extreme, point features give only sparse and non-robust information as there is no connectivity relationship. Line features are a very interesting compromise as they combine an important reduction of surface information with strong topological constraints as a line is an ordered list of points. In this paper we concentrate on the use of 'crest lines' introduced by [START_REF] Monga | From voxel to intrinsic surface features[END_REF] and developed by Thirion andGourdon (1995, 1996). Indeed they appear to be very good landmarks as they have been used successfully for rigid matching of 3-D medical images [START_REF] Ayache | Evaluating 3D registration of CT-scan images using crest lines[END_REF]. Moreover, they have a very strong anatomical meaning as we will see in the following.

Description of crest lines

Crest lines are defined by differential geometry parameters: let k 1 be the principal curvature with maximal curvature in absolute value and t 1 its associated principal direction, a point P belongs to a crest line when k 1 is maximal in the direction of t 1 which can be written as the zero crossings of the criterion e 1 = ∇ k 1 • t 1 (see Figure 3).

We find in Thirion andGourdon (1995, 1996) an original approach to computing crest lines on an isosurface defined by I (x, y, z) = I 0 , where I (x, y, z) is the intensity of the voxel localized at (x, y, z). It is based on the implicit representation of surfaces that leads to formulae which involve the differentials of the 3-D image up to third order: ∂ I /∂ x, ∂ 2 I /∂ x 2 , ∂ 3 I /∂ x 3 . They are calculated by Gaussian convolutions. Then the 'marching lines' algorithm follows the segments of crest lines that were extracted at each voxel by application of the two criteria e 1 = 0 and I = I 0 to create lines.

Due to their definition, the crest lines follow the salient lines of a surface. We can verify this in Figure 4 where crest lines of skull F emphasize the mandible, the orbits, the cheekbones or the temples and also, inside the cranium, the sphenoid and temporal bones as well as the foramen magnum. The automatic extraction gives 548 lines composed of 19 933 points.

Anatomical relevance of crest lines

Salient structures are also used by doctors as anatomical landmarks. For example, the crest lines definition is very close to the 'ridge lines' described in [START_REF] Bookstein | A proposal for the apprehension of curving cranofacial form in three dimensions[END_REF] and [START_REF] Cutting | Applications of computer graphics to the evaluation and treatment of major craniofacial malformations[END_REF]. In Figure 5, we display on the same skull the crest lines (in grey) and the ridge lines (in black) which were extracted semi-manually under the supervision of an anatomist [START_REF] Dean | Homology and the criteria for building deformable templates[END_REF]. The two sets of lines are very close, showing that crest lines would have a strong anatomical significance. Nevertheless, some crest lines appear noisy due to the discretization of the image (e.g. on the top of the skull in Figure 4). This problem can be partially solved by filtering the crest lines by a hysteresis thresholding on the value of the maximal curvature to avoid small curved curves (see Figure 6). Moreover, crest lines do not always correspond to the topology expected by anatomists; for example, the orbital crest lines are not closed. Only an a priori model could add this constraint.

Furthermore, we would also like to emphasize that the scheme itself checks the stability of the feature and its validity. In stage 2, we search for the common features and in stage 4, we compute statistical information. If no feature is common to all the database sets or if their variabilities are too large, this type of feature must be rejected. Otherwise the features can be considered as having good anatomical significance. This is very important because we can imagine testing features defined only by very complex mathematical formulae and seeing whether they could be relevant to characterizing anatomical structures.

Other feature lines

Other feature 3-D lines have been used in medical image processing:

• Geodesic lines. These lines are the shortest path between two points on a surface. In [START_REF] Cutting | A spline-based approach for averaging three-dimensional curves and surfaces[END_REF], geodesic lines are used to complement ridge lines to represent the surface of the anatomical structure with greater precision. Nevertheless, geodesic lines do not have a real anatomical meaning. Moreover, the optimization scheme for computing the geodesics is computationally very expensive. • Medial axis. In 2-D, these lines can be defined as the collection of all the centres of the circles which fit just inside the anatomical structure boundary [START_REF] Blum | A transformation for extracting new descriptors of shape[END_REF]. They define the local symmetry axis and indicate how a biological form is put together out of geometrically simpler pieces [START_REF] Bookstein | A proposal for the apprehension of curving cranofacial form in three dimensions[END_REF]. The generalization in 3-D gives medial surfaces but 3-D lines can be obtained by the intersection with another surface. For example, in [START_REF] Székely | Mapping the human cerebral cortex using 3D medial manifolds[END_REF] and [START_REF] Näf | 3D Voronoi skeletons and their usage for the characterization and recognition of 3D organ shape[END_REF], the medial surfaces of the brain data complement are intersected with the smoothed cortical surface which creates 3-D lines following the sulci. Nevertheless, computing medial surfaces is very expensive in terms of memory and time. • Junction lines. The skeletonization of the anatomical structure by discrete topology methods gives surfaces.

Then, the junction lines can be detected by a topological classification algorithm [START_REF] Malandain | Topological segmentation of discrete surfaces[END_REF]. Such lines have been used to define lines following the sulci on the cortical surface [START_REF] Vidal | Squelettes et Outils de Topologie Discrète[END_REF].

Description of stage 1:

• Extract features in A, B, C, D, E, F.

FEATURE REGISTRATION

Introduction

Given two sets of lines A and B extracted from two different patient images (see Figure 7), we want a twofold result:

• Line to line correspondence. Which line L i of A corresponds to which line L j of B? This allows us to find the common lines to all the sets in stage 2. • Point to point correspondence. Which point of A corresponds to which point of B? We need to know the corresponding points over the different sets to compute the average positions of lines in stage 3 and to analyse their variabilities in stage 4.

As a matter of fact, little work has been done on the registration of 3-D curves. In [START_REF] Bastuscheck | Object recognition by three-dimensional curve matching[END_REF] and [START_REF] Schwartz | Identification of partially obscured objects in two and three dimensions by matching noisy characteristic curves[END_REF], the rigid matching algorithm uses fast Fourier transforms to determine the least-squares difference between sequences of points sampled at equal intervals along two piecewise linear approximations of 3-D curves. [START_REF] Mokhtarian | Multi-scale, torsion-based shape representations for space curves[END_REF] proposes to model a 3-D line by its torsion profile for different scales, the extrema of which are then matched. [START_REF] Guéziec | Smoothing and matching of 3Dspace curves[END_REF] improve on the method described in [START_REF] Kishon | 3-D curve matching using splines[END_REF]: lines are indexed according to their differential characteristics computed by an approximation by B-splines. With hash tables, it is then fast to retrieve a point with given differential parameters and to test the accuracy of a rigid transformation. [START_REF] Pajdla | Matching of 3-D curves using semi-differential invariants[END_REF] use a semi-differential invariant description requiring only first derivatives and one reference point. All these methods have been developed for rigid matching and cannot be generalized to the non-rigid case as they are based on using Euclidean invariants.

Moreover, the development of a registration algorithm must deal the complexity of the data. The sets of lines are:

• very different in orientation, number of lines, shape, topology and discretization. • very dense as they are constituted of several hundred lines and tens of thousands of points.

The registration algorithm 4.2.1. General overview

To overcome the difficulties of the registration task, we propose to use a heuristic algorithm based on an iterative scheme. Such ideas are not new: for example, [START_REF] Burr | A dynamic model for image registration[END_REF] introduced an iterative technique to update gradually local registration of two different images, where each feature at one location influences matching decisions made at other locations. This principle is developed in the 'iterative closest-point' algorithm introduced by [START_REF] Besl | A method for registration of 3-D shapes[END_REF] and concurrently by [START_REF] Zhang | Iterative point matching for registration of freeform curves and surfaces[END_REF]. It consists of iteratively applying rigid transformations, based on a local point matching with its closest neighbour, to the set A in order to superimpose it on the set B. Such a method has already been generalized to non-rigid registration of surfaces in [START_REF] Feldmar | Rigid, affine and locally affine registration of free-form surfaces[END_REF].

We have adapted the original algorithm to our problem by:

• Generalizing it to non-rigid transformations. We have modelled deformations between anatomical structures by affine, polynomial and spline functions. • Taking into account constraints inferred by lines. The order of points along the lines constrains point correspondences. Thus, we discard inconsistent point matchings and compute registration parameters which define line correspondences.

In the following, we are going to review each step of the adapted ICP algorithm (see Figure 8).

Point matching

At each iteration, all the points of A's lines are linked with their closest neighbour in B with respect to the Euclidean distance. This is possible by using an efficient data structure [START_REF] Zhang | Iterative point matching for registration of freeform curves and surfaces[END_REF]) called a 'k-d tree' [START_REF] Preparata | Computational Geometryan Introduction[END_REF] or by precomputing a distance map [START_REF] Cuchet | Registration in neurosurgery and neuroradiotherapy applications[END_REF] that integrates the coordinates of the closest points. The Euclidean distance could be extended to include differential parameters (e.g. curvature) as proposed in [START_REF] Feldmar | Rigid, affine and locally affine registration of free-form surfaces[END_REF]. This preliminary simple matching gives a first list of point pairs, M 1 . Notice that the closest-point process is not bijective: each point of A has one and only one correspondent on B, whereas some points of B may have either no correspondent on A or more than one.

Line matching

If we want to estimate whether two lines L i ∈ A and L j ∈ B are registered, we need to compute the proportion p j i of points of L i which are matched with points of L j and the proportion p j i of points of L j which are matched with points of L i . If p j i or p j i are larger than a given threshold, for example 50% (it then corresponds to more than half of the line points being matched), we can conclude that the two lines are registered.

However, computing the line registration parameters p j i and p j i is not simple due to the non-bijectivity of points matching as we can see in Figure 9. On the left-hand side, we compute p j i = 100% and p j i = 40% (5/5 points of L i are registered with 5/13 points of L j ). Moreover the matchings 2 and 3 join a portion of L i to a portion of L j which has been already registered with another portion of L i by the matchings 1 and 4. We call this event a cross-matching. On the righthand side, the important difference of sampling along L i and L j yields that p j i = 100% and p j i = 67% (10/10 points of L i are registered with only 4/6 points of L j ). We also notice numerous multiple matchings (e.g. matchings 3, 4 and 5). To address this problem, we introduce two additional constraints:

• Injectivity constraint. Each point of B is linked to at most one point of A. • Monotonicity constraint. The ordering of corresponding points on L i and L j must be the same. In particular, this implies that the same portion of L j cannot be matched to two different portions of L i . Such a condition has also been described in [START_REF] Geiger | Matching elastic contours[END_REF].

To impose these constraints, we sort the matching between points of L i and of L j according to their distance. The closer two matched points are, the more likely the correspondence is. We begin at the most likely matched point P 0 of L i and we follow the line in the two directions. When we meet another point, P k (that is the kth one on L i after P 0 ), we look at its correspondent corres(P k ). If corres(P k ) does not belong to L j (we should not forget that we deal with hundreds of lines in B), we stop the propagation in this direction. Otherwise, if corres(P k ) (that belongs to L j ) has already been marked, it means that we are creating a cross or a multiple matching, which we prevent by discarding the current matching and stopping the propagation. If corres(P k ) has not been marked, we keep the matching (P k , corres(P k )), we mark all the points of the portion of L j , ]corres(P k-1 ), corres(P k )] and we continue.

When the process is terminated, we begin again with the most likely matched point of L i that has not already been met. In this way, we obtain consistent point correspondences. The algorithm is very fast as the complexity is proportional to the number of points of L i . To compute correctly the registration parameters between L i and L j , we have to follow the lines L j and to consider as virtually matched, the points which are localized between two points really matched with points of L i . Thus, in Figure 10, we now have on the left-hand side p j i = 60% and p j i = 77% (3/5 points of L i are registered with 10/13 points of L j ) and on the right-hand side p j i = 100% and p j i = 100% (10/10 points of L i are registered with 6/6 points of L j ). Based on this new matching, we obtain a second list of matched point pairs, M 2 , which is a consistent subset of M 1 .

Moreover, we can use the parameters p j i and p i j to make the registration algorithm more reliable. The idea is to consider that only matched points which belong to lines which are registered with a given threshold ( p j i > threshold or p j i > threshold) and to take only them into account in the following. We can then raise this threshold as iterations are performed, to improve only registration of lines already quite well registered. With this selection, we obtain a new list of matched point pairs, M 3 which can be considered as a reliable subset of M 2 .

Transformation computation

Based on the list of matched points M 3 = (P k , P k ), we can compute a transformation T of a given type, by minimizing the least-squares criterion:

P k ∈M 3 d 2 (T (P k ), P k )
where d is the Euclidean distance between two 3-D points.

What types of transformation can we use? In fact, we have to address the four following problems:

• Physical modelling. Others have tried to define some transformation classes to model inter-patient deformations as, for example, piecewise affine functions in Talairach and Tournoux (1988), but they remain much too simple to be very accurate. • Computation complexity. A transformation that can be written mathematically as a linear function of its coefficients is easy to compute by a least-squares method. • Regularity control. With a very general class of transformations, we can deform anything into anything. So, we need to have the possibility of controlling regularity by introducing constraints. • Topology conservation. The topology of the anatomical structure must be preserved by the transformation. In particular, we have to avoid self-intersection of the surface [START_REF] Christensen | Topological properties of smooth anatomic maps[END_REF].

According to these four criteria, we choose to use the following transformation types (see Figure 11):

• Rigid transformations at the beginning of the registration to align the two sets of lines. The least-squares computation can be performed by several methods [START_REF] Arun | Least-squares fitting of two 3-D point sets[END_REF][START_REF] Horn | Closed form solutions of absolute orientation using unit quaternions[END_REF]. • Affine transformations to retrieve the scalar differences between the two sets.

• Spline transformations to model more local and complex deformations. They have been used widely in 3-D medical image processing (e.g. [START_REF] Declerck | Automatic registration and alignment on a template of cardiac stess & rest SPECT images[END_REF][START_REF] Szeliski | Matching 3-D anatomical surfaces with non-rigid deformations using octree-splines[END_REF]. [START_REF] Bookstein | Principal warps: thin-plate splines and the decomposition of deformations[END_REF][START_REF] Bookstein | Landmark methods for forms without landmarks: morphometrics of group differences in outline shape[END_REF] proposes using a thin-plate spline interpolating function. Nevertheless, interpolation is relevant when the matched points of M 3 are totally reliable and distributed regularly (for example, with a few points being located manually).

In our case, these points are not totally reliable due to possible mismatches of the registration algorithm and are sparse in a few compact areas as they belong to lines. So, we prefer to use a spline approximation function which is regular enough to minimize the influence of an erroneous matched point [START_REF] Declerck | Automatic retrieval of anatomical structures in 3D medical images[END_REF]. The coordinate functions of T , (u, v, w), are then computed by a 3-D tensor product of B-spline basis functions. For instance, for u:

u(x, y, z) = n x -1 i=0 n y -1 j=0 n z -1 k=0 α i jk B x i,K (x) B y j,K (y) B z k,K (z)
with the following notation: n x is the number of control points in the x direction, which sets the accuracy of the approximation (eight in our experiments). α i jk is the 3-D mesh of the control points abscissae. These parameters define the transformation. B x i,K is the ith B-spline basis function; its order is K . The B x i,K generate the vector space of piecewise K -order polynomials. u is then a piecewise K th degree polynomial in each variable x, y and z. For their regularity properties, we choose cubic B-splines in our experiments (K = 3) with a regular grid of 3-D knots.

For a given number of control points and a set of B-spline basis functions, u is completely defined by the α i jk . They are calculated by minimizing a criterion computed with the set M 3 of matched points. In fact, the criterion splits into two parts: J (u) = J position (u) + J smooth (u).

• Position term. For each data point P k , u(P k ) must be as close as possible to P k . We choose a least-squares criterion:

J x position (u) = N k=1 u(x k , y k , z k ) -x k 2 .
Similar equations apply to J y position (v) and J z position (w) with y k and z k respectively.

• Smoothing term. B-splines have intrinsic smoothness properties, but these may be insufficient. We choose a second-order Tikhonov stabilizer: it measures how far from an affine transformation the deformation is: where ρ s is a weight coefficient that tunes the importance of smoothing.

J x smooth (u) = ρ s R 3 u 2 xx + u 2 yy + u 2 zz + 2u 2
J x is a positive quadratic function of the α i jk variables. To find the coefficients which minimize J x , we derive its expression with respect to all the α i jk : this yields n x × n y × n z linear equations. Rearranging these equations, we get a sparse, symmetric and positive linear system. We solve three systems (one for each coordinate) to completely estimate T .

We also have tried to use quadratic transformations which model some anatomical deformations of the brain according to [START_REF] Greitz | A computerized brain atlas: construction, anatomical content and some applications[END_REF], but spline functions appeared to be more general and more convenient to control the regularity.

Transformation implementation and iteration

The transformation T is then applied to A, bringing A closer to B and we iterate the process by modifying at each step two parameters: • the threshold on the registration parameter threshold used in the 'lines matching' step. • the type of the transformation T .

Presently, we use a constant iteration scheme: 30 iterations where the variable threshold is incremented from 0% to 50% and which consists of 10 rigid transformations, then 10 affine, and finally, 10 spline transformations. For spline functions, the smoothing parameter ρ s decreases from 10.0 (very rigid) to 1.0 (very deformable). We plan to use an adaptive scheme that automatically modifies the number of iterations as well as the number of each transformation type. But we would need to find a good criterion to evaluate the accuracy of the registration at each step (maybe based on the evolution of the mean-squared distance between matched points). At the end of the iterations, we obtain two results: the registered lines by thresholding the registration parameters p j i and p i j to the value 50% and the matched points belonging to M 3 .

In Figures 12131415, we show the result of the registration at different iterations (original position, after rigid, rigid+affine, and rigid + affine + spline) for two simplified sets of lines: we can see the registered lines and the matched points, in their deformed positions on the left, and in their original positions on the right.

Despite its simplicity and generality, this algorithm appears to be quite robust and is quite insensitive to discretization differences [for more details, see [START_REF] Subsol | Construction Automatique d'atlas Anatomiques Morphométriques à Partir d'images Médicales Tridimensionnelles[END_REF]]. Nevertheless, it is not symmetric in the sense that the registration of A towards B gives different results (small in general) from the ones obtained with the registration of B towards A.

A real example

In the following example, we register C towards B. The whole registration takes about 10 min on a DEC-Alpha workstation, 166 MHz. Notice that 80% of the CPU time is required to find the closest point. This step is already optimized with a k-d tree structure and an algorithm of complexity O(n 2/3 ) in the worst case, if n is the number of points stored in the tree [START_REF] Preparata | Computational Geometryan Introduction[END_REF]. In Figure 16 we can see on the left that the superimposition of the two sets is relatively accurate and we can conclude that line and point correspondences are reasonably correct. On the right, the two sets are in their original Table 1. The number of matched points and statistical parameters about the distribution of distances between matched points at the beginning of the registration (Begin), after rigid (Rigid ), rigid + affine (Affine), and rigid + affine + spline (Spline). These last values are to be compared with the set diameter which is around 200.0 mm. position which shows the extent and the complexity of the deformation between inter-patient anatomical structures.

Nb

In Table 1, we show the number of matched points and some statistical parameters about the distribution of distances between matched points at the beginning of the registration (Begin), after rigid (Rigid), rigid + affine (Affine) and rigid + affine + spline (Spline). All the values must be compared with the diameter of the anatomical structure that is about 200.0 mm. Whereas the number of matched points increases very fast after rigid (Begin -Rigid = +22%) and spline transformations (Affine -Spline = +13%, Begin -Spline = +36%), it stays stable between the rigid and affine transformations (Rigid -Affine = -1%). Nevertheless, affine transformations are very useful as they diminish the mean distance a lot (Rigid -Affine = -9%). During the whole process, the mean distance has diminished by 55% while the matched points increased by 36%. The standard deviation also decreases by 42%: more points are closer, which is confirmed by the evolution of the median distance: -61%.

Anatomical relevance of the registration

The previous figures give us only an evaluation of the quality of the superimposition of the two sets of lines. We assume then that a 'good' superimposition involves accurate registration results. But what about the real anatomical relevance of the registration? Thirion et al. (1996) present a technique to cross-validate different non-rigid matching techniques. The overall aim is to determine whether different methods, developed independently, give mutually coherent image superimposition results. In particular, a study was performed to compare three deformable techniques to superimpose skull images of different patients: the first method is based on ridge lines relying on the manual identification of anthropometric landmarks (see Section 3) [START_REF] Cutting | A spline-based approach for averaging three-dimensional curves and surfaces[END_REF][START_REF] Dean | Homology and the criteria for building deformable templates[END_REF], the second technique is the one described in the present paper and the last one is based on intensity (Thirion, 1995). The conclusion is that the three methods give mutually coherent results, with an average difference for feature location of 3-4 mm where the skull is highly curved. In smoother places, this average accuracy is reduced to 6-9 mm. It proves, in particular, that our registration technique is consistent with the superimposition based on ridge lines supervised by anatomists. Moreover, it confirms that crest lines are good landmarks as results are more stable around anatomically salient lines.

COMMON FEATURE IDENTIFICATION

The registration graph

With the registration algorithm, we can find the line correspondences between two sets, for example, A and B. More generally, we can find line correspondences between all the sets of the database (A → C, C → B etc.). A line that will be common to all the patients of the database will be a line that has a correspondent in all the data sets. So to find such lines, we need to build an overall representation of the correspondences between all the lines of all the sets.

We choose to use a graph representation where nodes represent a line of a set and the links mean 'is registered with' (see Figure 17). Notice that the links are oriented as the registration algorithm is not symmetric. To overcome this asymmetry, we perform the registration between two sets X and Y in both directions (X → Y and Y → X ). We then keep only the symmetric links (solid arrows) which we assume to correspond to robust line matchings.

Finding common lines

Sets of corresponding lines can be modelled as the connected components of the registration graph which are easy to compute by a classic propagation algorithm [for example, described in [START_REF] Cormen | Introduction to Algorithms[END_REF]] (see Figure 18). Since we want lines which are common to all the data, we keep only the connected subgraphs which contain at least one line of each set. In Figure 18, subgraph 3 is not taken into account as it does not contain a line from set C. All these subsets of common lines form the structure of the atlas.

In fact, to reduce the complexity of the atlas building, we do not always perform all the registrations between all the sets, an operation of complexity O(n 2 ) where n is the number of sets in the database. We can use a circular permutation that reduces the complexity to O(n):

A → B, B → C . . . Z → A.
Nevertheless, performing n registrations instead of n 2 reduces the number of links and may diminish the number of common line subsets. So in this case, we can then accept the subgraphs which include lines not from all the data sets but from a high proportion, for example, 80%.

Application to automatic labelling

If we are able to associate a label (for example, manually) to a common line of a data set (for example, 'mandible' associated to the line 123 of set A in the graph of Figure 18), we can propagate it to the subgraph, and then to all the corresponding common lines. This application is very useful for visualizing the results of registration and common feature identification.

Thus, by building the registration graph for the six skulls of the database, we find 63 subsets of common lines. We have represented in Figure 19 the lines of skulls B and C which are common. As we have labelled skull A, we can propagate the labelling to skulls B and C. Thus, we recognize in highlighted lines, the mandible (LMB+RMB), the nose (NOS), the orbits (LOR and ROR), the cheekbones (LCB and RCB) the temples (LTP and RTP) and the foramen magnum (FOR).

We also notice that the two parts (left and right) of the mandible have been merged as the same common line labelled RMB + LMB. This is due to the fact that in one data set, the mandible was in one part which has been registered to the two parts of the mandible of B and C, clustering then the two halves in the same common subset.

Automatic labelling can be generalized to automatic extraction of a part of the anatomical structure. Patches of the skull isosurface with points which are within a given distance of the lines identified by automatic labelling can be extracted. Thus, in Figure 20, skull B is automatically decomposed: we can identify the left and right orbits, the nose, the mandible, the left sphenoid bone and the foramen magnum.

FEATURE AVERAGE

Introduction

In this stage we wish to find the average positions of the features constituting the atlas, i.e. to average the sets of 3-D lines defining each common subset. We choose a common line L i of one data set (e.g. A). Thanks to the results of the previous stage 'common feature identification', we know the corresponding line(s) in the other data sets which we call L i (B) . . . L i (F). We can then compute the correspondences between the points of L i and those of L i (B) . . . L i (F), perform the average of corresponding points and reconstruct an average line. Nevertheless, in order to average the positions, we need to align all the data within the same frame.

Aligning the data in a reference frame

As a reference frame, we can choose, for example, the frame given by one data set (e.g. A). With the registration algorithm, we find the list of matched points between the reference data set and the other data sets. Then, we can easily compute by a least-squares criterion the rigid transformations to align the entire data set in the reference frame. But, as emphasized by [START_REF] David | Déformations ontogénétiques et évolutives des organismes: l'approche par la méthode des points homologues[END_REF], in ontogenetic and evolutive shape transformation studies we should not take into account differences of position, orientation and size, since these cannot be considered as true morphological differences. So, we compute by a least-squares criterion not only the global rigid transformation but a similarity that is the composition of a translation, a rotation and an isotropic scaling. After applying the similarity transformations to the data sets, all the subsets of common lines are in the same referenced frame and the residual deformations between lines are really meaningful morphometric differences (see Figure 21).

Finding the point correspondences between common lines

With the registration algorithm, we can find the list of matched points between the common line L i of A and the corresponding common lines L i (A) . . . L i (F) which are now aligned in the reference frame.

Nevertheless, some points of L i do not have any correspondent point on the corresponding lines. This is due to the line constraints (injectivity and monotonicity) implemented in the registration algorithms which discard some point matching. We notice this in Figure 22 where some points of the bottom of the mandible are not linked with points of L i (C). In order to find the whole deformation of L i , we have to find the missing correspondences. For each unmatched point P of L i , we compute the position of its potential correspondent point by a linear interpolation given by

--→ P P = d 1 --→ P 2 P 2 + d 2 --→ P 1 P 1 d 1 + d 2
where P 1 and P 2 are the two closest neighbours of P which have correspondent points denoted as P 1 and P 2 and d 1 and d 2 are the distances along L i between P and P 1 and P and P 2 .

After interpolation of the missing correspondent points, we obtain the deformation D i between L i and the corresponding common lines which are the vectors field given by the matched points (P i , P i ) (see Figure 23).

Smoothing point correspondences 6.4.1. Presentation of the problem

The deformation D i may still be quite irregular if the corresponding common lines L i , L i (B) . . . L i (F) are not smooth. This roughness is mainly due to the image discretization. If we average the D i directly, we may obtain a very irregular average deformation field D avg . To avoid this problem, we propose to filter the deformations D i by a low-pass filter, assuming that high-frequency oscillations are anatomically meaningless. To implement this filter, we propose to use modal analysis introduced in image processing in [START_REF] Pentland | Closed-form solutions for physically based shape modeling and recognition[END_REF] and developed by [START_REF] Nastar | Frequency-based nonrigid motion analysis: application to four dimensional medical images[END_REF].

A filtering method: modal analysis

Given a deformation field D X ,Y i composed of the n 3-D vectors (D X ,Y i [0], D X ,Y i [1] . . . D X ,Y i [n -1]
), modal analysis proposes to express the field in a modal basis using the following formulae [START_REF] Nastar | Frequency-based nonrigid motion analysis: application to four dimensional medical images[END_REF] (for the x coordinate):

d X ,Y i [k] x = n-1 p=0 D X ,Y i [ p] x • φ p [k]
where

φ p [k] = cos( pπ(2k+1)/2n) n-1 j=0 cos 2 ( pπ(2 j+1)/2n) -1/2 . Reciprocally, we have D X ,Y i [k] x = n-1 p=0 d X ,Y i [ p] x • φ p [k]. The n parameters d X ,Y i [k] x (respectively d X ,Y i [k]
[k] = 1 √ n d X,Y i [0] x = 1 √ n n-1 p=0 D X ,Y i [ p] x .
The other modes correspond to deformations of increasing complexity that leave the centre of line mass fixed. In Figure 24, we can see the effect of modes 1, 2 and 3 applied successively to the same mandibular line of A. The larger the mode number is, the more complex the deformation is.

Smoothing the deformation field

What is particularly interesting with modal analysis, as in Fourier analysis, is that we can approximate a deformation by taking into account only the first modes. Truncating the spectrum allows one to discard high-frequency deformations. Of course, the notion of frequency is only meaningful if the distance between points of the lines is constant, i.e. the points are with a constant amplitude to the mandibular line A (in grey, the original line; in black, the deformed line). The larger the mode number, the more complex the deformation.

uniformly sampled along the lines. As crest lines extraction does not verify this assumption, before modal analysis, we move the points along all the data lines in order to make their distances constant with a very simple algorithm described in [START_REF] Subsol | Construction Automatique d'atlas Anatomiques Morphométriques à Partir d'images Médicales Tridimensionnelles[END_REF]. A more sophisticated method, based on an approximation using spline curves, can be found in [START_REF] Guéziec | Smoothing and matching of 3Dspace curves[END_REF].

As we have noticed in Figure 24, the mode p introduces p sinusoids (in fact, p sinusoids for each coordinate axis), and so, can be associated with a wavelength of n/ p points. In our example, crest line segments are extracted in voxels of size 1.0×1.0×1.5 mm 3 . So we can assume that the segment length is at most around 1 mm. If we want to study details on a scale of 1 cm, we have to take into account deformations involving up to 10 points, with a wavelength of 10 points. This leads to the equation n/ p = 10 or p = n/10. Thus, in the following, we will keep the first 10% of the modes of the deformations spectra.

Computing the average lines

For each deformation field D A,X i , we compute the modal spectrum S A,X i that we truncate to the first 10% of the modes. Then we perform the average of the truncated spectra by averaging each amplitude independently. We integrate in this average the spectrum D A,A i that is null in order to take into account the influence of A in the atlas building. Based on the average spectrum S avg i we reconstruct the average deformation D avg i that we apply to the line L i to find the average common line L avg i .

In Figure 25, we present the average of common lines corresponding to the left orbit. Left, we can see the lines in their original positions. Right, we have discarded the differences of position, orientation and size by aligning all the data in the reference frame and we visualize the average line in black.

In Figure 26, we display some common lines with their average in black. In Figure 27, we present the set of all average common lines L avg i that constitute the geometry of the atlas. Notice their lateral symmetry.

Nevertheless, we could believe that the choice of the reference frame (A in our example) may influence the result of the averaging process. To study the effect, we have built an atlas successively from the six reference sets A . . . F. Then, we have aligned all the atlases in the same frame (the A one) by applying similitudes. In Figure 28, we notice how the six sets of average common lines are very close, which proves that the choice of the reference set has very little effect.

Obtaining an average surface

By registering the lines of the atlas and of the reference set A, we obtain a list of matched points. We have seen that we can compute a space spline transformation by a least-squares criterion that approximates these matched points [START_REF] Declerck | Automatic retrieval of anatomical structures in 3D medical images[END_REF]. If we apply this transformation to the surface of the reference set A, we obtain a surface representation of the atlas, based on average lines, which is presented in Figure 29.

We notice three things:

• The atlas is symmetric, whereas the averaging process is independent for each subset of common lines. Thus, it shows that the algorithm is consistent.

• Whereas the reference set A was very dolichocephalic (longer than wider head with a narrow face), the five other skulls are more brachycephalic (rounder heads).

The atlas is more brachycephalic, proving that the averaging process has correctly taken into account the data of the other skulls.

• Our atlas is visually very similar to the one presented in [START_REF] Cutting | A spline-based approach for averaging three-dimensional curves and surfaces[END_REF], which was created under the supervision of an anatomist.

FEATURE DEFORMATION ANALYSIS

Some previous work

In the previous stage, we have computed the average position of each common feature. We now have to estimate the variabilities of their shape.

Shape analysis was first based on the study of relative parameters such as the distance between two points or the angles between three points [START_REF] Abbot | Application and Comparison of Techniques for Threedimensional Analysis of Craniofacial Anomalies[END_REF]. Then a 'new morphometry' appeared a decade ago. According to [START_REF] Rohlf | A revolution in morphometrics[END_REF], it can be defined by the introduction of a tridimensional function fitting point relations that allows one to find the most meaningful parameters to analyse a shape and to build a taxonomy. The most famous method is based on thin-plate splines and is described in [START_REF] Bookstein | Principal warps: thin-plate splines and the decomposition of deformations[END_REF][START_REF] Bookstein | Landmark methods for forms without landmarks: morphometrics of group differences in outline shape[END_REF].

In our case, the differences between lines are modelled as the deformation fields D X ,Y i . We have to decompose them into meaningful principal deformations or modes. [START_REF] Martin | Shape analysis of brain structures using physical and experimental modes[END_REF] propose three main categories of decomposition:

• Predetermined modes. They are defined by the user as shearing, bending, tapering or pinching introduced in [START_REF] Barr | Global and local deformations of solid primitives[END_REF]. These modes have a concrete significance but they are limited and can only describe a rather simple deformation for simple structures. As anatomical structure deformations are very complex, we cannot use them. • Mathematical modes. Their formulae are only based on the geometry and the topology of the structure. We find in this category, Fourier decomposition [START_REF] Renaud | Fourier analysis applied to Stephanomys (Rodentia, Muridae) molars: nonprogressive evolutionary pattern in a gradual lineage[END_REF], the thin-plate splines method and the modal analysis we have previously used. The problem is that the modes have no anatomical or experimental significance. So they can be quite inappropriate for defining biological shapes with a few parameters. Moreover, when we truncate modal spectra, we discard high-frequency deformations that could be very important anatomically. • Experimental modes. In contrast to mathematical modes, their definition is based on the set of all the deformations D X ,Y i , which gives them an experimental validity. In particular, the principal components analysis used in [START_REF] Hill | Model-based interpretation of 3D medical images[END_REF], [START_REF] Cootes | The use of active shape models for locating structures in medical images[END_REF], [START_REF] Martin | Characterization of Neuropathological Shape Deformations[END_REF] and [START_REF] Székely | Segmentation of 2-D and 3-D objects from MRI volume data using constrained elastic deformations of flexible Fourier contour and surface models[END_REF] allows one to find a basis of modes, the importance of which can be perfectly quantified and sorted. Nevertheless, this method needs a very important training set (several tens of elements) in order to obtain a number of modes (for N elements, we can compute N -1 modes) that allows one to describe precisely a deformation that does not belong to the training set.

In our case, we have at the moment, only a few samples in the database, so we have preferred to use mathematical modes and we choose modal analysis in consistency with the previous stage. Moreover, modal analysis has already been used to study anatomical structures such as the mitral valve in [START_REF] Nastar | Frequency-based nonrigid motion analysis: application to four dimensional medical images[END_REF]. So, we have defined a very simple shape distance based on modal amplitudes.

A simple shape distance

We replace A by the lines of the atlas that we will call ATLAS in the following and we apply the same procedure as in the previous stage: we register the common lines L i of the atlas and those of the data in order to align all of them in the frame of ATLAS and to find the deformation fields D ATLAS,A i , D ATLAS,B i . . . D ATLAS,F i . We then deduce the spectra of these deformations and for each mode j, we compute the mean amplitude for the three axes, di [ j] x , di [ j] y , di [ j] z (which are very near to zero due to the averaging stage) and the associated standard deviations, σ i

[ j] x , σ i [ j] y and σ i [ j] z .
These statistical parameters allow us to define a distance between the shapes of two lines. Given a set of lines, X , we align it according to the frame of ATLAS by registering the lines and by computing with the leastsquares criterion a similitude that superimposes the best one. In this frame we then obtain the deformation D ATLAS,X i . We compute the modal spectrum that gives the amplitudes

d ATLAS,X i [0], d ATLAS,X i [1] . . . d ATLAS,X i [n -1] for the ith line.
For each mode, we can compare the amplitudes of the deformation towards X with respect to the amplitudes of those towards the different data by the amplitude distance defined (for the x-axis):

dist amp (D ATLAS,X i )[ j] x = d ATLAS,X i [ j] x -di [ j] x σ i [ j] x .
Large values of dist amp (D ATLAS,X i ) [ j] x allow us to find which modes characterize 'abnormal' deformations. If we assume that the distribution of the amplitude is a Gaussian law, we can associate the amplitude distance to a normalized centred Gaussian law. A value of dist amp (D ATLAS i , X ) [ j] x larger than 2.0 then indicates an 'abnormality' probability of 95%. This assumes that the modes are uncorrelated and can be studied independently. In fact, anatomical deformations are very complex and must be modelled by the combination of several modes, mixing the three axes. Nevertheless, we content ourselves with using this very simple amplitude distance that gives encouraging preliminary results.

Of course, we use this distance only for the 10% first modes of the deformation which are the only meaningful ones after the smoothing applied in the previous stage. In practice, we limit the study to the very first modes (4 or 5).

A CRANIOFACIAL APPLICATION: STUDY OF MAXILLARY DEFORMATION

In this section, we use the atlas in order to study a skull a S affected by a significant maxillary hypoplasia (see Figure 30).

Automatic labelling

By registering the crest lines of the atlas with those of S, we are able to label the latter, identifying the mandibular line. By taking points of the surface which are close to these lines, we automatically extract the mandible of S (see Figure 31).

Normalized registration

With the pairs of matched points found by the registration between the atlas and S, we compute the rigid and homothetic a Data from the Naturhistorisches Museum in Vienna. transformations which best superimpose (in a least-squares sense) the two skulls. In this way, we are able to contrast S and the atlas and emphasize the deformations of the mandible which appears laterally 'too wide' and vertically 'stretched'. One could imagine the potential use of such a dual display to plan a surgical procedure.

Quantitative shape analysis

Let us now analyse quantitatively the deformations of the mandible. With the method described in the previous section, we compute the first five modes of the S mandible deformations according to the average position given by the atlas. Figure 34. The automatic comparison between S and the atlas allows us to reposition the mandible in order to obtain a more 'normal' shape. We notice that a rotation and translation make the jaw less curved (left, the moved mandible is in mesh representation; middle, the original skull; right, after the proposed repositioning.

Then, for each of them, we compute the amplitude distance which are listed in Table 2. When the value is larger than 2.0 that corresponds to a variation of more than two standard deviations around the average value, the mode is considered as 'abnormal'. This is the case for the first and second x-mode, and the fourth z-mode. These corresponding deformations could be considered as typical of the mandibular abnormality (see Figure 32).

Application to computer-aided diagnosis

In order to visualize these abnormal basic deformations, we deform the atlas mandible according to these three modes (with an amplitude multiplied by 3 to exaggerate the corresponding deformations). A qualitative analysis of Figure 33 shows:

• the first x-mode can be associated to the breadth of the mandible.

• the effect of the second x-mode is not symmetrical. In fact, the mandible appears slightly skew and this mode can be associated to vertical twist.

• the fourth z-mode can be associated to lateral curvature of the mandible.

So, with only three automatically detected parameters, we could define and estimate the severity of the 3-D deformation of the mandible.

By registering the mandible crest line of S and the atlas, we are able to move the S mandible towards a more 'normal' postion. In this way, we could simulate the craniofacial surgery procedure of cutting some pieces of bone [START_REF] Marchac | New aspects of craniofacial surgery[END_REF] in order to reposition them (see Figure 34).

CONCLUSION

In this paper we have described a scheme for building automatically a morphometric anatomical atlas from 3-D medical images. We have shown how such an atlas could be used by presenting a number of preliminary experiments on a skull affected by a severe mandible deformation.

We would like to stress the fact that the described scheme is general and can be applied to other anatomical structures. We have already created an atlas of crest lines of the brain based on 10 different patient MR images and we have used it to make a specific study of the deformations of the cerebral ventricles Subsol et al. (1996a[START_REF] Subsol | Automatic analysis of cerebral atrophy[END_REF].

In future work, we plan to improve each stage of the scheme:

• Feature extraction. The extraction of crest lines could be made more robust and precise by using a multi-scale algorithm [START_REF] Fidrich | Following feature lines across scale[END_REF]. We have also begun to use 3-D skeleton lines extracted by mathematical morphology operators in order to characterize the sulci in a brain atlas. Moreover, we could mix line and point features as proposed in Thirion (1996b) to integrate accurate point features with robust line features.

• Common feature identification. At present, the common features are entire lines. We plan to detect common portions of lines. Moreover, we do not use all the information given by the matching graph that would allow us to obtain more robust common features. For example, two subgraphs of common features which contain many lines and which are joined by only one link could be considered distinct, the link being assumed to be an artefact.

• Average position. We plan to compare our average method with one developed by morphometricians [START_REF] Dean | The Middle Pleistocene Homo Erectus/Homo Sapiens Transition: New Evidence from Space Curve Statistics[END_REF]. Moreover, we wish to compare the filtering results of modal and Fourier analysis.

• Variability analysis. This is surely the weakest stage because it is based on the assumption of non-correlation of modes. We could use multivariate hypothesis testing in order to take into account the correlations. We also plan to use principal component analysis to decompose line correspondences into uncorrelated basic deformations, the importance of which will be completely quantified.

In parallel, we started to test the anatomical validity of the scheme by testing it on larger databases and by validating the results with specialists in skull and brain anatomy.

Finally, we are developing new applications for the skull atlas: a child skull growth study (Subsol et al., 1996b), skull evolution between prehistoric and contemporary man [START_REF] Subsol | Construction Automatique d'atlas Anatomiques Morphométriques à Partir d'images Médicales Tridimensionnelles[END_REF], facial reconstruction [START_REF] Quatrehomme | A fully three-dimensional method for facial reconstruction based on deformable models[END_REF] and sex assessment from the skull.

DESCRIPTION OF THE VIDEO

The accompanying video demonstrates the building of the atlas and its use for the study of a patient. The video runs for 1 min 50 s.

1. Automatically building of a 3-D skull atlas (50 s).

Non-rigid registration of crest lines extracted on two skulls (Section 4 of the paper).

Common crest line identification and labelling (Section 5 of the paper).

Average crest lines and surface of the skull (Section 6 of the paper).

2. Study of a patient (1 min) (Section 8 of the paper).

Presentation of a skull affected by a maxillary deformation.

Extraction of crest lines.

Registration with the atlas (Subsection 8.2 of the paper) for a qualitative study.

Morphometric comparison of crest lines of the mandible (in red, the atlas; in blue the patient) by modal analysis (Subsection 8.3 of the paper).

Figure 1 .

 1 Figure 1. The skulls A, B, C, D, E, F (left to right, top to bottom).These skulls were segmented from high-resolution CT scan images (acquired from two different devices, the voxel size is about 1 × 1 × 1.5 mm 3 ) by classic mathematical morphology and thresholding tools. We notice an important diversity in the skull size, orientation and shape.

.Figure 2 .

 2 Figure 2. A general scheme to build automatically computerized morphometrical anatomical atlases.

Figure 3 .

 3 Figure 3. Differential characteristics of a surface and the definition of a crest line.

Figure 4 .

 4 Figure 4. Crest lines of the skull F . Notice the inside lines emphasizing the sphenoid and temporal bones and also the foramen magnum.

Figure 5 .

 5 Figure 5. Comparison of crest lines (in grey) and ridge lines (in black) which were extracted semi-manually under the supervision of an anatomist. Their superimposition shows that crest lines have a strong anatomical significance.

Figure 6 .

 6 Figure 6. Crest lines of the skull F filtered by a hysteresis thresholding on the value of maximal curvature. Thus, we have discarded all the noisy lines which lay on the forehead and the top of the skull.

Figure 7 .

 7 Figure 7. Two sets of lines to be registered: the left-hand set is composed of 591 lines and 19 302 points; the right-hand one is composed of 583 lines and 19 368 points. We notice the variations in the shape, the number and the topology of lines.

Figure 8 .

 8 Figure 8. The feature registration algorithm.

Figure 9 .

 9 Figure 9. Computing registration parameters is not obvious due to the non-bijectivity of points matching.

Figure 10 .

 10 Figure 10. After discarding non-consistent matched points, we can compute the line registration parameters correctly.

Figure 11 .

 11 Figure 11. Example of application of transformations to a regular mesh displayed upper left: rigid (upper right) for position and orientation differences, affine (lower left) for scaling differences, and spline (lower right) for local and complex differences.

Figure 12 .

 12 Figure 12. The two simplified sets of crest lines in their original positions. The left-hand set will be deformed toward the right-hand one in order to find line and point correspondences.

Figure 13 .

 13 Figure 13. After the rigid transformations, the two sets of lines are aligned. We see only the registered lines and the matched points that are linked. Left, the sets are in their deformed positions. Right, they are in their original positions.

Figure 14 .

 14 Figure 14. After rigid and affine transformations, there is no longer any global size difference. Nevertheless, the orbits are still not aligned. We only see the registered lines and the matched points that are linked. Left, the sets are in their deformed positions. Right, they are in their original positions.

Figure 15 .

 15 Figure15. At the end of the iterations, after rigid, affine and spline transformations, the superimposition is very accurate and gives very accurate results for registered lines and matched points. We see only the registered lines and the matched points that are linked. Left, the sets are in their deformed positions. Right, they are in their original positions.

Figure 16 .

 16 Figure 16. Registration of C towards B. Left, we see the deformed set C with B. The matched points are linked. Notice how the two sets are reasonably superimposed. Right, C is in its original position. It allows us to estimate the extent of the deformation between the two sets.

Figure 17 .

 17 Figure 17. The registration graph: each node is a line of a set and the oriented link represents the relation 'is registered with'.

Figure 18 .

 18 Figure 18. The connected components of the registration graph define subsets of corresponding lines of different data sets. If these subgraphs contain at least one line of each data set, they define a subset of common lines.

Figure 19 .

 19 Figure19. The structure of the atlas is displayed for skulls B and C. Some subsets of common lines have been labelled automatically and highlighted: the mandible (bottom and top), the nose, the orbits, the cheekbones, the temples, the foramen magnum and the sphenoid and temporal bones.

Figure 20 .

 20 Figure 20. Automatic extraction of parts of the skull: on the left, the orbits, the nose and the mandible; on the right the left sphenoid bone and the foramen magnum.

•

  Register A and B, B and C, C and D, D and E, E and F , F and A, in both directions. Find the corresponding lines. • Build the registration graph. • Extract the connected components which contain at least one line from each data set. This gives the structure of the atlas.

Figure 21 .

 21 Figure 21. Left, the original data; right, data aligned in the frame defined by A. After application of similitudes, there are no longer global differences of position, orientation and scaling between the data; only meaningful morphometric deformations remain.

Figure 22 .

 22 Figure22. By using the registration algorithm, we can find the corresponding points between L i (in grey) and L i (C) (in black). But some points of L i remain without correspondent points because of the line topological constraints of the algorithm.

Figure 23 .

 23 Figure 23. After interpolation of the missing correspondent points, we obtain the deformation D i between L i and L i (C). We can check the result by applying D i to L i and comparing the obtained line (in black) with L i (C) which is displayed in the previous figure.

  y and d X ,Y i [k] z ) are the amplitudes (for the axes x, y and z) corresponding to the fundamental deformations φ p [i], called the modes. The set of all the amplitudes is called the spectrum of the deformation. The mode 0 represents the translation ∀k, φ 0

Figure 24 .

 24 Figure24. Application of mode 1 (top), 2 (middle) and 3 (bottom) with a constant amplitude to the mandibular line A (in grey, the original line; in black, the deformed line). The larger the mode number, the more complex the deformation.

Figure 25 .

 25 Figure 25. Left, the common lines corresponding to the left orbit in their original positions. Right, the average left orbit is in black among all the data aligned in the reference frame.

Figure 26 .

 26 Figure 26. Some common lines of the six skulls and in black the average common lines constituting the atlas.

Figure 27 .

 27 Figure 27. The set of all average common lines constituting the atlas. Notice their lateral symmetry.

Figure 28 .

 28 Figure 28. The six average models built from the six reference sets are very similar.

Figure 29 .

 29 Figure 29. The skull atlas obtained by the general scheme described in this paper.

Figure 30 .

 30 Figure 30. The skull S with a significant maxillary hypoplasia.

•

  Register ATLAS with respectively A, B, C, D, E and F . • Thanks to the registration result, align the data A, B, C, D, E and F in the reference frame defined ATLAS and remove similarity transformations. • For each common line of the L i of A: mode j, compute the mean value di [ j] and the standard deviation σ i [ j] of the amplitudes. These statistical parameters will define an amplitude distance on the modes, dist amp .

Figure 31 .

 31 Figure 31. The mandible of S automatically labelled and extracted by registration with the atlas (in black, the mandibular crest lines).

Figure 32 .

 32 Figure 32. Rigid and scale registration of the atlas (solid) with S (transparent) emphasizing the deformations of the mandible.

Figure 33 .

 33 Figure 33. The three 'abnormal' basic deformations of the S mandible. Left, the first x-mode which quantifies the breadth of the mandible; middle, the second x-mode which represents the twist of the mandible; right, the fourth z-mode which characterizes the lateral curvature.

Table 2 .

 2 Amplitude

		Mode 0 Mode 1 Mode 2 Mode 3 Mode 4
	x	0.124	2.442	2.476	0.745	0.353
	y	0.906	0.759	1.601	1.734	0.920
	z	1.806	0.896	1.062	1.017	2.267

distances for the first five modes of the mandible deformation between S and the atlas. According to the high values of these distances, the first and second x-mode and the fourth z-mode are considered 'abnormal'.

ACKNOWLEDGEMENTS

We thank General-Electric Medical Systems Europe, Dr Bruce Latimer from the Cleveland Museum of Natural History, Dr Court Cutting from New York University Medical Center, Dr David Dean from Case Western Reserve University of Cleveland and Dr André Guéziec for the CT scan data of the skulls. We would also like to thank Dr Janet Bertot, Professor Mike Brady, Dr David Dean, Dr Jérôme Declerck and Dr Jacques Feldmar for their substantial help. We are grateful to the reviewers for their useful comments and advice.