注意
转到末尾下载完整的示例代码。或通过 Binder 在浏览器中运行此示例
跟踪金属合金的凝固过程#
在此示例中，我们识别并跟踪正在凝固的镍基合金中的固液 (S-L) 界面。跟踪凝固过程中的变化可以计算凝固速度。这对于表征样品的凝固结构非常重要，并将用于为金属增材制造研究提供信息。图像序列由先进非铁结构合金中心 (CANFSA) 在阿贡国家实验室 (ANL) 的高级光子源 (APS) 使用同步加速器 X 射线照相术获得。此分析首次在会议上提出[1]。
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import plotly.express as px
import plotly.io
from skimage import filters, measure, restoration
from skimage.data import nickel_solidification
image_sequence = nickel_solidification()
y0, y1, x0, x1 = 0, 180, 100, 330
image_sequence = image_sequence[:, y0:y1, x0:x1]
print(f'shape: {image_sequence.shape}')
shape: (11, 180, 230)
该数据集是一个包含 11 帧（时间点）的 2D 图像堆栈。我们通过一个工作流程来可视化和分析它，其中第一个图像处理步骤在整个三维数据集（即跨空间和时间）上执行，以便优先去除局部瞬时噪声，而不是物理特征（例如，气泡、飞溅等），这些特征在前后帧中大致位于相同的位置。
fig = px.imshow(
image_sequence,
animation_frame=0,
binary_string=True,
labels={'animation_frame': 'time point'},
)
plotly.io.show(fig)
计算图像增量#
让我们首先应用高斯低通滤波器，以平滑图像并减少噪声。接下来，我们计算图像增量，即两个连续帧之间差异的序列。为此，我们将以倒数第二帧结束的图像序列减去以第二帧开始的图像序列。
smoothed = filters.gaussian(image_sequence)
image_deltas = smoothed[1:, :, :] - smoothed[:-1, :, :]
fig = px.imshow(
image_deltas,
animation_frame=0,
binary_string=True,
labels={'animation_frame': 'time point'},
)
plotly.io.show(fig)
裁剪最低和最高强度#
现在，我们计算 image_deltas 的第 5 个和第 95 个百分位强度：我们希望裁剪低于第 5 个百分位强度和高于第 95 个百分位强度的强度，同时还将强度值重新缩放到 [0, 1]。
p_low, p_high = np.percentile(image_deltas, [5, 95])
clipped = image_deltas - p_low
clipped[clipped < 0.0] = 0.0
clipped = clipped / p_high
clipped[clipped > 1.0] = 1.0
fig = px.imshow(
clipped,
animation_frame=0,
binary_string=True,
labels={'animation_frame': 'time point'},
)
plotly.io.show(fig)
反转和去噪#
我们反转 clipped 图像，以便最高强度区域将包括我们感兴趣跟踪的区域（即，S-L 界面）。然后，我们应用全变分去噪滤波器以减少界面之外的噪声。
inverted = 1 - clipped
denoised = restoration.denoise_tv_chambolle(inverted, weight=0.15)
fig = px.imshow(
denoised,
animation_frame=0,
binary_string=True,
labels={'animation_frame': 'time point'},
)
plotly.io.show(fig)
二值化#
我们的下一步是创建二值图像，将图像分为前景和背景：我们希望 S-L 界面成为每个二值图像前景中最突出的特征，以便最终将其与图像的其余部分分开。
我们需要一个阈值 thresh_val 来创建我们的二值图像 binarized。可以手动设置此值，但我们将使用 scikit-image 的 filters 子模块中的自动最小阈值方法（还有其他方法可能更适合不同的应用）。
thresh_val = filters.threshold_minimum(denoised)
binarized = denoised > thresh_val
fig = px.imshow(
binarized,
animation_frame=0,
binary_string=True,
labels={'animation_frame': 'time point'},
)
plotly.io.show(fig)
选择最大区域#
在我们的二值图像中，S-L 界面显示为连接像素的最大区域。对于工作流程的此步骤，我们将分别在每个 2D 图像上进行操作，而不是在整个 3D 数据集上进行操作，因为我们对每个区域的单个时间点感兴趣。
我们在二值图像上应用 skimage.measure.label()，以便每个区域都有自己的标签。然后，我们通过计算区域属性（包括 area 属性）并按 area 值进行排序来选择每个图像中最大的区域。函数 skimage.measure.regionprops_table() 返回一个区域属性表，该表可以读入 Pandas DataFrame。首先，让我们考虑此工作流程阶段的第一个图像增量 binarized[0, :, :]。
labeled_0 = measure.label(binarized[0, :, :])
props_0 = measure.regionprops_table(labeled_0, properties=('label', 'area', 'bbox'))
props_0_df = pd.DataFrame(props_0)
props_0_df = props_0_df.sort_values('area', ascending=False)
# Show top five rows
props_0_df.head()
因此，我们可以通过将其标签与上述（排序）表的第一行中的标签匹配来选择最大的区域。让我们可视化它及其红色边界框 (bbox)。
我们可以看到，通过将相同的边界框叠加到第 0 个原始图像上，该框的下限如何与 S-L 界面的底部对齐。此边界框是根据第 0 个和第 1 个图像之间的图像增量计算得出的，但该框最底部的区域对应于界面较早时间（第 0 个图像）的位置，因为界面正在向上移动。
fig = px.imshow(image_sequence[0, :, :], binary_string=True)
fig.add_shape(type='rect', x0=minc, y0=minr, x1=maxc, y1=maxr, line=dict(color='Red'))
plotly.io.show(fig)
现在，我们已准备好对序列中的所有图像增量执行此选择。我们还将存储 bbox 信息，这将用于跟踪 S-L 界面的位置。
largest_region = np.empty_like(binarized)
bboxes = []
for i in range(binarized.shape[0]):
labeled = measure.label(binarized[i, :, :])
props = measure.regionprops_table(labeled, properties=('label', 'area', 'bbox'))
props_df = pd.DataFrame(props)
props_df = props_df.sort_values('area', ascending=False)
largest_region[i, :, :] = labeled == props_df.iloc[0]['label']
bboxes.append([props_df.iloc[0][f'bbox-{i}'] for i in range(4)])
fig = px.imshow(
largest_region,
animation_frame=0,
binary_string=True,
labels={'animation_frame': 'time point'},
)
plotly.io.show(fig)
绘制界面位置随时间的变化#
此分析的最后一步是绘制固液界面随时间变化的位置。这可以通过绘制 maxr（边界框中的第三个元素）随时间变化来实现，因为该值显示了界面底部的 y 位置。此实验中的像素大小为 1.93 微米，帧率为每秒 80,000 帧，因此这些值用于将像素和图像编号转换为物理单位。我们通过将线性多项式拟合到散点图来计算平均凝固速度。速度是一阶系数。
ums_per_pixel = 1.93
fps = 80000
interface_y_um = [ums_per_pixel * bbox[2] for bbox in bboxes]
time_us = 1e6 / fps * np.arange(len(interface_y_um))
fig, ax = plt.subplots(dpi=100)
ax.scatter(time_us, interface_y_um)
c0, c1 = np.polynomial.polynomial.polyfit(time_us, interface_y_um, 1)
ax.plot(time_us, c1 * time_us + c0, label=f'Velocity: {abs(round(c1, 3))} m/s')
ax.set_title('S-L interface location vs. time')
ax.set_ylabel(r'Location ($\mu$m)')
ax.set_xlabel(r'Time ($\mu$s)')
plt.show()
脚本总运行时间：（0 分 4.572 秒）
