mirror of
https://github.com/wassname/scikit-image.git
synced 2026-07-02 01:37:54 +08:00
ENH: Add faster numpy hough version as well as a Cython version. Add appropriate setup.py files to build the Cython module. Fix hough tests.
Conflicts: scikits/image/setup.py scikits/image/transform/hough_transform.py scikits/image/transform/tests/test_hough_transform.py
This commit is contained in:
Stefan van der Walt
@@ -10,6 +10,7 @@ def configuration(parent_package='', top_path=None):
|
||||
config.add_subpackage('io')
|
||||
config.add_subpackage('morphology')
|
||||
config.add_subpackage('filter')
|
||||
config.add_subpackage('transform')
|
||||
|
||||
def add_test_directories(arg, dirname, fnames):
|
||||
if dirname.split(os.path.sep)[-1] == 'tests':
|
||||
|
||||
@@ -0,0 +1,63 @@
|
||||
cimport cython
|
||||
|
||||
import numpy as np
|
||||
cimport numpy as np
|
||||
|
||||
np.import_array()
|
||||
|
||||
|
||||
cdef extern from "math.h":
|
||||
double sqrt(double)
|
||||
double ceil(double)
|
||||
double round(double)
|
||||
|
||||
|
||||
cdef double PI_2 = 1.5707963267948966
|
||||
cdef double NEG_PI_2 = -PI_2
|
||||
|
||||
|
||||
@cython.boundscheck(False)
|
||||
def _hough(np.ndarray img, np.ndarray[ndim=1, dtype=np.double_t] theta=None):
|
||||
|
||||
if img.ndim != 2:
|
||||
raise ValueError('The input image must be 2D.')
|
||||
|
||||
# Compute the array of angles and their sine and cosine
|
||||
cdef np.ndarray[ndim=1, dtype=np.double_t] ctheta
|
||||
cdef np.ndarray[ndim=1, dtype=np.double_t] stheta
|
||||
|
||||
if theta is None:
|
||||
theta = np.linspace(PI_2, NEG_PI_2, 180)
|
||||
|
||||
ctheta = np.cos(theta)
|
||||
stheta = np.sin(theta)
|
||||
|
||||
# compute the bins and allocate the output array
|
||||
cdef np.ndarray[ndim=2, dtype=np.uint64_t] out
|
||||
cdef np.ndarray[ndim=1, dtype=np.double_t] bins
|
||||
cdef int max_distance, offset
|
||||
|
||||
max_distance = 2 * <int>ceil((sqrt(img.shape[0] * img.shape[0] +
|
||||
img.shape[1] * img.shape[1])))
|
||||
out = np.zeros((max_distance, theta.shape[0]), dtype=np.uint64)
|
||||
bins = np.linspace(-max_distance / 2.0, max_distance / 2.0, max_distance)
|
||||
offset = max_distance / 2
|
||||
|
||||
# compute the nonzero indexes
|
||||
cdef np.ndarray[ndim=1, dtype=np.int32_t] x_idxs, y_idxs
|
||||
y_idxs, x_idxs = np.PyArray_Nonzero(img)
|
||||
|
||||
# finally, run the transform
|
||||
cdef int nidxs, nthetas, i, j, x, y, out_idx
|
||||
nidxs = y_idxs.shape[0] # x and y are the same shape
|
||||
nthetas = theta.shape[0]
|
||||
for i in range(nidxs):
|
||||
x = x_idxs[i]
|
||||
y = y_idxs[i]
|
||||
for j in range(nthetas):
|
||||
out_idx = <int>round((ctheta[j] * x + stheta[j] * y)) + offset
|
||||
out[out_idx, j] += 1
|
||||
|
||||
return out, theta, bins
|
||||
|
||||
|
||||
@@ -1,52 +1,82 @@
|
||||
## Copyright (C) 2005 Stefan van der Walt <stefan@sun.ac.za>
|
||||
##
|
||||
## 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.
|
||||
##
|
||||
## 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.
|
||||
|
||||
__all__ = ['hough']
|
||||
|
||||
from itertools import izip
|
||||
import numpy as np
|
||||
|
||||
itype = np.uint16 # See ticket 225
|
||||
|
||||
def hough(img, angles=None):
|
||||
def _hough(img, theta=None):
|
||||
if img.ndim != 2:
|
||||
raise ValueError('The input image must be 2-D')
|
||||
|
||||
if theta is None:
|
||||
theta = np.linspace(-np.pi / 2, np.pi / 2, 180)
|
||||
|
||||
# compute the vertical bins (the distances)
|
||||
d = np.ceil(np.hypot(*img.shape))
|
||||
nr_bins = 2 * d
|
||||
bins = np.linspace(-d, d, nr_bins)
|
||||
|
||||
# allocate the output image
|
||||
out = np.zeros((nr_bins, len(theta)), dtype=np.uint64)
|
||||
|
||||
# precompute the sin and cos of the angles
|
||||
cos_theta = np.cos(theta)
|
||||
sin_theta = np.sin(theta)
|
||||
|
||||
# find the indices of the non-zero values in
|
||||
# the input image
|
||||
y, x = np.nonzero(img)
|
||||
|
||||
# x and y can be large, so we can't just broadcast to 2D
|
||||
# arrays as we may run out of memory. Instead we process
|
||||
# one vertical slice at a time.
|
||||
for i, (cT, sT) in enumerate(izip(cos_theta, sin_theta)):
|
||||
|
||||
# compute the base distances
|
||||
distances = x * cT + y * sT
|
||||
|
||||
# round the distances to the nearest integer
|
||||
# and shift them to a nonzero bin
|
||||
shifted = np.round(distances) - bins[0]
|
||||
|
||||
# cast the shifted values to ints to use as indices
|
||||
indices = shifted.astype(np.int)
|
||||
|
||||
# use bin count to accumulate the coefficients
|
||||
bincount = np.bincount(indices)
|
||||
|
||||
# finally assign the proper values to the out array
|
||||
out[:len(bincount), i] = bincount
|
||||
|
||||
return out, theta, bins
|
||||
|
||||
|
||||
# try to import and use the faster Cython version if it exists
|
||||
try:
|
||||
from ._hough_transform import _hough
|
||||
except ImportError:
|
||||
pass
|
||||
|
||||
|
||||
def hough(img, theta=None):
|
||||
"""Perform a straight line Hough transform.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
img : (M, N) bool ndarray
|
||||
Thresholded input image.
|
||||
angles : ndarray or list
|
||||
Angles at which to compute the transform.
|
||||
img : (M, N) ndarray
|
||||
Input image with nonzero values representing edges.
|
||||
theta :1D ndarray, dtype=double
|
||||
Angles at which to compute the transform, in radians.
|
||||
Defaults to -pi/2 - pi/2
|
||||
|
||||
Returns
|
||||
-------
|
||||
H : 2-D ndarray
|
||||
Hough transform coefficients.
|
||||
H : 2-D ndarray, uint64
|
||||
Hough transform accumulator.
|
||||
distances : ndarray
|
||||
Distance values.
|
||||
angles : ndarray
|
||||
Angle values.
|
||||
theta : ndarray
|
||||
Angles at which the transform was computed.
|
||||
|
||||
Examples
|
||||
--------
|
||||
@@ -62,7 +92,7 @@ def hough(img, angles=None):
|
||||
|
||||
Apply the Hough transform:
|
||||
|
||||
>>> out, angles, d = houghtf(img)
|
||||
>>> out, angles, d = hough(img)
|
||||
|
||||
Plot the results:
|
||||
|
||||
@@ -73,29 +103,4 @@ def hough(img, angles=None):
|
||||
>>> plt.show()
|
||||
|
||||
"""
|
||||
if img.ndim != 2:
|
||||
raise ValueError("Input must be a two-dimensional array")
|
||||
|
||||
img = img.astype(bool)
|
||||
|
||||
if angles is None:
|
||||
angles = np.linspace(-90,90,180)
|
||||
|
||||
theta = angles / 180. * np.pi
|
||||
d = np.ceil(np.hypot(*img.shape))
|
||||
nr_bins = 2*d - 1
|
||||
bins = np.linspace(-d, d, nr_bins)
|
||||
out = np.zeros((nr_bins, len(theta)), dtype=itype)
|
||||
|
||||
rows, cols = img.shape
|
||||
x,y = np.mgrid[:rows, :cols]
|
||||
|
||||
for i, (cT, sT) in enumerate(zip(np.cos(theta), np.sin(theta))):
|
||||
rho = np.round_(cT * x[img] + sT * y[img]) - bins[0] + 1
|
||||
rho = rho.astype(itype)
|
||||
rho[(rho < 0) | (rho > nr_bins)] = 0
|
||||
bc = np.bincount(rho.flat)[1:]
|
||||
out[:len(bc), i] = bc
|
||||
|
||||
return out, angles, bins
|
||||
|
||||
return _hough(img, theta)
|
||||
|
||||
@@ -0,0 +1,33 @@
|
||||
#!/usr/bin/env python
|
||||
|
||||
import os
|
||||
import shutil
|
||||
import hashlib
|
||||
|
||||
from scikits.image._build import cython
|
||||
|
||||
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('transform', parent_package, top_path)
|
||||
config.add_data_dir('tests')
|
||||
|
||||
cython(['_hough_transform.pyx'], working_path=base_path)
|
||||
|
||||
config.add_extension('_hough_transform', sources=['_hough_transform.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
|
||||
return config
|
||||
|
||||
if __name__ == '__main__':
|
||||
from numpy.distutils.core import setup
|
||||
setup(maintainer = 'Scikits.Image Developers',
|
||||
author = 'Scikits.Image Developers',
|
||||
maintainer_email = 'scikits-image@googlegroups.com',
|
||||
description = 'Transforms',
|
||||
url = 'http://stefanv.github.com/scikits.image/',
|
||||
license = 'SciPy License (BSD Style)',
|
||||
**(configuration(top_path='').todict())
|
||||
)
|
||||
@@ -5,17 +5,18 @@ from scikits.image.transform import *
|
||||
|
||||
def test_hough():
|
||||
# Generate a test image
|
||||
img = np.zeros((100, 150), dtype=bool)
|
||||
img[30, :] = 1
|
||||
img[:, 65] = 1
|
||||
img[35:45, 35:50] = 1
|
||||
for i in range(90):
|
||||
img[i, i] = 1
|
||||
img = np.zeros((100, 100), dtype=int)
|
||||
for i in range(25, 75):
|
||||
img[100 - i, i] = 1
|
||||
|
||||
out, angles, d = hough(img)
|
||||
|
||||
assert_equal(out.max(), 100)
|
||||
assert_equal(len(angles), 180)
|
||||
y, x = np.where(out == out.max())
|
||||
dist = d[y[0]]
|
||||
theta = angles[x[0]]
|
||||
|
||||
assert_equal(dist > 70, dist < 72)
|
||||
assert_equal(theta > 0.78, theta < 0.79)
|
||||
|
||||
def test_hough_angles():
|
||||
img = np.zeros((10, 10))
|
||||
|
||||
Reference in New Issue
Block a user