diff --git a/skimage/rank/_core16.pxd b/skimage/rank/_core16.pxd index 1a3194de..d00ef37e 100644 --- a/skimage/rank/_core16.pxd +++ b/skimage/rank/_core16.pxd @@ -1,18 +1,4 @@ -""" to compile this use: ->>> python setup.py build_ext --inplace - -to generate html report use: ->>> cython -a core16.pxd -""" - -#cython: cdivision=True -#cython: boundscheck=False -#cython: nonecheck=False -#cython: wraparound=False - -import numpy as np cimport numpy as np -from libc.stdlib cimport malloc, free #--------------------------------------------------------------------------- # 16 bit core kernel receives extra information about data bitdepth @@ -23,261 +9,4 @@ np.ndarray[np.uint16_t, ndim=2] image, np.ndarray[np.uint8_t, ndim=2] selem, np.ndarray[np.uint8_t, ndim=2] mask, np.ndarray[np.uint16_t, ndim=2] out, -char shift_x, char shift_y,int bitdepth): - """ Main loop, this function computes the histogram for each image point - - data is uint8 - - result is uint8 casted - """ - - cdef int rows = image.shape[0] - cdef int cols = image.shape[1] - cdef int srows = selem.shape[0] - cdef int scols = selem.shape[1] - - cdef int centre_r = int(selem.shape[0] / 2) + shift_y - cdef int centre_c = int(selem.shape[1] / 2) + shift_x - - # check that structuring element center is inside the element bounding box - assert centre_r >= 0 - assert centre_c >= 0 - assert centre_r < srows - assert centre_c < scols - assert bitdepth in range(2,13) - - maxbin_list = [0,0,4,8,16,32,64,128,256,512,1024,2048,4096] - midbin_list = [0,0,2,4,8,16,32,64,128,256,512,1024,2048] - - - #set maxbin and midbin - cdef int maxbin=maxbin_list[bitdepth],midbin=midbin_list[bitdepth] - - assert (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 - - # number of element in each attack border - cdef int num_se_n, num_se_s, num_se_e, num_se_w - - # the current local histogram distribution - cdef int* histo = malloc(maxbin * sizeof(int)) - - # these lists contain the relative pixel row and column for each of the 4 attack borders - # east, west, north and south - # e.g. se_e_r lists the rows of the east structuring element border - - cdef int* se_e_r = 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 - - num_se_n = num_se_s = num_se_e = num_se_w = 0 - - for r in range(srows): - for c in range(scols): - if t_e[r,c]: - se_e_r[num_se_e] = r - centre_r - se_e_c[num_se_e] = c - centre_c - num_se_e += 1 - if t_w[r,c]: - se_w_r[num_se_w] = r - centre_r - se_w_c[num_se_w] = c - centre_c - num_se_w += 1 - if t_n[r,c]: - se_n_r[num_se_n] = r - centre_r - se_n_c[num_se_n] = c - centre_c - num_se_n += 1 - if t_s[r,c]: - se_s_r[num_se_s] = r - centre_r - se_s_c[num_se_s] = c - centre_c - num_se_s += 1 - - # initial population and histogram - for i in range(maxbin): - histo[i] = 0 - - pop = 0 - - for r in range(srows): - for c in range(scols): - rr = r - cc = c - if selem[r, c]: - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] += 1 - pop += 1. - - r = 0 - c = 0 - # kernel ------------------------------------------- - out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c], - bitdepth,maxbin,midbin) - # kernel ------------------------------------------- - - # main loop - r = 0 - for even_row in range(0,rows,2): - # ---> west to east - for c in range(1,cols): - for s in range(num_se_e): - rr = r + se_e_r[s] + centre_r - cc = c + se_e_c[s] + centre_c - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] += 1 - pop += 1. - for s in range(num_se_w): - rr = r + se_w_r[s] + centre_r - cc = c + se_w_c[s] + centre_c - 1 - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] -= 1 - pop -= 1. - - # kernel ------------------------------------------- - out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c], - bitdepth,maxbin,midbin) - # kernel ------------------------------------------- - - r += 1 # pass to the next row - if r>=rows: - break - - # ---> north to south - for s in range(num_se_s): - rr = r + se_s_r[s] + centre_r - cc = c + se_s_c[s] + centre_c - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] += 1 - pop += 1. - for s in range(num_se_n): - rr = r + se_n_r[s] + centre_r - 1 - cc = c + se_n_c[s] + centre_c - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] -= 1 - pop -= 1. - - # kernel ------------------------------------------- - out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c], - bitdepth,maxbin,midbin) - # kernel ------------------------------------------- - - # ---> east to west - for c in range(cols-2,-1,-1): - for s in range(num_se_w): - rr = r + se_w_r[s] + centre_r - cc = c + se_w_c[s] + centre_c - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] += 1 - pop += 1. - for s in range(num_se_e): - rr = r + se_e_r[s] + centre_r - cc = c + se_e_c[s] + centre_c + 1 - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] -= 1 - pop -= 1. - - # kernel ------------------------------------------- - out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c], - bitdepth,maxbin,midbin) - # kernel ------------------------------------------- - - r += 1 # pass to the next row - if r>=rows: - break - - # ---> north to south - for s in range(num_se_s): - rr = r + se_s_r[s] + centre_r - cc = c + se_s_c[s] + centre_c - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] += 1 - pop += 1. - for s in range(num_se_n): - rr = r + se_n_r[s] + centre_r - 1 - cc = c + se_n_c[s] + centre_c - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] -= 1 - pop -= 1. - - # kernel ------------------------------------------- - out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c], - bitdepth,maxbin,midbin) - # kernel ------------------------------------------- - - # release memory allocated by malloc - - free(se_e_r) - free(se_e_c) - free(se_w_r) - free(se_w_c) - free(se_n_r) - free(se_n_c) - free(se_s_r) - free(se_s_c) - - free(histo) - - return out +char shift_x, char shift_y,int bitdepth) \ No newline at end of file diff --git a/skimage/rank/_core16.pyx b/skimage/rank/_core16.pyx new file mode 100644 index 00000000..1a3194de --- /dev/null +++ b/skimage/rank/_core16.pyx @@ -0,0 +1,283 @@ +""" to compile this use: +>>> python setup.py build_ext --inplace + +to generate html report use: +>>> cython -a core16.pxd +""" + +#cython: cdivision=True +#cython: boundscheck=False +#cython: nonecheck=False +#cython: wraparound=False + +import numpy as np +cimport numpy as np +from libc.stdlib cimport malloc, free + +#--------------------------------------------------------------------------- +# 16 bit core kernel receives extra information about data bitdepth +#--------------------------------------------------------------------------- + +cdef inline _core16(np.uint16_t kernel(int*, float, np.uint16_t, int ,int,int ), +np.ndarray[np.uint16_t, ndim=2] image, +np.ndarray[np.uint8_t, ndim=2] selem, +np.ndarray[np.uint8_t, ndim=2] mask, +np.ndarray[np.uint16_t, ndim=2] out, +char shift_x, char shift_y,int bitdepth): + """ Main loop, this function computes the histogram for each image point + - data is uint8 + - result is uint8 casted + """ + + cdef int rows = image.shape[0] + cdef int cols = image.shape[1] + cdef int srows = selem.shape[0] + cdef int scols = selem.shape[1] + + cdef int centre_r = int(selem.shape[0] / 2) + shift_y + cdef int centre_c = int(selem.shape[1] / 2) + shift_x + + # check that structuring element center is inside the element bounding box + assert centre_r >= 0 + assert centre_c >= 0 + assert centre_r < srows + assert centre_c < scols + assert bitdepth in range(2,13) + + maxbin_list = [0,0,4,8,16,32,64,128,256,512,1024,2048,4096] + midbin_list = [0,0,2,4,8,16,32,64,128,256,512,1024,2048] + + + #set maxbin and midbin + cdef int maxbin=maxbin_list[bitdepth],midbin=midbin_list[bitdepth] + + assert (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 + + # number of element in each attack border + cdef int num_se_n, num_se_s, num_se_e, num_se_w + + # the current local histogram distribution + cdef int* histo = malloc(maxbin * sizeof(int)) + + # these lists contain the relative pixel row and column for each of the 4 attack borders + # east, west, north and south + # e.g. se_e_r lists the rows of the east structuring element border + + cdef int* se_e_r = 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 + + num_se_n = num_se_s = num_se_e = num_se_w = 0 + + for r in range(srows): + for c in range(scols): + if t_e[r,c]: + se_e_r[num_se_e] = r - centre_r + se_e_c[num_se_e] = c - centre_c + num_se_e += 1 + if t_w[r,c]: + se_w_r[num_se_w] = r - centre_r + se_w_c[num_se_w] = c - centre_c + num_se_w += 1 + if t_n[r,c]: + se_n_r[num_se_n] = r - centre_r + se_n_c[num_se_n] = c - centre_c + num_se_n += 1 + if t_s[r,c]: + se_s_r[num_se_s] = r - centre_r + se_s_c[num_se_s] = c - centre_c + num_se_s += 1 + + # initial population and histogram + for i in range(maxbin): + histo[i] = 0 + + pop = 0 + + for r in range(srows): + for c in range(scols): + rr = r + cc = c + if selem[r, c]: + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] += 1 + pop += 1. + + r = 0 + c = 0 + # kernel ------------------------------------------- + out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c], + bitdepth,maxbin,midbin) + # kernel ------------------------------------------- + + # main loop + r = 0 + for even_row in range(0,rows,2): + # ---> west to east + for c in range(1,cols): + for s in range(num_se_e): + rr = r + se_e_r[s] + centre_r + cc = c + se_e_c[s] + centre_c + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] += 1 + pop += 1. + for s in range(num_se_w): + rr = r + se_w_r[s] + centre_r + cc = c + se_w_c[s] + centre_c - 1 + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] -= 1 + pop -= 1. + + # kernel ------------------------------------------- + out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c], + bitdepth,maxbin,midbin) + # kernel ------------------------------------------- + + r += 1 # pass to the next row + if r>=rows: + break + + # ---> north to south + for s in range(num_se_s): + rr = r + se_s_r[s] + centre_r + cc = c + se_s_c[s] + centre_c + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] += 1 + pop += 1. + for s in range(num_se_n): + rr = r + se_n_r[s] + centre_r - 1 + cc = c + se_n_c[s] + centre_c + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] -= 1 + pop -= 1. + + # kernel ------------------------------------------- + out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c], + bitdepth,maxbin,midbin) + # kernel ------------------------------------------- + + # ---> east to west + for c in range(cols-2,-1,-1): + for s in range(num_se_w): + rr = r + se_w_r[s] + centre_r + cc = c + se_w_c[s] + centre_c + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] += 1 + pop += 1. + for s in range(num_se_e): + rr = r + se_e_r[s] + centre_r + cc = c + se_e_c[s] + centre_c + 1 + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] -= 1 + pop -= 1. + + # kernel ------------------------------------------- + out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c], + bitdepth,maxbin,midbin) + # kernel ------------------------------------------- + + r += 1 # pass to the next row + if r>=rows: + break + + # ---> north to south + for s in range(num_se_s): + rr = r + se_s_r[s] + centre_r + cc = c + se_s_c[s] + centre_c + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] += 1 + pop += 1. + for s in range(num_se_n): + rr = r + se_n_r[s] + centre_r - 1 + cc = c + se_n_c[s] + centre_c + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] -= 1 + pop -= 1. + + # kernel ------------------------------------------- + out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c], + bitdepth,maxbin,midbin) + # kernel ------------------------------------------- + + # release memory allocated by malloc + + free(se_e_r) + free(se_e_c) + free(se_w_r) + free(se_w_c) + free(se_n_r) + free(se_n_c) + free(se_s_r) + free(se_s_c) + + free(histo) + + return out diff --git a/skimage/rank/_core16b.pxd b/skimage/rank/_core16b.pxd index 163662b5..b972ae9a 100644 --- a/skimage/rank/_core16b.pxd +++ b/skimage/rank/_core16b.pxd @@ -1,18 +1,4 @@ -""" to compile this use: ->>> python setup.py build_ext --inplace - -to generate html report use: ->>> cython -a core16b.pxd -""" - -#cython: cdivision=True -#cython: boundscheck=False -#cython: nonecheck=False -#cython: wraparound=False - -import numpy as np cimport numpy as np -from libc.stdlib cimport malloc, free #--------------------------------------------------------------------------- # 16 bit core kernel receives extra information about data bitdepth and bilateral interval @@ -23,262 +9,4 @@ np.ndarray[np.uint16_t, ndim=2] image, np.ndarray[np.uint8_t, ndim=2] selem, np.ndarray[np.uint8_t, ndim=2] mask, np.ndarray[np.uint16_t, ndim=2] out, -char shift_x, char shift_y,int bitdepth, int s0, int s1): - """ Main loop, this function computes the histogram for each image point - - data is uint8 - - result is uint8 casted - - only pixel inside [s0,s1] centered on g are taken into account - """ - - cdef int rows = image.shape[0] - cdef int cols = image.shape[1] - cdef int srows = selem.shape[0] - cdef int scols = selem.shape[1] - - cdef int centre_r = int(selem.shape[0] / 2) + shift_y - cdef int centre_c = int(selem.shape[1] / 2) + shift_x - - # check that structuring element center is inside the element bounding box - assert centre_r >= 0 - assert centre_c >= 0 - assert centre_r < srows - assert centre_c < scols - assert bitdepth in range(2,13) - - maxbin_list = [0,0,4,8,16,32,64,128,256,512,1024,2048,4096] - midbin_list = [0,0,2,4,8,16,32,64,128,256,512,1024,2048] - - - #set maxbin and midbin - cdef int maxbin=maxbin_list[bitdepth],midbin=midbin_list[bitdepth] - - assert (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 - - # number of element in each attack border - cdef int num_se_n, num_se_s, num_se_e, num_se_w - - # the current local histogram distribution - cdef int* histo = malloc(maxbin * sizeof(int)) - - # these lists contain the relative pixel row and column for each of the 4 attack borders - # east, west, north and south - # e.g. se_e_r lists the rows of the east structuring element border - - cdef int* se_e_r = 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 - - num_se_n = num_se_s = num_se_e = num_se_w = 0 - - for r in range(srows): - for c in range(scols): - if t_e[r,c]: - se_e_r[num_se_e] = r - centre_r - se_e_c[num_se_e] = c - centre_c - num_se_e += 1 - if t_w[r,c]: - se_w_r[num_se_w] = r - centre_r - se_w_c[num_se_w] = c - centre_c - num_se_w += 1 - if t_n[r,c]: - se_n_r[num_se_n] = r - centre_r - se_n_c[num_se_n] = c - centre_c - num_se_n += 1 - if t_s[r,c]: - se_s_r[num_se_s] = r - centre_r - se_s_c[num_se_s] = c - centre_c - num_se_s += 1 - - # initial population and histogram - for i in range(maxbin): - histo[i] = 0 - - pop = 0 - - for r in range(srows): - for c in range(scols): - rr = r - cc = c - if selem[r, c]: - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] += 1 - pop += 1. - - r = 0 - c = 0 - # kernel ------------------------------------------- - out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c], - bitdepth,maxbin,midbin,s0,s1) - # kernel ------------------------------------------- - - # main loop - r = 0 - for even_row in range(0,rows,2): - # ---> west to east - for c in range(1,cols): - for s in range(num_se_e): - rr = r + se_e_r[s] + centre_r - cc = c + se_e_c[s] + centre_c - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] += 1 - pop += 1. - for s in range(num_se_w): - rr = r + se_w_r[s] + centre_r - cc = c + se_w_c[s] + centre_c - 1 - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] -= 1 - pop -= 1. - - # kernel ------------------------------------------- - out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c], - bitdepth,maxbin,midbin,s0,s1) - # kernel ------------------------------------------- - - r += 1 # pass to the next row - if r>=rows: - break - - # ---> north to south - for s in range(num_se_s): - rr = r + se_s_r[s] + centre_r - cc = c + se_s_c[s] + centre_c - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] += 1 - pop += 1. - for s in range(num_se_n): - rr = r + se_n_r[s] + centre_r - 1 - cc = c + se_n_c[s] + centre_c - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] -= 1 - pop -= 1. - - # kernel ------------------------------------------- - out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c], - bitdepth,maxbin,midbin,s0,s1) - # kernel ------------------------------------------- - - # ---> east to west - for c in range(cols-2,-1,-1): - for s in range(num_se_w): - rr = r + se_w_r[s] + centre_r - cc = c + se_w_c[s] + centre_c - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] += 1 - pop += 1. - for s in range(num_se_e): - rr = r + se_e_r[s] + centre_r - cc = c + se_e_c[s] + centre_c + 1 - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] -= 1 - pop -= 1. - - # kernel ------------------------------------------- - out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c], - bitdepth,maxbin,midbin,s0,s1) - # kernel ------------------------------------------- - - r += 1 # pass to the next row - if r>=rows: - break - - # ---> north to south - for s in range(num_se_s): - rr = r + se_s_r[s] + centre_r - cc = c + se_s_c[s] + centre_c - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] += 1 - pop += 1. - for s in range(num_se_n): - rr = r + se_n_r[s] + centre_r - 1 - cc = c + se_n_c[s] + centre_c - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] -= 1 - pop -= 1. - - # kernel ------------------------------------------- - out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c], - bitdepth,maxbin,midbin,s0,s1) - # kernel ------------------------------------------- - - # release memory allocated by malloc - - free(se_e_r) - free(se_e_c) - free(se_w_r) - free(se_w_c) - free(se_n_r) - free(se_n_c) - free(se_s_r) - free(se_s_c) - - free(histo) - - return out +char shift_x, char shift_y,int bitdepth, int s0, int s1) \ No newline at end of file diff --git a/skimage/rank/_core16b.pyx b/skimage/rank/_core16b.pyx new file mode 100644 index 00000000..163662b5 --- /dev/null +++ b/skimage/rank/_core16b.pyx @@ -0,0 +1,284 @@ +""" to compile this use: +>>> python setup.py build_ext --inplace + +to generate html report use: +>>> cython -a core16b.pxd +""" + +#cython: cdivision=True +#cython: boundscheck=False +#cython: nonecheck=False +#cython: wraparound=False + +import numpy as np +cimport numpy as np +from libc.stdlib cimport malloc, free + +#--------------------------------------------------------------------------- +# 16 bit core kernel receives extra information about data bitdepth and bilateral interval +#--------------------------------------------------------------------------- + +cdef inline _core16b(np.uint16_t kernel(int*, float, np.uint16_t, int ,int,int,int,int), +np.ndarray[np.uint16_t, ndim=2] image, +np.ndarray[np.uint8_t, ndim=2] selem, +np.ndarray[np.uint8_t, ndim=2] mask, +np.ndarray[np.uint16_t, ndim=2] out, +char shift_x, char shift_y,int bitdepth, int s0, int s1): + """ Main loop, this function computes the histogram for each image point + - data is uint8 + - result is uint8 casted + - only pixel inside [s0,s1] centered on g are taken into account + """ + + cdef int rows = image.shape[0] + cdef int cols = image.shape[1] + cdef int srows = selem.shape[0] + cdef int scols = selem.shape[1] + + cdef int centre_r = int(selem.shape[0] / 2) + shift_y + cdef int centre_c = int(selem.shape[1] / 2) + shift_x + + # check that structuring element center is inside the element bounding box + assert centre_r >= 0 + assert centre_c >= 0 + assert centre_r < srows + assert centre_c < scols + assert bitdepth in range(2,13) + + maxbin_list = [0,0,4,8,16,32,64,128,256,512,1024,2048,4096] + midbin_list = [0,0,2,4,8,16,32,64,128,256,512,1024,2048] + + + #set maxbin and midbin + cdef int maxbin=maxbin_list[bitdepth],midbin=midbin_list[bitdepth] + + assert (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 + + # number of element in each attack border + cdef int num_se_n, num_se_s, num_se_e, num_se_w + + # the current local histogram distribution + cdef int* histo = malloc(maxbin * sizeof(int)) + + # these lists contain the relative pixel row and column for each of the 4 attack borders + # east, west, north and south + # e.g. se_e_r lists the rows of the east structuring element border + + cdef int* se_e_r = 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 + + num_se_n = num_se_s = num_se_e = num_se_w = 0 + + for r in range(srows): + for c in range(scols): + if t_e[r,c]: + se_e_r[num_se_e] = r - centre_r + se_e_c[num_se_e] = c - centre_c + num_se_e += 1 + if t_w[r,c]: + se_w_r[num_se_w] = r - centre_r + se_w_c[num_se_w] = c - centre_c + num_se_w += 1 + if t_n[r,c]: + se_n_r[num_se_n] = r - centre_r + se_n_c[num_se_n] = c - centre_c + num_se_n += 1 + if t_s[r,c]: + se_s_r[num_se_s] = r - centre_r + se_s_c[num_se_s] = c - centre_c + num_se_s += 1 + + # initial population and histogram + for i in range(maxbin): + histo[i] = 0 + + pop = 0 + + for r in range(srows): + for c in range(scols): + rr = r + cc = c + if selem[r, c]: + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] += 1 + pop += 1. + + r = 0 + c = 0 + # kernel ------------------------------------------- + out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c], + bitdepth,maxbin,midbin,s0,s1) + # kernel ------------------------------------------- + + # main loop + r = 0 + for even_row in range(0,rows,2): + # ---> west to east + for c in range(1,cols): + for s in range(num_se_e): + rr = r + se_e_r[s] + centre_r + cc = c + se_e_c[s] + centre_c + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] += 1 + pop += 1. + for s in range(num_se_w): + rr = r + se_w_r[s] + centre_r + cc = c + se_w_c[s] + centre_c - 1 + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] -= 1 + pop -= 1. + + # kernel ------------------------------------------- + out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c], + bitdepth,maxbin,midbin,s0,s1) + # kernel ------------------------------------------- + + r += 1 # pass to the next row + if r>=rows: + break + + # ---> north to south + for s in range(num_se_s): + rr = r + se_s_r[s] + centre_r + cc = c + se_s_c[s] + centre_c + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] += 1 + pop += 1. + for s in range(num_se_n): + rr = r + se_n_r[s] + centre_r - 1 + cc = c + se_n_c[s] + centre_c + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] -= 1 + pop -= 1. + + # kernel ------------------------------------------- + out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c], + bitdepth,maxbin,midbin,s0,s1) + # kernel ------------------------------------------- + + # ---> east to west + for c in range(cols-2,-1,-1): + for s in range(num_se_w): + rr = r + se_w_r[s] + centre_r + cc = c + se_w_c[s] + centre_c + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] += 1 + pop += 1. + for s in range(num_se_e): + rr = r + se_e_r[s] + centre_r + cc = c + se_e_c[s] + centre_c + 1 + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] -= 1 + pop -= 1. + + # kernel ------------------------------------------- + out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c], + bitdepth,maxbin,midbin,s0,s1) + # kernel ------------------------------------------- + + r += 1 # pass to the next row + if r>=rows: + break + + # ---> north to south + for s in range(num_se_s): + rr = r + se_s_r[s] + centre_r + cc = c + se_s_c[s] + centre_c + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] += 1 + pop += 1. + for s in range(num_se_n): + rr = r + se_n_r[s] + centre_r - 1 + cc = c + se_n_c[s] + centre_c + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] -= 1 + pop -= 1. + + # kernel ------------------------------------------- + out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c], + bitdepth,maxbin,midbin,s0,s1) + # kernel ------------------------------------------- + + # release memory allocated by malloc + + free(se_e_r) + free(se_e_c) + free(se_w_r) + free(se_w_c) + free(se_n_r) + free(se_n_c) + free(se_s_r) + free(se_s_c) + + free(histo) + + return out diff --git a/skimage/rank/_core16p.pxd b/skimage/rank/_core16p.pxd index 9e992b6b..d162540c 100644 --- a/skimage/rank/_core16p.pxd +++ b/skimage/rank/_core16p.pxd @@ -1,15 +1,8 @@ -#cython: cdivision=True -#cython: boundscheck=False -#cython: nonecheck=False -#cython: wraparound=False - -import numpy as np cimport numpy as np -from libc.stdlib cimport malloc, free # generic cdef functions -cdef inline int int_max(int a, int b): return a if a >= b else b -cdef inline int int_min(int a, int b): return a if a <= b else b +cdef inline int int_max(int a, int b) +cdef inline int int_min(int a, int b) #--------------------------------------------------------------------------- # 16 bit core kernel receives extra information about data inferior and superior percentiles @@ -20,263 +13,4 @@ np.ndarray[np.uint16_t, ndim=2] image, np.ndarray[np.uint8_t, ndim=2] selem, np.ndarray[np.uint8_t, ndim=2] mask, np.ndarray[np.uint16_t, ndim=2] out, -char shift_x, char shift_y,int bitdepth, float p0, float p1): - """ Main loop, this function computes the histogram for each image point - - data is uint16 - - result is uint16 casted - """ - - cdef int rows = image.shape[0] - cdef int cols = image.shape[1] - cdef int srows = selem.shape[0] - cdef int scols = selem.shape[1] - - cdef int centre_r = int(selem.shape[0] / 2) + shift_y - cdef int centre_c = int(selem.shape[1] / 2) + shift_x - - # check that structuring element center is inside the element bounding box - assert centre_r >= 0 - assert centre_c >= 0 - assert centre_r < srows - assert centre_c < scols - - assert bitdepth in range(2,13) - - maxbin_list = [0,0,4,8,16,32,64,128,256,512,1024,2048,4096] - midbin_list = [0,0,2,4,8,16,32,64,128,256,512,1024,2048] - - #set maxbin and midbin - cdef int maxbin=maxbin_list[bitdepth],midbin=midbin_list[bitdepth] - - assert (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 - - # number of element in each attack border - cdef int num_se_n, num_se_s, num_se_e, num_se_w - - # the current local histogram distribution - cdef int* histo = malloc(maxbin * sizeof(int)) - - # these lists contain the relative pixel row and column for each of the 4 attack borders - # east, west, north and south - # e.g. se_e_r lists the rows of the east structuring element border - - cdef int* se_e_r = 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 - - num_se_n = num_se_s = num_se_e = num_se_w = 0 - - for r in range(srows): - for c in range(scols): - if t_e[r,c]: - se_e_r[num_se_e] = r - centre_r - se_e_c[num_se_e] = c - centre_c - num_se_e += 1 - if t_w[r,c]: - se_w_r[num_se_w] = r - centre_r - se_w_c[num_se_w] = c - centre_c - num_se_w += 1 - if t_n[r,c]: - se_n_r[num_se_n] = r - centre_r - se_n_c[num_se_n] = c - centre_c - num_se_n += 1 - if t_s[r,c]: - se_s_r[num_se_s] = r - centre_r - se_s_c[num_se_s] = c - centre_c - num_se_s += 1 - - # initial population and histogram - for i in range(maxbin): - histo[i] = 0 - - pop = 0 - - for r in range(srows): - for c in range(scols): - rr = r - cc = c - if selem[r, c]: - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] += 1 - pop += 1. - - r = 0 - c = 0 - # kernel ------------------------------------------- - out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c], - bitdepth,maxbin,midbin,p0,p1) - # kernel ------------------------------------------- - - # main loop - r = 0 - for even_row in range(0,rows,2): - # ---> west to east - for c in range(1,cols): - for s in range(num_se_e): - rr = r + se_e_r[s] + centre_r - cc = c + se_e_c[s] + centre_c - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] += 1 - pop += 1. - for s in range(num_se_w): - rr = r + se_w_r[s] + centre_r - cc = c + se_w_c[s] + centre_c - 1 - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] -= 1 - pop -= 1. - - # kernel ------------------------------------------- - out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c], - bitdepth,maxbin,midbin,p0,p1) - # kernel ------------------------------------------- - - r += 1 # pass to the next row - if r>=rows: - break - - # ---> north to south - for s in range(num_se_s): - rr = r + se_s_r[s] + centre_r - cc = c + se_s_c[s] + centre_c - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] += 1 - pop += 1. - for s in range(num_se_n): - rr = r + se_n_r[s] + centre_r - 1 - cc = c + se_n_c[s] + centre_c - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] -= 1 - pop -= 1. - - # kernel ------------------------------------------- - out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c], - bitdepth,maxbin,midbin,p0,p1) - # kernel ------------------------------------------- - - # ---> east to west - for c in range(cols-2,-1,-1): - for s in range(num_se_w): - rr = r + se_w_r[s] + centre_r - cc = c + se_w_c[s] + centre_c - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] += 1 - pop += 1. - for s in range(num_se_e): - rr = r + se_e_r[s] + centre_r - cc = c + se_e_c[s] + centre_c + 1 - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] -= 1 - pop -= 1. - - # kernel ------------------------------------------- - out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c], - bitdepth,maxbin,midbin,p0,p1) - # kernel ------------------------------------------- - - r += 1 # pass to the next row - if r>=rows: - break - - # ---> north to south - for s in range(num_se_s): - rr = r + se_s_r[s] + centre_r - cc = c + se_s_c[s] + centre_c - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] += 1 - pop += 1. - for s in range(num_se_n): - rr = r + se_n_r[s] + centre_r - 1 - cc = c + se_n_c[s] + centre_c - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] -= 1 - pop -= 1. - - # kernel ------------------------------------------- - out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c], - bitdepth,maxbin,midbin,p0,p1) - # kernel ------------------------------------------- - - # release memory allocated by malloc - - free(se_e_r) - free(se_e_c) - free(se_w_r) - free(se_w_c) - free(se_n_r) - free(se_n_c) - free(se_s_r) - free(se_s_c) - - free(histo) - - return out - - +char shift_x, char shift_y,int bitdepth, float p0, float p1) \ No newline at end of file diff --git a/skimage/rank/_core16p.pyx b/skimage/rank/_core16p.pyx new file mode 100644 index 00000000..9e992b6b --- /dev/null +++ b/skimage/rank/_core16p.pyx @@ -0,0 +1,282 @@ +#cython: cdivision=True +#cython: boundscheck=False +#cython: nonecheck=False +#cython: wraparound=False + +import numpy as np +cimport numpy as np +from libc.stdlib cimport malloc, free + +# generic cdef functions +cdef inline int int_max(int a, int b): return a if a >= b else b +cdef inline int int_min(int a, int b): return a if a <= b else b + +#--------------------------------------------------------------------------- +# 16 bit core kernel receives extra information about data inferior and superior percentiles +#--------------------------------------------------------------------------- + +cdef inline _core16p(np.uint16_t kernel(int*, float, np.uint16_t,int,int,int, float, float), +np.ndarray[np.uint16_t, ndim=2] image, +np.ndarray[np.uint8_t, ndim=2] selem, +np.ndarray[np.uint8_t, ndim=2] mask, +np.ndarray[np.uint16_t, ndim=2] out, +char shift_x, char shift_y,int bitdepth, float p0, float p1): + """ Main loop, this function computes the histogram for each image point + - data is uint16 + - result is uint16 casted + """ + + cdef int rows = image.shape[0] + cdef int cols = image.shape[1] + cdef int srows = selem.shape[0] + cdef int scols = selem.shape[1] + + cdef int centre_r = int(selem.shape[0] / 2) + shift_y + cdef int centre_c = int(selem.shape[1] / 2) + shift_x + + # check that structuring element center is inside the element bounding box + assert centre_r >= 0 + assert centre_c >= 0 + assert centre_r < srows + assert centre_c < scols + + assert bitdepth in range(2,13) + + maxbin_list = [0,0,4,8,16,32,64,128,256,512,1024,2048,4096] + midbin_list = [0,0,2,4,8,16,32,64,128,256,512,1024,2048] + + #set maxbin and midbin + cdef int maxbin=maxbin_list[bitdepth],midbin=midbin_list[bitdepth] + + assert (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 + + # number of element in each attack border + cdef int num_se_n, num_se_s, num_se_e, num_se_w + + # the current local histogram distribution + cdef int* histo = malloc(maxbin * sizeof(int)) + + # these lists contain the relative pixel row and column for each of the 4 attack borders + # east, west, north and south + # e.g. se_e_r lists the rows of the east structuring element border + + cdef int* se_e_r = 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 + + num_se_n = num_se_s = num_se_e = num_se_w = 0 + + for r in range(srows): + for c in range(scols): + if t_e[r,c]: + se_e_r[num_se_e] = r - centre_r + se_e_c[num_se_e] = c - centre_c + num_se_e += 1 + if t_w[r,c]: + se_w_r[num_se_w] = r - centre_r + se_w_c[num_se_w] = c - centre_c + num_se_w += 1 + if t_n[r,c]: + se_n_r[num_se_n] = r - centre_r + se_n_c[num_se_n] = c - centre_c + num_se_n += 1 + if t_s[r,c]: + se_s_r[num_se_s] = r - centre_r + se_s_c[num_se_s] = c - centre_c + num_se_s += 1 + + # initial population and histogram + for i in range(maxbin): + histo[i] = 0 + + pop = 0 + + for r in range(srows): + for c in range(scols): + rr = r + cc = c + if selem[r, c]: + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] += 1 + pop += 1. + + r = 0 + c = 0 + # kernel ------------------------------------------- + out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c], + bitdepth,maxbin,midbin,p0,p1) + # kernel ------------------------------------------- + + # main loop + r = 0 + for even_row in range(0,rows,2): + # ---> west to east + for c in range(1,cols): + for s in range(num_se_e): + rr = r + se_e_r[s] + centre_r + cc = c + se_e_c[s] + centre_c + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] += 1 + pop += 1. + for s in range(num_se_w): + rr = r + se_w_r[s] + centre_r + cc = c + se_w_c[s] + centre_c - 1 + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] -= 1 + pop -= 1. + + # kernel ------------------------------------------- + out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c], + bitdepth,maxbin,midbin,p0,p1) + # kernel ------------------------------------------- + + r += 1 # pass to the next row + if r>=rows: + break + + # ---> north to south + for s in range(num_se_s): + rr = r + se_s_r[s] + centre_r + cc = c + se_s_c[s] + centre_c + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] += 1 + pop += 1. + for s in range(num_se_n): + rr = r + se_n_r[s] + centre_r - 1 + cc = c + se_n_c[s] + centre_c + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] -= 1 + pop -= 1. + + # kernel ------------------------------------------- + out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c], + bitdepth,maxbin,midbin,p0,p1) + # kernel ------------------------------------------- + + # ---> east to west + for c in range(cols-2,-1,-1): + for s in range(num_se_w): + rr = r + se_w_r[s] + centre_r + cc = c + se_w_c[s] + centre_c + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] += 1 + pop += 1. + for s in range(num_se_e): + rr = r + se_e_r[s] + centre_r + cc = c + se_e_c[s] + centre_c + 1 + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] -= 1 + pop -= 1. + + # kernel ------------------------------------------- + out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c], + bitdepth,maxbin,midbin,p0,p1) + # kernel ------------------------------------------- + + r += 1 # pass to the next row + if r>=rows: + break + + # ---> north to south + for s in range(num_se_s): + rr = r + se_s_r[s] + centre_r + cc = c + se_s_c[s] + centre_c + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] += 1 + pop += 1. + for s in range(num_se_n): + rr = r + se_n_r[s] + centre_r - 1 + cc = c + se_n_c[s] + centre_c + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] -= 1 + pop -= 1. + + # kernel ------------------------------------------- + out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c], + bitdepth,maxbin,midbin,p0,p1) + # kernel ------------------------------------------- + + # release memory allocated by malloc + + free(se_e_r) + free(se_e_c) + free(se_w_r) + free(se_w_c) + free(se_n_r) + free(se_n_c) + free(se_s_r) + free(se_s_c) + + free(histo) + + return out + + diff --git a/skimage/rank/_core8.pxd b/skimage/rank/_core8.pxd index d24297cb..aa3fe527 100644 --- a/skimage/rank/_core8.pxd +++ b/skimage/rank/_core8.pxd @@ -1,18 +1,4 @@ -""" to compile this use: ->>> python setup.py build_ext --inplace - -to generate html report use: ->>> cython -a core8.pxd -""" - -#cython: cdivision=True -#cython: boundscheck=False -#cython: nonecheck=False -#cython: wraparound=False - -import numpy as np cimport numpy as np -from libc.stdlib cimport malloc, free #--------------------------------------------------------------------------- # 8 bit core kernel @@ -23,247 +9,4 @@ np.ndarray[np.uint8_t, ndim=2] image, np.ndarray[np.uint8_t, ndim=2] selem, np.ndarray[np.uint8_t, ndim=2] mask, np.ndarray[np.uint8_t, ndim=2] out, -char shift_x, char shift_y): - """ Main loop, this function computes the histogram for each image point - - data is uint8 - - result is uint8 casted - """ - - cdef int rows = image.shape[0] - cdef int cols = image.shape[1] - cdef int srows = selem.shape[0] - cdef int scols = selem.shape[1] - - cdef int centre_r = int(selem.shape[0] / 2) + shift_y - cdef int centre_c = int(selem.shape[1] / 2) + shift_x - - # check that structuring element center is inside the element bounding box - assert centre_r >= 0 - assert centre_c >= 0 - assert centre_r < srows - assert centre_c < scols - - image = np.ascontiguousarray(image) - - if mask is None: - mask = np.ones((rows, cols), dtype=np.uint8) - else: - mask = np.ascontiguousarray(mask) - - if out is None: - out = np.zeros((rows, cols), dtype=np.uint8) - else: - out = np.ascontiguousarray(out) - - # create extended image and mask - cdef int erows = rows+srows-1 - cdef int ecols = cols+scols-1 - - cdef np.ndarray emask = np.zeros((erows, ecols), dtype=np.uint8) - cdef np.ndarray eimage = np.zeros((erows, ecols), dtype=np.uint8) - - eimage[centre_r:rows+centre_r,centre_c:cols+centre_c] = image - emask[centre_r:rows+centre_r,centre_c:cols+centre_c] = mask - - mask = np.ascontiguousarray(mask) - - # define pointers to the data - cdef np.uint8_t* eimage_data = 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 - - # number of element in each attack border - cdef int num_se_n, num_se_s, num_se_e, num_se_w - - # the current local histogram distribution - cdef int* histo = malloc(256 * sizeof(int)) - - # these lists contain the relative pixel row and column for each of the 4 attack borders - # east, west, north and south - # e.g. se_e_r lists the rows of the east structuring element border - - cdef int* se_e_r = 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 - - num_se_n = num_se_s = num_se_e = num_se_w = 0 - - for r in range(srows): - for c in range(scols): - if t_e[r,c]: - se_e_r[num_se_e] = r - centre_r - se_e_c[num_se_e] = c - centre_c - num_se_e += 1 - if t_w[r,c]: - se_w_r[num_se_w] = r - centre_r - se_w_c[num_se_w] = c - centre_c - num_se_w += 1 - if t_n[r,c]: - se_n_r[num_se_n] = r - centre_r - se_n_c[num_se_n] = c - centre_c - num_se_n += 1 - if t_s[r,c]: - se_s_r[num_se_s] = r - centre_r - se_s_c[num_se_s] = c - centre_c - num_se_s += 1 - - # initial population and histogram - for i in range(256): - histo[i] = 0 - - pop = 0 - - for r in range(srows): - for c in range(scols): - rr = r - cc = c - if selem[r, c]: - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] += 1 - pop += 1. - - r = 0 - c = 0 - # kernel ------------------------------------------- - out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c]) - # kernel ------------------------------------------- - - # main loop - r = 0 - for even_row in range(0,rows,2): - # ---> west to east - for c in range(1,cols): - for s in range(num_se_e): - rr = r + se_e_r[s] + centre_r - cc = c + se_e_c[s] + centre_c - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] += 1 - pop += 1. - for s in range(num_se_w): - rr = r + se_w_r[s] + centre_r - cc = c + se_w_c[s] + centre_c - 1 - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] -= 1 - pop -= 1. - - # kernel ------------------------------------------- - out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c]) - # kernel ------------------------------------------- - - r += 1 # pass to the next row - if r>=rows: - break - - # ---> north to south - for s in range(num_se_s): - rr = r + se_s_r[s] + centre_r - cc = c + se_s_c[s] + centre_c - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] += 1 - pop += 1. - for s in range(num_se_n): - rr = r + se_n_r[s] + centre_r - 1 - cc = c + se_n_c[s] + centre_c - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] -= 1 - pop -= 1. - - # kernel ------------------------------------------- - out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c]) - # kernel ------------------------------------------- - - # ---> east to west - for c in range(cols-2,-1,-1): - for s in range(num_se_w): - rr = r + se_w_r[s] + centre_r - cc = c + se_w_c[s] + centre_c - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] += 1 - pop += 1. - for s in range(num_se_e): - rr = r + se_e_r[s] + centre_r - cc = c + se_e_c[s] + centre_c + 1 - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] -= 1 - pop -= 1. - - # kernel ------------------------------------------- - out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c]) - # kernel ------------------------------------------- - - r += 1 # pass to the next row - if r>=rows: - break - - # ---> north to south - for s in range(num_se_s): - rr = r + se_s_r[s] + centre_r - cc = c + se_s_c[s] + centre_c - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] += 1 - pop += 1. - for s in range(num_se_n): - rr = r + se_n_r[s] + centre_r - 1 - cc = c + se_n_c[s] + centre_c - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] -= 1 - pop -= 1. - - # kernel ------------------------------------------- - out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c]) - # kernel ------------------------------------------- - - # release memory allocated by malloc - - free(se_e_r) - free(se_e_c) - free(se_w_r) - free(se_w_c) - free(se_n_r) - free(se_n_c) - free(se_s_r) - free(se_s_c) - - free(histo) - - return out - +char shift_x, char shift_y) diff --git a/skimage/rank/_core8.pyx b/skimage/rank/_core8.pyx new file mode 100644 index 00000000..d24297cb --- /dev/null +++ b/skimage/rank/_core8.pyx @@ -0,0 +1,269 @@ +""" to compile this use: +>>> python setup.py build_ext --inplace + +to generate html report use: +>>> cython -a core8.pxd +""" + +#cython: cdivision=True +#cython: boundscheck=False +#cython: nonecheck=False +#cython: wraparound=False + +import numpy as np +cimport numpy as np +from libc.stdlib cimport malloc, free + +#--------------------------------------------------------------------------- +# 8 bit core kernel +#--------------------------------------------------------------------------- + +cdef inline _core8(np.uint8_t kernel(int*, float, np.uint8_t), +np.ndarray[np.uint8_t, ndim=2] image, +np.ndarray[np.uint8_t, ndim=2] selem, +np.ndarray[np.uint8_t, ndim=2] mask, +np.ndarray[np.uint8_t, ndim=2] out, +char shift_x, char shift_y): + """ Main loop, this function computes the histogram for each image point + - data is uint8 + - result is uint8 casted + """ + + cdef int rows = image.shape[0] + cdef int cols = image.shape[1] + cdef int srows = selem.shape[0] + cdef int scols = selem.shape[1] + + cdef int centre_r = int(selem.shape[0] / 2) + shift_y + cdef int centre_c = int(selem.shape[1] / 2) + shift_x + + # check that structuring element center is inside the element bounding box + assert centre_r >= 0 + assert centre_c >= 0 + assert centre_r < srows + assert centre_c < scols + + image = np.ascontiguousarray(image) + + if mask is None: + mask = np.ones((rows, cols), dtype=np.uint8) + else: + mask = np.ascontiguousarray(mask) + + if out is None: + out = np.zeros((rows, cols), dtype=np.uint8) + else: + out = np.ascontiguousarray(out) + + # create extended image and mask + cdef int erows = rows+srows-1 + cdef int ecols = cols+scols-1 + + cdef np.ndarray emask = np.zeros((erows, ecols), dtype=np.uint8) + cdef np.ndarray eimage = np.zeros((erows, ecols), dtype=np.uint8) + + eimage[centre_r:rows+centre_r,centre_c:cols+centre_c] = image + emask[centre_r:rows+centre_r,centre_c:cols+centre_c] = mask + + mask = np.ascontiguousarray(mask) + + # define pointers to the data + cdef np.uint8_t* eimage_data = 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 + + # number of element in each attack border + cdef int num_se_n, num_se_s, num_se_e, num_se_w + + # the current local histogram distribution + cdef int* histo = malloc(256 * sizeof(int)) + + # these lists contain the relative pixel row and column for each of the 4 attack borders + # east, west, north and south + # e.g. se_e_r lists the rows of the east structuring element border + + cdef int* se_e_r = 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 + + num_se_n = num_se_s = num_se_e = num_se_w = 0 + + for r in range(srows): + for c in range(scols): + if t_e[r,c]: + se_e_r[num_se_e] = r - centre_r + se_e_c[num_se_e] = c - centre_c + num_se_e += 1 + if t_w[r,c]: + se_w_r[num_se_w] = r - centre_r + se_w_c[num_se_w] = c - centre_c + num_se_w += 1 + if t_n[r,c]: + se_n_r[num_se_n] = r - centre_r + se_n_c[num_se_n] = c - centre_c + num_se_n += 1 + if t_s[r,c]: + se_s_r[num_se_s] = r - centre_r + se_s_c[num_se_s] = c - centre_c + num_se_s += 1 + + # initial population and histogram + for i in range(256): + histo[i] = 0 + + pop = 0 + + for r in range(srows): + for c in range(scols): + rr = r + cc = c + if selem[r, c]: + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] += 1 + pop += 1. + + r = 0 + c = 0 + # kernel ------------------------------------------- + out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c]) + # kernel ------------------------------------------- + + # main loop + r = 0 + for even_row in range(0,rows,2): + # ---> west to east + for c in range(1,cols): + for s in range(num_se_e): + rr = r + se_e_r[s] + centre_r + cc = c + se_e_c[s] + centre_c + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] += 1 + pop += 1. + for s in range(num_se_w): + rr = r + se_w_r[s] + centre_r + cc = c + se_w_c[s] + centre_c - 1 + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] -= 1 + pop -= 1. + + # kernel ------------------------------------------- + out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c]) + # kernel ------------------------------------------- + + r += 1 # pass to the next row + if r>=rows: + break + + # ---> north to south + for s in range(num_se_s): + rr = r + se_s_r[s] + centre_r + cc = c + se_s_c[s] + centre_c + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] += 1 + pop += 1. + for s in range(num_se_n): + rr = r + se_n_r[s] + centre_r - 1 + cc = c + se_n_c[s] + centre_c + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] -= 1 + pop -= 1. + + # kernel ------------------------------------------- + out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c]) + # kernel ------------------------------------------- + + # ---> east to west + for c in range(cols-2,-1,-1): + for s in range(num_se_w): + rr = r + se_w_r[s] + centre_r + cc = c + se_w_c[s] + centre_c + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] += 1 + pop += 1. + for s in range(num_se_e): + rr = r + se_e_r[s] + centre_r + cc = c + se_e_c[s] + centre_c + 1 + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] -= 1 + pop -= 1. + + # kernel ------------------------------------------- + out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c]) + # kernel ------------------------------------------- + + r += 1 # pass to the next row + if r>=rows: + break + + # ---> north to south + for s in range(num_se_s): + rr = r + se_s_r[s] + centre_r + cc = c + se_s_c[s] + centre_c + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] += 1 + pop += 1. + for s in range(num_se_n): + rr = r + se_n_r[s] + centre_r - 1 + cc = c + se_n_c[s] + centre_c + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] -= 1 + pop -= 1. + + # kernel ------------------------------------------- + out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c]) + # kernel ------------------------------------------- + + # release memory allocated by malloc + + free(se_e_r) + free(se_e_c) + free(se_w_r) + free(se_w_c) + free(se_n_r) + free(se_n_c) + free(se_s_r) + free(se_s_c) + + free(histo) + + return out + diff --git a/skimage/rank/_core8p.pxd b/skimage/rank/_core8p.pxd index b1adac70..8d3181e8 100644 --- a/skimage/rank/_core8p.pxd +++ b/skimage/rank/_core8p.pxd @@ -1,15 +1,8 @@ -#cython: cdivision=True -#cython: boundscheck=False -#cython: nonecheck=False -#cython: wraparound=False - -import numpy as np cimport numpy as np -from libc.stdlib cimport malloc, free # generic cdef functions -cdef inline np.uint8_t uint8_max(np.uint8_t a, np.uint8_t b): return a if a >= b else b -cdef inline np.uint8_t uint8_min(np.uint8_t a, np.uint8_t b): return a if a <= b else b +cdef inline np.uint8_t uint8_max(np.uint8_t a, np.uint8_t b) +cdef inline np.uint8_t uint8_min(np.uint8_t a, np.uint8_t b) #--------------------------------------------------------------------------- # 8 bit core kernel receives extra information about data inferior and superior percentiles @@ -20,247 +13,5 @@ np.ndarray[np.uint8_t, ndim=2] image, np.ndarray[np.uint8_t, ndim=2] selem, np.ndarray[np.uint8_t, ndim=2] mask, np.ndarray[np.uint8_t, ndim=2] out, -char shift_x, char shift_y, float p0, float p1): - """ Main loop, this function computes the histogram for each image point - - data is uint8 - - result is uint8 casted - """ - - cdef int rows = image.shape[0] - cdef int cols = image.shape[1] - cdef int srows = selem.shape[0] - cdef int scols = selem.shape[1] - - cdef int centre_r = int(selem.shape[0] / 2) + shift_y - cdef int centre_c = int(selem.shape[1] / 2) + shift_x - - # check that structuring element center is inside the element bounding box - assert centre_r >= 0 - assert centre_c >= 0 - assert centre_r < srows - assert centre_c < scols - - image = np.ascontiguousarray(image) - - if mask is None: - mask = np.ones((rows, cols), dtype=np.uint8) - else: - mask = np.ascontiguousarray(mask) - - if out is None: - out = np.zeros((rows, cols), dtype=np.uint8) - else: - out = np.ascontiguousarray(out) - - # create extended image and mask - cdef int erows = rows+srows-1 - cdef int ecols = cols+scols-1 - - cdef np.ndarray emask = np.zeros((erows, ecols), dtype=np.uint8) - cdef np.ndarray eimage = np.zeros((erows, ecols), dtype=np.uint8) - - eimage[centre_r:rows+centre_r,centre_c:cols+centre_c] = image - emask[centre_r:rows+centre_r,centre_c:cols+centre_c] = mask - - mask = np.ascontiguousarray(mask) - - # define pointers to the data - cdef np.uint8_t* eimage_data = 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 - - # number of element in each attack border - cdef int num_se_n, num_se_s, num_se_e, num_se_w - - # the current local histogram distribution - cdef int* histo = malloc(256 * sizeof(int)) - - # these lists contain the relative pixel row and column for each of the 4 attack borders - # east, west, north and south - # e.g. se_e_r lists the rows of the east structuring element border - - cdef int* se_e_r = 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 - - num_se_n = num_se_s = num_se_e = num_se_w = 0 - - for r in range(srows): - for c in range(scols): - if t_e[r,c]: - se_e_r[num_se_e] = r - centre_r - se_e_c[num_se_e] = c - centre_c - num_se_e += 1 - if t_w[r,c]: - se_w_r[num_se_w] = r - centre_r - se_w_c[num_se_w] = c - centre_c - num_se_w += 1 - if t_n[r,c]: - se_n_r[num_se_n] = r - centre_r - se_n_c[num_se_n] = c - centre_c - num_se_n += 1 - if t_s[r,c]: - se_s_r[num_se_s] = r - centre_r - se_s_c[num_se_s] = c - centre_c - num_se_s += 1 - - # initial population and histogram - for i in range(256): - histo[i] = 0 - - pop = 0 - - for r in range(srows): - for c in range(scols): - rr = r - cc = c - if selem[r, c]: - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] += 1 - pop += 1. - - r = 0 - c = 0 - # kernel ------------------------------------------- - out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],p0,p1) - # kernel ------------------------------------------- - - # main loop - r = 0 - for even_row in range(0,rows,2): - # ---> west to east - for c in range(1,cols): - for s in range(num_se_e): - rr = r + se_e_r[s] + centre_r - cc = c + se_e_c[s] + centre_c - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] += 1 - pop += 1. - for s in range(num_se_w): - rr = r + se_w_r[s] + centre_r - cc = c + se_w_c[s] + centre_c - 1 - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] -= 1 - pop -= 1. - - # kernel ------------------------------------------- - out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],p0,p1) - # kernel ------------------------------------------- - - r += 1 # pass to the next row - if r>=rows: - break - - # ---> north to south - for s in range(num_se_s): - rr = r + se_s_r[s] + centre_r - cc = c + se_s_c[s] + centre_c - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] += 1 - pop += 1. - for s in range(num_se_n): - rr = r + se_n_r[s] + centre_r - 1 - cc = c + se_n_c[s] + centre_c - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] -= 1 - pop -= 1. - - # kernel ------------------------------------------- - out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],p0,p1) - # kernel ------------------------------------------- - - # ---> east to west - for c in range(cols-2,-1,-1): - for s in range(num_se_w): - rr = r + se_w_r[s] + centre_r - cc = c + se_w_c[s] + centre_c - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] += 1 - pop += 1. - for s in range(num_se_e): - rr = r + se_e_r[s] + centre_r - cc = c + se_e_c[s] + centre_c + 1 - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] -= 1 - pop -= 1. - - # kernel ------------------------------------------- - out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],p0,p1) - # kernel ------------------------------------------- - - r += 1 # pass to the next row - if r>=rows: - break - - # ---> north to south - for s in range(num_se_s): - rr = r + se_s_r[s] + centre_r - cc = c + se_s_c[s] + centre_c - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] += 1 - pop += 1. - for s in range(num_se_n): - rr = r + se_n_r[s] + centre_r - 1 - cc = c + se_n_c[s] + centre_c - if emask_data[rr * ecols + cc]: - value = eimage_data[rr * ecols + cc] - histo[value] -= 1 - pop -= 1. - - # kernel ------------------------------------------- - out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],p0,p1) - # kernel ------------------------------------------- - - # release memory allocated by malloc - - free(se_e_r) - free(se_e_c) - free(se_w_r) - free(se_w_c) - free(se_n_r) - free(se_n_c) - free(se_s_r) - free(se_s_c) - - free(histo) - - return out +char shift_x, char shift_y, float p0, float p1) diff --git a/skimage/rank/_core8p.pyx b/skimage/rank/_core8p.pyx new file mode 100644 index 00000000..b1adac70 --- /dev/null +++ b/skimage/rank/_core8p.pyx @@ -0,0 +1,266 @@ +#cython: cdivision=True +#cython: boundscheck=False +#cython: nonecheck=False +#cython: wraparound=False + +import numpy as np +cimport numpy as np +from libc.stdlib cimport malloc, free + +# generic cdef functions +cdef inline np.uint8_t uint8_max(np.uint8_t a, np.uint8_t b): return a if a >= b else b +cdef inline np.uint8_t uint8_min(np.uint8_t a, np.uint8_t b): return a if a <= b else b + +#--------------------------------------------------------------------------- +# 8 bit core kernel receives extra information about data inferior and superior percentiles +#--------------------------------------------------------------------------- + +cdef inline _core8p(np.uint8_t kernel(int*, float, np.uint8_t, float, float), +np.ndarray[np.uint8_t, ndim=2] image, +np.ndarray[np.uint8_t, ndim=2] selem, +np.ndarray[np.uint8_t, ndim=2] mask, +np.ndarray[np.uint8_t, ndim=2] out, +char shift_x, char shift_y, float p0, float p1): + """ Main loop, this function computes the histogram for each image point + - data is uint8 + - result is uint8 casted + """ + + cdef int rows = image.shape[0] + cdef int cols = image.shape[1] + cdef int srows = selem.shape[0] + cdef int scols = selem.shape[1] + + cdef int centre_r = int(selem.shape[0] / 2) + shift_y + cdef int centre_c = int(selem.shape[1] / 2) + shift_x + + # check that structuring element center is inside the element bounding box + assert centre_r >= 0 + assert centre_c >= 0 + assert centre_r < srows + assert centre_c < scols + + image = np.ascontiguousarray(image) + + if mask is None: + mask = np.ones((rows, cols), dtype=np.uint8) + else: + mask = np.ascontiguousarray(mask) + + if out is None: + out = np.zeros((rows, cols), dtype=np.uint8) + else: + out = np.ascontiguousarray(out) + + # create extended image and mask + cdef int erows = rows+srows-1 + cdef int ecols = cols+scols-1 + + cdef np.ndarray emask = np.zeros((erows, ecols), dtype=np.uint8) + cdef np.ndarray eimage = np.zeros((erows, ecols), dtype=np.uint8) + + eimage[centre_r:rows+centre_r,centre_c:cols+centre_c] = image + emask[centre_r:rows+centre_r,centre_c:cols+centre_c] = mask + + mask = np.ascontiguousarray(mask) + + # define pointers to the data + cdef np.uint8_t* eimage_data = 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 + + # number of element in each attack border + cdef int num_se_n, num_se_s, num_se_e, num_se_w + + # the current local histogram distribution + cdef int* histo = malloc(256 * sizeof(int)) + + # these lists contain the relative pixel row and column for each of the 4 attack borders + # east, west, north and south + # e.g. se_e_r lists the rows of the east structuring element border + + cdef int* se_e_r = 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 + + num_se_n = num_se_s = num_se_e = num_se_w = 0 + + for r in range(srows): + for c in range(scols): + if t_e[r,c]: + se_e_r[num_se_e] = r - centre_r + se_e_c[num_se_e] = c - centre_c + num_se_e += 1 + if t_w[r,c]: + se_w_r[num_se_w] = r - centre_r + se_w_c[num_se_w] = c - centre_c + num_se_w += 1 + if t_n[r,c]: + se_n_r[num_se_n] = r - centre_r + se_n_c[num_se_n] = c - centre_c + num_se_n += 1 + if t_s[r,c]: + se_s_r[num_se_s] = r - centre_r + se_s_c[num_se_s] = c - centre_c + num_se_s += 1 + + # initial population and histogram + for i in range(256): + histo[i] = 0 + + pop = 0 + + for r in range(srows): + for c in range(scols): + rr = r + cc = c + if selem[r, c]: + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] += 1 + pop += 1. + + r = 0 + c = 0 + # kernel ------------------------------------------- + out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],p0,p1) + # kernel ------------------------------------------- + + # main loop + r = 0 + for even_row in range(0,rows,2): + # ---> west to east + for c in range(1,cols): + for s in range(num_se_e): + rr = r + se_e_r[s] + centre_r + cc = c + se_e_c[s] + centre_c + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] += 1 + pop += 1. + for s in range(num_se_w): + rr = r + se_w_r[s] + centre_r + cc = c + se_w_c[s] + centre_c - 1 + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] -= 1 + pop -= 1. + + # kernel ------------------------------------------- + out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],p0,p1) + # kernel ------------------------------------------- + + r += 1 # pass to the next row + if r>=rows: + break + + # ---> north to south + for s in range(num_se_s): + rr = r + se_s_r[s] + centre_r + cc = c + se_s_c[s] + centre_c + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] += 1 + pop += 1. + for s in range(num_se_n): + rr = r + se_n_r[s] + centre_r - 1 + cc = c + se_n_c[s] + centre_c + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] -= 1 + pop -= 1. + + # kernel ------------------------------------------- + out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],p0,p1) + # kernel ------------------------------------------- + + # ---> east to west + for c in range(cols-2,-1,-1): + for s in range(num_se_w): + rr = r + se_w_r[s] + centre_r + cc = c + se_w_c[s] + centre_c + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] += 1 + pop += 1. + for s in range(num_se_e): + rr = r + se_e_r[s] + centre_r + cc = c + se_e_c[s] + centre_c + 1 + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] -= 1 + pop -= 1. + + # kernel ------------------------------------------- + out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],p0,p1) + # kernel ------------------------------------------- + + r += 1 # pass to the next row + if r>=rows: + break + + # ---> north to south + for s in range(num_se_s): + rr = r + se_s_r[s] + centre_r + cc = c + se_s_c[s] + centre_c + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] += 1 + pop += 1. + for s in range(num_se_n): + rr = r + se_n_r[s] + centre_r - 1 + cc = c + se_n_c[s] + centre_c + if emask_data[rr * ecols + cc]: + value = eimage_data[rr * ecols + cc] + histo[value] -= 1 + pop -= 1. + + # kernel ------------------------------------------- + out_data[r * cols + c] = kernel(histo,pop,eimage_data[(r+centre_r) * ecols + c + centre_c],p0,p1) + # kernel ------------------------------------------- + + # release memory allocated by malloc + + free(se_e_r) + free(se_e_c) + free(se_w_r) + free(se_w_c) + free(se_n_r) + free(se_n_c) + free(se_s_r) + free(se_s_c) + + free(histo) + + return out + diff --git a/skimage/rank/setup.py b/skimage/rank/setup.py index efe23515..20a59979 100644 --- a/skimage/rank/setup.py +++ b/skimage/rank/setup.py @@ -5,6 +5,8 @@ from skimage._build import cython base_path = os.path.abspath(os.path.dirname(__file__)) +import sys +sys.path.append('.') def configuration(parent_package='', top_path=None): from numpy.distutils.misc_util import Configuration, get_numpy_include_dirs @@ -12,12 +14,28 @@ def configuration(parent_package='', top_path=None): config = Configuration('rank', parent_package, top_path) # config.add_data_dir('tests') + + cython(['_core8.pyx'], working_path=base_path) + cython(['_core8p.pyx'], working_path=base_path) + cython(['_core16.pyx'], working_path=base_path) + cython(['_core16p.pyx'], working_path=base_path) + cython(['_core16b.pyx'], working_path=base_path) cython(['_crank8.pyx'], working_path=base_path) cython(['_crank8_percentiles.pyx'], working_path=base_path) cython(['_crank16.pyx'], working_path=base_path) cython(['_crank16_percentiles.pyx'], working_path=base_path) cython(['_crank16_bilateral.pyx'], working_path=base_path) + config.add_extension('_core8', sources=['_core8.c'], + include_dirs=[get_numpy_include_dirs()]) + config.add_extension('_core8p', sources=['_core8p.c'], + include_dirs=[get_numpy_include_dirs()]) + config.add_extension('_core16', sources=['_core16.c'], + include_dirs=[get_numpy_include_dirs()]) + config.add_extension('_core16p', sources=['_core16p.c'], + include_dirs=[get_numpy_include_dirs()]) + config.add_extension('_core16b', sources=['_core16b.c'], + include_dirs=[get_numpy_include_dirs()]) config.add_extension('_crank8', sources=['_crank8.c'], include_dirs=[get_numpy_include_dirs()]) config.add_extension('_crank8_percentiles', sources=['_crank8_percentiles.c'],