openpiv

¿Qué hace openpiv?

Descripción general

OpenPIV (Velocimetría por Imagen de Partículas Abierta) analiza el flujo de fluidos a partir de pares de imágenes PIV. Cubre el preprocesamiento, la correlación cruzada, la validación de vectores, el reemplazo de valores atípicos, el suavizado y el escalado a unidades físicas.

Todo lo que sigue está verificado contra openpiv 0.25.4. La API cambia entre versiones: comprueba inspect.signature() antes de confiar en un fragmento de código para una versión diferente.

Cuándo usarlo

Usa esta habilidad cuando trabajes con pares de imágenes PIV experimentales o de visualización de flujo: midiendo campos de velocidad 2D, ajustando parámetros de la ventana de interrogación, validando vectores o derivando vorticidad, tasa de deformación y estadísticas de turbulencia. Para simular flujo en lugar de medirlo, usa una habilidad de CFD.

Inicio rápido

Instala OpenPIV:

uv pip install openpiv

# Fíjalo cuando el análisis necesite ser reproducible: esta es la versión contra la que se verificó cada fragmento a continuación.
uv pip install "openpiv==0.25.4"

Ejecuta el análisis PIV en un par de imágenes:

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

# Correlación cruzada. Devuelve (u, v, s2n) siempre que sig2noise_method no sea 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",   # requerido 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 es un array booleano: True marca un vector espurio.
flags = validation.sig2noise_val(s2n, threshold=1.05)
u, v = filters.replace_outliers(u, v, flags, method="localmean", max_iter=3, kernel_size=2)

# Escala a unidades físicas, luego invierte coordenadas para graficar.
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)

O usa la CLI incluida, que envuelve exactamente ese flujo de trabajo:

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

Conceptos centrales

Fundamentos de PIV

La Velocimetría por Imagen de Partículas es un método óptico para medir la velocidad de un fluido rastreando partículas trazadoras iluminadas entre dos imágenes.

Flujo del proceso:

1. Captura un par de imágenes (frame_a, frame_b) separadas por un tiempo conocido dt.
2. Divide las imágenes en ventanas de interrogación.
3. Aplica correlación cruzada a las ventanas coincidentes para encontrar el desplazamiento pico.
4. Valida vectores (relación señal-ruido, rango global, mediana local).
5. Reemplaza vectores espurios con valores interpolados.
6. Escala los desplazamientos en píxeles a unidades físicas.

Parámetros de la ventana de interrogación

window_size — ventana de correlación en píxeles (típicamente 16–128). Ventanas más grandes dan mejor correlación pero resolución espacial más gruesa.

overlap — píxeles compartidos entre ventanas adyacentes (típicamente 50–75% de window_size). Un solapamiento mayor aumenta la densidad de vectores y el costo, pero los vectores adyacentes se correlacionan en lugar de ser independientes.

search_area_size — la ventana que se busca en el segundo cuadro. Debe ser ≥ window_size; unos pocos píxeles más grandes acomodan desplazamientos mayores. Combina un área de búsqueda extendida con correlation_method="linear" — el "circular" por defecto se basa en el envolvente de la FFT y aliasea desplazamientos grandes en pequeños. Consulta references/advanced_algorithms.md.

Reglas generales: mantén el mayor desplazamiento por debajo de aproximadamente un cuarto de window_size, y apunta a 5–10 partículas por ventana.

Relación señal-ruido

s2n mide qué tan distinto es el pico de correlación. sig2noise_method controla cómo se calcula: "peak2mean" (el valor por defecto de la función) o "peak2peak". Ambos están en escalas diferentes, por lo que un umbral ajustado para uno no tiene sentido para el otro. Los umbrales típicos de peak2peak son 1.05–1.3.

flags = validation.sig2noise_val(s2n, threshold=1.05)
# flags es bool: True == espurio. `~flags` selecciona los vectores buenos.

Operaciones comunes

Enmascaramiento dinámico

El enmascaramiento vive en openpiv.preprocess, no en un módulo openpiv.masking. Devuelve una tupla (imagen, máscara) y espera una imagen de punto flotante.

from openpiv import preprocess

# method="edges" para objetos oscuros de bordes definidos; "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
)

Alimenta la imagen devuelta al paso de correlación: ya tiene la región enmascarada puesta a cero. No multipliques el cuadro original por mask: el enmascaramiento ya está aplicado, y para method="edges" la máscara se devuelve como uint8 0/255 en lugar de booleano, por lo que multiplicar reescala la imagen por 255.

Procesamiento de múltiples pasos

El procesamiento de múltiples pasos (deformación de ventana) vive en openpiv.windef, controlado por una clase de datos PIVSettings. pyprocess no tiene un punto de entrada para múltiples pasos.

import numpy as np
from openpiv import scaling, windef

settings = windef.PIVSettings()
settings.windowsizes = (64, 32, 16)   # una entrada por paso, decreciente (este es también el valor por defecto)
settings.overlap = (32, 16, 8)        # misma longitud que windowsizes
settings.num_iterations = 3           # número de pasos a ejecutar realmente
settings.sig2noise_threshold = 1.05

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

# La salida está en PÍXELES POR CUADRO – convierte tú mismo. scaling.uniform solo divide
# por scaling_factor, así que aplica dt por separado.
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 ya valida, reemplaza valores atípicos, llena los NaN restantes con ceros y llama a transform_coordinates – no repitas esos pasos.

Trampa de unidades: PIVSettings tiene campos dt y scaling_factor, pero windef nunca los usa – first_pass llama a extended_search_area_piv sin dt, por lo que toda la cadena de múltiples pasos trabaja en píxeles por cuadro. Establecer settings.dt = 0.02 no cambia nada en los valores devueltos. Convierte después del hecho, como se muestra arriba.

Para tener control sobre pasos individuales, windef.first_pass y windef.multipass_img_deform son los bloques de construcción de nivel inferior.

Validación y posprocesamiento

Métodos de validación

Cada validador devuelve un array booleano donde True marca un vector espurio.

# Relación señal-ruido
flags = validation.sig2noise_val(s2n, threshold=1.05)

# Rango global – toma tuplas (min, max), posicionalmente o como u_thresholds/v_thresholds.
flags = validation.global_val(u, v, (-300, 300), (-300, 300))

# Mediana local – u_threshold y v_threshold son REQUERIDOS; size es la mitad del ancho del vecindario.
flags = validation.local_median_val(u, v, u_threshold=30.0, v_threshold=30.0, size=1)

# Combínalos con OR booleano (no np.maximum – estos son arrays booleanos).
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)
)

Establece estos umbrales en las unidades de u y v, no en píxeles por cuadro. extended_search_area_piv divide por dt, por lo que con dt=0.02 un desplazamiento de 3 px/cuadro llega como 150 px/s. Los umbrales anteriores se adaptan a ese caso; la cifra (-30, 30) que la literatura PIV y PIVSettings.min_max_u_disp usan es un límite de px/cuadro, y aplicarlo a la salida en px/s rechaza todo el campo. Valida antes de escalar, o escala los umbrales por 1/dt también.

Reemplazo de valores atípicos

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

method acepta "localmean", "disk" o "distance" – y solo esos tres. Un nombre no reconocido no se rechaza; pasa a un kernel de todos ceros y devuelve silenciosamente un campo inútil. Ten en cuenta que el reemplazo llena las posiciones marcadas con valores interpolados – si luego los sobrescribes con NaN, el reemplazo fue en vano. Elige uno u otro:

# Mantén los vectores marcados fuera del análisis por completo, en lugar de interpolar.
u = np.where(flags, np.nan, u)
v = np.where(flags, np.nan, v)

Suavizado

El suavizado es openpiv.smoothn.smoothn; no hay un módulo openpiv.smooth. Devuelve una tupla cuyo primer elemento es el campo suavizado, y no acepta entrada NaN.

from openpiv.smoothn import smoothn

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

Visualización

Graficación de campos vectoriales

display_vector_field lee un archivo de vectores guardado y llama a plt.show() internamente, por lo que selecciona un backend no interactivo para ejecuciones por lotes.

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,   # mismo factor usado en scaling.uniform, para mapear de vuelta a la imagen
    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)

Visualización 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, "Magnitud de velocidad", "viridis"),
    (axes[1], u, "Velocidad U", "RdBu_r"),
    (axes[2], v, "Velocidad 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)

Funciones de análisis

scripts/analyze.py agrupa estas funciones contra un params.npz escrito por runner.py. Infiere el espaciado de la cuadrícula física a partir de las coordenadas guardadas, por lo que las derivadas resultan por unidad de longitud:

import sys
sys.path.insert(0, "skills/openpiv/scripts")
from analyze import PIVAnalyzer

piv = PIVAnalyzer("results/params.npz")
vorticidad = piv.compute_vorticity()          # dv/dx - du/dy
exx, eyy, exy = piv.compute_strain()
estadisticas = piv.compute_statistics()       # u_mean, v_mean, rms_u, rms_v, tke
piv.plot_vector_field(save_path="quiver.png")

Las formas independientes, si prefieres calcularlas en línea:

Vorticidad

def compute_vorticity(u, v, dx=1.0, dy=None):
    """Vorticidad fuera del plano dv/dx - du/dy. Pasa el espaciado de la cuadrícula física, no 1.0."""
    dy = dx if dy is None else dy
    return np.gradient(v, dx, axis=1) - np.gradient(u, dy, axis=0)

El espaciado de la cuadrícula es (window_size - overlap) / scaling_factor en unidades físicas, por lo que dejar dx=1.0 produce vorticidad por celda de la cuadrícula, no por unidad de longitud.

Convención de signos: runner.py termina con transform_coordinates, que reetiqueta la cuadrícula en un marco y-hacia-arriba diestro pero deja las filas en orden de imagen, por lo que la y guardada disminuye a medida que el índice de fila aumenta. Las formas independientes anteriores asumen lo contrario, por lo que en un campo de params.npz devuelven -du/dy e invierten el signo de la vorticidad y la deformación cortante – niega las derivadas del axis=0, o usa PIVAnalyzer, que lee la orientación de las coordenadas guardadas.

Tasa de deformación

def compute_strain(u, v, dx=1.0, dy=None):
    """Devuelve (exx, eyy, exy) del tensor de tasa de deformación 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)

Estadísticas de turbulencia

def compute_statistics(u, v):
    """Estadísticas espaciales de un solo cuadro. NO es descomposición 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),
    }

Advertencia: restar la media espacial de un cuadro mide la varianza espacial, que equivale a la intensidad turbulenta solo para un campo homogéneo. La descomposición genuina de Reynolds necesita un conjunto de pares de imágenes: promedia sobre el eje temporal, luego resta ese campo medio de cada realización.

Uso de la CLI

# Ejecución básica
python skills/openpiv/scripts/runner.py \
    --image img1.bmp --image img2.bmp --output_dir results --verbose

# Parámetros ajustados con enmascaramiento 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

Opciones de la CLI

Opción	Valor por defecto	Descripción
--image	requerido	Archivo de imagen; especifica exactamente dos veces para el par
--output_dir	results	Directorio de salida (se crea si no existe)
--window_size	32	Tamaño de la ventana de interrogación (px)
--overlap	12	Solapamiento de la ventana (px)
--search_area	38	Tamaño del área de búsqueda (px), debe ser ≥ --window_size
--dt	0.02	Tiempo entre cuadros (s)
--scaling	96.52	Factor de escalado, píxeles por unidad física (ej. px/mm)
--threshold	1.05	Umbral de relación señal-ruido peak2peak
--mask	none	none o dynamic (openpiv.preprocess.dynamic_masking)
--mask_method	intensity	edges o intensity, usado solo si --mask dynamic
--filter_size	7	Tamaño del filtro para el enmascaramiento
--mask_threshold	0.005	Umbral para el enmascaramiento
--verbose	False	Salida detallada

Referencias

• OpenPIV Documentation
• references/advanced_algorithms.md – Correlación extendida, deformación de ventana y teoría de multipaso
• Raffel et al., Particle Image Velocimetry: A Practical Guide