diff --git a/scikits/image/setup.py b/scikits/image/setup.py index d0f08285..b1831d4e 100644 --- a/scikits/image/setup.py +++ b/scikits/image/setup.py @@ -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': diff --git a/scikits/image/transform/_hough_transform.pyx b/scikits/image/transform/_hough_transform.pyx new file mode 100644 index 00000000..3373c0b3 --- /dev/null +++ b/scikits/image/transform/_hough_transform.pyx @@ -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 * 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 = round((ctheta[j] * x + stheta[j] * y)) + offset + out[out_idx, j] += 1 + + return out, theta, bins + + diff --git a/scikits/image/transform/hough_transform.py b/scikits/image/transform/hough_transform.py index 2342c2f0..3c9af50f 100644 --- a/scikits/image/transform/hough_transform.py +++ b/scikits/image/transform/hough_transform.py @@ -1,52 +1,82 @@ -## Copyright (C) 2005 Stefan van der Walt -## -## 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) diff --git a/scikits/image/transform/setup.py b/scikits/image/transform/setup.py new file mode 100644 index 00000000..5e377186 --- /dev/null +++ b/scikits/image/transform/setup.py @@ -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()) + ) diff --git a/scikits/image/transform/tests/test_hough_transform.py b/scikits/image/transform/tests/test_hough_transform.py index a8c31d58..bd9612b7 100644 --- a/scikits/image/transform/tests/test_hough_transform.py +++ b/scikits/image/transform/tests/test_hough_transform.py @@ -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))