openpiv 是做什麼的？
概述
OpenPIV（開放粒子圖像測速）可從 PIV 圖像對中分析流體流動。它涵蓋了預處理、互相關、向量驗證、異常值替換、平滑以及縮放到物理單位。
以下所有內容均已針對 openpiv 0.25.4 版本進行驗證。API 在不同版本之間會有所變動——在使用不同版本的程式碼片段之前，請務必使用 inspect.signature() 進行檢查。
使用時機
在處理實驗 PIV 或流場可視化圖像對時使用此技能：測量 2D 速度場、調整詢問窗口參數、驗證向量，或推導渦量、應變率和湍流統計數據。如果是為了模擬流動而非測量它，請改用 CFD 技能。
快速入門
安裝 OpenPIV：
uv pip install openpiv
# 當分析需要可重現時，請鎖定版本——以下所有程式碼片段均針對此版本進行檢查。
uv pip install "openpiv==0.25.4"
在圖像對上執行 PIV 分析：
import numpy as np
from openpiv import tools, pyprocess, validation, filters, scaling
frame_a = tools.imread("image_a.bmp")
frame_b = tools.imread("image_b.bmp")
# 互相關。當 sig2noise_method 不為 None 時，返回 (u, v, s2n)。
u, v, s2n = pyprocess.extended_search_area_piv(
frame_a.astype(np.int32),
frame_b.astype(np.int32),
window_size=32,
overlap=12,
dt=0.02,
search_area_size=38,
correlation_method="linear", # 當 search_area_size > window_size 時必須使用
sig2noise_method="peak2peak",
)
x, y = pyprocess.get_coordinates(
image_size=frame_a.shape,
search_area_size=38,
overlap=12,
)
# flags 是一個布林陣列：True 表示虛假向量。
flags = validation.sig2noise_val(s2n, threshold=1.05)
u, v = filters.replace_outliers(u, v, flags, method="localmean", max_iter=3, kernel_size=2)
# 縮放到物理單位，然後翻轉以符合圖像座標進行繪圖。
x, y, u, v = scaling.uniform(x, y, u, v, scaling_factor=96.52)
x, y, u, v = tools.transform_coordinates(x, y, u, v)
tools.save("vectors.txt", x, y, u, v, flags)
或者使用內建的 CLI，它包裝了上述處理流程：
python skills/openpiv/scripts/runner.py \
--image frame_a.bmp --image frame_b.bmp --output_dir results --verbose
核心概念
PIV 基礎
粒子圖像測速是一種光學方法，通過追蹤兩張圖像之間被照亮的示蹤粒子來測量流體速度。
處理流程：
- 擷取間隔已知時間
dt的圖像對（frame_a、frame_b）。 - 將圖像分割成詢問窗口。
- 對匹配窗口執行互相關以找到峰值位移。
- 驗證向量（訊噪比、全局範圍、局部中位數）。
- 用插值替換虛假向量。
- 將像素位移縮放到物理單位。
詢問窗口參數
window_size — 以像素為單位的相關窗口（典型值為 16–128）。較大的窗口提供更好的相關性，但空間分辨率較低。
overlap — 相鄰窗口之間共享的像素數（典型值為 window_size 的 50–75%）。較高的重疊率會提高向量密度和計算成本，但相鄰向量會變得相關而非獨立。
search_area_size — 在第二幀中搜索的窗口。必須 ≥ window_size；大幾個像素可以容納較大的位移。將擴展搜索區域與 correlation_method="linear" 配對使用——預設的 "circular" 依賴於 FFT 環繞，會將較大的位移別名為較小的位移。請參閱 references/advanced_algorithms.md。
經驗法則：將最大位移保持在 window_size 的四分之一以下，並力求每個窗口中有 5–10 個粒子。
訊噪比
s2n 衡量相關峰值的顯著程度。sig2noise_method 控制其計算方式——"peak2mean"（函數預設值）或 "peak2peak"。這兩種方法使用的尺度不同，因此針對一種方法調整的閾值對另一種方法毫無意義。典型的 peak2peak 閾值範圍為 1.05–1.3。
flags = validation.sig2noise_val(s2n, threshold=1.05)
# flags 為布林值：True == 虛假。`~flags` 選擇良好的向量。
常見操作
動態遮罩
遮罩功能位於 openpiv.preprocess 中，而非 openpiv.masking 模組。它返回一個 (image, mask) 元組，並需要一個浮點數圖像。
from openpiv import preprocess
# method="edges" 用於暗色、邊緣銳利的物體；"intensity" 用於高對比度物體。
frame_a_masked, mask_a = preprocess.dynamic_masking(
frame_a.astype(np.float64), method="intensity", filter_size=7, threshold=0.005
)
frame_b_masked, mask_b = preprocess.dynamic_masking(
frame_b.astype(np.float64), method="intensity", filter_size=7, threshold=0.005
)
將返回的圖像輸入到相關步驟中——它已經將遮罩區域歸零。不要將原始幀與 mask 相乘：遮罩已經應用，而且對於 method="edges"，遮罩以 uint8 0/255 形式返回而非布林值，因此相乘會將圖像縮放 255 倍。
多級處理
多級處理（窗口變形）位於 openpiv.windef 中，由 PIVSettings 資料類別驅動。pyprocess 沒有多級處理的入口點。
import numpy as np
from openpiv import scaling, windef
settings = windef.PIVSettings()
settings.windowsizes = (64, 32, 16) # 每個級別一個條目，遞減（這也是預設值）
settings.overlap = (32, 16, 8) # 與 windowsizes 長度相同
settings.num_iterations = 3 # 實際運行的級別數量
settings.sig2noise_threshold = 1.05
x, y, u, v, flags = windef.simple_multipass(
frame_a.astype(np.int32), frame_b.astype(np.int32), settings
)
# 輸出單位為像素/幀——自行轉換。scaling.uniform 僅除以 scaling_factor，因此需要單獨應用 dt。
dt = 0.02
x, y, u, v = scaling.uniform(x, y, u, v, scaling_factor=96.52)
u, v = u / dt, v / dt
simple_multipass 已經執行了驗證、異常值替換、用零填充剩餘的 NaN，並調用了 transform_coordinates——不要重複這些步驟。
單位陷阱： PIVSettings 有 dt 和 scaling_factor 欄位，但 windef 從不使用它們——first_pass 調用 extended_search_area_piv 時不帶 dt，因此整個多級鏈以像素/幀為單位工作。設置 settings.dt = 0.02 不會改變返回的值。請按上述方法在事後進行轉換。
要精細控制各個級別，windef.first_pass 和 windef.multipass_img_deform 是更底層的構建模組。
驗證與後處理
驗證方法
每個驗證器都返回一個布林陣列，其中 True 標記了虛假向量。
# 訊噪比
flags = validation.sig2noise_val(s2n, threshold=1.05)
# 全局範圍——接受 (min, max) 元組，按位置或作為 u_thresholds/v_thresholds 傳遞。
flags = validation.global_val(u, v, (-300, 300), (-300, 300))
# 局部中位數——u_threshold 和 v_threshold 是必需的；size 是鄰域半寬。
flags = validation.local_median_val(u, v, u_threshold=30.0, v_threshold=30.0, size=1)
# 使用布林 OR 組合（不建議使用 np.maximum——這些是布林陣列）。
flags = (
validation.sig2noise_val(s2n, threshold=1.05)
| validation.global_val(u, v, (-300, 300), (-300, 300))
| validation.local_median_val(u, v, u_threshold=30.0, v_threshold=30.0)
)
請以 u 和 v 的單位設置這些閾值，而非像素/幀。
extended_search_area_piv 會除以 dt，因此當 dt=0.02 時，一個 3 像素/幀的位移會變成 150 像素/秒。上述閾值適用於這種情況；PIV 文獻和 PIVSettings.min_max_u_disp 中使用的 (-30, 30) 是一個像素/幀的限制，將其應用於像素/秒的輸出會拒絕整個場。請在縮放之前進行驗證，或者也將閾值乘以 1/dt。
異常值替換
u, v = filters.replace_outliers(
u, v, flags, method="localmean", max_iter=3, tol=1e-3, kernel_size=2
)
method 接受 "localmean"、"disk" 或 "distance"——且僅限這三種。無法識別的名稱不會被拒絕；它會回退到一個全零的內核，並默默地返回一個無用的場。請注意，替換填充了標記的位置並使用插值——如果你隨後將它們覆蓋為 NaN，則替換是浪費的。請選擇其中一個：
# 讓標記的向量完全排除在分析之外，而不是對它們進行插值。
u = np.where(flags, np.nan, u)
v = np.where(flags, np.nan, v)
平滑
平滑功能是 openpiv.smoothn.smoothn；沒有 openpiv.smooth 模組。它返回一個元組，其第一個元素是平滑後的場，並且它不接受 NaN 輸入。
from openpiv.smoothn import smoothn
u_smooth, *_ = smoothn(np.nan_to_num(u), s=0.5) # s: 值越大越平滑
v_smooth, *_ = smoothn(np.nan_to_num(v), s=0.5)
u_smooth = np.asarray(u_smooth)
可視化
向量場繪圖
display_vector_field 讀取一個已保存的向量文件，並在內部調用 plt.show()，因此對於批次運行請選擇非互動式後端。
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from openpiv import tools
fig, ax = plt.subplots(figsize=(8, 8))
tools.display_vector_field(
"vectors.txt",
ax=ax,
scaling_factor=96.52, # 與 scaling.uniform 中使用的係數相同，用於映射回圖像
scale=50,
width=0.0035,
on_img=True,
image_name="frame_a.bmp",
)
fig.savefig("vector_field.png", dpi=150, bbox_inches="tight")
plt.close(fig)
自定義可視化
import numpy as np
import matplotlib.pyplot as plt
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
mag = np.sqrt(u**2 + v**2)
for ax, field, title, cmap in [
(axes[0], mag, "速度大小", "viridis"),
(axes[1], u, "U 速度", "RdBu_r"),
(axes[2], v, "V 速度", "RdBu_r"),
]:
im = ax.imshow(field, cmap=cmap)
ax.set_title(title)
plt.colorbar(im, ax=ax)
fig.tight_layout()
fig.savefig("velocity_components.png")
plt.close(fig)
分析函數
scripts/analyze.py 將這些功能與 runner.py 寫入的 params.npz 打包在一起。它從保存的座標推斷物理網格間距，因此導數以單位長度輸出：
import sys
sys.path.insert(0, "skills/openpiv/scripts")
from analyze import PIVAnalyzer
piv = PIVAnalyzer("results/params.npz")
vorticity = piv.compute_vorticity() # dv/dx - du/dy
exx, eyy, exy = piv.compute_strain()
stats = piv.compute_statistics() # u_mean, v_mean, rms_u, rms_v, tke
piv.plot_vector_field(save_path="quiver.png")
如果你更願意內聯計算它們，這是獨立形式：
渦量
def compute_vorticity(u, v, dx=1.0, dy=None):
"""平面外渦量 dv/dx - du/dy。請傳遞物理網格間距，而不是 1.0。"""
dy = dx if dy is None else dy
return np.gradient(v, dx, axis=1) - np.gradient(u, dy, axis=0)
網格間距在物理單位中是 (window_size - overlap) / scaling_factor，因此保留 dx=1.0 會得到每個網格單元的渦量，而不是單位長度的渦量。
符號約定： runner.py 以 transform_coordinates 結束，該函數將網格重新標記為右手 y 向上框架，但保留圖像順序的行，因此保存的 y 隨著行索引的增加而減小。上述獨立形式假設相反，因此在 params.npz 場上，它們返回 -du/dy 並翻轉渦量和剪下應變的符號——請對 axis=0 導數取反，或者使用讀取保存座標方向的 PIVAnalyzer。
應變率
def compute_strain(u, v, dx=1.0, dy=None):
"""返回 2D 應變率張量的 (exx, eyy, exy)。"""
dy = dx if dy is None else dy
du_dx = np.gradient(u, dx, axis=1)
du_dy = np.gradient(u, dy, axis=0)
dv_dx = np.gradient(v, dx, axis=1)
dv_dy = np.gradient(v, dy, axis=0)
return du_dx, dv_dy, 0.5 * (du_dy + dv_dx)
湍流統計數據
def compute_statistics(u, v):
"""單幀空間統計數據。非雷諾分解。"""
u_prime = u - np.nanmean(u)
v_prime = v - np.nanmean(v)
rms_u, rms_v = np.nanstd(u_prime), np.nanstd(v_prime)
return {
"u_mean": np.nanmean(u),
"v_mean": np.nanmean(v),
"rms_u": rms_u,
"rms_v": rms_v,
"tke": 0.5 * (rms_u**2 + rms_v**2),
}
注意事項： 減去單個幀的空間平均值會測量空間方差，這僅在均勻場中才等於湍流強度。真正的雷諾分解需要一組圖像對：對時間軸進行平均，然後從每個實現中減去該平均場。
CLI 使用
# 基本運行
python skills/openpiv/scripts/runner.py \
--image img1.bmp --image img2.bmp --output_dir results --verbose
# 帶動態遮罩的調整參數
python skills/openpiv/scripts/runner.py \
--image frame_a.bmp \
--image frame_b.bmp \
--output_dir results \
--window_size 32 \
--overlap 12 \
--search_area 38 \
--dt 0.02 \
--scaling 96.52 \
--threshold 1.05 \
--mask dynamic \
--mask_method intensity \
--verbose
CLI 選項
| 選項 | 預設值 | 說明 |
|---|---|---|
--image | 必需 | 圖像文件；為配對指定恰好兩次 |
--output_dir | results | 輸出目錄（若不存在則創建） |
--window_size | 32 | 詢問窗口大小（像素） |
--overlap | 12 | 窗口重疊（像素） |
--search_area | 38 | 搜索區域大小（像素），必須 ≥ --window_size |
--dt | 0.02 | 幀之間的時間間隔（秒） |
--scaling | 96.52 | 縮放因子，每物理單位的像素數（例如 像素/毫米） |
--threshold | 1.05 | peak2peak 訊噪比閾值 |
--mask | none | none 或 dynamic（openpiv.preprocess.dynamic_masking） |
--mask_method | intensity | edges 或 intensity，在...上使用 |
...
[為翻譯而截斷]