diff --git a/RELEASE.txt b/RELEASE.txt index 80bf4cab..d29143aa 100644 --- a/RELEASE.txt +++ b/RELEASE.txt @@ -1,6 +1,8 @@ How to make a new release of ``skimage`` ======================================== +- Check ``TODO.txt`` for any outstanding tasks. + - Update release notes. - To show a list contributors, run ``doc/release/contributors.sh ``, diff --git a/TODO.txt b/TODO.txt new file mode 100644 index 00000000..2067d036 --- /dev/null +++ b/TODO.txt @@ -0,0 +1,14 @@ +Version 0.10 +------------ +* Remove deprecated functions: + - ``skimage.filter.rank.*`` +* Remove deprecated parameter ``epsilon`` of ``skimage.viewer.LineProfile`` + +Version 0.9 +----------- +* Remove deprecated functions + - ``skimage.filter.denoise_tv_chambolle`` + - ``skimage.morphology.is_local_maximum`` + - ``skimage.transform.hough`` + - ``skimage.transform.probabilistic_hough`` + - ``skimage.transform.hough_peaks`` diff --git a/bento.info b/bento.info index 9e3bd79b..648445d6 100644 --- a/bento.info +++ b/bento.info @@ -126,33 +126,18 @@ Library: Extension: skimage._shared.geometry Sources: skimage/_shared/geometry.pyx - Extension: skimage.filter.rank._core16 + Extension: skimage.filter.rank.generic_cy Sources: - skimage/filter/rank/_core16.pyx - Extension: skimage.filter.rank._crank8 + skimage/filter/rank/generic_cy.pyx + Extension: skimage.filter.rank.percentile_cy Sources: - skimage/filter/rank/_crank8.pyx - Extension: skimage.filter.rank._crank16 + skimage/filter/rank/percentile_cy.pyx + Extension: skimage.filter.rank.core_cy Sources: - skimage/filter/rank/_crank16.pyx - Extension: skimage.filter.rank._core8 + skimage/filter/rank/core_cy.pyx + Extension: skimage.filter.rank.bilateral_cy Sources: - skimage/filter/rank/_core8.pyx - Extension: skimage.filter.rank.bilateral_rank - Sources: - skimage/filter/rank/bilateral_rank.pyx - Extension: skimage.filter.rank._crank16_percentiles - Sources: - skimage/filter/rank/_crank16_percentiles.pyx - Extension: skimage.filter.rank.percentile_rank - Sources: - skimage/filter/rank/percentile_rank.pyx - Extension: skimage.filter.rank._crank8_percentiles - Sources: - skimage/filter/rank/_crank8_percentiles.pyx - Extension: skimage.filter.rank._crank16_bilateral - Sources: - skimage/filter/rank/_crank16_bilateral.pyx + skimage/filter/rank/bilateral_cy.pyx Executable: skivi Module: skimage.scripts.skivi diff --git a/doc/examples/applications/plot_rank_filters.py b/doc/examples/applications/plot_rank_filters.py index 31284226..fb69cd40 100644 --- a/doc/examples/applications/plot_rank_filters.py +++ b/doc/examples/applications/plot_rank_filters.py @@ -3,11 +3,11 @@ Rank filters ============ -Rank filters are non-linear filters using the local greylevels ordering to +Rank filters are non-linear filters using the local gray-level 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. +local gray-level histogram is computed on the neighborhood of a pixel (defined +by a 2-D 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: @@ -26,11 +26,9 @@ Rank filters can be used for several purposes such as: 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`. +In this example, we will see how to filter a gray-level image using some of the +linear and non-linear filters available in skimage. We use the `camera` image +from `skimage.data` for all comparisons. .. [1] Pierre Soille, On morphological operators based on rank filters, Pattern Recognition 35 (2002) 527-535. @@ -42,16 +40,16 @@ import matplotlib.pyplot as plt from skimage import data -ima = data.camera() -hist = np.histogram(ima, bins=np.arange(0, 256)) +noisy_image = data.camera() +hist = np.histogram(noisy_image, 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.imshow(noisy_image, 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') +plt.title('Histogram of grey values') """ @@ -65,50 +63,56 @@ randomly set to 0. The **median** filter is applied to remove the noise. .. note:: - there are different implementations of median filter : + 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]) +noise = np.random.random(noisy_image.shape) +noisy_image = data.camera() +noisy_image[noise > 0.99] = 255 +noisy_image[noise < 0.01] = 0 + +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.imshow(noisy_image, vmin=0, vmax=255) +plt.title('Noisy image') +plt.axis('off') + plt.subplot(2, 2, 2) -plt.imshow(lo, cmap=plt.cm.gray, vmin=0, vmax=255) -plt.xlabel('median $r=1$') +plt.imshow(median(noisy_image, disk(1)), vmin=0, vmax=255) +plt.title('Median $r=1$') +plt.axis('off') + plt.subplot(2, 2, 3) -plt.imshow(hi, cmap=plt.cm.gray, vmin=0, vmax=255) -plt.xlabel('median $r=5$') +plt.imshow(median(noisy_image, disk(5)), vmin=0, vmax=255) +plt.title('Median $r=5$') +plt.axis('off') + plt.subplot(2, 2, 4) -plt.imshow(ext, cmap=plt.cm.gray, vmin=0, vmax=255) -plt.xlabel('median $r=20$') +plt.imshow(median(noisy_image, disk(20)), vmin=0, vmax=255) +plt.title('Median $r=20$') +plt.axis('off') """ .. 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. +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 bigger sizes are filtered as well, such as the camera tripod. The +median filter is often used for noise removal because borders are preserved and +e.g. salt and pepper noise typically does not distort the gray-level. Image smoothing ================ -The example hereunder shows how a local **mean** smoothes the camera man image. +The example hereunder shows how a local **mean** filter smooths the camera man +image. """ @@ -116,13 +120,17 @@ from skimage.filter.rank import mean fig = plt.figure(figsize=[10, 7]) -loc_mean = mean(nima, disk(10)) +loc_mean = mean(noisy_image, disk(10)) + plt.subplot(1, 2, 1) -plt.imshow(ima, cmap=plt.cm.gray, vmin=0, vmax=255) -plt.xlabel('original') +plt.imshow(noisy_image, vmin=0, vmax=255) +plt.title('Original') +plt.axis('off') + plt.subplot(1, 2, 2) -plt.imshow(loc_mean, cmap=plt.cm.gray, vmin=0, vmax=255) -plt.xlabel('local mean $r=10$') +plt.imshow(loc_mean, vmin=0, vmax=255) +plt.title('Local mean $r=10$') +plt.axis('off') """ @@ -130,35 +138,42 @@ plt.xlabel('local mean $r=10$') 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 +that restricts the local neighborhood to pixel having a gray-level similar to the central one. .. note:: - a different implementation is available for color images in + A different implementation is available for color images in `skimage.filter.denoise_bilateral`. """ from skimage.filter.rank import bilateral_mean -ima = data.camera() +noisy_image = data.camera() selem = disk(10) -bilat = bilateral_mean(ima.astype(np.uint16), disk(20), s0=10, s1=10) +bilat = bilateral_mean(noisy_image.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.imshow(noisy_image, cmap=plt.cm.gray) +plt.title('Original') +plt.axis('off') + plt.subplot(2, 2, 3) plt.imshow(bilat, cmap=plt.cm.gray) -plt.xlabel('bilateral mean') +plt.title('Bilateral mean') +plt.axis('off') + plt.subplot(2, 2, 2) -plt.imshow(ima[200:350, 350:450], cmap=plt.cm.gray) +plt.imshow(noisy_image[200:350, 350:450], cmap=plt.cm.gray) +plt.axis('off') + plt.subplot(2, 2, 4) plt.imshow(bilat[200:350, 350:450], cmap=plt.cm.gray) +plt.axis('off') """ @@ -175,7 +190,7 @@ 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. +equalization emphasizes every local gray-level variations. .. [2] http://en.wikipedia.org/wiki/Histogram_equalization .. [3] http://en.wikipedia.org/wiki/Adaptive_histogram_equalization @@ -185,101 +200,113 @@ equalization emphasizes every local greylevel variations. from skimage import exposure from skimage.filter import rank -ima = data.camera() +noisy_image = data.camera() + # equalize globally and locally -glob = exposure.equalize(ima) * 255 -loc = rank.equalize(ima, disk(20)) +glob = exposure.equalize(noisy_image) * 255 +loc = rank.equalize(noisy_image, disk(20)) # extract histogram for each image -hist = np.histogram(ima, bins=np.arange(0, 256)) +hist = np.histogram(noisy_image, 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.imshow(noisy_image, interpolation='nearest') plt.axis('off') + plt.subplot(322) plt.plot(hist[1][:-1], hist[0], lw=2) -plt.title('histogram of grey values') +plt.title('Histogram of gray values') + plt.subplot(323) -plt.imshow(glob, cmap=plt.cm.gray, interpolation='nearest') +plt.imshow(glob, 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.title('Histogram of gray values') + plt.subplot(325) -plt.imshow(loc, cmap=plt.cm.gray, interpolation='nearest') +plt.imshow(loc, 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') +plt.title('Histogram of gray 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. +Another way to maximize the number of gray-levels used for an image is to apply +a local auto-leveling, i.e. the gray-value of a pixel is proportionally +remapped between local minimum and local maximum. -The following example shows how local autolevel enhances the camara man picture. +The following example shows how local auto-level enhances the camara man +picture. """ from skimage.filter.rank import autolevel -ima = data.camera() +noisy_image = data.camera() selem = disk(10) -auto = autolevel(ima.astype(np.uint16), disk(20)) +auto = autolevel(noisy_image.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.imshow(noisy_image, cmap=plt.cm.gray) +plt.title('Original') +plt.axis('off') + plt.subplot(1, 2, 2) plt.imshow(auto, cmap=plt.cm.gray) -plt.xlabel('local autolevel') +plt.title('Local autolevel') +plt.axis('off') """ .. 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. +This filter is very sensitive to local outliers, see the little white spot in +the left part of the sky. 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 auto-level 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 auto-level +result. """ -from skimage.filter.rank import percentile_autolevel +from skimage.filter.rank import autolevel_percentile 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) +loc_perc_autolevel0 = autolevel_percentile(image, selem=selem, p0=.00, p1=1.0) +loc_perc_autolevel1 = autolevel_percentile(image, selem=selem, p0=.01, p1=.99) +loc_perc_autolevel2 = autolevel_percentile(image, selem=selem, p0=.05, p1=.95) +loc_perc_autolevel3 = autolevel_percentile(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') +ax0.set_title('Original / auto-level') ax1.imshow( np.hstack((loc_perc_autolevel0, loc_perc_autolevel1)), vmin=0, vmax=255) -ax1.set_title('percentile autolevel 0%,1%') +ax1.set_title('Percentile auto-level 0%,1%') ax2.imshow( np.hstack((loc_perc_autolevel2, loc_perc_autolevel3)), vmin=0, vmax=255) -ax2.set_title('percentile autolevel 5% and 10%') +ax2.set_title('Percentile auto-level 5% and 10%') for ax in axes: ax.axis('off') @@ -289,29 +316,35 @@ for ax in axes: .. 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. +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 +from skimage.filter.rank import enhance_contrast -ima = data.camera() +noisy_image = data.camera() -enh = morph_contr_enh(ima, disk(5)) +enh = enhance_contrast(noisy_image, 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.imshow(noisy_image, cmap=plt.cm.gray) +plt.title('Original') +plt.axis('off') + plt.subplot(2, 2, 3) plt.imshow(enh, cmap=plt.cm.gray) -plt.xlabel('local morphlogical contrast enhancement') +plt.title('Local morphological contrast enhancement') +plt.axis('off') + plt.subplot(2, 2, 2) -plt.imshow(ima[200:350, 350:450], cmap=plt.cm.gray) +plt.imshow(noisy_image[200:350, 350:450], cmap=plt.cm.gray) +plt.axis('off') + plt.subplot(2, 2, 4) plt.imshow(enh[200:350, 350:450], cmap=plt.cm.gray) +plt.axis('off') """ @@ -322,24 +355,30 @@ percentile *p0* and *p1* instead of the local minimum and maximum. """ -from skimage.filter.rank import percentile_morph_contr_enh +from skimage.filter.rank import enhance_contrast_percentile -ima = data.camera() +noisy_image = data.camera() -penh = percentile_morph_contr_enh(ima, disk(5), p0=.1, p1=.9) +penh = enhance_contrast_percentile(noisy_image, 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.imshow(noisy_image, cmap=plt.cm.gray) +plt.title('Original') +plt.axis('off') + plt.subplot(2, 2, 3) plt.imshow(penh, cmap=plt.cm.gray) -plt.xlabel('local percentile morphlogical\n contrast enhancement') +plt.title('Local percentile morphological\n contrast enhancement') +plt.axis('off') + plt.subplot(2, 2, 2) -plt.imshow(ima[200:350, 350:450], cmap=plt.cm.gray) +plt.imshow(noisy_image[200:350, 350:450], cmap=plt.cm.gray) +plt.axis('off') + plt.subplot(2, 2, 4) plt.imshow(penh[200:350, 350:450], cmap=plt.cm.gray) +plt.axis('off') """ @@ -348,18 +387,18 @@ plt.imshow(penh[200:350, 350:450], cmap=plt.cm.gray) 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 Otsu threshold [1]_ method can be applied locally using the local gray- +level 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`. + Local is much slower than global thresholding. A function for global Otsu + thresholding can be found in : `skimage.filter.threshold_otsu`. .. [4] http://en.wikipedia.org/wiki/Otsu's_method @@ -382,27 +421,35 @@ 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.title('Original') plt.colorbar() +plt.axis('off') + plt.subplot(2, 2, 2) plt.imshow(t_loc_otsu, cmap=plt.cm.gray) -plt.xlabel('local Otsu ($radius=%d$)' % radius) +plt.title('Local Otsu ($r=%d$)' % radius) plt.colorbar() +plt.axis('off') + plt.subplot(2, 2, 3) plt.imshow(p8 >= t_loc_otsu, cmap=plt.cm.gray) -plt.xlabel('original>=local Otsu' % t_glob_otsu) +plt.title('Original >= local Otsu' % t_glob_otsu) +plt.axis('off') + plt.subplot(2, 2, 4) plt.imshow(glob_otsu, cmap=plt.cm.gray) -plt.xlabel('global Otsu ($t=%d$)' % t_glob_otsu) +plt.title('Global Otsu ($t=%d$)' % t_glob_otsu) +plt.axis('off') """ .. image:: PLOT2RST.current_figure -The following example shows how local Otsu's threshold handles a global level -shift applied to a synthetic image . +The following example shows how local Otsu thresholding handles a global level +shift applied to a synthetic image. """ @@ -413,13 +460,18 @@ 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.title('Original') +plt.axis('off') + plt.subplot(1, 2, 2) plt.imshow(m >= t, interpolation='nearest') -plt.xlabel('local Otsu ($radius=%d$)' % radius) +plt.title('Local Otsu ($r=%d$)' % radius) +plt.axis('off') """ @@ -428,7 +480,7 @@ plt.xlabel('local Otsu ($radius=%d$)' % radius) Image morphology ================ -Local maximum and local minimum are the base operators for greylevel +Local maximum and local minimum are the base operators for gray-level morphology. .. note:: @@ -436,33 +488,41 @@ morphology. `skimage.dilate` and `skimage.erode` are equivalent filters (see below for comparison). -Here is an example of the classical morphological greylevel filters: opening, +Here is an example of the classical morphological gray-level filters: opening, closing and morphological gradient. """ from skimage.filter.rank import maximum, minimum, gradient -ima = data.camera() +noisy_image = data.camera() -closing = maximum(minimum(ima, disk(5)), disk(5)) -opening = minimum(maximum(ima, disk(5)), disk(5)) -grad = gradient(ima, disk(5)) +closing = maximum(minimum(noisy_image, disk(5)), disk(5)) +opening = minimum(maximum(noisy_image, disk(5)), disk(5)) +grad = gradient(noisy_image, 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.imshow(noisy_image, cmap=plt.cm.gray) +plt.title('Original') +plt.axis('off') + plt.subplot(2, 2, 2) plt.imshow(closing, cmap=plt.cm.gray) -plt.xlabel('greylevel closing') +plt.title('Gray-level closing') +plt.axis('off') + plt.subplot(2, 2, 3) plt.imshow(opening, cmap=plt.cm.gray) -plt.xlabel('greylevel opening') +plt.title('Gray-level opening') +plt.axis('off') + plt.subplot(2, 2, 4) plt.imshow(grad, cmap=plt.cm.gray) -plt.xlabel('morphological gradient') +plt.title('Morphological gradient') +plt.axis('off') """ @@ -471,13 +531,14 @@ plt.xlabel('morphological gradient') Feature extraction =================== -Local histogram can be exploited to compute local entropy, which is related to +Local histograms 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 +filter returns the minimum number of bits needed to encode local gray-level 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. +`skimage.rank.entropy` returns the local entropy on a given structuring +element. The following example shows applies this filter on 8- and 16-bit +images. .. note:: @@ -492,47 +553,36 @@ 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 +image = data.camera() -ent8 = entropy(a8, disk(5)) # pixel value contain 10x the local entropy -ent16 = entropy(a16, disk(5)) # pixel value contain 1000x the local entropy +plt.figure(figsize=(10, 4)) -# 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.subplot(1, 2, 1) +plt.imshow(image, cmap=plt.cm.gray) +plt.title('Image') plt.colorbar() +plt.axis('off') -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.subplot(1, 2, 2) +plt.imshow(entropy(image, disk(5)), cmap=plt.cm.jet) +plt.title('Entropy') plt.colorbar() +plt.axis('off') """ .. 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. +The central part of the `skimage.rank` filters is build on a sliding window +that updates the local gray-level 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. + +In the following we compare the performance of different implementations +available in `skimage`. """ @@ -583,10 +633,10 @@ def ndi_med(image, n): Comparison between -* `rank.maximum` -* `cmorph.dilate` +* `filter.rank.maximum` +* `morphology.dilate` -on increasing structuring element size +on increasing structuring element size: """ @@ -603,18 +653,18 @@ for r in e_range: rec = np.asarray(rec) plt.figure() -plt.title('increasing element size') -plt.ylabel('time (ms)') -plt.xlabel('element radius') +plt.title('Performance with respect to element size') +plt.ylabel('Time (ms)') +plt.title('Element radius') plt.plot(e_range, rec) -plt.legend(['crank.maximum', 'cmorph.dilate']) +plt.legend(['filter.rank.maximum', 'morphology.dilate']) """ -and increasing image size - .. image:: PLOT2RST.current_figure +and increasing image size: + """ r = 9 @@ -623,7 +673,7 @@ 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') + a = (np.random.random((s, s)) * 256).astype(np.uint8) (rc, ms_rc) = cr_max(a, elem) (rcm, ms_rcm) = cm_dil(a, elem) rec.append((ms_rc, ms_rcm)) @@ -631,11 +681,11 @@ for s in s_range: rec = np.asarray(rec) plt.figure() -plt.title('increasing image size') -plt.ylabel('time (ms)') -plt.xlabel('image size') +plt.title('Performance with respect to image size') +plt.ylabel('Time (ms)') +plt.title('Image size') plt.plot(s_range, rec) -plt.legend(['crank.maximum', 'cmorph.dilate']) +plt.legend(['filter.rank.maximum', 'morphology.dilate']) """ @@ -644,11 +694,11 @@ plt.legend(['crank.maximum', 'cmorph.dilate']) Comparison between: -* `rank.median` -* `ctmf.median_filter` -* `ndimage.percentile` +* `filter.rank.median` +* `filter.median_filter` +* `scipy.ndimage.percentile` -on increasing structuring element size +on increasing structuring element size: """ @@ -666,27 +716,29 @@ for r in e_range: rec = np.asarray(rec) plt.figure() -plt.title('increasing element size') +plt.title('Performance with respect to element size') plt.plot(e_range, rec) -plt.legend(['rank.median', 'ctmf.median_filter', 'ndimage.percentile']) -plt.ylabel('time (ms)') -plt.xlabel('element radius') +plt.legend(['filter.rank.median', 'filter.median_filter', + 'scipy.ndimage.percentile']) +plt.ylabel('Time (ms)') +plt.title('Element radius') """ .. image:: PLOT2RST.current_figure -comparison of outcome of the three methods +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') +plt.title('filter.rank.median vs filtermedian_filter vs scipy.ndimage.percentile') +plt.axis('off') """ .. image:: PLOT2RST.current_figure -and increasing image size +and increasing image size: """ @@ -696,7 +748,7 @@ 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') + a = (np.random.random((s, s)) * 256).astype(np.uint8) (rc, ms_rc) = cr_med(a, elem) rctmf, ms_rctmf = ctmf_med(a, r) rndi, ms_ndi = ndi_med(a, r) @@ -705,11 +757,12 @@ for s in s_range: rec = np.asarray(rec) plt.figure() -plt.title('increasing image size') +plt.title('Performance with respect to image size') plt.plot(s_range, rec) -plt.legend(['rank.median', 'ctmf.median_filter', 'ndimage.percentile']) -plt.ylabel('time (ms)') -plt.xlabel('image size') +plt.legend(['filter.rank.median', 'filter.median_filter', + 'scipy.ndimage.percentile']) +plt.ylabel('Time (ms)') +plt.title('Image size') """ .. image:: PLOT2RST.current_figure diff --git a/doc/examples/plot_entropy.py b/doc/examples/plot_entropy.py index db230d8b..f5001e31 100644 --- a/doc/examples/plot_entropy.py +++ b/doc/examples/plot_entropy.py @@ -3,6 +3,9 @@ Entropy ======= +Image entropy is a quantity which is used to describe the amount of information +coded in an image. + """ import numpy as np import matplotlib.pyplot as plt @@ -13,33 +16,18 @@ from skimage.morphology import disk from skimage.util import img_as_ubyte -# defining a 8- and a 16-bit test images -a8 = img_as_ubyte(data.camera()) -a16 = a8.astype(np.uint16) * 4 +image = img_as_ubyte(data.camera()) -ent8 = entropy(a8, disk(5)) # pixel value contain 10x the local entropy -ent16 = entropy(a16, disk(5)) # pixel value contain 1000x the local entropy +fig, (ax0, ax1) = plt.subplots(ncols=2, figsize=(10, 4)) -# display results -plt.figure(figsize=(10, 10)) +img0 = ax0.imshow(image, cmap=plt.cm.gray) +ax0.set_title('Image') +ax0.axis('off') +plt.colorbar(img0, ax=ax0) -plt.subplot(2,2,1) -plt.imshow(a8, cmap=plt.cm.gray) -plt.xlabel('8-bit image') -plt.colorbar() +img1 = ax1.imshow(entropy(image, disk(5)), cmap=plt.cm.jet) +ax1.set_title('Entropy') +ax1.axis('off') +plt.colorbar(img1, ax=ax1) -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_rank_mean.py b/doc/examples/plot_rank_mean.py index e23beb18..6f16c440 100644 --- a/doc/examples/plot_rank_mean.py +++ b/doc/examples/plot_rank_mean.py @@ -6,9 +6,9 @@ Mean filters This example compares the following mean filters of the rank filter package: * **local mean**: all pixels belonging to the structuring element to compute - average gray level + average gray level. * **percentile mean**: only use values between percentiles p0 and p1 - (here 10% and 90%) + (here 10% and 90%). * **bilateral mean**: only use pixels of the structuring element having a gray level situated inside g-s0 and g+s1 (here g-500 and g+500) @@ -23,23 +23,30 @@ import matplotlib.pyplot as plt from skimage import data from skimage.morphology import disk -import skimage.filter.rank as rank +from skimage.filter import rank -a16 = (data.coins()).astype(np.uint16) * 16 + +image = (data.coins()).astype(np.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) +percentile_result = rank.mean_percentile(image, selem=selem, p0=.1, p1=.9) +bilateral_result = rank.mean_bilateral(image, selem=selem, s0=500, s1=500) +normal_result = rank.mean(image, selem=selem) -# display results -fig, axes = plt.subplots(nrows=3, figsize=(15, 10)) + +fig, axes = plt.subplots(nrows=3, figsize=(8, 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') +ax0.imshow(np.hstack((image, percentile_result))) +ax0.set_title('Percentile mean') +ax0.axis('off') + +ax1.imshow(np.hstack((image, bilateral_result))) +ax1.set_title('Bilateral mean') +ax1.axis('off') + +ax2.imshow(np.hstack((image, normal_result))) +ax2.set_title('Local mean') +ax2.axis('off') + plt.show() diff --git a/skimage/filter/rank/README.rst b/skimage/filter/rank/README.rst index cdf8205c..e5c5a9ad 100644 --- a/skimage/filter/rank/README.rst +++ b/skimage/filter/rank/README.rst @@ -23,10 +23,10 @@ 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. +images and 2 to 16-bit for 16-bit images depending on the maximum value of the +image. -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. +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 index deceaade..cfd034f1 100644 --- a/skimage/filter/rank/__init__.py +++ b/skimage/filter/rank/__init__.py @@ -1,36 +1,69 @@ -from ._rank import (autolevel, bottomhat, equalize, gradient, maximum, mean, - meansubtraction, median, minimum, modal, morph_contr_enh, - pop, threshold, tophat, noise_filter, entropy, otsu) -from .percentile_rank import (percentile_autolevel, percentile_gradient, - percentile_mean, percentile_mean_subtraction, - percentile_morph_contr_enh, percentile, - percentile_pop, percentile_threshold) -from .bilateral_rank import bilateral_mean, bilateral_pop +from .generic import (autolevel, bottomhat, equalize, gradient, maximum, mean, + subtract_mean, median, minimum, modal, enhance_contrast, + pop, threshold, tophat, noise_filter, entropy, otsu) +from .percentile import (autolevel_percentile, gradient_percentile, + mean_percentile, subtract_mean_percentile, + enhance_contrast_percentile, percentile, + pop_percentile, threshold_percentile) +from .bilateral import mean_bilateral, pop_bilateral + +from skimage._shared.utils import deprecated + + +percentile_autolevel = deprecated('autolevel_percentile')(autolevel_percentile) + +percentile_gradient = deprecated('gradient_percentile')(gradient_percentile) + +percentile_mean = deprecated('mean_percentile')(mean_percentile) +bilateral_mean = deprecated('mean_bilateral')(mean_bilateral) + +meansubtraction = deprecated('subtract_mean')(subtract_mean) +percentile_mean_subtraction = deprecated('subtract_mean_percentile')\ + (subtract_mean_percentile) + +morph_contr_enh = deprecated('enhance_contrast')(enhance_contrast) +percentile_morph_contr_enh = deprecated('enhance_contrast_percentile')\ + (enhance_contrast_percentile) + +percentile_pop = deprecated('pop_percentile')(pop_percentile) +bilateral_pop = deprecated('pop_bilateral')(pop_bilateral) + +percentile_threshold = deprecated('threshold_percentile')(threshold_percentile) __all__ = ['autolevel', + 'autolevel_percentile', 'bottomhat', 'equalize', 'gradient', + 'gradient_percentile', 'maximum', 'mean', - 'meansubtraction', + 'mean_percentile', + 'mean_bilateral', + 'subtract_mean', + 'subtract_mean_percentile', 'median', 'minimum', 'modal', - 'morph_contr_enh', + 'enhance_contrast', + 'enhance_contrast_percentile', 'pop', + 'pop_percentile', + 'pop_bilateral', 'threshold', + 'threshold_percentile', 'tophat', 'noise_filter', 'entropy', - 'otsu', + 'otsu' + 'percentile', + # Deprecated 'percentile_autolevel', 'percentile_gradient', 'percentile_mean', 'percentile_mean_subtraction', 'percentile_morph_contr_enh', - 'percentile', 'percentile_pop', 'percentile_threshold', 'bilateral_mean', diff --git a/skimage/filter/rank/_core16.pxd b/skimage/filter/rank/_core16.pxd deleted file mode 100644 index 5586aea1..00000000 --- a/skimage/filter/rank/_core16.pxd +++ /dev/null @@ -1,20 +0,0 @@ -cimport numpy as cnp - - -ctypedef cnp.uint16_t dtype_t - - -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(dtype_t kernel(Py_ssize_t *, float, dtype_t, - Py_ssize_t, Py_ssize_t, Py_ssize_t, float, - float, Py_ssize_t, Py_ssize_t), - cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[cnp.uint8_t, ndim=2] selem, - cnp.ndarray[cnp.uint8_t, ndim=2] mask, - cnp.ndarray[dtype_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 deleted file mode 100644 index 0c7a7a82..00000000 --- a/skimage/filter/rank/_core16.pyx +++ /dev/null @@ -1,255 +0,0 @@ -#cython: cdivision=True -#cython: boundscheck=False -#cython: nonecheck=False -#cython: wraparound=False - -import numpy as np - -cimport numpy as cnp -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, - dtype_t value): - histo[value] += 1 - pop[0] += 1 - - -cdef inline void histogram_decrement(Py_ssize_t * histo, float * pop, - dtype_t value): - histo[value] -= 1 - pop[0] -= 1 - - -cdef void _core16(dtype_t kernel(Py_ssize_t *, float, dtype_t, - Py_ssize_t, Py_ssize_t, Py_ssize_t, float, - float, Py_ssize_t, Py_ssize_t), - cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[cnp.uint8_t, ndim=2] selem, - cnp.ndarray[cnp.uint8_t, ndim=2] mask, - cnp.ndarray[dtype_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 dtype_t * out_data = out.data - cdef dtype_t * image_data = image.data - cdef cnp.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) < 0 - - t = np.hstack((np.zeros((selem.shape[0], 1)), selem)) - t_w = np.diff(t, axis=1) > 0 - - t = np.vstack((selem, np.zeros((1, selem.shape[1])))) - t_s = np.diff(t, axis=0) < 0 - - t = np.vstack((np.zeros((1, selem.shape[1])), selem)) - t_n = np.diff(t, axis=0) > 0 - - 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 deleted file mode 100644 index d3b6d8c2..00000000 --- a/skimage/filter/rank/_core8.pxd +++ /dev/null @@ -1,25 +0,0 @@ -cimport numpy as cnp - - -ctypedef cnp.uint8_t dtype_t - - -cdef dtype_t uint8_max(dtype_t a, dtype_t b) -cdef dtype_t uint8_min(dtype_t a, dtype_t b) - - -cdef dtype_t is_in_mask(Py_ssize_t rows, Py_ssize_t cols, - Py_ssize_t r, Py_ssize_t c, - dtype_t * mask) - - -# 8-bit core kernel receives extra information about data inferior and superior -# percentiles -cdef void _core8(dtype_t kernel(Py_ssize_t *, float, dtype_t, float, - float, Py_ssize_t, Py_ssize_t), - cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[dtype_t, ndim=2] selem, - cnp.ndarray[dtype_t, ndim=2] mask, - cnp.ndarray[dtype_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/_crank16.pyx b/skimage/filter/rank/_crank16.pyx deleted file mode 100644 index e704afe0..00000000 --- a/skimage/filter/rank/_crank16.pyx +++ /dev/null @@ -1,422 +0,0 @@ -#cython: cdivision=True -#cython: boundscheck=False -#cython: nonecheck=False -#cython: wraparound=False - -cimport numpy as cnp -from libc.math cimport log -from skimage.filter.rank._core16 cimport _core16 - - -# ----------------------------------------------------------------- -# kernels uint16 take extra parameter for defining the bitdepth -# ----------------------------------------------------------------- - - -ctypedef cnp.uint16_t dtype_t - - -cdef inline dtype_t kernel_autolevel(Py_ssize_t * histo, float pop, - dtype_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 (1. * (maxbin - 1) * (g - imin) / delta) - else: - return (imax - imin) - - -cdef inline dtype_t kernel_bottomhat(Py_ssize_t * histo, float pop, - dtype_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 (g - i) - else: - return (0) - -cdef inline dtype_t kernel_equalize(Py_ssize_t * histo, float pop, - dtype_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 (((maxbin - 1) * sum) / pop) - else: - return (0) - - -cdef inline dtype_t kernel_gradient(Py_ssize_t * histo, float pop, - dtype_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 (imax - imin) - else: - return (0) - - -cdef inline dtype_t kernel_maximum(Py_ssize_t * histo, float pop, - dtype_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 (i) - - return (0) - - -cdef inline dtype_t kernel_mean(Py_ssize_t * histo, float pop, - dtype_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 (mean / pop) - else: - return (0) - - -cdef inline dtype_t kernel_meansubtraction(Py_ssize_t * histo, - float pop, - dtype_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 ((g - mean / pop) / 2. + (midbin - 1)) - else: - return (0) - - -cdef inline dtype_t kernel_median(Py_ssize_t * histo, float pop, - dtype_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 (i) - else: - return (0) - - -cdef inline dtype_t kernel_minimum(Py_ssize_t * histo, float pop, - dtype_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 (i) - else: - return (0) - - -cdef inline dtype_t kernel_modal(Py_ssize_t * histo, float pop, - dtype_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 (imax) - else: - return (0) - - -cdef inline dtype_t kernel_morph_contr_enh(Py_ssize_t * histo, - float pop, - dtype_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 (imax) - else: - return (imin) - else: - return (0) - - -cdef inline dtype_t kernel_pop(Py_ssize_t * histo, float pop, - dtype_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 (pop) - - -cdef inline dtype_t kernel_threshold(Py_ssize_t * histo, float pop, - dtype_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 (g > (mean / pop)) - else: - return (0) - - -cdef inline dtype_t kernel_tophat(Py_ssize_t * histo, float pop, - dtype_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 (i - g) - else: - return (0) - -cdef inline dtype_t kernel_entropy(Py_ssize_t * histo, float pop, - dtype_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 * log(p) / 0.6931471805599453 - - return e * 1000 - else: - return (0) - -# ----------------------------------------------------------------- -# python wrappers -# ----------------------------------------------------------------- - - -def autolevel(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[cnp.uint8_t, ndim=2] selem, - cnp.ndarray[cnp.uint8_t, ndim=2] mask=None, - cnp.ndarray[dtype_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(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[cnp.uint8_t, ndim=2] selem, - cnp.ndarray[cnp.uint8_t, ndim=2] mask=None, - cnp.ndarray[dtype_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(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[cnp.uint8_t, ndim=2] selem, - cnp.ndarray[cnp.uint8_t, ndim=2] mask=None, - cnp.ndarray[dtype_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(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[cnp.uint8_t, ndim=2] selem, - cnp.ndarray[cnp.uint8_t, ndim=2] mask=None, - cnp.ndarray[dtype_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(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[cnp.uint8_t, ndim=2] selem, - cnp.ndarray[cnp.uint8_t, ndim=2] mask=None, - cnp.ndarray[dtype_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(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[cnp.uint8_t, ndim=2] selem, - cnp.ndarray[cnp.uint8_t, ndim=2] mask=None, - cnp.ndarray[dtype_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 meansubtraction(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[cnp.uint8_t, ndim=2] selem, - cnp.ndarray[cnp.uint8_t, ndim=2] mask=None, - cnp.ndarray[dtype_t, ndim=2] out=None, - char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8): - _core16(kernel_meansubtraction, image, selem, mask, out, shift_x, shift_y, - bitdepth, 0, 0, 0, 0) - - -def median(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[cnp.uint8_t, ndim=2] selem, - cnp.ndarray[cnp.uint8_t, ndim=2] mask=None, - cnp.ndarray[dtype_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(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[cnp.uint8_t, ndim=2] selem, - cnp.ndarray[cnp.uint8_t, ndim=2] mask=None, - cnp.ndarray[dtype_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(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[cnp.uint8_t, ndim=2] selem, - cnp.ndarray[cnp.uint8_t, ndim=2] mask=None, - cnp.ndarray[dtype_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(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[cnp.uint8_t, ndim=2] selem, - cnp.ndarray[cnp.uint8_t, ndim=2] mask=None, - cnp.ndarray[dtype_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(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[cnp.uint8_t, ndim=2] selem, - cnp.ndarray[cnp.uint8_t, ndim=2] mask=None, - cnp.ndarray[dtype_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(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[cnp.uint8_t, ndim=2] selem, - cnp.ndarray[cnp.uint8_t, ndim=2] mask=None, - cnp.ndarray[dtype_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(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[cnp.uint8_t, ndim=2] selem, - cnp.ndarray[cnp.uint8_t, ndim=2] mask=None, - cnp.ndarray[dtype_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(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[cnp.uint8_t, ndim=2] selem, - cnp.ndarray[cnp.uint8_t, ndim=2] mask=None, - cnp.ndarray[dtype_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 deleted file mode 100644 index e431e42b..00000000 --- a/skimage/filter/rank/_crank16_bilateral.pyx +++ /dev/null @@ -1,82 +0,0 @@ -#cython: cdivision=True -#cython: boundscheck=False -#cython: nonecheck=False -#cython: wraparound=False - -cimport numpy as cnp -from skimage.filter.rank._core16 cimport _core16 - - -# ----------------------------------------------------------------- -# kernels uint16 take extra parameter for defining the bitdepth -# ----------------------------------------------------------------- - - -ctypedef cnp.uint16_t dtype_t - - -cdef inline dtype_t kernel_mean(Py_ssize_t * histo, float pop, - dtype_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 (mean / bilat_pop) - else: - return (0) - else: - return (0) - - -cdef inline dtype_t kernel_pop(Py_ssize_t * histo, float pop, - dtype_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 (bilat_pop) - else: - return (0) - - -# ----------------------------------------------------------------- -# python wrappers -# ----------------------------------------------------------------- - - -def mean(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[cnp.uint8_t, ndim=2] selem, - cnp.ndarray[cnp.uint8_t, ndim=2] mask=None, - cnp.ndarray[dtype_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(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[cnp.uint8_t, ndim=2] selem, - cnp.ndarray[cnp.uint8_t, ndim=2] mask=None, - cnp.ndarray[dtype_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 deleted file mode 100644 index f4a4c9b2..00000000 --- a/skimage/filter/rank/_crank16_percentiles.pyx +++ /dev/null @@ -1,330 +0,0 @@ -#cython: cdivision=True -#cython: boundscheck=False -#cython: nonecheck=False -#cython: wraparound=False - -cimport numpy as cnp -from skimage.filter.rank._core16 cimport _core16, int_min, int_max - - -# ----------------------------------------------------------------- -# kernels uint16 (SOFT version using percentiles) -# ----------------------------------------------------------------- - - -ctypedef cnp.uint16_t dtype_t - - -cdef inline dtype_t kernel_autolevel(Py_ssize_t * histo, float pop, - dtype_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 (1.0 * (maxbin - 1) - * (int_min(int_max(imin, g), imax) - - imin) / delta) - else: - return (imax - imin) - else: - return (0) - - -cdef inline dtype_t kernel_gradient(Py_ssize_t * histo, float pop, - dtype_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 (imax - imin) - else: - return (0) - - -cdef inline dtype_t kernel_mean(Py_ssize_t * histo, float pop, - dtype_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 (1.0 * mean / n) - else: - return (0) - else: - return (0) - - -cdef inline dtype_t kernel_mean_subtraction(Py_ssize_t * histo, - float pop, - dtype_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 ((g - (mean / n)) * .5 + midbin) - else: - return (0) - else: - return (0) - - -cdef inline dtype_t kernel_morph_contr_enh(Py_ssize_t * histo, - float pop, - dtype_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 imax - if g < imin: - return imin - if imax - g < g - imin: - return imax - else: - return imin - else: - return (0) - - -cdef inline dtype_t kernel_percentile(Py_ssize_t * histo, float pop, - dtype_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 (i) - else: - return (0) - - -cdef inline dtype_t kernel_pop(Py_ssize_t * histo, float pop, - dtype_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 (n) - else: - return (0) - - -cdef inline dtype_t kernel_threshold(Py_ssize_t * histo, float pop, - dtype_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 ((maxbin - 1) * (g >= i)) - else: - return (0) - - -# ----------------------------------------------------------------- -# python wrappers -# ----------------------------------------------------------------- - - -def autolevel(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[cnp.uint8_t, ndim=2] selem, - cnp.ndarray[cnp.uint8_t, ndim=2] mask=None, - cnp.ndarray[dtype_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(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[cnp.uint8_t, ndim=2] selem, - cnp.ndarray[cnp.uint8_t, ndim=2] mask=None, - cnp.ndarray[dtype_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(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[cnp.uint8_t, ndim=2] selem, - cnp.ndarray[cnp.uint8_t, ndim=2] mask=None, - cnp.ndarray[dtype_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_subtraction(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[cnp.uint8_t, ndim=2] selem, - cnp.ndarray[cnp.uint8_t, ndim=2] mask=None, - cnp.ndarray[dtype_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_subtraction, image, selem, mask, out, shift_x, shift_y, - bitdepth, p0, p1, 0, 0) - - -def morph_contr_enh(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[cnp.uint8_t, ndim=2] selem, - cnp.ndarray[cnp.uint8_t, ndim=2] mask=None, - cnp.ndarray[dtype_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(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[cnp.uint8_t, ndim=2] selem, - cnp.ndarray[cnp.uint8_t, ndim=2] mask=None, - cnp.ndarray[dtype_t, ndim=2] out=None, - char shift_x=0, char shift_y=0, int bitdepth=8, - float p0=0.): - """return p0 percentile - """ - _core16(kernel_percentile, image, selem, mask, out, shift_x, shift_y, - bitdepth, p0, .0, 0, 0) - - -def pop(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[cnp.uint8_t, ndim=2] selem, - cnp.ndarray[cnp.uint8_t, ndim=2] mask=None, - cnp.ndarray[dtype_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(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[cnp.uint8_t, ndim=2] selem, - cnp.ndarray[cnp.uint8_t, ndim=2] mask=None, - cnp.ndarray[dtype_t, ndim=2] out=None, - char shift_x=0, char shift_y=0, int bitdepth=8, - float p0=0.): - """return (maxbin-1) if g > percentile p0 - """ - _core16(kernel_threshold, image, selem, mask, out, shift_x, shift_y, - bitdepth, p0, 0., 0, 0) diff --git a/skimage/filter/rank/_crank8.pyx b/skimage/filter/rank/_crank8.pyx deleted file mode 100644 index da511790..00000000 --- a/skimage/filter/rank/_crank8.pyx +++ /dev/null @@ -1,483 +0,0 @@ -#cython: cdivision=True -#cython: boundscheck=False -#cython: nonecheck=False -#cython: wraparound=False - -cimport numpy as cnp -from libc.math cimport log -from skimage.filter.rank._core8 cimport _core8 - - -# ----------------------------------------------------------------- -# kernels uint8 -# ----------------------------------------------------------------- - - -ctypedef cnp.uint8_t dtype_t - - -cdef inline dtype_t kernel_autolevel(Py_ssize_t * histo, float pop, - dtype_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 (255. * (g - imin) / delta) - else: - return (imax - imin) - else: - return (0) - - -cdef inline dtype_t kernel_bottomhat(Py_ssize_t * histo, float pop, - dtype_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 (g - i) - else: - return (0) - - -cdef inline dtype_t kernel_equalize(Py_ssize_t * histo, float pop, - dtype_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 ((255 * sum) / pop) - else: - return (0) - - -cdef inline dtype_t kernel_gradient(Py_ssize_t * histo, float pop, - dtype_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 (imax - imin) - else: - return (0) - - -cdef inline dtype_t kernel_maximum(Py_ssize_t * histo, float pop, - dtype_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 (i) - else: - return (0) - - -cdef inline dtype_t kernel_mean(Py_ssize_t * histo, float pop, - dtype_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 (mean / pop) - else: - return (0) - - -cdef inline dtype_t kernel_meansubtraction(Py_ssize_t * histo, float pop, - dtype_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 ((g - mean / pop) / 2. + 127) - else: - return (0) - - -cdef inline dtype_t kernel_median(Py_ssize_t * histo, float pop, - dtype_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 (i) - else: - return (0) - - -cdef inline dtype_t kernel_minimum(Py_ssize_t * histo, float pop, - dtype_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 (i) - else: - return (0) - - -cdef inline dtype_t kernel_modal(Py_ssize_t * histo, float pop, - dtype_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 (imax) - else: - return (0) - - -cdef inline dtype_t kernel_morph_contr_enh(Py_ssize_t * histo, float pop, - dtype_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 (imax) - else: - return (imin) - else: - return (0) - - -cdef inline dtype_t kernel_pop(Py_ssize_t * histo, float pop, - dtype_t g, float p0, float p1, - Py_ssize_t s0, Py_ssize_t s1): - - return (pop) - - -cdef inline dtype_t kernel_threshold(Py_ssize_t * histo, float pop, - dtype_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 (g > (mean / pop)) - else: - return (0) - - -cdef inline dtype_t kernel_tophat(Py_ssize_t * histo, float pop, - dtype_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 (i - g) - else: - return (0) - -cdef inline dtype_t kernel_noise_filter(Py_ssize_t * histo, float pop, - dtype_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 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 (i - g) - else: - return min_i - - -cdef inline dtype_t kernel_entropy(Py_ssize_t * histo, float pop, - dtype_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 * log(p) / 0.6931471805599453 - - return e * 10 - else: - return (0) - -cdef inline dtype_t kernel_otsu(Py_ssize_t * histo, float pop, dtype_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 (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 max_i - - -# ----------------------------------------------------------------- -# python wrappers -# used only internally -# ----------------------------------------------------------------- - - -def autolevel(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[dtype_t, ndim=2] selem, - cnp.ndarray[dtype_t, ndim=2] mask=None, - cnp.ndarray[dtype_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(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[dtype_t, ndim=2] selem, - cnp.ndarray[dtype_t, ndim=2] mask=None, - cnp.ndarray[dtype_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(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[dtype_t, ndim=2] selem, - cnp.ndarray[dtype_t, ndim=2] mask=None, - cnp.ndarray[dtype_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(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[dtype_t, ndim=2] selem, - cnp.ndarray[dtype_t, ndim=2] mask=None, - cnp.ndarray[dtype_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(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[dtype_t, ndim=2] selem, - cnp.ndarray[dtype_t, ndim=2] mask=None, - cnp.ndarray[dtype_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(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[dtype_t, ndim=2] selem, - cnp.ndarray[dtype_t, ndim=2] mask=None, - cnp.ndarray[dtype_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 meansubtraction(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[dtype_t, ndim=2] selem, - cnp.ndarray[dtype_t, ndim=2] mask=None, - cnp.ndarray[dtype_t, ndim=2] out=None, - char shift_x=0, char shift_y=0): - _core8(kernel_meansubtraction, image, selem, mask, out, shift_x, shift_y, - 0, 0, 0, 0) - - -def median(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[dtype_t, ndim=2] selem, - cnp.ndarray[dtype_t, ndim=2] mask=None, - cnp.ndarray[dtype_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(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[dtype_t, ndim=2] selem, - cnp.ndarray[dtype_t, ndim=2] mask=None, - cnp.ndarray[dtype_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(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[dtype_t, ndim=2] selem, - cnp.ndarray[dtype_t, ndim=2] mask=None, - cnp.ndarray[dtype_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(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[dtype_t, ndim=2] selem, - cnp.ndarray[dtype_t, ndim=2] mask=None, - cnp.ndarray[dtype_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(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[dtype_t, ndim=2] selem, - cnp.ndarray[dtype_t, ndim=2] mask=None, - cnp.ndarray[dtype_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(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[dtype_t, ndim=2] selem, - cnp.ndarray[dtype_t, ndim=2] mask=None, - cnp.ndarray[dtype_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(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[dtype_t, ndim=2] selem, - cnp.ndarray[dtype_t, ndim=2] mask=None, - cnp.ndarray[dtype_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(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[dtype_t, ndim=2] selem, - cnp.ndarray[dtype_t, ndim=2] mask=None, - cnp.ndarray[dtype_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(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[dtype_t, ndim=2] selem, - cnp.ndarray[dtype_t, ndim=2] mask=None, - cnp.ndarray[dtype_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(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[dtype_t, ndim=2] selem, - cnp.ndarray[dtype_t, ndim=2] mask=None, - cnp.ndarray[dtype_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 deleted file mode 100644 index 8e5cee9c..00000000 --- a/skimage/filter/rank/_crank8_percentiles.pyx +++ /dev/null @@ -1,294 +0,0 @@ -#cython: cdivision=True -#cython: boundscheck=False -#cython: nonecheck=False -#cython: wraparound=False - -cimport numpy as cnp -from skimage.filter.rank._core8 cimport _core8, uint8_max, uint8_min - - -# ----------------------------------------------------------------- -# kernels uint8 (SOFT version using percentiles) -# ----------------------------------------------------------------- - - -ctypedef cnp.uint8_t dtype_t - - -cdef inline dtype_t kernel_autolevel(Py_ssize_t * histo, float pop, - dtype_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 (255 * (uint8_min(uint8_max(imin, g), imax) - - imin) / delta) - else: - return (imax - imin) - else: - return (128) - - -cdef inline dtype_t kernel_gradient(Py_ssize_t * histo, float pop, - dtype_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 (imax - imin) - else: - return (0) - - -cdef inline dtype_t kernel_mean(Py_ssize_t * histo, float pop, - dtype_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 (1.0 * mean / n) - else: - return (0) - else: - return (0) - - -cdef inline dtype_t kernel_mean_subtraction(Py_ssize_t * histo, - float pop, - dtype_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 ((g - (mean / n)) * .5 + 127) - else: - return (0) - else: - return (0) - - -cdef inline dtype_t kernel_morph_contr_enh(Py_ssize_t * histo, - float pop, - dtype_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 imax - if g < imin: - return imin - if imax - g < g - imin: - return imax - else: - return imin - else: - return (0) - - -cdef inline dtype_t kernel_percentile(Py_ssize_t * histo, float pop, - dtype_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 (i) - else: - return (0) - - -cdef inline dtype_t kernel_pop(Py_ssize_t * histo, float pop, - dtype_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 (n) - else: - return (0) - - -cdef inline dtype_t kernel_threshold(Py_ssize_t * histo, float pop, - dtype_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 (255 * (g >= i)) - else: - return (0) - - -# ----------------------------------------------------------------- -# python wrappers -# ----------------------------------------------------------------- - - -def autolevel(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[dtype_t, ndim=2] selem, - cnp.ndarray[dtype_t, ndim=2] mask=None, - cnp.ndarray[dtype_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(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[dtype_t, ndim=2] selem, - cnp.ndarray[dtype_t, ndim=2] mask=None, - cnp.ndarray[dtype_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(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[dtype_t, ndim=2] selem, - cnp.ndarray[dtype_t, ndim=2] mask=None, - cnp.ndarray[dtype_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_subtraction(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[dtype_t, ndim=2] selem, - cnp.ndarray[dtype_t, ndim=2] mask=None, - cnp.ndarray[dtype_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_subtraction, image, selem, mask, out, shift_x, shift_y, - p0, p1, 0, 0) - - -def morph_contr_enh(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[dtype_t, ndim=2] selem, - cnp.ndarray[dtype_t, ndim=2] mask=None, - cnp.ndarray[dtype_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(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[dtype_t, ndim=2] selem, - cnp.ndarray[dtype_t, ndim=2] mask=None, - cnp.ndarray[dtype_t, ndim=2] out=None, - char shift_x=0, char shift_y=0, float p0=0.): - """return p0 percentile - """ - _core8(kernel_percentile, image, selem, mask, out, shift_x, shift_y, - p0, 0., 0, 0) - - -def pop(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[dtype_t, ndim=2] selem, - cnp.ndarray[dtype_t, ndim=2] mask=None, - cnp.ndarray[dtype_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(cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[dtype_t, ndim=2] selem, - cnp.ndarray[dtype_t, ndim=2] mask=None, - cnp.ndarray[dtype_t, ndim=2] out=None, - char shift_x=0, char shift_y=0, float p0=0.): - """return 255 if g > percentile p0 - """ - _core8(kernel_threshold, image, selem, mask, out, shift_x, shift_y, p0, 0., - 0, 0) diff --git a/skimage/filter/rank/_rank.py b/skimage/filter/rank/_rank.py deleted file mode 100644 index ea58ad86..00000000 --- a/skimage/filter/rank/_rank.py +++ /dev/null @@ -1,773 +0,0 @@ -"""The local histogram is computed using a sliding window similar to the method -described in [1]_. - -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. - -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. - -""" - -import numpy as np -from skimage import img_as_ubyte, img_as_uint -from skimage.filter.rank import _crank8, _crank16 -from skimage.filter.rank.generic import find_bitdepth - - -__all__ = ['autolevel', 'bottomhat', 'equalize', 'gradient', 'maximum', 'mean', - 'meansubtraction', '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 > 0) - 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.") - - is_8bit = image.dtype in (np.uint8, np.int8) - - if func8 is not None and (is_8bit or func16 is None): - out = _apply8(func8, image, selem, out, mask, shift_x, shift_y) - else: - image = img_as_uint(image) - if out is None: - out = np.zeros(image.shape, dtype=np.uint16) - bitdepth = find_bitdepth(image) - if bitdepth > 11: - image = image >> 4 - bitdepth = find_bitdepth(image) - func16(image, selem, shift_x=shift_x, shift_y=shift_y, mask=mask, - bitdepth=bitdepth + 1, out=out) - - return out - - -def _apply8(func8, image, selem, out, mask, shift_x, shift_y): - if out is None: - out = np.zeros(image.shape, dtype=np.uint8) - image = img_as_ubyte(image) - func8(image, selem, shift_x=shift_x, shift_y=shift_y, - mask=mask, out=out) - 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 meansubtraction(image, selem, out=None, mask=None, shift_x=False, - shift_y=False): - """Return image subtracted 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 meansubtraction. - - """ - - return _apply(_crank8.meansubtraction, _crank16.meansubtraction, 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 [1]_ 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 - ---------- - .. [1] 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). - 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 - Otsu's threshold values - - References - ---------- - .. [otsu] http://en.wikipedia.org/wiki/Otsu's_method - - Notes - ----- - * input image are 8-bit only - - Examples - -------- - >>> # Local entropy - >>> from skimage import data - >>> from skimage.filter.rank import otsu - >>> from skimage.morphology import disk - >>> # defining a 8-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/bilateral_rank.pyx b/skimage/filter/rank/bilateral.py similarity index 53% rename from skimage/filter/rank/bilateral_rank.pyx rename to skimage/filter/rank/bilateral.py index e2a1fcf3..f1b10fec 100644 --- a/skimage/filter/rank/bilateral_rank.pyx +++ b/skimage/filter/rank/bilateral.py @@ -3,19 +3,16 @@ The local histogram is computed using a sliding window similar to the method described in [1]_. -Input image must be 16-bit with a value < 4096 (i.e. 12 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 +* 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. +Result image is 8-/16-bit or double with respect to the input image and the +rank filter operation. References ---------- @@ -28,50 +25,27 @@ References 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 + +from . import bilateral_cy +from .generic import _handle_input -__all__ = ['bilateral_mean', 'bilateral_pop'] +__all__ = ['mean_bilateral', 'pop_bilateral'] -def _apply(func8, func16, image, selem, out, mask, shift_x, shift_y, s0, s1): - selem = img_as_ubyte(selem > 0) - image = np.ascontiguousarray(image) +def _apply(func, image, selem, out, mask, shift_x, shift_y, s0, s1, + out_dtype=None): - if mask is None: - mask = np.ones(image.shape, dtype=np.uint8) - else: - mask = np.ascontiguousarray(mask) - mask = img_as_ubyte(mask) + image, selem, out, mask, max_bin = _handle_input(image, selem, out, mask, + out_dtype) - 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.") + func(image, selem, shift_x=shift_x, shift_y=shift_y, mask=mask, + out=out, max_bin=max_bin, s0=s0, s1=s1) return out -def bilateral_mean(image, selem, out=None, mask=None, shift_x=False, +def mean_bilateral(image, selem, out=None, mask=None, shift_x=False, shift_y=False, s0=10, s1=10): """Apply a flat kernel bilateral filter. @@ -81,43 +55,38 @@ def bilateral_mean(image, selem, out=None, mask=None, shift_x=False, 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] + 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 (uint16). As the algorithm uses max. 12bit histogram, - an exception will be raised if image has a value > 4095 + image : ndarray (uint8, uint16) + Image array. selem : ndarray The neighborhood expressed as a 2-D array of 1's and 0's. - out : ndarray + out : ndarray (same dtype as input) If None, a new array will be allocated. - mask : ndarray (uint8) + mask : ndarray 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) + 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. + Define the [s0, s1] interval around the greyvalue of the center pixel + to be considered for computing the value. Returns ------- - out : uint16 array - The result of the local bilateral mean. + out : ndarray (same dtype as input image) + Output image. See also -------- - skimage.filter.denoise_bilateral() for a gaussian bilateral filter. - - Notes - ----- - - * input image are 16-bit only + skimage.filter.denoise_bilateral for a gaussian bilateral filter. Examples -------- @@ -128,13 +97,14 @@ def bilateral_mean(image, selem, out=None, mask=None, shift_x=False, >>> ima = data.camera().astype(np.uint16) >>> # 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, + return _apply(bilateral_cy._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, +def pop_bilateral(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 @@ -142,32 +112,27 @@ def bilateral_pop(image, selem, out=None, mask=None, shift_x=False, Parameters ---------- - image : ndarray - Image array (uint16). As the algorithm uses max. 12bit histogram, - an exception will be raised if image has a value > 4095 + image : ndarray (uint8, uint16) + Image array. selem : ndarray The neighborhood expressed as a 2-D array of 1's and 0's. - out : ndarray + out : ndarray (same dtype as input) If None, a new array will be allocated. - mask : ndarray (uint8) + mask : ndarray 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) + 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. + Define the [s0, s1] interval around the greyvalue of the center pixel + to be considered for computing the value. Returns ------- - out : uint16 array - the local number of pixels inside the bilateral neighborhood - - Notes - ----- - - * input image are 16-bit only + out : ndarray (same dtype as input image) + Output image. Examples -------- @@ -175,10 +140,10 @@ def bilateral_pop(image, selem, out=None, mask=None, shift_x=False, >>> from skimage.morphology import square >>> import skimage.filter.rank as rank >>> ima16 = 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.uint16) + ... [0, 1, 1, 1, 0], + ... [0, 1, 1, 1, 0], + ... [0, 1, 1, 1, 0], + ... [0, 0, 0, 0, 0]], dtype=np.uint16) >>> rank.bilateral_pop(ima16, square(3), s0=10,s1=10) array([[3, 4, 3, 4, 3], [4, 4, 6, 4, 4], @@ -188,5 +153,5 @@ def bilateral_pop(image, selem, out=None, mask=None, shift_x=False, """ - return _apply(None, _crank16_bilateral.pop, image, selem, out=out, + return _apply(bilateral_cy._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/bilateral_cy.pyx b/skimage/filter/rank/bilateral_cy.pyx new file mode 100644 index 00000000..de2b3e53 --- /dev/null +++ b/skimage/filter/rank/bilateral_cy.pyx @@ -0,0 +1,70 @@ +#cython: cdivision=True +#cython: boundscheck=False +#cython: nonecheck=False +#cython: wraparound=False + +cimport numpy as cnp +from libc.math cimport log + +from .core_cy cimport dtype_t, dtype_t_out, _core + + +cdef inline double _kernel_mean(Py_ssize_t* histo, double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): + + cdef Py_ssize_t i + cdef Py_ssize_t bilat_pop = 0 + cdef Py_ssize_t mean = 0 + + if pop: + for i in range(max_bin): + if (g > (i - s0)) and (g < (i + s1)): + bilat_pop += histo[i] + mean += histo[i] * i + if bilat_pop: + return mean / bilat_pop + else: + return 0 + else: + return 0 + + +cdef inline double _kernel_pop(Py_ssize_t* histo, double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): + + cdef Py_ssize_t i + cdef Py_ssize_t bilat_pop = 0 + + if pop: + for i in range(max_bin): + if (g > (i - s0)) and (g < (i + s1)): + bilat_pop += histo[i] + return bilat_pop + else: + return 0 + + +def _mean(dtype_t[:, ::1] image, + char[:, ::1] selem, + char[:, ::1] mask, + dtype_t_out[:, ::1] out, + char shift_x, char shift_y, Py_ssize_t s0, Py_ssize_t s1, + Py_ssize_t max_bin): + + _core(_kernel_mean[dtype_t], image, selem, mask, out, + shift_x, shift_y, 0, 0, s0, s1, max_bin) + + +def _pop(dtype_t[:, ::1] image, + char[:, ::1] selem, + char[:, ::1] mask, + dtype_t_out[:, ::1] out, + char shift_x, char shift_y, Py_ssize_t s0, Py_ssize_t s1, + Py_ssize_t max_bin): + + _core(_kernel_pop[dtype_t], image, selem, mask, out, + shift_x, shift_y, 0, 0, s0, s1, max_bin) diff --git a/skimage/filter/rank/core_cy.pxd b/skimage/filter/rank/core_cy.pxd new file mode 100644 index 00000000..2e97e50a --- /dev/null +++ b/skimage/filter/rank/core_cy.pxd @@ -0,0 +1,28 @@ +from numpy cimport uint8_t, uint16_t, double_t + + +ctypedef fused dtype_t: + uint8_t + uint16_t + +ctypedef fused dtype_t_out: + uint8_t + uint16_t + double_t + + +cdef dtype_t _max(dtype_t a, dtype_t b) +cdef dtype_t _min(dtype_t a, dtype_t b) + + +cdef void _core(double kernel(Py_ssize_t*, double, dtype_t, + Py_ssize_t, Py_ssize_t, double, + double, Py_ssize_t, Py_ssize_t), + dtype_t[:, ::1] image, + char[:, ::1] selem, + char[:, ::1] mask, + dtype_t_out[:, ::1] out, + char shift_x, char shift_y, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1, + Py_ssize_t max_bin) except * diff --git a/skimage/filter/rank/_core8.pyx b/skimage/filter/rank/core_cy.pyx similarity index 51% rename from skimage/filter/rank/_core8.pyx rename to skimage/filter/rank/core_cy.pyx index 79ae9bbf..02c2c8d0 100644 --- a/skimage/filter/rank/_core8.pyx +++ b/skimage/filter/rank/core_cy.pyx @@ -9,29 +9,29 @@ cimport numpy as cnp from libc.stdlib cimport malloc, free -cdef inline dtype_t uint8_max(dtype_t a, dtype_t b): +cdef inline dtype_t _max(dtype_t a, dtype_t b): return a if a >= b else b -cdef inline dtype_t uint8_min(dtype_t a, dtype_t b): +cdef inline dtype_t _min(dtype_t a, dtype_t b): return a if a <= b else b -cdef inline void histogram_increment(Py_ssize_t * histo, float * pop, +cdef inline void histogram_increment(Py_ssize_t* histo, double* pop, dtype_t value): histo[value] += 1 pop[0] += 1 -cdef inline void histogram_decrement(Py_ssize_t * histo, float * pop, +cdef inline void histogram_decrement(Py_ssize_t* histo, double* pop, dtype_t value): histo[value] -= 1 pop[0] -= 1 -cdef inline dtype_t is_in_mask(Py_ssize_t rows, Py_ssize_t cols, - Py_ssize_t r, Py_ssize_t c, - dtype_t * mask): +cdef inline char is_in_mask(Py_ssize_t rows, Py_ssize_t cols, + Py_ssize_t r, Py_ssize_t c, + char* 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 @@ -42,14 +42,17 @@ cdef inline dtype_t is_in_mask(Py_ssize_t rows, Py_ssize_t cols, return 0 -cdef void _core8(dtype_t kernel(Py_ssize_t *, float, dtype_t, float, - float, Py_ssize_t, Py_ssize_t), - cnp.ndarray[dtype_t, ndim=2] image, - cnp.ndarray[dtype_t, ndim=2] selem, - cnp.ndarray[dtype_t, ndim=2] mask, - cnp.ndarray[dtype_t, ndim=2] out, - char shift_x, char shift_y, float p0, float p1, - Py_ssize_t s0, Py_ssize_t s1) except *: +cdef void _core(double kernel(Py_ssize_t*, double, dtype_t, + Py_ssize_t, Py_ssize_t, double, + double, Py_ssize_t, Py_ssize_t), + dtype_t[:, ::1] image, + char[:, ::1] selem, + char[:, ::1] mask, + dtype_t_out[:, ::1] out, + char shift_x, char shift_y, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1, + Py_ssize_t max_bin) except *: """Compute histogram for each pixel neighborhood, apply kernel function and use kernel function return value for output image. """ @@ -59,8 +62,8 @@ cdef void _core8(dtype_t kernel(Py_ssize_t *, float, dtype_t, float, 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 + cdef Py_ssize_t centre_r = (selem.shape[0] / 2) + shift_y + cdef Py_ssize_t centre_c = (selem.shape[1] / 2) + shift_x # check that structuring element center is inside the element bounding box assert centre_r >= 0 @@ -68,54 +71,56 @@ cdef void _core8(dtype_t kernel(Py_ssize_t *, float, dtype_t, float, assert centre_r < srows assert centre_c < scols - # define pointers to the data + # add 1 to ensure maximum value is included in histogram -> range(max_bin) + max_bin += 1 - cdef dtype_t * out_data = out.data - cdef dtype_t * image_data = image.data - cdef dtype_t * mask_data = mask.data + cdef Py_ssize_t mid_bin = max_bin / 2 + + # define pointers to the data + cdef char* mask_data = &mask[0, 0] # 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 + # number of pixels actually inside the neighborhood (double) + cdef double pop = 0 # the current local histogram distribution - cdef Py_ssize_t * histo = malloc(256 * sizeof(Py_ssize_t)) + cdef Py_ssize_t* histo = malloc(max_bin * sizeof(Py_ssize_t)) + for i in range(max_bin): + histo[i] = 0 # 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 + 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 + num_se_n = num_se_s = num_se_e = num_se_w = 0 + + 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) < 0 + cdef char[:, :] t_e = (np.diff(t, axis=1) < 0).view(np.uint8) t = np.hstack((np.zeros((selem.shape[0], 1)), selem)) - t_w = np.diff(t, axis=1) > 0 + cdef char[:, :] t_w = (np.diff(t, axis=1) > 0).view(np.uint8) t = np.vstack((selem, np.zeros((1, selem.shape[1])))) - t_s = np.diff(t, axis=0) < 0 + cdef char[:, :] t_s = (np.diff(t, axis=0) < 0).view(np.uint8) t = np.vstack((np.zeros((1, selem.shape[1])), selem)) - t_n = np.diff(t, axis=0) > 0 - - num_se_n = num_se_s = num_se_e = num_se_w = 0 + cdef char[:, :] t_n = (np.diff(t, axis=0) > 0).view(np.uint8) for r in range(srows): for c in range(scols): @@ -136,92 +141,41 @@ cdef void _core8(dtype_t kernel(Py_ssize_t *, float, dtype_t, float, 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]) + histogram_increment(histo, &pop, image[rr, cc]) r = 0 c = 0 - # kernel ------------------------------------------------------------------- - out_data[r * cols + c] = kernel(histo, pop, image_data[r * cols + c], + out[r, c] = kernel(histo, pop, image[r, c], max_bin, mid_bin, 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]) + histogram_increment(histo, &pop, image[rr, 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]) + histogram_decrement(histo, &pop, image[rr, cc]) - # kernel ----------------------------------------------------------- - out_data[r * cols + c] = \ - kernel(histo, pop, image_data[r * cols + c], p0, p1, s0, s1) - # kernel ----------------------------------------------------------- + out[r, c] = kernel(histo, pop, image[r, c], + max_bin, mid_bin, p0, p1, s0, s1) - 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 + r += 1 # pass to the next row if r >= rows: break @@ -230,21 +184,55 @@ cdef void _core8(dtype_t kernel(Py_ssize_t *, float, dtype_t, float, 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]) + histogram_increment(histo, &pop, image[rr, 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]) + histogram_decrement(histo, &pop, image[rr, cc]) - # kernel --------------------------------------------------------------- - out_data[r * cols + c] = kernel(histo, pop, image_data[r * cols + c], - p0, p1, s0, s1) - # kernel --------------------------------------------------------------- + out[r, c] = kernel(histo, pop, image[r, c], + max_bin, mid_bin, p0, p1, s0, s1) + + # ---> 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[rr, 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[rr, cc]) + + out[r, c] = kernel(histo, pop, image[r, c], + max_bin, mid_bin, p0, p1, s0, s1) + + 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[rr, 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[rr, cc]) + + out[r, c] = kernel(histo, pop, image[r, c], + max_bin, mid_bin, p0, p1, s0, s1) # release memory allocated by malloc - free(se_e_r) free(se_e_c) free(se_w_r) @@ -253,5 +241,4 @@ cdef void _core8(dtype_t kernel(Py_ssize_t *, float, dtype_t, float, free(se_n_c) free(se_s_r) free(se_s_c) - free(histo) diff --git a/skimage/filter/rank/generic.py b/skimage/filter/rank/generic.py index 94fc3130..0bc8ca50 100644 --- a/skimage/filter/rank/generic.py +++ b/skimage/filter/rank/generic.py @@ -1,11 +1,740 @@ +"""The local histogram is computed using a sliding window similar to the method +described in [1]_. + +Input image can be 8-bit or 16-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-/16-bit or double with respect to the input image and the +rank filter operation. + +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. + +""" + +import warnings import numpy as np +from skimage import img_as_ubyte + +from . import generic_cy -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)) +__all__ = ['autolevel', 'bottomhat', 'equalize', 'gradient', 'maximum', 'mean', + 'subtract_mean', 'median', 'minimum', 'modal', 'enhance_contrast', + 'pop', 'threshold', 'tophat', 'noise_filter', 'entropy', 'otsu'] + + +def _handle_input(image, selem, out, mask, out_dtype=None): + + if image.dtype not in (np.uint8, np.uint16): + image = img_as_ubyte(image) + + selem = np.ascontiguousarray(img_as_ubyte(selem > 0)) + image = np.ascontiguousarray(image) + + if mask is None: + mask = np.ones(image.shape, dtype=np.uint8) else: - return 1 + mask = img_as_ubyte(mask) + mask = np.ascontiguousarray(mask) + + if out is None: + if out_dtype is None: + out_dtype = image.dtype + out = np.empty_like(image, dtype=out_dtype) + + if image is out: + raise NotImplementedError("Cannot perform rank operation in place.") + + is_8bit = image.dtype in (np.uint8, np.int8) + + if is_8bit: + max_bin = 255 + else: + max_bin = max(4, image.max()) + + bitdepth = int(np.log2(max_bin)) + if bitdepth > 10: + warnings.warn("Bitdepth of %d may result in bad rank filter " + "performance due to large number of bins." % bitdepth) + + return image, selem, out, mask, max_bin + + +def _apply(func, image, selem, out, mask, shift_x, shift_y, out_dtype=None): + + image, selem, out, mask, max_bin = _handle_input(image, selem, out, mask, + out_dtype) + + func(image, selem, shift_x=shift_x, shift_y=shift_y, mask=mask, + out=out, max_bin=max_bin) + + return out + + +def autolevel(image, selem, out=None, mask=None, shift_x=False, shift_y=False): + """Autolevel image using local histogram. + + Parameters + ---------- + image : ndarray (uint8, uint16) + Image array. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray (same dtype as input) + If None, a new array will be allocated. + mask : ndarray + 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 : ndarray (same dtype as input image) + Output image. + + 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(generic_cy._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 (uint8, uint16) + Image array. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray (same dtype as input) + If None, a new array will be allocated. + mask : ndarray + 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 + ------- + bottomhat : ndarray (same dtype as input image) + The result of the local bottomhat. + + """ + + return _apply(generic_cy._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 (uint8, uint16) + Image array. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray (same dtype as input) + If None, a new array will be allocated. + mask : ndarray + 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 : ndarray (same dtype as input image) + Output image. + + 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(generic_cy._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 (uint8, uint16) + Image array. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray (same dtype as input) + If None, a new array will be allocated. + mask : ndarray + 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 : ndarray (same dtype as input image) + Output image. + + """ + + return _apply(generic_cy._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 (uint8, uint16) + Image array. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray (same dtype as input) + If None, a new array will be allocated. + mask : ndarray + 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 : ndarray (same dtype as input image) + Output image. + + See also + -------- + skimage.morphology.dilation + + Note + ---- + * the lower algorithm complexity makes the rank.maximum() more efficient for + larger images and structuring elements + + """ + + return _apply(generic_cy._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 (uint8, uint16) + Image array. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray (same dtype as input) + If None, a new array will be allocated. + mask : ndarray + 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 : ndarray (same dtype as input image) + Output image. + + 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(generic_cy._mean, image, selem, out=out, + mask=mask, shift_x=shift_x, shift_y=shift_y) + + +def subtract_mean(image, selem, out=None, mask=None, shift_x=False, + shift_y=False): + """Return image subtracted from its local mean. + + Parameters + ---------- + image : ndarray (uint8, uint16) + Image array. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray (same dtype as input) + If None, a new array will be allocated. + mask : ndarray + 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 : ndarray (same dtype as input image) + Output image. + + """ + + return _apply(generic_cy._subtract_mean, 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 (uint8, uint16) + Image array. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray (same dtype as input) + If None, a new array will be allocated. + mask : ndarray + 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 : ndarray (same dtype as input image) + Output image. + + 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(generic_cy._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 (uint8, uint16) + Image array. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray (same dtype as input) + If None, a new array will be allocated. + mask : ndarray + 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 : ndarray (same dtype as input image) + Output image. + + See also + -------- + skimage.morphology.erosion + + Note + ---- + * the lower algorithm complexity makes the rank.minimum() more efficient + for larger images and structuring elements + + """ + + return _apply(generic_cy._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 (uint8, uint16) + Image array. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray (same dtype as input) + If None, a new array will be allocated. + mask : ndarray + 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 : ndarray (same dtype as input image) + Output image. + + """ + + return _apply(generic_cy._modal, image, selem, + out=out, mask=mask, shift_x=shift_x, shift_y=shift_y) + + +def enhance_contrast(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 (uint8, uint16) + Image array. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray (same dtype as input) + If None, a new array will be allocated. + mask : ndarray + 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 + Output image. + out : ndarray (same dtype as input image) + The result of the local enhance_contrast. + + Examples + -------- + >>> from skimage import data + >>> from skimage.morphology import disk + >>> from skimage.filter.rank import enhance_contrast + >>> # Load test image + >>> ima = data.camera() + >>> # Local mean + >>> avg = enhance_contrast(ima, disk(20)) + + """ + + return _apply(generic_cy._enhance_contrast, 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 (uint8, uint16) + Image array. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray (same dtype as input) + If None, a new array will be allocated. + mask : ndarray + 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 : ndarray (same dtype as input image) + Output image. + + 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(generic_cy._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 (uint8, uint16) + Image array. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray (same dtype as input) + If None, a new array will be allocated. + mask : ndarray + 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 : ndarray (same dtype as input image) + Output image. + + 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(generic_cy._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 (uint8, uint16) + Image array. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray (same dtype as input) + If None, a new array will be allocated. + mask : ndarray + 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 : ndarray (same dtype as input image) + Output image. + + """ + + return _apply(generic_cy._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 (uint8, uint16) + Image array. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray (same dtype as input) + If None, a new array will be allocated. + mask : ndarray + 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 : ndarray (same dtype as input image) + Output image. + + """ + + # 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(generic_cy._noise_filter, 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 [1]_ 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 (uint8, uint16) + Image array. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray (same dtype as input) + If None, a new array will be allocated. + mask : ndarray + 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 : ndarray (double) + Output image. + + References + ---------- + .. [1] 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 + >>> a8 = data.camera() + >>> ent8 = entropy(a8, disk(5)) + + """ + + return _apply(generic_cy._entropy, image, selem, + out=out, mask=mask, shift_x=shift_x, shift_y=shift_y, + out_dtype=np.double) + + +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). + 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 + 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 : ndarray (same dtype as input image) + Output image. + + 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-bit test images + >>> a8 = data.camera() + >>> loc_otsu = otsu(a8, disk(5)) + >>> thresh_image = a8 >= loc_otsu + + """ + + return _apply(generic_cy._otsu, image, selem, out=out, + mask=mask, shift_x=shift_x, shift_y=shift_y) diff --git a/skimage/filter/rank/generic_cy.pyx b/skimage/filter/rank/generic_cy.pyx new file mode 100644 index 00000000..dcf6e361 --- /dev/null +++ b/skimage/filter/rank/generic_cy.pyx @@ -0,0 +1,506 @@ +#cython: cdivision=True +#cython: boundscheck=False +#cython: nonecheck=False +#cython: wraparound=False + +cimport numpy as cnp +from libc.math cimport log + +from .core_cy cimport dtype_t, dtype_t_out, _core + + +cdef inline double _kernel_autolevel(Py_ssize_t* histo, double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): + + cdef Py_ssize_t i, imin, imax, delta + + if pop: + for i in range(max_bin - 1, -1, -1): + if histo[i]: + imax = i + break + for i in range(max_bin): + if histo[i]: + imin = i + break + delta = imax - imin + if delta > 0: + return (max_bin - 1) * (g - imin) / delta + else: + return 0 + else: + return 0 + + +cdef inline double _kernel_bottomhat(Py_ssize_t* histo, double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): + + cdef Py_ssize_t i + + if pop: + for i in range(max_bin): + if histo[i]: + break + return g - i + else: + return 0 + + +cdef inline double _kernel_equalize(Py_ssize_t* histo, double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): + + cdef Py_ssize_t i + cdef Py_ssize_t sum = 0 + + if pop: + for i in range(max_bin): + sum += histo[i] + if i >= g: + break + return ((max_bin - 1) * sum) / pop + else: + return 0 + + +cdef inline double _kernel_gradient(Py_ssize_t* histo, double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): + + cdef Py_ssize_t i, imin, imax + + if pop: + for i in range(max_bin - 1, -1, -1): + if histo[i]: + imax = i + break + for i in range(max_bin): + if histo[i]: + imin = i + break + return imax - imin + else: + return 0 + + +cdef inline double _kernel_maximum(Py_ssize_t* histo, double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): + + cdef Py_ssize_t i + + if pop: + for i in range(max_bin - 1, -1, -1): + if histo[i]: + return i + else: + return 0 + + +cdef inline double _kernel_mean(Py_ssize_t* histo, double pop,dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): + + cdef Py_ssize_t i + cdef Py_ssize_t mean = 0 + + if pop: + for i in range(max_bin): + mean += histo[i] * i + return mean / pop + else: + return 0 + + +cdef inline double _kernel_subtract_mean(Py_ssize_t* histo, double pop, + dtype_t g, + Py_ssize_t max_bin, + Py_ssize_t mid_bin, double p0, + double p1, Py_ssize_t s0, + Py_ssize_t s1): + + cdef Py_ssize_t i + cdef Py_ssize_t mean = 0 + + if pop: + for i in range(max_bin): + mean += histo[i] * i + return (g - mean / pop) / 2. + 127 + else: + return 0 + + +cdef inline double _kernel_median(Py_ssize_t* histo, double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): + + cdef Py_ssize_t i + cdef double sum = pop / 2.0 + + if pop: + for i in range(max_bin): + if histo[i]: + sum -= histo[i] + if sum < 0: + return i + else: + return 0 + + +cdef inline double _kernel_minimum(Py_ssize_t* histo, double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): + + cdef Py_ssize_t i + + if pop: + for i in range(max_bin): + if histo[i]: + return i + else: + return 0 + + +cdef inline double _kernel_modal(Py_ssize_t* histo, double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): + + cdef Py_ssize_t hmax = 0, imax = 0 + + if pop: + for i in range(max_bin): + if histo[i] > hmax: + hmax = histo[i] + imax = i + return imax + else: + return 0 + + +cdef inline double _kernel_enhance_contrast(Py_ssize_t* histo, double pop, + dtype_t g, + Py_ssize_t max_bin, + Py_ssize_t mid_bin, double p0, + double p1, Py_ssize_t s0, + Py_ssize_t s1): + + cdef Py_ssize_t i, imin, imax + + if pop: + for i in range(max_bin - 1, -1, -1): + if histo[i]: + imax = i + break + for i in range(max_bin): + if histo[i]: + imin = i + break + if imax - g < g - imin: + return imax + else: + return imin + else: + return 0 + + +cdef inline double _kernel_pop(Py_ssize_t* histo, double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): + + return pop + + +cdef inline double _kernel_threshold(Py_ssize_t* histo, double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): + + cdef Py_ssize_t i + cdef Py_ssize_t mean = 0 + + if pop: + for i in range(max_bin): + mean += histo[i] * i + return g > (mean / pop) + else: + return 0 + + +cdef inline double _kernel_tophat(Py_ssize_t* histo, double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): + + cdef Py_ssize_t i + + if pop: + for i in range(max_bin - 1, -1, -1): + if histo[i]: + break + return i - g + else: + return 0 + + +cdef inline double _kernel_noise_filter(Py_ssize_t* histo, double pop, + dtype_t g, Py_ssize_t max_bin, + Py_ssize_t mid_bin, double p0, + double 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 0 + + for i in range(g, -1, -1): + if histo[i]: + break + min_i = g - i + for i in range(g, max_bin): + if histo[i]: + break + if i - g < min_i: + return i - g + else: + return min_i + + +cdef inline double _kernel_entropy(Py_ssize_t* histo, double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): + cdef Py_ssize_t i + cdef double e, p + + if pop: + e = 0. + for i in range(max_bin): + p = histo[i] / pop + if p > 0: + e -= p * log(p) / 0.6931471805599453 + return e + else: + return 0 + + +cdef inline double _kernel_otsu(Py_ssize_t* histo, double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): + cdef Py_ssize_t i + cdef Py_ssize_t max_i + cdef double P, mu1, mu2, q1, new_q1, sigma_b, max_sigma_b + cdef double mu = 0. + + # compute local mean + if pop: + for i in range(max_bin): + mu += histo[i] * i + mu = mu / pop + else: + return 0 + + # maximizing the between class variance + max_i = 0 + q1 = histo[0] / pop + mu1 = 0. + max_sigma_b = 0. + + for i in range(1, max_bin): + 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 max_i + + +def _autolevel(dtype_t[:, ::1] image, + char[:, ::1] selem, + char[:, ::1] mask, + dtype_t_out[:, ::1] out, + char shift_x, char shift_y, Py_ssize_t max_bin): + + _core(_kernel_autolevel[dtype_t], image, selem, mask, out, + shift_x, shift_y, 0, 0, 0, 0, max_bin) + + +def _bottomhat(dtype_t[:, ::1] image, + char[:, ::1] selem, + char[:, ::1] mask, + dtype_t_out[:, ::1] out, + char shift_x, char shift_y, Py_ssize_t max_bin): + + _core(_kernel_bottomhat[dtype_t], image, selem, mask, out, + shift_x, shift_y, 0, 0, 0, 0, max_bin) + + +def _equalize(dtype_t[:, ::1] image, + char[:, ::1] selem, + char[:, ::1] mask, + dtype_t_out[:, ::1] out, + char shift_x, char shift_y, Py_ssize_t max_bin): + + _core(_kernel_equalize[dtype_t], image, selem, mask, out, + shift_x, shift_y, 0, 0, 0, 0, max_bin) + + +def _gradient(dtype_t[:, ::1] image, + char[:, ::1] selem, + char[:, ::1] mask, + dtype_t_out[:, ::1] out, + char shift_x, char shift_y, Py_ssize_t max_bin): + + _core(_kernel_gradient[dtype_t], image, selem, mask, out, + shift_x, shift_y, 0, 0, 0, 0, max_bin) + + +def _maximum(dtype_t[:, ::1] image, + char[:, ::1] selem, + char[:, ::1] mask, + dtype_t_out[:, ::1] out, + char shift_x, char shift_y, Py_ssize_t max_bin): + + _core(_kernel_maximum[dtype_t], image, selem, mask, out, + shift_x, shift_y, 0, 0, 0, 0, max_bin) + + +def _mean(dtype_t[:, ::1] image, + char[:, ::1] selem, + char[:, ::1] mask, + dtype_t_out[:, ::1] out, + char shift_x, char shift_y, Py_ssize_t max_bin): + + _core(_kernel_mean[dtype_t], image, selem, mask, out, + shift_x, shift_y, 0, 0, 0, 0, max_bin) + + +def _subtract_mean(dtype_t[:, ::1] image, + char[:, ::1] selem, + char[:, ::1] mask, + dtype_t_out[:, ::1] out, + char shift_x, char shift_y, Py_ssize_t max_bin): + + _core(_kernel_subtract_mean[dtype_t], image, selem, mask, + out, shift_x, shift_y, 0, 0, 0, 0, max_bin) + + +def _median(dtype_t[:, ::1] image, + char[:, ::1] selem, + char[:, ::1] mask, + dtype_t_out[:, ::1] out, + char shift_x, char shift_y, Py_ssize_t max_bin): + + _core(_kernel_median[dtype_t], image, selem, mask, out, + shift_x, shift_y, 0, 0, 0, 0, max_bin) + + +def _minimum(dtype_t[:, ::1] image, + char[:, ::1] selem, + char[:, ::1] mask, + dtype_t_out[:, ::1] out, + char shift_x, char shift_y, Py_ssize_t max_bin): + + _core(_kernel_minimum[dtype_t], image, selem, mask, out, + shift_x, shift_y, 0, 0, 0, 0, max_bin) + + +def _enhance_contrast(dtype_t[:, ::1] image, + char[:, ::1] selem, + char[:, ::1] mask, + dtype_t_out[:, ::1] out, + char shift_x, char shift_y, Py_ssize_t max_bin): + + _core(_kernel_enhance_contrast[dtype_t], image, selem, mask, + out, shift_x, shift_y, 0, 0, 0, 0, max_bin) + + +def _modal(dtype_t[:, ::1] image, + char[:, ::1] selem, + char[:, ::1] mask, + dtype_t_out[:, ::1] out, + char shift_x, char shift_y, Py_ssize_t max_bin): + + _core(_kernel_modal[dtype_t], image, selem, mask, out, + shift_x, shift_y, 0, 0, 0, 0, max_bin) + + +def _pop(dtype_t[:, ::1] image, + char[:, ::1] selem, + char[:, ::1] mask, + dtype_t_out[:, ::1] out, + char shift_x, char shift_y, Py_ssize_t max_bin): + + _core(_kernel_pop[dtype_t], image, selem, mask, out, + shift_x, shift_y, 0, 0, 0, 0, max_bin) + + +def _threshold(dtype_t[:, ::1] image, + char[:, ::1] selem, + char[:, ::1] mask, + dtype_t_out[:, ::1] out, + char shift_x, char shift_y, Py_ssize_t max_bin): + + _core(_kernel_threshold[dtype_t], image, selem, mask, out, + shift_x, shift_y, 0, 0, 0, 0, max_bin) + + +def _tophat(dtype_t[:, ::1] image, + char[:, ::1] selem, + char[:, ::1] mask, + dtype_t_out[:, ::1] out, + char shift_x, char shift_y, Py_ssize_t max_bin): + + _core(_kernel_tophat[dtype_t], image, selem, mask, out, + shift_x, shift_y, 0, 0, 0, 0, max_bin) + + +def _noise_filter(dtype_t[:, ::1] image, + char[:, ::1] selem, + char[:, ::1] mask, + dtype_t_out[:, ::1] out, + char shift_x, char shift_y, Py_ssize_t max_bin): + + _core(_kernel_noise_filter[dtype_t], image, selem, mask, out, + shift_x, shift_y, 0, 0, 0, 0, max_bin) + + +def _entropy(dtype_t[:, ::1] image, + char[:, ::1] selem, + char[:, ::1] mask, + dtype_t_out[:, ::1] out, + char shift_x, char shift_y, Py_ssize_t max_bin): + + _core(_kernel_entropy[dtype_t], image, selem, mask, out, + shift_x, shift_y, 0, 0, 0, 0, max_bin) + + +def _otsu(dtype_t[:, ::1] image, + char[:, ::1] selem, + char[:, ::1] mask, + dtype_t_out[:, ::1] out, + char shift_x, char shift_y, Py_ssize_t max_bin): + + _core(_kernel_otsu[dtype_t], image, selem, mask, out, + shift_x, shift_y, 0, 0, 0, 0, max_bin) diff --git a/skimage/filter/rank/percentile_rank.pyx b/skimage/filter/rank/percentile.py similarity index 53% rename from skimage/filter/rank/percentile_rank.pyx rename to skimage/filter/rank/percentile.py index 704f53c2..ff3b1559 100644 --- a/skimage/filter/rank/percentile_rank.pyx +++ b/skimage/filter/rank/percentile.py @@ -1,17 +1,17 @@ """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] +``autolevel_percentile`` 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]_. -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. +Input image can be 8-bit or 16-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. +Result image is 8-/16-bit or double with respect to the input image and the +rank filter operation. References ---------- @@ -23,55 +23,31 @@ References """ 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 + +from . import percentile_cy +from .generic import _handle_input -__all__ = ['percentile_autolevel', 'percentile_gradient', - 'percentile_mean', 'percentile_mean_subtraction', - 'percentile_morph_contr_enh', 'percentile', 'percentile_pop', - 'percentile_threshold'] +__all__ = ['autolevel_percentile', 'gradient_percentile', + 'mean_percentile', 'subtract_mean_percentile', + 'enhance_contrast_percentile', 'percentile', 'pop_percentile', + 'threshold_percentile'] -def _apply(func8, func16, image, selem, out, mask, shift_x, shift_y, p0, p1): - selem = img_as_ubyte(selem > 0) - image = np.ascontiguousarray(image) +def _apply(func, image, selem, out, mask, shift_x, shift_y, p0, p1, + out_dtype=None): - if mask is None: - mask = np.ones(image.shape, dtype=np.uint8) - else: - mask = np.ascontiguousarray(mask) - mask = img_as_ubyte(mask) + image, selem, out, mask, max_bin = _handle_input(image, selem, out, mask, + out_dtype) - 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.") + func(image, selem, shift_x=shift_x, shift_y=shift_y, mask=mask, + out=out, max_bin=max_bin, p0=p0, p1=p1) return out -def percentile_autolevel(image, selem, out=None, mask=None, shift_x=False, - shift_y=False, p0=.0, p1=1.): +def autolevel_percentile(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 @@ -79,15 +55,13 @@ def percentile_autolevel(image, selem, out=None, mask=None, shift_x=False, 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. + image : ndarray (uint8, uint16) + Image array. selem : ndarray The neighborhood expressed as a 2-D array of 1's and 0's. - out : ndarray + out : ndarray (same dtype as input) If None, a new array will be allocated. - mask : ndarray (uint8) + mask : ndarray 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 @@ -100,59 +74,56 @@ def percentile_autolevel(image, selem, out=None, mask=None, shift_x=False, Returns ------- - local autolevel : uint8 array or uint16 - The result of the local autolevel. + out : ndarray (same dtype as input image) + Output image. """ - 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, + return _apply(percentile_cy._autolevel, 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.): +def gradient_percentile(image, selem, out=None, mask=None, shift_x=False, + shift_y=False, p0=0, p1=1): + """Return greyscale local gradient of an image. + + gradient is computed on the given structuring element. Only + levels between percentiles [p0, p1] are used. + + Parameters + ---------- + image : ndarray (uint8, uint16) + Image array. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray (same dtype as input) + If None, a new array will be allocated. + mask : ndarray + 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 + ------- + out : ndarray (same dtype as input image) + Output image. + + """ + + return _apply(percentile_cy._gradient, + image, selem, out=out, mask=mask, shift_x=shift_x, + shift_y=shift_y, p0=p0, p1=p1) + + +def mean_percentile(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 @@ -160,15 +131,13 @@ def percentile_mean(image, selem, out=None, mask=None, shift_x=False, 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. + image : ndarray (uint8, uint16) + Image array. selem : ndarray The neighborhood expressed as a 2-D array of 1's and 0's. - out : ndarray + out : ndarray (same dtype as input) If None, a new array will be allocated. - mask : ndarray (uint8) + mask : ndarray 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 @@ -181,34 +150,32 @@ def percentile_mean(image, selem, out=None, mask=None, shift_x=False, Returns ------- - local mean : uint8 array or uint16 - The result of the local mean. + out : ndarray (same dtype as input image) + Output image. """ - return _apply(_crank8_percentiles.mean, _crank16_percentiles.mean, + return _apply(percentile_cy._mean, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y, p0=p0, p1=p1) -def percentile_mean_subtraction(image, selem, out=None, mask=None, - shift_x=False, shift_y=False, p0=.0, p1=1.): - """Return greyscale local mean_subtraction of an image. +def subtract_mean_percentile(image, selem, out=None, mask=None, + shift_x=False, shift_y=False, p0=0, p1=1): + """Return greyscale local subtract_mean of an image. - mean_subtraction is computed on the given structuring element. Only levels + subtract_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. + image : ndarray (uint8, uint16) + Image array. selem : ndarray The neighborhood expressed as a 2-D array of 1's and 0's. - out : ndarray + out : ndarray (same dtype as input) If None, a new array will be allocated. - mask : ndarray (uint8) + mask : ndarray 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 @@ -221,36 +188,32 @@ def percentile_mean_subtraction(image, selem, out=None, mask=None, Returns ------- - local mean_subtraction : uint8 array or uint16 - The result of the local mean_subtraction. + out : ndarray (same dtype as input image) + Output image. """ - return _apply(_crank8_percentiles.mean_subtraction, - _crank16_percentiles.mean_subtraction, + return _apply(percentile_cy._subtract_mean, 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. +def enhance_contrast_percentile(image, selem, out=None, mask=None, + shift_x=False, shift_y=False, p0=0, p1=1): + """Return greyscale local enhance_contrast of an image. - morph_contr_enh is computed on the given structuring element. Only levels + enhance_contrast 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. + image : ndarray (uint8, uint16) + Image array. selem : ndarray The neighborhood expressed as a 2-D array of 1's and 0's. - out : ndarray + out : ndarray (same dtype as input) If None, a new array will be allocated. - mask : ndarray (uint8) + mask : ndarray 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 @@ -263,19 +226,18 @@ def percentile_morph_contr_enh( Returns ------- - local morph_contr_enh : uint8 array or uint16 - The result of the local morph_contr_enh. + out : ndarray (same dtype as input image) + Output image. """ - return _apply(_crank8_percentiles.morph_contr_enh, - _crank16_percentiles.morph_contr_enh, + return _apply(percentile_cy._enhance_contrast, 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): + p0=0): """Return greyscale local percentile of an image. percentile is computed on the given structuring element. Returns the value @@ -283,15 +245,13 @@ def percentile(image, selem, out=None, mask=None, shift_x=False, shift_y=False, 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. + image : ndarray (uint8, uint16) + Image array. selem : ndarray The neighborhood expressed as a 2-D array of 1's and 0's. - out : ndarray + out : ndarray (same dtype as input) If None, a new array will be allocated. - mask : ndarray (uint8) + mask : ndarray 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 @@ -303,19 +263,18 @@ def percentile(image, selem, out=None, mask=None, shift_x=False, shift_y=False, Returns ------- - local percentile : uint8 array or uint16 - The result of the local percentile. + out : ndarray (same dtype as input image) + Output image. """ - return _apply(_crank8_percentiles.percentile, - _crank16_percentiles.percentile, + return _apply(percentile_cy._percentile, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y, p0=p0, p1=0.) -def percentile_pop(image, selem, out=None, mask=None, shift_x=False, - shift_y=False, p0=.0, p1=1.): +def pop_percentile(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 @@ -323,15 +282,13 @@ def percentile_pop(image, selem, out=None, mask=None, shift_x=False, 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. + image : ndarray (uint8, uint16) + Image array. selem : ndarray The neighborhood expressed as a 2-D array of 1's and 0's. - out : ndarray + out : ndarray (same dtype as input) If None, a new array will be allocated. - mask : ndarray (uint8) + mask : ndarray 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 @@ -344,18 +301,18 @@ def percentile_pop(image, selem, out=None, mask=None, shift_x=False, Returns ------- - local pop : uint8 array or uint16 - The result of the local pop. + out : ndarray (same dtype as input image) + Output image. """ - return _apply(_crank8_percentiles.pop, _crank16_percentiles.pop, + return _apply(percentile_cy._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): +def threshold_percentile(image, selem, out=None, mask=None, shift_x=False, + shift_y=False, p0=0): """Return greyscale local threshold of an image. threshold is computed on the given structuring element. Returns @@ -365,15 +322,13 @@ def percentile_threshold(image, selem, out=None, mask=None, shift_x=False, 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. + image : ndarray (uint8, uint16) + Image array. selem : ndarray The neighborhood expressed as a 2-D array of 1's and 0's. - out : ndarray + out : ndarray (same dtype as input) If None, a new array will be allocated. - mask : ndarray (uint8) + mask : ndarray 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 @@ -383,14 +338,13 @@ def percentile_threshold(image, selem, out=None, mask=None, shift_x=False, p0 : float in [0, ..., 1] Set the percentile value. - Returns - ------- - local threshold : uint8 array or uint16 + out : ndarray (same dtype as input image) + Output image. + local threshold : ndarray (same dtype as input) 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=0.) + return _apply(percentile_cy._threshold, + image, selem, out=out, mask=mask, shift_x=shift_x, + shift_y=shift_y, p0=p0, p1=0) diff --git a/skimage/filter/rank/percentile_cy.pyx b/skimage/filter/rank/percentile_cy.pyx new file mode 100644 index 00000000..e951a76e --- /dev/null +++ b/skimage/filter/rank/percentile_cy.pyx @@ -0,0 +1,301 @@ +#cython: cdivision=True +#cython: boundscheck=False +#cython: nonecheck=False +#cython: wraparound=False + +cimport numpy as cnp +from .core_cy cimport dtype_t, dtype_t_out, _core, _min, _max + + +cdef inline double _kernel_autolevel(Py_ssize_t* histo, double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): + + cdef Py_ssize_t i, imin, imax, sum, delta + + if pop: + sum = 0 + p1 = 1.0 - p1 + for i in range(max_bin): + sum += histo[i] + if sum > p0 * pop: + imin = i + break + sum = 0 + for i in range(max_bin - 1, -1, -1): + sum += histo[i] + if sum > p1 * pop: + imax = i + break + + delta = imax - imin + if delta > 0: + return (max_bin - 1) * (_min(_max(imin, g), imax) + - imin) / delta + else: + return imax - imin + else: + return 0 + + +cdef inline double _kernel_gradient(Py_ssize_t* histo, double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): + + cdef Py_ssize_t i, imin, imax, sum, delta + + if pop: + sum = 0 + p1 = 1.0 - p1 + for i in range(max_bin): + sum += histo[i] + if sum >= p0 * pop: + imin = i + break + sum = 0 + for i in range(max_bin - 1, -1, -1): + sum += histo[i] + if sum >= p1 * pop: + imax = i + break + + return imax - imin + else: + return 0 + + +cdef inline double _kernel_mean(Py_ssize_t* histo, double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): + + cdef Py_ssize_t i, sum, mean, n + + if pop: + sum = 0 + mean = 0 + n = 0 + for i in range(max_bin): + sum += histo[i] + if (sum >= p0 * pop) and (sum <= p1 * pop): + n += histo[i] + mean += histo[i] * i + + if n > 0: + return mean / n + else: + return 0 + else: + return 0 + + +cdef inline double _kernel_subtract_mean(Py_ssize_t* histo, double pop, + dtype_t g, + Py_ssize_t max_bin, + Py_ssize_t mid_bin, double p0, + double p1, Py_ssize_t s0, + Py_ssize_t s1): + + cdef Py_ssize_t i, sum, mean, n + + if pop: + sum = 0 + mean = 0 + n = 0 + for i in range(max_bin): + sum += histo[i] + if (sum >= p0 * pop) and (sum <= p1 * pop): + n += histo[i] + mean += histo[i] * i + if n > 0: + return (g - (mean / n)) * .5 + mid_bin + else: + return 0 + else: + return 0 + + +cdef inline double _kernel_enhance_contrast(Py_ssize_t* histo, double pop, + dtype_t g, + Py_ssize_t max_bin, + Py_ssize_t mid_bin, double p0, + double p1, Py_ssize_t s0, + Py_ssize_t s1): + + cdef Py_ssize_t i, imin, imax, sum, delta + + if pop: + sum = 0 + p1 = 1.0 - p1 + for i in range(max_bin): + sum += histo[i] + if sum > p0 * pop: + imin = i + break + sum = 0 + for i in range(max_bin - 1, -1, -1): + sum += histo[i] + if sum > p1 * pop: + imax = i + break + if g > imax: + return imax + if g < imin: + return imin + if imax - g < g - imin: + return imax + else: + return imin + else: + return 0 + + +cdef inline double _kernel_percentile(Py_ssize_t* histo, double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): + + cdef Py_ssize_t i + cdef Py_ssize_t sum = 0 + + if pop: + if p0 == 1: # make sure p0 = 1 returns the maximum filter + for i in range(max_bin - 1, -1, -1): + if histo[i]: + break + else: + for i in range(max_bin): + sum += histo[i] + if sum > p0 * pop: + break + return i + else: + return 0 + + +cdef inline double _kernel_pop(Py_ssize_t* histo, double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): + + cdef Py_ssize_t i, sum, n + + if pop: + sum = 0 + n = 0 + for i in range(max_bin): + sum += histo[i] + if (sum >= p0 * pop) and (sum <= p1 * pop): + n += histo[i] + return n + else: + return 0 + + +cdef inline double _kernel_threshold(Py_ssize_t* histo, double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1): + + cdef int i + cdef Py_ssize_t sum = 0 + + if pop: + for i in range(max_bin): + sum += histo[i] + if sum >= p0 * pop: + break + + return (max_bin - 1) * (g >= i) + else: + return 0 + + +def _autolevel(dtype_t[:, ::1] image, + char[:, ::1] selem, + char[:, ::1] mask, + dtype_t_out[:, ::1] out, + char shift_x, char shift_y, double p0, double p1, + Py_ssize_t max_bin): + + _core(_kernel_autolevel[dtype_t], image, selem, mask, out, + shift_x, shift_y, p0, p1, 0, 0, max_bin) + + +def _gradient(dtype_t[:, ::1] image, + char[:, ::1] selem, + char[:, ::1] mask, + dtype_t_out[:, ::1] out, + char shift_x, char shift_y, double p0, double p1, + Py_ssize_t max_bin): + + _core(_kernel_gradient[dtype_t], image, selem, mask, out, + shift_x, shift_y, p0, p1, 0, 0, max_bin) + + +def _mean(dtype_t[:, ::1] image, + char[:, ::1] selem, + char[:, ::1] mask, + dtype_t_out[:, ::1] out, + char shift_x, char shift_y, double p0, double p1, + Py_ssize_t max_bin): + + _core(_kernel_mean[dtype_t], image, selem, mask, out, + shift_x, shift_y, p0, p1, 0, 0, max_bin) + + +def _subtract_mean(dtype_t[:, ::1] image, + char[:, ::1] selem, + char[:, ::1] mask, + dtype_t_out[:, ::1] out, + char shift_x, char shift_y, double p0, double p1, + Py_ssize_t max_bin): + + _core(_kernel_subtract_mean[dtype_t], image, selem, mask, + out, shift_x, shift_y, p0, p1, 0, 0, max_bin) + + +def _enhance_contrast(dtype_t[:, ::1] image, + char[:, ::1] selem, + char[:, ::1] mask, + dtype_t_out[:, ::1] out, + char shift_x, char shift_y, double p0, double p1, + Py_ssize_t max_bin): + + _core(_kernel_enhance_contrast[dtype_t], image, selem, mask, + out, shift_x, shift_y, p0, p1, 0, 0, max_bin) + + +def _percentile(dtype_t[:, ::1] image, + char[:, ::1] selem, + char[:, ::1] mask, + dtype_t_out[:, ::1] out, + char shift_x, char shift_y, double p0, double p1, + Py_ssize_t max_bin): + + _core(_kernel_percentile[dtype_t], image, selem, mask, out, + shift_x, shift_y, p0, 1, 0, 0, max_bin) + + +def _pop(dtype_t[:, ::1] image, + char[:, ::1] selem, + char[:, ::1] mask, + dtype_t_out[:, ::1] out, + char shift_x, char shift_y, double p0, double p1, + Py_ssize_t max_bin): + + _core(_kernel_pop[dtype_t], image, selem, mask, out, + shift_x, shift_y, p0, p1, 0, 0, max_bin) + + +def _threshold(dtype_t[:, ::1] image, + char[:, ::1] selem, + char[:, ::1] mask, + dtype_t_out[:, ::1] out, + char shift_x, char shift_y, double p0, double p1, + Py_ssize_t max_bin): + + _core(_kernel_threshold[dtype_t], image, selem, mask, out, + shift_x, shift_y, p0, 1, 0, 0, max_bin) diff --git a/skimage/filter/rank/tests/test_rank.py b/skimage/filter/rank/tests/test_rank.py index d39386c6..1d573b81 100644 --- a/skimage/filter/rank/tests/test_rank.py +++ b/skimage/filter/rank/tests/test_rank.py @@ -33,10 +33,10 @@ def test_random_sizes(): shift_x=+1, shift_y=+1) assert_array_equal(image16.shape, out16.shape) - rank.percentile_mean(image=image16, mask=mask, out=out16, + rank.mean_percentile(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, + rank.mean_percentile(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) @@ -51,7 +51,7 @@ def test_compare_with_cmorph_dilate(): 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) + cm = cmorph._dilate(image=image, selem=elem) assert_array_equal(out, cm) @@ -65,7 +65,7 @@ def test_compare_with_cmorph_erode(): 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) + cm = cmorph._erode(image=image, selem=elem) assert_array_equal(out, cm) @@ -78,7 +78,7 @@ def test_bitdepth(): 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, + r = rank.mean_percentile(image=image, selem=elem, mask=mask, out=out, shift_x=0, shift_y=0, p0=.1, p1=.9) @@ -130,17 +130,6 @@ def test_structuring_element8(): 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 @@ -167,36 +156,34 @@ def test_compare_autolevels(): selem = disk(20) loc_autolevel = rank.autolevel(image, selem=selem) - loc_perc_autolevel = rank.percentile_autolevel(image, selem=selem, + loc_perc_autolevel = rank.autolevel_percentile(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 + # 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, + loc_perc_autolevel = rank.autolevel_percentile(image, selem=selem, p0=.0, p1=1.) assert_array_equal(loc_autolevel, loc_perc_autolevel) -def test_compare_uint_vs_float(): - # filters applied on 8-bit image ore 16-bit image (having only real 8-bit of - # dynamic) should be identical +def test_compare_ubyte_vs_float(): # Create signed int8 image that and convert it to uint8 - image_uint = img_as_uint(data.camera()) + image_uint = img_as_ubyte(data.camera()[:50, :50]) image_float = img_as_float(image_uint) methods = ['autolevel', 'bottomhat', 'equalize', 'gradient', 'threshold', - 'meansubtraction', 'morph_contr_enh', 'pop', 'tophat'] + 'subtract_mean', 'enhance_contrast', 'pop', 'tophat'] for method in methods: func = getattr(rank, method) @@ -206,8 +193,8 @@ def test_compare_uint_vs_float(): def test_compare_8bit_unsigned_vs_signed(): - # filters applied on 8-bit image ore 16-bit image (having only real 8-bit of - # dynamic) should be identical + # filters applied on 8-bit image ore 16-bit image (having only real 8-bit + # of dynamic) should be identical # Create signed int8 image that and convert it to uint8 image = img_as_ubyte(data.camera()) @@ -218,8 +205,8 @@ def test_compare_8bit_unsigned_vs_signed(): assert_array_equal(image_u, img_as_ubyte(image_s)) methods = ['autolevel', 'bottomhat', 'equalize', 'gradient', 'maximum', - 'mean', 'meansubtraction', 'median', 'minimum', 'modal', - 'morph_contr_enh', 'pop', 'threshold', 'tophat'] + 'mean', 'subtract_mean', 'median', 'minimum', 'modal', + 'enhance_contrast', 'pop', 'threshold', 'tophat'] for method in methods: func = getattr(rank, method) @@ -229,16 +216,16 @@ def test_compare_8bit_unsigned_vs_signed(): 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 + # filters applied on 8-bit image ore 16-bit image (having only real 8-bit + # of dynamic) should be identical image8 = util.img_as_ubyte(data.camera()) image16 = image8.astype(np.uint16) assert_array_equal(image8, image16) methods = ['autolevel', 'bottomhat', 'equalize', 'gradient', 'maximum', - 'mean', 'meansubtraction', 'median', 'minimum', 'modal', - 'morph_contr_enh', 'pop', 'threshold', 'tophat'] + 'mean', 'subtract_mean', 'median', 'minimum', 'modal', + 'enhance_contrast', 'pop', 'threshold', 'tophat'] for method in methods: func = getattr(rank, method) @@ -340,7 +327,8 @@ def test_smallest_selem16(): def test_empty_selem(): - # check that min, max and mean returns zeros if structuring element is empty + # 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) @@ -367,13 +355,11 @@ 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 = 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)) + 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) @@ -385,37 +371,41 @@ def test_entropy(): 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) + assert(np.max(rank.entropy(data, selem)) == 1) # 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) + assert(np.max(rank.entropy(data, selem)) == 2) # 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) + assert(np.max(rank.entropy(data, selem)) == 3) # 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) + assert(np.max(rank.entropy(data, selem)) == 4) # 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) + assert(np.max(rank.entropy(data, selem)) == 6) # 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) + assert(np.max(rank.entropy(data, selem)) == 8) # 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) + assert(np.max(rank.entropy(data, selem)) == 12) + + # make sure output is of dtype double + out = rank.entropy(data, np.ones((16, 16), dtype=np.uint8)) + assert out.dtype == np.double def test_selem_dtypes(): @@ -433,10 +423,81 @@ def test_selem_dtypes(): rank.mean(image=image, selem=elem, out=out, mask=mask, shift_x=0, shift_y=0) assert_array_equal(image, out) - rank.percentile_mean(image=image, selem=elem, out=out, mask=mask, + rank.mean_percentile(image=image, selem=elem, out=out, mask=mask, shift_x=0, shift_y=0) assert_array_equal(image, out) +def test_16bit(): + image = np.zeros((21, 21), dtype=np.uint16) + selem = np.ones((3, 3), dtype=np.uint8) + + for bitdepth in range(17): + value = 2 ** bitdepth - 1 + image[10, 10] = value + assert rank.minimum(image, selem)[10, 10] == 0 + assert rank.maximum(image, selem)[10, 10] == value + assert rank.mean(image, selem)[10, 10] == value / selem.size + + +def test_bilateral(): + image = np.zeros((21, 21), dtype=np.uint16) + selem = np.ones((3, 3), dtype=np.uint8) + + image[10, 10] = 1000 + image[10, 11] = 1010 + image[10, 9] = 900 + + assert rank.mean_bilateral(image, selem, s0=1, s1=1)[10, 10] == 1000 + assert rank.pop_bilateral(image, selem, s0=1, s1=1)[10, 10] == 1 + assert rank.mean_bilateral(image, selem, s0=11, s1=11)[10, 10] == 1005 + assert rank.pop_bilateral(image, selem, s0=11, s1=11)[10, 10] == 2 + + +def test_percentile_min(): + # check that percentile p0 = 0 is identical to local min + img = data.camera() + img16 = img.astype(np.uint16) + selem = disk(15) + # check for 8bit + img_p0 = rank.percentile(img, selem=selem, p0=0) + img_min = rank.minimum(img, selem=selem) + assert_array_equal(img_p0, img_min) + # check for 16bit + img_p0 = rank.percentile(img16, selem=selem, p0=0) + img_min = rank.minimum(img16, selem=selem) + assert_array_equal(img_p0, img_min) + + +def test_percentile_max(): + # check that percentile p0 = 1 is identical to local max + img = data.camera() + img16 = img.astype(np.uint16) + selem = disk(15) + # check for 8bit + img_p0 = rank.percentile(img, selem=selem, p0=1.) + img_max = rank.maximum(img, selem=selem) + assert_array_equal(img_p0, img_max) + # check for 16bit + img_p0 = rank.percentile(img16, selem=selem, p0=1.) + img_max = rank.maximum(img16, selem=selem) + assert_array_equal(img_p0, img_max) + + +def test_percentile_median(): + # check that percentile p0 = 0.5 is identical to local median + img = data.camera() + img16 = img.astype(np.uint16) + selem = disk(15) + # check for 8bit + img_p0 = rank.percentile(img, selem=selem, p0=.5) + img_max = rank.median(img, selem=selem) + assert_array_equal(img_p0, img_max) + # check for 16bit + img_p0 = rank.percentile(img16, selem=selem, p0=.5) + img_max = rank.median(img16, selem=selem) + assert_array_equal(img_p0, img_max) + + if __name__ == "__main__": run_module_suite() diff --git a/skimage/filter/setup.py b/skimage/filter/setup.py index c70730f0..33ad97df 100644 --- a/skimage/filter/setup.py +++ b/skimage/filter/setup.py @@ -14,46 +14,29 @@ def configuration(parent_package='', top_path=None): cython(['_ctmf.pyx'], working_path=base_path) cython(['_denoise_cy.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/percentile_rank.pyx'], working_path=base_path) - cython(['rank/bilateral_rank.pyx'], working_path=base_path) + cython(['rank/core_cy.pyx'], working_path=base_path) + cython(['rank/generic_cy.pyx'], working_path=base_path) + cython(['rank/percentile_cy.pyx'], working_path=base_path) + cython(['rank/bilateral_cy.pyx'], working_path=base_path) config.add_extension('_ctmf', sources=['_ctmf.c'], include_dirs=[get_numpy_include_dirs()]) config.add_extension('_denoise_cy', sources=['_denoise_cy.c'], include_dirs=[get_numpy_include_dirs(), '../_shared']) - config.add_extension('rank._core8', sources=['rank/_core8.c'], + config.add_extension('rank.core_cy', sources=['rank/core_cy.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'], + config.add_extension('rank.generic_cy', sources=['rank/generic_cy.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'], + 'rank.percentile_cy', sources=['rank/percentile_cy.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.percentile_rank', sources=['rank/percentile_rank.c'], - include_dirs=[get_numpy_include_dirs()]) - config.add_extension( - 'rank.bilateral_rank', sources=['rank/bilateral_rank.c'], + 'rank.bilateral_cy', sources=['rank/bilateral_cy.c'], include_dirs=[get_numpy_include_dirs()]) return config + if __name__ == '__main__': from numpy.distutils.core import setup setup(maintainer='scikit-image Developers',