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

(1)

Pattern Recognition ∎ (∎∎∎∎) ∎∎∎–∎∎∎

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

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

European Institute for Molecular Imaging, University of Münster, Germany

Cluster of Excellence EXC 1003, Cells in Motion, CiM, Münster, Germany

Department of Biomedical Imaging and Radiological Science, China Medical University, Taiwan

a r t i c l e i n f o Article history: Received 19 August 2013 Received in revised form 21 July 2014 Accepted 24 July 2014 Keywords: Generalized median Median surface

(2)

Graph search

Parameter handling

a b s t r a c t

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.

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 AU 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 algo-rithms 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: xjiang@uni-muenster.de (X. Jiang).

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 terrain-like 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.

(3)

Please cite this article as: Z. Wu, et al., Exact solution to median surface problem using 3D graph search and application to parameter space exploration, Pattern Recognition (2014), http://dx.doi.org/10.1016/j.patcog.2014.07.019i

(4)

(5)

(6)

(7)

(8)

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 position-sensitive weighting of the input surfaces. In Section 6 a scheme is presented for

(9)

experimental results to demonstrate the usability of median sur-face 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 representa-tion space U and a distance function dðp; qÞ to measure the dissimilarity between any two objects p; qAU. The generalized median p AU of S minimizes the sum of distances to all objects from S, i.e.

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

pAU_{q AS}

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

¼ arg min K 1 S ∑ d s S iÞ

(10)

s AUS

ð ;

ð Þ

i ¼ 1

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:

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

p AS qAS

The set median may serve as an approximative solution for the generalized median or an initialization for iterative generalized median computation. This is justified by the fact ∑q ASdðp^ ; qÞ= ∑qASdðp; qÞ r2_ 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 construc-tion) [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 compu-tational 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.

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

MN

ð2Þ

(11)

x ¼ 1 y ¼ 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 position-sensitive 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 compu-tation 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.

Please cite this article as: Z. Wu, et al., Exact solution to median surface problem using 3D graph search and application to parameter space exploration, Pattern Recognition (2014),

http://dx.doi.org/10.1016/j.patcog.2014.07.019i

3

Z. Wu et al. / Pattern Recognition ∎ (∎∎∎∎) ∎∎∎–∎∎∎

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 segmenta-tion result and the corresponding ground truth.

(12)

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:

¼ arg min K S ∑ d s S Þ s AUS ð ;

(13)

i i ¼ 1 ¼ arg min K M N

(14)

∑ ∑ ∑ w xy _ ρ s x y S ið x y ÞÞ s AUS ð ð ; Þ; ; i ¼ 1 x ¼ 1 y ¼ 1 ¼ arg min M

(15)

N K ∑ ∑ w _ ∑ ρ z ¼ s x y Þ; S x y ÞÞ s AUS xy ð ð ; ið ; x ¼ 1 y ¼ 1

(16)

i ¼ 1 cxyz ¼ arg min M N | fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} ∑ ∑ c xyz s AUS

(17)

x ¼ 1 y ¼ 1

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 para-meter 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 indepen-dently 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 explora-tion problem for image segmentation based on some baseline

(18)

segmentation algorithm. In this sense our approach is on a meta-level above the

segmentation level. The median surface computa-tion 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

|fflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflffl}CðsÞ

¼ arg min CðsÞ ð3Þ

s AUS

A candidate solution is a surface s AUS

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 41).

Note that Eq. (3) can be viewed as a first order MRF minimum energy problem, i.e.,

(19)

v AV ðu;vÞ AE

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; XvÞ is convex and Pðu; v;

Xu; XvÞ 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] corre-sponds 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

(20)

4 Z. Wu et al. / Pattern Recognition ∎ (∎∎∎∎) ∎∎∎–∎∎∎

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 c xy;z _ 1 ; z 1 ð4Þ Wðx; y; zÞ ¼ ( cxyz;_ z 4 1 ¼

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 41) has a direct edge to the node V ðx; y; z _ 1Þ below it, i.e.,

Ea ¼ f oVðx; y; zÞ; Vðx; y; z _ 1Þ 4 j z 41g

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 interposi-tion edges are defined by

8 f oV x; y; z ; V x 1; y; max 0; z _ xÞÞ 4 j xA 1; …; M _ 1g; z AZg [ E r oVð_{x; y; z} Þ ; V ð xþ 1; y; max ð 0; z x ÞÞ 4

(21)

xAf_{2; …; M ; z AZ} g [ > f ð Þ ð _ ð _ y 4 j f g ¼ > oV x; y; z ; V x; y 1; max 0; z y A 1; …; N 1 ; z AZ >

(22)

> < f ð Þ ð þ ð _ ÞÞ j f _ g g [ >

f oVðx; y; zÞ; Vðx; y _ 1; maxð0; z _ yÞÞ 4

j y Af2; …; Ng; z AZg

(23)

>

(24)

:

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 rNð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ðn2mcÞ, 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

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.

(25)

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 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 con-struction 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Þ ASi can be regarded as a sample

of a random variable Fxy, which is assumed to

form a Gaussian distribution. As an illustra-tion,

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Þ ASi, 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 illu-strated 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 correspond-ing 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

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

(26)

5

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 per-image basis, which may be approached, but cannot be achieved by segmentation algorithms in general. In all our tests the compar-ison between two surfaces, e.g. a segmented surface and a ground truth, was done by computing the average L1

(27)

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 best-performing 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 perfor-mance 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 position-sensitive weighting and presented a weighting scheme. A test series was performed to experimentally confirm the superior performance of weighting. Again, the median surfaces are com-pared 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 gener-ated 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

(28)

(29)

(30)

(31)

2.8 GHz

with 16G memory. The program was written in C

(Boykov–Kolmogorov algorithm) and Matlab (other parts). The

(32)

resolution of the subpart of a video, in which the near or the far

(33)

(34)

graph search algorithm is quite fast.

(35)

(36)

Fig. 3. Sonographic video. (a) Along with the time axis the 2D images build a 3D

Our approach to deal with the parameter problem is in line

with the ensemble paradigm which is a powerful boosting

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

technique in machine learning and pattern recognition [2,3]: by

referred to the web version of this paper.)

(37)

classi-fication can be achieved. In our case an ensemble of K surfaces is

Table 1

averaged to achieve an “optimal” segmentation. This principle has

Statistics of the 100 segmented surfaces: [minimum, maximum], mean, standard

been successfully validated for 2D contour detection [4] and image

deviation (unit: pixels).

segmentation [21,22] before. In the current work this ensemble

Video

#Images Near wall

(38)

Far wall

approach is demonstrated for parameter space exploration in the

(39)

1 86

[0.46, 1.00], 0.78, 0.17 [0.51, 0.71], 0.62, 0.05

In Section 4 we have discussed several parameter exploration

2 86

[0.41, 1.35], 0.82, 0.33 [0.47, 0.75], 0.61, 0.08

methods [27–30]. However, their common critical disadvantage is

3 111

[0.89, 2.82], 2.00, 0.73 [0.63, 0.96], 0.73, 0.09

(40)

4 73

[0.45, 1.76], 1.08, 0.56 [0.45, 1.45], 0.75, 0.22

(41)

in many situations. Our intention in the application part of this

Table 2

(42)

at least, without assuming the ground truth. Due to the very

Comparison with the ground truth (unit: pixels).

different nature of our unsupervised ensemble approach and other

Video

(43)

Near wall

Far wall

supervised methods we do not perform direct comparison with

them. However, we have shown that our median segmentation

(44)

Median vs. GT 0.63 0.57 results reach the quality of best-performing parameter setting BP vs. GT 0.61

(45)

0.55

(46)

2

Median vs. GT

0.48

0.50

cannot achieve the performance of best-performing parameter

BP vs. GT

0.46

0.49

setting, this result is an indirect proof of the power of our

3

Median vs. GT

(47)

0.67

approach compared to the existing supervised methods.

In our tests the involved dynamic programming approach from

BP vs. GT

0.96

(48)

[38] has two parameters only. In the case of more parameters the

4

Median vs. GT

0.48

0.54

number of sampled parameter settings certainly increases. One

BP vs. GT

(49)

0.51

can either fix some non-critical parameters and sample a small

(50)

(non-the studied parameter interval. The performance is

given in critical)

parameters relatively sparsely. Another solution is to

compute the K input surfaces in a parallel manner. It is important

Table 4. In the first test scenario similar behavior as in Table 2

to emphasize that this is a general issue related to many, if not all,

can be observed. The median computation basically achieves the

parameter learning methods from the literature.

(51)

The application of our median surface computation intends to

surprisingly, with the increasing percentage of non-optimal

para-alleviate the parameter problem. It is therefore important to make

meter settings, the performance of median computation slightly

sure that our approach does not introduce new parameter(s). This

drops. But even in the last test scenario the performance drop is

is the

case since all steps

described are parameter-free. For

not dramatic, indicating the potential of robustness by using the

(52)

median surface approach in dealing with the segmentation

para-on functipara-on fitting. Note that the cpara-ontinuity values x and

y in

meter problem.

(53)

com-The robustness is mainly due to the use of L1 norm. In robust

putation algorithm. Indeed, any baseline surface detection

algo-statistics [26] it is known that L1 is robust and can tolerate up to

rithm, which delivers the input for our median surface 50% outliers. Any other robust dissimilarity function ρ with a

(54)

computation, needs these values as well. They mean the degree

possibly

high breakdown point will

make the median surface

of surface continuity a user wishes to achieve and thus have to be

(55)

specified by the user, most likely in

an

application-dependent

Since our approach is an ensemble one, the ensemble quality

(56)

(57)

reasonable parameter subspace for sampling. In practice, this is

mostly not difficult to achieve. Together with a robust dissimilarity 8. Multiple surfaces median computation

(58)

tolerated by the median surface computation and an overall good Due to the imperfections of medical imaging techniques,

quality of median surface can be reached.

insufficient or even deficient image-derived information may be

7

(59)

(60)

Table 3

Definition 2. Two

coupled

terrain-like

surfaces contain two

Comparison weighted vs. unweighted (unit: pixels).

terrain-like single surfaces: f i : X _ Y-Z with X ¼ f1; 2; …; Mg,

(61)

Video Variant Near wall Far wall

Y ¼ f1; 2; …; Ng, and Z ¼ f1; 2; …; Lg, i¼ 1,2. The pre-defined

cou-pling constraints δl and δu specify, respectively, the minimum and

1 Weighted 0.63 0.57

(62)

maximum distances between the two non-intersecting surfaces,

Unweighted 0.67 0.57

i.e. δl rjf 1ðx; yÞ _ f 2ðx; yÞj rδu.

(63)

This definition can be easily extended to the case of K (42)

Unweighted 0.54 0.53

(64)

for each pair of neighboring surfaces fi and fj the related

con-Unweighted 1.08 0.68 straints δlij and δuij. 4 Weighted 0.48 0.54

(65)

coupled surfaces. In this case the individual graphs Gi ¼ ðV i; Ei; Unweighted 0.52 0.53 WiÞ;

(66)

available. This insufficiency can be alleviated by using clues from G ¼ ðV1 [ ⋯ [ VK ; E1 [ ⋯ [ EK [ Es; W1 [ ⋯ [ WK Þ

The extra intersurface edges in Es model the relations between

other mutually related boundaries or surfaces. Indeed,

co-surfaces. As done in [25], we specify the construction of Es for

detection of multiple coupled surfaces frequently yields superior

(67)

performance compared to the common single-surface detection

approaches [31,25,38–40]. Thus, it is essential to extend the single two surfaces. Supposing f1 is above f2, then

(68)

terrain-like surface median problem to the case of multiple Es x; y; z ; V2 x; y; z δu 4 z Zδu g δ l[ coupled surfaces. ¼ 8f oV2 ð_{x; y; z}Þ ; V1 ð_{x; y; z}_ δl Þ 4 jz oL

(69)

> f ð Þ ð _ Þ j _ g [ 8.1. Algorithmic extension < oV1 ð 0; 0; δl Þ ; V 2 0; 0; 0 Þ 4 g

(70)

> f

(71)

Note:that the formulation above is given for non-intersecting

We start with the definition of multiple surfaces (defined in the surfaces. The needed extension for intersecting surfaces can be

same way as in [25]).

found in [25], also the details of computing the K optimal surfaces.

(72)

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

Table 4

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

Boykov–Kolmogorov algorithm (Oðn2mcÞ) increases at least by a factor of K3 (because of K _ n nodes and K _ m plus extra intersur-face edges in the enlarged graph).

8.2. Experimental results Comparison type Near wall Far wall Video 1 50 best þ 10 worst Median vs. GT 0.57 0.55 Bp vs. GT 0.61 0.55 40 best þ 20 worst Median vs. GT 0.56 0.55 BP vs. GT 0.61 0.55 30 best þ 30 worst Median vs. GT 0.68 0.58 BP vs. GT 0.61 0.55 Video 2 50 best þ 10 worst Median vs. GT 0.44 0.49 BP vs. GT 0.46 0.49 40 best þ 20 worst Median vs. GT 0.44 0.49 BP vs. GT 0.46

(73)

0.49 30 best þ 30 worst Median vs. GT 0.51 0.50 BP vs. GT 0.46 0.49

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 experi-mental 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 seg-mentation 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

The 3D graph search algorithm in [25] designed for terrain-like

Video 3

surface detection has several variants of extensions which further

50 best þ 10 worst Median vs. GT 0.97 0.66

broaden the usability of median surface computation. The

exten-BP vs. GT 0.96 0.66

(74)

sion to the case of multiple surfaces in already discussed in Section

8. In this section we briefly discuss three additional extensions and

BP vs. GT 0.96 0.66

another issue related to computational cost.

30 best þ 30 worst Median vs. GT 1.10 0.67 9.1. Tube-like surfaces BP vs. GT 0.96 0.66 Video 4

(75)

Sometimes, the desired surface is required to be wraparound

along some direction (say, the x- or y-axis). For instance, this is a

BP vs. GT 0.46 0.51

common situation in medical image analysis to detect tube-like

(cylindrical) surfaces. In this case a tubular surface is first unfolded

(76)

BP vs. GT 0.46 0.51 30 best þ 30 worst Median vs. GT 0.48 0.53

[41]. Then, the boundary of the tubular surface in the original

image corresponds to a terrain-like surface to be detected in the

BP vs. GT 0.46 0.51

(77)

be applied. Moreover, the first and last rows along the unfolding

Table 5

plane should satisfy the continuity constraints as well, which is

guaranteed by a modified graph construction. Median surface

Computational time (seconds).

computation is directly applicable in this case.

Video #images Near wall Far wall Resolution Time Resolution

(78)

time 1 86 45 _ 121 3.2 49 _ 121 4.0 2 86 47 _ 163 6.3 46 _ 163 6.8 3 111 72 _ 188 21.8 56 _ 188 19.1 4 73 40 _ 126 2.4 45 _ 126 7.5 9.2. Closed surfaces

Besides the tube-like surfaces, closed surfaces are also fre-quently 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

(79)

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

9

Table 6

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

Video Comparison type Near wall

(80)

Intima Adventitia Intima Adventitia 1 Median vs. GT 0.63 0.51 0.57 0.47 BP vs. GT 0.61 0.48 0.55 0.47 2 Median vs. GT 0.48 0.80 0.50 0.81 BP vs. GT 0.46 0.80 0.49 0.80 3 Median vs. GT 0.98 1.43 0.67 0.76 BP vs. GT

(81)

0.96 1.46 0.66 0.75 4 Median vs. GT 0.48 0.74 0.54 0.53 BP vs. GT 0.46 0.69 0.51 0.50 Table 7

(82)

video #images Near wall Far wall Resolution Time Resolution Time 1 86 45 _ 121 93.44 49 _ 121 96.62 2 86 47 _ 163 172.79 46 _ 163 164.89 3 111 72 _ 188 405.76 56 _ 188 379.08 4 73 40 _ 126 74.36 45 _ 126 71.59

(83)

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 extrac-tion problems. Although it is quite natural for the variational methods such as active contours and level set approaches, impos-ing 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 single-surface 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 41 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 con-cretized 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 explora-tion 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 alter-nate 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 meta-level 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

(84)

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 pre-processing 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 high-level 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.

(85)

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

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.

L. Rokach, Ensemble-based classifiers, Artif. Intell. Rev. 33 (1–2) (2010) 1–39.

Z.-H. Zhou, Ensemble Methods: Foundations and Algorithms, CRC Press, Boca Raton, London, New York, 2012.

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.

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.

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 Recogni-tion, Lecture Notes in Computer Science, vol. 7877, Springer, Berlin, Heidel-berg, 2013, pp. 234–243.

D. Gusfield, Algorithms on Strings, Trees, and Sequences: Computer Science and Computational Biology, Cambridge University Press, Cambridge, UK, 1997.

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.

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.

R.I. Hartley, J. Trumpf, Y. Dai, H. Li, Rotation averaging, Int. J. Comput. Vis. 103

(3) (2013) 267–305.

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.

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.

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.

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.

S. Vega-Pons, J. Ruiz-Shulcloper, A survey of clustering ensemble algorithms, Int. J. Pattern Recognit. Artif. Intell. 25 (3) (2011) 337–372.

H. Alizadeh, B. Minaei-Bidgoli, H. Parvin, Optimizing fuzzy cluster ensemble in string representation, Int. J. Pattern Recognit. Artif. Intell. 27 (2) (2013) 1350005.

A. Feragen, S. Hauberg, M. Nielsen, F. Lauze, Means in spaces of tree-like shapes, in: Proceedings of ICCV, 2011, pp. 736–746.

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.

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. 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.

(86)

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.

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.

L. Franek, X. Jiang, Evolutionary weighted mean based framework for general-ized 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.

L. Franek, X. Jiang, Ensemble clustering by means of clustering embedding in vector spaces, Pattern Recognit. 47 (2) (2014) 833–842.

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.

C. Stewart, Robust parameter estimation in computer vision, SIAM Rev. 41 (3) (1999) 513–537.

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.

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.

L. Franek, X. Jiang, Alternating scheme for supervised parameter learning with application to image segmentation, in: Proceedings of CAIP, 2011, pp. 118–125.

L. Franek, X. Jiang, Orthogonal design of experiments for parameter learning in image segmentation, Signal Process. 93 (6) (2013) 1694–1704.

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.

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.

D. Greig, B. Porteous, A. Seheult, Exact maximum a posteriori estimation for binary images, J. R. Stat. Soc. Ser. B: Methodol. (1989) 271–279.

D. Hochbaum, An efficient algorithm for image segmentation, Markov random fields and related problems, J. ACM 48 (4) (2001) 686–701.

H. Ishikawa, D. Geiger, Segmentation by grouping junctions, in: Proceedings of CVPR, 1998, pp. 125–131.

H. Ishikawa, Exact optimization for Markov random fields with convex priors, IEEE Trans. Pattern Anal. Mach. Intell. 25 (10) (2003) 1333–1336.

X. Wu, D.Z. Chen, Optimal net surface problems with applications, in: Automata, Languages and Programming, Springer, Berlin, Heidelberg, 2002, pp. 1029–1042.

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.

A. Delong, Y. Boykov, Globally optimal segmentation of multi-region objects, in: Proceedings of ICCV, 2009, pp. 285–292.

F.R. Schmidt, Y. Boykov, Hausdorff distance constraint for multi-surface segmentation, in: Proceedings of ECCV, Springer, Berlin, Heidelberg, 2012, pp. 598–611.

M. Sonka, V. Hlavac, R. Boyle, Image Processing, Analysis, and Machine Vision, Cengage Learning, Boston, 2007. 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.

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.

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.

L. Xu, B. Stojkovic, H. Ding, Q. Song, X. Wu, M. Sonka, J. Xu, Faster segmenta-tion 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.

(87)

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.

11

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.

(88)

(89)