openpiv

O que openpiv faz?

Visão Geral

OpenPIV (Open Particle Image Velocimetry) analisa o fluxo de fluidos a partir de pares de imagens PIV. Ele cobre pré-processamento, correlação cruzada, validação de vetores, substituição de outliers, suavização e escalonamento para unidades físicas.

Tudo abaixo é verificado em relação ao openpiv 0.25.4. A API muda entre versões — verifique inspect.signature() antes de confiar em um trecho contra uma versão diferente.

Quando usar

Use este skill ao trabalhar com pares de imagens experimentais de PIV ou visualização de fluxo: medindo campos de velocidade 2D, ajustando parâmetros de janela de interrogação, validando vetores, ou derivando vorticidade, taxa de deformação e estatísticas de turbulência. Para simular fluxo em vez de medi-lo, use um skill de CFD.

Início Rápido

Instale o OpenPIV:

uv pip install openpiv

# Fixe a versão quando a análise precisar ser reproduzível -- esta é a versão na qual cada
# trecho abaixo foi verificado.
uv pip install "openpiv==0.25.4"

Execute a análise PIV em um par de imagens:

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

# Correlação cruzada. Retorna (u, v, s2n) sempre que sig2noise_method não for 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",   # necessário para 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 é um array booleano: True marca um vetor espúrio.
flags = validation.sig2noise_val(s2n, threshold=1.05)
u, v = filters.replace_outliers(u, v, flags, method="localmean", max_iter=3, kernel_size=2)

# Escalone para unidades físicas e inverta para coordenadas da imagem para plotagem.
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)

Ou use a CLI incluída, que encapsula exatamente esse pipeline:

python skills/openpiv/scripts/runner.py \
    --image frame_a.bmp --image frame_b.bmp --output_dir results --verbose

Conceitos Principais

Fundamentos do PIV

A Velocimetria por Imagem de Partículas é um método óptico para medir a velocidade de fluidos rastreando partículas traçadoras iluminadas entre duas imagens.

Fluxo do processo:

1. Capture um par de imagens (frame_a, frame_b) separadas por um tempo conhecido dt.
2. Divida as imagens em janelas de interrogação.
3. Correlacione cruzadamente as janelas correspondentes para encontrar o deslocamento de pico.
4. Valide os vetores (relação sinal-ruído, faixa global, mediana local).
5. Substitua vetores espúrios com valores interpolados.
6. Escalone os deslocamentos em pixels para unidades físicas.

Parâmetros da Janela de Interrogação

window_size — janela de correlação em pixels (tipicamente 16–128). Janelas maiores fornecem melhor correlação, mas resolução espacial mais grosseira.

overlap — pixels compartilhados entre janelas adjacentes (tipicamente 50–75% de window_size). Maior sobreposição aumenta a densidade de vetores e o custo, mas vetores adjacentes se tornam correlacionados em vez de independentes.

search_area_size — a janela pesquisada no segundo quadro. Deve ser ≥ window_size; alguns pixels a mais acomodam deslocamentos maiores. Combine uma área de busca estendida com correlation_method="linear" — o padrão "circular" depende da envoltória da FFT e aliasa grandes deslocamentos em pequenos. Veja references/advanced_algorithms.md.

Regras práticas: mantenha o maior deslocamento abaixo de cerca de um quarto de window_size e busque 5–10 partículas por janela.

Relação Sinal-Ruído

s2n mede quão distinto é o pico de correlação. sig2noise_method controla como é calculado — "peak2mean" (o padrão da função) ou "peak2peak". Os dois estão em escalas diferentes, então um limiar ajustado para um é sem sentido para o outro. Limiares típicos de peak2peak são 1.05–1.3.

flags = validation.sig2noise_val(s2n, threshold=1.05)
# flags é bool: True == espúrio. `~flags` seleciona os vetores bons.

Operações Comuns

Mascaramento Dinâmico

O mascaramento reside em openpiv.preprocess, não em um módulo openpiv.masking. Ele retorna uma tupla (image, mask) e espera uma imagem float.

from openpiv import preprocess

# method="edges" para objetos escuros e de bordas nítidas; "intensity" para objetos de alto contraste.
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
)

Alimente a imagem retornada na etapa de correlação — ela já tem a região mascarada zerada. Não multiplique o quadro original por mask: o mascaramento já foi aplicado, e para method="edges" a máscara retorna como uint8 0/255 em vez de booleano, então multiplicar redimensiona a imagem por 255.

Processamento Multi-Passagem

O processamento multi-passagem (deformação de janela) reside em openpiv.windef, dirigido por um dataclass PIVSettings. pyprocess não tem ponto de entrada multi-passagem.

import numpy as np
from openpiv import scaling, windef

settings = windef.PIVSettings()
settings.windowsizes = (64, 32, 16)   # uma entrada por passagem, decrescente (este também é o padrão)
settings.overlap = (32, 16, 8)        # mesmo comprimento que windowsizes
settings.num_iterations = 3           # número de passagens para realmente executar
settings.sig2noise_threshold = 1.05

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

# A saída está em PIXELS POR QUADRO -- converta você mesmo. scaling.uniform apenas divide
# por scaling_factor, então aplique dt separadamente.
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 já valida, substitui outliers, preenche NaNs restantes com zeros e chama transform_coordinates — não repita esses passos.

Armadilha de unidades: PIVSettings tem campos dt e scaling_factor, mas windef nunca usa nenhum dos dois — first_pass chama extended_search_area_piv sem dt, então toda a cadeia multi-passagem funciona em pixels por quadro. Definir settings.dt = 0.02 não muda nada nos valores retornados. Converta depois do fato, como acima.

Para controle sobre passagens individuais, windef.first_pass e windef.multipass_img_deform são os blocos de construção de nível inferior.

Validação e Pós-Processamento

Métodos de Validação

Cada validador retorna um array booleano onde True marca um vetor espúrio.

# Relação sinal-ruído
flags = validation.sig2noise_val(s2n, threshold=1.05)

# Faixa global -- recebe tuplas (min, max), posicionalmente ou como u_thresholds/v_thresholds.
flags = validation.global_val(u, v, (-300, 300), (-300, 300))

# Mediana local -- u_threshold e v_threshold são OBRIGATÓRIOS; size é a metade da largura da vizinhança.
flags = validation.local_median_val(u, v, u_threshold=30.0, v_threshold=30.0, size=1)

# Combine com OR booleano (não np.maximum -- estes são arrays bool).
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)
)

Defina esses limiares nas unidades de u e v, não em pixels por quadro. extended_search_area_piv divide por dt, então com dt=0.02 um deslocamento de 3 px/quadro chega como 150 px/s. Os limiares acima se adequam a esse caso; o valor (-30, 30) que a literatura PIV e PIVSettings.min_max_u_disp usam é um limite em px/quadro, e aplicá-lo à saída em px/s rejeita todo o campo. Valide antes de escalonar, ou escale os limiares por 1/dt também.

Substituição de Outliers

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

method aceita "localmean", "disk" ou "distance" — e apenas esses três. Um nome não reconhecido não é rejeitado; ele cai em um kernel todo-zero e silenciosamente retorna um campo inútil. Note que a substituição preenche as posições sinalizadas com valores interpolados — se você então sobrescrevê-los com NaN, a substituição foi desperdiçada. Escolha um ou outro:

# Mantenha vetores sinalizados fora da análise completamente, em vez de interpolá-los.
u = np.where(flags, np.nan, u)
v = np.where(flags, np.nan, v)

Suavização

A suavização é openpiv.smoothn.smoothn; não há um módulo openpiv.smooth. Ela retorna uma tupla cujo primeiro elemento é o campo suavizado, e não aceita entrada com NaN.

from openpiv.smoothn import smoothn

u_smooth, *_ = smoothn(np.nan_to_num(u), s=0.5)  # s: maior == mais suave
v_smooth, *_ = smoothn(np.nan_to_num(v), s=0.5)
u_smooth = np.asarray(u_smooth)

Visualização

Plotagem de Campos de Vetores

display_vector_field lê um arquivo de vetores salvo e chama plt.show() internamente, então selecione um backend não interativo para execuções em lote.

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,   # mesmo fator usado em scaling.uniform, para mapear de volta na imagem
    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)

Visualização Personalizada

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, "Magnitude da Velocidade", "viridis"),
    (axes[1], u, "Velocidade U", "RdBu_r"),
    (axes[2], v, "Velocidade 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)

Funções de Análise

scripts/analyze.py agrupa estas contra um params.npz escrito por runner.py. Ele infere o espaçamento da grade física a partir das coordenadas salvas, para que as derivadas saiam por unidade de comprimento:

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

As formas independentes, se você preferir calculá-las inline:

Vorticidade

def compute_vorticity(u, v, dx=1.0, dy=None):
    """Vorticidade fora do plano dv/dx - du/dy. Passe o espaçamento da grade física, não 1.0."""
    dy = dx if dy is None else dy
    return np.gradient(v, dx, axis=1) - np.gradient(u, dy, axis=0)

O espaçamento da grade é (window_size - overlap) / scaling_factor em unidades físicas, então deixar dx=1.0 produz vorticidade por célula da grade, não por unidade de comprimento.

Convenção de sinal: runner.py termina com transform_coordinates, que rotula a grade em um referencial destro com y para cima, mas deixa as linhas na ordem da imagem, então o y salvo diminui conforme o índice da linha cresce. As formas independentes acima assumem o oposto, então em um campo params.npz elas retornam -du/dy e invertem o sinal da vorticidade e da deformação cisalhante — negue as derivadas em axis=0, ou use PIVAnalyzer, que lê a orientação das coordenadas salvas.

Taxa de Deformação

def compute_strain(u, v, dx=1.0, dy=None):
    """Retorna (exx, eyy, exy) do tensor de taxa de deformação 2D."""
    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)

Estatísticas de Turbulência

def compute_statistics(u, v):
    """Estatísticas espaciais de quadro único. NÃO é decomposição de Reynolds."""
    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),
    }

Ressalva: subtrair a média espacial de um quadro mede a variância espacial, que é igual à intensidade turbulenta apenas para um campo homogêneo. A decomposição genuína de Reynolds precisa de um conjunto de pares de imagens: faça a média sobre o eixo do tempo, depois subtraia esse campo médio de cada realização.

Uso da CLI

# Execução básica
python skills/openpiv/scripts/runner.py \
    --image img1.bmp --image img2.bmp --output_dir results --verbose

# Parâmetros ajustados com mascaramento dinâmico
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

Opções da CLI

Opção	Padrão	Descrição
--image	obrigatório	Arquivo de imagem; especifique exatamente duas vezes para o par
--output_dir	results	Diretório de saída (criado se ausente)
--window_size	32	Tamanho da janela de interrogação (px)
--overlap	12	Sobreposição da janela (px)
--search_area	38	Tamanho da área de busca (px), deve ser ≥ --window_size
--dt	0.02	Tempo entre quadros (s)
--scaling	96.52	Fator de escalonamento, pixels por unidade física (ex.: px/mm)
--threshold	1.05	Limiar de relação sinal-ruído peak2peak
--mask	none	none ou dynamic (openpiv.preprocess.dynamic_masking)
--mask_method	intensity	edges ou intensity, usado em conjunto com --mask dynamic
--verbose	False	Saída detalhada

O runner salva todos os resultados (vetores, mapa de vorticidade, mapa de estatísticas) em --output_dir. Ele espera que as imagens estejam em formato BMP com tons de cinza de 8 bits. Para formatos ou níveis de cinza diferentes, converta antes.