mirror of
https://github.com/wassname/scikit-image.git
synced 2026-07-18 12:40:14 +08:00
Improve long rank filter example
This commit is contained in:
Johannes Schönberger
@@ -3,11 +3,11 @@
|
||||
Rank filters
|
||||
============
|
||||
|
||||
Rank filters are non-linear filters using the local greylevels ordering to
|
||||
Rank filters are non-linear filters using the local gray-level ordering to
|
||||
compute the filtered value. This ensemble of filters share a common base: the
|
||||
local grey-level histogram extraction computed on the neighborhood of a pixel
|
||||
(defined by a 2D structuring element). If the filtered value is taken as the
|
||||
middle value of the histogram, we get the classical median filter.
|
||||
local gray-level histogram is computed on the neighborhood of a pixel (defined
|
||||
by a 2-D structuring element). If the filtered value is taken as the middle
|
||||
value of the histogram, we get the classical median filter.
|
||||
|
||||
Rank filters can be used for several purposes such as:
|
||||
|
||||
@@ -26,11 +26,9 @@ Rank filters can be used for several purposes such as:
|
||||
Some well known filters are specific cases of rank filters [1]_ e.g.
|
||||
morphological dilation, morphological erosion, median filters.
|
||||
|
||||
The different implementation availables in `skimage` are compared.
|
||||
|
||||
In this example, we will see how to filter a greylevel image using some of the
|
||||
linear and non-linear filters availables in skimage. We use the `camera`
|
||||
image from `skimage.data`.
|
||||
In this example, we will see how to filter a gray-level image using some of the
|
||||
linear and non-linear filters available in skimage. We use the `camera` image
|
||||
from `skimage.data` for all comparisons.
|
||||
|
||||
.. [1] Pierre Soille, On morphological operators based on rank filters, Pattern
|
||||
Recognition 35 (2002) 527-535.
|
||||
@@ -42,16 +40,16 @@ import matplotlib.pyplot as plt
|
||||
|
||||
from skimage import data
|
||||
|
||||
ima = data.camera()
|
||||
hist = np.histogram(ima, bins=np.arange(0, 256))
|
||||
noisy_image = data.camera()
|
||||
hist = np.histogram(noisy_image, bins=np.arange(0, 256))
|
||||
|
||||
plt.figure(figsize=(8, 3))
|
||||
plt.subplot(1, 2, 1)
|
||||
plt.imshow(ima, cmap=plt.cm.gray, interpolation='nearest')
|
||||
plt.imshow(noisy_image, interpolation='nearest')
|
||||
plt.axis('off')
|
||||
plt.subplot(1, 2, 2)
|
||||
plt.plot(hist[1][:-1], hist[0], lw=2)
|
||||
plt.title('histogram of grey values')
|
||||
plt.title('Histogram of grey values')
|
||||
|
||||
"""
|
||||
|
||||
@@ -65,50 +63,56 @@ randomly set to 0. The **median** filter is applied to remove the noise.
|
||||
|
||||
.. note::
|
||||
|
||||
there are different implementations of median filter :
|
||||
There are different implementations of median filter:
|
||||
`skimage.filter.median_filter` and `skimage.filter.rank.median`
|
||||
|
||||
"""
|
||||
|
||||
noise = np.random.random(ima.shape)
|
||||
nima = data.camera()
|
||||
nima[noise > 0.99] = 255
|
||||
nima[noise < 0.01] = 0
|
||||
|
||||
from skimage.filter.rank import median
|
||||
from skimage.morphology import disk
|
||||
|
||||
fig = plt.figure(figsize=[10, 7])
|
||||
noise = np.random.random(noisy_image.shape)
|
||||
noisy_image = data.camera()
|
||||
noisy_image[noise > 0.99] = 255
|
||||
noisy_image[noise < 0.01] = 0
|
||||
|
||||
fig = plt.figure(figsize=(10, 7))
|
||||
|
||||
lo = median(nima, disk(1))
|
||||
hi = median(nima, disk(5))
|
||||
ext = median(nima, disk(20))
|
||||
plt.subplot(2, 2, 1)
|
||||
plt.imshow(nima, cmap=plt.cm.gray, vmin=0, vmax=255)
|
||||
plt.xlabel('noised image')
|
||||
plt.imshow(noisy_image, vmin=0, vmax=255)
|
||||
plt.title('Noisy image')
|
||||
plt.axis('off')
|
||||
|
||||
plt.subplot(2, 2, 2)
|
||||
plt.imshow(lo, cmap=plt.cm.gray, vmin=0, vmax=255)
|
||||
plt.xlabel('median $r=1$')
|
||||
plt.imshow(median(noisy_image, disk(1)), vmin=0, vmax=255)
|
||||
plt.title('Median $r=1$')
|
||||
plt.axis('off')
|
||||
|
||||
plt.subplot(2, 2, 3)
|
||||
plt.imshow(hi, cmap=plt.cm.gray, vmin=0, vmax=255)
|
||||
plt.xlabel('median $r=5$')
|
||||
plt.imshow(median(noisy_image, disk(5)), vmin=0, vmax=255)
|
||||
plt.title('Median $r=5$')
|
||||
plt.axis('off')
|
||||
|
||||
plt.subplot(2, 2, 4)
|
||||
plt.imshow(ext, cmap=plt.cm.gray, vmin=0, vmax=255)
|
||||
plt.xlabel('median $r=20$')
|
||||
plt.imshow(median(noisy_image, disk(20)), vmin=0, vmax=255)
|
||||
plt.title('Median $r=20$')
|
||||
plt.axis('off')
|
||||
|
||||
"""
|
||||
|
||||
.. image:: PLOT2RST.current_figure
|
||||
|
||||
The added noise is efficiently removed, as the image defaults are small (1 pixel
|
||||
wide), a small filter radius is sufficient. As the radius is increasing, objects
|
||||
with a bigger size are filtered as well, such as the camera tripod. The median
|
||||
filter is commonly used for noise removal because borders are preserved.
|
||||
The added noise is efficiently removed, as the image defaults are small (1
|
||||
pixel wide), a small filter radius is sufficient. As the radius is increasing,
|
||||
objects with bigger sizes are filtered as well, such as the camera tripod. The
|
||||
median filter is often used for noise removal because borders are preserved and
|
||||
e.g. salt and pepper noise typically does not distort the gray-level.
|
||||
|
||||
Image smoothing
|
||||
================
|
||||
|
||||
The example hereunder shows how a local **mean** smoothes the camera man image.
|
||||
The example hereunder shows how a local **mean** filter smooths the camera man
|
||||
image.
|
||||
|
||||
"""
|
||||
|
||||
@@ -116,13 +120,17 @@ from skimage.filter.rank import mean
|
||||
|
||||
fig = plt.figure(figsize=[10, 7])
|
||||
|
||||
loc_mean = mean(nima, disk(10))
|
||||
loc_mean = mean(noisy_image, disk(10))
|
||||
|
||||
plt.subplot(1, 2, 1)
|
||||
plt.imshow(ima, cmap=plt.cm.gray, vmin=0, vmax=255)
|
||||
plt.xlabel('original')
|
||||
plt.imshow(noisy_image, vmin=0, vmax=255)
|
||||
plt.title('Original')
|
||||
plt.axis('off')
|
||||
|
||||
plt.subplot(1, 2, 2)
|
||||
plt.imshow(loc_mean, cmap=plt.cm.gray, vmin=0, vmax=255)
|
||||
plt.xlabel('local mean $r=10$')
|
||||
plt.imshow(loc_mean, vmin=0, vmax=255)
|
||||
plt.title('Local mean $r=10$')
|
||||
plt.axis('off')
|
||||
|
||||
"""
|
||||
|
||||
@@ -130,35 +138,42 @@ plt.xlabel('local mean $r=10$')
|
||||
|
||||
One may be interested in smoothing an image while preserving important borders
|
||||
(median filters already achieved this), here we use the **bilateral** filter
|
||||
that restricts the local neighborhood to pixel having a greylevel similar to
|
||||
that restricts the local neighborhood to pixel having a gray-level similar to
|
||||
the central one.
|
||||
|
||||
.. note::
|
||||
|
||||
a different implementation is available for color images in
|
||||
A different implementation is available for color images in
|
||||
`skimage.filter.denoise_bilateral`.
|
||||
|
||||
"""
|
||||
|
||||
from skimage.filter.rank import bilateral_mean
|
||||
|
||||
ima = data.camera()
|
||||
noisy_image = data.camera()
|
||||
selem = disk(10)
|
||||
|
||||
bilat = bilateral_mean(ima.astype(np.uint16), disk(20), s0=10, s1=10)
|
||||
bilat = bilateral_mean(noisy_image.astype(np.uint16), disk(20), s0=10, s1=10)
|
||||
|
||||
# display results
|
||||
fig = plt.figure(figsize=[10, 7])
|
||||
|
||||
plt.subplot(2, 2, 1)
|
||||
plt.imshow(ima, cmap=plt.cm.gray)
|
||||
plt.xlabel('original')
|
||||
plt.imshow(noisy_image, cmap=plt.cm.gray)
|
||||
plt.title('Original')
|
||||
plt.axis('off')
|
||||
|
||||
plt.subplot(2, 2, 3)
|
||||
plt.imshow(bilat, cmap=plt.cm.gray)
|
||||
plt.xlabel('bilateral mean')
|
||||
plt.title('Bilateral mean')
|
||||
plt.axis('off')
|
||||
|
||||
plt.subplot(2, 2, 2)
|
||||
plt.imshow(ima[200:350, 350:450], cmap=plt.cm.gray)
|
||||
plt.imshow(noisy_image[200:350, 350:450], cmap=plt.cm.gray)
|
||||
plt.axis('off')
|
||||
|
||||
plt.subplot(2, 2, 4)
|
||||
plt.imshow(bilat[200:350, 350:450], cmap=plt.cm.gray)
|
||||
plt.axis('off')
|
||||
|
||||
"""
|
||||
|
||||
@@ -175,7 +190,7 @@ We compare here how the global histogram equalization is applied locally.
|
||||
|
||||
The equalized image [2]_ has a roughly linear cumulative distribution function
|
||||
for each pixel neighborhood. The local version [3]_ of the histogram
|
||||
equalization emphasizes every local greylevel variations.
|
||||
equalization emphasizes every local gray-level variations.
|
||||
|
||||
.. [2] http://en.wikipedia.org/wiki/Histogram_equalization
|
||||
.. [3] http://en.wikipedia.org/wiki/Adaptive_histogram_equalization
|
||||
@@ -185,74 +200,86 @@ equalization emphasizes every local greylevel variations.
|
||||
from skimage import exposure
|
||||
from skimage.filter import rank
|
||||
|
||||
ima = data.camera()
|
||||
noisy_image = data.camera()
|
||||
|
||||
# equalize globally and locally
|
||||
glob = exposure.equalize(ima) * 255
|
||||
loc = rank.equalize(ima, disk(20))
|
||||
glob = exposure.equalize(noisy_image) * 255
|
||||
loc = rank.equalize(noisy_image, disk(20))
|
||||
|
||||
# extract histogram for each image
|
||||
hist = np.histogram(ima, bins=np.arange(0, 256))
|
||||
hist = np.histogram(noisy_image, bins=np.arange(0, 256))
|
||||
glob_hist = np.histogram(glob, bins=np.arange(0, 256))
|
||||
loc_hist = np.histogram(loc, bins=np.arange(0, 256))
|
||||
|
||||
plt.figure(figsize=(10, 10))
|
||||
|
||||
plt.subplot(321)
|
||||
plt.imshow(ima, cmap=plt.cm.gray, interpolation='nearest')
|
||||
plt.imshow(noisy_image, interpolation='nearest')
|
||||
plt.axis('off')
|
||||
|
||||
plt.subplot(322)
|
||||
plt.plot(hist[1][:-1], hist[0], lw=2)
|
||||
plt.title('histogram of grey values')
|
||||
plt.title('Histogram of gray values')
|
||||
|
||||
plt.subplot(323)
|
||||
plt.imshow(glob, cmap=plt.cm.gray, interpolation='nearest')
|
||||
plt.imshow(glob, interpolation='nearest')
|
||||
plt.axis('off')
|
||||
|
||||
plt.subplot(324)
|
||||
plt.plot(glob_hist[1][:-1], glob_hist[0], lw=2)
|
||||
plt.title('histogram of grey values')
|
||||
plt.title('Histogram of gray values')
|
||||
|
||||
plt.subplot(325)
|
||||
plt.imshow(loc, cmap=plt.cm.gray, interpolation='nearest')
|
||||
plt.imshow(loc, interpolation='nearest')
|
||||
plt.axis('off')
|
||||
|
||||
plt.subplot(326)
|
||||
plt.plot(loc_hist[1][:-1], loc_hist[0], lw=2)
|
||||
plt.title('histogram of grey values')
|
||||
plt.title('Histogram of gray values')
|
||||
|
||||
"""
|
||||
|
||||
.. image:: PLOT2RST.current_figure
|
||||
|
||||
another way to maximize the number of greylevels used for an image is to apply
|
||||
a local autoleveling, i.e. here a pixel greylevel is proportionally remapped
|
||||
between local minimum and local maximum.
|
||||
Another way to maximize the number of gray-levels used for an image is to apply
|
||||
a local auto-leveling, i.e. the gray-value of a pixel is proportionally
|
||||
remapped between local minimum and local maximum.
|
||||
|
||||
The following example shows how local autolevel enhances the camara man picture.
|
||||
The following example shows how local auto-level enhances the camara man
|
||||
picture.
|
||||
|
||||
"""
|
||||
|
||||
from skimage.filter.rank import autolevel
|
||||
|
||||
ima = data.camera()
|
||||
noisy_image = data.camera()
|
||||
selem = disk(10)
|
||||
|
||||
auto = autolevel(ima.astype(np.uint16), disk(20))
|
||||
auto = autolevel(noisy_image.astype(np.uint16), disk(20))
|
||||
|
||||
# display results
|
||||
fig = plt.figure(figsize=[10, 7])
|
||||
|
||||
plt.subplot(1, 2, 1)
|
||||
plt.imshow(ima, cmap=plt.cm.gray)
|
||||
plt.xlabel('original')
|
||||
plt.imshow(noisy_image, cmap=plt.cm.gray)
|
||||
plt.title('Original')
|
||||
plt.axis('off')
|
||||
|
||||
plt.subplot(1, 2, 2)
|
||||
plt.imshow(auto, cmap=plt.cm.gray)
|
||||
plt.xlabel('local autolevel')
|
||||
plt.title('Local autolevel')
|
||||
plt.axis('off')
|
||||
|
||||
"""
|
||||
|
||||
.. image:: PLOT2RST.current_figure
|
||||
|
||||
This filter is very sensitive to local outlayers, see the little white spot in
|
||||
the sky left part. This is due to a local maximum which is very high comparing
|
||||
to the rest of the neighborhood. One can moderate this using the percentile
|
||||
version of the autolevel filter which uses given percentiles (one inferior,
|
||||
one superior) in place of local minimum and maximum. The example below
|
||||
illustrates how the percentile parameters influence the local autolevel result.
|
||||
This filter is very sensitive to local outliers, see the little white spot in
|
||||
the left part of the sky. This is due to a local maximum which is very high
|
||||
comparing to the rest of the neighborhood. One can moderate this using the
|
||||
percentile version of the auto-level filter which uses given percentiles (one
|
||||
inferior, one superior) in place of local minimum and maximum. The example
|
||||
below illustrates how the percentile parameters influence the local auto-level
|
||||
result.
|
||||
|
||||
"""
|
||||
|
||||
@@ -272,14 +299,14 @@ ax0, ax1, ax2 = axes
|
||||
plt.gray()
|
||||
|
||||
ax0.imshow(np.hstack((image, loc_autolevel)))
|
||||
ax0.set_title('original / autolevel')
|
||||
ax0.set_title('Original / auto-level')
|
||||
|
||||
ax1.imshow(
|
||||
np.hstack((loc_perc_autolevel0, loc_perc_autolevel1)), vmin=0, vmax=255)
|
||||
ax1.set_title('percentile autolevel 0%,1%')
|
||||
ax1.set_title('Percentile auto-level 0%,1%')
|
||||
ax2.imshow(
|
||||
np.hstack((loc_perc_autolevel2, loc_perc_autolevel3)), vmin=0, vmax=255)
|
||||
ax2.set_title('percentile autolevel 5% and 10%')
|
||||
ax2.set_title('Percentile auto-level 5% and 10%')
|
||||
|
||||
for ax in axes:
|
||||
ax.axis('off')
|
||||
@@ -289,29 +316,35 @@ for ax in axes:
|
||||
.. image:: PLOT2RST.current_figure
|
||||
|
||||
The morphological contrast enhancement filter replaces the central pixel by the
|
||||
local maximum if the original pixel value is closest to local maximum, otherwise
|
||||
by the minimum local.
|
||||
local maximum if the original pixel value is closest to local maximum,
|
||||
otherwise by the minimum local.
|
||||
|
||||
"""
|
||||
|
||||
from skimage.filter.rank import morph_contr_enh
|
||||
|
||||
ima = data.camera()
|
||||
noisy_image = data.camera()
|
||||
|
||||
enh = morph_contr_enh(ima, disk(5))
|
||||
enh = morph_contr_enh(noisy_image, disk(5))
|
||||
|
||||
# display results
|
||||
fig = plt.figure(figsize=[10, 7])
|
||||
plt.subplot(2, 2, 1)
|
||||
plt.imshow(ima, cmap=plt.cm.gray)
|
||||
plt.xlabel('original')
|
||||
plt.imshow(noisy_image, cmap=plt.cm.gray)
|
||||
plt.title('Original')
|
||||
plt.axis('off')
|
||||
|
||||
plt.subplot(2, 2, 3)
|
||||
plt.imshow(enh, cmap=plt.cm.gray)
|
||||
plt.xlabel('local morphlogical contrast enhancement')
|
||||
plt.title('Local morphological contrast enhancement')
|
||||
plt.axis('off')
|
||||
|
||||
plt.subplot(2, 2, 2)
|
||||
plt.imshow(ima[200:350, 350:450], cmap=plt.cm.gray)
|
||||
plt.imshow(noisy_image[200:350, 350:450], cmap=plt.cm.gray)
|
||||
plt.axis('off')
|
||||
|
||||
plt.subplot(2, 2, 4)
|
||||
plt.imshow(enh[200:350, 350:450], cmap=plt.cm.gray)
|
||||
plt.axis('off')
|
||||
|
||||
"""
|
||||
|
||||
@@ -324,22 +357,28 @@ percentile *p0* and *p1* instead of the local minimum and maximum.
|
||||
|
||||
from skimage.filter.rank import percentile_morph_contr_enh
|
||||
|
||||
ima = data.camera()
|
||||
noisy_image = data.camera()
|
||||
|
||||
penh = percentile_morph_contr_enh(ima, disk(5), p0=.1, p1=.9)
|
||||
penh = percentile_morph_contr_enh(noisy_image, disk(5), p0=.1, p1=.9)
|
||||
|
||||
# display results
|
||||
fig = plt.figure(figsize=[10, 7])
|
||||
plt.subplot(2, 2, 1)
|
||||
plt.imshow(ima, cmap=plt.cm.gray)
|
||||
plt.xlabel('original')
|
||||
plt.imshow(noisy_image, cmap=plt.cm.gray)
|
||||
plt.title('Original')
|
||||
plt.axis('off')
|
||||
|
||||
plt.subplot(2, 2, 3)
|
||||
plt.imshow(penh, cmap=plt.cm.gray)
|
||||
plt.xlabel('local percentile morphlogical\n contrast enhancement')
|
||||
plt.title('Local percentile morphological\n contrast enhancement')
|
||||
plt.axis('off')
|
||||
|
||||
plt.subplot(2, 2, 2)
|
||||
plt.imshow(ima[200:350, 350:450], cmap=plt.cm.gray)
|
||||
plt.imshow(noisy_image[200:350, 350:450], cmap=plt.cm.gray)
|
||||
plt.axis('off')
|
||||
|
||||
plt.subplot(2, 2, 4)
|
||||
plt.imshow(penh[200:350, 350:450], cmap=plt.cm.gray)
|
||||
plt.axis('off')
|
||||
|
||||
"""
|
||||
|
||||
@@ -348,18 +387,18 @@ plt.imshow(penh[200:350, 350:450], cmap=plt.cm.gray)
|
||||
Image threshold
|
||||
===============
|
||||
|
||||
The Otsu's threshold [1]_ method can be applied locally using the local
|
||||
greylevel distribution. In the example below, for each pixel, an "optimal"
|
||||
threshold is determined by maximizing the variance between two classes of pixels
|
||||
of the local neighborhood defined by a structuring element.
|
||||
The Otsu threshold [1]_ method can be applied locally using the local gray-
|
||||
level distribution. In the example below, for each pixel, an "optimal"
|
||||
threshold is determined by maximizing the variance between two classes of
|
||||
pixels of the local neighborhood defined by a structuring element.
|
||||
|
||||
The example compares the local threshold with the global threshold
|
||||
`skimage.filter.threshold_otsu`.
|
||||
|
||||
.. note::
|
||||
|
||||
Local thresholding is much slower than global one. There exists a function
|
||||
for global Otsu thresholding: `skimage.filter.threshold_otsu`.
|
||||
Local is much slower than global thresholding. A function for global Otsu
|
||||
thresholding can be found in : `skimage.filter.threshold_otsu`.
|
||||
|
||||
.. [4] http://en.wikipedia.org/wiki/Otsu's_method
|
||||
|
||||
@@ -382,27 +421,35 @@ t_glob_otsu = threshold_otsu(p8)
|
||||
glob_otsu = p8 >= t_glob_otsu
|
||||
|
||||
plt.figure()
|
||||
|
||||
plt.subplot(2, 2, 1)
|
||||
plt.imshow(p8, cmap=plt.cm.gray)
|
||||
plt.xlabel('original')
|
||||
plt.title('Original')
|
||||
plt.colorbar()
|
||||
plt.axis('off')
|
||||
|
||||
plt.subplot(2, 2, 2)
|
||||
plt.imshow(t_loc_otsu, cmap=plt.cm.gray)
|
||||
plt.xlabel('local Otsu ($radius=%d$)' % radius)
|
||||
plt.title('Local Otsu ($r=%d$)' % radius)
|
||||
plt.colorbar()
|
||||
plt.axis('off')
|
||||
|
||||
plt.subplot(2, 2, 3)
|
||||
plt.imshow(p8 >= t_loc_otsu, cmap=plt.cm.gray)
|
||||
plt.xlabel('original>=local Otsu' % t_glob_otsu)
|
||||
plt.title('Original >= local Otsu' % t_glob_otsu)
|
||||
plt.axis('off')
|
||||
|
||||
plt.subplot(2, 2, 4)
|
||||
plt.imshow(glob_otsu, cmap=plt.cm.gray)
|
||||
plt.xlabel('global Otsu ($t=%d$)' % t_glob_otsu)
|
||||
plt.title('Global Otsu ($t=%d$)' % t_glob_otsu)
|
||||
plt.axis('off')
|
||||
|
||||
"""
|
||||
|
||||
.. image:: PLOT2RST.current_figure
|
||||
|
||||
The following example shows how local Otsu's threshold handles a global level
|
||||
shift applied to a synthetic image .
|
||||
The following example shows how local Otsu thresholding handles a global level
|
||||
shift applied to a synthetic image.
|
||||
|
||||
"""
|
||||
|
||||
@@ -413,13 +460,18 @@ m = (np.tile(x, (n, 1)) * np.linspace(0.1, 1, n) * 128 + 128).astype(np.uint8)
|
||||
|
||||
radius = 10
|
||||
t = rank.otsu(m, disk(radius))
|
||||
|
||||
plt.figure()
|
||||
|
||||
plt.subplot(1, 2, 1)
|
||||
plt.imshow(m)
|
||||
plt.xlabel('original')
|
||||
plt.title('Original')
|
||||
plt.axis('off')
|
||||
|
||||
plt.subplot(1, 2, 2)
|
||||
plt.imshow(m >= t, interpolation='nearest')
|
||||
plt.xlabel('local Otsu ($radius=%d$)' % radius)
|
||||
plt.title('Local Otsu ($r=%d$)' % radius)
|
||||
plt.axis('off')
|
||||
|
||||
"""
|
||||
|
||||
@@ -428,7 +480,7 @@ plt.xlabel('local Otsu ($radius=%d$)' % radius)
|
||||
Image morphology
|
||||
================
|
||||
|
||||
Local maximum and local minimum are the base operators for greylevel
|
||||
Local maximum and local minimum are the base operators for gray-level
|
||||
morphology.
|
||||
|
||||
.. note::
|
||||
@@ -436,33 +488,41 @@ morphology.
|
||||
`skimage.dilate` and `skimage.erode` are equivalent filters (see below for
|
||||
comparison).
|
||||
|
||||
Here is an example of the classical morphological greylevel filters: opening,
|
||||
Here is an example of the classical morphological gray-level filters: opening,
|
||||
closing and morphological gradient.
|
||||
|
||||
"""
|
||||
|
||||
from skimage.filter.rank import maximum, minimum, gradient
|
||||
|
||||
ima = data.camera()
|
||||
noisy_image = data.camera()
|
||||
|
||||
closing = maximum(minimum(ima, disk(5)), disk(5))
|
||||
opening = minimum(maximum(ima, disk(5)), disk(5))
|
||||
grad = gradient(ima, disk(5))
|
||||
closing = maximum(minimum(noisy_image, disk(5)), disk(5))
|
||||
opening = minimum(maximum(noisy_image, disk(5)), disk(5))
|
||||
grad = gradient(noisy_image, disk(5))
|
||||
|
||||
# display results
|
||||
fig = plt.figure(figsize=[10, 7])
|
||||
|
||||
plt.subplot(2, 2, 1)
|
||||
plt.imshow(ima, cmap=plt.cm.gray)
|
||||
plt.xlabel('original')
|
||||
plt.imshow(noisy_image, cmap=plt.cm.gray)
|
||||
plt.title('Original')
|
||||
plt.axis('off')
|
||||
|
||||
plt.subplot(2, 2, 2)
|
||||
plt.imshow(closing, cmap=plt.cm.gray)
|
||||
plt.xlabel('greylevel closing')
|
||||
plt.title('Gray-level closing')
|
||||
plt.axis('off')
|
||||
|
||||
plt.subplot(2, 2, 3)
|
||||
plt.imshow(opening, cmap=plt.cm.gray)
|
||||
plt.xlabel('greylevel opening')
|
||||
plt.title('Gray-level opening')
|
||||
plt.axis('off')
|
||||
|
||||
plt.subplot(2, 2, 4)
|
||||
plt.imshow(grad, cmap=plt.cm.gray)
|
||||
plt.xlabel('morphological gradient')
|
||||
plt.title('Morphological gradient')
|
||||
plt.axis('off')
|
||||
|
||||
"""
|
||||
|
||||
@@ -471,13 +531,14 @@ plt.xlabel('morphological gradient')
|
||||
Feature extraction
|
||||
===================
|
||||
|
||||
Local histogram can be exploited to compute local entropy, which is related to
|
||||
Local histograms can be exploited to compute local entropy, which is related to
|
||||
the local image complexity. Entropy is computed using base 2 logarithm i.e. the
|
||||
filter returns the minimum number of bits needed to encode local greylevel
|
||||
filter returns the minimum number of bits needed to encode local gray-level
|
||||
distribution.
|
||||
|
||||
`skimage.rank.entropy` returns local entropy on a given structuring element.
|
||||
The following example shows this filter applied on 8- and 16- bit images.
|
||||
`skimage.rank.entropy` returns the local entropy on a given structuring
|
||||
element. The following example shows applies this filter on 8- and 16-bit
|
||||
images.
|
||||
|
||||
.. note::
|
||||
|
||||
@@ -492,47 +553,36 @@ from skimage.morphology import disk
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
# defining a 8- and a 16-bit test images
|
||||
a8 = data.camera()
|
||||
a16 = data.camera().astype(np.uint16) * 4
|
||||
image = data.camera()
|
||||
|
||||
ent8 = entropy(a8, disk(5)) # pixel value contain 10x the local entropy
|
||||
ent16 = entropy(a16, disk(5)) # pixel value contain 1000x the local entropy
|
||||
plt.figure(figsize=(10, 4))
|
||||
|
||||
# display results
|
||||
plt.figure(figsize=(10, 10))
|
||||
|
||||
plt.subplot(2, 2, 1)
|
||||
plt.imshow(a8, cmap=plt.cm.gray)
|
||||
plt.xlabel('8-bit image')
|
||||
plt.subplot(1, 2, 1)
|
||||
plt.imshow(image, cmap=plt.cm.gray)
|
||||
plt.title('Image')
|
||||
plt.colorbar()
|
||||
plt.axis('off')
|
||||
|
||||
plt.subplot(2, 2, 2)
|
||||
plt.imshow(ent8, cmap=plt.cm.jet)
|
||||
plt.xlabel('entropy*10')
|
||||
plt.colorbar()
|
||||
|
||||
plt.subplot(2, 2, 3)
|
||||
plt.imshow(a16, cmap=plt.cm.gray)
|
||||
plt.xlabel('16-bit image')
|
||||
plt.colorbar()
|
||||
|
||||
plt.subplot(2, 2, 4)
|
||||
plt.imshow(ent16, cmap=plt.cm.jet)
|
||||
plt.xlabel('entropy*1000')
|
||||
plt.subplot(1, 2, 2)
|
||||
plt.imshow(entropy(image, disk(5)), cmap=plt.cm.jet)
|
||||
plt.title('Entropy')
|
||||
plt.colorbar()
|
||||
plt.axis('off')
|
||||
|
||||
"""
|
||||
|
||||
.. image:: PLOT2RST.current_figure
|
||||
|
||||
Implementation
|
||||
================
|
||||
==============
|
||||
|
||||
The central part of the `skimage.rank` filters is build on a sliding window that
|
||||
update local greylevel histogram. This approach limits the algorithm complexity
|
||||
to O(n) where n is the number of image pixels. The complexity is also limited
|
||||
with respect to the structuring element size.
|
||||
The central part of the `skimage.rank` filters is build on a sliding window
|
||||
that updates the local gray-level histogram. This approach limits the algorithm
|
||||
complexity to O(n) where n is the number of image pixels. The complexity is
|
||||
also limited with respect to the structuring element size.
|
||||
|
||||
In the following we compare the performance of different implementations
|
||||
available in `skimage`.
|
||||
|
||||
"""
|
||||
|
||||
@@ -583,10 +633,10 @@ def ndi_med(image, n):
|
||||
|
||||
Comparison between
|
||||
|
||||
* `rank.maximum`
|
||||
* `cmorph.dilate`
|
||||
* `filter.rank.maximum`
|
||||
* `morphology.dilate`
|
||||
|
||||
on increasing structuring element size
|
||||
on increasing structuring element size:
|
||||
|
||||
"""
|
||||
|
||||
@@ -603,18 +653,18 @@ for r in e_range:
|
||||
rec = np.asarray(rec)
|
||||
|
||||
plt.figure()
|
||||
plt.title('increasing element size')
|
||||
plt.ylabel('time (ms)')
|
||||
plt.xlabel('element radius')
|
||||
plt.title('Performance with respect to element size')
|
||||
plt.ylabel('Time (ms)')
|
||||
plt.title('Element radius')
|
||||
plt.plot(e_range, rec)
|
||||
plt.legend(['crank.maximum', 'cmorph.dilate'])
|
||||
plt.legend(['filter.rank.maximum', 'morphology.dilate'])
|
||||
|
||||
"""
|
||||
|
||||
and increasing image size
|
||||
|
||||
.. image:: PLOT2RST.current_figure
|
||||
|
||||
and increasing image size:
|
||||
|
||||
"""
|
||||
|
||||
r = 9
|
||||
@@ -623,7 +673,7 @@ elem = disk(r + 1)
|
||||
rec = []
|
||||
s_range = range(100, 1000, 100)
|
||||
for s in s_range:
|
||||
a = (np.random.random((s, s)) * 256).astype('uint8')
|
||||
a = (np.random.random((s, s)) * 256).astype(np.uint8)
|
||||
(rc, ms_rc) = cr_max(a, elem)
|
||||
(rcm, ms_rcm) = cm_dil(a, elem)
|
||||
rec.append((ms_rc, ms_rcm))
|
||||
@@ -631,11 +681,11 @@ for s in s_range:
|
||||
rec = np.asarray(rec)
|
||||
|
||||
plt.figure()
|
||||
plt.title('increasing image size')
|
||||
plt.ylabel('time (ms)')
|
||||
plt.xlabel('image size')
|
||||
plt.title('Performance with respect to image size')
|
||||
plt.ylabel('Time (ms)')
|
||||
plt.title('Image size')
|
||||
plt.plot(s_range, rec)
|
||||
plt.legend(['crank.maximum', 'cmorph.dilate'])
|
||||
plt.legend(['filter.rank.maximum', 'morphology.dilate'])
|
||||
|
||||
|
||||
"""
|
||||
@@ -644,11 +694,11 @@ plt.legend(['crank.maximum', 'cmorph.dilate'])
|
||||
|
||||
Comparison between:
|
||||
|
||||
* `rank.median`
|
||||
* `ctmf.median_filter`
|
||||
* `ndimage.percentile`
|
||||
* `filter.rank.median`
|
||||
* `filter.median_filter`
|
||||
* `scipy.ndimage.percentile`
|
||||
|
||||
on increasing structuring element size
|
||||
on increasing structuring element size:
|
||||
|
||||
"""
|
||||
|
||||
@@ -666,27 +716,29 @@ for r in e_range:
|
||||
rec = np.asarray(rec)
|
||||
|
||||
plt.figure()
|
||||
plt.title('increasing element size')
|
||||
plt.title('Performance with respect to element size')
|
||||
plt.plot(e_range, rec)
|
||||
plt.legend(['rank.median', 'ctmf.median_filter', 'ndimage.percentile'])
|
||||
plt.ylabel('time (ms)')
|
||||
plt.xlabel('element radius')
|
||||
plt.legend(['filter.rank.median', 'filter.median_filter',
|
||||
'scipy.ndimage.percentile'])
|
||||
plt.ylabel('Time (ms)')
|
||||
plt.title('Element radius')
|
||||
|
||||
"""
|
||||
.. image:: PLOT2RST.current_figure
|
||||
|
||||
comparison of outcome of the three methods
|
||||
Comparison of outcome of the three methods:
|
||||
|
||||
"""
|
||||
|
||||
plt.figure()
|
||||
plt.imshow(np.hstack((rc, rctmf, rndi)))
|
||||
plt.xlabel('rank.median vs ctmf.median_filter vs ndimage.percentile')
|
||||
plt.title('filter.rank.median vs filtermedian_filter vs scipy.ndimage.percentile')
|
||||
plt.axis('off')
|
||||
|
||||
"""
|
||||
.. image:: PLOT2RST.current_figure
|
||||
|
||||
and increasing image size
|
||||
and increasing image size:
|
||||
|
||||
"""
|
||||
|
||||
@@ -696,7 +748,7 @@ elem = disk(r + 1)
|
||||
rec = []
|
||||
s_range = [100, 200, 500, 1000]
|
||||
for s in s_range:
|
||||
a = (np.random.random((s, s)) * 256).astype('uint8')
|
||||
a = (np.random.random((s, s)) * 256).astype(np.uint8)
|
||||
(rc, ms_rc) = cr_med(a, elem)
|
||||
rctmf, ms_rctmf = ctmf_med(a, r)
|
||||
rndi, ms_ndi = ndi_med(a, r)
|
||||
@@ -705,11 +757,12 @@ for s in s_range:
|
||||
rec = np.asarray(rec)
|
||||
|
||||
plt.figure()
|
||||
plt.title('increasing image size')
|
||||
plt.title('Performance with respect to image size')
|
||||
plt.plot(s_range, rec)
|
||||
plt.legend(['rank.median', 'ctmf.median_filter', 'ndimage.percentile'])
|
||||
plt.ylabel('time (ms)')
|
||||
plt.xlabel('image size')
|
||||
plt.legend(['filter.rank.median', 'filter.median_filter',
|
||||
'scipy.ndimage.percentile'])
|
||||
plt.ylabel('Time (ms)')
|
||||
plt.title('Image size')
|
||||
|
||||
"""
|
||||
.. image:: PLOT2RST.current_figure
|
||||
|
||||
@@ -20,7 +20,6 @@ image = img_as_ubyte(data.camera())
|
||||
|
||||
fig, (ax0, ax1) = plt.subplots(ncols=2, figsize=(10, 4))
|
||||
|
||||
|
||||
img0 = ax0.imshow(image, cmap=plt.cm.gray)
|
||||
ax0.set_title('Image')
|
||||
ax0.axis('off')
|
||||
|
||||
Reference in New Issue
Block a user