Image matching - Cylinder wake

Description

Image matching using skimage.transform.warp [source].

Setup

Packages

[1]:
%reload_ext autoreload
%autoreload 2
[2]:
import matplotlib.pyplot as plt
import numpy as np
import os

import pivuq

Load images

[3]:
parent_path = "./data/cylinder_wake/"
image_pair = np.array(
    [plt.imread(os.path.join(parent_path + ipath)).astype("float") for ipath in ["frameA.tif", "frameB.tif"]]
)

Plot raw images

[4]:
fig, axes = plt.subplots(ncols=3, sharex=True, sharey=True, figsize=(15, 5))

for i, ax in enumerate(axes[:2]):
    im = ax.imshow(image_pair[i], vmax=2**14)
    fig.colorbar(im, ax=ax)
    ax.set(title=f"$I_{i}$")

ax = axes[-1]
im = ax.imshow(np.abs(image_pair[1] - image_pair[0]), vmin=0, vmax=2**14)
fig.colorbar(im, ax=ax)
ax.set(title="$|I_1 - I_0|$")

for ax in axes:
    ax.set(xlim=(50, 400), ylim=(0, 600), xlabel="$j$", ylabel="$i$")
../_images/examples_02_image_matching_cylinder_8_0.png

Image Matching (PIV)

Load reference velocity

[5]:
data = np.loadtxt(os.path.join(parent_path + "vectors_PIV.dat"), skiprows=3).T

I, J = 75, 42
X = np.reshape(data[0], (I, J)) - 1  # zero-index
Y = np.reshape(data[1], (I, J)) - 1
U = np.stack((np.reshape(data[2], (I, J)), np.reshape(data[3], (I, J))))

Plot: reference fields

[6]:
fig, axes = plt.subplots(ncols=3, sharex=True, sharey=True, figsize=(15, 5))

for i, (ax, var) in enumerate(zip(axes[:2], ["U", "V"])):
    im = ax.pcolormesh(X, Y, U[i])
    fig.colorbar(im, ax=ax)
    ax.set(title=f"${var}$")

ax = axes[-1]
im = ax.pcolormesh(X, Y, np.linalg.norm(U, axis=0), vmax=20)
fig.colorbar(im, ax=ax)
ax.set(title="Magnitude")

for ax in axes:
    ax.set(xlim=(0, 450), ylim=(0, 800), xlabel="$j$", ylabel="$i$")
../_images/examples_02_image_matching_cylinder_13_0.png

Image matching

[7]:
%%time
warped_image_pair = pivuq.warp(
    image_pair,
    U,
    velocity_upsample_kind="linear",
    direction="center",
    nsteps=1,
    order=1,
)
CPU times: user 201 ms, sys: 21.8 ms, total: 222 ms
Wall time: 223 ms

Plot: warped image pairs

[8]:
fig, axes = plt.subplots(ncols=4, sharex=True, sharey=True, figsize=(20, 5))

for i, ax in enumerate(axes[:2]):
    im = ax.imshow(warped_image_pair[0], vmax=2**14)
    fig.colorbar(im, ax=ax)
    ax.set(title=f"$I_{i}$")

ax = axes[-2]
im = ax.imshow(np.abs(image_pair[1] - image_pair[0]), vmax=2**14)
fig.colorbar(im, ax=ax)
ax.set(title="$|I_1 - I_0|$")


ax = axes[-1]
im = ax.imshow(np.abs(warped_image_pair[1] - warped_image_pair[0]), vmax=2**14)
fig.colorbar(im, ax=ax)
ax.set(title=f"$|\hat{{I}}_1 - \hat{{I}}_0|$, sum = {np.abs(warped_image_pair[1] - warped_image_pair[0]).sum():g}")

for ax in axes:
    ax.set(xlim=(0, 450), ylim=(0, 800), xlabel="$j$", ylabel="$i$")
../_images/examples_02_image_matching_cylinder_17_0.png

Image Matching (Optical Flow: pixel-level)

Load reference velocity

[9]:
data = np.load(os.path.join(parent_path + "vectors_OF.npz"))

X = data["X"]  # already zero-indexed
Y = data["Y"]
U = data["U"]  # u, v

Plot: reference fields

[10]:
fig, axes = plt.subplots(ncols=3, sharex=True, sharey=True, figsize=(15, 5))

for i, (ax, var) in enumerate(zip(axes[:2], ["U", "V"])):
    im = ax.pcolormesh(X, Y, U[i])
    fig.colorbar(im, ax=ax)
    ax.set(title=f"${var}$")

ax = axes[-1]
im = ax.pcolormesh(X, Y, np.linalg.norm(U, axis=0), vmax=20)
fig.colorbar(im, ax=ax)
ax.set(title="Magnitude")

for ax in axes:
    ax.set(xlim=(0, 450), ylim=(0, 800), xlabel="$j$", ylabel="$i$")
../_images/examples_02_image_matching_cylinder_22_0.png

Image matching

[11]:
%%time
warped_image_pair = pivuq.warp(
    image_pair,
    U,
    velocity_upsample_kind="linear",
    direction="center",
    nsteps=1,
    order=1,
)
CPU times: user 57.8 ms, sys: 928 µs, total: 58.8 ms
Wall time: 58.5 ms

Plot: warped image pairs

[12]:
fig, axes = plt.subplots(ncols=4, sharex=True, sharey=True, figsize=(20, 5))

for i, ax in enumerate(axes[:2]):
    im = ax.imshow(warped_image_pair[0], vmax=2**14)
    fig.colorbar(im, ax=ax)
    ax.set(title=f"$I_{i}$")

ax = axes[-2]
im = ax.imshow(np.abs(image_pair[1] - image_pair[0]), vmax=2**14)
fig.colorbar(im, ax=ax)
ax.set(title="$|I_1 - I_0|$")


ax = axes[-1]
im = ax.imshow(np.abs(warped_image_pair[1] - warped_image_pair[0]), vmax=2**14)
fig.colorbar(im, ax=ax)
ax.set(title=f"$|\hat{{I}}_1 - \hat{{I}}_0|$, sum = {np.abs(warped_image_pair[1] - warped_image_pair[0]).sum():g}")

for ax in axes:
    ax.set(xlim=(0, 450), ylim=(0, 800), xlabel="$j$", ylabel="$i$")
../_images/examples_02_image_matching_cylinder_26_0.png