diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index c0816717..007f267a 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -28,7 +28,7 @@ Code to generate skimage logo. - Zachary Pincus - Tracing of low cost paths, FreeImage I/O plugin + Tracing of low cost paths, FreeImage I/O plugin, iso-contours - Almar Klein Binary heap class for graph algorithms diff --git a/doc/examples/plot_contours.py b/doc/examples/plot_contours.py new file mode 100644 index 00000000..020d3dfc --- /dev/null +++ b/doc/examples/plot_contours.py @@ -0,0 +1,42 @@ +""" +=============== +Contour finding +=============== + +``skimage.measure.find_contours`` uses a marching squares method to find +constant valued contours in an image. Array values are linearly interpolated +to provide better precision of the output contours. Contours which intersect +the image edge are open; all others are closed. + +The `marching squares algorithm +`__ is a special case of +the marching cubes algorithm (Lorensen, William and Harvey E. Cline. Marching +Cubes: A High Resolution 3D Surface Construction Algorithm. Computer Graphics +(SIGGRAPH 87 Proceedings) 21(4) July 1987, p. 163-170). + +""" + +from skimage import data +from skimage import measure + +import numpy as np +import matplotlib.pyplot as plt + +# Construct some test data +x, y = np.ogrid[-np.pi:np.pi:100j, -np.pi:np.pi:100j] +r = np.sin(np.exp((np.sin(x)**3 + np.cos(y)**2))) + +# Find contours at a constant value of 0.8 +contours = measure.find_contours(r, 0.8) + +# Display the image and plot all contours found +plt.imshow(r, interpolation='nearest') + +for n, contour in enumerate(contours): + plt.plot(contour[:, 1], contour[:, 0], linewidth=2) + +plt.axis('image') +plt.xticks([]) +plt.yticks([]) +plt.show() + diff --git a/skimage/measure/__init__.py b/skimage/measure/__init__.py new file mode 100755 index 00000000..a9eb1526 --- /dev/null +++ b/skimage/measure/__init__.py @@ -0,0 +1 @@ +from find_contours import find_contours diff --git a/skimage/measure/_find_contours.pyx b/skimage/measure/_find_contours.pyx new file mode 100644 index 00000000..17702f95 --- /dev/null +++ b/skimage/measure/_find_contours.pyx @@ -0,0 +1,186 @@ +# -*- python -*- +# cython: cdivision=True + +import numpy as np +cimport numpy as np + +np.import_array() + +cdef inline double _get_fraction(double from_value, double to_value, + double level): + if (to_value == from_value): + return 0 + return ((level - from_value) / (to_value - from_value)) + + +def iterate_and_store(np.ndarray[double, ndim=2, mode='c'] array, + double level, int vertex_connect_high): + """Iterate across the given array in a marching-squares fashion, + looking for segments that cross 'level'. If such a segment is + found, its coordinates are added to a growing list of segments, + which is returned by the function. if vertex_connect_high is + nonzero, high-values pixels are considered to be face+vertex + connected into objects; otherwise low-valued pixels are. + + """ + if array.shape[0] < 2 or array.shape[1] < 2: + raise ValueError("Input array must be at least 2x2.") + + cdef list arc_list = [] + cdef int n + + # The plan is to iterate a 2x2 square across the input array. This means + # that the upper-left corner of the square needs to iterate across a + # sub-array that's one-less-large in each direction (so that the square + # never steps out of bounds). The square is represented by four pointers: + # ul, ur, ll, and lr (for 'upper left', etc.). We also maintain the current + # 2D coordinates for the position of the upper-left pointer. Note that we + # ensured that the array is of type 'double' and is C-contiguous (last + # index varies the fastest). + + # Current coords start at 0,0. + cdef int[2] coords + coords[0] = 0 + coords[1] = 0 + + # Calculate the number of iterations we'll need + cdef int num_square_steps = (array.shape[0] - 1) * (array.shape[1] - 1) + + cdef unsigned char square_case = 0 + cdef tuple top, bottom, left, right + cdef double ul, ur, ll, lr + cdef int r0, r1, c0, c1 + + for n in range(num_square_steps): + # There are sixteen different possible square types, diagramed below. + # A + indicates that the vertex is above the contour value, and a - + # indicates that the vertex is below or equal to the contour value. + # The vertices of each square are: + # ul ur + # ll lr + # and can be treated as a binary value with the bits in that order. Thus + # each square case can be numbered: + # 0-- 1+- 2-+ 3++ 4-- 5+- 6-+ 7++ + # -- -- -- -- +- +- +- +- + # + # 8-- 9+- 10-+ 11++ 12-- 13+- 14-+ 15++ + # -+ -+ -+ -+ ++ ++ ++ ++ + # + # The position of the line segment that cuts through (or + # doesn't, in case 0 and 15) each square is clear, except in + # cases 6 and 9. In this case, where the segments are placed + # is determined by vertex_connect_high. If + # vertex_connect_high is false, then lines like \\ are drawn + # through square 6, and lines like // are drawn through square + # 9. Otherwise, the situation is reversed. + # Finally, recall that we draw the lines so that (moving from tail to + # head) the lower-valued pixels are on the left of the line. So, for + # example, case 1 entails a line slanting from the middle of the top of + # the square to the middle of the left side of the square. + + r0, c0 = coords[0], coords[1] + r1, c1 = r0 + 1, c0 + 1 + + ul = array[r0, c0] + ur = array[r0, c1] + ll = array[r1, c0] + lr = array[r1, c1] + + # now in advance the coords indices + if coords[1] < array.shape[1] - 2: + coords[1] += 1 + else: + coords[0] += 1 + coords[1] = 0 + + + square_case = 0 + if (ul > level): square_case += 1 + if (ur > level): square_case += 2 + if (ll > level): square_case += 4 + if (lr > level): square_case += 8 + + if (square_case != 0 and square_case != 15): + # only do anything if there's a line passing through the + # square. Cases 0 and 15 are entirely below/above the contour. + + top = r0, c0 + _get_fraction(ul, ur, level) + bottom = r1, c0 + _get_fraction(ll, lr, level) + left = r0 + _get_fraction(ul, ll, level), c0 + right = r0 + _get_fraction(ur, lr, level), c1 + + if (square_case == 1): + # top to left + arc_list.append(top) + arc_list.append(left) + elif (square_case == 2): + # right to top + arc_list.append(right) + arc_list.append(top) + elif (square_case == 3): + # right to left + arc_list.append(right) + arc_list.append(left) + elif (square_case == 4): + # left to bottom + arc_list.append(left) + arc_list.append(bottom) + elif (square_case == 5): + # top to bottom + arc_list.append(top) + arc_list.append(bottom) + elif (square_case == 6): + if vertex_connect_high: + arc_list.append(left) + arc_list.append(top) + + arc_list.append(right) + arc_list.append(bottom) + else: + arc_list.append(right) + arc_list.append(top) + arc_list.append(left) + arc_list.append(bottom) + elif (square_case == 7): + # right to bottom + arc_list.append(right) + arc_list.append(bottom) + elif (square_case == 8): + # bottom to right + arc_list.append(bottom) + arc_list.append(right) + elif (square_case == 9): + if vertex_connect_high: + arc_list.append(top) + arc_list.append(right) + + arc_list.append(bottom) + arc_list.append(left) + else: + arc_list.append(top) + arc_list.append(left) + + arc_list.append(bottom) + arc_list.append(right) + elif (square_case == 10): + # bottom to top + arc_list.append(bottom) + arc_list.append(top) + elif (square_case == 11): + # bottom to left + arc_list.append(bottom) + arc_list.append(left) + elif (square_case == 12): + # lef to right + arc_list.append(left) + arc_list.append(right) + elif (square_case == 13): + # top to right + arc_list.append(top) + arc_list.append(right) + elif (square_case == 14): + # left to top + arc_list.append(left) + arc_list.append(top) + + return arc_list diff --git a/skimage/measure/find_contours.py b/skimage/measure/find_contours.py new file mode 100755 index 00000000..582c1b52 --- /dev/null +++ b/skimage/measure/find_contours.py @@ -0,0 +1,194 @@ +import numpy as np +import _find_contours + +from collections import deque + +_param_options = ('high', 'low') + +def find_contours(array, level, + fully_connected='low', positive_orientation='low'): + """Find iso-valued contours in a 2D array for a given level value. + + Uses the "marching squares" method to compute a the iso-valued contours of + the input 2D array for a particular level value. Array values are linearly + interpolated to provide better precision for the output contours. + + Parameters + ---------- + array : 2D ndarray of double + Input data in which to find contours. + level : float + Value along which to find contours in the array. + fully_connected : str, {'low', 'high'} + Indicates whether array elements below the given level value are to be + considered fully-connected (and hence elements above the value will + only be face connected), or vice-versa. (See notes below for details.) + positive_orientation : either 'low' or 'high' + Indicates whether the output contours will produce positively-oriented + polygons around islands of low- or high-valued elements. If 'low' then + contours will wind counter- clockwise around elements below the + iso-value. Alternately, this means that low-valued elements are always + on the left of the contour. (See below for details.) + + Returns + ------- + contours : list of (n,2)-ndarrays + Each contour is an ndarray of shape ``(n, 2)``, + consisting of n ``(x, y)`` coordinates along the contour. + + Notes + ----- + The marching squares algorithm is a special case of the marching cubes + algorithm [1]_. A simple explanation is available here:: + + http://www.essi.fr/~lingrand/MarchingCubes/algo.html + + There is a single ambiguous case in the marching squares algorithm: when + a given ``2 x 2``-element square has two high-valued and two low-valued + elements, each pair diagonally adjacent. (Where high- and low-valued is + with respect to the contour value sought.) In this case, either the + high-valued elements can be 'connected together' via a thin isthmus that + separates the low-valued elements, or vice-versa. When elements are + connected together across a diagonal, they are considered 'fully + connected' (also known as 'face+vertex-connected' or '8-connected'). Only + high-valued or low-valued elements can be fully-connected, the other set + will be considered as 'face-connected' or '4-connected'. By default, + low-valued elements are considered fully-connected; this can be altered + with the 'fully_connected' parameter. + + Output contours are not guaranteed to be closed: contours which intersect + the array edge will be left open. All other contours will be closed. (The + closed-ness of a contours can be tested by checking whether the beginning + point is the same as the end point.) + + Contours are oriented. By default, array values lower than the contour + value are to the left of the contour and values greater than the contour + value are to the right. This means that contours will wind + counter-clockwise (i.e. in 'positive orientation') around islands of + low-valued pixels. This behavior can be altered with the + 'positive_orientation' parameter. + + The order of the contours in the output list is determined by the position + of the smallest ``x,y`` (in lexicographical order) coordinate in the + contour. This is a side-effect of how the input array is traversed, but + can be relied upon. + + .. warning:: + + Array coordinates/values are assumed to refer to the *center* of the + array element. Take a simple example input: ``[0, 1]``. The interpolated + position of 0.5 in this array is midway between the 0-element (at + ``x=0``) and the 1-element (at ``x=1``), and thus would fall at + ``x=0.5``. + + This means that to find reasonable contours, it is best to find contours + midway between the expected "light" and "dark" values. In particular, + given a binarized array, *do not* choose to find contours at the low or high + value of the array. This will often yield degenerate contours, especially + around structures that are a single array element wide. Instead choose + a middle value, as above. + + References + ---------- + .. [1] Lorensen, William and Harvey E. Cline. Marching Cubes: A High + Resolution 3D Surface Construction Algorithm. Computer Graphics + (SIGGRAPH 87 Proceedings) 21(4) July 1987, p. 163-170). + + """ + array = np.asarray(array, dtype=np.double) + if array.ndim != 2: + raise RuntimeError('Only 2D arrays are supported.') + level = float(level) + if (fully_connected not in _param_options or + positive_orientation not in _param_options): + raise ValueError('Parameters "fully_connected" and' + ' "positive_orientation" must be either "high" or "low".') + point_list = _find_contours.iterate_and_store(array, level, + fully_connected == 'high') + contours = _assemble_contours(_take_2(point_list)) + if positive_orientation == 'high': + contours = [c[::-1] for c in contours] + return contours + +def _take_2(seq): + iterator = iter(seq) + while(True): + n1 = iterator.next() + n2 = iterator.next() + yield (n1, n2) + +def _assemble_contours(points_iterator): + current_index = 0 + contours = {} + starts = {} + ends = {} + for from_point, to_point in points_iterator: + # Ignore degenerate segments. + # This happens when (and only when) one vertex of the square is + # exactly the contour level, and the rest are above or below. + # This degnerate vertex will be picked up later by neighboring squares. + if from_point == to_point: continue + + tail_data = starts.get(to_point) + head_data = ends.get(from_point) + + if tail_data is not None and head_data is not None: + tail, tail_num = tail_data + head, head_num = head_data + # We need to connect these two contours. + if tail is head: + # We need to closed a contour. + # Add the end point, and remove the contour from the + # 'starts' and 'ends' dicts. + head.append(to_point) + del starts[to_point] + del ends[from_point] + else: # tail is not head + # We need to join two distinct contours. + # We want to keep the first contour segment created, so that + # the final contours are ordered left->right, top->bottom. + if tail_num > head_num: + # tail was created second. Append tail to head. + head.extend(tail) + # remove all traces of tail: + del starts[to_point] + del ends[tail[-1]] + del contours[tail_num] + # remove the old end of head and add the new end. + del ends[from_point] + ends[head[-1]] = (head, head_num) + else: # tail_num <= head_num + # head was created second. Prepend head to tail. + tail.extendleft(reversed(head)) + # remove all traces of head: + del starts[head[0]] + del ends[from_point] + del contours[head_num] + # remove the old start of tail and add the new start. + del starts[to_point] + starts[tail[0]] = (tail, tail_num) + elif tail_data is None and head_data is None: + # we need to add a new contour + current_index += 1 + new_num = current_index + new_contour = deque((from_point, to_point)) + contours[new_num] = new_contour + starts[from_point] = (new_contour, new_num) + ends[to_point] = (new_contour, new_num) + elif tail_data is not None and head_data is None: + tail, tail_num = tail_data + # We've found a single contour to which the new segment should be + # prepended. + tail.appendleft(from_point) + del starts[to_point] + starts[from_point] = (tail, tail_num) + elif tail_data is None and head_data is not None: + head, head_num = head_data + # We've found a single contour to which the new segment should be + # appended + head.append(to_point) + del ends[from_point] + ends[to_point] = (head, head_num) + # end iteration over from_ and to_ points + + return [np.array(contour) for (num, contour) in sorted(contours.items())] diff --git a/skimage/measure/setup.py b/skimage/measure/setup.py new file mode 100644 index 00000000..9fc5ca2d --- /dev/null +++ b/skimage/measure/setup.py @@ -0,0 +1,29 @@ +#!/usr/bin/env python + +from skimage._build import cython + +import os +base_path = os.path.abspath(os.path.dirname(__file__)) + +def configuration(parent_package='', top_path=None): + from numpy.distutils.misc_util import Configuration, get_numpy_include_dirs + + config = Configuration('measure', parent_package, top_path) + config.add_data_dir('tests') + + cython(['_find_contours.pyx'], working_path=base_path) + + config.add_extension('_find_contours', sources=['_find_contours.c'], + include_dirs=[get_numpy_include_dirs()]) + + return config + +if __name__ == '__main__': + from numpy.distutils.core import setup + setup(maintainer = 'scikits.image Developers', + maintainer_email = 'scikits-image@googlegroups.com', + description = 'Graph-based Image-processing Algorithms', + url = 'https://github.com/scikits-image/scikits.image', + license = 'Modified BSD', + **(configuration(top_path='').todict()) + ) diff --git a/skimage/measure/tests/test_find_contours.py b/skimage/measure/tests/test_find_contours.py new file mode 100644 index 00000000..8d705878 --- /dev/null +++ b/skimage/measure/tests/test_find_contours.py @@ -0,0 +1,66 @@ +import numpy as np +from numpy.testing import * + +from skimage.measure import find_contours + +a = np.ones((8,8), dtype=np.float32) +a[1:-1, 1] = 0 +a[1, 1:-1] = 0 + +## array([[ 1., 1., 1., 1., 1., 1., 1., 1.], +## [ 1., 0., 0., 0., 0., 0., 0., 1.], +## [ 1., 0., 1., 1., 1., 1., 1., 1.], +## [ 1., 0., 1., 1., 1., 1., 1., 1.], +## [ 1., 0., 1., 1., 1., 1., 1., 1.], +## [ 1., 0., 1., 1., 1., 1., 1., 1.], +## [ 1., 0., 1., 1., 1., 1., 1., 1.], +## [ 1., 1., 1., 1., 1., 1., 1., 1.]], dtype=float32) + +x,y = np.mgrid[-1:1:5j,-1:1:5j] +r = np.sqrt(x**2 + y**2) + +def test_binary(): + contours = find_contours(a, 0.5) + assert len(contours) == 1 + assert_array_equal(contours[0], + [[ 6. , 1.5], + [ 5. , 1.5], + [ 4. , 1.5], + [ 3. , 1.5], + [ 2. , 1.5], + [ 1.5, 2. ], + [ 1.5, 3. ], + [ 1.5, 4. ], + [ 1.5, 5. ], + [ 1.5, 6. ], + [ 1. , 6.5], + [ 0.5, 6. ], + [ 0.5, 5. ], + [ 0.5, 4. ], + [ 0.5, 3. ], + [ 0.5, 2. ], + [ 0.5, 1. ], + [ 1. , 0.5], + [ 2. , 0.5], + [ 3. , 0.5], + [ 4. , 0.5], + [ 5. , 0.5], + [ 6. , 0.5], + [ 6.5, 1. ], + [ 6. , 1.5]]) + +def test_float(): + contours = find_contours(r, 0.5) + assert len(contours) == 1 + assert_array_equal(contours[0], + [[ 2., 3.], + [ 1., 2.], + [ 2., 1.], + [ 3., 2.], + [ 2., 3.]]) + + + +if __name__ == '__main__': + from numpy.testing import run_module_suite + run_module_suite() diff --git a/skimage/setup.py b/skimage/setup.py index 35a04f6b..f4b09a56 100644 --- a/skimage/setup.py +++ b/skimage/setup.py @@ -15,6 +15,7 @@ def configuration(parent_package='', top_path=None): config.add_subpackage('color') config.add_subpackage('draw') config.add_subpackage('feature') + config.add_subpackage('measure') def add_test_directories(arg, dirname, fnames): if dirname.split(os.path.sep)[-1] == 'tests':