openpiv 是做什么的？
概述
OpenPIV（开源粒子图像测速）用于分析 PIV 图像对中的流体流动。它涵盖了预处理、互相关、矢量验证、异常矢量替换、平滑处理以及物理单位缩放。
以下所有内容均针对 openpiv 0.25.4 版本验证。不同版本之间 API 可能发生变化——在信任某一代码片段之前，请先核对 inspect.signature()。
使用时机
在处理实验 PIV 或流动可视化图像对时使用此技能：测量二维速度场、调整查询窗口参数、验证矢量，或推导涡量、应变率和湍流统计量。如果需要模拟流动而非测量，请使用 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)
或使用内置的命令行工具，它封装了与上述完全相同的流程：
python skills/openpiv/scripts/runner.py \
--image frame_a.bmp --image frame_b.bmp --output_dir results --verbose
核心概念
PIV 基础
粒子图像测速是一种通过跟踪两个图像之间被照亮的示踪粒子来测量流体速度的光学方法。
处理流程：
- 捕获一对图像（
frame_a,frame_b），间隔已知时间dt。 - 将图像划分为查询窗口。
- 对匹配的窗口进行互相关，找到峰值位移。
- 验证矢量（信噪比、全局范围、局部中值）。
- 用插值替换异常矢量。
- 将像素位移缩放到物理单位。
查询窗口参数
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
)
# 输出单位为每帧像素（PIXELS PER FRAME）——需要自行转换。scaling.uniform 仅除以缩放因子，因此需单独应用 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 填充为 0，并调用 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)
# 全局范围——接受 (最小值, 最大值) 元组，位置参数或通过 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):
"""返回二维应变率张量的 (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，用于动态掩膜 |