The line- and block-like structures extraction via ingenious snake

Pattern Recognition Letters 112 (2018) 324–331

Contents lists available at ScienceDirect

Pattern Recognition Letters journal homepage: www.elsevier.com/locate/patrec

The line- and block-like structures extraction via ingenious snake Shoujun Zhou a,1,∗, Baolin Li a,1, Yuanquan Wang b,∗∗, Cheng Wang a, Tiexiang Wen a, Na Li a a b

Shenzhen Institutes of Advanced Technology, Chinese Academy of Sciences, 1068 Xueyuan Avenue, Shenzhen 518055, China School of Computer Science, Hebei University of Technology, No. 5340 Xiping Avenue, Beichen District, Tianjin 300401, China

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 6 March 2018 Available online 16 August 2018

Active contour model (ACM) plays an important role in computer vision and medical image analysis. The traditional ACMs were employed to extract closed contours of objects. While simultaneous extraction of line- and block-like objects, such as boundary contours, centerlines, as well as their topological relationship, remains open so far. Therefore, a novel ACM named "Ingenious Snake" is proposed to adaptively extract the feature curves. The proposed ingenious snake includes the following steps: 1) In the preprocessing, the line- and block-like structures are classified with k-means clustering, following up with morphological operation and enhancement. The gradient vector flow (GVF) field is then acquired from the resultant ridge feature map. 2) For the automatic initialization, the ridge-points are extracted by using the local phase measurement of GVF field, then the two-category of object ridgelines are obtained fast. 3) Finally, the contour deformation and curves evolvement are implemented with a management strategy. The resultant contours and centerlines well characterize the objects of interest. In the experiments, we compare the existing initialization methods and the adaptive extraction with a series of phantoms and testing images. The visual and quantitative assessments of structure extraction are satisfying in terms of effectiveness and accuracy. © 2018 Elsevier B.V. All rights reserved.

1. Introduction Simultaneous extraction of line and block structure as well as their topological relationship is important especially for medical image-guided intervention surgery-planning, which remains open so far. To our known, some related works used the state-of-the-art methods to extract the edge and contour, such as DeepContour [1], DeepEdge [2], and Holistically-Nested edge detection (HED) [3], where novel technologies involved deep convolutional neural networks (CNN), contour partitioning and model parameters fitting, multi-scale bifurcated deep network for top-down contour detection, and an end-to-end edge detection system to learn the type of rich hierarchical features. These data-driven methods have higher precision than the traditional methods (e.g., Sobel, Canny) in complicate edge detection, given that mass data processing be used to understand individual data. As the model-driven methods, Level Set also developed several efficient techniques such as piecewise-smooth model with approximately explicit solutions [4],

∗

Corresponding author. Co-corresponding author. E-mail addresses: [email protected] (S. [email protected] (Y. Wang). 1 Co-first authors. ∗∗

Zhou),

https://doi.org/10.1016/j.patrec.2018.08.018 0167-8655/© 2018 Elsevier B.V. All rights reserved.

[email protected]

(B.

Li),

active contours driven by regularized gradient flux flows [5], intensity inhomogeneity approach [6], simultaneous image segmentation and bias correction [7]. Another model-driven method is the active contour model (ACM, namely, Snake). Since the early debut of Snake [8], there has been a furry of works devoted to its initialization, construction and optimization of the external energy terms [2–20]. The open snake (or named open curve contour model) was early proposed by Wong et al. [12] and employed in [13,14] to extract the actin-network centerline. Different from close snake, the internal force term at the ports of open snake are cancelled, and the stretch forces are exerted on both ends. Generally, the initialization effects on calculation time and iterative approximation process, which is crucial especially for open snake. Many works strived for initializing ideal object curves and contour. Considering the slow evolvement process of the state-of-theart open snake, its initialization is more important than that of the close one. The representative researches are in [14–17]. For example, the initialization of the center of divergence (CoD) method [15] was attained by defining the GVF divergence central points, dealing with many initialization contours, and furthermore deleting the pseudo contours. In the force field segmentation (FFS) method [16], the initial contours were generally away from the object. Thus, it was faced with time-consuming iteration and approaching. For the Poisson inverse gradient (PIG) method [17], a few initialized models were specified in advance and followed with

S. Zhou et al. / Pattern Recognition Letters 112 (2018) 324–331

325

Fig. 1. Overview of the proposed Ingenious Snake for line- and block-like objects extraction. Two categories of objects are separated apart from the background, whose contours are extracted respectively. Symbols I, f, F, P, R and S denote the original image, the ridge feature map, the normalized GVF field, the ridge-point set, the ridgeline set, and the candidate contour set, respectively. The final extraction results of contour and centerline are grouped respectively.

solving Poisson equation to estimate the potential external energy field. Considering evolvement efficiency, open snake initialization is more important than the close snake’s, for which the internal seed points in the interesting objects were generally used [13,14] and the automatic initialization was realized by seed points extraction only along the natural coordinate axis, thus resulted in scattered and discontinuous initialization curves, especially in low contrast regions. In short, most of the relevant methods have a common deficiency, namely, topology connection between the region edges and tubular structure centerlines, which is important for path navigation, especially for the medical AI-assisted interventional surgery and the control of unmanned aerial vehicle. In this paper, a new method is proposed to fast and adaptively extract the two-category objects. The contributions are as follows: 1) Preprocessing of ambiguous objects uses the k-means clustering, following up with morphological operator and multi-scale filtering to preliminarily separate out the line- and block-like objects, further to obtain their ridge feature maps (RFM). 2) Based on the RFM, contours initialization is proceeding with isotropic measure (i.e., the local phase characteristic of the GVF field) and ridge-points tracking. The resultant curves possess sufficient ridge information. 3) A novel active contour model, namely, ingenious snake, is proposed to extract open- or closed-ended contours relating to line- and block-like objects respectively, where an effective management strategy is used for the snake evolvement.

2. Method A diagram of the proposed method is shown in Fig. 1. The ingenious snake performs the deformation and evolvement. First, the RFM is obtained by preprocessing the original image. Then the ridgelines are obtained by the automatic initialization algorithm. Finally, the contours of interesting is obtained by using the management strategy during evolvement.

2.1. Ingenious snake As a special application of ACM, ingenious snake is expressed as r (s ) = (x(s ), y(s )) where s represents the normalized arc-length with values s ∈ [0, 1]. The total energy Esnake consists of the internal energy Eint and the external energy Eext . The internal energy is like

the traditional snake definition:

 Eint =

1 0



 α (s )|rs (s )|2 + β (s )|rss (s )|2 ds

(1)

Where rs (s) and rss (s) stand for the first and second derivatives of the snake point (named "snaxel" in the follow-up contents) r(s) on the curve. The coefficients of elasticity and stiffness of each point on the deformation curve are expressed by α and β respectively, and different values of α and β can be obtained for the snaxels. Traditionally, the two coefficients are set as constant for the close snake contour [8] and be zero at the both tips for the open snake curve [12,13]. To express both the open curve and the closed contour, the external energy term of ingenious snake Eext is com posed of image term Eimg and stretch force term Estr :

 Eext =

1 0

kimg · Eimg (r (s ) ) + ζ · kstr · Estr (r (s ) )ds

(2)

The weight parameters kimg and kstr are respectively used to control the scale of the two external forces. In Eq. (2), the logical variable ζ , with ζ = 0 or 1, differentiates it from the traditional ACMs. Automatic identification of the port-type (namely, closedor open-ended curves) is as follows. If ζ = 0, Eq. (2) is reduced to describe the evolvement of closed contours. The form is identical to its traditional one, where the image energy kimg · Eimg (r(s)) acts as the external energy of the contour. The corresponding external force is constructed based on GVF field, that is, −∇ IGV F , as the similar one in [13]. If ζ = 1, Eq. (2) denotes an open snake model, which the stretch forces begin to take effect. The stretch forces are applied to both ends of the open curve and point to the tangent direction of the front-points (namely, the two end points). According to [18], the size of the stretching force is defined as Fstr (r (s ) ) = (I (r (s ) ) − Imean )/(I f − Ib )), where I(r(s)) denotes the pixel value covered by a snaxel r(s), while If , Ib , Imean denote the local intensities of the foreground, the background, and their mean, respectively. The differential of the stretch force is expressed as follows:

⎧ r s s( ) ⎪ ⎪ ⎨− |rs (s )| · Fstr (r (s )), s = 0, rs ( s ) ∇ Estr (r (s )) = · Fstr (r (s ) ), s = 1, ⎪ ⎪ | ⎩ rs (s )| 0 otherwise.

(3)

326

S. Zhou et al. / Pattern Recognition Letters 112 (2018) 324–331

the state-of-the-art technologies such as data-driven methods [1– 3] and model-driven methods [4–7]. Considering the limits of mass data samples and the problem of ambiguous object, we use the classical k-means clustering for the moment, following up with morphological operator and multi-scale filtering to obtain RFM. The process is specified as follows:

Fig. 2. Structural ridge enhancement: (a) the original image of sapling structure; (b) initial mask M is obtained by K-means clustering; (c) block-mask Mblo ; (d) linemask Mlin ; (e) the f1 from line-mask; (f) the f0 from block-mask.

When initial curves are deforming and evolving using the Eqs. (1)–(2), the ingenious snake activates an adaptive coefficient setting of α , β and the logical variable ζ in terms of their porttype. Let Ri be the i-th initial curve, Ri = {X, ζ }i , where X is the snaxel sequence composed of N points xn (with n = 1, . . . , N) and ζ is the logic variable of the curve port-type. For the constant coefficients α and β of the two front-points in Eq. (1), we modulate their values by using the window functions ω(n), namely



ω (n ) =

0, n = 1, 2, N − 1, N 1, others

(4)

Thus, the α , β in Eq. (1) are expressed as:

α (n ) = α0 · ω (n ) · ζ + α0 · (1 − ζ ) β (n ) = β0 · ω (n ) · ζ + β0 · (1 − ζ )

(6)

where X is the coordinate vector of the snaxels sequence ( x1 . . . ,xN )T with xi = (xi , yi ), 1 ≤ i ≤ N. At the t-th iteration, Xt can be computed by the Euler-Lagrangian equation −1

X = (A + γ I ) t

 =M◦k = Mk k Mblo ( r r) r

 kimg ∂ Eimg (X t−1 ) kstr ∂ Estr (X t−1 ) t−1 · γX − −ζ · ∂X ∂X (7)

where A is composed of coefficients α and β , γ is used to control the step size, and I is the unit matrix. In condition of ζ = 0, it can approach to the object and make the Eq. (7) convergence with few iteration-steps (i.e., nmax = 3). In condition of ζ = 1, it continues evolving until that the stretch force disappears or the displacement of its two front-points is less than a low limit (i.e., rˆmin ). 2.2. Preprocessing The original image I contains the region boundaries and the line structures, see Fig. 2, the preprocessing process is used to enhance the ridge characteristics of the two categories of objects. The enhanced image f0 and f1 (namely, the RFMs of the two categories) reflect ridge characteristics of boundaries and medialness w.r.t. block- and line-like structures. In general, the edge can be extracted using the Sobel operator or the Canny operator [20], or

(8)

where opening is the dilation of the erosion of the M by the structuring element kr ,  and  denote erosion and dilation, respectively. The erosion operation removes structures that are smaller than kr and dilation operation restores the shape of remaining  is generated by subtracting the mask structures. The line-mask Mlin  from mask M, namely M = M − M . Furthermore, an intact Mblo lin blo line-mask Mlin (e.g., Fig. 2(d)) is obtained by removing the small  . noise blocks with maximizing region connection in the mask Mlin Finally, we can obtain an intact mask Mblo (e.g., Fig. 2(c)) by subtracting the line object Mlin from the original mask M, namely

Mblo = M − Mlin (5)

where n = 1, . . . , N. For the external energy of Eq. (2), the discrete external energy Eext (X):

Eext (X ) = kimg · Eimg (X ) + ζ · (kstr · Estr (X )

2.2.1. Generation of category mask See Fig. 2 for an example, we use K-means clustering and morphological operators to extract the category-mask. Firstly, based on pixel gray distribution of the background and the foreground, Kmeans clustering algorithm aggregates the image pixels into two categories, in which the foreground with the smaller cluster center involves three objects (namely, leaves, root, and stem, see Fig. 2(b)), thus an initial mask M is obtained and contains these objects. Secondly, two morphological operators with erosion and dilation kernel kr are successively performed to change the mask. As the result, the thin line structure (stem) can be eliminated from the mask, while the block structures (leaves and root) are shrunken and regain their region. See the result in Fig. 2(c). The process runs by the morphology opening operation as follows.

(9)

Note, the Mlin and Mblo are the resultant line- and block-mask regions, denoted by f1 and f0 respectively. 2.2.2. Ridge enhancement The block-region boundary and line-structure center can be enhanced respectively. The resultant RFM is used to construct the external force field for ACM. The RFM of boundary is enhanced by using Canny edge extraction operator following with Gaussian kernel Gσ filter; while the RFM of line-structure center is enhanced by using the multi-scale filter. As shown in Fig. 2, the subfigures (e) and (f) are the f1 and f0 from (c) and (d), respectively. The resultant RFMs will be used for the initialization according to the following isotropic ridge-points extracting and tracking. 2.3. Automated initialization The GVF field is normalized and denoted by F(x, y), which can be obtained by solving a Euler iteration diffusion equation [9]. Firstly, the GVF field is obtained through computing the RFM. Then, the ridge-points are obtained by the isotropic direction measure of (i, j ) (i, j ) the GVF field. Given the two-dimensional vector Fx and Fy at the point p(i, j), it will be defined as a ridge-point if they satisfy one of the following conditions



i−1, j ) (i+1, j ) i−1, j ) (i+1, j ) i, j Fx( Fx + Fy( Fy ≤ 0 i f Fx( ) > 0 (i, j−1 ) (i, j+1 ) (i, j−1 ) (i, j+1 ) (i, j ) Fx Fx + Fy Fy ≤ 0 i f Fy > 0

(10) (i, j )

Image-marginal noise point can be reduced by making Fx (i, j ) and Fy greater than zero. For example, in Fig. 3, the 2D coronary angiography X-ray image is shown in (a), subfigure (b) shows

S. Zhou et al. / Pattern Recognition Letters 112 (2018) 324–331

327

Algorithm 1 The pseudo code of the ingenious snake.

Fig. 3. An example of ridge-point extraction from a coronary radiography image: (a) original image; (b) the GVF field is calculated and displayed overlapping with the RFM; (c) ridge-point extraction results; (d) isolated point cluster.

•

Fig. 4. Generating a ridgeline using ridge-point tracking, where two black points are the front-points, where the candidates locate in the fan shaped region.

the GVF field overlapping with the RFM, wherein the local amplification shows that the GVF vector points to the vessel center. The extraction of isotropic ridge-points via Eq. (10) is shown in (c) where most of the ridge-points are located at the vessel center. Few isolated points are derived from artifacts and local noises, as shown in Fig. 3(d) enlarged from Fig. 3(c). The isolated points can be automatically deleted according to a criterion that a ridge-point xi has no neighbour ridge-points xj with a prior range rˆmin < d|x j −xi | < rˆmax that is an observation from the minimal and maximal radius of vessel objects. Next, we provide a fast algorithm to automate a ridge-points tracking process which uses a local measure on the distances and the connecting vectors from the current point to the neighboring ridge-points. The steps are detailed as follows. •

•

Starting a tracking process using seed point: See Fig. 4, a seed point p0 is randomly generated from the ridge-point set P, when the neighboring two points p−1 and p1 are detected with the angle of the two direction vectors being no less than a big obtuse angle (e.g., 2π /3). Otherwise, the seed point is removed from the ridge-point set. Ridge-point tracking: Given seed point p0 , the tracking process continues to search the candidates along the path of opposite directions. Suppose that vt and v−t represent two reverse tracking directions, pt and p−t represent two front-points of the tracking. The tracking process sequentially update p ± t ± 1 , v ± t ± 1 to form a ridgeline in X. In current tracking step, the ridge-points around the path in the last step should be deleted

(namely, the deleted ridge-points are the ones whose distance to the ridgelines is less than rˆmax ), thereby updating the ridgepoints set P. Termination criteria: The tracking process will stop at the termination criteria, they are: (1) no candidates can be used to generate the front-points pt and p−t ; (2) after a length of tracking, the front-distance dɛ between the left and right front-points is less than a low-limit ɛ0 , namely dɛ < ɛ0 (usually ε0 = rˆmin ); (3) there is not seed points to start a new tracking when the whole ridge set P is empty, then the tracking process is finished.

Given L is the total number, we use R = {R1 , . . . , RL } representing the set of the ridgelines. Each ridgeline Ri = {X, ζ }i contains the initial coordinate sequence X = [x1 , x2 , x3 , . . .] and curves’ porttype ζ . When the termination terms fall into the second criteria, a closed contour is determined and labelled with ζ =0 in terms of dɛ < ɛ0 ; otherwise, if the termination terms fall into the first criteria, the curve will be determined as an open curve labelled with ζ =1. The resulting ridgelines are used as the initial contours for ingenious snake in the Section 2.1. 2.4. Management strategy of curves evolution Ridgelines contain closed contours and open curves. They are evolving by ingenious snake to attain convergence. The pseudocode of Ingenious Snake evolution is shown in the Algorithm 1. During the evolvement, the interactions such as curve merging, overlapping, and crossover would be managed. Meanwhile, inhomogeneous curves and contour (namely, the line structures versus the block ones) should be discriminatingly characterized. Therefore, a management strategy is used for the evolvement process. Firstly, all the closed contours R0 deforms in regular succession. Then the open curves R1 will be evolved one by one, wherein curves overlap and intersection (among open curves and closed contours) might occurs during the open curves’ evolvement or after their convergence. To facilitate the expression, the t-th step iterating evolvement of n-th curve is denoted by Stn , whose frontpoints are xt1 and xtN ; another converged curve is denoted by Sm with m = n. Evolvement maybe generate a connection between Stn and Sm . For connecting between two open curves, our program

328

S. Zhou et al. / Pattern Recognition Letters 112 (2018) 324–331

Fig. 5. The original images which include single-category and two-categories objects: (a) CT spine image with size of 225 × 415; (b) the coronary X-ray angiography with size of 400 × 400; (c) the fundus image with size of 328 × 329; (d) sapling image with size of 234 × 284; (e) ant phantom with size of 295 × 429; (f) abdomen MR angiography with size of 238 × 240.

permit Stn merging with one part of Sm to produce a complete branch, which is favourable to form a topological structure. Note, there are four connecting types for snake curves merging from the left or right side of current Sn to the left or right part of object Sm , where the corresponding operation is based on following two measures, named the natural smooth measures T1 and T2. −−→ T1: Given the candidate xˆ j ∈ Sm and the vector xti xˆ j connects t the front points xi to the snaxels xˆ j , the current snake curve Snt merges with the object curve Sm in condition that the intersection −−→ −−→ angle vti · xti xˆ j /d j ≤ π2 and the distances dtj = |xti xˆ j | ≤ rˆmax . T2: An optimum connecting point xˆ tj (i.e., the merging point) is −−→ located at max(|vm (xˆ j ) · xti xˆ j /dtj | ) where vm (xˆ j ) is the axial vector T1

at xˆ j in Sm , thus it indicates that the optimal connecting point xˆ tj is the one displaying the most consistent direction between front −−→ vector vti and candidates’ vector xti xˆ j . Note that our program prevents the Snake connection between the line- and block-structure. After all the evolvements finished, there comes into being lot of redundancies, i.e., curve-overlapping. To remove them, curve-overlapping check between Sn and Sm is conducted by the brute-force strategy according to [13], where Euclidean distances from each snaxels in Stn to Sm were calculated. 3. Experiment To validate the proposed method on extracting various object with single- or two-category contour. A set of 2D phantoms and medical images were employed in the experiments. Part of these data included the interesting objects such as block-like spine in CT image, line-like vessel in X-ray coronary and fundus angiographies; see Fig. 5(a)–(c). Another part of data included line- and block-like objects such as sapling phantom, ant phantom, and the aortaventralis magnetic resonance (MR) angiography; see Fig. 5(d)–(f). The proposed method is compared with the center of divergence (CoD) method [15], force field segmentation (FFS) method [16], Poisson inverse gradient (PIG) method [17] and stretching open active contours (SOACs) [14]. The traditional CoD, FFS, and PIG method were used to extract the close contour, while SOACs was used to extract the open contour. We will validate the ingenious snake extract both the line- and block-like objects with the close and open contours, based on visual and quantitative evaluation. Ground truths of Phantoms were provided manually for quantitative evaluation. The workstation is configurated as follows, Intel (R) Core (TM) i76700 CPU @ 3.41 GHz, RAM 32 G, operating system Win10 and programming language MATLAB r2014a. 3.1. Parameter setting In the experiments, the parameters of the traditional methods (CoD, FFS, PIG and SOACs) were set to the default according to [14∼17]. The main parameters were set for the proposed method,

Fig. 6. Single-category contour extraction from Fig. 5(a). From up-side to downside, show the results of initialization and convergence; from left-side to rightside, show the contour extraction results by the CoD, FFS, PIG, SOACs and proposed method.

see Table 1, which involved preprocessing, deforming, and evolvement management. The parameters for deformation and evolvement were set uniformly for the same type of datasets. With respect to the twocategory contour extraction from line- and block-like objects in Fig. 5(d)–(f), the centers of K-means clustering were set according to the number of structures and backgrounds. The morphological kernel kr uses a circular structure element with the scale approximate to rˆmax of the line structure. Meanwhile, the automated initialization also uses the scale parameters rˆmin and rˆmax in terms of the width range of interesting ridges, which was used to delete the isolated ridge-points. The low-limit of fronts’ distance ɛ0 is uniformly set with rˆmin . In the curves evolvement, parameters α , β , γ , kimg , kstr and nmax are also assigned uniformly. Meanwhile, curves are resampled to maintain point spacing strictly with unit pixel. In addition, curve-overlapping check was conducted by the bruteforce strategy [9], which involved calculating the Euclidean distance from each point of a curve to all the other converged curves. If the maximal distance to a converged curve was smaller than a threshold Dmin = rˆmin , the curve-overlap was detected. 3.2. Single-category contour extraction We compared different methods to extract single-category contour from block-like objects, such as the spinal boundaries in Fig. 5(a), whose extraction results are shown in Fig. 6. The quantitative comparison of these methods is shown in Table 2. During the automatic initialization, the CoD, FFS and proposed method produce 60, 210 and 7 closed contours, respectively; the PIG method preset 7 closed contours for initialization; the SOACs method generated 168 initial open curves. For iterative evolvement, the proposed method can adaptively evolve the close- and openended initial contours. While the CoD, FFS and PIG methods only performed the deformation of the close-ended contours. The proposed method was faster than the CoD and FFS method, for example in Table 2, the iterations and time-consuming of ours was 21 and 3.72 s, respectively. Although the PIG method is slightly faster than ours, it requires a priori number of the initial contours. The SOACs proceeds only with open snake evolvement. The iterations and time-consuming of it is 131 and 4.46 s for SOACs. The proposed method performed well in terms of efficiency and effectiveness. We also used the SOACs method and the proposed method to extract central lines in a set of coronary angiography image; see

S. Zhou et al. / Pattern Recognition Letters 112 (2018) 324–331

329

Table 1 Parameters settings of the proposed method for the six images shown in Fig. 5. The intensities of all images are scaled to be in the range [0,1]. Parameters settings

Fig. 5(a)

Fig. 5(b)

Fig. 5(c)

Fig. 5(d)

Fig. 5(e)

Fig. 5(f)

Num. of categories kr rˆmin rˆmax Evolv. Manag. Deformation

— — 2.0 4.5 ε0 = rˆmin , Dmin = rˆmin , nmax = 3 α = 0.1, β = 0.5, γ = 2, kimg = 2, kstr = 3

— — 2.0 4.5

— — 1.5 3.5

2 5∗ 5 2.0 5.0

2 5∗ 5 2.0 5.0

3 5∗ 5 3.0 8.0

Table 2 Comparing the computational efficiency of spinal boundaries extraction from CT spine image in Fig. 5(a). Method

Evolu. type

Init. quantity

Total iter. num.

Time consu. (s): ini./evo./total

Ave. MSE (pixels)

CoD FFS PIG SOACs Proposed

close close close open adaptive

60 210 7 168 7

336 450 75 131 21

2.12/14.94/17.06 8.66/12.44/21.10 1.29/2.29/3.58 3.21/1.25/4.46 3.01/0.71/3.72

13.2 2.7 3.1 2.3 2.2

Fig. 8. Two-category contour extraction from Fig. 5(d) and (e): (a1, a2) the original images, (b1, b2) initial curves and contours of the line- and block-like objects, (c1, c2) the evolvement results by the ingenious snake.

Fig. 7. Single-category contour extraction from Fig. 5(b),(c) by SOAC and the proposed method. From up-side to down-side, show the results of initialization and convergence.

the left four subfigures of Fig. 7, the SOAC’s and the proposed method generated 646 and 72 initial ridgelines, and eventually converged to 171 and 24 contours, respectively. In addition, complete central lines were obtained by the proposed method, according to their convergence results. For the fundus image in the right four subfigures of Fig. 7, the SOACs and the proposed method generated 440 and 83 initial ridgelines, and eventually converged to 101 and 59 contours, respectively. The proposed method greatly restrained the noise and thus avoided many fragments, according to the initialization in the four subfigures upside of Fig. 7. Furthermore, the vascular centerlines were extracted with more complete topology by the proposed model than that by SOAC’s. 3.3. Two-category contour extraction A set of images with line- and block-like structures, e.g., in Fig. 5(d)–(f), were used in the experiments to validate the capability of ingenious snake in two-category contour extraction, which cannot be realized by using the existing methods. As shown in Fig. 8, the subfigures (a1) and (a2) are the phantoms with Gaussian noise. The initial ridgelines were generated by automatic ridgepoint tracking and shown by overlapping on the phantoms, see the subfigures (b1) and (b2). We show their evolvement results by the ingenious snake in subfigures (c1) and (c2). These results quantitatively described their topology pattern and relationship. The extraction results were used for quantitative analysis. For the centerlines and boundary contours of interest, their ground truths were manually labeled. For the segmentation of the block

regions, i.e., the leaves from the left and right side, root block, and ant body, the resultant Dice coefficient is 0.9751, 0.9751, 0.9510, and 0.9882, respectively. To evaluate the precision of the extracted centerlines (i.e., the sapling branches and ant legs) and contours (i.e., left and right side, root block, and ant body) in Fig. 8, we refer to the document [9] and use the Vertex error (VE) and Hausdorff distance (HD) to evaluate the extraction results. Assume Sgt represents the center line gold standard, Sc represents the extracted centerline, and VE is defined as:

dV (Sc , Sgt ) =

1 min ||xc − xgt xgt ∈Sgt 2|Sc | x ∈S

c c



1



×

+ min ||xgt − xc

2|Sgt | xgt ∈Sgt xc ∈Sc

(11)

The HD definitions of Sgt and Sc are defined as:





dH (Sc , Sgt ) = max max min

xc inSc xgt ∈Sgt

||xc − xgt ||, max min ||xgt − xc || xgt inSgt xc ∈Sc

(12) For the contour extraction results by the ingenious snake in Fig. 8, the VEs and HDs were on average 0.528 and 2.304 pixels respectively, see Table 3. In comparison, SOACs method also obtained relevant VEs and HDs, i.e., on average 0.588 and 2.453 pixels respectively. In short, all the traditional methods could not produce an intact object extraction from line- and block-like objects. We used different methods to extract the two-category contour in the MR image of Fig. 9. In comparison, both the organ’s boundary and vessels’ central lines are the interesting contours. The CoD method initialized 42 small circles and distributed randomly inside the objects and background. The resultant convergence results could not well express either the vessels or the organ. The FFS method produced the initialization results outside the objects and converged at the boundaries of vessels and part of the organ, thus

330

S. Zhou et al. / Pattern Recognition Letters 112 (2018) 324–331 Table 3 The VEs and HDs of the extraction results from image in Fig. 8(c1) and (c2) are compared between the SOACs’ and the Ingenious snake. Method and phantoms in Fig. 8 SOACs [9]

(c1) (c2)

Ingenious snake

(c1) (c2)

Objects

Contour category

VEs

Mean VEs

HDs

Mean HDs

branches leaves & root ant legs ant body

open close open close Mean: open close open close Mean:

0.639 0.732 0.501 0.477

0.686

3.108 2.291 2.200 2.391

2.610

branches leaves & root ant legs ant body

Fig. 9. Two-category contour extraction from Fig. 5(e): The subfigures from left to right show the results from CoD, FFS, PIG, SOACs and proposed method. The subfigures in the first and second rows show the results of the initialization and the corresponding convergent results respectively.

it couldn’t distinguish between the two objects. The PIG’s only produced initialization results within the organ, where completeness description was lacked for both objects. The SOACs generated 212 initial ridgelines and eventually converged to 42 contours, where discontinuous and inaccurate results occurred at organ and vessels. The proposed method obtained the two-category contour at the same time. 4. Discussion and conclusion To extract the ambiguous object from the appointed datasets, model- and data-driven methods have their respective characteristics and advantages. The former concerns parametric settings and man-machine interaction, while the latter closely relates homogeneous mass data and the problem of ambiguous object and sample features. We proposed a novel parametric ACM, i.e., Ingenious Snake, to extract the line- and block-like objects together. While the traditional ACMs can only extract either of them, thus fail to describe the relationship between spatial topology and region, e.g., the structures in Fig. 5(d)–(e) and the topology relationship of vessels and organs in Fig. 5(f). For the resultant initialization, deformation, and evolvement management, there existed great difference between the traditional ACMs and the proposed method. Owing to the length limitation of the article, just six representative images were shown in Fig. 5. Visual and quantitative evaluation of the experimental data gave the proposed method best proofs. Table 3 showed that the VEs and HDs from ingenious snake were lower than that from SOACs. The evolution effect powerfully validated the ingenious snake is superior to SOACs. Besides well possessing the accuracy, the ingenious snake can efficiently extract the line- and block-like structure, further to describe the topology relationship between the two categories of them. These advantages own to pre-classification, automatic initialization, and the evolvement management strategy. Note, in Fig. 7, it is difficult to

0.580 0.605 0.498 0.425

0.489 0.588 0.593 0.462 0.528

2.732 2.136 2.190 2.155

2.296 2.453 2.434 2.173 2.304

thoroughly suppress the background noise, which produced finite noise ridge-points and the pseudo vessels ridgelines. To effectively compare the coronary vessels extraction in Fig. 7, our program shielded the ridgelines unconnected to the root (a manual label point) of the coronary artery tree. This operation is mainly to describe the completeness of the topological structures extraction, not to show the characteristics of the Ingenious snake. The parametric settings include the default Snake viscoelastic α , β and the update speed γ , tension forces’ kimg and kstr along radial and axial directions respectively. Without loss of generality, the relation 1 ≤ γ ≤ kimg ≤ kstr ≤ 3 can be adaptive to all the experimental images. Assumed that sampling snaxels interval be unit pixel, the above-mentioned parameter settings are effective. Another affect would be related to the sizes of image and objects, only concerned with the iteration step length and evolvement time. The minimal and maximal radius values rmin and rmax describe a prior range of the size of vessel-of-interest in the dataset including vascular structures. We may as well assign the lower limit by rˆmin = 1 (near the rmin ). The upper limit rˆmax (near the rmax ) is assigned with an empirical value according to visual observation based on multi-scale filtering with scale rˆmin ≤ σi ≤ rˆmax . Given the resultant RFM clearly include all the vessel-of-interest, the rˆmin and rˆmax are applicable. Because the ACM can be integrated with the results from other popular methods, our next work is to combine the ingenious snake with these state-of-the-art methods, to quantitively describe the structure topologies as well as their function relations in complicate image context. For the relevant application in AI-assisted intervention surgical plan and navigation, we will also solve the problems associating with the three-dimensional topology structure extraction, where deep learning methods will be greatly considered in research. Conflict of interest Where deep learning methods with respect to the topological contours extraction will be greatly considered in research. Acknowledgements This work is supported by the following grants: National High Technology Research and Development Program of China (National 863 Program) (No. 2015AA043203), National Science Foundation Program of China (No. 61471349), Basic Discipline Layout Project of Shenzhen City (No. JCYJ20150731154850923), and Shenzhen Engineering Laboratory for Key Technologies on Intervention Diagnosis and Treatment Integration. References [1] W. Shen, X. Wang, Y. Wang, X. Bai, Z. Zhang, Deepcontour: A deep convolutional feature learned by positivesharing loss for contour detection draft version, IEEE Conf. Comp. Vision Pattern Recog. (2015) 3982–3991.

S. Zhou et al. / Pattern Recognition Letters 112 (2018) 324–331 [2] G. Bertasius, J. Shi, L. Torresani, Deepedge: A multi-scale bifurcated deep network for top-down contour detection, IEEE Conf. Comp. Vision Pattern Recog. (2015) 4380–4389. [3] S. Xie, Z. Tu., Holistically-nested edge detection, in: IEEE International Conference on Computer Vision, 2015, pp. 1395–1403. [4] H. Song, Y. Zheng, K. Zhang, An efficient algorithm for piecewise-smooth model with approximately explicit solutions, Electron. Lett. 53 (4) (2017) 233–235. [5] H. Song, Active contours driven by regularised gradient flux flows for image segmentation, Electron. Lett. 50 (14) (2014) 992–994. [6] K. Zhang, L. Zhang, K.-M. Lam, D. Zhang, A level set approach to image segmentation with intensity inhomogeneity, IEEE Trans. Cybernet. 46 (2) (2016) 546–557. [7] K. Zhang, Q. Liu, H. Song, X. Li, A variational approach to simultaneous image segmentation and bias correction, IEEE Trans. Cybernet. 45 (8) (2015) 1426–1437. [8] M. Kass, A. Witkin, D. Terzopoulos, Snakes: active contour models, Int. J. Comput. Vision 1 (4) (1988) 321–331. [9] D.W. Ren, W.M. Zuo, X.F. Zhao, Fast gradient vector flow computation based on augmented Lagrangian method, Pattern Recognit. Lett. 34 (2) (2013) 219–225. [10] S. Zhang, J. Zhou, Centerline extraction for image segmentation using gradient and direction vector flow active contours, J. Signal Inform. Process. 4 (2013) 407–413. [11] Y. Wu, Y. Wang, Y. Jia, Adaptive diffusion flow active contours for image segmentation, Comp. Vison Image Understand. 117 (10) (2013) 1421–1435.

331

[12] Y.Y. Wong, P.C. Yuen, C.S. Tong, Segmented snake for contour detection, Pattern Recognit. 31 (11) (1998) 1669–1679. [13] T. Xu, D. Vavylonis, X. Huang, 3D actin network centerline extraction with multiple active contours, Med. Image Anal. 18 (2014) 272–284. [14] T. Xu, D. Vavylonis, F.C. Tsai, SOAX: a software for quantification of 3D biopolymer networks, Sci. Rep. 5 (March 2015) 9081, doi:10.1038/srep0908. [15] X. Ge, T. Jie, An automatic active contour model for multiple objects, in: proceedings of IEEE International Conference on Pattern Recognition, 2002, pp. 881–884. [16] C. Li, J. Liu, M.D. Fox, Segmentation of external force field for automatic initialization and splitting of snakes, Pattern Recognit. 38 (11) (2005) 1947–1960. [17] B. Li, S.T. Acton, Automatic active model initialization via Poisson inverse gradient, IEEE Trans. Image Process. 17 (8) (2008) 1406–1420. [18] Y. Wang, A. Narayanaswamy, C. Tsai, B. Roysam, A broadly applicable 3-D neuron tracing method based on open-curve snake, Neuroinformatics 9 (2) (2011) 193–217. [19] M. Abozahhad, R.R. Gharieb, S.M. Ahmed, Edge detection with a preprocessing approach, J. Signal Inform. Process. 5 (4) (2014) 123–134. [20] S. Chack, P. Sharma, An improved region based active contour model for medical image segmentation, Intern. J. Signal Process. Image Process. Pattern Recog. 8 (2015).