金属合金の凝固を追跡#

この例では、凝固中のニッケルベース合金の固液（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 画像スタックです。まず、画像処理の最初のステップを 3 次元データセット全体（つまり空間と時間全体）で実行するワークフローで可視化および分析します。これにより、形体がほぼ一定の位置（フレームからフレームへ）にある物理的特徴（例：気泡、はね）の除去に比べて、局在的で一時的なノイズの除去が優先されます。

fig = px.imshow(
    image_sequence,
    animation_frame=0,
    binary_string=True,
    labels={'animation_frame': 'time point'},
)
plotly.io.show(fig)

画像のデルタ値の計算#

最初に、ガウス低域通過フィルタを適用して画像を平滑にし、ノイズを軽減しましょう。次に、画像のデルタ値、つまり 2 つの連続したフレーム間の差のシーケンスを計算します。これを行うには、2 番目のフレームで始まる画像シーケンスから、2 番目に最後のフレームで終わる画像シーケンスを引きます。

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)

反転とノイズ低減#

追跡する領域（つまり、S-L インターフェイス）を含む領域に強度の最も高い領域が含まれるように、clipped 画像を反転します。次に、インターフェイスを超えるノイズを低減するために、合計変動ノイズ低減フィルタを適用します。

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)

2 値化#

次のステップは、バイナリ画像を作成して画像を前景と背景に分割することです。S-L インターフェイスをバイナリ画像の各前景にある最も目立つ特徴にすることで、最終的に画像の残りの部分から分離できるようになります。

バイナリ画像（binarized ) を作成するには、しきい値 thresh_val が必要です。これは手動で設定できますが、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 インターフェイスは連結されたピクセルの最大の領域として表示されます。ワークフローのこのステップでは、領域を時間の単一時点として扱う必要があるため、3 次元データセット全体ではなく、個々の 2 次元画像に対して処理を行います。

2 値画像に 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()

	ラベル	エリア	バウンディングボックス 0	バウンディングボックス 1	バウンディングボックス 2	バウンディングボックス 3
1	2	417.0	60	83	91	144
198	199	235.0	141	141	179	165
11	12	136.0	62	164	87	180
9	10	122.0	61	208	79	229
8	9	114.0	61	183	83	198

上記の（ソートされた）テーブルの最初の行にあるものとラベルを照合することにより、最大のリージョンを選択することができます。赤色のバウンディングボックス（bbox）と共に視覚化します。

largest_region_0 = labeled_0 == props_0_df.iloc[0]['label']
minr, minc, maxr, maxc = (props_0_df.iloc[0][f'bbox-{i}'] for i in range(4))
fig = px.imshow(largest_region_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)

同じバウンディングボックスを 0 番目の未加工画像に重ねることにより、ボックスの下限が S-L インターフェイスの下部にどのように揃っているかを視覚化できます。このバウンディングボックスは 0 番目の画像と 1 番目の画像の画像デルタから計算されましたが、ボックスの一番下のリージョンは、インターフェイスが上向きに移動しているため、時間的により前のインターフェイスの位置に対応します。

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 インターフェイスの位置を追跡するために使用される bbox 情報も格納します。

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 の 3 番目の要素）を時間に対してプロットすることで実現できます。この値はインターフェイスの下端部の y 位置を示すためです。この実験でのピクセルサイズは 1.93 ミクロンで、フレームレートは 1 秒あたり 80,000 フレームでした。したがって、これらの値はピクセルと画像番号を物理単位に変換するために使用されます。散布図に 1 次多項式をフィッティングすることで、平均凝固速度を計算します。速度は 1 次係数です。

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

ギャラリーは Sphinx-Gallery によって生成されました。