秩過濾器

秩過濾器是使用局部灰階排序來計算過濾值的非線性濾波器。這組濾波器有一個共同的基礎：局部灰階直方圖是在像素的鄰域 (由 2D 結構元素定義) 上計算的。如果將過濾值取為直方圖的中間值，則會得到經典的中值濾波器。

秩過濾器可用於多種目的，例如

• 影像品質增強，例如影像平滑、銳化

• 影像預處理，例如降噪、對比度增強

• 特徵提取，例如邊界偵測、隔離點偵測

• 影像後處理，例如小物件移除、物件分組、輪廓平滑

一些著名的濾波器 (例如，形態學膨脹和形態學侵蝕) 是秩過濾器的特定案例 [1]。

在此範例中，我們將看到如何使用 skimage 中提供的一些線性濾波器和非線性濾波器來過濾灰階影像。我們使用來自 skimage.data 的 camera 影像進行所有比較。

[1]

Pierre Soille, On morphological operators based on rank filters, Pattern Recognition 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.8/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()