Shape Deformation: SVM Regression and Application to Medical Image Segmentation Song Wang, Weiyu Zhu and Zhi-Pei Liang Department of ECE and Beckman Institute University of Illinois at Urbana-Champaign, IL 61801, U.S.A.
[email protected],
[email protected],
[email protected] Abstract This paper presents a novel landmark-based shape deformation method. This method effectively solves two problems inherent in landmark-based shape deformation: (a) identification of landmark points from a given input image, and (b) regularized deformation of the shape of an object defined in a template. The second problem is solved using a new constrained support vector machine (SVM) regression technique, in which a thin-plate kernel is utilized to provide non-rigid shape deformations. This method offers several advantages over existing landmark-based methods. First, it has a unique capability to detect and use multiple candidate landmark points in an input image to improve landmark detection. Second, it can handle the case of missing landmarks, which often arises in dealing with occluded images. We have applied the proposed method to extract the scalp contours from brain cryosection images with very encouraging results.
1 Introduction Image segmentation is an important task in image understanding and computer vision. Traditional methods are usually based on edge detection, region merging/splitting, statistical pixel classification or hybrid methods [19]. Although these methods are capable of partitioning an image into regions of connected pixels according to image intensity distributions, more often than not, the segmented regions do not correspond well to true objects. This limitation is particularly serious for medical image segmentation because medical images are often very noisy and the structures to be identified have large “internal” intensity variations. An effective way to solve this problem is to incorporate the known shape information of the desired objects into the segmentation process. Generally, the shape of an object is described by the boundary contour (or several closed contours) of the objects, which can be represented by a sequence of sampled landmark points. A number of shape-based segmentation methods have
been proposed. For example, the active contour (or snakes) [8, 1, 17, 4] model constrains the desired shape to be sufficiently continuous, smooth and differentiable during the deformation process. Although the original active contour cannot process multiple contours, the problem can be solved using the recent geodesic contour method [3]. More advanced shape models include the point distribution model (PDM) by Cootes et al.[5], which can learn shape variations from a set of segmented training images, each containing a set of landmarks to define the shape. PDM calculates the covariance matrix of these landmarks and then represents the shape variations along the directions of the most significant eigenvectors of the covariance matrix. However, it is a tedious job to build the training set manually and in some cases it is also difficult to identify the set of landmark points for all the shapes of interest. Similarly to PDM, Staib et al. [13] decompose the shape into items with different frequencies. Then the shape knowledge is learned by studying the distribution of each Fourier coefficient. The result is used to constrain the shape deformation. Jain [7] presents a similar approach which parameterizes the shape and uses the distribution of the parameters as the constraint for shape deformation. Along the same line, Leventon et al. [10] incorporate statistical shape information into geodesic active contours to segment medical images. The work presented in the paper is more closely related to the deformable models developed by Rueckert et al. [12] and Zhong et al. [18]. Similar to the active contour, the energy function of those models consists of two terms accounting for external and internal energy respectively. The external energy is a potential field describing the edge and region matching information between the template and the input image. The internal energy represents the shape deviation from the prior. However, it is not easy to determine the appropriate weight for the internal and external energy. Additionally, this energy function is nonlinear and includes too many unknowns to be well handled by gradient-based algorithms. This paper presents a new landmark-based shape deformation method which can effectively avoid the above ques-
tions. At first, edge- or/and region-based methods are applied to detect candidate landmark points in the input image within the neighboring area around each landmark in the template. Multiple candidates may be detected for one landmark when it is necessary. Then the template is deformed using a special support vector machine (SVM) regression technique. Based on the fitting result, a set of better landmarks is identified from the detected candidates. This process is repeated until the final result is obtained. Usually, a couple of iterations is sufficient. The rest of the paper is organized as follows. In Section 2, we describe the proposed shape-based deformation method. Section 3 discusses the selection of the regularization parameter and how to handle missing landmarks. Section 4 presents representative results for medical image segmentation, followed by the conclusion in Section 5.
2 The Proposed Method 2.1 Problem formulation For simplicity, consider the case in which the template contains only a single closed contour. Extension to the case of multiple contours is straightforward. The contour is represented as a series of ordered landmark points , where are the coordinates of the -th landmark. For each landmark , the proposed method first identifies a set of possible corresponding in the inlandmark points put image, where . Then the problem is solved in two steps:
based on the acquired ; . This process is repeated until !convergence. There are three main subproblems de in this approach:3 detection of landmark candidates ! in the input image, termination of the regression model and updating based on the regularized fitting results. Solution to these problems are discussed below.
2.2 Detection of landmark candidates
] 8 3
Given , we assume that in an input image falls in a circular area of radius centered at . According to the principle of anisotropic diffusions [11], non-rigid biologic shape deformation between the template and the input image can be done effectively by moving each landmark along the normal direction of the shape. This means that we can search for the corresponding only along the normal directions of as shown in Fig.1. Denote as the
73
^
n1 v 1
n2 v2
nn
v1 v2
vn vn
! #&%"' 2$#& %'#&% )'(* , #&%' +- /.0 1 Figure 1. Search for landmarks along the normal direction of the initial shape. 3 843 583 3 63 : from 1. Identify the best landmark point 73 8 9 83 unit vector along the normal direction at . Combining the landmark set ! such that 3 can the above assumptions, the corresponding landmark located in or near the real object boundary in the input ` H ]a^ ]cb be image; eL d ] obtained ] N . from the line segment _ 3 2. Deform the prior shape to match while keeping Because of noise and structural complexity, 83 from _ sometimes . In this the general shape characteristics of . it is difficult to accurately identify f #&%' (W case, we can extract several candidates in ! The template deformation is described asB 7 a C regression + Z 0 . [ 1 < > = 0 ? A @ B problem of finding a transform ; that 83 and the proposed algorithm will determine an optimal from ! . minimizes For different application methods can + D 5 3 8HJIKML g+ Zproblems, .0 : 1 . various be used to locate ! Generally, they can 1 FE8G ; ;ON (1) be categorized into two classes. detection methods. Since the desired shape is the QPF P&ARTS isKMaL lossR9function, IURVS is a regulariza- 1. Edge where G contour of an object in the input image, edge points S ; N tion parameter and is a regularization functional. in will be collected to construct . Any edge de_ ! W X < 5 Y 2 + Z 0 . [ \ 1 ; In this setting, ; is the detection algorithm can be applied here. A simplest apsired object shape in the input image. proach is to calculate the gradient vector for each pixel Generally, it is a combinatorial problem to determine the along _ . If its amplitude is very large (according 73 to a best . In this paper, we use an EM-like (Expectation Minpre-selected threshold), it will be included in ! . imization) two-step algorithm to solve it iteratively. Specif ically, we perform the shape-regularized fitting to an initial 2. Region matching methods. For any pixel along _ , 73 to get 73 an optimal ; and then update from detected if its neighboring area has similar features with the
!
neighboring area centered at in the template image, it will be included in . A number of matching metrics can be used for this problem, which include normalized cross-correlation, squared differences in pixel intensities, measures based on optical flow models and mutual information metrics [14]. Generally, region matching is often more robust than edge detection to process noisy images containing complex objects. However, region matching methods usually require the template landmarks to be very accurate, which is not necessary for edge detection methods. Selection of is related to the similarity between the template and the input image. For segmenting 3D medical images slice by slice, this is usually determined by the sampling interval between two neighboring slices. The larger the interval, the larger the . In practice, should also be adaptively modified based on the size of . If there are too many candidates in , then needs to be reduced, and vice versa. It is also possible that no landmark candidate is identified for some . This problem will be discussed in Section 3.2.
]
!
]
]
]!
2.3 SVM regression model for deformation
the insensitive margin parameter , landmark points can be identified as support vectors or non-support vectors. This result can be used to improve . This step is discussed in detail in Section 2.4. The -insensitive function is defined by
W3
nO 3 d ; 5 )n SnY3 d 5 rnd if nO8 3 d ; 5 rnA
;
else cS , this is just the modulo loss function. When
The only remaining component of the regression model is the additional constraints that should be placed on the deQ & P sired deformation ; . Obviously, that the desired ; lie in the normal direction of wethe hope shape centered at ,
; 5 m H ^ + Z.0 [ 1 (3) 0 B where W3 b is. detected along the normal directions of , Since we have 3 HJ ^ Q +- /.0 : 1 (4) B + / h . [ 1 where that : 1 , theNoting ; = Q?h b and n ^ /nv+- "+ Zare.0 known. regression i.e.,
73 from ! Q estimated landmarks + /.hGiven [ 1\a set, thisof section describes a regression model for problem can be rewritten as (1) and the corresponding algorithm. Generally, we expect Problem 1: the model to have the following properties: (a) 3 ; 5iregression P + is a good approximation for ; (b) ; is a geometD $ / ; 1 FE d HJIKL ;ON rically homologous mapping (an important feature in bio K L logic deformation); (c) ;ON represents biologic deKML ;N the 5i ; (d) well formation between and ; is invariant un9l Q0 Q0hQ . subject to constraints (3), where derQPF the affine or rigid transforms; and (e) The loss function P& is robust against noise andKML outliers; The thin-plate kernel is conditionally positive definite G From and the linear functions form the null space of the resultregularization theory, ;N can be defined as a norm 9 > = 0 ? ing , which must be of the form [2] ; in a reproducing kernel Hilbert space (or subspace) which =\k0kH*6 OH*6 OXHJ¡ ¡ FE[¢ j m ) can be uniquely represented by a positive definite (or con5 k 8 l 3 ?:5Mk£$Hy£Y ¤HJ£Y OXH E ~- j m )Y (5) ditionally positive definite) kernel function j . For [ F s t u medical image segmentation, we propose to use the thin5k 83lmonYpdq83rn nOvdi83rn , where + m plate kernel [2] j KL ;ON , also To solve Problem , we first fix and minimize ; named the bending energy, is given by with respect to ; under constraints (3). This is~-equivalent to } { > : = \ H | { 5 0 ? Q ~ x seeking minimizing subject to ; ¥ ¥ z KML ;N Vw$wyx >{|=: H*{}l?hQ~-~-
¦ / 6 ¦ > m
) | H 0 x (2) the constraint that ; must map to ^ , + Z.0 [ 1 . It is a typical zx thin-plate interpolation problem ` h a / 6 6 Q ` > £ a / Y £ a Z£Y HJ.h Hk } { Q & P m k andª the parameters , § ¨ where 583 5 ) k9. nY83 d 5 )n © ¢ ¢ ¢ h and « ª~0 /~- - : /~-h in (5), A popular loss function is G ; ; . can be calculated by solving the following equation: In this case, the model becomes the traditional thin-plate ¬®¯ ¯ ¬ © « ¬É\ ÉM : Introduce the Lagrange multipliers Éh ¦ , È ¦ ʦ É$¦ ¦ a ɦ ¦ ɦ h , Ë ÌÈ 5Í: Í 8 Íh ® Í: Í : Íh , the Lagrangian function for and Ë Problem Î is ÏÐÑ Q¿: ¿ ¦ È È ¦ Ë Ë ¦ m Ò DE ÀH À¦ )8H .+ Á ³ Hyà m » ³ Hyà m H¤ ´ H*à m » ´ Hyà k Ä d DFE É/
H*À:d6HJ Dd ɦ /
H À¦ HU6:dUMd D 5ÍÀ[H ͦ À¦ ) FE FE subject to
Ò ÓZÔÕ
¿¦
where . Imposing the first-order necessary conditions (setting the derivatives of the Lagrangian with respect to , and to zero) yields
and
Ï
UÖ È ¦ d È MdØ× É HUÍ ¦ Ò + Z.0 [ 1 É ¦ H Í ¦ Ò + Z.0 [ 1 É É ¦ Í Í RJS + Z.0 [ 1
¿
(9)
Ö Ã » à yH à » à z × ÖUà » ³ H*à » ´ O Substituting these results into Problem Î gives the following equivalent formulation Problem 4: Ù Ù ¾ ÏЦ È È ¦ k .+ È d È ¦ ÖU È d È ¦ H¤5×WH*Ú\ È d È ¦ 8HyÛ È H È ¦ (10) where
S JÉ É ¦ Ò Ü+- /.h [ 1 ÚÜ / : Z and Û¤Ü
:
. where Problem Ý is a simple quadratic programming problem which can be After the ¦ solved using standard algorithms. optimal È È are obtained, the desired in Problem Î can be calculated from (9). Finally, we apply (6) and (5) to get the desired ; . 73 2.4 Updating The next step to improve the W3 infromour ! iterative Q Þ+algorithm /.h [ 1 isbased estimation of on the results from the shape fitting step. To do so, we invoke the Kuhn-Tucker condition for Problem Ý , which states É 5 d d
dØÀ ¦ mS + Z.0 : 1 É ¦ 5 Ò d H
H À mS + Z.0 : 1 dÉ)rÀ¦ MS + Z.0 : 1 Ò d¶É¦ ) ÀMS + Z.0 [ 1$ É PmÉ ¦ É ¦ 2 Sh "+- /.h [ 1 . The input It is obvious that É data for which or have values are called is nonzero support vectors. Otherwise, a non-support vector. As shown if is a non-support vector, then 83 inliesFig.2,
-insensitive the current inside the margin centered 5i ,8i.e., V d
at the resulting ; . From (10), it 3 will have no effects also easy to see that on the final subject to
regression and has no contribution to the bending energy. In
ni
t( vi )
ε-insensitive margin
t( vj )
vi vi
vj
3
%
3 Discussion
and non-support
I
3.1 Selection of
this case, the current can be regarded as being accurate enough with the permissible error margin . On the other hand, any support vector lying outside the margin need to be improved if possible. Based on this consideration, the following simple exchanging strategy is adopted to update .
3
, we leave the current 3 un 2. For support vector , if another candidate in ! lies in 83 .the
-insensitive margin, then we adopt it as the new This strategy will reduces the bending energy if
is fixed. The remaining problem is the selection of
. Inspired by the area-shrinking technique used in SOM (self-organizing mapping), we reduce
gradually in the iteration process until it approaches zero. This is quite reasonable since the ac5i is usually quired ; inaccurate of the 3 at theW3 abeginning iteration for some mis-selected . Thus large
may increase the chance of including the true inside the5iinsensitive margin while eliminating bad outliers. As ; becomes more and more accurate, we can reduce
to improve the resolution of the final results, which is similar to what is 1. For non-support vector changed;
done in multiscale image analysis [9].
2.5 Summary of the algorithm The proposed algorithm is summarized below.
! o+ Z.0 [ 1 3 from ! Ð+- /.h [ 1 and ini2. Randomly select tialize
; 3. Perform SVM regression using kernel. 5i anda thin-plate Find the deformed shape ; label the support vectors and non-support vectors; W3 by replacing 4. Update with better Q +- /.0the support : 1 . If vectors candidates in ! no one is updated, reduce
. Return to step Î . If no one is updated and 1. Based on the template shape , detect a landmark candidate set in the input image along the normal direction of .
Ò
This section discusses two related issues in practical application of the proposed algorithm: (a) how to select the regularization parameter (or which bounds the and ), and (b) how to deal with the case when some is empty.
Ȧ
vj
Figure 2. Support vector vector .
S , stop the algorithm and return current ; 5i as
the2desired shape.
nj
I
È!
II
The regularization parameter plays an important role in the deformation. The optimal is chosen for our problem by minimizing the ordinary cross-validation (OCV) function
ßà
+ D Ik 1 FE On 3 d ;Õá â/ã 5 )n + Õ where ;á âZã is the minimizer of Problem except that the µ th Õ landmark point is left out, i.e., ;á â/ã minimizes + 1 ED rE ä nY 3 d ; 5 n HJIKL ;ON (11) â subject to constraints (3). This is obtained using the principle of “leave-out-one” for cross validation. In 83 is alsocross-values this case, the set as an unknown which is allowed â the normal direction ^ . Similar to vary freely along to the analysis in Section 2.3, solution to this problem can be obtained by solving the quadratic programming in Problem Ý together with two additional constraints É É ¦ S â toâ be a non-support vector and In this case, are fixed 3 â affect the results. Based on the acquired È â (so ¦ ) will not and È â , we can apply (9), (6),and (5) to get and ; . WahbaI [16] has conducted a thorough study on how to estimate for the scalar function approximation splines and also extended the OCV to generalized cross validation (GCV). I The main principle can also be adapted to select optimal for our SVM-based shape deformation. 3.2 No landmark candidate
+q