ipdp - classes and functions to facilitate the coding of watershed algorithms

```001. from numpy import array, ravel
002.
003. # constants
004. # neighbourhood 2D
005. N4 = array([[-1,0],[0,1],[1,0],[0,-1]])
006. N8 = array([[-1,0],[0,1],[1,0],[0,-1],[-1,-1],[-1,1],[1,1],[1,-1]])
007.
008. # neighbourhood 3D
009. N6 = array([[-1,0,0],[0,1,0],[1,0,0],[0,-1,0],[0,0,1],[0,0,-1]])
010.
011.
012. # values
013. inf = 1e400
014. BORDER = inf
015.
016. def isBorder(v):
017.     """ Tests if the value is a border value """
018.     return v == BORDER
019.
020. # common classes for the algorithms
021. class wsImage():
022.     """ Class for storing the images """
023.
024.     def __init__(self, array):
025.         """ Constructor for storing the N-D input image """
026.         self.input = array
027.         self.work = None
028.         self.label = None
029.         self.output = None
030.
031.     def begin(self, offsets):
032.         """ Prepare the image for processing """
033.         from numpy import zeros, ravel
034.
035.         if len(self.input.shape) != offsets.shape[1]:
036.             raise Exception("Image shape does not fit offsets dimensions")
037.
038.         # initialise the neighbour
039.         self.neighbour = wsNeighbour(offsets)
040.         # pad the input image
042.         # store the padded shape
043.         self.workshape = self.work.shape
044.         # ravel the padded image (make 1-D)
045.         self.work = ravel(self.work)
046.         # make a zeroed copy of it
047.         self.label = zeros(self.work.shape)
048.         # initialise the output
049.         self.output = None
050.         # initialise the shape of the image
051.         self.neighbour.setImageShape(self.workshape)
052.         self.neighbour.setImage(self.work)
053.         # initialise the domain object
054.         D = wsDomain(self.work.size)
055.         D.setImage(self.work)
056.
057.         # returns the neighbourhood relation, the working image, the label
058.         # image and the domain of the image
059.         return self.neighbour.N, self.work, self.label, D
060.
061.
063.         from numpy import zeros
064.         """ Pads the N-D image with the BORDER constant as necessary for
065.         containing all the offsets from the list
066.         """
067.         # generate the newshape by iterating through the original and
068.         # adding the extra space needed at each dimension
069.         newshape = tuple(map(lambda orig, d: orig + (d-1), self.input.shape, self.neighbour.shape))
070.         # generate the slicing list
071.         slicing = map(lambda orig, d: slice((d-1)/2, (d-1)/2 + orig), self.input.shape, self.neighbour.shape)
072.         # create the padded image
073.         workimage = zeros(newshape)
074.         workimage[:] = BORDER
075.         workimage[slicing] = self.input
076.         return workimage
077.
078.     def _crop(self):
079.         return self.crop(self.label)
080.
081.     def crop(self, x):
082.         """ Reshape and crops the N-D label image to the original size (same as self.input) """
083.         from numpy import reshape
084.         # generate the slicing list
085.         slicing = map(lambda orig, d: slice((d-1)/2, (d-1)/2 + orig), self.input.shape, self.neighbour.shape)
086.         # reshape the label image to the original shape
087.         temp = reshape(x, self.workshape)
088.         # crop the temp image
089.         return temp[slicing]
090.
091.     def end(self):
092.         return self._crop().astype('int32')
093.
094.     def makeWorkCopy(self, default=0):
095.         """ Make a copy of the work image filled with the value and type of parameter default """
096.         copied = self.work.copy()
097.         copied = copied.astype(type(default))
098.         copied.fill(default)
099.         return copied
100.
101.     def getCoordinate(self, p):
102.         """ Get the coordinates in (x,y) form for index p"""
103.         return (p / self.workshape[1], p % self.workshape[1])
104.
105.
106. class wsNeighbour():
107.     """ Class for neighbourhood processing """
108.
109.     def __init__(self, offsets):
110.         """ Constructor for the list of offsets in N-D (neighbours)
111.             offsets must be a m x N matrix, where m is the number of
112.             offsets (neighbours) and N is the dimensions of the image"""
113.         self.offsets = array(offsets)
114.         self._shape()
115.         self.s = None
116.
117.     def _shape(self):
118.         """ Calculates the shape of the offsets """
119.         N = self.offsets.shape[1]
120.         self.shape = []
121.         for i in range(N):
122.             dmax = max(self.offsets[:,i])
123.             dmin = min(self.offsets[:,i])
124.             if abs(dmax) != abs(dmin):
125.                 raise Exception("Offsets must be symmetrical")
126.             d = dmax - dmin + 1
127.             # make the dimension always odd
128.             if d % 2 == 0:
129.                 d += 1
130.             self.shape.append(d)
131.         self.shape = tuple(self.shape)
132.
133.     def setImageShape(self, imshape):
134.         """ Set the image shape and calculates the offsets in 1-D """
135.         self.s = imshape
136.         if len(self.s) != self.offsets.shape[1]:
137.             raise Exception("Image shape does not fit offsets dimensions")
138.
139.         # calculate the offsets in 1-D
140.         # the process occurs like this:
141.         # each offset is multiplied by the multiplication of the values of
142.         # the next components of the shape of the image and summed:
143.         # example:
144.         # shape: (3, 10, 15)
145.         # offset: [1, 1, 2]
146.         # offset in 1-D: (1 * 10 * 15) + (1 * 15) + (2)
147.         #
148.         # of course, the offsets must follow the order of the shape
149.         # (Nth-D, ..., 3rd-D, 2nd-D, 1st-D), that is usually
150.         # (time, channel, row, column) or in grayscale images (time, row, column)
151.         # or simple 2-D images (row, column)
152.
153.         # LONG VERSION
154.         # self.roffsets = []
155.         # for offset in self.offsets:
156.             # roffset = 0
157.             # for i in range(len(offset)):
158.                 # n = offset[i]
159.                 # roffset += n * reduce(lambda x,y: x*y, self.s[(i+1):], 1)
160.             # self.roffsets.append(roffset)
161.
162.         # SHORT VERSION (using map and reduce - I like this one better)
163.         self.roffsets = map(
164.             lambda offset: sum(
165.                 map(lambda n, i:
166.                         n * reduce(lambda x,y: x*y, self.s[(i+1):], 1),
167.                     offset, range(len(offset))
168.                     )
169.                 ),
170.             self.offsets)
171.
172.     def setImage(self, im):
173.         """ Set the working image to query for border values on neighbourhood calculation """
174.         self.im = im
175.
176.     def N(self, pixel):
177.         """ Returns the list of indexes of neighbours of pixel in 1-D """
178.         if not self.s:
179.             raise Exception("Set the image shape first!")
180.
181.         # calculate the coordinates of the neighbours based on the offsets in 1-D
182.         n = map(lambda c: c + pixel, self.roffsets)
183.         r = list()
184.
185.         for i in n:
186.             if isBorder(self.im[i]):
187.                 continue
188.
189.             r.append(i)
190.
191.         return r
192.
193.     def query(self, c):
194.         """ Lookup on the roffsets for the index of the offset c """
195.         for i in range(len(self.roffsets)):
196.             if self.roffsets[i] == c:
197.                 return i
198.         return -1
199.
201.         """ Adds the offset of the desired index to the value p """
202.         if index < 0 or index >= len(self.roffsets):
203.             return None
204.         else:
205.             c = self.roffsets[index]
206.             return p + c
207.
208. class wsDomain:
209.     def __init__(self, length):
210.         self.inner = xrange(length)
211.         self.count = 0
212.
213.     def next(self):
214.         if self.count >= len(self.inner):
215.             self.count = 0
216.             raise StopIteration
217.
218.         while isBorder(self.im[self.inner[self.count]]):
219.             self.count += 1
220.             if self.count >= len(self.inner):
221.                 self.count = 0
222.                 raise StopIteration
223.
224.         c = self.count
225.         self.count = c + 1
226.         return self.inner[c]
227.
228.     def __getitem__(self, item):
229.         return self.inner[item]
230.
231.     def __iter__(self):
232.         return self
233.
234.     def __len__(self):
235.         return len(self.inner)
236.
237.     def setImage(self, im):
238.         self.im = im
239.
240. class wsQueue():
241.     """ Simple queue class for abstracting list methods """
242.
243.     def __init__(self):
244.         self.array = list()
245.
246.     def push(self, a):
247.         """ Pushes an element to the end of the queue """
248.         self.array.append(a)
249.
250.     def pop(self):
251.         """ Pops the first element of the queue """
252.         return self.array.pop(0)
253.
254.     def empty(self):
255.         """ Returns whether the queue is empty or not """
256.         return len(self.array) == 0
257.
258.     def clear(self):
259.         """ Clear the list """
260.         self.array = list()
261.
262. class wsStack():
263.     """ Simple stack class for abstracting list methods """
264.
265.     def __init__(self):
266.         self.array = list()
267.
268.     def push(self, a):
269.         """ Pushes an element to the top of the stack """
270.         self.array.append(a)
271.
272.     def pop(self):
273.         """ Pops the top element of the stack """
274.         return self.array.pop()
275.
276.     def empty(self):
277.         """ Returns whether the stack is empty or not """
278.         return len(self.array) == 0
279.
280.     def clear(self):
281.         """ Clear the list """
282.         self.array = list()
283.
284. class wsHeapQueue():
285.     """ Priority queue class for abstracting list methods with FIFO policy """
286.
287.     def __init__(self):
288.         self.queue = dict()
289.
290.     def push(self, a, c):
291.         """ Pushes an element to the queue """
292.         if self.queue.has_key(c):
293.             self.queue[c].append(a)
294.         else:
295.             self.queue[c] = [a]
296.
297.     def pop(self):
298.         """ Pops the first element of the queue """
299.         key = min(self.queue.keys())
300.         element = self.queue[key].pop(0)
301.         if len(self.queue[key]) == 0:
302.             self.queue.pop(key)
303.         return element
304.
305.     def empty(self):
306.         """ Returns whether the queue is empty or not """
307.         return len(self.queue) == 0
308.
309.     def clear(self):
310.         """ Clear the queue """
311.         self.queue = dict()
312.
313.     def remove(self, a, c):
314.         """ Remove the element a at cost c """
315.         self.queue[c].remove(a)
316.
317.     def contains(self, a, c):
318.         """ Verifies if the queue contains element a at cost c """
319.         for x in self.queue[c]:
320.             if x == a:
321.                 return True
322.         return False
323.
324. class wsRandHeapQueue():
325.     """ Priority queue class for abstracting list methods without FIFO policy """
326.
327.     class wsPixel():
328.
329.         def __init__(self, p, l):
330.             self.p = p
331.             self.l = l
332.
333.         def __lt__(self, other):
334.             return self.l < other.l
335.
336.         def __le__(self, other):
337.             return self.l <= other.l
338.
339.         def __eq__(self, other):
340.             return self.l == other.l and self.p == other.p
341.
342.         def __ne__(self, other):
343.             return not (self.l == other.l and self.p == other.p)
344.
345.         def __gt__(self, other):
346.             return self.l > other.l
347.
348.         def __ge__(self, other):
349.             return self.l >= other.l
350.
351.     def __init__(self):
352.         self.queue = []
353.
354.     def push(self, a, c):
355.         from heapq import heappush
356.         """ Pushes an element to the queue """
357.         px = wsRandHeapQueue.wsPixel(a, c)
358.         heappush(self.queue, px)
359.
360.     def pop(self):
361.         from heapq import heappop
362.         """ Pops the first element of the queue """
363.         px = heappop(self.queue)
364.         return px.p
365.
366.     def empty(self):
367.         """ Returns whether the queue is empty or not """
368.         return len(self.queue) == 0
369.
370.     def clear(self):
371.         """ Clear the queue """
372.         self.queue = []
373.
374.     def remove(self, a, c):
375.         from heapq import heapify
376.         """ Remove the element a at cost c """
377.         self.queue.remove(wsRandHeapQueue.wsPixel(a,c))
378.         heapify(self.queue)
379.
380.     def contains(self, a, c):
381.         """ Verifies if the queue contains element a at cost c """
382.         for px in self.queue:
383.             if px.p == a and px.l == c:
384.                 return True
385.         return False
386.
387.
388. def findMinima(im, N, D):
389.
390.     UNVISITED = 0
391.     PENDING = 1
392.     VISITED = 2
393.
394.     work = im.copy()
395.     work[:] = UNVISITED
396.
397.     qPending = wsQueue()
398.     qMinima = wsQueue()
399.
400.     minima = list()
401.
402.     count = 1
403.
404.     for p in D:
405.         if work[p] != UNVISITED:
406.             continue
407.
408.         qPending.push(p)
409.         work[p] = PENDING
410.
411.         is_minima = True
412.         qMinima.clear()
413.
414.         while not qPending.empty():
415.
416.             q = qPending.pop()
417.             qMinima.push(q)
418.             work[q] = VISITED
419.
420.             for u in N(q):
421.
422.                 if im[u] == im[q] and work[u] == UNVISITED:
423.                     work[u] = PENDING
424.                     qPending.push(u)
425.                 elif im[u] < im[q]:
426.                     is_minima = False # keep flooding the plateau
427.
428.
429.
430.         if is_minima:
431.             count += 1
432.             seed = list()
433.             while not qMinima.empty():
434.                 q = qMinima.pop()
435.                 seed.append(q)
436.             minima.append(seed)
437.
438.     return minima
439.
440. def lowerComplete(im, offsets):
441.
442.     # initialise variables
443.     ws = wsImage(im)
444.     N, im, lc, D = ws.begin(offsets)
445.
446.     FICTITIOUS_PIXEL = -1
447.     queue = wsQueue()
448.
449.     lc[:] = BORDER
450.
451.     for p in D:
452.
453.         lc[p] = 0
454.         for q in N(p):
455.             if im[q] < im[p]:
456.                 queue.push(p)
457.                 lc[p] = -1
458.                 break
459.
460.     cur_dist = 1
461.
462.     queue.push(FICTITIOUS_PIXEL)
463.
464.     while not queue.empty():
465.         p = queue.pop()
466.         if p == FICTITIOUS_PIXEL:
467.             if not queue.empty():
468.                 queue.push(FICTITIOUS_PIXEL)
469.                 cur_dist += 1
470.         else:
471.             lc[p] = cur_dist
472.             for q in N(p):
473.
474.                 if im[q] == im[p] and lc[q] == 0:
475.                     queue.push(q)
476.                     lc[q] = -1
477.
478.     for p in D:
479.
480.         if lc[p] == 0:
481.             lc[p] = cur_dist * im[p]
482.         else:
483.             lc[p] = cur_dist * im[p] + lc[p] - 1
484.
485.     return ws.end()
486.
487. def se2offset(se):
488.     from numpy import array
489.     from ia870 import iaseshow
490.
491.     se = iaseshow(se, 'NORMAL')
492.
493.     offset = []
494.     for i in range(se.shape[0]):
495.         for j in range(se.shape[1]):
496.             if se[i,j] == True and (i-se.shape[0]/2 != 0 or j-se.shape[1]/2 != 0):
497.                 offset.append([i-se.shape[0]/2,j-se.shape[1]/2])
498.
499.     return array(offset)```
`1. pass`