等位函數法在影像切割之研究

全文

(1)國立交通大學資訊科學系碩士論文. 等位函數法在影像切割之研究 Research on Image Segmentation Using Level Set Methods. 研究生：何昌憲指導教授：薛元澤. 中華民國. 教授. 九十四. 年六月.

(2) 等位函數法在影像切割之研究 Research On Image Segmentation Using Level Set Methods. 研究生：何昌憲. Student：Chang-Xian Ho. 指導教授：薛元澤. Advisor：Yuang-Cheh Hsueh. 國立交通大學資訊科學系碩士論文. A Thesis Submitted to Institute of Computer and Information Science College of Electrical Engineering and Computer Science National Chiao Tung University in partial Fulfillment of the Requirements for the Degree of Master in. Computer and Information Science June 2005 Hsinchu, Taiwan, Republic of China. 中華民國九十四年六月.

(3) Pƒb¶Ê d~’5û˝. Nû`¤ : Œjè. çÞ : SõY. Å >¦×ç’mççÍ ( û˝F ) î=Ú. ¿b. ¡V, ÓOÚ7xíê, bP“ díTÜD@àÚ½b ÊrÖ@à5‡, d~’uø_½bí!…TÜ Ä¤ d~’AÑ dTÜ,íø_½b{æ …û˝í3bñíuû˝UàPƒb¶ (Level Set Methods) íG( (Active Contour) d~’¶ Pƒb¶uø_h˛íbMj¶, ªJ'ñqí[ý| (CÞ °v6xñqTÜR}«íÔ4 BbcÜ¸õdUàPƒb ¶íø<!…/½bí d~’¶, 1nFbÊõÒ«T,}Xƒí½æ£Ê d~’,íiÿõ. -i-.

(4) Research on Image Segmentation Using Level Set Methods. student : Chang-Xian Ho. Advisor : Dr. Yuang-Cheh Hsueh. Department (Institute) of Computer and Information Science National Chiao Tung University. ABSTRACT. In recent years, with the development of computer technology, digital image processing and its applications are more and more important. Image segmentation is an important and basic processing before other higher level applications. Image Segmentation has become one of the major subjects in image processing. The purpose of this thesis is to explore the image segmentation using Level Set Methods. Level Set Methods are novel numerical methods, they express curves or surfaces easily and have the convenient feature of dealing with PDE calculations. We implement several basic and important methods of image segmentation using level set methods. Then discuss the main problems encountered in practice and their pros. and cons. on image segmentation.. -ii-.

(5) Ð. á. ílb>áíuBíNû`¤ Œjè`¤, Ê¥sû˝FÞ®2k{“£ û˝,í2-NûD`, .céBo*Ö 25Dû˝5?‰, yéB¿… ñ}°çT0í−Ü, U…d)Jß‚êA, ã¤_,|y£íá< ?>áõðçÅ8² ¯Òkû˝,í£‡ ; °¢R+ pH ÒÆk {“,ínDû˝ ; ç!_Œ !ê ²– 6Ô NL p/U’8íõð kÅ¡— ¤ÕÔ>áÊÕø–J°Ü;í°¸£¤bí´-D2¥ |(, >áB|â=í‚£ðA, ÖVúBÖ:eKí< DÌ$Ì3í G|, UB?Êé íÞº2°ç ãø¤¹d.#Fb. SõY ãÐk Å >¦×ç’mçû˝F 2M¬Å ûý~. -iii-.

(6) CONTENTS. ABSTRACT(CHINESE) . . . . . . . . . . . . . . . . . . . . . . . . . .. i. ABSTRACT(ENGLISH) . . . . . . . . . . . . . . . . . . . . . . . . . .. ii. ACKNOWLEDGEMENT . . . . . . . . . . . . . . . . . . . . . . . . .. iii. CONTENTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. iv. LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . .. vi. 1. INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . .. 1. 1.1. Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 1. 1.2. Related Works . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 1. 1.3. Organization of Thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 2. LEVEL SET METHODS . . . . . . . . . . . . . . . . . . . . . . . .. 3. 2.1. Implicit Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 3. 2.2. Signed Distance Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 4. 2.3. Motion of Level Sets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 6. 2.3.1. External Velocity Field . . . . . . . . . . . . . . . . . . . . . . . . . .. 7. 2.3.2. Normal Direction with Mean Curvature . . . . . . . . . . . . . . . . .. 9. 2.3.3. Normal Direction with Constant Velocity Field . . . . . . . . . . . . .. 10. 2.3.4. Overall of Level Set Methods . . . . . . . . . . . . . . . . . . . . . .. 11. 2. 3. IMAGE SEGMENTATION USING LEVEL SET METHODS . . . . . . . . 14 3.1. Preview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 14. 3.2. Active Contours with Edge Stopping . . . . . . . . . . . . . . . . . . . . . . .. 15. -iv-.

(7) 3.3. Gradient Vector Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 17. 3.4. EdgeFlow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 19. 3.4.1. Useful Tools . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 20. 3.4.2. Intensity EdgeFlow . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 21. 3.4.3. Texture EdgeFlow . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 22. 3.4.4. EdgeFlow Vector . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 23. Region Based Image Segmentation . . . . . . . . . . . . . . . . . . . . . . . .. 26. 3.5.1. Description of model . . . . . . . . . . . . . . . . . . . . . . . . . . .. 26. 3.5.2. The Level Set Formulation of the Model . . . . . . . . . . . . . . . . .. 27. 3.5.3. Model of Vector-Valued Images . . . . . . . . . . . . . . . . . . . . .. 28. 3.5. 4. EXPERIMENT RESULTS . . . . . . . . . . . . . . . . . . . . . . . 30 4.1. 5. Vector Driven Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 30. 4.1.1. Curvature Term . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 32. 4.1.2. Edge Function Term . . . . . . . . . . . . . . . . . . . . . . . . . . .. 32. 4.1.3. Vector Field Term . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 35. 4.1.4. Force Balance Problem . . . . . . . . . . . . . . . . . . . . . . . . . .. 39. 4.2. Energy Balance Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 40. 4.3. Artificial Images . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 41. CONCLUSIONS AND FUTURE WORKS. . . . . . . . . . . . . . . . . 44. REFERENCES. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45. -v-.

(8) List of Figures. Fig. 2.1 parameter representation . . . . . . . . . . . . . . . . . . . . .. 4. Fig. 2.2 Implicit representation of x2 + y 2 = 1 . . . . . . . . . . . . . . . .. 4. Fig. 2.3 Motion of a curve using level set methods . . . . . . . . . . . . . .. 7. Fig. 2.4 Motion of a curve under external velocity field . . . . . . . . . . . .. 9. Fig. 2.5 Motion of a curve under curvature driven velocity . . . . . . . . . . . 10 Fig. 2.6 Motion of a curve under constant normal velocity . . . . . . . . . . . 11 Fig. 2.7 Topological changes of curves . . . . . . . . . . . . . . . . . . . 13. Fig. 3.1 Image Segmentation using level set methods . . . . . . . . . . . . . 15 Fig. 3.2 Edge function of Lenna . . . . . . . . . . . . . . . . . . . . . . 16 Fig. 3.3 Compare vector field of GAC with GVF of V-shape. . . . . . . . . . . 19 Fig. 3.4 Edge lies on direction θ from s . . . . . . . . . . . . . . . . . . . 21 Fig. 3.5 Filter responses in the Gabor filter dictionary . . . . . . . . . . . . . 23 Fig. 3.6 The procedure to generate EdgeFlow . . . . . . . . . . . . . . . . 25. Fig. 4.1 The procedure of image segmentation using level set methods . . . . . . 31 Fig. 4.2 The curve evolution under different κ weighting . . . . . . . . . . . . 32. -vi-.

(9) Fig. 4.3 Different control of generating edge function . . . . . . . . . . . . . 34 Fig. 4.4 Edge linking by threshold . . . . . . . . . . . . . . . . . . . . . 34 Fig. 4.5 Simplest segmentation by active curve . . . . . . . . . . . . . . . . 35 Fig. 4.6 Improvement by GVF . . . . . . . . . . . . . . . . . . . . . . 36 Fig. 4.7 Multi-Object curve stop at gaps . . . . . . . . . . . . . . . . . . 36 Fig. 4.8 Vector Field on the concave region . . . . . . . . . . . . . . . . . 37 Fig. 4.9 The curve evolution under different vector field . . . . . . . . . . . . 37 Fig. 4.10 U-shape Concave region . . . . . . . . . . . . . . . . . . . . . 38 Fig. 4.11 EdgeFlow of the Texture Box . . . . . . . . . . . . . . . . . . . 38 Fig. 4.12 Streamlines of EdgeFlow vector field . . . . . . . . . . . . . . . . 39 Fig. 4.13 The curve initialize outside the objects . . . . . . . . . . . . . . . 39 Fig. 4.14 Region based method initial in different location . . . . . . . . . . . 41 Fig. 4.15 Segment three objects by Gradient Vector Flow . . . . . . . . . . . . 42 Fig. 4.16 Segmentation of the 15% Gaussian noise image. . . . . . . . . . . . 42 Fig. 4.17 Segmentation of the 45% Gaussian noise image. . . . . . . . . . . . 42 Fig. 4.18 Segmentation of the 30% and 45% Gaussian noise image using region based method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43. -vii-.

(10) CHAPTER 1 INTRODUCTION 1.1 Motivation With the develop of computer and digital equipment, image processing has become one of major subjects in developing computer vision. Segmentation is an important technique used in image processing to identify objects in an image. Many different methods are used to partition an image into meaningful regions. One of these methods is to use an active contour to deform its shape and lactation according to the relative data from the image. The idea behind active contours for image segmentation is quite simple. We specify an initial guess for the contour, which is then moved by image driven forces to the boundaries of the desired objects. And recently, the level set methods are used to implement the active contours, based on various image analysis methods. Since the level set methods are novel and known by few people and its applications on the image processing are also developing, it is interesting to know how the method is and processing segmentation differ to well-known methods. Our motivation is to organize several basic and important methods of image segmentation using level set methods, and gives better understanding of their pros, cons, improvement between each other, and defects of these methods.. 1.2 Related Works There are two deformable models. The traditional one is the Parametric Active Contours (PACs), also called snakes, which trace the contours by discretizing them to finite number marker points. The PACs restricts the degree of topological adaptability, especially if the deformation involves splitting or merging of parts. In contrast, the Geometric Active Contours. -1-.

(11) (GACs) utilizes the level set methods introduced by Osher and Sethian [1][2] recently. The level set methods express the curve or surface as an implicit function, the topology changes naturally without any other dealing, and it is easy to calculate local properties. They give a smooth curve expression, and can be used in 3D models and simulation nature phenomena. There are three forces considered, the curvature, defined from the curve itself, are designed to keep the curve smooth during the deformation process, the normal direction constant force which are used to decide the expanding or shrinking of a curve and force from external felicity field, which is computed from the underlying image data, are defined to move the curve toward an object boundary or other desired features within the image.. 1.3 Organization of Thesis In chapter 2 we introduce the basic level set methods. We first introduce the basic idea of level set methods and the function of three forces. Then we discuss each force and there effect of curve evolution. In chapter 3 we introduce our main interest of using level set methods on image segmentations. We first introduce basic method from snake followed by fixed form in level set methods of Geodesic Active Contour (GAC) [3]. These two method generate small capture vector fields which point to the boundary of objects. The extension method is Gradient Vector Flow (GVF) [4], it improves the capture ability of the vectors. Different to previously methods which utilize gradients of an image, the Edgeflow, purposed in [5] is based on Gabor filter and Difference of Offset Gaussian to prediction the vector points to the boundary of objects. Finally, we introduce the region stability method [6]. In the chapter 4 we will test these methods on artificial images and to observe the noise effect on these methods.. -2-.

(12) CHAPTER 2 LEVEL SET METHODS. In this chapter we will introduce the background knowledge of level set methods introduced by Osher and Sethian [1, 2]. Level set methods are simple numerical methods used to solve various PDE problems. Recently, they have been used to many areas, including optimal design, CAD, 3D shape modelling, motion simulation, etc. We will introduce basic concepts of level sets in sections 1.1 and 1.2, then introduce evolving curve grounds on external velocity, mean curvature and normal direction velocity in section 1.3. Finally, we will discuss the advantages and disadvantage of using level set methods to implement curve evolution.. 2.1 Implicit Functions Consider a unit circle in two dimensions. We generally represent such a curve as a function x2 + y 2 = 1. A convenient way of approximating the representation is to parameterize the curve by discretizing it to a finite set of points s, s ∈ [0, 1] , such that each point x is on the curve, x = ~v (s) = (x(s), y(s)). The curve is denoted as C(s) : [0, 1] → R2 . Alternatively, the idea of level set methods is to represent the interface in Rn in one higherdimensional function. For example, the zero isocontour of function φ(x, y) = x2 + y 2 − 1, i.e. φ(x, y) = 0, also expresses the same unit circle as above. Note that φ is a three-dimensional function, and the unit circle is a two-dimensional curve. For an implicit function φ(x), here x = (x1 , . . . , xn ) ∈ Rn , define the interior region an open sets Ω− = { x | φ(x) < 0 }, and the exterior region Ω+ = { x | φ(x) > 0 }. The contour is defined by ∂Ω = { x | φ(x) = 0 }.. -3-.

(13) Fig. 2.1. Fig. 2.2. parameter representation. Implicit representation of x2 + y 2 = 1. 2.2 Signed Distance Function Given any interface Γ in Rn , a signed distance function is used to generate the underly implicit function in a different way. A distance function d(x), where x = (x1 , . . . , xn ), is defined as d(x) = min(|x − xI |) for all. xI ∈ ∂Ω,. (2.1). implying that d(x) = 0 if x ∈ Γ. A signed distance function is an implicit function φ with |φ(x)| = d(x). The φ(x) is defined. -4-.

(14) as φ(x) = −d(x) for. x ∈ Ω−. φ(x) = d(x). for. x ∈ Ω+. φ(x) = 0. for. x ∈ ∂Ω = Γ.. (2.2). In our example, we replace the implicit function φ(x, y) = x2 + y 2 − 1 with the signed p distance function φ(x, y) = x2 + y 2 − 1 in order to represent the unit circle ∂Ω = {x | |x| = 1}. Between the implicit function and signed distance function, there are several common properties and advantages of geometry operations. • Region definition The signed distance function gives the same boundary ∂Ω, interior region Ω− , and exterior region Ω+ , that φ(x, y) = x2 + y 2 − 1 does. • Local properties of the interface – Gradient and Normal Define the gradient as ∇φ =. ∂φ ∂φ , ∂x ∂y. .. (2.3). The gradient ∇φ is perpendicular to the isocontours, and points to the direction of increasing φ. This indicates that if x is a point at one isocontour, ∇φ(x) has the ~ at x. Thus the unit normal is same direction as normal N ~ = ∇φ . N |∇φ|. (2.4). Furthermore, the signed distance function has a new property that is |∇φ| = 1 for all x 6= 0.. -5-. (2.5).

(15) – Mean Curvature ~, The mean curvature of the interface is defined as the divergence of the normal N that is κ=∇·. ∇φ |∇φ|. ,. (2.6). so that κ > 0 for convex regions, κ < 0 for concave regions, and κ = 0 for a plane. • Unity and Binary operations For an implicit function φ(x), the complement of φ can be defined by φc (x) = −φ(x). If there are two implicit functions φ1 and φ2 , the union of the interior regions of φ1 and φ2 can be obtained by φ(x) = min(φ1 (x), φ2 (x)). Similarly, intersection of the interior regions of φ1 and φ2 can be obtained by φ(x) = max(φ1 (x), φ2 (x)). It’s easy to generate other operations of the interfaces like these do while using implicit functions. Fig. 2.3 is an example of using level set methods to express curve evolution. Fig. 2.3a is the p distance function φ(x) = x2 + y 2 − 1, Fig. 2.3b shows the contours of φ for each level set, we can see that the zero level set is a unit circle. If we suppose the unit circle evolves under the normal velocity of constant unit speed, after a time step 4t the φ(x) evolves as shown in Fig. 2.3c. Fig. 2.3d shows the contours of Fig. 2.3c, we have seen that now the zero level set is a circle with radius two. In fact, we can also let the φ be positive inside the circle and negative outside the circle, that will turn the cone upside down, and still gives the same zero level set of unit circle. At this p situation, φ(x) = 1 − x2 + y 2 .. 2.3 Motion of Level Sets The zero level set is generated to the corresponding curve that we are interested in. Then, evolving the φ(x, y) gives the curve evolution by tracking the zero level set at each time step.. -6-.

(16) (a) φ(x) at time t. (b) φ(x)=0 at time t. (c) φ(x) at time t + 4t. (d) φ(x)=0 at time t + 4t. Fig. 2.3. Motion of a curve using level set methods. The level set methods can also be used in higher dimensional hyper-surface, for example the three-dimensional surface and thus the φ is a four-dimensional function. Since our main goal is to use level set methods to implement curve evolution and image segmentation, we will only discuss two-dimensional curves and three-dimensional signed distance functions.. 2.3.1 External Velocity Field Let φ(x, y, t) denote φ(x, y) at time t and Γ(t) be the zero level set of φ(x, y, t). Suppose the external velocity of each point at Γ is given as V~ (x, y). Generating such external velocity filed V~ = (u, v), we wish to move Γ with this velocity. This gives the equation Γt = V~ (x, y).. (2.7). Since Γ(t) is the evolving zero level set of φ(x, y, t), this also gives φ(Γ(t), t) = 0. Taking. -7-.

(17) the derivative to t and applying chain rule gives φt + V~ · ∇φ = 0,. (2.8). this is the basic equation that will be used to introduce other kinds of velocity fields. Since we are working on the images, we can generate cartesian grids with size of images, then the numerical methods are used to solve the φ forward in time. A rather simple first-order accurate method is the forward Euler method given by φn+1 − φn ~ n + V · ∇φn = 0, 4t. (2.9). where 4t is the discretized time step, and tn+1 = tn + 4t, φn+1 =φ(x, y, tn+1 ). To write equation (2.9) in expanded form as φn+1 − φn + un φnx + v n φny = 0, 4t. (2.10). where (u, v) is the external velocity. We look at one single point (xi , yi ) with velocity (ui , vi ) as example, if ui > 0 the value of φ(xi , yi ) moves from left to right, and the method of characteristics tells us to look to the left of (xi , yi ) to determine what value of φ will land on the point (xi , yi ). Similarly, if ui < 0, we look to the right of (xi , yi ) to determine what value of φ will be. Clearly, backward difference −x +x Di,j φ should be used to approximate φx when ui > 0 and forward difference Di,j φ should be. used to approximate φx when ui < 0, where φi,j − φi−1,j 4t φi+1,j − φi,j = . 4t. −x Di,j = +x Di,j. (2.11) (2.12). The same is in approximating φny . This method of choosing an approximation to the spatial derivatives based on the sign of u and v is known as upwinding. There are many other methods can get second order or third order accurate in approximating φnx and φny , for example, Hamilton-Jacobin ENO, Hamilton-Jacobin WENO, and TVD RK (Total Variation Diminishing Runge-Kutta) method [7]. However, we use simple upwind method. -8-.

(18) and forward Euler in this thesis. The reason is that we are interesting in segmentation using curve evolution not in numerical accurate of simulating curve evolution. Fig. 2.4 shows the circle moving under the external velocity field V~ = (1, 0), this velocity field pushes the circle to move from left to right.. Fig. 2.4. Motion of a curve under external velocity field. 2.3.2 Normal Direction with Mean Curvature In the last subsection we discussed the motion of an interface in an external velocity filed. Now we consider the motion for a self-generated velocity. That is, the interface moves in normal direction with a velocity proportional to its curvature, which is generated from the φ itself. If ~ + Vt T~ , where N ~ is the normal vector and T~ is tangent vector to the we decompose V~ = Vn N interface, the level set equation (2.8) is equivalent to ~ + Vt T~ ) · ∇φ = 0, φt + (Vn N. (2.13). ~ and ∇φ point in the same direction, T~ · ∇φ = 0. Then equation (2.13) becomes since N ~ · ∇φ = 0. φt + Vn N. (2.14). Here Vn = −bκ, where b > 0 and using (2.4), finally we have φt − bκ |∇φ| = 0.. (2.15). The mean curvature has the following effects. At the convex region, the curvature κ > 0 such that the velocity in the normal direction Vn < 0. This results in shrinking curve to the direction opposite to normal and trying to smooth curve. Similarly, at the concave region the. -9-.

(19) Fig. 2.5. Motion of a curve under curvature driven velocity. curvature κ < 0 and Vn > 0, the curve expands as to smooth curve. This motivation lets any curve evolving into a circle and shrinking to disappear. Fig. 2.5 shows a star-shaped contour moving under curvature driven velocity. The tips, that are convex, move inward while the gaps, that are concave, move outward.. 2.3.3 Normal Direction with Constant Velocity Field In the last subsection the normal directional velocity of the interface is given by the mean curvature of the interface that determined by local geometric information. Now, we discuss the ~ . Just directly normal directional velocity of interface given by a constant velocity field V~ = aN replacing −bκ by a in the (2.15) yields the corresponding level set equation φt + a |∇φ| = 0.. (2.16). When a > 0 the interface moves in the normal direction, and when a < 0 the interface moves in the opposite direction. When a = 0, the interface stays still. Since φ is a signed distance function we have |∇φ| = 1, this reduces equation to φt = −a. By Forward Euler, φn+1 = φn − a4t, this equation gives a hint of the behavior of a curve.. -10-.

(20) (a) a > 0. (b) a < 0. Fig. 2.6. Motion of a curve under constant normal velocity. Suppose a > 0, if at some points x0 and x1 , which φ(x0 ) = 0 and φ(x1 ) = a4t, after a time step 4t, the φ(x0 ) = −a4t and φ(x1 ) = 0. This means that the zero level set moves to the normal direction with distance a4t. Similarly, if a < 0 the zero level set moves to the opposite direction with distance a4t. If we let our distance function φ be positive outside the circle and negative inside the circle, then the whole curve will expand outward for a > 0 and shrink inward for a < 0. Fig. 2.6a is an example of a star-shaped moving under the constant normal velocity a = 1 > 0, and Fig. 2.6b is under the constant normal velocity a = −1 < 0. Compare to the motion under the curvature driven velocity, where only the curves of the convex regions are shrinking, the motions of convex regions and concave regions under the constant normal velocity are either both shrinking or both expanding.. 2.3.4 Overall of Level Set Methods We have discussed the concept of using level set methods to express curves and their motions under three kinds of velocities. Combine these three kinds of velocities to get, φt − bκ |∇φ| + a |∇φ| + V~ · ∇φ = 0.. -11-. (2.17).

(21) Let F be the function of normal directional constant velocity of grid points. Weight each kind of velocity to yield the level set equation φt = ακ |∇φ| + βF |∇φ| − γ V~ · ∇φ.. (2.18). By Forward Euler method, we have φn+1 = φn + 4t(ακ |∇φ| + βF |∇φ| − γ V~ · ∇φ).. (2.19). The Courant-Friedreichs-Lewy condition (CFL condition) asserts that if want to have a stability approximation solution of (2.19), the numerical waves should propagate at least as fast as the physical waves. In other words, (ακ |∇φ|+βF |∇φ|−γ V~ ·∇φ) < 4x/4t over all grid points. Since we work with images, the grid size is set to the pixel size. Thus 4x = 1 pixel. Our time step must be chosen to satisfy 4t < 1/ max(ακ |∇φ| + βF |∇φ| − γ V~ · ∇φ). The first advantage of using level set methods is that we don’t need to care about the merging or separating of curves. The topological changes of curve evolution is not problematic as traditional curve tracking methods, which have to remove or add tracking points and reconnect if two curves are merging or separating. The level set methods automatically handle topological changes. Fig. 2.7a shows two single circles that are expanding to connect each other to generate a single curve. Fig. 2.7b shows two circles, the outside circle is shrinking while the inside expanding, then merge to a horn shape. The second advantage is that it is easy to calculate the local properties of the curve, such as normal and curvature discussed in Section 2.2. Third is the generality for arbitrary dimensions and no self-intersection problem which results in swallowtail in PACs. And the most disadvantage is its computational expensiveness.. -12-.

(22) (a) Merge of two expanding circles. (b) Merge of two circles which outside is shrinking and inside is expanding. Fig. 2.7. Topological changes of curves. -13-.

(23) CHAPTER 3 IMAGE SEGMENTATION USING LEVEL SET METHODS. In this chapter we will introduce our survey of image segmentation using level set methods. Section 3.1 gives a preview of classification and skeleton of image segmentation using level set methods. Section 3.2, 3.3 and 3.4 provide the detail of Geodesic Active Contour, Gradient Vector Flow, and EdgeFlow, respectively. Rather than previous edge-based methods, Section 3.5 introduce the region-based method to evolve contours via level set methods.. 3.1 Preview Image segmentation is a basic preprocessing of many image analysis applications. Most methods process image segmentation based on various attributes such as texture, color, and gray level intensity. Most previous approaches to image segmentation can be classified as follows: • Filtering-based methods : detect edges followed by edge linking. • Region growing and merging. • Global optimization based on energy functions and bayesian criteria. • Graph partitioning and clustering. • Curve evolution and active contour models. Being different to the other ways which give a fragment result needed further processing, curve evolution methods usually result in closed contours and give a whole complete object segments. The two different methods to represent and implement curve evolutions are parametric active contours (PACs) and geometric active contours (GACs). The major problem of using active contours to approach image segmentation is how to drive the curve to move to the boundaries of objects. In this thesis we will discuss several well-known. -14-.

(24) Image Segmentation Using Level Set Methods. Edge-based. Active Contour. Region-based. EdgeFlow. Mumford-Shah Model. Gradient Vector Flow (GVF) Fig. 3.1. Image Segmentation using level set methods. level set methods of image segmentation using active contours.. 3.2 Active Contours with Edge Stopping Caselles et al. [8] and Malladi et al. [9] are among the first to use level set methods to extract objects from an image in their parallel works. They suppose the curve is either expanding or shrinking and stop evolving at the object boundary. This gives the idea to design a constant velocity field in the normal direction such that the velocities are positive or negative when the curve is at the smooth region and zero while at the object boundary. Suppose C(s) : [0, 1] → R2 is a parametrization of a 2-D closed curve. The edge stopping is designed with the following edge function g(x, y) =. 1 , 1 + |∇Gσ (x, y) ∗ I(x, y)|. (3.1). where Gσ is a Gaussian function with standard deviation σ. In equation (3.1), it is obvious that at smooth regions we have |∇Gσ ∗ I| → 0 implies g → 1 and at edges we have |∇Gσ ∗ I| → ∞. -15-.

(25) implies g → 0. In fact, |∇Gσ ∗ I| as known is a simple method to find edges in an image. Another edge function which falls to zero faster on edges can be defined as: g(x, y) = e−|∇Gσ (x,y)∗I(x,y)|. Fig. 3.2. (3.2). g(x, y) = 1/(1 + |∇Gσ (x, y) ∗ I(x, y)|) of Lenna, normalize to 0 ∼ 255, where black means small value while white means large value.. Modify the curve evolution PDE as: ~, Ct = g(bκ + F )N. (3.3). the corresponding level set equation is φt = g(κ + F ) |∇φ| .. (3.4). The curve will expand if F is always positive and shrink if F is always negative. With the effects of g, the curve will slow down when progressing near to the edges, then stops while reaching the edges. The problem of this method is that in the numerical calculation if the curve propagates beyond the desired boundary, it continues propagating, there is no other mechanism to push the curve back to the object boundary.. -16-.

(26) To overcome this problem, Cassela et al. [3] introduce Geodesic Active Contour (GAC) method , which uses minimization of the energy function Z M in. 1. g(C) |C 0 (s)| ds. (3.5). 0. and finds the gradient descend equation ~ − (∇g · N ~ )N ~. Ct = g(κ + F )N. (3.6). The corresponding level set equation is φt = g(κ + F ) |∇φ| + ∇g · ∇φ.. (3.7). The term ∇g generates an external velocity vector field V~ = −∇g. The corresponding vector is almost zero at the smooth region, and points to the nearest edge at the neighborhoods of this edge, see Fig. 3.3a. As the contour moves beyond the edge, the vector pushes the contour back to the edge.. 3.3 Gradient Vector Flow Since the ∇g is almost zero at smooth region, the initial contour must be close to the true boundary or else it will evolve to ill result. In the papers [4] [10], Xu and Prince use an energy function to generate vector flow filed which is propagated from certain image generating vector field. Define f 1 (x, y) = |∇I(x, y)|2 , (3.8) 2. 2. f (x, y) = |∇Gσ (x, y) ∗ I(x, y)| . The field ∇f has vectors pointing toward the edges, but as mentioned above this kind of field has a narrow capture range. Xu an Prince generate vector filed ~v(x, y) = (u(x, y), v(x, y)),. -17-.

(27) called the gradient vector flow (GVF) field, from minimizing the energy functional ZZ ε=. µ(ux 2 + uy 2 + vx 2 + vy 2 ) + |∇f |2 |~v − ∇f |2 dx dy,. (3.9). where f is chosen as f 1 or f 2 and µ is a parameter governing the tradeoff between the first and the second terms. We see that when |∇f | is small, the vector is dominated by partial derivatives of the vector field, when |∇f | is large, the second term dominates the integrand, and is minimized by setting ~v = ∇f . This says that the minimum solution will hold the original vector at neighborhoods of edges and propagate the vector for the smooth region next to the edges. The calculus of variations gives the minimal solution of energy function in Euler equation forms. Given ZZ F (x, y, u, ux , uy )dx dy,. E[u] =. (3.10). D. the solution minimizing E[u] satisfies ∂F ∂ ∂F ∂ ∂F = + . ∂u ∂x ∂ux ∂y ∂uy. (3.11). Solving equation (3.9), the solution can be found from µ∇2 u − (u − fx )(fx 2 + fy 2 ) = 0 (3.12) 2. 2. 2. µ∇ v − (v − fy )(fx + fy ) = 0. Now, we turn to solve (3.12) by treating u and v as functions of time ut (x, y, t) = µ∇2 u(x, y, t) − (u(x, y, t) − fx (x, y)) · (fx (x, y)2 + fy (x, y)2 ) (3.13) 2. 2. 2. vt (x, y, t) = µ∇ v(x, y, t) − (v(x, y, t) − fy (x, y)) · (fx (x, y) + fy (x, y) ) where u(x, y, 0) and v(x, y, 0) are initialized as ∇f . In practice, the more iterations to solve (u, v) the more range vectors diffuse and the larger capture ability the result has. C. Xu and J. L. Prince use the GVF in the parametric curve model Ct (s, t) = αC 00 (s, t) + βC 0000 (s, t) + ~v.. -18-. (3.14).

(28) (a) −∇g. Fig. 3.3. (b) GVF v. Compare vector field of GAC with GVF of V-shape.. In this thesis we replace f with −∇g, in last subsection, to generate GVF, and use GVF as external velocity field. This results in modifying equation (3.7) as φt = g(κ + F ) |∇φ| − ~v · ∇φ.. (3.15). Since v has large capture rang, the curve evolution can be driven by field rather than constant speed F . The curve will shrink at somewhere and expand at others.. 3.4 EdgeFlow Previous active contour methods aim to identify an object by utilizing the local discontinuities such as edges. The weak connection between contours and images is only with edge function or vector flow generated from edge function. Obviously, it’s bad on nature images or texture images. As we will discuss in our experiment result, the segmentation is based on the founded edges. Alternatively, Wei et al. [11] generate the EdgeFlow from the intensity information of the pixel and its neighborhoods, then Sumengen et al. [5] use EdgeFlow in active contours. Decide the number of directions n and θ = {πi/n | i ∈ 0, . . . , n − 1}. Consider a location s in the image, define a tuple measurement F (s, θ) = [E(s, θ), P (s, θ), P (s, θ + π)], where. -19-.

(29) E(s, θ). edge energy at s along the orientation θ.. P (s, θ). probability of finding an image boundary in the direction θ from s.. P (s, θ + π). probability of finding an image boundary in the direction θ + π from s.. Sumengen et al. generate these measurements based on intensity and texture properties. For vector-valued image the intensity measurement can be calculated at each layer, RGB color image as example, we will have FR (s, θ), FG (s, θ), FB (s, θ) and FT EX (s, θ) for each direction θ from each location s. Finally, combining each orientation for each location s generates final EdgeFlow.. 3.4.1 Useful Tools We introduce several useful tools used to generate F (s, θ). The two-dimensional Gaussian function is defined as Gσ (x, y) = √. 1 2πσ. −(x2 + y 2 ) 2σ 2 e .. (3.16). The first derivative of Gaussian (GD) along the x-axis is GDσ (x, y) =. ∂Gσ (x, y) x = − 2 Gσ (x, y) ∂x σ. (3.17). and the difference of offset Gaussian (DOOG) along the x-axis with distance is defined as DOOGσ (x, y) = Gσ (x, y) − Gσ (x + d, y),. (3.18). where d is the offset between centers of two Gaussian kernels. To generate the derivative of Gaussian and difference of offset Gaussian along some orientation θ, we simply rotate the coordinates with −θ and calculate GD and DOOG. From linear algebra, we know that the new coordinates are x0 = x cos θ + y sin θ (3.19) 0. y = −x sin θ + y cos θ. -20-.

(30) thus GDσ,θ (x, y) = GDσ (x0 , y 0 ) (3.20) 0. 0. DOOGσ,θ (x, y) = DOOGσ (x , y ).. 3.4.2 Intensity EdgeFlow The GDs are generally used to determine if s lies near the edge and the edge lies in direction θ from s, i.e. the edge direction is around perpendicular to θ. If so, |GDσ,θ (x, y) ∗ I(x, y)| will have a large value. In fact, this is one dimensional derivative form different to |∇Gσ (x, y) ∗ I(x, y)|.. s. Fig. 3.4. •. θ. Edge lies on direction θ from s. For two orientation θ1 and θ2 , if |GDσ,θ1 (x, y) ∗ I(x, y)| > |GDσ,θ2 (x, y) ∗ I(x, y)|, then the edge may lies from s in direction θ1 than θ2 . The function E(s, θ) will be used to express the strength of vectors. Thus, we define edge energy as E(s, θ) = |GDσ,θ (x, y) ∗ I(x, y)| .. (3.21). Another prediction method utilities DOOGs, the simple idea is that if there is an edge lies within the distance d in direction θ from s, the intensity different between (x+d cos θ, y+d sin θ) and (x, y) should be large. The error in prediction is defined to be Error(s, θ) = |Iσ (x + d cos θ, y + d sin θ) − Iσ (x, y)| = |DOOGσ,θ (x, y) ∗ I(x, y)| .. -21-. (3.22).

(31) Therefore, the probabilities of EdgeFlow direction is P (s, θ) =. Error(s, θ) . Error(s, θ) + Error(s, θ + π). (3.23). Note that σ is the only parameter needed to control.. 3.4.3 Texture EdgeFlow The texture features are extracted based on a Gabor wavelet decomposition scheme [12]. A two dimensional Gabor function g(x, y) and its Fourier transform G(u, v) are 2   2 − 1 x + y + 2πjW x 2 2 1 2 σ σ x y g(x, y) = e , 2πσx σy 1 − G(u, v) = e 2. . (u − W )2 v 2 + 2 σu2 σv ,. (3.24). (3.25). where σu = 1/2πσx and σv = 1/2πσy . By dilations and rotations of g(x, y) through gmn (x, y) = a−m G(x0 , y 0 ), a > 1, m, n ∈ Z x0 = a−m (x cos ϕ + y sin ϕ). (3.26). y 0 = a−m (−x sin ϕ + y cos ϕ). this generates a class of function decompose frequency domain. Give Ul and Uh to be the lower and upper center frequencies of interest. Let K be the number of orientations and S the number of scales. Then design 1 Uh − S−1 ) Ul (a − 1)Uh √ σu = (a + 1) 2 ln 2 2 − 1 π σu (2 ln 2)2 σu2 2 σv = tan Uh − 2 ln 2 ln 2 − 2k Uh Uh2. a=(. (3.27). using (3.26) with W = Uh , m = 0, 1 . . . , S − 1, and k = 0, 1 . . . , K − 1. The total number of decomposed frequencies are N = S · K, see Fig. 3.5 as an example.. -22-.

(32) Fig. 3.5. Filter responses in the Gabor filter dictionary. The filter parameters used are Uh = 0.35, Ul = 0.1, K = 6 and S = 4.. Ma et al. define Ul = 1/4σ and Uh = 0.45, as in the intensity edges, the σ is the only parameter needed to control. Applying filters on the image gives the Gabor filtered images Oi (x, y) = gi (x, y) ∗ I(x, y) = mi (x, y) e[ϕi (x, y)]. (3.28). where 1 ≤ i ≤ N . For each location (x, y), now we have a texture feature vector Ψ(x, y) = [m1 (x, y), m2 (x, y), . . . , mN (x, y)]. (3.29). The texture edge energy is defined as E(s, θ) =. X. |GDσ,θ (x, y) ∗ mi (x, y)|. (3.30). |DOOGσ,θ (x, y) ∗ mi (x, y)|. (3.31). 1≤i≤N. and the error in prediction is defined as Error(s, θ) =. X 1≤i≤N. The probabilities P (s, θ) can be estimated using (3.23).. 3.4.4 EdgeFlow Vector For each location s, we had generated intensity and texture edge energy and probabilities in different orientation. The edge energy and probabilities can be combined together to form a. -23-.

(33) single EdgeFlow field. Consider E(s, θ) =. X. P (s, θ) =. X. Ea (s, θ) · ω(a). a∈A. (3.32) X. Pa (s, θ) · ω(a),. a∈A. ω(a) = 1. a∈A. where A = {intensity/color, texture}, the ω is the weighting coefficient with different image attribute a ∈ A. For example, EdgeFlow of the RGB color model image could be obtained with A = {red, green, blue, texture} and maybe let ω(texture) = 0.4 and ω(red) = ω(green) = ω(blue) = 0.2. After combine various edge energy and probabilities, there are only one edge energy and one probability in each orientation at location s. The following work is to deicide flow direction at each lactation s form {[E(s, θ), P (s, θ), P (s, θ + π)] | 0 ≤ θ < π} Define the continuous range of flow direction which maximize the sum of probabilities in a corresponding half plane (. ) X. Θ(s) = arg max θ. P (s, θ0 ) .. (3.33). θ≤θ0 <θ+π. The EdgeFlow vector is then defined by F~ (s) =. X. E(s, θ) · ejθ .. (3.34). Θ(s)≤θ<Θ(s)+π. The F~ (s) is a complex number with magnitude representing edge energy and angle representing the flow direction. In our implementation the direct flow vector is   X X V~ (s) =  E(s, θ) cos θ, E(s, θ) sin θ . Θ(s)≤θ<Θ(s)+π. (3.35). Θ(s)≤θ<Θ(s)+π. After this vector field is generated, the original method is followed by propagating vectors. The propagation ends and edges are made by checking the opposite direction vectors. If x or y component of two connect vectors change signs, the edges are marked. The boundaries are found by linking the edges. Segmentation is further processed with region merging. Sumengen et al. [5] use EdgeFlow as external vector field in front propagation. They utilize EdgeFlow to produce the normal direction constant velocity F as the edge function. The edge. -24-.

(34) Image. Intensity. Texture Gabor filter. channalG Eb (s, θi ). channalR Er (s, θi ) θ1. θ1. θ1. θ1 θ2. θ2 θn. g2 ∗ I. g1 ∗ I. channalB Eg (s, θi ). θ2. θ1. θ2. ··· .. .. gN ∗ I θ1 θ2. θ2 θn. θn. Etex (s, θ1 ). Etex (s, θ1 ). ···. Etex (s, θn ). .. . ω(r). θn. ω(g). ω(b). E(s, θ1 ). ω(t). E(s, θ2 ). ···. E(s, θn ). E(s, θ). Fig. 3.6. The procedure to generate EdgeFlow. function F is derived as follows F =. 1

(35)

(36) .

(37)

(38) 1 +

(39) V~

(40). (3.36). This is the reverse procedure of gradient vector flow method, which generates edge function first and then external velocity. We summarize the procedure as follows. -25-.

(41) 3.5 Region Based Image Segmentation Most methods of pushing the contour to the edge utilize vector flow field and stopping edge function, as we had surveyed prior, with level set active contour model (2.18). Another method is to generate level set functions based on minimization of an energy based-segmentation. Chan et al. [6] proposed a model using Mumford-Shah functional for segmentation [13].. 3.5.1 Description of model ¯ → R is the given image and C is the curve in Ω. The Mumford-Shah functional Suppose u0 : Ω for segmentation is F M S (u, C) =µ · Length(C) Z + λ |u0 (x, y) − u(x, y)|2 dx dy. (3.37). Ω. Z. |∇u(x, y)|2 dx dy,. + Ω\C. where µ and λ are positive parameters. The solution of image u is obtained by minimizing this functional and is formed by smooth regions Ri with sharp boundaries, denoted here by C. Mumford and Shah reduced F M S with restriction to piecewise constant functions u = ci where ci is constant on each connected component Ri of Ω\C. The model proposed by Chan et al. is a particular case of the Mumford and Shah’s model, which the Ω\C is decomposed according to only two piecewise components R1 = inside(C) and R2 = outside(C). The idea of Chan et al. is that assume the image u0 is formed by two regions of piecewiseconstant intensities, generally one object and one background. Then consider the fitting energy function Z. 2. Z. |u0 (x, y) − c1 | dx dy +. F1 (C) + F2 (C) = inside(C). |u0 (x, y) − c2 |2 dx dy, (3.38). outside(C). where c1 = average(inside(C)) and c2 = average(outside(C)). If C is fitting to the boundary. -26-.

(42) of the object, the F1 (C) + F2 (C) ≈ 0. Therefor in active contour model the curve C will minimize this function, and we add regularizing terms the length of C and the area inside C. This energy function is F (C, c1 , c2 ) =µ · Length(C) + ν · Area(inside(C)) Z + λ1 |u0 (x, y) − c1 |2 dx dy. (3.39). inside(C). Z. |u0 (x, y) − c2 |2 dx dy.. + λ2 outside(C). 3.5.2 The Level Set Formulation of the Model Represent the curve C with level set methods, here we inverse the φ to φ > 0 inside the curve and φ < 0 outside the curve. The function becomes: F (φ, c1 , c2 ) =µ · Length{φ = 0} + ν · Area{φ > 0} Z + λ1 |u0 (x, y) − c1 |2 dx dy. (3.40). φ>0. Z. |u0 (x, y) − c2 |2 dx dy.. + λ2 φ<0. Note integral terms in the function are over inside and outside regions of the curve respectively. Besides, after evolved with error of accuracy the zero level set lying between φ > 0 and φ < 0 may not exactly zero. In order to set the integral over the entire image domain Ω and measure the length of the curve we utilize the following two special functions: • Heaviside function  1 if x ≥ 0, H(x) = 0 if x ≤ 0.. (3.41). • Dirac measure δ δ(x) =. -27-. d H(x). dx. (3.42).

(43) Express the function F in Z. Z δ(φ) |∇φ| dx dy + ν. F (φ, c1 , c2 ) =µ Ω. H(φ)dx dy Ω. Z + λ1. |u0 (x, y) − c1 |2 H(φ)dx dy. (3.43). Ω. Z + λ2. |u0 (x, y) − c2 |2 (1 − H(φ))dx dy.. Ω. Minimizing F (φ, c1 , c2 ) using a gradient descent method, yields the associate Euler equation for φ ∇φ 2 2 φt = δ(φ) µ · div − ν − λ1 (u0 − c1 ) + λ2 (u0 − c2 ) . |∇φ|. (3.44). The additional length term is a regularizing term and has a scaling role. If µ is large, only lager objects are detected, while for small µ, much and smaller objects are also detected. Chan et al. set ν = 0 and λ1 = λ2 = 1. This is different to general level set active contour model. Actually, the term −ν − λ1 (u0 − c1 )2 + λ2 (u0 − c2 )2 can be seen as the constant speed in the normal direction but is dependent on the position and time.. 3.5.3 Model of Vector-Valued Images In the continuous paper of Chan et al. [14], they extended the method to vector-valued images. The model lets the curve minimize the energy function according to each channel u0,i of an image, with i = 1, . . . , N channels. + − − − Let c+ = (c+ 1 , . . . , cN ) and c = (c1 , . . . , cN ) be the vectors of the averages of inside and. outside curve on each channel respectively. The extension of model is F (φ, c+ , c− ) =µ · Length(C) Z N

(44) 2 1 X +

(45)

(46)

(47) dx dy + λi u0,i (x, y) − c+ i inside(C) N i=1 Z N

(48) 2 1 X −

(49)

(50)

(51) dx dy. + λi u0,i (x, y) − c− i N outside(C) i=1. -28-. (3.45).

(52) Rewriting it in level set form, we have −. +. Z δ(φ) |∇φ| dx dy. F (φ, c , c ) =µ Ω. Z + Ω. Z + Ω. N

(53) 2 1 X +

(54)

(55)

(56) H(φ)dx dy λi u0,i (x, y) − c+ i N i=1. (3.46). N

(57) 2 1 X −

(58)

(59)

(60) H(1 − φ)dx dy λi u0,i (x, y) − c− i N i=1. and the Euler equation for φ " φt = δ(φ) µ · div. . ∇φ |∇φ|. . # N N X 1 1 X + 2 2 − λ (u0,i − c+ λ− (u0,i − c− . i ) + i ) N i=1 i N i=1 i. (3.47). − The µ is regularizing term as above, λ+ i and λi are used to set different weighing according. each layer.. -29-.

(61) CHAPTER 4 EXPERIMENT RESULTS. In this chapter we will discuss the experiment results using the level set segmentation methods. First, we will discuss the effects of parameters when using level set methods on image segmentation. Then some artificial images will be the sample of image segmentation in order to compare the noise effects of these methods. In our experiment we can organize these methods from the view point of how they drive the curves. The first part is driven by generated vector field with edge function, this part including Geodesic Active Contour, Gradient Vector Flow and Edge Flow. The vector field V~ and edge function F are calculated before the curve evolution and fixed at all time. Another part is driven by energy balance which the curve evolution is according to the region stability which changes with time to time. We will discuss the improvement progress of these methods and observe the practical situation in our experiment.. 4.1 Vector Driven Methods The general level set equation of curve evolution is of the form φt = ακ |∇φ| + βF |∇φ| − γ V~ · ∇φ. We process the image segmentation with procedure as in Fig. 4.1 to generate the curve evolution segmentation. With the several methods we have introduced, each term can be summarized as • κ based on the curve itself, • F the edge function used to decide whether the curve is still expanding or shrinking or stops at the edges, • V~ the velocity field pushes the curve to move to the edges.. -30-.

(62) Image. Caculate Vector Flow V~. Initialization φ0. Edge Fun. F Caculate V~ and F dependent speed. Caculate curvature dependent speed. Update φt , t = t0 + n4t. !aa !! a ! ! Ct : {φt = 0}aaa No ! aa ! converges? aa !! aa!!!. Yes Stop Evoluation Fig. 4.1. The procedure of image segmentation using level set methods. The total amount of these three forces decide the curve evolution. Besides, α, β, and γ are weights to each term and σ is used in calculating Gaussian filtered image. In the experiment we manual set γ = 1 and α < 0.3 for each method. How to set a value for β is the most critical problem in curve evolution for segmentation. Too small β will make curve move slowly and may stick if there is no external velocity vector. Too large β will let curve move across the weak edges even there have external velocity pushing curve back to edges, but it is not enough to balance the force given by F . We first test the effect of each velocity as follows.. -31-.

(63) 4.1.1 Curvature Term In fact, the κ is set to −κ as we discussed in section 2.3.3. The κ keeps the curve evolving smoothly since it can suppress the high curvature location. In the segmentation, the κ also effects the sensitives of curve when moving pass the noise. Fig. 4.2 shows an example when a curve is moving across the noise, In this example, the size of each dot is about 10 pixels and β = 0.3. Fig. 4.2a is curve evolving with α = 0.1, Fig. 4.2b is shows the curve evolved with α = 0.5. The difference is that when the curve meets the noise with small α value, it will be stopped according the relative F ≈ 0, although κ velocity should make small expansion of the curve it is not large enough to overcome the influence of V~ which pushes the curve back to the boundary of dots, and with large α value, the κ term of velocity will be able to overcome the influence of V~ near the boundary.. (a) α = 0.1. (b) α = 0.5. Fig. 4.2. The curve evolution under different κ weighting. 4.1.2 Edge Function Term Next, we experiment the Edge Function term F = g, we find that this term controls the segmentation results in methods GAC and GVF. There are two variables that generate edge function and affect the curve evolution. First, the σ of the Gaussian smooth function decides the pixel rang of an edge to be found. Fig. 4.3b and Fig. 4.3c show σ = 0.75 and σ = 2, respectively. With small σ values the Gaussian smooth function is sensitive to the noises and the resultant. -32-.

(64) curve is hard to evolve because at most ranges F ≈ 0. We try to increase the α and β weights of κ and F , respectively, to overcome this situation in high noise images. We find that although this makes curve to easily pass through the noise in theory, but in practice, it does not work. Especially the new problem occurs in the hight curvature corner such as branch corner of a vessel, the high α value will let curve to expand and leave the corner. Use large σ value can blur the noise but it will also extend the small F area around the noises, in this situation, some weak but border noise will result extension area of F ≈ 0. We can see this at top left branch in Fig. 4.3c, there is a blurred edge across the vessel, when curve enter the narrow vessel it will stop. We use the following fixed edge function to nonlinearly enhance the F values. F =g=. 1 . 1 + a ∗ |∇Gσ ∗ I|. (4.1). Observe that F will have a larger value if 0 < a < 1. This means the F value of noise will not approximate zero, but this also increases the F value of real edges. We can use a threshold d to cutoff the F , if F < d we set F value equal to zero. And if a > 0 it can used to reduce the F in the weak edges in noiseless images, which have smooth areas but the edge between objects are not strong. We can chose a according to the original edge function to look if an image is filled with noise or filled with smooth regions with weak edges between them. The edge linking is another problem in F , if there are small gaps in the edge the curve will expand out from this gap. We use a threshold as above to link the edge, this method also set the F value of edges absolutely equal to zero. Because we work on the gray images, the maximum √ √ |Gσ ∗ I| is about 255 2. Thus the smallest F = 1/(1 + 255 2) ≈ 0.002, after several time steps, the curves will leave the edges. If we set the F values of edges to be zero, we can avoid this situation. But this improvement only works for gaps small enough. For large discontinuity, it is really the problem of edge detection. These problems are caused by that the edge functions are generated roughly from gradient. -33-.

(65) (a) σ = 0.75. (b) σ = 0.75. (c) σ = 2. (d) scale with a = 0.1. threshold d. (e) g(x, a) = 1/(1 + a |x|). Fig. 4.3. Different control of generating edge function. of gaussian smoothed image. The methods GAC and GFV design the curve to stop according to the edge function, that means where the final curve is, where the edges are found from gradient of gaussian smoothed image. Summarize the adjustment for F in the experiment • select suitable σ • select a to enhance edge function according the original edge function • select the threshold d to set F = 0 if F is less than the threshold. (a) No threshold. Fig. 4.4. (b) threshold with 0.5. (c) curve expand over the edge. Edge linking by threshold. -34-.

(66) 4.1.3 Vector Field Term The main difference between vector driven methods is the vector field generation. Before test the vector field generated by these methods, we explain the basic edge stopped method without external vector field. The basic method uses only κ and edge function F , and the initial curve is set to cover the objects or entirely inside an object. Then the curve expands or shrinks until reaches the edges. But after several time steps, the curves will leave the edges, as we have. (a). (b). Fig. 4.5. Simplest segmentation by active curve. discussed in last subsection. The improvement of GAC is to add the external vector field around the edges. This let the curve stop at the edges. We can see the same effect in the GVF too, since the GVF is extended from GAC. The drawback is that the curve still need to be initialized to cover the object or to be inside the object. The GVF propagates the vector field and emphasizes three improvements • the curve can be set across the object, • avoid the curve across the gap of edges, • curve can move to enter concave regions. We show this in Fig. 4.6, which the curve is initialized to across the object and is influenced by the vector field to surround the object without crossing top and bottom gaps. The firs improvement is because the vectors around the gaps will be generated after vector propagation and points to the opposite of edges as Fig. 4.6b shows. The second one can be used. -35-.

(67) (a). (b). Fig. 4.6. Improvement by GVF. to overcome the edge gaps problem as we have mentioned above. Unfortunately, the improvement only suits to the situation when there is only one object we want. In our experiment, we observe when the curve movement there are multi-object. The curve stops between the parallel corners or edges of objects.. Fig. 4.7. Multi-Object curve stop at gaps. The last improvement is that the curve can move into the concave regions of an object. Fig. 4.8 shows an example about this problem, the original −∇g in the concaves of w-shaped point horizontally and parallelly to left and right edges as shows in Fig. 4.8b. When the curve enters the concaves, it will be stopped. After propagation using GVF method, the vectors in the concave turn to point into concave as shown in Fig. 4.8c. Fig. 4.9 shows the curve evolution with −∇g and GVF respectively. But in our experiment we observe that this only works well if the concave region is open wide and reduce to narrow or almost the same width. If the concave region is from narrow to wide or with some bulge, the vectors after propagation will point outward the concave and the curve stays still. We can see this in Fig. 4.8d. Fig. 4.10 shows a clearer example, where vectors in the front part of concave point out. -36-.

(68) (a) w shape. (b) right concave of −∇g. Fig. 4.8. (c) right concave of GVF. (d) left concave of GVF. Vector Field on the concave region. (a) with −∇g. (b) with GV F. Fig. 4.9. The curve evolution under different vector field. and vectors in the rear area of concave point into concave. The GAC and GVF generate vector field based on gradient of edge function, i.e. F = g and V~ = −∇g, this means there are limits to the texture images. In contrast to generate vector field from edge function, EdgeFlow uses Gabor filter to give a better vector flow of texture images. Fig. 4.11 shows a texture box with similar intensity. In the experiment we utilizes the streamlines to compare the vector field generated by different parameters. The streamlines are the paths over which a dense number of free particles move under the influence of external forces when laced in the external force field. From the experiment, in the texture images, the σ has to set large, in order to avoid the disorderly vector groups which obstruct the curve evolution. We show this in the Fig. 4.12, with σ = 3 there are a little textures. -37-.

(69) (a) u shape. Fig. 4.10. (b) GVF. U-shape Concave region. in the inside box which are regarded as edges and results the vectors point to them. And we compare the EdgeFlows generated with intensity or texture only to look the effect of the prediction using Gabor filter. The most drawback of EdgeFlow is that it is very time consuming.. (a) Texture image. Box. (b) EdgeFlow of (a). (c) Segmentation with EdgeFlow. Fig. 4.11. EdgeFlow generate with 12 orientation, σ = 5, S = 5, K = 6, texture EdgeFlow only.. -38-.

(70) (a) σ = 5 texture only. Fig. 4.12. (b) σ = 3 texture only. (c) σ = 3 intensity only. Streamlines of EdgeFlow vector field. 4.1.4 Force Balance Problem We have tested the each vector driven method in terms of each force. In these methods we set the curve to cover, inside or across the objects. While the initial curve is far from the objects or edges, there are two methods to evolve curves. First method is to suppose that the vector field V~ is able to move the curves to the objects. But in the multi-object images or nature images, the vectors between the objects are disarray. And the curve will be ill evolved. Second is to set the curve expanding by F and will around the objects, this means in the smooth region the effect of ακ + βF is must larger than γ V~ , otherwise the curve will have ill-evolution as above. We show our experiment in Fig. 4.13. (a) Ill evolution. (b) Expanding around the objects. Fig. 4.13. The curve initialize outside the objects. But ακ + βF should also be less than γ V~ near the edges. Otherwise, the little values of ακ + βF will let the curve leave the edge after several time steps, the V~ is not able to push the curves back to the edges. And the curves will keep leave the edges.. -39-.

(71) F. F. V~. V~. to edge at edge leave edge. Edge. Summarize each parameter setting from experiment as follow • ακ + βF > γ V~ let curves expand to edges at smooth region. • ακ + βF small enough at edges avoid curve leave the edges. • ακ + βF < γ V~ let curves back to the edges when curves try to leave edges.. 4.2 Energy Balance Methods Instead of vector driven methods, we test on the region based methods. The region based equation is ". . φt = δ(φ) µ · div. ∇φ |∇φ|. . N N 1 X + 1 X − + 2 2 − λ (u0,i − ci ) + λ (u0,i − c− i ) N i=1 i N i=1 i. #. The algorithm of region based method is: Initialize φ0 , t = t0 for n less than maximum iterations Compute ci (φt ) ∇φ Compute div |∇φ| Update φt , t = t0 + n4t end The initial curve locations affect the segmented objects in region based method. Fig. 4.14 shows the segmentation progress of two initial curves. It is worthy to notice that the inside hole of the circle object is also segmented out by the curve. The time cost in region based method are much less than the vector driven methods, i.e. much fewer iterations of solving φ are needed.. -40-.

(72) (a). (b). Fig. 4.14. Region based method initial in different location. 4.3 Artificial Images We test the image segmentation on artificial images. In all tests the parameters are manually set to have best results. In Fig. 4.15 , we use GVF to segment the three objects shown. In this example, we also can see the curve merging and splitting. Then we add 15%, 30% and 45% Gaussian noise respectively, In the 15% Gaussian noise the GVF still can cut out the objects, as shown in Fig. 4.16, but with more noise it is hard to set a suitable parameter to find out the objects, since the methods are based on the edge and the edge function are noise sensitively. An ill-improvement is to set an initial curve which covers all objects then shrink it by a vector field. We use EdgeFlow to generate a under vector filed around edges. The segmented result with 45% gaussian noise added are shown in Fig. 4.17. With noise image, the region based method gives good results on Fig. 4.18a and Fig. 4.18b.. -41-.

(73) Fig. 4.15. Segment three objects by Gradient Vector Flow, σ = 0.5, µ = 0.1 with 60 iteration, α = 0.05, β = 0.2. Fig. 4.16. Added 15% Gaussian noise, σ = 2, µ = 0.1 with 60 iteration, α = 0.05, β = 0.3. Fig. 4.17. Added 45% Gaussian noise, using EdgeFlow, 8 orientation σ = 3, α = 0.1, β = 0.2. -42-.

(74) (a) Gaussian noise 30%. (b) Gaussian noise 45%. Fig. 4.18. Segmentation of the 30% and 45% Gaussian noise image using region based method. -43-.

(75) CHAPTER 5 CONCLUSIONS AND FUTURE WORKS. In this thesis, we survey the level set methods and their applications on image segmentation. We separate the conclusion discuss in two ways level set methods The level set methods are simpler to express the surfaces or curves, but they have the drawback of high time consuming on numerical PDE calculation. Since each implicit function φ can result in only one piece of segmentation, to get more segmentations we need more implicit functions to be initialized in different parts of an image. This also increases calculation time. Image segmentation The advantage we can find is its convenient for cutout the region after segmentation and the result of segmentation by curve is a complete region. The drawbacks of the level set methods in image segmentation are that: First, the level set is an active contour implementation method and the segmentation results are based on the methods of the image attribute analysis, for instance, edge detection uses gradient of Gaussian to smooth images, texture analysis utilizes Gabor filter. Second, the generated vector field and the driven curve are hardly to balance among the different velocities κ, F (Edge function) and V~ (External velocity field). The parameters need to be set differently case by case. Ideally, we hope we do not need to care about the parameters, α, β, of the level set curve evolution function and only need to care about the underlying designs of F and V~ . For this reason, the region based method seems to give an idea of obtaining another form of the level set function of curve evolution. Furthermore, to solve the parameter setting, designing a training system to get reasonable setting maybe workable.. -44-.

(76) REFERENCE. [1] S. Osher and J. A. Sethian, “Fronts propagating with curvature-dependent speed: Algorithms based on hamilton-jacobi formulations,” J. Comput. Phys., vol. 79, pp. 12–49, 1988. [2] J. A. Sethian, Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision and Materials Science.. Cam-. bridge University Press, 1988. [3] V. Caselles, R. Kimmel, and G. Sapiro, “Geodesic active contours,” Int. J. Comput. Vision, vol. 22, pp. 61–79, 1997. [4] C. Xu and J. L. Prince, “Gradient vector flow: A new external force for snakes,” in IEEE Pro. Conf. on Comp. Vis. Patt. Recog., 1997, pp. 66–71. [5] B. Sumengen, B. S. Manjunath, and C. Kenney, “Image segmentation using curve evolution,” in Proc. Asilomar Conf. Signals, Syst. Computers, vol. 2, 2001, pp. 1141–1145. [6] T. F. Chan and L. A. Vese, “An active contour model without edges,” in Scale-Space Theories in Computer Vision, 1999, pp. 141–151. [7] C. W. Shu and S. Osher, “Efficient implementation of essentially non-oscillatory shock capturing schems,” J. Comput. Phys., vol. 77, pp. 439–471, 1988. [8] V. Caselles, F. Catte, T. Coll, and F. Dibous, “A geometric model for active contours in image processing,” Numer. Math., vol. 66, no. 1, pp. 1–31, 1993. [9] R. Malladi, J. A. Sethian, and B. C. Vemuri, “Evolutionary fronts for topologyindependent shape modeling and recovery,” in Proc. European Conf. Computer Vision, 1994, pp. 3–13.. -45-.

(77) [10] C. Xu and J. L. Prince, “Snaks, shapes and gradient vector flow,” IEEE Trans. Image Processing, vol. 7, no. 3, pp. 359–369, 1998. [11] W. Y. Ma and B. S. Manjunath, “Edgeflow: A technique for boundary detection and image segmentation,” IEEE Trans. Image Processing, vol. 9, no. 8, pp. 1375–1388, 2000. [12] B. S. Manjunath and W. Y. Ma, “Texture features for browsing and retrieval of image data,” IEEE Trans. Pattern Anal. Machine Intell., vol. 18, no. 8, pp. 837–842, 1996. [13] D. Mumford and J. Shah, “Optimal approximation by piecewise smooth functions and associated variational problems,” Commun. Pure Appl. Math, vol. 42, pp. 577–685, 1989. [14] T. F. Chan, B. Y. Sanberg, and L. A. Vese, “Active contours without edges for vectorvalued images,” J. Visual Commun. and Image Rep., vol. 11, pp. 130–141, 2000. [15] T. F. Chan and L. A. Vese, “Active contours without edges,” IEEE Trans. Image Processing, vol. 10, no. 2, pp. 266–277, 2001. [16] S. Osher and R. Fedkiw, Level Set Methods and Dynamic Implicit Surfaces. Verlag, 2002.. -46-. Springer-.

(78)