Brain Structure Segmentation from MRI by Geometric Surface Flow

We present a method for semiautomatic segmentation of brain structures such as thalamus from MRI images based on the concept of geometric surface flow. Given an MRI image, the user can interactively initialize a seed model within region of interest. The model will then start to evolve by incorporating both boundary and region information following the principle of variational analysis. The deformation will stop when an equilibrium state is achieved. To overcome the low contrast of the original image data, a nonparametric kernel-based method is applied to simultaneously update the interior probability distribution during the model evolution. Our experiments on both 2D and 3D image data demonstrate that the new method is robust to image noise and inhomogeneity and will not leak from spurious edge gaps.


INTRODUCTION
Thalamus is the relay center for nerve impulses in the brain. It mediates communication among sensory, motor, and associative brain regions. Axons from almost every sensory system connect here as the last site before the information reaches the cerebral cortex. Information received from the diverse brain regions is passed on to the cortex through the thalamus. Anatomically, thalamus is the largest, most internal structure of the diencephalon consisting of dual lobe masses of gray matter. It is located at the rostral end of the midbrain on each side of the third ventricle. Each lobe is about 4 centimeters. Motor nuclei of the thalamus receive signals from the striatum and cerebellum and project into the motor and premotor areas of the cerebral cortex. The thalamus play a major role in the regulation of consciousness, alertness, arousal, and attention and is thus considered part of the limbic system.
Thalamus segmentation has become more and more essential for a wide range of clinical and research applications. For example, thalamus changes in terms of volume and intensity are involved in a large number of diseases, such as schizophrenia, Parkinson's disease and multiple sclerosis, and so forth. Manual segmentation is very labor intensive and the result is not reproducible. On the other hand, discrete methods such as thresholding or region growing are not reliable because of the low contrast and discontinuous edges in the MRI images of thalamus.
In this paper, we present a new semiautomatic framework for thalamus segmentation based on the concept of geometric surface flow. Unlike previous methods with edge information only, we apply a nonparametric kernel-based method that can simultaneously update the interior region statistics along with the boundary shape evolution. By integrating both boundary and region information, the new method is robust to image noise and inhomogeneity and will not leak from spurious edge gaps.

ACTIVE CONTOURS AND MEDICAL IMAGE SEGMENTATION
Segmenting structures from medical images and reconstructing a compact geometric representation of these structures is difficult due to the sheer size of the datasets and the complexity and variability of the anatomic shapes of interest. Furthermore, the shortcomings typical of sampled data, such as sampling artifacts, spatial aliasing, and noise, may cause the boundaries of structures to be indistinct and disconnected. The challenge is to extract boundary elements belonging to the same structure and integrate these elements into a coherent and consistent model of the structure. Among various segmentation techniques, active contours/deformable 2 International Journal of Biomedical Imaging models have been very successful since their invention in late 80's. The mathematical foundation of deformable models stems from the confluence of geometry, physics, and approximation theory. Terzopoulos pioneered the theory of continuous (multidimensional) deformable models based on Lagrangian dynamics [1] and formulated deformation energies for generalized splines with controlled continuity [1]. Kass et al. [2] introduced the active contour model or "snake," a deformable model which is essentially a 2D spline that minimizes an internal deformation energy subject to external forces derived from images. Later, a voluminous literature on deformable models appeared in computer vision, and especially in the field of medical image analysis [3][4][5][6][7][8][9]. For a collection of seminal works, see [10,11].

Geometric surface flow
Our new thalamus segmentation framework is based on the concept of geometric surface flow. The general formulation of the geometric surface flow is the following initial-value dynamical system of nonlinear PDEs: where F is the speed function, t is the time variable, k and k are the surface curvature and its derivative at the point p, s 0 ( p) is the initial surface, n is the surface normal vector.
Equation (1) can be either directly provided by the user, or more generally, obtained as a gradient descent flow of the Euler-Lagrange equation of certain energy functionals by the calculus of variations. In general, there are two approaches to numerically simulate PDEs such as (1): explicit Lagrangian approach or implicit level-set approach. Implicit level-set methods [3,6] are becoming popular recently mainly because of their ease in handling topology changes, however since all the computations are conducted in a higher dimension, the computational cost tends to be quite expensive. In this paper, we take the explicit Lagrangian approach to simulate the surface flow. In particular, we apply the new framework we recently proposed [12,13]. In our new framework, the geometry and topology of the model are always explicitly represented throughout the simulation process. To ensure the regularity of the model and the stability of the numerical integration process, powerful Laplacian tangential smoothing, along with commonly used mesh optimization techniques, are employed throughout the geometric deformation and topological variation process. In addition, a novel particlebased collision detection scheme is conducted to automatically handle topology changes during the deformation. The new framework for surface flow simulation is fast, simple, and accurate. More importantly, it allows the user to directly interact with the model during the deformation. For more details of our new framework, please refer to [12,13].

Thalamus segmentation
When applying the above geometric surface flow for thalamus segmentation, the speed function F in (1) is explicitly formulated as a linear combination of the following two terms: here F reg is the regularity term to maintain the smoothness of the model. F reg is usually a function of the curvature (e.g., mean curvature, Gaussian curvature). In this paper, we define F reg as the difference between the mean curvature H curr at the current position and the average mean curvature H avg of the whole model: Comparing with other commonly used regularity terms such as mean curvature or Gaussian curvature, (3) works much better in our experiment. This is probably because of the volume preserving property of (3), which can be considered as a shape prior term for the thalamus's oval shape. A more detailed explanation can be found in [14], where a similar term is used. F data is the term that interacts the model with the image data, and is often defined as a function of the edge information of the image data. However, since the thalamus has low contrast, using boundary information alone is not very reliable. It is more preferable to combine the edge information along with the region information as suggested by Huang et al. [15]. Hence, in this paper, we propose to employ both the edge information as well as the region information. In particular, we define F data as a linear combination of the following two terms: F edge is the speed function corresponding to the edge information in the image data that will attract the contour moving towards the edges of the image data and is defined as where Δ is the Laplacian operator, I is the image intensity function, and G σ * I is the smoothed intensity function by convoluting with a Gaussian filter with variance σ. Variance σ can be assigned according to the image resolution of the MRI dataset. In our experiment, σ is set as 1.0. F region is the speed function corresponding to the region information, and is defined as Greg Heckenberg et al.

3
where ∇ is the gradient operator, B(s) is the binary image created by the interior probability estimation of the current model, which will be explained in more detail in the following section (Sections 3.3.2 and 3.3.3). In a nutshell, F region will expand the contour outwards in regions that are considered compatible with the current segmented volume's statistics. Substituting (3), (4), (5), and (6) into (2), we obtained the following speed function for thalamus segmentation: Here, a 1 , a 2 , a 3 are the corresponding weighting coefficients. In our experiments, a 1 , a 2 , a 3 are set as 0.1, 1, 1, respectively.

Seed initialization
First, the user interactively selects a pixel/voxel inside the region of interest in the image data. Considering the thalamus has an oval shape, a circular contour centered at the pixel/voxel is then automatically created and serves as the initial seed model.

Interior statistics estimation
Then, the intensity probability density function of the interior regions enclosed by the seed contour is estimated. Specifically, we approximate the distribution by a Parzenwindow-function-based nonparametric method [16] because it is differentiable, more generic, and can represent complex multimodal intensity distributions. We choose the Gaussian kernel as the Parzen window function. Suppose the model M is placed on an image I, the volume of the image region bounded by the current model M is V , then the probability of a pixel's (voxel's) intensity value i being consistent with the model interior intensity can be derived as where σ is a constant that specifies the width of the Gaussian kernel and is set as 1.0 in our experiment. Since (8) is a simple integral, it can be calculated very efficiently as an incremental update, that is, only newly added voxels are calculated to update the integral value at each time step.

Binary image creation
Next, based on the interior probability density distribution of model M obtained in the previous step, the image intensity probability map P I of every pixel's (voxel's) intensity is obtained. Then a small threshold (e.g., the mean probability over the entire image domain) is applied on P I to produce a binary image B(s), in which pixels/voxels with probability higher than the threshold have value 1, zero otherwise.

Model evolution
Finally, the model will start to evolve according to (1) and (7). More specifically, the surface evolution process is numerically approximated using a simple, explicit iterative equation: When advancing the model, we must enforce a constraint on the size of the time step Δt. In particular, the time step Δt must satisfy the following stability condition: the velocity of change must be strictly restrained by the minimum detail in the system. In our system, this condition is where m e is the unit grid cell length of the image data and M F is the maximum magnitude of the speed F obtained by (7). Before each deformation step, we will calculate the velocity F at each vertex point and determine the maximum magnitude of the velocity M F . A proper time step can then be estimated from (10). At each deformation cycle, the model will loop from step 2 to step 4 until an equilibrium state is achieved, then the deformation stops and the geometry of the thalamus is extracted from the MRI image.

Experimental results
Some of the experimental results for 2D and 3D thalamus segmentations are shown in Figure 1 to Figure 4. In particular, an example of our 2D thalamus segmentation result is shown in Figures 1 and 2. Each of the two figures shows five snapshots of the model evolution process including the initial seed (a) and the final shape (e). Figure 1 shows the model superimposed with the original image, while Figure 2 shows the same five snapshots superimposed with the binary image created by the interior statistics estimation (Sections 3.3.2 and 3.3.3). The 3D thalamus segmentation result is shown in Figures 3 and 4. Figure 3 shows the four snapshots during the model evolution process. The four different views of the final 3D shape of the thalamus are shown in Figure 4.

CONCLUSION
We proposed a semiautomatic thalamus segmentation method based on an explicit simulation of geometric surface flow. To overcome the low contrast of the thalamus, we employ a nonparametric kernel-based statistics estimation that can incorporate both the boundary and region information in the model evolution process, so that the model is very robust to image noise and small gaps. In the future we would also like to apply our method for the more challenging task of thalamus nuclei segmentation from the newly available diffusion tensor images.

ACKNOWLEDGMENT
The brain MRI dataset used in this paper is kindly provided by Dr. Otto Muzik from the PET Imaging Center, School of Medicine, Wayne State University, Detroit, Michigan.