wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-07-17 11:32:45 +08:00
Code Issues Packages Projects Releases Wiki Activity

Initial integration of CellProfiler's watershed module

Browse Source
This commit is contained in:
emmanuelle
2011-10-05 21:13:46 +02:00
parent de38b5f93c
commit 656962cbc9
8 changed files with 1006 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files
scikits/image/filter/__init__.py
+1
View File
@@ -3,3 +3,4 @@ from ctmf import median_filter
from canny import canny
from edges import sobel, hsobel, vsobel, hprewitt, vprewitt, prewitt
from tv_denoise import tv_denoise
from rank_order import rank_order
scikits/image/morphology/__init__.py
+1
View File
@@ -1,3 +1,4 @@
from grey import *
from selem import *
from ccomp import label
from watershed import watershed, fast_watershed
scikits/image/morphology/_watershed.pyx
+99
View File
@@ -0,0 +1,99 @@
"""watershed.pyx - scithon implementation of guts of watershed
Originally part of CellProfiler, code licensed under both GPL and BSD licenses.
Website: http://www.cellprofiler.org
Copyright (c) 2003-2009 Massachusetts Institute of Technology
Copyright (c) 2009-2011 Broad Institute
All rights reserved.
Original author: Lee Kamentsky
"""
cdef extern from "numpy/arrayobject.h":
cdef void import_array()
import_array()
import numpy as np
cimport numpy as np
cimport cython
DTYPE_INT32 = np.int32
ctypedef np.int32_t DTYPE_INT32_t
DTYPE_BOOL = np.bool
ctypedef np.int8_t DTYPE_BOOL_t
include "heap_watershed.pxi"
@cython.boundscheck(False)
def watershed(np.ndarray[DTYPE_INT32_t,ndim=1,negative_indices=False, mode='c'] image,
np.ndarray[DTYPE_INT32_t,ndim=2,negative_indices=False, mode='c'] pq,
DTYPE_INT32_t age,
np.ndarray[DTYPE_INT32_t,ndim=2,negative_indices=False, mode='c'] structure,
DTYPE_INT32_t ndim,
np.ndarray[DTYPE_BOOL_t,ndim=1,negative_indices=False, mode='c'] mask,
np.ndarray[DTYPE_INT32_t,ndim=1,negative_indices=False, mode='c'] image_shape,
np.ndarray[DTYPE_INT32_t,ndim=1,negative_indices=False, mode='c'] output):
"""Do heavy lifting of watershed algorithm
image - the flattened image pixels, converted to rank-order
pq - the priority queue, starts with the marked pixels
the first element in each row is the image intensity
the second element is the age at entry into the queue
the third element is the index into the flattened image or labels
the remaining elements are the coordinates of the point
age - the next age to assign to a pixel
structure - a numpy int32 array containing the structuring elements
that define nearest neighbors. For each row, the first
element is the stride from the point to its neighbor
in a flattened array. The remaining elements are the
offsets from the point to its neighbor in the various
dimensions
ndim - # of dimensions in the image
mask - numpy boolean (char) array indicating which pixels to consider
and which to ignore. Also flattened.
image_shape - the dimensions of the image, for boundary checking,
a numpy array of np.int32
output - put the image labels in here
"""
cdef Heapitem elem
cdef Heapitem new_elem
cdef DTYPE_INT32_t nneighbors = structure.shape[0]
cdef DTYPE_INT32_t i = 0
cdef DTYPE_INT32_t index = 0
cdef DTYPE_INT32_t old_index = 0
cdef DTYPE_INT32_t max_index = image.shape[0]
cdef Heap *hp = <Heap *> heap_from_numpy2()
for i in range(pq.shape[0]):
elem.value = pq[i, 0]
elem.age = pq[i, 1]
elem.index = pq[i, 2]
heappush(hp, &elem)
while hp.items > 0:
#
# Pop off an item to work on
#
heappop(hp, &elem)
####################################################
# loop through each of the structuring elements
#
old_index = elem.index
for i in range(nneighbors):
# get the flattened address of the neighbor
index = structure[i,0]+old_index
if index < 0 or index >= max_index or output[index] or not mask[index]:
continue
new_elem.value = image[index]
new_elem.age = elem.age + 1
new_elem.index = index
age += 1
output[index] = output[old_index]
#
# Push the neighbor onto the heap to work on it later
#
heappush(hp, &new_elem)
heap_done(hp)
scikits/image/morphology/heap_general.pxi
+135
View File
@@ -0,0 +1,135 @@
"""
Originally part of CellProfiler, code licensed under both GPL and BSD
licenses.
Website: http://www.cellprofiler.org
Copyright (c) 2003-2009 Massachusetts Institute of Technology
Copyright (c) 2009-2011 Broad Institute
All rights reserved.
Original author: Lee Kamentsky
"""
cdef extern from "stdlib.h":
ctypedef unsigned long size_t
void free(void *ptr)
void *malloc(size_t size)
void *realloc(void *ptr, size_t size)
cdef struct Heap:
unsigned int items
unsigned int space
Heapitem *data
Heapitem **ptrs
cdef inline Heap *heap_from_numpy2():
cdef unsigned int k
cdef Heap *heap
heap = <Heap *> malloc(sizeof (Heap))
heap.items = 0
heap.space = 1000
heap.data = <Heapitem *> malloc(heap.space * sizeof(Heapitem))
heap.ptrs = <Heapitem **> malloc(heap.space * sizeof(Heapitem *))
for k in range(heap.space):
heap.ptrs[k] = heap.data + k
return heap
cdef inline void heap_done(Heap *heap):
free(heap.data)
free(heap.ptrs)
free(heap)
cdef inline void swap(unsigned int a, unsigned int b, Heap *h):
h.ptrs[a], h.ptrs[b] = h.ptrs[b], h.ptrs[a]
######################################################
# heappop - inlined
#
# pop an element off the heap, maintaining heap invariant
#
# Note: heap ordering is the same as python heapq, i.e., smallest first.
######################################################
cdef inline void heappop(Heap *heap,
Heapitem *dest):
cdef unsigned int i, smallest, l, r # heap indices
#
# Start by copying the first element to the destination
#
dest[0] = heap.ptrs[0][0]
heap.items -= 1
# if the heap is now empty, we can return, no need to fix heap.
if heap.items == 0:
return
#
# Move the last element in the heap to the first.
#
swap(0, heap.items, heap)
#
# Restore the heap invariant.
#
i = 0
smallest = i
while True:
# loop invariant here: smallest == i
# find smallest of (i, l, r), and swap it to i's position if necessary
l = i*2+1 #__left(i)
r = i*2+2 #__right(i)
if l < heap.items:
if smaller(heap.ptrs[l], heap.ptrs[i]):
smallest = l
if r < heap.items and smaller(heap.ptrs[r], heap.ptrs[smallest]):
smallest = r
else:
# this is unnecessary, but trims 0.04 out of 0.85 seconds...
break
# the element at i is smaller than either of its children, heap invariant restored.
if smallest == i:
break
# swap
swap(i, smallest, heap)
i = smallest
##################################################
# heappush - inlined
#
# push the element onto the heap, maintaining the heap invariant
#
# Note: heap ordering is the same as python heapq, i.e., smallest first.
##################################################
cdef inline void heappush(Heap *heap,
Heapitem *new_elem):
cdef unsigned int child = heap.items
cdef unsigned int parent
cdef unsigned int k
cdef Heapitem *new_data
# grow if necessary
if heap.items == heap.space:
heap.space = heap.space * 2
new_data = <Heapitem *> realloc(<void *> heap.data, <size_t> (heap.space * sizeof(Heapitem)))
heap.ptrs = <Heapitem **> realloc(<void *> heap.ptrs, <size_t> (heap.space * sizeof(Heapitem *)))
for k in range(heap.items):
heap.ptrs[k] = new_data + (heap.ptrs[k] - heap.data)
for k in range(heap.items, heap.space):
heap.ptrs[k] = new_data + k
heap.data = new_data
# insert new data at child
heap.ptrs[child][0] = new_elem[0]
heap.items += 1
# restore heap invariant, all parents <= children
while child>0:
parent = (child + 1) / 2 - 1 # __parent(i)
if smaller(heap.ptrs[child], heap.ptrs[parent]):
swap(parent, child, heap)
child = parent
else:
break
scikits/image/morphology/heap_watershed.pxi
+26
View File
@@ -0,0 +1,26 @@
"""
Originally part of CellProfiler, code licensed under both GPL and BSD
licenses.
Website: http://www.cellprofiler.org
Copyright (c) 2003-2009 Massachusetts Institute of Technology
Copyright (c) 2009-2011 Broad Institute
All rights reserved.
Original author: Lee Kamentsky
"""
import numpy as np
cimport numpy as np
cimport cython
cdef struct Heapitem:
np.int32_t value
np.int32_t age
np.int32_t index
cdef inline int smaller(Heapitem *a, Heapitem *b):
if a.value <> b.value:
return a.value < b.value
return a.age < b.age
include "heap_general.pxi"
scikits/image/morphology/setup.py
+3
View File
@@ -13,11 +13,14 @@ def configuration(parent_package='', top_path=None):
cython(['ccomp.pyx'], working_path=base_path)
cython(['cmorph.pyx'], working_path=base_path)
cython(['_watershed.pyx'], working_path=base_path)
config.add_extension('ccomp', sources=['ccomp.c'],
include_dirs=[get_numpy_include_dirs()])
config.add_extension('cmorph', sources=['cmorph.c'],
include_dirs=[get_numpy_include_dirs()])
config.add_extension('_watershed', sources=['_watershed.c'],
include_dirs=[get_numpy_include_dirs()])
return config
scikits/image/morphology/tests/test_watershed.py
+396
View File
@@ -0,0 +1,396 @@
"""test_watershed.py - tests the watershed function
Originally part of CellProfiler, code licensed under both GPL and BSD licenses.
Website: http://www.cellprofiler.org
Copyright (c) 2003-2009 Massachusetts Institute of Technology
Copyright (c) 2009-2011 Broad Institute
All rights reserved.
Original author: Lee Kamentsky
"""
#Portions of this test were taken from scipy's watershed test in test_ndimage.py
#
# Copyright (C) 2003-2005 Peter J. Verveer
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions
# are met:
#
# 1. Redistributions of source code must retain the above copyright
# notice, this list of conditions and the following disclaimer.
#
# 2. Redistributions in binary form must reproduce the above
# copyright notice, this list of conditions and the following
# disclaimer in the documentation and/or other materials provided
# with the distribution.
#
# 3. The name of the author may not be used to endorse or promote
# products derived from this software without specific prior
# written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS
# OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
# ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
# DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
# GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
# WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
import math
import time
import unittest
import numpy as np
import scipy.ndimage
from scikits.image.morphology.watershed import watershed,fast_watershed, \
_slow_watershed
eps = 1e-12
def diff(a, b):
if not isinstance(a, np.ndarray):
a = np.asarray(a)
if not isinstance(b, np.ndarray):
b = np.asarray(b)
if (0 in a.shape) and (0 in b.shape):
return 0.0
b[a==0] = 0
if (a.dtype in [np.complex64, np.complex128] or
b.dtype in [np.complex64, np.complex128]):
a = np.asarray(a, np.complex128)
b = np.asarray(b, np.complex128)
t = ((a.real - b.real)**2).sum() + ((a.imag - b.imag)**2).sum()
else:
a = np.asarray(a)
a = a.astype(np.float64)
b = np.asarray(b)
b = b.astype(np.float64)
t = ((a - b)**2).sum()
return math.sqrt(t)
class TestFastWatershed(unittest.TestCase):
eight = np.ones((3, 3),bool)
def test_watershed01(self):
"watershed 1"
data = np.array([[0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0],
[0, 1, 1, 1, 1, 1, 0],
[0, 1, 0, 0, 0, 1, 0],
[0, 1, 0, 0, 0, 1, 0],
[0, 1, 0, 0, 0, 1, 0],
[0, 1, 1, 1, 1, 1, 0],
[0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0]], np.uint8)
markers = np.array([[ -1, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 1, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0]],
np.int8)
out = watershed(data, markers,self.eight)
expected = np.array([[-1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1, -1],
[-1, 1, 1, 1, 1, 1, -1],
[-1, 1, 1, 1, 1, 1, -1],
[-1, 1, 1, 1, 1, 1, -1],
[-1, 1, 1, 1, 1, 1, -1],
[-1, 1, 1, 1, 1, 1, -1],
[-1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1, -1]])
error = diff(expected, out)
assert error < eps
out = _slow_watershed(data, markers, 8)
error = diff(expected, out)
assert error < eps
def test_watershed02(self):
"watershed 2"
data = np.array([[0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0],
[0, 1, 1, 1, 1, 1, 0],
[0, 1, 0, 0, 0, 1, 0],
[0, 1, 0, 0, 0, 1, 0],
[0, 1, 0, 0, 0, 1, 0],
[0, 1, 1, 1, 1, 1, 0],
[0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0]], np.uint8)
markers = np.array([[ -1, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 1, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0]],
np.int8)
out = fast_watershed(data, markers)
error = diff([[-1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1, -1],
[-1, -1, 1, 1, 1, -1, -1],
[-1, 1, 1, 1, 1, 1, -1],
[-1, 1, 1, 1, 1, 1, -1],
[-1, 1, 1, 1, 1, 1, -1],
[-1, -1, 1, 1, 1, -1, -1],
[-1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1, -1]], out)
self.failUnless(error < eps)
def test_watershed03(self):
"watershed 3"
data = np.array([[0, 0, 0, 0, 0, 0, 0],
[0, 1, 1, 1, 1, 1, 0],
[0, 1, 0, 1, 0, 1, 0],
[0, 1, 0, 1, 0, 1, 0],
[0, 1, 0, 1, 0, 1, 0],
[0, 1, 1, 1, 1, 1, 0],
[0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0]], np.uint8)
markers = np.array([[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 2, 0, 3, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, -1]],
np.int8)
out = fast_watershed(data, markers)
error = diff([[-1, -1, -1, -1, -1, -1, -1],
[-1, 0, 2, 0, 3, 0, -1],
[-1, 2, 2, 0, 3, 3, -1],
[-1, 2, 2, 0, 3, 3, -1],
[-1, 2, 2, 0, 3, 3, -1],
[-1, 0, 2, 0, 3, 0, -1],
[-1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1, -1]], out)
self.failUnless(error < eps)
def test_watershed04(self):
"watershed 4"
data = np.array([[0, 0, 0, 0, 0, 0, 0],
[0, 1, 1, 1, 1, 1, 0],
[0, 1, 0, 1, 0, 1, 0],
[0, 1, 0, 1, 0, 1, 0],
[0, 1, 0, 1, 0, 1, 0],
[0, 1, 1, 1, 1, 1, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0]], np.uint8)
markers = np.array([[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 2, 0, 3, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, -1]],
np.int8)
out = fast_watershed(data, markers, self.eight)
error = diff([[-1, -1, -1, -1, -1, -1, -1],
[-1, 2, 2, 0, 3, 3, -1],
[-1, 2, 2, 0, 3, 3, -1],
[-1, 2, 2, 0, 3, 3, -1],
[-1, 2, 2, 0, 3, 3, -1],
[-1, 2, 2, 0, 3, 3, -1],
[-1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1, -1]], out)
self.failUnless(error < eps)
def test_watershed05(self):
"watershed 5"
data = np.array([[0, 0, 0, 0, 0, 0, 0],
[0, 1, 1, 1, 1, 1, 0],
[0, 1, 0, 1, 0, 1, 0],
[0, 1, 0, 1, 0, 1, 0],
[0, 1, 0, 1, 0, 1, 0],
[0, 1, 1, 1, 1, 1, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0]], np.uint8)
markers = np.array([[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 3, 0, 2, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, -1]],
np.int8)
out = fast_watershed(data, markers, self.eight)
error = diff([[-1, -1, -1, -1, -1, -1, -1],
[-1, 3, 3, 0, 2, 2, -1],
[-1, 3, 3, 0, 2, 2, -1],
[-1, 3, 3, 0, 2, 2, -1],
[-1, 3, 3, 0, 2, 2, -1],
[-1, 3, 3, 0, 2, 2, -1],
[-1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1, -1]], out)
self.failUnless(error < eps)
def test_watershed06(self):
"watershed 6"
data = np.array([[0, 1, 0, 0, 0, 1, 0],
[0, 1, 0, 0, 0, 1, 0],
[0, 1, 0, 0, 0, 1, 0],
[0, 1, 1, 1, 1, 1, 0],
[0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0]], np.uint8)
markers = np.array([[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 1, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0],
[ -1, 0, 0, 0, 0, 0, 0]],
np.int8)
out = fast_watershed(data, markers, self.eight)
error = diff([[-1, 1, 1, 1, 1, 1, -1],
[-1, 1, 1, 1, 1, 1, -1],
[-1, 1, 1, 1, 1, 1, -1],
[-1, 1, 1, 1, 1, 1, -1],
[-1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1, -1]], out)
self.failUnless(error < eps)
def test_watershed07(self):
"A regression test of a competitive case that failed"
data = np.array([[255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255],
[255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255],
[255,255,255,255,255,204,204,204,204,204,204,255,255,255,255,255],
[255,255,255,204,204,183,153,153,153,153,183,204,204,255,255,255],
[255,255,204,183,153,141,111,103,103,111,141,153,183,204,255,255],
[255,255,204,153,111, 94, 72, 52, 52, 72, 94,111,153,204,255,255],
[255,255,204,153,111, 72, 39, 1, 1, 39, 72,111,153,204,255,255],
[255,255,204,183,141,111, 72, 39, 39, 72,111,141,183,204,255,255],
[255,255,255,204,183,141,111, 72, 72,111,141,183,204,255,255,255],
[255,255,255,255,204,183,141, 94, 94,141,183,204,255,255,255,255],
[255,255,255,255,255,204,153,103,103,153,204,255,255,255,255,255],
[255,255,255,255,204,183,141, 94, 94,141,183,204,255,255,255,255],
[255,255,255,204,183,141,111, 72, 72,111,141,183,204,255,255,255],
[255,255,204,183,141,111, 72, 39, 39, 72,111,141,183,204,255,255],
[255,255,204,153,111, 72, 39, 1, 1, 39, 72,111,153,204,255,255],
[255,255,204,153,111, 94, 72, 52, 52, 72, 94,111,153,204,255,255],
[255,255,204,183,153,141,111,103,103,111,141,153,183,204,255,255],
[255,255,255,204,204,183,153,153,153,153,183,204,204,255,255,255],
[255,255,255,255,255,204,204,204,204,204,204,255,255,255,255,255],
[255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255],
[255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255]])
mask = (data!=255)
markers = np.zeros(data.shape,int)
markers[6, 7] = 1
markers[14, 7] = 2
out = fast_watershed(data, markers, self.eight, mask=mask)
#
# The two objects should be the same size, except possibly for the
# border region
#
size1 = np.sum(out == 1)
size2 = np.sum(out == 2)
self.assertTrue(abs(size1-size2) <= 6)
def test_watershed08(self):
"The border pixels + an edge are all the same value"
data = np.array([[255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255],
[255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255],
[255,255,255,255,255,204,204,204,204,204,204,255,255,255,255,255],
[255,255,255,204,204,183,153,153,153,153,183,204,204,255,255,255],
[255,255,204,183,153,141,111,103,103,111,141,153,183,204,255,255],
[255,255,204,153,111, 94, 72, 52, 52, 72, 94,111,153,204,255,255],
[255,255,204,153,111, 72, 39, 1, 1, 39, 72,111,153,204,255,255],
[255,255,204,183,141,111, 72, 39, 39, 72,111,141,183,204,255,255],
[255,255,255,204,183,141,111, 72, 72,111,141,183,204,255,255,255],
[255,255,255,255,204,183,141, 94, 94,141,183,204,255,255,255,255],
[255,255,255,255,255,204,153,141,141,153,204,255,255,255,255,255],
[255,255,255,255,204,183,141, 94, 94,141,183,204,255,255,255,255],
[255,255,255,204,183,141,111, 72, 72,111,141,183,204,255,255,255],
[255,255,204,183,141,111, 72, 39, 39, 72,111,141,183,204,255,255],
[255,255,204,153,111, 72, 39, 1, 1, 39, 72,111,153,204,255,255],
[255,255,204,153,111, 94, 72, 52, 52, 72, 94,111,153,204,255,255],
[255,255,204,183,153,141,111,103,103,111,141,153,183,204,255,255],
[255,255,255,204,204,183,153,153,153,153,183,204,204,255,255,255],
[255,255,255,255,255,204,204,204,204,204,204,255,255,255,255,255],
[255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255],
[255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255]])
mask = (data != 255)
markers = np.zeros(data.shape,int)
markers[6,7] = 1
markers[14,7] = 2
out = fast_watershed(data, markers, self.eight, mask=mask)
#
# The two objects should be the same size, except possibly for the
# border region
#
size1 = np.sum(out == 1)
size2 = np.sum(out == 2)
self.assertTrue(abs(size1-size2) <= 6)
def test_watershed09(self):
"""Test on an image of reasonable size
This is here both for timing (does it take forever?) and to
ensure that the memory constraints are reasonable
"""
image = np.zeros((1000, 1000))
coords = np.random.uniform(0, 1000, (100, 2)).astype(int)
markers = np.zeros((1000, 1000),int)
idx = 1
for x,y in coords:
image[x,y] = 1
markers[x, y] = idx
idx += 1
image = scipy.ndimage.gaussian_filter(image, 4)
before = time.clock()
out = fast_watershed(image, markers, self.eight)
elapsed = time.clock() - before
print "Fast watershed ran a megapixel image in %f seconds"%(elapsed)
before = time.clock()
out = scipy.ndimage.watershed_ift(image.astype(np.uint16), markers, self.eight)
elapsed = time.clock() - before
print "Scipy watershed ran a megapixel image in %f seconds"%(elapsed)
scikits/image/morphology/watershed.py
+345
View File
@@ -0,0 +1,345 @@
"""watershed.py - watershed algorithm
This module implements a watershed algorithm that apportions pixels into
marked basins. The algorithm uses a priority queue to hold the pixels
with the metric for the priority queue being pixel value, then the time
of entry into the queue - this settles ties in favor of the closest marker.
Some ideas taken from
Soille, "Automated Basin Delineation from Digital Elevation Models Using
Mathematical Morphology", Signal Processing 20 (1990) 171-182.
The most important insight in the paper is that entry time onto the queue
solves two problems: a pixel should be assigned to the neighbor with the
largest gradient or, if there is no gradient, pixels on a plateau should
be split between markers on opposite sides.
Originally part of CellProfiler, code licensed under both GPL and BSD licenses.
Website: http://www.cellprofiler.org
Copyright (c) 2003-2009 Massachusetts Institute of Technology
Copyright (c) 2009-2011 Broad Institute
All rights reserved.
Original author: Lee Kamentsky
"""
from _heapq import heapify, heappush, heappop
import numpy as np
import scipy.ndimage
from ..filter import rank_order
import _watershed
def __get_strides_for_shape(shape):
"""Get the amount to multiply at each coord when converting to flat"""
lshape = list(shape)
lshape.reverse()
stride = [1]
for i in range(len(lshape) - 1):
stride.append(lshape[i] * stride[i])
stride.reverse()
return np.array(stride)
def fast_watershed(image, markers, connectivity=None, offset=None, mask=None):
"""
Return a matrix labeled using the watershed segmentation algorithm
Parameters
----------
image: ndarray (2-D, 3-D, ...) of integers
Data array where the lowest value points are labeled first.
markers: ndarray of the same shape as `image`
An array marking the basins with the values to be assigned in the
label matrix. Zero means not a marker. This array should be of an
integer type.
connectivity: ndarray, optional
An array with the same number of dimensions as `image` whose
non-zero elements indicate neighbors for connection.
Following the scipy convention, default is a one-connected array of
the dimension of the image.
offset: array_like of shape image.ndim, optional
offset of the connectivity (one offset per dimension)
mask: ndarray of bools or 0s and 1s, optional
Array of same shape as `image`. Only points at which mask == True
will be labeled.
Returns
-------
out: ndarray
A labeled matrix of the same type and shape as markers
Notes
-----
This function implements a watershed algorithm [1]_that apportions pixels
into marked basins. The algorithm uses a priority queue to hold the pixels
with the metric for the priority queue being pixel value, then the time of
entry into the queue - this settles ties in favor of the closest marker.
Some ideas taken from
Soille, "Automated Basin Delineation from Digital Elevation Models Using
Mathematical Morphology", Signal Processing 20 (1990) 171-182
The most important insight in the paper is that entry time onto the queue
solves two problems: a pixel should be assigned to the neighbor with the
largest gradient or, if there is no gradient, pixels on a plateau should
be split between markers on opposite sides.
This implementation converts all arguments to specific, lowest common
denominator types, then passes these to a C algorithm.
Markers can be determined manually, or automatically using for example
the local minima of the gradient of the image, or the local maxima of the
distance function to the background for separating overlapping objects
(see example).
References
----------
.. [1] http://en.wikipedia.org/wiki/Watershed_%28image_processing%29
.. [2] http://cmm.ensmp.fr/~beucher/wtshed.html
See also
--------
`regional_maximum`, `is_local_maximum`, `scipy.ndimage.label`
Examples
--------
The watershed algorithm is very useful to separate overlapping objects
>>> # Generate an initial image with two overlapping circles
>>> x, y = np.indices((80, 80))
>>> x1, y1, x2, y2 = 28, 28, 44, 52
>>> r1, r2 = 16, 20
>>> mask_circle1 = (x - x1)**2 + (y - y1)**2 < r1**2
>>> mask_circle2 = (x - x2)**2 + (y - y2)**2 < r2**2
>>> image = np.logical_or(mask_circle1, mask_circle2)
>>> # Now we want to separate the two objects in image
>>> # Generate the markers as local maxima of the distance
>>> # to the background
>>> from scipy import ndimage
>>> distance = ndimage.distance_transform_edt(image)
>>> local_maxi = regional_maximum(distance)
>>> markers = ndimage.label(local_maxi)[0]
>>> labels = watershed(-distance, markers, mask=image)
The algorithm works also for 3-D images, and can be used for example to
separate overlapping spheres.
"""
if connectivity == None:
c_connectivity = scipy.ndimage.generate_binary_structure(image.ndim, 1)
else:
c_connectivity = np.array(connectivity, bool)
if c_connectivity.ndim != image.ndim:
raise ValueError,"Connectivity dimension must be same as image"
if offset == None:
if any([x%2==0 for x in c_connectivity.shape]):
raise ValueError,"Connectivity array must have an unambiguous \
center"
#
# offset to center of connectivity array
#
offset = np.array(c_connectivity.shape) / 2
# pad the image, markers, and mask so that we can use the mask to
# keep from running off the edges
pads = offset
def pad(im):
new_im = np.zeros([i + 2*p for i, p in zip(im.shape, pads)], im.dtype)
new_im[[slice(p, -p, None) for p in pads]] = im
return new_im
if mask is not None:
mask = pad(mask)
else:
mask = pad(np.ones(image.shape, bool))
image = pad(image)
markers = pad(markers)
c_image = rank_order(image)[0].astype(np.int32)
c_markers = np.ascontiguousarray(markers, dtype=np.int32)
if c_markers.ndim != c_image.ndim:
raise ValueError,\
"markers (ndim=%d) must have same # of dimensions "\
"as image (ndim=%d)"%(c_markers.ndim, cimage.ndim)
if c_markers.shape != c_image.shape:
raise ValueError("image and markers must have the same shape")
if mask != None:
c_mask = np.ascontiguousarray(mask, dtype=bool)
if c_mask.ndim != c_markers.ndim:
raise ValueError, "mask must have same # of dimensions as image"
if c_markers.shape != c_mask.shape:
raise ValueError, "mask must have same shape as image"
c_markers[np.logical_not(mask)]=0
else:
c_mask = None
c_output = c_markers.copy()
#
# We pass a connectivity array that pre-calculates the stride for each
# neighbor.
#
# The result of this bit of code is an array with one row per
# point to be considered. The first column is the pre-computed stride
# and the second through last are the x,y...whatever offsets
# (to do bounds checking).
c = []
image_stride = __get_strides_for_shape(image.shape)
for i in range(np.product(c_connectivity.shape)):
multiplier = 1
offs = []
indexes = []
ignore = True
for j in range(len(c_connectivity.shape)):
elems = c_image.shape[j]
idx = (i / multiplier) % c_connectivity.shape[j]
off = idx - offset[j]
if off:
ignore = False
offs.append(off)
indexes.append(idx)
multiplier *= c_connectivity.shape[j]
if (not ignore) and c_connectivity.__getitem__(tuple(indexes)):
stride = np.dot(image_stride, np.array(offs))
offs.insert(0, stride)
c.append(offs)
c = np.array(c, np.int32)
pq, age = __heapify_markers(c_markers, c_image)
pq = np.ascontiguousarray(pq, dtype=np.int32)
if np.product(pq.shape) > 0:
# If nothing is labeled, the output is empty and we don't have to
# do anything
c_output = c_output.flatten()
if c_mask == None:
c_mask = np.ones(c_image.shape, np.int8).flatten()
else:
c_mask = c_mask.astype(np.int8).flatten()
_watershed.watershed(c_image.flatten(),
pq, age, c,
c_image.ndim,
c_mask,
np.array(c_image.shape,np.int32),
c_output)
c_output = c_output.reshape(c_image.shape)[[slice(1, -1, None)] *
image.ndim]
try:
return c_output.astype(markers.dtype)
except:
return c_output
watershed = fast_watershed
# ---------------------- deprecated ------------------------------
# Deprecate slower pure-Python code, that we keep only for
# pedagogical purposes
def __heapify_markers(markers,image):
"""Create a priority queue heap with the markers on it"""
stride = __get_strides_for_shape(image.shape)
coords = np.argwhere(markers != 0)
ncoords= coords.shape[0]
if ncoords > 0:
pixels = image[markers != 0]
age = np.array(range(ncoords))
offset = np.zeros(coords.shape[0], int)
for i in range(image.ndim):
offset = offset + stride[i] * coords[:,i]
pq = np.column_stack((pixels, age, offset, coords))
ordering = np.lexsort((age, pixels)) # pixels = top priority, age=second
pq = pq[ordering,:]
else:
pq = np.zeros((0, markers.ndim+3), int)
return (pq, ncoords)
def _slow_watershed(image, markers, connectivity=8, mask=None):
"""Return a matrix labeled using the watershed algorithm
Parameters
----------
image: 2-d ndarray of integers
a two-dimensional matrix where the lowest value points are
labeled first.
markers: 2-d ndarray of integers
a two-dimensional matrix marking the basins with the values
to be assigned in the label matrix. Zero means not a marker.
connectivity: {4, 8}, optional
either 4 for four-connected or 8 (default) for eight-connected
mask: 2-d ndarray of bools, optional
don't label points in the mask
Returns
-------
out: ndarray
A labeled matrix of the same type and shape as markers
Notes
-----
This function implements a watershed algorithm [1]_that apportions pixels
into marked basins. The algorithm uses a priority queue to hold the pixels
with the metric for the priority queue being pixel value, then the time of
entry into the queue - this settles ties in favor of the closest marker.
Some ideas taken from
Soille, "Automated Basin Delineation from Digital Elevation Models Using
Mathematical Morphology", Signal Processing 20 (1990) 171-182
The most important insight in the paper is that entry time onto the queue
solves two problems: a pixel should be assigned to the neighbor with the
largest gradient or, if there is no gradient, pixels on a plateau should
be split between markers on opposite sides.
This implementation converts all arguments to specific, lowest common
denominator types, then passes these to a C algorithm.
Markers can be determined manually, or automatically using for example
the local minima of the gradient of the image, or the local maxima of the
distance function to the background for separating overlapping objects.
"""
if connectivity not in (4, 8):
raise ValueError("Connectivity was %d: it should be either \
four or eight" %(connectivity))
image = np.array(image)
markers = np.array(markers)
labels = markers.copy()
max_x = markers.shape[0]
max_y = markers.shape[1]
if connectivity == 4:
connect_increments = ((1, 0), (0, 1), (-1, 0), (0, -1))
else:
connect_increments = ((1, 0), (1, 1), (0, 1), (-1, 1),
(-1, 0), (-1, -1), (0, -1), (1, -1))
pq, age = __heapify_markers(markers, image)
pq = pq.tolist()
#
# The second step pops a value off of the queue, then labels and pushes
# the neighbors
#
while len(pq):
pix_value, pix_age, ignore, pix_x, pix_y = heappop(pq)
pix_label = labels[pix_x, pix_y]
for xi, yi in connect_increments:
x = pix_x + xi
y = pix_y + yi
if x < 0 or y < 0 or x >= max_x or y >= max_y:
continue
if labels[x, y]:
continue
if mask != None and not mask[x, y]:
continue
# label the pixel
labels[x, y] = pix_label
# put the pixel onto the queue
heappush(pq, (image[x, y], age, 0, x, y))
age += 1
return labels