mirror of
https://github.com/wassname/scikit-image.git
synced 2026-07-18 12:40:14 +08:00
Merge pull request #75 from stefanv/convex_hull
ENH: Add binary convex hull computation.
This commit is contained in:
Stefan van der Walt
@@ -84,3 +84,5 @@
|
||||
- Nelle Varoquaux
|
||||
Renaming of the package to ``skimage``.
|
||||
|
||||
- W. Randolph Franklin
|
||||
Point in polygon test.
|
||||
|
||||
@@ -0,0 +1,35 @@
|
||||
"""
|
||||
===========
|
||||
Convex Hull
|
||||
===========
|
||||
|
||||
The convex hull of a binary image is the set of pixels included in the
|
||||
smallest convex polygon that surround all white pixels in the input.
|
||||
|
||||
In this example, we show how the input pixels (white) get filled in by the
|
||||
convex hull (white and grey).
|
||||
|
||||
A good overview of the algorithm is given on `Steve Eddin's blog
|
||||
<http://blogs.mathworks.com/steve/2011/10/04/binary-image-convex-hull-algorithm-notes/>`__.
|
||||
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
from skimage.morphology import convex_hull
|
||||
|
||||
image = np.array(
|
||||
[[0, 0, 0, 0, 0, 0, 0, 0, 0],
|
||||
[0, 0, 0, 0, 1, 0, 0, 0, 0],
|
||||
[0, 0, 0, 1, 0, 1, 0, 0, 0],
|
||||
[0, 0, 1, 0, 0, 0, 1, 0, 0],
|
||||
[0, 1, 0, 0, 0, 0, 0, 1, 0],
|
||||
[0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=float)
|
||||
|
||||
chull = convex_hull(image)
|
||||
image[chull] += 1.7
|
||||
image -= -1.7
|
||||
|
||||
plt.imshow(image, cmap=plt.cm.gray, interpolation='nearest')
|
||||
plt.show()
|
||||
@@ -3,3 +3,4 @@ from selem import *
|
||||
from .ccomp import label
|
||||
from watershed import watershed, is_local_maximum
|
||||
from skeletonize import skeletonize, medial_axis
|
||||
from .convex_hull import convex_hull_image
|
||||
|
||||
@@ -0,0 +1,59 @@
|
||||
# -*- python -*-
|
||||
|
||||
cimport numpy as np
|
||||
import numpy as np
|
||||
|
||||
def possible_hull(np.ndarray[dtype=np.uint8_t, ndim=2, mode="c"] img):
|
||||
"""Return positions of pixels that possibly belong to the convex hull.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
img : ndarray of bool
|
||||
Binary input image.
|
||||
|
||||
Returns
|
||||
-------
|
||||
coords : ndarray (N, 2)
|
||||
The ``(row, column)`` coordinates of all pixels that possibly belong to
|
||||
the convex hull.
|
||||
|
||||
"""
|
||||
cdef int i, j, k
|
||||
cdef unsigned int M, N
|
||||
|
||||
M = img.shape[0]
|
||||
N = img.shape[1]
|
||||
|
||||
# Output: M storage slots for left boundary pixels
|
||||
# N storage slots for top boundary pixels
|
||||
# M storage slots for right boundary pixels
|
||||
# N storage slots for bottom boundary pixels
|
||||
cdef np.ndarray[dtype=np.int_t, ndim=2] nonzero = \
|
||||
np.ones((2 * (M + N), 2), dtype=np.int)
|
||||
nonzero *= -1
|
||||
|
||||
k = 0
|
||||
for i in range(M):
|
||||
for j in range(N):
|
||||
if img[i, j] != 0:
|
||||
# Left check
|
||||
if nonzero[i, 1] == -1:
|
||||
nonzero[i, 0] = i
|
||||
nonzero[i, 1] = j
|
||||
|
||||
# Right check
|
||||
elif nonzero[M + N + i, 1] < j:
|
||||
nonzero[M + N + i, 0] = i
|
||||
nonzero[M + N + i, 1] = j
|
||||
|
||||
# Top check
|
||||
if nonzero[M + j, 1] == -1:
|
||||
nonzero[M + j, 0] = i
|
||||
nonzero[M + j, 1] = j
|
||||
|
||||
# Bottom check
|
||||
elif nonzero[2 * M + N + j, 0] < i:
|
||||
nonzero[2 * M + N + j, 0] = i
|
||||
nonzero[2 * M + N + j, 1] = j
|
||||
|
||||
return nonzero[nonzero[:, 0] != -1]
|
||||
@@ -0,0 +1,72 @@
|
||||
/* `pnpoly` is from
|
||||
|
||||
http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
|
||||
|
||||
Copyright (c) 1970-2003, Wm. Randolph Franklin
|
||||
|
||||
Permission is hereby granted, free of charge, to any person
|
||||
obtaining a copy of this software and associated documentation
|
||||
files (the "Software"), to deal in the Software without
|
||||
restriction, including without limitation the rights to use, copy,
|
||||
modify, merge, publish, distribute, sublicense, and/or sell copies
|
||||
of the Software, and to permit persons to whom the Software is
|
||||
furnished to do so, subject to the following conditions:
|
||||
|
||||
1. Redistributions of source code must retain the above copyright
|
||||
notice, this list of conditions and the following disclaimers.
|
||||
2. Redistributions in binary form must reproduce the above
|
||||
copyright notice in the documentation and/or other materials
|
||||
provided with the distribution.
|
||||
3. The name of W. Randolph Franklin may not be used to endorse or
|
||||
promote products derived from this Software without specific
|
||||
prior written permission.
|
||||
|
||||
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
||||
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
|
||||
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
||||
NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
|
||||
BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
|
||||
ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
|
||||
CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
|
||||
SOFTWARE. */
|
||||
|
||||
#ifdef __cplusplus
|
||||
extern "C" {
|
||||
#endif
|
||||
|
||||
unsigned char pnpoly(int nr_verts, double *xp, double *yp, double x, double y)
|
||||
{
|
||||
int i, j;
|
||||
unsigned char c = 0;
|
||||
for (i = 0, j = nr_verts-1; i < nr_verts; j = i++) {
|
||||
if ((((yp[i]<=y) && (y<yp[j])) ||
|
||||
((yp[j]<=y) && (y<yp[i]))) &&
|
||||
(x < (xp[j] - xp[i]) * (y - yp[i]) / (yp[j] - yp[i]) + xp[i]))
|
||||
|
||||
c = !c;
|
||||
}
|
||||
return c;
|
||||
}
|
||||
|
||||
void npnpoly(int nr_verts, double *xp, double *yp,
|
||||
int nr_points, double *x, double *y,
|
||||
unsigned char *result)
|
||||
/*
|
||||
* For N provided points, calculate whether they are in
|
||||
* the polygon defined by vertices *xp, *yp.
|
||||
*
|
||||
* nr_verts : number of vertices
|
||||
* *xp, *yp : x and y coordinates of vertices
|
||||
* nr_points : number of data points provided
|
||||
* *x, *y : data points
|
||||
*/
|
||||
{
|
||||
unsigned char n = 0;
|
||||
for (n = 0; n < nr_points; n++) {
|
||||
result[n] = pnpoly(nr_verts, xp, yp, x[n], y[n]);
|
||||
}
|
||||
}
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
@@ -0,0 +1,93 @@
|
||||
# -*- python -*-
|
||||
|
||||
cimport numpy as np
|
||||
import numpy as np
|
||||
|
||||
cdef extern from "_pnpoly.h":
|
||||
int pnpoly(int nr_verts, double *xp, double *yp,
|
||||
double x, double y)
|
||||
|
||||
void npnpoly(int nr_verts, double *xp, double *yp,
|
||||
int nr_points, double *x, double *y,
|
||||
unsigned char *result)
|
||||
|
||||
|
||||
def grid_points_inside_poly(shape, verts):
|
||||
"""Test whether points on a specified grid are inside a polygon.
|
||||
|
||||
For each ``(r, c)`` coordinate on a grid, i.e. ``(0, 0)``, ``(0, 1)`` etc.,
|
||||
test whether that point lies inside a polygon.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
shape : tuple (M, N)
|
||||
Shape of the grid.
|
||||
verts : (V, 2) array
|
||||
Specify the V vertices of the polygon, sorted either clockwise
|
||||
or anti-clockwise. The first point may (but does not need to be)
|
||||
duplicated.
|
||||
|
||||
Returns
|
||||
-------
|
||||
mask : (M, N) ndarray of bool
|
||||
True where the grid falls inside the polygon.
|
||||
|
||||
"""
|
||||
cdef np.ndarray[np.double_t, ndim=1, mode="c"] vx, vy
|
||||
verts = np.asarray(verts)
|
||||
|
||||
vx = verts[:, 0].astype(np.double)
|
||||
vy = verts[:, 1].astype(np.double)
|
||||
cdef int V = vx.shape[0]
|
||||
|
||||
cdef int M = shape[0]
|
||||
cdef int N = shape[1]
|
||||
cdef int m, n
|
||||
|
||||
cdef np.ndarray[dtype=np.uint8_t, ndim=2, mode="c"] out = \
|
||||
np.zeros((M, N), dtype=np.uint8)
|
||||
|
||||
for m in range(M):
|
||||
for n in range(N):
|
||||
out[m, n] = pnpoly(V, <double*>vx.data, <double*>vy.data, m, n)
|
||||
|
||||
return out.view(bool)
|
||||
|
||||
|
||||
def points_inside_poly(points, verts):
|
||||
"""Test whether points lie inside a polygon.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
points : (N, 2) array
|
||||
Input points, ``(x, y)``.
|
||||
verts : (M, 2) array
|
||||
Vertices of the polygon, sorted either clockwise or anti-clockwise.
|
||||
The first point may (but does not need to be) duplicated.
|
||||
|
||||
Returns
|
||||
-------
|
||||
mask : (N,) array of bool
|
||||
True if corresponding point is inside the polygon.
|
||||
|
||||
"""
|
||||
cdef np.ndarray[np.double_t, ndim=1, mode="c"] x, y, vx, vy
|
||||
|
||||
points = np.asarray(points)
|
||||
verts = np.asarray(verts)
|
||||
|
||||
x = points[:, 0].astype(np.double)
|
||||
y = points[:, 1].astype(np.double)
|
||||
|
||||
vx = verts[:, 0].astype(np.double)
|
||||
vy = verts[:, 1].astype(np.double)
|
||||
|
||||
cdef np.ndarray[np.uint8_t, ndim=1] out = \
|
||||
np.zeros(x.shape[0], dtype=np.uint8)
|
||||
|
||||
npnpoly(vx.shape[0], <double*>vx.data, <double*>vy.data,
|
||||
x.shape[0], <double*>x.data, <double*>y.data,
|
||||
<unsigned char*>out.data)
|
||||
|
||||
return out.astype(bool)
|
||||
|
||||
@@ -0,0 +1,65 @@
|
||||
__all__ = ['convex_hull_image']
|
||||
|
||||
import numpy as np
|
||||
from ._pnpoly import points_inside_poly, grid_points_inside_poly
|
||||
from ._convex_hull import possible_hull
|
||||
|
||||
def convex_hull_image(image):
|
||||
"""Compute the convex hull image of a binary image.
|
||||
|
||||
The convex hull is the set of pixels included in the smallest convex
|
||||
polygon that surround all white pixels in the input image.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Binary input image. This array is cast to bool before processing.
|
||||
|
||||
Returns
|
||||
-------
|
||||
hull : ndarray of uint8
|
||||
Binary image with pixels in convex hull set to 255.
|
||||
|
||||
References
|
||||
----------
|
||||
.. [1] http://blogs.mathworks.com/steve/2011/10/04/binary-image-convex-hull-algorithm-notes/
|
||||
|
||||
"""
|
||||
|
||||
image = image.astype(bool)
|
||||
|
||||
# Here we do an optimisation by choosing only pixels that are
|
||||
# the starting or ending pixel of a row or column. This vastly
|
||||
# limits the number of coordinates to examine for the virtual
|
||||
# hull.
|
||||
coords = possible_hull(image.astype(np.uint8))
|
||||
N = len(coords)
|
||||
|
||||
# Add a vertex for the middle of each pixel edge
|
||||
coords_corners = np.empty((N * 4, 2))
|
||||
for i, (x_offset, y_offset) in enumerate(zip((0, 0, -0.5, 0.5),
|
||||
(-0.5, 0.5, 0, 0))):
|
||||
coords_corners[i * N:(i + 1) * N] = coords + [x_offset, y_offset]
|
||||
|
||||
coords = coords_corners
|
||||
|
||||
try:
|
||||
from scipy.spatial import Delaunay
|
||||
except ImportError:
|
||||
raise ImportError('Could not import scipy.spatial, only available in '
|
||||
'scipy >= 0.9.')
|
||||
|
||||
# Find the convex hull
|
||||
chull = Delaunay(coords).convex_hull
|
||||
v = coords[np.unique(chull)]
|
||||
|
||||
# Sort vertices clock-wise
|
||||
v_centred = v - v.mean(axis=0)
|
||||
angles = np.arctan2(v_centred[:, 0], v_centred[:, 1])
|
||||
v = v[np.argsort(angles)]
|
||||
|
||||
# For each pixel coordinate, check whether that pixel
|
||||
# lies inside the convex hull
|
||||
mask = grid_points_inside_poly(image.shape[:2], v)
|
||||
|
||||
return mask
|
||||
@@ -15,6 +15,8 @@ def configuration(parent_package='', top_path=None):
|
||||
cython(['cmorph.pyx'], working_path=base_path)
|
||||
cython(['_watershed.pyx'], working_path=base_path)
|
||||
cython(['_skeletonize.pyx'], working_path=base_path)
|
||||
cython(['_pnpoly.pyx'], working_path=base_path)
|
||||
cython(['_convex_hull.pyx'], working_path=base_path)
|
||||
|
||||
config.add_extension('ccomp', sources=['ccomp.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
@@ -24,7 +26,10 @@ def configuration(parent_package='', top_path=None):
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
config.add_extension('_skeletonize', sources=['_skeletonize.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
|
||||
config.add_extension('_pnpoly', sources=['_pnpoly.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
config.add_extension('_convex_hull', sources=['_convex_hull.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
|
||||
return config
|
||||
|
||||
|
||||
@@ -0,0 +1,67 @@
|
||||
import numpy as np
|
||||
from numpy.testing import assert_array_equal
|
||||
from numpy.testing.decorators import skipif
|
||||
from skimage.morphology import convex_hull_image
|
||||
from skimage.morphology._convex_hull import possible_hull
|
||||
|
||||
try:
|
||||
import scipy.spatial
|
||||
scipy_spatial = True
|
||||
except ImportError:
|
||||
scipy_spatial = False
|
||||
|
||||
@skipif(not scipy_spatial)
|
||||
def test_basic():
|
||||
image = np.array(
|
||||
[[0, 0, 0, 0, 0, 0, 0, 0, 0],
|
||||
[0, 0, 0, 0, 1, 0, 0, 0, 0],
|
||||
[0, 0, 0, 1, 0, 1, 0, 0, 0],
|
||||
[0, 0, 1, 0, 0, 0, 1, 0, 0],
|
||||
[0, 1, 0, 0, 0, 0, 0, 1, 0],
|
||||
[0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=bool)
|
||||
|
||||
expected = np.array(
|
||||
[[0, 0, 0, 0, 0, 0, 0, 0, 0],
|
||||
[0, 0, 0, 0, 1, 0, 0, 0, 0],
|
||||
[0, 0, 0, 1, 1, 1, 0, 0, 0],
|
||||
[0, 0, 1, 1, 1, 1, 1, 0, 0],
|
||||
[0, 1, 1, 1, 1, 1, 1, 1, 0],
|
||||
[0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=bool)
|
||||
|
||||
assert_array_equal(convex_hull_image(image), expected)
|
||||
|
||||
@skipif(not scipy_spatial)
|
||||
def test_possible_hull():
|
||||
image = np.array(
|
||||
[[0, 0, 0, 0, 0, 0, 0, 0, 0],
|
||||
[0, 0, 0, 0, 1, 0, 0, 0, 0],
|
||||
[0, 0, 0, 1, 0, 1, 0, 0, 0],
|
||||
[0, 0, 1, 1, 1, 1, 1, 0, 0],
|
||||
[0, 1, 1, 1, 1, 1, 1, 1, 0],
|
||||
[0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=np.uint8)
|
||||
|
||||
expected = np.array([[1, 4],
|
||||
[2, 3],
|
||||
[3, 2],
|
||||
[4, 1],
|
||||
[4, 1],
|
||||
[3, 2],
|
||||
[2, 3],
|
||||
[1, 4],
|
||||
[2, 5],
|
||||
[3, 6],
|
||||
[4, 7],
|
||||
[2, 5],
|
||||
[3, 6],
|
||||
[4, 7],
|
||||
[4, 2],
|
||||
[4, 3],
|
||||
[4, 4],
|
||||
[4, 5],
|
||||
[4, 6]])
|
||||
|
||||
ph = possible_hull(image)
|
||||
assert_array_equal(ph, expected)
|
||||
|
||||
if __name__ == "__main__":
|
||||
np.testing.run_module_suite()
|
||||
@@ -0,0 +1,38 @@
|
||||
import numpy as np
|
||||
from numpy.testing import assert_array_equal
|
||||
|
||||
from skimage.morphology._pnpoly import points_inside_poly, \
|
||||
grid_points_inside_poly
|
||||
|
||||
class test_npnpoly():
|
||||
def test_square(self):
|
||||
v = np.array([[0, 0],
|
||||
[0, 1],
|
||||
[1, 1],
|
||||
[1, 0]])
|
||||
assert(points_inside_poly([[0.5, 0.5]], v)[0])
|
||||
assert(not points_inside_poly([[-0.1, 0.1]], v)[0])
|
||||
|
||||
def test_triangle(self):
|
||||
v = np.array([[0, 0],
|
||||
[1, 0],
|
||||
[0.5, 0.75]])
|
||||
assert(points_inside_poly([[0.5, 0.7]], v)[0])
|
||||
assert(not points_inside_poly([[0.5, 0.76]], v)[0])
|
||||
assert(not points_inside_poly([[0.7, 0.5]], v)[0])
|
||||
|
||||
def test_type(self):
|
||||
assert(points_inside_poly([[0, 0]], [[0, 0]]).dtype == np.bool)
|
||||
|
||||
def test_grid_points_inside_poly():
|
||||
v = np.array([[0, 0],
|
||||
[5, 0],
|
||||
[5, 5]])
|
||||
|
||||
expected = np.tril(np.ones((5, 5), dtype=bool))
|
||||
|
||||
assert_array_equal(grid_points_inside_poly((5, 5), v),
|
||||
expected)
|
||||
|
||||
if __name__ == "__main__":
|
||||
np.testing.run_module_suite()
|
||||
Reference in New Issue
Block a user