rectangular windows about openpiv-python HOT 14 CLOSED

Dear Alex,
Can you test these two functions? It should work right off the bat for rectangular windows of any shape and overlap.

def sliding_window_array(image, window_size = 64, overlap = 32):
    '''
    Basically, we have a 2d array and we want to perform cross-correlation
    over the interrogation windows. An approach could be to loop over the array
    but loops are expensive in python. So we create from the array a new array
    with three dimension, of size (n_windows, window_size, window_size), in
    which each slice, (along the first axis) is an interrogation window.
    '''
    if isinstance(window_size, tuple) == False and isinstance(window_size, list) == False:
        window_size = [window_size, window_size]
    if isinstance(overlap, tuple) == False and isinstance(overlap, list) == False:
        overlap = [overlap, overlap]

    x, y = get_rect_coordinates(image.shape, window_size, overlap, center_on_field = False)
    x = (x - window_size[1]//2).astype(int); y = (y - window_size[0]//2).astype(int)
    x, y = np.reshape(x, (-1,1,1)), np.reshape(y, (-1,1,1))

    win_x, win_y = np.meshgrid(np.arange(0, window_size[1]), np.arange(0, window_size[0]))
    win_x = win_x[np.newaxis,:,:] + x
    win_y = win_y[np.newaxis,:,:] + y
    windows = image[win_y, win_x]
    
    return windows

def get_rect_coordinates(frame_a, window_size, overlap, center_on_field = False):
    '''
    Rectangular grid wrapper of get_coordinates.
    '''
    if isinstance(window_size, tuple) == False and isinstance(window_size, list) == False:
        window_size = [window_size, window_size]
    if isinstance(overlap, tuple) == False and isinstance(overlap, list) == False:
        overlap = [overlap, overlap]
    _, y = get_coordinates(frame_a, window_size[0], overlap[0], center_on_field = False)
    x, _ = get_coordinates(frame_a, window_size[1], overlap[1], center_on_field = False)
    return np.meshgrid(x[0,:], y[:,0])

def get_coordinates(image_size, search_area_size, overlap, center_on_field = True):
    """Compute the x, y coordinates of the centers of the interrogation windows.
    the origin (0,0) is like in the image, top left corner
    positive x is an increasing column index from left to right
    positive y is increasing row index, from top to bottom


    Parameters
    ----------
    image_size: two elements tuple
        a two dimensional tuple for the pixel size of the image
        first element is number of rows, second element is
        the number of columns.

    search_area_size: int
        the size of the search area windows, sometimes it's equal to
        the interrogation window size in both frames A and B

    overlap: int = 0 (default is no overlap)
        the number of pixel by which two adjacent interrogation
        windows overlap.


    Returns
    -------
    x : 2d np.ndarray
        a two dimensional array containing the x coordinates of the
        interrogation window centers, in pixels.

    y : 2d np.ndarray
        a two dimensional array containing the y coordinates of the
        interrogation window centers, in pixels.

        Coordinate system 0,0 is at the top left corner, positive
        x to the right, positive y from top downwards, i.e.
        image coordinate system

    """

    # get shape of the resulting flow field
    field_shape = get_field_shape(image_size,
                                  search_area_size,
                                  overlap)

    # compute grid coordinates of the search area window centers
    # note the field_shape[1] (columns) for x
    x = (
        np.arange(field_shape[1]) * (search_area_size - overlap)
        + (search_area_size) / 2.0
    )
    # note the rows in field_shape[0]
    y = (
        np.arange(field_shape[0]) * (search_area_size - overlap)
        + (search_area_size) / 2.0
    )

    # moving coordinates further to the center, so that the points at the
    # extreme left/right or top/bottom
    # have the same distance to the window edges. For simplicity only integer
    # movements are allowed.
    if center_on_field == True:
        x += (
            image_size[1]
            - 1
            - ((field_shape[1] - 1) * (search_area_size - overlap) +
                (search_area_size - 1))
        ) // 2
        y += (
            image_size[0] - 1
            - ((field_shape[0] - 1) * (search_area_size - overlap) +
               (search_area_size - 1))
        ) // 2

        # the origin 0,0 is at top left
        # the units are pixels

    return np.meshgrid(x, y)

from openpiv-python.

Dear Alex,
Can you test these two functions? It should work right off the bat for rectangular windows of any shape and overlap.

def sliding_window_array(image, window_size = 64, overlap = 32):
    '''
    Basically, we have a 2d array and we want to perform cross-correlation
    over the interrogation windows. An approach could be to loop over the array
    but loops are expensive in python. So we create from the array a new array
    with three dimension, of size (n_windows, window_size, window_size), in
    which each slice, (along the first axis) is an interrogation window.
    '''
    if isinstance(window_size, tuple) == False and isinstance(window_size, list) == False:
        window_size = [window_size, window_size]
    if isinstance(overlap, tuple) == False and isinstance(overlap, list) == False:
        overlap = [overlap, overlap]

    x, y = get_rect_coordinates(image.shape, window_size, overlap, center_on_field = False)
    x = (x - window_size[1]//2).astype(int); y = (y - window_size[0]//2).astype(int)
    x, y = np.reshape(x, (-1,1,1)), np.reshape(y, (-1,1,1))

    win_x, win_y = np.meshgrid(np.arange(0, window_size[1]), np.arange(0, window_size[0]))
    win_x = win_x[np.newaxis,:,:] + x
    win_y = win_y[np.newaxis,:,:] + y
    windows = image[win_y, win_x]
    
    return windows

def get_rect_coordinates(frame_a, window_size, overlap, center_on_field = False):
    '''
    Rectangular grid wrapper of get_coordinates.
    '''
    if isinstance(window_size, tuple) == False and isinstance(window_size, list) == False:
        window_size = [window_size, window_size]
    if isinstance(overlap, tuple) == False and isinstance(overlap, list) == False:
        overlap = [overlap, overlap]
    _, y = get_coordinates(frame_a, window_size[0], overlap[0], center_on_field = False)
    x, _ = get_coordinates(frame_a, window_size[1], overlap[1], center_on_field = False)
    return np.meshgrid(x[0,:], y[:,0])

def get_coordinates(image_size, search_area_size, overlap, center_on_field = True):
    """Compute the x, y coordinates of the centers of the interrogation windows.
    the origin (0,0) is like in the image, top left corner
    positive x is an increasing column index from left to right
    positive y is increasing row index, from top to bottom


    Parameters
    ----------
    image_size: two elements tuple
        a two dimensional tuple for the pixel size of the image
        first element is number of rows, second element is
        the number of columns.

    search_area_size: int
        the size of the search area windows, sometimes it's equal to
        the interrogation window size in both frames A and B

    overlap: int = 0 (default is no overlap)
        the number of pixel by which two adjacent interrogation
        windows overlap.


    Returns
    -------
    x : 2d np.ndarray
        a two dimensional array containing the x coordinates of the
        interrogation window centers, in pixels.

    y : 2d np.ndarray
        a two dimensional array containing the y coordinates of the
        interrogation window centers, in pixels.

        Coordinate system 0,0 is at the top left corner, positive
        x to the right, positive y from top downwards, i.e.
        image coordinate system

    """

    # get shape of the resulting flow field
    field_shape = get_field_shape(image_size,
                                  search_area_size,
                                  overlap)

    # compute grid coordinates of the search area window centers
    # note the field_shape[1] (columns) for x
    x = (
        np.arange(field_shape[1]) * (search_area_size - overlap)
        + (search_area_size) / 2.0
    )
    # note the rows in field_shape[0]
    y = (
        np.arange(field_shape[0]) * (search_area_size - overlap)
        + (search_area_size) / 2.0
    )

    # moving coordinates further to the center, so that the points at the
    # extreme left/right or top/bottom
    # have the same distance to the window edges. For simplicity only integer
    # movements are allowed.
    if center_on_field == True:
        x += (
            image_size[1]
            - 1
            - ((field_shape[1] - 1) * (search_area_size - overlap) +
                (search_area_size - 1))
        ) // 2
        y += (
            image_size[0] - 1
            - ((field_shape[0] - 1) * (search_area_size - overlap) +
               (search_area_size - 1))
        ) // 2

        # the origin 0,0 is at top left
        # the units are pixels

    return np.meshgrid(x, y)

@ErichZimmer
please open a new pull request. For this function and the new notebook with the test.

please first create a clone or update (sync) a branch in your repo with the master, then add to this branch the new functions and features and submit a pull request. Thanks a lot for the great job.
it is very important that the tests and notebooks work - as someone might download a mismatched version.

from openpiv-python.

Additional use case
The measurement of boundary layer flow might also benefit from rectangular interrogation areas.

Additional example for the solution
JPIV also supports rectangular windows: https://eguvep.github.io/jpiv/settings_piv_interrogationwindow.html

Possible drawbacks
One drawback might be, that the performance of most FFT algorithms decrease, when the windows are not square and the length is not a power of two.

from openpiv-python.

@alexlib Is there a way we can incorporate rectangular windows through moving_window_array?

from openpiv-python.

from openpiv-python.

from openpiv-python.

I think I tried scikit-image view_as_windows function with some success before, but I'll try again soon. I think it produced a different array shape than numpy as_strided which caused some issues if I remember right.

from openpiv-python.

Okay, I got it to work, but it uses a lot of memory (more than NumPy's as_strided). However, it does seem to perform alright and seems to resolve some boundary layer issues. I'll play around with it and see what happens.
To get a rectangular grid for the x and y components, I simple performed get_coordinates and get_field_shape for each window size.

from openpiv-python.

Dear alex,
I created a very fast and low memory function for creating sub-windows from images.

def sliding_window_array(image, window_size = 64, overlap = 32):
    if isinstance(window_size, tuple) == False and isinstance(window_size, list) == False:
        window_size = [window_size, window_size]
    if isinstance(overlap, tuple) == False and isinstance(overlap, list) == False:
        overlap = [overlap, overlap]

    x, _ = get_coordinates(image.shape, window_size[0], overlap[0], center_on_field = False)
    _, y = get_coordinates(image.shape, window_size[1], overlap[1], center_on_field = False)
    x, y = np.meshgrid(x[0,:], y[:,0])
    
    x = (x - window//2).astype(int); y = (y - window//2).astype(int)
    x, y = np.reshape(x, (-1,1,1)), np.reshape(y, (-1,1,1))

    win_x, win_y = np.meshgrid(np.arange(0, window_size[0]), np.arange(0, window_size[1]))
    win_x = win_x[np.newaxis,:,:] + x
    win_y = win_y[np.newaxis,:,:] + y

    windows = image[win_y, win_x]
    return windows

So far, the results seem to be a copy of moving_window_array. Once I do further testing, I would recommend using this function due to its low memory usage.

Regards,
Erich

from openpiv-python.

I realized that I copy and pasted the wrong version of the function as this one has undefined variables. However, I'm not able to access GitHub with my computer right now (very slow internet in the single kilobyte/second range) so I'll have to edit it later. :(

from openpiv-python.

I realized that I copy and pasted the wrong version of the function as this one has undefined variables. However, I'm not able to access GitHub with my computer right now (very slow internet in the single kilobyte/second range) so I'll have to edit it later. :(

from openpiv-python.

I realized that I copy and pasted the wrong version of the function as this one has undefined variables. However, I'm not able to access GitHub with my computer right now (very slow internet in the single kilobyte/second range) so I'll have to edit it later. :(

from openpiv-python.

Here is a quick example I made.
https://github.com/OpenPIV/test_robustness/blob/main/notebooks/rectangular_windows_solution.ipynb

from openpiv-python.

FYI, the window/overlap size parameters are formatted to [y,x] so that functions like windef don't have to be modified and the axes match with the strided images. If we do add the above functions, I'll make sure to point this out in any future notebooks containing rectangular windows.

from openpiv-python.