Refactor rank filter package as combined implementation for 8- and 16-bit

This commit is contained in:
Johannes Schönberger
2013-07-12 23:16:49 +02:00
parent 9a36e3b270
commit 4de4053a9f
18 changed files with 1279 additions and 2273 deletions
+26 -67
View File
@@ -3,9 +3,6 @@
The local histogram is computed using a sliding window similar to the method
described in [1]_.
Input image must be 16-bit, the number of histogram bins is determined from the
maximum value present in the image.
The pixel neighborhood is defined by:
* the given structuring element
@@ -14,8 +11,6 @@ The pixel neighborhood is defined by:
The kernel is flat (i.e. each pixel belonging to the neighborhood contributes
equally).
Result image is 16-bit with respect to the input image.
References
----------
@@ -30,46 +25,19 @@ from skimage import img_as_ubyte
from ... import get_log
log = get_log()
from . import bilateral16_cy
from .generic import find_bitdepth
from . import bilateral_cy
from .generic import _handle_input
__all__ = ['bilateral_mean', 'bilateral_pop']
def _apply(func8, func16, image, selem, out, mask, shift_x, shift_y, s0, s1):
selem = img_as_ubyte(selem > 0)
image = np.ascontiguousarray(image)
def _apply(func, image, selem, out, mask, shift_x, shift_y, s0, s1):
if mask is None:
mask = np.ones(image.shape, dtype=np.uint8)
else:
mask = np.ascontiguousarray(mask)
mask = img_as_ubyte(mask)
image, selem, out, mask, max_bin = _handle_input(image, selem, out, mask)
if image is out:
raise NotImplementedError("Cannot perform rank operation in place.")
if image.dtype == np.uint8:
if func8 is None:
raise TypeError("Not implemented for uint8 image.")
if out is None:
out = np.zeros(image.shape, dtype=np.uint8)
func8(image, selem, shift_x=shift_x, shift_y=shift_y,
mask=mask, out=out, s0=s0, s1=s1)
elif image.dtype == np.uint16:
if func16 is None:
raise TypeError("Not implemented for uint16 image.")
if out is None:
out = np.zeros(image.shape, dtype=np.uint16)
bitdepth = find_bitdepth(image)
if bitdepth > 10:
log.warn("Bitdepth of %d may result in bad rank filter "
"performance." % bitdepth)
func16(image, selem, shift_x=shift_x, shift_y=shift_y, mask=mask,
bitdepth=bitdepth + 1, out=out, s0=s0, s1=s1)
else:
raise TypeError("Only uint8 and uint16 image supported.")
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
@@ -91,16 +59,16 @@ def bilateral_mean(image, selem, out=None, mask=None, shift_x=False,
Parameters
----------
image : ndarray (uint16)
Input image.
image : ndarray (uint8, uint16)
Image array.
selem : ndarray
The neighborhood expressed as a 2-D array of 1's and 0's.
out : ndarray
out : ndarray (same dtype as input)
If None, a new array will be allocated.
mask : ndarray (uint8)
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)
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).
@@ -110,17 +78,12 @@ def bilateral_mean(image, selem, out=None, mask=None, shift_x=False,
Returns
-------
out : uint16 array
out : ndarray (same dtype as input)
The result of the local bilateral mean.
See also
--------
skimage.filter.denoise_bilateral() for a gaussian bilateral filter.
Notes
-----
* input image are 16-bit only
skimage.filter.denoise_bilateral for a gaussian bilateral filter.
Examples
--------
@@ -131,9 +94,10 @@ def bilateral_mean(image, selem, out=None, mask=None, shift_x=False,
>>> ima = data.camera().astype(np.uint16)
>>> # bilateral filtering of cameraman image using a flat kernel
>>> bilat_ima = bilateral_mean(ima, disk(20), s0=10,s1=10)
"""
return _apply(None, bilateral16_cy.mean, image, selem, out=out,
return _apply(bilateral_cy._mean, image, selem, out=out,
mask=mask, shift_x=shift_x, shift_y=shift_y, s0=s0, s1=s1)
@@ -145,16 +109,16 @@ def bilateral_pop(image, selem, out=None, mask=None, shift_x=False,
Parameters
----------
image : ndarray (uint16)
Input image.
image : ndarray (uint8, uint16)
Image array.
selem : ndarray
The neighborhood expressed as a 2-D array of 1's and 0's.
out : ndarray
out : ndarray (same dtype as input)
If None, a new array will be allocated.
mask : ndarray (uint8)
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)
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).
@@ -164,24 +128,19 @@ def bilateral_pop(image, selem, out=None, mask=None, shift_x=False,
Returns
-------
out : uint16 array
out : ndarray (same dtype as input)
the local number of pixels inside the bilateral neighborhood
Notes
-----
* input image are 16-bit only
Examples
--------
>>> # Local mean
>>> from skimage.morphology import square
>>> import skimage.filter.rank as rank
>>> ima16 = 255 * 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.uint16)
... [0, 1, 1, 1, 0],
... [0, 1, 1, 1, 0],
... [0, 1, 1, 1, 0],
... [0, 0, 0, 0, 0]], dtype=np.uint16)
>>> rank.bilateral_pop(ima16, square(3), s0=10,s1=10)
array([[3, 4, 3, 4, 3],
[4, 4, 6, 4, 4],
@@ -191,5 +150,5 @@ def bilateral_pop(image, selem, out=None, mask=None, shift_x=False,
"""
return _apply(None, bilateral16_cy.pop, image, selem, out=out,
return _apply(bilateral_cy._pop, image, selem, out=out,
mask=mask, shift_x=shift_x, shift_y=shift_y, s0=s0, s1=s1)
-79
View File
@@ -1,79 +0,0 @@
#cython: cdivision=True
#cython: boundscheck=False
#cython: nonecheck=False
#cython: wraparound=False
cimport numpy as cnp
from .core16_cy cimport dtype_t, _core16
# -----------------------------------------------------------------
# kernels uint16 take extra parameter for defining the bitdepth
# -----------------------------------------------------------------
cdef inline dtype_t kernel_mean(Py_ssize_t* histo, float pop,
dtype_t g, Py_ssize_t bitdepth,
Py_ssize_t maxbin, Py_ssize_t midbin,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef int i, bilat_pop = 0
cdef float mean = 0.
if pop:
for i in range(maxbin):
if (g > (i - s0)) and (g < (i + s1)):
bilat_pop += histo[i]
mean += histo[i] * i
if bilat_pop:
return <dtype_t>(mean / bilat_pop)
else:
return <dtype_t>(0)
else:
return <dtype_t>(0)
cdef inline dtype_t kernel_pop(Py_ssize_t* histo, float pop,
dtype_t g, Py_ssize_t bitdepth,
Py_ssize_t maxbin, Py_ssize_t midbin,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef int i, bilat_pop = 0
if pop:
for i in range(maxbin):
if (g > (i - s0)) and (g < (i + s1)):
bilat_pop += histo[i]
return <dtype_t>(bilat_pop)
else:
return <dtype_t>(0)
# -----------------------------------------------------------------
# python wrappers
# -----------------------------------------------------------------
def mean(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0, int bitdepth=8, int s0=1, int s1=1):
"""average greylevel (clipped on uint8)
"""
_core16(kernel_mean, image, selem, mask, out, shift_x, shift_y,
bitdepth, 0., 0., s0, s1)
def pop(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0, int bitdepth=8, int s0=1, int s1=1):
"""returns the number of actual pixels of the structuring element inside
the mask
"""
_core16(kernel_pop, image, selem, mask, out, shift_x, shift_y,
bitdepth, .0, .0, s0, s1)
+80
View File
@@ -0,0 +1,80 @@
#cython: cdivision=True
#cython: boundscheck=False
#cython: nonecheck=False
#cython: wraparound=False
cimport numpy as cnp
from libc.math cimport log
from .core_cy cimport uint8_t, uint16_t, dtype_t, _core
cdef inline dtype_t _kernel_mean(Py_ssize_t* histo, float pop,
dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t i
cdef Py_ssize_t bilat_pop = 0
cdef Py_ssize_t mean = 0
if pop:
for i in range(max_bin):
if (g > (i - s0)) and (g < (i + s1)):
bilat_pop += histo[i]
mean += histo[i] * i
if bilat_pop:
return <dtype_t>(mean / bilat_pop)
else:
return <dtype_t>(0)
else:
return <dtype_t>(0)
cdef inline dtype_t _kernel_pop(Py_ssize_t* histo, float pop,
dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t i
cdef Py_ssize_t bilat_pop = 0
if pop:
for i in range(max_bin):
if (g > (i - s0)) and (g < (i + s1)):
bilat_pop += histo[i]
return <dtype_t>(bilat_pop)
else:
return <dtype_t>(0)
def _mean(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask,
dtype_t[:, ::1] out,
char shift_x, char shift_y, Py_ssize_t s0, Py_ssize_t s1,
Py_ssize_t max_bin):
if dtype_t is uint8_t:
_core[uint8_t](_kernel_mean[uint8_t], image, selem, mask, out,
shift_x, shift_y, 0, 0, s0, s1, max_bin)
elif dtype_t is uint16_t:
_core[uint16_t](_kernel_mean[uint16_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[:, ::1] out,
char shift_x, char shift_y, Py_ssize_t s0, Py_ssize_t s1,
Py_ssize_t max_bin):
if dtype_t is uint8_t:
_core[uint8_t](_kernel_pop[uint8_t], image, selem, mask, out,
shift_x, shift_y, 0, 0, s0, s1, max_bin)
elif dtype_t is uint16_t:
_core[uint16_t](_kernel_pop[uint16_t], image, selem, mask, out,
shift_x, shift_y, 0, 0, s0, s1, max_bin)
-20
View File
@@ -1,20 +0,0 @@
cimport numpy as cnp
ctypedef cnp.uint16_t dtype_t
cdef dtype_t uint16_max(dtype_t a, dtype_t b)
cdef dtype_t uint16_min(dtype_t a, dtype_t b)
# 16-bit core kernel receives extra information about data bitdepth
cdef void _core16(dtype_t kernel(Py_ssize_t*, float, dtype_t,
Py_ssize_t, Py_ssize_t, Py_ssize_t, float,
float, Py_ssize_t, Py_ssize_t),
dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask,
dtype_t[:, ::1] out,
char shift_x, char shift_y, Py_ssize_t bitdepth,
float p0, float p1, Py_ssize_t s0, Py_ssize_t s1) except *
-247
View File
@@ -1,247 +0,0 @@
#cython: cdivision=True
#cython: boundscheck=False
#cython: nonecheck=False
#cython: wraparound=False
import numpy as np
cimport numpy as cnp
from libc.stdlib cimport malloc, free
from .core8_cy cimport is_in_mask
cdef inline dtype_t uint16_max(dtype_t a, dtype_t b):
return a if a >= b else b
cdef inline dtype_t uint16_min(dtype_t a, dtype_t b):
return a if a <= b else b
cdef inline void histogram_increment(Py_ssize_t* histo, float* pop,
dtype_t value):
histo[value] += 1
pop[0] += 1
cdef inline void histogram_decrement(Py_ssize_t* histo, float* pop,
dtype_t value):
histo[value] -= 1
pop[0] -= 1
cdef void _core16(dtype_t kernel(Py_ssize_t*, float, dtype_t,
Py_ssize_t, Py_ssize_t, Py_ssize_t, float,
float, Py_ssize_t, Py_ssize_t),
dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask,
dtype_t[:, ::1] out,
char shift_x, char shift_y, Py_ssize_t bitdepth,
float p0, float p1, Py_ssize_t s0, Py_ssize_t s1) except *:
"""Compute histogram for each pixel neighborhood, apply kernel function and
use kernel function return value for output image.
"""
cdef Py_ssize_t rows = image.shape[0]
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 centre_r = int(selem.shape[0] / 2) + shift_y
cdef Py_ssize_t centre_c = int(selem.shape[1] / 2) + shift_x
# check that structuring element center is inside the element bounding box
assert centre_r >= 0
assert centre_c >= 0
assert centre_r < srows
assert centre_c < scols
maxbin_list = [0, 0, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096,
8192, 16384, 32768, 65536]
midbin_list = [int(m / 2) for m in maxbin_list]
# set maxbin and midbin
cdef Py_ssize_t maxbin = maxbin_list[bitdepth]
cdef Py_ssize_t midbin = midbin_list[bitdepth]
# define pointers to the data
cdef char* mask_data = &mask[0, 0]
# define local variable types
cdef Py_ssize_t r, c, rr, cc, s, value, local_max, i, even_row
# number of pixels actually inside the neighborhood (float)
cdef float pop
# allocate memory with malloc
cdef Py_ssize_t max_se = srows * scols
# number of element in each attack border
cdef Py_ssize_t num_se_n, num_se_s, num_se_e, num_se_w
# the current local histogram distribution
cdef Py_ssize_t* histo = <Py_ssize_t*>malloc(maxbin * sizeof(Py_ssize_t))
# these lists contain the relative pixel row and column for each of the 4
# attack borders east, west, north and south e.g. se_e_r lists the rows of
# the east structuring element border
cdef Py_ssize_t* se_e_r = <Py_ssize_t*>malloc(max_se * sizeof(Py_ssize_t))
cdef Py_ssize_t* se_e_c = <Py_ssize_t*>malloc(max_se * sizeof(Py_ssize_t))
cdef Py_ssize_t* se_w_r = <Py_ssize_t*>malloc(max_se * sizeof(Py_ssize_t))
cdef Py_ssize_t* se_w_c = <Py_ssize_t*>malloc(max_se * sizeof(Py_ssize_t))
cdef Py_ssize_t* se_n_r = <Py_ssize_t*>malloc(max_se * sizeof(Py_ssize_t))
cdef Py_ssize_t* se_n_c = <Py_ssize_t*>malloc(max_se * sizeof(Py_ssize_t))
cdef Py_ssize_t* se_s_r = <Py_ssize_t*>malloc(max_se * sizeof(Py_ssize_t))
cdef Py_ssize_t* se_s_c = <Py_ssize_t*>malloc(max_se * sizeof(Py_ssize_t))
# build attack and release borders by using difference along axis
t = np.hstack((selem, np.zeros((selem.shape[0], 1))))
cdef char[:, :] t_e = (np.diff(t, axis=1) < 0).view(np.uint8)
t = np.hstack((np.zeros((selem.shape[0], 1)), selem))
cdef char[:, :] t_w = (np.diff(t, axis=1) > 0).view(np.uint8)
t = np.vstack((selem, np.zeros((1, selem.shape[1]))))
cdef char[:, :] t_s = (np.diff(t, axis=0) < 0).view(np.uint8)
t = np.vstack((np.zeros((1, selem.shape[1])), selem))
cdef char[:, :] t_n = (np.diff(t, axis=0) > 0).view(np.uint8)
num_se_n = num_se_s = num_se_e = num_se_w = 0
for r in range(srows):
for c in range(scols):
if t_e[r, c]:
se_e_r[num_se_e] = r - centre_r
se_e_c[num_se_e] = c - centre_c
num_se_e += 1
if t_w[r, c]:
se_w_r[num_se_w] = r - centre_r
se_w_c[num_se_w] = c - centre_c
num_se_w += 1
if t_n[r, c]:
se_n_r[num_se_n] = r - centre_r
se_n_c[num_se_n] = c - centre_c
num_se_n += 1
if t_s[r, c]:
se_s_r[num_se_s] = r - centre_r
se_s_c[num_se_s] = c - centre_c
num_se_s += 1
# initial population and histogram
for i in range(maxbin):
histo[i] = 0
pop = 0
for r in range(srows):
for c in range(scols):
rr = r - centre_r
cc = c - centre_c
if selem[r, c]:
if is_in_mask(rows, cols, rr, cc, mask_data):
histogram_increment(histo, &pop, image[rr, cc])
r = 0
c = 0
# kernel -------------------------------------------
out[r, c] = kernel(histo, pop, image[r, c], bitdepth, maxbin, midbin,
p0, p1, s0, s1)
# kernel -------------------------------------------
# main loop
r = 0
for even_row in range(0, rows, 2):
# ---> west to east
for c in range(1, cols):
for s in range(num_se_e):
rr = r + se_e_r[s]
cc = c + se_e_c[s]
if is_in_mask(rows, cols, rr, cc, mask_data):
histogram_increment(histo, &pop, image[rr, cc])
for s in range(num_se_w):
rr = r + se_w_r[s]
cc = c + se_w_c[s] - 1
if is_in_mask(rows, cols, rr, cc, mask_data):
histogram_decrement(histo, &pop, image[rr, cc])
# kernel -------------------------------------------
out[r, c] = kernel(histo, pop, image[r, c], bitdepth, maxbin,
midbin, p0, p1, s0, s1)
# kernel -------------------------------------------
r += 1 # pass to the next row
if r >= rows:
break
# ---> north to south
for s in range(num_se_s):
rr = r + se_s_r[s]
cc = c + se_s_c[s]
if is_in_mask(rows, cols, rr, cc, mask_data):
histogram_increment(histo, &pop, image[rr, cc])
for s in range(num_se_n):
rr = r + se_n_r[s] - 1
cc = c + se_n_c[s]
if is_in_mask(rows, cols, rr, cc, mask_data):
histogram_decrement(histo, &pop, image[rr, cc])
# kernel -------------------------------------------
out[r, c] = kernel(histo, pop, image[r, c],
bitdepth, maxbin, midbin, p0, p1, s0, s1)
# kernel -------------------------------------------
# ---> east to west
for c in range(cols - 2, -1, -1):
for s in range(num_se_w):
rr = r + se_w_r[s]
cc = c + se_w_c[s]
if is_in_mask(rows, cols, rr, cc, mask_data):
histogram_increment(histo, &pop, image[rr, cc])
for s in range(num_se_e):
rr = r + se_e_r[s]
cc = c + se_e_c[s] + 1
if is_in_mask(rows, cols, rr, cc, mask_data):
histogram_decrement(histo, &pop, image[rr, cc])
# kernel -------------------------------------------
out[r, c] = kernel(histo, pop, image[r, c], bitdepth, maxbin,
midbin, p0, p1, s0, s1)
# kernel -------------------------------------------
r += 1 # pass to the next row
if r >= rows:
break
# ---> north to south
for s in range(num_se_s):
rr = r + se_s_r[s]
cc = c + se_s_c[s]
if is_in_mask(rows, cols, rr, cc, mask_data):
histogram_increment(histo, &pop, image[rr, cc])
for s in range(num_se_n):
rr = r + se_n_r[s] - 1
cc = c + se_n_c[s]
if is_in_mask(rows, cols, rr, cc, mask_data):
histogram_decrement(histo, &pop, image[rr, cc])
# kernel -------------------------------------------
out[r, c] = kernel(histo, pop, image[r, c], bitdepth, maxbin, midbin,
p0, p1, s0, s1)
# kernel -------------------------------------------
# release memory allocated by malloc
free(se_e_r)
free(se_e_c)
free(se_w_r)
free(se_w_c)
free(se_n_r)
free(se_n_c)
free(se_s_r)
free(se_s_c)
free(histo)
-25
View File
@@ -1,25 +0,0 @@
cimport numpy as cnp
ctypedef cnp.uint8_t dtype_t
cdef dtype_t uint8_max(dtype_t a, dtype_t b)
cdef dtype_t uint8_min(dtype_t a, dtype_t b)
cdef dtype_t is_in_mask(Py_ssize_t rows, Py_ssize_t cols,
Py_ssize_t r, Py_ssize_t c,
char* mask)
# 8-bit core kernel receives extra information about data inferior and superior
# percentiles
cdef void _core8(dtype_t kernel(Py_ssize_t*, float, dtype_t, float,
float, Py_ssize_t, Py_ssize_t),
dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask,
dtype_t[:, ::1] out,
char shift_x, char shift_y, float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1) except *
+23
View File
@@ -0,0 +1,23 @@
from numpy cimport uint8_t, uint16_t
ctypedef fused dtype_t:
uint8_t
uint16_t
cdef dtype_t _max(dtype_t a, dtype_t b)
cdef dtype_t _min(dtype_t a, dtype_t b)
cdef void _core(dtype_t kernel(Py_ssize_t*, float, dtype_t,
Py_ssize_t, Py_ssize_t, float,
float, Py_ssize_t, Py_ssize_t),
dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask,
dtype_t[:, ::1] out,
char shift_x, char shift_y,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1,
Py_ssize_t max_bin) except *
@@ -9,11 +9,11 @@ cimport numpy as cnp
from libc.stdlib cimport malloc, free
cdef inline dtype_t uint8_max(dtype_t a, dtype_t b):
cdef inline dtype_t _max(dtype_t a, dtype_t b):
return a if a >= b else b
cdef inline dtype_t uint8_min(dtype_t a, dtype_t b):
cdef inline dtype_t _min(dtype_t a, dtype_t b):
return a if a <= b else b
@@ -29,9 +29,9 @@ cdef inline void histogram_decrement(Py_ssize_t* histo, float* pop,
pop[0] -= 1
cdef inline dtype_t is_in_mask(Py_ssize_t rows, Py_ssize_t cols,
Py_ssize_t r, Py_ssize_t c,
char* mask):
cdef inline char is_in_mask(Py_ssize_t rows, Py_ssize_t cols,
Py_ssize_t r, Py_ssize_t c,
char* mask):
"""Check whether given coordinate is within image and mask is true."""
if r < 0 or r > rows - 1 or c < 0 or c > cols - 1:
return 0
@@ -42,14 +42,17 @@ cdef inline dtype_t is_in_mask(Py_ssize_t rows, Py_ssize_t cols,
return 0
cdef void _core8(dtype_t kernel(Py_ssize_t*, float, dtype_t, float,
float, Py_ssize_t, Py_ssize_t),
dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask,
dtype_t[:, ::1] out,
char shift_x, char shift_y, float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1) except *:
cdef void _core(dtype_t kernel(Py_ssize_t*, float, dtype_t,
Py_ssize_t, Py_ssize_t, float,
float, Py_ssize_t, Py_ssize_t),
dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask,
dtype_t[:, ::1] out,
char shift_x, char shift_y,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1,
Py_ssize_t max_bin) except *:
"""Compute histogram for each pixel neighborhood, apply kernel function and
use kernel function return value for output image.
"""
@@ -59,8 +62,8 @@ cdef void _core8(dtype_t kernel(Py_ssize_t*, float, dtype_t, float,
cdef Py_ssize_t srows = selem.shape[0]
cdef Py_ssize_t scols = selem.shape[1]
cdef Py_ssize_t centre_r = int(selem.shape[0] / 2) + shift_y
cdef Py_ssize_t centre_c = int(selem.shape[1] / 2) + shift_x
cdef Py_ssize_t centre_r = <Py_ssize_t>(selem.shape[0] / 2) + shift_y
cdef Py_ssize_t centre_c = <Py_ssize_t>(selem.shape[1] / 2) + shift_x
# check that structuring element center is inside the element bounding box
assert centre_r >= 0
@@ -68,26 +71,35 @@ cdef void _core8(dtype_t kernel(Py_ssize_t*, float, dtype_t, float,
assert centre_r < srows
assert centre_c < scols
# add 1 to ensure maximum value is included in histogram -> range(max_bin)
max_bin += 1
cdef Py_ssize_t mid_bin = max_bin / 2
# define pointers to the data
cdef char* mask_data = &mask[0, 0]
# define local variable types
cdef Py_ssize_t r, c, rr, cc, s, value, local_max, i, even_row
# number of pixels actually inside the neighborhood (float)
cdef float pop
# allocate memory with malloc
cdef Py_ssize_t max_se = srows * scols
# number of element in each attack border
cdef Py_ssize_t num_se_n, num_se_s, num_se_e, num_se_w
cdef float pop = 0
# the current local histogram distribution
cdef Py_ssize_t* histo = <Py_ssize_t*>malloc(256 * sizeof(Py_ssize_t))
cdef Py_ssize_t* histo = <Py_ssize_t*>malloc(max_bin * sizeof(Py_ssize_t))
for i in range(max_bin):
histo[i] = 0
# these lists contain the relative pixel row and column for each of the 4
# attack borders east, west, north and south e.g. se_e_r lists the rows of
# the east structuring element border
cdef Py_ssize_t max_se = srows * scols
# number of element in each attack border
cdef Py_ssize_t num_se_n, num_se_s, num_se_e, num_se_w
num_se_n = num_se_s = num_se_e = num_se_w = 0
cdef Py_ssize_t* se_e_r = <Py_ssize_t*>malloc(max_se * sizeof(Py_ssize_t))
cdef Py_ssize_t* se_e_c = <Py_ssize_t*>malloc(max_se * sizeof(Py_ssize_t))
cdef Py_ssize_t* se_w_r = <Py_ssize_t*>malloc(max_se * sizeof(Py_ssize_t))
@@ -110,8 +122,6 @@ cdef void _core8(dtype_t kernel(Py_ssize_t*, float, dtype_t, float,
t = np.vstack((np.zeros((1, selem.shape[1])), selem))
cdef char[:, :] t_n = (np.diff(t, axis=0) > 0).view(np.uint8)
num_se_n = num_se_s = num_se_e = num_se_w = 0
for r in range(srows):
for c in range(scols):
if t_e[r, c]:
@@ -131,13 +141,6 @@ cdef void _core8(dtype_t kernel(Py_ssize_t*, float, dtype_t, float,
se_s_c[num_se_s] = c - centre_c
num_se_s += 1
# initial population and histogram (kernel is centered on the first row and
# column)
for i in range(256):
histo[i] = 0
pop = 0
for r in range(srows):
for c in range(scols):
rr = r - centre_r
@@ -148,13 +151,13 @@ cdef void _core8(dtype_t kernel(Py_ssize_t*, float, dtype_t, float,
r = 0
c = 0
# kernel ------------------------------------------------------------------
out[r, c] = kernel(histo, pop, image[r, c], p0, p1, s0, s1)
# kernel ------------------------------------------------------------------
out[r, c] = kernel(histo, pop, image[r, c], max_bin, mid_bin,
p0, p1, s0, s1)
# main loop
r = 0
for even_row in range(0, rows, 2):
# ---> west to east
for c in range(1, cols):
for s in range(num_se_e):
@@ -169,9 +172,8 @@ cdef void _core8(dtype_t kernel(Py_ssize_t*, float, dtype_t, float,
if is_in_mask(rows, cols, rr, cc, mask_data):
histogram_decrement(histo, &pop, image[rr, cc])
# kernel ----------------------------------------------------------
out[r, c] = kernel(histo, pop, image[r, c], p0, p1, s0, s1)
# kernel ----------------------------------------------------------
out[r, c] = kernel(histo, pop, image[r, c], max_bin,
mid_bin, p0, p1, s0, s1)
r += 1 # pass to the next row
if r >= rows:
@@ -190,9 +192,8 @@ cdef void _core8(dtype_t kernel(Py_ssize_t*, float, dtype_t, float,
if is_in_mask(rows, cols, rr, cc, mask_data):
histogram_decrement(histo, &pop, image[rr, cc])
# kernel --------------------------------------------------------------
out[r, c] = kernel(histo, pop, image[r, c], p0, p1, s0, s1)
# kernel --------------------------------------------------------------
out[r, c] = kernel(histo, pop, image[r, c],
max_bin, mid_bin, p0, p1, s0, s1)
# ---> east to west
for c in range(cols - 2, -1, -1):
@@ -208,10 +209,8 @@ cdef void _core8(dtype_t kernel(Py_ssize_t*, float, dtype_t, float,
if is_in_mask(rows, cols, rr, cc, mask_data):
histogram_decrement(histo, &pop, image[rr, cc])
# kernel ----------------------------------------------------------
out[r, c] = kernel(
histo, pop, image[r, c], p0, p1, s0, s1)
# kernel ----------------------------------------------------------
out[r, c] = kernel(histo, pop, image[r, c], max_bin,
mid_bin, p0, p1, s0, s1)
r += 1 # pass to the next row
if r >= rows:
@@ -230,9 +229,8 @@ cdef void _core8(dtype_t kernel(Py_ssize_t*, float, dtype_t, float,
if is_in_mask(rows, cols, rr, cc, mask_data):
histogram_decrement(histo, &pop, image[rr, cc])
# kernel --------------------------------------------------------------
out[r, c] = kernel(histo, pop, image[r, c], p0, p1, s0, s1)
# kernel --------------------------------------------------------------
out[r, c] = kernel(histo, pop, image[r, c], max_bin, mid_bin,
p0, p1, s0, s1)
# release memory allocated by malloc
free(se_e_r)
@@ -243,5 +241,4 @@ cdef void _core8(dtype_t kernel(Py_ssize_t*, float, dtype_t, float,
free(se_n_c)
free(se_s_r)
free(se_s_c)
free(histo)
+123 -149
View File
@@ -20,7 +20,7 @@ from skimage import img_as_ubyte, img_as_uint
from ... import get_log
log = get_log()
from . import generic8_cy, generic16_cy
from . import generic_cy
__all__ = ['autolevel', 'bottomhat', 'equalize', 'gradient', 'maximum', 'mean',
@@ -28,56 +28,43 @@ __all__ = ['autolevel', 'bottomhat', 'equalize', 'gradient', 'maximum', 'mean',
'pop', 'threshold', 'tophat', 'noise_filter', 'entropy', 'otsu']
import numpy as np
def _handle_input(image, selem, out, mask):
if image.dtype not in (np.uint8, np.uint16):
image = img_as_ubyte(image)
def find_bitdepth(image):
"""returns the max bith depth of a uint16 image
"""
umax = np.max(image)
if umax > 2:
return int(np.log2(umax))
else:
return 1
def _apply(func8, func16, image, selem, out, mask, shift_x, shift_y):
selem = img_as_ubyte(selem > 0)
selem = np.ascontiguousarray(img_as_ubyte(selem > 0))
image = np.ascontiguousarray(image)
if mask is None:
mask = np.ones(image.shape, dtype=np.uint8)
else:
mask = np.ascontiguousarray(mask)
mask = img_as_ubyte(mask)
mask = np.ascontiguousarray(mask)
if out is None:
out = np.empty_like(image, dtype=image.dtype)
if image is out:
raise NotImplementedError("Cannot perform rank operation in place.")
is_8bit = image.dtype in (np.uint8, np.int8)
if func8 is not None and (is_8bit or func16 is None):
out = _apply8(func8, image, selem, out, mask, shift_x, shift_y)
if is_8bit:
max_bin = 255
else:
image = img_as_uint(image)
if out is None:
out = np.zeros(image.shape, dtype=np.uint16)
bitdepth = find_bitdepth(image)
if bitdepth > 10:
log.warn("Bitdepth of %d may result in bad rank filter "
"performance." % bitdepth)
func16(image, selem, shift_x=shift_x, shift_y=shift_y, mask=mask,
bitdepth=bitdepth + 1, out=out)
max_bin = max(4, image.max())
return out
return image, selem, out, mask, max_bin
def _apply8(func8, image, selem, out, mask, shift_x, shift_y):
if out is None:
out = np.zeros(image.shape, dtype=np.uint8)
image = img_as_ubyte(image)
func8(image, selem, shift_x=shift_x, shift_y=shift_y,
mask=mask, out=out)
def _apply(func, image, selem, out, mask, shift_x, shift_y):
image, selem, out, mask, max_bin = _handle_input(image, selem, out, mask)
func(image, selem, shift_x=shift_x, shift_y=shift_y, mask=mask,
out=out, max_bin=max_bin)
return out
@@ -86,13 +73,13 @@ def autolevel(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
Parameters
----------
image : ndarray
Image array (uint8 array or uint16).
image : ndarray (uint8, uint16)
Image array.
selem : ndarray
The neighborhood expressed as a 2-D array of 1's and 0's.
out : ndarray
out : ndarray (same dtype as input)
If None, a new array will be allocated.
mask : ndarray (uint8)
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
@@ -102,7 +89,7 @@ def autolevel(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
Returns
-------
out : uint8 array or uint16 array (same as input image)
out : ndarray (same dtype as input image)
The result of the local autolevel.
Examples
@@ -117,7 +104,7 @@ def autolevel(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
"""
return _apply(generic8_cy.autolevel, generic16_cy.autolevel, image, selem,
return _apply(generic_cy._autolevel, image, selem,
out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
@@ -126,13 +113,13 @@ def bottomhat(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
Parameters
----------
image : ndarray
Image array (uint8 array or uint16).
image : ndarray (uint8, uint16)
Image array.
selem : ndarray
The neighborhood expressed as a 2-D array of 1's and 0's.
out : ndarray
out : ndarray (same dtype as input)
If None, a new array will be allocated.
mask : ndarray (uint8)
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
@@ -142,12 +129,12 @@ def bottomhat(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
Returns
-------
local bottomhat : uint8 array or uint16 array depending on input image
bottomhat : ndarray (same dtype as input image)
The result of the local bottomhat.
"""
return _apply(generic8_cy.bottomhat, generic16_cy.bottomhat, image, selem,
return _apply(generic_cy._bottomhat, image, selem,
out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
@@ -156,13 +143,13 @@ def equalize(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
Parameters
----------
image : ndarray
Image array (uint8 array or uint16).
image : ndarray (uint8, uint16)
Image array.
selem : ndarray
The neighborhood expressed as a 2-D array of 1's and 0's.
out : ndarray
out : ndarray (same dtype as input)
If None, a new array will be allocated.
mask : ndarray (uint8)
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
@@ -172,7 +159,7 @@ def equalize(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
Returns
-------
out : uint8 array or uint16 array (same as input image)
out : ndarray (same dtype as input image)
The result of the local equalize.
Examples
@@ -187,7 +174,7 @@ def equalize(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
"""
return _apply(generic8_cy.equalize, generic16_cy.equalize, image, selem,
return _apply(generic_cy._equalize, image, selem,
out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
@@ -198,13 +185,13 @@ def gradient(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
Parameters
----------
image : ndarray
Image array (uint8 array or uint16).
image : ndarray (uint8, uint16)
Image array.
selem : ndarray
The neighborhood expressed as a 2-D array of 1's and 0's.
out : ndarray
out : ndarray (same dtype as input)
If None, a new array will be allocated.
mask : ndarray (uint8)
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
@@ -214,12 +201,12 @@ def gradient(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
Returns
-------
out : uint8 array or uint16 array (same as input image)
out : ndarray (same dtype as input image)
The local gradient.
"""
return _apply(generic8_cy.gradient, generic16_cy.gradient, image, selem,
return _apply(generic_cy._gradient, image, selem,
out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
@@ -229,13 +216,13 @@ def maximum(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
Parameters
----------
image : ndarray
Image array (uint8 array or uint16).
image : ndarray (uint8, uint16)
Image array.
selem : ndarray
The neighborhood expressed as a 2-D array of 1's and 0's.
out : ndarray
out : ndarray (same dtype as input)
If None, a new array will be allocated.
mask : ndarray (uint8)
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
@@ -245,7 +232,7 @@ def maximum(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
Returns
-------
out : uint8 array or uint16 array (same as input image)
out : ndarray (same dtype as input image)
The local maximum.
See also
@@ -254,13 +241,12 @@ def maximum(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
Note
----
* input image can be 8-bit or 16-bit with a value < 4096 (i.e. 12 bit)
* the lower algorithm complexity makes the rank.maximum() more efficient for
larger images and structuring elements
"""
return _apply(generic8_cy.maximum, generic16_cy.maximum, image, selem,
return _apply(generic_cy._maximum, image, selem,
out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
@@ -269,13 +255,13 @@ def mean(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
Parameters
----------
image : ndarray
Image array (uint8 array or uint16).
image : ndarray (uint8, uint16)
Image array.
selem : ndarray
The neighborhood expressed as a 2-D array of 1's and 0's.
out : ndarray
out : ndarray (same dtype as input)
If None, a new array will be allocated.
mask : ndarray (uint8)
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
@@ -285,7 +271,7 @@ def mean(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
Returns
-------
out : uint8 array or uint16 array (same as input image)
out : ndarray (same dtype as input image)
The local mean.
Examples
@@ -300,7 +286,7 @@ def mean(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
"""
return _apply(generic8_cy.mean, generic16_cy.mean, image, selem, out=out,
return _apply(generic_cy._mean, image, selem, out=out,
mask=mask, shift_x=shift_x, shift_y=shift_y)
@@ -310,13 +296,13 @@ def meansubtraction(image, selem, out=None, mask=None, shift_x=False,
Parameters
----------
image : ndarray
Image array (uint8 array or uint16).
image : ndarray (uint8, uint16)
Image array.
selem : ndarray
The neighborhood expressed as a 2-D array of 1's and 0's.
out : ndarray
out : ndarray (same dtype as input)
If None, a new array will be allocated.
mask : ndarray (uint8)
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
@@ -326,14 +312,13 @@ def meansubtraction(image, selem, out=None, mask=None, shift_x=False,
Returns
-------
out : uint8 array or uint16 array (same as input image)
out : ndarray (same dtype as input image)
The result of the local meansubtraction.
"""
return _apply(generic8_cy.meansubtraction, generic16_cy.meansubtraction,
image, selem, out=out, mask=mask, shift_x=shift_x,
shift_y=shift_y)
return _apply(generic_cy._meansubtraction, image, selem,
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):
@@ -341,13 +326,13 @@ def median(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
Parameters
----------
image : ndarray
Image array (uint8 array or uint16).
image : ndarray (uint8, uint16)
Image array.
selem : ndarray
The neighborhood expressed as a 2-D array of 1's and 0's.
out : ndarray
out : ndarray (same dtype as input)
If None, a new array will be allocated.
mask : ndarray (uint8)
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
@@ -357,7 +342,7 @@ def median(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
Returns
-------
out : uint8 array or uint16 array (same as input image)
out : ndarray (same dtype as input image)
The local median.
Examples
@@ -372,7 +357,7 @@ def median(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
"""
return _apply(generic8_cy.median, generic16_cy.median, image, selem,
return _apply(generic_cy._median, image, selem,
out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
@@ -381,13 +366,13 @@ def minimum(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
Parameters
----------
image : ndarray
Image array (uint8 array or uint16).
image : ndarray (uint8, uint16)
Image array.
selem : ndarray
The neighborhood expressed as a 2-D array of 1's and 0's.
out : ndarray
out : ndarray (same dtype as input)
If None, a new array will be allocated.
mask : ndarray (uint8)
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
@@ -397,7 +382,7 @@ def minimum(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
Returns
-------
out : uint8 array or uint16 array (same as input image)
out : ndarray (same dtype as input image)
The local minimum.
See also
@@ -406,13 +391,12 @@ def minimum(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
Note
----
* input image can be 8-bit or 16-bit with a value < 4096 (i.e. 12 bit)
* the lower algorithm complexity makes the rank.minimum() more efficient
for larger images and structuring elements
"""
return _apply(generic8_cy.minimum, generic16_cy.minimum, image, selem,
return _apply(generic_cy._minimum, image, selem,
out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
@@ -421,13 +405,13 @@ def modal(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
Parameters
----------
image : ndarray
Image array (uint8 array or uint16).
image : ndarray (uint8, uint16)
Image array.
selem : ndarray
The neighborhood expressed as a 2-D array of 1's and 0's.
out : ndarray
out : ndarray (same dtype as input)
If None, a new array will be allocated.
mask : ndarray (uint8)
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
@@ -437,12 +421,12 @@ def modal(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
Returns
-------
out : uint8 array or uint16 array (same as input image)
out : ndarray (same dtype as input image)
The local modal.
"""
return _apply(generic8_cy.modal, generic16_cy.modal, image, selem,
return _apply(generic_cy._modal, image, selem,
out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
@@ -454,13 +438,13 @@ def morph_contr_enh(image, selem, out=None, mask=None, shift_x=False,
Parameters
----------
image : ndarray
Image array (uint8 array or uint16).
image : ndarray (uint8, uint16)
Image array.
selem : ndarray
The neighborhood expressed as a 2-D array of 1's and 0's.
out : ndarray
out : ndarray (same dtype as input)
If None, a new array will be allocated.
mask : ndarray (uint8)
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
@@ -470,7 +454,7 @@ def morph_contr_enh(image, selem, out=None, mask=None, shift_x=False,
Returns
-------
out : uint8 array or uint16 array (same as input image)
out : ndarray (same dtype as input image)
The result of the local morph_contr_enh.
Examples
@@ -485,9 +469,8 @@ def morph_contr_enh(image, selem, out=None, mask=None, shift_x=False,
"""
return _apply(generic8_cy.morph_contr_enh, generic16_cy.morph_contr_enh,
image, selem, out=out, mask=mask, shift_x=shift_x,
shift_y=shift_y)
return _apply(generic_cy._morph_contr_enh, image, selem,
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):
@@ -496,13 +479,13 @@ def pop(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
Parameters
----------
image : ndarray
Image array (uint8 array or uint16).
image : ndarray (uint8, uint16)
Image array.
selem : ndarray
The neighborhood expressed as a 2-D array of 1's and 0's.
out : ndarray
out : ndarray (same dtype as input)
If None, a new array will be allocated.
mask : ndarray (uint8)
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
@@ -512,7 +495,7 @@ def pop(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
Returns
-------
out : uint8 array or uint16 array (same as input image)
out : ndarray (same dtype as input image)
The number of pixels belonging to the neighborhood.
Examples
@@ -534,7 +517,7 @@ def pop(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
"""
return _apply(generic8_cy.pop, generic16_cy.pop, image, selem, out=out,
return _apply(generic_cy._pop, image, selem, out=out,
mask=mask, shift_x=shift_x, shift_y=shift_y)
@@ -543,13 +526,13 @@ def threshold(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
Parameters
----------
image : ndarray
Image array (uint8 array or uint16).
image : ndarray (uint8, uint16)
Image array.
selem : ndarray
The neighborhood expressed as a 2-D array of 1's and 0's.
out : ndarray
out : ndarray (same dtype as input)
If None, a new array will be allocated.
mask : ndarray (uint8)
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
@@ -559,7 +542,7 @@ def threshold(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
Returns
-------
out : uint8 array or uint16 array (same as input image)
out : ndarray (same dtype as input image)
The result of the local threshold.
Examples
@@ -581,7 +564,7 @@ def threshold(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
"""
return _apply(generic8_cy.threshold, generic16_cy.threshold, image, selem,
return _apply(generic_cy._threshold, image, selem,
out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
@@ -590,13 +573,13 @@ def tophat(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
Parameters
----------
image : ndarray
Image array (uint8 array or uint16).
image : ndarray (uint8, uint16)
Image array.
selem : ndarray
The neighborhood expressed as a 2-D array of 1's and 0's.
out : ndarray
out : ndarray (same dtype as input)
If None, a new array will be allocated.
mask : ndarray (uint8)
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
@@ -606,12 +589,12 @@ def tophat(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
Returns
-------
out : uint8 array or uint16 array (same as input image)
out : ndarray (same dtype as input image)
The image tophat.
"""
return _apply(generic8_cy.tophat, generic16_cy.tophat, image, selem,
return _apply(generic_cy._tophat, image, selem,
out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
@@ -621,13 +604,13 @@ def noise_filter(image, selem, out=None, mask=None, shift_x=False,
Parameters
----------
image : ndarray
Image array (uint8 array or uint16).
image : ndarray (uint8, uint16)
Image array.
selem : ndarray
The neighborhood expressed as a 2-D array of 1's and 0's.
out : ndarray
out : ndarray (same dtype as input)
If None, a new array will be allocated.
mask : ndarray (uint8)
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
@@ -642,7 +625,7 @@ def noise_filter(image, selem, out=None, mask=None, shift_x=False,
Returns
-------
out : uint8 array or uint16 array (same as input image)
out : ndarray (same dtype as input image)
The image noise.
"""
@@ -654,7 +637,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(generic8_cy.noise_filter, None, image, selem_cpy, out=out,
return _apply(generic_cy._noise_filter, image, selem_cpy, out=out,
mask=mask, shift_x=shift_x, shift_y=shift_y)
@@ -665,13 +648,13 @@ def entropy(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
Parameters
----------
image : ndarray
Image array (uint8 array or uint16).
image : ndarray (uint8, uint16)
Image array.
selem : ndarray
The neighborhood expressed as a 2-D array of 1's and 0's.
out : ndarray
out : ndarray (same dtype as input)
If None, a new array will be allocated.
mask : ndarray (uint8)
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
@@ -681,8 +664,8 @@ def entropy(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
Returns
-------
out : uint8 array or uint16 array (same as input image)
Entropy x10 (uint8 images) and entropy x1000 (uint16 images)
out : ndarray (same dtype as input image)
Entropy of image.
References
----------
@@ -694,17 +677,12 @@ def entropy(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
>>> from skimage import data
>>> from skimage.filter.rank import entropy
>>> from skimage.morphology import disk
>>> # defining a 8- and a 16-bit test images
>>> a8 = data.camera()
>>> a16 = data.camera().astype(np.uint16) * 4
>>> # pixel values contain 10x the local entropy
>>> ent8 = entropy(a8, disk(5))
>>> # pixel values contain 1000x the local entropy
>>> ent16 = entropy(a16, disk(5))
"""
return _apply(generic8_cy.entropy, generic16_cy.entropy, image, selem,
return _apply(generic_cy._entropy, image, selem,
out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
@@ -719,7 +697,7 @@ def otsu(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
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 (uint8)
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
@@ -729,17 +707,13 @@ def otsu(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
Returns
-------
out : uint8 array
Otsu's threshold values
out : ndarray (same dtype as input image)
Otsu's threshold values.
References
----------
.. [otsu] http://en.wikipedia.org/wiki/Otsu's_method
Notes
-----
* input image are 8-bit only
Examples
--------
>>> # Local entropy
@@ -753,5 +727,5 @@ def otsu(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
"""
return _apply(generic8_cy.otsu, None, image, selem, out=out,
return _apply(generic_cy._otsu, image, selem, out=out,
mask=mask, shift_x=shift_x, shift_y=shift_y)
-418
View File
@@ -1,418 +0,0 @@
#cython: cdivision=True
#cython: boundscheck=False
#cython: nonecheck=False
#cython: wraparound=False
cimport numpy as cnp
from libc.math cimport log
from .core16_cy cimport dtype_t, _core16
# -----------------------------------------------------------------
# kernels uint16 take extra parameter for defining the bitdepth
# -----------------------------------------------------------------
cdef inline dtype_t kernel_autolevel(Py_ssize_t* histo, float pop,
dtype_t g, Py_ssize_t bitdepth,
Py_ssize_t maxbin, Py_ssize_t midbin,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t i, imin, imax, delta
if pop:
for i in range(maxbin - 1, -1, -1):
if histo[i]:
imax = i
break
for i in range(maxbin):
if histo[i]:
imin = i
break
delta = imax - imin
if delta > 0:
return <dtype_t>(1. * (maxbin - 1) * (g - imin) / delta)
else:
return <dtype_t>(imax - imin)
cdef inline dtype_t kernel_bottomhat(Py_ssize_t* histo, float pop,
dtype_t g, Py_ssize_t bitdepth,
Py_ssize_t maxbin, Py_ssize_t midbin,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t i
if pop:
for i in range(maxbin):
if histo[i]:
break
return <dtype_t>(g - i)
else:
return <dtype_t>(0)
cdef inline dtype_t kernel_equalize(Py_ssize_t* histo, float pop,
dtype_t g, Py_ssize_t bitdepth,
Py_ssize_t maxbin, Py_ssize_t midbin,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t i
cdef float sum = 0.
if pop:
for i in range(maxbin):
sum += histo[i]
if i >= g:
break
return <dtype_t>(((maxbin - 1) * sum) / pop)
else:
return <dtype_t>(0)
cdef inline dtype_t kernel_gradient(Py_ssize_t* histo, float pop,
dtype_t g, Py_ssize_t bitdepth,
Py_ssize_t maxbin, Py_ssize_t midbin,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t i, imin, imax
if pop:
for i in range(maxbin - 1, -1, -1):
if histo[i]:
imax = i
break
for i in range(maxbin):
if histo[i]:
imin = i
break
return <dtype_t>(imax - imin)
else:
return <dtype_t>(0)
cdef inline dtype_t kernel_maximum(Py_ssize_t* histo, float pop,
dtype_t g, Py_ssize_t bitdepth,
Py_ssize_t maxbin, Py_ssize_t midbin,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t i
if pop:
for i in range(maxbin - 1, -1, -1):
if histo[i]:
return <dtype_t>(i)
return <dtype_t>(0)
cdef inline dtype_t kernel_mean(Py_ssize_t* histo, float pop,
dtype_t g, Py_ssize_t bitdepth,
Py_ssize_t maxbin, Py_ssize_t midbin,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t i
cdef float mean = 0.
if pop:
for i in range(maxbin):
mean += histo[i] * i
return <dtype_t>(mean / pop)
else:
return <dtype_t>(0)
cdef inline dtype_t kernel_meansubtraction(Py_ssize_t* histo,
float pop,
dtype_t g,
Py_ssize_t bitdepth,
Py_ssize_t maxbin,
Py_ssize_t midbin,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t i
cdef float mean = 0.
if pop:
for i in range(maxbin):
mean += histo[i] * i
return <dtype_t>((g - mean / pop) / 2. + (midbin - 1))
else:
return <dtype_t>(0)
cdef inline dtype_t kernel_median(Py_ssize_t* histo, float pop,
dtype_t g, Py_ssize_t bitdepth,
Py_ssize_t maxbin, Py_ssize_t midbin,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t i
cdef float sum = pop / 2.0
if pop:
for i in range(maxbin):
if histo[i]:
sum -= histo[i]
if sum < 0:
return <dtype_t>(i)
else:
return <dtype_t>(0)
cdef inline dtype_t kernel_minimum(Py_ssize_t* histo, float pop,
dtype_t g, Py_ssize_t bitdepth,
Py_ssize_t maxbin, Py_ssize_t midbin,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t i
if pop:
for i in range(maxbin):
if histo[i]:
return <dtype_t>(i)
else:
return <dtype_t>(0)
cdef inline dtype_t kernel_modal(Py_ssize_t* histo, float pop,
dtype_t g, Py_ssize_t bitdepth,
Py_ssize_t maxbin, Py_ssize_t midbin,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t hmax = 0, imax = 0
if pop:
for i in range(maxbin):
if histo[i] > hmax:
hmax = histo[i]
imax = i
return <dtype_t>(imax)
else:
return <dtype_t>(0)
cdef inline dtype_t kernel_morph_contr_enh(Py_ssize_t* histo,
float pop,
dtype_t g,
Py_ssize_t bitdepth,
Py_ssize_t maxbin,
Py_ssize_t midbin,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t i, imin, imax
if pop:
for i in range(maxbin - 1, -1, -1):
if histo[i]:
imax = i
break
for i in range(maxbin):
if histo[i]:
imin = i
break
if imax - g < g - imin:
return <dtype_t>(imax)
else:
return <dtype_t>(imin)
else:
return <dtype_t>(0)
cdef inline dtype_t kernel_pop(Py_ssize_t* histo, float pop,
dtype_t g, Py_ssize_t bitdepth,
Py_ssize_t maxbin, Py_ssize_t midbin,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
return <dtype_t>(pop)
cdef inline dtype_t kernel_threshold(Py_ssize_t* histo, float pop,
dtype_t g, Py_ssize_t bitdepth,
Py_ssize_t maxbin, Py_ssize_t midbin,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t i
cdef float mean = 0.
if pop:
for i in range(maxbin):
mean += histo[i] * i
return <dtype_t>(g > (mean / pop))
else:
return <dtype_t>(0)
cdef inline dtype_t kernel_tophat(Py_ssize_t* histo, float pop,
dtype_t g, Py_ssize_t bitdepth,
Py_ssize_t maxbin, Py_ssize_t midbin,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t i
if pop:
for i in range(maxbin - 1, -1, -1):
if histo[i]:
break
return <dtype_t>(i - g)
else:
return <dtype_t>(0)
cdef inline dtype_t kernel_entropy(Py_ssize_t* histo, float pop,
dtype_t g, Py_ssize_t bitdepth,
Py_ssize_t maxbin, Py_ssize_t midbin,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t i
cdef float e, p
if pop:
e = 0.
for i in range(maxbin):
p = histo[i] / pop
if p > 0:
e -= p * log(p) / 0.6931471805599453
return <dtype_t>e * 1000
else:
return <dtype_t>(0)
# -----------------------------------------------------------------
# python wrappers
# -----------------------------------------------------------------
def autolevel(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
_core16(kernel_autolevel, image, selem, mask, out, shift_x, shift_y,
bitdepth, 0, 0, <Py_ssize_t>0, <Py_ssize_t>0)
def bottomhat(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
_core16(kernel_bottomhat, image, selem, mask, out, shift_x, shift_y,
bitdepth, 0, 0, <Py_ssize_t>0, <Py_ssize_t>0)
def equalize(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
_core16(kernel_equalize, image, selem, mask, out, shift_x, shift_y,
bitdepth, 0, 0, <Py_ssize_t>0, <Py_ssize_t>0)
def gradient(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
_core16(kernel_gradient, image, selem, mask, out, shift_x, shift_y,
bitdepth, 0, 0, <Py_ssize_t>0, <Py_ssize_t>0)
def maximum(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
_core16(kernel_maximum, image, selem, mask, out, shift_x, shift_y,
bitdepth, 0, 0, <Py_ssize_t>0, <Py_ssize_t>0)
def mean(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
_core16(kernel_mean, image, selem, mask, out, shift_x, shift_y,
bitdepth, 0, 0, <Py_ssize_t>0, <Py_ssize_t>0)
def meansubtraction(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
_core16(kernel_meansubtraction, image, selem, mask, out, shift_x, shift_y,
bitdepth, 0, 0, <Py_ssize_t>0, <Py_ssize_t>0)
def median(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
_core16(kernel_median, image, selem, mask, out, shift_x, shift_y,
bitdepth, 0, 0, <Py_ssize_t>0, <Py_ssize_t>0)
def minimum(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
_core16(kernel_minimum, image, selem, mask, out, shift_x, shift_y,
bitdepth, 0, 0, <Py_ssize_t>0, <Py_ssize_t>0)
def morph_contr_enh(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
_core16(kernel_morph_contr_enh, image, selem, mask, out, shift_x, shift_y,
bitdepth, 0, 0, <Py_ssize_t>0, <Py_ssize_t>0)
def modal(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
_core16(kernel_modal, image, selem, mask, out, shift_x, shift_y,
bitdepth, 0, 0, <Py_ssize_t>0, <Py_ssize_t>0)
def pop(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
_core16(kernel_pop, image, selem, mask, out, shift_x, shift_y,
bitdepth, 0, 0, <Py_ssize_t>0, <Py_ssize_t>0)
def threshold(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
_core16(kernel_threshold, image, selem, mask, out, shift_x, shift_y,
bitdepth, 0, 0, <Py_ssize_t>0, <Py_ssize_t>0)
def tophat(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
_core16(kernel_tophat, image, selem, mask, out, shift_x, shift_y,
bitdepth, 0, 0, <Py_ssize_t>0, <Py_ssize_t>0)
def entropy(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
_core16(kernel_entropy, image, selem, mask, out, shift_x, shift_y,
bitdepth, 0, 0, <Py_ssize_t>0, <Py_ssize_t>0)
-480
View File
@@ -1,480 +0,0 @@
#cython: cdivision=True
#cython: boundscheck=False
#cython: nonecheck=False
#cython: wraparound=False
cimport numpy as cnp
from libc.math cimport log
from .core8_cy cimport dtype_t, _core8
# -----------------------------------------------------------------
# kernels uint8
# -----------------------------------------------------------------
cdef inline dtype_t kernel_autolevel(Py_ssize_t* histo, float pop,
dtype_t g, float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t i, imin, imax, delta
if pop:
for i in range(255, -1, -1):
if histo[i]:
imax = i
break
for i in range(256):
if histo[i]:
imin = i
break
delta = imax - imin
if delta > 0:
return <dtype_t>(255. * (g - imin) / delta)
else:
return <dtype_t>(imax - imin)
else:
return <dtype_t>(0)
cdef inline dtype_t kernel_bottomhat(Py_ssize_t* histo, float pop,
dtype_t g, float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t i
if pop:
for i in range(256):
if histo[i]:
break
return <dtype_t>(g - i)
else:
return <dtype_t>(0)
cdef inline dtype_t kernel_equalize(Py_ssize_t* histo, float pop,
dtype_t g, float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t i
cdef float sum = 0.
if pop:
for i in range(256):
sum += histo[i]
if i >= g:
break
return <dtype_t>((255 * sum) / pop)
else:
return <dtype_t>(0)
cdef inline dtype_t kernel_gradient(Py_ssize_t* histo, float pop,
dtype_t g, float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t i, imin, imax
if pop:
for i in range(255, -1, -1):
if histo[i]:
imax = i
break
for i in range(256):
if histo[i]:
imin = i
break
return <dtype_t>(imax - imin)
else:
return <dtype_t>(0)
cdef inline dtype_t kernel_maximum(Py_ssize_t* histo, float pop,
dtype_t g, float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t i
if pop:
for i in range(255, -1, -1):
if histo[i]:
return <dtype_t>(i)
else:
return <dtype_t>(0)
cdef inline dtype_t kernel_mean(Py_ssize_t* histo, float pop,
dtype_t g, float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t i
cdef float mean = 0.
if pop:
for i in range(256):
mean += histo[i] * i
return <dtype_t>(mean / pop)
else:
return <dtype_t>(0)
cdef inline dtype_t kernel_meansubtraction(Py_ssize_t* histo, float pop,
dtype_t g, float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t i
cdef float mean = 0.
if pop:
for i in range(256):
mean += histo[i] * i
return <dtype_t>((g - mean / pop) / 2. + 127)
else:
return <dtype_t>(0)
cdef inline dtype_t kernel_median(Py_ssize_t* histo, float pop,
dtype_t g, float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t i
cdef float sum = pop / 2.0
if pop:
for i in range(256):
if histo[i]:
sum -= histo[i]
if sum < 0:
return <dtype_t>(i)
else:
return <dtype_t>(0)
cdef inline dtype_t kernel_minimum(Py_ssize_t* histo, float pop,
dtype_t g, float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t i
if pop:
for i in range(256):
if histo[i]:
return <dtype_t>(i)
else:
return <dtype_t>(0)
cdef inline dtype_t kernel_modal(Py_ssize_t* histo, float pop,
dtype_t g, float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t hmax = 0, imax = 0
if pop:
for i in range(256):
if histo[i] > hmax:
hmax = histo[i]
imax = i
return <dtype_t>(imax)
else:
return <dtype_t>(0)
cdef inline dtype_t kernel_morph_contr_enh(Py_ssize_t* histo, float pop,
dtype_t g, float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t i, imin, imax
if pop:
for i in range(255, -1, -1):
if histo[i]:
imax = i
break
for i in range(256):
if histo[i]:
imin = i
break
if imax - g < g - imin:
return <dtype_t>(imax)
else:
return <dtype_t>(imin)
else:
return <dtype_t>(0)
cdef inline dtype_t kernel_pop(Py_ssize_t* histo, float pop,
dtype_t g, float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
return <dtype_t>(pop)
cdef inline dtype_t kernel_threshold(Py_ssize_t* histo, float pop,
dtype_t g, float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t i
cdef float mean = 0.
if pop:
for i in range(256):
mean += histo[i] * i
return <dtype_t>(g > (mean / pop))
else:
return <dtype_t>(0)
cdef inline dtype_t kernel_tophat(Py_ssize_t* histo, float pop,
dtype_t g, float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t i
if pop:
for i in range(255, -1, -1):
if histo[i]:
break
return <dtype_t>(i - g)
else:
return <dtype_t>(0)
cdef inline dtype_t kernel_noise_filter(Py_ssize_t* histo, float pop,
dtype_t g, float p0, float 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 <dtype_t>0
for i in range(g, -1, -1):
if histo[i]:
break
min_i = g - i
for i in range(g, 256):
if histo[i]:
break
if i - g < min_i:
return <dtype_t>(i - g)
else:
return <dtype_t>min_i
cdef inline dtype_t kernel_entropy(Py_ssize_t* histo, float pop,
dtype_t g, float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t i
cdef float e, p
if pop:
e = 0.
for i in range(256):
p = histo[i] / pop
if p > 0:
e -= p * log(p) / 0.6931471805599453
return <dtype_t>e * 10
else:
return <dtype_t>(0)
cdef inline dtype_t kernel_otsu(Py_ssize_t* histo, float pop, dtype_t g,
float p0, float p1, Py_ssize_t s0,
Py_ssize_t s1):
cdef Py_ssize_t i
cdef Py_ssize_t max_i
cdef float P, mu1, mu2, q1, new_q1, sigma_b, max_sigma_b
cdef float mu = 0.
# compute local mean
if pop:
for i in range(256):
mu += histo[i] * i
mu = (mu / pop)
else:
return <dtype_t>(0)
# maximizing the between class variance
max_i = 0
q1 = histo[0] / pop
m1 = 0.
max_sigma_b = 0.
for i in range(1, 256):
P = histo[i] / pop
new_q1 = q1 + P
if new_q1 > 0:
mu1 = (q1 * mu1 + i * P) / new_q1
mu2 = (mu - new_q1 * mu1) / (1. - new_q1)
sigma_b = new_q1 * (1. - new_q1) * (mu1 - mu2) ** 2
if sigma_b > max_sigma_b:
max_sigma_b = sigma_b
max_i = i
q1 = new_q1
return <dtype_t>max_i
# -----------------------------------------------------------------
# python wrappers
# used only internally
# -----------------------------------------------------------------
def autolevel(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0):
_core8(kernel_autolevel, image, selem, mask, out, shift_x, shift_y,
0, 0, <Py_ssize_t>0, <Py_ssize_t>0)
def bottomhat(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0):
_core8(kernel_bottomhat, image, selem, mask, out, shift_x, shift_y,
0, 0, <Py_ssize_t>0, <Py_ssize_t>0)
def equalize(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0):
_core8(kernel_equalize, image, selem, mask, out, shift_x, shift_y,
0, 0, <Py_ssize_t>0, <Py_ssize_t>0)
def gradient(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0):
_core8(kernel_gradient, image, selem, mask, out, shift_x, shift_y,
0, 0, <Py_ssize_t>0, <Py_ssize_t>0)
def maximum(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0):
_core8(kernel_maximum, image, selem, mask, out, shift_x, shift_y,
0, 0, <Py_ssize_t>0, <Py_ssize_t>0)
def mean(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0):
_core8(kernel_mean, image, selem, mask, out, shift_x, shift_y,
0, 0, <Py_ssize_t>0, <Py_ssize_t>0)
def meansubtraction(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0):
_core8(kernel_meansubtraction, image, selem, mask, out, shift_x, shift_y,
0, 0, <Py_ssize_t>0, <Py_ssize_t>0)
def median(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0):
_core8(kernel_median, image, selem, mask, out, shift_x, shift_y,
0, 0, <Py_ssize_t>0, <Py_ssize_t>0)
def minimum(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0):
_core8(kernel_minimum, image, selem, mask, out, shift_x, shift_y,
0, 0, <Py_ssize_t>0, <Py_ssize_t>0)
def morph_contr_enh(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0):
_core8(kernel_morph_contr_enh, image, selem, mask, out, shift_x, shift_y,
0, 0, <Py_ssize_t>0, <Py_ssize_t>0)
def modal(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0):
_core8(kernel_modal, image, selem, mask, out, shift_x, shift_y,
0, 0, <Py_ssize_t>0, <Py_ssize_t>0)
def pop(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0):
_core8(kernel_pop, image, selem, mask, out, shift_x, shift_y,
0, 0, <Py_ssize_t>0, <Py_ssize_t>0)
def threshold(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0):
_core8(kernel_threshold, image, selem, mask, out, shift_x, shift_y, 0, 0,
<Py_ssize_t>0, <Py_ssize_t>0)
def tophat(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0):
_core8(kernel_tophat, image, selem, mask, out, shift_x, shift_y,
0, 0, <Py_ssize_t>0, <Py_ssize_t>0)
def noise_filter(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0):
_core8(kernel_noise_filter, image, selem, mask, out, shift_x, shift_y,
0, 0, <Py_ssize_t>0, <Py_ssize_t>0)
def entropy(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0):
_core8(kernel_entropy, image, selem, mask, out, shift_x, shift_y,
0, 0, <Py_ssize_t>0, <Py_ssize_t>0)
def otsu(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0):
_core8(kernel_otsu, image, selem, mask, out, shift_x, shift_y,
0, 0, <Py_ssize_t>0, <Py_ssize_t>0)
+576
View File
@@ -0,0 +1,576 @@
#cython: cdivision=True
#cython: boundscheck=False
#cython: nonecheck=False
#cython: wraparound=False
cimport numpy as cnp
from libc.math cimport log
from .core_cy cimport uint8_t, uint16_t, dtype_t, _core
cdef inline dtype_t _kernel_autolevel(Py_ssize_t* histo, float pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t i, imin, imax, delta
if pop:
for i in range(max_bin - 1, -1, -1):
if histo[i]:
imax = i
break
for i in range(max_bin):
if histo[i]:
imin = i
break
delta = imax - imin
if delta > 0:
return <dtype_t>(<float>(max_bin - 1) * (g - imin) / delta)
else:
return <dtype_t>(imax - imin)
else:
return <dtype_t>(0)
cdef inline dtype_t _kernel_bottomhat(Py_ssize_t* histo, float pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
float p0, float 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]:
break
return <dtype_t>(g - i)
else:
return <dtype_t>(0)
cdef inline dtype_t _kernel_equalize(Py_ssize_t* histo, float pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t i
cdef Py_ssize_t sum = 0
if pop:
for i in range(max_bin):
sum += histo[i]
if i >= g:
break
return <dtype_t>(((max_bin - 1) * sum) / pop)
else:
return <dtype_t>(0)
cdef inline dtype_t _kernel_gradient(Py_ssize_t* histo, float pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t i, imin, imax
if pop:
for i in range(max_bin - 1, -1, -1):
if histo[i]:
imax = i
break
for i in range(max_bin):
if histo[i]:
imin = i
break
return <dtype_t>(imax - imin)
else:
return <dtype_t>(0)
cdef inline dtype_t _kernel_maximum(Py_ssize_t* histo, float pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
float p0, float 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 <dtype_t>(i)
else:
return <dtype_t>(0)
cdef inline dtype_t _kernel_mean(Py_ssize_t* histo, float pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
float p0, float 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 <dtype_t>(mean / pop)
else:
return <dtype_t>(0)
cdef inline dtype_t _kernel_meansubtraction(Py_ssize_t* histo, float pop,
dtype_t g, Py_ssize_t max_bin,
Py_ssize_t mid_bin, float p0,
float 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 <dtype_t>((g - mean / pop) / 2. + 127)
else:
return <dtype_t>(0)
cdef inline dtype_t _kernel_median(Py_ssize_t* histo, float pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t i
cdef float sum = pop / 2.0
if pop:
for i in range(max_bin):
if histo[i]:
sum -= histo[i]
if sum < 0:
return <dtype_t>(i)
else:
return <dtype_t>(0)
cdef inline dtype_t _kernel_minimum(Py_ssize_t* histo, float pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
float p0, float 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 <dtype_t>(i)
else:
return <dtype_t>(0)
cdef inline dtype_t _kernel_modal(Py_ssize_t* histo, float pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t hmax = 0, imax = 0
if pop:
for i in range(max_bin):
if histo[i] > hmax:
hmax = histo[i]
imax = i
return <dtype_t>(imax)
else:
return <dtype_t>(0)
cdef inline dtype_t _kernel_morph_contr_enh(Py_ssize_t* histo, float pop,
dtype_t g, Py_ssize_t max_bin,
Py_ssize_t mid_bin, float p0,
float p1, Py_ssize_t s0,
Py_ssize_t s1):
cdef Py_ssize_t i, imin, imax
if pop:
for i in range(max_bin - 1, -1, -1):
if histo[i]:
imax = i
break
for i in range(max_bin):
if histo[i]:
imin = i
break
if imax - g < g - imin:
return <dtype_t>(imax)
else:
return <dtype_t>(imin)
else:
return <dtype_t>(0)
cdef inline dtype_t _kernel_pop(Py_ssize_t* histo, float pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
return <dtype_t>(pop)
cdef inline dtype_t _kernel_threshold(Py_ssize_t* histo, float pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
float p0, float 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 <dtype_t>(g > (mean / pop))
else:
return <dtype_t>(0)
cdef inline dtype_t _kernel_tophat(Py_ssize_t* histo, float pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
float p0, float 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]:
break
return <dtype_t>(i - g)
else:
return <dtype_t>(0)
cdef inline dtype_t _kernel_noise_filter(Py_ssize_t* histo, float pop,
dtype_t g, Py_ssize_t max_bin,
Py_ssize_t mid_bin, float p0, float 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 <dtype_t>0
for i in range(g, -1, -1):
if histo[i]:
break
min_i = g - i
for i in range(g, max_bin):
if histo[i]:
break
if i - g < min_i:
return <dtype_t>(i - g)
else:
return <dtype_t>min_i
cdef inline dtype_t _kernel_entropy(Py_ssize_t* histo, float pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t i
cdef float e, p
if pop:
e = 0.
for i in range(max_bin):
p = histo[i] / pop
if p > 0:
e -= p * log(p) / 0.6931471805599453
return <dtype_t>e
else:
return <dtype_t>(0)
cdef inline dtype_t _kernel_otsu(Py_ssize_t* histo, float pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t i
cdef Py_ssize_t max_i
cdef float P, mu1, mu2, q1, new_q1, sigma_b, max_sigma_b
cdef float mu = 0.
# compute local mean
if pop:
for i in range(max_bin):
mu += histo[i] * i
mu = (mu / pop)
else:
return <dtype_t>(0)
# maximizing the between class variance
max_i = 0
q1 = histo[0] / pop
m1 = 0.
max_sigma_b = 0.
for i in range(1, max_bin):
P = histo[i] / pop
new_q1 = q1 + P
if new_q1 > 0:
mu1 = (q1 * mu1 + i * P) / new_q1
mu2 = (mu - new_q1 * mu1) / (1. - new_q1)
sigma_b = new_q1 * (1. - new_q1) * (mu1 - mu2) ** 2
if sigma_b > max_sigma_b:
max_sigma_b = sigma_b
max_i = i
q1 = new_q1
return <dtype_t>max_i
def _autolevel(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask,
dtype_t[:, ::1] out,
char shift_x, char shift_y, Py_ssize_t max_bin):
if dtype_t is uint8_t:
_core[uint8_t](_kernel_autolevel[uint8_t], image, selem, mask, out,
shift_x, shift_y, 0, 0, 0, 0, max_bin)
elif dtype_t is uint16_t:
_core[uint16_t](_kernel_autolevel[uint16_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[:, ::1] out,
char shift_x, char shift_y, Py_ssize_t max_bin):
if dtype_t is uint8_t:
_core[uint8_t](_kernel_bottomhat[uint8_t], image, selem, mask, out,
shift_x, shift_y, 0, 0, 0, 0, max_bin)
elif dtype_t is uint16_t:
_core[uint16_t](_kernel_bottomhat[uint16_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[:, ::1] out,
char shift_x, char shift_y, Py_ssize_t max_bin):
if dtype_t is uint8_t:
_core[uint8_t](_kernel_equalize[uint8_t], image, selem, mask, out,
shift_x, shift_y, 0, 0, 0, 0, max_bin)
elif dtype_t is uint16_t:
_core[uint16_t](_kernel_equalize[uint16_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[:, ::1] out,
char shift_x, char shift_y, Py_ssize_t max_bin):
if dtype_t is uint8_t:
_core[uint8_t](_kernel_gradient[uint8_t], image, selem, mask, out,
shift_x, shift_y, 0, 0, 0, 0, max_bin)
elif dtype_t is uint16_t:
_core[uint16_t](_kernel_gradient[uint16_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[:, ::1] out,
char shift_x, char shift_y, Py_ssize_t max_bin):
if dtype_t is uint8_t:
_core[uint8_t](_kernel_maximum[uint8_t], image, selem, mask, out,
shift_x, shift_y, 0, 0, 0, 0, max_bin)
elif dtype_t is uint16_t:
_core[uint16_t](_kernel_maximum[uint16_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[:, ::1] out,
char shift_x, char shift_y, Py_ssize_t max_bin):
if dtype_t is uint8_t:
_core[uint8_t](_kernel_mean[uint8_t], image, selem, mask, out,
shift_x, shift_y, 0, 0, 0, 0, max_bin)
elif dtype_t is uint16_t:
_core[uint16_t](_kernel_mean[uint16_t], image, selem, mask, out,
shift_x, shift_y, 0, 0, 0, 0, max_bin)
def _meansubtraction(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask,
dtype_t[:, ::1] out,
char shift_x, char shift_y, Py_ssize_t max_bin):
if dtype_t is uint8_t:
_core[uint8_t](_kernel_meansubtraction[uint8_t], image, selem, mask,
out, shift_x, shift_y, 0, 0, 0, 0, max_bin)
elif dtype_t is uint16_t:
_core[uint16_t](_kernel_meansubtraction[uint16_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[:, ::1] out,
char shift_x, char shift_y, Py_ssize_t max_bin):
if dtype_t is uint8_t:
_core[uint8_t](_kernel_median[uint8_t], image, selem, mask, out,
shift_x, shift_y, 0, 0, 0, 0, max_bin)
elif dtype_t is uint16_t:
_core[uint16_t](_kernel_median[uint16_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[:, ::1] out,
char shift_x, char shift_y, Py_ssize_t max_bin):
if dtype_t is uint8_t:
_core[uint8_t](_kernel_minimum[uint8_t], image, selem, mask, out,
shift_x, shift_y, 0, 0, 0, 0, max_bin)
elif dtype_t is uint16_t:
_core[uint16_t](_kernel_minimum[uint16_t], image, selem, mask, out,
shift_x, shift_y, 0, 0, 0, 0, max_bin)
def _morph_contr_enh(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask,
dtype_t[:, ::1] out,
char shift_x, char shift_y, Py_ssize_t max_bin):
if dtype_t is uint8_t:
_core[uint8_t](_kernel_morph_contr_enh[uint8_t], image, selem, mask,
out, shift_x, shift_y, 0, 0, 0, 0, max_bin)
elif dtype_t is uint16_t:
_core[uint16_t](_kernel_morph_contr_enh[uint16_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[:, ::1] out,
char shift_x, char shift_y, Py_ssize_t max_bin):
if dtype_t is uint8_t:
_core[uint8_t](_kernel_modal[uint8_t], image, selem, mask, out,
shift_x, shift_y, 0, 0, 0, 0, max_bin)
elif dtype_t is uint16_t:
_core[uint16_t](_kernel_modal[uint16_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[:, ::1] out,
char shift_x, char shift_y, Py_ssize_t max_bin):
if dtype_t is uint8_t:
_core[uint8_t](_kernel_pop[uint8_t], image, selem, mask, out,
shift_x, shift_y, 0, 0, 0, 0, max_bin)
elif dtype_t is uint16_t:
_core[uint16_t](_kernel_pop[uint16_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[:, ::1] out,
char shift_x, char shift_y, Py_ssize_t max_bin):
if dtype_t is uint8_t:
_core[uint8_t](_kernel_threshold[uint8_t], image, selem, mask, out,
shift_x, shift_y, 0, 0, 0, 0, max_bin)
elif dtype_t is uint16_t:
_core[uint16_t](_kernel_threshold[uint16_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[:, ::1] out,
char shift_x, char shift_y, Py_ssize_t max_bin):
if dtype_t is uint8_t:
_core[uint8_t](_kernel_tophat[uint8_t], image, selem, mask, out,
shift_x, shift_y, 0, 0, 0, 0, max_bin)
elif dtype_t is uint16_t:
_core[uint16_t](_kernel_tophat[uint16_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[:, ::1] out,
char shift_x, char shift_y, Py_ssize_t max_bin):
if dtype_t is uint8_t:
_core[uint8_t](_kernel_noise_filter[uint8_t], image, selem, mask, out,
shift_x, shift_y, 0, 0, 0, 0, max_bin)
elif dtype_t is uint16_t:
_core[uint16_t](_kernel_noise_filter[uint16_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[:, ::1] out,
char shift_x, char shift_y, Py_ssize_t max_bin):
if dtype_t is uint8_t:
_core[uint8_t](_kernel_entropy[uint8_t], image, selem, mask, out,
shift_x, shift_y, 0, 0, 0, 0, max_bin)
elif dtype_t is uint16_t:
_core[uint16_t](_kernel_entropy[uint16_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[:, ::1] out,
char shift_x, char shift_y, Py_ssize_t max_bin):
if dtype_t is uint8_t:
_core[uint8_t](_kernel_otsu[uint8_t], image, selem, mask, out,
shift_x, shift_y, 0, 0, 0, 0, max_bin)
elif dtype_t is uint16_t:
_core[uint16_t](_kernel_otsu[uint16_t], image, selem, mask, out,
shift_x, shift_y, 0, 0, 0, 0, max_bin)
+63 -93
View File
@@ -26,8 +26,8 @@ from skimage import img_as_ubyte
from ... import get_log
log = get_log()
from . import percentile8_cy, percentile16_cy
from .generic import find_bitdepth
from . import percentile_cy
from .generic import _handle_input
__all__ = ['percentile_autolevel', 'percentile_gradient',
@@ -36,45 +36,18 @@ __all__ = ['percentile_autolevel', 'percentile_gradient',
'percentile_threshold']
def _apply(func8, func16, image, selem, out, mask, shift_x, shift_y, p0, p1):
selem = img_as_ubyte(selem > 0)
image = np.ascontiguousarray(image)
def _apply(func, image, selem, out, mask, shift_x, shift_y, p0, p1):
if mask is None:
mask = np.ones(image.shape, dtype=np.uint8)
else:
mask = np.ascontiguousarray(mask)
mask = img_as_ubyte(mask)
image, selem, out, mask, max_bin = _handle_input(image, selem, out, mask)
if image is out:
raise NotImplementedError("Cannot perform rank operation in place.")
if image.dtype == np.uint8:
if func8 is None:
raise TypeError("Not implemented for uint8 image.")
if out is None:
out = np.zeros(image.shape, dtype=np.uint8)
func8(image, selem, shift_x=shift_x, shift_y=shift_y,
mask=mask, out=out, p0=p0, p1=p1)
elif image.dtype == np.uint16:
if func16 is None:
raise TypeError("Not implemented for uint16 image.")
if out is None:
out = np.zeros(image.shape, dtype=np.uint16)
bitdepth = find_bitdepth(image)
if bitdepth > 10:
log.warn("Bitdepth of %d may result in bad rank filter "
"performance." % bitdepth)
func16(image, selem, shift_x=shift_x, shift_y=shift_y, mask=mask,
bitdepth=bitdepth + 1, out=out, p0=p0, p1=p1)
else:
raise TypeError("Only uint8 and uint16 image supported.")
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
def percentile_autolevel(image, selem, out=None, mask=None, shift_x=False,
shift_y=False, p0=.0, p1=1.):
shift_y=False, p0=0, p1=1):
"""Return greyscale local autolevel of an image.
Autolevel is computed on the given structuring element. Only levels between
@@ -82,13 +55,13 @@ def percentile_autolevel(image, selem, out=None, mask=None, shift_x=False,
Parameters
----------
image : ndarray
Image array (uint8 array or uint16).
image : ndarray (uint8, uint16)
Image array.
selem : ndarray
The neighborhood expressed as a 2-D array of 1's and 0's.
out : ndarray
out : ndarray (same dtype as input)
If None, a new array will be allocated.
mask : ndarray (uint8)
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
@@ -101,18 +74,18 @@ def percentile_autolevel(image, selem, out=None, mask=None, shift_x=False,
Returns
-------
local autolevel : uint8 array or uint16
local autolevel : ndarray (same dtype as input)
The result of the local autolevel.
"""
return _apply(percentile8_cy.autolevel, percentile16_cy.autolevel,
return _apply(percentile_cy._autolevel,
image, selem, out=out, mask=mask, shift_x=shift_x,
shift_y=shift_y, p0=p0, p1=p1)
def percentile_gradient(image, selem, out=None, mask=None, shift_x=False,
shift_y=False, p0=.0, p1=1.):
shift_y=False, p0=0, p1=1):
"""Return greyscale local percentile_gradient of an image.
percentile_gradient is computed on the given structuring element. Only
@@ -120,13 +93,13 @@ def percentile_gradient(image, selem, out=None, mask=None, shift_x=False,
Parameters
----------
image : ndarray
Image array (uint8 array or uint16).
image : ndarray (uint8, uint16)
Image array.
selem : ndarray
The neighborhood expressed as a 2-D array of 1's and 0's.
out : ndarray
out : ndarray (same dtype as input)
If None, a new array will be allocated.
mask : ndarray (uint8)
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
@@ -139,18 +112,18 @@ def percentile_gradient(image, selem, out=None, mask=None, shift_x=False,
Returns
-------
local percentile_gradient : uint8 array or uint16
local percentile_gradient : ndarray (same dtype as input)
The result of the local percentile_gradient.
"""
return _apply(percentile8_cy.gradient, percentile16_cy.gradient,
return _apply(percentile_cy._gradient,
image, selem, out=out, mask=mask, shift_x=shift_x,
shift_y=shift_y, p0=p0, p1=p1)
def percentile_mean(image, selem, out=None, mask=None, shift_x=False,
shift_y=False, p0=.0, p1=1.):
shift_y=False, p0=0, p1=1):
"""Return greyscale local mean of an image.
Mean is computed on the given structuring element. Only levels between
@@ -158,13 +131,13 @@ def percentile_mean(image, selem, out=None, mask=None, shift_x=False,
Parameters
----------
image : ndarray
Image array (uint8 array or uint16).
image : ndarray (uint8, uint16)
Image array.
selem : ndarray
The neighborhood expressed as a 2-D array of 1's and 0's.
out : ndarray
out : ndarray (same dtype as input)
If None, a new array will be allocated.
mask : ndarray (uint8)
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
@@ -177,18 +150,18 @@ def percentile_mean(image, selem, out=None, mask=None, shift_x=False,
Returns
-------
local mean : uint8 array or uint16
local mean : ndarray (same dtype as input)
The result of the local mean.
"""
return _apply(percentile8_cy.mean, percentile16_cy.mean,
return _apply(percentile_cy._mean,
image, selem, out=out, mask=mask, shift_x=shift_x,
shift_y=shift_y, p0=p0, p1=p1)
def percentile_mean_subtraction(image, selem, out=None, mask=None,
shift_x=False, shift_y=False, p0=.0, p1=1.):
shift_x=False, shift_y=False, p0=0, p1=1):
"""Return greyscale local mean_subtraction of an image.
mean_subtraction is computed on the given structuring element. Only levels
@@ -196,13 +169,13 @@ def percentile_mean_subtraction(image, selem, out=None, mask=None,
Parameters
----------
image : ndarray
Image array (uint8 array or uint16).
image : ndarray (uint8, uint16)
Image array.
selem : ndarray
The neighborhood expressed as a 2-D array of 1's and 0's.
out : ndarray
out : ndarray (same dtype as input)
If None, a new array will be allocated.
mask : ndarray (uint8)
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
@@ -215,19 +188,18 @@ def percentile_mean_subtraction(image, selem, out=None, mask=None,
Returns
-------
local mean_subtraction : uint8 array or uint16
local mean_subtraction : ndarray (same dtype as input)
The result of the local mean_subtraction.
"""
return _apply(percentile8_cy.mean_subtraction,
percentile16_cy.mean_subtraction,
return _apply(percentile_cy._mean_subtraction,
image, selem, out=out, mask=mask, shift_x=shift_x,
shift_y=shift_y, p0=p0, p1=p1)
def percentile_morph_contr_enh(image, selem, out=None, mask=None,
shift_x=False, shift_y=False, p0=.0, p1=1.):
shift_x=False, shift_y=False, p0=0, p1=1):
"""Return greyscale local morph_contr_enh of an image.
morph_contr_enh is computed on the given structuring element. Only levels
@@ -235,13 +207,13 @@ def percentile_morph_contr_enh(image, selem, out=None, mask=None,
Parameters
----------
image : ndarray
Image array (uint8 array or uint16).
image : ndarray (uint8, uint16)
Image array.
selem : ndarray
The neighborhood expressed as a 2-D array of 1's and 0's.
out : ndarray
out : ndarray (same dtype as input)
If None, a new array will be allocated.
mask : ndarray (uint8)
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
@@ -254,19 +226,18 @@ def percentile_morph_contr_enh(image, selem, out=None, mask=None,
Returns
-------
local morph_contr_enh : uint8 array or uint16
local morph_contr_enh : ndarray (same dtype as input)
The result of the local morph_contr_enh.
"""
return _apply(percentile8_cy.morph_contr_enh,
percentile16_cy.morph_contr_enh,
return _apply(percentile_cy._morph_contr_enh,
image, selem, out=out, mask=mask, shift_x=shift_x,
shift_y=shift_y, p0=p0, p1=p1)
def percentile(image, selem, out=None, mask=None, shift_x=False, shift_y=False,
p0=.0):
p0=0):
"""Return greyscale local percentile of an image.
percentile is computed on the given structuring element. Returns the value
@@ -274,13 +245,13 @@ def percentile(image, selem, out=None, mask=None, shift_x=False, shift_y=False,
Parameters
----------
image : ndarray
Image array (uint8 array or uint16).
image : ndarray (uint8, uint16)
Image array.
selem : ndarray
The neighborhood expressed as a 2-D array of 1's and 0's.
out : ndarray
out : ndarray (same dtype as input)
If None, a new array will be allocated.
mask : ndarray (uint8)
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
@@ -292,19 +263,18 @@ def percentile(image, selem, out=None, mask=None, shift_x=False, shift_y=False,
Returns
-------
local percentile : uint8 array or uint16
local percentile : ndarray (same dtype as input)
The result of the local percentile.
"""
return _apply(percentile8_cy.percentile,
percentile16_cy.percentile,
return _apply(percentile_cy._percentile,
image, selem, out=out, mask=mask, shift_x=shift_x,
shift_y=shift_y, p0=p0, p1=0.)
def percentile_pop(image, selem, out=None, mask=None, shift_x=False,
shift_y=False, p0=.0, p1=1.):
shift_y=False, p0=0, p1=1):
"""Return greyscale local pop of an image.
pop is computed on the given structuring element. Only levels between
@@ -312,13 +282,13 @@ def percentile_pop(image, selem, out=None, mask=None, shift_x=False,
Parameters
----------
image : ndarray
Image array (uint8 array or uint16).
image : ndarray (uint8, uint16)
Image array.
selem : ndarray
The neighborhood expressed as a 2-D array of 1's and 0's.
out : ndarray
out : ndarray (same dtype as input)
If None, a new array will be allocated.
mask : ndarray (uint8)
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
@@ -331,18 +301,18 @@ def percentile_pop(image, selem, out=None, mask=None, shift_x=False,
Returns
-------
local pop : uint8 array or uint16
local pop : ndarray (same dtype as input)
The result of the local pop.
"""
return _apply(percentile8_cy.pop, percentile16_cy.pop,
return _apply(percentile_cy._pop,
image, selem, out=out, mask=mask, shift_x=shift_x,
shift_y=shift_y, p0=p0, p1=p1)
def percentile_threshold(image, selem, out=None, mask=None, shift_x=False,
shift_y=False, p0=.0):
shift_y=False, p0=0):
"""Return greyscale local threshold of an image.
threshold is computed on the given structuring element. Returns
@@ -352,13 +322,13 @@ def percentile_threshold(image, selem, out=None, mask=None, shift_x=False,
Parameters
----------
image : ndarray
Image array (uint8 array or uint16).
image : ndarray (uint8, uint16)
Image array.
selem : ndarray
The neighborhood expressed as a 2-D array of 1's and 0's.
out : ndarray
out : ndarray (same dtype as input)
If None, a new array will be allocated.
mask : ndarray (uint8)
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
@@ -370,11 +340,11 @@ def percentile_threshold(image, selem, out=None, mask=None, shift_x=False,
Returns
-------
local threshold : uint8 array or uint16
local threshold : ndarray (same dtype as input)
The result of the local threshold.
"""
return _apply(percentile8_cy.threshold, percentile16_cy.threshold,
return _apply(percentile_cy._threshold,
image, selem, out=out, mask=mask, shift_x=shift_x,
shift_y=shift_y, p0=p0, p1=0.)
shift_y=shift_y, p0=p0, p1=0)
-326
View File
@@ -1,326 +0,0 @@
#cython: cdivision=True
#cython: boundscheck=False
#cython: nonecheck=False
#cython: wraparound=False
cimport numpy as cnp
from .core16_cy cimport dtype_t, _core16, uint16_min, uint16_max
# -----------------------------------------------------------------
# kernels uint16 (SOFT version using percentiles)
# -----------------------------------------------------------------
cdef inline dtype_t kernel_autolevel(Py_ssize_t* histo, float pop,
dtype_t g, Py_ssize_t bitdepth,
Py_ssize_t maxbin, Py_ssize_t midbin,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef int i, imin, imax, sum, delta
if pop:
sum = 0
p1 = 1.0 - p1
for i in range(maxbin):
sum += histo[i]
if sum > p0 * pop:
imin = i
break
sum = 0
for i in range(maxbin - 1, -1, -1):
sum += histo[i]
if sum > p1 * pop:
imax = i
break
delta = imax - imin
if delta > 0:
return <dtype_t>(1.0 * (maxbin - 1)
* (uint16_min(uint16_max(imin, g), imax)
- imin) / delta)
else:
return <dtype_t>(imax - imin)
else:
return <dtype_t>(0)
cdef inline dtype_t kernel_gradient(Py_ssize_t* histo, float pop,
dtype_t g, Py_ssize_t bitdepth,
Py_ssize_t maxbin, Py_ssize_t midbin,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef int i, imin, imax, sum, delta
if pop:
sum = 0
p1 = 1.0 - p1
for i in range(maxbin):
sum += histo[i]
if sum >= p0 * pop:
imin = i
break
sum = 0
for i in range((maxbin - 1), -1, -1):
sum += histo[i]
if sum >= p1 * pop:
imax = i
break
return <dtype_t>(imax - imin)
else:
return <dtype_t>(0)
cdef inline dtype_t kernel_mean(Py_ssize_t* histo, float pop,
dtype_t g, Py_ssize_t bitdepth,
Py_ssize_t maxbin, Py_ssize_t midbin,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef int i, sum, mean, n
if pop:
sum = 0
mean = 0
n = 0
for i in range(maxbin):
sum += histo[i]
if (sum >= p0 * pop) and (sum <= p1 * pop):
n += histo[i]
mean += histo[i] * i
if n > 0:
return <dtype_t>(1.0 * mean / n)
else:
return <dtype_t>(0)
else:
return <dtype_t>(0)
cdef inline dtype_t kernel_mean_subtraction(Py_ssize_t* histo,
float pop,
dtype_t g,
Py_ssize_t bitdepth,
Py_ssize_t maxbin,
Py_ssize_t midbin,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef int i, sum, mean, n
if pop:
sum = 0
mean = 0
n = 0
for i in range(maxbin):
sum += histo[i]
if (sum >= p0 * pop) and (sum <= p1 * pop):
n += histo[i]
mean += histo[i] * i
if n > 0:
return <dtype_t>((g - (mean / n)) * .5 + midbin)
else:
return <dtype_t>(0)
else:
return <dtype_t>(0)
cdef inline dtype_t kernel_morph_contr_enh(Py_ssize_t* histo,
float pop,
dtype_t g,
Py_ssize_t bitdepth,
Py_ssize_t maxbin,
Py_ssize_t midbin,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef int i, imin, imax, sum, delta
if pop:
sum = 0
p1 = 1.0 - p1
for i in range(maxbin):
sum += histo[i]
if sum > p0 * pop:
imin = i
break
sum = 0
for i in range((maxbin - 1), -1, -1):
sum += histo[i]
if sum > p1 * pop:
imax = i
break
if g > imax:
return <dtype_t>imax
if g < imin:
return <dtype_t>imin
if imax - g < g - imin:
return <dtype_t>imax
else:
return <dtype_t>imin
else:
return <dtype_t>(0)
cdef inline dtype_t kernel_percentile(Py_ssize_t* histo, float pop,
dtype_t g, Py_ssize_t bitdepth,
Py_ssize_t maxbin, Py_ssize_t midbin,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef int i
cdef float sum = 0.
if pop:
for i in range(maxbin):
sum += histo[i]
if sum >= p0 * pop:
break
return <dtype_t>(i)
else:
return <dtype_t>(0)
cdef inline dtype_t kernel_pop(Py_ssize_t* histo, float pop,
dtype_t g, Py_ssize_t bitdepth,
Py_ssize_t maxbin, Py_ssize_t midbin,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef int i, sum, n
if pop:
sum = 0
n = 0
for i in range(maxbin):
sum += histo[i]
if (sum >= p0 * pop) and (sum <= p1 * pop):
n += histo[i]
return <dtype_t>(n)
else:
return <dtype_t>(0)
cdef inline dtype_t kernel_threshold(Py_ssize_t* histo, float pop,
dtype_t g, Py_ssize_t bitdepth,
Py_ssize_t maxbin, Py_ssize_t midbin,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef int i
cdef float sum = 0.
if pop:
for i in range(maxbin):
sum += histo[i]
if sum >= p0 * pop:
break
return <dtype_t>((maxbin - 1) * (g >= i))
else:
return <dtype_t>(0)
# -----------------------------------------------------------------
# python wrappers
# -----------------------------------------------------------------
def autolevel(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0, int bitdepth=8,
float p0=0., float p1=0.):
"""bottom hat
"""
_core16(kernel_autolevel, image, selem, mask, out, shift_x, shift_y,
bitdepth, p0, p1, <Py_ssize_t>0, <Py_ssize_t>0)
def gradient(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0, int bitdepth=8,
float p0=0., float p1=0.):
"""return p0,p1 percentile gradient
"""
_core16(kernel_gradient, image, selem, mask, out, shift_x, shift_y,
bitdepth, p0, p1, <Py_ssize_t>0, <Py_ssize_t>0)
def mean(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0, int bitdepth=8,
float p0=0., float p1=0.):
"""return mean between [p0 and p1] percentiles
"""
_core16(kernel_mean, image, selem, mask, out, shift_x, shift_y,
bitdepth, p0, p1, <Py_ssize_t>0, <Py_ssize_t>0)
def mean_subtraction(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0, int bitdepth=8,
float p0=0., float p1=0.):
"""return original - mean between [p0 and p1] percentiles *.5 +127
"""
_core16(
kernel_mean_subtraction, image, selem, mask, out, shift_x, shift_y,
bitdepth, p0, p1, <Py_ssize_t>0, <Py_ssize_t>0)
def morph_contr_enh(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0, int bitdepth=8,
float p0=0., float p1=0.):
"""reforce contrast using percentiles
"""
_core16(kernel_morph_contr_enh, image, selem, mask, out, shift_x, shift_y,
bitdepth, p0, p1, <Py_ssize_t>0, <Py_ssize_t>0)
def percentile(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0, int bitdepth=8,
float p0=0.):
"""return p0 percentile
"""
_core16(kernel_percentile, image, selem, mask, out, shift_x, shift_y,
bitdepth, p0, .0, <Py_ssize_t>0, <Py_ssize_t>0)
def pop(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0, int bitdepth=8,
float p0=0., float p1=0.):
"""return nb of pixels between [p0 and p1]
"""
_core16(kernel_pop, image, selem, mask, out, shift_x, shift_y,
bitdepth, p0, p1, <Py_ssize_t>0, <Py_ssize_t>0)
def threshold(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0, int bitdepth=8,
float p0=0.):
"""return (maxbin-1) if g > percentile p0
"""
_core16(kernel_threshold, image, selem, mask, out, shift_x, shift_y,
bitdepth, p0, 0., <Py_ssize_t>0, <Py_ssize_t>0)
-291
View File
@@ -1,291 +0,0 @@
#cython: cdivision=True
#cython: boundscheck=False
#cython: nonecheck=False
#cython: wraparound=False
cimport numpy as cnp
from .core8_cy cimport dtype_t, _core8, uint8_max, uint8_min
# -----------------------------------------------------------------
# kernels uint8 (SOFT version using percentiles)
# -----------------------------------------------------------------
cdef inline dtype_t kernel_autolevel(Py_ssize_t* histo, float pop,
dtype_t g, float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef int i, imin, imax, sum, delta
if pop:
sum = 0
p1 = 1.0 - p1
imin = 0
imax = 255
for i in range(256):
sum += histo[i]
if sum > (p0 * pop):
imin = i
break
sum = 0
for i in range(255, -1, -1):
sum += histo[i]
if sum > (p1 * pop):
imax = i
break
delta = imax - imin
if delta > 0:
return <dtype_t>(255 * (uint8_min(uint8_max(imin, g), imax)
- imin) / delta)
else:
return <dtype_t>(imax - imin)
else:
return <dtype_t>(128)
cdef inline dtype_t kernel_gradient(Py_ssize_t* histo, float pop,
dtype_t g, float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef int i, imin, imax, sum, delta
if pop:
sum = 0
p1 = 1.0 - p1
for i in range(256):
sum += histo[i]
if sum >= p0 * pop:
imin = i
break
sum = 0
for i in range(255, -1, -1):
sum += histo[i]
if sum >= p1 * pop:
imax = i
break
return <dtype_t>(imax - imin)
else:
return <dtype_t>(0)
cdef inline dtype_t kernel_mean(Py_ssize_t* histo, float pop,
dtype_t g, float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef int i, sum, mean, n
if pop:
sum = 0
mean = 0
n = 0
for i in range(256):
sum += histo[i]
if (sum >= p0 * pop) and (sum <= p1 * pop):
n += histo[i]
mean += histo[i] * i
if n > 0:
return <dtype_t>(1.0 * mean / n)
else:
return <dtype_t>(0)
else:
return <dtype_t>(0)
cdef inline dtype_t kernel_mean_subtraction(Py_ssize_t* histo,
float pop,
dtype_t g,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef int i, sum, mean, n
if pop:
sum = 0
mean = 0
n = 0
for i in range(256):
sum += histo[i]
if (sum >= p0 * pop) and (sum <= p1 * pop):
n += histo[i]
mean += histo[i] * i
if n > 0:
return <dtype_t>((g - (mean / n)) * .5 + 127)
else:
return <dtype_t>(0)
else:
return <dtype_t>(0)
cdef inline dtype_t kernel_morph_contr_enh(Py_ssize_t* histo,
float pop,
dtype_t g, float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef int i, imin, imax, sum, delta
if pop:
sum = 0
p1 = 1.0 - p1
for i in range(256):
sum += histo[i]
if sum >= p0 * pop:
imin = i
break
sum = 0
for i in range(255, -1, -1):
sum += histo[i]
if sum >= p1 * pop:
imax = i
break
if g > imax:
return <dtype_t>imax
if g < imin:
return <dtype_t>imin
if imax - g < g - imin:
return <dtype_t>imax
else:
return <dtype_t>imin
else:
return <dtype_t>(0)
cdef inline dtype_t kernel_percentile(Py_ssize_t* histo, float pop,
dtype_t g, float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef int i
cdef float sum = 0.
if pop:
for i in range(256):
sum += histo[i]
if sum >= p0 * pop:
break
return <dtype_t>(i)
else:
return <dtype_t>(0)
cdef inline dtype_t kernel_pop(Py_ssize_t* histo, float pop,
dtype_t g, float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef int i, sum, n
if pop:
sum = 0
n = 0
for i in range(256):
sum += histo[i]
if (sum >= p0 * pop) and (sum <= p1 * pop):
n += histo[i]
return <dtype_t>(n)
else:
return <dtype_t>(0)
cdef inline dtype_t kernel_threshold(Py_ssize_t* histo, float pop,
dtype_t g, float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef int i
cdef float sum = 0.
if pop:
for i in range(256):
sum += histo[i]
if sum >= p0 * pop:
break
return <dtype_t>(255 * (g >= i))
else:
return <dtype_t>(0)
# -----------------------------------------------------------------
# python wrappers
# -----------------------------------------------------------------
def autolevel(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0, float p0=0., float p1=0.):
"""autolevel
"""
_core8(kernel_autolevel, image, selem, mask, out, shift_x, shift_y, p0, p1,
<Py_ssize_t>0, <Py_ssize_t>0)
def gradient(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0, float p0=0., float p1=0.):
"""return p0,p1 percentile gradient
"""
_core8(kernel_gradient, image, selem, mask, out, shift_x, shift_y, p0, p1,
<Py_ssize_t>0, <Py_ssize_t>0)
def mean(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0, float p0=0., float p1=0.):
"""return mean between [p0 and p1] percentiles
"""
_core8(kernel_mean, image, selem, mask, out, shift_x, shift_y, p0, p1,
<Py_ssize_t>0, <Py_ssize_t>0)
def mean_subtraction(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0, float p0=0., float p1=0.):
"""return original - mean between [p0 and p1] percentiles *.5 +127
"""
_core8(kernel_mean_subtraction, image, selem, mask, out, shift_x, shift_y,
p0, p1, <Py_ssize_t>0, <Py_ssize_t>0)
def morph_contr_enh(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0, float p0=0., float p1=0.):
"""reforce contrast using percentiles
"""
_core8(kernel_morph_contr_enh, image, selem, mask, out, shift_x, shift_y,
p0, p1, <Py_ssize_t>0, <Py_ssize_t>0)
def percentile(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0, float p0=0.):
"""return p0 percentile
"""
_core8(kernel_percentile, image, selem, mask, out, shift_x, shift_y,
p0, 0., <Py_ssize_t>0, <Py_ssize_t>0)
def pop(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0, float p0=0., float p1=0.):
"""return nb of pixels between [p0 and p1]
"""
_core8(kernel_pop, image, selem, mask, out, shift_x, shift_y, p0, p1,
<Py_ssize_t>0, <Py_ssize_t>0)
def threshold(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask=None,
dtype_t[:, ::1] out=None,
char shift_x=0, char shift_y=0, float p0=0.):
"""return 255 if g > percentile p0
"""
_core8(kernel_threshold, image, selem, mask, out, shift_x, shift_y, p0, 0.,
<Py_ssize_t>0, <Py_ssize_t>0)
+325
View File
@@ -0,0 +1,325 @@
#cython: cdivision=True
#cython: boundscheck=False
#cython: nonecheck=False
#cython: wraparound=False
cimport numpy as cnp
from .core_cy cimport uint8_t, uint16_t, dtype_t, _core, _min, _max
cdef inline dtype_t _kernel_autolevel(Py_ssize_t* histo, float pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t i, imin, imax, sum, delta
if pop:
sum = 0
p1 = 1.0 - p1
for i in range(max_bin - 1):
sum += histo[i]
if sum > p0 * pop:
imin = i
break
sum = 0
for i in range(max_bin - 1, -1, -1):
sum += histo[i]
if sum > p1 * pop:
imax = i
break
delta = imax - imin
if delta > 0:
return <dtype_t>(<float>(max_bin - 1) * (_min(_max(imin, g), imax)
- imin) / delta)
else:
return <dtype_t>(imax - imin)
else:
return <dtype_t>(0)
cdef inline dtype_t _kernel_gradient(Py_ssize_t* histo, float pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t i, imin, imax, sum, delta
if pop:
sum = 0
p1 = 1.0 - p1
for i in range(max_bin):
sum += histo[i]
if sum >= p0 * pop:
imin = i
break
sum = 0
for i in range((max_bin - 1), -1, -1):
sum += histo[i]
if sum >= p1 * pop:
imax = i
break
return <dtype_t>(imax - imin)
else:
return <dtype_t>(0)
cdef inline dtype_t _kernel_mean(Py_ssize_t* histo, float pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t i, sum, mean, n
if pop:
sum = 0
mean = 0
n = 0
for i in range(max_bin):
sum += histo[i]
if (sum >= p0 * pop) and (sum <= p1 * pop):
n += histo[i]
mean += histo[i] * i
if n > 0:
return <dtype_t>(mean / n)
else:
return <dtype_t>(0)
else:
return <dtype_t>(0)
cdef inline dtype_t _kernel_mean_subtraction(Py_ssize_t* histo, float pop,
dtype_t g, Py_ssize_t max_bin,
Py_ssize_t mid_bin, float p0,
float p1, Py_ssize_t s0,
Py_ssize_t s1):
cdef Py_ssize_t i, sum, mean, n
if pop:
sum = 0
mean = 0
n = 0
for i in range(max_bin):
sum += histo[i]
if (sum >= p0 * pop) and (sum <= p1 * pop):
n += histo[i]
mean += histo[i] * i
if n > 0:
return <dtype_t>((g - (mean / n)) * .5 + mid_bin)
else:
return <dtype_t>(0)
else:
return <dtype_t>(0)
cdef inline dtype_t _kernel_morph_contr_enh(Py_ssize_t* histo, float pop,
dtype_t g, Py_ssize_t max_bin,
Py_ssize_t mid_bin, float p0,
float p1, Py_ssize_t s0,
Py_ssize_t s1):
cdef Py_ssize_t i, imin, imax, sum, delta
if pop:
sum = 0
p1 = 1.0 - p1
for i in range(max_bin):
sum += histo[i]
if sum > p0 * pop:
imin = i
break
sum = 0
for i in range((max_bin - 1), -1, -1):
sum += histo[i]
if sum > p1 * pop:
imax = i
break
if g > imax:
return <dtype_t>imax
if g < imin:
return <dtype_t>imin
if imax - g < g - imin:
return <dtype_t>imax
else:
return <dtype_t>imin
else:
return <dtype_t>(0)
cdef inline dtype_t _kernel_percentile(Py_ssize_t* histo, float pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t i
cdef Py_ssize_t sum = 0
if pop:
for i in range(max_bin):
sum += histo[i]
if sum >= p0 * pop:
break
return <dtype_t>(i)
else:
return <dtype_t>(0)
cdef inline dtype_t _kernel_pop(Py_ssize_t* histo, float pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef Py_ssize_t i, sum, n
if pop:
sum = 0
n = 0
for i in range(max_bin):
sum += histo[i]
if (sum >= p0 * pop) and (sum <= p1 * pop):
n += histo[i]
return <dtype_t>(n)
else:
return <dtype_t>(0)
cdef inline dtype_t _kernel_threshold(Py_ssize_t* histo, float pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
float p0, float p1,
Py_ssize_t s0, Py_ssize_t s1):
cdef int i
cdef Py_ssize_t sum = 0
if pop:
for i in range(max_bin):
sum += histo[i]
if sum >= p0 * pop:
break
return <dtype_t>((max_bin - 1) * (g >= i))
else:
return <dtype_t>(0)
def _autolevel(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask,
dtype_t[:, ::1] out,
char shift_x, char shift_y, float p0, float p1,
Py_ssize_t max_bin):
if dtype_t is uint8_t:
_core[uint8_t](_kernel_autolevel[uint8_t], image, selem, mask, out,
shift_x, shift_y, p0, p1, 0, 0, max_bin)
elif dtype_t is uint16_t:
_core[uint16_t](_kernel_autolevel[uint16_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[:, ::1] out,
char shift_x, char shift_y, float p0, float p1,
Py_ssize_t max_bin):
if dtype_t is uint8_t:
_core[uint8_t](_kernel_gradient[uint8_t], image, selem, mask, out,
shift_x, shift_y, p0, p1, 0, 0, max_bin)
elif dtype_t is uint16_t:
_core[uint16_t](_kernel_gradient[uint16_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[:, ::1] out,
char shift_x, char shift_y, float p0, float p1,
Py_ssize_t max_bin):
if dtype_t is uint8_t:
_core[uint8_t](_kernel_mean[uint8_t], image, selem, mask, out,
shift_x, shift_y, p0, p1, 0, 0, max_bin)
elif dtype_t is uint16_t:
_core[uint16_t](_kernel_mean[uint16_t], image, selem, mask, out,
shift_x, shift_y, p0, p1, 0, 0, max_bin)
def _mean_subtraction(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask,
dtype_t[:, ::1] out,
char shift_x, char shift_y, float p0, float p1,
Py_ssize_t max_bin):
if dtype_t is uint8_t:
_core[uint8_t](_kernel_mean_subtraction[uint8_t], image, selem, mask,
out, shift_x, shift_y, p0, p1, 0, 0, max_bin)
elif dtype_t is uint16_t:
_core[uint16_t](_kernel_mean_subtraction[uint16_t], image, selem, mask,
out, shift_x, shift_y, p0, p1, 0, 0, max_bin)
def _morph_contr_enh(dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask,
dtype_t[:, ::1] out,
char shift_x, char shift_y, float p0, float p1,
Py_ssize_t max_bin):
if dtype_t is uint8_t:
_core[uint8_t](_kernel_morph_contr_enh[uint8_t], image, selem, mask,
out, shift_x, shift_y, p0, p1, 0, 0, max_bin)
elif dtype_t is uint16_t:
_core[uint16_t](_kernel_morph_contr_enh[uint16_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[:, ::1] out,
char shift_x, char shift_y, float p0, Py_ssize_t max_bin):
if dtype_t is uint8_t:
_core[uint8_t](_kernel_percentile[uint8_t], image, selem, mask, out,
shift_x, shift_y, p0, 1, 0, 0, max_bin)
elif dtype_t is uint16_t:
_core[uint16_t](_kernel_percentile[uint16_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[:, ::1] out,
char shift_x, char shift_y, float p0, float p1,
Py_ssize_t max_bin):
if dtype_t is uint8_t:
_core[uint8_t](_kernel_pop[uint8_t], image, selem, mask, out,
shift_x, shift_y, p0, p1, 0, 0, max_bin)
elif dtype_t is uint16_t:
_core[uint16_t](_kernel_pop[uint16_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[:, ::1] out,
char shift_x, char shift_y, float p0, Py_ssize_t max_bin):
if dtype_t is uint8_t:
_core[uint8_t](_kernel_threshold[uint8_t], image, selem, mask, out,
shift_x, shift_y, p0, 1, 0, 0, max_bin)
elif dtype_t is uint16_t:
_core[uint16_t](_kernel_threshold[uint16_t], image, selem, mask, out,
shift_x, shift_y, p0, 1, 0, 0, max_bin)
+9 -11
View File
@@ -176,12 +176,10 @@ def test_compare_autolevels_16bit():
assert_array_equal(loc_autolevel, loc_perc_autolevel)
def test_compare_uint_vs_float():
# filters applied on 8-bit image ore 16-bit image (having only real 8-bit of
# dynamic) should be identical
def test_compare_ubyte_vs_float():
# Create signed int8 image that and convert it to uint8
image_uint = img_as_uint(data.camera()[:50, :50])
image_uint = img_as_ubyte(data.camera()[:50, :50])
image_float = img_as_float(image_uint)
methods = ['autolevel', 'bottomhat', 'equalize', 'gradient', 'threshold',
@@ -372,37 +370,37 @@ def test_entropy():
selem = np.ones((16, 16), dtype=np.uint8)
# 1 bit per pixel
data = np.tile(np.asarray([0, 1]), (100, 100)).astype(np.uint8)
assert(np.max(rank.entropy(data, selem)) == 10)
assert(np.max(rank.entropy(data, selem)) == 1)
# 2 bit per pixel
data = np.tile(np.asarray([[0, 1], [2, 3]]), (10, 10)).astype(np.uint8)
assert(np.max(rank.entropy(data, selem)) == 20)
assert(np.max(rank.entropy(data, selem)) == 2)
# 3 bit per pixel
data = np.tile(
np.asarray([[0, 1, 2, 3], [4, 5, 6, 7]]), (10, 10)).astype(np.uint8)
assert(np.max(rank.entropy(data, selem)) == 30)
assert(np.max(rank.entropy(data, selem)) == 3)
# 4 bit per pixel
data = np.tile(
np.reshape(np.arange(16), (4, 4)), (10, 10)).astype(np.uint8)
assert(np.max(rank.entropy(data, selem)) == 40)
assert(np.max(rank.entropy(data, selem)) == 4)
# 6 bit per pixel
data = np.tile(
np.reshape(np.arange(64), (8, 8)), (10, 10)).astype(np.uint8)
assert(np.max(rank.entropy(data, selem)) == 60)
assert(np.max(rank.entropy(data, selem)) == 6)
# 8-bit per pixel
data = np.tile(
np.reshape(np.arange(256), (16, 16)), (10, 10)).astype(np.uint8)
assert(np.max(rank.entropy(data, selem)) == 80)
assert(np.max(rank.entropy(data, selem)) == 8)
# 12 bit per pixel
selem = np.ones((64, 64), dtype=np.uint8)
data = np.tile(
np.reshape(np.arange(4096), (64, 64)), (2, 2)).astype(np.uint16)
assert(np.max(rank.entropy(data, selem)) == 12000)
assert(np.max(rank.entropy(data, selem)) == 12)
def test_selem_dtypes():
+8 -18
View File
@@ -14,34 +14,24 @@ def configuration(parent_package='', top_path=None):
cython(['_ctmf.pyx'], working_path=base_path)
cython(['_denoise_cy.pyx'], working_path=base_path)
cython(['rank/core8_cy.pyx'], working_path=base_path)
cython(['rank/core16_cy.pyx'], working_path=base_path)
cython(['rank/generic8_cy.pyx'], working_path=base_path)
cython(['rank/percentile8_cy.pyx'], working_path=base_path)
cython(['rank/generic16_cy.pyx'], working_path=base_path)
cython(['rank/percentile16_cy.pyx'], working_path=base_path)
cython(['rank/bilateral16_cy.pyx'], working_path=base_path)
cython(['rank/core_cy.pyx'], working_path=base_path)
cython(['rank/generic_cy.pyx'], working_path=base_path)
cython(['rank/percentile_cy.pyx'], working_path=base_path)
cython(['rank/bilateral_cy.pyx'], working_path=base_path)
config.add_extension('_ctmf', sources=['_ctmf.c'],
include_dirs=[get_numpy_include_dirs()])
config.add_extension('_denoise_cy', sources=['_denoise_cy.c'],
include_dirs=[get_numpy_include_dirs(), '../_shared'])
config.add_extension('rank.core8_cy', sources=['rank/core8_cy.c'],
config.add_extension('rank.core_cy', sources=['rank/core_cy.c'],
include_dirs=[get_numpy_include_dirs()])
config.add_extension('rank.core16_cy', sources=['rank/core16_cy.c'],
include_dirs=[get_numpy_include_dirs()])
config.add_extension('rank.generic8_cy', sources=['rank/generic8_cy.c'],
config.add_extension('rank.generic_cy', sources=['rank/generic_cy.c'],
include_dirs=[get_numpy_include_dirs()])
config.add_extension(
'rank.percentile8_cy', sources=['rank/percentile8_cy.c'],
include_dirs=[get_numpy_include_dirs()])
config.add_extension('rank.generic16_cy', sources=['rank/generic16_cy.c'],
'rank.percentile_cy', sources=['rank/percentile_cy.c'],
include_dirs=[get_numpy_include_dirs()])
config.add_extension(
'rank.percentile16_cy', sources=['rank/percentile16_cy.c'],
include_dirs=[get_numpy_include_dirs()])
config.add_extension(
'rank.bilateral16_cy', sources=['rank/bilateral16_cy.c'],
'rank.bilateral_cy', sources=['rank/bilateral_cy.c'],
include_dirs=[get_numpy_include_dirs()])
return config