diff --git a/skimage/rank/README.rst b/skimage/rank/README.rst new file mode 100644 index 00000000..b15001e6 --- /dev/null +++ b/skimage/rank/README.rst @@ -0,0 +1,5 @@ +To use this to build your Cython file use the commandline options: + +.. sourcecode:: text + + $ python setup.py build_ext --inplace \ No newline at end of file diff --git a/skimage/rank/app.py b/skimage/rank/app.py new file mode 100644 index 00000000..acfc33ec --- /dev/null +++ b/skimage/rank/app.py @@ -0,0 +1,171 @@ +import numpy as np +from time import time +import matplotlib.pyplot as plt +from skimage import data + +from tools import log_timing,init_logger + +import crank +import crank16 +import crank_percentiles +import crank16_percentiles +from pyrankfilter import filter +from cmorph import dilate + + +@log_timing +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) + +def compare(): + """comparison between + - Cython maximum rankfilter implementation + - weaves maximum rankfilter implementation + - cmorph.dilate cython implementation + on increasing structuring element size and increasing image size + """ + a = (np.random.random((500,500))*256).astype('uint8') + + rec = [] + for r in range(1,20,1): + 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) + + plt.plot(rec) + plt.legend(['sliding cython','sliding weaves','cmorph']) + plt.figure() + plt.imshow(np.hstack((rc,rw,rcm))) + + r = 9 + elem = np.ones((r,r),dtype='uint8') + + rec = [] + 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) + + plt.figure() + plt.plot(rec) + plt.legend(['sliding cython','sliding weaves','cmorph']) + plt.figure() + plt.imshow(np.hstack((rc,rw,rcm))) + + plt.show() + +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 + + return True + + +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 + + 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') + a = (data.coins()).astype('uint16')*16 + selem = np.ones((20,20),dtype='uint8') +# 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) +# plt.imshow(f2) + plt.imshow(np.hstack((f1,f2))) + plt.colorbar() + plt.show() + + diff --git a/skimage/rank/cmorph.pyx b/skimage/rank/cmorph.pyx new file mode 100644 index 00000000..9b8b3a27 --- /dev/null +++ b/skimage/rank/cmorph.pyx @@ -0,0 +1,118 @@ +#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 + + +def dilate(np.ndarray[np.uint8_t, ndim=2] image, + np.ndarray[np.uint8_t, ndim=2] selem, + np.ndarray[np.uint8_t, ndim=2] out=None, + char shift_x=0, char shift_y=0): + + 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 + + image = np.ascontiguousarray(image) + if out is None: + out = np.zeros((rows, cols), dtype=np.uint8) + else: + out = np.ascontiguousarray(out) + + cdef np.uint8_t* out_data = out.data + cdef np.uint8_t* image_data = image.data + + cdef int r, c, rr, cc, s, value, local_max + + cdef int selem_num = np.sum(selem != 0) + cdef int* sr = malloc(selem_num * sizeof(int)) + cdef int* sc = malloc(selem_num * sizeof(int)) + + s = 0 + for r in range(srows): + for c in range(scols): + if selem[r, c] != 0: + sr[s] = r - centre_r + sc[s] = c - centre_c + s += 1 + + for r in range(rows): + for c in range(cols): + local_max = 0 + for s in range(selem_num): + rr = r + sr[s] + cc = c + sc[s] + if 0 <= rr < rows and 0 <= cc < cols: + value = image_data[rr * cols + cc] + if value > local_max: + local_max = value + + out_data[r * cols + c] = local_max + + free(sr) + free(sc) + + return out + + +def erode(np.ndarray[np.uint8_t, ndim=2] image, + np.ndarray[np.uint8_t, ndim=2] selem, + np.ndarray[np.uint8_t, ndim=2] out=None, + char shift_x=0, char shift_y=0): + + 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 + + image = np.ascontiguousarray(image) + if out is None: + out = np.zeros((rows, cols), dtype=np.uint8) + else: + out = np.ascontiguousarray(out) + + cdef np.uint8_t* out_data = out.data + cdef np.uint8_t* image_data = image.data + + cdef int r, c, rr, cc, s, value, local_min + + cdef int selem_num = np.sum(selem != 0) + cdef int* sr = malloc(selem_num * sizeof(int)) + cdef int* sc = malloc(selem_num * sizeof(int)) + + s = 0 + for r in range(srows): + for c in range(scols): + if selem[r, c] != 0: + sr[s] = r - centre_r + sc[s] = c - centre_c + s += 1 + + for r in range(rows): + for c in range(cols): + local_min = 255 + for s in range(selem_num): + rr = r + sr[s] + cc = c + sc[s] + if 0 <= rr < rows and 0 <= cc < cols: + value = image_data[rr * cols + cc] + if value < local_min: + local_min = value + + out_data[r * cols + c] = local_min + + free(sr) + free(sc) + + return out diff --git a/skimage/rank/core.pxd b/skimage/rank/core.pxd new file mode 100644 index 00000000..72facde3 --- /dev/null +++ b/skimage/rank/core.pxd @@ -0,0 +1,1054 @@ +""" to compile this use: +>>> python setup.py build_ext --inplace + +to generate html report use: +>>> cython -a core.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 + +#--------------------------------------------------------------------------- +# 16 bit core, kernel receive 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 + +#--------------------------------------------------------------------------- +# 8 bit core kernel receive 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 + +#--------------------------------------------------------------------------- +# 16 bit core kernel receive 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 \ No newline at end of file diff --git a/skimage/rank/crank.pyx b/skimage/rank/crank.pyx new file mode 100644 index 00000000..09048917 --- /dev/null +++ b/skimage/rank/crank.pyx @@ -0,0 +1,325 @@ +""" to compile this use: +>>> python setup.py build_ext --inplace + +to generate html report use: +>>> cython -a crank.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 core cimport rank8 + +# todo +# - manage float output, +# - manage different bit depth input +# - add auxiliary parameters (spectral_interval, infSup) + +# ----------------------------------------------------------------- +# kernels uint8 +# ----------------------------------------------------------------- + +cdef inline np.uint8_t kernel_autolevel(int* histo, float pop, np.uint8_t g): + cdef int i,imin,imax,delta + + if pop: + for i in range(255,-1,-1): + if histo[i]: + imax = i + break + for i in range(256): + if histo[i]: + imin = i + break + delta = imax-imin + if delta>0: + return (255.*(g-imin)/delta) + else: + return (imax-imin) + +cdef inline np.uint8_t kernel_bottomhat(int* histo, float pop, np.uint8_t g): + cdef int i + + for i in range(256): + if histo[i]: + break + + return (g-i) + + +cdef inline np.uint8_t kernel_egalise(int* histo, float pop, np.uint8_t g): + cdef int i + cdef float sum = 0. + + if pop: + for i in range(256): + sum += histo[i] + if i>=g: + break + + return ((255*sum)/pop) + else: + return (0) + +cdef inline np.uint8_t kernel_gradient(int* histo, float pop, np.uint8_t g): + cdef int i,imin,imax + + + if pop: + for i in range(255,-1,-1): + if histo[i]: + imax = i + break + for i in range(256): + if histo[i]: + imin = i + break + return (imax-imin) + else: + return (0) + +cdef inline np.uint8_t kernel_maximum(int* histo, float pop, np.uint8_t g): + cdef int i + + if pop: + for i in range(255,-1,-1): + if histo[i]: + return (i) + + return (0) + +cdef inline np.uint8_t kernel_mean(int* histo, float pop, np.uint8_t g): + cdef int i + cdef float mean = 0. + + if pop: + for i in range(256): + mean += histo[i]*i + return (mean/pop) + else: + return (0) + +cdef inline np.uint8_t kernel_meansubstraction(int* histo, float pop, np.uint8_t g): + cdef int i + cdef float mean = 0. + + if pop: + for i in range(256): + mean += histo[i]*i + return ((g-mean/pop)/2.+127) + else: + return (0) + +cdef inline np.uint8_t kernel_median(int* histo, float pop, np.uint8_t g): + cdef int i + cdef float sum = pop/2.0 + + if pop: + for i in range(256): + if histo[i]: + sum -= histo[i] + if sum<0: + return (i) + + return (0) + +cdef inline np.uint8_t kernel_minimum(int* histo, float pop, np.uint8_t g): + cdef int i + + if pop: + for i in range(256): + if histo[i]: + return (i) + + return (0) + +cdef inline np.uint8_t kernel_modal(int* histo, float pop, np.uint8_t g): + cdef int hmax=0,imax=0 + + if pop: + for i in range(256): + if histo[i]>hmax: + hmax = histo[i] + imax = i + return (imax) + + return (0) + +cdef inline np.uint8_t kernel_morph_contr_enh(int* histo, float pop, np.uint8_t g): + cdef int i,imin,imax + + if pop: + for i in range(255,-1,-1): + if histo[i]: + imax = i + break + for i in range(256): + if histo[i]: + imin = i + break + if imax-g < g-imin: + return (imax) + else: + return (imin) + else: + return (0) + +cdef inline np.uint8_t kernel_pop(int* histo, float pop, np.uint8_t g): + return (pop) + +cdef inline np.uint8_t kernel_threshold(int* histo, float pop, np.uint8_t g): + cdef int i + cdef float mean = 0. + + if pop: + for i in range(256): + mean += histo[i]*i + return (g>(mean/pop)) + else: + return (0) + +cdef inline np.uint8_t kernel_tophat(int* histo, float pop, np.uint8_t g): + cdef int i + + for i in range(255,-1,-1): + if histo[i]: + break + + return (i-g) + +# ----------------------------------------------------------------- +# python wrappers +# ----------------------------------------------------------------- +def autolevel(np.ndarray[np.uint8_t, ndim=2] image, + np.ndarray[np.uint8_t, ndim=2] selem, + np.ndarray[np.uint8_t, ndim=2] mask=None, + np.ndarray[np.uint8_t, ndim=2] out=None, + char shift_x=0, char shift_y=0): + """bottom hat + """ + return rank8(kernel_autolevel,image,selem,mask,out,shift_x,shift_y) + +def bottomhat(np.ndarray[np.uint8_t, ndim=2] image, + np.ndarray[np.uint8_t, ndim=2] selem, + np.ndarray[np.uint8_t, ndim=2] mask=None, + np.ndarray[np.uint8_t, ndim=2] out=None, + char shift_x=0, char shift_y=0): + """bottom hat + """ + return rank8(kernel_bottomhat,image,selem,mask,out,shift_x,shift_y) + +def egalise(np.ndarray[np.uint8_t, ndim=2] image, + np.ndarray[np.uint8_t, ndim=2] selem, + np.ndarray[np.uint8_t, ndim=2] mask=None, + np.ndarray[np.uint8_t, ndim=2] out=None, + char shift_x=0, char shift_y=0): + """local egalisation of the gray level + """ + return rank8(kernel_egalise,image,selem,mask,out,shift_x,shift_y) + +def gradient(np.ndarray[np.uint8_t, ndim=2] image, + np.ndarray[np.uint8_t, ndim=2] selem, + np.ndarray[np.uint8_t, ndim=2] mask=None, + np.ndarray[np.uint8_t, ndim=2] out=None, + char shift_x=0, char shift_y=0): + """local maximum - local minimum gray level + """ + return rank8(kernel_gradient,image,selem,mask,out,shift_x,shift_y) + +def maximum(np.ndarray[np.uint8_t, ndim=2] image, + np.ndarray[np.uint8_t, ndim=2] selem, + np.ndarray[np.uint8_t, ndim=2] mask=None, + np.ndarray[np.uint8_t, ndim=2] out=None, + char shift_x=0, char shift_y=0): + """local maximum gray level + """ + return rank8(kernel_maximum,image,selem,mask,out,shift_x,shift_y) + +def mean(np.ndarray[np.uint8_t, ndim=2] image, + np.ndarray[np.uint8_t, ndim=2] selem, + np.ndarray[np.uint8_t, ndim=2] mask=None, + np.ndarray[np.uint8_t, ndim=2] out=None, + char shift_x=0, char shift_y=0): + """average gray level (clipped on uint8) + """ + return rank8(kernel_mean,image,selem,mask,out,shift_x,shift_y) + +def meansubstraction(np.ndarray[np.uint8_t, ndim=2] image, + np.ndarray[np.uint8_t, ndim=2] selem, + np.ndarray[np.uint8_t, ndim=2] mask=None, + np.ndarray[np.uint8_t, ndim=2] out=None, + char shift_x=0, char shift_y=0): + """(g - average gray level)/2+127 (clipped on uint8) + """ + return rank8(kernel_meansubstraction,image,selem,mask,out,shift_x,shift_y) + +def median(np.ndarray[np.uint8_t, ndim=2] image, + np.ndarray[np.uint8_t, ndim=2] selem, + np.ndarray[np.uint8_t, ndim=2] mask=None, + np.ndarray[np.uint8_t, ndim=2] out=None, + char shift_x=0, char shift_y=0): + """local median + """ + return rank8(kernel_median,image,selem,mask,out,shift_x,shift_y) + +def minimum(np.ndarray[np.uint8_t, ndim=2] image, + np.ndarray[np.uint8_t, ndim=2] selem, + np.ndarray[np.uint8_t, ndim=2] mask=None, + np.ndarray[np.uint8_t, ndim=2] out=None, + char shift_x=0, char shift_y=0): + """local minimum gray level + """ + return rank8(kernel_minimum,image,selem,mask,out,shift_x,shift_y) + +def morph_contr_enh(np.ndarray[np.uint8_t, ndim=2] image, + np.ndarray[np.uint8_t, ndim=2] selem, + np.ndarray[np.uint8_t, ndim=2] mask=None, + np.ndarray[np.uint8_t, ndim=2] out=None, + char shift_x=0, char shift_y=0): + """morphological contrast enhancement + """ + return rank8(kernel_morph_contr_enh,image,selem,mask,out,shift_x,shift_y) + +def modal(np.ndarray[np.uint8_t, ndim=2] image, + np.ndarray[np.uint8_t, ndim=2] selem, + np.ndarray[np.uint8_t, ndim=2] mask=None, + np.ndarray[np.uint8_t, ndim=2] out=None, + char shift_x=0, char shift_y=0): + """local mode + """ + return rank8(kernel_modal,image,selem,mask,out,shift_x,shift_y) + +def pop(np.ndarray[np.uint8_t, ndim=2] image, + np.ndarray[np.uint8_t, ndim=2] selem, + np.ndarray[np.uint8_t, ndim=2] mask=None, + np.ndarray[np.uint8_t, ndim=2] out=None, + char shift_x=0, char shift_y=0): + """returns the number of actual pixels of the structuring element inside the mask + """ + return rank8(kernel_pop,image,selem,mask,out,shift_x,shift_y) + +def threshold(np.ndarray[np.uint8_t, ndim=2] image, + np.ndarray[np.uint8_t, ndim=2] selem, + np.ndarray[np.uint8_t, ndim=2] mask=None, + np.ndarray[np.uint8_t, ndim=2] out=None, + char shift_x=0, char shift_y=0): + """returns 255 if gray level higher than local mean, 0 else + """ + return rank8(kernel_threshold,image,selem,mask,out,shift_x,shift_y) + +def tophat(np.ndarray[np.uint8_t, ndim=2] image, + np.ndarray[np.uint8_t, ndim=2] selem, + np.ndarray[np.uint8_t, ndim=2] mask=None, + np.ndarray[np.uint8_t, ndim=2] out=None, + char shift_x=0, char shift_y=0): + """top hat + """ + return rank8(kernel_tophat,image,selem,mask,out,shift_x,shift_y) + diff --git a/skimage/rank/crank16.pyx b/skimage/rank/crank16.pyx new file mode 100644 index 00000000..78e5a2af --- /dev/null +++ b/skimage/rank/crank16.pyx @@ -0,0 +1,323 @@ +""" to compile this use: +>>> python setup.py build_ext --inplace + +to generate html report use: +>>> cython -a crank.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 core cimport rank16 + +# todo +# - manage float output, +# - manage different bit depth input +# - add auxiliary parameters (spectral_interval, infSup) + +# ----------------------------------------------------------------- +# 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): + cdef int i + cdef float mean = 0. + + if pop: + for i in range(maxbin): + mean += histo[i]*i + return (mean/pop) + 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): + return (pop) + +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): + """bottom hat + """ + return rank16(kernel_autolevel,image,selem,mask,out,shift_x,shift_y,bitdepth) + +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): + """bottom hat + """ + return rank16(kernel_bottomhat,image,selem,mask,out,shift_x,shift_y,bitdepth) + +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): + """local egalisation of the gray level + """ + return rank16(kernel_egalise,image,selem,mask,out,shift_x,shift_y,bitdepth) + +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): + """local maximum - local minimum gray level + """ + return rank16(kernel_gradient,image,selem,mask,out,shift_x,shift_y,bitdepth) + +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): + """local maximum gray level + """ + return rank16(kernel_maximum,image,selem,mask,out,shift_x,shift_y,bitdepth) + +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): + """average gray level (clipped on uint8) + """ + return rank16(kernel_mean,image,selem,mask,out,shift_x,shift_y,bitdepth) + +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): + """(g - average gray level)/2+midbin (clipped on uint8) + """ + return rank16(kernel_meansubstraction,image,selem,mask,out,shift_x,shift_y,bitdepth) + +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): + """local median + """ + return rank16(kernel_median,image,selem,mask,out,shift_x,shift_y,bitdepth) + +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): + """local minimum gray level + """ + return rank16(kernel_minimum,image,selem,mask,out,shift_x,shift_y,bitdepth) + +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): + """morphological contrast enhancement + """ + return rank16(kernel_morph_contr_enh,image,selem,mask,out,shift_x,shift_y,bitdepth) + +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): + """local mode + """ + return rank16(kernel_modal,image,selem,mask,out,shift_x,shift_y,bitdepth) + +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): + """returns the number of actual pixels of the structuring element inside the mask + """ + return rank16(kernel_pop,image,selem,mask,out,shift_x,shift_y,bitdepth) + +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): + """returns maxbin-1 if gray level higher than local mean, 0 else + """ + return rank16(kernel_threshold,image,selem,mask,out,shift_x,shift_y,bitdepth) + +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): + """top hat + """ + return rank16(kernel_tophat,image,selem,mask,out,shift_x,shift_y,bitdepth) diff --git a/skimage/rank/crank16_percentiles.pyx b/skimage/rank/crank16_percentiles.pyx new file mode 100644 index 00000000..b967c61a --- /dev/null +++ b/skimage/rank/crank16_percentiles.pyx @@ -0,0 +1,268 @@ +""" to compile this use: +>>> python setup.py build_ext --inplace + +to generate html report use: +>>> cython -a crank_percentiles.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 core cimport rank16_percentile + +# todo +# - manage float output, +# - manage different bit depth input +# - add auxiliary parameters (spectral_interval, infSup) + +# ----------------------------------------------------------------- +# kernels uint8 (SOFT version using percentiles) +# ----------------------------------------------------------------- + +cdef inline np.uint16_t kernel_autolevel(int* histo, float pop, np.uint16_t g,int bitdepth,int maxbin, int midbin, float p0, float p1): + cdef int i,imin,imax,sum,delta + + if pop: + sum = 0 + p1 = 1.0-p1 + for i in range(maxbin): + sum += histo[i] + if sum>=p0*pop: + imin = i + break + sum = 0 + for i in range(maxbin-1,-1,-1): + sum += histo[i] + if sum>=p1*pop: + imax = i + break + + delta = imax-imin + if g>imax: + return (maxbin-1) + if g(0) + if delta>0: + return ((maxbin-1)*1.*(g-imin)/delta) + else: + return (0) + else: + return (0) + + +cdef inline np.uint16_t kernel_gradient(int* histo, float pop, np.uint16_t g,int bitdepth,int maxbin, int midbin, float p0, float p1): + cdef int i,imin,imax,sum,delta + + if pop: + sum = 0 + p1 = 1.0-p1 + for i in range(maxbin): + sum += histo[i] + if sum>=p0*pop: + imin = i + break + sum = 0 + for i in range((maxbin-1),-1,-1): + sum += histo[i] + if sum>=p1*pop: + imax = i + break + + return (imax-imin) + else: + return (0) + + +cdef inline np.uint16_t kernel_mean(int* histo, float pop, np.uint16_t g,int bitdepth,int maxbin, int midbin, float p0, float p1): + cdef int i,sum,mean,n + + if pop: + sum = 0 + mean = 0 + n = 0 + for i in range(maxbin): + sum += histo[i] + if (sum>=p0*pop) and (sum<=p1*pop): + n += histo[i] + mean += histo[i]*i + + if n>0: + return (1.0*mean/n) + else: + return (0) + else: + return (0) + +cdef inline np.uint16_t kernel_mean_substraction(int* histo, float pop, np.uint16_t g,int bitdepth,int maxbin, int midbin, float p0, float p1): + cdef int i,sum,mean,n + + if pop: + sum = 0 + mean = 0 + n = 0 + for i in range(maxbin): + sum += histo[i] + if (sum>=p0*pop) and (sum<=p1*pop): + n += histo[i] + mean += histo[i]*i + if n>0: + return ((g-(mean/n))*.5+midbin) + else: + return (0) + else: + 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, float p0, float p1): + cdef int i,imin,imax,sum,delta + + if pop: + sum = 0 + p1 = 1.0-p1 + for i in range(maxbin): + sum += histo[i] + if sum>=p0*pop: + imin = i + break + sum = 0 + for i in range((maxbin-1),-1,-1): + sum += histo[i] + if sum>=p1*pop: + imax = i + break + if g>imax: + return imax + if gimin + if imax-g < g-imin: + return imax + else: + return imin + else: + return (0) + +cdef inline np.uint16_t kernel_percentile(int* histo, float pop, np.uint16_t g,int bitdepth,int maxbin, int midbin, float p0, float p1): + cdef int i + cdef float sum = 0. + + if pop: + for i in range(maxbin): + sum += histo[i] + if sum>=p0*pop: + break + + return (i) + else: + return (0) + +cdef inline np.uint16_t kernel_pop(int* histo, float pop, np.uint16_t g,int bitdepth,int maxbin, int midbin, float p0, float p1): + cdef int i,sum,n + + if pop: + sum = 0 + n = 0 + for i in range(maxbin): + sum += histo[i] + if (sum>=p0*pop) and (sum<=p1*pop): + n += histo[i] + return (n) + else: + return (0) + +cdef inline np.uint16_t kernel_threshold(int* histo, float pop, np.uint16_t g,int bitdepth,int maxbin, int midbin, float p0, float p1): + cdef int i + cdef float sum = 0. + + if pop: + for i in range(maxbin): + sum += histo[i] + if sum>=p0*pop: + break + + return ((maxbin-1)*(g>=i)) + else: + return (0) + +# ----------------------------------------------------------------- +# 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, float p0=0., float p1=0.): + """bottom hat + """ + return rank16_percentile(kernel_autolevel,image,selem,mask,out,shift_x,shift_y,bitdepth,p0,p1) + + +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, float p0=0., float p1=0.): + """return p0,p1 percentile gradient + """ + return rank16_percentile(kernel_gradient,image,selem,mask,out,shift_x,shift_y,bitdepth,p0,p1) + +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, float p0=0., float p1=0.): + """return mean between [p0 and p1] percentiles + """ + return rank16_percentile(kernel_mean,image,selem,mask,out,shift_x,shift_y,bitdepth,p0,p1) + +def mean_substraction(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, float p0=0., float p1=0.): + """return original - mean between [p0 and p1] percentiles *.5 +127 + """ + return rank16_percentile(kernel_mean_substraction,image,selem,mask,out,shift_x,shift_y,bitdepth,p0,p1) + +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, float p0=0., float p1=0.): + """reforce contrast using percentiles + """ + return rank16_percentile(kernel_morph_contr_enh,image,selem,mask,out,shift_x,shift_y,bitdepth,p0,p1) + + +def percentile(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, float p0=0., float p1=0.): + """return p0 percentile + """ + return rank16_percentile(kernel_percentile,image,selem,mask,out,shift_x,shift_y,bitdepth,p0,p1) + + +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, float p0=0., float p1=0.): + """return nb of pixels between [p0 and p1] + """ + return rank16_percentile(kernel_pop,image,selem,mask,out,shift_x,shift_y,bitdepth,p0,p1) + +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, float p0=0., float p1=0.): + """return (maxbin-1) if g > percentile p0 + """ + return rank16_percentile(kernel_threshold,image,selem,mask,out,shift_x,shift_y,bitdepth,p0,p1) diff --git a/skimage/rank/crank_percentiles.pyx b/skimage/rank/crank_percentiles.pyx new file mode 100644 index 00000000..51a9e01f --- /dev/null +++ b/skimage/rank/crank_percentiles.pyx @@ -0,0 +1,263 @@ +""" to compile this use: +>>> python setup.py build_ext --inplace + +to generate html report use: +>>> cython -a crank_percentiles.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 core cimport rank8_percentile + +# todo +# - manage float output, +# - manage different bit depth input +# - add auxiliary parameters (spectral_interval, infSup) + +# ----------------------------------------------------------------- +# kernels uint8 (SOFT version using percentiles) +# ----------------------------------------------------------------- + +cdef inline np.uint8_t kernel_autolevel(int* histo, float pop, np.uint8_t g, float p0, float p1): + cdef int i,imin,imax,sum,delta + + if pop: + sum = 0 + p1 = 1.0-p1 + for i in range(256): + sum += histo[i] + if sum>=p0*pop: + imin = i + break + sum = 0 + for i in range(255,-1,-1): + sum += histo[i] + if sum>=p1*pop: + imax = i + break + + delta = imax-imin + if delta>0: + return (255.*(g-imin)/delta) + else: + return (0) + else: + return (0) + + +cdef inline np.uint8_t kernel_gradient(int* histo, float pop, np.uint8_t g, float p0, float p1): + cdef int i,imin,imax,sum,delta + + if pop: + sum = 0 + p1 = 1.0-p1 + for i in range(256): + sum += histo[i] + if sum>=p0*pop: + imin = i + break + sum = 0 + for i in range(255,-1,-1): + sum += histo[i] + if sum>=p1*pop: + imax = i + break + + return (imax-imin) + else: + return (0) + + +cdef inline np.uint8_t kernel_mean(int* histo, float pop, np.uint8_t g, float p0, float p1): + cdef int i,sum,mean,n + + if pop: + sum = 0 + mean = 0 + n = 0 + for i in range(256): + sum += histo[i] + if (sum>=p0*pop) and (sum<=p1*pop): + n += histo[i] + mean += histo[i]*i + if n>0: + return (1.0*mean/n) + else: + return (0) + else: + return (0) + +cdef inline np.uint8_t kernel_mean_substraction(int* histo, float pop, np.uint8_t g, float p0, float p1): + cdef int i,sum,mean,n + + if pop: + sum = 0 + mean = 0 + n = 0 + for i in range(256): + sum += histo[i] + if (sum>=p0*pop) and (sum<=p1*pop): + n += histo[i] + mean += histo[i]*i + if n>0: + return ((g-(mean/n))*.5+127) + else: + return (0) + else: + return (0) + +cdef inline np.uint8_t kernel_morph_contr_enh(int* histo, float pop, np.uint8_t g, float p0, float p1): + cdef int i,imin,imax,sum,delta + + if pop: + sum = 0 + p1 = 1.0-p1 + for i in range(256): + sum += histo[i] + if sum>=p0*pop: + imin = i + break + sum = 0 + for i in range(255,-1,-1): + sum += histo[i] + if sum>=p1*pop: + imax = i + break + if g>imax: + return imax + if gimin + if imax-g < g-imin: + return imax + else: + return imin + else: + return (0) + +cdef inline np.uint8_t kernel_percentile(int* histo, float pop, np.uint8_t g, float p0, float p1): + cdef int i + cdef float sum = 0. + + if pop: + for i in range(256): + sum += histo[i] + if sum>=p0*pop: + break + + return (i) + else: + return (0) + +cdef inline np.uint8_t kernel_pop(int* histo, float pop, np.uint8_t g, float p0, float p1): + cdef int i,sum,n + + if pop: + sum = 0 + n = 0 + for i in range(256): + sum += histo[i] + if (sum>=p0*pop) and (sum<=p1*pop): + n += histo[i] + return (n) + else: + return (0) + +cdef inline np.uint8_t kernel_threshold(int* histo, float pop, np.uint8_t g, float p0, float p1): + cdef int i + cdef float sum = 0. + + if pop: + for i in range(256): + sum += histo[i] + if sum>=p0*pop: + break + + return (255*(g>=i)) + else: + return (0) + +# ----------------------------------------------------------------- +# python wrappers +# ----------------------------------------------------------------- +def autolevel(np.ndarray[np.uint8_t, ndim=2] image, + np.ndarray[np.uint8_t, ndim=2] selem, + np.ndarray[np.uint8_t, ndim=2] mask=None, + np.ndarray[np.uint8_t, ndim=2] out=None, + char shift_x=0, char shift_y=0, float p0=0., float p1=0.): + """bottom hat + """ + return rank8_percentile(kernel_autolevel,image,selem,mask,out,shift_x,shift_y,p0,p1) + + +def gradient(np.ndarray[np.uint8_t, ndim=2] image, + np.ndarray[np.uint8_t, ndim=2] selem, + np.ndarray[np.uint8_t, ndim=2] mask=None, + np.ndarray[np.uint8_t, ndim=2] out=None, + char shift_x=0, char shift_y=0, float p0=0., float p1=0.): + """return p0,p1 percentile gradient + """ + return rank8_percentile(kernel_gradient,image,selem,mask,out,shift_x,shift_y,p0,p1) + +def mean(np.ndarray[np.uint8_t, ndim=2] image, + np.ndarray[np.uint8_t, ndim=2] selem, + np.ndarray[np.uint8_t, ndim=2] mask=None, + np.ndarray[np.uint8_t, ndim=2] out=None, + char shift_x=0, char shift_y=0, float p0=0., float p1=0.): + """return mean between [p0 and p1] percentiles + """ + return rank8_percentile(kernel_mean,image,selem,mask,out,shift_x,shift_y,p0,p1) + +def mean_substraction(np.ndarray[np.uint8_t, ndim=2] image, + np.ndarray[np.uint8_t, ndim=2] selem, + np.ndarray[np.uint8_t, ndim=2] mask=None, + np.ndarray[np.uint8_t, ndim=2] out=None, + char shift_x=0, char shift_y=0, float p0=0., float p1=0.): + """return original - mean between [p0 and p1] percentiles *.5 +127 + """ + return rank8_percentile(kernel_mean_substraction,image,selem,mask,out,shift_x,shift_y,p0,p1) + +def morph_contr_enh(np.ndarray[np.uint8_t, ndim=2] image, + np.ndarray[np.uint8_t, ndim=2] selem, + np.ndarray[np.uint8_t, ndim=2] mask=None, + np.ndarray[np.uint8_t, ndim=2] out=None, + char shift_x=0, char shift_y=0, float p0=0., float p1=0.): + """reforce contrast using percentiles + """ + return rank8_percentile(kernel_morph_contr_enh,image,selem,mask,out,shift_x,shift_y,p0,p1) + + +def percentile(np.ndarray[np.uint8_t, ndim=2] image, + np.ndarray[np.uint8_t, ndim=2] selem, + np.ndarray[np.uint8_t, ndim=2] mask=None, + np.ndarray[np.uint8_t, ndim=2] out=None, + char shift_x=0, char shift_y=0, float p0=0., float p1=0.): + """return p0 percentile + """ + return rank8_percentile(kernel_percentile,image,selem,mask,out,shift_x,shift_y,p0,p1) + + +def pop(np.ndarray[np.uint8_t, ndim=2] image, + np.ndarray[np.uint8_t, ndim=2] selem, + np.ndarray[np.uint8_t, ndim=2] mask=None, + np.ndarray[np.uint8_t, ndim=2] out=None, + char shift_x=0, char shift_y=0, float p0=0., float p1=0.): + """return nb of pixels between [p0 and p1] + """ + return rank8_percentile(kernel_pop,image,selem,mask,out,shift_x,shift_y,p0,p1) + +def threshold(np.ndarray[np.uint8_t, ndim=2] image, + np.ndarray[np.uint8_t, ndim=2] selem, + np.ndarray[np.uint8_t, ndim=2] mask=None, + np.ndarray[np.uint8_t, ndim=2] out=None, + char shift_x=0, char shift_y=0, float p0=0., float p1=0.): + """return 255 if g > percentile p0 + """ + return rank8_percentile(kernel_threshold,image,selem,mask,out,shift_x,shift_y,p0,p1) diff --git a/skimage/rank/setup.py b/skimage/rank/setup.py new file mode 100644 index 00000000..6d18e27f --- /dev/null +++ b/skimage/rank/setup.py @@ -0,0 +1,15 @@ +import numpy as np + +from distutils.core import setup +from distutils.extension import Extension +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()]), + 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_percentiles", ["crank16_percentiles.pyx"], include_dirs=[np.get_include()])] +) \ No newline at end of file diff --git a/skimage/rank/tools.py b/skimage/rank/tools.py new file mode 100644 index 00000000..835a8884 --- /dev/null +++ b/skimage/rank/tools.py @@ -0,0 +1,52 @@ +__author__ = 'Olivier Debeir 2021' + +import logging +import time + +def init_logger(logfilename = 'myapp.log'): + """add logger capabilities + """ + FORMAT = '%(asctime)-15s %(processName)s %(process)d %(message)s' + logging.basicConfig(filename=logfilename,format=FORMAT,filemode='wt') + logger = logging.getLogger() + logger.setLevel(logging.DEBUG) + + # create console handler and set level to debug + ch = logging.StreamHandler() + ch.setLevel(logging.DEBUG) + + # add ch to logger + logger.addHandler(ch) + logger.info('start logging in %s' % logfilename) + return logger + + +logger = logging.getLogger() + +def log_timing(func): + + def wrapper(*arg): + log_timing.level += 1 + t1 = time.time() + res = func(*arg) + t2 = time.time() + ms = (t2-t1)*1000.0 + logger.info('%s%s took %0.3f ms' % (log_timing.level*'-',func.func_name, ms)) + log_timing.level -= 1 + return (res,ms) + + return wrapper + +log_timing.level = 0 + +def tumbnail_it(data): + """display image with its histogram + """ + h = np.histogram(data[:],100) + hn = 512*h[0]/np.max(h[0]) + + plt.subplot(1,2,1) + plt.imshow(ima8,interpolation='nearest',cmap=cm.gray) + plt.subplot(1,2,2) + plt.plot(hn) + plt.colorbar()