Merge pull request #348 from odebeir/rank_filters

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