openpiv

OpenPIV

概述

OpenPIV（Open Particle Image Velocimetry，开源粒子图像测速）可从 PIV 图像对中分析流体流动。它涵盖 预处理、互相关、矢量验证、异常值替换、平滑处理，以及到物理单位的尺度换算。

面向 openpiv 0.26.1。合成检查在 Python 3.13、NumPy 2.5.3、SciPy 1.18.1、 scikit-image 0.26.0 和 Matplotlib 3.11.2 上运行。未安装 Rust；请使用 backend="scipy" 走 经过测试的路径。参见审查证据。

何时使用

在处理实验 PIV 或流动可视化图像对时使用本技能：测量 2D 速度场、调整查询窗口参数、验证矢量，或推导涡量、 应变率和湍流统计量。若要模拟流动而非测量流动，请改用 CFD 技能。

快速开始

安装 OpenPIV：

uv pip install openpiv

# Pin it when the analysis needs to be reproducible -- this is the version every
# snippet below was checked against.
uv pip install "openpiv==0.26.1"

对图像对运行 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")

# Always returns (u, v, s2n); s2n is all NaN if sig2noise_method=None.
u, v, s2n = pyprocess.extended_search_area_piv(
    frame_a.astype(np.float32),
    frame_b.astype(np.float32),
    window_size=32,
    overlap=12,
    dt=0.02,
    search_area_size=38,
    correlation_method="linear",   # circular also supports extended search
    normalized_correlation=True,
    backend="scipy",
    sig2noise_method="peak2peak",
)

x, y = pyprocess.get_coordinates(
    image_size=frame_a.shape,
    search_area_size=38,
    overlap=12,
    center_on_field=False,         # matches sliding-window positions
)

# flags is a boolean array: True marks a spurious vector.
flags = validation.sig2noise_val(s2n, threshold=1.05) | ~np.isfinite(s2n)
flags |= ~np.isfinite(u) | ~np.isfinite(v)
u, v = filters.replace_outliers(u, v, flags, method="localmean", max_iter=3, kernel_size=2)

# Scale, then convert to right-handed image-boundary coordinates.
x, y, u, v = scaling.uniform(x, y, u, v, scaling_factor=96.52)
y, v = frame_a.shape[0] / 96.52 - y, -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 基础

粒子图像测速是一种光学方法，通过跟踪两幅图像之间被照亮的 示踪粒子来测量流体速度。

处理流程：

1. 采集相隔已知时间 dt 的一对图像（frame_a、frame_b）。
2. 将图像划分为查询窗口。
3. 对匹配的窗口进行互相关，以找到峰值位移。
4. 验证矢量（信噪比、全局范围、局部中位数）。
5. 用插值替换虚假矢量。
6. 将像素位移换算为物理单位。

查询窗口参数

window_size —— 以像素为单位的相关窗口（通常为 16–128）。窗口越大，相关性越好， 但空间分辨率越粗。

overlap —— 相邻窗口之间共享的像素数（通常为 window_size 的 50–75%）。更高的 重叠度会提高矢量密度和成本，但相邻矢量会变得相关，而非 独立。

search_area_size —— 在第二帧中搜索的窗口。必须 ≥ window_size；大几个 像素可容纳更大的位移。"circular" 和 "linear" 都支持扩展 搜索。CLI 为其选择零填充的 "linear" 并设置 normalized_correlation=True；这 并不能恢复任意位移或消除虚假峰值。网格步长为 search_area_size - overlap，因此重叠必须小于搜索区域。

经验法则：将最大位移保持在 window_size 的大约四分之一以下，并力争每个窗口 有 5–10 个粒子。

信噪比

s2n 衡量相关峰的显著程度。sig2noise_method 控制其计算方式—— "peak2mean"（函数默认值）或 "peak2peak"。两者的尺度不同，因此为其中一种 调校的阈值对另一种毫无意义。典型的 peak2peak 阈值为 1.05–1.3。

flags = validation.sig2noise_val(s2n, threshold=1.05)
# flags is bool: True == spurious. `~flags` selects the good vectors.

常见操作

动态掩码

掩码功能位于 openpiv.preprocess 中，返回 (image, mask)。使用浮点副本 以避免无符号减法的伪影。在分析前检查两帧上的掩码。

from openpiv import preprocess

# Intensity masking uses Otsu; the threshold argument is for edges.
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_a | mask_b 作为被排除的物理区域。CLI 在 查询中心对其采样、保存它，并在替换后将这些速度保持为 NaN。跨越物体的窗口 仍可能有偏差：为实验检查/膨胀掩码。在 0.26.1 中， method="edges" 使用 uint8 而非 bool 进行索引，可能损坏图像或失败。 CLI 拒绝使用它；请使用强度掩码或在自定义工作流中使用经过验证的布尔掩码。

多遍处理

多遍（窗口变形）位于 openpiv.windef 中，由 PIVSettings 数据类驱动。 pyprocess 没有多遍处理入口。

import numpy as np
from openpiv import scaling, windef

settings = windef.PIVSettings(
    windowsizes=(64, 32, 16), overlap=(32, 16, 8), num_iterations=3,
    backend="scipy", sig2noise_method="peak2peak", sig2noise_threshold=1.05,
)

x, y, u, v, flags = windef.simple_multipass(
    frame_a.astype(np.float32), frame_b.astype(np.float32), settings
)

# Output is in PIXELS PER FRAME -- convert yourself. scaling.uniform only divides
# by scaling_factor, so apply dt separately.
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 填充为零，并转换坐标。请保留它的这些标记：修复/补零后的输出并非独立测量结果。 始终传入 settings：它的无 settings 路径在 0.26.1 中有两个窗口尺寸，但有三轮迭代。

单位陷阱： simple_multipass 及其别名 multigrid_windef 会忽略 settings.dt 和 settings.scaling_factor；请按上文所述对其数组进行转换。批量调用的 windef.piv(settings) 会在保存前应用 这两者——请勿再对其文件进行缩放。这个简单封装并不会遵循所有批量 预处理/输出开关。

若想控制各个单独遍次，windef.first_pass 和 windef.multipass_img_deform 是 更底层的构建模块。

校验与后处理

校验方法

每个校验器都会返回一个布尔数组，其中 True 标记伪矢量。同时还要拒绝 非有限的 u、v 和 s2n：仅靠 sig2noise_val 并不会标记 NaN。

# Signal-to-noise
flags = validation.sig2noise_val(s2n, threshold=1.05)

# Global range -- takes (min, max) TUPLES, positionally or as u_thresholds/v_thresholds.
flags = validation.global_val(u, v, (-300, 300), (-300, 300))

# Local median -- u_threshold and v_threshold are REQUIRED; size is the neighbourhood half-width.
flags = validation.local_median_val(u, v, u_threshold=30.0, v_threshold=30.0, size=1)

# Combine with boolean OR (not np.maximum -- these are bool arrays).
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 px/帧的位移会以 150 px/s 的形式到达。上面的阈值适用于这种情况；PIV 文献和 PIVSettings.min_max_u_disp 所使用的 (-30, 30) 数值是一个 px/帧限制，若将其应用于 px/s 输出，会拒绝 整个矢量场。请以 dt=1 进行相关计算，以便在转换前校验位移，或者将 阈值调整为实际速度单位（同时包含时间和校准系数）。

异常值替换

u, v = filters.replace_outliers(
    u, v, flags, method="localmean", max_iter=3, tol=1e-3, kernel_size=2
)

method 接受 "localmean"、"disk" 或 "distance"；无效名称在填充缺失单元格时会引发 ValueError。 对于 "distance"，请使用 kernel_size>=2：在 0.26.1 中尺寸 1 会把所有邻居 权重截断为零并留下空洞。替换会填充被标记的 位置并写入插值结果——如果你随后用 NaN 覆盖它们，这次替换就 白费了。请二选一：

# Keep flagged vectors out of the analysis entirely, instead of interpolating them.
u = np.where(flags, np.nan, u)
v = np.where(flags, np.nan, v)

平滑

平滑是 openpiv.smoothn.smoothn；不存在 openpiv.smooth 模块。它返回一个元组， 其第一个元素是平滑后的场。NaN/Inf 表示缺失观测；请传入副本，因为 它可能会修改输入。先对空洞补零会错误地将它们视为已测量的零。

from openpiv.smoothn import smoothn

u_smooth, *_ = smoothn(np.asarray(u, dtype=float).copy(), s=0.5)  # larger == smoother
v_smooth, *_ = smoothn(np.asarray(v, dtype=float).copy(), 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,   # same factor used in scaling.uniform, to map back onto the image
    scale=50,
    width=0.0035,
    on_img=False,
)
image = tools.imread("frame_a.bmp")
height, width = image.shape
ax.imshow(image, cmap="gray", origin="upper", zorder=-1,
          extent=(0, width / 96.52, 0, height / 96.52))
fig.savefig("vector_field.png", dpi=150, bbox_inches="tight")
plt.close(fig)

请显式设置图像范围：该辅助函数的 on_img=True 会从最后一个矢量推断范围， 当窗口留下未使用边距时可能会拉伸叠加图层。命令行界面会使用实际图像边界进行绘图。

自定义可视化

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, "Velocity Magnitude", "viridis"),
    (axes[1], u, "U Velocity", "RdBu_r"),
    (axes[2], v, "V Velocity", "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()             # excludes flagged and masked vectors
piv.plot_vector_field(save_path="quiver.png")

如果你更愿意直接内联计算，也可以使用独立形式：

涡量

def compute_vorticity(u, v, dx=1.0, dy=None):
    """Out-of-plane vorticity dv/dx - du/dy. Pass the physical grid spacing, not 1.0."""
    dy = dx if dy is None else dy
    return np.gradient(v, dx, axis=1) - np.gradient(u, dy, axis=0)

对于单遍扩展搜索网格，间距为 (search_area_size - overlap) / scaling_factor，以物理单位表示；仅当两个窗口大小相同时，才简化为 (window_size - overlap) / scaling_factor。优先使用保存的 x 和 y 坐标之差，尤其是在多遍处理后。保留 dx=1.0 会得到每个网格单元的涡量，而非每单位长度的涡量。参见 OpenPIV 坐标生成。

符号约定： 命令行写入 y = image_height / scaling_factor - y 并对 v 取负，保持行按图像顺序排列。保存的 y 沿数组向下递减。独立形式需要带符号的 dy（此处为负），或使用 PIVAnalyzer，它会读取方向。正的 dy 会改变 y 导数的贡献，甚至可能抵消真实的旋转。

应变率

def compute_strain(u, v, dx=1.0, dy=None):
    """Return (exx, eyy, exy) of the 2D strain-rate tensor."""
    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):
    """Single-frame spatial statistics. NOT Reynolds decomposition."""
    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),
    }

注意事项： 减去单帧的空间均值衡量的是空间方差，只有在合理的均匀性/遍历性假设下才等于湍流强度。雷诺分解需要系综：从每个实现中减去其均值场。遗留的 tke 键只是双分量空间方差能量，并非完整的湍流动能。分析器默认排除被标记的向量；compute_statistics(include_interpolated=True) 会包含修复后的向量。

命令行用法

# Basic run
python skills/openpiv/scripts/runner.py \
    --image img1.bmp --image img2.bmp --output_dir results --verbose

# Tuned parameters with dynamic masking
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

命令行选项

选项	默认值	描述
--image	必需	图像文件；对于图像对必须恰好指定两次
--output_dir	results	输出目录（不存在时创建）
--window_size	32	查询窗口大小（px）
--overlap	12	搜索网格重叠（px），非负且小于搜索区域
--search_area	38	搜索区域大小（px），必须 ≥ --window_size
--dt	0.02	帧间时间（s）
--scaling	96.52	缩放因子，每物理单位的像素数（例如 px/mm）
--threshold	1.05	peak2peak 信噪比阈值
--mask	none	none 或 dynamic（openpiv.preprocess.dynamic_masking）
--mask_method	intensity	在 0.26.1 中仅 intensity 可用；edges 被拒绝
--backend	scipy	scipy、auto 或 rust；Rust 需要其扩展
--drop_invalid	关闭	将标记的向量置为 NaN，而不是保留插值
--verbose	关闭	打印进度消息

使用 OpenPIV 自带的图像对端到端验证安装：

python skills/openpiv/scripts/run_example.py --output_dir /tmp/openpiv-demo

输出文件

• vectors.txt — 制表符分隔，%.4e 格式，带有 # x y u v flags mask 注释头
• params.npz — x、y、u、v、flags、mask、s2n、计时/校准/窗口参数、请求的后端以及 OpenPIV 版本。坐标使用物理单位；速度使用单位/秒。
• vector_field.png — 绘制在第一帧上的向量场

# x	y	u	v	flags	mask
2.1757e-01	3.5226e+00	-6.2220e-02	-2.7081e+00	0.0000e+00	0.0000e+00
4.8695e-01	3.5226e+00	-3.1587e-01	-2.9800e+00	0.0000e+00	0.0000e+00

flags 为浮点数，0 表示有效，1 表示被标记。mask=1 标记被排除的对象区域，与 flags 无关。请与测量场一同保留两者。

最佳实践

参数选择

1. 窗口大小 — 32×32 适合大多数情况。64/128 用于在较粗分辨率下获得更好的相关性；16/24 用于更精细的分辨率，代价是噪声。
2. 重叠 — 窗口大小的 50–75%。
3. 阈值 — 提高它以拒绝更多向量；在切换 sig2noise_method 后务必重新调整。
4. 缩放因子 — 根据已知参考（如校准网格）进行校准，并保持单位一致（OpenPIV 的 test1 教程数据中的 96.52 是 px/mm）。

图像质量

• 粒子可见且分布均匀，每个查询窗口 5–10 个
• 无饱和或过曝区域
• 背景噪声最小；考虑在整个运行过程中进行背景减除

处理技巧

1. 从默认值开始，然后根据得到的向量场进行调优。
2. 检查 s2n 分布 —— 中位数偏低意味着相关性差，而不是阈值不佳。
3. 尽早可视化；明显的问题（向量均匀、边缘伪影）会立刻显现。
4. 对于存在大速度梯度或位移的流场，使用多遍（windef）。
5. 对反射和固体边界进行掩膜处理，而不是让它们生成向量。

资源

references/

• advanced_algorithms.md —— 相关性与亚像素方法、多遍窗口变形、 PIVSettings 字段、3D 与相分离模块

当需要详细的算法或设置信息时，加载该参考文档。