注意
前往結尾下載完整範例程式碼。或透過 Binder 在您的瀏覽器中執行此範例
追蹤金屬合金的凝固過程#
在本範例中，我們識別並追蹤正在凝固的鎳基合金中的固液 (S-L) 介面。隨著時間追蹤凝固過程可以計算凝固速度。這對於表徵樣品的凝固結構非常重要，並將用於為金屬的增材製造研究提供資訊。影像序列由先進非鐵結構合金中心 (CANFSA) 使用同步加速器 X 射線成像在阿貢國家實驗室 (ANL) 的先進光子源 (APS) 獲得。此分析首次在會議上發表 [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)
資料集是一個 2D 影像堆疊，包含 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。首先，讓我們考慮工作流程此階段的第一個影像增量，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 秒）
