备注
转到结尾 以下载完整示例代码。或者通过 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帧（时间点）的二维图像堆栈。我们在一个工作流程中对其进行可视化和分析，其中首先对整个三维数据集（即跨越空间和时间）执行图像处理步骤，以便优先去除局部、瞬态噪声，而不是物理特征（例如，气泡、飞溅等），这些物理特征在从一帧到下一帧中大致处于相同位置。
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 。首先，让我们考虑工作流这一阶段的第一个图像 delta，即 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)
我们现在准备好为序列中的所有图像增量执行此选择。我们还将存储边界框信息，这些信息将用于跟踪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.786 秒)
