From f87db0a1ecf53faf038b421ad9c2a69b1e240a08 Mon Sep 17 00:00:00 2001 From: Geoffrey French Date: Sun, 31 Aug 2014 15:12:01 +0100 Subject: [PATCH 01/13] Fixes a bug in _update_doc in skimage/io/__init__.py that will attempt to compute max of empty list if no plugins are found. --- skimage/io/__init__.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/skimage/io/__init__.py b/skimage/io/__init__.py index 144fef57..8f3b9bca 100644 --- a/skimage/io/__init__.py +++ b/skimage/io/__init__.py @@ -40,7 +40,8 @@ def _update_doc(doc): info_table = [(p, plugin_info(p).get('description', 'no description')) for p in available_plugins if not p == 'test'] - name_length = max([len(n) for (n, _) in info_table]) + name_length = max([len(n) for (n, _) in info_table]) if len(info_table) > 0 else 0 + description_length = WRAP_LEN - 1 - name_length column_lengths = [name_length, description_length] _format_plugin_info_table(info_table, column_lengths) From d63d89497bd9a6e3b924c3db155be4430347a584 Mon Sep 17 00:00:00 2001 From: Geoffrey French Date: Sun, 31 Aug 2014 16:04:38 +0100 Subject: [PATCH 02/13] Modified rank filters package so that the _core function in core.pyx outputs to a 3D image, permitting the generation of images with arbitrary size feature vector pixels. Implemented windowed_histogram that generates a windowed histogram of an image. --- skimage/filter/rank/__init__.py | 5 +- skimage/filter/rank/_percentile.py | 2 +- skimage/filter/rank/bilateral.py | 2 +- skimage/filter/rank/bilateral_cy.pyx | 55 ++-- skimage/filter/rank/core_cy.pxd | 8 +- skimage/filter/rank/core_cy.pyx | 28 +- skimage/filter/rank/generic.py | 99 +++++-- skimage/filter/rank/generic_cy.pyx | 377 ++++++++++++++------------ skimage/filter/rank/percentile_cy.pyx | 176 ++++++------ 9 files changed, 431 insertions(+), 321 deletions(-) diff --git a/skimage/filter/rank/__init__.py b/skimage/filter/rank/__init__.py index 8641b984..22b0e881 100644 --- a/skimage/filter/rank/__init__.py +++ b/skimage/filter/rank/__init__.py @@ -1,6 +1,6 @@ from .generic import (autolevel, bottomhat, equalize, gradient, maximum, mean, subtract_mean, median, minimum, modal, enhance_contrast, - pop, threshold, tophat, noise_filter, entropy, otsu, sum) + pop, threshold, tophat, noise_filter, entropy, otsu, sum, windowed_histogram) from ._percentile import (autolevel_percentile, gradient_percentile, mean_percentile, subtract_mean_percentile, enhance_contrast_percentile, percentile, @@ -37,4 +37,5 @@ __all__ = ['autolevel', 'noise_filter', 'entropy', 'otsu', - 'percentile'] + 'percentile', + 'windowed_histogram'] diff --git a/skimage/filter/rank/_percentile.py b/skimage/filter/rank/_percentile.py index 01dd0b49..9b93195b 100644 --- a/skimage/filter/rank/_percentile.py +++ b/skimage/filter/rank/_percentile.py @@ -43,7 +43,7 @@ def _apply(func, image, selem, out, mask, shift_x, shift_y, p0, p1, func(image, selem, shift_x=shift_x, shift_y=shift_y, mask=mask, out=out, max_bin=max_bin, p0=p0, p1=p1) - return out + return out.reshape(out.shape[:2]) def autolevel_percentile(image, selem, out=None, mask=None, shift_x=False, diff --git a/skimage/filter/rank/bilateral.py b/skimage/filter/rank/bilateral.py index d01680db..73147c7f 100644 --- a/skimage/filter/rank/bilateral.py +++ b/skimage/filter/rank/bilateral.py @@ -42,7 +42,7 @@ def _apply(func, image, selem, out, mask, shift_x, shift_y, s0, s1, func(image, selem, shift_x=shift_x, shift_y=shift_y, mask=mask, out=out, max_bin=max_bin, s0=s0, s1=s1) - return out + return out.reshape(out.shape[:2]) def mean_bilateral(image, selem, out=None, mask=None, shift_x=False, diff --git a/skimage/filter/rank/bilateral_cy.pyx b/skimage/filter/rank/bilateral_cy.pyx index 25a81766..0e232465 100644 --- a/skimage/filter/rank/bilateral_cy.pyx +++ b/skimage/filter/rank/bilateral_cy.pyx @@ -9,10 +9,11 @@ from libc.math cimport log from .core_cy cimport dtype_t, dtype_t_out, _core -cdef inline double _kernel_mean(Py_ssize_t* histo, double pop, dtype_t g, - Py_ssize_t max_bin, Py_ssize_t mid_bin, - double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): +cdef inline void _kernel_mean(dtype_t_out[:] out, Py_ssize_t* histo, + double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): cdef Py_ssize_t i cdef Py_ssize_t bilat_pop = 0 @@ -24,17 +25,18 @@ cdef inline double _kernel_mean(Py_ssize_t* histo, double pop, dtype_t g, bilat_pop += histo[i] mean += histo[i] * i if bilat_pop: - return mean / bilat_pop + out[0] = mean / bilat_pop else: - return 0 + out[0] = 0 else: - return 0 + out[0] = 0 -cdef inline double _kernel_pop(Py_ssize_t* histo, double pop, dtype_t g, - Py_ssize_t max_bin, Py_ssize_t mid_bin, - double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): +cdef inline void _kernel_pop(dtype_t_out[:] out, Py_ssize_t* histo, + double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): cdef Py_ssize_t i cdef Py_ssize_t bilat_pop = 0 @@ -43,14 +45,15 @@ cdef inline double _kernel_pop(Py_ssize_t* histo, double pop, dtype_t g, for i in range(max_bin): if (g > (i - s0)) and (g < (i + s1)): bilat_pop += histo[i] - return bilat_pop + out[0] = bilat_pop else: - return 0 + out[0] = 0 -cdef inline double _kernel_sum(Py_ssize_t* histo, double pop, dtype_t g, - Py_ssize_t max_bin, Py_ssize_t mid_bin, - double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): +cdef inline void _kernel_sum(dtype_t_out[:] out, Py_ssize_t* histo, + double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): cdef Py_ssize_t i cdef Py_ssize_t bilat_pop = 0 @@ -62,40 +65,40 @@ cdef inline double _kernel_sum(Py_ssize_t* histo, double pop, dtype_t g, bilat_pop += histo[i] sum += histo[i] * i if bilat_pop: - return sum + out[0] = sum else: - return 0 + out[0] = 0 else: - return 0 + out[0] = 0 def _mean(dtype_t[:, ::1] image, char[:, ::1] selem, char[:, ::1] mask, - dtype_t_out[:, ::1] out, + dtype_t_out[:, :, ::1] out, char shift_x, char shift_y, Py_ssize_t s0, Py_ssize_t s1, Py_ssize_t max_bin): - _core(_kernel_mean[dtype_t], image, selem, mask, out, + _core(_kernel_mean[dtype_t_out, dtype_t], image, selem, mask, out, shift_x, shift_y, 0, 0, s0, s1, max_bin) def _pop(dtype_t[:, ::1] image, char[:, ::1] selem, char[:, ::1] mask, - dtype_t_out[:, ::1] out, + dtype_t_out[:, :, ::1] out, char shift_x, char shift_y, Py_ssize_t s0, Py_ssize_t s1, Py_ssize_t max_bin): - _core(_kernel_pop[dtype_t], image, selem, mask, out, + _core(_kernel_pop[dtype_t_out, dtype_t], image, selem, mask, out, shift_x, shift_y, 0, 0, s0, s1, max_bin) def _sum(dtype_t[:, ::1] image, char[:, ::1] selem, char[:, ::1] mask, - dtype_t_out[:, ::1] out, + dtype_t_out[:, :, ::1] out, char shift_x, char shift_y, Py_ssize_t s0, Py_ssize_t s1, Py_ssize_t max_bin): - _core(_kernel_sum[dtype_t], image, selem, mask, out, + _core(_kernel_sum[dtype_t_out, dtype_t], image, selem, mask, out, shift_x, shift_y, 0, 0, s0, s1, max_bin) diff --git a/skimage/filter/rank/core_cy.pxd b/skimage/filter/rank/core_cy.pxd index 2e97e50a..7d23f12a 100644 --- a/skimage/filter/rank/core_cy.pxd +++ b/skimage/filter/rank/core_cy.pxd @@ -15,13 +15,13 @@ cdef dtype_t _max(dtype_t a, dtype_t b) cdef dtype_t _min(dtype_t a, dtype_t b) -cdef void _core(double kernel(Py_ssize_t*, double, dtype_t, - Py_ssize_t, Py_ssize_t, double, - double, Py_ssize_t, Py_ssize_t), +cdef void _core(void kernel(dtype_t_out[:] out, Py_ssize_t*, double, + dtype_t, Py_ssize_t, Py_ssize_t, double, + double, Py_ssize_t, Py_ssize_t), dtype_t[:, ::1] image, char[:, ::1] selem, char[:, ::1] mask, - dtype_t_out[:, ::1] out, + dtype_t_out[:, :, ::1] out, char shift_x, char shift_y, double p0, double p1, Py_ssize_t s0, Py_ssize_t s1, diff --git a/skimage/filter/rank/core_cy.pyx b/skimage/filter/rank/core_cy.pyx index 02c2c8d0..3c8c9e42 100644 --- a/skimage/filter/rank/core_cy.pyx +++ b/skimage/filter/rank/core_cy.pyx @@ -42,13 +42,13 @@ cdef inline char is_in_mask(Py_ssize_t rows, Py_ssize_t cols, return 0 -cdef void _core(double kernel(Py_ssize_t*, double, dtype_t, - Py_ssize_t, Py_ssize_t, double, - double, Py_ssize_t, Py_ssize_t), +cdef void _core(void kernel(dtype_t_out[:] out_feat, Py_ssize_t*, double, dtype_t, + Py_ssize_t, Py_ssize_t, double, + double, Py_ssize_t, Py_ssize_t), dtype_t[:, ::1] image, char[:, ::1] selem, char[:, ::1] mask, - dtype_t_out[:, ::1] out, + dtype_t_out[:, :, ::1] out, char shift_x, char shift_y, double p0, double p1, Py_ssize_t s0, Py_ssize_t s1, @@ -151,8 +151,8 @@ cdef void _core(double kernel(Py_ssize_t*, double, dtype_t, r = 0 c = 0 - out[r, c] = kernel(histo, pop, image[r, c], max_bin, mid_bin, - p0, p1, s0, s1) + kernel(out[r, c, :], histo, pop, image[r, c], max_bin, mid_bin, + p0, p1, s0, s1) # main loop r = 0 @@ -172,8 +172,8 @@ cdef void _core(double kernel(Py_ssize_t*, double, dtype_t, if is_in_mask(rows, cols, rr, cc, mask_data): histogram_decrement(histo, &pop, image[rr, cc]) - out[r, c] = kernel(histo, pop, image[r, c], - max_bin, mid_bin, p0, p1, s0, s1) + kernel(out[r, c, :], histo, pop, image[r, c], max_bin, mid_bin, + p0, p1, s0, s1) r += 1 # pass to the next row if r >= rows: @@ -192,8 +192,8 @@ cdef void _core(double kernel(Py_ssize_t*, double, dtype_t, if is_in_mask(rows, cols, rr, cc, mask_data): histogram_decrement(histo, &pop, image[rr, cc]) - out[r, c] = kernel(histo, pop, image[r, c], - max_bin, mid_bin, p0, p1, s0, s1) + kernel(out[r, c, :], histo, pop, image[r, c], max_bin, mid_bin, + p0, p1, s0, s1) # ---> east to west for c in range(cols - 2, -1, -1): @@ -209,8 +209,8 @@ cdef void _core(double kernel(Py_ssize_t*, double, dtype_t, if is_in_mask(rows, cols, rr, cc, mask_data): histogram_decrement(histo, &pop, image[rr, cc]) - out[r, c] = kernel(histo, pop, image[r, c], - max_bin, mid_bin, p0, p1, s0, s1) + kernel(out[r, c, :], histo, pop, image[r, c], max_bin, mid_bin, + p0, p1, s0, s1) r += 1 # pass to the next row if r >= rows: @@ -229,8 +229,8 @@ cdef void _core(double kernel(Py_ssize_t*, double, dtype_t, if is_in_mask(rows, cols, rr, cc, mask_data): histogram_decrement(histo, &pop, image[rr, cc]) - out[r, c] = kernel(histo, pop, image[r, c], - max_bin, mid_bin, p0, p1, s0, s1) + kernel(out[r, c, :], histo, pop, image[r, c], + max_bin, mid_bin, p0, p1, s0, s1) # release memory allocated by malloc free(se_e_r) diff --git a/skimage/filter/rank/generic.py b/skimage/filter/rank/generic.py index ccc1166e..07c592af 100644 --- a/skimage/filter/rank/generic.py +++ b/skimage/filter/rank/generic.py @@ -28,7 +28,7 @@ __all__ = ['autolevel', 'bottomhat', 'equalize', 'gradient', 'maximum', 'mean', 'pop', 'threshold', 'tophat', 'noise_filter', 'entropy', 'otsu'] -def _handle_input(image, selem, out, mask, out_dtype=None): +def _handle_input(image, selem, out, mask, out_dtype=None, pixel_size=1): if image.dtype not in (np.uint8, np.uint16): image = img_as_ubyte(image) @@ -45,7 +45,10 @@ def _handle_input(image, selem, out, mask, out_dtype=None): if out is None: if out_dtype is None: out_dtype = image.dtype - out = np.empty_like(image, dtype=out_dtype) + out = np.empty(image.shape+(pixel_size,), dtype=out_dtype) + else: + if len(out.shape) == 2: + out = out.reshape(out.shape+(pixel_size,)) if image is out: raise NotImplementedError("Cannot perform rank operation in place.") @@ -65,7 +68,7 @@ def _handle_input(image, selem, out, mask, out_dtype=None): return image, selem, out, mask, max_bin -def _apply(func, image, selem, out, mask, shift_x, shift_y, out_dtype=None): +def _apply_scalar_per_pixel(func, image, selem, out, mask, shift_x, shift_y, out_dtype=None): image, selem, out, mask, max_bin = _handle_input(image, selem, out, mask, out_dtype) @@ -73,6 +76,17 @@ def _apply(func, image, selem, out, mask, shift_x, shift_y, out_dtype=None): func(image, selem, shift_x=shift_x, shift_y=shift_y, mask=mask, out=out, max_bin=max_bin) + return out.reshape(out.shape[:2]) + + +def _apply_vector_per_pixel(func, image, selem, out, mask, shift_x, shift_y, out_dtype=None, pixel_size=1): + + image, selem, out, mask, max_bin = _handle_input(image, selem, out, mask, + out_dtype, pixel_size=pixel_size) + + func(image, selem, shift_x=shift_x, shift_y=shift_y, mask=mask, + out=out, max_bin=max_bin) + return out @@ -113,7 +127,7 @@ def autolevel(image, selem, out=None, mask=None, shift_x=False, shift_y=False): """ - return _apply(generic_cy._autolevel, image, selem, + return _apply_scalar_per_pixel(generic_cy._autolevel, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y) @@ -154,7 +168,7 @@ def bottomhat(image, selem, out=None, mask=None, shift_x=False, shift_y=False): """ - return _apply(generic_cy._bottomhat, image, selem, + return _apply_scalar_per_pixel(generic_cy._bottomhat, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y) @@ -192,7 +206,7 @@ def equalize(image, selem, out=None, mask=None, shift_x=False, shift_y=False): """ - return _apply(generic_cy._equalize, image, selem, + return _apply_scalar_per_pixel(generic_cy._equalize, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y) @@ -230,7 +244,7 @@ def gradient(image, selem, out=None, mask=None, shift_x=False, shift_y=False): """ - return _apply(generic_cy._gradient, image, selem, + return _apply_scalar_per_pixel(generic_cy._gradient, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y) @@ -277,7 +291,7 @@ def maximum(image, selem, out=None, mask=None, shift_x=False, shift_y=False): """ - return _apply(generic_cy._maximum, image, selem, + return _apply_scalar_per_pixel(generic_cy._maximum, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y) @@ -315,7 +329,7 @@ def mean(image, selem, out=None, mask=None, shift_x=False, shift_y=False): """ - return _apply(generic_cy._mean, image, selem, out=out, + return _apply_scalar_per_pixel(generic_cy._mean, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y) @@ -354,7 +368,7 @@ def subtract_mean(image, selem, out=None, mask=None, shift_x=False, """ - return _apply(generic_cy._subtract_mean, image, selem, + return _apply_scalar_per_pixel(generic_cy._subtract_mean, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y) @@ -392,7 +406,7 @@ def median(image, selem, out=None, mask=None, shift_x=False, shift_y=False): """ - return _apply(generic_cy._median, image, selem, + return _apply_scalar_per_pixel(generic_cy._median, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y) @@ -439,7 +453,7 @@ def minimum(image, selem, out=None, mask=None, shift_x=False, shift_y=False): """ - return _apply(generic_cy._minimum, image, selem, + return _apply_scalar_per_pixel(generic_cy._minimum, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y) @@ -479,7 +493,7 @@ def modal(image, selem, out=None, mask=None, shift_x=False, shift_y=False): """ - return _apply(generic_cy._modal, image, selem, + return _apply_scalar_per_pixel(generic_cy._modal, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y) @@ -522,7 +536,7 @@ def enhance_contrast(image, selem, out=None, mask=None, shift_x=False, """ - return _apply(generic_cy._enhance_contrast, image, selem, + return _apply_scalar_per_pixel(generic_cy._enhance_contrast, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y) @@ -571,7 +585,7 @@ def pop(image, selem, out=None, mask=None, shift_x=False, shift_y=False): """ - return _apply(generic_cy._pop, image, selem, out=out, + return _apply_scalar_per_pixel(generic_cy._pop, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y) @@ -620,7 +634,7 @@ def sum(image, selem, out=None, mask=None, shift_x=False, shift_y=False): """ - return _apply(generic_cy._sum, image, selem, out=out, + return _apply_scalar_per_pixel(generic_cy._sum, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y) @@ -669,7 +683,7 @@ def threshold(image, selem, out=None, mask=None, shift_x=False, shift_y=False): """ - return _apply(generic_cy._threshold, image, selem, + return _apply_scalar_per_pixel(generic_cy._threshold, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y) @@ -710,7 +724,7 @@ def tophat(image, selem, out=None, mask=None, shift_x=False, shift_y=False): """ - return _apply(generic_cy._tophat, image, selem, + return _apply_scalar_per_pixel(generic_cy._tophat, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y) @@ -761,7 +775,7 @@ def noise_filter(image, selem, out=None, mask=None, shift_x=False, selem_cpy = selem.copy() selem_cpy[centre_r, centre_c] = 0 - return _apply(generic_cy._noise_filter, image, selem_cpy, out=out, + return _apply_scalar_per_pixel(generic_cy._noise_filter, image, selem_cpy, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y) @@ -806,7 +820,7 @@ def entropy(image, selem, out=None, mask=None, shift_x=False, shift_y=False): """ - return _apply(generic_cy._entropy, image, selem, + return _apply_scalar_per_pixel(generic_cy._entropy, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y, out_dtype=np.double) @@ -850,5 +864,48 @@ def otsu(image, selem, out=None, mask=None, shift_x=False, shift_y=False): """ - return _apply(generic_cy._otsu, image, selem, out=out, + return _apply_scalar_per_pixel(generic_cy._otsu, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y) + + +def windowed_histogram(image, selem, out=None, mask=None, shift_x=False, shift_y=False): + """Sliding window histogram + + Parameters + ---------- + image : ndarray + Image array (uint8 array). + selem : 2-D array + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray + If None, a new array will be allocated. + mask : ndarray + Mask array that defines (>0) area of the image included in the local + neighborhood. If None, the complete image is used (default). + shift_x, shift_y : int + Offset added to the structuring element center point. Shift is bounded + to the structuring element sizes (center must be inside the given + structuring element). + + Returns + ------- + out : 3-D array (same dtype as input image) whose extra dimension + Output image. + + References + ---------- + .. [otsu] http://en.wikipedia.org/wiki/Otsu's_method + + Examples + -------- + >>> from skimage import data + >>> from skimage.filter.rank import windowed_histogram + >>> from skimage.morphology import disk + >>> img = data.camera() + >>> hist_img = windowed_histogram(img, disk(5)) + >>> thresh_image = img >= local_otsu + + """ + + return _apply_vector_per_pixel(generic_cy._windowed_hist, image, selem, out=out, + mask=mask, shift_x=shift_x, shift_y=shift_y, pixel_size=image.max()+1) diff --git a/skimage/filter/rank/generic_cy.pyx b/skimage/filter/rank/generic_cy.pyx index 1c26cf53..88766f44 100644 --- a/skimage/filter/rank/generic_cy.pyx +++ b/skimage/filter/rank/generic_cy.pyx @@ -9,10 +9,11 @@ from libc.math cimport log from .core_cy cimport dtype_t, dtype_t_out, _core -cdef inline double _kernel_autolevel(Py_ssize_t* histo, double pop, dtype_t g, - Py_ssize_t max_bin, Py_ssize_t mid_bin, - double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): +cdef inline void _kernel_autolevel(dtype_t_out[:] out, Py_ssize_t* histo, + double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): cdef Py_ssize_t i, imin, imax, delta @@ -27,17 +28,18 @@ cdef inline double _kernel_autolevel(Py_ssize_t* histo, double pop, dtype_t g, break delta = imax - imin if delta > 0: - return (max_bin - 1) * (g - imin) / delta + out[0] = (max_bin - 1) * (g - imin) / delta else: - return 0 + out[0] = 0 else: - return 0 + out[0] = 0 -cdef inline double _kernel_bottomhat(Py_ssize_t* histo, double pop, dtype_t g, - Py_ssize_t max_bin, Py_ssize_t mid_bin, - double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): +cdef inline void _kernel_bottomhat(dtype_t_out[:] out, Py_ssize_t* histo, + double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): cdef Py_ssize_t i @@ -45,15 +47,16 @@ cdef inline double _kernel_bottomhat(Py_ssize_t* histo, double pop, dtype_t g, for i in range(max_bin): if histo[i]: break - return g - i + out[0] = g - i else: - return 0 + out[0] = 0 -cdef inline double _kernel_equalize(Py_ssize_t* histo, double pop, dtype_t g, - Py_ssize_t max_bin, Py_ssize_t mid_bin, - double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): +cdef inline void _kernel_equalize(dtype_t_out[:] out, Py_ssize_t* histo, + double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): cdef Py_ssize_t i cdef Py_ssize_t sum = 0 @@ -63,15 +66,16 @@ cdef inline double _kernel_equalize(Py_ssize_t* histo, double pop, dtype_t g, sum += histo[i] if i >= g: break - return ((max_bin - 1) * sum) / pop + out[0] = (((max_bin - 1) * sum) / pop) else: - return 0 + out[0] = 0 -cdef inline double _kernel_gradient(Py_ssize_t* histo, double pop, dtype_t g, - Py_ssize_t max_bin, Py_ssize_t mid_bin, - double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): +cdef inline void _kernel_gradient(dtype_t_out[:] out, Py_ssize_t* histo, + double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): cdef Py_ssize_t i, imin, imax @@ -84,65 +88,68 @@ cdef inline double _kernel_gradient(Py_ssize_t* histo, double pop, dtype_t g, if histo[i]: imin = i break - return imax - imin + out[0] = (imax - imin) else: - return 0 + out[0] = 0 -cdef inline double _kernel_maximum(Py_ssize_t* histo, double pop, dtype_t g, - Py_ssize_t max_bin, Py_ssize_t mid_bin, - double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): +cdef inline void _kernel_maximum(dtype_t_out[:] out, Py_ssize_t* histo, + double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): cdef Py_ssize_t i if pop: for i in range(max_bin - 1, -1, -1): if histo[i]: - return i + out[0] = i + return else: - return 0 + out[0] = 0 -cdef inline double _kernel_mean(Py_ssize_t* histo, double pop,dtype_t g, +cdef inline void _kernel_mean(dtype_t_out[:] out, Py_ssize_t* histo, + double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): + + cdef Py_ssize_t i + cdef Py_ssize_t mean = 0 + + if pop: + for i in range(max_bin): + mean += histo[i] * i + out[0] = (mean / pop) + else: + out[0] = 0 + + +cdef inline void _kernel_subtract_mean(dtype_t_out[:] out, Py_ssize_t* histo, + double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): + + cdef Py_ssize_t i + cdef Py_ssize_t mean = 0 + + if pop: + for i in range(max_bin): + mean += histo[i] * i + out[0] = ((g - mean / pop) / 2. + 127) + else: + out[0] = 0 + + +cdef inline void _kernel_median(dtype_t_out[:] out, Py_ssize_t* histo, + double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, Py_ssize_t s0, Py_ssize_t s1): - cdef Py_ssize_t i - cdef Py_ssize_t mean = 0 - - if pop: - for i in range(max_bin): - mean += histo[i] * i - return mean / pop - else: - return 0 - - -cdef inline double _kernel_subtract_mean(Py_ssize_t* histo, double pop, - dtype_t g, - Py_ssize_t max_bin, - Py_ssize_t mid_bin, double p0, - double p1, Py_ssize_t s0, - Py_ssize_t s1): - - cdef Py_ssize_t i - cdef Py_ssize_t mean = 0 - - if pop: - for i in range(max_bin): - mean += histo[i] * i - return (g - mean / pop) / 2. + 127 - else: - return 0 - - -cdef inline double _kernel_median(Py_ssize_t* histo, double pop, dtype_t g, - Py_ssize_t max_bin, Py_ssize_t mid_bin, - double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): - cdef Py_ssize_t i cdef double sum = pop / 2.0 @@ -151,30 +158,34 @@ cdef inline double _kernel_median(Py_ssize_t* histo, double pop, dtype_t g, if histo[i]: sum -= histo[i] if sum < 0: - return i + out[0] = i + return else: - return 0 + out[0] = 0 -cdef inline double _kernel_minimum(Py_ssize_t* histo, double pop, dtype_t g, - Py_ssize_t max_bin, Py_ssize_t mid_bin, - double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): +cdef inline void _kernel_minimum(dtype_t_out[:] out, Py_ssize_t* histo, + double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): cdef Py_ssize_t i if pop: for i in range(max_bin): if histo[i]: - return i + out[0] = i + return else: - return 0 + out[0] = 0 -cdef inline double _kernel_modal(Py_ssize_t* histo, double pop, dtype_t g, - Py_ssize_t max_bin, Py_ssize_t mid_bin, - double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): +cdef inline void _kernel_modal(dtype_t_out[:] out, Py_ssize_t* histo, + double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): cdef Py_ssize_t hmax = 0, imax = 0 @@ -183,17 +194,19 @@ cdef inline double _kernel_modal(Py_ssize_t* histo, double pop, dtype_t g, if histo[i] > hmax: hmax = histo[i] imax = i - return imax + out[0] = imax else: - return 0 + out[0] = 0 -cdef inline double _kernel_enhance_contrast(Py_ssize_t* histo, double pop, - dtype_t g, - Py_ssize_t max_bin, - Py_ssize_t mid_bin, double p0, - double p1, Py_ssize_t s0, - Py_ssize_t s1): +cdef inline void _kernel_enhance_contrast(dtype_t_out[:] out, + Py_ssize_t* histo, + double pop, + dtype_t g, + Py_ssize_t max_bin, + Py_ssize_t mid_bin, double p0, + double p1, Py_ssize_t s0, + Py_ssize_t s1): cdef Py_ssize_t i, imin, imax @@ -207,25 +220,27 @@ cdef inline double _kernel_enhance_contrast(Py_ssize_t* histo, double pop, imin = i break if imax - g < g - imin: - return imax + out[0] = imax else: - return imin + out[0] = imin else: - return 0 + out[0] = 0 -cdef inline double _kernel_pop(Py_ssize_t* histo, double pop, dtype_t g, - Py_ssize_t max_bin, Py_ssize_t mid_bin, - double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): +cdef inline void _kernel_pop(dtype_t_out[:] out, Py_ssize_t* histo, + double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): - return pop + out[0] = pop -cdef inline double _kernel_sum(Py_ssize_t* histo, double pop,dtype_t g, - Py_ssize_t max_bin, Py_ssize_t mid_bin, - double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): +cdef inline void _kernel_sum(dtype_t_out[:] out, Py_ssize_t* histo, + double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): cdef Py_ssize_t i cdef Py_ssize_t sum = 0 @@ -233,15 +248,16 @@ cdef inline double _kernel_sum(Py_ssize_t* histo, double pop,dtype_t g, if pop: for i in range(max_bin): sum += histo[i] * i - return sum + out[0] = sum else: - return 0 + out[0] = 0 -cdef inline double _kernel_threshold(Py_ssize_t* histo, double pop, dtype_t g, - Py_ssize_t max_bin, Py_ssize_t mid_bin, - double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): +cdef inline void _kernel_threshold(dtype_t_out[:] out, Py_ssize_t* histo, + double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): cdef Py_ssize_t i cdef Py_ssize_t mean = 0 @@ -249,15 +265,16 @@ cdef inline double _kernel_threshold(Py_ssize_t* histo, double pop, dtype_t g, if pop: for i in range(max_bin): mean += histo[i] * i - return g > (mean / pop) + out[0] = (g > (mean / pop)) else: - return 0 + out[0] = 0 -cdef inline double _kernel_tophat(Py_ssize_t* histo, double pop, dtype_t g, - Py_ssize_t max_bin, Py_ssize_t mid_bin, - double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): +cdef inline void _kernel_tophat(dtype_t_out[:] out, Py_ssize_t* histo, + double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): cdef Py_ssize_t i @@ -265,23 +282,23 @@ cdef inline double _kernel_tophat(Py_ssize_t* histo, double pop, dtype_t g, for i in range(max_bin - 1, -1, -1): if histo[i]: break - return i - g + out[0] = (i - g) else: - return 0 + out[0] = 0 -cdef inline double _kernel_noise_filter(Py_ssize_t* histo, double pop, - dtype_t g, Py_ssize_t max_bin, - Py_ssize_t mid_bin, double p0, - double p1, Py_ssize_t s0, - Py_ssize_t s1): +cdef inline void _kernel_noise_filter(dtype_t_out[:] out, Py_ssize_t* histo, + double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): cdef Py_ssize_t i cdef Py_ssize_t min_i # early stop if at least one pixel of the neighborhood has the same g if histo[g] > 0: - return 0 + out[0] = 0 for i in range(g, -1, -1): if histo[i]: @@ -291,15 +308,16 @@ cdef inline double _kernel_noise_filter(Py_ssize_t* histo, double pop, if histo[i]: break if i - g < min_i: - return i - g + out[0] = (i - g) else: - return min_i + out[0] = min_i -cdef inline double _kernel_entropy(Py_ssize_t* histo, double pop, dtype_t g, - Py_ssize_t max_bin, Py_ssize_t mid_bin, - double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): +cdef inline void _kernel_entropy(dtype_t_out[:] out, Py_ssize_t* histo, + double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): cdef Py_ssize_t i cdef double e, p @@ -309,15 +327,16 @@ cdef inline double _kernel_entropy(Py_ssize_t* histo, double pop, dtype_t g, p = histo[i] / pop if p > 0: e -= p * log(p) / 0.6931471805599453 - return e + out[0] = e else: - return 0 + out[0] = 0 -cdef inline double _kernel_otsu(Py_ssize_t* histo, double pop, dtype_t g, - Py_ssize_t max_bin, Py_ssize_t mid_bin, - double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): +cdef inline void _kernel_otsu(dtype_t_out[:] out, Py_ssize_t* histo, + double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): cdef Py_ssize_t i cdef Py_ssize_t max_i cdef double P, mu1, mu2, q1, new_q1, sigma_b, max_sigma_b @@ -329,7 +348,7 @@ cdef inline double _kernel_otsu(Py_ssize_t* histo, double pop, dtype_t g, mu += histo[i] * i mu = mu / pop else: - return 0 + out[0] = 0 # maximizing the between class variance max_i = 0 @@ -349,183 +368,205 @@ cdef inline double _kernel_otsu(Py_ssize_t* histo, double pop, dtype_t g, max_i = i q1 = new_q1 - return max_i + out[0] = max_i + + +cdef inline void _kernel_win_hist(dtype_t_out[:] out, Py_ssize_t* histo, + double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): + cdef Py_ssize_t i + cdef Py_ssize_t max_i + for i in xrange(out.shape[0]): + out[i] = histo[i] + def _autolevel(dtype_t[:, ::1] image, char[:, ::1] selem, char[:, ::1] mask, - dtype_t_out[:, ::1] out, + dtype_t_out[:, :, ::1] out, char shift_x, char shift_y, Py_ssize_t max_bin): - _core(_kernel_autolevel[dtype_t], image, selem, mask, out, + _core(_kernel_autolevel[dtype_t_out, dtype_t], image, selem, mask, out, shift_x, shift_y, 0, 0, 0, 0, max_bin) def _bottomhat(dtype_t[:, ::1] image, char[:, ::1] selem, char[:, ::1] mask, - dtype_t_out[:, ::1] out, + dtype_t_out[:, :, ::1] out, char shift_x, char shift_y, Py_ssize_t max_bin): - _core(_kernel_bottomhat[dtype_t], image, selem, mask, out, + _core(_kernel_bottomhat[dtype_t_out, dtype_t], image, selem, mask, out, shift_x, shift_y, 0, 0, 0, 0, max_bin) def _equalize(dtype_t[:, ::1] image, char[:, ::1] selem, char[:, ::1] mask, - dtype_t_out[:, ::1] out, + dtype_t_out[:, :, ::1] out, char shift_x, char shift_y, Py_ssize_t max_bin): - _core(_kernel_equalize[dtype_t], image, selem, mask, out, + _core(_kernel_equalize[dtype_t_out, dtype_t], image, selem, mask, out, shift_x, shift_y, 0, 0, 0, 0, max_bin) def _gradient(dtype_t[:, ::1] image, char[:, ::1] selem, char[:, ::1] mask, - dtype_t_out[:, ::1] out, + dtype_t_out[:, :, ::1] out, char shift_x, char shift_y, Py_ssize_t max_bin): - _core(_kernel_gradient[dtype_t], image, selem, mask, out, + _core(_kernel_gradient[dtype_t_out, dtype_t], image, selem, mask, out, shift_x, shift_y, 0, 0, 0, 0, max_bin) def _maximum(dtype_t[:, ::1] image, char[:, ::1] selem, char[:, ::1] mask, - dtype_t_out[:, ::1] out, + dtype_t_out[:, :, ::1] out, char shift_x, char shift_y, Py_ssize_t max_bin): - _core(_kernel_maximum[dtype_t], image, selem, mask, out, + _core(_kernel_maximum[dtype_t_out, dtype_t], image, selem, mask, out, shift_x, shift_y, 0, 0, 0, 0, max_bin) def _mean(dtype_t[:, ::1] image, char[:, ::1] selem, char[:, ::1] mask, - dtype_t_out[:, ::1] out, + dtype_t_out[:, :, ::1] out, char shift_x, char shift_y, Py_ssize_t max_bin): - _core(_kernel_mean[dtype_t], image, selem, mask, out, + _core(_kernel_mean[dtype_t_out, dtype_t], image, selem, mask, out, shift_x, shift_y, 0, 0, 0, 0, max_bin) def _subtract_mean(dtype_t[:, ::1] image, char[:, ::1] selem, char[:, ::1] mask, - dtype_t_out[:, ::1] out, + dtype_t_out[:, :, ::1] out, char shift_x, char shift_y, Py_ssize_t max_bin): - _core(_kernel_subtract_mean[dtype_t], image, selem, mask, + _core(_kernel_subtract_mean[dtype_t_out, dtype_t], image, selem, mask, out, shift_x, shift_y, 0, 0, 0, 0, max_bin) def _median(dtype_t[:, ::1] image, char[:, ::1] selem, char[:, ::1] mask, - dtype_t_out[:, ::1] out, + dtype_t_out[:, :, ::1] out, char shift_x, char shift_y, Py_ssize_t max_bin): - _core(_kernel_median[dtype_t], image, selem, mask, out, + _core(_kernel_median[dtype_t_out, dtype_t], image, selem, mask, out, shift_x, shift_y, 0, 0, 0, 0, max_bin) def _minimum(dtype_t[:, ::1] image, char[:, ::1] selem, char[:, ::1] mask, - dtype_t_out[:, ::1] out, + dtype_t_out[:, :, ::1] out, char shift_x, char shift_y, Py_ssize_t max_bin): - _core(_kernel_minimum[dtype_t], image, selem, mask, out, + _core(_kernel_minimum[dtype_t_out, dtype_t], image, selem, mask, out, shift_x, shift_y, 0, 0, 0, 0, max_bin) def _enhance_contrast(dtype_t[:, ::1] image, char[:, ::1] selem, char[:, ::1] mask, - dtype_t_out[:, ::1] out, + dtype_t_out[:, :, ::1] out, char shift_x, char shift_y, Py_ssize_t max_bin): - _core(_kernel_enhance_contrast[dtype_t], image, selem, mask, + _core(_kernel_enhance_contrast[dtype_t_out, dtype_t], image, selem, mask, out, shift_x, shift_y, 0, 0, 0, 0, max_bin) def _modal(dtype_t[:, ::1] image, char[:, ::1] selem, char[:, ::1] mask, - dtype_t_out[:, ::1] out, + dtype_t_out[:, :, ::1] out, char shift_x, char shift_y, Py_ssize_t max_bin): - _core(_kernel_modal[dtype_t], image, selem, mask, out, + _core(_kernel_modal[dtype_t_out, dtype_t], image, selem, mask, out, shift_x, shift_y, 0, 0, 0, 0, max_bin) def _pop(dtype_t[:, ::1] image, char[:, ::1] selem, char[:, ::1] mask, - dtype_t_out[:, ::1] out, + dtype_t_out[:, :, ::1] out, char shift_x, char shift_y, Py_ssize_t max_bin): - _core(_kernel_pop[dtype_t], image, selem, mask, out, + _core(_kernel_pop[dtype_t_out, dtype_t], image, selem, mask, out, shift_x, shift_y, 0, 0, 0, 0, max_bin) def _sum(dtype_t[:, ::1] image, char[:, ::1] selem, char[:, ::1] mask, - dtype_t_out[:, ::1] out, + dtype_t_out[:, :, ::1] out, char shift_x, char shift_y, Py_ssize_t max_bin): - _core(_kernel_sum[dtype_t], image, selem, mask, + _core(_kernel_sum[dtype_t_out, dtype_t], image, selem, mask, out, shift_x, shift_y, 0, 0, 0, 0, max_bin) def _threshold(dtype_t[:, ::1] image, char[:, ::1] selem, char[:, ::1] mask, - dtype_t_out[:, ::1] out, + dtype_t_out[:, :, ::1] out, char shift_x, char shift_y, Py_ssize_t max_bin): - _core(_kernel_threshold[dtype_t], image, selem, mask, out, + _core(_kernel_threshold[dtype_t_out, dtype_t], image, selem, mask, out, shift_x, shift_y, 0, 0, 0, 0, max_bin) def _tophat(dtype_t[:, ::1] image, char[:, ::1] selem, char[:, ::1] mask, - dtype_t_out[:, ::1] out, + dtype_t_out[:, :, ::1] out, char shift_x, char shift_y, Py_ssize_t max_bin): - _core(_kernel_tophat[dtype_t], image, selem, mask, out, + _core(_kernel_tophat[dtype_t_out, dtype_t], image, selem, mask, out, shift_x, shift_y, 0, 0, 0, 0, max_bin) def _noise_filter(dtype_t[:, ::1] image, char[:, ::1] selem, char[:, ::1] mask, - dtype_t_out[:, ::1] out, + dtype_t_out[:, :, ::1] out, char shift_x, char shift_y, Py_ssize_t max_bin): - _core(_kernel_noise_filter[dtype_t], image, selem, mask, out, + _core(_kernel_noise_filter[dtype_t_out, dtype_t], image, selem, mask, out, shift_x, shift_y, 0, 0, 0, 0, max_bin) def _entropy(dtype_t[:, ::1] image, char[:, ::1] selem, char[:, ::1] mask, - dtype_t_out[:, ::1] out, + dtype_t_out[:, :, ::1] out, char shift_x, char shift_y, Py_ssize_t max_bin): - _core(_kernel_entropy[dtype_t], image, selem, mask, out, + _core(_kernel_entropy[dtype_t_out, dtype_t], image, selem, mask, out, shift_x, shift_y, 0, 0, 0, 0, max_bin) def _otsu(dtype_t[:, ::1] image, char[:, ::1] selem, char[:, ::1] mask, - dtype_t_out[:, ::1] out, + dtype_t_out[:, :, ::1] out, char shift_x, char shift_y, Py_ssize_t max_bin): - _core(_kernel_otsu[dtype_t], image, selem, mask, out, + _core(_kernel_otsu[dtype_t_out, dtype_t], image, selem, mask, out, + shift_x, shift_y, 0, 0, 0, 0, max_bin) + + +def _windowed_hist(dtype_t[:, ::1] image, + char[:, ::1] selem, + char[:, ::1] mask, + dtype_t_out[:, :, ::1] out, + char shift_x, char shift_y, Py_ssize_t max_bin): + + _core(_kernel_win_hist[dtype_t_out, dtype_t], image, selem, mask, out, shift_x, shift_y, 0, 0, 0, 0, max_bin) diff --git a/skimage/filter/rank/percentile_cy.pyx b/skimage/filter/rank/percentile_cy.pyx index 38d04b33..af33546d 100644 --- a/skimage/filter/rank/percentile_cy.pyx +++ b/skimage/filter/rank/percentile_cy.pyx @@ -7,10 +7,11 @@ cimport numpy as cnp from .core_cy cimport dtype_t, dtype_t_out, _core, _min, _max -cdef inline double _kernel_autolevel(Py_ssize_t* histo, double pop, dtype_t g, - Py_ssize_t max_bin, Py_ssize_t mid_bin, - double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): +cdef inline void _kernel_autolevel(dtype_t_out[:] out, Py_ssize_t* histo, + double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): cdef Py_ssize_t i, imin, imax, sum, delta @@ -31,18 +32,19 @@ cdef inline double _kernel_autolevel(Py_ssize_t* histo, double pop, dtype_t g, delta = imax - imin if delta > 0: - return (max_bin - 1) * (_min(_max(imin, g), imax) - - imin) / delta + out[0] = ((max_bin - 1) * (_min(_max(imin, g), imax) + - imin) / delta) else: - return imax - imin + out[0] = (imax - imin) else: - return 0 + out[0] = 0 -cdef inline double _kernel_gradient(Py_ssize_t* histo, double pop, dtype_t g, - Py_ssize_t max_bin, Py_ssize_t mid_bin, - double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): +cdef inline void _kernel_gradient(dtype_t_out[:] out, Py_ssize_t* histo, + double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): cdef Py_ssize_t i, imin, imax, sum, delta @@ -61,15 +63,16 @@ cdef inline double _kernel_gradient(Py_ssize_t* histo, double pop, dtype_t g, imax = i break - return imax - imin + out[0] = (imax - imin) else: - return 0 + out[0] = 0 -cdef inline double _kernel_mean(Py_ssize_t* histo, double pop, dtype_t g, - Py_ssize_t max_bin, Py_ssize_t mid_bin, - double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): +cdef inline void _kernel_mean(dtype_t_out[:] out, Py_ssize_t* histo, + double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): cdef Py_ssize_t i, sum, mean, n @@ -84,16 +87,17 @@ cdef inline double _kernel_mean(Py_ssize_t* histo, double pop, dtype_t g, mean += histo[i] * i if n > 0: - return mean / n + out[0] = (mean / n) else: - return 0 + out[0] = 0 else: - return 0 + out[0] = 0 -cdef inline double _kernel_sum(Py_ssize_t* histo, double pop, dtype_t g, - Py_ssize_t max_bin, Py_ssize_t mid_bin, - double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): +cdef inline void _kernel_sum(dtype_t_out[:] out, Py_ssize_t* histo, + double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): cdef Py_ssize_t i, sum, sum_g, n @@ -108,18 +112,18 @@ cdef inline double _kernel_sum(Py_ssize_t* histo, double pop, dtype_t g, sum_g += histo[i] * i if n > 0: - return sum_g + out[0] = (sum_g) else: - return 0 + out[0] = 0 else: - return 0 + out[0] = 0 -cdef inline double _kernel_subtract_mean(Py_ssize_t* histo, double pop, - dtype_t g, - Py_ssize_t max_bin, - Py_ssize_t mid_bin, double p0, - double p1, Py_ssize_t s0, - Py_ssize_t s1): +cdef inline void _kernel_subtract_mean(dtype_t_out[:] out, Py_ssize_t* histo, + double pop, dtype_t g, + Py_ssize_t max_bin, + Py_ssize_t mid_bin, double p0, + double p1, Py_ssize_t s0, + Py_ssize_t s1): cdef Py_ssize_t i, sum, mean, n @@ -133,19 +137,20 @@ cdef inline double _kernel_subtract_mean(Py_ssize_t* histo, double pop, n += histo[i] mean += histo[i] * i if n > 0: - return (g - (mean / n)) * .5 + mid_bin + out[0] = ((g - (mean / n)) * .5 + mid_bin) else: - return 0 + out[0] = 0 else: - return 0 + out[0] = 0 -cdef inline double _kernel_enhance_contrast(Py_ssize_t* histo, double pop, - dtype_t g, - Py_ssize_t max_bin, - Py_ssize_t mid_bin, double p0, - double p1, Py_ssize_t s0, - Py_ssize_t s1): +cdef inline void _kernel_enhance_contrast(dtype_t_out[:] out, + Py_ssize_t* histo, double pop, + dtype_t g, + Py_ssize_t max_bin, + Py_ssize_t mid_bin, double p0, + double p1, Py_ssize_t s0, + Py_ssize_t s1): cdef Py_ssize_t i, imin, imax, sum, delta @@ -164,21 +169,22 @@ cdef inline double _kernel_enhance_contrast(Py_ssize_t* histo, double pop, imax = i break if g > imax: - return imax + out[0] = imax if g < imin: - return imin + out[0] = imin if imax - g < g - imin: - return imax + out[0] = imax else: - return imin + out[0] = imin else: - return 0 + out[0] = 0 -cdef inline double _kernel_percentile(Py_ssize_t* histo, double pop, dtype_t g, - Py_ssize_t max_bin, Py_ssize_t mid_bin, - double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): +cdef inline void _kernel_percentile(dtype_t_out[:] out, Py_ssize_t* histo, + double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): cdef Py_ssize_t i cdef Py_ssize_t sum = 0 @@ -193,15 +199,16 @@ cdef inline double _kernel_percentile(Py_ssize_t* histo, double pop, dtype_t g, sum += histo[i] if sum > p0 * pop: break - return i + out[0] = i else: - return 0 + out[0] = 0 -cdef inline double _kernel_pop(Py_ssize_t* histo, double pop, dtype_t g, - Py_ssize_t max_bin, Py_ssize_t mid_bin, - double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): +cdef inline void _kernel_pop(dtype_t_out[:] out, Py_ssize_t* histo, + double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): cdef Py_ssize_t i, sum, n @@ -212,15 +219,16 @@ cdef inline double _kernel_pop(Py_ssize_t* histo, double pop, dtype_t g, sum += histo[i] if (sum >= p0 * pop) and (sum <= p1 * pop): n += histo[i] - return n + out[0] = n else: - return 0 + out[0] = 0 -cdef inline double _kernel_threshold(Py_ssize_t* histo, double pop, dtype_t g, - Py_ssize_t max_bin, Py_ssize_t mid_bin, - double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): +cdef inline void _kernel_threshold(dtype_t_out[:] out, Py_ssize_t* histo, + double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): cdef int i cdef Py_ssize_t sum = 0 @@ -231,103 +239,103 @@ cdef inline double _kernel_threshold(Py_ssize_t* histo, double pop, dtype_t g, if sum >= p0 * pop: break - return (max_bin - 1) * (g >= i) + out[0] = ((max_bin - 1) * (g >= i)) else: - return 0 + out[0] = 0 def _autolevel(dtype_t[:, ::1] image, char[:, ::1] selem, char[:, ::1] mask, - dtype_t_out[:, ::1] out, + dtype_t_out[:, :, ::1] out, char shift_x, char shift_y, double p0, double p1, Py_ssize_t max_bin): - _core(_kernel_autolevel[dtype_t], image, selem, mask, out, + _core(_kernel_autolevel[dtype_t_out, dtype_t], image, selem, mask, out, shift_x, shift_y, p0, p1, 0, 0, max_bin) def _gradient(dtype_t[:, ::1] image, char[:, ::1] selem, char[:, ::1] mask, - dtype_t_out[:, ::1] out, + dtype_t_out[:, :, ::1] out, char shift_x, char shift_y, double p0, double p1, Py_ssize_t max_bin): - _core(_kernel_gradient[dtype_t], image, selem, mask, out, + _core(_kernel_gradient[dtype_t_out, dtype_t], image, selem, mask, out, shift_x, shift_y, p0, p1, 0, 0, max_bin) def _mean(dtype_t[:, ::1] image, char[:, ::1] selem, char[:, ::1] mask, - dtype_t_out[:, ::1] out, + dtype_t_out[:, :, ::1] out, char shift_x, char shift_y, double p0, double p1, Py_ssize_t max_bin): - _core(_kernel_mean[dtype_t], image, selem, mask, out, + _core(_kernel_mean[dtype_t_out, dtype_t], image, selem, mask, out, shift_x, shift_y, p0, p1, 0, 0, max_bin) def _sum(dtype_t[:, ::1] image, char[:, ::1] selem, char[:, ::1] mask, - dtype_t_out[:, ::1] out, + dtype_t_out[:, :, ::1] out, char shift_x, char shift_y, double p0, double p1, Py_ssize_t max_bin): - _core(_kernel_sum[dtype_t], image, selem, mask, out, + _core(_kernel_sum[dtype_t_out, dtype_t], image, selem, mask, out, shift_x, shift_y, p0, p1, 0, 0, max_bin) def _subtract_mean(dtype_t[:, ::1] image, char[:, ::1] selem, char[:, ::1] mask, - dtype_t_out[:, ::1] out, + dtype_t_out[:, :, ::1] out, char shift_x, char shift_y, double p0, double p1, Py_ssize_t max_bin): - _core(_kernel_subtract_mean[dtype_t], image, selem, mask, + _core(_kernel_subtract_mean[dtype_t_out, dtype_t], image, selem, mask, out, shift_x, shift_y, p0, p1, 0, 0, max_bin) def _enhance_contrast(dtype_t[:, ::1] image, char[:, ::1] selem, char[:, ::1] mask, - dtype_t_out[:, ::1] out, + dtype_t_out[:, :, ::1] out, char shift_x, char shift_y, double p0, double p1, Py_ssize_t max_bin): - _core(_kernel_enhance_contrast[dtype_t], image, selem, mask, + _core(_kernel_enhance_contrast[dtype_t_out, dtype_t], image, selem, mask, out, shift_x, shift_y, p0, p1, 0, 0, max_bin) def _percentile(dtype_t[:, ::1] image, char[:, ::1] selem, char[:, ::1] mask, - dtype_t_out[:, ::1] out, + dtype_t_out[:, :, ::1] out, char shift_x, char shift_y, double p0, double p1, Py_ssize_t max_bin): - _core(_kernel_percentile[dtype_t], image, selem, mask, out, + _core(_kernel_percentile[dtype_t_out, dtype_t], image, selem, mask, out, shift_x, shift_y, p0, 1, 0, 0, max_bin) def _pop(dtype_t[:, ::1] image, char[:, ::1] selem, char[:, ::1] mask, - dtype_t_out[:, ::1] out, + dtype_t_out[:, :, ::1] out, char shift_x, char shift_y, double p0, double p1, Py_ssize_t max_bin): - _core(_kernel_pop[dtype_t], image, selem, mask, out, + _core(_kernel_pop[dtype_t_out, dtype_t], image, selem, mask, out, shift_x, shift_y, p0, p1, 0, 0, max_bin) def _threshold(dtype_t[:, ::1] image, char[:, ::1] selem, char[:, ::1] mask, - dtype_t_out[:, ::1] out, + dtype_t_out[:, :, ::1] out, char shift_x, char shift_y, double p0, double p1, Py_ssize_t max_bin): - _core(_kernel_threshold[dtype_t], image, selem, mask, out, + _core(_kernel_threshold[dtype_t_out, dtype_t], image, selem, mask, out, shift_x, shift_y, p0, 1, 0, 0, max_bin) From 7ebb2388d29e8fca1a1af4d0ddfea6c4ff00b3ca Mon Sep 17 00:00:00 2001 From: Geoffrey French Date: Sun, 31 Aug 2014 16:22:35 +0100 Subject: [PATCH 03/13] Fixed some test failures and added a test for windowed_histogram. --- skimage/filter/rank/generic.py | 6 +++--- skimage/filter/rank/tests/test_rank.py | 26 ++++++++++++++++++++++++++ 2 files changed, 29 insertions(+), 3 deletions(-) diff --git a/skimage/filter/rank/generic.py b/skimage/filter/rank/generic.py index 07c592af..487d34b5 100644 --- a/skimage/filter/rank/generic.py +++ b/skimage/filter/rank/generic.py @@ -42,6 +42,9 @@ def _handle_input(image, selem, out, mask, out_dtype=None, pixel_size=1): mask = img_as_ubyte(mask) mask = np.ascontiguousarray(mask) + if image is out: + raise NotImplementedError("Cannot perform rank operation in place.") + if out is None: if out_dtype is None: out_dtype = image.dtype @@ -50,9 +53,6 @@ def _handle_input(image, selem, out, mask, out_dtype=None, pixel_size=1): if len(out.shape) == 2: out = out.reshape(out.shape+(pixel_size,)) - if image is out: - raise NotImplementedError("Cannot perform rank operation in place.") - is_8bit = image.dtype in (np.uint8, np.int8) if is_8bit: diff --git a/skimage/filter/rank/tests/test_rank.py b/skimage/filter/rank/tests/test_rank.py index 5cbffd90..390a7790 100644 --- a/skimage/filter/rank/tests/test_rank.py +++ b/skimage/filter/rank/tests/test_rank.py @@ -542,6 +542,32 @@ def test_sum(): rank.sum_bilateral(image=image16, selem=elem, out=out16, mask=mask,s0=1000,s1=1000) assert_array_equal(r, out16) +def test_windowed_histogram(): + # check the number of valid pixels in the neighborhood + + image8 = np.array([[0, 0, 0, 0, 0], + [0, 1, 1, 1, 0], + [0, 1, 1, 1, 0], + [0, 1, 1, 1, 0], + [0, 0, 0, 0, 0]], dtype=np.uint8) + elem = np.ones((3, 3), dtype=np.uint8) + out8 = np.empty(image8.shape+(2,), dtype=np.uint8) + mask = np.ones(image8.shape, dtype=np.uint8) + + r0 = np.array([[3, 4, 3, 4, 3], + [4, 5, 3, 5, 4], + [3, 3, 0, 3, 3], + [4, 5, 3, 5, 4], + [3, 4, 3, 4, 3]], dtype=np.uint8) + r1 = np.array([[1, 2, 3, 2, 1], + [2, 4, 6, 4, 2], + [3, 6, 9, 6, 3], + [2, 4, 6, 4, 2], + [1, 2, 3, 2, 1]], dtype=np.uint8) + rank.windowed_histogram(image=image8, selem=elem, out=out8, mask=mask) + assert_array_equal(r0, out8[:,:,0]) + assert_array_equal(r1, out8[:,:,1]) + if __name__ == "__main__": run_module_suite() From da93619e59d3aac0ec4945ca5d4646beea770009 Mon Sep 17 00:00:00 2001 From: Geoffrey French Date: Sun, 31 Aug 2014 22:25:45 +0100 Subject: [PATCH 04/13] Fixed some issues with the new windowed_histogram function in filter.rank. It used to be able to output uint8 histogram, whose max pixel counts of 255 could easily overflow. Addressed by limiting the output type to float. Normalized histograms are now generated, otherwise behaviour is unpredictable at boundaries or when pixels are not permitted by a mask. The new optional n_bins parameter allows the caller to specify the size of the histogram generated. Having it fixed to image.max()+1 could result in feature vectors being shorter than desired just due to a value not being used in an image. An example has been added to the docs, that demonstrates the application of windows histograms in object matching; a single coin is extracted and found by chi squared histogram matching. --- doc/examples/plot_windowed_histogram.py | 129 ++++++++++++++++++++++++ skimage/filter/rank/generic.py | 31 ++++-- skimage/filter/rank/generic_cy.pyx | 10 +- skimage/filter/rank/tests/test_rank.py | 37 ++++--- 4 files changed, 182 insertions(+), 25 deletions(-) create mode 100644 doc/examples/plot_windowed_histogram.py diff --git a/doc/examples/plot_windowed_histogram.py b/doc/examples/plot_windowed_histogram.py new file mode 100644 index 00000000..7b0bcb3a --- /dev/null +++ b/doc/examples/plot_windowed_histogram.py @@ -0,0 +1,129 @@ +""" +======================== +Sliding window histogram +======================== + +This example extracts a single coin from the coins image and generates a +histogram of its greyscale values. + +It then computes a sliding window histogram of the complete image using +rank.windowed_histogram. The local histogram for the region surrounding +each pixel in the image is compared to that of the single coin, with +a similarity measure being computed and displayed. + +To demonstrate the rotational invariance of the technique, the same +test is performed on a version of the coins image rotated by 45 degrees. +""" +import numpy as np +import matplotlib +import matplotlib.pyplot as plt + +from skimage import data +from skimage.util.dtype import dtype_range +from skimage.util import img_as_ubyte +from skimage import exposure +from skimage.morphology import disk +from skimage.filter import rank +from skimage import transform + + +matplotlib.rcParams['font.size'] = 9 + + +def windowed_histogram_similarity(image, selem, reference_hist, n_bins): + # Compute normalized windowed histogram feature vector for each pixel + px_histograms = rank.windowed_histogram(image, selem, n_bins=n_bins) + + # Reshape coin histogram to (1,1,N) for broadcast when we want to use it in + # arithmetic operations with the windowed histograms fro the image + reference_hist = reference_hist.reshape((1,1) + reference_hist.shape) + + # Compute Chi squared distance metric: sum((X-Y)**2 / (X+Y); + # a measure of distance between histograms + X = px_histograms + Y = reference_hist + num = (X-Y)*(X-Y) + denom = X+Y + frac = num / denom + frac[denom==0] = 0 + chi_sqr = np.sum(frac, axis=2) * 0.5 + + # Generate a similarity measure. It needs to be low when distance is high. + # and high when distance is low; taking the reciprocal will do this. + # Chi squared will always be >= 0. Add small value to prevent divide by 0. + # Square the denominator to push low values toward 0; this makes the + # high similarity regions stand out in the figure created below; this + # us just done for aesthetics. + similarity = 1 / (chi_sqr + 1.0e-6)**2 + + return similarity + + +# Load the coins image +img = img_as_ubyte(data.coins()) +# img = img_as_ubyte(plt.imread('../../skimage/data/coins.png')) + +# Quantize to 16 levels of grayscale; this way the output image will have a +# 16-dimensional feature vector per pixel +quantized_img = img/16 + +# Select the coin from the 4th column, second row. +# Co-ordinate ordering: [x1,y1,x2,y2] +coin_coords = [184,100,228,148] # 44 x 44 region +coin = quantized_img[coin_coords[1]:coin_coords[3], coin_coords[0]:coin_coords[2]] + +# Compute coin histogram and normalize +coin_hist, _ = np.histogram(coin.flatten(), bins=16, range=(0,16)) +coin_hist = coin_hist.astype(float) / np.sum(coin_hist) + + +# Compute a disk shaped mask that will define the shape of our sliding window +# Example coin is ~44px across, so make a disk 61px wide (2*rad+1) to be big +# enough for other coins too. +selem = disk(30) + + +# Compute the similarity across the complete image +similarity = windowed_histogram_similarity(quantized_img, selem, coin_hist, + coin_hist.shape[0]) + +# Now try a rotated image +rotated_img = img_as_ubyte(transform.rotate(img, 45.0, resize=True)) +# Quantize to 16 levels as before +quantized_rotated_image = rotated_img/16 +# Similarity on rotated image +rotated_similarity = windowed_histogram_similarity(quantized_rotated_image, + selem, coin_hist, + coin_hist.shape[0]) + + + +# Plot it all +fig, axes = plt.subplots(nrows=5, figsize=(6, 18)) +ax0, ax1, ax2, ax3, ax4 = axes + +ax0.imshow(img, cmap='gray') +ax0.set_title('Original image') +ax0.axis('off') + +ax1.imshow(quantized_img, cmap='gray') +ax1.set_title('Quantized image') +ax1.axis('off') + +ax2.imshow(coin, cmap='gray') +ax2.set_title('Coin from 2nd row, 4th column') +ax2.axis('off') + +ax3.imshow(img, cmap='gray') +# While jet is not a great colormap, it makes the high similarity areas +# stand out +ax3.imshow(similarity, cmap='jet', alpha=0.5) +ax3.set_title('Original image with overlayed similarity') +ax3.axis('off') + +ax4.imshow(rotated_img, cmap='gray') +ax4.imshow(rotated_similarity, cmap='jet', alpha=0.5) +ax4.set_title('Rotated image with overlayed similarity') +ax4.axis('off') + +plt.show() diff --git a/skimage/filter/rank/generic.py b/skimage/filter/rank/generic.py index 487d34b5..2af2428e 100644 --- a/skimage/filter/rank/generic.py +++ b/skimage/filter/rank/generic.py @@ -868,8 +868,8 @@ def otsu(image, selem, out=None, mask=None, shift_x=False, shift_y=False): mask=mask, shift_x=shift_x, shift_y=shift_y) -def windowed_histogram(image, selem, out=None, mask=None, shift_x=False, shift_y=False): - """Sliding window histogram +def windowed_histogram(image, selem, out=None, mask=None, shift_x=False, shift_y=False, n_bins=None): + """Normalized sliding window histogram Parameters ---------- @@ -886,15 +886,19 @@ def windowed_histogram(image, selem, out=None, mask=None, shift_x=False, shift_y Offset added to the structuring element center point. Shift is bounded to the structuring element sizes (center must be inside the given structuring element). + n_bins : int or None + The number of histogram bins. Will default to image.max() + 1 if None + is passed. Returns ------- - out : 3-D array (same dtype as input image) whose extra dimension - Output image. - - References - ---------- - .. [otsu] http://en.wikipedia.org/wiki/Otsu's_method + out : 3-D array with float dtype of dimensions (H,W,N), where (H,W) are + the dimensions of the input image and N is n_bins or image.max()+1 + if no value is provided as a parameter. Effectively, each pixel + is an N-dimensional feature vector that is the histogram. + The sum of the elements in the feature vector will be 1, unless + no pixels in the window were covered by both selem and mask, in which + case all elements will be 0. Examples -------- @@ -903,9 +907,14 @@ def windowed_histogram(image, selem, out=None, mask=None, shift_x=False, shift_y >>> from skimage.morphology import disk >>> img = data.camera() >>> hist_img = windowed_histogram(img, disk(5)) - >>> thresh_image = img >= local_otsu """ - return _apply_vector_per_pixel(generic_cy._windowed_hist, image, selem, out=out, - mask=mask, shift_x=shift_x, shift_y=shift_y, pixel_size=image.max()+1) + if n_bins is None: + n_bins = image.max() + 1 + + return _apply_vector_per_pixel(generic_cy._windowed_hist, image, selem, + out=out, mask=mask, + shift_x=shift_x, shift_y=shift_y, + out_dtype=np.double, + pixel_size=n_bins) diff --git a/skimage/filter/rank/generic_cy.pyx b/skimage/filter/rank/generic_cy.pyx index 88766f44..bafd2432 100644 --- a/skimage/filter/rank/generic_cy.pyx +++ b/skimage/filter/rank/generic_cy.pyx @@ -378,8 +378,14 @@ cdef inline void _kernel_win_hist(dtype_t_out[:] out, Py_ssize_t* histo, Py_ssize_t s0, Py_ssize_t s1): cdef Py_ssize_t i cdef Py_ssize_t max_i - for i in xrange(out.shape[0]): - out[i] = histo[i] + cdef double scale + if pop: + scale = 1.0 / pop + for i in xrange(out.shape[0]): + out[i] = (histo[i] * scale) + else: + for i in xrange(out.shape[0]): + out[i] = 0 diff --git a/skimage/filter/rank/tests/test_rank.py b/skimage/filter/rank/tests/test_rank.py index 390a7790..48d02462 100644 --- a/skimage/filter/rank/tests/test_rank.py +++ b/skimage/filter/rank/tests/test_rank.py @@ -551,22 +551,35 @@ def test_windowed_histogram(): [0, 1, 1, 1, 0], [0, 0, 0, 0, 0]], dtype=np.uint8) elem = np.ones((3, 3), dtype=np.uint8) - out8 = np.empty(image8.shape+(2,), dtype=np.uint8) + outf = np.empty(image8.shape+(2,), dtype=float) mask = np.ones(image8.shape, dtype=np.uint8) + # Population so we can normalize the expected output while maintaining + # code readability + pop = np.array([[4, 6, 6, 6, 4], + [6, 9, 9, 9, 6], + [6, 9, 9, 9, 6], + [6, 9, 9, 9, 6], + [4, 6, 6, 6, 4]], dtype=float) + r0 = np.array([[3, 4, 3, 4, 3], - [4, 5, 3, 5, 4], - [3, 3, 0, 3, 3], - [4, 5, 3, 5, 4], - [3, 4, 3, 4, 3]], dtype=np.uint8) + [4, 5, 3, 5, 4], + [3, 3, 0, 3, 3], + [4, 5, 3, 5, 4], + [3, 4, 3, 4, 3]], dtype=float) / pop r1 = np.array([[1, 2, 3, 2, 1], - [2, 4, 6, 4, 2], - [3, 6, 9, 6, 3], - [2, 4, 6, 4, 2], - [1, 2, 3, 2, 1]], dtype=np.uint8) - rank.windowed_histogram(image=image8, selem=elem, out=out8, mask=mask) - assert_array_equal(r0, out8[:,:,0]) - assert_array_equal(r1, out8[:,:,1]) + [2, 4, 6, 4, 2], + [3, 6, 9, 6, 3], + [2, 4, 6, 4, 2], + [1, 2, 3, 2, 1]], dtype=float) / pop + rank.windowed_histogram(image=image8, selem=elem, out=outf, mask=mask) + assert_array_equal(r0, outf[:,:,0]) + assert_array_equal(r1, outf[:,:,1]) + + # Test n_bins parameter + larger_output = rank.windowed_histogram(image=image8, selem=elem, + mask=mask, n_bins=5) + assert larger_output.shape[2] == 5 if __name__ == "__main__": From 29763582f55733729e2f46d2fbc1c8c499286096 Mon Sep 17 00:00:00 2001 From: Geoffrey French Date: Mon, 1 Sep 2014 20:29:17 +0100 Subject: [PATCH 05/13] Fixed Python 3 error in plot_windowed_histogram caused by use of incorrect division operator. Also tweaked the similarity computation a little. --- doc/examples/plot_windowed_histogram.py | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/doc/examples/plot_windowed_histogram.py b/doc/examples/plot_windowed_histogram.py index 7b0bcb3a..cfa9345d 100644 --- a/doc/examples/plot_windowed_histogram.py +++ b/doc/examples/plot_windowed_histogram.py @@ -1,3 +1,4 @@ +from __future__ import division """ ======================== Sliding window histogram @@ -50,11 +51,8 @@ def windowed_histogram_similarity(image, selem, reference_hist, n_bins): # Generate a similarity measure. It needs to be low when distance is high. # and high when distance is low; taking the reciprocal will do this. - # Chi squared will always be >= 0. Add small value to prevent divide by 0. - # Square the denominator to push low values toward 0; this makes the - # high similarity regions stand out in the figure created below; this - # us just done for aesthetics. - similarity = 1 / (chi_sqr + 1.0e-6)**2 + # Chi squared will always be >= 0, add small value to prevent divide by 0. + similarity = 1 / (chi_sqr + 1.0e-4) return similarity @@ -65,7 +63,7 @@ img = img_as_ubyte(data.coins()) # Quantize to 16 levels of grayscale; this way the output image will have a # 16-dimensional feature vector per pixel -quantized_img = img/16 +quantized_img = img//16 # Select the coin from the 4th column, second row. # Co-ordinate ordering: [x1,y1,x2,y2] @@ -90,7 +88,7 @@ similarity = windowed_histogram_similarity(quantized_img, selem, coin_hist, # Now try a rotated image rotated_img = img_as_ubyte(transform.rotate(img, 45.0, resize=True)) # Quantize to 16 levels as before -quantized_rotated_image = rotated_img/16 +quantized_rotated_image = rotated_img//16 # Similarity on rotated image rotated_similarity = windowed_histogram_similarity(quantized_rotated_image, selem, coin_hist, From b9ac11da3f20dd30fe3c7bab30124df72c2fb152 Mon Sep 17 00:00:00 2001 From: Geoffrey French Date: Mon, 1 Sep 2014 22:39:41 +0100 Subject: [PATCH 06/13] Exchanged array output parameter for pointer and Py_ssize_t in kernel functions for slight speedup. --- skimage/filter/rank/bilateral_cy.pyx | 9 ++-- skimage/filter/rank/core_cy.pxd | 2 +- skimage/filter/rank/core_cy.pyx | 13 +++--- skimage/filter/rank/generic_cy.pyx | 61 ++++++++++++++++++--------- skimage/filter/rank/percentile_cy.pyx | 27 ++++++++---- 5 files changed, 72 insertions(+), 40 deletions(-) diff --git a/skimage/filter/rank/bilateral_cy.pyx b/skimage/filter/rank/bilateral_cy.pyx index 0e232465..20a97788 100644 --- a/skimage/filter/rank/bilateral_cy.pyx +++ b/skimage/filter/rank/bilateral_cy.pyx @@ -9,7 +9,8 @@ from libc.math cimport log from .core_cy cimport dtype_t, dtype_t_out, _core -cdef inline void _kernel_mean(dtype_t_out[:] out, Py_ssize_t* histo, +cdef inline void _kernel_mean(dtype_t_out* out, Py_ssize_t odepth, + Py_ssize_t* histo, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, @@ -32,7 +33,8 @@ cdef inline void _kernel_mean(dtype_t_out[:] out, Py_ssize_t* histo, out[0] = 0 -cdef inline void _kernel_pop(dtype_t_out[:] out, Py_ssize_t* histo, +cdef inline void _kernel_pop(dtype_t_out* out, Py_ssize_t odepth, + Py_ssize_t* histo, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, @@ -49,7 +51,8 @@ cdef inline void _kernel_pop(dtype_t_out[:] out, Py_ssize_t* histo, else: out[0] = 0 -cdef inline void _kernel_sum(dtype_t_out[:] out, Py_ssize_t* histo, +cdef inline void _kernel_sum(dtype_t_out* out, Py_ssize_t odepth, + Py_ssize_t* histo, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, diff --git a/skimage/filter/rank/core_cy.pxd b/skimage/filter/rank/core_cy.pxd index 7d23f12a..9627c0d9 100644 --- a/skimage/filter/rank/core_cy.pxd +++ b/skimage/filter/rank/core_cy.pxd @@ -15,7 +15,7 @@ cdef dtype_t _max(dtype_t a, dtype_t b) cdef dtype_t _min(dtype_t a, dtype_t b) -cdef void _core(void kernel(dtype_t_out[:] out, Py_ssize_t*, double, +cdef void _core(void kernel(dtype_t_out*, Py_ssize_t, Py_ssize_t*, double, dtype_t, Py_ssize_t, Py_ssize_t, double, double, Py_ssize_t, Py_ssize_t), dtype_t[:, ::1] image, diff --git a/skimage/filter/rank/core_cy.pyx b/skimage/filter/rank/core_cy.pyx index 3c8c9e42..c1216c0d 100644 --- a/skimage/filter/rank/core_cy.pyx +++ b/skimage/filter/rank/core_cy.pyx @@ -42,7 +42,7 @@ cdef inline char is_in_mask(Py_ssize_t rows, Py_ssize_t cols, return 0 -cdef void _core(void kernel(dtype_t_out[:] out_feat, Py_ssize_t*, double, dtype_t, +cdef void _core(void kernel(dtype_t_out*, Py_ssize_t, Py_ssize_t*, double, dtype_t, Py_ssize_t, Py_ssize_t, double, double, Py_ssize_t, Py_ssize_t), dtype_t[:, ::1] image, @@ -61,6 +61,7 @@ cdef void _core(void kernel(dtype_t_out[:] out_feat, Py_ssize_t*, double, dtype_ cdef Py_ssize_t cols = image.shape[1] cdef Py_ssize_t srows = selem.shape[0] cdef Py_ssize_t scols = selem.shape[1] + cdef Py_ssize_t odepth = out.shape[2] cdef Py_ssize_t centre_r = (selem.shape[0] / 2) + shift_y cdef Py_ssize_t centre_c = (selem.shape[1] / 2) + shift_x @@ -151,7 +152,7 @@ cdef void _core(void kernel(dtype_t_out[:] out_feat, Py_ssize_t*, double, dtype_ r = 0 c = 0 - kernel(out[r, c, :], histo, pop, image[r, c], max_bin, mid_bin, + kernel(&out[r, c, 0], odepth, histo, pop, image[r, c], max_bin, mid_bin, p0, p1, s0, s1) # main loop @@ -172,7 +173,7 @@ cdef void _core(void kernel(dtype_t_out[:] out_feat, Py_ssize_t*, double, dtype_ if is_in_mask(rows, cols, rr, cc, mask_data): histogram_decrement(histo, &pop, image[rr, cc]) - kernel(out[r, c, :], histo, pop, image[r, c], max_bin, mid_bin, + kernel(&out[r, c, 0], odepth, histo, pop, image[r, c], max_bin, mid_bin, p0, p1, s0, s1) r += 1 # pass to the next row @@ -192,7 +193,7 @@ cdef void _core(void kernel(dtype_t_out[:] out_feat, Py_ssize_t*, double, dtype_ if is_in_mask(rows, cols, rr, cc, mask_data): histogram_decrement(histo, &pop, image[rr, cc]) - kernel(out[r, c, :], histo, pop, image[r, c], max_bin, mid_bin, + kernel(&out[r, c, 0], odepth, histo, pop, image[r, c], max_bin, mid_bin, p0, p1, s0, s1) # ---> east to west @@ -209,7 +210,7 @@ cdef void _core(void kernel(dtype_t_out[:] out_feat, Py_ssize_t*, double, dtype_ if is_in_mask(rows, cols, rr, cc, mask_data): histogram_decrement(histo, &pop, image[rr, cc]) - kernel(out[r, c, :], histo, pop, image[r, c], max_bin, mid_bin, + kernel(&out[r, c, 0], odepth, histo, pop, image[r, c], max_bin, mid_bin, p0, p1, s0, s1) r += 1 # pass to the next row @@ -229,7 +230,7 @@ cdef void _core(void kernel(dtype_t_out[:] out_feat, Py_ssize_t*, double, dtype_ if is_in_mask(rows, cols, rr, cc, mask_data): histogram_decrement(histo, &pop, image[rr, cc]) - kernel(out[r, c, :], histo, pop, image[r, c], + kernel(&out[r, c, 0], odepth, histo, pop, image[r, c], max_bin, mid_bin, p0, p1, s0, s1) # release memory allocated by malloc diff --git a/skimage/filter/rank/generic_cy.pyx b/skimage/filter/rank/generic_cy.pyx index bafd2432..bc54855f 100644 --- a/skimage/filter/rank/generic_cy.pyx +++ b/skimage/filter/rank/generic_cy.pyx @@ -9,7 +9,8 @@ from libc.math cimport log from .core_cy cimport dtype_t, dtype_t_out, _core -cdef inline void _kernel_autolevel(dtype_t_out[:] out, Py_ssize_t* histo, +cdef inline void _kernel_autolevel(dtype_t_out* out, Py_ssize_t odepth, + Py_ssize_t* histo, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, @@ -35,7 +36,8 @@ cdef inline void _kernel_autolevel(dtype_t_out[:] out, Py_ssize_t* histo, out[0] = 0 -cdef inline void _kernel_bottomhat(dtype_t_out[:] out, Py_ssize_t* histo, +cdef inline void _kernel_bottomhat(dtype_t_out* out, Py_ssize_t odepth, + Py_ssize_t* histo, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, @@ -52,7 +54,8 @@ cdef inline void _kernel_bottomhat(dtype_t_out[:] out, Py_ssize_t* histo, out[0] = 0 -cdef inline void _kernel_equalize(dtype_t_out[:] out, Py_ssize_t* histo, +cdef inline void _kernel_equalize(dtype_t_out* out, Py_ssize_t odepth, + Py_ssize_t* histo, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, @@ -71,7 +74,8 @@ cdef inline void _kernel_equalize(dtype_t_out[:] out, Py_ssize_t* histo, out[0] = 0 -cdef inline void _kernel_gradient(dtype_t_out[:] out, Py_ssize_t* histo, +cdef inline void _kernel_gradient(dtype_t_out* out, Py_ssize_t odepth, + Py_ssize_t* histo, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, @@ -93,7 +97,8 @@ cdef inline void _kernel_gradient(dtype_t_out[:] out, Py_ssize_t* histo, out[0] = 0 -cdef inline void _kernel_maximum(dtype_t_out[:] out, Py_ssize_t* histo, +cdef inline void _kernel_maximum(dtype_t_out* out, Py_ssize_t odepth, + Py_ssize_t* histo, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, @@ -110,7 +115,8 @@ cdef inline void _kernel_maximum(dtype_t_out[:] out, Py_ssize_t* histo, out[0] = 0 -cdef inline void _kernel_mean(dtype_t_out[:] out, Py_ssize_t* histo, +cdef inline void _kernel_mean(dtype_t_out* out, Py_ssize_t odepth, + Py_ssize_t* histo, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, @@ -127,7 +133,8 @@ cdef inline void _kernel_mean(dtype_t_out[:] out, Py_ssize_t* histo, out[0] = 0 -cdef inline void _kernel_subtract_mean(dtype_t_out[:] out, Py_ssize_t* histo, +cdef inline void _kernel_subtract_mean(dtype_t_out* out, Py_ssize_t odepth, + Py_ssize_t* histo, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, @@ -144,7 +151,8 @@ cdef inline void _kernel_subtract_mean(dtype_t_out[:] out, Py_ssize_t* histo, out[0] = 0 -cdef inline void _kernel_median(dtype_t_out[:] out, Py_ssize_t* histo, +cdef inline void _kernel_median(dtype_t_out* out, Py_ssize_t odepth, + Py_ssize_t* histo, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, @@ -164,7 +172,8 @@ cdef inline void _kernel_median(dtype_t_out[:] out, Py_ssize_t* histo, out[0] = 0 -cdef inline void _kernel_minimum(dtype_t_out[:] out, Py_ssize_t* histo, +cdef inline void _kernel_minimum(dtype_t_out* out, Py_ssize_t odepth, + Py_ssize_t* histo, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, @@ -181,7 +190,8 @@ cdef inline void _kernel_minimum(dtype_t_out[:] out, Py_ssize_t* histo, out[0] = 0 -cdef inline void _kernel_modal(dtype_t_out[:] out, Py_ssize_t* histo, +cdef inline void _kernel_modal(dtype_t_out* out, Py_ssize_t odepth, + Py_ssize_t* histo, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, @@ -199,7 +209,8 @@ cdef inline void _kernel_modal(dtype_t_out[:] out, Py_ssize_t* histo, out[0] = 0 -cdef inline void _kernel_enhance_contrast(dtype_t_out[:] out, +cdef inline void _kernel_enhance_contrast(dtype_t_out* out, + Py_ssize_t odepth, Py_ssize_t* histo, double pop, dtype_t g, @@ -227,7 +238,8 @@ cdef inline void _kernel_enhance_contrast(dtype_t_out[:] out, out[0] = 0 -cdef inline void _kernel_pop(dtype_t_out[:] out, Py_ssize_t* histo, +cdef inline void _kernel_pop(dtype_t_out* out, Py_ssize_t odepth, + Py_ssize_t* histo, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, @@ -236,7 +248,8 @@ cdef inline void _kernel_pop(dtype_t_out[:] out, Py_ssize_t* histo, out[0] = pop -cdef inline void _kernel_sum(dtype_t_out[:] out, Py_ssize_t* histo, +cdef inline void _kernel_sum(dtype_t_out* out, Py_ssize_t odepth, + Py_ssize_t* histo, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, @@ -253,7 +266,8 @@ cdef inline void _kernel_sum(dtype_t_out[:] out, Py_ssize_t* histo, out[0] = 0 -cdef inline void _kernel_threshold(dtype_t_out[:] out, Py_ssize_t* histo, +cdef inline void _kernel_threshold(dtype_t_out* out, Py_ssize_t odepth, + Py_ssize_t* histo, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, @@ -270,7 +284,8 @@ cdef inline void _kernel_threshold(dtype_t_out[:] out, Py_ssize_t* histo, out[0] = 0 -cdef inline void _kernel_tophat(dtype_t_out[:] out, Py_ssize_t* histo, +cdef inline void _kernel_tophat(dtype_t_out* out, Py_ssize_t odepth, + Py_ssize_t* histo, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, @@ -287,7 +302,8 @@ cdef inline void _kernel_tophat(dtype_t_out[:] out, Py_ssize_t* histo, out[0] = 0 -cdef inline void _kernel_noise_filter(dtype_t_out[:] out, Py_ssize_t* histo, +cdef inline void _kernel_noise_filter(dtype_t_out* out, Py_ssize_t odepth, + Py_ssize_t* histo, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, @@ -313,7 +329,8 @@ cdef inline void _kernel_noise_filter(dtype_t_out[:] out, Py_ssize_t* histo, out[0] = min_i -cdef inline void _kernel_entropy(dtype_t_out[:] out, Py_ssize_t* histo, +cdef inline void _kernel_entropy(dtype_t_out* out, Py_ssize_t odepth, + Py_ssize_t* histo, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, @@ -332,7 +349,8 @@ cdef inline void _kernel_entropy(dtype_t_out[:] out, Py_ssize_t* histo, out[0] = 0 -cdef inline void _kernel_otsu(dtype_t_out[:] out, Py_ssize_t* histo, +cdef inline void _kernel_otsu(dtype_t_out* out, Py_ssize_t odepth, + Py_ssize_t* histo, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, @@ -371,7 +389,8 @@ cdef inline void _kernel_otsu(dtype_t_out[:] out, Py_ssize_t* histo, out[0] = max_i -cdef inline void _kernel_win_hist(dtype_t_out[:] out, Py_ssize_t* histo, +cdef inline void _kernel_win_hist(dtype_t_out* out, Py_ssize_t odepth, + Py_ssize_t* histo, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, @@ -381,10 +400,10 @@ cdef inline void _kernel_win_hist(dtype_t_out[:] out, Py_ssize_t* histo, cdef double scale if pop: scale = 1.0 / pop - for i in xrange(out.shape[0]): + for i in xrange(odepth): out[i] = (histo[i] * scale) else: - for i in xrange(out.shape[0]): + for i in xrange(odepth): out[i] = 0 diff --git a/skimage/filter/rank/percentile_cy.pyx b/skimage/filter/rank/percentile_cy.pyx index af33546d..6d54999f 100644 --- a/skimage/filter/rank/percentile_cy.pyx +++ b/skimage/filter/rank/percentile_cy.pyx @@ -7,7 +7,8 @@ cimport numpy as cnp from .core_cy cimport dtype_t, dtype_t_out, _core, _min, _max -cdef inline void _kernel_autolevel(dtype_t_out[:] out, Py_ssize_t* histo, +cdef inline void _kernel_autolevel(dtype_t_out* out, Py_ssize_t odepth, + Py_ssize_t* histo, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, @@ -40,7 +41,8 @@ cdef inline void _kernel_autolevel(dtype_t_out[:] out, Py_ssize_t* histo, out[0] = 0 -cdef inline void _kernel_gradient(dtype_t_out[:] out, Py_ssize_t* histo, +cdef inline void _kernel_gradient(dtype_t_out* out, Py_ssize_t odepth, + Py_ssize_t* histo, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, @@ -68,7 +70,8 @@ cdef inline void _kernel_gradient(dtype_t_out[:] out, Py_ssize_t* histo, out[0] = 0 -cdef inline void _kernel_mean(dtype_t_out[:] out, Py_ssize_t* histo, +cdef inline void _kernel_mean(dtype_t_out* out, Py_ssize_t odepth, + Py_ssize_t* histo, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, @@ -93,7 +96,8 @@ cdef inline void _kernel_mean(dtype_t_out[:] out, Py_ssize_t* histo, else: out[0] = 0 -cdef inline void _kernel_sum(dtype_t_out[:] out, Py_ssize_t* histo, +cdef inline void _kernel_sum(dtype_t_out* out, Py_ssize_t odepth, + Py_ssize_t* histo, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, @@ -118,7 +122,8 @@ cdef inline void _kernel_sum(dtype_t_out[:] out, Py_ssize_t* histo, else: out[0] = 0 -cdef inline void _kernel_subtract_mean(dtype_t_out[:] out, Py_ssize_t* histo, +cdef inline void _kernel_subtract_mean(dtype_t_out* out, Py_ssize_t odepth, + Py_ssize_t* histo, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, @@ -144,7 +149,8 @@ cdef inline void _kernel_subtract_mean(dtype_t_out[:] out, Py_ssize_t* histo, out[0] = 0 -cdef inline void _kernel_enhance_contrast(dtype_t_out[:] out, +cdef inline void _kernel_enhance_contrast(dtype_t_out* out, + Py_ssize_t odepth, Py_ssize_t* histo, double pop, dtype_t g, Py_ssize_t max_bin, @@ -180,7 +186,8 @@ cdef inline void _kernel_enhance_contrast(dtype_t_out[:] out, out[0] = 0 -cdef inline void _kernel_percentile(dtype_t_out[:] out, Py_ssize_t* histo, +cdef inline void _kernel_percentile(dtype_t_out* out, Py_ssize_t odepth, + Py_ssize_t* histo, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, @@ -204,7 +211,8 @@ cdef inline void _kernel_percentile(dtype_t_out[:] out, Py_ssize_t* histo, out[0] = 0 -cdef inline void _kernel_pop(dtype_t_out[:] out, Py_ssize_t* histo, +cdef inline void _kernel_pop(dtype_t_out* out, Py_ssize_t odepth, + Py_ssize_t* histo, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, @@ -224,7 +232,8 @@ cdef inline void _kernel_pop(dtype_t_out[:] out, Py_ssize_t* histo, out[0] = 0 -cdef inline void _kernel_threshold(dtype_t_out[:] out, Py_ssize_t* histo, +cdef inline void _kernel_threshold(dtype_t_out* out, Py_ssize_t odepth, + Py_ssize_t* histo, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, From 0936522064dd15daf856885b34fd6e513a55559e Mon Sep 17 00:00:00 2001 From: Geoffrey French Date: Mon, 1 Sep 2014 22:44:45 +0100 Subject: [PATCH 07/13] Improved plot_windowed_histogram doctoring, removed unnecessary imports and removed an internal testing line. --- doc/examples/plot_windowed_histogram.py | 15 ++++++--------- 1 file changed, 6 insertions(+), 9 deletions(-) diff --git a/doc/examples/plot_windowed_histogram.py b/doc/examples/plot_windowed_histogram.py index cfa9345d..6d8f9ed0 100644 --- a/doc/examples/plot_windowed_histogram.py +++ b/doc/examples/plot_windowed_histogram.py @@ -4,13 +4,13 @@ from __future__ import division Sliding window histogram ======================== -This example extracts a single coin from the coins image and generates a -histogram of its greyscale values. +This example extracts a single coin from the `skimage.data.coins` image and +generates a histogram of its greyscale values. It then computes a sliding window histogram of the complete image using -rank.windowed_histogram. The local histogram for the region surrounding -each pixel in the image is compared to that of the single coin, with -a similarity measure being computed and displayed. +`skimage.filter.rank.windowed_histogram`. The local histogram for the region +surrounding each pixel in the image is compared to that of the single coin, +with a similarity measure being computed and displayed. To demonstrate the rotational invariance of the technique, the same test is performed on a version of the coins image rotated by 45 degrees. @@ -20,9 +20,7 @@ import matplotlib import matplotlib.pyplot as plt from skimage import data -from skimage.util.dtype import dtype_range from skimage.util import img_as_ubyte -from skimage import exposure from skimage.morphology import disk from skimage.filter import rank from skimage import transform @@ -57,9 +55,8 @@ def windowed_histogram_similarity(image, selem, reference_hist, n_bins): return similarity -# Load the coins image +# Load the `skimage.data.coins` image img = img_as_ubyte(data.coins()) -# img = img_as_ubyte(plt.imread('../../skimage/data/coins.png')) # Quantize to 16 levels of grayscale; this way the output image will have a # 16-dimensional feature vector per pixel From 9f5a6fddbfe21c9d3ba4bb4caa8310c44967e9b5 Mon Sep 17 00:00:00 2001 From: Geoffrey French Date: Mon, 1 Sep 2014 22:46:12 +0100 Subject: [PATCH 08/13] Fixed comment in plot_windowed_histogram. --- doc/examples/plot_windowed_histogram.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/examples/plot_windowed_histogram.py b/doc/examples/plot_windowed_histogram.py index 6d8f9ed0..73cf82a7 100644 --- a/doc/examples/plot_windowed_histogram.py +++ b/doc/examples/plot_windowed_histogram.py @@ -37,7 +37,7 @@ def windowed_histogram_similarity(image, selem, reference_hist, n_bins): # arithmetic operations with the windowed histograms fro the image reference_hist = reference_hist.reshape((1,1) + reference_hist.shape) - # Compute Chi squared distance metric: sum((X-Y)**2 / (X+Y); + # Compute Chi squared distance metric: sum((X-Y)^2 / (X+Y)); # a measure of distance between histograms X = px_histograms Y = reference_hist From a90096555ba7c252bf5dabea8d5e98f7f1d60fc1 Mon Sep 17 00:00:00 2001 From: Geoffrey French Date: Mon, 1 Sep 2014 22:59:59 +0100 Subject: [PATCH 09/13] Docstring and comment improvements and fixes in plot_windowed_histogram. Readability improvement to skimage/io/__init__.py --- doc/examples/plot_windowed_histogram.py | 10 ++++++++-- skimage/io/__init__.py | 5 ++++- 2 files changed, 12 insertions(+), 3 deletions(-) diff --git a/doc/examples/plot_windowed_histogram.py b/doc/examples/plot_windowed_histogram.py index 73cf82a7..8e2c4384 100644 --- a/doc/examples/plot_windowed_histogram.py +++ b/doc/examples/plot_windowed_histogram.py @@ -12,6 +12,12 @@ It then computes a sliding window histogram of the complete image using surrounding each pixel in the image is compared to that of the single coin, with a similarity measure being computed and displayed. +The histogram of the single coin is computed using `numpy.histogram` on a +box shaped region surrounding the coin, while the sliding window histograms +are computed using a disc shaped structural element of a slightly different +size. This is done in aid of demonstrating that the technique still finds +similarity inspite of these differences. + To demonstrate the rotational invariance of the technique, the same test is performed on a version of the coins image rotated by 45 degrees. """ @@ -34,7 +40,7 @@ def windowed_histogram_similarity(image, selem, reference_hist, n_bins): px_histograms = rank.windowed_histogram(image, selem, n_bins=n_bins) # Reshape coin histogram to (1,1,N) for broadcast when we want to use it in - # arithmetic operations with the windowed histograms fro the image + # arithmetic operations with the windowed histograms from the image reference_hist = reference_hist.reshape((1,1) + reference_hist.shape) # Compute Chi squared distance metric: sum((X-Y)^2 / (X+Y)); @@ -47,7 +53,7 @@ def windowed_histogram_similarity(image, selem, reference_hist, n_bins): frac[denom==0] = 0 chi_sqr = np.sum(frac, axis=2) * 0.5 - # Generate a similarity measure. It needs to be low when distance is high. + # Generate a similarity measure. It needs to be low when distance is high # and high when distance is low; taking the reciprocal will do this. # Chi squared will always be >= 0, add small value to prevent divide by 0. similarity = 1 / (chi_sqr + 1.0e-4) diff --git a/skimage/io/__init__.py b/skimage/io/__init__.py index 8f3b9bca..85461216 100644 --- a/skimage/io/__init__.py +++ b/skimage/io/__init__.py @@ -40,7 +40,10 @@ def _update_doc(doc): info_table = [(p, plugin_info(p).get('description', 'no description')) for p in available_plugins if not p == 'test'] - name_length = max([len(n) for (n, _) in info_table]) if len(info_table) > 0 else 0 + if len(info_table) > 0: + name_length = max([len(n) for (n, _) in info_table]) + else: + name_length = 0 description_length = WRAP_LEN - 1 - name_length column_lengths = [name_length, description_length] From bf02a92ee399208874171fa251ee92e55cc61515 Mon Sep 17 00:00:00 2001 From: Geoffrey French Date: Mon, 1 Sep 2014 23:04:24 +0100 Subject: [PATCH 10/13] Small fix to windowed_histogram doctoring. --- skimage/filter/rank/generic.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/skimage/filter/rank/generic.py b/skimage/filter/rank/generic.py index 2af2428e..cec1624d 100644 --- a/skimage/filter/rank/generic.py +++ b/skimage/filter/rank/generic.py @@ -887,8 +887,8 @@ def windowed_histogram(image, selem, out=None, mask=None, shift_x=False, shift_y to the structuring element sizes (center must be inside the given structuring element). n_bins : int or None - The number of histogram bins. Will default to image.max() + 1 if None - is passed. + The number of histogram bins. Will default to `image.max() + 1` + if None is passed. Returns ------- From e6bda5accd0d9133e37eb1d9ed456c5191d8db19 Mon Sep 17 00:00:00 2001 From: Geoffrey French Date: Mon, 1 Sep 2014 23:41:03 +0100 Subject: [PATCH 11/13] Fix to skimage.filter.rank.windowed_histogram docstring. Better explanation of technique in plot_windowed_histogram example, along with (hopefully correct) citations. Relevant additions to release_dev.txt and CONTRIBUTORS.txt. --- CONTRIBUTORS.txt | 3 +++ doc/examples/plot_windowed_histogram.py | 30 +++++++++++++++++++------ doc/release/release_dev.txt | 1 + skimage/filter/rank/generic.py | 8 +++---- 4 files changed, 31 insertions(+), 11 deletions(-) diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index 30f432f4..b0a82519 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -185,3 +185,6 @@ - Adam Feuer PIL Image import and export improvements + +- Geoffrey French + skimage.filters.rank.windowed_histogram and plot_windowed_histogram example. \ No newline at end of file diff --git a/doc/examples/plot_windowed_histogram.py b/doc/examples/plot_windowed_histogram.py index 8e2c4384..2cc77139 100644 --- a/doc/examples/plot_windowed_histogram.py +++ b/doc/examples/plot_windowed_histogram.py @@ -4,22 +4,38 @@ from __future__ import division Sliding window histogram ======================== -This example extracts a single coin from the `skimage.data.coins` image and -generates a histogram of its greyscale values. +Histogram matching can be used for object detection in images [1]_. +This example extracts a single coin from the `skimage.data.coins` image +and uses histogram matching to attempt to locate it within the original +image. -It then computes a sliding window histogram of the complete image using -`skimage.filter.rank.windowed_histogram`. The local histogram for the region -surrounding each pixel in the image is compared to that of the single coin, -with a similarity measure being computed and displayed. +First, a box-shaped region of the image containing the target coin is +extracted and a histogram of its greyscale values is computed. + +Next, for each pixel in the test image, a histogram of the greyscale values +in a region of the image surrounding the pixel is computed. +`skimage.filter.rank.windowed_histogram` is used for this task, as it +employs an efficient sliding window based algorithm that is able to compute +these histograms quickly [2]_. +The local histogram for the region surrounding each pixel in the image is +compared to that of the single coin, with a similarity measure being +computed and displayed. The histogram of the single coin is computed using `numpy.histogram` on a box shaped region surrounding the coin, while the sliding window histograms are computed using a disc shaped structural element of a slightly different size. This is done in aid of demonstrating that the technique still finds -similarity inspite of these differences. +similarity in spite of these differences. To demonstrate the rotational invariance of the technique, the same test is performed on a version of the coins image rotated by 45 degrees. + +References +---------- +.. [1] Porikli, F. "Integral Histogram: A Fast Way to Extract Histograms + in Cartesian Spaces" CVPR, 2005. Vol. 1. IEEE, 2005 +.. [2] S.Perreault and P.Hebert. Median filtering in constant time. + Trans. Image Processing, 16(9):2389-2394, 2007. """ import numpy as np import matplotlib diff --git a/doc/release/release_dev.txt b/doc/release/release_dev.txt index cbeded27..77c8d102 100644 --- a/doc/release/release_dev.txt +++ b/doc/release/release_dev.txt @@ -20,6 +20,7 @@ Region Adjacency Graphs - Similarity RAGs (#1080) - Normalized Cut on RAGs (#1080) - RAG Drawing (#1087) +Sliding Windowed Histogram (#1127) Improvements ------------ diff --git a/skimage/filter/rank/generic.py b/skimage/filter/rank/generic.py index cec1624d..018d7c0b 100644 --- a/skimage/filter/rank/generic.py +++ b/skimage/filter/rank/generic.py @@ -895,10 +895,10 @@ def windowed_histogram(image, selem, out=None, mask=None, shift_x=False, shift_y out : 3-D array with float dtype of dimensions (H,W,N), where (H,W) are the dimensions of the input image and N is n_bins or image.max()+1 if no value is provided as a parameter. Effectively, each pixel - is an N-dimensional feature vector that is the histogram. - The sum of the elements in the feature vector will be 1, unless - no pixels in the window were covered by both selem and mask, in which - case all elements will be 0. + is a N-D feature vector that is the histogram. The sum of the + elements in the feature vector will be 1, unless no pixels in the + window were covered by both selem and mask, in which case all + elements will be 0. Examples -------- From 5342299572dc26b7bff56ec0999e1158860b5e2a Mon Sep 17 00:00:00 2001 From: Geoffrey French Date: Tue, 2 Sep 2014 19:16:41 +0100 Subject: [PATCH 12/13] PEP8 compliance and doc formatting fixes. --- doc/examples/plot_windowed_histogram.py | 14 ++--- skimage/filter/rank/bilateral.py | 3 +- skimage/filter/rank/bilateral_cy.pyx | 2 + skimage/filter/rank/core_cy.pyx | 16 ++--- skimage/filter/rank/generic.py | 84 ++++++++++++++++--------- skimage/filter/rank/generic_cy.pyx | 2 +- skimage/filter/rank/percentile_cy.pyx | 2 + 7 files changed, 75 insertions(+), 48 deletions(-) diff --git a/doc/examples/plot_windowed_histogram.py b/doc/examples/plot_windowed_histogram.py index 2cc77139..1c57721b 100644 --- a/doc/examples/plot_windowed_histogram.py +++ b/doc/examples/plot_windowed_histogram.py @@ -57,7 +57,7 @@ def windowed_histogram_similarity(image, selem, reference_hist, n_bins): # Reshape coin histogram to (1,1,N) for broadcast when we want to use it in # arithmetic operations with the windowed histograms from the image - reference_hist = reference_hist.reshape((1,1) + reference_hist.shape) + reference_hist = reference_hist.reshape((1, 1) + reference_hist.shape) # Compute Chi squared distance metric: sum((X-Y)^2 / (X+Y)); # a measure of distance between histograms @@ -66,7 +66,7 @@ def windowed_histogram_similarity(image, selem, reference_hist, n_bins): num = (X-Y)*(X-Y) denom = X+Y frac = num / denom - frac[denom==0] = 0 + frac[denom == 0] = 0 chi_sqr = np.sum(frac, axis=2) * 0.5 # Generate a similarity measure. It needs to be low when distance is high @@ -80,17 +80,18 @@ def windowed_histogram_similarity(image, selem, reference_hist, n_bins): # Load the `skimage.data.coins` image img = img_as_ubyte(data.coins()) -# Quantize to 16 levels of grayscale; this way the output image will have a +# Quantize to 16 levels of greyscale; this way the output image will have a # 16-dimensional feature vector per pixel quantized_img = img//16 # Select the coin from the 4th column, second row. # Co-ordinate ordering: [x1,y1,x2,y2] -coin_coords = [184,100,228,148] # 44 x 44 region -coin = quantized_img[coin_coords[1]:coin_coords[3], coin_coords[0]:coin_coords[2]] +coin_coords = [184, 100, 228, 148] # 44 x 44 region +coin = quantized_img[coin_coords[1]:coin_coords[3], + coin_coords[0]:coin_coords[2]] # Compute coin histogram and normalize -coin_hist, _ = np.histogram(coin.flatten(), bins=16, range=(0,16)) +coin_hist, _ = np.histogram(coin.flatten(), bins=16, range=(0, 16)) coin_hist = coin_hist.astype(float) / np.sum(coin_hist) @@ -114,7 +115,6 @@ rotated_similarity = windowed_histogram_similarity(quantized_rotated_image, coin_hist.shape[0]) - # Plot it all fig, axes = plt.subplots(nrows=5, figsize=(6, 18)) ax0, ax1, ax2, ax3, ax4 = axes diff --git a/skimage/filter/rank/bilateral.py b/skimage/filter/rank/bilateral.py index 73147c7f..aeb318d1 100644 --- a/skimage/filter/rank/bilateral.py +++ b/skimage/filter/rank/bilateral.py @@ -158,8 +158,9 @@ def pop_bilateral(image, selem, out=None, mask=None, shift_x=False, return _apply(bilateral_cy._pop, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y, s0=s0, s1=s1) + def sum_bilateral(image, selem, out=None, mask=None, shift_x=False, - shift_y=False, s0=10, s1=10): + shift_y=False, s0=10, s1=10): """Apply a flat kernel bilateral filter. This is an edge-preserving and noise reducing denoising filter. It averages diff --git a/skimage/filter/rank/bilateral_cy.pyx b/skimage/filter/rank/bilateral_cy.pyx index 20a97788..b4e22d3f 100644 --- a/skimage/filter/rank/bilateral_cy.pyx +++ b/skimage/filter/rank/bilateral_cy.pyx @@ -51,6 +51,7 @@ cdef inline void _kernel_pop(dtype_t_out* out, Py_ssize_t odepth, else: out[0] = 0 + cdef inline void _kernel_sum(dtype_t_out* out, Py_ssize_t odepth, Py_ssize_t* histo, double pop, dtype_t g, @@ -96,6 +97,7 @@ def _pop(dtype_t[:, ::1] image, _core(_kernel_pop[dtype_t_out, dtype_t], image, selem, mask, out, shift_x, shift_y, 0, 0, s0, s1, max_bin) + def _sum(dtype_t[:, ::1] image, char[:, ::1] selem, char[:, ::1] mask, diff --git a/skimage/filter/rank/core_cy.pyx b/skimage/filter/rank/core_cy.pyx index c1216c0d..4bd21df0 100644 --- a/skimage/filter/rank/core_cy.pyx +++ b/skimage/filter/rank/core_cy.pyx @@ -42,8 +42,8 @@ cdef inline char is_in_mask(Py_ssize_t rows, Py_ssize_t cols, return 0 -cdef void _core(void kernel(dtype_t_out*, Py_ssize_t, Py_ssize_t*, double, dtype_t, - Py_ssize_t, Py_ssize_t, double, +cdef void _core(void kernel(dtype_t_out*, Py_ssize_t, Py_ssize_t*, double, + dtype_t, Py_ssize_t, Py_ssize_t, double, double, Py_ssize_t, Py_ssize_t), dtype_t[:, ::1] image, char[:, ::1] selem, @@ -173,8 +173,8 @@ cdef void _core(void kernel(dtype_t_out*, Py_ssize_t, Py_ssize_t*, double, dtype if is_in_mask(rows, cols, rr, cc, mask_data): histogram_decrement(histo, &pop, image[rr, cc]) - kernel(&out[r, c, 0], odepth, histo, pop, image[r, c], max_bin, mid_bin, - p0, p1, s0, s1) + kernel(&out[r, c, 0], odepth, histo, pop, image[r, c], max_bin, + mid_bin, p0, p1, s0, s1) r += 1 # pass to the next row if r >= rows: @@ -193,8 +193,8 @@ cdef void _core(void kernel(dtype_t_out*, Py_ssize_t, Py_ssize_t*, double, dtype if is_in_mask(rows, cols, rr, cc, mask_data): histogram_decrement(histo, &pop, image[rr, cc]) - kernel(&out[r, c, 0], odepth, histo, pop, image[r, c], max_bin, mid_bin, - p0, p1, s0, s1) + kernel(&out[r, c, 0], odepth, histo, pop, image[r, c], max_bin, + mid_bin, p0, p1, s0, s1) # ---> east to west for c in range(cols - 2, -1, -1): @@ -210,8 +210,8 @@ cdef void _core(void kernel(dtype_t_out*, Py_ssize_t, Py_ssize_t*, double, dtype if is_in_mask(rows, cols, rr, cc, mask_data): histogram_decrement(histo, &pop, image[rr, cc]) - kernel(&out[r, c, 0], odepth, histo, pop, image[r, c], max_bin, mid_bin, - p0, p1, s0, s1) + kernel(&out[r, c, 0], odepth, histo, pop, image[r, c], max_bin, + mid_bin, p0, p1, s0, s1) r += 1 # pass to the next row if r >= rows: diff --git a/skimage/filter/rank/generic.py b/skimage/filter/rank/generic.py index 018d7c0b..f2e75051 100644 --- a/skimage/filter/rank/generic.py +++ b/skimage/filter/rank/generic.py @@ -68,7 +68,8 @@ def _handle_input(image, selem, out, mask, out_dtype=None, pixel_size=1): return image, selem, out, mask, max_bin -def _apply_scalar_per_pixel(func, image, selem, out, mask, shift_x, shift_y, out_dtype=None): +def _apply_scalar_per_pixel(func, image, selem, out, mask, shift_x, shift_y, + out_dtype=None): image, selem, out, mask, max_bin = _handle_input(image, selem, out, mask, out_dtype) @@ -79,10 +80,12 @@ def _apply_scalar_per_pixel(func, image, selem, out, mask, shift_x, shift_y, out return out.reshape(out.shape[:2]) -def _apply_vector_per_pixel(func, image, selem, out, mask, shift_x, shift_y, out_dtype=None, pixel_size=1): +def _apply_vector_per_pixel(func, image, selem, out, mask, shift_x, shift_y, + out_dtype=None, pixel_size=1): image, selem, out, mask, max_bin = _handle_input(image, selem, out, mask, - out_dtype, pixel_size=pixel_size) + out_dtype, + pixel_size=pixel_size) func(image, selem, shift_x=shift_x, shift_y=shift_y, mask=mask, out=out, max_bin=max_bin) @@ -128,7 +131,8 @@ def autolevel(image, selem, out=None, mask=None, shift_x=False, shift_y=False): """ return _apply_scalar_per_pixel(generic_cy._autolevel, image, selem, - out=out, mask=mask, shift_x=shift_x, shift_y=shift_y) + out=out, mask=mask, + shift_x=shift_x, shift_y=shift_y) def bottomhat(image, selem, out=None, mask=None, shift_x=False, shift_y=False): @@ -169,7 +173,8 @@ def bottomhat(image, selem, out=None, mask=None, shift_x=False, shift_y=False): """ return _apply_scalar_per_pixel(generic_cy._bottomhat, image, selem, - out=out, mask=mask, shift_x=shift_x, shift_y=shift_y) + out=out, mask=mask, + shift_x=shift_x, shift_y=shift_y) def equalize(image, selem, out=None, mask=None, shift_x=False, shift_y=False): @@ -207,7 +212,8 @@ def equalize(image, selem, out=None, mask=None, shift_x=False, shift_y=False): """ return _apply_scalar_per_pixel(generic_cy._equalize, image, selem, - out=out, mask=mask, shift_x=shift_x, shift_y=shift_y) + out=out, mask=mask, + shift_x=shift_x, shift_y=shift_y) def gradient(image, selem, out=None, mask=None, shift_x=False, shift_y=False): @@ -245,7 +251,8 @@ def gradient(image, selem, out=None, mask=None, shift_x=False, shift_y=False): """ return _apply_scalar_per_pixel(generic_cy._gradient, image, selem, - out=out, mask=mask, shift_x=shift_x, shift_y=shift_y) + out=out, mask=mask, + shift_x=shift_x, shift_y=shift_y) def maximum(image, selem, out=None, mask=None, shift_x=False, shift_y=False): @@ -292,7 +299,8 @@ def maximum(image, selem, out=None, mask=None, shift_x=False, shift_y=False): """ return _apply_scalar_per_pixel(generic_cy._maximum, image, selem, - out=out, mask=mask, shift_x=shift_x, shift_y=shift_y) + out=out, mask=mask, + shift_x=shift_x, shift_y=shift_y) def mean(image, selem, out=None, mask=None, shift_x=False, shift_y=False): @@ -330,7 +338,8 @@ def mean(image, selem, out=None, mask=None, shift_x=False, shift_y=False): """ return _apply_scalar_per_pixel(generic_cy._mean, image, selem, out=out, - mask=mask, shift_x=shift_x, shift_y=shift_y) + out=out, mask=mask, + shift_x=shift_x, shift_y=shift_y) def subtract_mean(image, selem, out=None, mask=None, shift_x=False, @@ -369,7 +378,8 @@ def subtract_mean(image, selem, out=None, mask=None, shift_x=False, """ return _apply_scalar_per_pixel(generic_cy._subtract_mean, image, selem, - out=out, mask=mask, shift_x=shift_x, shift_y=shift_y) + out=out, mask=mask, + shift_x=shift_x, shift_y=shift_y) def median(image, selem, out=None, mask=None, shift_x=False, shift_y=False): @@ -407,7 +417,8 @@ def median(image, selem, out=None, mask=None, shift_x=False, shift_y=False): """ return _apply_scalar_per_pixel(generic_cy._median, image, selem, - out=out, mask=mask, shift_x=shift_x, shift_y=shift_y) + out=out, mask=mask, + shift_x=shift_x, shift_y=shift_y) def minimum(image, selem, out=None, mask=None, shift_x=False, shift_y=False): @@ -454,7 +465,8 @@ def minimum(image, selem, out=None, mask=None, shift_x=False, shift_y=False): """ return _apply_scalar_per_pixel(generic_cy._minimum, image, selem, - out=out, mask=mask, shift_x=shift_x, shift_y=shift_y) + out=out, mask=mask, + shift_x=shift_x, shift_y=shift_y) def modal(image, selem, out=None, mask=None, shift_x=False, shift_y=False): @@ -494,7 +506,8 @@ def modal(image, selem, out=None, mask=None, shift_x=False, shift_y=False): """ return _apply_scalar_per_pixel(generic_cy._modal, image, selem, - out=out, mask=mask, shift_x=shift_x, shift_y=shift_y) + out=out, mask=mask, + shift_x=shift_x, shift_y=shift_y) def enhance_contrast(image, selem, out=None, mask=None, shift_x=False, @@ -537,7 +550,8 @@ def enhance_contrast(image, selem, out=None, mask=None, shift_x=False, """ return _apply_scalar_per_pixel(generic_cy._enhance_contrast, image, selem, - out=out, mask=mask, shift_x=shift_x, shift_y=shift_y) + out=out, mask=mask, + shift_x=shift_x, shift_y=shift_y) def pop(image, selem, out=None, mask=None, shift_x=False, shift_y=False): @@ -586,7 +600,8 @@ def pop(image, selem, out=None, mask=None, shift_x=False, shift_y=False): """ return _apply_scalar_per_pixel(generic_cy._pop, image, selem, out=out, - mask=mask, shift_x=shift_x, shift_y=shift_y) + mask=mask, shift_x=shift_x, + shift_y=shift_y) def sum(image, selem, out=None, mask=None, shift_x=False, shift_y=False): @@ -635,7 +650,8 @@ def sum(image, selem, out=None, mask=None, shift_x=False, shift_y=False): """ return _apply_scalar_per_pixel(generic_cy._sum, image, selem, out=out, - mask=mask, shift_x=shift_x, shift_y=shift_y) + mask=mask, shift_x=shift_x, + shift_y=shift_y) def threshold(image, selem, out=None, mask=None, shift_x=False, shift_y=False): @@ -684,7 +700,8 @@ def threshold(image, selem, out=None, mask=None, shift_x=False, shift_y=False): """ return _apply_scalar_per_pixel(generic_cy._threshold, image, selem, - out=out, mask=mask, shift_x=shift_x, shift_y=shift_y) + out=out, mask=mask, + shift_x=shift_x, shift_y=shift_y) def tophat(image, selem, out=None, mask=None, shift_x=False, shift_y=False): @@ -725,7 +742,8 @@ def tophat(image, selem, out=None, mask=None, shift_x=False, shift_y=False): """ return _apply_scalar_per_pixel(generic_cy._tophat, image, selem, - out=out, mask=mask, shift_x=shift_x, shift_y=shift_y) + out=out, mask=mask, + shift_x=shift_x, shift_y=shift_y) def noise_filter(image, selem, out=None, mask=None, shift_x=False, @@ -775,8 +793,9 @@ def noise_filter(image, selem, out=None, mask=None, shift_x=False, selem_cpy = selem.copy() selem_cpy[centre_r, centre_c] = 0 - return _apply_scalar_per_pixel(generic_cy._noise_filter, image, selem_cpy, out=out, - mask=mask, shift_x=shift_x, shift_y=shift_y) + return _apply_scalar_per_pixel(generic_cy._noise_filter, image, selem_cpy, + out=out, mask=mask, + shift_x=shift_x, shift_y=shift_y) def entropy(image, selem, out=None, mask=None, shift_x=False, shift_y=False): @@ -821,8 +840,9 @@ def entropy(image, selem, out=None, mask=None, shift_x=False, shift_y=False): """ return _apply_scalar_per_pixel(generic_cy._entropy, image, selem, - out=out, mask=mask, shift_x=shift_x, shift_y=shift_y, - out_dtype=np.double) + out=out, mask=mask, + shift_x=shift_x, shift_y=shift_y, + out_dtype=np.double) def otsu(image, selem, out=None, mask=None, shift_x=False, shift_y=False): @@ -865,10 +885,12 @@ def otsu(image, selem, out=None, mask=None, shift_x=False, shift_y=False): """ return _apply_scalar_per_pixel(generic_cy._otsu, image, selem, out=out, - mask=mask, shift_x=shift_x, shift_y=shift_y) + mask=mask, shift_x=shift_x, + shift_y=shift_y) -def windowed_histogram(image, selem, out=None, mask=None, shift_x=False, shift_y=False, n_bins=None): +def windowed_histogram(image, selem, out=None, mask=None, + shift_x=False, shift_y=False, n_bins=None): """Normalized sliding window histogram Parameters @@ -887,18 +909,18 @@ def windowed_histogram(image, selem, out=None, mask=None, shift_x=False, shift_y to the structuring element sizes (center must be inside the given structuring element). n_bins : int or None - The number of histogram bins. Will default to `image.max() + 1` + The number of histogram bins. Will default to ``image.max() + 1`` if None is passed. Returns ------- out : 3-D array with float dtype of dimensions (H,W,N), where (H,W) are - the dimensions of the input image and N is n_bins or image.max()+1 - if no value is provided as a parameter. Effectively, each pixel - is a N-D feature vector that is the histogram. The sum of the - elements in the feature vector will be 1, unless no pixels in the - window were covered by both selem and mask, in which case all - elements will be 0. + the dimensions of the input image and N is n_bins or + ``image.max() + 1`` if no value is provided as a parameter. + Effectively, each pixel is a N-D feature vector that is the histogram. + The sum of the elements in the feature vector will be 1, unless no + pixels in the window were covered by both selem and mask, in which + case all elements will be 0. Examples -------- diff --git a/skimage/filter/rank/generic_cy.pyx b/skimage/filter/rank/generic_cy.pyx index bc54855f..98c8cd8a 100644 --- a/skimage/filter/rank/generic_cy.pyx +++ b/skimage/filter/rank/generic_cy.pyx @@ -407,7 +407,6 @@ cdef inline void _kernel_win_hist(dtype_t_out* out, Py_ssize_t odepth, out[i] = 0 - def _autolevel(dtype_t[:, ::1] image, char[:, ::1] selem, char[:, ::1] mask, @@ -527,6 +526,7 @@ def _pop(dtype_t[:, ::1] image, _core(_kernel_pop[dtype_t_out, dtype_t], image, selem, mask, out, shift_x, shift_y, 0, 0, 0, 0, max_bin) + def _sum(dtype_t[:, ::1] image, char[:, ::1] selem, char[:, ::1] mask, diff --git a/skimage/filter/rank/percentile_cy.pyx b/skimage/filter/rank/percentile_cy.pyx index 6d54999f..5ff584d3 100644 --- a/skimage/filter/rank/percentile_cy.pyx +++ b/skimage/filter/rank/percentile_cy.pyx @@ -285,6 +285,7 @@ def _mean(dtype_t[:, ::1] image, _core(_kernel_mean[dtype_t_out, dtype_t], image, selem, mask, out, shift_x, shift_y, p0, p1, 0, 0, max_bin) + def _sum(dtype_t[:, ::1] image, char[:, ::1] selem, char[:, ::1] mask, @@ -295,6 +296,7 @@ def _sum(dtype_t[:, ::1] image, _core(_kernel_sum[dtype_t_out, dtype_t], image, selem, mask, out, shift_x, shift_y, p0, p1, 0, 0, max_bin) + def _subtract_mean(dtype_t[:, ::1] image, char[:, ::1] selem, char[:, ::1] mask, From c518f78ff72e6b4a65e3493d42db610ad615c3b9 Mon Sep 17 00:00:00 2001 From: Geoffrey French Date: Tue, 2 Sep 2014 20:17:43 +0100 Subject: [PATCH 13/13] Fixed mistake in generic.py --- skimage/filter/rank/generic.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/skimage/filter/rank/generic.py b/skimage/filter/rank/generic.py index f2e75051..a474d050 100644 --- a/skimage/filter/rank/generic.py +++ b/skimage/filter/rank/generic.py @@ -338,8 +338,7 @@ def mean(image, selem, out=None, mask=None, shift_x=False, shift_y=False): """ return _apply_scalar_per_pixel(generic_cy._mean, image, selem, out=out, - out=out, mask=mask, - shift_x=shift_x, shift_y=shift_y) + mask=mask, shift_x=shift_x, shift_y=shift_y) def subtract_mean(image, selem, out=None, mask=None, shift_x=False,