Merge pull request #626 from ahojnnes/rank-fixes

Refactor rank filters.
This commit is contained in:
Stefan van der Walt
2013-08-01 07:42:57 -07:00
27 changed files with 2380 additions and 3398 deletions
+2
View File
@@ -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 <commit>``,
+14
View File
@@ -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``
+8 -23
View File
@@ -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
+244 -191
View File
@@ -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
+13 -25
View File
@@ -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()
+22 -15
View File
@@ -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()
+6 -6
View File
@@ -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.
+45 -12
View File
@@ -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',
-20
View File
@@ -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 *
-255
View File
@@ -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 = <dtype_t * >out.data
cdef dtype_t * image_data = <dtype_t * >image.data
cdef cnp.uint8_t * mask_data = <cnp.uint8_t * >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 = <Py_ssize_t * >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 = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
cdef Py_ssize_t * se_e_c = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
cdef Py_ssize_t * se_w_r = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
cdef Py_ssize_t * se_w_c = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
cdef Py_ssize_t * se_n_r = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
cdef Py_ssize_t * se_n_c = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
cdef Py_ssize_t * se_s_r = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
cdef Py_ssize_t * se_s_c = <Py_ssize_t * >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)
-25
View File
@@ -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 *
-422
View File
@@ -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 <dtype_t>(1. * (maxbin - 1) * (g - imin) / delta)
else:
return <dtype_t>(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 <dtype_t>(g - i)
else:
return <dtype_t>(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 <dtype_t>(((maxbin - 1) * sum) / pop)
else:
return <dtype_t>(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 <dtype_t>(imax - imin)
else:
return <dtype_t>(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 <dtype_t>(i)
return <dtype_t>(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 <dtype_t>(mean / pop)
else:
return <dtype_t>(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 <dtype_t>((g - mean / pop) / 2. + (midbin - 1))
else:
return <dtype_t>(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 <dtype_t>(i)
else:
return <dtype_t>(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 <dtype_t>(i)
else:
return <dtype_t>(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 <dtype_t>(imax)
else:
return <dtype_t>(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 <dtype_t>(imax)
else:
return <dtype_t>(imin)
else:
return <dtype_t>(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 <dtype_t>(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 <dtype_t>(g > (mean / pop))
else:
return <dtype_t>(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 <dtype_t>(i - g)
else:
return <dtype_t>(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 <dtype_t>e * 1000
else:
return <dtype_t>(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, <Py_ssize_t>0, <Py_ssize_t>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, <Py_ssize_t>0, <Py_ssize_t>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, <Py_ssize_t>0, <Py_ssize_t>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, <Py_ssize_t>0, <Py_ssize_t>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, <Py_ssize_t>0, <Py_ssize_t>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, <Py_ssize_t>0, <Py_ssize_t>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, <Py_ssize_t>0, <Py_ssize_t>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, <Py_ssize_t>0, <Py_ssize_t>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, <Py_ssize_t>0, <Py_ssize_t>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, <Py_ssize_t>0, <Py_ssize_t>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, <Py_ssize_t>0, <Py_ssize_t>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, <Py_ssize_t>0, <Py_ssize_t>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, <Py_ssize_t>0, <Py_ssize_t>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, <Py_ssize_t>0, <Py_ssize_t>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, <Py_ssize_t>0, <Py_ssize_t>0)
@@ -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 <dtype_t>(mean / bilat_pop)
else:
return <dtype_t>(0)
else:
return <dtype_t>(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 <dtype_t>(bilat_pop)
else:
return <dtype_t>(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)
@@ -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 <dtype_t>(1.0 * (maxbin - 1)
* (int_min(int_max(imin, g), imax)
- imin) / delta)
else:
return <dtype_t>(imax - imin)
else:
return <dtype_t>(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 <dtype_t>(imax - imin)
else:
return <dtype_t>(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 <dtype_t>(1.0 * mean / n)
else:
return <dtype_t>(0)
else:
return <dtype_t>(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 <dtype_t>((g - (mean / n)) * .5 + midbin)
else:
return <dtype_t>(0)
else:
return <dtype_t>(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 <dtype_t>imax
if g < imin:
return <dtype_t>imin
if imax - g < g - imin:
return <dtype_t>imax
else:
return <dtype_t>imin
else:
return <dtype_t>(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 <dtype_t>(i)
else:
return <dtype_t>(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 <dtype_t>(n)
else:
return <dtype_t>(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 <dtype_t>((maxbin - 1) * (g >= i))
else:
return <dtype_t>(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, <Py_ssize_t>0, <Py_ssize_t>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, <Py_ssize_t>0, <Py_ssize_t>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, <Py_ssize_t>0, <Py_ssize_t>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, <Py_ssize_t>0, <Py_ssize_t>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, <Py_ssize_t>0, <Py_ssize_t>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, <Py_ssize_t>0, <Py_ssize_t>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, <Py_ssize_t>0, <Py_ssize_t>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., <Py_ssize_t>0, <Py_ssize_t>0)
-483
View File
@@ -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 <dtype_t>(255. * (g - imin) / delta)
else:
return <dtype_t>(imax - imin)
else:
return <dtype_t>(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 <dtype_t>(g - i)
else:
return <dtype_t>(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 <dtype_t>((255 * sum) / pop)
else:
return <dtype_t>(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 <dtype_t>(imax - imin)
else:
return <dtype_t>(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 <dtype_t>(i)
else:
return <dtype_t>(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 <dtype_t>(mean / pop)
else:
return <dtype_t>(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 <dtype_t>((g - mean / pop) / 2. + 127)
else:
return <dtype_t>(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 <dtype_t>(i)
else:
return <dtype_t>(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 <dtype_t>(i)
else:
return <dtype_t>(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 <dtype_t>(imax)
else:
return <dtype_t>(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 <dtype_t>(imax)
else:
return <dtype_t>(imin)
else:
return <dtype_t>(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 <dtype_t>(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 <dtype_t>(g > (mean / pop))
else:
return <dtype_t>(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 <dtype_t>(i - g)
else:
return <dtype_t>(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 <dtype_t>0
for i in range(g, -1, -1):
if histo[i]:
break
min_i = g - i
for i in range(g, 256):
if histo[i]:
break
if i - g < min_i:
return <dtype_t>(i - g)
else:
return <dtype_t>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 <dtype_t>e * 10
else:
return <dtype_t>(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 <dtype_t>(0)
# maximizing the between class variance
max_i = 0
q1 = histo[0] / pop
m1 = 0.
max_sigma_b = 0.
for i in range(1, 256):
P = histo[i] / pop
new_q1 = q1 + P
if new_q1 > 0:
mu1 = (q1 * mu1 + i * P) / new_q1
mu2 = (mu - new_q1 * mu1) / (1. - new_q1)
sigma_b = new_q1 * (1. - new_q1) * (mu1 - mu2) ** 2
if sigma_b > max_sigma_b:
max_sigma_b = sigma_b
max_i = i
q1 = new_q1
return <dtype_t>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, <Py_ssize_t>0, <Py_ssize_t>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, <Py_ssize_t>0, <Py_ssize_t>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, <Py_ssize_t>0, <Py_ssize_t>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, <Py_ssize_t>0, <Py_ssize_t>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, <Py_ssize_t>0, <Py_ssize_t>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, <Py_ssize_t>0, <Py_ssize_t>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, <Py_ssize_t>0, <Py_ssize_t>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, <Py_ssize_t>0, <Py_ssize_t>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, <Py_ssize_t>0, <Py_ssize_t>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, <Py_ssize_t>0, <Py_ssize_t>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, <Py_ssize_t>0, <Py_ssize_t>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, <Py_ssize_t>0, <Py_ssize_t>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,
<Py_ssize_t>0, <Py_ssize_t>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, <Py_ssize_t>0, <Py_ssize_t>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, <Py_ssize_t>0, <Py_ssize_t>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, <Py_ssize_t>0, <Py_ssize_t>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, <Py_ssize_t>0, <Py_ssize_t>0)
-294
View File
@@ -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 <dtype_t>(255 * (uint8_min(uint8_max(imin, g), imax)
- imin) / delta)
else:
return <dtype_t>(imax - imin)
else:
return <dtype_t>(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 <dtype_t>(imax - imin)
else:
return <dtype_t>(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 <dtype_t>(1.0 * mean / n)
else:
return <dtype_t>(0)
else:
return <dtype_t>(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 <dtype_t>((g - (mean / n)) * .5 + 127)
else:
return <dtype_t>(0)
else:
return <dtype_t>(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 <dtype_t>imax
if g < imin:
return <dtype_t>imin
if imax - g < g - imin:
return <dtype_t>imax
else:
return <dtype_t>imin
else:
return <dtype_t>(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 <dtype_t>(i)
else:
return <dtype_t>(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 <dtype_t>(n)
else:
return <dtype_t>(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 <dtype_t>(255 * (g >= i))
else:
return <dtype_t>(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,
<Py_ssize_t>0, <Py_ssize_t>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,
<Py_ssize_t>0, <Py_ssize_t>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,
<Py_ssize_t>0, <Py_ssize_t>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, <Py_ssize_t>0, <Py_ssize_t>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, <Py_ssize_t>0, <Py_ssize_t>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., <Py_ssize_t>0, <Py_ssize_t>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,
<Py_ssize_t>0, <Py_ssize_t>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.,
<Py_ssize_t>0, <Py_ssize_t>0)
-773
View File
@@ -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)
@@ -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)
+70
View File
@@ -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)
+28
View File
@@ -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 *
@@ -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 = <Py_ssize_t>(selem.shape[0] / 2) + shift_y
cdef Py_ssize_t centre_c = <Py_ssize_t>(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 = <dtype_t * >out.data
cdef dtype_t * image_data = <dtype_t * >image.data
cdef dtype_t * mask_data = <dtype_t * >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 = <Py_ssize_t * >malloc(256 * sizeof(Py_ssize_t))
cdef Py_ssize_t* histo = <Py_ssize_t*>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 = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
cdef Py_ssize_t * se_e_c = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
cdef Py_ssize_t * se_w_r = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
cdef Py_ssize_t * se_w_c = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
cdef Py_ssize_t * se_n_r = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
cdef Py_ssize_t * se_n_c = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
cdef Py_ssize_t * se_s_r = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
cdef Py_ssize_t * se_s_c = <Py_ssize_t * >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 = <Py_ssize_t*>malloc(max_se * sizeof(Py_ssize_t))
cdef Py_ssize_t* se_e_c = <Py_ssize_t*>malloc(max_se * sizeof(Py_ssize_t))
cdef Py_ssize_t* se_w_r = <Py_ssize_t*>malloc(max_se * sizeof(Py_ssize_t))
cdef Py_ssize_t* se_w_c = <Py_ssize_t*>malloc(max_se * sizeof(Py_ssize_t))
cdef Py_ssize_t* se_n_r = <Py_ssize_t*>malloc(max_se * sizeof(Py_ssize_t))
cdef Py_ssize_t* se_n_c = <Py_ssize_t*>malloc(max_se * sizeof(Py_ssize_t))
cdef Py_ssize_t* se_s_r = <Py_ssize_t*>malloc(max_se * sizeof(Py_ssize_t))
cdef Py_ssize_t* se_s_c = <Py_ssize_t*>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] = <dtype_t_out>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] = <dtype_t_out>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] = <dtype_t_out>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] = <dtype_t_out>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] = <dtype_t_out>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)
+736 -7
View File
@@ -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)
+506
View File
@@ -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 <double>(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)
@@ -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)
+301
View File
@@ -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 <double>(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)
+109 -48
View File
@@ -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()
+9 -26
View File
@@ -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',