Notebook

In [1]:

# %load test_windef.py
"""
Created on Fri Oct  4 14:33:21 2019

@author: Theo
"""

import matplotlib.pyplot as plt

import numpy as np
import openpiv.windef as windef
from test_process import create_pair, shift_u, shift_v, threshold
from openpiv.pyprocess import get_coordinates, extended_search_area_piv, get_field_shape
from openpiv import validation, filters
import scipy.ndimage as scn
# from openpiv.windef import frame_interpolation

from scipy.interpolate import RectBivariateSpline

In [2]:

frame_a, frame_b = create_pair(image_size=256)

# this test are created only to test the displacement evaluation of the
# function the validation methods are not tested here ant therefore
# are disabled.

In [3]:

# circular cross correlation
def test_first_pass_circ():
    """ test of the first pass """
    x, y, u, v, s2n = windef.first_pass(
        frame_a,
        frame_b,
        window_size=64,
        overlap=32,
        iterations=1,
        correlation_method="circular",
        subpixel_method="gaussian",
        do_sig2noise=True,
        sig2noise_method="peak2peak",
        sig2noise_mask=2,
    )
    # print("\n", x, y, u, v, s2n)
    plt.quiver(x,y,u,v)
    
    assert np.mean(np.abs(u - shift_u)) < threshold
    assert np.mean(np.abs(v - shift_v)) < threshold

In [4]:

test_first_pass_circ()

def multipass_img_deform( frame_a, frame_b, window_size, overlap, iterations, current_iteration, x_old, y_old, u_old, v_old, correlation_method="circular", subpixel_method="gaussian", do_sig2noise=False, sig2noise_method="peak2peak", sig2noise_mask=2, MinMaxU=(-100, 50), MinMaxV=(-50, 50), std_threshold=5, median_threshold=2, median_size=1, filter_method="localmean", max_filter_iteration=10, filter_kernel_size=2, interpolation_order=3, ): """ Multi pass of the PIV evaluation. This function does the PIV evaluation of the second and other passes. It returns the coordinates of the interrogation window centres, the displacement u, v for each interrogation window as well as the mask which indicates wether the displacement vector was interpolated or not. Parameters ---------- frame_a : 2d np.ndarray the first image frame_b : 2d np.ndarray the second image window_size : tuple of ints the size of the interrogation window overlap : tuple of ints the overlap of the interrogation window, e.g. window_size/2 x_old : 2d np.ndarray the x coordinates of the vector field of the previous pass y_old : 2d np.ndarray the y coordinates of the vector field of the previous pass u_old : 2d np.ndarray the u displacement of the vector field of the previous pass v_old : 2d np.ndarray the v displacement of the vector field of the previous pass subpixel_method: string the method used for the subpixel interpolation. one of the following methods to estimate subpixel location of the peak: 'centroid' [replaces default if correlation map is negative], 'gaussian' [default if correlation map is positive], 'parabolic' MinMaxU : two elements tuple sets the limits of the u displacment component Used for validation. MinMaxV : two elements tuple sets the limits of the v displacment component Used for validation. std_threshold : float sets the threshold for the std validation median_threshold : float sets the threshold for the median validation filter_method : string the method used to replace the non-valid vectors Methods: 'localmean', 'disk', 'distance', max_filter_iteration : int maximum of filter iterations to replace nans filter_kernel_size : int size of the kernel used for the filtering interpolation_order : int the order of the spline interpolation used for the image deformation Returns ------- x : 2d np.array array containg the x coordinates of the interrogation window centres y : 2d np.array array containg the y coordinates of the interrogation window centres u : 2d np.array array containing the u displacement for every interrogation window u : 2d np.array array containing the u displacement for every interrogation window mask : 2d np.array array containg the mask values (bool) which contains information if the vector was filtered """ x, y = get_coordinates(np.shape(frame_a), window_size, overlap) "calculate the y and y coordinates of the interrogation window centres" """The interpolation function dont like meshgrids as input. Hence, the edges must be extracted to provide the sufficient input. x_old and y_old are the coordinates of the old grid. x_int and y_int are the coordinates of the new grid""" import pdb # pdb.set_trace() print(f"Iteration {current_iteration}") y_old = y_old[:, 0] # y_old = y_old[::-1] x_old = x_old[0, :] print(x_old,y_old) y_int = y[:, 0] # y_int = y_int[::-1] x_int = x[0, :] print(x_int,y_int) # interpolating the displacements from the old grid onto the new grid # y befor x because of numpy works row major ip = RectBivariateSpline(y_old, x_old, u_old, kx=2, ky=2) u_pre = ip(y_int, x_int) ip2 = RectBivariateSpline(y_old, x_old, v_old, kx=2, ky=2) v_pre = ip2(y_int, x_int) # this one is doing the image deformation (see above) frame_b_deform = frame_interpolation( frame_b, x, y, u_pre, v_pre, interpolation_order=interpolation_order ) if do_sig2noise is True and \ current_iteration == iterations and \ iterations != 1: sig2noise_method = sig2noise_method else: sig2noise_method = None u, v, s2n = extended_search_area_piv( frame_a, frame_b_deform, window_size=window_size, overlap=overlap, search_area_size=window_size, width=sig2noise_mask, subpixel_method=subpixel_method, sig2noise_method=sig2noise_method ) shapes = np.array(get_field_shape(frame_a.shape, window_size, overlap)) u = u.reshape(shapes) v = v.reshape(shapes) s2n = s2n.reshape(shapes) # adding the recent displacment on to the displacment of the previous pass u += u_pre v -= v_pre # validation using gloabl limits and local median u, v, mask_g = validation.global_val(u, v, MinMaxU, MinMaxV) u, v, mask_s = validation.global_std(u, v, std_threshold=std_threshold) u, v, mask_m = validation.local_median_val( u, v, u_threshold=median_threshold, v_threshold=median_threshold, size=median_size, ) # adding masks to add the effect of alle the validations mask = mask_g + mask_m + mask_s # mask=np.zeros_like(u) # filter to replace the values that where marked by the validation if current_iteration != iterations: # filter to replace the values that where marked by the validation u, v = filters.replace_outliers( u, v, method=filter_method, max_iter=max_filter_iteration, kernel_size=filter_kernel_size, ) return x, y, u, v, s2n, mask

In [5]:

def test_multi_pass_circ():
    """ test fot the multipass """
    window_size = (128, 64, 32)
    overlap = (64, 32, 16)
    iterations = 3

    x, y, u, v, s2n = windef.first_pass(
        frame_a,
        frame_b,
        window_size[0],
        overlap[0],
        iterations,
        correlation_method="circular",
        subpixel_method="gaussian",
        do_sig2noise=True,
        sig2noise_method="peak2peak",
        sig2noise_mask=2,
    )
    u_old = u.copy()
    v_old = v.copy()
    plt.figure()
    plt.quiver(x,y,u,v,color='r')
    
    for i in range(2, iterations + 1):
        x, y, u, v, s2n, mask = windef.multipass_img_deform(
            frame_a,
            frame_b,
            window_size[i - 1],
            overlap[i - 1],
            iterations,
            i,
            x,
            y,
            u,
            v,
            correlation_method="circular",
            subpixel_method="gaussian",
            do_sig2noise=False,
            sig2noise_method="peak2peak",
            sig2noise_mask=2,
            MinMaxU=(-100, 50),
            MinMaxV=(-50, 50),
            std_threshold=1000000,
            median_threshold=200000,
            median_size=1,
            filter_method="localmean",
            max_filter_iteration=10,
            filter_kernel_size=2,
            interpolation_order=3,
        )
        plt.figure()
        plt.quiver(x,y,u,v,color='b')

    # print("\n", x, y, u, v, s2n)
    
    assert np.mean(np.abs(u - shift_u)) < threshold and np.any(u != u_old)
    assert np.mean(np.abs(v - shift_v)) < threshold and np.any(v != v_old)
    # the second condition is to check if the multipass is done.
    # It need's a little numerical inaccuracy.

In [6]:

test_multi_pass_circ()

In [7]:

# linear cross correlation
def test_first_pass_lin():
    """ test of the first pass """
    x, y, u, v, s2n = windef.first_pass(
        frame_a,
        frame_b,
        window_size=64,
        overlap=32,
        iterations=1,
        correlation_method="linear",
        subpixel_method="gaussian",
        do_sig2noise=True,
        sig2noise_method="peak2peak",
        sig2noise_mask=2,
    )
    print("\n", x, y, u, v, s2n)
    plt.quiver(x,y,u,v)
    assert np.mean(np.abs(u - shift_u)) < threshold
    assert np.mean(np.abs(v - shift_v)) < threshold


def test_multi_pass_lin():
    """ test fot the multipass """
    window_size = (128, 64, 32)
    overlap = (64, 32, 16)
    iterations = 3

    x, y, u, v, s2n = windef.first_pass(
        frame_a,
        frame_b,
        window_size[0],
        overlap[0],
        iterations,
        correlation_method="linear",
        subpixel_method="gaussian",
        do_sig2noise=True,
        sig2noise_method="peak2peak",
        sig2noise_mask=2,
    )
    u_old = u.copy()
    v_old = v.copy()
    i = 1
    for i in range(2, iterations + 1):
        x, y, u, v, sn, m = windef.multipass_img_deform(
            frame_a,
            frame_b,
            window_size[i - 1],
            overlap[i - 1],
            iterations,
            i,
            x,
            y,
            u,
            v,
            correlation_method="linear",
            subpixel_method="gaussian",
            do_sig2noise=False,
            sig2noise_method="peak2peak",
            sig2noise_mask=2,
            MinMaxU=(-100, 50),
            MinMaxV=(-50, 50),
            std_threshold=1000000,
            median_threshold=200000,
            median_size=1,
            filter_method="localmean",
            max_filter_iteration=10,
            filter_kernel_size=2,
            interpolation_order=3,
        )

    print("\n", x, y, u, v, s2n)
    plt.quiver(x,y,u,v)
    assert np.max(np.abs(u - shift_u)) < threshold and np.any(u != u_old)
    assert np.max(np.abs(v - shift_v)) < threshold and np.any(v != v_old)

    # the second condition is to check if the multipass is done.
    # It need's a little numerical inaccuracy.

In [8]:

test_first_pass_lin()

 [[ 31.5  63.5  95.5 127.5 159.5 191.5 223.5]
 [ 31.5  63.5  95.5 127.5 159.5 191.5 223.5]
 [ 31.5  63.5  95.5 127.5 159.5 191.5 223.5]
 [ 31.5  63.5  95.5 127.5 159.5 191.5 223.5]
 [ 31.5  63.5  95.5 127.5 159.5 191.5 223.5]
 [ 31.5  63.5  95.5 127.5 159.5 191.5 223.5]
 [ 31.5  63.5  95.5 127.5 159.5 191.5 223.5]] [[ 31.5  31.5  31.5  31.5  31.5  31.5  31.5]
 [ 63.5  63.5  63.5  63.5  63.5  63.5  63.5]
 [ 95.5  95.5  95.5  95.5  95.5  95.5  95.5]
 [127.5 127.5 127.5 127.5 127.5 127.5 127.5]
 [159.5 159.5 159.5 159.5 159.5 159.5 159.5]
 [191.5 191.5 191.5 191.5 191.5 191.5 191.5]
 [223.5 223.5 223.5 223.5 223.5 223.5 223.5]] [[-5.62150081 -5.63118854 -5.62711791 -5.61035926 -5.60842178 -5.64271644
  -5.6286937 ]
 [-5.61766025 -5.64184492 -5.62928086 -5.37469241 -5.61659489 -5.63686507
  -5.63101431]
 [-5.6254699  -5.65763585 -5.62191295 -5.61393824 -5.61914938 -5.60429813
  -5.64055421]
 [-5.62509928 -5.60967399 -5.59611104 -5.6095403  -5.64228957 -5.60900026
  -5.60216998]
 [-5.62244183 -5.62372048 -5.61031928 -5.61108526 -5.37292814 -5.64033133
  -5.39686542]
 [-5.38230552 -5.61752006 -5.62347095 -5.61228985 -5.61453877 -5.63714391
  -5.34771702]
 [-5.36775971 -5.6102013  -5.62860802 -5.61919778 -5.59542295 -5.61681
  -5.31502461]] [[3.24915689 3.28665092 3.26388789 3.32571495 3.27621682 3.25181558
  3.22616336]
 [3.24568837 3.29902151 3.29057884 3.23828797 3.2839164  3.23560294
  3.22271697]
 [3.26101256 3.25922452 3.26327762 3.26083888 3.26748495 3.24949335
  3.25512838]
 [3.29743577 3.27280519 3.27866063 3.26743418 3.2537556  3.244571
  3.25972992]
 [3.28027768 3.25960931 3.24293682 3.24000981 3.20112112 3.26138836
  3.26407557]
 [3.25283338 3.29824393 3.2690165  3.2359582  3.28863406 3.2807047
  3.26799498]
 [3.24410612 3.26092364 3.26917099 3.25335362 3.24524533 3.27838901
  3.24274377]] [[11.20926706 12.42393723 11.46690691 10.47058    12.76496112 11.29643274
  12.92448741]
 [11.95504236 11.52272746 11.74513278 11.25058796 13.85809624 11.49433842
  14.11509896]
 [10.96628904 12.35334914 12.71852201 10.85532175 11.42462637 12.0810409
  12.45498773]
 [11.76330316 11.30505963 11.39726804 12.69497315 12.52202166 12.20898187
  11.35453332]
 [11.75371364 12.72162247 12.87225487 12.37211671 12.89221981 12.71851731
  11.63151007]
 [12.53352641 13.05355922 13.79648244 12.97578834 11.66695102 14.67601484
  12.12166189]
 [10.88329178 12.26896229 12.06663699 12.84482765 13.29658085 11.83799976
  12.99743431]]