注意
转到末尾 下载完整的示例代码。或通过 Binder 在浏览器中运行此示例
Li 阈值分割#
1993 年，Li 和 Lee 提出了一种新的标准来找到“最佳”阈值，以区分图像的背景和前景[1]。他们提出，最小化前景与前景均值以及背景与背景均值之间的交叉熵，将在大多数情况下给出最佳阈值。
然而，直到 1998 年，找到此阈值的方法是尝试所有可能的阈值，然后选择交叉熵最小的阈值。在那时，Li 和 Tam 实现了一种新的迭代方法，通过使用交叉熵的斜率更快地找到最佳点[2]。这是 scikit-image 中实现的方法skimage.filters.threshold_li()。
在这里，我们演示了交叉熵及其 Li 迭代方法的优化。请注意，我们正在使用私有函数_cross_entropy，它不应在生产代码中使用，因为它可能会发生变化。
import numpy as np
import matplotlib.pyplot as plt
from skimage import data
from skimage import filters
from skimage.filters.thresholding import _cross_entropy
cell = data.cell()
camera = data.camera()
首先，让我们绘制skimage.data.camera()图像在所有可能阈值下的交叉熵。
thresholds = np.arange(np.min(camera) + 1.5, np.max(camera) - 1.5)
entropies = [_cross_entropy(camera, t) for t in thresholds]
optimal_camera_threshold = thresholds[np.argmin(entropies)]
fig, ax = plt.subplots(1, 3, figsize=(8, 3))
ax[0].imshow(camera, cmap='gray')
ax[0].set_title('image')
ax[0].set_axis_off()
ax[1].imshow(camera > optimal_camera_threshold, cmap='gray')
ax[1].set_title('thresholded')
ax[1].set_axis_off()
ax[2].plot(thresholds, entropies)
ax[2].set_xlabel('thresholds')
ax[2].set_ylabel('cross-entropy')
ax[2].vlines(
optimal_camera_threshold,
ymin=np.min(entropies) - 0.05 * np.ptp(entropies),
ymax=np.max(entropies) - 0.05 * np.ptp(entropies),
)
ax[2].set_title('optimal threshold')
fig.tight_layout()
print('The brute force optimal threshold is:', optimal_camera_threshold)
print('The computed optimal threshold is:', filters.threshold_li(camera))
plt.show()
The brute force optimal threshold is: 78.5
The computed optimal threshold is: 78.91288426606151
接下来，让我们使用iter_callback功能threshold_li来检查优化过程是如何发生的。
iter_thresholds = []
optimal_threshold = filters.threshold_li(camera, iter_callback=iter_thresholds.append)
iter_entropies = [_cross_entropy(camera, t) for t in iter_thresholds]
print('Only', len(iter_thresholds), 'thresholds examined.')
fig, ax = plt.subplots()
ax.plot(thresholds, entropies, label='all threshold entropies')
ax.plot(iter_thresholds, iter_entropies, label='optimization path')
ax.scatter(iter_thresholds, iter_entropies, c='C1')
ax.legend(loc='upper right')
plt.show()
Only 5 thresholds examined.
这显然比蛮力方法效率高得多。但是，在某些图像中，交叉熵不是凸的，这意味着只有一个最优点。在这种情况下，梯度下降可能会产生不是最优的阈值。在这个例子中，我们看到了优化过程的糟糕初始猜测如何导致糟糕的阈值选择。
iter_thresholds2 = []
opt_threshold2 = filters.threshold_li(
cell, initial_guess=64, iter_callback=iter_thresholds2.append
)
thresholds2 = np.arange(np.min(cell) + 1.5, np.max(cell) - 1.5)
entropies2 = [_cross_entropy(cell, t) for t in thresholds]
iter_entropies2 = [_cross_entropy(cell, t) for t in iter_thresholds2]
fig, ax = plt.subplots(1, 3, figsize=(8, 3))
ax[0].imshow(cell, cmap='magma')
ax[0].set_title('image')
ax[0].set_axis_off()
ax[1].imshow(cell > opt_threshold2, cmap='gray')
ax[1].set_title('thresholded')
ax[1].set_axis_off()
ax[2].plot(thresholds2, entropies2, label='all threshold entropies')
ax[2].plot(iter_thresholds2, iter_entropies2, label='optimization path')
ax[2].scatter(iter_thresholds2, iter_entropies2, c='C1')
ax[2].legend(loc='upper right')
plt.show()
在这张图片中，令人惊奇的是，默认的初始猜测，即图像平均值，实际上正好位于目标函数两个“谷”之间的峰值之上。不提供初始猜测，Li 的阈值分割方法根本不会执行任何操作！
iter_thresholds3 = []
opt_threshold3 = filters.threshold_li(cell, iter_callback=iter_thresholds3.append)
iter_entropies3 = [_cross_entropy(cell, t) for t in iter_thresholds3]
fig, ax = plt.subplots(1, 3, figsize=(8, 3))
ax[0].imshow(cell, cmap='magma')
ax[0].set_title('image')
ax[0].set_axis_off()
ax[1].imshow(cell > opt_threshold3, cmap='gray')
ax[1].set_title('thresholded')
ax[1].set_axis_off()
ax[2].plot(thresholds2, entropies2, label='all threshold entropies')
ax[2].plot(iter_thresholds3, iter_entropies3, label='optimization path')
ax[2].scatter(iter_thresholds3, iter_entropies3, c='C1')
ax[2].legend(loc='upper right')
plt.show()
为了了解正在发生的事情，让我们定义一个函数li_gradient，它复制 Li 方法的内部循环并返回从当前阈值到下一个阈值的变化。当此梯度为 0 时，我们处于所谓的驻点，Li 返回此值。当它为负时，下一个阈值猜测将更低，当它为正时，下一个猜测将更高。
在下图中，我们显示了初始猜测位于该熵峰值右侧时交叉熵和 Li 更新路径。我们叠加了阈值更新梯度，用threshold_li标记了 0 梯度线和默认初始猜测。
def li_gradient(image, t):
"""Find the threshold update at a given threshold."""
foreground = image > t
mean_fore = np.mean(image[foreground])
mean_back = np.mean(image[~foreground])
t_next = (mean_back - mean_fore) / (np.log(mean_back) - np.log(mean_fore))
dt = t_next - t
return dt
iter_thresholds4 = []
opt_threshold4 = filters.threshold_li(
cell, initial_guess=68, iter_callback=iter_thresholds4.append
)
iter_entropies4 = [_cross_entropy(cell, t) for t in iter_thresholds4]
print(len(iter_thresholds4), 'examined, optimum:', opt_threshold4)
gradients = [li_gradient(cell, t) for t in thresholds2]
fig, ax1 = plt.subplots()
ax1.plot(thresholds2, entropies2)
ax1.plot(iter_thresholds4, iter_entropies4)
ax1.scatter(iter_thresholds4, iter_entropies4, c='C1')
ax1.set_xlabel('threshold')
ax1.set_ylabel('cross entropy', color='C0')
ax1.tick_params(axis='y', labelcolor='C0')
ax2 = ax1.twinx()
ax2.plot(thresholds2, gradients, c='C3')
ax2.hlines(
[0], xmin=thresholds2[0], xmax=thresholds2[-1], colors='gray', linestyles='dashed'
)
ax2.vlines(
np.mean(cell),
ymin=np.min(gradients),
ymax=np.max(gradients),
colors='gray',
linestyles='dashed',
)
ax2.set_ylabel(r'$\Delta$(threshold)', color='C3')
ax2.tick_params(axis='y', labelcolor='C3')
fig.tight_layout()
plt.show()
8 examined, optimum: 111.68876119648344
除了允许用户提供数字作为初始猜测外，skimage.filters.threshold_li()可以接收一个函数，该函数根据图像强度进行猜测，就像numpy.mean()默认情况下所做的那样。当需要处理许多具有不同范围的图像时，这可能是一个不错的选择。
def quantile_95(image):
# you can use np.quantile(image, 0.95) if you have NumPy>=1.15
return np.percentile(image, 95)
iter_thresholds5 = []
opt_threshold5 = filters.threshold_li(
cell, initial_guess=quantile_95, iter_callback=iter_thresholds5.append
)
iter_entropies5 = [_cross_entropy(cell, t) for t in iter_thresholds5]
print(len(iter_thresholds5), 'examined, optimum:', opt_threshold5)
fig, ax1 = plt.subplots()
ax1.plot(thresholds2, entropies2)
ax1.plot(iter_thresholds5, iter_entropies5)
ax1.scatter(iter_thresholds5, iter_entropies5, c='C1')
ax1.set_xlabel('threshold')
ax1.set_ylabel('cross entropy', color='C0')
ax1.tick_params(axis='y', labelcolor='C0')
plt.show()
5 examined, optimum: 111.68876119648344
脚本的总运行时间：（0 分钟 4.998 秒）
