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
041.         self.work = self._pad()
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. 
062.     def _pad(self):
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. 
200.     def addOffset(self, p, index):
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