openpiv

Was macht openpiv?

Übersicht

OpenPIV (Open Particle Image Velocimetry) analysiert Fluidströmungen aus PIV-Bildpaaren. Es umfasst Vorverarbeitung, Kreuzkorrelation, Vektorvalidierung, Ausreißerersetzung, Glättung und Skalierung auf physikalische Einheiten.

Alles unten ist gegen openpiv 0.25.4 verifiziert. Die API ändert sich zwischen den Versionen – überprüfen Sie inspect.signature(), bevor Sie einem Code-Snippet für eine andere Version vertrauen.

Wann verwenden

Verwenden Sie dieses Skill bei der Arbeit mit experimentellen PIV- oder Strömungsvisualisierungs-Bildpaaren: Messen von 2D-Geschwindigkeitsfeldern, Anpassen von Abfragefensterparametern, Validieren von Vektoren oder Ableiten von Vortizität, Dehnungsrate und Turbulenzstatistiken. Verwenden Sie für die Simulation von Strömungen anstelle der Messung ein CFD-Skill.

Schnellstart

OpenPIV installieren:

uv pip install openpiv

# Fixieren, wenn die Analyse reproduzierbar sein soll – dies ist die Version, gegen die jedes
# unten stehende Snippet getestet wurde.
uv pip install "openpiv==0.25.4"

PIV-Analyse an einem Bildpaar durchführen:

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

# Kreuzkorrelation. Gibt (u, v, s2n) zurück, wenn sig2noise_method nicht None ist.
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",   # erforderlich für 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 ist ein boolesches Array: True markiert einen fehlerhaften Vektor.
flags = validation.sig2noise_val(s2n, threshold=1.05)
u, v = filters.replace_outliers(u, v, flags, method="localmean", max_iter=3, kernel_size=2)

# Auf physikalische Einheiten skalieren, dann für die Darstellung in Bildkoordinaten umwandeln.
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)

Oder die gebündelte CLI verwenden, die genau diese Pipeline beinhaltet:

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

Kernkonzepte

PIV-Grundlagen

Particle Image Velocimetry ist eine optische Methode zur Messung der Fluidgeschwindigkeit durch Verfolgung beleuchteter Tracer-Partikel zwischen zwei Bildern.

Prozessablauf:

1. Ein Bildpaar aufnehmen (frame_a, frame_b), getrennt durch eine bekannte Zeit dt.
2. Die Bilder in Abfragefenster unterteilen.
3. Passende Fenster kreuzkorrelieren, um die Spitzenverschiebung zu finden.
4. Vektoren validieren (Signal-Rausch-Verhältnis, globaler Bereich, lokaler Median).
5. Fehlerhafte Vektoren durch interpolierte Werte ersetzen.
6. Pixelverschiebungen auf physikalische Einheiten skalieren.

Parameter der Abfragefenster

window_size — Korrelationsfenster in Pixeln (typischerweise 16–128). Größere Fenster ergeben eine bessere Korrelation, aber eine geringere räumliche Auflösung.

overlap — Pixel, die zwischen benachbarten Fenstern geteilt werden (typischerweise 50–75 % von window_size). Höherer Überlapp erhöht die Vektordichte und die Kosten, aber benachbarte Vektoren werden korreliert statt unabhängig.

search_area_size — Das im zweiten Bild durchsuchte Fenster. Muss ≥ window_size sein; ein paar Pixel mehr erlauben größere Verschiebungen. Paaren Sie einen erweiterten Suchbereich mit correlation_method="linear" – die Standardeinstellung "circular" basiert auf FFT-Wrap-around und aliasiert große Verschiebungen zu kleinen. Siehe references/advanced_algorithms.md.

Faustregeln: Halten Sie die größte Verschiebung unter etwa einem Viertel von window_size und zielen Sie auf 5–10 Partikel pro Fenster.

Signal-Rausch-Verhältnis

s2n misst, wie ausgeprägt die Korrelationsspitze ist. sig2noise_method steuert, wie es berechnet wird – "peak2mean" (der Funktionsstandard) oder "peak2peak". Die beiden befinden sich auf unterschiedlichen Skalen, daher ist ein für eine Methode abgestimmter Schwellenwert für die andere bedeutungslos. Typische peak2peak-Schwellenwerte sind 1,05–1,3.

flags = validation.sig2noise_val(s2n, threshold=1.05)
# flags ist bool: True == fehlerhaft. `~flags` wählt die guten Vektoren aus.

Häufige Operationen

Dynamische Maskierung

Die Maskierung befindet sich in openpiv.preprocess, nicht in einem openpiv.masking-Modul. Sie gibt ein (image, mask)-Tupel zurück und erwartet ein Float-Bild.

from openpiv import preprocess

# method="edges" für dunkle, scharfkantige Objekte; "intensity" für kontrastreiche Objekte.
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
)

Das zurückgegebene Bild in den Korrelationsschritt einspeisen – es hat den maskierten Bereich bereits auf Null gesetzt. Das ursprüngliche Bild nicht mit mask multiplizieren: Die Maskierung ist bereits angewendet, und bei method="edges" wird die Maske als uint8 0/255 statt bool zurückgegeben, sodass die Multiplikation das Bild um 255 neu skaliert.

Mehrfachdurchgangsverarbeitung

Mehrfachdurchgang (Fensterdeformation) befindet sich in openpiv.windef, gesteuert durch eine PIVSettings-Dataclass. pyprocess hat keinen Mehrfachdurchgangs-Einstiegspunkt.

import numpy as np
from openpiv import scaling, windef

settings = windef.PIVSettings()
settings.windowsizes = (64, 32, 16)   # ein Eintrag pro Durchgang, abnehmend (dies ist auch der Standard)
settings.overlap = (32, 16, 8)        # gleiche Länge wie windowsizes
settings.num_iterations = 3           # Anzahl der tatsächlich auszuführenden Durchgänge
settings.sig2noise_threshold = 1.05

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

# Ausgabe in PIXEL PRO FRAME -- selbst umrechnen. scaling.uniform dividiert nur
# durch scaling_factor, also dt separat anwenden.
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 validiert bereits, ersetzt Ausreißer, füllt verbleibende NaN mit Nullen und ruft transform_coordinates auf – wiederholen Sie diese Schritte nicht.

Einheitenfalle: PIVSettings hat dt- und scaling_factor-Felder, aber windef verwendet keines von beiden – first_pass ruft extended_search_area_piv ohne dt auf, sodass die gesamte Mehrfachdurchgangskette in Pixeln pro Frame arbeitet. Das Setzen von settings.dt = 0.02 ändert nichts an den zurückgegebenen Werten. Konvertieren Sie wie oben nach der Tatsache.

Für die Kontrolle über einzelne Durchgänge sind windef.first_pass und windef.multipass_img_deform die Bausteine auf niedrigerer Ebene.

Validierung und Nachbearbeitung

Validierungsmethoden

Jeder Validator gibt ein boolesches Array zurück, in dem True einen fehlerhaften Vektor markiert.

# Signal-Rausch-Verhältnis
flags = validation.sig2noise_val(s2n, threshold=1.05)

# Globaler Bereich -- nimmt (min, max) TUPEL, positionell oder als u_thresholds/v_thresholds.
flags = validation.global_val(u, v, (-300, 300), (-300, 300))

# Lokaler Median -- u_threshold und v_threshold sind ERFORDERLICH; size ist die Nachbarschaftshalbbreite.
flags = validation.local_median_val(u, v, u_threshold=30.0, v_threshold=30.0, size=1)

# Kombinieren mit booleschem ODER (nicht np.maximum -- dies sind boolesche 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)
)

Setzen Sie diese Schwellenwerte in den Einheiten von u und v, nicht in Pixeln pro Frame. extended_search_area_piv dividiert durch dt, sodass bei dt=0.02 eine Verschiebung von 3 px/Frame als 150 px/s ankommt. Die obigen Schwellenwerte passen für diesen Fall; die Angabe (-30, 30), die in der PIV-Literatur und PIVSettings.min_max_u_disp verwendet wird, ist ein px/Frame-Limit, und ihre Anwendung auf die px/s-Ausgabe verwirft das gesamte Feld. Validieren Sie entweder vor der Skalierung oder skalieren Sie die Schwellenwerte ebenfalls um 1/dt.

Ausreißerersetzung

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

method akzeptiert "localmean", "disk" oder "distance" – und nur diese drei. Ein nicht erkannter Name wird nicht abgelehnt; er fällt auf einen Null-Kernel zurück und gibt stillschweigend ein nutzloses Feld zurück. Beachten Sie, dass die Ersetzung die markierten Positionen mit interpolierten Werten füllt – wenn Sie sie dann mit NaN überschreiben, war die Ersetzung vergeblich. Wählen Sie das eine oder das andere:

# Markierte Vektoren vollständig aus der Analyse heraushalten, anstatt sie zu interpolieren.
u = np.where(flags, np.nan, u)
v = np.where(flags, np.nan, v)

Glättung

Die Glättung erfolgt mit openpiv.smoothn.smoothn; es gibt kein openpiv.smooth-Modul. Es gibt ein Tupel zurück, dessen erstes Element das geglättete Feld ist, und es akzeptiert keine NaN-Eingabe.

from openpiv.smoothn import smoothn

u_smooth, *_ = smoothn(np.nan_to_num(u), s=0.5)  # s: größer == glatter
v_smooth, *_ = smoothn(np.nan_to_num(v), s=0.5)
u_smooth = np.asarray(u_smooth)

Visualisierung

Vektorfelddiagramm

display_vector_field liest eine gespeicherte Vektordatei und ruft intern plt.show() auf, wählen Sie daher für Batch-Läufe ein nicht-interaktives Backend.

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,   # gleicher Faktor wie in scaling.uniform, um auf das Bild zurückzuführen
    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)

Benutzerdefinierte Visualisierung

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, "Geschwindigkeitsbetrag", "viridis"),
    (axes[1], u, "U-Geschwindigkeit", "RdBu_r"),
    (axes[2], v, "V-Geschwindigkeit", "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)

Analysefunktionen

scripts/analyze.py bündelt diese gegen eine von runner.py geschriebene params.npz. Es ermittelt den physikalischen Gitterabstand aus den gespeicherten Koordinaten, sodass die Ableitungen pro Längeneinheit ausgegeben werden:

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

Die eigenständigen Formen, falls Sie lieber inline rechnen möchten:

Vortizität

def compute_vorticity(u, v, dx=1.0, dy=None):
    """Out-of-plane-Vortizität dv/dx - du/dy. Übergeben Sie den physikalischen Gitterabstand, nicht 1.0."""
    dy = dx if dy is None else dy
    return np.gradient(v, dx, axis=1) - np.gradient(u, dy, axis=0)

Der Gitterabstand ist (window_size - overlap) / scaling_factor in physikalischen Einheiten, sodass das Belassen von dx=1.0 die Vortizität pro Gitterzelle ergibt, nicht pro Längeneinheit.

Vorzeichenkonvention: runner.py endet mit transform_coordinates, welches das Gitter in einen rechtshändigen y-nach-oben-Rahmen umbenennt, die Zeilen aber in der Bildreihenfolge belässt, sodass das gespeicherte y abnimmt, während der Zeilenindex wächst. Die obigen eigenständigen Formen nehmen das Gegenteil an, daher liefern sie auf einem params.npz-Feld -du/dy zurück und kehren das Vorzeichen der Vortizität und der Scherdehnung um – negieren Sie die axis=0-Ableitungen oder verwenden Sie PIVAnalyzer, der die Ausrichtung von den gespeicherten Koordinaten abliest.

Dehnungsrate

def compute_strain(u, v, dx=1.0, dy=None):
    """Gibt (exx, eyy, exy) des 2D-Dehnungsratentensors zurück."""
    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)

Turbulenzstatistiken

def compute_statistics(u, v):
    """Räumliche Statistik eines Einzelbildes. KEINE Reynolds-Zerlegung."""
    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),
    }

Einschränkung: Das Subtrahieren des räumlichen Mittelwerts eines Bildes misst die räumliche Varianz, die nur für ein homogenes Feld gleich der turbulenten Intensität ist. Eine echte Reynolds-Zerlegung erfordert eine Sammlung von Bildpaaren: Mitteln Sie über die Zeitachse und subtrahieren Sie dann dieses mittlere Feld von jeder Realisierung.

CLI-Verwendung

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

# Abgestimmte Parameter mit dynamischer Maskierung
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-Optionen

Option	Standard	Beschreibung
--image	erforderlich	Bilddatei; genau zweimal für das Paar angeben
--output_dir	results	Ausgabeverzeichnis (wird erstellt, falls nicht vorhanden)
--window_size	32	Abfragefenstergröße (px)
--overlap	12	Fensterüberlapp (px)
--search_area	38	Suchbereichsgröße (px), muss ≥ --window_size sein
--dt	0.02	Zeit zwischen den Bildern (s)
--scaling	96.52	Skalierungsfaktor, Pixel pro physikalischer Einheit (z.B. px/mm)
--threshold	1.05	peak2peak Signal-Rausch-Schwellenwert
--mask	none	none oder dynamic (openpiv.preprocess.dynamic_masking)
--mask_method	intensity	edges oder intensity, verwendet in …

… [für die Übersetzung abgeschnitten]