Segmentation of the hippocampus from brain MRI using deformable contours

Computerized Medical Imaging and Graphics

Pergamon Computerized Medical Imaging and Graphics 22 (1998) 203–216

Segmentation of the hippocampus from brain MRI using deformable contours Amir Ghanei

a,b,c,

*, Hamid Soltanian-Zadeh

a,b,d

, Joe P. Windham

a

a

Department of Diagnostic Radiology and Medical Imaging, Henry Ford Hospital, Detroit, MI 48202, USA b Department of Electrical and Computer Engineering, University of Tehran, Tehran 14399, Iran c Department of Electrical Engineering and Computer Science, University of Michigan, Ann Arbor, MI 48105, USA d Department of Radiology, Case Western Reserve University, Cleveland, OH 44106, USA Received 22 May 1997; revised 13 April 1998

Abstract The application of a discrete dynamic contour model for segmentation of the hippocampus from brain MRI has been investigated. Solutions to several common problems of dynamic contours in this case and similar cases have been developed. A new method for extracting the discontinuous boundary of a structure with multiple edges near the structure has been developed. The method is based on detecting and following edges by external forces. The reliability of the final contour and the model stability have been improved by using a continuous mapping of the external energy and limiting movements of the contour. The problem of optimizing the internal force weight has been overcome by making it dependent on the amount of the external force. Finally, the results of applying the proposed algorithm, which implements the above modifications, to multiple applications have been evaluated. q 1998 Elsevier Science Ltd. All rights reserved. Keywords: Dynamic contours; Segmentation; Medical imaging; Hippocampus; MRI

1. Introduction One of the most important stages in medical image analysis is segmentation of the objects or definition of their contours. Medical images are relatively hard to segment as a result of several undesirable properties such as low signal-to-noise and contrast-to-noise ratios, and multiple and discontinuous edges. Manual segmentation is time consuming, is not reproducible, and is not accurate due to human error. Standard edge detection methods are not appropriate either, because of objects discontinuous boundaries and low contrast-to-noise ratios. Deformable models seem more appropriate than conventional segmentation and edge detection methods, due to their relative power in treating each structure as a unit object, their flexibility, and producing a closed contour. In this article, we focus on the application of dynamic contours to the segmentation of the hippocampus in brain MRI. The hippocampus is a brain structure whose change from a normal state may be related to some of the important * Corresponding author. Tel: 001 313 876 1428; fax: 001 313 876 7925; email: [email protected]

0895-6111/98/$19.00 q 1998 Elsevier Science Ltd. All rights reserved. PII: S0 89 5 -6 1 11 ( 98 ) 00 0 26 - 3

brain diseases, such as epilepsy and schizophrenia. Segmentation of this structure from brain MRI is necessary for the volume analysis required for detection of epilepsy from MRI which is preferable to EEG tests that take nearly a week and are invasive. After segmentation, the shape and volume of the hippocampus can be estimated and compared with the natural shape and volume [1–5]. The hippocampus has relatively low contrast, multiple edges, and also discontinuous edges in MRI, making it relatively hard for automatic segmentation [6]. Although better than standard segmentation techniques, current methods for dynamic contours (e.g. [7–10]) have several problems when faced with such conditions. We discuss these problems and develop appropriate solutions for them. As a result, the goal of our work is to develop a tool to assist the operator with definition and segmentation of objects like the hippocampus from medical images. Such a tool should reduce the required user interaction significantly, thereby saving time and improving accuracy and reproducibility. The application of the method is not restricted to medical images, it can be used in computer graphics and animation or object segmentation in other types of images.

204

A. Ghanei et al./Computerized Medical Imaging and Graphics 22 (1998) 203–216

2. Background Dynamic contour models have been introduced by Terzopoulos, et al. [7] and have been investigated and applied in many ways, e.g. Kass et al. [11], Hill et al. [12]. Terzopoulos’s model is based on minimizing the energy of spline segments that make the contour. The energy consists of an internal energy which is the bending energy of spline and an external energy which is calculated by integrating image features; a user defined constraint energy can also be added. A considerable amount of research has been done concerning different aspects of deformable models. Different authors have used different methods for solving differential equations related to the contour evolution [13]. The models that use differential methods for minimizing the contour energy (e.g. [7,11]) are computationally inefficient and can lead to numerical instability, due to the discrete nature of the image. In contrast, different dynamic programming methods have been used by several authors. Geriger et al. [14] and Amini et al. [15] have used dynamic programming for tracking and matching the deformable contours and minimizing its energy. Grzeszczuk and Levin [8] used simulated annealing for evolution of the contour, to make the model stable. Poon et al. [16] have also used simulated annealing. They have incorporated statistical region-based image features to improve the reliability of the algorithm and used it in multi-region segmentation. In Sebbahi’s model [17], the energy of eight neighboring pixels for each contour point is tested to find the local minimum of the energy. As a result, their model can be easily trapped in local minima in the case of multiple edges. Tsai et al. [18] have used a Hopfield network for minimizing the energy of the contour to provide faster convergence and retain better stability. The aforementioned methods are very slow and thus lack the important advantage of automatic segmentation over manual segmentation (speed). One of the problems associated with discrete contour models is clustering of the vertices. Ranganath [19] has introduced a solution to this problem by defining a constraint for the minimum distance between adjacent points on the contour. This constraint limits the flexibility of the model and it is hard to find a correct threshold. Cohen and Cohen [20] have used a finite element method for minimizing the contour energy. They have suggested a solution for instability and oscillation. Their method replaces image forces by a binary force, which weakens the contour behavior and makes it very sensitive to noise. They have also introduced a balloon model to make the snake inflate until stopped by image forces. Although their method helps to make the model somewhat independent of the initial contour, optimizing the value of the balloon force in the presence of external and internal forces is a hard task, if not impossible. Some of the authors used analytical or parametric models to represent a contour (see Refs [21,22]). Although this makes the model more compact and is useful for

applications such as motion analysis [23], structures with irregular shapes and sharp corners can not be easily represented by these models and may need a large number of terms and a considerable increase in calculations. In contrast, discrete models seem to be more flexible in fitting into edges with high curvature [24,17]. Christensen et al. [25–27] have developed a digital anatomical textbook model in which, the image of the patient’s brain is transformed into digital atlas domain by using a probabilistic transformation that has the physical properties of deformable elastic solids. Different regions of the patient’s brain and its deformation from the normal shape (for normal variations among humans only) will be obtained from the digital atlas accordingly. Leitner and Cinquin [9] have redefined internal forces and used the results for 3D segmentation. Malladi et al. [28] used level sets of a function to represent the contour. Kichenassamy et al. [29] have used gradient flows for the contour evolution. Many researchers have used deformable contours for a variety of applications (e.g. Refs [30–32]). Miller et al. [33–35] have introduced a geometric deformable model. Although their model remedies the problem of finite difference approach in minimizing the contour energy, it does not effectively solve the problem of shrinking and self-cutting. Recently, Lobregt and Viergever [36] have introduced a new discrete dynamic contour model. They have used a curvature based internal force to solve the above problems. However, their model is not well suited for segmentation of objects that have multiple edges near their boundaries. They have introduced a damping force for stabilizing the contour, but it is hard to find proper weighting for this force from image to image or within different areas of an image. Also, this force can not always make the contour stable without decreasing its accuracy. They have not addressed the problem of extracting a proper image feature for defining external forces. Our model follows the geometry introduced in [36], but we have redefined external forces by using an edge-tracking method. This method helps the contour to be absorbed to the correct edges. We have addressed the problems associated with optimizing the internal force weight, contour stability and extraction of image features for external energy calculations. Contributions of this work are six-fold: (1) a new approach in using the external force for deformation of the contour to produce better results near multiple and discontinuous edges; (2) introducing adaptive force weight for overcoming the problem of weight optimization; (3) limiting the contour movement in each iteration using an analytical function to obtain faster convergence and more reliable results; (4) using a continuous estimation of the image energy to generate more accurate and stable results; (5) comparing some standard edge detection operators for deriving external force; and (6) illustrating the performance of the new approach for the segmentation of the hippocampus from brain MRI.

A. Ghanei et al./Computerized Medical Imaging and Graphics 22 (1998) 203–216

3. Methods In this section, the original structure of the model and force fields are briefly described first. Then, improvements in the model are described and changes introduced into the deformation process to solve the problem of multiple edges are discussed. Next, adaptive internal force weights are described and a solution for the contour oscillation and instability is presented. Finally, the role of appropriate functions for extracting desired image features is presented. 3.1. Representation of a discrete contour Here, we will briefly describe the original model and refer the reader to Lobregt and Viergever [36] for a complete description. The model consists of a set of vertices connected to each other by edges (see Fig. 1). The position of vertex i is represented by vector p i and the edge segment between vertices i and i þ 1 by D i. Therefore, we may write: Di ¼ pi þ 1 ¹ pi

(1)

Vertices are indexed in a clockwise direction. The unit tangential vector at a segment is defined along the segment direction and at a vertex, it is considered to be in the direction of the average of unit tangential vectors of two adjacent segments, i.e. ti ¼

di ¹ 1 þ di kdi ¹ 1 þ di k

(2)

where t i is the unit tangential vector at vertex i and d i is the unit vector in the direction of D i and k·k is the Euclidean norm. The radial vector, r i, is obtained by rotating the tangential vector 908 clockwise. The movement of each vertex during the deformation is along its radial vector. In the model, elastic forces are based on the curvature of the contour. Curvature is defined at a vertex as the vector difference between the unit tangential vectors of two joining edge segments: ci ¼ di ¹ di ¹ 1

(3)

This vector is parallel to r i, therefore, c i can be represented by its magnitude, C i, and the vector r i as: c i ¼ Ci ri

(4)

Deformation process starts from an initial contour. Each vertex is moved according to the sum of internal, external, and a damping force. The internal force, f in,i, is computed from the curvature of the contour by convolution of the C i sequence with a discrete filter to prevent changes in parts of the contour with average curvature of zero: f in, i ¼ (Ci # ki )ri

205

to the image data. The damping force, f damp,i, is set proportional to the square of the vertex speed. The contour deformation is computed in discrete positions in time pi (t þ Dt) ¼ pi (t) þ vi (t)Dt

(6)

vi (t þ Dt) ¼ vi (t) þ ai (t)Dt

(7)

ai (t þ Dt) ¼ f i (t þ Dt)=mi

(8)

f i ¼ wex f ex, i þ win f in, i þ wdamp f damp, i

(9)

where a i, v i and m i are vertex acceleration, velocity and mass, respectively, and w ex, w in, w damp are weights for the external, internal, and damping forces, respectively. 3.2. Model improvements 3.2.1. Edge tracking When there are multiple edges near the boundaries of a desired object, such as the hippocampus, the problem is that the contour may get trapped at undesired edges. Even if the initial contour is near the boundary of the object, some effects like internal forces can push parts of the contour inside the gradient force field of an unwanted edge. In existing models, this will cause that part to be absorbed by the wrong edge. For example, if an object with a stronger edge exists near the desired object, its gradient field can easily absorb the contour. The problem also happens when a relatively strong part of the desired object edge pulls the contour from other sides of the object by its strong gradient field. Unfortunately, above cases occur quite often in CT and MRI images. On the other hand, to ensure that different initial contours, drawn at different times or by different operators, will feel the force field of the desired edge, we may use a sufficiently large window for computing the potential energy of image pixels or use multi-scale and blur functions. However, in cases like the hippocampus, these procedures can mix boundaries of the different neighboring structures and move the contour to an undesired boundary. Another property needed in cases like the hippocampus is that the contour should close the boundary in discontinuous parts of the object edge without expansion to the outside. These problems limit application of dynamic contour models. We have developed an edge-tracking method to

(5)

where the # sign denotes the convolution and k ¼ {…,0,0,-1/2,1,-1/2,0,0,…}. The external force, f ex,i, is equal to the inverse of the gradient of the image energy, E im, which is obtained by applying an appropriate operator

Fig. 1. The model consists of a polygon. The position of each vertex is represented by a vector p i and vertices are connected by segments D i. t i and r i are tangential and radial vectors, respectively.

206

A. Ghanei et al./Computerized Medical Imaging and Graphics 22 (1998) 203–216

solve these problems. The procedure is as follows. For each vertex like i, we search along r i to find an edge. A location on r i is interpreted as edge, if the external energy has a minimum at that location. More precisely, each point on r i can be represented by a parametric vector such as s a,i ¼ p i-ar i. We slide sa;i along r i by changing a. a is changed between -a max and a max, using Da as the step size. The value of a max determines the maximum allowed distance between initial polygon and the object edges, and the value of Da determines the maximum possible error from local minimum of energy. We have used Da ¼ 0.5 and a max ¼ 15 (both in pixel unit) in our experiments. In each step, we compare the energy of s a,i to its neighboring locations on r i. This search is continued until we reach a local minimum of the energy, where we interpret s a,i as edge (or when a exceeds a max). We name this point s min,i ¼ p i-a minr i at which: Eim (smin, i ) # Eim (smin, i 6 Dari )

(10)

where E im(x) is the image energy at point x, which will be discussed in Section 3.2.7. The external force pulls the vertex i towards the s min,i, and the amount of force is proportional to the distance of the vertex from s min,i: f ex, i ¼ b((smin, i ¹ pi )·ri )ri ¼ ¹ bamin ri

(11)

where b is a positive constant with dimension of kg s ¹2 and changes the length dimension, m, to the force dimension Kg m s ¹2. Using local minimum of the energy for finding the edge is more reliable than using a threshold, which can be hard to find a proper value for it. Here, we use another threshold which we name E t, for finding an edge. We consider the energy of s a,i for comparison with its neighbors only if it is less than E t. This approach ignores faint minima in the energy of image pixels due to noise. This is entirely different from using the former threshold value. As E t gets larger, it becomes less probable to ignore some edges during the search but we lose some of the advantages of canceling noise effects. The optimal value of E t can be obtained by comparing an area of the image that does not have any edges, with another area that contains edges. Once an optimal value is found for E t, it can be used for the same type of images and applications. In the case of multiple edges, we take advantage of the above algorithm. In this case, the contour will be absorbed by boundaries of the first object. Here, if the vertex reached to the edge, then s min,i will not change further and f ex,i will point to the edge, even if the vertex moves away due to other reasons. It should be noted that the described method requires the initial contour to be drawn closer to the desired object in the presence of multiple edges. However, this is a reasonable criterion, because the intended object must be distinguished somehow from neighboring objects for applying the model.

3.2.2. Balloon force The above edge-tracking method generates the flexibility to introduce a balloon force to the model. This option can be used to enable a small initial contour to be absorbed to the far boundaries of the object. The balloon force acts on the vertices for which no minimum of the energy is found along the search length. This force pushes the vertex outward: f ex, i ¼ ¹ bamax ri

(12)

The minus sign causes the force to be outward. In this method, the contour expands until it is stopped by an edge, so even an initial small contour can detect a large pattern. Matsumoto et al. [24] have designed a balloon model by introducing anchor points outside the contour and defining forces between the contour points and anchors. Their approach, however, uses a predetermined inflating force that is independent of the image and degrades the contour dynamics. It can also cause the contour to ignore weak edges. Other authors like Kucera and Martin [37] and Cohen and Cohen [20] have also introduced balloon models. Cohen have used the balloon model in the form of a constant outward force which inflates the contour, but our model has a basic difference to Cohen’s model. Here, the balloon force is not defined for all contour points equally, but only for those without any edge in vicinity. As a result, the problem of optimizing the amount of balloon force in the face of external forces does not exist in our model. In fact, these two forces will not exist for a vertex simultaneously, and the balloon force disappears as soon as the vertex gets close to an edge. If we remove the minus sign in the previous equation, the contour will shrink instead of inflate and can detect an object that is inside the initial contour. In the above inflating (shrinking) model, we limit our search for the edge to the outside (inside) of the polygon. In our applications, it was possible to satisfy this criterion when drawing the initial contour. It can be seen that, in this method, we do not need to use a large window size for computing the external energy. The problem of far edges from the initial contour can be overcome in this method by setting a large value for the search distance or using a balloon force. Another commonly used solution is a multi-scale approach, e.g. [36], but again the problem of undesired edges in multi-edge places will arise. 3.2.3. Virtual edges (forces) A virtual edge generates a virtual force, which is in contrast to the balloon force and thus they can not be used together. When there is no edge in the vertex vicinity, for example in the discontinuous parts of an object boundary, we can use virtual edges. Virtual edge can be the initial vertex location, or a user defined anchor point. We virtually assign the minimum energy to this location: smin, i ¼ sfix, i

(13)

where s fix,i is the virtual edge. The amount of force toward

A. Ghanei et al./Computerized Medical Imaging and Graphics 22 (1998) 203–216

this edge can be different from real edges. It can vary between zero (maximum freedom for vertex movement) and infinity (vertex fixed at its location). The amount of this force is related to the amount of deformation from the initial contour that is desirable for areas without an edge. Virtual edges are useful for the slices where the hippocampus has discontinuous boundaries. In fact it helps the contour to close the missing parts of the object edges using: internal forces; other parts of the object boundary (which are present); and the user drawn initial contour which can convey the human knowledge about the missing parts of the boundary. In the existing models, the last term (initial shape) makes no contribution, in a controlled manner, in estimating the missing edges, i.e. the user has no control on the final result by intentionally changing the initial contour. The following algorithm briefs the procedures described in Section 3.2.1, Section 3.2.2, Section 3.2.3: 1. 2. 3. 4. 5. 6. 7.

Set i ¼ 0. If i ¼ N, then exit. If virtual-force, set s a,i ¼ s fix,i. Go to step 7. Set a ¼ -a max. Set s a,i ¼ p i-ar i. If E im(s a,i) # E im(sa;i 6 Daðri Þ and Eim ðsa;i Þ # Et then, Edge is found. Set s min,i ¼ s a,i and a min ¼ a. Compute f ex,i from Eq. (11). Set i ¼ i þ 1 and go to step 2. 8. Otherwise, if a , a max, set a ¼ a þ Da and go to step 5. 9. If balloon model, go to step 7. 10. Set f ex,i ¼ 0. Set i ¼ i þ 1 and go to step 2.

3.2.4. Adaptive internal force weight One of the problems with deformable contours is that elastic or internal forces act equally on all contour points. Therefore, when there are large differences in image driven forces (external forces) along the structure boundaries, it causes non-uniform smoothness. This makes it hard to find a globally optimal internal force weight and obtain uniform smoothness and uniform matching of the contour to the object boundaries. If we use a large value for w in, there will be a large difference between the contour and the object edges in sharp corners. If we use a small value, we do not get enough smoothness and alliance of the contour points. To solve this problem, Liao and Medioni [38] have used a variable amount for elastic energy coefficient which is proportional to the ratio of the external energy to the internal energy. Their method cancels the effect of the internal energy and sets the amount of the elastic force proportional to the external force. Therefore, in their model, the local geometric shape of the model does not contribute to the acting forces. We use variable w in, but we do not use direct proportionality. Our adaptive value for w in is not related to inverse of the internal forces. In this manner, w in changes slightly and partially with the changes of f ex. The ratio of the effect of the local f ex to the local f in can be controlled with a parameter like

207

0 # l # 1 as follows: win, i ¼ (1 ¹ l þ l

kf ex, i k )w maxkf ex k 0

(14)

where w 0 is the initial value for w in. The above relation is more appropriate than a direct proportionality between w in and the external force, which cancels the effect of f in completely in places without f ex. At the extreme values of l we have: l ¼ 0 : win will not change and will be independent of f ex : l ¼ 1 : win will be directly proportional to f ex : It must be noted that none of the extreme values are desired, specially the later one, since it may cause the contour to be absorbed to weak edges. Between the extreme values, w in is partially dependent on f ex, but it also has a constant bias that eliminates the above undesired effect. 3.2.5. Limited movement The deformation process described in Section 3.1 may generate undesirable results because of large movements of the contour. For example, a large external force for a vertex generates a large displacement for it, pushing the vertex away from the desired object. We have observed that although it is an effective way to stabilize the contour, damping force of the original model can not totally solve the above problem. Also, it needs a lot of trial and error to find a proper value for w damp. Cohen and Cohen [20] used a binary external force which is either zero or a constant to remedy this problem. Therefore, their model is very noise sensitive. In our model, we overcome the above problem by limiting the maximum movement of the vertex. This can be done by a proper mapping of the vertex velocity to its displacement. Such a mapping should be linear for small values, preserving the linearity of the deformation process. We use an exponential function for this purpose: v :r (15) Dpi ¼ Lmax {1 ¹ exp( ¹ kvi kDt=Lmax )} i i kvi k where L max is a constant and determines the maximum movement of the vertex, Dp i is vertex displacement and v i is its speed, both along the radial vector. This function has an upper limit and is monotonic increasing. It is zero at origin and approximates Dtv i for small values of kv ik. It should be noted that the smaller the L max, the slower is the deformation and the larger the L max, the faster is the deformation, at the expense of less reliability. We have used L max ¼ 1 in our experiments which gives a good compromise between resolution and speed. 3.2.6. Energy estimation The final contour will oscillate if there are two states with minimum energy for the contour, due to dealing with discrete energy functions. In each state, there exists an attraction to the other state, causing the contour to oscillate

208

A. Ghanei et al./Computerized Medical Imaging and Graphics 22 (1998) 203–216

between these states. For solving this problem, we use an estimation procedure for the external energy inside each pixel. Although other authors may have used a similar solution, we have not seen any reference to it in the literature. Therefore, we present a short discussion here. We use an interpolation for the image energy to generate a continuous map. The energy assigned to a location is interpolated from the energy of the pixel at that location, as well as its eight neighboring pixels. It is similar to the bilinear interpolation, but it takes into account the energy of all the pixel neighbors, while retaining continuity at the pixel boundaries. This continuous mapping of energy produces a single state with minimum energy between the two previous minima. In addition to solving the oscillation problem, this approach makes the minimum movement of the vertex less than a pixel and eliminates digitization error. In location (x þ a,y þ b) where x and y are integer values and 0 # a, b , 1 (a ¼ b ¼ 1/2 corresponds to the pixel center), we define: Eim (x þ a, y þ b) ¼ (1 ¹ a)(1 ¹ b)Eim (x ¹ 1, y ¹ 1)

mentioned previously, the hippocampus has a relatively low contrast and low signal-to-noise ratio in MRI, requiring a better suited function for extracting the image features. Here, we explain two suitable edge detection operators for this purpose and compare them with gradient based operators. 3.2.7.1. Step expansion filter (SEF). SEF has been introduced by Raghunath and Ben-Arie [39]. It is an edge detector that uses expansion matching and restoration. It is optimal in the sense of a special figure named discriminative signal-to-noise-ratio [39]. Assuming the noise in the image can be modeled with white noise and the edge model is a step function, the impulse response of the SEF filter is: h(x) ¼ sgn(x)e ¹ 2lxl=j

(17)

where j 2 is the variance of the white noise and sgn(·) is the sign function.

þ 2(1 ¹ a)(1 ¹ lb ¹ 0:5l)Eim (x ¹ 1, y)

3.2.7.2. Canny edge detector (CED). Canny [40] showed that the ideal operator for maximizing the conventional signal-to-noise ratio in detecting a particular edge is the correlation with the edge model itself. Canny’s edge detector for step edges is well approximated by the derivative of the Gaussian mask. It can be written as:

þ M(1¹la ¹ 0:5l)(1 ¹ lb ¹ 0:5l)Eim (x, y)

h(x) ¼ xe ¹ x

þ 2a(1 ¹ lb ¹ 0:5l)E(x þ 1, y)

where d is the variance of the detector. The values of parameters d and j are dependent on the imaging system and procedures and parameters used. However, usually an acceptable estimation for them can be obtained by analyzing uniform regions of the image. These values can be estimated once for each kind of application and then be used in similar cases.

þ 2(1 ¹ la ¹ 0:5l)(1 ¹ b)Eim (x, y ¹ 1) þ a(1 ¹ b)Eim (x þ 1, y ¹ 1)

þ (1 ¹ a)bEim (x ¹ 1, y þ 1) þ 2(1 ¹ la ¹ 0:5l)bEim (x, y þ 1) þ abEim (x þ 1, y þ 1)

ð16Þ

The coefficients on the left hand side of each energy term represent the vertical and horizontal distances of the point from the center of the corresponding pixel. The value of M determines the contribution of the energy of the pixel that contains the point — pixel at (x,y) — to the total energy. It can easily be seen that with M ¼ 4 this function is continuous over the image (see Fig. 2). 3.2.7. Energy functions Using the gradient of the image gray levels as the external energy may be suitable for some applications, but as

Fig. 2. Continuity of energy estimation. The energy of a point like (x9,y9) on a pixel boundary can be interpolated as it is inside the pixel (x,y) or inside the pixel (x þ 1,y); both give the same result.

2

=2d

(18)

3.2.7.3. Other operators (gradient based). Although not as good as SEF or CED, gray level gradient in multiple directions can be used. It is a good operator in the sense of zero localization error, since it is a linear operator [41]. Also, it is relatively simple and fast. We have tested operators like the absolute difference of the pixel with its neighbors. In most cases, we obtained the best result when using the SEF function. Any of these functions can be applied along the two orthogonal directions — like (0,1) and (1,0) — and the sum of the absolute values of the results can be considered as an energy field for the edge feature over the image. This is a cityblock metric, which gives a satisfactory approximation to the Euclidean metric but it is much faster [42]. There are a number of noise filters that can be applied to the image (e.g. Refs [43–46]), before the edge operator is applied, to improve the performance of the edge detector. However, it should be noted that the optimal selection of filter parameters is crucial in obtaining smooth results while preserving low-contrast edges.

A. Ghanei et al./Computerized Medical Imaging and Graphics 22 (1998) 203–216

209

Fig. 3. Using the model on a simulation to show different aspects of edge-tracking. (a) Initial polygon; (b) final contour when free vertices (those not near an edge) are fixed; (c) final contour when free vertices are not fixed; and (d) result when balloon force is used.

4. Results We have implemented the described improvement strategies in the model and applied the result to MRI. Although our main aim was segmentation of the hippocampus from MRI, the application of the method is not limited to this case. The method is fast; it takes a few seconds to run on a SUN Sparcstation 20. In the remainder of this section, we describe the parameters we have used and illustrate the results of our experiments. Because of the numerous parameters, the degree of freedom in our model may seem too high, but most of the discussed parameters can be fixed for a wide range of applications, once they have been set by the designer. In our experiments, we set m i ¼ 256 and w damp ¼ -1. The value of a max is set to 15 from observations of several typical initial polygons drawn by different operators, and their maximum distance to the object boundaries. Da is set to 0.5 for maximum accuracy. A smaller value does practically nothing except slow down the model and a larger value reduces the resolution. We set b ¼ 1/a max which normalizes the length of f ex. These parameters usually do not need to be changed by the user. Other parameters which are typically changed by the user are: w ex, w in and l. Nevertheless, once a correct set of parameters has been found for an application, it can be used for all similar cases. Sometimes, adjustment may be required for some of the parameters. We set w ex ¼ 1, w in ¼ 2.5 and l ¼ 0.2 for all cases except for a few ones which will be noted. An option for the user is to select the type of external force used in the model. The proper selection should be clear by the descriptions given in Section 3.2.1– Section 3.2.3. For example, for an object with discontinuous boundaries, edge-tracking with virtual-force may be used, and

for a closed object, the balloon-force is more appropriate to save time in defining the initial contour. The deformation process is stopped automatically when the contour movements become small, or when the number of iterations exceeds a limit. Lobregt and Viergever [36] have defined the convergence when v i ; a i ; f i ; 0 for all i. We have found that this condition is usually very hard to reach. Therefore, we define the convergence criterion to be the maximum displacement of the contour points less than a threshold in three consecutive iterations. This threshold is experimentally set to be 20% of the maximum displacement of the first few iterations (in which the contour has the most movement because it is farthest from the object boundaries). In this method, the threshold value is independent of the image type and object boundaries. We integrated the this method into a user friendly software package which allows the user to load all the slices of a volume data and use the initial or final polygon of each slice as the initial polygon for the next slice(s). This polygon can be modified before it is used for the new slice. This can dramatically reduce the interaction required, because of the similarities between the adjacent intersections in volumetric data. The above flexibilities increase the advantages of the model over manual segmentation in regard to speed, easiness, and reproducibility. In segmenting six sets of hippocampus images, the total time spent for segmentation of each slice, including user interactions, was between 5 and 30 s (14 s on average). This is in contrast to the manual segmentation, which was between 40 s and 2 min (67 s on average). In the following, we illustrate the results of applying the model to several practical situations. Fig. 3 shows the results for a simulated image. The effect of fixing parts of the contour that is not near any edge

Fig. 4. A CT image of a phantom for which the balloon model is used to extract the object. (a) Initial polygon; (b) intermediate result of (a); (c) final result of (a) with the balloon force; and (d) final result on (SEF) energy distribution.

210

A. Ghanei et al./Computerized Medical Imaging and Graphics 22 (1998) 203–216

Fig. 5. An MRI image of the egg phantom used for the data presented in Table 1. Initial polygon: (a) on the image; (b) on the energy distribution (SEF). Final contour: (c) on the image; (d) on the energy distribution.

(virtual-edge) is shown in Fig. 3(b) and can be compared with Fig. 3(c) in which those vertices are not fixed. In Fig. 3(d), using the balloon force caused the contour to be absorbed to the far edges of the object. Fig. 4 shows the result for a CT image of an egg phantom. The balloon force is used to detect the closed boundary of the object, which is far from the initial polygon. The final contour on the energy distribution illustrates the accuracy of the model. Fig. 5 shows the result for an MRI image of another egg phantom. Both the initial and final contours have been shown on the image and on the energy distribution. This example has been used in testing the accuracy of the model, which will be described later. Examples of the hippocampus segmentation are shown in Figs 6 and 7. Fig. 6 shows the result of the algorithm for a raw (not filtered) image of the hippocampus. The noise power is relatively high and edges are weak. The initial polygon has four vertices. The number of vertices increases during deformation by resampling. In this example, we have used a SEF operator with j ¼ 5 to obtain the external energy which is shown in Fig. 6(b). Intermediate and final results are shown in Fig. 6(c) and (d), respectively. Fig. 7 shows

another slice of the hippocampus. In this slice multiple undesired edges exist near the hippocampus. Therefore, the window size for calculating the external energy is reduced to 3 3 3 to avoid undesired edges. This makes the estimated external energy noisier than the previous case. We use a nonlinear edge-preserving filter, specifically designed for MRI by Soltanian-Zadeh et al. [46], to reduce the noise level. The energy distribution is smoother than the previous example, but edges are very unclear. In both cases, the results have been acceptable to the collaborating neuroradiologist. Fig. 8 shows the significance of using the edge-tracking algorithm in the proposed method. In the original method (without edge-tracking), the contour gets trapped in undesired edges because there are several edges near the object boundaries, but in our model this effect has been overcome. Figs 9 and 10 are examples of the balloon approach. In Fig. 9 balloon force is used to let a small initial contour detect the object boundaries. The result is compared with that of drawing the initial polygon close to the object boundaries. Since the boundary is closed, the method is successful. As mentioned before, our balloon model is not sensitive to the amount of balloon force. Fig. 10 shows another example

Fig. 6. (a) An image of the hippocampus with an initial contour defined by four vertices. Noise power is high. (b) External energy distribution and the initial contour; (c) intermediate result of contour deformation; and (d) final result.

Fig. 7. (a) Image of a slice of hippocampus with the initial polygon, image is noise filtered but multiple edges are present; (b) initial polygon on energy distribution; (c) intermediate result; and (d) final result.

A. Ghanei et al./Computerized Medical Imaging and Graphics 22 (1998) 203–216

211

Fig. 8. Effect of using the edge-tracking in the model for the hippocampus segmentation: (a) initial contour; (b) initial polygon on energy distribution; (c) result without edge-tracking; and (d) final result of the method (with edge-tracking).

of using balloon force. For this example, we selected a slice of the hippocampus with discontinuous boundaries to illustrate the aforementioned limitation of using the balloon model (see Section 3.2.2). In this case, the contour expands outside of the object. In Fig. 11, the results of using different weights for the internal force are compared. With larger w in, we get a smoother result and somewhat compensate the effects of the noise on the detected boundary. This is achieved at the price of loosing some details [compare Fig. 11(c) and (d)]. Fig. 12 shows the advantage of using adaptive w in. An image of the hippocampus is shown and a deformable contour is used to segment the object and separate the extra part in the right side. We need to use a large w in for this purpose, but this significantly reduces the accuracy in other sides [Fig. 12(b)]. For increasing the accuracy and obtaining more precise edges, w in has been decreased, but the powerful edges in the right hand side (top and bottom) have absorbed the contour [Fig. 12(c)]. Using the adaptive value for w in, we can separate the undesired area from the

object and obtain the desired accuracy in other parts of the boundary [Fig. 12(d)]. Note that the undesired part has also been excluded in the initial polygon [Fig. 12(a)]. In Fig. 13, it is shown that estimation of energy and limitation of movement can prevent the model from jumping over desired edges and being absorbed by wrong edges. In this image, several structures are present near each other. This situation also causes instability of the contour. Applying proper estimation of the energy and limiting the contour movements have remedied the problem. In Fig. 14, different edge detector operators (external energy functions) are compared. With absolute difference, details are best preserved and edge localization error is not present. However, this operator is sensitive to noise and is not suitable for low contrast boundaries of the object. Similarly, the gradient operator is good in the sense of simplicity and edge localization, but it destroys details, especially in places with multiple edges like the hippocampus. SEF operator is better than the previous two operators in the sense of low sensitivity to noise and preserving details,

Fig. 9. Using the balloon force for detecting closed boundaries of an object that is far from the initial polygon: (a) initial polygon, which is drawn very small; (b) the result of using (a); (c) a typical initial polygon for the same case; and (d) the result of using (c). Comparison with (b) shows good reproducibility of the method.

Fig. 10. Effect of the balloon force for detecting discontinous boundaries: (a) Initial contour; (b) result using the balloon force, note that the contour has expanded outside the structure; (c) result without the balloon force.

212

A. Ghanei et al./Computerized Medical Imaging and Graphics 22 (1998) 203–216

Fig. 11. Effect of internal force weight on the result: (a) initial contour; (b) initial polygon in energy distribution; (c) final result of the method with w in ¼ 2; and (d) final result with w in ¼ 4.

given proper value for j is used. With a large j, this operator acts like the standard gradient operator over its window. For CED with a large value for d, we lose details but get more noise resistance. In Fig. 15, the result of the proposed method is compared with the manual contours drawn by a neuroradiologist on two hippocampus images from a brain MRI study. The hippocampus boundaries are shown using white dots for the manual results and using gray contours for the computer result. This Figure presents the best and worst cases in the sense of matching with the manual results. Minimum and maximum relative difference for the areas determined by the two methods have been 2.07% and 8.7%, respectively. The false positive error (i.e. pixels that have been segmented by the model but are not in the manual segmentation) was 4.57% and 12% for these two cases, while the false negative error (i.e. pixels that have not been segmented by the model, but labeled as hippocampus in the manual segmentation) was 2.5% and 3.3%, respectively. To show the versatility of the proposed method, we have applied it to two clinical CT images. Fig. 16 shows the results of the model, segmenting rectum from an image

with a large noise power and discontinuous object boundary. We have used the SEF edge detector with j ¼ 11 in this case. In Fig. 17 the algorithm is used to segment liver from a CT image. The boundaries are very unclear in some areas. In this case, the tissue in the lower part of the liver, which is part of the gallbladder and does not have any detectable boundary, can be excluded by increasing w in. Table 1 shows the results obtained from the application of the method to volumetry of an egg phantom. Volumes of egg white and egg yolk have also been measured physically, using water displacement (within specified accuracy). When using the deformable model, in order to prevent losing parts of the objects, we included partial volume regions. As such, the results of the model over-estimated the volumes, but are still within an acceptable range (,5%). This example shows an initial estimate for the accuracy of the model. Note that the slice thickness of 7.5 mm has caused a major part of the relative error. It has generated a relatively large partial volume region around the boundaries. This partial volume is shown for one of the slices in Fig. 5. Table 2 compares the result of the manual and automatic segmentation of the hippocampus from a set of MR images.

Fig. 12. Effect of using adaptive w in: (a) initial polygon; (b) result of (a) using a large win to exclude the right area from the contour (w in ¼ 3); (c) result of (a) using a smaller value for w in to obtain better boundary matching (w in ¼ 1); and (d) result using adaptive w in with l ¼ 0.5 and w 0 ¼ 3.

Fig. 13. Effect of estimation of energy and limitation of contour movements on the stability of the result: (a) initial contour; (b) result of (a) after a few iterations; (c) result of (a) after more iteration; and (d) final result of (a) using interpolation of energy and limitation of movement.

A. Ghanei et al./Computerized Medical Imaging and Graphics 22 (1998) 203–216

213

Fig. 14. Comparison between different edge operators described in Section 3.2.7: (a) maximum of absolute difference with neighboring pixels; (b) gradient amplitude of gray levels (window size ¼ 7 3 7 pixels); (c) SEF edge operator with j ¼ 0.8; (d) SEF with j ¼ 0.4; (e) CEF with d ¼ 0.7; and (f) CEF with d ¼ 2.5.

The hippocampus slices are segmented by a neuroradiologist. A 1 and A 2 are the area (in pixels) obtained by manual and automatic methods, respectively. The fourth column is the area covered by both of the methods. The relative amount of agreement and disagreement between the model and manual segmentation are shown by R c and R n, respectively, and for each slice they are computed as: Rc ¼ 2

A1 ∩ A 2 A1 þ A 2

(19)

Rn ¼ 2

A1 ∪ A2 ¹ A 1 ∩ A2 A1 þ A 2

(20)

where R c and R n are the ratio of the common pixels (pixels that marked as the hippocampus by both of methods), and non-common pixels (pixels that are marked by one method only) to the average area, respectively. It can be seen that there exists a good agreement between the results obtained by the two methods. Initial results for evaluating reproducibility of the model are shown in Table 3. A set of the hippocampus images have been segmented twice. The average of the relative differences between the results is within an acceptable range (less than 5%). In all of the above cases, we have used two-pixel resolution (maximum length of segments are two pixels). The execution time has been less than 3 s per slice, using a SUN Sparcstation 20.

5. Concluding remarks Solutions to some of the problems associated with deformable models in their application to medical image segmentation have been found and implemented. We developed a method, which can be used for the segmentation of the hippocampus from brain MRI for the purpose of detecting epilepsy. The method is applicable to other structures with low contrast and multiple and discontinuous edges. As examples, its use for segmentation of rectum and liver from CT images were shown. Initial tests on the program by two operators suggest it is reproducible, accurate, easy to use and useful. Although we have used adaptive weights, the model is still dependent on parameter values and the initial polygon. Future work can be directed towards making the model more independent of the parameters and the initial polygon. It can also be directed towards the extension of the model to 3D. Geometrical structure of the model has a good potential for extension to 3D. We are working on this subject and initial results are promising. Extension of the method to incorporate information from multiple images of the same location, normally acquired in MRI studies, will be a third direction for the future work.

Table 1 Accuracy of the model. Volume estimation in an egg phantom and comparison of the computer results and physical measurement Volume of Egg yolk Egg white Egg

Computer result Measured value (cm 3) (cm 3) 23.34 35.04 58.38

22.14 6 0.49 34.76 6 0.46 56.90 6 0.16

Relative error (%) 5.27 0.8 2.53

Fig. 15. Comparison between results of the algorithm and results of manual detection of boundaries by a neuroradiologist: (a) image with minimum area difference; and (b) image with maximum area difference between manual and computer results.

214

A. Ghanei et al./Computerized Medical Imaging and Graphics 22 (1998) 203–216

Fig. 16. Application of the method for segmentation of rectum from a CT image. Image is very noisy and object boundary is unclear and is discontinuous in its top portion. (a) Initial contour; (b) energy distribution with the initial contour; (c) intermediate result; and (d) final result.

Fig. 17. Application of the method for segmentation of liver from a CT image: (a) initial contour; (b) intermediate result; (c) final result; and (d) final result with a smaller w in.

6. Summary In this work, we have investigated the problem of segmenting the hippocampus from brain MRI. The segmentation of the hippocampus and measurement of its volume from brain MRI is a crucial step in diagnosis of some important mental diseases such as epilepsy. We have used an improved discrete deformable model for this purpose. We have introduced a new external force, which is based on searching for local minima of the image energy in the contour normal direction. This method solves the problem of contour deformation near multiple edges. With

the help of this definition, we have introduced a new balloon force, which is insensitive to the forces weights. It helps the contour to be absorbed to the far edges of the object, without needing multi-scale or blurring techniques, which can lead the contour to undesired edges. Reliability of the final contour and the model stability have been improved by using a continuous mapping of external energy and limiting movements of the contour. Also, the problem of optimizing the internal force weight has been overcome by making it dependent on the amount of external and internal force. Finally, capabilities and limitations of the model have been tested and illustrated using real clinical MRI studies, phantom studies, and CT images.

Table 2 Comparison between the result of manual segmnentation and using the purposed model on a set of hippocampus slices Slice

A1

A2

A1 ∩ A2

1 2 3 4 5 6 7 8 9 10 11 Total volume

306 251 174 180 195 198 268 443 499 353 306 3173

311 211 185 204 192 211 261 435 441 390 311 3152

315 251 201 206 201 218 268 475 508 399 325 3367

Rc

Rn 0.98 0.91 0.88 0.93 0.96. 0.93 0.99 0.92 0.92 0.93 0.95 0.94

0.04 0.18 0.24 0.14 0.08 0.14 0.02 0.16 0.16 0.14 0.1 0.12

A 1 and A 2 are the manual and model results, respectively. R c and R n show the relative agreement and disagreement between the two methods. The areas are in pixel unit.

Table 3 Reproducibility of the model. The results of the hippocampus segmentation by the model in two different trials (using different initial contours and parameters) are compared Slice location

Try no. 2

Try no. 1

12.2 14.9 17.5 22.8 25.5

216 185 400 394 643

203 174 406 401 611

Average error Average absolute error Standard deviation of error

Relative difference 6.2% 6.1% ¹1.5% ¹1.8% 5.1%

2.8% 4.1% 4.1%

A. Ghanei et al./Computerized Medical Imaging and Graphics 22 (1998) 203–216

Acknowledgements We would like to thank the following people for their help in preparing this work: Jeff Hasenau for his help in preparing data and images and testing the method, Dr William Sanders and Dr Kost Elisavich for manual segmentation of the hippocampus and sharing clinical information with us, and Lucie Bower for development of a user friendly interface for the method.

References [1] Cook MJ. Mesial temporal sclerosis and volumetric investigations. Acta Neurolog. Scand. Suppl., 1994;152:109–114. [2] Heinz R, Ferris N, Lee EK. et al. MR and positron emission tomography in the diagnosis of surgically correctable temporal lobe epilepsy. Am. J. Neuroradiol., 1994;15:1341–1348. [3] Jackson D, Cook MJ. Volume and T2 quantitation boost diagnostic accuracy. MR: News Mag. Magn. Reson. 32–45; Mar. 1993. [4] Haller J, Botteron K, Brunsden B. et al. Hippocampal MR volumetry . SPIE, 1994;2359:660–671. [5] Bogarts B, Lieberman JA. et al. Hippocampus-amygdala volumes and psychopathology in chronic schizophrenia. Biol. Psychiatr., 1993;33:236–246. [6] Clifford R, Jack M, Bentley D, Colleen K, Twomey A, Zinsmeister R. MR Imaging-based volume measurments of the hippocampal formation and anterior temporal lobe: validation studies. Radiology, 1990;176:205–209. [7] Terzopoulos D, Platt J, Barr A, Fleischer K. Elastically deformable models. Comput. Graph., 1987;21:205–214. [8] Grzeszczuk P, Levin DN. Brownian strings: segmenting images with stochastically deformable contours. SPIE, 1994;2359:72–89. [9] Leitner F, Cinquin P. From splines and snakes to snake splines. Workshop: Geometric Reasoning for Preception and Action 264–281; 1993. [10] Terzopoulos D, Vasilescu M. Adaptive surface reconstruction. SPIE, 1990;1383:257–264. [11] Kass M, Witkin A. Terzopoulos, D. Snakes: active contour models. Int. J. Comput. Vis. 321–331; 1988. [12] Hill A, Cootes TF, Taylor CJ, Lindley K. Medical image interpretation: a generic approach using deformable templates. Med. Informatics, 1994;19:47–59. [13] Davatzikos C, Bryan RN. Using a deformable surface to obtain a shape represntation of the cortex. IEEE Trans. Med. Imag., 1996;15(6):785–795. [14] Geiger D, Gupta A, Costa LA, Vlontzos J. Dynamic programming for detecting, tracking and matching deformable contours. IEEE Trans. Patt. Anal. Mach. Intell., 1995;17:294–302. [15] Amini AA, Weymouth TE, Anderson DJ. A parallel algorithm for determining two-dimensional object positions using incomplete information about their boundaries. Patt. Recog., 1989;22:21–28. [16] Poon S, Braun M, Fahrig R, Ginige A, Dorrel A. Segmentation of medical images using an active contour model incorporating regionbased image features. SPIE, 1993;2359:90–97. [17] Sebbahi A, Herment A, de Cesare A, Mousseaux A. Multimodality cardiovascular image segmentation using a deformable contour model. Comput. Med. Imag. Graph., 1997;21(2):79–89. [18] Tsai T, Sun YN, Chung PC. Minimising the energy of active contour model using a Hopfield network. IEEE Proc-E, 1993;140:297–303. [19] Ranganath, S. Analysis of the effects of snake parameters on contour extraction. Proc. Int. Conf. Automation, Robotics Comput. Vis. 4.5.1– 4.5.6; Sept. 1992.

215

[20] Cohen D, Cohen I. A finite element method applied to new active contour models and 3D reconstruction from cross sections. Proc. 3rd Int. IEEE Conf. Comput. Vis. 1990. [21] DeCarlo D, Metaxas D. Blended deformable models. IEEE Trans. PAMI, 1996;18(4):443–447. [22] Szekely G, Kelemen A. Christian Brechbuhler and Guido Gerig segmentation of 2-D and 3-D objects from MRI volume data using constrained elastic deformation of flexible Fourier contour and surface models. Med. Imag. Anal., 1996;1(1):19–34. [23] Park J, Metaxas D, Axel L. Analysis of left ventricular wall motion based on volumetric deformable models and MRI-SPAMM. Med. Imag. Anal., 1996;1(1):53–71. [24] Matsumoto S, Asato R, Okada T, Konishi J. Intracranial contour extraction with active contour models. JMIR, 1997;7:353–360. [25] Christensen GE, Joshi SC, Miller MI. Volumetric transformation of brain anatomy. IEEE Trans. Med. Imag., 1997;16(6):864–877. [26] Christensen GE, Rabbitt RD, Miller MI. 3D brain mapping using a deformable neuroanatomy. Phys. Med. Biol., 1994;39:609–618. [27] Miller MI, Christensen GE, Amit Y, Grenander U. Mathematical textbook of deformable neuroanatomies. Proc. Natl Acad. Sci. USA, 1993;90(24):11944–11948. [28] Malladi R, Sethian JA, Vemuri BC. Shape modeling with front propogation. IEEE Trans. Patt. Anal. Mach. Intell., 1995;17:158–177. [29] Kichenassamy S, Kumar A, Olver P, Tannenbaum A, Yezzi A. Gradient flows and geometric active contours. Proc. Int. Conf. Comput. Vis. ICCV ’95 810–815; 1995. [30] Ranganath S. Contour extraction from cardiac MRI studies using snakes. IEEE Trans. Med. Imag., 1995;14:328–338. [31] Sandor R, Leahy R. Automatic recognition of brain regions from magnetic resonance images. IEEE: Asilomar Conf. Signal, Syst. Comput, 172–177, 1992. [32] Davatzikos A, Prince JL. An active contour model for mapping the cortex. IEEE Trans. Med. Imag., 1995;14:65–81. [33] Miller V, Breen DE, Wozny MJ. Extracting geometric models through constraint minimization. Rensselear Polytechnic Institue Technical Report, No. 90024, 1990. [34] Miller V. On GDM’s: geometrically deformed models for the extraction of closed shapes from volume data. Rensselaer Polytechnic Institue, Technical Report, No. 90051, 1990. [35] Miller V, Breen DE, Lorensen WE, O’Bara RM, Wozny MJ. Geometrically deformed models: a method to extract closed geometric models from volume data. Comput. Graph., 1991;25: 217–226. [36] Lobregt S, Viergever MA. A discrete dynamic contour model. IEEE Trans. Med. Imag., 1995;14:12–24. [37] Kucera D, Martin RW. Segmentation of sequences of echocardiographic images using a simplified 3D active contour model with region-based external forces. Comput. Med. Imag. Graph., 1997;21(1):1–21. [38] Liao C, Medioni G. Surface approximation of a cloud of 3d points. Graph Models Imag. Proc., 1995;57:67–74. [39] Raghunath Rao K, Ben-Arie J. Optimal edge detection using expansion matching and restoration. IEEE Trans. Patt. Anal. Mach. Intell., 1994;16:1169–1182. [40] Canny J. A computational approach to edge detection. IEEE Trans. Patt. Anal. Mach. Intell., 1986;8:679–697. [41] Koplowitz GV, Greco V. On the edge location error for local maximum and zero-crossing edge detectors. IEEE Tran. Patt. Anal. Med. Imag., 1994;16:1207–1212. [42] Gonzalez RC. Digital image processing. Reading, MA: Addison-Wesley, 1992. [43] Shen J, Castsan S.. Towards the unification of band-limited derivative operators for edge detection. Signal Proc., 1993;31:103–109. [44] Sundaram R, Ersoy OK, Hansen D. Adaptive approach to edge detection. Optical Engineering, 1995;34:3271–3276. [45] van der Heijden F. Edge and line feature extraction based on covariance models. IEEE Trans. Patt. Anal. Mach. Intell., 1995;17:16–33.

216

A. Ghanei et al./Computerized Medical Imaging and Graphics 22 (1998) 203–216

[46] Soltanian-Zadeh H, Windham JP, Yagle AE. A multidimensional nonlinear edge-preserving filter for magnetic resonance image restoration. IEEE Trans. Imag. Proc., 1995;4:147–161. Amir Ghanei was born in Tehran, Iran in 1970. He received his B.Sc. degree in Electrical Engineering from Sharif University of Technology, Tehran, Iran, in 1994, and his M.Sc. degree in the same field from University of Tehran, Tehran, Iran, in 1996. From 1991 to 1993 he was a researcher at the Iran Telecommunication Research Center, Tehran, Iran. He was with the department of Industrial Automation of Kerman Tablo Co., Tehran, Iran, from 1993 to 1995, where he was working on industrial computers and software development for computer automation. Since 1995 he has been with Department of Diagnostic Radiology of Henry Ford Hospital, Detroit, MI, where he is currently researching in medical image processing. He is currently pursuing a Ph.D. in the field of Electrical Engineering: Signal Processing at the University of Michigan, Ann Arbor, MI. His research interests include computer automation, computer vision, medical imaging, and image processing.

Hamid Soltanian-Zadeh was born in Yazd, Iran in 1960. He received B.Sc. and M.Sc. degrees in Electrical Engineering: Electronics from the University of Tehran, Tehran, Iran, in 1986, and M.S.E. and Ph.D. degrees in Electrical Engineering: Systems from the University of Michigan, Ann Arbor, MI, in 1991 and 1992, respectively. From 1985 to 1986 he was a design engineer at Iran Telecommunication Research Center, Tehran, Iran. In 1987 he was a lecturer of Electrical Engineering at the University of Tehran. Since September 1988 he has been with the Department of Diagnostic Radiology at Henry Ford Hospital, Detroit, MI, where he is currently a Senior Staff Scientist. Since Feb. 1994, he has also been with the Department of Electrical and Computer Engineering of the University of Tehran, where he is currently an Associate Professor of Electrical and Biomedical Engineering. In addition, since Aug. 1997, he has been an Assistant Professor of Radiology at Case Western Reserve University, Cleveland, OH. His research interests include medical imaging, image reconstruction and processing, pattern recognition, and neural networks.

Joe P. Windham received his Ph.D. degree in Nuclear Science and Engineering from the University of Cincinnati, OH in 1972. He was at the Medical College of Ohio in Toledo, OH, from 1972 to 1984 where he was an Associate Professor and Head of the Division of Medical Physics of the Department of Radiology. He has been with the Department of Radiology of Henry Ford Hospital, Detroit, MI, since June 1984 where he is currently Head of the Division of Radiological Physics. He is certified by the American Board of Radiology in Radiological Physics. He is active in the American College of Radiology where he is a member of the commission on Physics and Radiation Safety, commission on Education, and commission on Radiation Oncology. He is chairman of the committee on Education under the commission on Physics and Radiation Safety. His research interests include image processing and quantitative analysis of medical images (particularly images obtained from magnetic resonance), computed tomography, and nuclear medicine.