使用滚球算法估计背景强度

滚球算法在曝光不均匀的情况下估计灰度图像的背景强度。它常用于生物医学图像处理，并由 Stanley R. Sternberg 于 1983 年首次提出 [1]。

该算法作为一种滤波器工作，并且非常直观。我们将图像视为一个表面，在其每个像素位置上堆叠了大小为单位的块。块的数量，以及因此的表面高度，由像素的强度决定。为了获得所需（像素）位置的背景强度，我们想象在所需位置将一个球体浸入表面下方。一旦它完全被块体覆盖，球体的顶点就决定了该位置的背景强度。然后，我们可以滚动这个球体在表面下方，以获得整个图像的背景值。

Scikit-image 实现了一个此滚球算法的通用版本，它允许您不仅可以使用球体，还可以使用任意形状作为内核，并且适用于 n 维 ndimages。这允许您直接过滤 RGB 图像或沿任何（或所有）空间维度过滤图像堆栈。

[1]

Sternberg，Stanley R.“生物医学图像处理。”计算机 1（1983）：22-34。DOI:10.1109/MC.1983.1654163

经典滚球

在 scikit-image 中，滚球算法假设您的背景强度较低（黑色），而特征强度较高（白色）。如果您的图像情况如此，您可以直接使用如下所示的滤波器

import matplotlib.pyplot as plt
import numpy as np
import pywt

from skimage import data, restoration, util


def plot_result(image, background):
    fig, ax = plt.subplots(nrows=1, ncols=3)

    ax[0].imshow(image, cmap='gray')
    ax[0].set_title('Original image')
    ax[0].axis('off')

    ax[1].imshow(background, cmap='gray')
    ax[1].set_title('Background')
    ax[1].axis('off')

    ax[2].imshow(image - background, cmap='gray')
    ax[2].set_title('Result')
    ax[2].axis('off')

    fig.tight_layout()


image = data.coins()

background = restoration.rolling_ball(image)

plot_result(image, background)
plt.show()

白色背景

如果您在明亮背景上有深色特征，则需要在将其传递到算法之前反转图像，然后反转结果。这可以通过以下方式实现：

image = data.page()
image_inverted = util.invert(image)

background_inverted = restoration.rolling_ball(image_inverted, radius=45)
filtered_image_inverted = image_inverted - background_inverted
filtered_image = util.invert(filtered_image_inverted)
background = util.invert(background_inverted)

fig, ax = plt.subplots(nrows=1, ncols=3)

ax[0].imshow(image, cmap='gray')
ax[0].set_title('Original image')
ax[0].axis('off')

ax[1].imshow(background, cmap='gray')
ax[1].set_title('Background')
ax[1].axis('off')

ax[2].imshow(filtered_image, cmap='gray')
ax[2].set_title('Result')
ax[2].axis('off')

fig.tight_layout()

plt.show()

小心不要在减去明亮背景时出现整数下溢。例如，此代码看起来是正确的，但可能会遇到下溢导致不希望出现的伪影。您可以在可视化的右上角看到这一点。

image = data.page()
image_inverted = util.invert(image)

background_inverted = restoration.rolling_ball(image_inverted, radius=45)
background = util.invert(background_inverted)
underflow_image = image - background  # integer underflow occurs here

# correct subtraction
correct_image = util.invert(image_inverted - background_inverted)

fig, ax = plt.subplots(nrows=1, ncols=2)

ax[0].imshow(underflow_image, cmap='gray')
ax[0].set_title('Background Removal with Underflow')
ax[0].axis('off')

ax[1].imshow(correct_image, cmap='gray')
ax[1].set_title('Correct Background Removal')
ax[1].axis('off')

fig.tight_layout()

plt.show()

图像数据类型

rolling_ball 可以处理除 np.uint8 之外的其他数据类型。您可以以相同的方式将它们传递到函数中。

image = data.coins()[:200, :200].astype(np.uint16)

background = restoration.rolling_ball(image, radius=70.5)
plot_result(image, background)
plt.show()

但是，如果您使用归一化到 [0, 1] 的浮点图像，则需要注意。在这种情况下，球体将比图像强度大得多，这可能导致意外的结果。

image = util.img_as_float(data.coins()[:200, :200])

background = restoration.rolling_ball(image, radius=70.5)
plot_result(image, background)
plt.show()

因为 radius=70.5 比图像的最大强度大得多，所以有效内核大小会大大减小，即，仅滚动图像中球体的一小部分（大约 radius=10）。您可以在下面的“高级形状”部分中找到此奇怪效果的再现。

为了获得预期的结果，您需要降低内核的强度。这可以通过使用 kernel 参数手动指定内核来完成。

注意：半径等于椭圆的半轴长度，即完整轴的一半。因此，内核形状乘以 2。

normalized_radius = 70.5 / 255
image = util.img_as_float(data.coins())
kernel = restoration.ellipsoid_kernel((70.5 * 2, 70.5 * 2), normalized_radius * 2)

background = restoration.rolling_ball(image, kernel=kernel)
plot_result(image, background)
plt.show()

高级形状

默认情况下，rolling_ball 使用球形内核（不出所料）。有时，这可能过于限制——如上例所示——因为强度维度与空间维度具有不同的比例，或者因为图像维度可能具有不同的含义——一个维度可能是图像堆栈中的堆栈计数器。

为了解决这个问题，rolling_ball 有一个 kernel 参数，允许您指定要使用的内核。内核必须与图像具有相同的维度（注意：维度，而不是形状）。为了帮助创建它，skimage 提供了两个默认内核。ball_kernel 指定球形内核，并用作默认内核。ellipsoid_kernel 指定椭球形内核。

image = data.coins()
kernel = restoration.ellipsoid_kernel((70.5 * 2, 70.5 * 2), 70.5 * 2)

background = restoration.rolling_ball(image, kernel=kernel)
plot_result(image, background)
plt.show()

您还可以使用 ellipsoid_kernel 重新创建先前的意外结果，并看到有效（空间）滤波器大小已减小。

image = data.coins()

kernel = restoration.ellipsoid_kernel((10 * 2, 10 * 2), 255 * 2)

background = restoration.rolling_ball(image, kernel=kernel)
plot_result(image, background)
plt.show()

更高维度

rolling_ball 的另一个特性是您可以直接将其应用于更高维度的图像，例如，在共聚焦显微镜期间获得的图像 z 堆栈。内核维度数必须与图像维度数匹配，因此内核形状现在是 3 维的。

image = data.cells3d()[:, 1, ...]
background = restoration.rolling_ball(
    image, kernel=restoration.ellipsoid_kernel((1, 21, 21), 0.1)
)

plot_result(image[30, ...], background[30, ...])
plt.show()

内核大小为 1 不会沿此轴进行滤波。换句话说，上述滤波器将分别应用于堆栈中的每个图像。

但是，您还可以通过指定除 1 之外的值来同时沿所有 3 个维度进行滤波。

image = data.cells3d()[:, 1, ...]
background = restoration.rolling_ball(
    image, kernel=restoration.ellipsoid_kernel((5, 21, 21), 0.1)
)

plot_result(image[30, ...], background[30, ...])
plt.show()

另一种可能性是沿平面轴（z 堆栈轴）过滤单个像素。

image = data.cells3d()[:, 1, ...]
background = restoration.rolling_ball(
    image, kernel=restoration.ellipsoid_kernel((100, 1, 1), 0.1)
)

plot_result(image[30, ...], background[30, ...])
plt.show()

1D 信号滤波

作为 rolling_ball 的 n 维特征的另一个示例，我们展示了一个 1D 数据的实现。在这里，我们感兴趣的是去除 ECG 波形的背景信号以检测突出峰值（高于局部基线的值）。使用较小的半径值可以去除更平滑的峰值。

x = pywt.data.ecg()
background = restoration.rolling_ball(x, radius=80)
background2 = restoration.rolling_ball(x, radius=10)
plt.figure()
plt.plot(x, label='original')
plt.plot(x - background, label='radius=80')
plt.plot(x - background2, label='radius=10')
plt.legend()
plt.show()