注意
转到末尾 下载完整的示例代码。或通过 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)
绘制界面位置随时间的变化#
本分析的最后一步是绘制 S-L 界面位置随时间的变化。这可以通过绘制 maxr（bbox 中的第三个元素）随时间的变化来实现，因为该值显示了界面的底部在 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 分钟 2.622 秒）
