From 71a54377940a3dd81e46d87ed27988ea9ed9b2ac Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Wed, 26 Nov 2014 17:03:50 -0500 Subject: [PATCH] Add fused type specialization for moments functions --- bento.info | 4 +- skimage/measure/{_moments.pyx => _moments.py} | 69 ++++--------------- skimage/measure/_moments_cy.pyx | 65 +++++++++++++++++ skimage/measure/setup.py | 4 +- 4 files changed, 81 insertions(+), 61 deletions(-) rename skimage/measure/{_moments.pyx => _moments.py} (67%) create mode 100644 skimage/measure/_moments_cy.pyx diff --git a/bento.info b/bento.info index 1ba9bbc3..6a8a0000 100644 --- a/bento.info +++ b/bento.info @@ -46,9 +46,9 @@ Library: Extension: skimage.measure._find_contours_cy Sources: skimage/measure/_find_contours_cy.pyx - Extension: skimage.measure._moments + Extension: skimage.measure._moments_cy Sources: - skimage/measure/_moments.pyx + skimage/measure/_moments_cy.pyx Extension: skimage.measure._marching_cubes_cy Sources: skimage/measure/_marching_cubes_cy.pyx diff --git a/skimage/measure/_moments.pyx b/skimage/measure/_moments.py similarity index 67% rename from skimage/measure/_moments.pyx rename to skimage/measure/_moments.py index 5a9bc11a..a1816acb 100644 --- a/skimage/measure/_moments.pyx +++ b/skimage/measure/_moments.py @@ -1,11 +1,8 @@ -#cython: cdivision=True -#cython: boundscheck=False -#cython: nonecheck=False -#cython: wraparound=False import numpy as np +from . import _moments_cy -cpdef moments(double[:, :] image, Py_ssize_t order=3): +def moments(image, order=3): """Calculate all raw image moments up to a certain order. The following properties can be calculated from raw image moments: @@ -17,7 +14,7 @@ cpdef moments(double[:, :] image, Py_ssize_t order=3): Parameters ---------- - image : 2D double array + image : 2D double or uint8 array Rasterized shape as image. order : int, optional Maximum order of moments. Default is 3. @@ -39,11 +36,10 @@ cpdef moments(double[:, :] image, Py_ssize_t order=3): .. [4] http://en.wikipedia.org/wiki/Image_moment """ - return moments_central(image, 0, 0, order) + return _moments_cy.moments_central(image, 0, 0, order) -cpdef moments_central(double[:, :] image, double cr, double cc, - Py_ssize_t order=3): +def moments_central(image, cr, cc, order=3): """Calculate all central image moments up to a certain order. Note that central moments are translation invariant but not scale and @@ -51,7 +47,7 @@ cpdef moments_central(double[:, :] image, double cr, double cc, Parameters ---------- - image : 2D double array + image : 2D double or uint8 array Rasterized shape as image. cr : double Center row coordinate. @@ -77,25 +73,11 @@ cpdef moments_central(double[:, :] image, double cr, double cc, .. [4] http://en.wikipedia.org/wiki/Image_moment """ - cdef Py_ssize_t p, q, r, c - cdef double[:, ::1] mu = np.zeros((order + 1, order + 1), dtype=np.double) - cdef double val, dr, dc, dcp, drq - for r in range(image.shape[0]): - dr = r - cr - for c in range(image.shape[1]): - dc = c - cc - val = image[r, c] - dcp = 1 - for p in range(order + 1): - drq = 1 - for q in range(order + 1): - mu[p, q] += val * drq * dcp - drq *= dr - dcp *= dc - return np.asarray(mu) + + return _moments_cy.moments_central(image, cr, cc, order) -def moments_normalized(double[:, :] mu, Py_ssize_t order=3): +def moments_normalized(mu, order=3): """Calculate all normalized central image moments up to a certain order. Note that normalized central moments are translation and scale invariant @@ -125,18 +107,10 @@ def moments_normalized(double[:, :] mu, Py_ssize_t order=3): .. [4] http://en.wikipedia.org/wiki/Image_moment """ - cdef Py_ssize_t p, q - cdef double[:, ::1] nu = np.zeros((order + 1, order + 1), dtype=np.double) - for p in range(order + 1): - for q in range(order + 1): - if p + q >= 2: - nu[p, q] = mu[p, q] / mu[0, 0] ** ((p + q) / 2 + 1) - else: - nu[p, q] = np.nan - return np.asarray(nu) + return _moments_cy.moments_normalized(mu.astype(np.double), order) -def moments_hu(double[:, :] nu): +def moments_hu(nu): """Calculate Hu's set of image moments. Note that this set of moments is proofed to be translation, scale and @@ -167,23 +141,4 @@ def moments_hu(double[:, :] nu): """ - cdef double[::1] hu = np.zeros((7, ), dtype=np.double) - cdef double t0 = nu[3, 0] + nu[1, 2] - cdef double t1 = nu[2, 1] + nu[0, 3] - cdef double q0 = t0 * t0 - cdef double q1 = t1 * t1 - cdef double n4 = 4 * nu[1, 1] - cdef double s = nu[2, 0] + nu[0, 2] - cdef double d = nu[2, 0] - nu[0, 2] - hu[0] = s - hu[1] = d * d + n4 * nu[1, 1] - hu[3] = q0 + q1 - hu[5] = d * (q0 - q1) + n4 * t0 * t1 - t0 *= q0 - 3 * q1 - t1 *= 3 * q0 - q1 - q0 = nu[3, 0]- 3 * nu[1, 2] - q1 = 3 * nu[2, 1] - nu[0, 3] - hu[2] = q0 * q0 + q1 * q1 - hu[4] = q0 * t0 + q1 * t1 - hu[6] = q1 * t0 - q0 * t1 - return np.asarray(hu) + return _moments_cy.moments_hu(nu.astype(np.double)) diff --git a/skimage/measure/_moments_cy.pyx b/skimage/measure/_moments_cy.pyx new file mode 100644 index 00000000..a1d0ebde --- /dev/null +++ b/skimage/measure/_moments_cy.pyx @@ -0,0 +1,65 @@ +#cython: cdivision=True +#cython: boundscheck=False +#cython: nonecheck=False +#cython: wraparound=False +import numpy as np +import cython + + +ctypedef fused image_t: + cython.uchar[:, :] + cython.double[:, :] + + +cpdef moments_central(image_t image, double cr, double cc, Py_ssize_t order): + cdef Py_ssize_t p, q, r, c + cdef double[:, ::1] mu = np.zeros((order + 1, order + 1), dtype=np.double) + cdef double val, dr, dc, dcp, drq + for r in range(image.shape[0]): + dr = r - cr + for c in range(image.shape[1]): + dc = c - cc + val = image[r, c] + dcp = 1 + for p in range(order + 1): + drq = 1 + for q in range(order + 1): + mu[p, q] += val * drq * dcp + drq *= dr + dcp *= dc + return np.asarray(mu) + + +def moments_normalized(double[:, :] mu, Py_ssize_t order=3): + cdef Py_ssize_t p, q + cdef double[:, ::1] nu = np.zeros((order + 1, order + 1), dtype=np.double) + for p in range(order + 1): + for q in range(order + 1): + if p + q >= 2: + nu[p, q] = mu[p, q] / mu[0, 0] ** ((p + q) / 2 + 1) + else: + nu[p, q] = np.nan + return np.asarray(nu) + + +def moments_hu(double[:, :] nu): + cdef double[::1] hu = np.zeros((7, ), dtype=np.double) + cdef double t0 = nu[3, 0] + nu[1, 2] + cdef double t1 = nu[2, 1] + nu[0, 3] + cdef double q0 = t0 * t0 + cdef double q1 = t1 * t1 + cdef double n4 = 4 * nu[1, 1] + cdef double s = nu[2, 0] + nu[0, 2] + cdef double d = nu[2, 0] - nu[0, 2] + hu[0] = s + hu[1] = d * d + n4 * nu[1, 1] + hu[3] = q0 + q1 + hu[5] = d * (q0 - q1) + n4 * t0 * t1 + t0 *= q0 - 3 * q1 + t1 *= 3 * q0 - q1 + q0 = nu[3, 0]- 3 * nu[1, 2] + q1 = 3 * nu[2, 1] - nu[0, 3] + hu[2] = q0 * q0 + q1 * q1 + hu[4] = q0 * t0 + q1 * t1 + hu[6] = q1 * t0 - q0 * t1 + return np.asarray(hu) diff --git a/skimage/measure/setup.py b/skimage/measure/setup.py index 9fce3cf0..ef92450a 100644 --- a/skimage/measure/setup.py +++ b/skimage/measure/setup.py @@ -14,7 +14,7 @@ def configuration(parent_package='', top_path=None): cython(['_ccomp.pyx'], working_path=base_path) cython(['_find_contours_cy.pyx'], working_path=base_path) - cython(['_moments.pyx'], working_path=base_path) + cython(['_moments_cy.pyx'], working_path=base_path) cython(['_marching_cubes_cy.pyx'], working_path=base_path) cython(['_pnpoly.pyx'], working_path=base_path) @@ -22,7 +22,7 @@ def configuration(parent_package='', top_path=None): include_dirs=[get_numpy_include_dirs()]) config.add_extension('_find_contours_cy', sources=['_find_contours_cy.c'], include_dirs=[get_numpy_include_dirs()]) - config.add_extension('_moments', sources=['_moments.c'], + config.add_extension('_moments_cy', sources=['_moments_cy.c'], include_dirs=[get_numpy_include_dirs()]) config.add_extension('_marching_cubes_cy', sources=['_marching_cubes_cy.c'],