Exact solution to median surface problem using 3D graph search and application to parameter space exploration

Pattern Recognition 48 (2015) 380–390

Contents lists available at ScienceDirect

Pattern Recognition journal homepage: www.elsevier.com/locate/pr

Exact solution to median surface problem using 3D graph search and application to parameter space exploration Zhengwang Wu a, Xiaoyi Jiang b,c,d,n , Nanning Zheng a, Yuehu Liu a, Dachuan Cheng e a

Institute of Artificial Intelligence and Robotics, Xi'an Jiaotong University, China Department of Mathematics and Computer Science, University of Münster, Einsteinstrasse 62, 48149 Münster, Germany c European Institute for Molecular Imaging, University of Münster, Germany d Cluster of Excellence EXC 1003, Cells in Motion, CiM, Münster, Germany e Department of Biomedical Imaging and Radiological Science, China Medical University, Taiwan b

art ic l e i nf o

a b s t r a c t

Article history: Received 19 August 2013 Received in revised form 21 July 2014 Accepted 24 July 2014 Available online 1 August 2014

In this paper we propose the generalized median surface problem and present an exact solution by applying a 3D graph search algorithm. To our knowledge this is a novel class of median computation problem which has not been considered before in the literature. In addition to the theoretic interest in median surface computation we also demonstrate its practical value by means of the task of parameter space exploration without ground truth, which is an effective means of dealing with the difficult parameter selection problem for image segmentation. We present a concrete application for artery boundary detection in sonography (ultrasound imaging). It will be shown that the median computation can not only avoid the parameter training, but also potentially achieve results comparable with trained parameters. The median-based approach can thus be a good alternate, particularly in situations with no ground truth available. & 2014 Elsevier Ltd. All rights reserved.

Keywords: Generalized median Median surface Graph search Parameter handling

1. Introduction Median computation is a useful concept in pattern recognition [1]. Given an object set S in space U, the generalized median is defined by x A U which minimizes the sum of distances to all objects in S. It can be considered as a good representative of the given set. As another motivation, median computation helps us to eliminate (smooth out) erroneous objects by averaging over all objects. From a general point of view, the median concept is motivated by well established results from consensus learning, e.g. multiple classifier combination: by averaging the results of several classifiers a more reliable classification can be achieved [2,3]. The median concept has been concretized to a lot of domains (see Section 2 for a brief overview). In particular, the 2D median contour problem has been investigated [4]. In this work we go a step further from 2D to 3D and consider the related 3D median surface problem. Numerous median computation problems and related algorithms have been developed for specific domains to integrate as much as possible domain-specific knowledge in order to obtain

n Corresponding author at: Department of Mathematics and Computer Science, University of Münster, Einsteinstrasse 62, 48149 Münster, Germany. Tel.: þ 49 251 8333759; fax: þ 49 251 8333755. E-mail address: [email protected] (X. Jiang).

http://dx.doi.org/10.1016/j.patcog.2014.07.019 0031-3203/& 2014 Elsevier Ltd. All rights reserved.

possibly exact solutions in an efficient way. For instance, the generalized median string problem based on the edit distance is shown to be NP-hard, but simplified histogram-based distances reduce the complexity to low-order polynomial time only [5]. For 2D contours dynamic programming can be used to determine the optimal median contour in a time linear to the image size [4]. In this work we will show that for the class of the so-called terrainlike surfaces (to be formally defined in Section 3) and some extensions, a 3D graph search algorithm can be applied to exactly and efficiently solve the median surface problem. The purpose of this paper is twofold. First, we study the median surface problem and solve it by applying a 3D graph search algorithm. To our knowledge this is a class of median computation problem which has not been considered before in the literature. In addition to the theoretic interest in median surface computation in its own right, we also demonstrate its practical value motivated by recent works on exploring the parameter space of segmentation algorithms without ground truth. Similar to multiple classifier systems, ensemble techniques should be developed to achieve the best possible segmentation result on a per-image basis. The remainder of this paper is structured as follows. In the next section we give a brief overview of the generalized median concept. In Section 3 the median surface problem under consideration is defined, which is further motivated by parameter exploration in the context of segmentation problems in Section 4.

Z. Wu et al. / Pattern Recognition 48 (2015) 380–390

The median surface problem will be exactly solved by applying a 3D graph search algorithm (Section 5). Our definition of median surface computation is very general and can cope with positionsensitive weighting of the input surfaces. In Section 6 a scheme is presented for determining the weighting factors. We report experimental results to demonstrate the usability of median surface computation in the context of segmentation parameter exploration in Section 7. Section 8 describes the extension from single surface to multiple surfaces. Additional extensions of our work are discussed in Section 9. Finally, some discussions in Section 10 conclude the paper. This paper is an extended version of the conference paper [6]. The generalized median concept is introduced with more details. We present more algorithmic details of our approach, in particular the weighting scheme in Section 6 and the multiple surface median computation in Section 8. The experimental work in Section 7 has been reorganized and substantially extended. Finally, various extensions are discussed in Section 9.

2. Generalized median concept Assume that we are given a set S of objects in some representation space U and a distance function dðp; qÞ to measure the dissimilarity between any two objects p; q A U. The generalized median p A U of S minimizes the sum of distances to all objects from S, i.e.

381

generalized median of strings [11,23], graphs [14], and clusterings [24], all being of NP-hard class.

3. Median surface problem We start our discussion with terrain-like (height-field) as specified in Definition 1 (defined in the same way as in [25]). Definition 1. A terrain-like surface is a function: f : X  Y-Z with X ¼ f1; 2; …; Mg, Y ¼ f1; 2; …; Ng, and Z ¼ f1; 2; …; Lg. In order to guarantee surface connectivity in 3D, an additional continuity constraint requires jf ðx þ 1; yÞ  f ðx; yÞj r Δx and jf ðx; yþ 1Þ  f ðx; yÞj r Δy for small positive constants Δx and Δy. It corresponds to a discrete version of the Monge patch in differential geometry. For the sake of simplicity we will use the term “surface” only in the following. This class of surfaces is very common in image analysis. In 3D biomedical volume datasets, an important task is to detect such terrain-like surfaces, possibly in an optimal manner. Besides, stacking 2D images along an additional axis, e.g. time, also results in 3D volume datasets. Then, contour detection in an image sequence can be recasted into a 3D surface detection problem. We assume a set S of K input surfaces fS1 ; S2 ; …; SK g and a distance function dðÞ which measures the dissimilarity of two surfaces. Then, the general median surface is defined by K

p ¼ arg min ∑ dðp; qÞ:

S ¼ arg min ∑ dðs; Si Þ

In general p is neither a member of S nor unique. A related concept is the so-called set median, which results from constraining the search to the given set S:

where US represents the space (universe) of all potential solutions, i.e. surfaces within the search volume X  Y  Z. The distance function is defined by

p^ ¼ arg min ∑ dðp; qÞ:

dðs; Si Þ ¼ ∑

pAU qAS

pAS qAS

The set median may serve as an approximative solution for the generalized median or an initialization for iterative generalized ^ qÞ= median computation. This is justified by the fact ∑q A S dðp; ∑q A S dðp; qÞ r 2  2=jSj if the triangle inequality applies to the distance function dðp; qÞ (see [7] for a proof).1 That is, the set median has a consensus error that is at most 2 2=jSj times the consensus error of the generalized median. The median concept has been concretized to many domains including vectors [8], tensors [9], rotations [10], strings [11], graphs [12–14], clusterings [15,16], tree-like shapes [17], time series [18], ranking (permutation) [19], images (for atlas construction) [20], and image segmentations [21,22]. In [4] the 2D median contour problem is investigated. In this work we extend it from 2D to 3D for studying the related median surface problem. We are not aware of the previous work on median computation for the class of surfaces under consideration. Our work thus enriches the generalized median related research. There exist only very few general frameworks for median computation. One such framework described in [14] is based on an embedding into the vector space. The median vector is computed by means of the Weiszfeld algorithm [8] and inversely transformed to the original space. Another general framework [23] computes the weighted mean of pairs of objects in an evolutionary scheme. Both frameworks are approximative only and therefore suitable for those median problems with inherently high computational complexity. Indeed, they have been applied to computing 1

The proof in [7] is given for strings and edit distance. But the proof itself does not use any string-specific property and thus holds for any distance function which satisfies the triangle inequality.

s A US i ¼ 1

M

N

∑ wxy  ρðsðx; yÞ; Si ðx; yÞÞ

ð1Þ

ð2Þ

x¼1y¼1

where ρ is a dissimilarity function for scalar values. Any function suitable for a certain application, e.g. the Minkowski distance Lp, can be used for this purpose. In particular, those from robust statistics [26] may help us to achieve improved performance against outliers in the input surface data. In the simplest case the weight wxy can be set to be constant for all ðx; yÞ positions. However, the larger the wxy is, the more the influence of the particular position and the input surface Si will be imposed on the final median surface. Generally, a positionsensitive weighting gives us more flexibility to incorporate problem-specific knowledge. For our segmentation parameter exploration, we will fully utilize this flexibility (see Section 6). Note that the case of constant weight corresponds to computing a mean surface, i.e. averaging at each point but with guaranteed surface continuity. A comparison of weighted and unweighted median computation will be presented in Section 7.3. The median surface problem defined above can be extended in various ways. We can consider multiple surfaces simultaneously instead of one single surface (see Section 8). In addition an extension to tube-like (cylindrical) surfaces and closed surfaces substantially broadens the applicability of median surface computation in practice (see Section 9).

4. Motivation One motivation of median surface computation is exploring segmentation parameter space without ground truth, i.e., to get the best segmentation result without the help of ground-truth for parameter training or selection.

382

Z. Wu et al. / Pattern Recognition 48 (2015) 380–390

Segmentation algorithms mostly have some parameters and their optimal values are not easy to find. Traditionally, a training image set with (manual) ground truth segmentation is assumed to be available. Then, a subspace of the parameter space is explored to find out the best parameter setting. For each candidate of parameter setting in this subspace, a performance measure is computed in the following way:

 Segment each image in the training set based on this parameter  

setting. Compute a performance measure by comparing the segmentation result and the corresponding ground truth. Compute the average performance measure over all images of the training set.

detection problem, which is solvable by a 3D graph search algorithm in low-order polynomial time. First, we reformulate Eq. (1) as follows: K

S ¼ arg min ∑ dðs; Si Þ s A US i ¼ 1 K

M

N

¼ arg min ∑ ∑ ∑ wxy  ρðsðx; yÞ; Si ðx; yÞÞ s A US i ¼ 1 x ¼ 1 y ¼ 1 M

N

M

N

K

¼ arg min ∑ ∑ wxy  ∑ ρðz ¼ sðx; yÞ; Si ðx; yÞÞ s A US x ¼ 1 y ¼ 1 i¼1 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} cxyz

¼ arg min ∑ ∑ cxyz s A US x ¼ 1 y ¼ 1

|fflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflffl}CðsÞ

The optimal parameter setting is selected as the one with the largest average performance measure. Since fully exploring the subspace can be very costly, several methods have been proposed to reduce the search space such as space subsampling [27], genetic search [28], alternating scheme [29], and orthogonal experiment design [30]. Although this general approach is reasonable and has been successfully practiced in several applications, its critical disadvantage is the need of ground truth segmentation. The manual generation of ground truth is always painful and thus a main barrier of wide use in many situations. Recently, a novel approach is proposed by applying the concept of generalized median for implicitly exploring the parameter space without the need of ground truth segmentation. We assume a reasonable subspace of the parameter space (i.e. a lower and a upper bound for each parameter). This subspace is sampled into a finite number M of parameter settings. Then, the segmentation procedure is run for all the M parameter settings and the generalized median of the M segmentation results is computed. The rationale here is that in line with the ensemble paradigm, the median result tends to be a good one within the explored parameter subspace, as already successfully demonstrated for 2D contour detection [4] and region segmentation [21,22]. Segmentation of surfaces is one of the most important problems in (biomedical) image analysis. But in most cases it involves parameter selection, which is however typically not easy to handle. Thus, median surface computation can directly help us to alleviate the parameter problem in 3D surface segmentation as well. As another situation of potential application we consider the segmentation of 2D images along the time axis. Many algorithms from the literature, e.g. [31], perform the segmentation independently on all images and thus cannot guarantee a continuous segmentation over time, which is naturally desired when working with sequences of images. If the parameter space exploration technique described above is applied to the 3D volumes formed by stacking all frame-wise segmentations along the time axis, we obtain a continuous temporal segmentation without any extra effort as a nice spin-off of handling the parameter problem. It is important to emphasize that we do not intend to solve the segmentation problem, but instead study the parameter exploration problem for image segmentation based on some baseline segmentation algorithm. In this sense our approach is on a metalevel above the segmentation level. The median surface computation is not only an interesting topic in its own right but also of substantial practical value. This motivates us to find an efficient way for exact median surface computation.

5. Exact computation by 3D graph search In this section we show that the median surface problem defined in Eq. (1) can be transformed into an optimal 3D surface

¼ arg min CðsÞ

ð3Þ

s A US

A candidate solution is a surface s A U S characterized by the z-value sðx; yÞ for each position (x,y) on the grid X  Y. Each point ðx; y; zÞ in the volume X  Y  Z is assigned a cost cxyz, which is determined by its deviations (in z-direction) from the K input surfaces Si ðx; yÞ and the position-specific weight wxy. Then, the goodness of a candidate solution surface s is evaluated by the measure C(s), i.e. summing up the costs of all positions. Therefore, the median surface is simply the optimal surface with minimal cost from the solution space US (consisting of all terrain-like surfaces within the volume X  Y  Z). This discussion immediately leads to the following new optimization problem. We first compute a cost cxyz for each point ðx; y; zÞ in the volume X  Y  Z. Then, the median surface is determined by finding the terrain-like surface within the volume with the minimal sum of costs. It is important to mention that this optimization problem cannot be solved by computing the optimal z-value for each of the M  N positions (x,y) independently (i.e. by enumerating all z-values from Z and minimizing cxyz). Doing it in this simple way, we will encounter the trouble of generating a discontinuous resultant surface in general. Only for some elementary cases (e.g. constant weight wxy and ρ ¼ L22 ) the straightforward position-wise optimization will deliver an optimal continuous resultant surface. But in general, a global optimization approach is needed. For the special case N ¼1 (i.e. the y-axis vanishes), the 3D optimal surface segmentation becomes a 2D optimal contour detection problem. This reduced problem is much easier to solve, e.g. by a highly efficient dynamic programming algorithm [4]. Unfortunately, there is no direct way of extending the dynamic programming solution to the general 3D case (N 4 1). Note that Eq. (3) can be viewed as a first order MRF minimum energy problem, i.e., EðXÞ ¼ ∑ Gðv; X v Þ þ ∑ Pðu; v; X u ; X v Þ vAV

ðu;vÞ A E

but without the prior (the second term) where X is a labeling function. Obtaining the exact solution of the first order MRF is a classical problem in computer vision [32]. The work [33] solved the binary variable first order MRF with specified priors. The work [34] discussed using the first order MRF to solve the image segmentation problems and gave the exact solution when Gðv; X v Þ is convex and Pðu; v; X u ; X v Þ is a linear function. In [35,36] a general first order MRF energy function with convex priors and linear order labels was solved with exact solution. In [37] an optimal net surface problem is proposed, which can be viewed as a concrete case of [35,36], because the base graph in [37] corresponds to the original graph in [35,36]. Later, in [25], the problem of extracting the optimal surface from a 3D volume is proposed, which is a special case of [37]. Fortunately, our problem can be viewed as the optimal surface search as formulated in [25] and we

Z. Wu et al. / Pattern Recognition 48 (2015) 380–390

thus use the algorithm in [25] to solve our optimal surface detection problem in Eq. (3). In the following we give a brief presentation of the most important steps of this algorithm and the readers are referred to [25] for further details. A node-weighted directed graph G ¼ ðV ; E; WÞ is constructed as follows. For each point ðx; y; zÞ in the volume X  Y  Z a corresponding node Vðx; y; zÞ is defined in G, whose weight Wðx; y; zÞ is assigned according to ( cxyz  cxy;z  1 ; z 4 1 Wðx; y; zÞ ¼ ð4Þ z¼1 cxyz ; where cxyz is the cost defined in Eq. (3). Two types of edges are defined for G: E ¼ Ea [ Er . The set Ea of intraposition edges models the connections within the same position (x,y). Each node Vðx; y; zÞ (z 4 1) has a direct edge to the node V ðx; y; z 1Þ below it, i.e., Ea ¼ f o Vðx; y; zÞ; V ðx; y; z 1Þ 4 j z 4 1g The set Er of interposition edges models the connections of adjacent positions. The node V ðx; y; zÞ is connected to the nodes Vðx þ 1; y; maxð0; z  Δx ÞÞ at adjacent position ðx þ1; yÞ, i.e. Δx below z; z Z Δx . Similar edges are drawn to the other three adjacent positions ðx  1; yÞ; ðx; y þ 1Þ, and ðx; y 1Þ. Overall, the interposition edges are defined by 8 f o Vðx; y; zÞ; > > > > < f o Vðx; y; zÞ; r E ¼ f o Vðx; y; zÞ; > > > > : f o Vðx; y; zÞ;

Vðx þ 1; y; maxð0; z  Δx ÞÞ 4

j x A f1; …; M  1g; z A Zg [

Vðx  1; y; maxð0; z  Δx ÞÞ 4 Vðx; y þ1; maxð0; z  Δy ÞÞ 4 Vðx; y 1; maxð0; z  Δy ÞÞ 4

j x A f2; …; Mg; z A Zg [ j y A f1; …; N 1g; z A Zg [ j y A f2; …; Ng; z A Zg

Given the constructed directed graph G, a closed set C is a subset of nodes such that all successors of any nodes in C are also contained in C. The cost of a closed set is the total cost of all its nodes. For any feasible surface N , the subset of nodes on or below N in G, namely C ¼ fV ðx; y; zÞ j z r N ðx; yÞg, forms a closed set in G. Due to the node cost defined in (4), the costs of N and C are clearly equal. In [25] it is shown that the original optimal surface detection problem is equivalent to finding a minimum nonempty closed set in G. This problem is well studied in graph theory and can be solved by computing a minimum s  t cut in a related graph Gst (see [25] for the details of constructing Gst from G). In our implementation we used the Boykov–Kolmogorov algorithm [32] to compute the minimum s  t cut. For a graph with n nodes and m edges, the theoretical worst-case time complexity for this algorithm is Oðn2 mcÞ, where c is the cost of the minimum cut. In summary our approach to median surface computation consists of the following steps:

 Given a set S of K surfaces fS1 ; S2 ; …; SK g, we construct a 3D search volume X  Y  Z. Each point ðx; y; zÞ in this volume is



383

assigned a cost cxyz, which is determined by its deviations (in z-direction) from the K input surfaces Si ðx; yÞ and the positionspecific weight wxy. Then, the goodness of a candidate solution surface s is measured by C(s) in Eq. (3), i.e. summing up the costs of all (x,y) positions. Thus, the median surface is simply the optimal surface with minimal cost in the search volume. Apply the 3D graph search algorithm in [25] to find the optimal surface with minimal cost in the search volume.

The first step transforms the median surface computation problem into one of the optimal surface detections, which is then solved by the 3D graph search algorithm in the second step. Note that our approach can guarantee the true median surface under the hard continuity constraint. The reason is that the median surface should also satisfy the continuity constraint and the search graph construction thus uses the continuity constraint. It is worth mentioning that the number of input surfaces K does not pose any problem (extra complexity) to our approach at all because the input surfaces are only used to calculate the voxels' cost in the search volume (the search volume itself is fixed). Afterwards, the complexity of median surface computation does not depend on this number anymore.

6. Weighting scheme The median surface definition depends on the weight wxy. Generally, a position-sensitive weighting gives us more flexibility to incorporate problem-specific knowledge. In this section we present a position-wise weighting scheme along this line. For any position (x,y) in the cost volume X  Y  Z, each zi ðx; yÞ A Si can be regarded as a sample of a random variable Fxy, which is assumed to form a Gaussian distribution. As an illustration, Fig. 1 shows the distribution of the z-values from all input surfaces for a particular position (x,y) and the estimated Gaussian density function. For the zi ðx; yÞ A Si , the corresponding wxy is then chosen to be the density at z ¼ sðx; yÞ. This weighting scheme has a simple interpretation. If some z-value appears frequently in the input surfaces, it is likely that the median surface should have a similar z-value at this position. The density-based weighting reflects this expectation and reinforces achieving the expected result. The importance of position-sensitive weighting is further illustrated in Fig. 2. The distributions for two positions ðx1 ; y1 Þ and ðx2 ; y2 Þ are considered. The green circle and the yellow triangle symbolize two different input surfaces. At position ðx1 ; y1 Þ the green surface should receive a high weight since the corresponding z-value plays a strong majority role. In contrast the yellow surface should be under-weighted. The situation at position ðx2 ; y2 Þ is exactly opposite. This illustrative example clearly demonstrates that a constant weight is generally not a good choice. In our experimental work we compare the constant weighting and the weighting scheme presented above, and the results reported in Section 7 confirm the superiority of position-sensitive weighting.

7. Experimental results In this section we demonstrate a practical use of median surface computation by applying the algorithm described above to segmentation parameter exploration in sonography. 7.1. Sonographic data and experimental settings Fig. 1. Weighting scheme: the distribution of z-value from all input surfaces for a particular position (x,y) and its estimated Gaussian density function.

The task considered here is to alleviate the painful parameter selection problem in the extraction of artery boundaries from

384

Z. Wu et al. / Pattern Recognition 48 (2015) 380–390

Fig. 2. Illustration of importance of position-sensitive weighting. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

sonographic videos. An artery has a near wall and a far wall,2 as illustrated in Fig. 3(b). Along with the time axis the 2D images can be regarded as a 3D volume, see Fig. 3(a). Four sonographic videos from patients were used in our experiments. For these videos, a ground truth of the arterial walls (golden standard) was labeled manually. In total we thus have 8 test cases (4 videos, near and far wall). The dynamic programming approach from [38] was applied to detect the two contours. This algorithm has two parameters. For both parameters a reasonable interval ½0:1; 1:9 was identified and sampled into 10 values (0.1, 0.3, …, 1.7, 1.9). In total we thus considered 100 different parameter settings. For each video, these parameter settings were used to generate 100 near wall surfaces and 100 far wall surfaces. Then, a median surface was computed from each of the two surface ensembles. The diversity of quality of the 100 segmented surfaces in each test case is documented by the statistics in Table 1. The quality is measured by comparing each surface with the corresponding GT based on L1. The 100 input surfaces thus lead to 100 performance measures (average deviation per (x,y) position). Each table item presents their value interval, mean, and standard deviation. For all 8 test cases the corresponding ensemble shows some diversity, indicating the varying suitability of the distinct parameter settings. It is this diversity that promises potential by means of fusion (median computation). For comparison purpose the best-performing one among the 100 parameter settings was determined by using all test videos and comparing with the ground truth. This is the ultimate optimal segmentation (subject to the explored parameter space) on a perimage basis, which may be approached, but cannot be achieved by segmentation algorithms in general. In all our tests the comparison between two surfaces, e.g. a segmented surface and a ground truth, was done by computing the average L1 deviation in z per (x,y) position. 7.2. Comparison with the ground truth We compared our median result with the ground truth. The average L1 deviation in z is shown in Table 2. Contrasting with the ensemble statistics in Table 1, the ensemble solution by means of median computation consistently achieves a quality level between the minimum and the mean, in many cases actually near the minimum. In addition, the results from the best-performing parameter setting (BP) were also compared with GT. Four examples of comparison from two different frames are given in Fig. 4 for 2 More precisely, these are the intimas of the near and the far wall, see Section 8.2.

illustration purpose. While in (b) and (d) the median and BP results are absolutely comparable, (a) and (c) illustrate superior median computation compared to BP. As can be seen in Table 2, our median segmentation results reach the quality of bestperforming parameter setting very well. Using our median surface algorithm thus can not only avoid the parameter training (which is only possible with existing ground truth), but also potentially achieve segmentation results similar to the best parameters. This fact is clearly due to the ensemble nature of the median surface computation. Due to the small number of videos we did not perform cross validation. Instead, the best-performing parameter setting was determined by using all available videos. This results in performance measures which usually cannot be fully reached when the optimal parameters are learned using training data and applied to unknown test data. Therefore, the comparison in Table 2 is based on an overoptimistic estimation of training-based performance, making the performance of median computation even more favorable. 7.3. Comparison with unweighted median computation In Section 6 we have argued the importance of a positionsensitive weighting and presented a weighting scheme. A test series was performed to experimentally confirm the superior performance of weighting. Again, the median surfaces are compared with GT based on L1. The average difference reported in Table 3 shows an almost consistent improvement of the weighting scheme over the unweighted median computation. The weighting scheme presented in Section 6 is still rather simple. Note that the performance evaluation in this work mainly follows the goal of demonstrating the usability of the median computation and with regard to this test series the potential of weighting. From an application point of view, it certainly remains to improve the weighting scheme and to maximize the benefit of position-sensitive weighting. But this is beyond the scope of our current work. 7.4. Robustness study We manually changed the population of input surfaces to simulate various test situations. Three test scenarios were generated by selecting 60 out of the 100 input surfaces:

 50 best and 10 worst input surfaces in terms of GT comparison.  40 best and 20 worst input surfaces.  30 best and 30 worst input surfaces. In particular, the last test scenario has a high percentage (50%) of parameter settings, which are relatively non-optimal in terms of

Z. Wu et al. / Pattern Recognition 48 (2015) 380–390

385

7.5. Computational time The experiments were performed on an Intel Xeon x5560 2.8 GHz with 16G memory. The program was written in C (Boykov–Kolmogorov algorithm) and Matlab (other parts). The computational time is listed in Table 5. Also shown is the resolution of the subpart of a video, in which the near or the far wall is searched for. The computational time increases with the volume size (image size and number of images). Generally, the graph search algorithm is quite fast. 7.6. Discussions Fig. 3. Sonographic video. (a) Along with the time axis the 2D images build a 3D volume. (b) The near (blue) and the far (red) wall of the artery in a single image. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

Table 1 Statistics of the 100 segmented surfaces: [minimum, maximum], mean, standard deviation (unit: pixels). Video

#Images

Near wall

Far wall

1 2 3 4

86 86 111 73

[0.46, 1.00], 0.78, 0.17 [0.41, 1.35], 0.82, 0.33 [0.89, 2.82], 2.00, 0.73 [0.45, 1.76], 1.08, 0.56

[0.51, 0.71], 0.62, 0.05 [0.47, 0.75], 0.61, 0.08 [0.63, 0.96], 0.73, 0.09 [0.45, 1.45], 0.75, 0.22

Table 2 Comparison with the ground truth (unit: pixels). Video

Comparison type

Near wall

Far wall

1

Median vs. GT BP vs. GT

0.63 0.61

0.57 0.55

2

Median vs. GT BP vs. GT

0.48 0.46

0.50 0.49

3

Median vs. GT BP vs. GT

0.98 0.96

0.67 0.66

4

Median vs. GT BP vs. GT

0.48 0.46

0.54 0.51

the studied parameter interval. The performance is given in Table 4. In the first test scenario similar behavior as in Table 2 can be observed. The median computation basically achieves the same performance as the best-performing parameter setting. Not surprisingly, with the increasing percentage of non-optimal parameter settings, the performance of median computation slightly drops. But even in the last test scenario the performance drop is not dramatic, indicating the potential of robustness by using the median surface approach in dealing with the segmentation parameter problem. The robustness is mainly due to the use of L1 norm. In robust statistics [26] it is known that L1 is robust and can tolerate up to 50% outliers. Any other robust dissimilarity function ρ with a possibly high breakdown point will make the median surface computation robust as well. Since our approach is an ensemble one, the ensemble quality certainly influences the final result. It is thus important to select a reasonable parameter subspace for sampling. In practice, this is mostly not difficult to achieve. Together with a robust dissimilarity function it can be expected that the “bad” samples will be well tolerated by the median surface computation and an overall good quality of median surface can be reached.

Our approach to deal with the parameter problem is in line with the ensemble paradigm which is a powerful boosting technique in machine learning and pattern recognition [2,3]: by averaging the results of several classifiers a more reliable classification can be achieved. In our case an ensemble of K surfaces is averaged to achieve an “optimal” segmentation. This principle has been successfully validated for 2D contour detection [4] and image segmentation [21,22] before. In the current work this ensemble approach is demonstrated for parameter space exploration in the context of 3D surface segmentation. In Section 4 we have discussed several parameter exploration methods [27–30]. However, their common critical disadvantage is the need of ground truth segmentation. The manual generation of ground truth is always painful and thus a main barrier of wide use in many situations. Our intention in the application part of this work is to overcome the difficulty with parameters, to some extent at least, without assuming the ground truth. Due to the very different nature of our unsupervised ensemble approach and other supervised methods we do not perform direct comparison with them. However, we have shown that our median segmentation results reach the quality of best-performing parameter setting very well. Since even ground truth based parameter learning cannot achieve the performance of best-performing parameter setting, this result is an indirect proof of the power of our approach compared to the existing supervised methods. In our tests the involved dynamic programming approach from [38] has two parameters only. In the case of more parameters the number of sampled parameter settings certainly increases. One can either fix some non-critical parameters and sample a small number of remaining critical parameters, or sample the (noncritical) parameters relatively sparsely. Another solution is to compute the K input surfaces in a parallel manner. It is important to emphasize that this is a general issue related to many, if not all, parameter learning methods from the literature. The application of our median surface computation intends to alleviate the parameter problem. It is therefore important to make sure that our approach does not introduce new parameter(s). This is the case since all steps described are parameter-free. For instance, the weighting scheme described in Section 6 is based on function fitting. Note that the continuity values Δx and Δy in Definition 1 are not parameters introduced by our median computation algorithm. Indeed, any baseline surface detection algorithm, which delivers the input for our median surface computation, needs these values as well. They mean the degree of surface continuity a user wishes to achieve and thus have to be specified by the user, most likely in an application-dependent manner. 8. Multiple surfaces median computation Due to the imperfections of medical imaging techniques, insufficient or even deficient image-derived information may be

386

Z. Wu et al. / Pattern Recognition 48 (2015) 380–390

Fig. 4. Examples of comparison of GT, BP, and median.

Table 3 Comparison weighted vs. unweighted (unit: pixels). Video

Variant

Near wall

Far wall

1

Weighted Unweighted

0.63 0.67

0.57 0.57

2

Weighted Unweighted

0.48 0.54

0.50 0.53

3

Weighted Unweighted

0.98 1.08

0.67 0.68

4

Weighted Unweighted

0.48 0.52

0.54 0.53

Definition 2. Two coupled terrain-like surfaces contain two terrain-like single surfaces: f i : X  Y-Z with X ¼ f1; 2; …; Mg, Y ¼ f1; 2; …; Ng, and Z ¼ f1; 2; …; Lg, i¼1,2. The pre-defined coupling constraints δl and δu specify, respectively, the minimum and maximum distances between the two non-intersecting surfaces, i.e. δl r jf 1 ðx; yÞ  f 2 ðx; yÞj rδu . This definition can be easily extended to the case of K ( 4 2) surfaces. Denoting the K surfaces by f 1 ; f 2 ; …; f K , we need to define for each pair of neighboring surfaces fi and fj the related constraints δijl and δiju. The 3D graph search algorithm [25] has an extension for K coupled surfaces. In this case the individual graphs Gi ¼ ðV i ; Ei ; W i Þ; i ¼ 1; …; K, are combined to build an overall search graph: G ¼ ðV 1 [ ⋯ [ V K ; E1 [ ⋯ [ EK [ Es ; W 1 [ ⋯ [ W K Þ

available. This insufficiency can be alleviated by using clues from other mutually related boundaries or surfaces. Indeed, codetection of multiple coupled surfaces frequently yields superior performance compared to the common single-surface detection approaches [31,25,38–40]. Thus, it is essential to extend the single terrain-like surface median problem to the case of multiple coupled surfaces. 8.1. Algorithmic extension We start with the definition of multiple surfaces (defined in the same way as in [25]).

The extra intersurface edges in Es model the relations between surfaces. As done in [25], we specify the construction of Es for K ¼2. The ideas can be easily generalized to handling more than two surfaces. Supposing f1 is above f2, then 8 > < f o V 1 ðx; y; zÞ; V 2 ðx; y; z  δu Þ 4 jz Zδu g [ s E ¼ f o V 2 ðx; y; zÞ; V 1 ðx; y; z  δl Þ 4 jz o L  δl g [ > : f o V ð0; 0; δ Þ; V ð0; 0; 0Þ 4 g 1

l

2

Note that the formulation above is given for non-intersecting surfaces. The needed extension for intersecting surfaces can be found in [25], also the details of computing the K optimal surfaces.

Z. Wu et al. / Pattern Recognition 48 (2015) 380–390

387

The K-surface median surface computation, however, may suffer from high demand of computation. Compared to the original single-surface algorithm, the worst-case time complexity of the

Boykov–Kolmogorov algorithm (Oðn2 mcÞ) increases at least by a factor of K3 (because of K  n nodes and K  m plus extra intersurface edges in the enlarged graph).

Table 4 Performance of three populations of input surfaces (unit: pixels).

8.2. Experimental results

Comparison type

Near wall

Far wall

Median vs. GT Bp vs. GT

0.57 0.61

0.55 0.55

40 best þ 20 worst

Median vs. GT BP vs. GT

0.56 0.61

0.55 0.55

30 best þ 30 worst

Median vs. GT BP vs. GT

0.68 0.61

0.58 0.55

Median vs. GT BP vs. GT

0.44 0.46

0.49 0.49

40 best þ 20 worst

Median vs. GT BP vs. GT

0.44 0.46

0.49 0.49

30 best þ 30 worst

Median vs. GT BP vs. GT

0.51 0.46

0.50 0.49

Video 1 50 best þ 10 worst

Video 2 50 best þ 10 worst

The experimental work for the single-surface case reported in Section 7 is restricted to the intima of near and far walls only. In fact the another so-called adventitia is also of interest, see Fig. 5 for an illustration. Using the same sonographic data and experimental settings as described in Section 7.1, we have performed experiments for multiple surface median computation (K ¼2). For the near wall, a pair of median intima and median adventitia are computed simultaneously and the same is also done for the far wall. Similar to Table 2, Table 6 lists a comparison with the GT and the best-performing parameter setting. Again, our median segmentation results reach the quality of best-performing parameter setting very well. The corresponding computational time is reported in Table 7. As expected, it is much higher than that for the single surface case.

9. Extensions

Median vs. GT BP vs. GT

0.97 0.96

0.66 0.66

40 best þ 20 worst

Median vs. GT BP vs. GT

0.96 0.96

0.66 0.66

The 3D graph search algorithm in [25] designed for terrain-like surface detection has several variants of extensions which further broaden the usability of median surface computation. The extension to the case of multiple surfaces in already discussed in Section 8. In this section we briefly discuss three additional extensions and another issue related to computational cost.

30 best þ 30 worst

Median vs. GT BP vs. GT

1.10 0.96

0.67 0.66

9.1. Tube-like surfaces

Median vs. GT BP vs. GT

0.47 0.46

0.51 0.51

40 best þ 20 worst

Median vs. GT BP vs. GT

0.47 0.46

0.51 0.51

30 best þ 30 worst

Median vs. GT BP vs. GT

0.48 0.46

0.53 0.51

Video 3 50 best þ 10 worst

Video 4 50 best þ 10 worst

Table 5 Computational time (seconds). Video

1 2 3 4

#images

86 86 111 73

Near wall

Sometimes, the desired surface is required to be wraparound along some direction (say, the x- or y-axis). For instance, this is a common situation in medical image analysis to detect tube-like (cylindrical) surfaces. In this case a tubular surface is first unfolded into a terrain-like surface using cylindrical coordinate transform [41]. Then, the boundary of the tubular surface in the original image corresponds to a terrain-like surface to be detected in the unfolded image and the 3D graph search algorithm [25] can thus be applied. Moreover, the first and last rows along the unfolding plane should satisfy the continuity constraints as well, which is guaranteed by a modified graph construction. Median surface computation is directly applicable in this case.

Far wall

Resolution

Time

Resolution

time

45  121 47  163 72  188 40  126

3.2 6.3 21.8 2.4

49  121 46  163 56  188 45  126

4.0 6.8 19.1 7.5

9.2. Closed surfaces Besides the tube-like surfaces, closed surfaces are also frequently encountered in practical situations, but cannot be directly transformed to a general volume like the tube-like surfaces. Fortunately, as we mentioned in Section 5, the optimization

Fig. 5. Illustration of intima and adventitia for near wall and far wall in artery ultrasound image (compare Fig. 3).

388

Z. Wu et al. / Pattern Recognition 48 (2015) 380–390

Table 6 Comparison with the ground truth (unit: pixels) in multiple surface setting. Video

Comparison type

Near wall

Far wall

Intima

Adventitia

Intima

Adventitia

1

Median vs. GT BP vs. GT

0.63 0.61

0.51 0.48

0.57 0.55

0.47 0.47

2

Median vs. GT BP vs. GT

0.48 0.46

0.80 0.80

0.50 0.49

0.81 0.80

3

Median vs. GT BP vs. GT

0.98 0.96

1.43 1.46

0.67 0.66

0.76 0.75

4

Median vs. GT BP vs. GT

0.48 0.46

0.74 0.69

0.54 0.51

0.53 0.50

Table 7 Computational time (seconds) in multiple surface setting. video

1 2 3 4

#images

86 86 111 73

Near wall

Far wall

Resolution

Time

Resolution

Time

45  121 47  163 72  188 40  126

93.44 172.79 405.76 74.36

49  121 46  163 56  188 45  126

96.62 164.89 379.08 71.59

problem can be regarded as a general MRF problem and both [36] and [37] provide ways for an exact computation. Some other closed surface works [42,43] are based on dedicated base graph generation and the algorithm in [37]. Therefore, the median surface computation is also applicable for closed surfaces.

9.3. Smoothness constraint Smoothness constraint is doubtless desirable for surface extraction problems. Although it is quite natural for the variational methods such as active contours and level set approaches, imposing smoothness constraint to graph search based algorithms is not trivial, especially for problems of boundary or surface extraction. In [44,45], the graph search algorithm from [37] is augmented by smoothness constraints by means of arc penalty. For our median surface computation these extended versions can be used if smoothness should be integrated into the detection.

9.4. Computational cost As discussed before, the K-surface median surface computation may suffer from high demand of computation. Even for the singlesurface scenario the computation demand may significantly increase when dealing with videos. A large number T of frames will cause both space and time problems. For typical image sizes in medical image analysis, for instance, a large T may result in memory overflow in practice. Even if sufficient memory can be allocated, the computation time can become too long to be useful for real applications. Both situations (K 4 1 and large T) call for a practical solution. To reduce the complexity one can partition a video into blocks and process each block by the algorithm from [25]. A separate handling of the blocks, however, cannot guarantee the consistency between two adjacent blocks. Therefore, suitable methods are needed to take care of smooth transition from one block to another. We are currently investigating various block-wise approximate variants.

10. Discussions and conclusion In this paper we have formulated the generalized median surface problem and presented its exact solution by means of a 3D graph search algorithm. The median concept has been concretized for many domains. This work can be understood as an extension of the simpler case of 2D median contour computation by means of dynamic programming [4] to 3D median surface. Since dynamic programming cannot be directly extended to 3D in an efficient manner, we resort to the 3D graph search algorithm [25]. In addition to terrain-like surfaces including single and multiple coupled surfaces, the median computation can also be applied to deal with tubular surfaces and closed surfaces (of both types). This work is motivated by the task of parameter space exploration without ground truth, which is an effective means of dealing with the difficult parameter problem and has been successfully applied to domains like 2D contour detection [4] and region segmentation [21,22]. Our median surface computation algorithm thus provides a useful tool for parameter exploration in 3D surface segmentation or 2D contour segmentation in a temporal context. A concrete application has been demonstrated on artery boundary detection in sonography, which confirmed the findings from the previous studies. That is, the median computation can not only avoid the parameter training, but also potentially achieve results comparable with trained parameters. Parameter training is only possible with the existing ground truth, which is not always available. The median-based approach can thus be a good alternate in case of no ground truth. One question that may be of interest to the readers is why we do not directly apply the graph search to solve the 3D surface segmentation problem. Our motivation in this work is not the segmentation itself, but exploring segmentation parameter space without ground truth. In this sense our approach is on a metalevel above the segmentation level and can be based on any baseline segmentation algorithm. Even if it is not our intention to directly apply the graph search algorithm for 3D segmentation, we also tested this option by applying the same preprocessing as in the baseline segmentation method used in this work [38]. For video 1, see Table 2, the performance measure (average L1 deviation in z) is 0.57 (near wall) and 0.53 (far wall), which is comparable with our ensemble solution (0.63 and 0.57, respectively). In contrast, video 3 leads to 1.74 (near all) and 1.44 (far wall), which is substantially worse than our ensemble solution (0.98 and 0.67, respectively). Overall, the achieved segmentation performance is not always satisfactory and can neither beat the baseline segmentation algorithm nor our ensemble approach. The reason for this observation may be the strong basic assumption of graph search that the desired surface must be the minimum-cost surface. Generally, this assumption poses substantial challenges to feature computation in the preprocessing step. In practice, high-level properties are thus typically adopted as constraints to relax this assumption which may not be easy to realize in the framework of graph search. This motivates further work on improved preprocessing that may enhance the performance of the graph search algorithm for 3D segmentation, but this is a research topic different from the main motivation of the current work, i.e. exploring segmentation parameter space without ground truth. Additional future work will be to study alternative weighting schemes, for example, weighting schemes that adopt more highlevel and prior properties such as curvature. Furthermore, we will extend our ensemble approach to more applications as motivated in Section 4. As an example, the median result could also serve for selecting the best preprocessing and feature from various potential options.

Z. Wu et al. / Pattern Recognition 48 (2015) 380–390

Conflict of interest None declared. Acknowledgments This work is supported by the State Key Program of National Nature Science Foundation of China, Grant no. 61231018. Xiaoyi Jiang was supported by the Deutsche Forschungsgemeinschaft (DFG): SFB 656 MoBil (project B3) and EXC 1003 Cells in Motion – Cluster of Excellence. References [1] X. Jiang, H. Bunke, Learning by generalized median concept, in: P.Wang (Ed.), Pattern Recognition and Machine Vision, River Publishers, Aalborg, 2010, pp. 1–16. [2] L. Rokach, Ensemble-based classifiers, Artif. Intell. Rev. 33 (1–2) (2010) 1–39. [3] Z.-H. Zhou, Ensemble Methods: Foundations and Algorithms, CRC Press, Boca Raton, London, New York, 2012. [4] P. Wattuya, X. Jiang, A class of generalized median contour problem with exact solution, in: Proceedings of the Joint IAPR International Workshop on SSPR and SPR, Lecture Notes in Computer Science, vol. 4109, Springer, Berlin, Heidelberg, 2006, pp. 109–117. [5] C. Solnon, J.-M. Jolion, Generalized vs set median strings for histogram-based distances: algorithms and classification results in the image domain, in: Proceedings of GbR, 2007, pp. 404–414. [6] Z. Wu, X. Jiang, N. Zheng, Y. Liu, D.-C. Cheng, Exact computation of median surfaces using optimal 3d graph search, in: Proceedings of 9th IAPR-TC-15 International Workshop on Graph-Based Representations in Pattern Recognition, Lecture Notes in Computer Science, vol. 7877, Springer, Berlin, Heidelberg, 2013, pp. 234–243. [7] D. Gusfield, Algorithms on Strings, Trees, and Sequences: Computer Science and Computational Biology, Cambridge University Press, Cambridge, UK, 1997. [8] E. Weiszfeld, F. Plastria, On the point for which the sum of the distances to ngiven points is minimum, Ann. Oper. Res. 167 (2009) 7–41. [9] M. Welk, J. Weickert, F. Becker, C. Schnörr, C. Feddern, B. Burgeth, Median and related local filters for tensor-valued images, Signal Process. 87 (2) (2007) 291–308. [10] R.I. Hartley, J. Trumpf, Y. Dai, H. Li, Rotation averaging, Int. J. Comput. Vis. 103 (3) (2013) 267–305. [11] X. Jiang, J. Wentker, M. Ferrer, Generalized median string computation by means of string embedding in vector spaces, Pattern Recognit. Lett. 33 (7) (2012) 842–852. [12] X. Jiang, A. Münger, H. Bunke, On median graphs: properties, algorithms, and applications, IEEE Trans. Pattern Anal. Mach. Intell. 23 (10) (2001) 1144–1151. [13] L. Mukherjee, V. Singh, J. Peng, J. Xu, M.J. Zeitz, R. Berezney, Generalized median graphs and applications, J. Comb. Optim. 17 (1) (2009) 21–44. [14] M. Ferrer, D. Karatzas, E. Valveny, I. Bardají, H. Bunke, A generic framework for median graph computation based on a recursive embedding approach, Comput. Vis. Image Underst. 115 (7) (2011) 919–928. [15] S. Vega-Pons, J. Ruiz-Shulcloper, A survey of clustering ensemble algorithms, Int. J. Pattern Recognit. Artif. Intell. 25 (3) (2011) 337–372. [16] H. Alizadeh, B. Minaei-Bidgoli, H. Parvin, Optimizing fuzzy cluster ensemble in string representation, Int. J. Pattern Recognit. Artif. Intell. 27 (2) (2013) 1350005. [17] A. Feragen, S. Hauberg, M. Nielsen, F. Lauze, Means in spaces of tree-like shapes, in: Proceedings of ICCV, 2011, pp. 736–746. [18] F. Petitjean, A. Ketterlin, P. Gançarski, A global averaging method for dynamic time warping, with applications to clustering, Pattern Recognit. 44 (3) (2011) 678–693. [19] S.C. Boulakia, A. Denise, S. Hamel, Using medians to generate consensus rankings for biological data, in: Proceedings of 23rd International Conference on Scientific and Statistical Database Management, 2011, pp. 73–90. [20] Y. Xie, J. Ho, B.C. Vemuri, Image atlas construction via intrinsic averaging on the manifold of images, in: Proceedings of CVPR, 2010, pp. 2933–2939.

389

[21] L. Franek, D.D. Abdala, S. Vega-Pons, X. Jiang, Image segmentation fusion using general ensemble clustering methods, in: Proceedings of Asian Conference on Computer Vision, vol. 4, 2010, pp. 373–384. [22] L. Franek, X. Jiang, P. Wattuya, Local instability problem of image segmentation algorithms: systematic study and an ensemble-based solution, Int. J. Pattern Recognit. Artif. Intell. 26 (5) (2012) 1265004. [23] L. Franek, X. Jiang, Evolutionary weighted mean based framework for generalized median computation with application to strings, in: Proceedings of Joint IAPR International Workshop on SSPR and SPR, Lecture Notes in Computer Science, vol. 7626, Springer, Berlin, Heidelberg, 2012, pp. 70–78. [24] L. Franek, X. Jiang, Ensemble clustering by means of clustering embedding in vector spaces, Pattern Recognit. 47 (2) (2014) 833–842. [25] K. Li, X. Wu, D. Chen, M. Sonka, Optimal surface segmentation in volumetric images – a graph-theoretic approach, IEEE Trans. Pattern Anal. Mach. Intell. 28 (1) (2006) 119–134. [26] C. Stewart, Robust parameter estimation in computer vision, SIAM Rev. 41 (3) (1999) 513–537. [27] J. Min, M.W. Powell, K.W. Bowyer, Automated performance evaluation of range image segmentation algorithms, IEEE Trans. Syst. Man Cybern. Part B 34 (1) (2004) 263–271. [28] G. Pignalberi, R. Cucchiara, L. Cinque, S. Levialdi, Tuning range image segmentation by genetic algorithm, EURASIP J. Adv. Signal Process. 2003 (8) (2003) 780–790. [29] L. Franek, X. Jiang, Alternating scheme for supervised parameter learning with application to image segmentation, in: Proceedings of CAIP, 2011, pp. 118–125. [30] L. Franek, X. Jiang, Orthogonal design of experiments for parameter learning in image segmentation, Signal Process. 93 (6) (2013) 1694–1704. [31] D. Cheng, X. Jiang, Detections of arterial wall in sonographic artery images using dual dynamic programming, IEEE Trans. Inf. Technol. Biomed. 12 (6) (2008) 792–799. [32] Y. Boykov, V. Kolmogorov, An experimental comparison of min-cut/max-flow algorithms for energy minimization in vision, IEEE Trans. Pattern Anal. Mach. Intell. 26 (9) (2004) 1124–1137. [33] D. Greig, B. Porteous, A. Seheult, Exact maximum a posteriori estimation for binary images, J. R. Stat. Soc. Ser. B: Methodol. (1989) 271–279. [34] D. Hochbaum, An efficient algorithm for image segmentation, Markov random fields and related problems, J. ACM 48 (4) (2001) 686–701. [35] H. Ishikawa, D. Geiger, Segmentation by grouping junctions, in: Proceedings of CVPR, 1998, pp. 125–131. [36] H. Ishikawa, Exact optimization for Markov random fields with convex priors, IEEE Trans. Pattern Anal. Mach. Intell. 25 (10) (2003) 1333–1336. [37] X. Wu, D.Z. Chen, Optimal net surface problems with applications, in: Automata, Languages and Programming, Springer, Berlin, Heidelberg, 2002, pp. 1029–1042. [38] D.-C. Cheng, A. Schmidt-Trucksss, C.-H. Liu, S.-H. Liu, Automated detection of the arterial inner walls of the common carotid artery based on dynamic B-mode signals, Sensors 10 (12) (2010) 10601–10619. [39] A. Delong, Y. Boykov, Globally optimal segmentation of multi-region objects, in: Proceedings of ICCV, 2009, pp. 285–292. [40] F.R. Schmidt, Y. Boykov, Hausdorff distance constraint for multi-surface segmentation, in: Proceedings of ECCV, Springer, Berlin, Heidelberg, 2012, pp. 598–611. [41] M. Sonka, V. Hlavac, R. Boyle, Image Processing, Analysis, and Machine Vision, Cengage Learning, Boston, 2007. [42] K. Li, S. Millington, X. Wu, D.Z. Chen, M. Sonka, Simultaneous segmentation of multiple closed surfaces using optimal graph searching, in: Proceedings of International Conference on Information Processing in Medical Imaging, Springer, Berlin, Heidelberg, 2005, pp. 406–417. [43] Y. Yin, X. Zhang, R. Williams, X. Wu, D.D. Anderson, M. Sonka, Logismoslayered optimal graph image segmentation of multiple objects and surfaces: cartilage segmentation in the knee joint, IEEE Trans. Med. Imaging 29 (12) (2010) 2023–2037. [44] L. Xu, B. Stojkovic, Y. Zhu, Q. Song, X. Wu, M. Sonka, J. Xu, Efficient algorithms for segmenting globally optimal and smooth multi-surfaces, in: Proceedings of 22th International Conference on Information Processing in Medical Imaging, 2011, pp. 208–220. [45] L. Xu, B. Stojkovic, H. Ding, Q. Song, X. Wu, M. Sonka, J. Xu, Faster segmentation algorithm for optical coherence tomography images with guaranteed smoothness, in: Proceedings of 2nd International Workshop on Machine Learning in Medical Imaging, 2011, pp. 308–316.

Zhengwang Wu received the B.S. degree in Computer Science from Xi'an Jiaotong University, Xi'an, China, in 2006. He is currently a Ph.D. student in Institute of Artificial Intelligence and Robotics, Xi'an Jiaotong University. His research interests include medical image processing, structural pattern recognition, and combinatory optimization.

Xiaoyi Jiang studied Computer Science at Peking University and received his Ph.D. and Venia Docendi (Habilitation) degree from University of Bern, Switzerland. He was an associate professor at Technical University of Berlin. Since 2002 he is a full professor of Computer Science at University of Münster, Germany. Currently, he is Editor-in-Chief of International Journal of Pattern Recognition and Artificial Intelligence and also serves on the advisory board and editorial board of several journals including Pattern Recognition, IEEE Transactions on Cybernetics, and Chinese Science Bulletin. He is a senior member of IEEE and fellow of IAPR.

390

Z. Wu et al. / Pattern Recognition 48 (2015) 380–390

Nanning Zheng graduated in 1975 from the Department of Electrical Engineering, Xi'an Jiaotong University, Xi'an, China, received the M.E. degree in information and control engineering from Xi'an Jiaotong University in 1981, and the Ph.D. degree in electrical engineering from Keio University, Japan, in 1985. He is currently a full professor and the Director of the Institute of Artificial Intelligence and Robotics at Xi'an Jiaotong University. His research interests include computer vision, pattern recognition, computational intelligence, image processing, and hardware implementation of intelligent systems. He is a member of the Chinese Academy Engineering and IEEE fellow.

Yuehu Liu received the B.S. and M.E. degrees in computer science at Xi'an Jiaotong University, China, in 1984 and 1989, respectively, and the Ph.D. degree in electrical engineering from Keio University, Japan, in 2000. He is currently a professor of the Institute of Artificial Intelligence and Robotics at Xi'an Jiaotong University. His research interests include computer vision, pattern recognition, computational intelligence, image processing.

Da-Chuan Cheng studied Electrical Engineering (B.S.) at Chinese Culture University, Taiwan, Biomedical Engineering (M.S. and Ph.D.) at National Cheng Kung University, Taiwan. He is an associate professor at the Department of Biomedical Imaging and Radiological Science of the China Medical University, Taiwan. His research interests include medical image processing and optimization.