Alexander V Tuzikov
email: tuzikov@mpen.bas-net.

Olivier Colliot
email: olivier.colliot@enst.fr

Isabelle Bloch
email: isabelle.bloch¥@enst.fr

Brain symmetry plane computation in MR images using inertia axes and optimization

Keywords: Symmetry plane, symmetry measure, MR brain image, mid-sagittal plane, inertia axes, optimization

Detection of the best symmetry plane in 3D images can be treated as a registration problem between the original and the reflected images. The registration is performed in a 3D space of parameters defining orientation and shift of a reflection plane. We use the normalized ¦ £ metric as the similarity measure between original and reflected images and investigate an algorithm for computation of the best symmetry plane. The algorithm computes first an initial position of the plane by analyzing principal inertia axes. We demonstrate on several MR brain images that the initial position is in the neighborhood of the global maximum. Therefore the downhill simplex method is further used for the computation of the best symmetry plane. The proposed algorithm was tested on simulated and real MR brain images.

Introduction

The mid-sagittal plane is defined as the plane separating both brain hemispheres. The automatic detection of this plane in brain images is a useful task. It can be used to reorient images, for example for further alignment in the Talairach reference frame [START_REF] Liu | Robust midsagittal plane extraction from normal and pathological 3D neuroradiology images[END_REF]. It can also be a first step in intra-subject registration [START_REF] Ardekani | Automatic detection of the mid-sagittal plane in 3D brain im-ages[END_REF] or serve as a basis for the study of dissymetries in the brain [START_REF] Thirion | Statistical Analysis of Normal and Abnormal Dissymmetry in Volumetric Medical Images[END_REF].

Brain images possess a high degree of plane symmetry although they are not exactly symmetrical. Medical studies (see, for example, [START_REF] Amunts | Interhemispheric asymmetry of the human motor cortex related to handedness and gender[END_REF] for some corresponding references) demonstrate brain anatomical asymmetries in sizes of the frontal and occipital lobes, in surface of the planum temporale and planum parietale, in the anterior speech region. The plane with respect to which the brain is the most symmetrical is called the brain symmetry plane. This plane is often considered as the first approximation of the mid-sagittal plane [START_REF] Kruggel | Alignment of magneticresonance brain datasets with the steretactical coordinate system[END_REF]. The main goal of this paper is to develop an algorithm for the symmetry plane computation in MR brain images.

There is a vast literature in mathematics, image processing and computer vision domains dealing with different kinds of symmetry (central, reflection, rotation, skew) of shapes and images. Let us point out, for example, some publications [START_REF] Marola | On the detection of the axes of symmetry of symmetric and almost symmetric planar images[END_REF][START_REF] Masuda | Detection of partial symmetry using correlation with rotated-reflected images[END_REF] devoted to the reflection symmetry. Most of practical algorithms developed for symmetry measurements are applied in the 2D case and only some of them can be extended to 3D [START_REF] Zabrodsky | Symmetry as a continuous feature[END_REF]. Principal axes of inertia were used in [START_REF] O'mara | Measuring bilateral symmetry in digital images[END_REF] to define the best plane of reflection symmetry. An algorithm for finding symmetry planes of 3D objects using extended Gaussian image representation was developed in [START_REF] Sun | 3-D symmetry detection using the extended Gaussian image[END_REF]. Octree representation of 3D objects is used in [START_REF] Minovic | Symmetry identification of a 3d object represented by octree[END_REF] to compute their symmetry degree. Recently several new algorithms were proposed for symmetry plane computation in 3D images of the brain [START_REF] Ardekani | Automatic detection of the mid-sagittal plane in 3D brain im-ages[END_REF][START_REF] Liu | Robust midsagittal plane extraction from normal and pathological 3D neuroradiology images[END_REF][START_REF] Prima | Computation of the Mid-Sagittal Plane in 3D Images of the Brain[END_REF][START_REF] Thirion | Statistical Analysis of Normal and Abnormal Dissymmetry in Volumetric Medical Images[END_REF].

The problem of finding symmetry plane and reflection symmetry measure can be treated as a rigid registration problem between the original and the reflected images. Therefore, one can use different similarity measures known in the literature [START_REF] Maintz | A survey of medical image registration[END_REF] (for example, correlation coefficient, correlation ratio, mutual information). Since the problem belongs to the class of monomodal registration, then the cor-relation coefficient or ¦ £ metric seem to be the most suitable similarity measures.

A method for computation of the brain symmetry plane using optimization technique is presented in [START_REF] Ardekani | Automatic detection of the mid-sagittal plane in 3D brain im-ages[END_REF]. A number of plane positions equally distributed in the space of plane orientations are used there as candidates for finding the initial plane position for the optimization procedure. We propose (Section 3) an algorithm that uses a single initial plane for the further optimization. This plane is computed from the ellipsoid of inertia and we show in Section 4 for MR brain images that it belongs to the region of the global maximum of the similarity measure under consideration. Therefore, starting the optimization from this plane gives more chances to achieve the global maximum and besides reduces substantially the complexity of the algorithm.

Symmetry measure definition

Let be a unit vector in

¡ £¢

, ¤ £ is the unit sphere of all possible directions in

¡ ¥¢

, and ¦ ¨ § © is a plane in ¡ £¢ orthogonal to the vector and passing at the distance from the coordinate origin. We denote by § © the reflection of image with respect to the plane ¦ § © :

§ © ! " # %$ '& )( 0& )1 2 43 5 ¥ § © 6 7$ '& 8( 0& 91 @ ) BA
An image is called reflection symmetrical if there ex- ists a reflection plane ¦ ¨ § © such that C § © " D3 E . We say in this case that ¦ ¨ § © is a plane of reflection symmetry for .

Suppose that one is able to compute the symmetry measure F ' § © " of the image with respect to an arbitrary re- flection plane ¦ § © . Then the problem we are dealing with in this paper is to find the plane of maximal symmetry degree of image and the corresponding value of symmetry measure F " :

F " 43 G IH 6P § Q !R @S B© T TQ VU XW F X § © " YA (1)
When comparing images, it is usually assumed that a distance or a similarity measure ` ¢ & 9 £ is defined on the set of images. In this case the symmetry measure F § © ! " can be defined as the similarity between images and § © 6 , i.e. F § © ! " 43 ` " a& § © ! " 8 . Different similarity measures are known in the literature for comparing images. In this report we use ¦ £ metric (¦ £ a& cb 43 ed fd hg ib Xd fd ) and the following symmetry measure:

F § © ! " 43 qp rg d sd tg V § © " Td sd £ u d fd £d fd £ A

Symmetry plane computation

To compute the reflection symmetry plane one needs to check all possible orientations and translations of the reflec-tion plane. It is an optimization problem in 3 dimensional parametric space (here v and w are the angles of the plane normal vector with the $ axis and $ x( plane, respectively, and is the distance to the plane from the coordinate ori- gin). This computation is time consuming. Therefore it seems useful to check first the most appropriate positions of reflection planes. These positions can be used as initial points for the further optimization.

The initial positions can be interpreted, for example, in terms of ellipsoid of inertia [START_REF] Goldstein | Classical Mechanics[END_REF]. The axes of the ellipsoid of inertia for a 3D body are defined by the eigenvectors of the corresponding covariation matrix.

It is known that if a body possesses a plane reflection symmetry then the symmetry plane passes through its center of mass and is orthogonal to some ellipsoid axis. Therefore for computation of a reflection symmetry measure in this case it is sufficient to check 3 planes only.

Suppose that the center of mass of the ellipsoid of inertia of image is y £ 3 %$ & )( & )1 and the eigenvectors are

¢ & £ & ¢ . Then one can compute G IH 6P ¤ F § @ 9© § 2 9 2 V " B& F § S © § S @ " Y& F § ! B© § Y @ V ¥ &
where #y denotes the inner product of vectors and y .

The best position of the reflection symmetry plane found from the ellipsoid of inertia can be further improved using optimization in the 3D space of reflection plane parameters. In our implementation we used the downhill simplex method [START_REF] Press | Numerical Recipes in C. 2nd Edition[END_REF]. Since the downhill simplex method is a local optimization method there is no guaranty in general of finding the global maximum. However, in Section 4 we demonstrate experimentally that one can obtain good solutions for MR brain images using this approach.

Algorithm 1

1. Computation of the initial symmetry plane from the ellipsoid of inertia.

2. Improving the orientation and translation of the symmetry plane and the symmetry measure computation using the optimization technique.

Experiment discussion

We investigated the graphs of F ' § © for different val- ues of and for 5 real MR normal brain images. These graphs have a similar form (see Fig. 1). Fig. 1(a) shows the profile of F § © for the brain image presented in The normal brain database contains simulated brain MR images for three modalities (T1-, T2-, and PD(proton-density)weighted) and a variety of slice thicknesses, noise levels, and levels of intensity non-uniformity [START_REF] Collins | Design and construction of a realistic digital brain phantom[END_REF][START_REF] Kwan | MRI simulation-based evaluation of image-processing and classification methods[END_REF]. We tested the algorithm for 54 images for each modality (with 3 slice thicknesses 1,3,5 mm, 6 levels of noise 0%,1%,3%,5%,7%,9% and 3 levels of intensity nonuniformity 0%,20%,40%). An additional advantage of using this database for symmetry analysis is that the images are aligned in the stereotactic space and the mid-sagittal plane coincides with ( 1 plane. For all these images the orientation of the best symmetry plane found by Algorithm 1 coincides with the mid-sagittal plane. The maximal absolute difference ¤ 9d w d for w is A ¥ u degree, the mean absolute value ¤ §¦ )d w d 3 A sp ©¨degree and the standard deviation c )d w d 3 A sp degree. The values for v are

¤ )d v d 3 A & ¤ §¦ 9d v d 3 A fp ¨and c )d v d 3 A
degrees, respectively. This confirms that the proposed algorithm is stable in the presence of noise and image non-uniformity. Algorithm 1 was also tested for 24 real MR normal brain images (in T1 modality) and for all of them we obtained good final results (see an example in Fig. 2). For 19 images we applied the process directly on original images (without any pre-processing or segmentation) and obtained good final results (see third and fourth rows of Fig. 2). For 5 others we compared the position of the symmetry plane found for the original image and for the corresponding segmented brain. The segmentation of the brain was obtained by a 3D mathematical morphology algorithm [START_REF] Mangin | From 3D magnetic resonance images to structural representations of the cortex topography using topology preserving deformations[END_REF]. The final results were similar. The mean value and the standard deviation of absolute differences for w and v were de- grees respectively. Here indices 1 and 2 correspond to original and segmented brain images. Therefore, it seems that the algorithm can be applied successfully for original MR images. For 4 images we compared the results obtained for a binary mask of the segmented brain and for the gray-level image itself. The final results were also simi-

¤ §¦ )d w ¢ g w £ d 3 A & c 9d w ¢ g w £ d 3 A and ¤ ¦ )d v ¢ g v £ d 3 A & c )d v ¢ g v £ d 3 A ! "
lar (¤ §¦ )d w £ g w ¢ d t3 A & c )d w £ g w ¢ d t3 A ¥ #¨and ¤ ¦ )d v £ g v ¢ d 3 A ẍp & £ c )d v £ g v ¢ d 3 A degrees as-
suming that index 3 corresponds to a brain mask image). It means that the brain shape defines the position of the symmetry plane with a high precision.

Conclusion

We proposed an algorithm for computation of the best symmetry plane and the plane symmetry measure for 3D images. The initial position of the plane is computed from the corresponding ellipsoid of inertia. This position is further improved using the downhill simplex method. This algorithm can be applied for MR brain images. We demonstrated on several simulated and real MR brain images that the initial position is in the neighborhood of the symmetry measure global maximum. Therefore a local optimization method allows to find efficiently the global maximum.

It was shown on synthetic MR brain images that the algorithm is stable to noise and intensity non-uniformity. Further work includes comparing the method to other approaches and studying its behavior in pathological cases.

A. Tuzikov was supported by a sabbatical grant at the Department of Signal and Image Processing (ENST) and is very grateful to the department for its hospitality. We would like to thank CHU La Pitié-Salpétrière (Paris), CHU La Timone (Marseille), INSERM 280 and Hôpital Neurologique (Lyon) for providing the images.

Fig. 2

 2 Fig. 2 for reflection planes passing through the image center of mass. The function values are computed for angles w & v g r u & ) u with an incremental step of A @ radian. The function graph has a regular form with a global

Figure 1 .

 1 Figure 1. The graphs of function F § © com- puted for different values of and .

Figure 2 .

 2 Figure 2. Coronal and axial slices of the segmented and non-segmented brain images: (first row) initial symmetry plane computed from the ellipsoid of inertia; (second row) the best symmetry plane found by Algorithm 1; (third row) initial plane; (fourth row) the best symmetry plane.


Olivier Colliot
email: olivier.colliot@enst.fr

Isabelle Bloch
email: isabelle.bloch@enst.fr

Alexander V Tuzikov
email: tuzikov@mpen.bas-net.

Characterization of approximate plane symmetries for 3D fuzzy objects

Keywords: fuzzy object, plane symmetry, symmetry measure, measure of comparison

We are interested in finding and characterizing the symmetry planes of fuzzy objects in 3D space. We introduce first a fuzzy symmetry measure which defines an object symmetry degree with respect to a given plane. It is computed by measuring the similarity between the original object and its reflection. The choice of an appropriate measure of comparison is based on the desired properties. In a second part, a method for finding the best symmetry planes of fuzzy objects is proposed. We then apply these results to the representation of directional relationships.

Introduction

Symmetry is an important property of objects and also a useful feature for their recognition. However, exact symmetry does not exist in real objects and one has to deal with approximate symmetries. There is a vast literature in mathematics, image processing and computer vision domains dealing with different kinds of symmetry (central, reflection, rotation, skew) of shapes and images. Many works on inexact symmetries quantify the degree of symmetry, using a symmetry measure or distance (see, for example, [START_REF] Marola | On the detection of the axes of symmetry of symmetric and almost symmetric planar images[END_REF][START_REF] O'mara | Measuring bilateral symmetry in digital images[END_REF][START_REF] Tuzikov | Symmetry measure computation for convex polyhedra[END_REF][START_REF] Zabrodsky | Symmetry as a continuous feature[END_REF]). Most results on symmetry degree evaluation are obtained for precisely defined objects. Zabrodsky et al. [START_REF] Zabrodsky | Symmetry as a continuous feature[END_REF] consider shapes with uncertainty, i.e. shapes for which the location of each point is given as a probability distribution.

In this paper we are interested in characterizing approximate symmetries of fuzzy 3D objects. To our knowledge, this problem has not been addressed before. Following a classical approach used for crisp shapes and images, we first introduce a fuzzy symmetry measure which characterizes an object symmetry degree with respect to a given plane (Section 2). For this we use a measure of comparison 1 between the object and its reflection. Various measures of comparison have been proposed in the literature for fuzzy sets. The choice of a measure that is appropriate to our problem is based on its properties (Section 3). In Sec-tion 4, we present an algorithm for finding the best symmetry plane of an object. In Section 5, the results of this paper are applied to the representation of directional spatial relationships.

2 Symmetry measure

Reflection of a fuzzy object

Let Π be a plane in the 3D space R 3 . Given a point x, we denote by e Π (x) its image under the reflection with respect to Π. e Π is a bijective transformation in R 3 . Therefore, one can define the reflection of a fuzzy set as follows.

Definition 2.1 The reflection of a fuzzy set A is a fuzzy set e Π (A) defined as:

µ e Π (A) (e Π (x)) = µ A (x) for every x ∈ S.
We denote by e u,d the reflection with respect to a plane Π u,d which is orthogonal to u and passing at the signed distance d from the origin. In spherical coordinates a unit vector u is defined by two angles β ∈]π, π] and α ∈ [-π/2, π/2] (see Fig. 1). As vectors u and -u define the same plane, we use β ∈ [0, π[, α ∈]π/2, π/2] and d ∈ R.

Symmetry measure

We want to define a symmetry degree of a fuzzy object with respect to a given plane Π. One option is to compare A and e Π (A). A symmetry measure σ A can be defined as a measure of comparison between the original object and its reflection:

σ A (Π) = S(A, e Π (A)),
where S is a measure of comparison between fuzzy objects. As before, we use notations

σ A (u, d) = σ A (α, β, d) = σ A (Π u,d ).
Various measures of comparison have been proposed in the literature. They possess different properties and the choice of a measure depends on the application and of the concept one wants to describe. Below we discuss the measures that can be used to define a symmetry measure.

Deriving symmetry measures from measures of comparison

First we present properties that can be used to distinguish measures of comparison and discuss which of them should be satisfied by a symmetry measure. Then we summarize which of these properties hold for different measures of comparison proposed in the literature and select some of them to define symmetry measures.

Desired properties of symmetry measures

Measures of similitude Bouchon-Meunier et al. [START_REF] Bouchon-Meunier | Towards general measures of comparison of objects[END_REF] have proposed a classification of measures of comparison between fuzzy sets, in particular M-measures of comparison which are derived from a fuzzy measure M .

Definition 3.1 [4] An M-measure of com- parison is a mapping S : F × F → [0, 1] such that S(A, B) = F S M (A ∩ B), M (B - A), M (A -B) for a given mapping F S : R + × R + × R + → [0, 1].
A particular class of measures of comparison is composed of measures of similitude.

Definition 3.2 [START_REF] Bouchon-Meunier | Towards general measures of comparison of objects[END_REF] An M-measure of similitude is an M-measure of comparison S such that F S (u, v, w) is non-decreasing in u, nonincreasing in v and w.

M-measures of similitude are well suited for describing symmetries: symmetry is stronger if the measure of intersection between the original object and its reflection increases, and it is weaker if the difference between them increases. Measures of similitude include measures of satisfiability and measures of resemblance.

Definition 3.3 [4]
1. An M-mesure of satisfiability is an Mmeasure of similitude such that

• F S (u, v, w) is independent of w; • F S (0, v, w) = 0, for all v, w (exclusivity); • F S (u, 0, w) = 1,
for all u = 0 (inclusion).

2. An M-measure of resemblance is an Mmeasure of similitude such that

• S is reflexive, i.e. S(A, A) = 1. • S is symmetrical, i.e. S(A, B) = S(B, A).
In our case M (Ae Π (A)) = M (e Π (A) -A), therefore measures of satisfiability are also measures of resemblance [START_REF] Bouchon-Meunier | Towards general measures of comparison of objects[END_REF]. Moreover, the exclusivity property entails that the symmetry measure is equal to zero when the plane passes outside the support of the object. The inclusion property, as well as reflexivity, entails that the symmetry degree is equal to 1 when the object coincides with its reflection, i.e. when there is an exact symmetry. The symmetry property implies that the symmetry measure for an object A with respect to a given plane Π is equal to the measure computed for its reflection e Π (A) 2 . Therefore, Mmeasures of satisfiability seem to be suitable for the definition of measures of symmetry.

Additionnal properties

Pappis [START_REF] Pappis | A comparative assessment of measures of similarity of fuzzy values[END_REF] proposes the following additional properties which are in fact the reverse implications of reflexivity and exclusivity:

S(A, B) = 1 ⇐⇒ A = B,
2 Proofs are omitted here for the sake of conciseness

S(A, B) = 0 ⇐⇒ supp(A) ∩ supp(B) = ∅.
The first property, also called separability for distances, expresses that the symmetry measure is equal to 1 if and only if there is an exact symmetry. The second one expresses that the symmetry measure equals zero if and only if the plane passes outside the support of the object.

Geometrical properties

Intuitively speaking a symmetry measure should be invariant with respect to translation, rotation and scaling. If S is invariant w.r.t. translation (resp. rotation) then so is σ. This is also true for scaling but as the scaling of a fuzzy set in the discrete case is not clearly defined, we will not consider it later on.

Definition 3.4

The symmetry measure of A with respect to Π is defined as:

σ A (Π) = S(A, e Π (A)),
where S is a measure of comparison with the following properties:

(P1) Symmetry: S(A, B) = S(B, A); (P2) Reflexivity: S(A, B) = 1 ⇐⇒ A = B;
(P3) S(A, B) = 0 if and only if the supports of A and B are disjoint;

(P4) S is invariant w.r.t. translation;

(P5) S is invariant w.r.t. rotation.

Other properties of measures of comparison considered for instance in [START_REF] Bouchon-Meunier | Towards general measures of comparison of objects[END_REF][START_REF] Liu | Entropy, distance measure and similarity measure of fuzzy sets and their relations[END_REF][START_REF] Pappis | A comparative assessment of measures of similarity of fuzzy values[END_REF] are either equivalent to these ones or not interesting for deriving symmetry measures.

Existing measures

We use here a classification of measures of comparison that is very similar to those used in [START_REF] Zwick | Measures of similarity among fuzzy concepts: A comparative analysis[END_REF] and [START_REF] Bloch | On fuzzy distances and their use in image processing under imprecision[END_REF].

Set-theoretic approach

Most of the measures discussed in this section have been derived from a general measure proposed by Tversky [START_REF] Tversky | Features of similarity[END_REF] and are based on combinations of µ A and µ B using t-norms and t-conorms. They satisfy (P1) as t-norms and t-conorms are commutative. The following measure has been used by several authors [START_REF] Bouchon-Meunier | Towards general measures of comparison of objects[END_REF][START_REF] Dubois | Fuzzy Sets and Systems, Theory and Applications[END_REF][START_REF] Pappis | A comparative assessment of measures of similarity of fuzzy values[END_REF][START_REF] Zwick | Measures of similarity among fuzzy concepts: A comparative analysis[END_REF] 3 :

S 1 (A, B) = x∈S (µ A (x), µ B (x)) x∈S ⊥(µ A (x), µ B (x)
) Property (P2) holds if and only if = min and ⊥ = max. Property (P3) is fulfilled for t-norms "minimum" and "product" but is not for "drastic" and "Lukasiewicz" ones [START_REF] Dubois | Fuzzy Sets and Systems, Theory and Applications[END_REF]. Properties (P4) and (P5) are fulfilled.

Wang [START_REF] Wang | New similarity measures on fuzzy sets and on elements[END_REF] proposed the following measure of comparison:

S 2 (A, B) = 1 |S| × x∈S (µ A (x), µ B (x)) ⊥(µ A (x), µ B (x))
with 0 0 = 1. S 2 satisfies (P2) if and only if = min and ⊥ = max but does not satisfy (P3).

However, it is easy to check that a modified version of S 2 defined as follows:

S 3 (A, B) = 1 |supp(A) ∪ supp(B)| × x∈supp(A)∪supp(B) (µ A (x), µ B (x)) ⊥(µ A (x), µ B (x))
satisfies property (P3) for t-norms "minimum" and "product". Properties (P4) and (P5) are also fulfilled.

Hyung et al. [START_REF] Hyung | Similarity measure between fuzzy sets and between elements[END_REF] proposed to use a measure of comparison defined as

S 4 (A, B) = max x∈S (µ A (x), µ B (x)).
S 4 satisfies property (P3) for "minimum" and "product" t-norms but does not satisfy (P2).

L p distance approach

In this section we use the L p distance between fuzzy sets A and B:

A -B p = x∈S |µ A (x) -µ B (x)| p 1 p
3 Here we deal with the finite discrete case. In the continuous case, the sum is replaced by an integral if it converges.

A -B ∞ = max x∈S (|µ A (x) -µ B (x)|).
Measures of comparison based on the L p distance have the following general form:

S(A, B) = 1 - A -B p K ,
where K is a normalization coefficient. It is easy to see that properties (P1), (P2), (P4) and (P5) are fulfilled for measures of this type.

For example, Wang [START_REF] Wang | New similarity measures on fuzzy sets and on elements[END_REF] and Bouchon-Meunier et al. [START_REF] Bouchon-Meunier | Towards general measures of comparison of objects[END_REF] proposed the following measure

S 5 (A, B) = 1 - A -B 1 |S| .
This measure does not satisfy property (P3).

The following measure of comparison proposed by Pappis [START_REF] Pappis | A comparative assessment of measures of similarity of fuzzy values[END_REF]

S(A, B) = 1 -x∈S |µ A (x) -µ B (x)| x∈S [µ A (x) + µ B (x)] can be generalized as S 6 (A, B) = 1 - A -B p x∈S µ A (x) p + µ B (x) p 1 p
Measure S 6 satisfies property (P3).

Pappis [START_REF] Pappis | A comparative assessment of measures of similarity of fuzzy values[END_REF] also proposed to use the L ∞ distance

S 7 (A, B) = 1 -A -B ∞ S 7
does not satisfy property (P3). However, when A and B are normalized fuzzy sets and their supports are disjoint S 7 (A, B) = 0. But the converse implication is still false.

It is easy to verify that the measure of comparison (β > 0)

S 8 (A, B) = e -β A-B p
proposed in [START_REF] Bouchon-Meunier | Towards general measures of comparison of objects[END_REF] does not satisfy property (P3) either.

Correlation coefficient

Gerstenkorn [START_REF] Gerstenkorn | Correlation of intuitionistic fuzzy sets[END_REF] introduced a correlation coefficient between fuzzy sets:

S 9 (A, B) = C(A, B) T (A)T (B) ,
where

C(A, B) = x∈S [µ A (x)µ B (x) + (1 -µ A (x))(1 -µ B (x))]
and

T (A) = x∈S µ A (x) 2 + (1 -µ A (x)) 2 .
Measure S 9 (A, B) satisfies properties (P1) and (P2) but does not satisfy (P3).

Chosen symmetry measures

The measures of comparison S 1 , S 3 (for tnorm "minimum") and S 6 satisfy properties (P1)-(P5). Therefore we define three symmetry measures:

σ 1,A (Π) = x∈S min(µ A (x), µ e Π (A) (x)) x∈S max(µ A (x), µ e Π (A) (x)) σ 2,A (Π) = 1 |supp(A) ∪ supp(e Π (A))| × x∈supp(A)∪supp(e Π (A)) min(µ A (x), µ e Π (A) (x)) max(µ A (x), µ e Π (A) (x)) σ 3,A (Π) = 1 - A -e Π (A) p x∈S µ A (x) p + µ e Π (A) (x) p 1 p .

Finding symmetry planes of a fuzzy object

It is of interest to locate automatically symmetry planes of a given object. In this section we show how a symmetry measure can be used for this purpose. Figure 3 shows another 2D fuzzy object. In the direction defined by β = 0, the maximum of σ 1,A is obtained for d = 10 which is the position of the symmetry plane of the alpha-cut of level 0.5. This result fits well with intuition. It can be useful to perform some preliminary smoothing on σ A to exclude some nonsignificant local maxima. They can appear due to the object itself or due to artefacts generated by discretization. These artefacts can appear, for example, when steps on α or β become small. We do not study in this paper the influence of discretization on the symmetry measure.

Efficient search of symmetry planes

Whereas the computation of σ A for a sufficiently small step is feasible in the 2D case, it is a far too expensive operation for 3D objects.

In many cases, one only wants to locate the best symmetry plane of an object corresponding to the largest symmetry measure value.

We propose a method that expresses the problem of finding the best symmetry plane as an optimization problem in the parametric space ]π/2, π/2] × [0, π[×R.

The optimization procedure needs a starting point. We suggest to use the ellipsoid of inertia to define candidates for this starting point. The ellipsoid of inertia has already been used in [START_REF] O'mara | Measuring bilateral symmetry in digital images[END_REF] to define the symmetry plane of an object. Here it is only taken as an initialization. The directions of axes are defined as the eigenvectors of the covariation matrix: Here m pqr defines a central moment of order

p + q + r m pqr (A) = S µ A (x, y, z)(x -x c ) p (y -y c ) q (z -z c ) r ,
where c = (x c , y c , z c ) is the object center of mass. If a 3D object possesses an exact plane of symmetry it passes through its center of mass and is orthogonal to one of the ellipsoid axes. Let us denote by u 1 , u 2 and u 3 the eigenvectors of the covariation matrix. We consider then three planes orthogonal to these vectors and passing through the center of mass:

Π 1 = Π u 1 ,u 1 •c , Π 2 = Π u 2 ,u 2 •c and Π 3 = Π u 3 ,u 3 •c . Our initial symmetry plane Π i maximizes the symmetry measure, i.e. σ A (Π i ) = max{σ A (Π 1 ), σ A (Π 2 ), σ A (Π 3 )}.
This is only possible when the eigenvectors are different. Otherwise, one gets an ellipsoid of revolution.

Then, the best symmetry plane is found using an optimization technique in the parametric space ]π/2, π/2] × [0, π[×R. We use the Nelder-Mead downhill simplex method [START_REF] Press | Numerical Recipes in C. 2nd Edition[END_REF] which was also used in [START_REF] Ardekani | Automatic detection of the midsagittal plane in 3D brain images[END_REF] but with a different initialization. This method is often used when one does not know the function derivatives. It is accurate and robust under a good starting point. However, it is a local optimization method and, in general, one has no guaranty to find the global maximum.

Figure 4 shows the symmetry plane found by this method. The image is a fuzzy segmentation of the lateral ventricles in an MR image. The symmetry measure σ 1,A of the object with respect to this plane is 0.73. We also applied this method on gray-level images (but with a different symmetry measure) to compute the brain symmetry plane in MR images [START_REF] Tuzikov | Brain symmetry plane computation in MR images using inertia axes and optimization[END_REF]. It has shown very good results.

Application to the representation of directional spatial relationships

Spatial relationships can be very useful for scene recognition and interpretation. Here we deal with directional relationships such as to the left of, above etc. There are basically two ways to define directional relationships: with respect either to an extrinsic frame or to an instrinsic frame which is defined by the reference object. The latter case often occurs when the reference object presents an approximate reflectional symmetry. The frame can then be partially defined by the symmetry plane of the reference object. For example, in the common cases of a human body or a car, relation-ships left and right refer to directions which are orthogonal to the symmetry plane of this object. A typical example can be found in neuroanatomy: the brain is an approximately symmetrical scene in which directions are defined with respect to an intrinsic frame. This frame is partially defined by the mid-saggital plane which approximately corresponds to the brain symmetry plane.

The results of this paper can be easily integrated in the representation proposed in [START_REF] Bloch | Fuzzy relative position between objects in image processing: a morphological approach[END_REF].

Let us recall the principle of this method (the relationship is defined between a reference object R and a target object A):

1. Definition of a fuzzy landscape around the reference object. This landscape is a fuzzy subset of S such that the membership value of each point corresponds to the degree of satisfaction of the spatial relation.

Comparison of the object

A to the fuzzy landscape. This is done using a fuzzy pattern matching approach.

It is shown that the fuzzy landscape associated to the direction defined by the unit vector u α,β (α ∈]π/2, π/2] and β ∈ [0, 2π[) corresponds to the fuzzy dilation of R by the following structuring element B α,β : The algorithm presented in Section 4 can be used directly to define directional relationships. The plane orientation is defined by a normal vector u which can be used directly for the computation of the fuzzy landscape.

µ B α,β (x) = f (arccos v x • u α,β v
When not only the reference object but the whole scene is approximately reflection symmetrical, the scene is often better described by relationships outside and inside than left and right. Indeed, left and right would have to

Conclusion

In this paper, we studied approximate plane symmetries of fuzzy objects. Three symmetry measures were derived from measures of comparison. They show good properties. Using these measures, we also proposed an algorithm for finding the best symmetry plane of a 3D fuzzy object. This method was applied to the definition of spatial relationships in symmetrical scenes.

We are now working on the use of symmetry measures to represent imprecisely located symmetry planes as fuzzy subsets of the parametric space representing orientations and directions. It would also be interesting to use symmetry measures as features for scene recognition. They could be used, for example, to define symmetry attributes in a fuzzy attributed graph [START_REF] Perchant | Semantic spatial fuzzy attribute design for graph modeling[END_REF].

Figure 1 :

 1 Figure 1: Angles α and β define a unit vector u We also use notation e α,β,d instead of e u,d in 3D, and e β,d = e 0,β,d in 2D.

4. 1 Figure 2 :

 12 Figure 2: (a) A 2D fuzzy set A. (b) σ 1,A . (c) σ 1,A for d = 0. (d) σ 1,A for β = 0.

Figure 3 :

 3 Figure 3: (a) A 2D fuzzy set A. (b) σ 1,A . (c) σ 1,A for d = 0. (d) σ 1,A for β = 0.

m

  200 m 110 m 101 m 110 m 020 m 011 m 101 m 011 m 002  

Figure 4 :

 4 Figure 4: Left, one slice of a 3D object with its symmetry plane. Right, 3D renderings of the alpha-cut of level 0.5 of this object (with and without the symmetry plane).

  x ) for all x ∈ S, and µ B α,β (o) = 1, where f is a decreasing function on [0, π] e.g. f (θ) = max[0, 1-(2/π)θ], o is the center of the structuring element and v x is a vector from o to x.

Figure 5 Figure 5 :

 55 Figure 5 shows an example of such fuzzy landscapes.

We prefer to use the expression "measure of comparison" as in[START_REF] Bouchon-Meunier | Towards general measures of comparison of objects[END_REF] instead of "similarity measure", since different authors assume different properties for the notion of similarity measure.

Acknowledgments

A. Tuzikov was supported by a sabbatical grant at the Department of Signal and Image Processing (ENST) and is very grateful to the department for its hospitality.


O Camara
email: camara@enst.fr

O Colliot

G Delso

I Bloch

3D nonlinear PET-CT image registration algorithm with constrained Free-Form Deformations

This paper presents a 3D nonlinear PET-CT image registration method guided by a B-Spline Free-Form Deformations (FFD) model, dedicated to thoracic and abdominal regions. It is divided into two stages: one FFD-based registration of structures that can be identified in both images; and a whole-image intensity registration step constrained by the FFD computed during the first step. Different similarity criteria have been adopted for both stages: Root Mean Square (RMS) to register recognized structures and Normalized Mutual Information (NMI) for optimizing the whole-image intensity stage. Structure segmentation is performed according to a hierarchical procedure, where the extraction of a given structure is driven by information derived from a simpler one. This information is composed of spatial constraints and expressed by the means of regions of interest, in which a 3D simplex mesh deformable model based method is applied. The results have been very positively evaluated by three medical experts.

Introduction

In oncology, the joint use of anatomical and functional images is increasing, which can be explained both by the development of acquisition devices and methods, and by the complementarity between such modalities. On the functional side, the development of Positron Emission Tomography (PET) acquisitions with fluorine 18 fluorodeoxyglucose (18-FDG) tracer gives access to a rich information for diagnosis and therapeutic follow-up of both primary and metastatic cancers. But the anatomical information is reduced, making difficult to localize the tumors with a high precision with respect to the organs. This anatomical information is provided by Computed Tomography (CT) or Magnetic Resonance Imaging (MRI) and allows the physician to get an accurate localization of the lesions, as well as size and shape measures. Unfortunately it does not provide suf-ficient knowledge about the lesion malignancy. Therefore, combining information from these two modalities would have a significant impact on improving medical decisions for diagnosis, therapy and treatment planning [START_REF] Wagner | Fused Image Tomography: An Integrating Force[END_REF]. Such a combination calls for a registration step in order to achieve a good correspondence between the images and the structures they contain.

The aim of this paper is to propose a contribution to the registration of CT and PET images for thoracic and abdominal applications. For such problems, rigid registration is not sufficient since it does not account for the strong deformations existing between both images, and it is necessary to develop nonrigid registration methods. These deformations, induced by the specificities of the acquisitions, are mainly due to the different acquisition protocols involved and the elastic nature of the imaged regions.

We first present a brief review of nonlinear registration algorithms. We then discuss the segmentation procedure and describe the registration method for segmented structures and for grey-level intensities. Finally, some results and a discussion of the technique are presented.

Registration framework

The goal of image registration is to determine the transformation that maps the information contained in one image into its anatomical correspondence in the other. Image registration methods can be divided according to several criteria: area of interest, nature of the transformation, similarity measure and optimization method. An exhaustive review of registration methods can be found in [START_REF] Maintz | A Survey of Medical Image Registration[END_REF].

Several nonlinear transformations can be found in the image processing literature. B-Spline Free Form Deformations (FFD) is a parametric model which provides a flexible nonlinear transformation due to the fact that no assumptions on the images or structures to register are made. This model have been successfully used in different medical imaging applications, such as pre-and post contrast MR mammogram registration [START_REF] Rueckert | Nonrigid Registration Using Free-Form Deformations: Applications to Breast MR Images[END_REF], brain registration [START_REF] Hartkens | Using points and surfaces to improve voxel-based non-rigid registration[END_REF] or cardiac segmentation [START_REF] Lotjonen | Segmentation of MR images using deformable models: Application to cardiac images[END_REF]. In general, they perform a linear transformation before the nonlinear phase to get closer enough to the final solution. Anyway, most of these applications work with monomodality data, where relationship between corresponding intensities in both images to register is simpler than in multimodality applications. This fact is aggravated when functional images are considered due to the large amount of noise and artifacts. Nevertheless, Mattes et al. [START_REF] Mattes | Pet-ct image registration in the chest using free-form deformations[END_REF] have used the FFD model to register chest transmission PET-CT images, applying a hierarchical and multiresolution scheme to avoid local minima and to alleviate the need for accurate initialization. One problem of this method is that it considers a functional relation between transmission and emission PET image acquisitions. Besides, authors point out that results are not satisfactory in the regions with larger deformations, such as at the diaphragm and the abdomen.

One way to avoid this problem is to constrain these deformations to avoid the convergence towards local minima and to reduce the computational cost of the algorithm. These constraints led us to propose a methodology divided into an initialization phase registering segmented structures from both images, and a second registration based on the whole-image intensities, refining the previous phase. Transformations in both steps are modeled by means of Free Form Deformations (FFD), governed by a grid of several control points for each dimension. The segmentation step is achieved by a 3D simplex mesh deformable model over regions of interest in a hierarchical identification procedure.

Structure segmentation

Overview of the procedure

The initial registration step requires the segmentation of some thoracic and abdominal structures in both images. This first registration will then be refined using grey-level information. Possible segmentation errors will not be propagated to the final result, as the fine registration step will be able to correct them.

Different levels of difficulty in the segmentation of structures suggest the use of a hierarchical procedure: the extraction of a given structure will be driven by information derived from a simpler one. This information can be composed of spatial constraints inferred from the previously segmented structures and be expressed by the means of Regions Of Interest (ROI) in which the search for new structures will take place. The structures to be segmented for further registration are lungs, kidneys and liver (in this order). Although they will not be used in the registration process, the skin and skeleton are also extracted as first steps of the Lungs, kidneys and liver are treated in two different stages: a first stage is composed of automatic thresholding and mathematical morphology operations in the ROI defined by previously segmented objects. The second stage consists in refining the result using a 3D deformable model. Skin and skeleton are segmented using only the first stage.

First stage: rough segmentation

To constrain the segmentation, a ROI is defined using spatial relationships with respect to other structures. These relationships include directions (for example the liver is below the lungs) and exclusion constraints (previously obtained structures are subtracted from the ROI so that no pair of objects overlap). The ROI for each structure are defined as follows (see Figure 1 for examples of ROI): skin and skeleton: they constitute the first steps of the procedure in CT and therefore no ROI is used. To locate the pelvis, we compute the area of the bounding box of the skeleton on each axial slice. The first significant decrease of this size (starting from the bottom) gives an indication of the top of the pelvis; lungs: in CT, the ROI is derived from the skin. In PET, we segment the lungs on the transmission image and dilate them to produce a ROI in the emission image; kidneys: in CT images, the region is bounded using the chest dimensions we have learned from the segmentation of the skeleton. An upper bound in the z axis is derived from the lungs: the ROI is defined below a line linking the lower-left limit of the right lung and the lower-right limit of the left lung (a line is drawn on each coronal slice). A lower bound can be computed using the top of the pelvis extracted from the skele- ton. In PET images, as the skeleton is not available, we only use the upper bound derived from the lungs; liver: the ROI is the same as above except that it excludes the kidneys. The ROI has been particularly useful to separate it from the heart and the kidneys.

Within this defined region, we then perform the following pipeline: k-means automatic thresholding, binary erosion, selection of connected components, binary dilation and 3D hole filling. In the case of the lungs and the liver, we select the biggest connected components while for the kidneys the two most symmetrical components with respect to the body symmetry plane are extracted using an algorithm proposed in [START_REF] Colliot | Characterization of approximate plane symmetries for 3D fuzzy objects[END_REF].

Second stage: refinement using a 3D deformable model

The first stage cannot be considered as a final segmentation. The main problem is the lack of regularization term. This has proven to be a problem in the case of subtle structures which could be detected in one modality but not in the other, thus introducing a difference that the registration process would wrongly interpret as a deformation and would try to compensate. A 3D deformable model has been implemented to overcome this problem. Deformable models are curves or surfaces defined within an image that evolve under constraints computed from the image data and regularity constraints.

We chose to implement a discrete model based on simplex meshes (introduced by Delingette [START_REF] Delingette | General object reconstruction based on simplex meshes[END_REF]). As a good initialization is very useful to achieve a fast convergence of the model, an initial surface is derived from the first segmentation stage. The segmentation obtained at the first stage can be eroded to ensure that the initialization is inside the object. Then, it is transformed to a triangulation using an isosurface algorithm [START_REF] Lorensen | Marching cube, a high resolution 3D surface reconstruction algorithm[END_REF]. It is decimated and converted to a simplex mesh using the dual operation.

The evolution of the deformable surface is described by the following dynamic force equation:

¡ ¡ £¢ ¥¤ §¦ © ¦ ! #" $ % (1)
where ¦ !& ' is the internal force that specifies the regularity of the surface and ¦ (" $ is the external force that drives the surface towards image edges. The chosen internal force is:

¦ ©& ' )¤ 10 32 54 76 98 2 @4 A B2 54 %

where 0 and 8

respectively control the surface tension (prevent it from stretching) and rigidity (prevent it from bending) and 2 4 is the Laplacian operator.

In our case, the external force is not only derived from image edges but also constrains the deformable model to stay in the ROI. It can be written as a linear combination:

¦ #" $ ¤ 1C ED GF ¦ ©H PI Q (3)
where D is a Gradient Vector Flow (GVF) (introduced by Xu et al. [START_REF] Xu | Gradient vector flow: A new external force for snakes[END_REF]), ¦ !H PI Q is a force attached to the ROI and C and F are weighting parameters. A GVF field is computed by diffusion of the gradient vector of a given edge map.

The edge map is derived from the gradient after performing an anisotropic diffusion. Anisotropic diffusion [START_REF] Gerig | Nonlinear anisotropic filtering of MRI data[END_REF] is an efficient way to remove noise in homogeneous regions while preserving and even enhancing edges, which proves to be very helpful for noisy structures such as the liver in PET images.

The second term of the external force is used to prevent the deformable model from going outside the ROI. ¦ H PI Q is a distance potential force and it can be written as follows:

¦ H PI Q ¡ ¤ 6 2 ¢ ¡ £ 2 ¢ ¡ £ (4)
where ¢ is a distance map to the ROI (the force is switchedoff inside the ROI). It should be noted that we also use the ROI as a mask on the GVF and thus the GVF is equal to zero outside the ROI. Finally, we also use the ROI as a mask on the obtained segmentation to ensure that no objects are overlapping. Segmentation results are shown in Figure 2.

Nonlinear registration procedure 4.1 Overview

The nonlinear registration procedure is separated into two stages: the first step consists of the registration between segmented structures in both images; and the second step, initialized with the precedently computed transformation, performs a refined registration between the whole intensity images.

Deformation model

A nonlinear transformation based on B-Spline Free Form Deformations (FFD), introduced by Sederberg et al. [START_REF] Sederberg | Free form deformation of solid geometric models[END_REF], has been chosen to establish the correspondence between images. The choice of this method over other more constrained parametric models is due to the great variability of the structures in our application. So, we preferred the flexibility that FFDs get from the fact that no assumptions on the structures are made. On the other hand, the speed requirements of the system make FFDs preferable to other more realistic and time consuming deformation frameworks, such as elastic or fluid models.

In this technique, deformations of the object volume are achieved by tuning an underlying mesh of control points. The number of control points would define the locality of the deformations allowed by the FFD model. On the other hand, convergence times will notably increase with a more densely populated grid. Thus, a trade-off concerning these two aspects must be taken, our choice being a mesh of 10 control points for dimension.

FFD model has been used for both stages of the registration procedure. This implementation allows us to easily integrate the two phases into the same framework, structures registration being considered as just an extra step of the multiresolution chain used in most voxel-based techniques. Thus, fine registration procedure will start with the grid found in the structure registration phase, which provides an initial transformation very close to the final solution, at least in the neighborhood of the segmented structures.

The optimization procedure is based on an iterative gradient descent technique over the entire grid of control points. At each iteration, we compute a local gradient estimation for each control point by finite differences. Furthermore, a local spring force regularization term has been added to prevent the nodes from intersecting, which could lead to unwanted alterations of the structure topology.

Structure registration

The aim of this step is to provide an initialization to the grey-level registration as close as possible to the desired final result. This transformation will constrain the search of the global solution that will undergo the next stage. Therefore, this registration phase can be seen as a multiresolution step, filtering out of all the data but the main structures, then forwarding the result as an initial estimate to a higher level where finer detail will be considered.

Before performing nonlinear registration between anatomical features, we compute an initial approximation of the transformation between both images. This approximation is composed of a rigid motion, an independent scaling along the three axes and the elimination of the parts of the volumes that have no interest for our application.

The FFD framework implies tuning the control points of the grid to minimize a given similarity criterion. The choice of this criterion is straightforward in our case, as we are working with segmented images with a linear intensity relation. Thus, the Root Mean Square (RMS) difference of corresponding pixel grey levels, summed across the whole volume, will be used to determine the optimal deformation parameters.

An overlap measure consisting of the quotient between intersection and union among structures (being equal to 1 if total superimposition is achieved) has been used to evaluate CT and nonlinear registered PET recognized features. We have obtained a value of 0.658 for the linear phase and a value of 0.903 for the nonlinear phase, clearly improving structure registration results, as can be seen in Fig. 3.

Fine registration

Transformation produced by the initialization step is not necessarily valid for those regions away from the segmented structures, so the computation of their displacement has to be done by this fine registration stage. Another objective of this stage is the correction of the errors that may have been introduced by the structure segmentation procedure, taking advantage of the image grey level information.

As the deformation framework has not changed, the algorithm is essentially the same as the one used in the structure registration stage. Nevertheless, the fact of working with the whole image intensity levels forces us to change the similarity criterion to maximize. The choice of a similarity measure is strongly related to the imaging modalities to be registered. A particularly complex situation arises when the intensity distributions of two different modalities do not follow a functional relation, such as in our case, implicating PET and CT images. Mutual Information (MI), a criterion proposed by Viola [START_REF] Viola | Alignment by maximization of mutual information[END_REF], has been demonstrated to be a powerful tool for multimodal image registration with a nonlinear intensity relation. MI expresses how much information from an image I is contained in another image J. Therefore, Mutual Information will be maximal if the images are geometrically aligned. NMI is a variant of MI, introduced by Studholme [START_REF] Studholme | An overlap invariant entropy measure of 3D medical image alignment[END_REF] to prevent the actual amount of image overlap from affecting the measure. Its computation requires an estimation of the marginal and joint probability distributions from both images. We use a frequency-based approximation, ¢¡ 1¤ ¤£ ¢¡ ¦¥ ¦£ , where §¡ is the estimated probability of having an intensity i in one image and j in the other, £ ¢¡ being the number of voxels with these inten- sities and n the total number of voxels.

Results and conclusions

We tested our method on a set of 15 CT and PET scans of the thoracic and abdominal regions. Some slices of the superimposition between CT and registered PET image can be seen in Figure 4 (All slices can be found at: http://www.tsi.enst.fr/ ocamara/Evaluation3.html). We have devised an evaluation protocol which has allowed a group of three clinicians (Dr.Foehrenbach, Dr. Rigo and Dr. Marchandise) to assess the registration results using an online procedure. The measure generated provides an error less than 1cm over the most relevant structures (lungs, liver, kidneys, heart), which conforms to the objective of the application (errors smaller than PET resolution images).

We have verified that the registration is better achieved around segmented structures which have been already registered in the initialization stage, independently of its belonging to thoracic or abdominal zones. One exception is the errors found in the stomach because of its important movements between two acquisitions and the fact that no constraints have been imposed on it, thus NMI registration has not converged towards the correct registration. The fact of constraining the FFD deformation by means of an initialization stage has speed up the convergence of the fine registration step, the computational cost being reduced over 80%.

As a conclusion, the results presented in this paper indicate that our method can provide an useful tool for data analysis in thoracic and abdominal oncology applications. The nonrigidity in the imaged regions is effectively modeled by means of a Free Form Deformation (FFD), and satisfactory registration results can be obtained by minimizing a Normalized Mutual Information criterion, given a good enough initialization. A progressive segmentation method has been proposed to provide such an initialization, integrating it as a first step of the multiresolution procedure.

Further work will focus on allowing the initialization phase to assign a weight to the nodes of the FFD according to the mechanical properties of the underlying tissue.

Figure 1 :

 1 Figure 1: ROI used for the segmentation of the kidneys (brighter area). Left: CT image. Right: PET image.

Figure 2 :

 2 Figure 2: Top-left: Axial slice of GVF computed on a CT image (left) and detail (right). Bottom-left: Example of a simplex mesh. Middle: 3D rendering of segmented structures in CT image superimposed on a slice. Right: 3D rendering of segmented structures in PET image superimposed on a slice.

Figure 3 :

 3 Figure 3: Example of structure registration. Left: 3D rendering of PET segmented structures after rigid registration with slices from CT original image. Right: 3D rendering of PET segmented structures after nonlinear structure registration with slices from CT original image, illustrating the extent of the deformation.

Figure 4 :

 4 Figure 4: An axial and a coronal slice of registration results. Left: CT original image. Center: Registered PET image. Right: Superimposition of both images. All slices can be found at: http://www.tsi.enst.fr/ ocamara/Evaluation3.html

Acknowledgments

The authors would like to thank Dr. Hervé Foehrenbach, Dr. Pierre Rigo, Dr. Xavier Marchandise, Yves Martelli, the research team of S.H.F.J. in Orsay and the members of Segami Corporation for their contribution to this project. This work was partially supported by the French Ministry for Research (grant number 01B0267).