OpenPIV

Que fait OpenPIV ?

Présentation

OpenPIV (Open Particle Image Velocimetry) analyse les écoulements de fluides à partir de paires d'images PIV. Il couvre le prétraitement, la corrélation croisée, la validation des vecteurs, le remplacement des valeurs aberrantes, le lissage et la mise à l'échelle en unités physiques.

Tout ce qui suit est vérifié pour openpiv 0.25.4. L'API change entre les versions — vérifiez inspect.signature() avant de faire confiance à un extrait avec une version différente.

Quand l'utiliser

Utilisez cette compétence lorsque vous travaillez avec des paires d'images PIV expérimentales ou de visualisation d'écoulement : mesurer des champs de vélocité 2D, ajuster les paramètres de fenêtre d'interrogation, valider des vecteurs ou dériver la vorticité, le taux de déformation et les statistiques de turbulence. Pour simuler un écoulement plutôt que de le mesurer, utilisez plutôt une compétence CFD.

Démarrage rapide

Installez OpenPIV :

uv pip install openpiv

# Figez la version lorsque l'analyse doit être reproductible -- c'est la version avec laquelle
# chaque extrait ci-dessous a été vérifié.
uv pip install "openpiv==0.25.4"

Exécutez une analyse PIV sur une paire d'images :

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

# Corrélation croisée. Renvoie (u, v, s2n) dès que sig2noise_method n'est pas 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",   # requis pour 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 est un tableau booléen : True marque un vecteur aberrant.
flags = validation.sig2noise_val(s2n, threshold=1.05)
u, v = filters.replace_outliers(u, v, flags, method="localmean", max_iter=3, kernel_size=2)

# Mise à l'échelle en unités physiques, puis basculement en coordonnées image pour le tracé.
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 utilisez l'interface en ligne de commande intégrée, qui encapsule exactement ce pipeline :

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

Concepts fondamentaux

Principes de base PIV

La vélocimétrie par images de particules est une méthode optique pour mesurer la vitesse d'un fluide en suivant les particules traceuses éclairées entre deux images.

Déroulement du processus :

1. Capturez une paire d'images (frame_a, frame_b) séparées par un temps connu dt.
2. Divisez les images en fenêtres d'interrogation.
3. Corrélez croisement les fenêtres correspondantes pour trouver le déplacement du pic.
4. Validez les vecteurs (rapport signal/bruit, plage globale, médiane locale).
5. Remplacez les vecteurs aberrants par des valeurs interpolées.
6. Mettez à l'échelle les déplacements en pixels en unités physiques.

Paramètres de la fenêtre d'interrogation

window_size — fenêtre de corrélation en pixels (généralement 16–128). Des fenêtres plus grandes donnent une meilleure corrélation mais une résolution spatiale plus grossière.

overlap — pixels partagés entre les fenêtres adjacentes (généralement 50–75 % de window_size). Un chevauchement plus élevé augmente la densité des vecteurs et le coût, mais les vecteurs adjacents deviennent corrélés plutôt qu'indépendants.

search_area_size — la fenêtre recherchée dans la deuxième image. Doit être ≥ window_size ; quelques pixels supplémentaires permettent des déplacements plus importants. Associez une zone de recherche étendue avec correlation_method="linear" — la valeur par défaut "circular" repose sur le repliement FFT et alias les grands déplacements en petits. Voir references/advanced_algorithms.md.

Règles empiriques : maintenez le plus grand déplacement à environ un quart de window_size, et visez 5 à 10 particules par fenêtre.

Rapport signal/bruit

s2n mesure la netteté du pic de corrélation. sig2noise_method contrôle la façon dont il est calculé — "peak2mean" (la valeur par défaut de la fonction) ou "peak2peak". Les deux sont sur des échelles différentes, donc un seuil ajusté pour l'un n'a pas de sens pour l'autre. Les seuils typiques pour peak2peak sont de 1,05 à 1,3.

flags = validation.sig2noise_val(s2n, threshold=1.05)
# flags est bool : True == aberrant. `~flags` sélectionne les bons vecteurs.

Opérations courantes

Masquage dynamique

Le masquage se trouve dans openpiv.preprocess, pas dans un module openpiv.masking. Il renvoie un tuple (image, mask) et attend une image flottante.

from openpiv import preprocess

# method="edges" pour les objets sombres aux bords nets ; "intensity" pour les objets à fort 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
)

Transmettez l'image renvoyée à l'étape de corrélation — elle a déjà la zone masquée mise à zéro. Ne multipliez pas l'image originale par mask : le masquage est déjà appliqué, et pour method="edges", le masque revient sous forme uint8 0/255 plutôt que booléen, donc la multiplication remet à l'échelle l'image par 255.

Traitement multi-passes

Le multi-passes (déformation de fenêtre) se trouve dans openpiv.windef, piloté par une classe de données PIVSettings. pyprocess n'a pas de point d'entrée multi-passes.

import numpy as np
from openpiv import scaling, windef

settings = windef.PIVSettings()
settings.windowsizes = (64, 32, 16)   # une entrée par passe, décroissante (c'est aussi la valeur par défaut)
settings.overlap = (32, 16, 8)        # même longueur que windowsizes
settings.num_iterations = 3           # nombre de passes à exécuter réellement
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 sortie est en PIXELS PAR TRAME -- convertissez vous-même. scaling.uniform ne fait que diviser
# par scaling_factor, donc appliquez dt séparément.
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 valide déjà, remplace les valeurs aberrantes, remplit les NaN restants avec des zéros et appelle transform_coordinates — ne répétez pas ces étapes.

Piège d'unités : PIVSettings a des champs dt et scaling_factor, mais windef ne les utilise jamais — first_pass appelle extended_search_area_piv sans dt, donc toute la chaîne multi-passes travaille en pixels par trame. Définir settings.dt = 0.02 ne change rien aux valeurs renvoyées. Convertissez après coup, comme ci-dessus.

Pour un contrôle sur les passes individuelles, windef.first_pass et windef.multipass_img_deform sont les blocs de construction de plus bas niveau.

Validation et post-traitement

Méthodes de validation

Chaque validateur renvoie un tableau booléen où True marque un vecteur aberrant.

# Rapport signal/bruit
flags = validation.sig2noise_val(s2n, threshold=1.05)

# Plage globale -- prend des TUPLES (min, max), positionnels ou sous forme de u_thresholds/v_thresholds.
flags = validation.global_val(u, v, (-300, 300), (-300, 300))

# Médiane locale -- u_threshold et v_threshold sont OBLIGATOIRES ; size est la demi-largeur du voisinage.
flags = validation.local_median_val(u, v, u_threshold=30.0, v_threshold=30.0, size=1)

# Combinez avec OU booléen (pas np.maximum -- ce sont des tableaux booléens).
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)
)

Définissez ces seuils dans les unités de u et v, pas en pixels par trame. extended_search_area_piv divise par dt, donc avec dt=0.02 un déplacement de 3 px/trame arrive à 150 px/s. Les seuils ci-dessus conviennent à ce cas ; la valeur (-30, 30) que la littérature PIV et PIVSettings.min_max_u_disp utilisent est une limite en px/trame, et l'appliquer à une sortie en px/s rejette la totalité du champ. Soit validez avant la mise à l'échelle, soit mettez également à l'échelle les seuils par 1/dt.

Remplacement des valeurs aberrantes

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

method accepte "localmean", "disk" ou "distance" — et seulement ces trois-là. Un nom non reconnu n'est pas rejeté ; il passe à un noyau tout zéro et renvoie silencieusement un champ inutile. Notez que le remplacement remplit les positions marquées avec des valeurs interpolées — si vous les écrasez ensuite avec NaN, le remplacement était inutile. Choisissez l'un ou l'autre :

# Gardez les vecteurs marqués hors de l'analyse complètement, plutôt que de les interpoler.
u = np.where(flags, np.nan, u)
v = np.where(flags, np.nan, v)

Lissage

Le lissage est openpiv.smoothn.smoothn ; il n'y a pas de module openpiv.smooth. Il renvoie un tuple dont le premier élément est le champ lissé, et il n'accepte pas d'entrée NaN.

from openpiv.smoothn import smoothn

u_smooth, *_ = smoothn(np.nan_to_num(u), s=0.5)  # s : plus grand == plus lisse
v_smooth, *_ = smoothn(np.nan_to_num(v), s=0.5)
u_smooth = np.asarray(u_smooth)

Visualisation

Tracé du champ vectoriel

display_vector_field lit un fichier de vecteurs sauvegardé et appelle plt.show() en interne, donc sélectionnez un backend non interactif pour les exécutions par lots.

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,   # même facteur utilisé dans scaling.uniform, pour mapper de nouveau sur l'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)

Visualisation personnalisée

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 de la vitesse", "viridis"),
    (axes[1], u, "Vitesse U", "RdBu_r"),
    (axes[2], v, "Vitesse 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)

Fonctions d'analyse

scripts/analyze.py regroupe ces fonctions autour d'un params.npz écrit par runner.py. Il déduit l'espacement spatial physique des coordonnées sauvegardées, donc les dérivées sont exprimées par unité de longueur :

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

Les formes autonomes, si vous préférez les calculer en ligne :

Vorticité

def compute_vorticity(u, v, dx=1.0, dy=None):
    """Vorticité hors plan dv/dx - du/dy. Passez l'espacement spatial physique, pas 1.0."""
    dy = dx if dy is None else dy
    return np.gradient(v, dx, axis=1) - np.gradient(u, dy, axis=0)

L'espacement spatial est (window_size - overlap) / scaling_factor en unités physiques, donc laisser dx=1.0 donne une vorticité par cellule de grille, pas par unité de longueur.

Convention de signe : runner.py se termine par transform_coordinates, qui réétiquette la grille en un repère y vers le haut droitier mais laisse les lignes dans l'ordre image, donc le y sauvegardé décroît à mesure que l'indice de ligne augmente. Les formes autonomes ci-dessus supposent l'inverse, donc sur un champ params.npz elles renvoient -du/dy et inversent le signe de la vorticité et du cisaillement — annulez les dérivées selon axis=0, ou utilisez PIVAnalyzer, qui lit l'orientation depuis les coordonnées sauvegardées.

Taux de déformation

def compute_strain(u, v, dx=1.0, dy=None):
    """Renvoie (exx, eyy, exy) du tenseur de taux de déformation 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)

Statistiques de turbulence

def compute_statistics(u, v):
    """Statistiques spatiales d'une seule trame. PAS une décomposition 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),
    }

Mise en garde : soustraire la moyenne spatiale d'une trame mesure la variance spatiale, qui est égale à l'intensité turbulente uniquement pour un champ homogène. Une véritable décomposition de Reynolds nécessite un ensemble de paires d'images : faites la moyenne sur l'axe temporel, puis soustrayez ce champ moyen à chaque réalisation.

Utilisation en ligne de commande

# Exécution de base
python skills/openpiv/scripts/runner.py \
    --image img1.bmp --image img2.bmp --output_dir results --verbose

# Paramètres ajustés avec masquage dynamique
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

Options de l'interface en ligne de commande

Option	Par défaut	Description
--image	requis	Fichier image ; spécifiez exactement deux fois pour la paire
--output_dir	results	Répertoire de sortie (créé s'il n'existe pas)
--window_size	32	Taille de la fenêtre d'interrogation (px)
--overlap	12	Chevauchement des fenêtres (px)
--search_area	38	Taille de la zone de recherche (px), doit être ≥ --window_size
--dt	0.02	Temps entre les images (s)
--scaling	96.52	Facteur d'échelle, pixels par unité physique (ex. px/mm)
--threshold	1.05	Seuil de rapport signal/bruit peak2peak
--mask	none	none ou dynamic (openpiv.preprocess.dynamic_masking)
--mask_method	intensity	edges ou intensity, utilisé lorsque --mask=dynamic
--verbose	False	Sortie détaillée

Script PIVAnalyzer

PIVAnalyzer dans scripts/analyze.py lit le fichier params.npz produit par runner.py et fournit des méthodes pour les analyses en aval. Il détecte automatiquement l'orientation de la grille à partir des coordonnées sauvegardées.

from analyze import PIVAnalyzer

piv = PIVAnalyzer("results/params.npz")

# Vorticité
vort = piv.compute_vorticity()

# Taux de déformation
exx, eyy, exy = piv.compute_strain()

# Statistiques (moyennes, RMS, énergie cinétique turbulente)
stats = piv.compute_statistics()
print(f"Énergie cinétique turbulente : {stats['tke']:.3f}")

# Tracé
display.plot_vector_field(save_path="quiver.png")

Scripts supplémentaires

Les scripts listés sous skills/openpiv/scripts/ incluent un exécuteur autonome et des wrappers d'analyse. Consultez le README de ce répertoire pour les options avancées.