等级过滤器

等级过滤器是使用局部灰度级排序来计算滤波值的非线性过滤器。这组滤镜共享一个共同的基础：局部灰度级直方图是在像素的邻域(由2D结构元素定义)上计算的。如果将滤波值作为直方图的中间值，则得到经典的中值滤波。

等级过滤器可用于多种目的，例如：

• 图像质量增强，例如图像平滑、锐化

• 图像预处理，如降噪、对比度增强

• 特征提取，例如边界检测、孤立点检测

• 图像后处理，如小目标去除、目标分组、轮廓平滑

一些众所周知的过滤器(例如，形态膨胀和形态侵蚀)是等级过滤器的具体情况 1.

在本例中，我们将了解如何使用skImage中提供的一些线性和非线性滤镜来过滤灰度图像。我们使用 camera 图片来自 skimage.data 用于所有的比较。

1

Pierre Soille，关于基于等级过滤器的形态算子，模式识别35(2002)527-535， DOI:10.1016/S0031-3203(01)00047-4

import numpy as np
import matplotlib.pyplot as plt

from skimage.util import img_as_ubyte
from skimage import data
from skimage.exposure import histogram

noisy_image = img_as_ubyte(data.camera())
hist, hist_centers = histogram(noisy_image)

fig, ax = plt.subplots(ncols=2, figsize=(10, 5))

ax[0].imshow(noisy_image, cmap=plt.cm.gray)
ax[0].axis('off')

ax[1].plot(hist_centers, hist, lw=2)
ax[1].set_title('Gray-level histogram')

plt.tight_layout()

降噪

图像中添加了一些噪声：1%的像素随机设置为255，1%的像素随机设置为0。这个 中位数 采用滤波法去除噪声。

from skimage.filters.rank import median
from skimage.morphology import disk, ball

noise = np.random.random(noisy_image.shape)
noisy_image = img_as_ubyte(data.camera())
noisy_image[noise > 0.99] = 255
noisy_image[noise < 0.01] = 0

fig, axes = plt.subplots(2, 2, figsize=(10, 10), sharex=True, sharey=True)
ax = axes.ravel()

ax[0].imshow(noisy_image, vmin=0, vmax=255, cmap=plt.cm.gray)
ax[0].set_title('Noisy image')

ax[1].imshow(median(noisy_image, disk(1)), vmin=0, vmax=255, cmap=plt.cm.gray)
ax[1].set_title('Median $r=1$')

ax[2].imshow(median(noisy_image, disk(5)), vmin=0, vmax=255, cmap=plt.cm.gray)
ax[2].set_title('Median $r=5$')

ax[3].imshow(median(noisy_image, disk(20)), vmin=0, vmax=255, cmap=plt.cm.gray)
ax[3].set_title('Median $r=20$')

for a in ax:
    a.axis('off')

plt.tight_layout()

由于图像的缺省值很小(1-像素宽)，因此可以有效地去除添加的噪声，较小的过滤器半径就足够了。随着半径的增加，较大尺寸的对象也会被过滤，例如摄影机三脚架。中值滤波通常用于去噪，因为它保留了边界。例如，考虑仅位于整个图像的几个像素上的噪声，就像椒盐噪声的情况一样 2: 中值过滤器将忽略噪声像素，因为它们将显示为异常值；因此，与移动平均过滤器所做的相反，它不会显著改变一组局部像素的中值。

2

https://en.wikipedia.org/wiki/Salt-and-pepper_noise

图像平滑

下面的示例显示了一个本地 mean 滤镜使摄像师图像变得平滑。

from skimage.filters.rank import mean

loc_mean = mean(noisy_image, disk(10))

fig, ax = plt.subplots(ncols=2, figsize=(10, 5), sharex=True, sharey=True)

ax[0].imshow(noisy_image, vmin=0, vmax=255, cmap=plt.cm.gray)
ax[0].set_title('Original')

ax[1].imshow(loc_mean, vmin=0, vmax=255, cmap=plt.cm.gray)
ax[1].set_title('Local mean $r=10$')

for a in ax:
    a.axis('off')

plt.tight_layout()

人们可能对在保持重要边界的同时平滑图像感兴趣(中值过滤器已经做到了这一点)。在这里，我们使用 双边 滤镜将本地邻域限制为灰度级与中心灰度级相似的像素。

备注

中的彩色图像有不同的实现方式 skimage.restoration.denoise_bilateral() 。

from skimage.filters.rank import mean_bilateral

noisy_image = img_as_ubyte(data.camera())

bilat = mean_bilateral(noisy_image.astype(np.uint16), disk(20), s0=10, s1=10)

fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(10, 10),
                         sharex='row', sharey='row')
ax = axes.ravel()

ax[0].imshow(noisy_image, cmap=plt.cm.gray)
ax[0].set_title('Original')

ax[1].imshow(bilat, cmap=plt.cm.gray)
ax[1].set_title('Bilateral mean')

ax[2].imshow(noisy_image[100:250, 350:450], cmap=plt.cm.gray)

ax[3].imshow(bilat[100:250, 350:450], cmap=plt.cm.gray)

for a in ax:
    a.axis('off')

plt.tight_layout()

人们可以看到，图像的大部分连续部分(例如天空)被平滑，而其他细节被保留。

对比度增强

我们在这里比较全局直方图均衡化是如何局部应用的。

均衡后的图像 3 对于每个像素邻域具有大致线性的累积分布函数。本地版本 4 直方图均衡化强调每个局部灰度的变化。

3

https://en.wikipedia.org/wiki/Histogram_equalization

4

https://en.wikipedia.org/wiki/Adaptive_histogram_equalization

from skimage import exposure
from skimage.filters import rank

noisy_image = img_as_ubyte(data.camera())

# equalize globally and locally
glob = exposure.equalize_hist(noisy_image) * 255
loc = rank.equalize(noisy_image, disk(20))

# extract histogram for each image
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))

fig, axes = plt.subplots(nrows=3, ncols=2, figsize=(12, 12))
ax = axes.ravel()

ax[0].imshow(noisy_image, cmap=plt.cm.gray)
ax[0].axis('off')

ax[1].plot(hist[1][:-1], hist[0], lw=2)
ax[1].set_title('Histogram of gray values')

ax[2].imshow(glob, cmap=plt.cm.gray)
ax[2].axis('off')

ax[3].plot(glob_hist[1][:-1], glob_hist[0], lw=2)
ax[3].set_title('Histogram of gray values')

ax[4].imshow(loc, cmap=plt.cm.gray)
ax[4].axis('off')

ax[5].plot(loc_hist[1][:-1], loc_hist[0], lw=2)
ax[5].set_title('Histogram of gray values')

plt.tight_layout()

最大化用于图像的灰度级数量的另一种方法是应用局部自动平整，即像素的灰度值在局部最小值和局部最大值之间按比例重新映射。

以下示例显示了本地自动调平如何增强Camara MAN画面。

from skimage.filters.rank import autolevel

noisy_image = img_as_ubyte(data.camera())

auto = autolevel(noisy_image.astype(np.uint16), disk(20))

fig, ax = plt.subplots(ncols=2, figsize=(10, 5), sharex=True, sharey=True)

ax[0].imshow(noisy_image, cmap=plt.cm.gray)
ax[0].set_title('Original')

ax[1].imshow(auto, cmap=plt.cm.gray)
ax[1].set_title('Local autolevel')

for a in ax:
    a.axis('off')

plt.tight_layout()

该滤波器对局部异常值非常敏感。人们可以使用自动级别过滤器的百分位数版本来调节这一点，该过滤器使用给定的百分位数(一个较低的，一个较高的)来代替局部最小值和最大值。下面的示例说明了百分位数参数如何影响本地自动级别调整结果。

from skimage.filters.rank import autolevel_percentile

image = data.camera()

selem = disk(20)
loc_autolevel = autolevel(image, selem=selem)
loc_perc_autolevel0 = autolevel_percentile(image, selem=selem, p0=.01, p1=.99)
loc_perc_autolevel1 = autolevel_percentile(image, selem=selem, p0=.05, p1=.95)
loc_perc_autolevel2 = autolevel_percentile(image, selem=selem, p0=.1, p1=.9)
loc_perc_autolevel3 = autolevel_percentile(image, selem=selem, p0=.15, p1=.85)

fig, axes = plt.subplots(nrows=3, ncols=2, figsize=(10, 10),
                         sharex=True, sharey=True)
ax = axes.ravel()

title_list = ['Original',
              'auto_level',
              'auto-level 1%',
              'auto-level 5%',
              'auto-level 10%',
              'auto-level 15%']
image_list = [image,
              loc_autolevel,
              loc_perc_autolevel0,
              loc_perc_autolevel1,
              loc_perc_autolevel2,
              loc_perc_autolevel3]

for i in range(0, len(image_list)):
    ax[i].imshow(image_list[i], cmap=plt.cm.gray, vmin=0, vmax=255)
    ax[i].set_title(title_list[i])
    ax[i].axis('off')

plt.tight_layout()

输出：

/scikit-image/doc/examples/applications/plot_rank_filters.py:260: FutureWarning:

`selem` is a deprecated argument name for `autolevel`. It will be removed in version 1.0. Please use `footprint` instead.

/scikit-image/doc/examples/applications/plot_rank_filters.py:261: FutureWarning:

`selem` is a deprecated argument name for `autolevel_percentile`. It will be removed in version 1.0. Please use `footprint` instead.

/scikit-image/doc/examples/applications/plot_rank_filters.py:262: FutureWarning:

`selem` is a deprecated argument name for `autolevel_percentile`. It will be removed in version 1.0. Please use `footprint` instead.

/scikit-image/doc/examples/applications/plot_rank_filters.py:263: FutureWarning:

`selem` is a deprecated argument name for `autolevel_percentile`. It will be removed in version 1.0. Please use `footprint` instead.

/scikit-image/doc/examples/applications/plot_rank_filters.py:264: FutureWarning:

`selem` is a deprecated argument name for `autolevel_percentile`. It will be removed in version 1.0. Please use `footprint` instead.

如果原始像素值最接近局部最大值，则形态对比度增强过滤器用局部最大值替换中心像素，否则用最小局部最小值替换中心像素。

from skimage.filters.rank import enhance_contrast

noisy_image = img_as_ubyte(data.camera())

enh = enhance_contrast(noisy_image, disk(5))

fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(10, 10),
                         sharex='row', sharey='row')
ax = axes.ravel()

ax[0].imshow(noisy_image, cmap=plt.cm.gray)
ax[0].set_title('Original')

ax[1].imshow(enh, cmap=plt.cm.gray)
ax[1].set_title('Local morphological contrast enhancement')

ax[2].imshow(noisy_image[100:250, 350:450], cmap=plt.cm.gray)

ax[3].imshow(enh[100:250, 350:450], cmap=plt.cm.gray)

for a in ax:
    a.axis('off')

plt.tight_layout()

局部形态对比度增强的百分位数使用百分位数 p0 和 p1 而不是局部最小值和最大值。

from skimage.filters.rank import enhance_contrast_percentile

noisy_image = img_as_ubyte(data.camera())

penh = enhance_contrast_percentile(noisy_image, disk(5), p0=.1, p1=.9)

fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(10, 10),
                         sharex='row', sharey='row')
ax = axes.ravel()

ax[0].imshow(noisy_image, cmap=plt.cm.gray)
ax[0].set_title('Original')

ax[1].imshow(penh, cmap=plt.cm.gray)
ax[1].set_title('Local percentile morphological\n contrast enhancement')

ax[2].imshow(noisy_image[100:250, 350:450], cmap=plt.cm.gray)

ax[3].imshow(penh[100:250, 350:450], cmap=plt.cm.gray)

for a in ax:
    a.axis('off')

plt.tight_layout()

图像阈值

大津门限法 5 可以使用局部灰度级分布局部应用。在下面的例子中，对于每个像素，通过使由结构元素定义的局部邻域的两类像素之间的方差最大化来确定“最佳”阈值。

这些算法既可以用于2D图像，也可以用于3D图像。

该示例将局部阈值处理与全局阈值处理进行了比较，后者由 skimage.filters.threshold_otsu() 。请注意，前者比后者慢得多。

5

https://en.wikipedia.org/wiki/Otsu's_method

from skimage.filters.rank import otsu
from skimage.filters import threshold_otsu
from skimage import exposure

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

fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(12, 12),
                         sharex=True, sharey=True)
ax = axes.ravel()

fig.colorbar(ax[0].imshow(p8, cmap=plt.cm.gray), ax=ax[0])
ax[0].set_title('Original')

fig.colorbar(ax[1].imshow(t_loc_otsu, cmap=plt.cm.gray), ax=ax[1])
ax[1].set_title('Local Otsu ($r=%d$)' % radius)

ax[2].imshow(p8 >= t_loc_otsu, cmap=plt.cm.gray)
ax[2].set_title('Original >= local Otsu' % t_glob_otsu)

ax[3].imshow(glob_otsu, cmap=plt.cm.gray)
ax[3].set_title('Global Otsu ($t=%d$)' % t_glob_otsu)

for a in ax:
    a.axis('off')

plt.tight_layout()

下面的示例执行相同的比较，这次使用的是3D图像。

brain = exposure.rescale_intensity(data.brain().astype(float))

radius = 5
neighborhood = ball(radius)

# t_loc_otsu is an image
t_loc_otsu = rank.otsu(brain, neighborhood)
loc_otsu = brain >= t_loc_otsu

# t_glob_otsu is a scalar
t_glob_otsu = threshold_otsu(brain)
glob_otsu = brain >= t_glob_otsu

fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(12, 12),
                         sharex=True, sharey=True)
ax = axes.ravel()

slice_index = 3

fig.colorbar(ax[0].imshow(brain[slice_index], cmap=plt.cm.gray), ax=ax[0])
ax[0].set_title('Original')

fig.colorbar(ax[1].imshow(t_loc_otsu[slice_index], cmap=plt.cm.gray), ax=ax[1])
ax[1].set_title('Local Otsu ($r=%d$)' % radius)

ax[2].imshow(brain[slice_index] >= t_loc_otsu[slice_index], cmap=plt.cm.gray)
ax[2].set_title('Original >= local Otsu' % t_glob_otsu)

ax[3].imshow(glob_otsu[slice_index], cmap=plt.cm.gray)
ax[3].set_title('Global Otsu ($t=%d$)' % t_glob_otsu)

for a in ax:
    a.axis('off')

fig.tight_layout()

输出：

/vpy/lib/python3.9/site-packages/scikit_image-0.20.0.dev0-py3.9-linux-x86_64.egg/skimage/filters/rank/generic.py:282: UserWarning:

Possible precision loss converting image of type float64 to uint8 as required by rank filters. Convert manually using skimage.util.img_as_ubyte to silence this warning.

以下示例显示了局部OTSU阈值处理如何处理应用于合成图像的全局级别移动。

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))

fig, ax = plt.subplots(ncols=2, figsize=(10, 5),
                       sharex=True, sharey=True)

ax[0].imshow(m, cmap=plt.cm.gray)
ax[0].set_title('Original')

ax[1].imshow(m >= t, cmap=plt.cm.gray)
ax[1].set_title('Local Otsu ($r=%d$)' % radius)

for a in ax:
    a.axis('off')

plt.tight_layout()

图像形态

局部极大值和局部极小值是灰度形态的基本算子。

这里有一个经典的形态灰度滤镜的例子：开启、闭合和形态梯度。

from skimage.filters.rank import maximum, minimum, gradient

noisy_image = img_as_ubyte(data.camera())

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, axes = plt.subplots(nrows=2, ncols=2, figsize=(10, 10),
                         sharex=True, sharey=True)
ax = axes.ravel()

ax[0].imshow(noisy_image, cmap=plt.cm.gray)
ax[0].set_title('Original')

ax[1].imshow(closing, cmap=plt.cm.gray)
ax[1].set_title('Gray-level closing')

ax[2].imshow(opening, cmap=plt.cm.gray)
ax[2].set_title('Gray-level opening')

ax[3].imshow(grad, cmap=plt.cm.gray)
ax[3].set_title('Morphological gradient')

for a in ax:
    a.axis('off')

plt.tight_layout()

特征提取

局部直方图可以用来计算与局部图像复杂度有关的局部熵。使用以2为底的对数计算熵，即，过滤器返回编码局部灰度级分布所需的最小比特数。

skimage.filters.rank.entropy() 返回给定结构元素的局部熵。下面的示例将此滤镜应用于8位和16位图像。

备注

为了更好地使用可用的图像位，该函数为8位图像返回10倍的熵，为16位图像返回1000倍的熵。

from skimage import data
from skimage.filters.rank import entropy
from skimage.morphology import disk
import numpy as np
import matplotlib.pyplot as plt

image = data.camera()

fig, ax = plt.subplots(ncols=2, figsize=(12, 6), sharex=True, sharey=True)

fig.colorbar(ax[0].imshow(image, cmap=plt.cm.gray), ax=ax[0])
ax[0].set_title('Image')

fig.colorbar(ax[1].imshow(entropy(image, disk(5)), cmap=plt.cm.gray), ax=ax[1])
ax[1].set_title('Entropy')

for a in ax:
    a.axis('off')

plt.tight_layout()

实施

城市的中心部分 skimage.filters.rank 滤镜建立在更新局部灰度级直方图的滑动窗口上。该方法将算法复杂度限制在O(N)，其中n是图像像素数。关于结构化元素的大小，复杂性也是有限的。

在下面的内容中，我们将比较 skimage 。

from time import time

from scipy.ndimage import percentile_filter
from skimage.morphology import dilation
from skimage.filters.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 ndi_med(image, n):
    return percentile_filter(image, 50, size=n * 2 - 1)

比较两种情况

• skimage.filters.rank.maximum

• skimage.morphology.dilation

关于增加结构元素大小：

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)

fig, ax = plt.subplots(figsize=(10, 10), sharey=True)
ax.set_title('Performance with respect to element size')
ax.set_ylabel('Time (ms)')
ax.set_xlabel('Element radius')
ax.plot(e_range, rec)
ax.legend(['filters.rank.maximum', 'morphology.dilate'])

plt.tight_layout()

输出：

/scikit-image/doc/examples/applications/plot_rank_filters.py:588: FutureWarning:

`selem` is a deprecated argument name for `maximum`. It will be removed in version 1.0. Please use `footprint` instead.

/scikit-image/doc/examples/applications/plot_rank_filters.py:593: FutureWarning:

`selem` is a deprecated argument name for `dilation`. It will be removed in version 1.0. Please use `footprint` instead.

以及增加图像大小：

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(np.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)

fig, ax = plt.subplots()
ax.set_title('Performance with respect to image size')
ax.set_ylabel('Time (ms)')
ax.set_xlabel('Image size')
ax.plot(s_range, rec)
ax.legend(['filters.rank.maximum', 'morphology.dilate'])

plt.tight_layout()

输出：

/scikit-image/doc/examples/applications/plot_rank_filters.py:588: FutureWarning:

`selem` is a deprecated argument name for `maximum`. It will be removed in version 1.0. Please use `footprint` instead.

/scikit-image/doc/examples/applications/plot_rank_filters.py:593: FutureWarning:

`selem` is a deprecated argument name for `dilation`. It will be removed in version 1.0. Please use `footprint` instead.

比较：

• skimage.filters.rank.median

• scipy.ndimage.percentile_filter

关于增加结构元素大小：

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)
    rndi, ms_ndi = ndi_med(a, r)
    rec.append((ms_rc, ms_ndi))

rec = np.asarray(rec)

fig, ax = plt.subplots()
ax.set_title('Performance with respect to element size')
ax.plot(e_range, rec)
ax.legend(['filters.rank.median', 'scipy.ndimage.percentile'])
ax.set_ylabel('Time (ms)')
ax.set_xlabel('Element radius')

输出：

/scikit-image/doc/examples/applications/plot_rank_filters.py:583: FutureWarning:

`selem` is a deprecated argument name for `median`. It will be removed in version 1.0. Please use `footprint` instead.


Text(0.5, 23.52222222222222, 'Element radius')

两种方法结果的比较：

fig, ax = plt.subplots(ncols=2, figsize=(10, 5), sharex=True, sharey=True)

ax[0].set_title('filters.rank.median')
ax[0].imshow(rc, cmap=plt.cm.gray)

ax[1].set_title('scipy.ndimage.percentile')
ax[1].imshow(rndi, cmap=plt.cm.gray)

for a in ax:
    a.axis('off')

plt.tight_layout()

关于增加图像大小：

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(np.uint8)
    (rc, ms_rc) = cr_med(a, elem)
    rndi, ms_ndi = ndi_med(a, r)
    rec.append((ms_rc, ms_ndi))

rec = np.asarray(rec)

fig, ax = plt.subplots()
ax.set_title('Performance with respect to image size')
ax.plot(s_range, rec)
ax.legend(['filters.rank.median', 'scipy.ndimage.percentile'])
ax.set_ylabel('Time (ms)')
ax.set_xlabel('Image size')

plt.tight_layout()

plt.show()

输出：

/scikit-image/doc/examples/applications/plot_rank_filters.py:583: FutureWarning:

`selem` is a deprecated argument name for `median`. It will be removed in version 1.0. Please use `footprint` instead.