wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-06-29 17:05:09 +08:00
Code Issues Packages Projects Releases Wiki Activity

adjust code for new project

Browse Source
This commit is contained in:
odebeir
2012-10-04 09:13:33 +02:00
parent 26f9186bc9
commit 413db93d09
13 changed files with 1832 additions and 108 deletions
Download Patch File Download Diff File Expand all files Collapse all files
skimage/rank/app.log
+1
View File
@@ -0,0 +1 @@
2012-10-04 09:12:19,785 MainProcess 5265 start logging in app.log
skimage/rank/app.py
+92 -81
View File
@@ -1,3 +1,5 @@
import unittest
import numpy as np
from time import time
import matplotlib.pyplot as plt
@@ -9,7 +11,7 @@ import crank
import crank16
import crank_percentiles
import crank16_percentiles
from pyrankfilter import filter
import crank16_bilateral
from cmorph import dilate
@@ -17,10 +19,6 @@ from cmorph import dilate
def c_max(image,selem):
return crank.maximum(image=image,selem = selem)
@log_timing
def w_max(image,selem):
return filter.maximum(image,struct_elem = selem)
@log_timing
def cm_max(image,selem):
return dilate(image=image,selem = selem)
@@ -39,10 +37,8 @@ def compare():
elem = np.ones((r,r),dtype='uint8')
# elem = (np.random.random((r,r))>.5).astype('uint8')
(rc,ms_rc) = c_max(a,elem)
(rw, ms_rw) = w_max(a,elem)
(rcm,ms_rcm) = cm_max(a,elem)
rec.append((ms_rc,ms_rw,ms_rcm))
assert (rc==rw).all()
assert (rc==rcm).all()
rec = np.asarray(rec)
@@ -59,10 +55,8 @@ def compare():
for s in range(100,1000,100):
a = (np.random.random((s,s))*256).astype('uint8')
(rc,ms_rc) = c_max(a,elem)
(rw, ms_rw) = w_max(a,elem)
(rcm,ms_rcm) = cm_max(a,elem)
rec.append((ms_rc,ms_rw,ms_rcm))
assert (rc==rw).all()
assert (rc==rcm).all()
rec = np.asarray(rec)
@@ -71,89 +65,105 @@ def compare():
plt.plot(rec)
plt.legend(['sliding cython','sliding weaves','cmorph'])
plt.figure()
plt.imshow(np.hstack((rc,rw,rcm)))
plt.imshow(np.hstack((rc,rcm)))
plt.show()
class TestSequenceFunctions(unittest.TestCase):
def test_image_size():
"""try several image sizes to check bounds conditions
"""
niter = 10
elem = np.asarray([[1,1,1],[1,1,1],[1,1,1]],dtype='uint8')
for m,n in np.random.random_integers(1,100,size=(10,2)):
a = np.ones((m,n),dtype='uint8')
r = crank.mean(image=a,selem = elem,shift_x=0,shift_y=0)
assert a.shape == r.shape
r = crank.mean(image=a,selem = elem,shift_x=0,shift_y=-1)
assert a.shape == r.shape
r = crank.mean(image=a,selem = elem,shift_x=0,shift_y=+1)
assert a.shape == r.shape
r = crank.mean(image=a,selem = elem,shift_x=-1,shift_y=0)
assert a.shape == r.shape
r = crank.mean(image=a,selem = elem,shift_x=+1,shift_y=0)
assert a.shape == r.shape
r = crank.mean(image=a,selem = elem,shift_x=-1,shift_y=-1)
assert a.shape == r.shape
r = crank.mean(image=a,selem = elem,shift_x=+1,shift_y=+1)
assert a.shape == r.shape
def setUp(self):
pass
return True
def test_random_sizes(self):
# make sure the size is not a problem
niter = 10
elem = np.asarray([[1,1,1],[1,1,1],[1,1,1]],dtype='uint8')
for m,n in np.random.random_integers(1,100,size=(10,2)):
a8 = np.ones((m,n),dtype='uint8')
r = crank.mean(image=a8,selem = elem,shift_x=0,shift_y=0)
self.assertTrue(a8.shape == r.shape)
r = crank.mean(image=a8,selem = elem,shift_x=+1,shift_y=+1)
self.assertTrue(a8.shape == r.shape)
for m,n in np.random.random_integers(1,100,size=(10,2)):
a16 = np.ones((m,n),dtype='uint16')
r = crank16.mean(image=a16,selem = elem,shift_x=0,shift_y=0)
self.assertTrue(a16.shape == r.shape)
r = crank16.mean(image=a16,selem = elem,shift_x=+1,shift_y=+1)
self.assertTrue(a16.shape == r.shape)
for m,n in np.random.random_integers(1,100,size=(10,2)):
a16 = np.ones((m,n),dtype='uint16')
r = crank16_percentiles.mean(image=a16,selem = elem,shift_x=0,shift_y=0,p0=.1,p1=.9)
self.assertTrue(a16.shape == r.shape)
r = crank16_percentiles.mean(image=a16,selem = elem,shift_x=+1,shift_y=+1,p0=.1,p1=.9)
self.assertTrue(a16.shape == r.shape)
def test_compare_with_cmorph(self):
#compare the result of maximum filter with dilate
a = (np.random.random((500,500))*256).astype('uint8')
for r in range(1,20,1):
elem = np.ones((r,r),dtype='uint8')
# elem = (np.random.random((r,r))>.5).astype('uint8')
rc = crank.maximum(image=a,selem = elem)
cm = dilate(image=a,selem = elem)
self.assertTrue((rc==cm).all())
def test_bitdepth(self):
elem = np.ones((3,3),dtype='uint8')
a16 = np.ones((100,100),dtype='uint16')*255
r = crank16_percentiles.mean(image=a16,selem = elem,shift_x=0,shift_y=0,p0=.1,p1=.9,bitdepth=8)
a16 = np.ones((100,100),dtype='uint16')*255*2
r = crank16_percentiles.mean(image=a16,selem = elem,shift_x=0,shift_y=0,p0=.1,p1=.9,bitdepth=9)
a16 = np.ones((100,100),dtype='uint16')*255*4
r = crank16_percentiles.mean(image=a16,selem = elem,shift_x=0,shift_y=0,p0=.1,p1=.9,bitdepth=10)
a16 = np.ones((100,100),dtype='uint16')*255*8
r = crank16_percentiles.mean(image=a16,selem = elem,shift_x=0,shift_y=0,p0=.1,p1=.9,bitdepth=11)
a16 = np.ones((100,100),dtype='uint16')*255*16
r = crank16_percentiles.mean(image=a16,selem = elem,shift_x=0,shift_y=0,p0=.1,p1=.9,bitdepth=12)
def test_population(self):
a = np.zeros((5,5),dtype='uint8')
elem = np.ones((3,3),dtype='uint8')
p = crank.pop(image=a,selem = elem)
r = np.asarray([[4, 6, 6, 6, 4],
[6, 9, 9, 9, 6],
[6, 9, 9, 9, 6],
[6, 9, 9, 9, 6],
[4, 6, 6, 6, 4]])
np.testing.assert_array_equal(r,p)
def test_structuring_element(self):
a = np.zeros((6,6),dtype='uint8')
a[2,2] = 255
elem = np.asarray([[1,1,0],[1,1,1],[0,0,1]],dtype='uint8')
f = crank.maximum(image=a,selem = elem,shift_x=1,shift_y=1)
r = np.asarray([[ 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0],
[ 0, 0, 255, 0, 0, 0],
[ 0, 0, 255, 255, 255, 0],
[ 0, 0, 0, 255, 255, 0],
[ 0, 0, 0, 0, 0, 0]])
np.testing.assert_array_equal(r,f)
@unittest.expectedFailure
def test_fail_on_bitdepth(self):
# should fail because data bitdepth is too high for the function
a16 = np.ones((100,100),dtype='uint16')*255
elem = np.ones((3,3),dtype='uint8')
f = crank16_percentiles.mean(image=a16,selem = elem,shift_x=0,shift_y=0,p0=.1,p1=.9,bitdepth=4)
if __name__ == '__main__':
logger = init_logger('app.log')
a = np.zeros((10,10),dtype='uint8')
a[2,2] = 255
# a[2,3] = 255
# a[2,4] = 255
print a
# unittest.main()
suite = unittest.TestLoader().loadTestsFromTestCase(TestSequenceFunctions)
unittest.TextTestRunner(verbosity=2).run(suite)
mask = np.ones_like(a)
# mask[:3,:3] = 0
# elem = np.asarray([[0,1,0],[1,1,1],[0,1,0]],dtype='uint8')
elem = np.asarray([[1,1,0],[1,1,1],[0,0,1]],dtype='uint8')
niter = 1
t0 = time()
for iter in range(niter):
r = crank.mean(image=a,selem = elem,shift_x=0,shift_y=0,mask = mask)
p = crank.pop(image=a,selem = elem,shift_x=0,shift_y=0,mask = mask)
t1 = time()
print '%f msec'%(t1-t0)
print 'cython mean'
print r
print p
t0 = time()
for iter in range(niter):
r = filter.mean(a,struct_elem = elem,struct_elem_center=(1,1),mask = mask)
t1 = time()
print '%f msec'%(t1-t0)
print 'filter.mean:'
print r
print a
r = crank.maximum(image=a,selem = elem,shift_x=0,shift_y=0,mask = mask)
print r
r = crank.gradient(image=r,selem = elem,shift_x=0,shift_y=0,mask = mask)
print r
im = np.zeros((10,10),dtype='uint8')
im[2:6,2:6] = 255
elem = np.asarray([[1,1,1],[1,1,1],[1,1,1]],dtype='uint8')
f = crank.gradient(image=im,selem = elem)
print f
f = crank.egalise(image=im,selem = elem)
print f
# compare()
# test_image_size()
# a = (data.coins()).astype('uint8')
a8 = (data.coins()).astype('uint8')
@@ -162,9 +172,10 @@ if __name__ == '__main__':
# f1 = filter.soft_gradient(a,struct_elem = selem,bitDepth=8,infSup=[.1,.9])
# f2 = crank16.bottomhat(a,selem = selem,bitdepth=12)
f1 = crank_percentiles.mean(a8,selem = selem,p0=.1,p1=.9)
f2 = crank16_percentiles.mean(a,selem = selem,bitdepth=12,p0=.1,p1=.9)
# f2 = crank16_percentiles.mean(a,selem = selem,bitdepth=12,p0=.1,p1=.9)
f2 = crank16_bilateral.mean(a,selem = selem,bitdepth=12,s0=500,s1=500)
# plt.imshow(f2)
plt.imshow(np.hstack((f1,f2)))
plt.imshow(np.hstack((a,f2)))
plt.colorbar()
plt.show()
skimage/rank/core16.pxd
+285
View File
@@ -0,0 +1,285 @@
""" 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
# 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 bitdepth
#---------------------------------------------------------------------------
cdef inline rank16(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
eimage = np.ascontiguousarray(eimage)
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
cdef int n_se_n, n_se_s, n_se_e, n_se_w
cdef int selem_num = np.sum(selem != 0)
cdef int* sr = <int*>malloc(selem_num * sizeof(int))
cdef int* sc = <int*>malloc(selem_num * sizeof(int))
cdef int* histo = <int*>malloc(maxbin * sizeof(int))
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
n_se_n = n_se_s = n_se_e = n_se_w = 0
for r in range(srows):
for c in range(scols):
if t_e[r,c]:
se_e_r[n_se_e] = r - centre_r
se_e_c[n_se_e] = c - centre_c
n_se_e += 1
if t_w[r,c]:
se_w_r[n_se_w] = r - centre_r
se_w_c[n_se_w] = c - centre_c
n_se_w += 1
if t_n[r,c]:
se_n_r[n_se_n] = r - centre_r
se_n_c[n_se_n] = c - centre_c
n_se_n += 1
if t_s[r,c]:
se_s_r[n_se_s] = r - centre_r
se_s_c[n_se_s] = c - centre_c
n_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(n_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(n_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(n_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(n_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(n_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(n_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(n_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(n_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(sr)
free(sc)
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
+286
View File
@@ -0,0 +1,286 @@
""" 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
# 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 bitdepth and bilateral interval
#---------------------------------------------------------------------------
cdef inline rank16b(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
eimage = np.ascontiguousarray(eimage)
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
cdef int n_se_n, n_se_s, n_se_e, n_se_w
cdef int selem_num = np.sum(selem != 0)
cdef int* sr = <int*>malloc(selem_num * sizeof(int))
cdef int* sc = <int*>malloc(selem_num * sizeof(int))
cdef int* histo = <int*>malloc(maxbin * sizeof(int))
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
n_se_n = n_se_s = n_se_e = n_se_w = 0
for r in range(srows):
for c in range(scols):
if t_e[r,c]:
se_e_r[n_se_e] = r - centre_r
se_e_c[n_se_e] = c - centre_c
n_se_e += 1
if t_w[r,c]:
se_w_r[n_se_w] = r - centre_r
se_w_c[n_se_w] = c - centre_c
n_se_w += 1
if t_n[r,c]:
se_n_r[n_se_n] = r - centre_r
se_n_c[n_se_n] = c - centre_c
n_se_n += 1
if t_s[r,c]:
se_s_r[n_se_s] = r - centre_r
se_s_c[n_se_s] = c - centre_c
n_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(n_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(n_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(n_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(n_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(n_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(n_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(n_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(n_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(sr)
free(sc)
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
+287
View File
@@ -0,0 +1,287 @@
""" to compile this use:
>>> python setup.py build_ext --inplace
to generate html report use:
>>> cython -a core16p.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
# 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 rank16_percentile(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
eimage = np.ascontiguousarray(eimage)
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
cdef int n_se_n, n_se_s, n_se_e, n_se_w
cdef int selem_num = np.sum(selem != 0)
cdef int* sr = <int*>malloc(selem_num * sizeof(int))
cdef int* sc = <int*>malloc(selem_num * sizeof(int))
cdef int* histo = <int*>malloc(maxbin * sizeof(int))
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
n_se_n = n_se_s = n_se_e = n_se_w = 0
for r in range(srows):
for c in range(scols):
if t_e[r,c]:
se_e_r[n_se_e] = r - centre_r
se_e_c[n_se_e] = c - centre_c
n_se_e += 1
if t_w[r,c]:
se_w_r[n_se_w] = r - centre_r
se_w_c[n_se_w] = c - centre_c
n_se_w += 1
if t_n[r,c]:
se_n_r[n_se_n] = r - centre_r
se_n_c[n_se_n] = c - centre_c
n_se_n += 1
if t_s[r,c]:
se_s_r[n_se_s] = r - centre_r
se_s_c[n_se_s] = c - centre_c
n_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(n_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(n_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(n_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(n_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(n_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(n_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(n_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(n_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(sr)
free(sc)
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
+271
View File
@@ -0,0 +1,271 @@
""" 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
# 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
#---------------------------------------------------------------------------
# 8 bit core kernel
#---------------------------------------------------------------------------
cdef inline rank8(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
eimage = np.ascontiguousarray(eimage)
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
cdef int n_se_n, n_se_s, n_se_e, n_se_w
cdef int selem_num = np.sum(selem != 0)
cdef int* sr = <int*>malloc(selem_num * sizeof(int))
cdef int* sc = <int*>malloc(selem_num * sizeof(int))
cdef int* histo = <int*>malloc(256 * sizeof(int))
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
n_se_n = n_se_s = n_se_e = n_se_w = 0
for r in range(srows):
for c in range(scols):
if t_e[r,c]:
se_e_r[n_se_e] = r - centre_r
se_e_c[n_se_e] = c - centre_c
n_se_e += 1
if t_w[r,c]:
se_w_r[n_se_w] = r - centre_r
se_w_c[n_se_w] = c - centre_c
n_se_w += 1
if t_n[r,c]:
se_n_r[n_se_n] = r - centre_r
se_n_c[n_se_n] = c - centre_c
n_se_n += 1
if t_s[r,c]:
se_s_r[n_se_s] = r - centre_r
se_s_c[n_se_s] = c - centre_c
n_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(n_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(n_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(n_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(n_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(n_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(n_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(n_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(n_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(sr)
free(sc)
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
+271
View File
@@ -0,0 +1,271 @@
""" to compile this use:
>>> python setup.py build_ext --inplace
to generate html report use:
>>> cython -a core8p.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
# 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
#---------------------------------------------------------------------------
# 8 bit core kernel receives extra information about data inferior and superior percentiles
#---------------------------------------------------------------------------
cdef inline rank8_percentile(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
eimage = np.ascontiguousarray(eimage)
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
cdef int n_se_n, n_se_s, n_se_e, n_se_w
cdef int selem_num = np.sum(selem != 0)
cdef int* sr = <int*>malloc(selem_num * sizeof(int))
cdef int* sc = <int*>malloc(selem_num * sizeof(int))
cdef int* histo = <int*>malloc(256 * sizeof(int))
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
n_se_n = n_se_s = n_se_e = n_se_w = 0
for r in range(srows):
for c in range(scols):
if t_e[r,c]:
se_e_r[n_se_e] = r - centre_r
se_e_c[n_se_e] = c - centre_c
n_se_e += 1
if t_w[r,c]:
se_w_r[n_se_w] = r - centre_r
se_w_c[n_se_w] = c - centre_c
n_se_w += 1
if t_n[r,c]:
se_n_r[n_se_n] = r - centre_r
se_n_c[n_se_n] = c - centre_c
n_se_n += 1
if t_s[r,c]:
se_s_r[n_se_s] = r - centre_r
se_s_c[n_se_s] = c - centre_c
n_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(n_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(n_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(n_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(n_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(n_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(n_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(n_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(n_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(sr)
free(sc)
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/crank.pyx
+1 -6
View File
@@ -15,12 +15,7 @@ import numpy as np
cimport numpy as np
# import main loop
from core cimport rank8
# todo
# - manage float output,
# - manage different bit depth input
# - add auxiliary parameters (spectral_interval, infSup)
from core8 cimport rank8
# -----------------------------------------------------------------
# kernels uint8
skimage/rank/crank16.pyx
+2 -7
View File
@@ -2,7 +2,7 @@
>>> python setup.py build_ext --inplace
to generate html report use:
>>> cython -a crank.pxd
>>> cython -a crank16.pxd
"""
@@ -15,12 +15,7 @@ import numpy as np
cimport numpy as np
# import main loop
from core cimport rank16
# todo
# - manage float output,
# - manage different bit depth input
# - add auxiliary parameters (spectral_interval, infSup)
from core16 cimport rank16
# -----------------------------------------------------------------
# kernels uint16 take extra parameter for defining the bitdepth
skimage/rank/crank16_bilateral.pyx
+332
View File
@@ -0,0 +1,332 @@
""" to compile this use:
>>> python setup.py build_ext --inplace
to generate html report use:
>>> cython -a crank16.pxd
"""
#cython: cdivision=True
#cython: boundscheck=False
#cython: nonecheck=False
#cython: wraparound=False
import numpy as np
cimport numpy as np
# import main loop
from core16b cimport rank16b
# -----------------------------------------------------------------
# kernels uint16 take extra parameter for defining the bitdepth
# -----------------------------------------------------------------
#cdef inline np.uint16_t kernel_autolevel(int* histo, float pop, np.uint16_t g,int bitdepth,int maxbin, int midbin):
# cdef int 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 <np.uint16_t>(maxbin*1.*(g-imin)/delta)
# else:
# return <np.uint16_t>(imax-imin)
#
#cdef inline np.uint16_t kernel_bottomhat(int* histo, float pop, np.uint16_t g,int bitdepth,int maxbin, int midbin):
# cdef int i
#
# for i in range(maxbin):
# if histo[i]:
# break
#
# return <np.uint16_t>(g-i)
#
#
#cdef inline np.uint16_t kernel_egalise(int* histo, float pop, np.uint16_t g,int bitdepth,int maxbin, int midbin):
# cdef int i
# cdef float sum = 0.
#
# if pop:
# for i in range(maxbin):
# sum += histo[i]
# if i>=g:
# break
#
# return <np.uint16_t>((maxbin*1.*sum)/pop)
# else:
# return <np.uint16_t>(0)
#
#cdef inline np.uint16_t kernel_gradient(int* histo, float pop, np.uint16_t g,int bitdepth,int maxbin, int midbin):
# cdef int 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 <np.uint16_t>(imax-imin)
# else:
# return <np.uint16_t>(0)
#
#cdef inline np.uint16_t kernel_maximum(int* histo, float pop, np.uint16_t g,int bitdepth,int maxbin, int midbin):
# cdef int i
#
# if pop:
# for i in range(maxbin-1,-1,-1):
# if histo[i]:
# return <np.uint16_t>(i)
#
# return <np.uint16_t>(0)
cdef inline np.uint16_t kernel_mean(int* histo, float pop, np.uint16_t g,int bitdepth,int maxbin, int midbin, int s0, int 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 <np.uint16_t>(mean/bilat_pop)
else:
return <np.uint16_t>(0)
else:
return <np.uint16_t>(0)
#cdef inline np.uint16_t kernel_meansubstraction(int* histo, float pop, np.uint16_t g,int bitdepth,int maxbin, int midbin):
# cdef int i
# cdef float mean = 0.
#
# if pop:
# for i in range(maxbin):
# mean += histo[i]*i
# return <np.uint16_t>((g-mean/pop)/2.+midbin)
# else:
# return <np.uint16_t>(0)
#
#cdef inline np.uint16_t kernel_median(int* histo, float pop, np.uint16_t g,int bitdepth,int maxbin, int midbin):
# cdef int i
# cdef float sum = pop/2.0
#
# if pop:
# for i in range(maxbin):
# if histo[i]:
# sum -= histo[i]
# if sum<0:
# return <np.uint16_t>(i)
#
# return <np.uint16_t>(0)
#
#cdef inline np.uint16_t kernel_minimum(int* histo, float pop, np.uint16_t g,int bitdepth,int maxbin, int midbin):
# cdef int i
#
# if pop:
# for i in range(maxbin):
# if histo[i]:
# return <np.uint16_t>(i)
#
# return <np.uint16_t>(0)
#
#cdef inline np.uint16_t kernel_modal(int* histo, float pop, np.uint16_t g,int bitdepth,int maxbin, int midbin):
# cdef int hmax=0,imax=0
#
# if pop:
# for i in range(maxbin):
# if histo[i]>hmax:
# hmax = histo[i]
# imax = i
# return <np.uint16_t>(imax)
#
# return <np.uint16_t>(0)
#
#cdef inline np.uint16_t kernel_morph_contr_enh(int* histo, float pop, np.uint16_t g,int bitdepth,int maxbin, int midbin):
# cdef int 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 <np.uint16_t>(imax)
# else:
# return <np.uint16_t>(imin)
# else:
# return <np.uint16_t>(0)
#
cdef inline np.uint16_t kernel_pop(int* histo, float pop, np.uint16_t g,int bitdepth,int maxbin, int midbin, int s0, int 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 <np.uint16_t>(bilat_pop)
else:
return <np.uint16_t>(0)
#
#cdef inline np.uint16_t kernel_threshold(int* histo, float pop, np.uint16_t g,int bitdepth,int maxbin, int midbin):
# cdef int i
# cdef float mean = 0.
#
# if pop:
# for i in range(maxbin):
# mean += histo[i]*i
# return <np.uint16_t>(g>(mean/pop))
# else:
# return <np.uint16_t>(0)
#
#cdef inline np.uint16_t kernel_tophat(int* histo, float pop, np.uint16_t g,int bitdepth,int maxbin, int midbin):
# cdef int i
#
# for i in range(maxbin-1,-1,-1):
# if histo[i]:
# break
#
# return <np.uint16_t>(i-g)
# -----------------------------------------------------------------
# python wrappers
# -----------------------------------------------------------------
#def autolevel(np.ndarray[np.uint16_t, ndim=2] image,
# np.ndarray[np.uint8_t, ndim=2] selem,
# np.ndarray[np.uint8_t, ndim=2] mask=None,
# np.ndarray[np.uint16_t, ndim=2] out=None,
# char shift_x=0, char shift_y=0, int bitdepth=8, int s0=1, int s1=1):
# """bottom hat
# """
# return rank16b(kernel_autolevel,image,selem,mask,out,shift_x,shift_y,bitdepth,s0,s1)
#
#def bottomhat(np.ndarray[np.uint16_t, ndim=2] image,
# np.ndarray[np.uint8_t, ndim=2] selem,
# np.ndarray[np.uint8_t, ndim=2] mask=None,
# np.ndarray[np.uint16_t, ndim=2] out=None,
# char shift_x=0, char shift_y=0, int bitdepth=8, int s0=1, int s1=1):
# """bottom hat
# """
# return rank16b(kernel_bottomhat,image,selem,mask,out,shift_x,shift_y,bitdepth,s0,s1)
#
#def egalise(np.ndarray[np.uint16_t, ndim=2] image,
# np.ndarray[np.uint8_t, ndim=2] selem,
# np.ndarray[np.uint8_t, ndim=2] mask=None,
# np.ndarray[np.uint16_t, ndim=2] out=None,
# char shift_x=0, char shift_y=0, int bitdepth=8, int s0=1, int s1=1):
# """local egalisation of the gray level
# """
# return rank16b(kernel_egalise,image,selem,mask,out,shift_x,shift_y,bitdepth,s0,s1)
#
#def gradient(np.ndarray[np.uint16_t, ndim=2] image,
# np.ndarray[np.uint8_t, ndim=2] selem,
# np.ndarray[np.uint8_t, ndim=2] mask=None,
# np.ndarray[np.uint16_t, ndim=2] out=None,
# char shift_x=0, char shift_y=0, int bitdepth=8, int s0=1, int s1=1):
# """local maximum - local minimum gray level
# """
# return rank16b(kernel_gradient,image,selem,mask,out,shift_x,shift_y,bitdepth,s0,s1)
#
#def maximum(np.ndarray[np.uint16_t, ndim=2] image,
# np.ndarray[np.uint8_t, ndim=2] selem,
# np.ndarray[np.uint8_t, ndim=2] mask=None,
# np.ndarray[np.uint16_t, ndim=2] out=None,
# char shift_x=0, char shift_y=0, int bitdepth=8, int s0=1, int s1=1):
# """local maximum gray level
# """
# return rank16b(kernel_maximum,image,selem,mask,out,shift_x,shift_y,bitdepth,s0,s1)
def mean(np.ndarray[np.uint16_t, ndim=2] image,
np.ndarray[np.uint8_t, ndim=2] selem,
np.ndarray[np.uint8_t, ndim=2] mask=None,
np.ndarray[np.uint16_t, ndim=2] out=None,
char shift_x=0, char shift_y=0, int bitdepth=8, int s0=1, int s1=1):
"""average gray level (clipped on uint8)
"""
return rank16b(kernel_mean,image,selem,mask,out,shift_x,shift_y,bitdepth,s0,s1)
#def meansubstraction(np.ndarray[np.uint16_t, ndim=2] image,
# np.ndarray[np.uint8_t, ndim=2] selem,
# np.ndarray[np.uint8_t, ndim=2] mask=None,
# np.ndarray[np.uint16_t, ndim=2] out=None,
# char shift_x=0, char shift_y=0, int bitdepth=8, int s0=1, int s1=1):
# """(g - average gray level)/2+midbin (clipped on uint8)
# """
# return rank16b(kernel_meansubstraction,image,selem,mask,out,shift_x,shift_y,bitdepth,s0,s1)
#
#def median(np.ndarray[np.uint16_t, ndim=2] image,
# np.ndarray[np.uint8_t, ndim=2] selem,
# np.ndarray[np.uint8_t, ndim=2] mask=None,
# np.ndarray[np.uint16_t, ndim=2] out=None,
# char shift_x=0, char shift_y=0, int bitdepth=8, int s0=1, int s1=1):
# """local median
# """
# return rank16b(kernel_median,image,selem,mask,out,shift_x,shift_y,bitdepth,s0,s1)
#
#def minimum(np.ndarray[np.uint16_t, ndim=2] image,
# np.ndarray[np.uint8_t, ndim=2] selem,
# np.ndarray[np.uint8_t, ndim=2] mask=None,
# np.ndarray[np.uint16_t, ndim=2] out=None,
# char shift_x=0, char shift_y=0, int bitdepth=8, int s0=1, int s1=1):
# """local minimum gray level
# """
# return rank16b(kernel_minimum,image,selem,mask,out,shift_x,shift_y,bitdepth,s0,s1)
#
#def morph_contr_enh(np.ndarray[np.uint16_t, ndim=2] image,
# np.ndarray[np.uint8_t, ndim=2] selem,
# np.ndarray[np.uint8_t, ndim=2] mask=None,
# np.ndarray[np.uint16_t, ndim=2] out=None,
# char shift_x=0, char shift_y=0, int bitdepth=8, int s0=1, int s1=1):
# """morphological contrast enhancement
# """
# return rank16b(kernel_morph_contr_enh,image,selem,mask,out,shift_x,shift_y,bitdepth,s0,s1)
#
#def modal(np.ndarray[np.uint16_t, ndim=2] image,
# np.ndarray[np.uint8_t, ndim=2] selem,
# np.ndarray[np.uint8_t, ndim=2] mask=None,
# np.ndarray[np.uint16_t, ndim=2] out=None,
# char shift_x=0, char shift_y=0, int bitdepth=8, int s0=1, int s1=1):
# """local mode
# """
# return rank16b(kernel_modal,image,selem,mask,out,shift_x,shift_y,bitdepth,s0,s1)
#
def pop(np.ndarray[np.uint16_t, ndim=2] image,
np.ndarray[np.uint8_t, ndim=2] selem,
np.ndarray[np.uint8_t, ndim=2] mask=None,
np.ndarray[np.uint16_t, ndim=2] 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
"""
return rank16b(kernel_pop,image,selem,mask,out,shift_x,shift_y,bitdepth,s0,s1)
#def threshold(np.ndarray[np.uint16_t, ndim=2] image,
# np.ndarray[np.uint8_t, ndim=2] selem,
# np.ndarray[np.uint8_t, ndim=2] mask=None,
# np.ndarray[np.uint16_t, ndim=2] out=None,
# char shift_x=0, char shift_y=0, int bitdepth=8, int s0=1, int s1=1):
# """returns maxbin-1 if gray level higher than local mean, 0 else
# """
# return rank16b(kernel_threshold,image,selem,mask,out,shift_x,shift_y,bitdepth,s0,s1)
#
#def tophat(np.ndarray[np.uint16_t, ndim=2] image,
# np.ndarray[np.uint8_t, ndim=2] selem,
# np.ndarray[np.uint8_t, ndim=2] mask=None,
# np.ndarray[np.uint16_t, ndim=2] out=None,
# char shift_x=0, char shift_y=0, int bitdepth=8, int s0=1, int s1=1):
# """top hat
# """
# return rank16b(kernel_tophat,image,selem,mask,out,shift_x,shift_y,bitdepth,s0,s1)
skimage/rank/crank16_percentiles.pyx
+1 -6
View File
@@ -15,12 +15,7 @@ import numpy as np
cimport numpy as np
# import main loop
from core cimport rank16_percentile
# todo
# - manage float output,
# - manage different bit depth input
# - add auxiliary parameters (spectral_interval, infSup)
from core16p cimport rank16_percentile
# -----------------------------------------------------------------
# kernels uint8 (SOFT version using percentiles)
skimage/rank/crank_percentiles.pyx
+1 -6
View File
@@ -15,12 +15,7 @@ import numpy as np
cimport numpy as np
# import main loop
from core cimport rank8_percentile
# todo
# - manage float output,
# - manage different bit depth input
# - add auxiliary parameters (spectral_interval, infSup)
from core8p cimport rank8_percentile
# -----------------------------------------------------------------
# kernels uint8 (SOFT version using percentiles)
skimage/rank/setup.py
+2 -2
View File
@@ -6,10 +6,10 @@ from Cython.Distutils import build_ext
setup(
cmdclass = {'build_ext': build_ext},
ext_modules = [Extension("helloworld", ["helloworld.pyx"]),
Extension("cmorph", ["cmorph.pyx"], include_dirs=[np.get_include()]),
ext_modules = [Extension("cmorph", ["cmorph.pyx"], include_dirs=[np.get_include()]),
Extension("crank", ["crank.pyx"], include_dirs=[np.get_include()]),
Extension("crank_percentiles", ["crank_percentiles.pyx"], include_dirs=[np.get_include()]),
Extension("crank16", ["crank16.pyx"], include_dirs=[np.get_include()]),
Extension("crank16_bilateral", ["crank16_bilateral.pyx"], include_dirs=[np.get_include()]),
Extension("crank16_percentiles", ["crank16_percentiles.pyx"], include_dirs=[np.get_include()])]
)