shortestPaths

Synopsis

Prototype:
y = shortestPaths(f, offsets)
Description:
The shortestPaths function implements the watershed transform by calculating shortest paths for each pixel. This algorithm attempts to visit every pixel the least times as possible. That is accomplished using queues for ordering processing. This algorithm was proposed by Osma-Ruiz et al. [OsmaRuizPR2007]
Definition:
LC-WT
Exploration:
Depth-First

Algorithm

Algoritmo 1: Shortest Paths

References

[OsmaRuizPR2007]V. Osma-Ruiz, J. I. Godino-Llorente, N. Sáenz-Lechón, and P. Gómez-Vilda, ``An improved watershed algorithm based on efficient computation of shortest paths,'' Pattern Recognition, vol. 40, no. 3, pp. 1078–1090, 2007.

Source Code

  1 from ipdp.common import *
  2 
  3 #constants
  4 UNVISITED = -32767
  5 PENDING = -32766
  6 
  7 def shortestPaths(im, offsets):
  8 
  9 
 10     # initialise variables
 11     ws = wsImage(im)
 12     N, im, lab, D = ws.begin(offsets)
 13     lab[:] = UNVISITED
 14     basins = 1
 15 
 16     qPending = wsQueue()
 17     qEdge = wsQueue()
 18     qInner = wsQueue()
 19     qDescending = wsQueue()
 20 
 21     def arrow(p,q):
 22         # finds the offset
 23         c = q - p
 24         # query the neighbourhood for the offset
 25         idx = ws.neighbour.query(c)
 26         if idx == -1:
 27             raise Exception("Neighbour not found")
 28         return - idx # inverted signal
 29 
 30     def pointed(p, direction):
 31         # add the offset of the previously detected direction
 32         direction = int(direction)
 33         q = ws.neighbour.addOffset(p, -direction) # inverted signal
 34         if not q:
 35             raise Exception("Invalid direction")
 36         return q
 37 
 38     # start looping
 39     for p in D:
 40 
 41         if lab[p] != UNVISITED:
 42             continue
 43 
 44         minvalue = min(im[N(p)])
 45         minpx = None
 46 
 47         for q in N(p):
 48 
 49             if im[q] == im[p]:
 50                 if qPending.empty():
 51                     lab[p] = PENDING
 52                     qPending.push(p)
 53                 lab[q] = PENDING
 54                 qPending.push(q)
 55             elif im[q] < im[p] and im[q] == minvalue:
 56                 minpx = q
 57 
 58         if not qPending.empty():
 59             while not qPending.empty():
 60                 q = qPending.pop()
 61                 if q != p:
 62 
 63                     minvalue = min(im[N(q)])
 64                     minpx = None
 65 
 66                     for u in N(q):
 67 
 68                         if im[u] == im[p]:
 69                             if lab[u] == UNVISITED:
 70                                 lab[u] = PENDING
 71                                 qPending.push(u)
 72                         elif im[u] < im[q] and im[u] == minvalue:
 73                             minpx = u
 74 
 75                 if not minpx is None:
 76                     lab[q] = arrow(q, minpx)
 77                     qEdge.push(q)
 78                 else:
 79                     qInner.push(q)
 80 
 81             if not qEdge.empty():
 82                 if not qInner.empty():
 83                     while not qEdge.empty():
 84                         q = qEdge.pop()
 85                         for u in N(q):
 86 
 87                             if im[u] == im[q] and lab[u] == PENDING:
 88                                 lab[u] = arrow(u,q)
 89                                 qEdge.push(u)
 90                 else:
 91                     qEdge.clear()
 92                 qInner.clear()
 93             else:
 94                 while not qInner.empty():
 95                     q = qInner.pop()
 96                     lab[q] = basins
 97                 basins += 1
 98         else:
 99             if minpx is None:
100                 lab[p] = basins
101                 basins += 1
102             else:
103                 lab[p] = arrow(p,minpx)
104 
105 
106     for p in D:
107 
108         q = p
109         while lab[q] <= 0:
110             qDescending.push(q)
111             q = pointed(q, lab[q])
112 
113         while not qDescending.empty():
114             u = qDescending.pop()
115             lab[u] = lab[q]
116 
117     # crop and return
118     return ws.end()