diff --git a/skimage/rank/_core16.pyx b/skimage/rank/_core16.pyx index 2e55aa8b..edef264b 100644 --- a/skimage/rank/_core16.pyx +++ b/skimage/rank/_core16.pyx @@ -22,6 +22,14 @@ from libc.stdlib cimport malloc, free 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 void histogram_increment(Py_ssize_t* histo,float *pop,np.uint16_t value): + histo[value] += 1 + pop[0] += 1. + +cdef inline void histogram_decrement(Py_ssize_t* histo,float *pop,np.uint16_t value): + histo[value] -= 1 + pop[0] -= 1. + cdef inline np.uint8_t is_in_mask(Py_ssize_t rows, Py_ssize_t cols,Py_ssize_t r, Py_ssize_t c,np.uint8_t* mask): """ returns 1 if given(r,c) coordinate are within the image frame ([0-rows],[0-cols]) and inside the given mask @@ -79,6 +87,9 @@ cdef inline _core16(np.uint16_t kernel(Py_ssize_t*, float, np.uint16_t,Py_ssize_ else: mask = np.ascontiguousarray(mask) + if image is out: + raise NotImplementedError("Cannot perform rank operation in place.") + if out is None: out = np.zeros((rows, cols), dtype=np.uint16) else: @@ -165,9 +176,7 @@ cdef inline _core16(np.uint16_t kernel(Py_ssize_t*, float, np.uint16_t,Py_ssize_ cc = c - centre_c if selem[r, c]: if is_in_mask(rows,cols,rr,cc,mask_data): - value = image_data[rr * cols + cc] - histo[value] += 1 - pop += 1. + histogram_increment(histo,&pop,image_data[rr * cols + cc]) r = 0 c = 0 @@ -185,16 +194,13 @@ cdef inline _core16(np.uint16_t kernel(Py_ssize_t*, float, np.uint16_t,Py_ssize_ rr = r + se_e_r[s] cc = c + se_e_c[s] if is_in_mask(rows,cols,rr,cc,mask_data): - value = image_data[rr * cols + cc] - histo[value] += 1 - pop += 1. + histogram_increment(histo,&pop,image_data[rr * cols + cc]) + for s in range(num_se_w): rr = r + se_w_r[s] cc = c + se_w_c[s] - 1 if is_in_mask(rows,cols,rr,cc,mask_data): - value = image_data[rr * cols + cc] - histo[value] -= 1 - pop -= 1. + histogram_decrement(histo,&pop,image_data[rr * cols + cc]) # kernel ------------------------------------------- out_data[r * cols + c] = kernel(histo,pop,image_data[r * cols + c ], @@ -210,16 +216,13 @@ cdef inline _core16(np.uint16_t kernel(Py_ssize_t*, float, np.uint16_t,Py_ssize_ rr = r + se_s_r[s] cc = c + se_s_c[s] if is_in_mask(rows,cols,rr,cc,mask_data): - value = image_data[rr * cols + cc] - histo[value] += 1 - pop += 1. + histogram_increment(histo,&pop,image_data[rr * cols + cc]) + for s in range(num_se_n): rr = r + se_n_r[s] - 1 cc = c + se_n_c[s] if is_in_mask(rows,cols,rr,cc,mask_data): - value = image_data[rr * cols + cc] - histo[value] -= 1 - pop -= 1. + histogram_decrement(histo,&pop,image_data[rr * cols + cc]) # kernel ------------------------------------------- out_data[r * cols + c] = kernel(histo,pop,image_data[r * cols + c], @@ -232,16 +235,13 @@ cdef inline _core16(np.uint16_t kernel(Py_ssize_t*, float, np.uint16_t,Py_ssize_ rr = r + se_w_r[s] cc = c + se_w_c[s] if is_in_mask(rows,cols,rr,cc,mask_data): - value = image_data[rr * cols + cc] - histo[value] += 1 - pop += 1. + histogram_increment(histo,&pop,image_data[rr * cols + cc]) + for s in range(num_se_e): rr = r + se_e_r[s] cc = c + se_e_c[s] + 1 if is_in_mask(rows,cols,rr,cc,mask_data): - value = image_data[rr * cols + cc] - histo[value] -= 1 - pop -= 1. + histogram_decrement(histo,&pop,image_data[rr * cols + cc]) # kernel ------------------------------------------- out_data[r * cols + c] = kernel(histo,pop,image_data[r * cols + c ], @@ -257,16 +257,13 @@ cdef inline _core16(np.uint16_t kernel(Py_ssize_t*, float, np.uint16_t,Py_ssize_ rr = r + se_s_r[s] cc = c + se_s_c[s] if is_in_mask(rows,cols,rr,cc,mask_data): - value = image_data[rr * cols + cc] - histo[value] += 1 - pop += 1. + histogram_increment(histo,&pop,image_data[rr * cols + cc]) + for s in range(num_se_n): rr = r + se_n_r[s] - 1 cc = c + se_n_c[s] if is_in_mask(rows,cols,rr,cc,mask_data): - value = image_data[rr * cols + cc] - histo[value] -= 1 - pop -= 1. + histogram_decrement(histo,&pop,image_data[rr * cols + cc]) # kernel ------------------------------------------- out_data[r * cols + c] = kernel(histo,pop,image_data[r * cols + c ], diff --git a/skimage/rank/_core8.pyx b/skimage/rank/_core8.pyx index 19e09aec..86c40a6d 100644 --- a/skimage/rank/_core8.pyx +++ b/skimage/rank/_core8.pyx @@ -76,6 +76,9 @@ cdef inline _core8(np.uint8_t kernel(Py_ssize_t*, float, np.uint8_t, float, floa else: mask = np.ascontiguousarray(mask) + if image is out: + raise NotImplementedError("Cannot perform rank operation in place.") + if out is None: out = np.zeros((rows, cols), dtype=np.uint8) else: @@ -183,7 +186,7 @@ cdef inline _core8(np.uint8_t kernel(Py_ssize_t*, float, np.uint8_t, float, floa cc = c + se_e_c[s] if is_in_mask(rows,cols,rr,cc,mask_data): histogram_increment(histo,&pop,image_data[rr * cols + cc]) - + for s in range(num_se_w): rr = r + se_w_r[s] cc = c + se_w_c[s] - 1 diff --git a/skimage/rank/tests/test_suite.py b/skimage/rank/tests/test_suite.py index 9ed39956..d0e2c319 100644 --- a/skimage/rank/tests/test_suite.py +++ b/skimage/rank/tests/test_suite.py @@ -16,6 +16,7 @@ class TestSequenceFunctions(unittest.TestCase): def test_random_sizes(self): # make sure the size is not a problem + niter = 10 elem = np.asarray([[1,1,1],[1,1,1],[1,1,1]],dtype='uint8') for m,n in np.random.random_integers(1,100,size=(10,2)): @@ -39,8 +40,9 @@ class TestSequenceFunctions(unittest.TestCase): r = _crank16_percentiles.mean(image=a16,selem = elem,shift_x=+1,shift_y=+1,p0=.1,p1=.9) self.assertTrue(a16.shape == r.shape) - def test_compare_with_cmorph(self): + def test_compare_with_cmorph_dilate(self): #compare the result of maximum filter with dilate + a = (np.random.random((500,500))*256).astype('uint8') for r in range(1,20,1): @@ -50,7 +52,21 @@ class TestSequenceFunctions(unittest.TestCase): cm = cmorph.dilate(image=a,selem = elem) self.assertTrue((rc==cm).all()) + def test_compare_with_cmorph_erode(self): + #compare the result of maximum filter with erode + + a = (np.random.random((500,500))*256).astype('uint8') + + for r in range(1,20,1): + elem = np.ones((r,r),dtype='uint8') + # elem = (np.random.random((r,r))>.5).astype('uint8') + rc = _crank8.minimum(image=a,selem = elem) + cm = cmorph.erode(image=a,selem = elem) + self.assertTrue((rc==cm).all()) + def test_bitdepth(self): + # test the different bit depth for rank16 + elem = np.ones((3,3),dtype='uint8') a16 = np.ones((100,100),dtype='uint16')*255 r = _crank16_percentiles.mean(image=a16,selem = elem,shift_x=0,shift_y=0,p0=.1,p1=.9,bitdepth=8) @@ -64,6 +80,8 @@ class TestSequenceFunctions(unittest.TestCase): r = _crank16_percentiles.mean(image=a16,selem = elem,shift_x=0,shift_y=0,p0=.1,p1=.9,bitdepth=12) def test_population(self): + # check the number of valid pixels in the neighborhood + a = np.zeros((5,5),dtype='uint8') elem = np.ones((3,3),dtype='uint8') p = _crank8.pop(image=a,selem = elem) @@ -75,6 +93,8 @@ class TestSequenceFunctions(unittest.TestCase): np.testing.assert_array_equal(r,p) def test_structuring_element(self): + # check the output for a custom structuring element + a = np.zeros((6,6),dtype='uint8') a[2,2] = 255 elem = np.asarray([[1,1,0],[1,1,1],[0,0,1]],dtype='uint8') @@ -91,11 +111,37 @@ class TestSequenceFunctions(unittest.TestCase): @unittest.expectedFailure def test_fail_on_bitdepth(self): # should fail because data bitdepth is too high for the function + a16 = np.ones((100,100),dtype='uint16')*255 elem = np.ones((3,3),dtype='uint8') f = _crank16_percentiles.mean(image=a16,selem = elem,shift_x=0,shift_y=0,p0=.1,p1=.9,bitdepth=4) + def test_output(self): + #check rank function with external OUT output array + + selem = disk(20) + a = (np.random.random((500,500))*256).astype('uint8') + out = np.zeros_like(a) + f1 = rank.mean(a,selem,out=out) + f2 = rank.mean(a,selem) + np.testing.assert_array_equal(f1,f2) + np.testing.assert_array_equal(out,f2) + + @unittest.expectedFailure + def test_inplace_output(self): + #rank filters are not supposed to filter inplace + + selem = disk(20) + a = (np.random.random((500,500))*256).astype('uint8') + out = a + f = rank.mean(a,selem,out=out) + np.testing.assert_array_equal(f,out) + + def test_compare_autolevels(self): + # compare autolevel and percentile autolevel with p0=0.0 and p1=1.0 + # should returns the same arrays + image = data.camera() selem = disk(20) @@ -105,6 +151,9 @@ class TestSequenceFunctions(unittest.TestCase): assert (loc_autolevel==loc_perc_autolevel).all() def test_compare_autolevels_16bit(self): + # compare autolevel(16bit) and percentile autolevel(16bit) with p0=0.0 and p1=1.0 + # should returns the same arrays + image = data.camera().astype(np.uint16)*4 selem = disk(20) @@ -114,7 +163,9 @@ class TestSequenceFunctions(unittest.TestCase): assert (loc_autolevel==loc_perc_autolevel).all() def test_compare_8bit_vs_16bit(self): - # filters applied on 8bit image ore 16bit image (having only real 8bit of dynamic) should be identical + # filters applied on 8bit image ore 16bit image (having only real 8bit of dynamic) + # should be identical + i8 = data.camera() i16 = i8.astype(np.uint16) assert (i8==i16).all()