秩序滤波器

秩序滤波器是非线性滤波器，它使用局部灰度级排序来计算滤波后的值。 这些滤波器集合共享一个共同的基础：局部灰度级直方图是在像素的邻域（由 2D 结构元素定义）上计算的。 如果滤波后的值取自直方图的中间值，我们将得到经典的中值滤波器。

秩序滤波器可用于多种目的，例如

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

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

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

• 图像后处理，例如去除小物体、物体分组、轮廓平滑

一些众所周知的滤波器（例如形态学膨胀和形态学腐蚀）是秩序滤波器的特例 [1].

在本示例中，我们将看到如何使用 skimage 中提供的一些线性滤波器和非线性滤波器来过滤灰度图像。 我们使用 skimage.data 中的 camera 图像进行所有比较。

[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].set_axis_off()

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

fig.tight_layout()

噪声去除

在图像中添加了一些噪声：1% 的像素随机设置为 255，1% 随机设置为 0。 应用 **中值** 滤波器来去除噪声。

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

rng = np.random.default_rng()
noise = rng.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.set_axis_off()

fig.tight_layout()

添加的噪声被有效地去除，因为图像默认值很小（1 像素宽），较小的滤波器半径就足够了。 随着半径的增加，尺寸更大的物体也会被过滤掉，例如相机三脚架。 中值滤波器通常用于降噪，因为它可以保留边界。 例如，考虑仅位于整个图像中几个像素上的噪声，就像椒盐噪声的情况一样 [2]：中值滤波器会忽略噪声像素，因为它们将显示为异常值； 因此，它不会显着改变一组局部像素的中值，这与移动平均滤波器所做的事情不同。

[2]

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

图像平滑

以下示例显示了局部 **均值** 滤波器如何平滑摄像师图像。

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

fig.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.set_axis_off()

fig.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].set_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].set_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].set_axis_off()

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

fig.tight_layout()

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

以下示例显示了局部自动调平如何增强摄像师图片。

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

fig.tight_layout()

此滤波器对局部异常值非常敏感。 您可以使用自动调平滤波器的百分位数版本来缓解这种情况，该版本使用给定的百分位数（一个较低，一个较高）来代替局部最小值和最大值。 以下示例说明了百分位数参数如何影响局部自动调平的结果。

from skimage.filters.rank import autolevel_percentile

image = data.camera()

footprint = disk(20)
loc_autolevel = autolevel(image, footprint=footprint)
loc_perc_autolevel0 = autolevel_percentile(image, footprint=footprint, p0=0.01, p1=0.99)
loc_perc_autolevel1 = autolevel_percentile(image, footprint=footprint, p0=0.05, p1=0.95)
loc_perc_autolevel2 = autolevel_percentile(image, footprint=footprint, p0=0.1, p1=0.9)
loc_perc_autolevel3 = autolevel_percentile(image, footprint=footprint, p0=0.15, p1=0.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].set_axis_off()

fig.tight_layout()

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

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

fig.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=0.1, p1=0.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.set_axis_off()

fig.tight_layout()

图像阈值

Otsu 阈值方法 [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
footprint = disk(radius)

# t_loc_otsu is an image
t_loc_otsu = otsu(p8, footprint)
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(f'Local Otsu ($r={radius}$)')

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

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

for a in ax:
    a.set_axis_off()

fig.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(f'Local Otsu ($r={radius}$)')

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

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

for a in ax:
    a.set_axis_off()

fig.tight_layout()

/opt/hostedtoolcache/Python/3.12.3/x64/lib/python3.12/site-packages/skimage/filters/rank/generic.py:353: 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(f'Local Otsu ($r={radius}$)')

for a in ax:
    a.set_axis_off()

fig.tight_layout()

图像形态学

局部最大值和局部最小值是灰度级形态学的基础算子。

以下是一个经典的形态学灰度级滤波器的示例：开运算、闭运算和形态学梯度。

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

noisy_image = img_as_ubyte(data.camera())

opening = maximum(minimum(noisy_image, disk(5)), disk(5))
closing = 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.set_axis_off()

fig.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.set_axis_off()

fig.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, footprint):
    return median(image=image, footprint=footprint)


@exec_and_timeit
def cr_max(image, footprint):
    return maximum(image=image, footprint=footprint)


@exec_and_timeit
def cm_dil(image, footprint):
    return dilation(image=image, footprint=footprint)


@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'])

fig.tight_layout()

以及图像大小的增加

r = 9
elem = disk(r + 1)

rec = []
s_range = range(100, 1000, 100)
for s in s_range:
    a = (rng.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'])

fig.tight_layout()

比较

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

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

fig.tight_layout()

随着图像大小的增加

r = 9
elem = disk(r + 1)

rec = []
s_range = [100, 200, 500, 1000]
for s in s_range:
    a = (rng.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')

fig.tight_layout()

plt.show()