Lennard-Jones force field for geometric active contour

ARTICLE IN PRESS Signal Processing 90 (2010) 1249–1266

Contents lists available at ScienceDirect

Signal Processing journal homepage: www.elsevier.com/locate/sigpro

Lennard-Jones force field for geometric active contour Zhenglong Li a,, Qingshan Liu a,b,, Hanqing Lu a, Dimitris N. Metaxas b a b

National Laboratory of Pattern Recognition, Institute of Automation, Chinese Academy of Sciences, P.O. Box 2728, Beijing 100190, PR China Department of Computer Sciences, Rutgers, The State University of New Jersey, 110 Frelinghuysen Road, Piscataway, NJ 08854-8019, USA

a r t i c l e in fo

abstract

Article history: Received 5 September 2008 Received in revised form 28 June 2009 Accepted 12 October 2009 Available online 21 October 2009

This paper presents a new geometric active contour (GAC) model based on Lennard-Jones (L-J) force field, which is inspired by the theory of intermolecular interaction. Different from conventional gradient based GAC models, the proposed model does not rely on any pre-computed edge map. We take each pixel of image as a particle, and design an L-J force field for GAC model according to interaction between pixels. We introduce a parameter of distance regularization to make the force tunable, and define an energy function for the L-J function to integrate various image features efficiently. A switch parameter c generates two different characteristics for the L-J force field: in the case of c ¼ 0, the force vector flows bi-directionally converge to boundaries, while in the case of ca0 a morphological effect is formed. To quantitatively evaluate the proposed force field, we present two criteria: directional uniformity and average amplitude, and compare it with the popular GVF. In addition, we also prove that the so-called electrostatic model is actually a special case of the proposed model. & 2009 Elsevier B.V. All rights reserved.

Keywords: Image segmentation Geometric active contour Lennard-Jones force

1. Introduction Active contours have been widely used in the communities of image processing. From their mathematical representations, active contours can be classified into two categories: parametric active contours (PAC) [1,2] and geometric active contours [3,4]. PACs employ the explicit Lagrangian formulation, and it is difficult for them to handle topological changes of curves, i.e., curve split and merging, during evolution. GACs use an implicit Eulerian formulation which can be solved by level set method [5]. Compared with PACs, GACs can deal with topological changes automatically. For this advantage, we will focus on the latter in this paper. GACs were originally introduced in [3,4,6], respectively. The main idea is to use an image-related force field to

 Corresponding authors at: National Laboratory of Pattern Recognition, Institute of Automation, Chinese Academy of Sciences, P.O. Box 2728, Beijing 100190, P.R. China. E-mail addresses: [email protected] (Z. Li), [email protected] (Q. Liu).

0165-1684/$ - see front matter & 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.sigpro.2009.10.008

drive implicitly represented curves to positions of interests. Denote curve as C ¼ ½x; yT , and image as I . The curve evolution is described by a partial differential equation as follows:

Ct ¼ ðgðjrI ðCÞjÞðkðCÞ þ bÞ þ F Img  NÞN;

ð1Þ

where F Img ¼ rgðjrI ðCÞjÞ is an edge gradient deduced force field. r is the gradient operator, and the symbol ‘‘  ’’ stands for the inner product operator. In (1), k is the curvature regularization term, and N is the inward unit normal to C. The function gðÞ is a monotonically decreasing function such that g : ½0; þ1Þ/Rþ , with gð0Þ ¼ 1 and gðþ1Þ ¼ 0; gðÞ serves as an indicator of boundaries to trap the moving curve. The constant term b pushes the curve from farther feature-lacking region to positions of interests [7]. Note that an improper b will be prone to resulting in edge leakage or sensitivity to noise. Numerically, the solution of equation (1) can be obtained by the level set method [8].

ARTICLE IN PRESS 1250

Z. Li et al. / Signal Processing 90 (2010) 1249–1266

According to [8], the curve evolution of (1) can be expressed in the form of the net force F : Ct ¼ FN. The net force F is composed of three terms: F ¼ FProp þ FCurv þ FAdvec , where FProp , FCurv and FAdvec are propagation expansion speed, regularization, and advection term, respectively. Furthermore let C be embedded in a level set function f, i.e., fðCÞ ¼ 0, and (1) can be rewritten as

ft þ FProp jrfj þ FCurv jrfj þ FAdvec jrfj ¼ 0;

ð2Þ

where FProp ¼ b is the constant propagation expansion term (or the balloon force), FCurv ¼ k is a geometric regularization to ensure smoothness of C, and FAdvec ¼ ½u; vT  N the image-related advection term. The geometric properties of C, such as normal and curvature, can be represented by f [8], i.e., N ¼ rf=jrfj and k ¼ div rf=jrfj. Therefore (2) can be reformulated as

ft þ FProp jrfj þ FCurv jrfj þ ½u; vT  rf ¼ 0:

ð3Þ

In (2), FAdvec is coupled with the image I by the imagerelated force field ½u; vT , and it plays an important role in GAC. A popular and intuitive choice for ½u; vT is the gradient, i.e., ½u; vT ¼ rgðjrI jÞ. However, such a gradient based force field has two issues: (1) its support range is locally bounded; (2) it is only based on two-dimensional grayscale gradient information, and much important image information is ignored, such as color and texture etc., so it cannot deal with some complicated image cases well. To circumvent the limited support range of the gradient related force field, the gradient vector flow (GVF) [9] was successfully introduced into GAC [10–12]. The GVF is performed on a pre-computed edge map, and it essentially extends the support range of the gradient by a generalized diffusion process. It means the pre-computed edge map has great influence on the performance of the GVF. In [13,14], the electrostatic model was proposed for active contour model. However, it still relies on the precomputed edge map to generate the desirable force field as the GVF does. In order to overcome the insufficiency of the twodimensional gradient information, some techniques have been proposed. A so-called ‘‘generalized gradient’’ [15] defined in RN was integrated into GAC in [10,12,16]. As well as the two-dimensional gradient, the ‘‘generalized gradient’’ has the problem of the myopic support range too. Ref. [12] used EdgeFlow [17] to incorporate multicues for GAC. EdgeFlow is in essence an edge gradient based method, so it is sensitive to noise. Ref. [10] employed the mean shift clustering [18] as a pre-process to integrate the region information (mainly color). In [19,20], the statistical information is employed by using Gaussian modeling. However, there is no unified framework to integrate various image information in GAC so far. In the paper, we propose a new GAC based on the L-J force field to handle the above issues. In the theory of intermolecular interaction, there exist attractive and repulsive interactions between pairwise particles. The attraction predominates when two particles are in a moderate distance, while the repulsive effect becomes dominant if they reach each other too close. This is the

physics essence of the L-J force. Inspired by this physical principle, we regard each pixel in an image as a particle, and then an L-J force field is designed in image domain according to the interactions between the pixels. In order to integrate various image information into GAC, we define an energy function to replace the energy constant, which is embodied by a distance metric to the features of pairwise pixels. Moreover, different emphasis is put on different pairwise pixels by this energy function. We also introduce a parameter of distance regularization to make the force tunable (cf. Eq. (7)). A switch parameter c makes the L-J force field have two different characteristics, and it makes the proposed GAC model applicable for different applications. To quantitatively evaluate the proposed force field, we compare it with the popular GVF based on two criteria, i.e., directional uniformity of force field and average amplitude (cf. Section 4.1). To our best knowledge, we are the first one to evaluate the force filed in GAC quantitatively. In addition, the so-called electrostatic model [13,14] is actually a special case of the proposed model. We conducted the experiments on various images including medical images, natural images, and text line detection in handwritten documents, and the experimental results showed that the proposed method has an encouraging performance. The rest of the paper is organized as follows. The concept of the L-J force is introduced in Section 2. In Section 3, we first present the L-J force field in the image domain, and then we analyze its characteristics. Finally, we address how to integrate the L-J force field into GAC. Section 4 discusses the quantitative criteria for the force field evaluation, and the comparison between the proposed method and the GVF; we also prove that the electrostatic model is actually a special case of the proposed L-J force model in this section. Section 5 reports the experimental results, and Section 6 concludes the paper. 2. Introduction of Lennard-Jones force The L-J potential ULJ is a physical concept to approximately describe interactions between pairwise particles, such as atoms or molecules [21,22]. Given a system composed of two particles fe1 ; e2 g located at positions fs1 ; s2 g, respectively, the L-J potential equation is expressed as   m  n  s s ; ð4Þ  c2 ULJ ðs1 ; s2 Þ ¼ e c1 jrj jrj where jrj ¼ js1  s2 j is the interatomic separation (spacing distance); s is the particle diameter, and c1 , c2 , m and n ðm4nÞ are constants in which the constant e represents the amplitude of potential energy (or well depth). The L-J force is computed as following:   smþ1 snþ1 ð5Þ F LJ ðs1 ; s2 Þ ¼ rr ULJ ðs1 ; s2 Þ ¼ e c mþ2  nþ2 r; jrj jrj where e ¼ enc2 and c ¼ mc1 =nc2 . It is the negative gradient of ULJ w.r.t. r. The L-J force is composed of two different parts: the attractive force and the repulsive force. The attractive force, represented by ðsnþ1 =jrjnþ2 Þr, is

ARTICLE IN PRESS Z. Li et al. / Signal Processing 90 (2010) 1249–1266

the L-J force was exerted on the particles to keep the elastic connections, yet allowing for the rearrangement of elements over large deformations. In this section, we propose it to improve the performance of GAC. We first present the L-J force field in image domain and analyze its characteristics. Then we discuss how to integrate image information in the L-J force field. Finally we address how to design the L-J force field driven GAC.

Interaction Force (in nN)

3 Lennard−Jones Force Attractive Effect

2

Repulsive Effect

1 r0

0

1251

−1

3.1. L-J force field in image domain −2 −3 0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

Distance (in nm) Fig. 1. L-J force between the probe tip of scanning probe microscopy and testing sample (the ordinate is rescaled for display) [23].

collectively called van der Waals forces, or dispersion force. The repulsive force is corresponding to cðsmþ1 =jrjmþ2 Þr, and it is called Pauli repulsion, which is caused by overlapping of electronic clouds surrounding the particles. When two particles are at a moderate distance, the attractive force is dominant, and they mildly attract each other; when the separation between the particles is short, they strongly repel each other due to the leading repulsive force. To illustrate the L-J force, we give a real-world example: the interaction between the probe tip of scanning probe microscopy and testing sample (tipsample)1 [23]. Fig. 1 shows the plot of tip-sample L-J force. The x-axis represents the tip-sample distance, and the y-axis denotes the force amplitude. The sign of amplitude represents the direction of the force: the plus sign for the repulsive force and the negative for the attractive force. The repulsive force (dot-dashed red curve) is dominant when the distance is small, while the attractive force (dashed blue curve) gets predominant as the distance increasing. In the figure there are two critical points: the equilibrium and the extremum points. The shape of the L-J force curve is determined by the two points. The equilibrium point r0 is the zero-crossing point of the plot, which splits the x-axis into two parts: the repulsion-predominant part and the attraction-predominant part, i.e., F LJ jjrj¼r0 ¼ 0, and r0 ¼ c1=ðmnÞ s. The extremum point r1 makes (5) reach the minimum value at r1 , i.e., ðd=drÞF LJ jjrj¼r1 ¼ 0, and r1 ¼ ððn þ 1Þ= ðm þ 1ÞÞ1=ðnmÞ r0 . 3. L-J force field and the proposed GAC model The L-J force was successfully used by Tonnesen and Szeliski et al. [24,25] in the field of computer graphics for surface modeling and reconstruction. They modeled fluids and deformable solids as a dynamic particle system, and 1 For the scanning probe microscopy, probe tip and testing sample can be approximately treated as particles for their comparable scale to molecule or atom.

Denote the image I as fF ; Sg, where F ¼ fe1 ; e2 ; . . . ; eM g is the feature set for pixels, and S ¼ fs1 ; s2 ; . . . ; sM g are the pixel positions on the image plane. I is considered as a grid system where each pixel at sðÞ is a particle, and each pixel interacts with all the other pixels. Because the interaction between two particles becomes trivial when the distance between them is large enough, we just take account of the interaction in a local neighborhood of si , X  F LJ ðsi ; sj Þ; ð6Þ F LJimg ðsi ; N i Þ ¼ sj 2N i

F LJ ðsi ; sj Þ

where is the pairwise interaction between two pixels, and neighborhood N i ¼ fsj : jsj  si joConstant and jaig. F LJ ðsi ; sj Þ is modified from the original L-J force function (5) in two aspects. One is that we replace the energy constant e with an energy function, which is defined on F as eðei ; ej Þ ¼ expðD2 ðei ; ej Þ=kÞ, where Dðei ; ej Þ is a distance metric (in the paper the Euclidean metric is used simply), and the positive constant k controls the interaction intensity between pixel si and sj . By the energy function eðÞ, various image information can be easily integrated besides the conventional two-dimensional grayscale gradient information that is popular in the previous works. A small k means much intense interaction and more details are preserved; a large k will result in strong smoothing effect and coarse result will be obtained. The details are discussed in Section 5. L-J force in image domain is derived from the concept of L-J potential in physics. The essence of L-J force in image domain is the energy function eðei ; ej Þ that reflects simulated interaction between pixels. Any proper energy function that reflects this simulated interaction property can be used in the proposed method. For example we can design a new energy function just by replacing the Euclidean metric with other proper metric in (5), and a good option is the metric proposed in [26]. Another is that we introduce a control parameter d to the distance jrj as a regularization and control parameter, !"  mþ1 D2 ðei ; ej Þ s c F LJ ðsi ; sj Þ ¼ exp  k jrj þ d  nþ1 # s r  ; ð7Þ jrj þ d jrj where r ¼ si  sj . Given the parameters m and n, the support effect of the force is approximately proportional to d. A large d means large support response and high average amplitude (cf. Eq. (13) in Section 4.1.2), and vice versa. We will analyze its performance in Section 4.2.

ARTICLE IN PRESS 1252

Z. Li et al. / Signal Processing 90 (2010) 1249–1266

According to the inverse-square law in physics, we choose m ¼ 2, n ¼ 1, and s ¼ 1 for simplicity in the paper. We take the parameter c as a zero-nonzero switch, which controls two different characteristics of the force field. When c ¼ 0, only the attractive part takes effect. In this case, the vector flows bi-directionally converge to the boundaries. If ca0, the combination of the attractive effect and the repulsive effect produces equilibrium adjacent to the boundary with a distance. In this situation, the morphological dilation effect is achieved by the equilibrium positions. We will use these two different characteristics for different applications (see details in Section 5). In the following, we take two toy examples, the line edge and the step edge, to show these characteristics of the L-J force intuitively. 3.1.1. Case of c ¼ 0 This case means that we only use the attractive term in (7), ! D2 ðei ; ej Þ 1 r : ð8Þ F Attr ¼ exp  k ðjrj þ dÞ2 jrj Fig. 2(a) shows the attraction performance of L-J force field on the line edge. The top of Fig. 2(a) is the vector flows of the line edge map, and the amplitude response of L-J force field is shown at the bottom. The sign of amplitude response curve represents the direction of vector flows: the negative stands for direction from right to left, and the positive for the opposite direction. In Fig. 2(a), the vector flows bi-directionally converge to the line edge, i.e., the L-J force forms a head-to-head vector flows across the line edge. The performance of L-J force field on the step edge is illustrated in Fig. 2(b). The vector flows and the amplitude responses of the step edge are similar to those of the line edge in Fig. 2(a). It implies that the L-J force field not only forms the vector flows bidirectionally converging to the boundaries in the line-like edge map, but also it can achieve the same effect to the step-like edge map or piecewise-content image.2 3.1.2. Case of ca0 The case of ca0 means both the attractive force and the repulsive force are taken into account. Fig. 2(c) illustrates the performance of L-J force field on the line edge map in this case, where the top shows the vector flows of the line edge and the bottom is the amplitude response curve. Across the line edge, the vector flows converge to the equilibrium positions at both sides of the line edge with a distance to the line edge. The equilibriums correspond to the zero-crossing points in the amplitude response curve. This phenomenon can be used to form the morphological dilation effect. The L-J force field of the step edge is shown in Fig. 2(d), and it is similar to Fig. 2(c): for either the step edge or the line edge the L-J force field can achieve the same morphological dilation effect.

3.2. L-J force driven GAC Our purpose is to design a GAC model based on the L-J force field. An intuitive way is to replace F Img in (1) with the L-J force F LJimg , i.e., Ct ¼ ðkðCÞ þ b þ F LJimg  NÞN, where k and b are same as in (1). However, the constant propagation expansion term b usually cannot respond properly to the local situation. At the boundary positions whose responses are weaker than b, i.e., the weak edges, the evolving curve will be divergent. If we drop out the balloon force or set it with a small value improperly, the curve evolution will show the phenomenon of conglutination at those positions where F  N  0, or where the amplitude of the local force is too weak to drive C. To deal with the above issues, we use an adaptive propagation expansion term b according to the local force field [11]: b ðsÞ ¼ sign ðFðsÞ  NÞgðjFðsÞ  NjÞ;

ð9Þ



where sign ðxÞ is a sign function to indicate the local direction change (inflation or deflation) of the curve: sign ðxÞ ¼ 1 when xZb, otherwise sign ðxÞ ¼ 1. b is a constant to balance the expansion and the contraction effect. And gðÞ is a monotonically decreasing function that adaptively controls the curve evolution according to the local vector field. Therefore the proposed L-J force driven GAC model is given as Ct ¼ ðkðCÞ þ b ðCÞ þ F LJimg  NÞN:

ð10Þ

U shape is often used to test the performance of the GAC model for the special case with deep concave shape. In [9], GVF showed good property for U shape. Fig. 3 shows the results of U shape by the proposed method and GVF. It can be found that the proposed L-J force driven GAC deals with the U shape as well as GVF. Fig. 4 shows the performances with different expansion term schemes. Fig. 4(c) is the result with b ¼ 0. The curve cannot go into the deep narrow gaps, because the vectors near the boundaries are nearly perpendicular to the boundaries and the width of the gap is very narrow so that there are no vectors digging into the depth of the gap. Fig. 4(d) shows the result using a constant expansion term ðb ¼ 0:3Þ. Since b is relatively larger than some edge responses, some weak edges of matches are leaked; a constant expansion term is not suitable for the situation where the edge responses are non-uniform. Fig. 4(b) shows the result using the local adaptive expansion term. The whole boundaries of matches are well segmented. 4. Analysis of vector flow field In this section, we first discuss how to evaluate the force field in GAC model quantitatively. Then we will give a comparison between the L-J force field and GVF. At last, we will prove that the electrostatic model [13,20] is actually a special case of the L-J force field model. 4.1. Two quantitative criteria for the force field evaluation

2

By selecting proper DðÞ for feature set F , image can be decomposed into approximately ‘‘stepwise’’ components in the sense of homogeneity of metric.

To our best knowledge, all the previous works evaluated the force field used in GAC by observing the

ARTICLE IN PRESS Z. Li et al. / Signal Processing 90 (2010) 1249–1266

1 0.8 0.6 0.4 0.2 0 −0.2 −0.4 −0.6 −0.8 −1

1253

1 0.8 0.6 0.4 0.2 0 −0.2 −0.4 −0.6 −0.8 −1 0

10

20

30

40

50

60

0

10

20

30

40

50

60

1 0.8 0.6 0.4 0.2 0 −0.2 −0.4 −0.6 −0.8 −1

1 0.8 0.6 0.4 0.2 0 −0.2 −0.4 −0.6 −0.8 −1 0

10

20

30

40

50

60

0

10

20

30

40

50

60

Fig. 2. L-J force field of the line edge and the step edge. Top row: the vector fields; bottom row: amplitude response curves. (a) The case of c ¼ 0 on the line edge. (b) The case of c ¼ 0 on the step edge. (c) The case of ca0 on the line edge ðc ¼ 5Þ. (d) The case of ca0 on the step edge ðc ¼ 20Þ.

final segmentation results qualitatively, and no objective measurement was given. In this paper, we present two quantitative criteria to evaluate the force field used, i.e., directional uniformity and average amplitude. Our design is based on the fact that the direction and the amplitude can determine one point in the two-dimensional force field. A good force field should have the locally uniform directions except the positions around boundaries, and it should be enough in the amplitude even at those positions that are far from the boundaries.

4.1.1. Directional uniformity of the force field The moving direction of the curve is consistent with the direction of the force vectors. Thus, a desirable force field should have head-to-head force vector flows around the boundaries, i.e., the vector flows ‘‘collided’’ at the boundaries to trap the curve. In those regions far from the boundaries, the vector flows should be smoothly and avoid being collided together, i.e., no sinks or curls. We use divergence to indicate where new flow is generated (source) or current flow is annihilated (sink) at a specified

ARTICLE IN PRESS 1254

Z. Li et al. / Signal Processing 90 (2010) 1249–1266

Fig. 3. U-shape segmentation. (a) Result by GVF. (b) Result by L-J model.

Fig. 4. Results by various expansion terms (balloon force). (a) The initial curve. (b) The result using local adaptive expansion term (cf. Eq. (9)). (c) Conglutination with expansion term b ¼ 0. (d) Edge leakage caused by a constant expansion term.

point, and to give the change rate quantitatively. Since we only focus on the directional information of the force field in this case, we define the directional uniformity of the force field at one point s in the force field F as   FðsÞ FðsÞ MðsÞ ¼ neg div div ; ð11Þ jFðsÞj jFðsÞj where the function negðÞ is to preserve the responses of the collided vector flows. If xo0, negðÞ ¼ 1; otherwise, it is equal to 0. The collided vector flows mean that the vector flows ‘‘disappear’’ at the position, i.e., sinks. We only focus on the collided vector flows, for those departed (tail-to-tail) vector flows have no effect on the convergence of the curve. div FðsÞ=jFðsÞj has a strong response at the position of vector collision or large direction torsion, whereas it produces small response at those positions with smooth vector flows. Based on all the points’ uniformity, MðsÞ’s, we can draw the directional uniformity map of the force field. MðsÞ measures the local directional uniformity at one point, and the directional uniformity map gives an intuitive illustration of the directional uniformity of the force field. To overview the directional uniformity of the force field quantitatively, we use the information entropy of directional uniformity map (EDU) to give a quantitative measure, X ð12Þ HEDU ðMÞ ¼  PrðMðsÞÞ log2 PrðMðsÞÞ; s2S

where MðÞ is the directional uniformity map, and PrðÞ is the probability mass function.

4.1.2. Average amplitude of the force field How to evaluate the support range of the force field quantitatively is another issue. The amplitude response of the force to the boundary is a good indicator of the capture range. A usual option is bandwidth, which means the length between the peak response (at the boundary) and the response at the position far from the boundary.3 For example, 10 dB-bandwidth is the length between the peak response and 10% of the peak value. Based on the concept of bandwidth, we define the average amplitude to evaluate the force field, which is the ratio between the integral of the force amplitude (the area under the response curve) over an interval and the length of this interval: AðFÞ ¼

1 L

Z

L

jFj dx;

ð13Þ

0

where L is the length between the peak value and the cutoff amplitude, e.g., 10% of the peak value. 3 For the two-dimensional force field, the bandwidth is defined along the direction perpendicular to the boundary. In experiments, we use the line edge as the benchmark input.

ARTICLE IN PRESS Z. Li et al. / Signal Processing 90 (2010) 1249–1266

1255

1 0.9 Amplitude of response vectors

0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0

10

20

30

40

50

60

Distance (in pixel) to edge. Fig. 5. Comparison between the L-J force and the GVF force responses to the line edge with different d and m. Solid black curves: the responses of the L-J force. Dashed red curves: the responses of the GVF. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Table 1 The 10 dB-bandwidth and the average amplitude.

4.2. Comparison with GVF With the help of the defined criteria, we compare the proposed L-J force field with the GVF from the aspects of robustness and support range in this subsection. For a fair comparison, we take only the attractive part of the L-J force into account (cf. (8)). Before we do comparison between the L-J force field and GVF, we briefly review the idea of GVF. The essence of GVF is a generalized diffusion process as follows:      u u ¼ jrf j2  rf ; mr2 ð14Þ v v

d

0.1

0.5

0.9

1.3

L-J force Bandwidth Average amplitude

13.853 0.086

19.752 0.119

25.267 0.153

32.941 0.188

0.14

0.18

0.22

10.288 0.087

11.152 0.098

12.128 0.107

m

0.1

GVF Bandwidth Average amplitude

9.607 0.075

2

where r is the Laplacian operator,4 and f is an edge function, ½u; vT is the vector for diffusion. In (14), the gradient of the edge function, rf ¼ ½@f =@x; @f =@yT , can be replaced by any edge deduced vectors. The term r2 ½uv corresponds to the diffusion term. The right-hand side of (14) factorizes into two factors: the first factor controls the diffusion intensity, while the second term, ð½uv  rf Þ, guarantees ½u; vT not to deviate far from the original rf , i.e., the fidelity term. When jrf j2  0, the whole equation is dominated by the diffusion term. In (14), the GVF is coupled with the image via the gradient of edge function rf . To solve (14) numerically, u, v are treated as the functions of time, i.e., u ¼ uðx; y; tÞ, v ¼ vðx; y; tÞ (cf. (16a) and (16b) in [9] for the detailed discretization scheme). According to the Courant–Friedrichs–Lewy condition (CFL condition) [27], the time step Dt should satisfy 2

4 Here the Laplacian operator r is redefined for vector calculus: for 2 2 2 x ¼ ½x1 ; . . . ; xn T , r x ¼ ½r x1 ; . . . ; r xn T .

DtrDxDy=4m to ensure the convergence, where Dx and Dy are the spatial step (usually Dx ¼ Dy). In [9], m is interpreted as a regularization term to balance the diffusion effect and the fidelity. Moreover, [9] suggested to set a large m to intensely smooth the resulting vector field under noisy environment. Actually the effect of m is very limited (for details see Fig. 5 and Table 1). Following the CFL condition, there is moDxDy=4Dt ¼ 1=4Dt (if set Dx ¼ Dy ¼ 1). When Dt ¼ 1, m is upper bounded by the CFL condition of 0.25. An empirical value for m was 0.2 in [9]. Therefore, it is not practical for GVF to enhance the robustness and to extend the support range by increasing m without limitation. An alternative is to decrease the spatial step and the time step to get a large m, and not to violate the CFL condition. However, high spatial resolution will incur image interpolation and numerical problems, while high time resolution will slow down the numerical solving process.

ARTICLE IN PRESS 1256

Z. Li et al. / Signal Processing 90 (2010) 1249–1266

In addition, the result by GVF is also influenced by the number of iterations. When the number of iterations is large, resulting vector flows can be smooth. It is possible to improve the vector flows by increasing the number of iterations, although the number of iterations in practice is often set empirically. In contrast with GVF, L-J force field is computed in one pass. 4.2.1. Support range The purpose of GVF is to enlarge the support range of the gradient force. As we mentioned above, it is practically controlled by the parameter m in (14). In the L-J force, the regularization constant d in (8) plays a similar role to the parameter m in the GVF. Without loss of generality, we take the line edge as example to investigate the influence of m and d on GVF and the L-J force, respectively. Fig. 5 shows the amplitude responses with different m and d, where the solid curves and the dashed curves represent the response amplitudes of the L-J force and GVF, respectively. Table 1 lists the values of m and d that are used in Fig. 5. According to the CFL condition, m for the GVF has to be limited in a small range, while d in the proposed method has no such constraint. In Fig. 5, the amplitude responses of the GVF decay much rapidly compared to those of the L-J force. It means that the support range of the L-J force field is larger than that of the GVF force field. Table 1 lists the 10 dB-bandwidth and the average amplitudes of two force fields with different d and m, respectively. We can see that the support range of the L-J force gets enlarged notably with the increase of d, whereas m only influence the support range in a very limited manner. 4.2.2. Robustness As mentioned before, the GVF is performed on the precomputed edge map f ðÞ, therefore the quality of the edge map will have much influence on the force field. Different from GVF, the L-J force field is directly computed from the image data, and it is more robust than the GVF. We take two real medical images for example to compare the robustness. The implementation of the GVF refers to [9], including the extraction of the pre-computed edge map, i.e., a Gaussian filter of standard variance 1.5 is used first, and then the sided difference operator is performed to extract the gradient map. The L-J force field here is using grayscale feature, and the parameter setting is given in Section 5. Fig. 6(a) is the original CT image of human lung. The left-most subfigure of Fig. 6(b) shows the amplitude map of the L-J force field, and the right-most is the directional uniformity map of the L-J force field. The middle subfigure illustrates the zoomed force field of the blue color masked region in the amplitude map. From the zoomed force field, we can see clearly that the vector flows converge to the boundaries bi-directionally, and the force field is very smooth. Only along the boundaries, there exist the strong responses. In the cavity of lung, the white fiber-like tissues do not produce the strong non-uniformity responses, and the cavity part of the L-J force field is clear and smooth, due to the aggregate effect of the L-J force field.

Correspondingly, Fig. 6(c) lists the GVF’s amplitude map, the zoomed force field, and the directional uniformity map. Compared to the L-J force field, the force field is not so smooth as the L-J force field, even though the GVF vector flows also converge to the boundaries bidirectionally. At the bottom-right part of the GVF force field, the directions of the force vectors are not uniform, and some curls are generated. Non-uniformities also appear in the cavity part in the directional uniformity map. The GVF’s energy is intensely concentrated near edges in a narrow region, whereas the L-J force field’s energy is distributed extensively. The EDU of the L-J force field is 2.4930, while the EDU of the GVF force field is 3.0895. Fig. 7 shows the results of the MRI image of human left ventricle. Fig. 7(b) shows the L-J force field, and Fig. 7(c) is for the GVF force field. Because the imaging mechanics of MRI, the image of ventricle is not sharp and clear as usual optical images. In spite of the influence of noise, the L-J force field outperforms the GVF. On this image, we also investigate the performance of the GVF force field when the CFL condition is violated. The results are shown in Fig. 7(d), we can see that the force field becomes very disordered in this case. The EDUs of the L-J force field, the convergent GVF force field, and the divergent GVF force field are 2.9649, 3.9720, and 4.0217, respectively. In Figs. 8 and 9 we give the comparisons between the L-J force driven GAC and GVF snakes. Fig. 8 shows the comparison on segmentation of the CT image of human lung. Fig. 8(a) gives the result by the proposed method, and Fig. 8(b) by GVF snakes as comparison. In Fig. 8(a) the curve evolution process is shown from left to right and then from top to bottom. The initial curve is not totally in the cavities (the black parts) of the lung: the parts of the curve are outside of or inside of the cavities. The curves are attracted to the boundaries wherever the parts of the curves are outside of or inside of the boundaries: the inside parts of the curve inflate to the boundaries, while the outside parts deflate to the boundaries. Finally the initial curve splits into two curves convergent to the boundaries of the cavities. In contrast with the proposed method, GVF snakes are easily influenced by the white fiber-like tissues. In the right cavity the contour is attracted to the white tissues, and there exists the same phenomenon in the left cavity. That is because the vector field formed by GVF is not smooth as the L-J method. Moreover GVF snakes cannot handle the topological changes: it cannot split to correspond to the two cavities. Fig. 9 shows segmentation of the MRI image of left human ventricle. Due to imaging characteristics of MRI, the ventricle image is noisy and not sharp or clear as its optical counterpart. The contrast between the foreground and the background is very low, and on the top part of the ventricle there is a weak edge. Fig. 9(a) is the result by L-J force driven GAC. We use a low k to suppress the noise influence for L-J force field, and the contours converge well to the boundary, and the weak edge is correctly extracted. Fig. 9(b) shows the result by GVF snakes. The boundary of the ventricle is not well extracted.

ARTICLE IN PRESS Z. Li et al. / Signal Processing 90 (2010) 1249–1266

1257

HEDU = 2.4930

HEDU = 3.0895 Fig. 6. The L-J force field and the GVF force field of the CT image of human lung. (a) The original CT image. (b) The L-J force field: (from left to right) the amplitude map, the zoomed part, and the directional uniformity map. (c) The GVF force field: the amplitude map, the zoomed part, and the directional uniformity map (m ¼ 0:2o0:25 satisfying the CFL condition).

4.3. Relation to electrostatic model In [13,14] the electrostatic models were proposed to produce the force field. Similar to the GVF, it is based on a pre-computed edge map, but it regards each point of edge as an electric charge. Let F denotes the charge quantity (1 is for edges and 0 is for background), and S still stands for the position of each pixel. Then according to Coulomb’s law, the electrostatic force between unit charge on si and one edge point sj is given as (assume a unit charge is placed at si ) F Cou ðs1 ; s2 Þ ¼ kC

r ; jrj3

ð15Þ

where kC is the Coulomb constant, and r ¼ s1  s2 . The total electrostatic force from all the edge points P to si is calculated as F Cou ðsi ; BÞ ¼ sj 2B F Cou ðsi ; sj Þ

where B ¼ fsi : si is an edge pointg. Consider (8), F Attr ¼ eðei ; ej Þð1=ðjrj þ dÞ2 Þr=jrj, and let eðei ; ej Þ ¼ kCou , d ¼ 0, the L-J force FLJ is degenerated to the electrostatic force F Cou . Therefore the electrostatic force model is actually a special case of the proposed L-J force field model. 5. Experiments We test the proposed L-J force driven GAC on several kinds of images. The experiments can be classified into four groups: (1) the experimental data are grayscale images (mainly medical images and one blob image). We use only the grayscale information in this group, i.e., the feature set F is the image grayscale set; (2) we investigate to integrate the color feature and the texture feature, respectively; (3) we evaluate the performance of GAC integrating with the statistical information, where we

ARTICLE IN PRESS 1258

Z. Li et al. / Signal Processing 90 (2010) 1249–1266

HEDU = 2.9649

HEDU = 3.9720

HEDU = 4.0217 Fig. 7. The L-J force field and the GVF force field of the MRI image of left human ventricle. (a) The original image. (b), (c) and (d) are for L-J force field, convergent GVF, and divergent GVF, respectively: the first column is the amplitude map, the second is the zoomed force field, and the last column is the directionally uniformity map. In (c) m ¼ 0:2o0:25 satisfies the CFL condition. In (d) a divergent result is given: m ¼ 0:340:25 violating the CFL condition.

design a strategy to integrate the statistical information into F ; (4) finally, we show the morphological dilation effect of the proposed method by the applications of text

line extraction in document layout analysis. The former three group experiments employ the attractive part of the L-J force field, and the last group experiment considers the

ARTICLE IN PRESS Z. Li et al. / Signal Processing 90 (2010) 1249–1266

1259

Fig. 8. Segmentation of the CT image of human lung. (a) Result by L-J force driven GAC. (b) Result by GVF snakes.

combination effect of attraction and repulsion of the L-J force. In (7) and (8), the control parameter d equals 2 for all the experiments, but the value of parameter k is set according to the image content. Under noisy environment, a large k can suppress the interaction to emphasize the smoothing of the aggregate effect. When the contrast between foreground and background is low, a low k is needed to encourage the interaction to include details as much as possible. Therefore, we set a large k as 2 to suppress the noise in the MRI image of ventricle (Fig. 7), the blob image (Fig. 11), and the X-ray image of the Cygnus loop (Fig. 15). For the B-mode ultrasonic image (Fig. 10) with low contrast and the squirrel image in Fig. 13 with rich details, k is given a small value, say 0.1 and 0.5, respectively. For the other images, k ¼ 1. The size of the neighborhood N is another parameter, and it is set according to the foreground size empirically, i.e., the

radius of N is between the half size of the object and double the size of the object.

5.1. Group one In this group, we test the proposed method on medical image5 and microscopy image. The type of medical image is B-mode ultrasonic image. F is composed of the image grayscale values. Fig. 10 shows the curve evolution for a B-mode ultrasonic image of human ventricle. Due to the characteristics of ultrasound imaging, there are a lot of speckle noises in the image. A standard anisotropic 5 For other types of medical images, i.e., magnetic resonance imaging (MRI) and computer tomography (CT), see Figs. 8 and 9 in Section 4.

ARTICLE IN PRESS 1260

Z. Li et al. / Signal Processing 90 (2010) 1249–1266

Fig. 9. Segmentation of left ventricle. (a) Result by L-J force driven GAC. (b) Result by GVF snakes.

Fig. 10. Segmentation of the B-mode ultrasonic image of ventricle.

ARTICLE IN PRESS Z. Li et al. / Signal Processing 90 (2010) 1249–1266

1261

Fig. 11. Segmentation of the blob image.

Fig. 12. Segmentation of the horse image using color feature.

diffusion [28] in practice is first used to suppress the speckle noise, and then the L-J force field is computed. Although it is well-known that the B-mode ultrasonic image is not easy to obtain a good segmentation result, the proposed method ends up with a satisfactory ventricle segmentation result. The curve evolution on a microscopy image is shown in Fig. 11. The whole image is very blurry, and there are several cell blobs distributed in the image. Given an initial curve randomly, the L-J force driven GAC automatically handles the typological change of the curve, and finally succeeds in extracting the boundaries of the cells.

and a flower image. On the horse image and the squirrel image, we simply use the color features. We use the CIE L a b color space, for it is preferred over RGB space for visual uniformity [29]. Each pixel ei is represented by a 3vector of fL ; a ; b g components. The flower image has a textural background, in which three red flowers distributed over the white-border leaves. Only by color, it is difficult to distinguish the flowers from the textural background. For this image, we use the texture features. The six-dimensional Carson’s texture descriptor6 [29] is used. Fig. 12 shows the segmentation result of the horse color image. The horse is basically rufous except tail, and

5.2. Group two In this group, we test the proposed method on three color images including a horse image, a squirrel image,

6 The complete Carson’s texture descriptor is eight-dimensional, and the last two-dimensional feature is position vector s to guarantee the clustering coherence. In the experiment, we drop out these two terms.

ARTICLE IN PRESS 1262

Z. Li et al. / Signal Processing 90 (2010) 1249–1266

Fig. 13. Segmentation of the squirrel image.

Fig. 14. Flower segmentation with textural background.

the grass is green. Due to the color distinction, the horse body is well segmented, while the horse tail is preserved. Segmentation of the squirrel image is shown in Fig. 13. The tail of the squirrel is hairy and puffy, so it is hard for the conventional gradient-like methods to get a good performance. The proposed method finely extracts the complete contour of the squirrel. Fig. 14 shows the results of the flower image, which has the rich textural background. We use the Carson’s texture descriptor to distinguish the flowers from the textural background. It is found that the red flowers are extracted successfully from the textural leaves. Moreover, the topological change happens in the textural background. 5.3. Group three In this group experiment, we aim at integrating the statistical information into the L-J force field. Inspired by [30,31], we assume each pixel with a binary label fFrg; Bkgg standing for the foreground and the background, respectively. For each pixel, the log-likelihoods of labeling the foreground and the background are computed, respectively, log PrðI jFrgÞ, log PrðI jBkgÞ. Then, we use these two log-likelihoods to construct F as F ¼ log PrðI jFrgÞ  log PrðI jBkgÞ:

ð16Þ

Several methods can be used to compute the loglikelihood of the foreground and the background. An easy option is using Gaussian mixture model (GMM) to get the likelihoods [19,32]. Some unsupervised techniques such as GMM with EM algorithm [29], Dirichlet process [33], factor analysis [34] can do the job too, although these techniques are relatively complicated. In this paper, we investigate the unsupervised method on the Cygnus loop image and the supervised method on the natural images. 5.3.1. Cygnus loop image segmentation Fig. 15 is the gamma-band image of the Cygnus loop. On this image, we test the unsupervised method to integrate the statistical information. We first compute the histogram of the image, and then an automatic Otsu thresholding is used to separate the histogram into two parts7: one is for the estimate of PrðI jFrgÞ and the other for PrðI jBkgÞ. Fig. 15 shows the curve evolution results. The initial curve is partly outside of the outer loop of the Cygnus. Since the L-J force vectors converge to the boundaries bi-directionally, the curve is attracted to 7 In fact the image is composed of three parts: the outer loop, the inner part of Cygnus and the black background. In the experiment, the grayscale range [0,1] (the whole range of grayscale is [0,255]) is assumed as black background.

ARTICLE IN PRESS Z. Li et al. / Signal Processing 90 (2010) 1249–1266

1263

Fig. 15. Unsupervised segmentation of the Cygnus loop using histogram information.

Fig. 16. Semi-supervised image segmentation of zebra using statistical information. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

the boundaries wherever it is outside of or inside of the Cygnus loop. 5.3.2. Natural images with complex textures Here, we will show the experiments on the natural images with complex textures, in which we use a semisupervised method to compute the likelihoods. We model the foreground and background as single Gaussian distribution,8 GðtjRFrg ; lFrg Þ and GðtjRBkg ; lBkg Þ, respectively, where t is the six-dimensional Carson’s texture 8 In the paper we use single Gaussian to model the foreground and background for simplicity. In practice if background contains multiple different textures, single Gaussian assumption is insufficient. At this situation the background can be modeled using Gaussian mixture model to accommodate the varieties brought by the multiple different textures.

descriptor. Similar to [31], some pixels in the foreground and background are labeled for learning the parameters of Gaussian, i.e., the covariance matrix R and the mean l. Then for each pixel the foreground and background energy, i.e., the likelihoods, are calculated, respectively. In Fig. 16, two green strokes marked on the zebra body and grass ground, respectively, are used to estimate R and l. The feature set F is computed according to (16). Although the zebra has the rich textures. The proposed method accomplishes extracting the whole contour of the zebra well. Fig. 17 shows the segmentation of a tiger in the weed. Two green marked regions, respectively, on the tiger body and the background are used to estimate the parameters of the Gaussian models. We can see that the tiger is perfectly extracted from the clutter background. Fig. 18 shows the segmentation results of

ARTICLE IN PRESS 1264

Z. Li et al. / Signal Processing 90 (2010) 1249–1266

the leaves from the background of grapes. The background is textural, and the color difference between the leaves and the grapes is not distinct. By integrating the statistical information, the proposed GAC well segments the leaves from the grapes.

5.4. Text line detection As mentioned in Section 3.1, the proposed method can produce the morphological dilation effect when ca0. Although [35] (cf. Chapter 4.16) presented to use GAC to

Fig. 17. Semi-supervised image segmentation of tiger with natural textures using statistical information.

Fig. 18. Semi-supervised image segmentation using statistical information.

Fig. 19. Morphological dilation effect used for text line extraction in text layout analysis.

ARTICLE IN PRESS Z. Li et al. / Signal Processing 90 (2010) 1249–1266

1265

Fig. 20. Text line detection in the Chinese handwritten documents.

get the effect of the morphological dilation, our method is more robust and easy to implement. We demonstrate this characteristic with the experiments of the text line detection. In document layout analysis, the text line detection is an important pre-process step for optical character recognition (OCR). Especially for handwritten document analysis, the text lines are usually curvilinear and irregular, so it is difficult to deal with text line extraction in handwritten documents. In [36], GAC was introduced to extract text lines, but GAC was actually used as a postprocess and was followed by a conventional text line detection method. We use the proposed GAC to detect the text lines. The effect of dilation is controlled by the equilibrium position r0 . Because r0 is actually determined by the parameter c (r0 ¼ c in this paper), c is equivalent to the structure element of morphological image processing. For example, a constant c represents a structure element of the disk with radius c. In text line detection, the horizontal characters are desired to be glued into one line, while the neighboring characters under different lines should be separated as widely as possible. So at the horizontal direction the dilation effect should be stronger than at the vertical direction. We design an elliptical ‘‘structure element’’ for the proposed method, and c is defined as a polar elliptical equation with

center at si , ab cðsi ; sj Þ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; a2 ð1  cos2 yÞ þ b2 cos2 y

ð17Þ

where y ¼ angleðsj  si Þ, and a, b are the major axis, the minor axis, respectively. Instead of computing angle that involves inverse trigonometric function, we directly compute cosy ¼ sj  si =jsj  si j  ½10. We set a as the double width of the character and b as half the height of the character. We test our method on one English handwritten document and two Chinese handwritten documents. Fig. 19 shows the result on the English handwritten document. Figs. 20(a) and (b) report the results of two Chinese handwritten documents. All these results are very successful. 6. Conclusions In this paper, we proposed the L-J force driven GAC model, which is inspired by the theory of intermolecular interactions. Different from the previous works that are dependent on the pre-computed edge map, the proposed method is directly based on the original image. Moreover, it could integrate various image information effectively. The L-J force field has two different characteristics: (1) the attractive force can form head-to-head vector flows

ARTICLE IN PRESS 1266

Z. Li et al. / Signal Processing 90 (2010) 1249–1266

converging at the boundaries, and it is suitable for image segmentation; (2) the morphological dilation effect can be obtained if the repulsive force is taken into account too. To evaluate the performance of the proposed force field, we compared it with the popular GVF model with two quantitative criteria: the directional uniformity of the vector field and the average amplitude. To our best knowledge, we are the first one to evaluate the force filed in GAC model quantitatively. In addition, we proved that the electrostatic model is actually a special case of our model. We tested the proposed method on various images, and the experimental results are encouraging.

Acknowledgments The authors would like to thank the anonymous referees for their constructive comments. This work is funded by research grants from Natural Sciences Foundation of China under grant No. 60605004, 60833006, 60675003 and 60835002. References [1] M. Kass, A. Witkin, D. Terzopoulos, Snakes: active contour models, Int. J. Comput. Vision 1 (4) (1988) 321–331. [2] A.L. Yuille, D.S. Cohen, P.W. Hallinan, Feature extraction from faces using deformable templates, in: IEEE Proceedings of Conference on Computer Vision and Pattern Recognition, 1989, pp. 104–109. [3] V. Caselles, R. Kimmel, G. Sapiro, Geodesic active contours, Int. J. Comput. Vision 22 (1) (1997) 61–79. [4] R. Malladi, J.A. Sethian, B.C. Vemuri, Shape modeling with front propagation: a level set approach, IEEE Trans. Pattern Anal. Mach. Intell. 17 (2) (1995) 158–175. [5] S. Osher, J.A. Sethian, Fronts propagating with curvature dependent speed: algorithms based on Hamilton–Jacobi formulations, Journal of Computational Physics 79 (1988) 12–49. [6] S. Kichenassamy, A. Kumar, P. Olver, A. Tannenbaum, A. Yezzi, Gradient flows and geometric active contour models, in: IEEE Proceedings of the International Conference on Computer Vision, 1995, pp. 810–815. [7] L.D. Cohen, On active contour models and balloons, Comput. Vision Graphics Image Process. Image Understanding 53 (2) (1991) 211–218. [8] J.A. Sethian, Level Set Methods and Fast Marching Methods, second ed., Cambridge University Press, Cambridge, 2002. [9] C. Xu, J.L. Prince, Snakes, shapes, and gradient vector flow, IEEE Trans. Image Process. 7 (3) (1998) 359–369. [10] X. Xie, M. Mirmehdi, RAGS: region-aided geometric snake, IEEE Trans. Image Process. 13 (5) (2004) 640–652. [11] N. Paragios, O. Mellina-Gottardo, V. Ramesh, Gradient vector flow fast geometric active contours, IEEE Trans. Pattern Anal. Mach. Intell. 26 (3) (2004) 402–407. [12] B. Sumengen, B.S. Manjunath, Edgeflow-driven variational image segmentation: theory and performance evaluation, Technical Report, Department of Electrical and Computer Engineering, University of California, Santa Barbara, August 2005. [13] D. Yuan, S. Lu, Simulated static electric field (SSEF) snakes for deformable models, in: Proceedings of the International Conference on Pattern Recognition, vol. 1, 2002, pp. 83–86.

[14] A.C. Jalba, M.H. Wilkinson, J.B. Roerdink, CPM: a deformable model for shape recovery and segmentation based on charged particles, IEEE Trans. Pattern Anal. Mach. Intell. 26 (10) (2004) 1320–1335. [15] S.D. Zenzo, A note on the gradient of a multi-image, Comput. Vision Graphics Image Process. Graph Mod. Image Process. 33 (1) (1986) 116–125. [16] G. Sapiro, Vector (self) snakes: a geometric framework for color, texture, and multiscale image segmentation, in: IEEE Proceedings of the International Conference on Image Processing, vol. 1, 1996, pp. 817–820. [17] W. Ma, B.S. Manjunath, EdgeFlow: a technique for boundary detection and image segmentation, IEEE Trans. Image Process. 9 (8) (2000) 1375–1388. [18] D. Comaniciu, P. Meer, Mean shift analysis and applications, in: IEEE Proceedings of the International Conference on Computer Vision, vol. 2, 1999, pp. 1197–1203. [19] N. Paragios, R. Deriche, Geodesic active region for supervised texture segmentation, in: IEEE Proceedings of the International Conference on Computer Vision, vol. 2, 1999, pp. 926–932. [20] D. Cremers, M. Rousson, R. Deriche, A review of statistical approaches to level set segmentation: integrating color, texture, motion and shape, Int. J. Comput. Vision 72 (2) (2007) 195–215. [21] A.J. Stone, The Theory of Intermolecular Forces, Oxford University Press, USA, 1997. [22] J.E. Lennard-Jones, Cohesion, Proc. Phys. Soc. 43 (1931) 461–482. [23] D. Sarid, Exploring Scanning Probe Microscopy with MATHEMATICA, second ed, Wiley-Interscience, New York, 2007. [24] D. Tonnesen, Dynamically coupled particle systems for geometric modeling, reconstruction, and animation, Ph.D. Thesis, University of Toronto, 1998. [25] R. Szeliski, D. Tonnesen, D. Terzopoulos, Modeling surfaces of arbitrary topology with dynamic particles, in: IEEE Proceedings of the Conference on Computer Vision and Pattern Recognition, 1993, pp. 82–87. [26] E.P. Xing, A.Y. Ng, M.I. Jordan, S. Russell, Distance metric learning, with application to clustering with side-information, in: Advances in Neural Information Processing Systems, 2003, pp. 505–512. [27] J.D. Anderson Jr., Computational Fluid Dynamics: The Basics with Applications, McGraw-Hill, New York, 1995. [28] P. Perona, J. Malik, Scale-space and edge detection using anisotropic diffusion, IEEE Trans. Pattern Anal. Mach. Intell. 12 (7) (1990) 629–639. [29] C. Carson, S. Belongie, H. Greenspan, J. Malik, Blobworld: image segmentation using expectation–maximization and its application to image querying, IEEE Trans. Pattern Anal. Mach. Intell. 24 (8) (2002) 1026–1038. [30] B. Sumengen, B.S. Manjunath, Graph partitioning active contours (GPAC) for image segmentation, IEEE Trans. Pattern Anal. Mach. Intell. 28 (4) (2005) 509–521. [31] Y.Y. Boykov, M.-P. Jolly, Interactive graph cuts for optimal boundary and region segmentation of objects in N-D images, in: IEEE Proceedings of the International Conference on Computer Vision, 2001, pp. 105–112. [32] N. Paragios, R. Deriche, Geodesic active regions and level set methods for supervised texture segmentation, Int. J. Comput. Vision 46 (3) (2002) 223–247. [33] C.E. Rasmussen, The infinite Gaussian mixture model, in: Advances in Neural Information Processing Systems, vol. 12, 2000, pp. 554–560. [34] M.E. Tipping, C.M. Bishop, Probabilistic principal component analysis, Technical Report, Neural Computing Research Group, Aston University, Birmingham, UK, Technical Report NCRG/97/010, September 1997. [35] A. Bovik (Ed.), Handbook of Image and Video Processing, second ed., Academic Press, New York, 2005. [36] Y. Li, Y. Zheng, D. Doermann, S. Jaeger, Script-independent text line segmentation in freestyle handwritten documents, IEEE Trans. Pattern Anal. Mach. Intell. 30 (8) (2007) 1313–1329.