A simple boundary reinforcement technique for segmentation without prior

A simple boundary reinforcement technique for segmentation without prior

Accepted Manuscript A simple boundary reinforcement technique for segmentation without prior Nóirín Duggan, Egil Bae, Edward Jones, Martin Glavin, Lum...

2MB Sizes 3 Downloads 71 Views

Accepted Manuscript A simple boundary reinforcement technique for segmentation without prior Nóirín Duggan, Egil Bae, Edward Jones, Martin Glavin, Luminita Vese PII: DOI: Reference:

S0167-8655(14)00132-9 http://dx.doi.org/10.1016/j.patrec.2014.04.014 PATREC 6003

To appear in:

Pattern Recognition Letters

Please cite this article as: Duggan, N., Bae, E., Jones, E., Glavin, M., Vese, L., A simple boundary reinforcement technique for segmentation without prior, Pattern Recognition Letters (2014), doi: http://dx.doi.org/10.1016/ j.patrec.2014.04.014

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

A simple boundary reinforcement technique for segmentation without prior

1

2

3 4

N´oir´ın Duggana,b,1,, Egil Baeb , Edward Jonesa , Martin Glavina , Luminita Veseb a

8

Electrical & Electronic Engineering, National University of Ireland, Galway, Galway, Ireland b Department of Mathematics, University of California Los Angeles, 405 Hilgard Avenue, Los Angeles, CA 90095-1555

9

Abstract

5 6 7

10

Accurate boundary detection is a critical step of the image segmentation

11

process. While most edge detectors rely on the presence of strong intensity

12

gradients, this criteria can limit robustness in many real world cases. In

13

this work we propose a scheme which makes use of a combination of both

14

global and local segmentation methods to capture the boundary of target

15

objects in low contrast images. This approach has several advantages: The

16

globally convex segmentation scheme is immune from initial conditions and

17

is easily adapted to the data. The addition of a segmentation scheme based

18

on local curve evolution produces a solution which is shown to help preserve

19

topology between the initial and target shape, a property lacking in globally

20

convex segmentation schemes. Experimental results show that the proposed

21

method achieves enhanced performance compared to classical data-driven

22

segmentation schemes proposed in the literature.

23

Keywords: image segmentation, convex optimization, variational methods,

24

edge detection, echocardiography images

Email address: [email protected] (N´oir´ın Duggan ) Preprint submitted to Elsevier

April 24, 2014

25

1. Introduction

26

Image segmentation is arguably one of the most important tasks in com-

27

puter vision and has been an active area of research over the last three

28

decades. Efforts directed at this task have included variational approaches,

29

statistical and more recently combinatorial methods. Many of the most suc-

30

cessful segmentation methods have been formulated as energy minimization

31

problems in which the key object detection criteria are incorporated into the

32

energy functional.

33

One of the first proposals to cast image segmentation as an energy min-

34

imization problem was the approach of Kass et al. [11] who computed the

35

segmentation of a given image by evolving curves in the direction of the

36

negative energy gradient using appropriate partial differential equations. In

37

this approach, boundaries are detected using the strength of the gradient

38

at each pixel. In [11] the curve is represented explicitly, which can lead to

39

some drawbacks, in particular in order to avoid self-intersection and overlap

40

of contour points, reparameterization is required. Similarly, sophisticated

41

reparameterization schemes are needed to handle topological changes, which

42

are necessary for segmenting multiple objects or objects with an unknown

43

topology. To overcome the need for reparameterization, techniques based on implicit curve evolution theory [15, 3] allow for motion based on geometric measures such as unit normal and curvature. To obtain a new length constraint which is independent of parameterization, Caselles et al. [3] and Kichenassamy et al. [12] simultaneously proposed the implicit geodesic active contour, the

2

energy functional of which is given by: ∫ 1 E(C) = g(|∇I(C(q))|)|C ′ (q)|dq

(1)

0

where g =

1 with p = 1 or 2 1 + β|∇I|p

44

If we let Ω denote the image domain, where Ω ⊂ RN then in the above

45

equation, the image is represented by, I : Ω → R, C(q) : [0, 1] → Ω represents

46

the curve, while β controls the sharpness of the detected edges. In the level set

47

method of Osher and Sethian [15], the curve is represented implicitly as the

48

zero level line of some embedding function ϕ : Ω → R: C = {x ∈ Ω | ϕ(x) =

49

0} where ϕ(x) ≤ 0 in the interior of the curve and ϕ(x) > 0 in the exterior.

50

The level set method evolves a curve by updating the level set function at

51

fixed coordinates through time, rather than tracking a curve through time, as

52

in a parametric setting. In an implicit formulation, both topological changes

53

and also extension to higher dimensions are handled naturally.

54

Another energy minimization approach for detecting edges in an image

55

is the Mumford-Shah model [14], which in its most general form seeks both

56

an edge set and an approximation of the image which is smooth everywhere

57

except across the edge set. There has been a considerable amount of research

58

on computing minimizers by numerical algorithms. In [5, 23] the authors

59

presented a level-set formulation of the piecewise constant variant of the

60

Mumford-Shah model. Considering an image with two regions, the object to

61

be segmented is denoted S and the background denoted Ω\S, the authors

62

proposed the following model: ∫ ∫ 2 min |I(x) − c1 | dx + |I(x) − c2 |2 dx + ν |∂S| . S,c1 ,c2

Ω\S

S

3

(2)

The last term of the energy functional is the length of the boundary of S weighted by the parameter ν, while the first terms represent the data-fidelity where c1 and c2 are two scalars that attempt to approximate the image in the interior and exterior region. In order to find a minimizer by a numerical algorithm, it was proposed in [5, 23] to represent the partition with a level set function, resulting in the following problem: (∫ min |I(x, y) − c1 |2 H(ϕ)dx ϕ,c1 ,c2 ∫Ω + |I(x, y) − c2 |2 (1 − H(ϕ))dx Ω ) ∫ +ν |∇H(ϕ)|dx

(3)

Ω 63

where H is the Heaviside function, defined such that H(z) = 1 if z ≥

64

1 and H(z) = 0 if z < 0 and is used to select either the interior or exte-

65

rior of the curve. This work is known as Active Contours without Edges

66

(ACWE) or Chan-Vese (CV) method. More recently, efficient convex opti-

67

mization algorithms have been developed which result in global minimizers,

68

e.g. the Split Bregman algorithm [9] and continuous max-flow [26, 27].

69

To increase robustness, schemes have been proposed which combine region

70

and edge terms [16, 1]. In [16], an energy minimization approach for tracking

71

objects was proposed, where regional and edge information were integrated

72

in separate terms of an energy functional, in addition to a motion based term.

73

In order to solve the resulting minimization problem, the level set method

74

was applied. The quality of segmentation often depends on the quality of the

75

detected edges. Bresson et al. [1] proposed a method which combines the

76

data fidelity term of the ACWE model and the edge attraction term of the

77

geodesic active contour (GAC) model in a global minimization framework. 4

78

Another line of research attempts to minimize the full Mumford-Shah

79

model with curves that contain endpoints in order to detect edges. The

80

Mumford-Shah model and the GAC model were extended to include more

81

general edge sets in [22, 21]. In [24], the authors represented the edges by

82

a binary function, leading to a non-convex minimization problem closely in

83

spirit to the phase field representation. A convex relaxation of the piecewise

84

smooth Mumford-Shah model was proposed in [18], by lifting the the problem

85

to a higher dimensional space. Although computationally expensive, this

86

approach has the advantage of converging to a close approximation of a global

87

minimizer, without getting stuck in a local minimizer. In [2] a generalization

88

of the Mumford-Shah model was proposed with an extra term which enforced

89

closeness between the gradient of the reconstructed image and the original

90

image. The motivation behind this model is to better reconstruct edges with

91

low contrast changes. A numerical algorithm based on the level set method

92

was proposed. Closely related is also the diffusion equation of Perona and

93

Malik [17], which converges to a local minimizer of the Mumford-Shah model

94

if an appropriate data fitting term is added.

95

In Rajpoot et. al [19] an intensity invariant, real-time method was pro-

96

posed to extract boundary information by analyzing the monogenic signal

97

[7] . The proposed filter uses a local phase based method to extract edge

98

information. The method produces a finer edge extraction than gradient de-

99

tectors; however, due to its intensity invariance, tends to over-segment the

100

image. Another approach for generating feature detectors proposed for ap-

101

plication to echocardiography, is the method of Mulet-Parada et. al [13],

102

a 2D+T method, in which it was demonstrated that a local phase-based

5

103

approach produces more accurate edge results than conventional gradient

104

magnitude methods.

105

The problem we addressed in the current work was how to increase seg-

106

mentation robustness using a purely data driven approach. Incorporating

107

prior information into a segmentation scheme is of course an advantage if

108

such information exists, however this is not always the case; with this in

109

mind, in the current work we sought to develop a method for images with

110

weak boundaries by looking at new ways to incorporate pre-existing informa-

111

tion. In the preliminary conference paper [6], a combination technique was

112

presented using Geometric Split Bregman method with a topology preserv-

113

ing level set technique for application to echocardiography with the similar

114

objective of detecting low-contrast boundaries. The method described here

115

similarly proposes a sequential model using a global segmentation scheme,

116

however in this work we make use of the alternate formulation using continu-

117

ous maximum flow method while using the Geodesic Active Contour method

118

to capture the final boundary. Similarly to [24, 2], we also apply a variation

119

of the Mumford Shah model to improve the edge detection step. While in

120

[24] the goal was to improve the ‘standard’ canny edge model, our proposal is

121

to delineate structures with low contrast boundaries by reinforcing the weak

122

edges.

123

To summarize our contributions, in this paper we propose:

124

1. A new technique for generating an edge detector function which exploits

125

the natural properties of the Mumford-Shah functional to effectively

126

reinforce weak boundaries.

127

2. A combination of the proposed edge detector with the GAC method 6

128

and demonstrate its applicability to images with weak edges, such as

129

echocardiography medical images.

130

The remainder of the paper is organized as follows: In Section 2, we

131

describe the proposed segmentation model. Section 3 contains experimental

132

results on both medical and non-medical real images. Section 4 presents the

133

discussion of the results and Section 5 concludes the paper.

134

2. Methodology

135

In this paper we propose a sequential model which addresses some of the

136

weaknesses of these previous methods. The type of images focused on in this

137

work often have weak edges and boundaries. Therefore, an edge detection

138

method based directly on the image gradients is not expected to work well.

139

In order to obtain a clearer edge map, we first approximate the image by a

140

piecewise constant function, by efficiently computing a global minimizer to

141

the model (2) with the Continuous Maximum Flow (CMF) algorithm [26, 27]

142

Then in order to obtain a single curve which captures the boundary of the

143

object, we apply the geodesic active contour (GAC) method where the edge

144

attraction term is constructed from the edge map in the first step.

145

2.1. Computation of Edge Detector from Global Segmentation output

146

The first step of our method is to extract a rough edge map by using the

147

ACWE model with two regions. In recent work, efficient algorithms have

148

been proposed for computing global minima to this model. In [4] it was

149

shown that (2) can be exactly minimized via the convex problem

7

∫ min ϕ(x)∈[0,1]

|I(x) − c1 |2 ϕ(x) + |I(x) − c2 |2 (1 − ϕ(x))dx Ω ∫ +ν |∇ϕ(x)|dx .

(4)

Ω 150

It was shown that if ϕ∗ is a minimizer of (4) and t ∈ (0, 1] is any threshold

151

level, the partition S = {x ∈ Ω : ϕ(x) ≥ t}, Ω\S = {x ∈ Ω : ϕ(x) < t} is a

152

153

global minimizer to the model (2). The binary function   1 , ϕ(x) ≥ t ϕt (x) := ,  0 , ϕ(x) < t

(5)

is the characteristic function of the region S.

154

A very efficient algorithm was proposed for solving the problem (4) in

155

[26, 27]. The basic idea is to derive an augmented Lagragian algorithm

156

based on the following dual problem of (4): ∫ sup ps (x) dx ps ,pt ,p

(6)



subject to |p(x)|2 ≤ ν,

∀x ∈ Ω;

(7)

ps (x) ≤ |I(x) − c1 |2 ,

∀x ∈ Ω;

(8)

pt (x) ≤ |I(x) − c2 |2 ,

∀x ∈ Ω;

(9)

on ∂Ω.

(10)

a.e. x ∈ Ω.

(11)

p · n = 0, div p(x) − ps (x) + pt (x) = 0, 157

where the dual variables are scalar and vector functions: ps , pt : Ω 7→ R and

158

p : Ω 7→ RN , where N is the dimension of the image domain Ω. The dual 8

159

problem can be interpreted as a maximum flow problem over a continuous

160

domain, therefore the following algorithm was referred to as a continuous

161

max-flow algorithm (CMF). In order to solve the maximization problem, a Lagrange multiplier ϕ was introduced for the constraint (11), and the augmented Lagrangian functional was formulated as follows:



Lc (ps , pt , p, ϕ) :=

ps dx ∫ ( ) + ϕ div p − ps + pt dx Ω



− 162

163

164

where for a function a, ∥a∥2 =

∫ Ω

c ∥div p − ps + pt ∥2 2

(12)

|a(x)|2 dx. Assume the problem has been

discretized such that ps , pt and ϕ are defined for each pixel in the discrete ∫ image domain Ω and , ∇ and div = −∇∗ are some discrete integration

165

and differential operators. In all experiments in this paper, we have used

166

a mimetic discretization of the differential operators [10]. An augmented

167

lagrangian method could be applied by alternatively maximizing Lc for the

168

dual variables ps , pt , p with constraints (7)-(10) and update the lagrange mul-

169

tiplier ϕ as follows:

170

171

Initialize p1s , p1t , p1 and ϕ1 . For k = 1, ... until convergence: •

(∫ pk+1 := s

arg max ps (x)≤|I(x)−c1 |2 ∀x∈Ω

ps dx Ω

2 c − ps − pt k − div pk + ϕk /c 2 9

)

which can easily be computed pointwise in closed form.

172



2 c pk+1 := arg max − div p − ps k+1 + pt k − ϕk /c , ∥p∥∞ ≤ν 2 173

where ∥p∥∞ = supx∈Ω |p(x)|2 . This problem can either be solved itera-

174

tively or approximately in one step via a simple linearization [27]. In

175

our implementation we used the linearization. • pk+1 := t

arg max pt (x)≤|I(x)−c2 |2 ∀x∈Ω



2 c

pt − ps k+1 + div pk+1 − ϕk /c 2

This problem can also easily be computed in closed form pointwise.

176

• ϕk+1 = ϕk − c (div pk+1 − pk+1 + pk+1 ); s t 177

• Set k = k + 1 and repeat.

178

The output ϕ at convergence will be a solution to (4) and one can obtain

179

a partition which solves (2) by the thresholding procedure described above.

180

More details on implementation of the algorithm can be found in the paper

181

[27].

182

2.2. Boundary Localization through local curve evolution This initial segmentation corresponds to a simplified binary representation of the image, in which 2 prominent regions are identified and represented by binary functions with an intensity of either 1 or 0. However, the segmentation result obtained in this first step is not expected to be a good 10

representation of the object we wish to detect. Globally optimal solutions of the region based model (2) often have a complicated topology, which do not represent a simple connected object with a closed boundary. In the second step, we make use of the edge of the binary image from section 2.1 to create a more prominent edge attraction force in the geodesic active contour model defined in equation (1). The evolution equation of the GAC model is given by ∂C = (κg − ⟨∇g, N ⟩)N ∂t 183

(13)

where κ is the Euclidean curvature, and N is the unit inward normal. In our proposed method, the edge map g, is generated using the clean binary output of step 1 using: g(∇ICM F ) =

184

1 1 + |G ∗ ∇ICM F |

(14)

where G represents a Gaussian Kernel.

185

Compared to the globally optimal output of Section 2.1, the final out-

186

put of the combined approach is expected to have a much simpler topology,

187

consisting in most cases of one connected component.

188

The GAC method was implemented using the AOS (Additive Operator

189

Splitting) scheme for solving PDE equations in the form of ut = div(g ∗ u)).

190

The details of the implementation can be found in [25].

191

To summarize our approach, we first create a binary approximation of

192

the original image and extract the major edge set using CMF method. This

193

information is used to yield a better edge attraction force in the GAC method,

194

which locates the object boundary. 11

195

In the following section, the experimental results of the proposed method

196

are presented.

197

3. Experiments

198

To test the method a sample combination of medical and non-medical real

199

world data was selected which naturally exhibit low contrast boundaries. In

200

the case of the medical data, a small test sample of 5 frames each were cho-

201

sen from 5 low quality echocardiography datasets. The other test images

202

consist of an image of a galaxy and also an image of a swarm of birds. The

203

method was validated by comparing it to two classical data driven segmenta-

204

tion models: ACWE approach and the standard GAC model. In the medical

205

image test, we also compared the method to the local-phase based feature

206

detection approach of [19, 20], this method having been designed specifically

207

for echocardiography. For each test case and dataset, parameters controlling

208

length as well as the initial estimate of regional intensity were chosen based

209

on the data, parameters were similarly optimized for each of the comparison

210

methods. For the experiments on medical data, we kept the intensity esti-

211

mates, c1 and c2 from equation (2), constant over the course of the energy

212

minimization (minimizing with respect to length). For the experiments, c1

213

was chosen to be in the range 16e-3 to 5e-2, with c2 in the range between

214

0.12 to 0.28, while best results were obtained with the length parameter ν,

215

between 1e − 3 and 0.5.

216

For the ACWE method the following smooth weight values were used

217

[15, 2.03, 1.23, 18], while the contrast parameters were set to [950, 318.75, 78.75, 1200].

218

For the standard GAC method the length parameter was set to 10 through12

219

out, while the expansion term was chosen to be between 0.6 and 0.9. For

220

the local-phase method, the following values were chosen for the wavelength

221

parameter: [42, 7, 16, 28]. Each method was manually initialized by plac-

222

ing a contour in each frame. Sensitivity to variance in parameter choice is

223

presented in Table 2.

224

For the medical data, quantitative results showing the similarity between

225

manual delineations carried out by 2 expert observers and each of the meth-

226

ods is presented in Table 1. A sample of the final results obtained on the

227

echocardiography data, together with those obtained with the GAC and

228

ACWE and local-phase based methods is shown in Figures 4 and 5 while

229

results obtained from non-medical data are shown in Figures 1 and 2. Figure

230

3 shows an example of the CMF output and its corresponding edge map.

231

The average CPU time required to compute the global segmentation (step

232

1) was found to be in the order of 0.6 seconds, while the full curve evolution

233

method took on average 9.7 seconds. In total, the proposed method took

234

an average of 10.5 seconds. If we compare this with the reference methods,

235

the Local Phase Method took on average 10.1 seconds, the Geodesic method

236

14.7 seconds and the ACWE method took 5 seconds. The experiments were

237

carried out both in C and in Matlab, on a 64-bit, Intel Core 2 Dual CPU

238

3.00GHz processor with 4 GB of installed memory. The average image size

239

was 376x376. The parallel processing capability of the processor was not

240

used.

241

Table 1 presents the average accuracy results for each dataset. With

242

respect to the Hausdorff measurement, the CMF + GAC obtained an average

243

Hausdorff distance of 4.48 mm from the reference contour, which was an

13

244

improvement on each of the other methods, with the ACWE, GAC and

245

local-phase method obtaining averages of 6.29 mm, 4.7 mm and 6.77 mm

246

respectively. The second table presents Mean Absolute Distance results.

247

In this case, the CMF + GAC achieved an average distance of 1.31 mm

248

compared to 1.9 mm, 1.84 mm and 1.9 mm respectively for the ACWE,

249

GAC and local-phase methods. The third table presents the result of the

250

dice overlap metric. When comparing the overlap percentage between the

251

CMF + GAC method and the manually delineated results it achieved a result

252

of 91% compared to 86% achieved by the ACWE method, 88% achieved by

253

the GAC method and 87% achieved by the local-phase method.

254

Qualitative results on the medical data are presented in Figures 4 and 5.

255

Figure 4 shows the final result obtained with each method, while Figure 5

256

shows the similarity of each method to the reference delineation.

257

Starting with Row 1 in Figures 4 and 5, with respect to the GAC method,

258

it can be seen (for example in Figure 4 (R1,C2)), that the absence of a sharp

259

contrast in the upper right corner causes the method to over segment the

260

chamber in this region. The local-phase based method performs compara-

261

tively well, however over segments the chamber in the upper left region. The

262

ACWE method produces a better result, however slightly under segments the

263

chamber in this same region. The proposed CMF + GAC method achieves

264

a result visually quite similar to the manual delineation.

265

In Row 2 of Figures (4,5), the ACWE obtains the most accurate result

266

(see Figure 5 (R2,C2)). The absence of a sharp boundary disrupts the GAC

267

method while the local-phase method over segments the chamber. The pro-

268

posed CMF + GAC method obtains a mostly accurate result, however under

14

269

segments the chamber in the lower left region.

270

Row 3 presents an example of an abnormally shaped chamber as well as

271

limited window width in the upper right region. As in the previous examples,

272

the GAC and local-phase based methods tend to over-segment (see Figures

273

5 (R3,C2-C3) while the ACWE method, together with the CMF + GAC

274

approach capture the true shape more accurately. The CMF + GAC method

275

obtains a greater degree of similarity in the upper region of the chamber.

276

The frame in Row 4 presents an example of the lack of intensity homo-

277

geneity often found in echocardiographic data. The ACWE method (Figure

278

4 (R4,C4) and Figure 5 (R4,C4)) finds the cleanest part of the chamber,

279

however, as highlighted by the manual delineation, this does not represent

280

the true boundary. The results obtained by both the local-phase as well as

281

the CMF + GAC methods are visually quite similar, the phase based method

282

achieves a smoother boundary. The GAC method slightly under-segments

283

the upper chamber region.

284

Figure 1 presents the results of testing the CMF + GAC method, together

285

with the ACWE and GAC method on an image of a galaxy. The results show

286

that the GAC method finds the strongest boundary therefore missing the

287

lower gradient regions, while the ACWE method captures the bright star-

288

like points dotted around the image, which are not part of the galaxy, it can

289

be observed that the CMF + GAC results achieves the best approximation

290

of the boundary.

291

Figure 2 gives an example of birds swarming. In this example, the final

292

results of both the GAC and the CMF + GAC methods, are quite simi-

293

lar, both giving the approximate boundary of the swarm; while the ACWE

15

294

method partitions the image into visually the most intensity homogeneous

295

regions, consequently missing, in many areas, the actual boundary. Both the

296

GAC and CMF + GAC methods achieve comparable results, with a smoother

297

boundary than ACWE.

298

Table 2 gives the average for all datasets of the percentage error due to

299

change in the parameters used for each method. The CMF + GAC method

300

shows robust performance with respect to the length parameter, however

301

changing the initial estimates for inner and outer region intensity produces

302

a significant fluctuation in performance, this is mirrored, though to a lesser

303

degree, in the percentage error due to the image or contrast parameter in the

304

ACWE methods. Among the details which are highlighted in this table is

305

the lack of stability of the GAC method with this data; shifting the length

306

parameter by 5% for example, produces an error of almost 44% in the final

307

result.

308

4. Discussion

309

It can be observed from the results of the experiments on non-medical

310

data, that the CMF + GAC approach achieves either similar performance

311

to both the ACWE and GAC methods or enhanced performance in certain

312

cases; for example due to the lack of sharp contrast as illustrated in Figure

313

1, causing the GAC approach to fail, or when the predominant intensity

314

densities do not represent the actual boundary as in Figure 2, which disrupts

315

the ACWE method. The CMF + GAC approach achieves enhanced results

316

compared to either method in these instances. The qualitative results for the

317

medical data in Figures 4 to 5 largely mirror this observation. 16

318

With respect to the quantitative results in Table 1, for each metric the

319

proposed approach achieved the best overall result compared to each refer-

320

ence method; achieving an average Hausdorff metric of 4.48mm, an average

321

Mean Absolute Distance of 1.31mm and an average Dice coefficient of 0.91.

322

Comparing these results to the next most accurate method, which was the

323

GAC method, we see that for the Hausdorff metric, this represented a 5%

324

improvement on the GAC method, a 40% improvement as measured by Mean

325

Absolute Distance, while for the Dice coefficient, an improvement of 3% was

326

gained. With respect to individual datasets, we can observe that in 4 out of

327

5 cases the proposed approach achieved the lowest Hausdorff Distance, for

328

the Mean Absolute Distance this was the case in 3 out of 5 cases, while again

329

in 4 out 5 cases the proposed method achieved the best, or joint best result

330

with respect to the Dice coefficient.

331

The results also highlight some limitations of the model, such as in Figure

332

5 (R2,C5). In this example, due to the extremely low contrast boundary in

333

the left region of the chamber, the length parameter of the CMF method

334

was strongly relaxed, which has the positive effect of reinforcing the bound-

335

ary but also preserving noise (as exhibited in Figures 3(e),(f)). This effect

336

disrupted the CMF + GAC approach in this instance. This result is also

337

connected to the fluctuation in performance highlighted in sensitivity results

338

in Table 2. Due to the lack of clear distinction between the intensity char-

339

acteristics of the regions (i.e. the inner chamber and the inner wall of the

340

myocardium), choosing data fidelity parameters which are robust with re-

341

spect to variation is a challenging task. The sensitivity results in Table 2 for

342

each of the CMF + GAC, ACWE and GAC approaches reflect this finding.

17

343

These results highlight the limitations of data driven approaches in general,

344

in handling particularly low quality data. Notwithstanding this observed

345

fluctuation, however; and considering the results as a whole, it can be ob-

346

served that the CMF + GAC approach gives the most robust response with

347

respect to accuracy. Overall the sensitivity results indicate that the local-

348

phase approach is the most robust with respect to parameter variance, while

349

this approach achieved the least accuracy in terms of similarity to manual

350

delineation.

351

When considering the quantitative results in Table 1, it can be seen that

352

using the level-set ACWE method on the whole, while less accurate than

353

the CMF + GAC approach, is capable of achieving reasonable results on

354

low quality data. With this in mind, it could be asked what advantage does

355

solving the Mumford-Shah functional using the CMF approach have over a

356

traditional level set formulation?

357

In the traditional ACWE approach, tuning all algorithmic parameters

358

can be a challenging and rather unintuitive process, for example it needs to

359

be decided how often and accurately the level set function should be initial-

360

ized, how the time step sizes should be chosen etc. Choosing the correct

361

parameters for the CMF approach essentially means selecting the regulariza-

362

tion parameter ν and the penalty parameter c in the Augmented Lagrangian

363

method. The algorithm converges reliably for a wide range of c, (in all of

364

our experiments, this value remained unchanged). In addition, the CMF ap-

365

proach shows good robustness with respect to parameter ν (compare rows 2

366

to 3) under CMF + GAC approach, to rows 2-3 of ACWE method, in Table

367

2. The CMF approach is guaranteed to converge to a global minimizer, in

18

368

contrast to the level set approach which may get stuck in an inferior local

369

minimum.

370

The combination of the CMF approach with the GAC method also has

371

some interesting properties with respect to topology. Although the method

372

contains no explicit topology preserving constraint, by relaxing the length

373

constraint for the CMF approach, it becomes possible to favor the generation

374

of closed regions, this is shown for example in Figure 3(e). Figure 4 (R1,C3),

375

however, shows an example where this is impossible. In this example, the top

376

left portion of the chamber is outside of the viewing window, as such there

377

is not enough information for the CMF method to ‘reinforce’ the boundary,

378

however, by increasing the stiffness properties of the GAC contour, it is still

379

possible to maintain a closed topology. The reason for this is that the GAC

380

method converges to a local minimizer with simple topology, which in this

381

application is a better solution than the global minimizer. In this way the

382

combination of local and global segmentation methods works to preserve

383

topology.

384

While it has been shown that the proposed method is capable of accu-

385

rately segmenting low contrast images; it should be noted that the described

386

method to generate edge sets (step 1) is not proposed as a general edge detec-

387

tor approach, the detection of extremely fine structures, for example, would

388

prove problematic for the model. Instead the method has been designed for

389

the partition of structures with weak boundaries, as exemplified in a vari-

390

ety of real-world applications, a sample of which we have highlighted in the

391

current paper.

392

Finally, when observing the results it can be seen that global features

19

393

play a more significant role in the results achieved by the CMF + GAC

394

method, than local edge features, this is a key aspect of the method. While

395

the overall objective of the scheme is accurate boundary estimation, this goal

396

is achieved by integrating regional information into the edge detector process

397

using the CMF method, it is this regional information which makes boundary

398

reinforcement possible.

399

5. Conclusion

400

In this work we have presented a method to generate edge detectors which

401

allow for the effective segmentation of regions bounded by weak edges and

402

have demonstrated its performance on a sample of real-world medical and

403

non-medical data. It was shown that compared to classical image driven

404

approaches, the proposed method is capable of more accurately detecting

405

boundaries from low contrast images.

406

The method we have proposed increases robustness in two ways - first by

407

incorporating regional information into the process of edge map generation,

408

weak boundaries are effectively reinforced, secondly the edge map output

409

which corresponds to a 2 region global segmentation result, gives a more

410

robust solution compared to traditional purely variational based schemes.

411

Overall, we believe that both the qualitative and quantitative results

412

clearly demonstrate the improved performance of the CMF + GAC approach

413

over purely regional or intensity based methods as exemplified in the GAC

414

and ACWE methods.

20

(a)

(b)

(c)

(d)

(e)

Figure 1: Sombrero galaxy image [8]. The blue glow corresponds to the combined light of billions of old stars and the orange ring marks the sites of young stars: (a) original image, (b) initial phi,(c) Geodesic Result,(d) ACWE result, (e) CMF + GAC result. The results show that the GAC method finds the strongest boundary therefore missing the lower gradient regions (which also form part of the galaxy), while both the ACWE and CMF + GAC methods achieve accurate results. Parameters used for the CMF + GAC method: c1 = 0.18; c2 = 0.09; ν = 5e − 1, expansion term = 1, length parameter = 10. Courtesy NASA/JPL-Caltech

21

(a)

(b)

(c)

(d)

(e)

Figure 2: Birds; (a) original, (b) initial phi, (c) Geodesic Result, (d) ACWE result, (e) CMF + GAC result. In this example, both the GAC and the CMF + GAC methods achieve good results while the ACWE method partitions the image into visually the most intensity homogeneous regions missing, in many areas, the actual boundary. Parameters CMF + GAC method: initial c1 = 0.02; c2 = 0.7, ν = 0.01 ; expansion term = 0.6, length parameter = 10; Photograph by Lucas Felzmann. Used with permission.

22

Table 1: Quantitative results computed using Hausdorff, Mean Absolute Distance and Dice metrics for medical data. Experiments were conducted on a sample of 55 frames across 5 datasets. Represented here is the average distance for each metric computed over a test set of frames from each dataset

HD LAX1 LAX2 LAX3 LAX4 LAX6 Average % diff Our proposed Method 4.6

5.28

3.34 4.69 4.49

4.48

ACWE 6.78

9.91

3.37

5.91

5.47

6.29

40.4

Standard geodesic 5.31

3.8

4.08

5.64

4.66

4.7

4.9

Using Phase derived edge maps 6.62

5.48

5.02

8.22

8.52

6.77

51.1

MAD LAX1 LAX2 LAX3 LAX4 LAX6 Average % diff Our proposed Method 1.17

1.51

0.94

ACWE 1.82

3.46

0.89

1.9

1.45

1.9

45

Standard geodesic 2.06

1.39

1.81

1.85

2.11

1.84

40.5

1.6

1.98

2.26

1.95

1.9

45

Using Phase derived edge maps 1.72

1.54 1.41

1.31

Dice LAX1 LAX2 LAX3 LAX4 LAX6 Average % diff Our proposed Method 0.93

0.92

0.93

0.88

0.9

0.91

ACWE 0.87

0.75

0.93

0.84

0.9

0.86

5

Standard geodesic 0.89

0.93

0.87

0.87

0.86

0.88

3

0.9

0.85

0.83

0.87

0.87

4

Using Phase derived edge maps

0.9

23

(a) Original Image(b)

Output

of

(c) Edge map

of

(f) Edge map

CMF method

(d) Original Image(e)

Output

CMF method

(g) CMF

Output

of(h) Corresponding

method,Edge map

increased ν

Figure 3: Example output of CMF method for 2 datasets; (a) and (d) original frames; (b), (e) and (g) corresponding output of CMF method, using values of ν = 0.5, ν = 1e − 02 and ν = 0.5 respectively; (c), (f) and (h) edge map generated from binary segmentation. Figures (e) and (g) show the difference produced by varying ν,(g) shows a cleaner result however (e) more accurately represents the underlying (closed) structure. The output of the CMF method is only shown to illustrate each step of the proposed method and is not part of the final output. To save space, the CMF result is only shown for a subset of the experiments.

24

R1

R2

R3

R4 C1

C2

C3

C4

C5

Figure 4: Results achieved by each method on four sample frames (Rows R1 R4). Each frame shows an example of boundaries characterized by weak edges. C1: The original frame. C2: Standard GAC method. C3: Local-phase based method. C4: ACWE method. C5: Results obtained using the proposed method.

25

R1

R2

R3

R4 C1

C2

C3

C4

C5

Figure 5: Results showing the similarity between the proposed method and the expert manual delineations. C1 − C5: from left to right, the GAC, local-phase method, ACWE and the CMF+GAC method combined with the manual delineation. It can be observed that the standard GAC method, as well as the Local-Phase Based method tend to oversegment the chamber. C4: shows the results of the ACWE method, which on average gives a good result, however it can be seen in R1, R3 and R4; that the proposed method achieves the greater degree of similarity to the reference contour.

26

Table 2: Percentage error due to each parameter. Average taken over all frames Proposed Method

HD Average Std vary length parameter by +/- 5% 0.92 1.12 vary length parameter by +/- 10% 1.74 1.94 vary length parameter by +/- 20% 6.33 11.57 vary c1 or c2 by +/- 5% 5.59 9.47 vary c1 or c2 by +/- 10% 12.12 22.97 vary c1 or c2 by +/- 20% 28.71 44.09 Expansion Parameter (tuned for 1 dataset) vary expansion parameter by +/-5% 6.39 6.39 vary expansion parameter by +/-10% 14.5 1.72 vary expansion parameter by +/-20% 74.69 60.21 ACWE vary length parameter by +/-5% vary length parameter by +/-10% vary length parameter by +/-20% vary image parameter by +/-5% vary image parameter by +/-10% vary image parameter by +/-20% GAC

MAD Average 0.16 0.38 2.68 3.36 6.41 13.64

Dice Std Average Std 0.15 0.01 0.01 0.37 0.03 0.03 5.83 0.26 0.61 4.48 0.93 1.35 8.99 1.34 1.75 16.51 2.37 2.48

2.99 1.51 12.13 10.65 28.28 25.74

0.41 1.32 5.99

0 0.91 5.5

HD Average 4.19 11.71 11.03 15.96 20.15 31.57

MAD Dice Std Average Std Average Std 9.38 3.12 6.02 1.12 2.96 19.02 5.55 6.47 1.52 2.3 8.48 6.12 4.76 1.59 2.23 20.53 7.09 8.64 1.95 3.18 28.63 7.62 8.55 2.06 2.86 33.87 16.6 15.49 3.86 4.08

HD Average 25.64 48.92 102.89 43.78 44.48 47.74

Std 36.9 35.57 56.99 61.76 61.26 60.45

MAD Average 21.09 47.1 87.07 34.42 35.53 38.8

Dice Std Average Std 30.5 4.68 7.69 28.44 10.37 7.86 50.38 13.98 10.32 44.15 6.49 7.21 43.12 5.54 7.07 40.23 6.39 6.58

HD MAD Average Std Average vary wavelength parameter by +/-5% 6.63 5.06 2.55 vary wavelength parameter by +/-10% 7.5 5.25 3.9 vary wavelength parameter by +/-20% 10.31 5.28 5.95 Expansion Parameter (tuned for 3 datasets) vary expansion parameter by +/-5% 0 0 0.01 vary expansion parameter by +/-10% 2.18 2.36 1.24 vary expansion parameter by +/-20% 11.74 10.81 6.85

Dice Std Average Std 1.75 0.41 0.37 2.51 0.68 0.55 3.03 0.91 0.59

vary expansion parameter by +/-5% vary expansion parameter by +/-10% vary expansion parameter by +/-20% vary contour parameter by +/-5% vary contour parameter by +/-10% vary contour parameter by +/-20% Phase

27

0.01 1.12 4.42

0 0.23 1.27

0 0.2 0.87

415

Acknowledgements

416

The authors would especially like to thank the authors of [19] for pro-

417

viding us with the implementation of the Local Phase method as well as

418

the anonymous reviewers for their valuable comments and suggestions. This

419

work was supported by the National Science Foundation grant NSF DMS-

420

1217239, UC Lab Fees Research grant 12-LR-23660, Norwegian Research

421

Council eVita project 214889 and the Irish Research Council.

422

References

423

[1] Bresson, X., Esedoglu, S., Vandergheynst, P., Thiran, J.P., Osher, S.,

424

2007. Fast global minimization of the active contour/snake model. J

425

Math Imaging Vis 28, 151–167.

426

427

428

429

[2] Bui, T.D., Gao, S., Zhang, Q., 2005. A generalized mumford-shah model for roof-edge detection, in: ICIP (2), pp. 1214–1217. [3] Caselles, V., Kimmel, R., Sapiro, G., 1997. Geodesic active contours. Int J Comput Vision 22, 61–79.

430

[4] Chan, T., Esedoglu, S., Nikolova, M., 2006. Algorithms for finding global

431

minimizers of image segmentation and denoising models. Siam J Appl

432

Math , 17 pp.

433

434

[5] Chan, T., Vese, L., 2001. Active contours without edges. IEEE T Image Process 10, 266–277.

435

[6] Duggan, N., Schaeffer, H., Guyader, C.L., Jones, E., Glavin, M., Vese,

436

L., 2013. Boundary detection in echocardiography using a split bregman 28

437

edge detector and a topology-preserving level set approach, in: Interna-

438

tional Symposium on Biomedical Imaging (ISBI), IEEE.

439

440

[7] Felsberg, M., Sommer, G., 2001. The monogenic signal. IEEE Transactions on Signal Processing 49, 3136–3144.

441

[8] Ford, H.C., Hui, X., Ciardullo, R., Jacoby, G.H., Freeman, K.C., 1996.

442

The stellar halo of M104. I. A survey for planetary nebulae and the

443

planetary nebula luminosity function distance. Astrophysical Journal

444

458, 455–466. doi:doi:10.1086/176828.

445

[9] Goldstein, T., Bresson, X., Osher, S., 2010. Geometric applications of

446

the split bregman method: Segmentation and surface reconstruction. J.

447

Sci. Comput. 45, 272–293.

448

[10] Hyman, J.M., Shashkov, M.J., 1997. Natural discretizations for the

449

divergence, gradient, and curl on logically rectangular grids. Comput.

450

Math. Appl. 33, 81–104.

451

452

[11] Kass, M., Witkin, A., Terzopoulos, D., 1987. Snakes: Active contour models. Int J Comput Vision , 321 –331.

453

[12] Kichenassamy, S., Kumar, A., Olver, P., Tannenbaum, A., Yezzi, A.,

454

1996. Conformal curvature flows: from phase transitions to active vision.

455

Arch. Rat. Mech. Anal. 134, 275–301.

456

457

[13] Mulet-Parada, M., Noble, J., 2000. 2d+t acoustic boundary detection in echocardiography. Medical Image Analyis 4, 21–30.

29

458

[14] Mumford,

D.,

Shah,

J.,

1989.

Optimal approximations by

459

piecewise smooth functions and associated variational problems.

460

Communications on Pure and Applied Mathematics 42, 577–685.

461

doi:10.1002/cpa.3160420503.

462

[15] Osher, S., Sethian, J.A., 1988. Fronts propagating with curvature-

463

dependent speed: algorithms based on hamilton-jacobi formulations. J

464

Comput Phys 79, 12–49.

465

[16] Paragios, N., Deriche, R., 1999. Unifying boundary and region-based

466

information for geodesic active tracking, in: IEEE Conference on Com-

467

puter Vision and Pattern Recognition.

468

[17] Perona, P., Malik, J., 1990.

Scale-space and edge detection using

469

anisotropic diffusion. IEEE Trans. Pattern Anal. Mach. Intell. 12, 629–

470

639.

471

[18] Pock, T., Cremers, D., Bischof, H., Chambolle, A., 2009. An algorithm

472

for minimizing the mumford-shah functional, in: ICCV, pp. 1133–1140.

473

[19] Rajpoot, K., Grau, V., Noble, J., 2009a. Local-phase based 3d bound-

474

ary detection using monogenic signal and its application to real-time

475

3-d echocardiography images, in: IEEE International Symposium on

476

Biomedical Imaging, pp. 783 –786.

477

[20] Rajpoot, K., Noble, J.A., Grau, V., Szmigielski, C., Becher, H., 2009b.

478

Image-driven cardiac left ventricle segmentation for the evaluation of

479

multiview fused real-time 3-dimensional echocardiography images, in:

480

MICCAI, pp. 893–900. 30

481

[21] Schaeffer, H., 2012. Active arcs and contours. UCLA CAM Report 12-

482

54. URL: ftp://ftp.math.ucla.edu/pub/camreport/cam12-54.pdf.

483

[22] Schaeffer, H., Vese, L., 2013. Active contours with free endpoints. Jour-

484

nal of Mathematical Imaging and Vision , 1–17.

485

[23] Vese, L., Chan, T., 2002. A multiphase level set framework for image

486

segmentation using the mumford and shah model. Int J Comput Vision

487

50, 271–293.

488

[24] Wang, L.L., Shi, Y., Tai, X.C., 2012. Robust edge detection using

489

mumford-shah model and binary level set method, in: Bruckstein, A.,

490

Haar Romeny, B., Bronstein, A., Bronstein, M. (Eds.), Scale Space and

491

Variational Methods in Computer Vision. volume 6667 of LNCS, pp.

492

291–301.

493

[25] Weickert, J., Romeny, B.M.T.H., Viergever, M.A., 1998. Efficient and

494

reliable schemes for nonlinear diffusion filtering. IEEE Transactions on

495

Image Processing 7, 398–410.

496

[26] Yuan, J., Bae, E., Tai, X.C., 2010. A study on continuous max-flow

497

and min-cut approaches, in: Computer Vision and Pattern Recognition

498

(CVPR), 2010 IEEE Conference on, pp. 2217–2224.

499

[27] Yuan, J., Bae, E., Tai, X.C., Boykov, Y., 2014. A spatially continuous

500

max-flow and min-cut framework for binary labeling problems. Nu-

501

merische Mathematik 126, 559–587.

31

Highlights: An method for segmentation of low contrast images is proposed The global optimal output of the CV model is used to create the edge attraction force in the GAC method. Experimental results demonstrate improved performance over existing data driven methods