使用填充修复斑点角膜图像

光学相干断层扫描 (OCT) 是一种非侵入性成像技术，眼科医生使用它来拍摄患者眼睛后部的图像 [1]。在执行 OCT 时，灰尘可能会粘在设备的参考镜上，导致图像上出现黑点。问题在于这些污点覆盖了体内组织的区域，从而隐藏了感兴趣的数据。我们的目标是根据其边界附近的像素恢复（重建）隐藏区域。

本教程改编自 Jules Scholler [2] 分享的应用程序。这些图像由 Viacheslav Mazlin 采集（请参阅 skimage.data.palisades_of_vogt()）。

[1]

David Turbert，由 Ninel Z Gregori, MD (2023) 审查 什么是光学相干断层扫描？，美国眼科学会。

[2]

Jules Scholler (2019) “使用填充进行图像去噪” https://www.jscholler.com/2019-02-28-remove-dots/

import matplotlib.pyplot as plt
import numpy as np
import plotly.io
import plotly.express as px

import skimage as ski

我们在这里使用的数据集是人体体内组织的图像序列（电影！）。具体来说，它显示了给定角膜样本的Vogt 嵴柱。

加载图像数据

image_seq = ski.data.palisades_of_vogt()

print(f'number of dimensions: {image_seq.ndim}')
print(f'shape: {image_seq.shape}')
print(f'dtype: {image_seq.dtype}')

number of dimensions: 3
shape: (60, 1440, 1440)
dtype: uint16

数据集是一个具有 60 帧（时间点）和 2 个空间维度的图像堆栈。让我们通过每 6 个时间点采样可视化 10 帧：我们可以看到照明的一些变化。我们利用 Plotly 的 imshow 函数中的 animation_frame 参数。作为旁注，当 binary_string 参数设置为 True 时，图像表示为灰度。

fig = px.imshow(
    image_seq[::6, :, :],
    animation_frame=0,
    binary_string=True,
    labels={'animation_frame': '6-step time point'},
    title='Sample of in-vivo human cornea',
)
fig.update_layout(autosize=False, minreducedwidth=250, minreducedheight=250)
plotly.io.show(fig)

按时间聚合

首先，我们想要检测丢失数据的那些污点。用技术术语来说，我们想要分割污点（对于序列中的所有帧）。与实际数据（信号）不同，污点不会在帧之间移动；它们是静止的。因此，我们首先计算图像序列的时间聚合。我们将使用中值图像来分割污点，后者相对于背景（模糊信号）脱颖而出。此外，为了了解（移动）数据，让我们计算方差。

image_med = np.median(image_seq, axis=0)
image_var = np.var(image_seq, axis=0)

assert image_var.shape == image_med.shape

print(f'shape: {image_med.shape}')

fig, ax = plt.subplots(ncols=2, figsize=(12, 6))

ax[0].imshow(image_med, cmap='gray')
ax[0].set_title('Image median over time')
ax[1].imshow(image_var, cmap='gray')
ax[1].set_title('Image variance over time')

fig.tight_layout()

shape: (1440, 1440)

使用局部阈值化

为了分割污点，我们使用阈值化。我们正在处理的图像照明不均匀，这会导致前景和背景（绝对）强度的空间变化。例如，一个区域的平均背景强度可能与另一个（遥远的）区域不同。因此，计算图像中不同的阈值值更合适，每个区域一个。这称为自适应（或局部）阈值化，而不是通常的阈值化过程，该过程对图像中的所有像素使用单个（全局）阈值。

当从 filters 模块调用 threshold_local 函数时，我们可以更改默认的邻域大小（block_size），即照明变化的典型大小（像素数），以及 offset（偏移邻域的加权平均值）。让我们尝试两个不同的 block_size 值。

thresh_1 = ski.filters.threshold_local(image_med, block_size=21, offset=15)
thresh_2 = ski.filters.threshold_local(image_med, block_size=43, offset=15)

mask_1 = image_med < thresh_1
mask_2 = image_med < thresh_2

让我们定义一个方便的函数来并排显示两个图，以便我们更容易地比较它们。

def plot_comparison(plot1, plot2, title1, title2):
    fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(12, 6), sharex=True, sharey=True)
    ax1.imshow(plot1, cmap='gray')
    ax1.set_title(title1)
    ax2.imshow(plot2, cmap='gray')
    ax2.set_title(title2)


plot_comparison(mask_1, mask_2, "block_size = 21", "block_size = 43")

在第二个掩模中，即使用较大的 block_size 值生成的掩模中，“污点”似乎更加明显。我们注意到，将偏移参数的值从其默认值 0 更改会产生更均匀的背景，让感兴趣的对象更明显地突出显示。请注意，切换参数值可以让我们更深入地了解所使用的方法，这通常可以使我们更接近所需的結果。

thresh_0 = ski.filters.threshold_local(image_med, block_size=43)

mask_0 = image_med < thresh_0

plot_comparison(mask_0, mask_2, "No offset", "Offset = 15")

去除细粒度特征

我们使用形态学滤波器来锐化掩模并专注于污点。两个基本的形态学算子是膨胀和腐蚀，其中膨胀（分别腐蚀）将像素设置为由结构元素（脚印）定义的邻域的最亮（分别最暗）值。

在这里，我们使用 morphology 模块中的 diamond 函数创建一个菱形脚印。腐蚀后跟膨胀称为开运算。首先，我们应用开运算滤波器，以去除小物体和细线，同时保留较大物体的形状和大小。

footprint = ski.morphology.diamond(3)
mask_open = ski.morphology.opening(mask_2, footprint)
plot_comparison(mask_2, mask_open, "mask before", "after opening")

由于“打开”图像始于腐蚀操作，因此小于结构元素的明亮区域已被删除。应用开运算滤波器时，调整脚印参数可以控制删除特征的细粒度。例如，如果我们在上面使用footprint = ski.morphology.diamond(1)，我们只会看到更小的特征被过滤掉，因此在掩模中保留更多斑点。相反，如果我们使用相同半径的圆盘形脚印，即footprint = ski.morphology.disk(3)，则会过滤掉更多细粒度的特征，因为圆盘的面积大于菱形的面积。

接下来，我们可以通过应用膨胀滤波器来使检测到的区域变宽。

mask_dilate = ski.morphology.dilation(mask_open, footprint)
plot_comparison(mask_open, mask_dilate, "Before", "After dilation")

膨胀会扩大明亮区域并缩小黑暗区域。请注意，白色斑点确实变厚了。

单独修复每个帧

我们现在可以对每个帧应用修复。为此，我们使用来自restoration模块的函数inpaint_biharmonic。它实现了基于双调和方程的算法。此函数接受两个数组作为输入：要恢复的图像和与我们要修复的区域相对应的掩模（具有相同的形状）。

image_seq_inpainted = np.zeros(image_seq.shape)

for i in range(image_seq.shape[0]):
    image_seq_inpainted[i] = ski.restoration.inpaint_biharmonic(
        image_seq[i], mask_dilate
    )

让我们可视化一个恢复的图像，其中污点已被修复。首先，我们找到污点的轮廓（好吧，掩模的轮廓），以便我们可以在恢复的图像上绘制它们。

contours = ski.measure.find_contours(mask_dilate)

# Gather all (row, column) coordinates of the contours
x = []
y = []
for contour in contours:
    x.append(contour[:, 0])
    y.append(contour[:, 1])
# Note that the following one-liner is equivalent to the above:
# x, y = zip(*((contour[:, 0], contour[:, 1]) for contour in contours))

# Flatten the coordinates
x_flat = np.concatenate(x).ravel().round().astype(int)
y_flat = np.concatenate(y).ravel().round().astype(int)
# Create mask of these contours
contour_mask = np.zeros(mask_dilate.shape, dtype=bool)
contour_mask[x_flat, y_flat] = 1
# Pick one frame
sample_result = image_seq_inpainted[12]
# Normalize it (so intensity values range [0, 1])
sample_result /= sample_result.max()

我们使用来自color模块的函数label2rgb，使用透明度（alpha 参数）将恢复的图像与分割的斑点叠加。

color_contours = ski.color.label2rgb(
    contour_mask, image=sample_result, alpha=0.4, bg_color=(1, 1, 1)
)

fig, ax = plt.subplots(figsize=(6, 6))

ax.imshow(color_contours)
ax.set_title('Segmented spots over restored image')

fig.tight_layout()

请注意，位于 (x, y) ~ (719, 1237) 的污点很突出；理想情况下，它应该被分割并修复。我们可以看到，我们“丢失”了它，因为在去除细粒度特征时，使用了开运算处理步骤。

plt.show()