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

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

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

Scikit-image 实现了一个通用的滚动球算法版本，它不仅允许你使用球体，还可以使用任意形状作为核，并且适用于 n 维图像。这使你可以直接过滤 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 参数手动指定内核来完成的。

注意：半径等于椭圆的一个半轴的长度，即全轴的 一半 。因此，核形状乘以二。

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 的另一个特点是你可以直接将其应用于更高维度的图像，例如，通过共聚焦显微镜获得的图像堆栈。核的维度数量必须与图像的维度匹配，因此核的形状现在是三维的。

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的值，同时沿着所有三个维度进行过滤。

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-stack轴）过滤单个像素。

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