diff --git a/skimage/rank/app.log b/skimage/rank/app.log new file mode 100644 index 00000000..7eae20f5 --- /dev/null +++ b/skimage/rank/app.log @@ -0,0 +1 @@ +2012-10-04 09:12:19,785 MainProcess 5265 start logging in app.log diff --git a/skimage/rank/app.py b/skimage/rank/app.py index acfc33ec..0cf24408 100644 --- a/skimage/rank/app.py +++ b/skimage/rank/app.py @@ -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() diff --git a/skimage/rank/core16.pxd b/skimage/rank/core16.pxd new file mode 100644 index 00000000..9973025f --- /dev/null +++ b/skimage/rank/core16.pxd @@ -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 (imageeimage.data + cdef np.uint8_t* emask_data = emask.data + + cdef np.uint16_t* out_data = out.data + cdef np.uint16_t* image_data = image.data + cdef np.uint8_t* mask_data = 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 = malloc(selem_num * sizeof(int)) + cdef int* sc = malloc(selem_num * sizeof(int)) + cdef int* histo = malloc(maxbin * sizeof(int)) + cdef int* se_e_r = malloc(max_se * sizeof(int)) + cdef int* se_e_c = malloc(max_se * sizeof(int)) + cdef int* se_w_r = malloc(max_se * sizeof(int)) + cdef int* se_w_c = malloc(max_se * sizeof(int)) + cdef int* se_n_r = malloc(max_se * sizeof(int)) + cdef int* se_n_c = malloc(max_se * sizeof(int)) + cdef int* se_s_r = malloc(max_se * sizeof(int)) + cdef int* se_s_c = 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 diff --git a/skimage/rank/core16b.pxd b/skimage/rank/core16b.pxd new file mode 100644 index 00000000..38832c11 --- /dev/null +++ b/skimage/rank/core16b.pxd @@ -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 (imageeimage.data + cdef np.uint8_t* emask_data = emask.data + + cdef np.uint16_t* out_data = out.data + cdef np.uint16_t* image_data = image.data + cdef np.uint8_t* mask_data = 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 = malloc(selem_num * sizeof(int)) + cdef int* sc = malloc(selem_num * sizeof(int)) + cdef int* histo = malloc(maxbin * sizeof(int)) + cdef int* se_e_r = malloc(max_se * sizeof(int)) + cdef int* se_e_c = malloc(max_se * sizeof(int)) + cdef int* se_w_r = malloc(max_se * sizeof(int)) + cdef int* se_w_c = malloc(max_se * sizeof(int)) + cdef int* se_n_r = malloc(max_se * sizeof(int)) + cdef int* se_n_c = malloc(max_se * sizeof(int)) + cdef int* se_s_r = malloc(max_se * sizeof(int)) + cdef int* se_s_c = 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 diff --git a/skimage/rank/core16p.pxd b/skimage/rank/core16p.pxd new file mode 100644 index 00000000..dbf51c54 --- /dev/null +++ b/skimage/rank/core16p.pxd @@ -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 (imageeimage.data + cdef np.uint8_t* emask_data = emask.data + + cdef np.uint16_t* out_data = out.data + cdef np.uint16_t* image_data = image.data + cdef np.uint8_t* mask_data = 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 = malloc(selem_num * sizeof(int)) + cdef int* sc = malloc(selem_num * sizeof(int)) + cdef int* histo = malloc(maxbin * sizeof(int)) + cdef int* se_e_r = malloc(max_se * sizeof(int)) + cdef int* se_e_c = malloc(max_se * sizeof(int)) + cdef int* se_w_r = malloc(max_se * sizeof(int)) + cdef int* se_w_c = malloc(max_se * sizeof(int)) + cdef int* se_n_r = malloc(max_se * sizeof(int)) + cdef int* se_n_c = malloc(max_se * sizeof(int)) + cdef int* se_s_r = malloc(max_se * sizeof(int)) + cdef int* se_s_c = 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 + + diff --git a/skimage/rank/core8.pxd b/skimage/rank/core8.pxd new file mode 100644 index 00000000..0c901c8b --- /dev/null +++ b/skimage/rank/core8.pxd @@ -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 = eimage.data + cdef np.uint8_t* emask_data = emask.data + + cdef np.uint8_t* out_data = out.data + cdef np.uint8_t* image_data = image.data + cdef np.uint8_t* mask_data = 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 = malloc(selem_num * sizeof(int)) + cdef int* sc = malloc(selem_num * sizeof(int)) + cdef int* histo = malloc(256 * sizeof(int)) + cdef int* se_e_r = malloc(max_se * sizeof(int)) + cdef int* se_e_c = malloc(max_se * sizeof(int)) + cdef int* se_w_r = malloc(max_se * sizeof(int)) + cdef int* se_w_c = malloc(max_se * sizeof(int)) + cdef int* se_n_r = malloc(max_se * sizeof(int)) + cdef int* se_n_c = malloc(max_se * sizeof(int)) + cdef int* se_s_r = malloc(max_se * sizeof(int)) + cdef int* se_s_c = 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 + diff --git a/skimage/rank/core8p.pxd b/skimage/rank/core8p.pxd new file mode 100644 index 00000000..25e3ffcf --- /dev/null +++ b/skimage/rank/core8p.pxd @@ -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 = eimage.data + cdef np.uint8_t* emask_data = emask.data + + cdef np.uint8_t* out_data = out.data + cdef np.uint8_t* image_data = image.data + cdef np.uint8_t* mask_data = 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 = malloc(selem_num * sizeof(int)) + cdef int* sc = malloc(selem_num * sizeof(int)) + cdef int* histo = malloc(256 * sizeof(int)) + cdef int* se_e_r = malloc(max_se * sizeof(int)) + cdef int* se_e_c = malloc(max_se * sizeof(int)) + cdef int* se_w_r = malloc(max_se * sizeof(int)) + cdef int* se_w_c = malloc(max_se * sizeof(int)) + cdef int* se_n_r = malloc(max_se * sizeof(int)) + cdef int* se_n_c = malloc(max_se * sizeof(int)) + cdef int* se_s_r = malloc(max_se * sizeof(int)) + cdef int* se_s_c = 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 + diff --git a/skimage/rank/crank.pyx b/skimage/rank/crank.pyx index 09048917..59016eed 100644 --- a/skimage/rank/crank.pyx +++ b/skimage/rank/crank.pyx @@ -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 diff --git a/skimage/rank/crank16.pyx b/skimage/rank/crank16.pyx index 78e5a2af..cd0bacd4 100644 --- a/skimage/rank/crank16.pyx +++ b/skimage/rank/crank16.pyx @@ -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 diff --git a/skimage/rank/crank16_bilateral.pyx b/skimage/rank/crank16_bilateral.pyx new file mode 100644 index 00000000..313b83d6 --- /dev/null +++ b/skimage/rank/crank16_bilateral.pyx @@ -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 (maxbin*1.*(g-imin)/delta) +# else: +# return (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 (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 ((maxbin*1.*sum)/pop) +# else: +# return (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 (imax-imin) +# else: +# return (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 (i) +# +# return (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 (mean/bilat_pop) + else: + return (0) + else: + return (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 ((g-mean/pop)/2.+midbin) +# else: +# return (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 (i) +# +# return (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 (i) +# +# return (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 (imax) +# +# return (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 (imax) +# else: +# return (imin) +# else: +# return (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 (bilat_pop) + else: + return (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 (g>(mean/pop)) +# else: +# return (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 (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) diff --git a/skimage/rank/crank16_percentiles.pyx b/skimage/rank/crank16_percentiles.pyx index b967c61a..7756bc19 100644 --- a/skimage/rank/crank16_percentiles.pyx +++ b/skimage/rank/crank16_percentiles.pyx @@ -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) diff --git a/skimage/rank/crank_percentiles.pyx b/skimage/rank/crank_percentiles.pyx index 51a9e01f..d4d6312f 100644 --- a/skimage/rank/crank_percentiles.pyx +++ b/skimage/rank/crank_percentiles.pyx @@ -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) diff --git a/skimage/rank/setup.py b/skimage/rank/setup.py index 6d18e27f..e9af9446 100644 --- a/skimage/rank/setup.py +++ b/skimage/rank/setup.py @@ -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()])] ) \ No newline at end of file