From fbf40b0ae7d55ad2479ccecc61e657bab213cbb4 Mon Sep 17 00:00:00 2001 From: Vighnesh Birodkar Date: Sat, 15 Mar 2014 13:43:10 +0530 Subject: [PATCH] corrections to cython file, switched dtype to double --- skimage/feature/_hessian_det_appx.pyx | 50 +++++++++++++++------------ skimage/feature/blob.py | 10 +++--- skimage/feature/tests/test_blob.py | 2 +- 3 files changed, 34 insertions(+), 28 deletions(-) diff --git a/skimage/feature/_hessian_det_appx.pyx b/skimage/feature/_hessian_det_appx.pyx index f620c35b..358a4973 100644 --- a/skimage/feature/_hessian_det_appx.pyx +++ b/skimage/feature/_hessian_det_appx.pyx @@ -1,10 +1,14 @@ +# cython: cdivision=True +# cython: boundscheck=False +# cython: nonecheck=False +# cython: wraparound=False import numpy as np -cimport numpy as np +cimport numpy as cnp from skimage.transform import integral_image, integrate from skimage import util -cdef inline int _clip(np.int_t x, np.int_t low, np.int_t high): +cdef inline Py_ssize_t _clip(Py_ssize_t x, Py_ssize_t low, Py_ssize_t high): """Clips coordinate between high and low. This method was created so that `hessian_det_appx` does not have to make @@ -33,8 +37,9 @@ cdef inline int _clip(np.int_t x, np.int_t low, np.int_t high): return x -cdef inline int _integ(np.int_t[:, :] img, np.int_t r, np.int_t c, - np.int_t rl, np.int_t cl): +cdef inline cnp.double_t _integ( + cnp.double_t[:, :] img, Py_ssize_t r, Py_ssize_t c, + Py_ssize_t rl, Py_ssize_t cl): """Integrate over the integral image in the given window This method was created so that `hessian_det_appx` does not have to make @@ -66,14 +71,14 @@ cdef inline int _integ(np.int_t[:, :] img, np.int_t r, np.int_t c, r2 = _clip(r + rl, 0, img.shape[0] - 1) c2 = _clip(c + cl, 0, img.shape[1] - 1) - cdef np.int_t ans = img[r, c] + img[r2, c2] - img[r, c2] - img[r2, c] + cdef cnp.double_t ans = img[r, c] + img[r2, c2] - img[r, c2] - img[r2, c] if (ans < 0): return 0 return ans -def _hessian_det_appx(np.ndarray[np.int_t, ndim=2] image, float sigma): +def _hessian_det_appx(cnp.double_t[:, :] img, float sigma): """Computes the approximate Hessian Determinant over an image. This method uses box filters over integral images to compute the @@ -81,7 +86,7 @@ def _hessian_det_appx(np.ndarray[np.int_t, ndim=2] image, float sigma): Parameters ---------- - image : array + img : array The integral image over which to compute Hessian Determinant. sigma : float Standard deviation used for the Gaussian kernel, used for the Hessian @@ -107,19 +112,18 @@ def _hessian_det_appx(np.ndarray[np.int_t, ndim=2] image, float sigma): determinant. """ - cdef np.int_t[:, :] img = image - cdef int size = int(3 * sigma) - cdef int height = img.shape[0] - cdef int width = img.shape[1] - cdef int r, c - cdef int s2 = (size - 1) / 2 - cdef int s3 = size / 3 - cdef int l = size/3 - cdef int w = size - cdef int b = (size - 1)/2 - cdef int mid, side - zeros = np.zeros_like(img) - cdef np.ndarray[np.float_t, ndim = 2] out = zeros.astype(np.float) + cdef Py_ssize_t size = int(3 * sigma) + cdef Py_ssize_t height = img.shape[0] + cdef Py_ssize_t width = img.shape[1] + cdef Py_ssize_t r, c + cdef Py_ssize_t s2 = (size - 1) / 2 + cdef Py_ssize_t s3 = size / 3 + cdef Py_ssize_t l = size / 3 + cdef Py_ssize_t w = size + cdef Py_ssize_t b = (size - 1) / 2 + cdef cnp.double_t mid, side, tl, tr, bl, br + cdef cnp.double_t[:, ::1] out = np.zeros_like(img, dtype=np.double) + cdef cnp.double_t w_i = 1.0 / size / size cdef float dxx, dyy, dxy @@ -135,19 +139,19 @@ def _hessian_det_appx(np.ndarray[np.int_t, ndim=2] image, float sigma): tr = _integ(img, r + 1, c - s3, s3, s3) # top right dxy = bl + tr - tl - br - dxy = -dxy / w / w + dxy = -dxy * w_i mid = _integ(img, r - s3 + 1, c - s2, 2 * s3 - 1, w) # middle box side = _integ(img, r - s3 + 1, c - s3 / 2, 2 * s3 - 1, s3) # sides dxx = mid - 3 * side - dxx = -dxx / w / w + dxx = -dxx * w_i mid = _integ(img, r - s2, c - s3 + 1, w, 2 * s3 - 1) side = _integ(img, r - s3 / 2, c - s3 + 1, s3, 2 * s3 - 1) dyy = mid - 3 * side - dyy = -dyy / w / w + dyy = -dyy * w_i out[r, c] = (dxx * dyy - 0.81 * (dxy * dxy)) diff --git a/skimage/feature/blob.py b/skimage/feature/blob.py index 3d602b3a..7c6a5649 100644 --- a/skimage/feature/blob.py +++ b/skimage/feature/blob.py @@ -4,7 +4,7 @@ import itertools as itt import math from math import sqrt, hypot, log from numpy import arccos -from skimage.util import img_as_float, img_as_ubyte +from skimage.util import img_as_float from .peak import peak_local_max from ._hessian_det_appx import _hessian_det_appx from skimage.transform import integral_image @@ -302,7 +302,7 @@ def blob_log(image, min_sigma=1, max_sigma=50, num_sigma=10, threshold=.2, return _prune_blobs(local_maxima, overlap) -def blob_doh(image, min_sigma=1, max_sigma=30, num_sigma=10, threshold=500, +def blob_doh(image, min_sigma=1, max_sigma=30, num_sigma=10, threshold=0.01, overlap=.5, log_scale=False): """Finds blobs in the given grayscale image. @@ -367,12 +367,14 @@ def blob_doh(image, min_sigma=1, max_sigma=30, num_sigma=10, threshold=500, [192, 212, 23], [193, 275, 23], [195, 100, 23], + [197, 44, 20], [197, 153, 20], [260, 173, 30], [262, 243, 23], [265, 113, 23], [270, 363, 30]]) + Notes ----- The radius of each blob is approximately `sigma`. @@ -387,8 +389,8 @@ def blob_doh(image, min_sigma=1, max_sigma=30, num_sigma=10, threshold=500, if image.ndim != 2: raise ValueError("'image' must be a grayscale ") - image = img_as_ubyte(image).astype(np.uint8) - image = integral_image(image).astype(np.int) + image = img_as_float(image) + image = integral_image(image) if log_scale: start, stop = log(min_sigma, 10), log(max_sigma, 10) diff --git a/skimage/feature/tests/test_blob.py b/skimage/feature/tests/test_blob.py index 4166ff2a..10167214 100644 --- a/skimage/feature/tests/test_blob.py +++ b/skimage/feature/tests/test_blob.py @@ -101,7 +101,7 @@ def test_blob_doh(): min_sigma=1, max_sigma=60, num_sigma=10, - threshold=1000) + threshold=.05) radius = lambda x: x[2] s = sorted(blobs, key=radius)