openpiv

OpenPIV

概述

OpenPIV（Open Particle Image Velocimetry）通过 PIV 图像对分析流体流动。它涵盖预处理、互相关、矢量验证、异常值替换、平滑处理，以及缩放到物理单位。

以下内容均已针对 openpiv 0.25.4 验证。不同版本的 API 可能有所变化；使用针对其他版本的代码片段前，请先检查 inspect.signature()。

适用场景

在处理实验 PIV 或流动可视化图像对时，可使用此技能：测量二维速度场、调整询问窗口参数、验证矢量，或推导涡量、应变率和湍流统计数据。如果要模拟流动而非测量流动，请改用 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.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")

# Cross-correlate. Returns (u, v, s2n) whenever sig2noise_method is not None.
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",   # required for 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 is a boolean array: True marks a spurious vector.
flags = validation.sig2noise_val(s2n, threshold=1.05)
u, v = filters.replace_outliers(u, v, flags, method="localmean", max_iter=3, kernel_size=2)

# Scale to physical units, then flip to image coordinates for plotting.
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 基础知识

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

处理流程：

1. 获取一对图像（frame_a、frame_b），两者间隔已知时间 dt。
2. 将图像划分为询问窗口。
3. 对匹配的窗口进行互相关，以找到峰值位移。
4. 验证矢量（信噪比、全局范围、局部中值）。
5. 使用插值后的值替换异常矢量。
6. 将像素位移缩放为物理单位。

询问窗口参数

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 is bool: True == spurious. `~flags` selects the good vectors.

常见操作

动态掩膜

掩膜功能位于 openpiv.preprocess 中，不在 openpiv.masking 模块中。它返回一个 (image, mask) 元组，并要求输入为浮点图像。

from openpiv import preprocess

# method="edges" for dark, sharp-edged objects; "intensity" for high-contrast objects.
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)   # one entry per pass, decreasing (this is also the default)
settings.overlap = (32, 16, 8)        # same length as windowsizes
settings.num_iterations = 3           # number of passes to actually run
settings.sig2noise_threshold = 1.05

x, y, u, v, flags = windef.simple_multipass(
    frame_a.astype(np.int32), frame_b.astype(np.int32), 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 填充为零，并调用 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 表示异常向量。

# 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/frame 的位移会变为 150 px/s。上述阈值适用于这种情况；PIV 文献和 PIVSettings.min_max_u_disp 使用的 (-30, 30) 数值是 px/frame 限值，用于 px/s 输出时会拒绝整个场。请在缩放前进行验证，或将阈值也乘以 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，替换就白费了。二者选其一：

# 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 输入。

from openpiv.smoothn import smoothn

u_smooth, *_ = smoothn(np.nan_to_num(u), s=0.5)  # s: larger == smoother
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,   # same factor used in scaling.uniform, to map back onto the image
    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, "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()             # 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):
    """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)

网格间距在物理单位下为 (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):
    """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),
    }

**注意：**减去单帧的空间均值测量的是空间方差，只有对于均匀场，它才等于湍流强度。真正的雷诺分解需要使用一组图像对：先沿时间轴求平均，再从每个实现中减去该均值场。

CLI 用法

# 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

CLI 选项

选项	默认值	描述
--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	edges 或 intensity，仅与 --mask dynamic 搭配使用
--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 数组的 NumPy 归档文件
• 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 表示被标记的矢量。

最佳实践

参数选择

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 和相位分离模块

需要详细了解算法或设置时，请加载此参考资料。

相关技能

• scvelo：使用 scVelo 对具有剪接和未剪接层的单细胞 RNA-seq AnnData 对象执行 RNA 速度分析（来自 velocyto、ST…）