Computation of the mid-sagittal plane in diffusion tensor MR brain images

TLDR: This plane is estimated as the one that best superposes the two hemispheres of the brain by reflection symmetry via the automated minimisation of a correlation-type global criterion over the tensor image via the NEWUOA algorithm in a multiresolution framework.

Abstract: We propose a method for the automated computation of the mid-sagittal plane of the brain in diffusion tensor MR images. We estimate this plane as the one that best superposes the two hemispheres of the brain by reflection symmetry. This is done via the automated minimisation of a correlation-type global criterion over the tensor image. The minimisation is performed using the NEWUOA algorithm in a multiresolution framework. We validate our algorithm on synthetic diffusion tensor MR images. We quantitatively compare this computed plane with similar planes obtained from scalar diffusion images (such as FA and ADC maps) and from the B0 image (that is, without diffusion sensitisation). Finally, we show some results on real diffusion tensor MR images.

Figures

Figure 2. Robustness and accuracy of the algorithm: δ = max(δi) measures its capture range, whereas = max( i) measures its accuracy.

Figure 1. Synthetic experiments. In a real DTI (matrix size = 256×256, field of view = 22.4cm, number of slices = 60, slice thickness = 2mm), the mid-sagittal plane of the brain is manually aligned with the centre of the image grid (a). Half of the image is then removed and replaced by the flipped version of the other half (b). In this new image, the central plane of the grid is a perfect symmetry plane for the brain. This image is then deformed by a rigid transformation R1 that moves each point p to p′ following: p′ = R1(p) = Rx(φx)Ry(φy)Rz(φz)p+ t, where Rx (resp. Ry and Rz) is a rotation around the x- (resp. y- and z-) axis and t = (tx, ty, tz) is a translation (c). We then apply our algorithm to compute the transformation R2 correcting the tilt (d). The three parameters φx, ty and tz leave the symmetry plane unchanged; they are set to 0. We lead 400 different experiments, with random parameters φy, φz and tx between ±20◦ and ±20 mm respectively.

Figure 4. Results on real DTI data. Two patients with Parkinson’s disease (left: original image and right: realigned image). The color coding is: red for left-right, green for antero-posterior and blue for top-bottom. It is applied to the principal diffusion of direction (eigenvector of the diffusion tensor with the largest eigenvalue), with the brightness proportional to the fractional anisotropy (FA).

Figure 3. Plot of the robustness/accuracy results for the 400 random experiments: horizontal axis = δ and vertical axis = ; the algorithm never fails when δ < 43 mm. The first error occurs for an experiment where δ = 43.82 mm, which corresponds to (tx, φy, φz) = (-18.26, 18.45, 2.68).

Full text

HAL Id: inserm-00140454
https://www.hal.inserm.fr/inserm-00140454
Submitted on 3 May 2011
HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entic 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 diusion de documents
scientiques 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.
Computation of the mid-sagittal plane in diusion
tensor MR brain images
Sylvain Prima, Nicolas Wiest-Daesslé
To cite this version:
Sylvain Prima, Nicolas Wiest-Daesslé. Computation of the mid-sagittal plane in diusion tensor
MR brain images. Medical Imaging 2007: Image Processing, Feb 2007, San Diego, United States.
pp.65121I, �10.1117/12.709467�. �inserm-00140454�

Computation of the mid-sagittal plane
in diffusion tensor MR brain images
Sylvain Prima
a,b,c,d
and Nicolas Wiest-Daessl´e
a,b,c,d
a
INRIA, VisAGeS U746 Unit/Project, IRISA, Campus de Beaulieu, F-35042 Rennes, France
b
University of Rennes I, IRISA, Campus de Beaulieu, F-35042 Rennes, France
c
CNRS, IRISA UMR 6074, Campus de Beaulieu, F-35042 Rennes, France
d
INSERM, VisAGeS U746 Unit/Project, IRISA, Campus de Beaulieu, F-35042 Rennes, France
ABSTRACT
We propose a method for the automated computation of the mid-sagittal plane of the brain in diffusion tensor
MR images. We estimate this plane as the one that best superposes the two hemispheres of the brain by reflection
symmetry. This is done via the automated minimisation of a correlation-type global criterion over the tensor
image. The minimisation is performed using the NEWUOA algorithm in a multiresolution framework. We
validate our algorithm on synthetic diffusion tensor MR images. We quantitatively compare this computed plane
with similar planes obtained from scalar diffusion images (such as FA and ADC maps) and from the B0 image
(that is, without diffusion sensitisation). Finally, we show some results on real diffusion tensor MR images.
Keywords: Registration, optimisation methods, diffusion tensor imaging, brain symmetry, mid-sagittal plane.
1. INTRODUCTION
The human brain displays an approximately bilateral symmetry with respect to the mid-sagittal plane (MSP).
This gross symmetry is often partially hidden in medical images, when the head is scanned in a tilted position.
In this situation, homologous bilateral structures (ventricles, grey nuclei, etc.) do not appear on the same 2D
axial or coronal slices, which can lead to misinterpretation of the images. Having a coherent display of medical
images with respect to the cerebral MSP can be useful for many applications.
First, it allows to remove the inter-hemispheric differences only due to the head tilt. After correction of the
tilt, the remaining normal or abnormal asymmetries can be assessed qualitatively or quantitatively, which can
be relevant for many pathologies. Second, the MSP itself is a crucial landmark for further analyses, such as
morphometry of the corpus callosum, delineation of the AC-PC line, or global quantification of brain symmetries
and asymmetries. Third, knowing the position of this plane is useful in registration tasks, where such an a priori
knowledge about the images to be matched can improve the robustness, accuracy and computation speed of the
methods.
Many algorithms have been developed for the automated computation of the cerebral mid-sagittal plane in
medical images. Two different approaches have been mainly used. Some works define the mid-sagittal plane
as the plane that best matches the inter-hemispheric fissure.
1–4
Such techniques are only applicable when this
fissure and its characteristic features are clearly visible on images, limiting their application to MR and CT
images. A more general approach is to view the mid-sagittal plane as the one that best separates the brain into
two similar parts. Many methods have been developed based on this idea, dealing with either anatomical (MR
and CT) or functional images (PET, SPECT).
5–12
Diffusion Tensor Magnetic Resonance Imaging (DT-MRI or DTI for short) complements conventional, anatom-
ical, MR imaging in that it adds specific structural information about the brain. By measuring diffusion properties
of water molecules in the brain, it allows the indirect visualisation of the microscopic organisation of tissues,
and more especially the white matter. There is an increasing interest in how DTI conveys brain symmetry and
asymmetry.
13–20
Further author information: (Send correspondence to Sylvain Prima)
Sylvain Prima: E-mail: sprima@irisa.fr, Telephone: 33 2 99 84 73 59
Medical Imaging 2007: Image Processing, edited by Josien P. W. Pluim, Joseph M. Reinhardt,
Proc. of SPIE Vol. 6512, 65121I, (2007) · 1605-7422/07/$18 · doi: 10.1117/12.709467
Proc. of SPIE Vol. 6512 65121I-1

In this paper, we propose a method for the automated computation of the mid-sagittal plane of the brain
in DTI. We estimate this plane as the one that best superposes the two hemispheres of the brain by reflection
symmetry. This is done via the automated minimisation of a correlation-type global criterion over the image
data, i.e. diffusion tensors. The minimisation is performed using the NEWUOA algorithm in a multiresolution
framework.
Such an algorithm has the same potential applications as the methods developed for anatomical and functional
images. More specifically, it has a direct application on DTI visualisation, which is often a critical task due to the
complex, non-scalar, nature of the data under study. For example, the classical red-green-blue colormap applied
to the principal direction of diffusion for easier display of DTI (red for left-right, green for antero-posterior and
blue for top-bottom)
21
is of poor interest if the mid-sagittal plane (and thus the left-right direction) is not aligned
with the image grid.
We detail our method in Section 2. We make a quantitative evaluation of its robustness and accuracy
properties, and present some results on real data in Section 3. We conclude and give some perspectives in
Section 4.
2. METHOD
2.1. Data acquisition
Images were obtained on a whole-body Philips 3T MR scanner. A single-shot echo-planar sequence combined
with coil sensitivity encoding (SENSE) was used to acquire the diffusion-weighted MR images with a spin-echo
Stejskal-Tanner sequence (TE/TR(ms) = 56/7991, flip angle = 90
◦
).
22
For each subject, an image without
diffusion sensitisation was acquired (b value = 0 sec/mm
2
), together with diffusion-weighted MR images (b value
= 800 sec/mm
2
) with diffusion-sensitising gradients applied in 16 different directions. The characteristics of the
diffusion-weighted MR images are: matrix size = 256 × 256, pixel size (mm) = 0.875 × 0.875, FOV (mm) =
224 × 224 × 120, number of slices = 60, slice thickness (mm) = 2.
2.2. Pre-processing
The pre-processing stage includes the following steps:
• Correction of distortions. Diffusion-weighted images acquired with echo-planar sequences are signifi-
cantly distorted, mainly due to the eddy currents induced by the large diffusion gradients. To deal with
this artefact, each diffusion-weighted image was registered to the image without diffusion sensitisation. The
alignment is based on the maximisation of a similarity measure on image intensities (the local correlation)
and a transformation model based on the imaging physics.
23
• Computation of a brain mask as the intersection of individual brain masks computed on each diffusion-
weighted MR image by histogram analysis.
24
• Estimation of the diffusion tensor by multiple linear regression on a voxel-by-voxel basis.
25
• Diagonalisation of the diffusion tensor by classical linear algebra (computation of the eigenvalues and
eigenvectors).
• Computation of different quantitative MR diffusion parameters and related 3D maps: Appar-
ent Diffusion Coefficient (ADC), Fractional Anisotropy (FA), Relative Anisotropy (RA), direction of the
eigenvector with the largest eigenvalue.
Proc. of SPIE Vol. 6512 65121I-2

2.3. Computation of the mid-sagittal plane
2.3.1. Some notations
• I is a diffusion tensor image, defined on a regular grid G
• For the sake of simplicity and without loss of generality we assume that G is isotropic with voxel size 1
mm
3
• v is a voxel on the grid G, and its value I(v) is a tensor (symmetric positive definite (SPD) 3x3 matrix)
• a(b) is the value of the function a for argument b, whereas ab is the product of matrices a and b;inthe
following, v is considered as either the argument of a function or as a column vector, depending on the
context
• S
P
is the reflection symmetry with respect to the plane P
• M
T
is the transpose of matrix M
• x, y and z are the left-right, antero-posterior and top-bottom axes, respectively
2.4. Formulation of the problem
If the brain was perfectly symmetrical, there would exist a symmetry plane P superposing each voxel v of the
grid with its anatomical counterpart. If we note v

= S
P
(v), it means that the tensors I(v)andI(v

) would be
mirror-symmetrical with respect to P , which writes:
I(v)=S
P
(I(v

))
In real life, the brain is only grossly symmetrical, and such an ideal symmetry plane does not exist. Instead,
we define the approximate symmetry (or mid-sagittal) plane P as the plane minimising the average distance
between the tensor I(v) and the reflection of I(v

) with respect to P , where v

= S
P
(v) is the voxel homologous
to v in the contralateral hemisphere. This average distance would be null for a perfect symmetry plane. The
optimisation problem can then be written as:
˜
P = arg min
P
1
card(G
P
)

v ∈ G
P
d (I(v),S
P
(I(v

))) , where:
• v

= S
P
(v)
• G
P
is the overlapping area between the two sets of voxels to compare: G
P
= G ∩ S
P
(G)
• d(., .) is a metric on tensors
2.5. Implementation details
2.5.1. Parameterisation of P
Any plane in IR
3
can be characterised by a normal unit vector n and its distance d to the origin. The unit vector
n is characterised by its angles α and β in a spherical coordinate system, where:
• α is the latitude angle (the equator is in the xy-plane)
• β is the longitude angle (the prime meridian is in the xz-plane)
Proc. of SPIE Vol. 6512 65121I-3

The coordinates of n are:
n
T
=

cos(β)cos(α)sin(β)cos(α)sin(α)

If we note H = Id − 2nn
T
(Id is the 3x3 identity matrix) and t =2dn, then it can be easily shown that the
reflection of the voxel v and of the tensor I(v) with respect to P are equal to, respectively:
S
P
(v)=Hv + t and S
P
(I(v)) = HI(v)H
T
= HI(v)H,asH = H
T
2.5.2. Interpolation scheme
For an arbitrary plane P , the voxel v

= S
P
(v) does not necessarily coincide with a grid point of G.Thus,an
interpolation is needed to estimate the tensor value I(v

) at this non-grid point. A natural choice would be to
interpolate in log-Euclidean or affine-invariant frameworks, which would preserve the constraint of symmetric
positive definiteness.
26–29
This unfortunately leads to a high computational cost and a prohibitive computation
time for the overall algorithm. Our choice is to simply use the nearest neighbor interpolation instead, which also
guarantees that the interpolated value is a SPD matrix.
2.5.3. Distance function
The distance function d is simply based on the Frobenius norm. For two tensors C and D, d is defined as:
d(C, D)=




3

i=1
3

j=1
(C
ij
− D
ij
)
2
2.5.4. Optimisation
We choose the algorithm NEWUOA (NEW Unconstrained Optimisation Algorithm) to solve the optimisation
problem. This algorithm, originally proposed by M.J.D. Powell,
30
has been shown to be significantly more
accurate, more robust and faster than other deterministic, derivative-free algorithms such as Nelder-Mead’s
downhill simplex or Powell’s direction set algorithms for intensity-based rigid-body image registration.
31
The
three parameters of interest are α, β and d (that define the plane P and thus the associated reflection symmetry
S
P
). The algorithm is initialised with a plane at the centre of the image grid, orthogonal to the left-right axis.
A multiresolution (simple decimation) approach is implemented to increase the robustness, accuracy and
speed of the algorithm. A subsampled image is built by taking only one voxel out of f
x
(resp. f
y
, f
z
)inthe
x-(resp. y-, z-) direction in the original image. A first solution is found at this low resolution, and is used to
initialise the algorithm at a higher resolution, and so on. In this coarse-to-fine scheme, typically two or three
levels can be used, depending on the resolution of the original image.
2.5.5. Realignment of the plane
Once the MSP is estimated, it is realigned at the centre of the image grid for improved image display.
7
Such
a realignment is performed by applying the rigid transformation R to both the anatomical and the tensor
information.
32
The rigid transformation R is computed as:
R =(S
K
◦ S
P
)
−1/2
where S
P
is the reflection symmetry with respect to the estimated mid-sagittal plane and S
K
is the reflection
symmetry with respect to K, the plane in the middle of the image grid. Let r be the rotation part of R.Each
tensor I(v) of the image I is rotated following: R(I(v)) = rI(v)r
T
. The log-Euclidean interpolation
28
is used for
computation of the realigned DTI. The overall computation time for a two-level multiresolution scheme is about
13 minutes for a typical DT image (matrix size 256 × 256 and 60 slices) on a standard PC (OS Linux), with an
Intel Xeon 4 CPU at 2.8GHz, 2 GBytes of RAM.
Proc. of SPIE Vol. 6512 65121I-4

Frequently asked questions

Q1. What contributions have the authors mentioned in the paper "Computation of the mid-sagittal plane in diffusion tensor mr brain images" ?

The authors propose a method for the automated computation of the mid-sagittal plane of the brain in diffusion tensor MR images. Finally, the authors show some results on real diffusion tensor MR images. 

Q2. What are the future works mentioned in the paper "Computation of the mid-sagittal plane in diffusion tensor mr brain images" ?

In this paper, the authors have proposed a method for the automated computation of the mid-sagittal plane of the brain in diffusion tensor MR images. In the future, the authors plan to experiment other implementation choices and compare them with the present method.