使用修复技术恢复斑点角膜图像

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

本教程改编自Jules Scholler [2] 分享的应用程序。图片由Viacheslav Mazlin获取（参见 skimage.data.palisades_of_vogt()）。

[1]

David Turbert，由Ninel Z Gregori，MD审阅（2023）`什么是光学相干断层扫描？ <https://www.aao.org/eye-health/treatments/what-is-optical-coherence-tomography>`_，美国眼科学会。

[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个空间维度的图像堆栈。让我们通过每隔六个时间点采样来可视化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 值生成的掩码中，“污点”似乎更加明显。我们注意到，将偏移参数的值从默认的零值增加，会得到一个更均匀的背景，使感兴趣的物体更加显眼。请注意，切换参数值可以让我们更深入地理解所使用的方法，这通常会使我们更接近预期的结果。

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