# Evaluation of the symmetry plane in 3D MR brain images

Alexander V. Tuzikov, Olivier Colliot, Isabelle Bloch

► To cite this version:

Alexander V. Tuzikov, Olivier Colliot, Isabelle Bloch. Evaluation of the symmetry plane in 3D MR brain images. Pattern Recognition Letters, 2003, 24 (14), pp.2219-2233. 10.1016/S0167-8655(03)00049-7. hal-01227827

HAL Id: hal-01227827

https://inria.hal.science/hal-01227827

Submitted on 13 Jan 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Evaluation of the symmetry plane in 3D MR
brain images

Alexander V. Tuzikov<sup>a</sup>, Olivier Colliot<sup>b</sup>, Isabelle Bloch<sup>b,*</sup>

<sup>a</sup> Research-Engineering Center of Informational Technologies, National Academy of Sciences of Belarus, Akademicheskaja 25, 220072 Minsk, Belarus

<sup>b</sup>Ecole Nationale Supérieure des Télécommunications, Département TSI, CNRS URA 820, 46 rue Barrault, 75634 Paris Cedex 13, France

**Abstract**

Although the brain is not perfectly symmetrical with respect to the mid-sagittal plane, the automatic detection of this plane and of the degree of symmetry is of interest for many anatomical and functional studies. We propose a method for detecting the best symmetry plane in 3D MR brain images. We express this problem as a registration problem and compute a degree of similarity between the image and its reflection with respect to a plane. The best plane is then obtained by maximizing the similarity measure. This optimization is performed using the downhill simplex method and is initialized by a plane obtained from principal inertia axes, which proves to be close to the global optimum. This is demonstrated on several MR brain images. The proposed algorithm is then successfully tested on simulated and real 3D MR brain images. We also investigated the influence of the optimization procedure control parameters on the computation speed and result precision. Preliminary results obtained on CT and SPECT images suggest that the method can be extended to other modalities.

*Key words:* Symmetry plane, symmetry measure, MR brain image, mid-sagittal plane, inertia axes, optimization

# 1 Introduction

The mid-sagittal plane is defined as the plane that best separates both brain hemispheres. The automatic detection of this plane in brain images such as magnetic resonance images

<footnote>* Corresponding author. Tel: +33 1 45 81 75 85, Fax: +33 1 45 81 37 94
Email addresses: tuzikov@mpen.bas-net.by (Alexander V. Tuzikov),
Olivier.Colliot@enst.fr (Olivier Colliot), Isabelle.Bloch@enst.fr (Isabelle Bloch).</footnote>
(MRI) is a useful task. It can be used to reorient images, for example for further alignment in the Talairach reference frame [1]. It can also be a first step in intra-subject registration [2] or serve as a basis for the study of dissymmetries in the brain [3].

Brain images possess a high degree of plane symmetry although they are not exactly symmetrical. Medical studies (see, for example, [4] for some corresponding references) demonstrate brain anatomical asymmetries in size of the frontal and occipital lobes, in surface of the planum temporale and planum parietale, in the anterior speech region. Symmetry is a geometric feature. Since real brains and their images do not exhibit exact symmetry, it is of interest to measure the degree of symmetry in them. It can be computed using the notion of symmetry measure [5] or a distance function used for comparing objects or images. A plane corresponding to the maximal value of the symmetry measure is called a symmetry plane. This plane is often considered as the first approximation of the mid-sagittal plane [6]. The constraint of obtaining a plane is useful for some of the foreseen applications, in particular for computing spatial relationships with this plane as a reference. Although the true separation between both hemispheres can usually be detected in MRI and CT images, this is not possible on PET images for instance. Estimating the plane to be found as a symmetry plane allows to use the same algorithm for several modalities, as will be shown in our experiments. 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 [7–12] 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 [13]. Principal axes of inertia were used in [14] to define the best symmetry plane. An algorithm for finding symmetry planes of 3D objects using extended Gaussian image representation was developed in [15]. Octree representation of 3D objects is used in [16] to compute their symmetry degree.

Recently several new algorithms were proposed for symmetry plane computation in 3D images of the brain [2, 1607, 17, 1630].

The cross-correlation between original and reflected images is used in [2] as a symmetry measure of the image with respect to a considered plane. Preliminary image filtering, smoothing and size reduction are used there to reduce the computational complexity of finding the best symmetry plane. The optimization procedure includes the evaluation step for finding the initial plane position that is further improved using a downhill simplex optimization method. The evaluation step computes the symmetry measure for a number of plane positions equally distributed in the space of plane orientations.

A two step algorithm for computation of the symmetry plane in 3D brain images is proposed in [3]. It is assumed that some initial plane is already given. This initial plane can be either the middle plane of the image, generally a good starting point for the optimization procedure, or the plane computed based on principal inertia axes, if the brain is too tilted. The algorithm
allows automatically to reorient and re-center this plane. At the first step a point-to-point correspondence is established between two hemispheres using the demons algorithm. Then this correspondence is used for finding the new position of the plane minimizing the least squares criterion. The authors propose also a method to estimate the 3D dissymmetry field and apply it to the following problems: the study of the normal dissymmetry within a given population; the comparison of the dissymmetry between two populations; the detection of the significant abnormal dissymmetries of a patient with respect to a reference population. An improved version of the symmetry plane computation algorithm is presented in [17]. In the first step, the demons algorithm is replaced with a block matching. In the second step, the authors use a robust least trimmed squares criterion. Finally, the whole process is iterated. As it is reported the algorithm is efficient and achieves a good accuracy for anatomical and functional images. The method also works for pathological, highly asymmetrical brains.

A robust algorithm was developed in [1] for finding the symmetry plane for 3D normal and pathological neuroimages. The algorithm sequentially deals with preprocessed 2D slices and uses a cross-correlation criterion for finding lines of symmetry for every slice. The slice preprocessing includes smoothing, subsampling and edge detection. The global symmetry plane parameters are computed using a robust estimation of the information obtained from image slices processing. The authors compared the developed algorithm with one based on the maximization of mutual information registration, and reported the superior performance of the developed algorithm under tested conditions of noise and bias fields. However, it is not clear whether the authors applied a special preprocessing of the initial 3D images while using maximization of a mutual information criterion.

The problem of finding symmetry plane and symmetry measure can be treated as a con-
strained rigid registration problem between original and reflected images (with respect to
some fixed plane). Different similarity measures known in the literature [18] can be classified
depending on the type of registration: for monomodal registration one can use, for exam-
ple, the correlation coefficient [19], normalized cross-correlation, sum of absolute intensity
differences, sum of squared intensity differences, stochastic sign change [20]; for multimodal
registration one can use the correlation ratio [21], mutual information [22] or normalized
mutual information [23]. Our problem is a particular case of monomodal registration: both
images differ by a geometric transformation only and therefore have the same contrast. In the
following we use the normalized sum of squared intensity differences as a similarity measure.

The paper is organized as follows. In Section 2 we discuss the problem of plane symmetry evaluation in 3D images and the symmetry measure used. In Section 3 the problem is reformulated in terms of a constrained registration between the original image and its reflection with respect to a fixed plane. Then in Section 4 we propose an algorithm for the computation of the brain symmetry plane. The algorithm is very close to the one reported in [2]. The main difference is in the choice of the initial plane for the further optimization. In Section 5 we provide the results of investigating the profiles of the symmetry measure for several real MR brain images. The main observation is that these profiles have a regular form and that a plane computed from the ellipsoid of inertia belongs to the region of the global maximum of the similarity measure under consideration. Therefore, starting the optimization from this
plane with a local technique gives more chances to achieve the global maximum and besides reduces substantially the complexity of the algorithm. We also discuss in this section the results of processing simulated and real MR brain images, the properties of the proposed algorithm and the choice of its control parameters. The algorithm works directly on 3D images and not slice by slice as the one proposed in [1]. The algorithm is a natural extension of earlier methods that used the best plane obtained from principal axes of inertia as the symmetry plane (see for example [14]). The orientation and position of this plane is improved in our method by maximizing the correspondence between all voxels from its left and right parts.

## 2 Symmetry measure

In this paper we deal with 3D gray-level images, i.e. non-negative functions $f: F \subset \mathbb{R}^3 \to V$. Here $V$ is a subset of $\mathbb{R}_+$ or $\mathbb{N}$.

Given a point $\mathbf{b} = (b_1, b_2, b_3) \in \mathbb{R}^3$, we denote by $f_{\mathbf{b}}$ or by $f + \mathbf{b}$ the translation of image $f$ at point $\mathbf{b}$, i.e. $f_{\mathbf{b}}(x, y, z) = f(x - b_1, y - b_2, z - b_3)$.

Let $\mathbf{u}$ be a unit vector in $\mathbb{R}^3$, $S^2$ the unit sphere of all possible directions in $\mathbb{R}^3$, $S_+^2$ the hemisphere of the unit sphere containing points with non-negative $x$ coordinate, and $\Pi_{\mathbf{u},d}$ the plane in $\mathbb{R}^3$ orthogonal to the vector $\mathbf{u}$ and passing at the signed distance $d$ from the coordinate origin. We denote by $e_{\mathbf{u},d}(f)$ the reflection of image $f$ with respect to the plane $\Pi_{\mathbf{u},d}$:

$$e_{\mathbf{u},d}(f)(x, y, z) = f(e_{\mathbf{u},d}(x, y, z)).$$

**Definition 1** An image $f$ is called plane symmetrical if there exists a symmetry plane $\Pi_{\mathbf{u},d}$ such that $e_{\mathbf{u},d}(f) = f$. We say in this case that $\Pi_{\mathbf{u},d}$ is a plane of symmetry for $f$.

Note that images can have several symmetry planes. In practice we often deal with images that are almost plane symmetrical but do not possess the exact symmetry. In this case the notion of symmetry measure is very useful to access the degree of symmetry in images. Suppose that one is able to compute the symmetry measure $\mu_{\mathbf{u},d}(f)$ of the image $f$ with respect to an arbitrary plane $\Pi_{\mathbf{u},d}$. Then the problem we are dealing with in this paper is to find a plane of maximal symmetry measure of image $f$ and the corresponding value of the symmetry measure $\mu(f)$. It can be formulated in an equivalent form as (the search space can be restricted to $(S_+^2, \mathbb{R})$ or $(S^2, \mathbb{R}_+)$ since $\Pi_{\mathbf{u},d}$ and $\Pi_{-\mathbf{u},-d}$ define the same plane):

$$\mu(f) = \max_{\mathbf{u} \in S^2, d \in \mathbb{R}_+} \mu_{\mathbf{u},d}(f) = \max_{\mathbf{u} \in S_+^2, d \in \mathbb{R}} \mu_{\mathbf{u},d}(f). \qquad (1)$$

In the right part of the formula we use a hemisphere (as in our implementation) and therefore we use the signed distance between the symmetry plane and the coordinate origin. We use everywhere in this paper max instead of sup and assume that the maximum is achieved (which in fact is true under some continuity constraints on $\mu_{\mathbf{u},d}$).
When comparing images it is usually assumed that a distance or a similarity measure $\sigma(f_1, f_2)$ is defined on the set of images. In this case the symmetry measure $\mu_{u,d}(f)$ can be defined as the similarity between images $f$ and $e_{u,d}(f)$, i.e.:

$$ \mu_{u,d}(f) = \sigma(f, e_{u,d}(f)). \qquad (2) $$

One can see that $e_{u,d}(f)$ is a translation of $e_{u,0}(f)$ at the point $2du = (2du_1, 2du_2, 2du_3)$, i.e.:

$$ e_{u,d}(f) = e_{u,0}(f) + 2du. \qquad (3) $$

Therefore, one has:

$$ \mu_{u,d}(f) = \sigma(f, e_{u,0}(f) + 2du). $$

It is clear that for the computation of the symmetry measure $\mu(f)$ one needs to check all possible orientations and translations of the symmetry plane. It is an optimization problem in a three-dimensional parametric space. Here two oriented angles $\alpha$ and $\beta$ define a normal vector (see Fig. 1) to the symmetry plane and one parameter is needed for the signed distance to this plane. As $\mathbf{u}$ vary in $S_+^2$, $\alpha$ and $\beta$ will vary in $[-\pi/2, \pi/2)$. Later in Section 3 we discuss also how the problem of computing the symmetry measure and the corresponding plane can be reformulated in terms of a constrained registration between the original image and its reflection with respect to a fixed plane.

<figure><img src="image_2.png" /><figcaption>Fig. 1. The angles $\alpha$ and $\beta$ define a normal vector to the symmetry plane.</figcaption></figure>

Given two images $f$ and $g$, the $l_2$ metric is defined as follows:

$$ l_2(f, g) = ||f - g|| = \sqrt{\iiint_{\mathbb{R}^3} (f(x, y, z) - g(x, y, z))^2 dxdydz.} $$

We use the following symmetry measure in accordance with our view of the problem as a monomodal registration problem:

$$ \mu_{u,d}(f) = 1 - \frac{||f - e_{u,d}(f)||^2}{2||f||^2}. \qquad (4) $$
It holds $0 \le \mu_{u,d}(f) \le 1$ and $\mu_{u,d}(f) = 1$ if and only if $\Pi_{u,d}$ is a plane of symmetry for $f$. One can show that this measure coincides with the normalized cross-correlation.

### 3 Symmetry analysis in terms of registration

In this section we show that the computation of the symmetry plane is equivalent to computing a constrained rigid transformation between $f$ and its reflection with respect to a fixed plane. Denote by $r_{u,\alpha}$ a rotation in $\mathbb{R}^3$ about the oriented axis directed along vector $u$ over an angle $\alpha$ (counter-clockwise direction).

Given two unit vectors $u_1, u_2$ from the unit sphere $S^2$, denote by $\alpha(u_1, u_2)$ the angle in $[0, \pi]$ between them. Following the outline of [24] the composition of plane symmetries can be expressed as a rotation in the following way:

$$e_{u_2,0}e_{u_1,0} = r_{u, 2\alpha(u_1,u_2)}, \qquad (5)$$

where $u = \frac{u_1 \times u_2}{||u_1 \times u_2||}$ and $u_1 \times u_2$ denotes the outer product of vectors $u_1$ and $u_2$ (see Fig. 2).

<figure><img src="image_3.png" /><figcaption>Fig. 2. Vectors $u, u_1, u_2$ and angle $\alpha(u_1, u_2)$ used in formula (5).</figcaption></figure>

Since reflections $e_{u_2,d}$ and $e_{u_2,0}$ differ by a shift $2d\mathbf{u}_2$ along the vector $\mathbf{u}_2$ (3), then one gets:

$$e_{u_2,d} = r_{u, 2\alpha(u_1,u_2)}e_{u_1,0} + 2d\mathbf{u}_2.$$

Let us take for simplicity the vector $u_1 = (0, 0, 1)$ and denote by $f'$ the reflection of $f$ with respect to the $xy$ plane, i.e. $f' = e_{u_1,0}(f)$. Denote also by $S^1$ the great circle of $S^2$ orthogonal to $u_1$ (i.e. the equator). It is clear that in this case the vector $u$ always belongs to $S^1$. We will call the rotations about axes passing through the equator as feasible rotations. Therefore, every reflection $e_{u_2,d}$ of image $f$ can be represented as a feasible rotation of image $f'$ followed by a shift along vector $\mathbf{u}_2 = r_{u,\alpha}\mathbf{u}_1$ at the distance $2d$, i.e.,

$$e_{u_2,d}(f) = r_{u,2\alpha}(f') + 2d(r_{u,\alpha}\mathbf{u}_1).$$
To make $\mathbf{u}_2$ vary over $S^2$ is equivalent to make $\mathbf{u}$ vary over $S^1$ and $\alpha$ vary over $[0, \pi]$. Now the problem of symmetry measure calculation defined by equation (1) can be written as follows in terms of feasible rotations:

$$ \mu(f) = \max_{\mathbf{u} \in S^1, \alpha \in [0, \pi], d \in \mathbb{R}_+} \sigma(f, r_{\mathbf{u},2\alpha}(f') + 2d(r_{\mathbf{u},\alpha}\mathbf{u}_1)). \quad (6) $$

Therefore, one needs to find a rotation of image $f'$ about an axis belonging to the $xy$ plane followed by a translation that maximizes the similarity measure with $f$. Suppose that the maximum in (6) is achieved for $(\mathbf{u}^*, \alpha^*, d^*)$ then the corresponding symmetry plane is defined by the normal vector $\mathbf{u}_2 = r_{\mathbf{u}^*, \alpha^*} \mathbf{u}_1$ and the distance $d^*$ from the coordinate origin. Note that a similar idea to find a symmetry axis in 2D slices of a 3D image was used in [1].

We would like also to note that it can be shown as in [24] that every feasible rotation followed by a translation in a specified direction defines a reflection. The direction of the translation is defined by the rotation axis and angle.

In this paper we implemented the direct maximization of the symmetry measure given by equation (1). Equation (6) provides another option to compute the symmetry measure. However, in this paper it has more a theoretical than a practical value and it is presented here just to demonstrate another parameterization for computing symmetry measures.

## 4 Symmetry plane computation

The computation of symmetry measure given by equation (1) is a time consuming non convex optimization problem. Therefore it seems useful to check first the most appropriate positions of symmetry planes. These positions can be used as initial points for the further optimization. In this section we first give the details concerning the computation of initial points. Then we discuss the optimization method used and, finally, formulate an algorithm proposed for the computation of symmetry measure and plane.

The initial positions can be interpreted, for example, in terms of ellipsoid of inertia [25]. In this case we consider an image $f$ as a 3D dimensional body $F$ with a mass distribution $f(x, y, z)$ at every point $(x, y, z) \in F$.

The axes of the ellipsoid of inertia for a 3D body are defined by the eigenvectors of the corresponding covariance matrix:

$$ \begin{pmatrix} m_{200} & m_{110} & m_{101} \\ m_{110} & m_{020} & m_{011} \\ m_{101} & m_{011} & m_{002} \end{pmatrix} $$
Here $m_{pqr}$ defines a central moment of the order $p+q+r$:

$$m_{pqr}(f) = \iiint_{F} f(x, y, z)(x - x_{c})^{p}(y - y_{c})^{q}(z - z_{c})^{r}dxdydz,$$

where the integral is taken for the volume of the body and $(x_c, y_c, z_c)$ are the coordinates of the center of mass.

It is known that if a body possesses a plane symmetry then the symmetry plane passes through its center of mass orthogonal to some ellipsoid axis. Therefore, for the computation of a symmetry measure in this case it is sufficient to check 3 planes only. However, this is true when the eigenvectors are different. Otherwise, one gets an ellipsoid of revolution and its axes are not uniquely defined.

Suppose that the center of mass of the ellipsoid of inertia of image $f$ is $\mathbf{b}_c = (x_c, y_c, z_c)$ and that $\mathbf{u}_1, \mathbf{u}_2, \mathbf{u}_3$ are unit eigenvectors taken in each of the three distinct eigenspaces. Then by computing the symmetry measure for the three major planes of the ellipsoid of inertia one gets the following value $\mu'(f)$ which is a lower bound for the symmetry measure $\mu(f)$:

$$\mu'(f) = \max\{\mu_{\mathbf{u}_1, \mathbf{u}_1 \cdot \mathbf{b}_c}(f), \mu_{\mathbf{u}_2, \mathbf{u}_2 \cdot \mathbf{b}_c}(f), \mu_{\mathbf{u}_3, \mathbf{u}_3 \cdot \mathbf{b}_c}(f)\}, \quad (7)$$

where $\mathbf{u} \cdot \mathbf{b}$ denotes the inner product of vectors $\mathbf{u}$ and $\mathbf{b}$.

The plane passing through $\mathbf{b}_c$ and orthogonal to one $\mathbf{u}_i$, $i \in \{1, 2, 3\}$, such that $\mu'(f) = \mu_{\mathbf{u}_i, \mathbf{u}_i \cdot \mathbf{b}_c}(f)$ is chosen as an initial symmetry plane. Its position and orientation are further improved using the optimization in the 3D space of symmetry plane parameters. In our implementation we used the *downhill simplex method* [26]. The method is convenient in cases when one does not have the function derivatives. It requires only function evaluations. As it is reported in the literature [27] the method is quite accurate and robust under good starting points. For the optimization in a 3D space the method needs 4 non-coplanar points as an input. We used 4 points that define a tetrahedron with one vertex at the initial point. Since the downhill simplex method is a local optimization method there is no guarantee in general of finding the global maximum. However, in Section 5 we demonstrate experimentally that one can obtain good solutions for MR brain images using this approach.

The proposed algorithm consists of two steps.

**Algorithm for symmetry plane and measure computation:**

(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.
# 5 Experiments and discussion

## 5.1 Evaluation of the initialization

We investigated the graphs of $\mu_{u,d}(f)$ for different values of $u$ and $d$ for 5 real MR normal brain images. These graphs have a similar form and one of them is shown in Fig. 3. Fig. 3(a) shows the profile of $\mu_{u,d}(f)$ for the image presented in Fig. 5 for symmetry planes passing through the image center of mass. It is assumed that the coordinate origin is in the image center of mass. However, we do not align this coordinate system with the principal axes of ellipsoid of inertia, it is only a translation of the original one. Parameters $\alpha, \beta, d$ defining the orientation and position of the symmetry plane are given in this coordinate system. The function values are computed for angles $\alpha, \beta \in [-\pi/2, \pi/2)$ with an incremental step 0.05. The function graph has a regular form with a maximum in the neighborhood of the initial point. One can see in Fig. 3(c) that the initial position is in the neighborhood of this maximum. Fig. 3(b) shows the profile of $\mu_{u,d}(f)$ for angles $\alpha, \beta \in [0, 0.1]$ and a neighborhood of the center of mass $d \in [d_c - 5, d_c + 5]$. One can see that the maximum of the symmetry measure is in a neighborhood of the image center of mass. Summarizing we can conclude that in the tested MR brain images the initial position of the plane obtained from the ellipsoid of inertia belongs to the region of the global maximum of the symmetry measure. It allows us to go up to the function global maximum using the local optimization technique.

## 5.2 Evaluation on simulated MRI

The algorithm was tested on a public available 3-D simulated brain database from the McConnell Brain Imaging Center at McGill University (http://www.bic.mni.mcgill.ca/brainweb/). 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 [28–30]. We tested the algorithm on 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 non-uniformity 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 the $yz$ plane. For all these images the orientation and position of the best symmetry plane found by the algorithm coincide with the $yz$ plane. Under the assumption that for these images the mid-sagittal plane and the symmetry plane coincide, this shows the success of the algorithm in finding the symmetry plane. The maximal absolute difference $md(|\alpha|)$ for $\alpha$ is 0.42 degree, the mean absolute value $mv(|\alpha|) = 0.16$ degree and the standard deviation $std(|\alpha|) = 0.13$ degree. The values for $\beta$ equal $md(|\beta|) = 0.58$, $mv(|\beta|) = 0.16$ and $std(|\beta|) = 0.08$ degrees, respectively. This confirms that the proposed algorithm is stable in the presence of noise and image non-uniformity. In Fig. 4 we demonstrate as an example the best symmetry planes found for images T1-1-3-20, T2-1-7-40 and PD-3-1-0. Here the
<figure><img src="image_4.png" /><figcaption>Fig. 3. The graphs of function $\mu_{u,d}(f)$ computed for different values of $u$ and $d$ for the image shown in Fig. 5: (a) symmetry planes pass through the image center of mass ($d=0$) and $\alpha, \beta \in [-\pi/2, \pi/2)$ with an incremental step 0.05; (b) symmetry planes pass in the neighborhood of 5 voxels of the image center of mass and $\alpha, \beta \in [-0.05, 0.05]$ with an incremental step 0.005, i.e. in our case $\alpha$ and $\beta$ are in the neighborhood of the global maximum; to simplify the visualization every pair $(\alpha, \beta)$ is shown as one point on the axis and for every fixed $\alpha$ values of $\beta$ are running from -0.05 to 0.05; (c) Isolevels of the graph plotted in (a); + and * indicate the initial and final orientation (position in the $\alpha, \beta$ space) of the symmetry plane obtained after the optimization step, respectively.</figcaption></figure>

first part of the notation denotes the image modality, the second one - the slice thickness, the third one - the level of noise and the last one - the level of intensity non-uniformity.

## 5.3 Evaluation on real images

The algorithm was also tested on 24 real MR normal brain images (in T1 modality) and for all of them we got visually good final results (see Fig. 5 and 6). We applied the procedure
<figure><img src="image_5.png" /><figcaption>Fig. 4. Symmetry planes found for 3D simulated brain images from the McGill data base (coronal and axial slices): (a),(b) image T1-1-3-20; (c),(d) image T2-1-7-40; (e),(f) image PD-3-1-0.</figcaption></figure>
directly on the original image (without any preprocessing or segmentation) and obtained visually good final results (see one of them in Fig. 7). For 5 images we compared the position of the symmetry plane found for the original image and for the corresponding segmented brain. The segmentation of the brain from outer layers was obtained by a 3D mathematical morphology algorithm [31]. The final results were similar. The mean value and the standard deviation of absolute differences for α and β were mv(|α₁ − α₂|) = 0.5, std(|α₁ − α₂|) = 0.39 and mv(|β₁ − β₂|) = 0.76, std(|β₁ − β₂|) = 0.73 degrees 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 6 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 similar (mv(|α₂ − α₃|) = 0.59, std(|α₂ − α₃|) = 0.46 and mv(|β₂ − β₃|) = 0.61, std(|β₂ − β₃|) = 0.58 degrees assuming that index 3 corresponds to a brain mask image). It means that the brain shape defines the orientation of the symmetry plane with a high precision. In routine applications of the algorithm, we used directly the gray-level images without any preprocessing. However, one can use a binary mask for finding the symmetry plane for brain images with a pathology when there is a substantial asymmetry in gray-level values.

<figure><img src="image_6.png" /><figcaption>Fig. 5. Comparison of the best symmetry plane found for segmented and non-segmented brains: (a),(b) the initial symmetry plane of the segmented brain image computed from the ellipsoid of inertia; (c),(d) the best symmetry plane of the segmented brain image computed by the algorithm. (e),(f) the initial symmetry plane of the original image computed from the ellipsoid of inertia; (g),(h) the best symmetry plane of the original image computed by the algorithm.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 6. Best symmetry plane found by the algorithm for a segmented brain image. 3D images have been visualized on the Anatomist software (www.anatomist.info), developed at S.H.F.J, Orsay.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Several coronal and axial slices of the best symmetry plane found by the algorithm for a non-segmented brain image.</figcaption></figure>

5.4 Evaluation on pathological MR images and on images from other modalities

Additionally, we tested our procedure on 7 MR images with different kinds of pathologies. Four of them are presented in Fig. 8. The results are visually good. It would be interesting to
perform additional tests with pathologies of increasing size in order to evaluate the robustness
of the algorithm with respect to asymmetries. On the first image (Fig. 8 (a)), the automatic
initialization failed and the result was obtained thanks to a manual initialization using the
middle plane of the image. We also tested the algorithm for other modalities than MRI.
Fig. 9 presents results obtained on a CT scan and on four SPECT images. Good results were
obtained in both cases.

<figure><img src="image_9.png" /><figcaption>Fig. 8. Results on pathological MR images. (a) Thalamic haematoma. (b) Resection of the anterior pole of the temporal lobe. (c) Paramedian frontal meningioma. (d) Paramedian tumoral lesion.</figcaption></figure>

5.5 *Study of the optimization procedure*

We investigated the influence of the optimization procedure control parameters on the speed of the algorithm and the final result precision. The precision is defined by the symmetry measure value and the speed by the number of iterations in the optimization procedure. The control parameters of the downhill simplex algorithm include its stopping condition and the size of the input simplex (tetrahedron) defined in (α, β, d) space [26]. Our experiments were performed on 24 real images.

Let us denote by $\mu_1 \ge \mu_2 \ge \mu_3 \ge \mu_4$ the symmetry measure values corresponding to the simplex vertices. Every simplex vertex defines a symmetry plane and the symmetry measure is computed for this plane. The following termination criterion was used: the symmetry measure value for the final output simplex is fractionally smaller than some tolerance $f_{tol}$, i.e. $2(\mu_1 - \mu_4)/(|\mu_1| + |\mu_4|) < f_{tol}$. It is a standard criterion for the downhill simplex algorithm.

Let us denote the size parameters of the input simplex by $\Delta\alpha$, $\Delta\beta$, $\Delta d$ which define the
possible absolute difference in the simplex vertex coordinates and denote $(\alpha_1, \beta_1, d_1)$ the
parameters corresponding to the initial plane. We investigated three different inputs:
<figure><img src="image_10.png" /><figcaption>Fig. 9. Results on images from other modalities. (a) CT. (b) DaTSCAN SPECT of a patient with Parkinson's Disease. (c) ⁹⁹m-Tc-Bicisate Perfusion SPECT of a patient with a paramedian frontal meningioma. (d) ⁹⁹m-Tc-HMPAO Perfusion SPECT of a patient with a right parietal ischemia (e) Thallium-201 SPECT of a patient with a paramedian tumoral lesion.</figcaption></figure>

* a single simplex with the initial position at one of its vertices (its vertices have the following coordinates $M_1(\alpha_1, \beta_1, d_1)$, $M_2(\alpha_1 + \Delta\alpha, \beta_1, d_1)$, $M_3(\alpha_1, \beta_1 + \Delta\beta, d_1)$ and $M_4(\alpha_1, \beta_1, d_1 + \Delta d)$);

* a single simplex covering the neighborhood of the initial position (its vertices have the following coordinates $M_1(\alpha_1 + \Delta\alpha, \beta_1, d_1 - \Delta d)$, $M_2(\alpha_1 - \Delta\alpha, \beta_1 + \Delta\beta, d_1 - \Delta d)$, $M_3(\alpha_1 - \Delta\alpha, \beta_1 - \Delta\beta, d_1 - \Delta d)$ and $M_4(\alpha_1, \beta_1, d_1 + \Delta d)$);

* 8 simplexes covering the neighborhood of the initial position (the simplexes contain the initial position in one common vertex, one simplex is from the first input and the others
are generated in a symmetrical way).

We have not observed substantial differences between the first and the second types of input simplexes, but we prefer to use the first one, since in this case one has a guarantee to obtain a solution that is not worse than the initial plane. The third type of input is naturally more time-consuming, and in our experiments it did not result in a significant improvement of the symmetry measure in comparison to the first one. The experiment results were compared for all images with the references computed with parameters $f_{tol} = 10^{-4}$, $\Delta\alpha = \Delta\beta = 17.2$ degrees and $\Delta d = 10$ voxels that were visually good. We obtained, as it was expected, that the number of iterations decreases with increasing $f_{tol}$. Table 1 and Fig. 10 show that, even in the worst case, the use of $f_{tol} = 10^{-2}$ instead of $10^{-4}$ leads to no significant difference for the final symmetry plane (the maximal difference in this case is 0.5 degrees for angles $\alpha$ and $\beta$ and 0.0026 for the symmetry measure value). Therefore, our default value $f_{tol} = 10^{-2}$ is a compromise between the number of iterations and final result precision.

<caption>Table 1</caption>

Influence of $f_{tol}$ on the number of iterations and result precision. Here 0 corresponds to the reference symmetry plane obtained for $f_{tol} = 10^{-4}$ and $|\mu - \mu_0|$ denotes the absolute difference between the symmetry measures computed for the reference plane and planes obtained for other values of $f_{tol}$. Mean and maximal differences are displayed in the table for the dataset of 24 images.

<table><thead><tr><th></th><th></th><th>10<sup>-4</sup></th><th>5 &times; 10<sup>-4</sup></th><th>10<sup>-3</sup></th><th>5 &times; 10<sup>-3</sup></th><th>10<sup>-2</sup></th><th>5 &times; 10<sup>-2</sup></th><th>10<sup>-1</sup></th></tr></thead><tbody><tr><td colspan="2">Mean number of iterations</td><td>100.3</td><td>78.5</td><td>68.5</td><td>50.8</td><td>46.2</td><td>31.2</td><td>26</td></tr><tr><td rowspan="2">|μ - μ<sub>0</sub>|</td><td>mean</td><td>0</td><td>0.00001</td><td>0.0001</td><td>0.0006</td><td>0.0009</td><td>0.0052</td><td>0.0099</td></tr><tr><td>max</td><td>0</td><td>0.0002</td><td>0.0007</td><td>0.0026</td><td>0.0026</td><td>0.0177</td><td>0.0416</td></tr><tr><td rowspan="2">|α - α<sub>0</sub>|</td><td>mean</td><td>0</td><td>0.0088</td><td>0.0170</td><td>0.0716</td><td>0.0980</td><td>0.4547</td><td>0.4719</td></tr><tr><td>max</td><td>0</td><td>0.0394</td><td>0.0741</td><td>0.3403</td><td>0.3403</td><td>1.6210</td><td>2.0016</td></tr><tr><td rowspan="2">|β - β<sub>0</sub>|</td><td>mean</td><td>0</td><td>0.0042</td><td>0.0198</td><td>0.0925</td><td>0.1271</td><td>0.4534</td><td>0.7843</td></tr><tr><td>max</td><td>0</td><td>0.0159</td><td>0.1951</td><td>0.2624</td><td>0.5327</td><td>1.7288</td><td>1.7288</td></tr></tbody></table>

The absolute difference for angles $\alpha$ and $\beta$ between the initial and final planes was in our experiments always smaller than 5 degrees and the difference for $d$ was always smaller than 1. The mean values were respectively 1.4 for $\alpha$, 1.8 for $\beta$ and 0.45 for $d$. Therefore, the best symmetry planes found by the algorithm are close to the initial plane. We tested the algorithm for the following values $\Delta\alpha = \Delta\beta = 17.2, 11.5, 5.7, 2.9, 1.7, 0.6, 0.3$ degrees, $\Delta d = 10, 5, 2, 1$ and $f_{tol} = 10^{-4}$. The mean number of iterations for different sizes of the input simplex is shown in Fig. 11. The minimal number of iterations is obtained for $\Delta\alpha = \Delta\beta = 2.9$ degrees and $\Delta d = 1$. For smaller or higher values, the number of iterations increases. It is likely that the downhill simplex algorithm has an optimal number of iterations for a given size of the input simplex which should be related to the difference between the initial point and the position of the global maximum. This means that the previous values approximately correspond to the *distance* between the symmetry plane and the initial plane. For these
<figure><img src="image_11.png" /><figcaption>Fig. 10. Results of the optimization procedure for different values of $f_{tol}$ (other parameters are $\Delta\alpha = \Delta\beta = 0.3$ radians, $\Delta d = 10$ and index 0 corresponds to the reference symmetry plane): (a) mean number of iterations; (b) mean and maximum values of the difference for symmetry measure with respect to the reference; (c),(d) mean and maximum values of the absolute difference for $\alpha$ and $\beta$.</figcaption></figure>

parameters, the symmetry measure value is close to the reference experiment (the maximal difference is 0.004).

Therefore, we decided to use the following control parameters: $f_{tol} = 10^{-2}$, $\Delta\alpha = \Delta\beta = 0.05$ radians (about 2.9 degrees) and $\Delta d = 1$. For these parameters the symmetry measure never decreased of more than 0.0054 from the reference result. The results were visually as good as those for the references. The complete procedure performs about 30 iterations and takes about 3 minutes for a 128x128x62 image on a SUN Ultra5 350 MHz workstation.
<figure><img src="image_12.png" /><figcaption>Fig. 11. The mean number of iterations in the optimization procedure for $f_{tol} = 10^{-4}$ and different values of $\Delta\alpha, \Delta\beta, \Delta d$.</figcaption></figure>

# 6 Conclusion

We proposed an algorithm for computing the best symmetry plane and the plane symmetry measure for 3D images. The initial position of the plane is computed from the ellipsoid of inertia of the image. 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 with respect to noise and intensity non-uniformity. However, it is sensitive to the position of the initial plane used in the optimization procedure.

An obvious limit of the proposed approach, shared with the other approaches of the literature, is that if the head is only partially contained in the volume (in particular if a large part of one hemisphere is missing), then the method will not provide good results, since there is no more symmetry in the image.

The algorithm achieved good results on all real MR images of our experiments. We also performed some tests on pathological MR images which were satisfying. However, it is likely that the symmetry plane will deviate from the mid-sagittal plane in case of highly asymmetrical pathologies and it would be interesting to test quantitatively the robustness with respect to asymmetrical lesions. Although no extensive experiments have been performed, preliminary results obtained on CT and SPECT images suggest that the method can be extended to other modalities. We also investigated the control parameters of the optimization procedure to get a compromise between the algorithm speed and the result precision.

Further extensions of this work may include the evaluation of the influence of various simi-
larity measures on the position of the best symmetry plane. We also plan to use this plane
as a feature for structure recognition in MR images, for example, using its orientation to
define directional relationships between objects [32,33].

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. We would like
to thank CHU La Pitié-Salpêtrière (Paris), CHU La Timone (Marseille), Hôpital Militaire du
Val-de-Grâce (Paris), Service Hospitalier Frédéric Joliot (Orsay) INSERM 280 and Hôpital
Neurologique (Lyon) for providing the images.

We would like to thank very much the reviewers for suggestions that have greatly improved
the paper.

References

[1] Y. Liu, R. Collins, W. Rothfus, Robust midsagittal plane extraction from normal and pathological 3D neuroradiology images, IEEE Transactions on Medical Imaging 20 (3) (2001) 175-192.

[2] B. Ardekani, J. Kershaw, M. Braun, I. Kanno, Automatic detection of the mid-sagittal plane in 3D brain images, IEEE Transactions on Medical Imaging 16 (6) (1997) 947-952.

[3] J.-P. Thirion, S. Prima, G. Subsol, N. Roberts, Statistical analysis of normal and abnormal dissymmetry in volumetric medical images, Medical Image Analysis 4 (2) (2000) 111-121.

[4] K. Amunts, L. Jäncke, H. Mohlberg, H. Steinmetz, K. Zilles, Interhemispheric asymmetry of the human motor cortex related to handedness and gender, Neuropsychologia 38 (2000) 304-312.

[5] B. Grünbaum, Measures of symmetry for convex sets, in: Proc. Sympos. Pure Math., Vol. 7, Providence, USA, 1963, pp. 233–270.

[6] F. Kruggel, D. von Cramon, Alignment of magnetic-resonance brain datasets with the stereotactical coordinate system, Medical Image Analysis 3 (2) (1999) 175–185.

[7] M. J. Attalah, On symmetry detection, IEEE Transaction on Computer 34 (1985) 663–666.

[8] B. A. deValcourt, Measures of axial symmetry for ovals, Bull. Amer. Math. Soc. 72 (2) (1966) 289-290.

[9] G. Marola, On the detection of the axes of symmetry of symmetric and almost symmetric planar images, IEEE Transactions on Pattern Analysis and Machine Intelligence 11 (1989) 104–108.

[10] T. Masuda, K. Yamamoto, H. Yamada, Detection of partial symmetry using correlation with rotated-reflected images, Pattern Recognition 26 (1993) 1245–1253.
[11] C. Sun, Symmetry detection using gradient information, Pattern Recognition Letters 16 (9) (1995) 987-996.

[12] P. J. van Otterloo, A contour-oriented approach to digital shape analysis, Ph.D. thesis, Delft University of Technology, Delft, The Netherlands (1988).

[13] H. Zabrodsky, S. Peleg, D. Avnir, Symmetry as a continuous feature, IEEE Transactions on Pattern Analysis and Machine Intelligence 17 (1995) 1154-1166.

[14] D. O'Mara, R. Owens, Measuring bilateral symmetry in digital images, in: Proceedings of TENCON'96 IEEE Conference: Digital Signal Processing Applications, Vol. 1, 1996, pp. 151-156.

[15] C. Sun, J. Sherrah, 3-D symmetry detection using the extended Gaussian image, IEEE Transactions on Pattern Analysis and Machine Intelligence 19 (2) (1997) 164-168.

[16] P. Minovic, S. Ishikawa, K. Kato, Symmetry identification of a 3D object represented by octree, IEEE Transactions on Pattern Analysis and Machine Intelligence 15 (5) (1993) 507-514.

[17] S. Prima, S. Ourselin, N. Ayache, Computation of the mid-sagittal plane in 3D brain images, IEEE Transactions on Medical Imaging 21 (2) (2002) 122-138.

[18] J. Maintz, M. Viergever, A survey of medical image registration, Medical Image Analysis 2 (1) (1998) 1-36.

[19] L. Brown, A survey of image registration techniques, ACM Computing Surveys 24 (4) (1992) 325-375.

[20] A. Venot, J. Lebruchec, J. Roucayrol, A new class of similarity measures for robust image registration, Computer Vision, Graphics and Image Processing 28 (1984) 176-184.

[21] A. Roche, G. Malandain, N. Ayache, Unifying maximum likelihood approaches in medical image registration, Journal of Imaging Systems and Technology 11 (1) (2000) 71-80.

[22] F. Maes, A. Collignon, D. Vandermeulen, G. Marchal, P. Suetens, Multimodality image registration by maximisation of mutual information, IEEE Transactions on Medical Imaging 16 (2) (1997) 187-198.

[23] C. Studholme, D. Hill, D. Hawkes, An overlap invariant entropy measure of 3D medical image alignment, Pattern Recognition 32 (1999) 71-86.

[24] A. Tuzikov, S. Sheynin, Symmetry measure computation for convex polyhedra, Journal of Mathematical Imaging and Vision 16 (1) (2002) 41-56.

[25] H. Goldstein, Classical Mechanics, Addison-Wesley, Reading, MA, 1950.

[26] W. Press, S. Teukolsky, W. Vetterling, B. Flannery, Numerical Recipes in C, 2nd Edition. Cambridge University Press, Cambridge, 1992.

[27] J. Bernon, V. Boudousq, J. Rohmer, M. Fourcade, M. Zanca, M. Rossi, D. Mariano-Coulart, A comparative study of Powell's and Downhill Simplex algorithms for a fast multimodal surface matching in brain imaging, Computerized Medical Imaging and Graphics 25 (2001) 287-297.
[28] D. Collins, A. Zijdenbos, V. Kollokian, J. Sled, N. Kabani, C. Holmes, A. Evans, Design and construction of a realistic digital brain phantom, IEEE Transactions on Medical Imaging 17 (3) (1998) 463–468.

[29] R.-S. Kwan, A. Evans, G. Pike, MRI simulation-based evaluation of image-processing and classification methods, IEEE Transactions on Medical Imaging 18 (11) (1999) 1085–1097.

[30] R.-S. Kwan, A. Evans, G. Pike, An extensible MRI simulator for post-processing evaluation, in: Visualization in Biomedical Computing, Vol. 1131 of Lecture Notes in Computer Science, Springer-Verlag, 1996, pp. 135–140.

[31] J. F. Mangin, V. Frouin, I. Bloch, J. Regis, J. Lopez-Krahe, From 3D magnetic resonance images to structural representations of the cortex topography using topology preserving deformations, Journal of Mathematical Imaging and Vision 5 (1995) 297–318.

[32] I. Bloch, Fuzzy relative position between objects in image processing: a morphological approach, IEEE Transactions on Pattern Analysis and Machine Intelligence 21 (7) (1999) 657–664.

[33] O. Colliot, I. Bloch, A. Tuzikov, Characterization of approximate plane symmetries for 3D fuzzy objects, in: Information Processing and Management of Uncertainty IPMU, Vol. 3, Annecy, France, 2002, pp. 1749–1756.