berge

Synopsis

Prototype:
y = berge(f, offsets)
Description:
The berge function implements the watershed transform by the Berge's [Berge1962] algorithm for shortest path forests. This algorithm requires the input image f has no plateaus. This algorithm was proposed by Meyer [MeyerSP1994] and corrected by Roerdink and Meijster [RoerdinkMeijsterFI2000] to generate the watershed pixels.
Definition:
TD-WT
Exploration:
Any order

Algorithm

Algoritmo 1: Berge

References

[Berge1962]C. Berge, The Theory of graphs and its applications, 1st ed. New York: John Wiley & Sons, Inc., 1962.
[MeyerSP1994]F. Meyer, ``Topographic distance and watershed lines,'' Signal Processing, vol. 38, no. 1, pp. 113–125, 1994.
[RoerdinkMeijsterFI2000]J. B. T. M. Roerdink and A. Meijster, ``The watershed transform: definitions, algorithms and parallelization strategies,'' Fundam. Inf., vol. 41, no. 1-2, pp. 187–228, 2000.

Source Code

 1 from ipdp.common import *
 2 
 3 # constants
 4 MASK = -2
 5 WSHED = 0
 6 
 7 def berge(im, offsets):
 8 
 9     # initialise variables
10     lc = lowerComplete(im, offsets)
11     ws = wsImage(lc)
12     N, im, lab, D = ws.begin(offsets)
13 
14     # find minima
15     M = findMinima(im, N, D)
16 
17     def cost(p, q):
18         LSp = im[p] - min(im[N(p)])
19         LSq = im[q] - min(im[N(q)])
20         if im[p] > im[q]:
21             return LSp
22         elif im[q] > im[p]:
23             return LSq
24         else:
25             return (LSp + LSq)/2.0
26 
27     def propagate(p, Q, cmp):
28         s = True
29         for q in Q:
30             if not cmp(q, p):
31                 continue
32 
33             c = cost(p,q)
34             if dist[p] + c < dist[q]:
35                 dist[q] = dist[p] + c
36                 lab[q] = lab[p]
37                 s = False
38             elif dist[p] + c == dist[q] and lab[q] != lab[p] and lab[q] != WSHED:
39                 lab[q] = WSHED
40                 s = False
41 
42         return s
43 
44     lab[:] = MASK
45     dist = ws.makeWorkCopy(inf)
46     stable = False
47 
48     D = list(D)
49 
50     for m in xrange(len(M)):
51         for p in M[m]:
52             lab[p] = m+1
53             dist[p] = im[p]
54 
55     while not stable:
56         stable = True
57         for p in D:
58             stable = stable and propagate(p, N(p), lambda u,v: u > v)
59 
60         for p in reversed(D):
61             stable = stable and propagate(p, N(p), lambda u,v: u < v)
62 
63     return ws.end()