diff --git a/.gitignore b/.gitignore index 1bd4352c..2623ea74 100644 --- a/.gitignore +++ b/.gitignore @@ -22,4 +22,5 @@ doc/source/auto_examples/images/plot_*.png doc/source/auto_examples/images/thumb doc/source/auto_examples/applications/ doc/source/_static/random.js - +.idea/ +*.log diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index c73b14e2..8ba62e0b 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -117,3 +117,6 @@ - Petter Strandmark Perimeter calculation in regionprops. + +- Olivier Debeir + Rank filters (8- and 16-bits) using sliding window. \ No newline at end of file diff --git a/doc/__init__.py b/doc/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/doc/examples/applications/plot_rank_filters.py b/doc/examples/applications/plot_rank_filters.py new file mode 100644 index 00000000..c182671d --- /dev/null +++ b/doc/examples/applications/plot_rank_filters.py @@ -0,0 +1,719 @@ +""" +============ +Rank filters +============ + +Rank filters are non-linear filters using the local greylevels ordering to +compute the filtered value. This ensemble of filters share a common base: the +local grey-level histogram extraction computed on the neighborhood of a pixel +(defined by a 2D structuring element). If the filtered value is taken as the +middle value of the histogram, we get the classical median filter. + +Rank filters can be used for several purposes such as: + +* image quality enhancement + e.g. image smoothing, sharpening + +* image pre-processing + e.g. noise reduction, contrast enhancement + +* feature extraction + e.g. border detection, isolated point detection + +* post-processing + e.g. small object removal, object grouping, contour smoothing + +Some well known filters are specific cases of rank filters [1]_ e.g. +morphological dilation, morphological erosion, median filters. + +The different implementation availables in `skimage` are compared. + +In this example, we will see how to filter a greylevel image using some of the +linear and non-linear filters availables in skimage. We use the `camera` +image from `skimage.data`. + +.. [1] Pierre Soille, On morphological operators based on rank filters, Pattern + Recognition 35 (2002) 527-535. + +""" + +import numpy as np +import matplotlib.pyplot as plt + +from skimage import data + +ima = data.camera() +hist = np.histogram(ima, bins=np.arange(0, 256)) + +plt.figure(figsize=(8, 3)) +plt.subplot(1, 2, 1) +plt.imshow(ima, cmap=plt.cm.gray, interpolation='nearest') +plt.axis('off') +plt.subplot(1, 2, 2) +plt.plot(hist[1][:-1], hist[0], lw=2) +plt.title('histogram of grey values') + +""" + +.. image:: PLOT2RST.current_figure + +Noise removal +============= + +Some noise is added to the image, 1% of pixels are randomly set to 255, 1% are +randomly set to 0. The **median** filter is applied to remove the noise. + +.. note:: + + there are different implementations of median filter : + `skimage.filter.median_filter` and `skimage.filter.rank.median` + +""" + +noise = np.random.random(ima.shape) +nima = data.camera() +nima[noise > 0.99] = 255 +nima[noise < 0.01] = 0 + +from skimage.filter.rank import median +from skimage.morphology import disk + +fig = plt.figure(figsize=[10, 7]) + +lo = median(nima, disk(1)) +hi = median(nima, disk(5)) +ext = median(nima, disk(20)) +plt.subplot(2, 2, 1) +plt.imshow(nima, cmap=plt.cm.gray, vmin=0, vmax=255) +plt.xlabel('noised image') +plt.subplot(2, 2, 2) +plt.imshow(lo, cmap=plt.cm.gray, vmin=0, vmax=255) +plt.xlabel('median $r=1$') +plt.subplot(2, 2, 3) +plt.imshow(hi, cmap=plt.cm.gray, vmin=0, vmax=255) +plt.xlabel('median $r=5$') +plt.subplot(2, 2, 4) +plt.imshow(ext, cmap=plt.cm.gray, vmin=0, vmax=255) +plt.xlabel('median $r=20$') + +""" + +.. image:: PLOT2RST.current_figure + +The added noise is efficiently removed, as the image defaults are small (1 pixel +wide), a small filter radius is sufficient. As the radius is increasing, objects +with a bigger size are filtered as well, such as the camera tripod. The median +filter is commonly used for noise removal because borders are preserved. + +Image smoothing +================ + +The example hereunder shows how a local **mean** smoothes the camera man image. + +""" + +from skimage.filter.rank import mean + +fig = plt.figure(figsize=[10, 7]) + +loc_mean = mean(nima, disk(10)) +plt.subplot(1, 2, 1) +plt.imshow(ima, cmap=plt.cm.gray, vmin=0, vmax=255) +plt.xlabel('original') +plt.subplot(1, 2, 2) +plt.imshow(loc_mean, cmap=plt.cm.gray, vmin=0, vmax=255) +plt.xlabel('local mean $r=10$') + +""" + +.. image:: PLOT2RST.current_figure + +One may be interested in smoothing an image while preserving important borders +(median filters already achieved this), here we use the **bilateral** filter +that restricts the local neighborhood to pixel having a greylevel similar to +the central one. + +.. note:: + + a different implementation is available for color images in + `skimage.filter.denoise_bilateral`. + +""" + +from skimage.filter.rank import bilateral_mean + +ima = data.camera() +selem = disk(10) + +bilat = bilateral_mean(ima.astype(np.uint16), disk(20), s0=10, s1=10) + +# display results +fig = plt.figure(figsize=[10, 7]) +plt.subplot(2, 2, 1) +plt.imshow(ima, cmap=plt.cm.gray) +plt.xlabel('original') +plt.subplot(2, 2, 3) +plt.imshow(bilat, cmap=plt.cm.gray) +plt.xlabel('bilateral mean') +plt.subplot(2, 2, 2) +plt.imshow(ima[200:350, 350:450], cmap=plt.cm.gray) +plt.subplot(2, 2, 4) +plt.imshow(bilat[200:350, 350:450], cmap=plt.cm.gray) + +""" + +.. image:: PLOT2RST.current_figure + +One can see that the large continuous part of the image (e.g. sky) is smoothed +whereas other details are preserved. + + +Contrast enhancement +==================== + +We compare here how the global histogram equalization is applied locally. + +The equalized image [2]_ has a roughly linear cumulative distribution function +for each pixel neighborhood. The local version [3]_ of the histogram +equalization emphasizes every local greylevel variations. + +.. [2] http://en.wikipedia.org/wiki/Histogram_equalization +.. [3] http://en.wikipedia.org/wiki/Adaptive_histogram_equalization + +""" + +from skimage import exposure +from skimage.filter import rank + +ima = data.camera() +# equalize globally and locally +glob = exposure.equalize(ima) * 255 +loc = rank.equalize(ima, disk(20)) + +# extract histogram for each image +hist = np.histogram(ima, bins=np.arange(0, 256)) +glob_hist = np.histogram(glob, bins=np.arange(0, 256)) +loc_hist = np.histogram(loc, bins=np.arange(0, 256)) + +plt.figure(figsize=(10, 10)) +plt.subplot(321) +plt.imshow(ima, cmap=plt.cm.gray, interpolation='nearest') +plt.axis('off') +plt.subplot(322) +plt.plot(hist[1][:-1], hist[0], lw=2) +plt.title('histogram of grey values') +plt.subplot(323) +plt.imshow(glob, cmap=plt.cm.gray, interpolation='nearest') +plt.axis('off') +plt.subplot(324) +plt.plot(glob_hist[1][:-1], glob_hist[0], lw=2) +plt.title('histogram of grey values') +plt.subplot(325) +plt.imshow(loc, cmap=plt.cm.gray, interpolation='nearest') +plt.axis('off') +plt.subplot(326) +plt.plot(loc_hist[1][:-1], loc_hist[0], lw=2) +plt.title('histogram of grey values') + +""" + +.. image:: PLOT2RST.current_figure + +another way to maximize the number of greylevels used for an image is to apply +a local autoleveling, i.e. here a pixel greylevel is proportionally remapped +between local minimum and local maximum. + +The following example shows how local autolevel enhances the camara man picture. + +""" + +from skimage.filter.rank import autolevel + +ima = data.camera() +selem = disk(10) + +auto = autolevel(ima.astype(np.uint16), disk(20)) + +# display results +fig = plt.figure(figsize=[10, 7]) +plt.subplot(1, 2, 1) +plt.imshow(ima, cmap=plt.cm.gray) +plt.xlabel('original') +plt.subplot(1, 2, 2) +plt.imshow(auto, cmap=plt.cm.gray) +plt.xlabel('local autolevel') + +""" + +.. image:: PLOT2RST.current_figure + +This filter is very sensitive to local outlayers, see the little white spot in +the sky left part. This is due to a local maximum which is very high comparing +to the rest of the neighborhood. One can moderate this using the percentile +version of the autolevel filter which uses given percentiles (one inferior, +one superior) in place of local minimum and maximum. The example below +illustrates how the percentile parameters influence the local autolevel result. + +""" + +from skimage.filter.rank import percentile_autolevel + +image = data.camera() + +selem = disk(20) +loc_autolevel = autolevel(image, selem=selem) +loc_perc_autolevel0 = percentile_autolevel(image, selem=selem, p0=.00, p1=1.0) +loc_perc_autolevel1 = percentile_autolevel(image, selem=selem, p0=.01, p1=.99) +loc_perc_autolevel2 = percentile_autolevel(image, selem=selem, p0=.05, p1=.95) +loc_perc_autolevel3 = percentile_autolevel(image, selem=selem, p0=.1, p1=.9) + +fig, axes = plt.subplots(nrows=3, figsize=(7, 8)) +ax0, ax1, ax2 = axes +plt.gray() + +ax0.imshow(np.hstack((image, loc_autolevel))) +ax0.set_title('original / autolevel') + +ax1.imshow( + np.hstack((loc_perc_autolevel0, loc_perc_autolevel1)), vmin=0, vmax=255) +ax1.set_title('percentile autolevel 0%,1%') +ax2.imshow( + np.hstack((loc_perc_autolevel2, loc_perc_autolevel3)), vmin=0, vmax=255) +ax2.set_title('percentile autolevel 5% and 10%') + +for ax in axes: + ax.axis('off') + +""" + +.. image:: PLOT2RST.current_figure + +The morphological contrast enhancement filter replaces the central pixel by the +local maximum if the original pixel value is closest to local maximum, otherwise +by the minimum local. + +""" + +from skimage.filter.rank import morph_contr_enh + +ima = data.camera() + +enh = morph_contr_enh(ima, disk(5)) + +# display results +fig = plt.figure(figsize=[10, 7]) +plt.subplot(2, 2, 1) +plt.imshow(ima, cmap=plt.cm.gray) +plt.xlabel('original') +plt.subplot(2, 2, 3) +plt.imshow(enh, cmap=plt.cm.gray) +plt.xlabel('local morphlogical contrast enhancement') +plt.subplot(2, 2, 2) +plt.imshow(ima[200:350, 350:450], cmap=plt.cm.gray) +plt.subplot(2, 2, 4) +plt.imshow(enh[200:350, 350:450], cmap=plt.cm.gray) + +""" + +.. image:: PLOT2RST.current_figure + +The percentile version of the local morphological contrast enhancement uses +percentile *p0* and *p1* instead of the local minimum and maximum. + +""" + +from skimage.filter.rank import percentile_morph_contr_enh + +ima = data.camera() + +penh = percentile_morph_contr_enh(ima, disk(5), p0=.1, p1=.9) + +# display results +fig = plt.figure(figsize=[10, 7]) +plt.subplot(2, 2, 1) +plt.imshow(ima, cmap=plt.cm.gray) +plt.xlabel('original') +plt.subplot(2, 2, 3) +plt.imshow(penh, cmap=plt.cm.gray) +plt.xlabel('local percentile morphlogical\n contrast enhancement') +plt.subplot(2, 2, 2) +plt.imshow(ima[200:350, 350:450], cmap=plt.cm.gray) +plt.subplot(2, 2, 4) +plt.imshow(penh[200:350, 350:450], cmap=plt.cm.gray) + +""" + +.. image:: PLOT2RST.current_figure + +Image threshold +=============== + +The Otsu's threshold [1]_ method can be applied locally using the local +greylevel distribution. In the example below, for each pixel, an "optimal" +threshold is determined by maximizing the variance between two classes of pixels +of the local neighborhood defined by a structuring element. + +The example compares the local threshold with the global threshold +`skimage.filter.threshold_otsu`. + +.. note:: + + Local thresholding is much slower than global one. There exists a function + for global Otsu thresholding: `skimage.filter.threshold_otsu`. + +.. [1] http://en.wikipedia.org/wiki/Otsu's_method + +""" + +from skimage.filter.rank import otsu +from skimage.filter import threshold_otsu + +p8 = data.page() + +radius = 10 +selem = disk(radius) + +# t_loc_otsu is an image +t_loc_otsu = otsu(p8, selem) +loc_otsu = p8 >= t_loc_otsu + +# t_glob_otsu is a scalar +t_glob_otsu = threshold_otsu(p8) +glob_otsu = p8 >= t_glob_otsu + +plt.figure() +plt.subplot(2, 2, 1) +plt.imshow(p8, cmap=plt.cm.gray) +plt.xlabel('original') +plt.colorbar() +plt.subplot(2, 2, 2) +plt.imshow(t_loc_otsu, cmap=plt.cm.gray) +plt.xlabel('local Otsu ($radius=%d$)' % radius) +plt.colorbar() +plt.subplot(2, 2, 3) +plt.imshow(p8 >= t_loc_otsu, cmap=plt.cm.gray) +plt.xlabel('original>=local Otsu' % t_glob_otsu) +plt.subplot(2, 2, 4) +plt.imshow(glob_otsu, cmap=plt.cm.gray) +plt.xlabel('global Otsu ($t=%d$)' % t_glob_otsu) + +""" + +.. image:: PLOT2RST.current_figure + +The following example shows how local Otsu's threshold handles a global level +shift applied to a synthetic image . + +""" + +n = 100 +theta = np.linspace(0, 10 * np.pi, n) +x = np.sin(theta) +m = (np.tile(x, (n, 1)) * np.linspace(0.1, 1, n) * 128 + 128).astype(np.uint8) + +radius = 10 +t = rank.otsu(m, disk(radius)) +plt.figure() +plt.subplot(1, 2, 1) +plt.imshow(m) +plt.xlabel('original') +plt.subplot(1, 2, 2) +plt.imshow(m >= t, interpolation='nearest') +plt.xlabel('local Otsu ($radius=%d$)' % radius) + +""" + +.. image:: PLOT2RST.current_figure + +Image morphology +================ + +Local maximum and local minimum are the base operators for greylevel +morphology. + +.. note:: + + `skimage.dilate` and `skimage.erode` are equivalent filters (see below for + comparison). + +Here is an example of the classical morphological greylevel filters: opening, +closing and morphological gradient. + +""" + +from skimage.filter.rank import maximum, minimum, gradient + +ima = data.camera() + +closing = maximum(minimum(ima, disk(5)), disk(5)) +opening = minimum(maximum(ima, disk(5)), disk(5)) +grad = gradient(ima, disk(5)) + +# display results +fig = plt.figure(figsize=[10, 7]) +plt.subplot(2, 2, 1) +plt.imshow(ima, cmap=plt.cm.gray) +plt.xlabel('original') +plt.subplot(2, 2, 2) +plt.imshow(closing, cmap=plt.cm.gray) +plt.xlabel('greylevel closing') +plt.subplot(2, 2, 3) +plt.imshow(opening, cmap=plt.cm.gray) +plt.xlabel('greylevel opening') +plt.subplot(2, 2, 4) +plt.imshow(grad, cmap=plt.cm.gray) +plt.xlabel('morphological gradient') + +""" + +.. image:: PLOT2RST.current_figure + +Feature extraction +=================== + +Local histogram can be exploited to compute local entropy, which is related to +the local image complexity. Entropy is computed using base 2 logarithm i.e. the +filter returns the minimum number of bits needed to encode local greylevel +distribution. + +`skimage.rank.entropy` returns local entropy on a given structuring element. +The following example shows this filter applied on 8- and 16- bit images. + +.. note:: + + to better use the available image bit, the function returns 10x entropy for + 8-bit images and 1000x entropy for 16-bit images. + +""" + +from skimage import data +from skimage.filter.rank import entropy +from skimage.morphology import disk +import numpy as np +import matplotlib.pyplot as plt + +# defining a 8- and a 16-bit test images +a8 = data.camera() +a16 = data.camera().astype(np.uint16) * 4 + +ent8 = entropy(a8, disk(5)) # pixel value contain 10x the local entropy +ent16 = entropy(a16, disk(5)) # pixel value contain 1000x the local entropy + +# display results +plt.figure(figsize=(10, 10)) + +plt.subplot(2, 2, 1) +plt.imshow(a8, cmap=plt.cm.gray) +plt.xlabel('8-bit image') +plt.colorbar() + +plt.subplot(2, 2, 2) +plt.imshow(ent8, cmap=plt.cm.jet) +plt.xlabel('entropy*10') +plt.colorbar() + +plt.subplot(2, 2, 3) +plt.imshow(a16, cmap=plt.cm.gray) +plt.xlabel('16-bit image') +plt.colorbar() + +plt.subplot(2, 2, 4) +plt.imshow(ent16, cmap=plt.cm.jet) +plt.xlabel('entropy*1000') +plt.colorbar() + +""" + +.. image:: PLOT2RST.current_figure + +Implementation +================ + +The central part of the `skimage.rank` filters is build on a sliding window that +update local greylevel histogram. This approach limits the algorithm complexity +to O(n) where n is the number of image pixels. The complexity is also limited +with respect to the structuring element size. + +""" + +from time import time + +from scipy.ndimage.filters import percentile_filter +from skimage.morphology import dilation +from skimage.filter import median_filter +from skimage.filter.rank import median, maximum + + +def exec_and_timeit(func): + """Decorator that returns both function results and execution time.""" + def wrapper(*arg): + t1 = time() + res = func(*arg) + t2 = time() + ms = (t2 - t1) * 1000.0 + return (res, ms) + return wrapper + + +@exec_and_timeit +def cr_med(image, selem): + return median(image=image, selem=selem) + + +@exec_and_timeit +def cr_max(image, selem): + return maximum(image=image, selem=selem) + + +@exec_and_timeit +def cm_dil(image, selem): + return dilation(image=image, selem=selem) + + +@exec_and_timeit +def ctmf_med(image, radius): + return median_filter(image=image, radius=radius) + + +@exec_and_timeit +def ndi_med(image, n): + return percentile_filter(image, 50, size=n * 2 - 1) + +""" + +Comparison between + +* `rank.maximum` +* `cmorph.dilate` + +on increasing structuring element size + +""" + +a = data.camera() + +rec = [] +e_range = range(1, 20, 2) +for r in e_range: + elem = disk(r + 1) + rc, ms_rc = cr_max(a, elem) + rcm, ms_rcm = cm_dil(a, elem) + rec.append((ms_rc, ms_rcm)) + +rec = np.asarray(rec) + +plt.figure() +plt.title('increasing element size') +plt.ylabel('time (ms)') +plt.xlabel('element radius') +plt.plot(e_range, rec) +plt.legend(['crank.maximum', 'cmorph.dilate']) + +""" + +and increasing image size + +.. image:: PLOT2RST.current_figure + +""" + +r = 9 +elem = disk(r + 1) + +rec = [] +s_range = range(100, 1000, 100) +for s in s_range: + a = (np.random.random((s, s)) * 256).astype('uint8') + (rc, ms_rc) = cr_max(a, elem) + (rcm, ms_rcm) = cm_dil(a, elem) + rec.append((ms_rc, ms_rcm)) + +rec = np.asarray(rec) + +plt.figure() +plt.title('increasing image size') +plt.ylabel('time (ms)') +plt.xlabel('image size') +plt.plot(s_range, rec) +plt.legend(['crank.maximum', 'cmorph.dilate']) + + +""" + +.. image:: PLOT2RST.current_figure + +Comparison between: + +* `rank.median` +* `ctmf.median_filter` +* `ndimage.percentile` + +on increasing structuring element size + +""" + +a = data.camera() + +rec = [] +e_range = range(2, 30, 4) +for r in e_range: + elem = disk(r + 1) + rc, ms_rc = cr_med(a, elem) + rctmf, ms_rctmf = ctmf_med(a, r) + rndi, ms_ndi = ndi_med(a, r) + rec.append((ms_rc, ms_rctmf, ms_ndi)) + +rec = np.asarray(rec) + +plt.figure() +plt.title('increasing element size') +plt.plot(e_range, rec) +plt.legend(['rank.median', 'ctmf.median_filter', 'ndimage.percentile']) +plt.ylabel('time (ms)') +plt.xlabel('element radius') + +""" +.. image:: PLOT2RST.current_figure + +comparison of outcome of the three methods + +""" + +plt.figure() +plt.imshow(np.hstack((rc, rctmf, rndi))) +plt.xlabel('rank.median vs ctmf.median_filter vs ndimage.percentile') + +""" +.. image:: PLOT2RST.current_figure + +and increasing image size + +""" + +r = 9 +elem = disk(r + 1) + +rec = [] +s_range = [100, 200, 500, 1000] +for s in s_range: + a = (np.random.random((s, s)) * 256).astype('uint8') + (rc, ms_rc) = cr_med(a, elem) + rctmf, ms_rctmf = ctmf_med(a, r) + rndi, ms_ndi = ndi_med(a, r) + rec.append((ms_rc, ms_rctmf, ms_ndi)) + +rec = np.asarray(rec) + +plt.figure() +plt.title('increasing image size') +plt.plot(s_range, rec) +plt.legend(['rank.median', 'ctmf.median_filter', 'ndimage.percentile']) +plt.ylabel('time (ms)') +plt.xlabel('image size') + +""" +.. image:: PLOT2RST.current_figure + +""" + +plt.show() diff --git a/doc/examples/plot_16bitbilateral.py b/doc/examples/plot_16bitbilateral.py new file mode 100644 index 00000000..473fcadd --- /dev/null +++ b/doc/examples/plot_16bitbilateral.py @@ -0,0 +1,47 @@ +""" +============================== +Bilateral mean +============================== +This example compares + +* local mean +* percentile mean +* bilateral mean + +build on the local histogram distribution +local mean uses all pixels belonging to the structuring element to compute average gray level, +percentile mean uses only values between percentiles p0 and p1 (here 10% and 90%), +whereas bilateral mean uses only pixels of the structuring element having a gray level situated inside +g-s0 and g+s1 (here g-500 and g+500). +The filters are applied on a 16 bit image (actual bitdepth is 12bit). + +Percentile and usual mean give here similar results, these filters smooth the complete image (background and details). +Bilateral mean exhibits a high filtering rate for continuous area (i.e. background) while image higher frequencies +remains untouched. + +""" +import numpy as np +import matplotlib.pyplot as plt + +from skimage import data +from skimage.morphology import disk +import skimage.filter.rank as rank + +a16 = (data.coins()).astype('uint16') * 16 +selem = disk(20) + +f1 = rank.percentile_mean(a16, selem=selem, p0=.1, p1=.9) +f2 = rank.bilateral_mean(a16, selem=selem, s0=500, s1=500) +f3 = rank.mean(a16, selem=selem) + +# display results +fig, axes = plt.subplots(nrows=3, figsize=(15, 10)) +ax0, ax1, ax2 = axes + +ax0.imshow(np.hstack((a16, f1))) +ax0.set_title('percentile mean') +ax1.imshow(np.hstack((a16, f2))) +ax1.set_title('bilateral mean') +ax2.imshow(np.hstack((a16, f3))) +ax2.set_title('local mean') +plt.show() diff --git a/doc/examples/plot_entropy.py b/doc/examples/plot_entropy.py new file mode 100644 index 00000000..f019d79c --- /dev/null +++ b/doc/examples/plot_entropy.py @@ -0,0 +1,44 @@ +""" +=================== +Entropy +=================== + + +""" +from skimage import data +from skimage.filter.rank import entropy +from skimage.morphology import disk +import numpy as np +import matplotlib.pyplot as plt + +# defining a 8- and a 16-bit test images +a8 = data.camera() +a16 = data.camera().astype(np.uint16)*4 + +ent8 = entropy(a8,disk(5)) # pixel value contain 10x the local entropy +ent16 = entropy(a16,disk(5)) # pixel value contain 1000x the local entropy + +# display results +plt.figure(figsize=(10, 10)) + +plt.subplot(2,2,1) +plt.imshow(a8, cmap=plt.cm.gray) +plt.xlabel('8-bit image') +plt.colorbar() + +plt.subplot(2,2,2) +plt.imshow(ent8, cmap=plt.cm.jet) +plt.xlabel('entropy*10') +plt.colorbar() + +plt.subplot(2,2,3) +plt.imshow(a16, cmap=plt.cm.gray) +plt.xlabel('16-bit image') +plt.colorbar() + +plt.subplot(2,2,4) +plt.imshow(ent16, cmap=plt.cm.jet) +plt.xlabel('entropy*1000') +plt.colorbar() +plt.show() + diff --git a/doc/examples/plot_local_equalize.py b/doc/examples/plot_local_equalize.py new file mode 100644 index 00000000..bc458505 --- /dev/null +++ b/doc/examples/plot_local_equalize.py @@ -0,0 +1,84 @@ +""" +=============================== +Local Histogram Equalization +=============================== + +This examples enhances an image with low contrast, using a method called +*local histogram equalization*, which "spreads out the most frequent intensity +values" in an image . +The equalized image [1]_ has a roughly linear cumulative distribution function for each pixel neighborhood. +The local version [2]_ of the histogram equalization emphasized every local graylevel variations. + +.. [1] http://en.wikipedia.org/wiki/Histogram_equalization +.. [2] http://en.wikipedia.org/wiki/Adaptive_histogram_equalization + +""" + +from skimage import data +from skimage.util.dtype import dtype_range +from skimage import exposure +from skimage.morphology import disk + +import matplotlib.pyplot as plt + +import numpy as np +from skimage.filter import rank + + +def plot_img_and_hist(img, axes, bins=256): + """Plot an image along with its histogram and cumulative histogram. + + """ + ax_img, ax_hist = axes + ax_cdf = ax_hist.twinx() + + # Display image + ax_img.imshow(img, cmap=plt.cm.gray) + ax_img.set_axis_off() + + # Display histogram + ax_hist.hist(img.ravel(), bins=bins) + ax_hist.ticklabel_format(axis='y', style='scientific', scilimits=(0, 0)) + ax_hist.set_xlabel('Pixel intensity') + + xmin, xmax = dtype_range[img.dtype.type] + ax_hist.set_xlim(xmin, xmax) + + # Display cumulative distribution + img_cdf, bins = exposure.cumulative_distribution(img, bins) + ax_cdf.plot(bins, img_cdf, 'r') + + return ax_img, ax_hist, ax_cdf + + +# Load an example image +img = data.moon() + +# Contrast stretching +p2 = np.percentile(img, 2) +p98 = np.percentile(img, 98) +img_rescale = exposure.equalize(img) + +# Equalization +selem = disk(30) +img_eq = rank.equalize(img, selem=selem) + + +# Display results +f, axes = plt.subplots(2, 3, figsize=(8, 4)) + +ax_img, ax_hist, ax_cdf = plot_img_and_hist(img, axes[:, 0]) +ax_img.set_title('Low contrast image') +ax_hist.set_ylabel('Number of pixels') + +ax_img, ax_hist, ax_cdf = plot_img_and_hist(img_rescale, axes[:, 1]) +ax_img.set_title('Global equalise') + +ax_img, ax_hist, ax_cdf = plot_img_and_hist(img_eq, axes[:, 2]) +ax_img.set_title('Local equalize') +ax_cdf.set_ylabel('Fraction of total intensity') + + +# prevent overlap of y-axis labels +plt.subplots_adjust(wspace=0.4) +plt.show() diff --git a/doc/examples/plot_local_otsu.py b/doc/examples/plot_local_otsu.py new file mode 100644 index 00000000..968ce6e1 --- /dev/null +++ b/doc/examples/plot_local_otsu.py @@ -0,0 +1,49 @@ +""" +===================== +Local Otsu Threshold +===================== +This example shows how Otsu's threshold [1]_ method can be applied locally. +For each pixel, an "optimal" threshold is determined by maximizing the variance between two classes of pixels +of the local neighborhood defined by a structuring element. + +The example compares the local threshold with the global threshold. + +.. note: local threshold is much slower than global one. + +.. [1] http://en.wikipedia.org/wiki/Otsu's_method + +""" +import matplotlib.pyplot as plt + +from skimage import data +from skimage.morphology.selem import disk +import skimage.filter.rank as rank +from skimage.filter import threshold_otsu + + +p8 = data.page() + +radius = 10 +selem = disk(radius) + +loc_otsu = rank.otsu(p8, selem) +t_glob_otsu = threshold_otsu(p8) +glob_otsu = p8 >= t_glob_otsu + + +plt.figure() +plt.subplot(2, 2, 1) +plt.imshow(p8, cmap=plt.cm.gray) +plt.xlabel('original') +plt.colorbar() +plt.subplot(2, 2, 2) +plt.imshow(loc_otsu, cmap=plt.cm.gray) +plt.xlabel('local Otsu ($radius=%d$)' % radius) +plt.colorbar() +plt.subplot(2, 2, 3) +plt.imshow(p8 >= loc_otsu, cmap=plt.cm.gray) +plt.xlabel('original>=local Otsu' % t_glob_otsu) +plt.subplot(2, 2, 4) +plt.imshow(glob_otsu, cmap=plt.cm.gray) +plt.xlabel('global Otsu ($t=%d$)' % t_glob_otsu) +plt.show() diff --git a/doc/examples/plot_marked_watershed.py b/doc/examples/plot_marked_watershed.py new file mode 100644 index 00000000..e97280c7 --- /dev/null +++ b/doc/examples/plot_marked_watershed.py @@ -0,0 +1,54 @@ +""" +================================ +Markers for watershed transform +================================ + +The watershed is a classical algorithm used for **segmentation**, that +is, for separating different objects in an image. + +Here a marker image is build from the region of low gradient inside the image. + +See Wikipedia_ for more details on the algorithm. + +.. _Wikipedia: http://en.wikipedia.org/wiki/Watershed_(image_processing) + +""" + +from scipy import ndimage +import matplotlib.pyplot as plt +from skimage.morphology import watershed, disk +from skimage import data + +# original data +from skimage.filter import rank + +image = data.camera() + +# denoise image +denoised = rank.median(image, disk(2)) + +# find continuous region (low gradient) --> markers +markers = rank.gradient(denoised, disk(5)) < 10 +markers = ndimage.label(markers)[0] + +#local gradient +gradient = rank.gradient(denoised, disk(2)) + +# process the watershed +labels = watershed(gradient, markers) + +# display results +fig, axes = plt.subplots(ncols=4, figsize=(8, 2.7)) +ax0, ax1, ax2, ax3 = axes + +ax0.imshow(image, cmap=plt.cm.gray, interpolation='nearest') +ax1.imshow(gradient, cmap=plt.cm.spectral, interpolation='nearest') +ax2.imshow(markers, cmap=plt.cm.spectral, interpolation='nearest') +ax3.imshow(image, cmap=plt.cm.gray, interpolation='nearest') +ax3.imshow(labels, cmap=plt.cm.spectral, interpolation='nearest', alpha=.7) + +for ax in axes: + ax.axis('off') + +plt.subplots_adjust(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0, right=1) +plt.show() diff --git a/doc/examples/plot_watershed.py b/doc/examples/plot_watershed.py index a1cd18cf..50857b2c 100644 --- a/doc/examples/plot_watershed.py +++ b/doc/examples/plot_watershed.py @@ -26,7 +26,7 @@ See Wikipedia_ for more details on the algorithm. """ import numpy as np -from scipy import ndimage + import matplotlib.pyplot as plt from skimage.morphology import watershed, is_local_maximum diff --git a/skimage/filter/rank/.gitignore b/skimage/filter/rank/.gitignore new file mode 100644 index 00000000..08a24d5c --- /dev/null +++ b/skimage/filter/rank/.gitignore @@ -0,0 +1 @@ +demo/ \ No newline at end of file diff --git a/skimage/filter/rank/README.rst b/skimage/filter/rank/README.rst new file mode 100644 index 00000000..cdf8205c --- /dev/null +++ b/skimage/filter/rank/README.rst @@ -0,0 +1,32 @@ +To do +----- + +* add simple examples, adapt documentation on existing examples +* add/check existing doc +* adapting tests for each type of filter + +General remarks +--------------- + +Basically these filters compute local histogram for each pixel. A histogram is +built using a moving window in order to limit redundant computation. The path +followed by the moving window is given hereunder + + ...-----------------------\ +/--------------------------/ +\-------------------------- ... + +We compare cmorph.dilate to this histogram based method to show how +computational costs increase with respect to image size or structuring element +size. This implementation gives better results for large structuring elements. + +The local histogram is updated at each pixel as the structuring element window +moves by, i.e. only those pixels entering and leaving the structuring element +update the local histogram. The histogram size is 8-bit (256 bins) for 8-bit +images and 2 to 12-bit (up to 4096 bins) for 16-bit images depending on the +maximum value of the image. Pixel values higher than 4095 raise a ValueError. + +The filter is applied up to the image border, the neighboorhood used is adjusted +accordingly. The user may provide a mask image (same size as input image) where +non zero values are the part of the image participating in the histogram +computation. By default the entire image is filtered. diff --git a/skimage/filter/rank/__init__.py b/skimage/filter/rank/__init__.py new file mode 100644 index 00000000..30d936db --- /dev/null +++ b/skimage/filter/rank/__init__.py @@ -0,0 +1,3 @@ +from .rank import * +from .percentile_rank import * +from .bilateral_rank import * diff --git a/skimage/filter/rank/_core16.pxd b/skimage/filter/rank/_core16.pxd new file mode 100644 index 00000000..5f7de9df --- /dev/null +++ b/skimage/filter/rank/_core16.pxd @@ -0,0 +1,17 @@ +cimport numpy as np + + +cdef int int_max(int a, int b) +cdef int int_min(int a, int b) + + +# 16-bit core kernel receives extra information about data bitdepth +cdef void _core16(np.uint16_t kernel(Py_ssize_t *, float, np.uint16_t, + Py_ssize_t, Py_ssize_t, Py_ssize_t, float, + float, Py_ssize_t, Py_ssize_t), + 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, Py_ssize_t bitdepth, + float p0, float p1, Py_ssize_t s0, Py_ssize_t s1) except * diff --git a/skimage/filter/rank/_core16.pyx b/skimage/filter/rank/_core16.pyx new file mode 100644 index 00000000..aa959fd0 --- /dev/null +++ b/skimage/filter/rank/_core16.pyx @@ -0,0 +1,254 @@ +#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 +from _core8 cimport is_in_mask + + +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 void _core16(np.uint16_t kernel(Py_ssize_t *, float, np.uint16_t, + Py_ssize_t, Py_ssize_t, Py_ssize_t, float, + float, Py_ssize_t, Py_ssize_t), + 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, Py_ssize_t bitdepth, + float p0, float p1, Py_ssize_t s0, Py_ssize_t s1) except *: + """Compute histogram for each pixel neighborhood, apply kernel function and + use kernel function return value for output image. + """ + + cdef Py_ssize_t rows = image.shape[0] + cdef Py_ssize_t cols = image.shape[1] + cdef Py_ssize_t srows = selem.shape[0] + cdef Py_ssize_t scols = selem.shape[1] + + cdef Py_ssize_t centre_r = int(selem.shape[0] / 2) + shift_y + cdef Py_ssize_t 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 Py_ssize_t maxbin = maxbin_list[bitdepth] + cdef Py_ssize_t midbin = midbin_list[bitdepth] + + assert (image < maxbin).all() + + # define pointers to the 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 Py_ssize_t r, c, rr, cc, s, value, local_max, i, even_row + # number of pixels actually inside the neighborhood (float) + cdef float pop + + # allocate memory with malloc + cdef Py_ssize_t max_se = srows * scols + + # number of element in each attack border + cdef Py_ssize_t num_se_n, num_se_s, num_se_e, num_se_w + + # the current local histogram distribution + cdef Py_ssize_t * histo = malloc(maxbin * sizeof(Py_ssize_t)) + + # 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 Py_ssize_t * se_e_r = malloc(max_se * sizeof(Py_ssize_t)) + cdef Py_ssize_t * se_e_c = malloc(max_se * sizeof(Py_ssize_t)) + cdef Py_ssize_t * se_w_r = malloc(max_se * sizeof(Py_ssize_t)) + cdef Py_ssize_t * se_w_c = malloc(max_se * sizeof(Py_ssize_t)) + cdef Py_ssize_t * se_n_r = malloc(max_se * sizeof(Py_ssize_t)) + cdef Py_ssize_t * se_n_c = malloc(max_se * sizeof(Py_ssize_t)) + cdef Py_ssize_t * se_s_r = malloc(max_se * sizeof(Py_ssize_t)) + cdef Py_ssize_t * se_s_c = malloc(max_se * sizeof(Py_ssize_t)) + + # 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 - centre_r + cc = c - centre_c + if selem[r, c]: + if is_in_mask(rows, cols, rr, cc, mask_data): + histogram_increment(histo, &pop, image_data[rr * cols + cc]) + + r = 0 + c = 0 + # kernel ------------------------------------------- + out_data[r * cols + c] = kernel(histo, pop, image_data[r * cols + c], + bitdepth, maxbin, midbin, p0, p1, 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] + 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 + if is_in_mask(rows, cols, rr, cc, mask_data): + histogram_decrement(histo, &pop, image_data[rr * cols + cc]) + + # kernel ------------------------------------------- + out_data[r * cols + c] = kernel( + histo, pop, image_data[r * cols + c], + bitdepth, maxbin, midbin, p0, p1, 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] + cc = c + se_s_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_n): + rr = r + se_n_r[s] - 1 + cc = c + se_n_c[s] + if is_in_mask(rows, cols, rr, cc, mask_data): + histogram_decrement(histo, &pop, image_data[rr * cols + cc]) + + # kernel ------------------------------------------- + out_data[r * cols + c] = kernel(histo, pop, image_data[r * cols + c], + bitdepth, maxbin, midbin, p0, p1, 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] + cc = c + se_w_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_e): + rr = r + se_e_r[s] + cc = c + se_e_c[s] + 1 + if is_in_mask(rows, cols, rr, cc, mask_data): + histogram_decrement(histo, &pop, image_data[rr * cols + cc]) + + # kernel ------------------------------------------- + out_data[r * cols + c] = kernel( + histo, pop, image_data[r * cols + c], + bitdepth, maxbin, midbin, p0, p1, 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] + cc = c + se_s_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_n): + rr = r + se_n_r[s] - 1 + cc = c + se_n_c[s] + if is_in_mask(rows, cols, rr, cc, mask_data): + histogram_decrement(histo, &pop, image_data[rr * cols + cc]) + + # kernel ------------------------------------------- + out_data[r * cols + c] = kernel(histo, pop, image_data[r * cols + c], + bitdepth, maxbin, midbin, p0, p1, 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) diff --git a/skimage/filter/rank/_core8.pxd b/skimage/filter/rank/_core8.pxd new file mode 100644 index 00000000..38236b87 --- /dev/null +++ b/skimage/filter/rank/_core8.pxd @@ -0,0 +1,22 @@ +cimport numpy as np + + +cdef np.uint8_t uint8_max(np.uint8_t a, np.uint8_t b) +cdef np.uint8_t uint8_min(np.uint8_t a, np.uint8_t b) + + +cdef 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) + + +# 8-bit core kernel receives extra information about data inferior and superior +# percentiles +cdef void _core8(np.uint8_t kernel(Py_ssize_t *, float, np.uint8_t, float, + float, Py_ssize_t, Py_ssize_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, float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1) except * diff --git a/skimage/filter/rank/_core8.pyx b/skimage/filter/rank/_core8.pyx new file mode 100644 index 00000000..9f4bc693 --- /dev/null +++ b/skimage/filter/rank/_core8.pyx @@ -0,0 +1,256 @@ +#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 + + +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 void histogram_increment(Py_ssize_t * histo, float * pop, + np.uint8_t value): + histo[value] += 1 + pop[0] += 1 + + +cdef inline void histogram_decrement(Py_ssize_t * histo, float * pop, + np.uint8_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): + """Check whether given coordinate is within image and mask is true.""" + if r < 0 or r > rows - 1 or c < 0 or c > cols - 1: + return 0 + else: + if mask[r * cols + c]: + return 1 + else: + return 0 + + +cdef void _core8(np.uint8_t kernel(Py_ssize_t *, float, np.uint8_t, float, + float, Py_ssize_t, Py_ssize_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, float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1) except *: + """Compute histogram for each pixel neighborhood, apply kernel function and + use kernel function return value for output image. + """ + + cdef Py_ssize_t rows = image.shape[0] + cdef Py_ssize_t cols = image.shape[1] + cdef Py_ssize_t srows = selem.shape[0] + cdef Py_ssize_t scols = selem.shape[1] + + cdef Py_ssize_t centre_r = int(selem.shape[0] / 2) + shift_y + cdef Py_ssize_t 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 + + # define pointers to the 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 Py_ssize_t r, c, rr, cc, s, value, local_max, i, even_row + + # number of pixels actually inside the neighborhood (float) + cdef float pop + + # allocate memory with malloc + cdef Py_ssize_t max_se = srows * scols + + # number of element in each attack border + cdef Py_ssize_t num_se_n, num_se_s, num_se_e, num_se_w + + # the current local histogram distribution + cdef Py_ssize_t * histo = malloc(256 * sizeof(Py_ssize_t)) + + # 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 Py_ssize_t * se_e_r = malloc(max_se * sizeof(Py_ssize_t)) + cdef Py_ssize_t * se_e_c = malloc(max_se * sizeof(Py_ssize_t)) + cdef Py_ssize_t * se_w_r = malloc(max_se * sizeof(Py_ssize_t)) + cdef Py_ssize_t * se_w_c = malloc(max_se * sizeof(Py_ssize_t)) + cdef Py_ssize_t * se_n_r = malloc(max_se * sizeof(Py_ssize_t)) + cdef Py_ssize_t * se_n_c = malloc(max_se * sizeof(Py_ssize_t)) + cdef Py_ssize_t * se_s_r = malloc(max_se * sizeof(Py_ssize_t)) + cdef Py_ssize_t * se_s_c = malloc(max_se * sizeof(Py_ssize_t)) + + # 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 (kernel is centered on the first row and + # column) + for i in range(256): + histo[i] = 0 + + pop = 0 + + for r in range(srows): + for c in range(scols): + rr = r - centre_r + cc = c - centre_c + if selem[r, c]: + if is_in_mask(rows, cols, rr, cc, mask_data): + histogram_increment(histo, &pop, image_data[rr * cols + cc]) + + r = 0 + c = 0 + # kernel ------------------------------------------------------------------- + out_data[r * cols + c] = kernel(histo, pop, image_data[r * cols + c], + p0, p1, 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] + 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 + if is_in_mask(rows, cols, rr, cc, mask_data): + histogram_decrement(histo, &pop, image_data[rr * cols + cc]) + + # kernel ----------------------------------------------------------- + out_data[r * cols + c] = \ + kernel(histo, pop, image_data[r * cols + c], p0, p1, 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] + cc = c + se_s_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_n): + rr = r + se_n_r[s] - 1 + cc = c + se_n_c[s] + if is_in_mask(rows, cols, rr, cc, mask_data): + histogram_decrement(histo, &pop, image_data[rr * cols + cc]) + + # kernel --------------------------------------------------------------- + out_data[r * cols + c] = kernel(histo, pop, image_data[r * cols + c], + p0, p1, 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] + cc = c + se_w_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_e): + rr = r + se_e_r[s] + cc = c + se_e_c[s] + 1 + if is_in_mask(rows, cols, rr, cc, mask_data): + histogram_decrement(histo, &pop, image_data[rr * cols + cc]) + + # kernel ----------------------------------------------------------- + out_data[r * cols + c] = kernel( + histo, pop, image_data[r * cols + c], p0, p1, 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] + cc = c + se_s_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_n): + rr = r + se_n_r[s] - 1 + cc = c + se_n_c[s] + if is_in_mask(rows, cols, rr, cc, mask_data): + histogram_decrement(histo, &pop, image_data[rr * cols + cc]) + + # kernel --------------------------------------------------------------- + out_data[r * cols + c] = kernel(histo, pop, image_data[r * cols + c], + p0, p1, 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) diff --git a/skimage/filter/rank/_crank16.pyx b/skimage/filter/rank/_crank16.pyx new file mode 100644 index 00000000..e81bf81f --- /dev/null +++ b/skimage/filter/rank/_crank16.pyx @@ -0,0 +1,420 @@ +#cython: cdivision=True +#cython: boundscheck=False +#cython: nonecheck=False +#cython: wraparound=False + +import numpy as np +cimport numpy as np +from libc.math cimport log2 +from skimage.filter.rank._core16 cimport _core16 + + +# ----------------------------------------------------------------- +# kernels uint16 take extra parameter for defining the bitdepth +# ----------------------------------------------------------------- + + +cdef inline np.uint16_t kernel_autolevel(Py_ssize_t * histo, float pop, + np.uint16_t g, Py_ssize_t bitdepth, + Py_ssize_t maxbin, Py_ssize_t midbin, + float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + cdef Py_ssize_t i, imin, imax, delta + + if pop: + for i in range(maxbin - 1, -1, -1): + if histo[i]: + imax = i + break + for i in range(maxbin): + if histo[i]: + imin = i + break + delta = imax - imin + if delta > 0: + return < np.uint16_t > (1. * (maxbin - 1) * (g - imin) / delta) + else: + return < np.uint16_t > (imax - imin) + + +cdef inline np.uint16_t kernel_bottomhat(Py_ssize_t * histo, float pop, + np.uint16_t g, Py_ssize_t bitdepth, + Py_ssize_t maxbin, Py_ssize_t midbin, + float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + cdef Py_ssize_t i + + if pop: + for i in range(maxbin): + if histo[i]: + break + + return < np.uint16_t > (g - i) + else: + return < np.uint16_t > (0) + +cdef inline np.uint16_t kernel_equalize(Py_ssize_t * histo, float pop, + np.uint16_t g, Py_ssize_t bitdepth, + Py_ssize_t maxbin, Py_ssize_t midbin, + float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + cdef Py_ssize_t i + cdef float sum = 0. + + if pop: + for i in range(maxbin): + sum += histo[i] + if i >= g: + break + + return < np.uint16_t > (((maxbin - 1) * sum) / pop) + else: + return < np.uint16_t > (0) + + +cdef inline np.uint16_t kernel_gradient(Py_ssize_t * histo, float pop, + np.uint16_t g, Py_ssize_t bitdepth, + Py_ssize_t maxbin, Py_ssize_t midbin, + float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + cdef Py_ssize_t i, imin, imax + + if pop: + for i in range(maxbin - 1, -1, -1): + if histo[i]: + imax = i + break + for i in range(maxbin): + if histo[i]: + imin = i + break + return < np.uint16_t > (imax - imin) + else: + return < np.uint16_t > (0) + + +cdef inline np.uint16_t kernel_maximum(Py_ssize_t * histo, float pop, + np.uint16_t g, Py_ssize_t bitdepth, + Py_ssize_t maxbin, Py_ssize_t midbin, + float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + cdef Py_ssize_t i + + if pop: + for i in range(maxbin - 1, -1, -1): + if histo[i]: + return < np.uint16_t > (i) + + return < np.uint16_t > (0) + + +cdef inline np.uint16_t kernel_mean(Py_ssize_t * histo, float pop, + np.uint16_t g, Py_ssize_t bitdepth, + Py_ssize_t maxbin, Py_ssize_t midbin, + float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + cdef Py_ssize_t i + cdef float mean = 0. + + if pop: + for i in range(maxbin): + mean += histo[i] * i + return < np.uint16_t > (mean / pop) + else: + return < np.uint16_t > (0) + + +cdef inline np.uint16_t kernel_meansubstraction(Py_ssize_t * histo, + float pop, + np.uint16_t g, + Py_ssize_t bitdepth, + Py_ssize_t maxbin, + Py_ssize_t midbin, + float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + cdef Py_ssize_t i + cdef float mean = 0. + + if pop: + for i in range(maxbin): + mean += histo[i] * i + return < np.uint16_t > ((g - mean / pop) / 2. + (midbin - 1)) + else: + return < np.uint16_t > (0) + + +cdef inline np.uint16_t kernel_median(Py_ssize_t * histo, float pop, + np.uint16_t g, Py_ssize_t bitdepth, + Py_ssize_t maxbin, Py_ssize_t midbin, + float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + cdef Py_ssize_t i + cdef float sum = pop / 2.0 + + if pop: + for i in range(maxbin): + if histo[i]: + sum -= histo[i] + if sum < 0: + return < np.uint16_t > (i) + else: + return < np.uint16_t > (0) + + +cdef inline np.uint16_t kernel_minimum(Py_ssize_t * histo, float pop, + np.uint16_t g, Py_ssize_t bitdepth, + Py_ssize_t maxbin, Py_ssize_t midbin, + float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + cdef Py_ssize_t i + + if pop: + for i in range(maxbin): + if histo[i]: + return < np.uint16_t > (i) + else: + return < np.uint16_t > (0) + + +cdef inline np.uint16_t kernel_modal(Py_ssize_t * histo, float pop, + np.uint16_t g, Py_ssize_t bitdepth, + Py_ssize_t maxbin, Py_ssize_t midbin, + float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + cdef Py_ssize_t hmax = 0, imax = 0 + + if pop: + for i in range(maxbin): + if histo[i] > hmax: + hmax = histo[i] + imax = i + return < np.uint16_t > (imax) + else: + return < np.uint16_t > (0) + + +cdef inline np.uint16_t kernel_morph_contr_enh(Py_ssize_t * histo, + float pop, + np.uint16_t g, + Py_ssize_t bitdepth, + Py_ssize_t maxbin, + Py_ssize_t midbin, + float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + cdef Py_ssize_t i, imin, imax + + if pop: + for i in range(maxbin - 1, -1, -1): + if histo[i]: + imax = i + break + for i in range(maxbin): + if histo[i]: + imin = i + break + if imax - g < g - imin: + return < np.uint16_t > (imax) + else: + return < np.uint16_t > (imin) + else: + return < np.uint16_t > (0) + + +cdef inline np.uint16_t kernel_pop(Py_ssize_t * histo, float pop, + np.uint16_t g, Py_ssize_t bitdepth, + Py_ssize_t maxbin, Py_ssize_t midbin, + float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + return < np.uint16_t > (pop) + + +cdef inline np.uint16_t kernel_threshold(Py_ssize_t * histo, float pop, + np.uint16_t g, Py_ssize_t bitdepth, + Py_ssize_t maxbin, Py_ssize_t midbin, + float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + cdef Py_ssize_t i + cdef float mean = 0. + + if pop: + for i in range(maxbin): + mean += histo[i] * i + return < np.uint16_t > (g > (mean / pop)) + else: + return < np.uint16_t > (0) + + +cdef inline np.uint16_t kernel_tophat(Py_ssize_t * histo, float pop, + np.uint16_t g, Py_ssize_t bitdepth, + Py_ssize_t maxbin, Py_ssize_t midbin, + float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + cdef Py_ssize_t i + + if pop: + for i in range(maxbin - 1, -1, -1): + if histo[i]: + break + + return < np.uint16_t > (i - g) + else: + return < np.uint16_t > (0) + +cdef inline np.uint16_t kernel_entropy(Py_ssize_t * histo, float pop, + np.uint16_t g, Py_ssize_t bitdepth, + Py_ssize_t maxbin, Py_ssize_t midbin, + float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + cdef Py_ssize_t i + cdef float e, p + + if pop: + e = 0. + + for i in range(maxbin): + p = histo[i] / pop + if p > 0: + e -= p * log2(p) + + return < np.uint16_t > e * 1000 + else: + return < np.uint16_t > (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, Py_ssize_t bitdepth=8): + _core16(kernel_autolevel, image, selem, mask, out, shift_x, shift_y, + bitdepth, 0, 0, 0, 0) + + +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, Py_ssize_t bitdepth=8): + _core16(kernel_bottomhat, image, selem, mask, out, shift_x, shift_y, + bitdepth, 0, 0, 0, 0) + + +def equalize(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, Py_ssize_t bitdepth=8): + _core16(kernel_equalize, image, selem, mask, out, shift_x, shift_y, + bitdepth, 0, 0, 0, 0) + + +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, Py_ssize_t bitdepth=8): + _core16(kernel_gradient, image, selem, mask, out, shift_x, shift_y, + bitdepth, 0, 0, 0, 0) + + +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, Py_ssize_t bitdepth=8): + _core16(kernel_maximum, image, selem, mask, out, shift_x, shift_y, + bitdepth, 0, 0, 0, 0) + + +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, Py_ssize_t bitdepth=8): + _core16(kernel_mean, image, selem, mask, out, shift_x, shift_y, + bitdepth, 0, 0, 0, 0) + + +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, Py_ssize_t bitdepth=8): + _core16(kernel_meansubstraction, image, selem, mask, out, shift_x, shift_y, + bitdepth, 0, 0, 0, 0) + + +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, Py_ssize_t bitdepth=8): + _core16(kernel_median, image, selem, mask, out, shift_x, shift_y, + bitdepth, 0, 0, 0, 0) + + +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, Py_ssize_t bitdepth=8): + _core16(kernel_minimum, image, selem, mask, out, shift_x, shift_y, + bitdepth, 0, 0, 0, 0) + + +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, Py_ssize_t bitdepth=8): + _core16(kernel_morph_contr_enh, image, selem, mask, out, shift_x, shift_y, + bitdepth, 0, 0, 0, 0) + + +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, Py_ssize_t bitdepth=8): + _core16(kernel_modal, image, selem, mask, out, shift_x, shift_y, + bitdepth, 0, 0, 0, 0) + + +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, Py_ssize_t bitdepth=8): + _core16(kernel_pop, image, selem, mask, out, shift_x, shift_y, + bitdepth, 0, 0, 0, 0) + + +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, Py_ssize_t bitdepth=8): + _core16(kernel_threshold, image, selem, mask, out, shift_x, shift_y, + bitdepth, 0, 0, 0, 0) + + +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, Py_ssize_t bitdepth=8): + _core16(kernel_tophat, image, selem, mask, out, shift_x, shift_y, + bitdepth, 0, 0, 0, 0) + + +def entropy(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, Py_ssize_t bitdepth=8): + _core16(kernel_entropy, image, selem, mask, out, shift_x, shift_y, + bitdepth, 0, 0, 0, 0) diff --git a/skimage/filter/rank/_crank16_bilateral.pyx b/skimage/filter/rank/_crank16_bilateral.pyx new file mode 100644 index 00000000..c71ccc5c --- /dev/null +++ b/skimage/filter/rank/_crank16_bilateral.pyx @@ -0,0 +1,80 @@ +#cython: cdivision=True +#cython: boundscheck=False +#cython: nonecheck=False +#cython: wraparound=False + +import numpy as np +cimport numpy as np +from skimage.filter.rank._core16 cimport _core16 + + +# ----------------------------------------------------------------- +# kernels uint16 take extra parameter for defining the bitdepth +# ----------------------------------------------------------------- + + +cdef inline np.uint16_t kernel_mean(Py_ssize_t * histo, float pop, + np.uint16_t g, Py_ssize_t bitdepth, + Py_ssize_t maxbin, Py_ssize_t midbin, + float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + + cdef int i, bilat_pop = 0 + cdef float mean = 0. + + if pop: + for i in range(maxbin): + if (g > (i - s0)) and (g < (i + s1)): + bilat_pop += histo[i] + mean += histo[i] * i + if bilat_pop: + return < np.uint16_t > (mean / bilat_pop) + else: + return < np.uint16_t > (0) + else: + return < np.uint16_t > (0) + + +cdef inline np.uint16_t kernel_pop(Py_ssize_t * histo, float pop, + np.uint16_t g, Py_ssize_t bitdepth, + Py_ssize_t maxbin, Py_ssize_t midbin, + float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + + cdef int i, bilat_pop = 0 + + if pop: + for i in range(maxbin): + if (g > (i - s0)) and (g < (i + s1)): + bilat_pop += histo[i] + return < np.uint16_t > (bilat_pop) + else: + return < np.uint16_t > (0) + + +# ----------------------------------------------------------------- +# python wrappers +# ----------------------------------------------------------------- + + +def mean(np.ndarray[np.uint16_t, ndim=2] image, + np.ndarray[np.uint8_t, ndim=2] selem, + np.ndarray[np.uint8_t, ndim=2] mask=None, + np.ndarray[np.uint16_t, ndim=2] out=None, + char shift_x=0, char shift_y=0, int bitdepth=8, int s0=1, int s1=1): + """average greylevel (clipped on uint8) + """ + _core16(kernel_mean, image, selem, mask, out, shift_x, shift_y, + bitdepth, 0., 0., s0, s1) + + +def pop(np.ndarray[np.uint16_t, ndim=2] image, + np.ndarray[np.uint8_t, ndim=2] selem, + np.ndarray[np.uint8_t, ndim=2] mask=None, + np.ndarray[np.uint16_t, ndim=2] out=None, + char shift_x=0, char shift_y=0, int bitdepth=8, int s0=1, int s1=1): + """returns the number of actual pixels of the structuring element inside + the mask + """ + _core16(kernel_pop, image, selem, mask, out, shift_x, shift_y, + bitdepth, .0, .0, s0, s1) diff --git a/skimage/filter/rank/_crank16_percentiles.pyx b/skimage/filter/rank/_crank16_percentiles.pyx new file mode 100644 index 00000000..220d2386 --- /dev/null +++ b/skimage/filter/rank/_crank16_percentiles.pyx @@ -0,0 +1,327 @@ +#cython: cdivision=True +#cython: boundscheck=False +#cython: nonecheck=False +#cython: wraparound=False + +import numpy as np +cimport numpy as np +from skimage.filter.rank._core16 cimport _core16, int_min, int_max + + +# ----------------------------------------------------------------- +# kernels uint16 (SOFT version using percentiles) +# ----------------------------------------------------------------- + + +cdef inline np.uint16_t kernel_autolevel(Py_ssize_t * histo, float pop, + np.uint16_t g, Py_ssize_t bitdepth, + Py_ssize_t maxbin, Py_ssize_t midbin, + float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + + 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 delta > 0: + return < np.uint16_t > (1.0 * (maxbin - 1) + * (int_min(int_max(imin, g), imax) - imin) / delta) + else: + return < np.uint16_t > (imax - imin) + else: + return < np.uint16_t > (0) + + +cdef inline np.uint16_t kernel_gradient(Py_ssize_t * histo, float pop, + np.uint16_t g, Py_ssize_t bitdepth, + Py_ssize_t maxbin, Py_ssize_t midbin, + float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + + 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 < np.uint16_t > (imax - imin) + else: + return < np.uint16_t > (0) + + +cdef inline np.uint16_t kernel_mean(Py_ssize_t * histo, float pop, + np.uint16_t g, Py_ssize_t bitdepth, + Py_ssize_t maxbin, Py_ssize_t midbin, + float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + + 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 < np.uint16_t > (1.0 * mean / n) + else: + return < np.uint16_t > (0) + else: + return < np.uint16_t > (0) + + +cdef inline np.uint16_t kernel_mean_substraction(Py_ssize_t * histo, + float pop, + np.uint16_t g, + Py_ssize_t bitdepth, + Py_ssize_t maxbin, + Py_ssize_t midbin, + float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + + 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 < np.uint16_t > ((g - (mean / n)) * .5 + midbin) + else: + return < np.uint16_t > (0) + else: + return < np.uint16_t > (0) + + +cdef inline np.uint16_t kernel_morph_contr_enh(Py_ssize_t * histo, + float pop, + np.uint16_t g, + Py_ssize_t bitdepth, + Py_ssize_t maxbin, + Py_ssize_t midbin, + float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + + 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 < np.uint16_t > imax + if g < imin: + return < np.uint16_t > imin + if imax - g < g - imin: + return < np.uint16_t > imax + else: + return < np.uint16_t > imin + else: + return < np.uint16_t > (0) + + +cdef inline np.uint16_t kernel_percentile(Py_ssize_t * histo, float pop, + np.uint16_t g, Py_ssize_t bitdepth, + Py_ssize_t maxbin, Py_ssize_t midbin, + float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + + cdef int i + cdef float sum = 0. + + if pop: + for i in range(maxbin): + sum += histo[i] + if sum >= p0 * pop: + break + + return < np.uint16_t > (i) + else: + return < np.uint16_t > (0) + + +cdef inline np.uint16_t kernel_pop(Py_ssize_t * histo, float pop, + np.uint16_t g, Py_ssize_t bitdepth, + Py_ssize_t maxbin, Py_ssize_t midbin, + float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + + 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 < np.uint16_t > (n) + else: + return < np.uint16_t > (0) + + +cdef inline np.uint16_t kernel_threshold(Py_ssize_t * histo, float pop, + np.uint16_t g, Py_ssize_t bitdepth, + Py_ssize_t maxbin, Py_ssize_t midbin, + float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + + cdef int i + cdef float sum = 0. + + if pop: + for i in range(maxbin): + sum += histo[i] + if sum >= p0 * pop: + break + + return < np.uint16_t > ((maxbin - 1) * (g >= i)) + else: + return < np.uint16_t > (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 + """ + _core16(kernel_autolevel, image, selem, mask, out, shift_x, shift_y, + bitdepth, p0, p1, 0, 0) + + +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 + """ + _core16(kernel_gradient, image, selem, mask, out, shift_x, shift_y, + bitdepth, p0, p1, 0, 0) + + +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 + """ + _core16(kernel_mean, image, selem, mask, out, shift_x, shift_y, + bitdepth, p0, p1, 0, 0) + + +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 + """ + _core16( + kernel_mean_substraction, image, selem, mask, out, shift_x, shift_y, + bitdepth, p0, p1, 0, 0) + + +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 + """ + _core16(kernel_morph_contr_enh, image, selem, mask, out, shift_x, shift_y, + bitdepth, p0, p1, 0, 0) + + +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 + """ + _core16(kernel_percentile, image, selem, mask, out, shift_x, shift_y, + bitdepth, p0, p1, 0, 0) + + +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] + """ + _core16(kernel_pop, image, selem, mask, out, shift_x, shift_y, + bitdepth, p0, p1, 0, 0) + + +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 + """ + _core16(kernel_threshold, image, selem, mask, out, shift_x, shift_y, + bitdepth, p0, p1, 0, 0) diff --git a/skimage/filter/rank/_crank8.pyx b/skimage/filter/rank/_crank8.pyx new file mode 100644 index 00000000..8bff9703 --- /dev/null +++ b/skimage/filter/rank/_crank8.pyx @@ -0,0 +1,481 @@ +#cython: cdivision=True +#cython: boundscheck=False +#cython: nonecheck=False +#cython: wraparound=False + +import numpy as np +cimport numpy as np +from libc.math cimport log2 +from skimage.filter.rank._core8 cimport _core8 + + +# ----------------------------------------------------------------- +# kernels uint8 +# ----------------------------------------------------------------- + + +cdef inline np.uint8_t kernel_autolevel(Py_ssize_t * histo, float pop, + np.uint8_t g, float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + + cdef Py_ssize_t 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 < np.uint8_t > (255. * (g - imin) / delta) + else: + return < np.uint8_t > (imax - imin) + else: + return < np.uint8_t > (0) + + +cdef inline np.uint8_t kernel_bottomhat(Py_ssize_t * histo, float pop, + np.uint8_t g, float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + + cdef Py_ssize_t i + + if pop: + for i in range(256): + if histo[i]: + break + + return < np.uint8_t > (g - i) + else: + return < np.uint8_t > (0) + + +cdef inline np.uint8_t kernel_equalize(Py_ssize_t * histo, float pop, + np.uint8_t g, float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + + cdef Py_ssize_t i + cdef float sum = 0. + + if pop: + for i in range(256): + sum += histo[i] + if i >= g: + break + + return < np.uint8_t > ((255 * sum) / pop) + else: + return < np.uint8_t > (0) + + +cdef inline np.uint8_t kernel_gradient(Py_ssize_t * histo, float pop, + np.uint8_t g, float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + + cdef Py_ssize_t 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 < np.uint8_t > (imax - imin) + else: + return < np.uint8_t > (0) + + +cdef inline np.uint8_t kernel_maximum(Py_ssize_t * histo, float pop, + np.uint8_t g, float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + + cdef Py_ssize_t i + + if pop: + for i in range(255, -1, -1): + if histo[i]: + return < np.uint8_t > (i) + else: + return < np.uint8_t > (0) + + +cdef inline np.uint8_t kernel_mean(Py_ssize_t * histo, float pop, + np.uint8_t g, float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + + cdef Py_ssize_t i + cdef float mean = 0. + + if pop: + for i in range(256): + mean += histo[i] * i + return < np.uint8_t > (mean / pop) + else: + return < np.uint8_t > (0) + + +cdef inline np.uint8_t kernel_meansubstraction(Py_ssize_t * histo, float pop, + np.uint8_t g, float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + + cdef Py_ssize_t i + cdef float mean = 0. + + if pop: + for i in range(256): + mean += histo[i] * i + return < np.uint8_t > ((g - mean / pop) / 2. + 127) + else: + return < np.uint8_t > (0) + + +cdef inline np.uint8_t kernel_median(Py_ssize_t * histo, float pop, + np.uint8_t g, float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + + cdef Py_ssize_t i + cdef float sum = pop / 2.0 + + if pop: + for i in range(256): + if histo[i]: + sum -= histo[i] + if sum < 0: + return < np.uint8_t > (i) + else: + return < np.uint8_t > (0) + + +cdef inline np.uint8_t kernel_minimum(Py_ssize_t * histo, float pop, + np.uint8_t g, float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + + cdef Py_ssize_t i + + if pop: + for i in range(256): + if histo[i]: + return < np.uint8_t > (i) + else: + return < np.uint8_t > (0) + + +cdef inline np.uint8_t kernel_modal(Py_ssize_t * histo, float pop, + np.uint8_t g, float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + + cdef Py_ssize_t hmax = 0, imax = 0 + + if pop: + for i in range(256): + if histo[i] > hmax: + hmax = histo[i] + imax = i + return < np.uint8_t > (imax) + else: + return < np.uint8_t > (0) + + +cdef inline np.uint8_t kernel_morph_contr_enh(Py_ssize_t * histo, float pop, + np.uint8_t g, float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + + cdef Py_ssize_t 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 < np.uint8_t > (imax) + else: + return < np.uint8_t > (imin) + else: + return < np.uint8_t > (0) + + +cdef inline np.uint8_t kernel_pop(Py_ssize_t * histo, float pop, + np.uint8_t g, float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + + return < np.uint8_t > (pop) + + +cdef inline np.uint8_t kernel_threshold(Py_ssize_t * histo, float pop, + np.uint8_t g, float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + + cdef Py_ssize_t i + cdef float mean = 0. + + if pop: + for i in range(256): + mean += histo[i] * i + return < np.uint8_t > (g > (mean / pop)) + else: + return < np.uint8_t > (0) + + +cdef inline np.uint8_t kernel_tophat(Py_ssize_t * histo, float pop, + np.uint8_t g, float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + + cdef Py_ssize_t i + + if pop: + for i in range(255, -1, -1): + if histo[i]: + break + + return < np.uint8_t > (i - g) + else: + return < np.uint8_t > (0) + +cdef inline np.uint8_t kernel_noise_filter(Py_ssize_t * histo, float pop, + np.uint8_t g, float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + + cdef Py_ssize_t i + cdef Py_ssize_t min_i + + # early stop if at least one pixel of the neighborhood has the same g + if histo[g] > 0: + return < np.uint8_t > 0 + + for i in range(g, -1, -1): + if histo[i]: + break + min_i = g - i + for i in range(g, 256): + if histo[i]: + break + if i - g < min_i: + return < np.uint8_t > (i - g) + else: + return < np.uint8_t > min_i + + +cdef inline np.uint8_t kernel_entropy(Py_ssize_t * histo, float pop, + np.uint8_t g, float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + cdef Py_ssize_t i + cdef float e, p + + if pop: + e = 0. + + for i in range(256): + p = histo[i] / pop + if p > 0: + e -= p * log2(p) + + return < np.uint8_t > e * 10 + else: + return < np.uint8_t > (0) + +cdef inline np.uint8_t kernel_otsu(Py_ssize_t * histo, float pop, np.uint8_t g, + float p0, float p1, Py_ssize_t s0, + Py_ssize_t s1): + cdef Py_ssize_t i + cdef Py_ssize_t max_i + cdef float P, mu1, mu2, q1, new_q1, sigma_b, max_sigma_b + cdef float mu = 0. + + # compute local mean + if pop: + for i in range(256): + mu += histo[i] * i + mu = (mu / pop) + else: + return < np.uint8_t > (0) + + # maximizing the between class variance + max_i = 0 + q1 = histo[0] / pop + m1 = 0. + max_sigma_b = 0. + + for i in range(1, 256): + P = histo[i] / pop + new_q1 = q1 + P + if new_q1 > 0: + mu1 = (q1 * mu1 + i * P) / new_q1 + mu2 = (mu - new_q1 * mu1) / (1. - new_q1) + sigma_b = new_q1 * (1. - new_q1) * (mu1 - mu2) ** 2 + if sigma_b > max_sigma_b: + max_sigma_b = sigma_b + max_i = i + q1 = new_q1 + + return < np.uint8_t > max_i + + +# ----------------------------------------------------------------- +# python wrappers +# used only internally +# ----------------------------------------------------------------- + + +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): + _core8(kernel_autolevel, image, selem, mask, out, shift_x, shift_y, + 0, 0, 0, 0) + + +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): + _core8(kernel_bottomhat, image, selem, mask, out, shift_x, shift_y, + 0, 0, 0, 0) + + +def equalize(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): + _core8(kernel_equalize, image, selem, mask, out, shift_x, shift_y, + 0, 0, 0, 0) + + +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): + _core8(kernel_gradient, image, selem, mask, out, shift_x, shift_y, + 0, 0, 0, 0) + + +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): + _core8(kernel_maximum, image, selem, mask, out, shift_x, shift_y, + 0, 0, 0, 0) + + +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): + _core8(kernel_mean, image, selem, mask, out, shift_x, shift_y, + 0, 0, 0, 0) + + +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): + _core8(kernel_meansubstraction, image, selem, mask, out, shift_x, shift_y, + 0, 0, 0, 0) + + +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): + _core8(kernel_median, image, selem, mask, out, shift_x, shift_y, + 0, 0, 0, 0) + + +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): + _core8(kernel_minimum, image, selem, mask, out, shift_x, shift_y, + 0, 0, 0, 0) + + +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): + _core8(kernel_morph_contr_enh, image, selem, mask, out, shift_x, shift_y, + 0, 0, 0, 0) + + +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): + _core8(kernel_modal, image, selem, mask, out, shift_x, shift_y, + 0, 0, 0, 0) + + +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): + _core8(kernel_pop, image, selem, mask, out, shift_x, shift_y, + 0, 0, 0, 0) + + +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): + _core8(kernel_threshold, image, selem, mask, out, shift_x, shift_y, 0, 0, + 0, 0) + + +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): + _core8(kernel_tophat, image, selem, mask, out, shift_x, shift_y, + 0, 0, 0, 0) + + +def noise_filter(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): + _core8(kernel_noise_filter, image, selem, mask, out, shift_x, shift_y, + 0, 0, 0, 0) + + +def entropy(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): + _core8(kernel_entropy, image, selem, mask, out, shift_x, shift_y, + 0, 0, 0, 0) + + +def otsu(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): + _core8(kernel_otsu, image, selem, mask, out, shift_x, shift_y, + 0, 0, 0, 0) diff --git a/skimage/filter/rank/_crank8_percentiles.pyx b/skimage/filter/rank/_crank8_percentiles.pyx new file mode 100644 index 00000000..31b00742 --- /dev/null +++ b/skimage/filter/rank/_crank8_percentiles.pyx @@ -0,0 +1,292 @@ +#cython: cdivision=True +#cython: boundscheck=False +#cython: nonecheck=False +#cython: wraparound=False + +import numpy as np +cimport numpy as np +from skimage.filter.rank._core8 cimport _core8, uint8_max, uint8_min + + +# ----------------------------------------------------------------- +# kernels uint8 (SOFT version using percentiles) +# ----------------------------------------------------------------- + + +cdef inline np.uint8_t kernel_autolevel(Py_ssize_t * histo, float pop, + np.uint8_t g, float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + cdef int i, imin, imax, sum, delta + + if pop: + sum = 0 + p1 = 1.0 - p1 + imin = 0 + imax = 255 + + 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 < np.uint8_t > (255 + * (uint8_min(uint8_max(imin, g), imax) - imin) / delta) + else: + return < np.uint8_t > (imax - imin) + else: + return < np.uint8_t > (128) + + +cdef inline np.uint8_t kernel_gradient(Py_ssize_t * histo, float pop, + np.uint8_t g, float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + 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 < np.uint8_t > (imax - imin) + else: + return < np.uint8_t > (0) + + +cdef inline np.uint8_t kernel_mean(Py_ssize_t * histo, float pop, + np.uint8_t g, float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + 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 < np.uint8_t > (1.0 * mean / n) + else: + return < np.uint8_t > (0) + else: + return < np.uint8_t > (0) + + +cdef inline np.uint8_t kernel_mean_substraction(Py_ssize_t * histo, + float pop, + np.uint8_t g, + float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + 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 < np.uint8_t > ((g - (mean / n)) * .5 + 127) + else: + return < np.uint8_t > (0) + else: + return < np.uint8_t > (0) + + +cdef inline np.uint8_t kernel_morph_contr_enh(Py_ssize_t * histo, + float pop, + np.uint8_t g, float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + 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 < np.uint8_t > imax + if g < imin: + return < np.uint8_t > imin + if imax - g < g - imin: + return < np.uint8_t > imax + else: + return < np.uint8_t > imin + else: + return < np.uint8_t > (0) + + +cdef inline np.uint8_t kernel_percentile(Py_ssize_t * histo, float pop, + np.uint8_t g, float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + cdef int i + cdef float sum = 0. + + if pop: + for i in range(256): + sum += histo[i] + if sum >= p0 * pop: + break + + return < np.uint8_t > (i) + else: + return < np.uint8_t > (0) + + +cdef inline np.uint8_t kernel_pop(Py_ssize_t * histo, float pop, + np.uint8_t g, float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + 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 < np.uint8_t > (n) + else: + return < np.uint8_t > (0) + + +cdef inline np.uint8_t kernel_threshold(Py_ssize_t * histo, float pop, + np.uint8_t g, float p0, float p1, + Py_ssize_t s0, Py_ssize_t s1): + cdef int i + cdef float sum = 0. + + if pop: + for i in range(256): + sum += histo[i] + if sum >= p0 * pop: + break + + return < np.uint8_t > (255 * (g >= i)) + else: + return < np.uint8_t > (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.): + """autolevel + """ + _core8(kernel_autolevel, image, selem, mask, out, shift_x, shift_y, p0, p1, + 0, 0) + + +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 + """ + _core8(kernel_gradient, image, selem, mask, out, shift_x, shift_y, p0, p1, + 0, 0) + + +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 + """ + _core8(kernel_mean, image, selem, mask, out, shift_x, shift_y, p0, p1, + 0, 0) + + +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 + """ + _core8(kernel_mean_substraction, image, selem, mask, out, shift_x, shift_y, + p0, p1, 0, 0) + + +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 + """ + _core8(kernel_morph_contr_enh, image, selem, mask, out, shift_x, shift_y, + p0, p1, 0, 0) + + +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 + """ + _core8(kernel_percentile, image, selem, mask, out, shift_x, shift_y, + p0, p1, 0, 0) + + +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] + """ + _core8(kernel_pop, image, selem, mask, out, shift_x, shift_y, p0, p1, + 0, 0) + + +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 + """ + _core8(kernel_threshold, image, selem, mask, out, shift_x, shift_y, p0, p1, + 0, 0) diff --git a/skimage/filter/rank/bilateral_rank.pyx b/skimage/filter/rank/bilateral_rank.pyx new file mode 100644 index 00000000..0d3cdd88 --- /dev/null +++ b/skimage/filter/rank/bilateral_rank.pyx @@ -0,0 +1,188 @@ +"""Approximate bilateral rank filter for local (custom kernel) mean. + +The local histogram is computed using a sliding window similar to the method +described in: + +.. [1] Reference: Huang, T. ,Yang, G. ; Tang, G.. "A fast two-dimensional + median filtering algorithm", IEEE Transactions on Acoustics, Speech and + Signal Processing, Feb 1979. Volume: 27 , Issue: 1, Page(s): 13 - 18. + +Input image can be 8-bit or 16-bit with a value < 4096 (i.e. 12 bit), 8-bit +images are casted in 16-bit the number of histogram bins is determined from the +maximum value present in the image. + +The pixel neighborhood is defined by: + +* the given structuring element +* an interval [g-s0,g+s1] in greylevel around g the processed pixel greylevel + +The kernel is flat (i.e. each pixel belonging to the neighborhood contributes +equally). + +Result image is 16-bit with respect to the input image. + +""" + +import numpy as np +from skimage import img_as_ubyte +from skimage.filter.rank import _crank16_bilateral +from skimage.filter.rank.generic import find_bitdepth + + +__all__ = ['bilateral_mean', 'bilateral_pop'] + + +def _apply(func8, func16, image, selem, out, mask, shift_x, shift_y, s0, s1): + selem = img_as_ubyte(selem) + image = np.ascontiguousarray(image) + + if mask is None: + mask = np.ones(image.shape, dtype=np.uint8) + else: + mask = np.ascontiguousarray(mask) + mask = img_as_ubyte(mask) + + if image is out: + raise NotImplementedError("Cannot perform rank operation in place.") + + if image.dtype == np.uint8: + if func8 is None: + raise TypeError("Not implemented for uint8 image.") + if out is None: + out = np.zeros(image.shape, dtype=np.uint8) + func8(image, selem, shift_x=shift_x, shift_y=shift_y, + mask=mask, out=out, s0=s0, s1=s1) + elif image.dtype == np.uint16: + if func16 is None: + raise TypeError("Not implemented for uint16 image.") + if out is None: + out = np.zeros(image.shape, dtype=np.uint16) + bitdepth = find_bitdepth(image) + if bitdepth > 11: + raise ValueError("Only uint16 <4096 image (12bit) supported.") + func16(image, selem, shift_x=shift_x, shift_y=shift_y, mask=mask, + bitdepth=bitdepth + 1, out=out, s0=s0, s1=s1) + else: + raise TypeError("Only uint8 and uint16 image supported.") + + return out + + +def bilateral_mean(image, selem, out=None, mask=None, shift_x=False, + shift_y=False, s0=10, s1=10): + """Apply a flat kernel bilateral filter. + + This is an edge-preserving and noise reducing denoising filter. It averages + pixels based on their spatial closeness and radiometric similarity. + + Spatial closeness is measured by considering only the local pixel + neighborhood given by a structuring element (selem). + + Radiometric similarity is defined by the greylevel interval [g-s0,g+s1] + where g is the current pixel greylevel. Only pixels belonging to the + structuring element AND having a greylevel inside this interval are + averaged. Return greyscale local bilateral_mean of an image. + + Parameters + ---------- + image : ndarray + Image array (uint8 array or uint16). If image is uint16, as the + algorithm uses max. 12bit histogram, an exception will be raised if + image has a value > 4095 + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray + If None, a new array will be allocated. + mask : ndarray (uint8) + Mask array that defines (>0) area of the image included in the local + neighborhood. If None, the complete image is used (default). + shift_x, shift_y : (int) + Offset added to the structuring element center point. Shift is bounded + to the structuring element sizes (center must be inside the given + structuring element). + s0, s1 : int + define the [s0, s1] interval to be considered for computing the value. + + Returns + ------- + out : uint16 array (uint8 image are casted to uint16) + The result of the local bilateral mean. + + See also + -------- + skimage.filter.denoise_bilateral() for a gaussian bilateral filter. + + Notes + ----- + + * input image can be 8-bit or 16-bit with a value < 4096 (i.e. 12 bit) + + * 8-bit images are casted in 16-bit + + Examples + -------- + >>> from skimage import data + >>> from skimage.morphology import disk + >>> from skimage.filter.rank import bilateral_mean + >>> # Load test image + >>> ima = data.camera() + >>> # bilateral filtering of cameraman image using a flat kernel + >>> bilat_ima = bilateral_mean(ima, disk(20), s0=10,s1=10) + """ + + return _apply(None, _crank16_bilateral.mean, image, selem, out=out, + mask=mask, shift_x=shift_x, shift_y=shift_y, s0=s0, s1=s1) + + +def bilateral_pop(image, selem, out=None, mask=None, shift_x=False, + shift_y=False, s0=10, s1=10): + """Return the number (population) of pixels actually inside the bilateral + neighborhood, i.e. being inside the structuring element AND having a gray + level inside the interval [g-s0, g+s1]. + + Parameters + ---------- + image : ndarray + Image array (uint8 array or uint16). If image is uint16, as the + algorithm uses max. 12bit histogram, an exception will be raised if + image has a value > 4095 + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray + If None, a new array will be allocated. + mask : ndarray (uint8) + Mask array that defines (>0) area of the image included in the local + neighborhood. If None, the complete image is used (default). + shift_x, shift_y : (int) + Offset added to the structuring element center point. Shift is bounded + to the structuring element sizes (center must be inside the given + structuring element). + s0, s1 : int + define the [s0, s1] interval to be considered for computing the value. + + Returns + ------- + out : uint16 array (uint8 image are casted to uint16) + the local number of pixels inside the bilateral neighborhood + + Examples + -------- + >>> # Local mean + >>> from skimage.morphology import square + >>> import skimage.filter.rank as rank + >>> ima8 = 255 * np.array([[0, 0, 0, 0, 0], + ... [0, 1, 1, 1, 0], + ... [0, 1, 1, 1, 0], + ... [0, 1, 1, 1, 0], + ... [0, 0, 0, 0, 0]], dtype=np.uint8) + >>> rank.bilateral_pop(ima8, square(3), s0=10,s1=10) + array([[3, 4, 3, 4, 3], + [4, 4, 6, 4, 4], + [3, 6, 9, 6, 3], + [4, 4, 6, 4, 4], + [3, 4, 3, 4, 3]], dtype=uint16) + + """ + + return _apply(None, _crank16_bilateral.pop, image, selem, out=out, + mask=mask, shift_x=shift_x, shift_y=shift_y, s0=s0, s1=s1) diff --git a/skimage/filter/rank/generic.py b/skimage/filter/rank/generic.py new file mode 100644 index 00000000..94fc3130 --- /dev/null +++ b/skimage/filter/rank/generic.py @@ -0,0 +1,11 @@ +import numpy as np + + +def find_bitdepth(image): + """returns the max bith depth of a uint16 image + """ + umax = np.max(image) + if umax > 2: + return int(np.log2(umax)) + else: + return 1 diff --git a/skimage/filter/rank/percentile_rank.pyx b/skimage/filter/rank/percentile_rank.pyx new file mode 100644 index 00000000..62cc4576 --- /dev/null +++ b/skimage/filter/rank/percentile_rank.pyx @@ -0,0 +1,396 @@ +"""Inferior and superior ranks, provided by the user, are passed to the kernel +function to provide a softer version of the rank filters. E.g. +percentile_autolevel will stretch image levels between percentile [p0, p1] +instead of using [min, max]. It means that isolated bright or dark pixels will +not produce halos. + +The local histogram is computed using a sliding window similar to the method +described in [1]. + +References +========== + +.. [1] Huang, T. ,Yang, G. ; Tang, G.. "A fast two-dimensional + median filtering algorithm", IEEE Transactions on Acoustics, Speech and + Signal Processing, Feb 1979. Volume: 27 , Issue: 1, Page(s): 13 - 18. + +Input image can be 8-bit or 16-bit with a value < 4096 (i.e. 12 bit), for 16-bit +input images, the number of histogram bins is determined from the maximum value +present in the image. + +Result image is 8 or 16-bit with respect to the input image. + +""" + +import numpy as np +from skimage import img_as_ubyte +from skimage.filter.rank.generic import find_bitdepth +from skimage.filter.rank import _crank16_percentiles, _crank8_percentiles + + +__all__ = ['percentile_autolevel', 'percentile_gradient', + 'percentile_mean', 'percentile_mean_substraction', + 'percentile_morph_contr_enh', 'percentile', 'percentile_pop', + 'percentile_threshold'] + + +def _apply(func8, func16, image, selem, out, mask, shift_x, shift_y, p0, p1): + selem = img_as_ubyte(selem) + image = np.ascontiguousarray(image) + + if mask is None: + mask = np.ones(image.shape, dtype=np.uint8) + else: + mask = np.ascontiguousarray(mask) + mask = img_as_ubyte(mask) + + if image is out: + raise NotImplementedError("Cannot perform rank operation in place.") + + if image.dtype == np.uint8: + if func8 is None: + raise TypeError("Not implemented for uint8 image.") + if out is None: + out = np.zeros(image.shape, dtype=np.uint8) + func8(image, selem, shift_x=shift_x, shift_y=shift_y, + mask=mask, out=out, p0=p0, p1=p1) + elif image.dtype == np.uint16: + if func16 is None: + raise TypeError("Not implemented for uint16 image.") + if out is None: + out = np.zeros(image.shape, dtype=np.uint16) + bitdepth = find_bitdepth(image) + if bitdepth > 11: + raise ValueError("Only uint16 <4096 image (12bit) supported.") + func16(image, selem, shift_x=shift_x, shift_y=shift_y, mask=mask, + bitdepth=bitdepth + 1, out=out, p0=p0, p1=p1) + else: + raise TypeError("Only uint8 and uint16 image supported.") + + return out + + +def percentile_autolevel(image, selem, out=None, mask=None, shift_x=False, + shift_y=False, p0=.0, p1=1.): + """Return greyscale local autolevel of an image. + + Autolevel is computed on the given structuring element. Only levels between + percentiles [p0, p1] are used. + + Parameters + ---------- + image : ndarray + Image array (uint8 array or uint16). If image is uint16, as the + algorithm uses max. 12bit histogram, an exception will be raised if + image has a value > 4095. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray + If None, a new array will be allocated. + mask : ndarray (uint8) + Mask array that defines (>0) area of the image included in the local + neighborhood. If None, the complete image is used (default). + shift_x, shift_y : int + Offset added to the structuring element center point. Shift is bounded + to the structuring element sizes (center must be inside the given + structuring element). + p0, p1 : float in [0, ..., 1] + Define the [p0, p1] percentile interval to be considered for computing + the value. + + Returns + ------- + local autolevel : uint8 array or uint16 + The result of the local autolevel. + + """ + + return _apply( + _crank8_percentiles.autolevel, _crank16_percentiles.autolevel, + image, selem, out=out, mask=mask, shift_x=shift_x, + shift_y=shift_y, p0=p0, p1=p1) + + +def percentile_gradient(image, selem, out=None, mask=None, shift_x=False, + shift_y=False, p0=.0, p1=1.): + """Return greyscale local percentile_gradient of an image. + + percentile_gradient is computed on the given structuring element. Only + levels between percentiles [p0, p1] are used. + + Parameters + ---------- + image : ndarray + Image array (uint8 array or uint16). If image is uint16, as the + algorithm uses max. 12bit histogram, an exception will be raised if + image has a value > 4095. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray + If None, a new array will be allocated. + mask : ndarray (uint8) + Mask array that defines (>0) area of the image included in the local + neighborhood. If None, the complete image is used (default). + shift_x, shift_y : int + Offset added to the structuring element center point. Shift is bounded + to the structuring element sizes (center must be inside the given + structuring element). + p0, p1 : float in [0, ..., 1] + Define the [p0, p1] percentile interval to be considered for computing + the value. + + Returns + ------- + local percentile_gradient : uint8 array or uint16 + The result of the local percentile_gradient. + + """ + + return _apply(_crank8_percentiles.gradient, _crank16_percentiles.gradient, + image, selem, out=out, mask=mask, shift_x=shift_x, + shift_y=shift_y, p0=p0, p1=p1) + + +def percentile_mean(image, selem, out=None, mask=None, shift_x=False, + shift_y=False, p0=.0, p1=1.): + """Return greyscale local mean of an image. + + Mean is computed on the given structuring element. Only levels between + percentiles [p0, p1] are used. + + Parameters + ---------- + image : ndarray + Image array (uint8 array or uint16). If image is uint16, as the + algorithm uses max. 12bit histogram, an exception will be raised if + image has a value > 4095. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray + If None, a new array will be allocated. + mask : ndarray (uint8) + Mask array that defines (>0) area of the image included in the local + neighborhood. If None, the complete image is used (default). + shift_x, shift_y : int + Offset added to the structuring element center point. Shift is bounded + to the structuring element sizes (center must be inside the given + structuring element). + p0, p1 : float in [0, ..., 1] + Define the [p0, p1] percentile interval to be considered for computing + the value. + + Returns + ------- + local mean : uint8 array or uint16 + The result of the local mean. + + """ + + return _apply(_crank8_percentiles.mean, _crank16_percentiles.mean, + image, selem, out=out, mask=mask, shift_x=shift_x, + shift_y=shift_y, p0=p0, p1=p1) + + +def percentile_mean_substraction(image, selem, out=None, mask=None, + shift_x=False, shift_y=False, p0=.0, p1=1.): + """Return greyscale local mean_substraction of an image. + + mean_substraction is computed on the given structuring element. Only levels + between percentiles [p0, p1] are used. + + Parameters + ---------- + image : ndarray + Image array (uint8 array or uint16). If image is uint16, as the + algorithm uses max. 12bit histogram, an exception will be raised if + image has a value > 4095. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray + If None, a new array will be allocated. + mask : ndarray (uint8) + Mask array that defines (>0) area of the image included in the local + neighborhood. If None, the complete image is used (default). + shift_x, shift_y : int + Offset added to the structuring element center point. Shift is bounded + to the structuring element sizes (center must be inside the given + structuring element). + p0, p1 : float in [0, ..., 1] + Define the [p0, p1] percentile interval to be considered for computing + the value. + + Returns + ------- + local mean_substraction : uint8 array or uint16 + The result of the local mean_substraction. + + """ + + return _apply(_crank8_percentiles.mean_substraction, + _crank16_percentiles.mean_substraction, + image, selem, out=out, mask=mask, shift_x=shift_x, + shift_y=shift_y, p0=p0, p1=p1) + + +def percentile_morph_contr_enh( + image, selem, out=None, mask=None, shift_x=False, + shift_y=False, p0=.0, p1=1.): + """Return greyscale local morph_contr_enh of an image. + + morph_contr_enh is computed on the given structuring element. Only levels + between percentiles [p0, p1] are used. + + Parameters + ---------- + image : ndarray + Image array (uint8 array or uint16). If image is uint16, as the + algorithm uses max. 12bit histogram, an exception will be raised if + image has a value > 4095. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray + If None, a new array will be allocated. + mask : ndarray (uint8) + Mask array that defines (>0) area of the image included in the local + neighborhood. If None, the complete image is used (default). + shift_x, shift_y : int + Offset added to the structuring element center point. Shift is bounded + to the structuring element sizes (center must be inside the given + structuring element). + p0, p1 : float in [0, ..., 1] + Define the [p0, p1] percentile interval to be considered for computing + the value. + + Returns + ------- + local morph_contr_enh : uint8 array or uint16 + The result of the local morph_contr_enh. + + """ + + return _apply(_crank8_percentiles.morph_contr_enh, + _crank16_percentiles.morph_contr_enh, + image, selem, out=out, mask=mask, shift_x=shift_x, + shift_y=shift_y, p0=p0, p1=p1) + + +def percentile(image, selem, out=None, mask=None, shift_x=False, shift_y=False, + p0=.0, p1=1.): + """Return greyscale local percentile of an image. + + percentile is computed on the given structuring element. Only levels between + percentiles [p0, p1] are used. + + Parameters + ---------- + image : ndarray + Image array (uint8 array or uint16). If image is uint16, as the + algorithm uses max. 12bit histogram, an exception will be raised if + image has a value > 4095. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray + If None, a new array will be allocated. + mask : ndarray (uint8) + Mask array that defines (>0) area of the image included in the local + neighborhood. If None, the complete image is used (default). + shift_x, shift_y : int + Offset added to the structuring element center point. Shift is bounded + to the structuring element sizes (center must be inside the given + structuring element). + p0, p1 : float in [0, ..., 1] + Define the [p0, p1] percentile interval to be considered for computing + the value. + + Returns + ------- + local percentile : uint8 array or uint16 + The result of the local percentile. + + """ + + return _apply(_crank8_percentiles.percentile, + _crank16_percentiles.percentile, + image, selem, out=out, mask=mask, shift_x=shift_x, + shift_y=shift_y, p0=p0, p1=p1) + + +def percentile_pop(image, selem, out=None, mask=None, shift_x=False, + shift_y=False, p0=.0, p1=1.): + """Return greyscale local pop of an image. + + pop is computed on the given structuring element. Only levels between + percentiles [p0, p1] are used. + + Parameters + ---------- + image : ndarray + Image array (uint8 array or uint16). If image is uint16, as the + algorithm uses max. 12bit histogram, an exception will be raised if + image has a value > 4095. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray + If None, a new array will be allocated. + mask : ndarray (uint8) + Mask array that defines (>0) area of the image included in the local + neighborhood. If None, the complete image is used (default). + shift_x, shift_y : int + Offset added to the structuring element center point. Shift is bounded + to the structuring element sizes (center must be inside the given + structuring element). + p0, p1 : float in [0, ..., 1] + Define the [p0, p1] percentile interval to be considered for computing + the value. + + Returns + ------- + local pop : uint8 array or uint16 + The result of the local pop. + + """ + + return _apply(_crank8_percentiles.pop, _crank16_percentiles.pop, + image, selem, out=out, mask=mask, shift_x=shift_x, + shift_y=shift_y, p0=p0, p1=p1) + + +def percentile_threshold(image, selem, out=None, mask=None, shift_x=False, + shift_y=False, p0=.0, p1=1.): + """Return greyscale local threshold of an image. + + threshold is computed on the given structuring element. Only levels between + percentiles [p0, p1] are used. + + Parameters + ---------- + image : ndarray + Image array (uint8 array or uint16). If image is uint16, as the + algorithm uses max. 12bit histogram, an exception will be raised if + image has a value > 4095. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray + If None, a new array will be allocated. + mask : ndarray (uint8) + Mask array that defines (>0) area of the image included in the local + neighborhood. If None, the complete image is used (default). + shift_x, shift_y : int + Offset added to the structuring element center point. Shift is bounded + to the structuring element sizes (center must be inside the given + structuring element). + p0, p1 : float in [0, ..., 1] + Define the [p0, p1] percentile interval to be considered for computing + the value. + + Returns + ------- + local threshold : uint8 array or uint16 + The result of the local threshold. + + """ + + return _apply( + _crank8_percentiles.threshold, _crank16_percentiles.threshold, + image, selem, out=out, mask=mask, shift_x=shift_x, + shift_y=shift_y, p0=p0, p1=p1) diff --git a/skimage/filter/rank/rank.pyx b/skimage/filter/rank/rank.pyx new file mode 100644 index 00000000..be078ddd --- /dev/null +++ b/skimage/filter/rank/rank.pyx @@ -0,0 +1,764 @@ +"""The local histogram is computed using a sliding window similar to the method +described in: + +.. [1] Reference: Huang, T. ,Yang, G. ; Tang, G.. "A fast two-dimensional + median filtering algorithm", IEEE Transactions on Acoustics, Speech and + Signal Processing, Feb 1979. Volume: 27 , Issue: 1, Page(s): 13 - 18. + +Input image can be 8-bit or 16-bit with a value < 4096 (i.e. 12 bit), for 16-bit +input images, the number of histogram bins is determined from the maximum value +present in the image. + +Result image is 8 or 16-bit with respect to the input image. + +""" + +import numpy as np +from skimage import img_as_ubyte +from skimage.filter.rank import _crank8, _crank16 +from skimage.filter.rank.generic import find_bitdepth + + +__all__ = ['autolevel', 'bottomhat', 'equalize', 'gradient', 'maximum', 'mean', + 'meansubstraction', 'median', 'minimum', 'modal', 'morph_contr_enh', + 'pop', 'threshold', 'tophat', 'noise_filter', 'entropy', 'otsu'] + + +def _apply(func8, func16, image, selem, out, mask, shift_x, shift_y): + selem = img_as_ubyte(selem) + image = np.ascontiguousarray(image) + + if mask is None: + mask = np.ones(image.shape, dtype=np.uint8) + else: + mask = np.ascontiguousarray(mask) + mask = img_as_ubyte(mask) + + if image is out: + raise NotImplementedError("Cannot perform rank operation in place.") + + if image.dtype == np.uint8: + if func8 is None: + raise TypeError("Not implemented for uint8 image.") + if out is None: + out = np.zeros(image.shape, dtype=np.uint8) + func8(image, selem, shift_x=shift_x, shift_y=shift_y, + mask=mask, out=out) + elif image.dtype == np.uint16: + if func16 is None: + raise TypeError("Not implemented for uint16 image.") + if out is None: + out = np.zeros(image.shape, dtype=np.uint16) + bitdepth = find_bitdepth(image) + if bitdepth > 11: + raise ValueError("Only uint16 <4096 image (12bit) supported.") + func16(image, selem, shift_x=shift_x, shift_y=shift_y, mask=mask, + bitdepth=bitdepth + 1, out=out) + else: + raise TypeError("Only uint8 and uint16 image supported.") + + return out + + +def autolevel(image, selem, out=None, mask=None, shift_x=False, shift_y=False): + """Autolevel image using local histogram. + + Parameters + ---------- + image : ndarray + Image array (uint8 array or uint16). If image is uint16, the algorithm + uses max. 12bit histogram, an exception will be raised if image has a + value > 4095. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray + If None, a new array will be allocated. + mask : ndarray (uint8) + Mask array that defines (>0) area of the image included in the local + neighborhood. If None, the complete image is used (default). + shift_x, shift_y : int + Offset added to the structuring element center point. Shift is bounded + to the structuring element sizes (center must be inside the given + structuring element). + + Returns + ------- + out : uint8 array or uint16 array (same as input image) + The result of the local autolevel. + + Examples + -------- + >>> from skimage import data + >>> from skimage.morphology import disk + >>> from skimage.filter.rank import autolevel + >>> # Load test image + >>> ima = data.camera() + >>> # Stretch image contrast locally + >>> auto = autolevel(ima, disk(20)) + + """ + + return _apply(_crank8.autolevel, _crank16.autolevel, image, selem, out=out, + mask=mask, shift_x=shift_x, shift_y=shift_y) + + +def bottomhat(image, selem, out=None, mask=None, shift_x=False, shift_y=False): + """Returns greyscale local bottomhat of an image. + + Parameters + ---------- + image : ndarray + Image array (uint8 array or uint16). If image is uint16, the algorithm + uses max. 12bit histogram, an exception will be raised if image has a + value > 4095. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray + If None, a new array will be allocated. + mask : ndarray (uint8) + Mask array that defines (>0) area of the image included in the local + neighborhood. If None, the complete image is used (default). + shift_x, shift_y : int + Offset added to the structuring element center point. Shift is bounded + to the structuring element sizes (center must be inside the given + structuring element). + + Returns + ------- + local bottomhat : uint8 array or uint16 array depending on input image + The result of the local bottomhat. + + """ + + return _apply(_crank8.bottomhat, _crank16.bottomhat, image, selem, out=out, + mask=mask, shift_x=shift_x, shift_y=shift_y) + + +def equalize(image, selem, out=None, mask=None, shift_x=False, shift_y=False): + """Equalize image using local histogram. + + Parameters + ---------- + image : ndarray + Image array (uint8 array or uint16). If image is uint16, the algorithm + uses max. 12bit histogram, an exception will be raised if image has a + value > 4095. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray + If None, a new array will be allocated. + mask : ndarray (uint8) + Mask array that defines (>0) area of the image included in the local + neighborhood. If None, the complete image is used (default). + shift_x, shift_y : int + Offset added to the structuring element center point. Shift is bounded + to the structuring element sizes (center must be inside the given + structuring element). + + Returns + ------- + out : uint8 array or uint16 array (same as input image) + The result of the local equalize. + + Examples + -------- + >>> from skimage import data + >>> from skimage.morphology import disk + >>> from skimage.filter.rank import equalize + >>> # Load test image + >>> ima = data.camera() + >>> # Local equalization + >>> equ = equalize(ima, disk(20)) + + """ + + return _apply(_crank8.equalize, _crank16.equalize, image, selem, out=out, + mask=mask, shift_x=shift_x, shift_y=shift_y) + + +def gradient(image, selem, out=None, mask=None, shift_x=False, shift_y=False): + """Return greyscale local gradient of an image (i.e. local maximum - local + minimum). + + + Parameters + ---------- + image : ndarray + Image array (uint8 array or uint16). If image is uint16, the algorithm + uses max. 12bit histogram, an exception will be raised if image has a + value > 4095. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray + If None, a new array will be allocated. + mask : ndarray (uint8) + Mask array that defines (>0) area of the image included in the local + neighborhood. If None, the complete image is used (default). + shift_x, shift_y : int + Offset added to the structuring element center point. Shift is bounded + to the structuring element sizes (center must be inside the given + structuring element). + + Returns + ------- + out : uint8 array or uint16 array (same as input image) + The local gradient. + + """ + + return _apply(_crank8.gradient, _crank16.gradient, image, selem, out=out, + mask=mask, shift_x=shift_x, shift_y=shift_y) + + +def maximum(image, selem, out=None, mask=None, shift_x=False, shift_y=False): + """Return greyscale local maximum of an image. + + + Parameters + ---------- + image : ndarray + Image array (uint8 array or uint16). If image is uint16, the algorithm + uses max. 12bit histogram, an exception will be raised if image has a + value > 4095. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray + If None, a new array will be allocated. + mask : ndarray (uint8) + Mask array that defines (>0) area of the image included in the local + neighborhood. If None, the complete image is used (default). + shift_x, shift_y : int + Offset added to the structuring element center point. Shift is bounded + to the structuring element sizes (center must be inside the given + structuring element). + + Returns + ------- + out : uint8 array or uint16 array (same as input image) + The local maximum. + + See also + -------- + skimage.morphology.dilation + + Note + ---- + * input image can be 8-bit or 16-bit with a value < 4096 (i.e. 12 bit) + * the lower algorithm complexity makes the rank.maximum() more efficient for + larger images and structuring elements + + """ + + return _apply(_crank8.maximum, _crank16.maximum, image, selem, out=out, + mask=mask, shift_x=shift_x, shift_y=shift_y) + + +def mean(image, selem, out=None, mask=None, shift_x=False, shift_y=False): + """Return greyscale local mean of an image. + + Parameters + ---------- + image : ndarray + Image array (uint8 array or uint16). If image is uint16, the algorithm + uses max. 12bit histogram, an exception will be raised if image has a + value > 4095. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray + If None, a new array will be allocated. + mask : ndarray (uint8) + Mask array that defines (>0) area of the image included in the local + neighborhood. If None, the complete image is used (default). + shift_x, shift_y : int + Offset added to the structuring element center point. Shift is bounded + to the structuring element sizes (center must be inside the given + structuring element). + + Returns + ------- + out : uint8 array or uint16 array (same as input image) + The local mean. + + Examples + -------- + >>> from skimage import data + >>> from skimage.morphology import disk + >>> from skimage.filter.rank import mean + >>> # Load test image + >>> ima = data.camera() + >>> # Local mean + >>> avg = mean(ima, disk(20)) + + """ + + return _apply(_crank8.mean, _crank16.mean, image, selem, out=out, + mask=mask, shift_x=shift_x, shift_y=shift_y) + + +def meansubstraction(image, selem, out=None, mask=None, shift_x=False, + shift_y=False): + """Return image substracted from its local mean. + + Parameters + ---------- + image : ndarray + Image array (uint8 array or uint16). If image is uint16, the algorithm + uses max. 12bit histogram, an exception will be raised if image has a + value > 4095. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray + If None, a new array will be allocated. + mask : ndarray (uint8) + Mask array that defines (>0) area of the image included in the local + neighborhood. If None, the complete image is used (default). + shift_x, shift_y : int + Offset added to the structuring element center point. Shift is bounded + to the structuring element sizes (center must be inside the given + structuring element). + + Returns + ------- + out : uint8 array or uint16 array (same as input image) + The result of the local meansubstraction. + + """ + + return _apply(_crank8.meansubstraction, _crank16.meansubstraction, image, + selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y) + + +def median(image, selem, out=None, mask=None, shift_x=False, shift_y=False): + """Return greyscale local median of an image. + + Parameters + ---------- + image : ndarray + Image array (uint8 array or uint16). If image is uint16, the algorithm + uses max. 12bit histogram, an exception will be raised if image has a + value > 4095. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray + If None, a new array will be allocated. + mask : ndarray (uint8) + Mask array that defines (>0) area of the image included in the local + neighborhood. If None, the complete image is used (default). + shift_x, shift_y : int + Offset added to the structuring element center point. Shift is bounded + to the structuring element sizes (center must be inside the given + structuring element). + + Returns + ------- + out : uint8 array or uint16 array (same as input image) + The local median. + + Examples + -------- + >>> from skimage import data + >>> from skimage.morphology import disk + >>> from skimage.filter.rank import median + >>> # Load test image + >>> ima = data.camera() + >>> # Local mean + >>> avg = median(ima, disk(20)) + + """ + + return _apply(_crank8.median, _crank16.median, image, selem, out=out, + mask=mask, shift_x=shift_x, shift_y=shift_y) + + +def minimum(image, selem, out=None, mask=None, shift_x=False, shift_y=False): + """Return greyscale local minimum of an image. + + Parameters + ---------- + image : ndarray + Image array (uint8 array or uint16). If image is uint16, the algorithm + uses max. 12bit histogram, an exception will be raised if image has a + value > 4095. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray + If None, a new array will be allocated. + mask : ndarray (uint8) + Mask array that defines (>0) area of the image included in the local + neighborhood. If None, the complete image is used (default). + shift_x, shift_y : int + Offset added to the structuring element center point. Shift is bounded + to the structuring element sizes (center must be inside the given + structuring element). + + Returns + ------- + out : uint8 array or uint16 array (same as input image) + The local minimum. + + See also + -------- + skimage.morphology.erosion + + Note + ---- + * input image can be 8-bit or 16-bit with a value < 4096 (i.e. 12 bit) + * the lower algorithm complexity makes the rank.minimum() more efficient + for larger images and structuring elements + + """ + + return _apply(_crank8.minimum, _crank16.minimum, image, selem, out=out, + mask=mask, shift_x=shift_x, shift_y=shift_y) + + +def modal(image, selem, out=None, mask=None, shift_x=False, shift_y=False): + """Return greyscale local mode of an image. + + Parameters + ---------- + image : ndarray + Image array (uint8 array or uint16). If image is uint16, the algorithm + uses max. 12bit histogram, an exception will be raised if image has a + value > 4095. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray + If None, a new array will be allocated. + mask : ndarray (uint8) + Mask array that defines (>0) area of the image included in the local + neighborhood. If None, the complete image is used (default). + shift_x, shift_y : int + Offset added to the structuring element center point. Shift is bounded + to the structuring element sizes (center must be inside the given + structuring element). + + Returns + ------- + out : uint8 array or uint16 array (same as input image) + The local modal. + + """ + + return _apply(_crank8.modal, _crank16.modal, image, selem, out=out, + mask=mask, shift_x=shift_x, shift_y=shift_y) + + +def morph_contr_enh(image, selem, out=None, mask=None, shift_x=False, + shift_y=False): + """Enhance an image replacing each pixel by the local maximum if pixel + greylevel is closest to maximimum than local minimum OR local minimum + otherwise. + + Parameters + ---------- + image : ndarray + Image array (uint8 array or uint16). If image is uint16, the algorithm + uses max. 12bit histogram, an exception will be raised if image has a + value > 4095. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray + If None, a new array will be allocated. + mask : ndarray (uint8) + Mask array that defines (>0) area of the image included in the local + neighborhood. If None, the complete image is used (default). + shift_x, shift_y : int + Offset added to the structuring element center point. Shift is bounded + to the structuring element sizes (center must be inside the given + structuring element). + + Returns + ------- + out : uint8 array or uint16 array (same as input image) + The result of the local morph_contr_enh. + + Examples + -------- + >>> from skimage import data + >>> from skimage.morphology import disk + >>> from skimage.filter.rank import morph_contr_enh + >>> # Load test image + >>> ima = data.camera() + >>> # Local mean + >>> avg = morph_contr_enh(ima, disk(20)) + + """ + + return _apply(_crank8.morph_contr_enh, _crank16.morph_contr_enh, image, + selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y) + + +def pop(image, selem, out=None, mask=None, shift_x=False, shift_y=False): + """Return the number (population) of pixels actually inside the + neighborhood. + + Parameters + ---------- + image : ndarray + Image array (uint8 array or uint16). If image is uint16, the algorithm + uses max. 12bit histogram, an exception will be raised if image has a + value > 4095. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray + If None, a new array will be allocated. + mask : ndarray (uint8) + Mask array that defines (>0) area of the image included in the local + neighborhood. If None, the complete image is used (default). + shift_x, shift_y : int + Offset added to the structuring element center point. Shift is bounded + to the structuring element sizes (center must be inside the given + structuring element). + + Returns + ------- + out : uint8 array or uint16 array (same as input image) + The number of pixels belonging to the neighborhood. + + Examples + -------- + >>> # Local mean + >>> from skimage.morphology import square + >>> import skimage.filter.rank as rank + >>> ima = 255 * np.array([[0, 0, 0, 0, 0], + ... [0, 1, 1, 1, 0], + ... [0, 1, 1, 1, 0], + ... [0, 1, 1, 1, 0], + ... [0, 0, 0, 0, 0]], dtype=np.uint8) + >>> rank.pop(ima, square(3)) + array([[4, 6, 6, 6, 4], + [6, 9, 9, 9, 6], + [6, 9, 9, 9, 6], + [6, 9, 9, 9, 6], + [4, 6, 6, 6, 4]], dtype=uint8) + + """ + + return _apply(_crank8.pop, _crank16.pop, image, selem, out=out, + mask=mask, shift_x=shift_x, shift_y=shift_y) + + +def threshold(image, selem, out=None, mask=None, shift_x=False, shift_y=False): + """Return greyscale local threshold of an image. + + Parameters + ---------- + image : ndarray + Image array (uint8 array or uint16). If image is uint16, the algorithm + uses max. 12bit histogram, an exception will be raised if image has a + value > 4095. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray + If None, a new array will be allocated. + mask : ndarray (uint8) + Mask array that defines (>0) area of the image included in the local + neighborhood. If None, the complete image is used (default). + shift_x, shift_y : int + Offset added to the structuring element center point. Shift is bounded + to the structuring element sizes (center must be inside the given + structuring element). + + Returns + ------- + out : uint8 array or uint16 array (same as input image) + The result of the local threshold. + + Examples + -------- + >>> # Local threshold + >>> from skimage.morphology import square + >>> from skimage.filter.rank import threshold + >>> ima = 255 * np.array([[0, 0, 0, 0, 0], + ... [0, 1, 1, 1, 0], + ... [0, 1, 1, 1, 0], + ... [0, 1, 1, 1, 0], + ... [0, 0, 0, 0, 0]], dtype=np.uint8) + >>> threshold(ima, square(3)) + array([[0, 0, 0, 0, 0], + [0, 1, 1, 1, 0], + [0, 1, 0, 1, 0], + [0, 1, 1, 1, 0], + [0, 0, 0, 0, 0]], dtype=uint8) + + """ + + return _apply(_crank8.threshold, _crank16.threshold, image, selem, out=out, + mask=mask, shift_x=shift_x, shift_y=shift_y) + + +def tophat(image, selem, out=None, mask=None, shift_x=False, shift_y=False): + """Return greyscale local tophat of an image. + + Parameters + ---------- + image : ndarray + Image array (uint8 array or uint16). If image is uint16, the algorithm + uses max. 12bit histogram, an exception will be raised if image has a + value > 4095. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray + If None, a new array will be allocated. + mask : ndarray (uint8) + Mask array that defines (>0) area of the image included in the local + neighborhood. If None, the complete image is used (default). + shift_x, shift_y : int + Offset added to the structuring element center point. Shift is bounded + to the structuring element sizes (center must be inside the given + structuring element). + + Returns + ------- + out : uint8 array or uint16 array (same as input image) + The image tophat. + + """ + + return _apply(_crank8.tophat, _crank16.tophat, image, selem, out=out, + mask=mask, shift_x=shift_x, shift_y=shift_y) + + +def noise_filter(image, selem, out=None, mask=None, shift_x=False, + shift_y=False): + """Returns the noise feature as described in [Hashimoto12]_ + + Parameters + ---------- + image : ndarray + Image array (uint8 array or uint16). If image is uint16, the algorithm + uses max. 12bit histogram, an exception will be raised if image has a + value > 4095. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray + If None, a new array will be allocated. + mask : ndarray (uint8) + Mask array that defines (>0) area of the image included in the local + neighborhood. If None, the complete image is used (default). + shift_x, shift_y : int + Offset added to the structuring element center point. Shift is bounded + to the structuring element sizes (center must be inside the given + structuring element). + + References + ---------- + .. [Hashimoto12] N. Hashimoto et al. Referenceless image quality evaluation + for whole slide imaging. J Pathol Inform 2012;3:9. + + Returns + ------- + out : uint8 array or uint16 array (same as input image) + The image noise . + + """ + + # ensure that the central pixel in the structuring element is empty + centre_r = int(selem.shape[0] / 2) + shift_y + centre_c = int(selem.shape[1] / 2) + shift_x + # make a local copy + selem_cpy = selem.copy() + selem_cpy[centre_r, centre_c] = 0 + + return _apply(_crank8.noise_filter, None, image, selem_cpy, out=out, + mask=mask, shift_x=shift_x, shift_y=shift_y) + + +def entropy(image, selem, out=None, mask=None, shift_x=False, shift_y=False): + """Returns the entropy [wiki_entropy]_ computed locally. Entropy is computed + using base 2 logarithm i.e. the filter returns the minimum number of + bits needed to encode local greylevel distribution. + + Parameters + ---------- + image : ndarray + Image array (uint8 array or uint16). If image is uint16, the algorithm + uses max. 12bit histogram, an exception will be raised if image has a + value > 4095. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray + If None, a new array will be allocated. + mask : ndarray (uint8) + Mask array that defines (>0) area of the image included in the local + neighborhood. If None, the complete image is used (default). + shift_x, shift_y : int + Offset added to the structuring element center point. Shift is bounded + to the structuring element sizes (center must be inside the given + structuring element). + + Returns + ------- + out : uint8 array or uint16 array (same as input image) + entropy x10 (uint8 images) and entropy x1000 (uint16 images) + + References + ---------- + .. [wiki_entropy] http://en.wikipedia.org/wiki/Entropy_(information_theory) + + Examples + -------- + >>> # Local entropy + >>> from skimage import data + >>> from skimage.filter.rank import entropy + >>> from skimage.morphology import disk + >>> # defining a 8- and a 16-bit test images + >>> a8 = data.camera() + >>> a16 = data.camera().astype(np.uint16) * 4 + >>> # pixel values contain 10x the local entropy + >>> ent8 = entropy(a8, disk(5)) + >>> # pixel values contain 1000x the local entropy + >>> ent16 = entropy(a16, disk(5)) + + """ + + return _apply(_crank8.entropy, _crank16.entropy, image, selem, out=out, + mask=mask, shift_x=shift_x, shift_y=shift_y) + + +def otsu(image, selem, out=None, mask=None, shift_x=False, shift_y=False): + """Returns the Otsu's threshold value for each pixel. + + Parameters + ---------- + image : ndarray + Image array (uint8 array or uint16). If image is uint16, the algorithm + uses max. 12bit histogram, an exception will be raised if image has a + value > 4095. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray + If None, a new array will be allocated. + mask : ndarray (uint8) + Mask array that defines (>0) area of the image included in the local + neighborhood. If None, the complete image is used (default). + shift_x, shift_y : int + Offset added to the structuring element center point. Shift is bounded + to the structuring element sizes (center must be inside the given + structuring element). + + Returns + ------- + out : uint8 array or uint16 array (same as input image) + Otsu's threshold values + + References + ---------- + .. [otsu] http://en.wikipedia.org/wiki/Otsu's_method + + Examples + -------- + >>> # Local entropy + >>> from skimage import data + >>> from skimage.filter.rank import otsu + >>> from skimage.morphology import disk + >>> # defining a 8- and a 16-bit test images + >>> a8 = data.camera() + >>> loc_otsu = otsu(a8, disk(5)) + >>> thresh_image = a8 >= loc_otsu + + """ + + return _apply(_crank8.otsu, None, image, selem, out=out, + mask=mask, shift_x=shift_x, shift_y=shift_y) diff --git a/skimage/filter/rank/tests/test_rank.py b/skimage/filter/rank/tests/test_rank.py new file mode 100644 index 00000000..43ff030f --- /dev/null +++ b/skimage/filter/rank/tests/test_rank.py @@ -0,0 +1,380 @@ +import numpy as np +from numpy.testing import run_module_suite, assert_array_equal, assert_raises + +from skimage import data +from skimage.morphology import cmorph, disk +from skimage.filter import rank + + +def test_random_sizes(): + # make sure the size is not a problem + + niter = 10 + elem = np.array([[1, 1, 1], [1, 1, 1], [1, 1, 1]], dtype=np.uint8) + for m, n in np.random.random_integers(1, 100, size=(10, 2)): + mask = np.ones((m, n), dtype=np.uint8) + + image8 = np.ones((m, n), dtype=np.uint8) + out8 = np.empty_like(image8) + rank.mean(image=image8, selem=elem, mask=mask, out=out8, + shift_x=0, shift_y=0) + assert_array_equal(image8.shape, out8.shape) + rank.mean(image=image8, selem=elem, mask=mask, out=out8, + shift_x=+1, shift_y=+1) + assert_array_equal(image8.shape, out8.shape) + + image16 = np.ones((m, n), dtype=np.uint16) + out16 = np.empty_like(image8, dtype=np.uint16) + rank.mean(image=image16, selem=elem, mask=mask, out=out16, + shift_x=0, shift_y=0) + assert_array_equal(image16.shape, out16.shape) + rank.mean(image=image16, selem=elem, mask=mask, out=out16, + shift_x=+1, shift_y=+1) + assert_array_equal(image16.shape, out16.shape) + + rank.percentile_mean(image=image16, mask=mask, out=out16, + selem=elem, shift_x=0, shift_y=0, p0=.1, p1=.9) + assert_array_equal(image16.shape, out16.shape) + rank.percentile_mean(image=image16, mask=mask, out=out16, + selem=elem, shift_x=+1, shift_y=+1, p0=.1, p1=.9) + assert_array_equal(image16.shape, out16.shape) + + +def test_compare_with_cmorph_dilate(): + # compare the result of maximum filter with dilate + + image = (np.random.random((100, 100)) * 256).astype(np.uint8) + out = np.empty_like(image) + mask = np.ones(image.shape, dtype=np.uint8) + + for r in range(1, 20, 1): + elem = np.ones((r, r), dtype=np.uint8) + rank.maximum(image=image, selem=elem, out=out, mask=mask) + cm = cmorph.dilate(image=image, selem=elem) + assert_array_equal(out, cm) + + +def test_compare_with_cmorph_erode(): + # compare the result of maximum filter with erode + + image = (np.random.random((100, 100)) * 256).astype(np.uint8) + out = np.empty_like(image) + mask = np.ones(image.shape, dtype=np.uint8) + + for r in range(1, 20, 1): + elem = np.ones((r, r), dtype=np.uint8) + rank.minimum(image=image, selem=elem, out=out, mask=mask) + cm = cmorph.erode(image=image, selem=elem) + assert_array_equal(out, cm) + + +def test_bitdepth(): + # test the different bit depth for rank16 + + elem = np.ones((3, 3), dtype=np.uint8) + out = np.empty((100, 100), dtype=np.uint16) + mask = np.ones((100, 100), dtype=np.uint8) + + for i in range(5): + image = np.ones((100, 100), dtype=np.uint16) * 255 * 2 ** i + r = rank.percentile_mean(image=image, selem=elem, mask=mask, + out=out, shift_x=0, shift_y=0, p0=.1, p1=.9) + + +def test_population(): + # check the number of valid pixels in the neighborhood + + image = np.zeros((5, 5), dtype=np.uint8) + elem = np.ones((3, 3), dtype=np.uint8) + out = np.empty_like(image) + mask = np.ones(image.shape, dtype=np.uint8) + + rank.pop(image=image, selem=elem, out=out, mask=mask) + r = np.array([[4, 6, 6, 6, 4], + [6, 9, 9, 9, 6], + [6, 9, 9, 9, 6], + [6, 9, 9, 9, 6], + [4, 6, 6, 6, 4]]) + assert_array_equal(r, out) + + +def test_structuring_element8(): + # check the output for a custom structuring element + + r = np.array([[0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0], + [0, 0, 255, 0, 0, 0], + [0, 0, 255, 255, 255, 0], + [0, 0, 0, 255, 255, 0], + [0, 0, 0, 0, 0, 0]]) + + # 8-bit + image = np.zeros((6, 6), dtype=np.uint8) + image[2, 2] = 255 + elem = np.asarray([[1, 1, 0], [1, 1, 1], [0, 0, 1]], dtype=np.uint8) + out = np.empty_like(image) + mask = np.ones(image.shape, dtype=np.uint8) + + rank.maximum(image=image, selem=elem, out=out, mask=mask, + shift_x=1, shift_y=1) + assert_array_equal(r, out) + + # 16-bit + image = np.zeros((6, 6), dtype=np.uint16) + image[2, 2] = 255 + out = np.empty_like(image) + + rank.maximum(image=image, selem=elem, out=out, mask=mask, + shift_x=1, shift_y=1) + assert_array_equal(r, out) + + +def test_fail_on_bitdepth(): + # should fail because data bitdepth is too high for the function + + image = np.ones((100, 100), dtype=np.uint16) * 2 ** 12 + elem = np.ones((3, 3), dtype=np.uint8) + out = np.empty_like(image) + mask = np.ones(image.shape, dtype=np.uint8) + assert_raises(ValueError, rank.percentile_mean, image=image, + selem=elem, out=out, mask=mask, shift_x=0, shift_y=0) + + +def test_pass_on_bitdepth(): + # should pass because data bitdepth is not too high for the function + + image = np.ones((100, 100), dtype=np.uint16) * 2 ** 11 + elem = np.ones((3, 3), dtype=np.uint8) + out = np.empty_like(image) + mask = np.ones(image.shape, dtype=np.uint8) + + +def test_inplace_output(): + # rank filters are not supposed to filter inplace + + selem = disk(20) + image = (np.random.random((500, 500)) * 256).astype(np.uint8) + out = image + assert_raises(NotImplementedError, rank.mean, image, selem, out=out) + + +def test_compare_autolevels(): + # compare autolevel and percentile autolevel with p0=0.0 and p1=1.0 + # should returns the same arrays + + image = data.camera() + + selem = disk(20) + loc_autolevel = rank.autolevel(image, selem=selem) + loc_perc_autolevel = rank.percentile_autolevel(image, selem=selem, + p0=.0, p1=1.) + + assert_array_equal(loc_autolevel, loc_perc_autolevel) + + +def test_compare_autolevels_16bit(): + # compare autolevel(16-bit) and percentile autolevel(16-bit) with p0=0.0 and + # p1=1.0 should returns the same arrays + + image = data.camera().astype(np.uint16) * 4 + + selem = disk(20) + loc_autolevel = rank.autolevel(image, selem=selem) + loc_perc_autolevel = rank.percentile_autolevel(image, selem=selem, + p0=.0, p1=1.) + + assert_array_equal(loc_autolevel, loc_perc_autolevel) + + +def test_compare_8bit_vs_16bit(): + # filters applied on 8-bit image ore 16-bit image (having only real 8-bit of + # dynamic) should be identical + + image8 = data.camera() + image16 = image8.astype(np.uint16) + assert_array_equal(image8, image16) + + methods = ['autolevel', 'bottomhat', 'equalize', 'gradient', 'maximum', + 'mean', 'meansubstraction', 'median', 'minimum', 'modal', + 'morph_contr_enh', 'pop', 'threshold', 'tophat'] + + for method in methods: + func = getattr(rank, method) + f8 = func(image8, disk(3)) + f16 = func(image16, disk(3)) + assert_array_equal(f8, f16) + + +def test_trivial_selem8(): + # check that min, max and mean returns identity if structuring element + # contains only central pixel + + image = np.zeros((5, 5), dtype=np.uint8) + out = np.zeros_like(image) + mask = np.ones_like(image, dtype=np.uint8) + image[2, 2] = 255 + image[2, 3] = 128 + image[1, 2] = 16 + + elem = np.array([[0, 0, 0], [0, 1, 0], [0, 0, 0]], dtype=np.uint8) + rank.mean(image=image, selem=elem, out=out, mask=mask, + shift_x=0, shift_y=0) + assert_array_equal(image, out) + rank.minimum(image=image, selem=elem, out=out, mask=mask, + shift_x=0, shift_y=0) + assert_array_equal(image, out) + rank.maximum(image=image, selem=elem, out=out, mask=mask, + shift_x=0, shift_y=0) + assert_array_equal(image, out) + + +def test_trivial_selem16(): + # check that min, max and mean returns identity if structuring element + # contains only central pixel + + image = np.zeros((5, 5), dtype=np.uint16) + out = np.zeros_like(image) + mask = np.ones_like(image, dtype=np.uint8) + image[2, 2] = 255 + image[2, 3] = 128 + image[1, 2] = 16 + + elem = np.array([[0, 0, 0], [0, 1, 0], [0, 0, 0]], dtype=np.uint8) + rank.mean(image=image, selem=elem, out=out, mask=mask, + shift_x=0, shift_y=0) + assert_array_equal(image, out) + rank.minimum(image=image, selem=elem, out=out, mask=mask, + shift_x=0, shift_y=0) + assert_array_equal(image, out) + rank.maximum(image=image, selem=elem, out=out, mask=mask, + shift_x=0, shift_y=0) + assert_array_equal(image, out) + + +def test_smallest_selem8(): + # check that min, max and mean returns identity if structuring element + # contains only central pixel + + image = np.zeros((5, 5), dtype=np.uint8) + out = np.zeros_like(image) + mask = np.ones_like(image, dtype=np.uint8) + image[2, 2] = 255 + image[2, 3] = 128 + image[1, 2] = 16 + + elem = np.array([[1]], dtype=np.uint8) + rank.mean(image=image, selem=elem, out=out, mask=mask, + shift_x=0, shift_y=0) + assert_array_equal(image, out) + rank.minimum(image=image, selem=elem, out=out, mask=mask, + shift_x=0, shift_y=0) + assert_array_equal(image, out) + rank.maximum(image=image, selem=elem, out=out, mask=mask, + shift_x=0, shift_y=0) + assert_array_equal(image, out) + + +def test_smallest_selem16(): + # check that min, max and mean returns identity if structuring element + # contains only central pixel + + image = np.zeros((5, 5), dtype=np.uint16) + out = np.zeros_like(image) + mask = np.ones_like(image, dtype=np.uint8) + image[2, 2] = 255 + image[2, 3] = 128 + image[1, 2] = 16 + + elem = np.array([[1]], dtype=np.uint8) + rank.mean(image=image, selem=elem, out=out, mask=mask, + shift_x=0, shift_y=0) + assert_array_equal(image, out) + rank.minimum(image=image, selem=elem, out=out, mask=mask, + shift_x=0, shift_y=0) + assert_array_equal(image, out) + rank.maximum(image=image, selem=elem, out=out, mask=mask, + shift_x=0, shift_y=0) + assert_array_equal(image, out) + + +def test_empty_selem(): + # check that min, max and mean returns zeros if structuring element is empty + + image = np.zeros((5, 5), dtype=np.uint16) + out = np.zeros_like(image) + mask = np.ones_like(image, dtype=np.uint8) + res = np.zeros_like(image) + image[2, 2] = 255 + image[2, 3] = 128 + image[1, 2] = 16 + + elem = np.array([[0, 0, 0], [0, 0, 0]], dtype=np.uint8) + + rank.mean(image=image, selem=elem, out=out, mask=mask, + shift_x=0, shift_y=0) + assert_array_equal(res, out) + rank.minimum(image=image, selem=elem, out=out, mask=mask, + shift_x=0, shift_y=0) + assert_array_equal(res, out) + rank.maximum(image=image, selem=elem, out=out, mask=mask, + shift_x=0, shift_y=0) + assert_array_equal(res, out) + + +def test_otsu(): + # test the local Otsu segmentation on a synthetic image + # (left to right ramp * sinus) + + test = np.tile( + [128, 145, 103, 127, 165, 83, 127, 185, 63, 127, 205, 43, + 127, 225, 23, 127], + (16, 1)) + test = test.astype(np.uint8) + res = np.tile([1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1], + (16, 1)) + selem = np.ones((6, 6), dtype=np.uint8) + th = 1 * (test >= rank.otsu(test, selem)) + assert_array_equal(th, res) + + +def test_entropy(): + # verify that entropy is coherent with bitdepth of the input data + + selem = np.ones((16, 16), dtype=np.uint8) + # 1 bit per pixel + data = np.tile(np.asarray([0, 1]), (100, 100)).astype(np.uint8) + assert(np.max(rank.entropy(data, selem)) == 10) + + # 2 bit per pixel + data = np.tile(np.asarray([[0, 1], [2, 3]]), (10, 10)).astype(np.uint8) + assert(np.max(rank.entropy(data, selem)) == 20) + + # 3 bit per pixel + data = np.tile( + np.asarray([[0, 1, 2, 3], [4, 5, 6, 7]]), (10, 10)).astype(np.uint8) + assert(np.max(rank.entropy(data, selem)) == 30) + + # 4 bit per pixel + data = np.tile( + np.reshape(np.arange(16), (4, 4)), (10, 10)).astype(np.uint8) + assert(np.max(rank.entropy(data, selem)) == 40) + + # 6 bit per pixel + data = np.tile( + np.reshape(np.arange(64), (8, 8)), (10, 10)).astype(np.uint8) + assert(np.max(rank.entropy(data, selem)) == 60) + + # 8-bit per pixel + data = np.tile( + np.reshape(np.arange(256), (16, 16)), (10, 10)).astype(np.uint8) + assert(np.max(rank.entropy(data, selem)) == 80) + + # 12 bit per pixel + selem = np.ones((64, 64), dtype=np.uint8) + data = np.tile( + np.reshape(np.arange(4096), (64, 64)), (2, 2)).astype(np.uint16) + assert(np.max(rank.entropy(data, selem)) == 12000) + + +if __name__ == "__main__": + run_module_suite() diff --git a/skimage/filter/setup.py b/skimage/filter/setup.py index b996055b..56c1e9e5 100644 --- a/skimage/filter/setup.py +++ b/skimage/filter/setup.py @@ -14,11 +14,47 @@ def configuration(parent_package='', top_path=None): cython(['_ctmf.pyx'], working_path=base_path) cython(['_denoise.pyx'], working_path=base_path) + cython(['rank/_core8.pyx'], working_path=base_path) + cython(['rank/_core16.pyx'], working_path=base_path) + cython(['rank/_crank8.pyx'], working_path=base_path) + cython(['rank/_crank8_percentiles.pyx'], working_path=base_path) + cython(['rank/_crank16.pyx'], working_path=base_path) + cython(['rank/_crank16_percentiles.pyx'], working_path=base_path) + cython(['rank/_crank16_bilateral.pyx'], working_path=base_path) + cython(['rank/rank.pyx'], working_path=base_path) + cython(['rank/percentile_rank.pyx'], working_path=base_path) + cython(['rank/bilateral_rank.pyx'], working_path=base_path) config.add_extension('_ctmf', sources=['_ctmf.c'], - include_dirs=[get_numpy_include_dirs()]) + include_dirs=[get_numpy_include_dirs()]) config.add_extension('_denoise', sources=['_denoise.c'], - include_dirs=[get_numpy_include_dirs(), '../_shared']) + include_dirs=[get_numpy_include_dirs(), '../_shared']) + config.add_extension('rank/_core8', sources=['rank/_core8.c'], + include_dirs=[get_numpy_include_dirs()]) + config.add_extension('rank/_core16', sources=['rank/_core16.c'], + include_dirs=[get_numpy_include_dirs()]) + config.add_extension('rank/_crank8', sources=['rank/_crank8.c'], + include_dirs=[get_numpy_include_dirs()]) + config.add_extension( + 'rank/_crank8_percentiles', sources=['rank/_crank8_percentiles.c'], + include_dirs=[get_numpy_include_dirs()]) + config.add_extension('rank/_crank16', sources=['rank/_crank16.c'], + include_dirs=[get_numpy_include_dirs()]) + config.add_extension( + 'rank/_crank16_percentiles', sources=['rank/_crank16_percentiles.c'], + include_dirs=[get_numpy_include_dirs()]) + config.add_extension( + 'rank/_crank16_bilateral', sources=['rank/_crank16_bilateral.c'], + include_dirs=[get_numpy_include_dirs()]) + config.add_extension( + 'rank/rank', sources=['rank/rank.c'], + include_dirs=[get_numpy_include_dirs()]) + config.add_extension( + 'rank/percentile_rank', sources=['rank/percentile_rank.c'], + include_dirs=[get_numpy_include_dirs()]) + config.add_extension( + 'rank/bilateral_rank', sources=['rank/bilateral_rank.c'], + include_dirs=[get_numpy_include_dirs()]) return config