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