Generalized ellipsoids and anisotropic filtering for segmentation improvement in 3D medical imaging

Generalized Ellipsoids and Anisotropic filtering for Segmentation Improvement in 3D Medical Imaging1 R. Dosil, X. M. Pardo Dept. Electrónica e Computación Universidade de Santiago de Compostela, Spain. [email protected], [email protected] Deformable models have demonstrated to be very useful techniques for image segmentation. However, they present several weak points. Two of the main problems with deformable models are the following: (1) results are often dependent on the initial model location, and (2) the generation of image potentials is very sensitive to noise. Modeling and preprocessing methods presented in this paper contribute to solve these problems. We propose an initialization tool to obtain a good approximation to global shape and location of a given object into a 3D image. We also introduce a novel technique for corner preserving anisotropic diffusion filtering to improve contrast and corner measures. This is useful for both guiding initialization (global shape) and subsequent deformation for fine tuning (local shape). Keywords: registration, deformable models, segmentation, anisotropic diffusion, surface patch saliency, 3D medical images. Corresponding author: Raquel Dosil Lago Dept. Electrónica e Computación Universidade de Santiago de Compostela Monte da Condesa, Campus Sur 15782, Santiago de Compostela, SPAIN e-mail: [email protected] Fax: +34 981 599412 1 This work was supported by Spanish Government and Xunta de Galicia by projects TIC2000-0399-C02-02 and PGIDT99PXI20606B respectively Generalized Ellipsoids and Anisotropic filtering for Segmentation Improvement in 3D Medical Imaging1 R. Dosil, X. M. Pardo Dept. Electrónica e Computación Universidade de Santiago de Compostela, Spain. [email protected], [email protected] Abstract. Deformable models have demonstrated to be very useful techniques for image segmentation. However, they present several weak points. Two of the main problems with deformable models are the following: (1) results are often dependent on the initial model location, and (2) the generation of image potentials is very sensitive to noise. Modeling and preprocessing methods presented in this paper contribute to solve these problems. We propose an initialization tool to obtain a good approximation to global shape and location of a given object into a 3D image. We also introduce a novel technique for corner preserving anisotropic diffusion filtering to improve contrast and corner measures. This is useful for both guiding initialization (global shape) and subsequent deformation for fine tuning (local shape). Keywords: registration, deformable models, segmentation, anisotropic diffusion, surface patch saliency, 3D medical images. 1. Introduction A deformable model [13] is an energy minimization method, where the energy functional is defined in terms of intrinsic shape attributes (internal energy) and desired 1 This work was supported by Spanish Government and Xunta de Galicia by projects TIC2000-0399-C02-02 and PGIDT99PXI20606B respectively image features (external potential), such as gradient and curvature. The external potential originates forces that attract the model to specific image features while the internal energy causes stress forces that try to maintain model continuity and smoothness. When forces are balanced, the model reaches equilibrium, and the geometric deformation finishes. Therefore, the model deforms itself from its initial location to approach the nearest energy minimum, which maybe does not correspond to the target object surface. Deformable matching can efficiently deal with small and local shape changes, but fails if global misalignment is too large. It is very important that the starting model location and the object boundary are near enough, so a good initialization method should be valuable. When deformable models are applied to 3D medical data, two kinds of segmentation are possible: 2 ½ D (slice by slice) and 3D. Initialization is simpler for 2 ½ D segmentation, since it is applied just to the first slice. Each slice is initialized with the result of the previous one. Manual edition of initial 3D surfaces is very laborious, and the automatic or semiautomatic initialization is usually more complex than the 2D counterpart. As a counterweight to the low robustness due to the usual low accuracy in initialization, deformable surface models have the power to ensure smoothness and coherence in 3D shapes. In deformable model literature, we can find different approaches to cope with the initialization problem in 3D. Some approaches are based on the incorporation of balloon forces to overcome potential minima. McInerney and Terzopoulos [17], for example, used a balloon model to initialize the ventricle tracking in MRI data. The main difficulty with this scheme is concerned with the inherent trade-off in the choice of the inflation/deflation force. Other approaches are based on the manual or semiautomatic selection of anchor points. Among them is the imposing of interactive constraints in the form of springs and volcanos [13], and more recently the method of Neuenschwander et al. [21], that allows the user to fix a set of seed points and their normal vectors, which cannot be changed during deformation. Some authors propose multiple initialization through several seed models. During the deformation process, several initial models will merge and the superfluous ones should be removed. This solution was used, among others, by Tek and Kimia [31], and Leonardis et al. [15]. Important decisions have to do with: (1) the number and location of initial seed models, (2) the stopping criteria of the seed growing process, and (3) choosing one result among the final models. Fully automatic initialization can be achieved by matching the object in the image with a prototype, as done by Bajcsy and Kovacic [1] who used brain atlases in their specificpurpose initialization techniques. There are several problems with the deformable atlas approach. The technique is sensitive to initial positioning of the atlas and the presence of neighboring features may cause matching problems. One solution is to use image preprocessing in conjunction with the deformable atlas; Sandor and Leahy [26] used this approach. Several researchers, as Cootes et al. [6], and Staib et al. [30], augmented snake-like models with prior information about typical mean shapes and normal variations. A number or researchers have incorporated knowledge of object shape using deformable shape templates. Among them, superquadrics have gained popularity in medical image research [12, 18, 32]. The segmentation of human organs from CT or MR images is a good example where model-based reconstruction can be applied. The model-based reconstruction problem could be stated in two (no necessarily disjoint) phases [7]: registration and free form deformation. On the one hand, registration describes a transformation with far less degrees of freedom than the free form deformations. Therefore, their ability to represent shape variations is less important than the free form deformation. On the other hand, because of their restricted degrees of freedom, they tend to be more robust than free form deformations. Therefore, the first is better for describing global shape and location and the second is better in detecting fine details. In this work, we propose a fully automatic initialization method (registration) based on matching 3D data with models that comprise global shape and high level feature information. The main idea is to perform the initialization phase not taking into account the relevance of individual point features, but the properties and global saliency of connected point features (surface patches in 3D). The method includes the construction of a priori models from test images. Both shape model construction from test images and matching with a new image are achieved by means of fitting a parametric surface to a cloud of points belonging to the object boundary. This parametric surface is represented by a generalized ellipsoid scheme. Boundary points are extracted from the image, grouped in patches and finally selected according to high level properties. In this way, we try to extend to 3D the ideas developed in a previous paper [22] where a priori knowledge on contour segments was used in order to obtain an improved initialization of 2D deformable models. Extraction of boundary points is not a critical step because the relevance of the point feature depends on all the neighbors in the same surface patch, and matching with global shape models is very robust. However, accuracy in boundary detection has certain influence in the size of detected surface patches and contributes to eliminate noise influence. Noise elimination and boundary detection can be achieved simultaneously using the derivative of Gaussian filter, but this presents several drawbacks, such as false negatives and dislocation of gradient maxima at great scales, and false positives at low scales. Nonlinear diffusion methods [33] permit noise elimination while preserve meaningful structures. They simulate a heat-spreading phenomenon in which diffusivity depends on local properties of the image. Perona and Malik [23] introduced a gradient dependent diffusion coefficient to stop diffusion at boundary points, eliminating boundary blurring but also maintaining noise at boundary points. Weickert [34] proposed a nonlinear anisotropic diffusion method to smooth surfaces just along the tangent plane at boundary points, reducing noise also at boundaries without blurring. The shortcoming of Weickert’s approach to anisotropic diffusion is that it causes a rounding effect on corners. Anisotropic diffusion has been frequently used in medical image processing [11, 14] [27]. In particular, Krissian et al. [14] have already applied it to surface extraction successfully. They developed a directional anisotropic diffusion technique to enhance surface extraction on vessel images, obtaining better results than the ones provided by Gaussian filtering. In their approach, diffusivity in the maximum curvature direction is annulated to reduce the rounding problem. This is made at the expense of noise reduction in the aforementioned direction. The second contribution of the present work is to introduce a corner-preserving anisotropic diffusion method, as a preprocessing tool to improve gradient and corner measure and detection. This is achieved by defining curvature dependent diffusion coefficients. The incorporation of the preprocessing step represents an incoming for both initialization, improving surface detection, and deformation too, since false minima elimination and correct placement of boundaries and corners improves the external potential definition. Moreover, the curvature based approach proposed in this paper leads to a better location of gradient potential minima and enhancement of the curvature based potential, as rounding problem is avoided. The paper is organized as follows. In next section, an overview of the complete initialization methodology is described. In section 3, the generalized ellipsoid model is defined. Section 4 presents the optimization technique used for model construction and matching. In section 5 all preprocessing steps are described. These are, image denoising, boundary surface patches detection and patch classification. Finally, in section 6 efficacy of the whole process is studied, paying special attention to its robustness in the presence of noise and loss of information caused by contrast variations over surfaces. 2. Initialization methodology The goal of 3D reconstruction methods is to obtain a detailed description of objects present in a volume data. To this end, deformable models are a good choice, because they possess great flexibility and guarantee certain smoothness and continuity properties. The problem with deformable models is that they interact locally with image features. It is necessary to define a good starting geometric configuration to ensure the success of the deformation process. The initial surface must be close to the boundary of the object of interest, which implies to determine position, orientation and global structure of that object in the 3D domain. Figure 1 illustrates the initialization process proposed in this paper. The attainmento of a high level description of certain object is separated in two stages: 1) modeling of global shape of the object class and 2) matching the class model with the instance object. The first stage is performed off-line over a set of sample shapes of the same object class (figure 1, step 1). The resulting class model is called the a priori model. The sample images are segmented manually to avoid typical automatic segmentation problems and then a surface prototype is extracted from them. Afterwards, global shape of the object surface prototype is determined by fitting a mathematical model to it. To represent global shape we have chosen the generalized ellipsoid model, also called superquadric. Once an a priori model is available, position, orientation and scale of an instance object in a new image is determined by a matching technique. To perform matching it is necessary to extract some feature from the image that describes the object surface properly and to find the transformations that lead to a correspondence between image features and model surface. In this work surface points are used as image features (figure 1, step 2). The extraction of such low level features from gradient information is not robust. Again, the use of prior knowledge is needed to distinguish the target object surface points from others present in the image as, for example, points belonging to structures other that the one under study or noise artifacts. To this end, points are not considered individually, but they are grouped in patches. Resulting patches are characterized by some of their average properties as, for example, area, contrast or shape descriptors, and these are used in a selection process to exclude undesired patches (figure 1, step 3). The starting surface for deformation is obtained after matching between model and object (figure 1, step 4). This is done by finding the parameters of the rigid transformation that minimizes an error function related to the distances from selected image points to the parametric surface. In short, what is presented here is a method that takes advantage from both bottom-up and top-down processes. Initialization is accomplished by obtaining higher and higher level descriptions of image contents: from volume points with associated gray levels, to boundary point representation, from boundary points to surface patches, featured by local and global descriptors that allow discriminating between desired and spurious patches, and from surface patches to a global surface model with the help of prior knowledge. The next stage, not described here, would walk the inverse way. Starting from the coarse representation of the surface resulting from initialization, the deformation process introduces local degrees of freedom to reach a detailed description of the surface object. Thus, the conflicting goals of high robustness and high resolution can be achieved. 3. Global shape models Both boundary detection and calculus of shape descriptors, based on curvature measures, require the computation of directional derivatives of gray level values. The derivative of Gaussian operator performs differentiation and smoothing simultaneously [19], but it modifies gradient maxima position, dislocating surfaces. Moreover, small structures can be eliminated. Nonlinear anisotropic filters offer better results. In this work, a corner preserving anisotropic filter has been developed to smooth 3D images without alter gradient and curvature values neither misplacing boundaries nor rounding corners. After image denoising, derivatives are approximated by central differences. Next subsections describe the main steps in feature extraction: anisotropic smoothing, boundary detection, and surface patch selection. 5.1 Anisotropic diffusion Let us introduce the general method for diffusion filtering in 3D proposed by Weickert [34]. Let Ω be a 3D image domain and ∂Ω its boundary. Given an image I(x, y, z), its filtered version u(x, y, z, t) is obtained by the next expression, with reflecting boundary conditions ( ) ( ) ( ) () ()() () ( ) () ∞×Ω∂ Ω ∞×Ω =∇ == ∇∇∇=∂ ,0 ,0 on on on 0,,,, ,,0,,, ,,,,,, ntzyxuC zyxItzyxu tzyxuuCtzyxu t (15) where n is the outer normal, 〈·,·〉 is the inner product and subscripts stand for partial derivatives. In the isotropic case, diffusion coefficient C is a scalar magnitude. Usually, it is a decreasing function of ||∇u|| with values belonging to the interval [0, 1]. In this way, diffusion is stopped in the presence of boundaries. Previous expression is often related to the energy minimization formulation. The energy functional E(u) is defined as the integral over the image of a potential Φ(||∇u||). Equation (15) is obtained by applying the gradient descent method to minimize the energy functional. Both approaches are related by () ( ) uuuC ∇∇Φ=∇ ' (16) Using this relation and operating with equation (15), next expression for nonlinear isotropic diffusion can be obtained ( ) ( ) ξξξξ uuuuut−∆∇Φ+Φ= ''' (17) where uξξ stands for the second derivative of u in the normal direction ξ. A possible choice for the potential function is the one proposed by Green [10] Φ(s)=( α 2/2)log cosh(s) (18) where α represents the gradient threshold at which diffusivity stops growing. Other approaches were stated by Perona and Malik [23] and Charbonnier et al. [5] among others. A diffusion process is called anisotropic when diffusivity takes different values {λ1, λ 2 , λ3} in different directions {e1, e2, e3} of space. As a result, diffusivity is a tensorial magnitude, and then, the flux vector C∇u is not parallel to the gradient direction. In the reference frame defined by the basis {e1, e2, e3}, diffusivity turns into a diagonal tensor D = diag(λ1, λ2, λ3). The expression for the diffusion tensor C in a general frame is C = TDT T where T is the matrix formed by the column basis vectors. Another way of constructing an anisotropic filter is generalizing equation (17), allowing independent diffusion coefficients for each term, as done by Krissian et al. [14]. The alternative expression is attained by taking the tensorial approach and setting the diffusivity eigenvalues such that λ i = Φi (||∇u||)/||∇u||, with i = 1, 2, 3, and applying some mathematical relations, obtaining ( ) ( ) 2211 321 '''' ηηηηξξ uuuuuut∇Φ+∇Φ+Φ= (19) Using this approach to filter an image for a given time t is faster than using the tensorial approach, because the second derivatives of u in the extreme curvature directions can be computed directly without calculating the Hessian eigenvectors, using the next relation 21 2211 kuukuu ∇−=∇−= ηηηη (20) where k1 and k2 are the maximum and minimum curvatures respectively. A typical scheme for anisotropic diffusion is constructed by setting diffusivities associated to tangent directions to constant values, while taking a decreasing function of ||∇u|| in the gradient direction. Thus, diffusion is stopped in the presence of boundaries only in the normal direction but not in the tangent plane, eliminating noise also at surfaces. 5.1.1 Corner preserving diffusion Smoothing in the tangent directions lessens curvature values as the system evolves. This effect eliminates noise by flattening surfaces, but also alters the shape of objects, eliminating small details and rounding corners. One way to prevent rounding is avoiding diffusion in the maximum curvature direction. In this way, the highest curvature value is not modified but, consequently, noise is maintained in the correspondent direction too. A better choice would be reducing diffusion only in the presence of corners, while keeping it at flat regions or surfaces where the tangent vectors variation is smooth. In this work, we propose such a diffusion method. To this aim, diffusivity in the maximum curvature direction has been modified in relation to the isotropic approach to introduce a dependency on a certain corner measure c. This diffusion coefficient is expected to be high in the absence of corners and decrease as corner measure grows. This behavior can be modeled with the Green function presented for contrast preserving filtering, just by changing the boundary detector by a corner detector and the gradient threshold by a corner threshold β . Therefore, diffusion coefficients in equation (19) are () ( ) 1'tanh'cosh'' 32 2 1=∇Φ=∇Φ∇=Φ −uccuu ββα (21) The corner detector is related to the local principal curvatures of the image. Curvature measures the variation of the tangent vector with the arc length in some direction, but it cannot be considered as a corner detector itself, because high curvature values can be originated by noise. To distinguish between real features and noise artifacts, curvature is usually multiplied by some power of the gradient magnitude. Different approaches for a corner detector with these characteristics are studied and compared by Sporring et al. [29]. Among the various possibilities discussed there, here it has been selected next max kuc ⋅∇= (22) As a result, an image point is considered a corner when both its gradient modulus and its maximum curvature have high values. Higher order powers of the gradient modulus can reject corners from structures with low contrast. The selection of the threshold values α and β is crucial in the accuracy of the obtained results, since they establish whether a feature must be smoothed or preserved. To obtain automatic parameter estimation, a tool from robust statistics [3, 24] is employed. The median absolute deviation about the median, MAD, is taken as a measure of the robust scale σ e of some magnitude. If medianI is the median of some magnitude computed from all points belonging to image I, then the robust scale for gradient is () ( ) [ ] 6745.0medianmedian6745.0MAD III IIe ∇−∇=∇= σ (23) Constant 0.6745 is the MAD of a zero-mean normal distribution with unit variance. Scale σ e is the contrast value at which flux must stop growing. If the stopping criterion is to reach a certain fraction x of the asymptotic limit of the isotropic flux function f∞, parameter α can be related to σ e by ∞ ⋅ = = ∇ fxIf e),||(|| α σ . For the Green function f∞ = α , so the threshold parameter is ( ) x eatanh σ α = (24) Here, it has been taken x = tanh(1) = 0.7619, so that α = σ e. The same estimation can be done for β , computing MADI (c). 5.2 Boundary detection Once the image is smoothed, directional derivatives can be computed using a central finite differences scheme. An image point is considered a surface point if it is a local maximum of the gradient modulus. Monga and Benayoun [19] accomplish gradient maxima detection by comparing the modulus magnitude of each point r only with values correspondent to previous and next points in the gradient direction ∇u(r), represented by r+ and r– and calculated by ( ) ( ) rrrr uu ∇∇±= ±/ (25) When ( ) ( ) ( ) { } _ ,rrr uuu ∇∇>∇ + (26) r is a gradient maximum. If r+ or r– do not coincide with an image position, gradient modulus is approximated by trilinear interpolation in a vicinity of r. 5.3 Surface patches selection After gradient maxima detection, image contents are described by a set of candidate boundary points. Now, it is necessary to determine what points are to be used in the fitting process. To this end, a collection of surface points is not an appropriate description of image objects, since points do not carry information about which object they belong to and what is the global aspect of that object. Simple thresholding of individual point attributes, as gradient modulus or curvatures, can cause discontinuities on relevant surfaces due to local fluctuations. Furthermore, not only structures under study are recovered, but also undesired surfaces or noise artifacts can appear. For those reasons, boundary point representation of image objects is replaced by a higher level description. Extending the idea pointed by Pardo and Cabello [22] from 2D to 3D, gradient maxima are grouped in connected components to obtain surface patches. Global information about objects can be extracted from this new description as, for example, area, contrast or shape. These global features are used to determine the value one unique label for each point that indicates whether it belongs to the object surface or not. In a more general case this label may represent the degree in which that point can be said to belong to the surface object. This is, each point r i belonging to certain surface patch P j is characterized by a label value determined by a labeling function L such that L (r i ) = L ( P j ), ∀r i ∈ P j , or what is the same, all points in a patch have the same label value and this value is determined from the contributions of all individual points. Labeling function dependency on surface patch attributes is determined according to prior knowledge about image contents. L must show great values for patches with global descriptor values similar to the ones expected a priori for that object class and low values for any other patch. Therefore, an expression for L must be constructed from those a priori descriptor values for each object class. To accomplish grouping, a 26-connectivity criterion is used. Many of the points detected as gradient maxima do not correspond to surfaces of the desired objects in the image, but they are noise artifacts. When grouping boundary points to construct surface patches, those points must not be considered. Simple thresholding is not a good technique to distinguish noise artifacts from real surface points, as contrast, in general, is not uniform over the object surface. Here, hysteresis thresholding is used. Hysteresis involves defining two threshold levels. The lowest threshold level determines whether a point belongs to a surface. If this level is chosen appropriately, surface fragmentation is reduced. In addition, each connected component must have, at least, a number n of points with gradient greater or equal to the other threshold level. This is supposed to exclude surface patches originated by noise. Here, minimum number of points over the highest threshold level has been taken n = 1. Our goal is to construct a representation of the image that emphasizes salient locations. We seek to associate a measure of saliency, denoted by the labeling function, to each surface patch. A property that seems to play an important role in boundary saliency is the combination of size, global and/or local (smoothness) shape, and contrast. A labeling function that would account for our working examples is one that favors long, smooth shape and high gradient surface patches. In our proposal, smoothness is related to curvature type or curvature variations. The exact formulation of L can be adapted to the specific application domain. Here, working hypotheses are that true surfaces have higher area and higher average gradient level than noise artifacts. In addition, they can be useful to discriminate among various anatomical structures with known properties. For example, cortical bone tissue in CT images is characterized by its high contrast in relation to muscle and trabecular bone tissues. Relative sizes of objects are also known in general. Shape descriptors are used to discriminate anatomical structures with different morphologies. In this work Gaussian curvature K and mean curvature H are used as shape descriptors. () −+−+ ⋅ = + =kkKkkH 2 (27) where k+ and k– are the extreme curvature values at each point. Using their values, next classification of points can be made: H > 0 H = 0 H < 0 K > 0 concave elliptic − convex elliptic K = 0 concave cylindrical plane convex cylindrical K < 0 concave hyperbolic saddle point convex hyperbolic Surface patch descriptors can be obtained from point descriptors by averaging their values. However, the error committed in the calculus of K is the product of the errors correspondent to k+ and k–. As curvature values are very sensitive to noise, results obtained for K are not very reliable. Better results should be obtained by averaging extreme curvatures and then using resulting values in equation (27), at least for the Gaussian curvature sign, despite mean(k+ · k–) ≠ mean(k+)·mean(k–). Average Gaussian and mean curvature signs classify surfaces in a very rough manner. Their utility is limited to simple objects. If maximum and minimum curvature signs vary strongly over the surface, average measures are no longer descriptive of surface shape. Let us think, for example, in the shape of a vertebra. In these cases, the curvature variations along the surface patch can be used as a measure of smoothness. Labeling function must depend on those patch attributes in such a manner that target surfaces have high label values while the remainder patches have low values. Labels can be used in several ways to decide how the fitting process is to be done. A general method involves considering labels as weight factors in the calculus of the error function, modulating the contribution of each surface point to the total error. Then, the error function can be redefined in the next way () () ( ) ∑ = =N i ii DLE 1 22 ,qrrq (28) Then, if a point belongs to a salient patch, its contribution to the accumulated distance is decisive in the result. As label value decreases, the influence of all points on the patch is reduced. 6. Results Performance of the whole methodology depends on three aspects: the capability of the chosen model to represent objects in an appropriate manner in certain field of application, the efficiency of the optimization method and the effectiveness of the preprocessing technique in surface patches extraction. In this section, the initialization method is tested taking all this points into account. 6.1 Modeling with superquadrics In section 4.1, several error metrics were presented. Each of them represents a fitness measure in certain metric. To make a comparison between results obtained with each one, it is necessary to establish a metric-independent quality criterion. The procedure is as follows. A synthetic superquadric surface, illustrated in figure 4, is designed with vector parameter q, representing it as a point cloud. Afterwards, parameters are estimated optimizing three different error functions: D2, D3 and D4 correspondent to equations (11), (12) and (13) respectively. Then, parameters q’ estimated with each error function are compared with real values to determine which metric to use. The AG scheme employed to fit the surface has the following characteristics. Parameter vector has been represented with Gray code of 16 bits. Population size has been set to 80 individual for modeling and 40 for matching. The error function is mapped by linear rank before selection of candidates to be used by genetic operators. Selection method is probabilistic tournament. Genetic operators applied are two-point crossover and mutation with probabilities pc = 0.8 and pm = 0.1 respectively. The best 20% individuals are reproduced in next iteration of the algorithm. To compare real and estimated parameters, the discrepancy d between the two surfaces has been defined. It represents the Euclidean distance between the two parameter vectors. To assign the same weight to each parameter in the distance measure, each component qi is scaled with a factor α i. Scale factors are determined heuristically and are related to the range of variation of each parameter. Then, d is () ( ) ∑−= i iii qqd 22 ', α q'q (29) Parameters of the synthetic superquadric are shown in the first column of table 1. The second column contains the weight factors for each parameter. The remainder columns show parameter vectors estimated with three different error functions. They have been () () ( ) ( )()     <><> = otherwise0 and,0,if1 minmax 2/12 thjthjjj j kPkkPkPHlaPg PL (30) where a is the patch area, g(Pj) is the gradient modulus averaged over all patch points and l is a threshold value estimated from typical noise, muscle and bone gradient and area values. Here it has been taken l = 1e5. The use of this labeling function as a weight factor in the fitting process is equivalent to the elimination of less salient patches. Selected points, presented in forth row of figure 19, contribute equally to the error function. The surface is not complete due to attenuation of the tibia intensity level in the knee region. The matching process has been realized using the geometric model obtained in section 6.1. Results in fifth row of figure 19 show that initial model is near the surface to be modeled. As this model has been obtained from the same image by manual segmentation, this provides a tool to measure the correctness of the rigid transformation parameter estimation. Results in table 5 show parameters q obtained in the model construction phase and parameters q’ obtained by matching the model to the processed image. If the disparity measure between both parameter vectors is computed, the result is d = 0.06839, which a good result for a relatively complex surface. MRI Images are characterized by their high signal to noise ratio (SNR). Fist row of figure 20 shows an example of an MRI image of the aorta artery. It can be seen that it is very noisy and the average contrast is low. In addition, other structures are present in the image besides the aorta. The application of anisotropic diffusion of scale σ = 5 enhances the MRI image (figure 20, 2nd row), given that noise is almost completely eliminated but boundaries are not blurred. Threshold parameters are α = 2.7 and β = 0.42. Because of the low contrast of the image, high threshold must be also low. Here, it has been taken t1 = 40. Low threshold must be relatively high, t2 = 25, to avoid connection of distinct surfaces. Connected components are shown in third row of figure 20. Again, thresholding is now enough to obtain desired surface and patches selection is required. The same criterion used in the previous example is used here to extract aorta surface, which is also cylindrical, but concave now. () () ( ) ( )()     <>>> = otherwise0 and,0,if1 minmax 2/12 thjthjjj j kPkkPkPHlaPg PL (31) For this kind of MRI images, it is taken l = 1e4. After selection of patches, results in forth row of figure 20 are obtained. Patches that do not belong to the aorta have been eliminated, but part of the surface has been lost during preprocessing. A priori model has not been extracted from test images this time. A simple cylinder model is enough to represent roughly the artery shape. Results of the matching process are shown in fifth row of figure 20. 7. Conclusions Deformable models are very reliable modeling techniques. They force continuity and smoothness in segmentation but they have a local field of activity, so starting configuration determines the success of the procedure. In this paper, an automatic initialization tool to guide deformation is presented. To this end, a complete methodology has been developed to obtain a description of global shape, location, orientation and size of an object from a 3D image. Object identification is accomplished by combining high level information from images and introducing a priori knowledge. A methodology, based on a priori knowledge about surface features, has been designed to isolate surface patches belonging to the target object. This method permits to eliminate structures that correspond to noise or other objects present in the image. Modeling with superquadrics does the remaining work. Surface patches are used only to determine rigid transformation parameters, and fine tuning of learned geometric features. The principal advantage of this technique is that initialization is highly automated and it provides a good approximation between surfaces of the desired object and its model, ensuring proximity to the correct energy minimum in a deformable model scheme. Generality is also important. A priori models can be easily constructed in any application domain, so that a detailed anatomical atlas is not necessary. Description with implicit surfaces simplifies the estimation of matching error, since it is not necessary to determine correspondences between model points and image points. The discretization and bounding of the solution space allows using a GA to find the global optimum without an excessive time cost. In future works we expect to supply the technique with mechanisms to detect different components of multipart or branched objects. Thus, each part can be initialized with a different model in a hierarchical scheme. The other contribution in this work is the corner preserving anisotropic filtering. Curvature dependent diffusion coefficients have been designed, so that diffusion is stopped at corner points in the tangent direction correspondent to the maximum curvature level. Therefore, noise is eliminated at every image region, inter-region points included, without dislocate boundaries. This fact implies an important improving for surface detection, since shape is now more reliable. Surface merging also is avoided. Furthermore, this processing technique is expected to improve the measures for the external energy of the deformable model. Results presented show that both gradient and corner detectors offer better responses after anisotropic filtering in comparison with Gaussian blurring. The effect of the definition of the threshold parameters is very relevant for this application. The automatic setting of these parameters provides a useful tool to find a compromise between smoothing and boundary or corner preserving. However, there are cases in which the definition of global thresholds is insufficient. Variations on the background intensity level, attenuation of the features intensity or the presence of structures of different intensities can be the reason of important loss of information. To solve the problem, the threshold parameters should be estimated locally on a vicinity of each point. This is a proposal for future works, where viability of the approach must be studied in terms of computational efficiency. Bibliography [1] Bajcsy, R and Kovacic, S, Multiresolution elastic matching, Computer Vision, Graphics, and Image Processing, 46 (1989) 1-21. [2] Bardinet, E, Cohen, L D, and Ayache, N, A parametric deformable model to fit unstructured 3D data, report nº 2617 (INRIA, Sophia-Antipolis, France, 1995). [3] Black, M, Sapiro, G, Marimont, D and Heeger, D, Robust anisotropic diffusion, Trans. on image processing, 7(3) (1998) 421-432. [4] Bolle, R M and Vemuri, B C, On three-dimensional surface reconstruction methods, IEEE. Trans on Pattern Analysis and Machine Intelligence, 13(1) (1991) 1-13. [5] Charbonnier, P, Aubert, G, Blanc-Feraud, M and Barlaud, M, Two deterministic half-quadratic regularization algorithms for computed imaging, IEEE Int. Conf. On Image Proc., (Austin, Texas, 1994) 168-172. [6] Cootes, T, Hill, A, Taylor, C, and Haslam, J, The use of active shape models for locating structures in medical images, in: Barrett, H H and Gmitro, A F, eds. Information Processing Medical Imaging, LNCS 687, (Springer-Verlag, Berlin, 1993), 33-47. [7] Delingette, H, Montagnat, J, General deformable model approach for model-based reconstruction, Proc. of the IEEE International Workshop on Model-Bases 3D Image Analysis, (Bombay, India, 1998). [8] Dosil, R, Elipsoides xeneralizados activos: aplicación á segmentación de imaxes médicas 3D, MA Thesis, (Universidade de Santiago de Compostela, Spain, 2000). [9] Goldberg, D E, Genetic algorithms in search, optimization and machine learning, (Addison-Wesley, 1999). [10] Green, P J, Bayesian reconstruction from emission tomography data using a modified EM algorithm, IEEE Trans. on Medical Imaging, 9 (1990), 84-93. [11] Gerig, G, Kübler, O, Kikinis, R and Jolesz, A, Nonlinear anisotropic filtering of MRI data, IEEE Trans. on Medical Imaging, 11(2) (1992) 221-231. [12] Jaklic, A, Leonardis, A and Solina, F, Segmentation and Recovery of Superquadrics, Computational imaging and vision, Vol. 20, Kluwer, Dordrecth, 2000. [13] Kass, M, Witkin, A, Terzopoulos, D, Snakes: active contour models, Int. Journal of Computer Vision, 1(4) (1988) 321-331. [14] Krissian, K, Malandain, G and Ayache, N, Directional anisotropic diffusion applied to segmentation of vessels in 3D images, report nº 3064, (INRIA, SophiaAntipolis, France, 1995). [15] Leonardis, A, Jaklic, A, Solina, F, Superquadrics for segmenting and modeling range data, IEEE. Trans on Pattern Analysis and Machine Intelligence, 19(11) (1997) 1289-1295. [16] Linares, P, Rodríguez, L and Montilla, G, Genetic algorithms fitting of deformable superquadrics applied to left ventricle visualization, Computers in Cardiology, 25 (1998) 657-660. [17] McInerney, T and Terzopoulos, D, A finite-element model for 3D shape reconstruction and nonrigid motion tracking, ICCV’93, (1993) 518-523. [18] McInerney, T and Terzopoulos, D, Deformable models in medical image analysis: a survey, Medical Image Analysis, 1(2) (1996). [19] Monga, O, and Benayoun, S, Using partial derivatives of 3D images to extract typical surface features, report nº1599 (INRIA, Rocquencourt, France, 1989). [20] Monga, O, Deriche, R, Malandain, G and Cocquerez, J-P, Recursive filtering and edge closing: two primary tools for 3D edge detection, report nº 1103 (INRIA, Rocquencourt, France, 1989). [21] Neuenschwander, W, Fua, P, Székely, G and Kübler, O, Velcro surfaces: fast initialization of deformable models, Computer Vision and Image Understanding, 65(2) (1997) 237-245. [22] Pardo, X M and Cabello, D, Biomedical active segmentation guided by edge saliency, Pattern Recognition Letters, 21 (2000) 559-572. [23] Perona, P and Malik, J, Scale-space and edge detection using anisotropic diffusion, IEEE. Trans on Pattern Analysis and Machine Intelligence, 12(7) (1990). [24] Rousseeuw, P J and Leroy, A M, Robust regression and outlier detection, (New York: Wiley, 1987). [25] Saint-Marc, P, Chen, J-S and Medioni, G, Adaptive smoothing: a general tool for early vision, IEEE Trans. on Pattern Analysis and Machine Intelligence, 13(6) (1991) 514-529. [26] Sandor, S and Leahy, R, Surface-based labeling of cortical anatomy using a deformable atlas, Medical Image, 16(1) (1997) 41-54. [27] Solé, F, Ngan, S-C, Sapiro, G, Hu, X and López, A, Anisotropic 2D and 3D averaging of fMRI signals, IEEE Trans. on Medical Imaging, 20(2) (2001) 86-93. [28] Solina, F and Bajcsy, R, Recovery of parametric models from range images: the case for superquadrics with global and local deformations, IEEE Trans. on PAMI, 12(2) (1990) 131-146. [29] Sporring, J, Nielsen, M, Weickert, J and Olsen, O F, A note on differential corner measures, Proc. 14th Int. Conf. Pattern Recognition, IEEE Computer Society Press, Los Alamitos, 1 (1998) 652-654. [30] Staib, L H, Chakraborty, A and Duncan, J S, An integrated approach for locating neuroanatomical structure from MRI, Int. Journal of Pattern Recognition and Artificial Intelligence, 11(8) (1997) 1247-1269. [31] Tek, H and Kimia, B B, Volumetric segmentation of medical images by threedimensional bubbles, Computer Vision and Image Understanding, 65(2) (1997) 246-258. [32] Terzopoulos, D, and Metaxas, D, Dynamic 3D models with local and global deformations: deformable superquadrics, IEEE Trans. on Pattern Analysis and Machine Intelligence, 13(7) (1991) 703-714. [33] Weickert, J, A review of nonlinear diffusion filtering, in: ter Haar Romeny, B, Florack, L, Koenderink, J and Viergever, M, eds., Scale-Space Theory in Computer Vision, Lecture Notes in Comp. Science, 1252, (Springer, Berlin, 1997) 3-28. [34] Weickert, J, Scale-space properties of nonlinear diffusion filtering with a diffusion tensor, report nº 110, (Laboratory of Technomathematics, University of Kaiserslautern, Germany, 1994). [35] Whaite, P and Ferrie, F P, From uncertainty to visual exploration, IEEE Trans. on Pattern Analysis and Machine Intelligence, 13(10) (1991) 1038-1049. Captions Figure 1. Sequence of steps to achieve the starting configuration for the deformable model. Figure 2. Superquadric surfaces for different ε1 and ε2 values. Figure 3. Global deformations. Figure 4. Synthetic superquadric used to test fitting method. Figure 5. Tibia prototype to be modeled. Figure 6. Tibia model obtained with error function D4. Figure 7. Tibia model obtained with error function D3. Figure 8. (a) Original test image, consisting of a cube (1) a cylinder (2) and a sphere (3). It has been smoothed with (b) Gaussian filter (c) corner preserving anisotropic diffusion, both with scale σ = 6. Slices represent plane z = 30. Figure 9. Gradient at slice z = 30 of (a) Gaussian filtering (b) anisotropic filtering. Figure 10. Boundaries detected at slice z = 30 with (a) Gaussian filter (b) anisotropic filter. Figure 11. Corner measure at slice z = 30 of (a) Gaussian filtering (b) anisotropic filtering. Figure 12. Displacement of corners with filtering of scale σ = 6. Corners are labeled by quadrant. Figure 13. Dislocation of sphere boundary through diffusion time. Figure 14. Extreme curvatures estimation at different scales of the Gaussian filter. Figure 15. Extreme curvatures estimation at different scales of anisotropic diffusion. Figure 16. Different sections of a synthetic image to test patch classification. Figure 17. Different selections of surface patches. (a) Non-planar | kmax | > kth (b) Planar | kmax | < kth (c) Concave H > kth (d) Convex H < −kth (e) Elliptical | kmax | > kth & | kmin | > kth & K > 0 (f) Hyperbolic | kmax | > kth & | kmin | > kth & K < 0 (g) Cylindrical | kmax | > kth & | kmin | < kth (h) Selection of individual cylindrical points Results are presented for slice z = 25 except in (h), which shows z = 40. Figure 18. 1st row: test image with Gaussian noise of variance σ n = 5 and variation of the background color along y axis. 2nd row: Image smoothed with anisotropic diffusion with scale σ = 5. 3rd row: Selection of cylindrical patches. 4th row: cylinder model –resulting from matching with extracted patch –superimposed to original image. Figure 19. 1st row: Slices of the original tibia CT image. 2nd row: image filtered with anisotropic diffusion of scale σ = 5. 3rd row: convex cylindrical patches selected with u1 = 20 and u2 = 50. 4th row: initial model superimposed to original image. Figure 20. 1st row: Slices of the original aorta MRI image. 2nd row: image filtered with anisotropic diffusion of scale σ = 5. 3rd row: concave cylindrical patches selected with u1 = 25 and u2 = 40. 4th row: initial model superimposed to original image. Table 1. Columns from left to right: synthetic superquadric parameters, correspondent weight factors, average estimated parameters qi’ obtained with different error functions. Figure 14. Figure 15. y = 50 z = 25 x = 70 x = 150 x = 213 Figure 16. Table 3. Patch k1 k2 H K Plane 0.001219 −0.006986 −0.002883 −0.000009 Plane 0.001622 −0.007492 −0.002935 −0.000012 Convex hyperbola 0.010841 −0.036092 −0.012626 −0.000391 Concave hyperbola 0.040352 −0.015881 0.012236 −0.000641 Convex cylinder −0.000721 −0.038698 −0.019710 0.000028 Concave cylinder 0.052256 −0.009572 0.021342 −0.000500 Convex sphere −0.042745 −0.047356 −0.045050 0.002024 Concave sphere 0.074956 0.069251 0.072103 0.005191 (a) (b) (c) (d) (e) (f) (g) (h) Figure 17. z = 50 x = 50 y = 50 Figure 18. q i α i q i’ gauss q i’ anisotropic a 0 1 1 9.758e−1 1.069e+1 β 0 0.159 −1.203e−1 1.452e−2 t 1 50 0.02 5.011e+1 5.032e+1 t 2 50 0.02 5.013e+1 5.203e+1 t 3 50 0.02 4.924e+1 5.122e+1 Table 4. z = 50 z = 140 x = 98 y = 112 Figure 19. Table 5 q i q i’ a 0 1.00e+0 9.54e−1 α 5.31e−1 3.88e−1 β −9.10e−2 −1.08e−1 γ −1.27e−1 1.03e−1 t 1 9.65e+1 9.77e+1 t 2 1.04e+2 1.04e+2 t 3 8.73e+1 8.78e+1 z = t3 x = t1 y = t2 Figure 20.