wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-07-02 17:40:31 +08:00
Code Issues Packages Projects Releases Wiki Activity

fix pxd-pyx confusion

Browse Source
This commit is contained in:
odebeir
2012-10-13 09:47:25 +02:00
parent 894cb13f50
commit e83a456ebc
11 changed files with 1411 additions and 1324 deletions
Download Patch File Download Diff File Expand all files Collapse all files
skimage/rank/_core16.pxd
+1 -272
View File
@@ -1,18 +1,4 @@
""" to compile this use:
>>> python setup.py build_ext --inplace
to generate html report use:
>>> cython -a core16.pxd
"""
#cython: cdivision=True
#cython: boundscheck=False
#cython: nonecheck=False
#cython: wraparound=False
import numpy as np
cimport numpy as np
from libc.stdlib cimport malloc, free
#---------------------------------------------------------------------------
# 16 bit core kernel receives extra information about data bitdepth
@@ -23,261 +9,4 @@ np.ndarray[np.uint16_t, ndim=2] image,
np.ndarray[np.uint8_t, ndim=2] selem,
np.ndarray[np.uint8_t, ndim=2] mask,
np.ndarray[np.uint16_t, ndim=2] out,
char shift_x, char shift_y,int bitdepth):
""" Main loop, this function computes the histogram for each image point
- data is uint8
- result is uint8 casted
"""
cdef int rows = image.shape[0]
cdef int cols = image.shape[1]
cdef int srows = selem.shape[0]
cdef int scols = selem.shape[1]
cdef int centre_r = int(selem.shape[0] / 2) + shift_y
cdef int 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
assert bitdepth in range(2,13)
maxbin_list = [0,0,4,8,16,32,64,128,256,512,1024,2048,4096]
midbin_list = [0,0,2,4,8,16,32,64,128,256,512,1024,2048]
#set maxbin and midbin
cdef int maxbin=maxbin_list[bitdepth],midbin=midbin_list[bitdepth]
assert (image<maxbin).all()
image = np.ascontiguousarray(image)
if mask is None:
mask = np.ones((rows, cols), dtype=np.uint8)
else:
mask = np.ascontiguousarray(mask)
if out is None:
out = np.zeros((rows, cols), dtype=np.uint16)
else:
out = np.ascontiguousarray(out)
# create extended image and mask
cdef int erows = rows+srows-1
cdef int ecols = cols+scols-1
cdef np.ndarray emask = np.zeros((erows, ecols), dtype=np.uint8)
cdef np.ndarray eimage = np.zeros((erows, ecols), dtype=np.uint16)
eimage[centre_r:rows+centre_r,centre_c:cols+centre_c] = image
emask[centre_r:rows+centre_r,centre_c:cols+centre_c] = mask
mask = np.ascontiguousarray(mask)
# define pointers to the data
cdef np.uint16_t* eimage_data = <np.uint16_t*>eimage.data
cdef np.uint8_t* emask_data = <np.uint8_t*>emask.data
cdef np.uint16_t* out_data = <np.uint16_t*>out.data
cdef np.uint16_t* image_data = <np.uint16_t*>image.data
cdef np.uint8_t* mask_data = <np.uint8_t*>mask.data
# define local variable types
cdef int r, c, rr, cc, s, value, local_max, i, even_row
cdef float pop # number of pixels actually inside the neighborhood (float)
# allocate memory with malloc
cdef int max_se = srows*scols
# number of element in each attack border
cdef int num_se_n, num_se_s, num_se_e, num_se_w
# the current local histogram distribution
cdef int* histo = <int*>malloc(maxbin * sizeof(int))
# 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 int* se_e_r = <int*>malloc(max_se * sizeof(int))
cdef int* se_e_c = <int*>malloc(max_se * sizeof(int))
cdef int* se_w_r = <int*>malloc(max_se * sizeof(int))
cdef int* se_w_c = <int*>malloc(max_se * sizeof(int))
cdef int* se_n_r = <int*>malloc(max_se * sizeof(int))
cdef int* se_n_c = <int*>malloc(max_se * sizeof(int))
cdef int* se_s_r = <int*>malloc(max_se * sizeof(int))
cdef int* se_s_c = <int*>malloc(max_se * sizeof(int))
# build attack and release borders
# by using difference along axis
t = np.hstack((selem,np.zeros((selem.shape[0],1))))
t_e = np.diff(t,axis=1)==-1
t = np.hstack((np.zeros((selem.shape[0],1)),selem))
t_w = np.diff(t,axis=1)==1
t = np.vstack((selem,np.zeros((1,selem.shape[1]))))
t_s = np.diff(t,axis=0)==-1
t = np.vstack((np.zeros((1,selem.shape[1])),selem))
t_n = np.diff(t,axis=0)==1
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
cc = c
if selem[r, c]:
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
r = 0
c = 0
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],
bitdepth,maxbin,midbin)
# 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] + centre_r
cc = c + se_e_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
for s in range(num_se_w):
rr = r + se_w_r[s] + centre_r
cc = c + se_w_c[s] + centre_c - 1
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] -= 1
pop -= 1.
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],
bitdepth,maxbin,midbin)
# 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] + centre_r
cc = c + se_s_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
for s in range(num_se_n):
rr = r + se_n_r[s] + centre_r - 1
cc = c + se_n_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] -= 1
pop -= 1.
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],
bitdepth,maxbin,midbin)
# 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] + centre_r
cc = c + se_w_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
for s in range(num_se_e):
rr = r + se_e_r[s] + centre_r
cc = c + se_e_c[s] + centre_c + 1
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] -= 1
pop -= 1.
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],
bitdepth,maxbin,midbin)
# 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] + centre_r
cc = c + se_s_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
for s in range(num_se_n):
rr = r + se_n_r[s] + centre_r - 1
cc = c + se_n_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] -= 1
pop -= 1.
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],
bitdepth,maxbin,midbin)
# 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)
return out
char shift_x, char shift_y,int bitdepth)
skimage/rank/_core16.pyx
+283
View File
@@ -0,0 +1,283 @@
""" to compile this use:
>>> python setup.py build_ext --inplace
to generate html report use:
>>> cython -a core16.pxd
"""
#cython: cdivision=True
#cython: boundscheck=False
#cython: nonecheck=False
#cython: wraparound=False
import numpy as np
cimport numpy as np
from libc.stdlib cimport malloc, free
#---------------------------------------------------------------------------
# 16 bit core kernel receives extra information about data bitdepth
#---------------------------------------------------------------------------
cdef inline _core16(np.uint16_t kernel(int*, float, np.uint16_t, int ,int,int ),
np.ndarray[np.uint16_t, ndim=2] image,
np.ndarray[np.uint8_t, ndim=2] selem,
np.ndarray[np.uint8_t, ndim=2] mask,
np.ndarray[np.uint16_t, ndim=2] out,
char shift_x, char shift_y,int bitdepth):
""" Main loop, this function computes the histogram for each image point
- data is uint8
- result is uint8 casted
"""
cdef int rows = image.shape[0]
cdef int cols = image.shape[1]
cdef int srows = selem.shape[0]
cdef int scols = selem.shape[1]
cdef int centre_r = int(selem.shape[0] / 2) + shift_y
cdef int 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
assert bitdepth in range(2,13)
maxbin_list = [0,0,4,8,16,32,64,128,256,512,1024,2048,4096]
midbin_list = [0,0,2,4,8,16,32,64,128,256,512,1024,2048]
#set maxbin and midbin
cdef int maxbin=maxbin_list[bitdepth],midbin=midbin_list[bitdepth]
assert (image<maxbin).all()
image = np.ascontiguousarray(image)
if mask is None:
mask = np.ones((rows, cols), dtype=np.uint8)
else:
mask = np.ascontiguousarray(mask)
if out is None:
out = np.zeros((rows, cols), dtype=np.uint16)
else:
out = np.ascontiguousarray(out)
# create extended image and mask
cdef int erows = rows+srows-1
cdef int ecols = cols+scols-1
cdef np.ndarray emask = np.zeros((erows, ecols), dtype=np.uint8)
cdef np.ndarray eimage = np.zeros((erows, ecols), dtype=np.uint16)
eimage[centre_r:rows+centre_r,centre_c:cols+centre_c] = image
emask[centre_r:rows+centre_r,centre_c:cols+centre_c] = mask
mask = np.ascontiguousarray(mask)
# define pointers to the data
cdef np.uint16_t* eimage_data = <np.uint16_t*>eimage.data
cdef np.uint8_t* emask_data = <np.uint8_t*>emask.data
cdef np.uint16_t* out_data = <np.uint16_t*>out.data
cdef np.uint16_t* image_data = <np.uint16_t*>image.data
cdef np.uint8_t* mask_data = <np.uint8_t*>mask.data
# define local variable types
cdef int r, c, rr, cc, s, value, local_max, i, even_row
cdef float pop # number of pixels actually inside the neighborhood (float)
# allocate memory with malloc
cdef int max_se = srows*scols
# number of element in each attack border
cdef int num_se_n, num_se_s, num_se_e, num_se_w
# the current local histogram distribution
cdef int* histo = <int*>malloc(maxbin * sizeof(int))
# 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 int* se_e_r = <int*>malloc(max_se * sizeof(int))
cdef int* se_e_c = <int*>malloc(max_se * sizeof(int))
cdef int* se_w_r = <int*>malloc(max_se * sizeof(int))
cdef int* se_w_c = <int*>malloc(max_se * sizeof(int))
cdef int* se_n_r = <int*>malloc(max_se * sizeof(int))
cdef int* se_n_c = <int*>malloc(max_se * sizeof(int))
cdef int* se_s_r = <int*>malloc(max_se * sizeof(int))
cdef int* se_s_c = <int*>malloc(max_se * sizeof(int))
# build attack and release borders
# by using difference along axis
t = np.hstack((selem,np.zeros((selem.shape[0],1))))
t_e = np.diff(t,axis=1)==-1
t = np.hstack((np.zeros((selem.shape[0],1)),selem))
t_w = np.diff(t,axis=1)==1
t = np.vstack((selem,np.zeros((1,selem.shape[1]))))
t_s = np.diff(t,axis=0)==-1
t = np.vstack((np.zeros((1,selem.shape[1])),selem))
t_n = np.diff(t,axis=0)==1
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
cc = c
if selem[r, c]:
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
r = 0
c = 0
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],
bitdepth,maxbin,midbin)
# 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] + centre_r
cc = c + se_e_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
for s in range(num_se_w):
rr = r + se_w_r[s] + centre_r
cc = c + se_w_c[s] + centre_c - 1
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] -= 1
pop -= 1.
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],
bitdepth,maxbin,midbin)
# 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] + centre_r
cc = c + se_s_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
for s in range(num_se_n):
rr = r + se_n_r[s] + centre_r - 1
cc = c + se_n_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] -= 1
pop -= 1.
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],
bitdepth,maxbin,midbin)
# 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] + centre_r
cc = c + se_w_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
for s in range(num_se_e):
rr = r + se_e_r[s] + centre_r
cc = c + se_e_c[s] + centre_c + 1
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] -= 1
pop -= 1.
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],
bitdepth,maxbin,midbin)
# 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] + centre_r
cc = c + se_s_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
for s in range(num_se_n):
rr = r + se_n_r[s] + centre_r - 1
cc = c + se_n_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] -= 1
pop -= 1.
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],
bitdepth,maxbin,midbin)
# 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)
return out
skimage/rank/_core16b.pxd
+1 -273
View File
@@ -1,18 +1,4 @@
""" to compile this use:
>>> python setup.py build_ext --inplace
to generate html report use:
>>> cython -a core16b.pxd
"""
#cython: cdivision=True
#cython: boundscheck=False
#cython: nonecheck=False
#cython: wraparound=False
import numpy as np
cimport numpy as np
from libc.stdlib cimport malloc, free
#---------------------------------------------------------------------------
# 16 bit core kernel receives extra information about data bitdepth and bilateral interval
@@ -23,262 +9,4 @@ np.ndarray[np.uint16_t, ndim=2] image,
np.ndarray[np.uint8_t, ndim=2] selem,
np.ndarray[np.uint8_t, ndim=2] mask,
np.ndarray[np.uint16_t, ndim=2] out,
char shift_x, char shift_y,int bitdepth, int s0, int s1):
""" Main loop, this function computes the histogram for each image point
- data is uint8
- result is uint8 casted
- only pixel inside [s0,s1] centered on g are taken into account
"""
cdef int rows = image.shape[0]
cdef int cols = image.shape[1]
cdef int srows = selem.shape[0]
cdef int scols = selem.shape[1]
cdef int centre_r = int(selem.shape[0] / 2) + shift_y
cdef int 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
assert bitdepth in range(2,13)
maxbin_list = [0,0,4,8,16,32,64,128,256,512,1024,2048,4096]
midbin_list = [0,0,2,4,8,16,32,64,128,256,512,1024,2048]
#set maxbin and midbin
cdef int maxbin=maxbin_list[bitdepth],midbin=midbin_list[bitdepth]
assert (image<maxbin).all()
image = np.ascontiguousarray(image)
if mask is None:
mask = np.ones((rows, cols), dtype=np.uint8)
else:
mask = np.ascontiguousarray(mask)
if out is None:
out = np.zeros((rows, cols), dtype=np.uint16)
else:
out = np.ascontiguousarray(out)
# create extended image and mask
cdef int erows = rows+srows-1
cdef int ecols = cols+scols-1
cdef np.ndarray emask = np.zeros((erows, ecols), dtype=np.uint8)
cdef np.ndarray eimage = np.zeros((erows, ecols), dtype=np.uint16)
eimage[centre_r:rows+centre_r,centre_c:cols+centre_c] = image
emask[centre_r:rows+centre_r,centre_c:cols+centre_c] = mask
mask = np.ascontiguousarray(mask)
# define pointers to the data
cdef np.uint16_t* eimage_data = <np.uint16_t*>eimage.data
cdef np.uint8_t* emask_data = <np.uint8_t*>emask.data
cdef np.uint16_t* out_data = <np.uint16_t*>out.data
cdef np.uint16_t* image_data = <np.uint16_t*>image.data
cdef np.uint8_t* mask_data = <np.uint8_t*>mask.data
# define local variable types
cdef int r, c, rr, cc, s, value, local_max, i, even_row
cdef float pop # number of pixels actually inside the neighborhood (float)
# allocate memory with malloc
cdef int max_se = srows*scols
# number of element in each attack border
cdef int num_se_n, num_se_s, num_se_e, num_se_w
# the current local histogram distribution
cdef int* histo = <int*>malloc(maxbin * sizeof(int))
# 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 int* se_e_r = <int*>malloc(max_se * sizeof(int))
cdef int* se_e_c = <int*>malloc(max_se * sizeof(int))
cdef int* se_w_r = <int*>malloc(max_se * sizeof(int))
cdef int* se_w_c = <int*>malloc(max_se * sizeof(int))
cdef int* se_n_r = <int*>malloc(max_se * sizeof(int))
cdef int* se_n_c = <int*>malloc(max_se * sizeof(int))
cdef int* se_s_r = <int*>malloc(max_se * sizeof(int))
cdef int* se_s_c = <int*>malloc(max_se * sizeof(int))
# build attack and release borders
# by using difference along axis
t = np.hstack((selem,np.zeros((selem.shape[0],1))))
t_e = np.diff(t,axis=1)==-1
t = np.hstack((np.zeros((selem.shape[0],1)),selem))
t_w = np.diff(t,axis=1)==1
t = np.vstack((selem,np.zeros((1,selem.shape[1]))))
t_s = np.diff(t,axis=0)==-1
t = np.vstack((np.zeros((1,selem.shape[1])),selem))
t_n = np.diff(t,axis=0)==1
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
cc = c
if selem[r, c]:
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
r = 0
c = 0
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],
bitdepth,maxbin,midbin,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] + centre_r
cc = c + se_e_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
for s in range(num_se_w):
rr = r + se_w_r[s] + centre_r
cc = c + se_w_c[s] + centre_c - 1
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] -= 1
pop -= 1.
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],
bitdepth,maxbin,midbin,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] + centre_r
cc = c + se_s_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
for s in range(num_se_n):
rr = r + se_n_r[s] + centre_r - 1
cc = c + se_n_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] -= 1
pop -= 1.
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],
bitdepth,maxbin,midbin,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] + centre_r
cc = c + se_w_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
for s in range(num_se_e):
rr = r + se_e_r[s] + centre_r
cc = c + se_e_c[s] + centre_c + 1
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] -= 1
pop -= 1.
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],
bitdepth,maxbin,midbin,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] + centre_r
cc = c + se_s_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
for s in range(num_se_n):
rr = r + se_n_r[s] + centre_r - 1
cc = c + se_n_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] -= 1
pop -= 1.
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],
bitdepth,maxbin,midbin,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)
return out
char shift_x, char shift_y,int bitdepth, int s0, int s1)
skimage/rank/_core16b.pyx
+284
View File
@@ -0,0 +1,284 @@
""" to compile this use:
>>> python setup.py build_ext --inplace
to generate html report use:
>>> cython -a core16b.pxd
"""
#cython: cdivision=True
#cython: boundscheck=False
#cython: nonecheck=False
#cython: wraparound=False
import numpy as np
cimport numpy as np
from libc.stdlib cimport malloc, free
#---------------------------------------------------------------------------
# 16 bit core kernel receives extra information about data bitdepth and bilateral interval
#---------------------------------------------------------------------------
cdef inline _core16b(np.uint16_t kernel(int*, float, np.uint16_t, int ,int,int,int,int),
np.ndarray[np.uint16_t, ndim=2] image,
np.ndarray[np.uint8_t, ndim=2] selem,
np.ndarray[np.uint8_t, ndim=2] mask,
np.ndarray[np.uint16_t, ndim=2] out,
char shift_x, char shift_y,int bitdepth, int s0, int s1):
""" Main loop, this function computes the histogram for each image point
- data is uint8
- result is uint8 casted
- only pixel inside [s0,s1] centered on g are taken into account
"""
cdef int rows = image.shape[0]
cdef int cols = image.shape[1]
cdef int srows = selem.shape[0]
cdef int scols = selem.shape[1]
cdef int centre_r = int(selem.shape[0] / 2) + shift_y
cdef int 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
assert bitdepth in range(2,13)
maxbin_list = [0,0,4,8,16,32,64,128,256,512,1024,2048,4096]
midbin_list = [0,0,2,4,8,16,32,64,128,256,512,1024,2048]
#set maxbin and midbin
cdef int maxbin=maxbin_list[bitdepth],midbin=midbin_list[bitdepth]
assert (image<maxbin).all()
image = np.ascontiguousarray(image)
if mask is None:
mask = np.ones((rows, cols), dtype=np.uint8)
else:
mask = np.ascontiguousarray(mask)
if out is None:
out = np.zeros((rows, cols), dtype=np.uint16)
else:
out = np.ascontiguousarray(out)
# create extended image and mask
cdef int erows = rows+srows-1
cdef int ecols = cols+scols-1
cdef np.ndarray emask = np.zeros((erows, ecols), dtype=np.uint8)
cdef np.ndarray eimage = np.zeros((erows, ecols), dtype=np.uint16)
eimage[centre_r:rows+centre_r,centre_c:cols+centre_c] = image
emask[centre_r:rows+centre_r,centre_c:cols+centre_c] = mask
mask = np.ascontiguousarray(mask)
# define pointers to the data
cdef np.uint16_t* eimage_data = <np.uint16_t*>eimage.data
cdef np.uint8_t* emask_data = <np.uint8_t*>emask.data
cdef np.uint16_t* out_data = <np.uint16_t*>out.data
cdef np.uint16_t* image_data = <np.uint16_t*>image.data
cdef np.uint8_t* mask_data = <np.uint8_t*>mask.data
# define local variable types
cdef int r, c, rr, cc, s, value, local_max, i, even_row
cdef float pop # number of pixels actually inside the neighborhood (float)
# allocate memory with malloc
cdef int max_se = srows*scols
# number of element in each attack border
cdef int num_se_n, num_se_s, num_se_e, num_se_w
# the current local histogram distribution
cdef int* histo = <int*>malloc(maxbin * sizeof(int))
# 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 int* se_e_r = <int*>malloc(max_se * sizeof(int))
cdef int* se_e_c = <int*>malloc(max_se * sizeof(int))
cdef int* se_w_r = <int*>malloc(max_se * sizeof(int))
cdef int* se_w_c = <int*>malloc(max_se * sizeof(int))
cdef int* se_n_r = <int*>malloc(max_se * sizeof(int))
cdef int* se_n_c = <int*>malloc(max_se * sizeof(int))
cdef int* se_s_r = <int*>malloc(max_se * sizeof(int))
cdef int* se_s_c = <int*>malloc(max_se * sizeof(int))
# build attack and release borders
# by using difference along axis
t = np.hstack((selem,np.zeros((selem.shape[0],1))))
t_e = np.diff(t,axis=1)==-1
t = np.hstack((np.zeros((selem.shape[0],1)),selem))
t_w = np.diff(t,axis=1)==1
t = np.vstack((selem,np.zeros((1,selem.shape[1]))))
t_s = np.diff(t,axis=0)==-1
t = np.vstack((np.zeros((1,selem.shape[1])),selem))
t_n = np.diff(t,axis=0)==1
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
cc = c
if selem[r, c]:
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
r = 0
c = 0
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],
bitdepth,maxbin,midbin,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] + centre_r
cc = c + se_e_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
for s in range(num_se_w):
rr = r + se_w_r[s] + centre_r
cc = c + se_w_c[s] + centre_c - 1
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] -= 1
pop -= 1.
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],
bitdepth,maxbin,midbin,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] + centre_r
cc = c + se_s_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
for s in range(num_se_n):
rr = r + se_n_r[s] + centre_r - 1
cc = c + se_n_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] -= 1
pop -= 1.
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],
bitdepth,maxbin,midbin,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] + centre_r
cc = c + se_w_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
for s in range(num_se_e):
rr = r + se_e_r[s] + centre_r
cc = c + se_e_c[s] + centre_c + 1
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] -= 1
pop -= 1.
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],
bitdepth,maxbin,midbin,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] + centre_r
cc = c + se_s_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
for s in range(num_se_n):
rr = r + se_n_r[s] + centre_r - 1
cc = c + se_n_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] -= 1
pop -= 1.
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],
bitdepth,maxbin,midbin,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)
return out
skimage/rank/_core16p.pxd
+3 -269
View File
@@ -1,15 +1,8 @@
#cython: cdivision=True
#cython: boundscheck=False
#cython: nonecheck=False
#cython: wraparound=False
import numpy as np
cimport numpy as np
from libc.stdlib cimport malloc, free
# generic cdef functions
cdef inline int int_max(int a, int b): return a if a >= b else b
cdef inline int int_min(int a, int b): return a if a <= b else b
cdef inline int int_max(int a, int b)
cdef inline int int_min(int a, int b)
#---------------------------------------------------------------------------
# 16 bit core kernel receives extra information about data inferior and superior percentiles
@@ -20,263 +13,4 @@ np.ndarray[np.uint16_t, ndim=2] image,
np.ndarray[np.uint8_t, ndim=2] selem,
np.ndarray[np.uint8_t, ndim=2] mask,
np.ndarray[np.uint16_t, ndim=2] out,
char shift_x, char shift_y,int bitdepth, float p0, float p1):
""" Main loop, this function computes the histogram for each image point
- data is uint16
- result is uint16 casted
"""
cdef int rows = image.shape[0]
cdef int cols = image.shape[1]
cdef int srows = selem.shape[0]
cdef int scols = selem.shape[1]
cdef int centre_r = int(selem.shape[0] / 2) + shift_y
cdef int 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
assert bitdepth in range(2,13)
maxbin_list = [0,0,4,8,16,32,64,128,256,512,1024,2048,4096]
midbin_list = [0,0,2,4,8,16,32,64,128,256,512,1024,2048]
#set maxbin and midbin
cdef int maxbin=maxbin_list[bitdepth],midbin=midbin_list[bitdepth]
assert (image<maxbin).all()
image = np.ascontiguousarray(image)
if mask is None:
mask = np.ones((rows, cols), dtype=np.uint8)
else:
mask = np.ascontiguousarray(mask)
if out is None:
out = np.zeros((rows, cols), dtype=np.uint16)
else:
out = np.ascontiguousarray(out)
# create extended image and mask
cdef int erows = rows+srows-1
cdef int ecols = cols+scols-1
cdef np.ndarray emask = np.zeros((erows, ecols), dtype=np.uint8)
cdef np.ndarray eimage = np.zeros((erows, ecols), dtype=np.uint16)
eimage[centre_r:rows+centre_r,centre_c:cols+centre_c] = image
emask[centre_r:rows+centre_r,centre_c:cols+centre_c] = mask
mask = np.ascontiguousarray(mask)
# define pointers to the data
cdef np.uint16_t* eimage_data = <np.uint16_t*>eimage.data
cdef np.uint8_t* emask_data = <np.uint8_t*>emask.data
cdef np.uint16_t* out_data = <np.uint16_t*>out.data
cdef np.uint16_t* image_data = <np.uint16_t*>image.data
cdef np.uint8_t* mask_data = <np.uint8_t*>mask.data
# define local variable types
cdef int r, c, rr, cc, s, value, local_max, i, even_row
cdef float pop # number of pixels actually inside the neighborhood (float)
# allocate memory with malloc
cdef int max_se = srows*scols
# number of element in each attack border
cdef int num_se_n, num_se_s, num_se_e, num_se_w
# the current local histogram distribution
cdef int* histo = <int*>malloc(maxbin * sizeof(int))
# 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 int* se_e_r = <int*>malloc(max_se * sizeof(int))
cdef int* se_e_c = <int*>malloc(max_se * sizeof(int))
cdef int* se_w_r = <int*>malloc(max_se * sizeof(int))
cdef int* se_w_c = <int*>malloc(max_se * sizeof(int))
cdef int* se_n_r = <int*>malloc(max_se * sizeof(int))
cdef int* se_n_c = <int*>malloc(max_se * sizeof(int))
cdef int* se_s_r = <int*>malloc(max_se * sizeof(int))
cdef int* se_s_c = <int*>malloc(max_se * sizeof(int))
# build attack and release borders
# by using difference along axis
t = np.hstack((selem,np.zeros((selem.shape[0],1))))
t_e = np.diff(t,axis=1)==-1
t = np.hstack((np.zeros((selem.shape[0],1)),selem))
t_w = np.diff(t,axis=1)==1
t = np.vstack((selem,np.zeros((1,selem.shape[1]))))
t_s = np.diff(t,axis=0)==-1
t = np.vstack((np.zeros((1,selem.shape[1])),selem))
t_n = np.diff(t,axis=0)==1
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
cc = c
if selem[r, c]:
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
r = 0
c = 0
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],
bitdepth,maxbin,midbin,p0,p1)
# 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] + centre_r
cc = c + se_e_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
for s in range(num_se_w):
rr = r + se_w_r[s] + centre_r
cc = c + se_w_c[s] + centre_c - 1
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] -= 1
pop -= 1.
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],
bitdepth,maxbin,midbin,p0,p1)
# 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] + centre_r
cc = c + se_s_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
for s in range(num_se_n):
rr = r + se_n_r[s] + centre_r - 1
cc = c + se_n_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] -= 1
pop -= 1.
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],
bitdepth,maxbin,midbin,p0,p1)
# 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] + centre_r
cc = c + se_w_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
for s in range(num_se_e):
rr = r + se_e_r[s] + centre_r
cc = c + se_e_c[s] + centre_c + 1
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] -= 1
pop -= 1.
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],
bitdepth,maxbin,midbin,p0,p1)
# 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] + centre_r
cc = c + se_s_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
for s in range(num_se_n):
rr = r + se_n_r[s] + centre_r - 1
cc = c + se_n_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] -= 1
pop -= 1.
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],
bitdepth,maxbin,midbin,p0,p1)
# 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)
return out
char shift_x, char shift_y,int bitdepth, float p0, float p1)
skimage/rank/_core16p.pyx
+282
View File
@@ -0,0 +1,282 @@
#cython: cdivision=True
#cython: boundscheck=False
#cython: nonecheck=False
#cython: wraparound=False
import numpy as np
cimport numpy as np
from libc.stdlib cimport malloc, free
# generic cdef functions
cdef inline int int_max(int a, int b): return a if a >= b else b
cdef inline int int_min(int a, int b): return a if a <= b else b
#---------------------------------------------------------------------------
# 16 bit core kernel receives extra information about data inferior and superior percentiles
#---------------------------------------------------------------------------
cdef inline _core16p(np.uint16_t kernel(int*, float, np.uint16_t,int,int,int, float, float),
np.ndarray[np.uint16_t, ndim=2] image,
np.ndarray[np.uint8_t, ndim=2] selem,
np.ndarray[np.uint8_t, ndim=2] mask,
np.ndarray[np.uint16_t, ndim=2] out,
char shift_x, char shift_y,int bitdepth, float p0, float p1):
""" Main loop, this function computes the histogram for each image point
- data is uint16
- result is uint16 casted
"""
cdef int rows = image.shape[0]
cdef int cols = image.shape[1]
cdef int srows = selem.shape[0]
cdef int scols = selem.shape[1]
cdef int centre_r = int(selem.shape[0] / 2) + shift_y
cdef int 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
assert bitdepth in range(2,13)
maxbin_list = [0,0,4,8,16,32,64,128,256,512,1024,2048,4096]
midbin_list = [0,0,2,4,8,16,32,64,128,256,512,1024,2048]
#set maxbin and midbin
cdef int maxbin=maxbin_list[bitdepth],midbin=midbin_list[bitdepth]
assert (image<maxbin).all()
image = np.ascontiguousarray(image)
if mask is None:
mask = np.ones((rows, cols), dtype=np.uint8)
else:
mask = np.ascontiguousarray(mask)
if out is None:
out = np.zeros((rows, cols), dtype=np.uint16)
else:
out = np.ascontiguousarray(out)
# create extended image and mask
cdef int erows = rows+srows-1
cdef int ecols = cols+scols-1
cdef np.ndarray emask = np.zeros((erows, ecols), dtype=np.uint8)
cdef np.ndarray eimage = np.zeros((erows, ecols), dtype=np.uint16)
eimage[centre_r:rows+centre_r,centre_c:cols+centre_c] = image
emask[centre_r:rows+centre_r,centre_c:cols+centre_c] = mask
mask = np.ascontiguousarray(mask)
# define pointers to the data
cdef np.uint16_t* eimage_data = <np.uint16_t*>eimage.data
cdef np.uint8_t* emask_data = <np.uint8_t*>emask.data
cdef np.uint16_t* out_data = <np.uint16_t*>out.data
cdef np.uint16_t* image_data = <np.uint16_t*>image.data
cdef np.uint8_t* mask_data = <np.uint8_t*>mask.data
# define local variable types
cdef int r, c, rr, cc, s, value, local_max, i, even_row
cdef float pop # number of pixels actually inside the neighborhood (float)
# allocate memory with malloc
cdef int max_se = srows*scols
# number of element in each attack border
cdef int num_se_n, num_se_s, num_se_e, num_se_w
# the current local histogram distribution
cdef int* histo = <int*>malloc(maxbin * sizeof(int))
# 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 int* se_e_r = <int*>malloc(max_se * sizeof(int))
cdef int* se_e_c = <int*>malloc(max_se * sizeof(int))
cdef int* se_w_r = <int*>malloc(max_se * sizeof(int))
cdef int* se_w_c = <int*>malloc(max_se * sizeof(int))
cdef int* se_n_r = <int*>malloc(max_se * sizeof(int))
cdef int* se_n_c = <int*>malloc(max_se * sizeof(int))
cdef int* se_s_r = <int*>malloc(max_se * sizeof(int))
cdef int* se_s_c = <int*>malloc(max_se * sizeof(int))
# build attack and release borders
# by using difference along axis
t = np.hstack((selem,np.zeros((selem.shape[0],1))))
t_e = np.diff(t,axis=1)==-1
t = np.hstack((np.zeros((selem.shape[0],1)),selem))
t_w = np.diff(t,axis=1)==1
t = np.vstack((selem,np.zeros((1,selem.shape[1]))))
t_s = np.diff(t,axis=0)==-1
t = np.vstack((np.zeros((1,selem.shape[1])),selem))
t_n = np.diff(t,axis=0)==1
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
cc = c
if selem[r, c]:
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
r = 0
c = 0
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],
bitdepth,maxbin,midbin,p0,p1)
# 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] + centre_r
cc = c + se_e_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
for s in range(num_se_w):
rr = r + se_w_r[s] + centre_r
cc = c + se_w_c[s] + centre_c - 1
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] -= 1
pop -= 1.
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],
bitdepth,maxbin,midbin,p0,p1)
# 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] + centre_r
cc = c + se_s_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
for s in range(num_se_n):
rr = r + se_n_r[s] + centre_r - 1
cc = c + se_n_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] -= 1
pop -= 1.
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],
bitdepth,maxbin,midbin,p0,p1)
# 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] + centre_r
cc = c + se_w_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
for s in range(num_se_e):
rr = r + se_e_r[s] + centre_r
cc = c + se_e_c[s] + centre_c + 1
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] -= 1
pop -= 1.
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],
bitdepth,maxbin,midbin,p0,p1)
# 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] + centre_r
cc = c + se_s_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
for s in range(num_se_n):
rr = r + se_n_r[s] + centre_r - 1
cc = c + se_n_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] -= 1
pop -= 1.
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],
bitdepth,maxbin,midbin,p0,p1)
# 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)
return out
skimage/rank/_core8.pxd
+1 -258
View File
@@ -1,18 +1,4 @@
""" to compile this use:
>>> python setup.py build_ext --inplace
to generate html report use:
>>> cython -a core8.pxd
"""
#cython: cdivision=True
#cython: boundscheck=False
#cython: nonecheck=False
#cython: wraparound=False
import numpy as np
cimport numpy as np
from libc.stdlib cimport malloc, free
#---------------------------------------------------------------------------
# 8 bit core kernel
@@ -23,247 +9,4 @@ np.ndarray[np.uint8_t, ndim=2] image,
np.ndarray[np.uint8_t, ndim=2] selem,
np.ndarray[np.uint8_t, ndim=2] mask,
np.ndarray[np.uint8_t, ndim=2] out,
char shift_x, char shift_y):
""" Main loop, this function computes the histogram for each image point
- data is uint8
- result is uint8 casted
"""
cdef int rows = image.shape[0]
cdef int cols = image.shape[1]
cdef int srows = selem.shape[0]
cdef int scols = selem.shape[1]
cdef int centre_r = int(selem.shape[0] / 2) + shift_y
cdef int 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
image = np.ascontiguousarray(image)
if mask is None:
mask = np.ones((rows, cols), dtype=np.uint8)
else:
mask = np.ascontiguousarray(mask)
if out is None:
out = np.zeros((rows, cols), dtype=np.uint8)
else:
out = np.ascontiguousarray(out)
# create extended image and mask
cdef int erows = rows+srows-1
cdef int ecols = cols+scols-1
cdef np.ndarray emask = np.zeros((erows, ecols), dtype=np.uint8)
cdef np.ndarray eimage = np.zeros((erows, ecols), dtype=np.uint8)
eimage[centre_r:rows+centre_r,centre_c:cols+centre_c] = image
emask[centre_r:rows+centre_r,centre_c:cols+centre_c] = mask
mask = np.ascontiguousarray(mask)
# define pointers to the data
cdef np.uint8_t* eimage_data = <np.uint8_t*>eimage.data
cdef np.uint8_t* emask_data = <np.uint8_t*>emask.data
cdef np.uint8_t* out_data = <np.uint8_t*>out.data
cdef np.uint8_t* image_data = <np.uint8_t*>image.data
cdef np.uint8_t* mask_data = <np.uint8_t*>mask.data
# define local variable types
cdef int r, c, rr, cc, s, value, local_max, i, even_row
cdef float pop # number of pixels actually inside the neighborhood (float)
# allocate memory with malloc
cdef int max_se = srows*scols
# number of element in each attack border
cdef int num_se_n, num_se_s, num_se_e, num_se_w
# the current local histogram distribution
cdef int* histo = <int*>malloc(256 * sizeof(int))
# 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 int* se_e_r = <int*>malloc(max_se * sizeof(int))
cdef int* se_e_c = <int*>malloc(max_se * sizeof(int))
cdef int* se_w_r = <int*>malloc(max_se * sizeof(int))
cdef int* se_w_c = <int*>malloc(max_se * sizeof(int))
cdef int* se_n_r = <int*>malloc(max_se * sizeof(int))
cdef int* se_n_c = <int*>malloc(max_se * sizeof(int))
cdef int* se_s_r = <int*>malloc(max_se * sizeof(int))
cdef int* se_s_c = <int*>malloc(max_se * sizeof(int))
# build attack and release borders
# by using difference along axis
t = np.hstack((selem,np.zeros((selem.shape[0],1))))
t_e = np.diff(t,axis=1)==-1
t = np.hstack((np.zeros((selem.shape[0],1)),selem))
t_w = np.diff(t,axis=1)==1
t = np.vstack((selem,np.zeros((1,selem.shape[1]))))
t_s = np.diff(t,axis=0)==-1
t = np.vstack((np.zeros((1,selem.shape[1])),selem))
t_n = np.diff(t,axis=0)==1
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(256):
histo[i] = 0
pop = 0
for r in range(srows):
for c in range(scols):
rr = r
cc = c
if selem[r, c]:
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
r = 0
c = 0
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c])
# 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] + centre_r
cc = c + se_e_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
for s in range(num_se_w):
rr = r + se_w_r[s] + centre_r
cc = c + se_w_c[s] + centre_c - 1
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] -= 1
pop -= 1.
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c])
# 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] + centre_r
cc = c + se_s_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
for s in range(num_se_n):
rr = r + se_n_r[s] + centre_r - 1
cc = c + se_n_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] -= 1
pop -= 1.
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c])
# 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] + centre_r
cc = c + se_w_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
for s in range(num_se_e):
rr = r + se_e_r[s] + centre_r
cc = c + se_e_c[s] + centre_c + 1
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] -= 1
pop -= 1.
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c])
# 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] + centre_r
cc = c + se_s_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
for s in range(num_se_n):
rr = r + se_n_r[s] + centre_r - 1
cc = c + se_n_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] -= 1
pop -= 1.
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c])
# 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)
return out
char shift_x, char shift_y)
skimage/rank/_core8.pyx
+269
View File
@@ -0,0 +1,269 @@
""" to compile this use:
>>> python setup.py build_ext --inplace
to generate html report use:
>>> cython -a core8.pxd
"""
#cython: cdivision=True
#cython: boundscheck=False
#cython: nonecheck=False
#cython: wraparound=False
import numpy as np
cimport numpy as np
from libc.stdlib cimport malloc, free
#---------------------------------------------------------------------------
# 8 bit core kernel
#---------------------------------------------------------------------------
cdef inline _core8(np.uint8_t kernel(int*, float, np.uint8_t),
np.ndarray[np.uint8_t, ndim=2] image,
np.ndarray[np.uint8_t, ndim=2] selem,
np.ndarray[np.uint8_t, ndim=2] mask,
np.ndarray[np.uint8_t, ndim=2] out,
char shift_x, char shift_y):
""" Main loop, this function computes the histogram for each image point
- data is uint8
- result is uint8 casted
"""
cdef int rows = image.shape[0]
cdef int cols = image.shape[1]
cdef int srows = selem.shape[0]
cdef int scols = selem.shape[1]
cdef int centre_r = int(selem.shape[0] / 2) + shift_y
cdef int 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
image = np.ascontiguousarray(image)
if mask is None:
mask = np.ones((rows, cols), dtype=np.uint8)
else:
mask = np.ascontiguousarray(mask)
if out is None:
out = np.zeros((rows, cols), dtype=np.uint8)
else:
out = np.ascontiguousarray(out)
# create extended image and mask
cdef int erows = rows+srows-1
cdef int ecols = cols+scols-1
cdef np.ndarray emask = np.zeros((erows, ecols), dtype=np.uint8)
cdef np.ndarray eimage = np.zeros((erows, ecols), dtype=np.uint8)
eimage[centre_r:rows+centre_r,centre_c:cols+centre_c] = image
emask[centre_r:rows+centre_r,centre_c:cols+centre_c] = mask
mask = np.ascontiguousarray(mask)
# define pointers to the data
cdef np.uint8_t* eimage_data = <np.uint8_t*>eimage.data
cdef np.uint8_t* emask_data = <np.uint8_t*>emask.data
cdef np.uint8_t* out_data = <np.uint8_t*>out.data
cdef np.uint8_t* image_data = <np.uint8_t*>image.data
cdef np.uint8_t* mask_data = <np.uint8_t*>mask.data
# define local variable types
cdef int r, c, rr, cc, s, value, local_max, i, even_row
cdef float pop # number of pixels actually inside the neighborhood (float)
# allocate memory with malloc
cdef int max_se = srows*scols
# number of element in each attack border
cdef int num_se_n, num_se_s, num_se_e, num_se_w
# the current local histogram distribution
cdef int* histo = <int*>malloc(256 * sizeof(int))
# 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 int* se_e_r = <int*>malloc(max_se * sizeof(int))
cdef int* se_e_c = <int*>malloc(max_se * sizeof(int))
cdef int* se_w_r = <int*>malloc(max_se * sizeof(int))
cdef int* se_w_c = <int*>malloc(max_se * sizeof(int))
cdef int* se_n_r = <int*>malloc(max_se * sizeof(int))
cdef int* se_n_c = <int*>malloc(max_se * sizeof(int))
cdef int* se_s_r = <int*>malloc(max_se * sizeof(int))
cdef int* se_s_c = <int*>malloc(max_se * sizeof(int))
# build attack and release borders
# by using difference along axis
t = np.hstack((selem,np.zeros((selem.shape[0],1))))
t_e = np.diff(t,axis=1)==-1
t = np.hstack((np.zeros((selem.shape[0],1)),selem))
t_w = np.diff(t,axis=1)==1
t = np.vstack((selem,np.zeros((1,selem.shape[1]))))
t_s = np.diff(t,axis=0)==-1
t = np.vstack((np.zeros((1,selem.shape[1])),selem))
t_n = np.diff(t,axis=0)==1
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(256):
histo[i] = 0
pop = 0
for r in range(srows):
for c in range(scols):
rr = r
cc = c
if selem[r, c]:
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
r = 0
c = 0
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c])
# 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] + centre_r
cc = c + se_e_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
for s in range(num_se_w):
rr = r + se_w_r[s] + centre_r
cc = c + se_w_c[s] + centre_c - 1
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] -= 1
pop -= 1.
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c])
# 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] + centre_r
cc = c + se_s_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
for s in range(num_se_n):
rr = r + se_n_r[s] + centre_r - 1
cc = c + se_n_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] -= 1
pop -= 1.
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c])
# 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] + centre_r
cc = c + se_w_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
for s in range(num_se_e):
rr = r + se_e_r[s] + centre_r
cc = c + se_e_c[s] + centre_c + 1
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] -= 1
pop -= 1.
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c])
# 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] + centre_r
cc = c + se_s_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
for s in range(num_se_n):
rr = r + se_n_r[s] + centre_r - 1
cc = c + se_n_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] -= 1
pop -= 1.
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c])
# 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)
return out
skimage/rank/_core8p.pxd
+3 -252
View File
@@ -1,15 +1,8 @@
#cython: cdivision=True
#cython: boundscheck=False
#cython: nonecheck=False
#cython: wraparound=False
import numpy as np
cimport numpy as np
from libc.stdlib cimport malloc, free
# generic cdef functions
cdef inline np.uint8_t uint8_max(np.uint8_t a, np.uint8_t b): return a if a >= b else b
cdef inline np.uint8_t uint8_min(np.uint8_t a, np.uint8_t b): return a if a <= b else b
cdef inline np.uint8_t uint8_max(np.uint8_t a, np.uint8_t b)
cdef inline np.uint8_t uint8_min(np.uint8_t a, np.uint8_t b)
#---------------------------------------------------------------------------
# 8 bit core kernel receives extra information about data inferior and superior percentiles
@@ -20,247 +13,5 @@ np.ndarray[np.uint8_t, ndim=2] image,
np.ndarray[np.uint8_t, ndim=2] selem,
np.ndarray[np.uint8_t, ndim=2] mask,
np.ndarray[np.uint8_t, ndim=2] out,
char shift_x, char shift_y, float p0, float p1):
""" Main loop, this function computes the histogram for each image point
- data is uint8
- result is uint8 casted
"""
cdef int rows = image.shape[0]
cdef int cols = image.shape[1]
cdef int srows = selem.shape[0]
cdef int scols = selem.shape[1]
cdef int centre_r = int(selem.shape[0] / 2) + shift_y
cdef int 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
image = np.ascontiguousarray(image)
if mask is None:
mask = np.ones((rows, cols), dtype=np.uint8)
else:
mask = np.ascontiguousarray(mask)
if out is None:
out = np.zeros((rows, cols), dtype=np.uint8)
else:
out = np.ascontiguousarray(out)
# create extended image and mask
cdef int erows = rows+srows-1
cdef int ecols = cols+scols-1
cdef np.ndarray emask = np.zeros((erows, ecols), dtype=np.uint8)
cdef np.ndarray eimage = np.zeros((erows, ecols), dtype=np.uint8)
eimage[centre_r:rows+centre_r,centre_c:cols+centre_c] = image
emask[centre_r:rows+centre_r,centre_c:cols+centre_c] = mask
mask = np.ascontiguousarray(mask)
# define pointers to the data
cdef np.uint8_t* eimage_data = <np.uint8_t*>eimage.data
cdef np.uint8_t* emask_data = <np.uint8_t*>emask.data
cdef np.uint8_t* out_data = <np.uint8_t*>out.data
cdef np.uint8_t* image_data = <np.uint8_t*>image.data
cdef np.uint8_t* mask_data = <np.uint8_t*>mask.data
# define local variable types
cdef int r, c, rr, cc, s, value, local_max, i, even_row
cdef float pop # number of pixels actually inside the neighborhood (float)
# allocate memory with malloc
cdef int max_se = srows*scols
# number of element in each attack border
cdef int num_se_n, num_se_s, num_se_e, num_se_w
# the current local histogram distribution
cdef int* histo = <int*>malloc(256 * sizeof(int))
# 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 int* se_e_r = <int*>malloc(max_se * sizeof(int))
cdef int* se_e_c = <int*>malloc(max_se * sizeof(int))
cdef int* se_w_r = <int*>malloc(max_se * sizeof(int))
cdef int* se_w_c = <int*>malloc(max_se * sizeof(int))
cdef int* se_n_r = <int*>malloc(max_se * sizeof(int))
cdef int* se_n_c = <int*>malloc(max_se * sizeof(int))
cdef int* se_s_r = <int*>malloc(max_se * sizeof(int))
cdef int* se_s_c = <int*>malloc(max_se * sizeof(int))
# build attack and release borders
# by using difference along axis
t = np.hstack((selem,np.zeros((selem.shape[0],1))))
t_e = np.diff(t,axis=1)==-1
t = np.hstack((np.zeros((selem.shape[0],1)),selem))
t_w = np.diff(t,axis=1)==1
t = np.vstack((selem,np.zeros((1,selem.shape[1]))))
t_s = np.diff(t,axis=0)==-1
t = np.vstack((np.zeros((1,selem.shape[1])),selem))
t_n = np.diff(t,axis=0)==1
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(256):
histo[i] = 0
pop = 0
for r in range(srows):
for c in range(scols):
rr = r
cc = c
if selem[r, c]:
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
r = 0
c = 0
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],p0,p1)
# 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] + centre_r
cc = c + se_e_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
for s in range(num_se_w):
rr = r + se_w_r[s] + centre_r
cc = c + se_w_c[s] + centre_c - 1
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] -= 1
pop -= 1.
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],p0,p1)
# 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] + centre_r
cc = c + se_s_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
for s in range(num_se_n):
rr = r + se_n_r[s] + centre_r - 1
cc = c + se_n_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] -= 1
pop -= 1.
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],p0,p1)
# 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] + centre_r
cc = c + se_w_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
for s in range(num_se_e):
rr = r + se_e_r[s] + centre_r
cc = c + se_e_c[s] + centre_c + 1
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] -= 1
pop -= 1.
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],p0,p1)
# 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] + centre_r
cc = c + se_s_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
for s in range(num_se_n):
rr = r + se_n_r[s] + centre_r - 1
cc = c + se_n_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] -= 1
pop -= 1.
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],p0,p1)
# 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)
return out
char shift_x, char shift_y, float p0, float p1)
skimage/rank/_core8p.pyx
+266
View File
@@ -0,0 +1,266 @@
#cython: cdivision=True
#cython: boundscheck=False
#cython: nonecheck=False
#cython: wraparound=False
import numpy as np
cimport numpy as np
from libc.stdlib cimport malloc, free
# generic cdef functions
cdef inline np.uint8_t uint8_max(np.uint8_t a, np.uint8_t b): return a if a >= b else b
cdef inline np.uint8_t uint8_min(np.uint8_t a, np.uint8_t b): return a if a <= b else b
#---------------------------------------------------------------------------
# 8 bit core kernel receives extra information about data inferior and superior percentiles
#---------------------------------------------------------------------------
cdef inline _core8p(np.uint8_t kernel(int*, float, np.uint8_t, float, float),
np.ndarray[np.uint8_t, ndim=2] image,
np.ndarray[np.uint8_t, ndim=2] selem,
np.ndarray[np.uint8_t, ndim=2] mask,
np.ndarray[np.uint8_t, ndim=2] out,
char shift_x, char shift_y, float p0, float p1):
""" Main loop, this function computes the histogram for each image point
- data is uint8
- result is uint8 casted
"""
cdef int rows = image.shape[0]
cdef int cols = image.shape[1]
cdef int srows = selem.shape[0]
cdef int scols = selem.shape[1]
cdef int centre_r = int(selem.shape[0] / 2) + shift_y
cdef int 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
image = np.ascontiguousarray(image)
if mask is None:
mask = np.ones((rows, cols), dtype=np.uint8)
else:
mask = np.ascontiguousarray(mask)
if out is None:
out = np.zeros((rows, cols), dtype=np.uint8)
else:
out = np.ascontiguousarray(out)
# create extended image and mask
cdef int erows = rows+srows-1
cdef int ecols = cols+scols-1
cdef np.ndarray emask = np.zeros((erows, ecols), dtype=np.uint8)
cdef np.ndarray eimage = np.zeros((erows, ecols), dtype=np.uint8)
eimage[centre_r:rows+centre_r,centre_c:cols+centre_c] = image
emask[centre_r:rows+centre_r,centre_c:cols+centre_c] = mask
mask = np.ascontiguousarray(mask)
# define pointers to the data
cdef np.uint8_t* eimage_data = <np.uint8_t*>eimage.data
cdef np.uint8_t* emask_data = <np.uint8_t*>emask.data
cdef np.uint8_t* out_data = <np.uint8_t*>out.data
cdef np.uint8_t* image_data = <np.uint8_t*>image.data
cdef np.uint8_t* mask_data = <np.uint8_t*>mask.data
# define local variable types
cdef int r, c, rr, cc, s, value, local_max, i, even_row
cdef float pop # number of pixels actually inside the neighborhood (float)
# allocate memory with malloc
cdef int max_se = srows*scols
# number of element in each attack border
cdef int num_se_n, num_se_s, num_se_e, num_se_w
# the current local histogram distribution
cdef int* histo = <int*>malloc(256 * sizeof(int))
# 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 int* se_e_r = <int*>malloc(max_se * sizeof(int))
cdef int* se_e_c = <int*>malloc(max_se * sizeof(int))
cdef int* se_w_r = <int*>malloc(max_se * sizeof(int))
cdef int* se_w_c = <int*>malloc(max_se * sizeof(int))
cdef int* se_n_r = <int*>malloc(max_se * sizeof(int))
cdef int* se_n_c = <int*>malloc(max_se * sizeof(int))
cdef int* se_s_r = <int*>malloc(max_se * sizeof(int))
cdef int* se_s_c = <int*>malloc(max_se * sizeof(int))
# build attack and release borders
# by using difference along axis
t = np.hstack((selem,np.zeros((selem.shape[0],1))))
t_e = np.diff(t,axis=1)==-1
t = np.hstack((np.zeros((selem.shape[0],1)),selem))
t_w = np.diff(t,axis=1)==1
t = np.vstack((selem,np.zeros((1,selem.shape[1]))))
t_s = np.diff(t,axis=0)==-1
t = np.vstack((np.zeros((1,selem.shape[1])),selem))
t_n = np.diff(t,axis=0)==1
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(256):
histo[i] = 0
pop = 0
for r in range(srows):
for c in range(scols):
rr = r
cc = c
if selem[r, c]:
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
r = 0
c = 0
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],p0,p1)
# 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] + centre_r
cc = c + se_e_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
for s in range(num_se_w):
rr = r + se_w_r[s] + centre_r
cc = c + se_w_c[s] + centre_c - 1
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] -= 1
pop -= 1.
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],p0,p1)
# 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] + centre_r
cc = c + se_s_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
for s in range(num_se_n):
rr = r + se_n_r[s] + centre_r - 1
cc = c + se_n_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] -= 1
pop -= 1.
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],p0,p1)
# 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] + centre_r
cc = c + se_w_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
for s in range(num_se_e):
rr = r + se_e_r[s] + centre_r
cc = c + se_e_c[s] + centre_c + 1
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] -= 1
pop -= 1.
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],p0,p1)
# 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] + centre_r
cc = c + se_s_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] += 1
pop += 1.
for s in range(num_se_n):
rr = r + se_n_r[s] + centre_r - 1
cc = c + se_n_c[s] + centre_c
if emask_data[rr * ecols + cc]:
value = eimage_data[rr * ecols + cc]
histo[value] -= 1
pop -= 1.
# kernel -------------------------------------------
out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],p0,p1)
# 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)
return out
skimage/rank/setup.py
+18
View File
@@ -5,6 +5,8 @@ from skimage._build import cython
base_path = os.path.abspath(os.path.dirname(__file__))
import sys
sys.path.append('.')
def configuration(parent_package='', top_path=None):
from numpy.distutils.misc_util import Configuration, get_numpy_include_dirs
@@ -12,12 +14,28 @@ def configuration(parent_package='', top_path=None):
config = Configuration('rank', parent_package, top_path)
# config.add_data_dir('tests')
cython(['_core8.pyx'], working_path=base_path)
cython(['_core8p.pyx'], working_path=base_path)
cython(['_core16.pyx'], working_path=base_path)
cython(['_core16p.pyx'], working_path=base_path)
cython(['_core16b.pyx'], working_path=base_path)
cython(['_crank8.pyx'], working_path=base_path)
cython(['_crank8_percentiles.pyx'], working_path=base_path)
cython(['_crank16.pyx'], working_path=base_path)
cython(['_crank16_percentiles.pyx'], working_path=base_path)
cython(['_crank16_bilateral.pyx'], working_path=base_path)
config.add_extension('_core8', sources=['_core8.c'],
include_dirs=[get_numpy_include_dirs()])
config.add_extension('_core8p', sources=['_core8p.c'],
include_dirs=[get_numpy_include_dirs()])
config.add_extension('_core16', sources=['_core16.c'],
include_dirs=[get_numpy_include_dirs()])
config.add_extension('_core16p', sources=['_core16p.c'],
include_dirs=[get_numpy_include_dirs()])
config.add_extension('_core16b', sources=['_core16b.c'],
include_dirs=[get_numpy_include_dirs()])
config.add_extension('_crank8', sources=['_crank8.c'],
include_dirs=[get_numpy_include_dirs()])
config.add_extension('_crank8_percentiles', sources=['_crank8_percentiles.c'],