Suggestion for an alternative fast correlation method

[bschumac]

Currently working with the code for another project and saw that you have a function moving_window_array in the pyprocess.py which is unused.
Do you plan with it in the future?

[alexlib]

HI @bschumac - we do not use it apparently - do you know maybe where it can fit?

[bschumac]

There are other methods of correlating the windows which were used earlier in the PIV process. My suggestion is to implement the Greyscale Difference Method which is quite different to the cross correlation method. Futher info is found in
Kaga A. et al. (1992): Application of a fast algorithm for pattern tracking on airflow measurement. In: Proc 6th international symposium on flow visualization pp 853-857.

[alexlib]

There are other methods of correlating the windows which were used earlier in the PIV process. My suggestion is to implement the Greyscale Difference Method which is quite different to the cross correlation method. Futher info is found in
Kaga A. et al. (1992): Application of a fast algorithm for pattern tracking on airflow measurement. In: Proc 6th international symposium on flow visualization pp 853-857.

@bschumac This is a great idea. Feel free to fork the repository, add this subroutine, run tests on it and submit a pull request. if you need some help, send me a link to your code in your fork and I'll also take a look.

[bschumac]

There are other methods of correlating the windows which were used earlier in the PIV process. My suggestion is to implement the Greyscale Difference Method which is quite different to the cross correlation method. Futher info is found in
Kaga A. et al. (1992): Application of a fast algorithm for pattern tracking on airflow measurement. In: Proc 6th international symposium on flow visualization pp 853-857.

@bschumac This is a great idea. Feel free to fork the repository, add this subroutine, run tests on it and submit a pull request. if you need some help, send me a link to your code in your fork and I'll also take a look.

I can do that. I am using the code slightly different tho. Therefore it will need some time before I can integrate this method. Feel free to do it yourself if you have some time.
Cheers from NZ

[alexlib]

Any luck @bschumac ?

[ErichZimmer]

@bschumac
Is the process similar to minimum squared (quadratic) differences? If so, then it can be written as the following;

def mqd(image_a, image_b, correlation_method = 'circular'):
    '''
    FFT accelarated minimum quadratic differences correlation.
    MQD, at its simplest, can be calculated by corr[x,y] = (im1[x,y] - im2[x-u, y-v])^2.
    However, looping through the interoogation windows is expensive, so FFTs
    are used to speed up the process. Using FOIL, we can break down the equation
    to im1^2  - 2(im1*im2) + im2^2, making its application much easier.
    
    Parameters
    ----------
    image_a : 3d np.ndarray, first dimension is the number of windows,
        and two last dimensions are interrogation windows of the first image

    image_b : similar

    correlation_method : string
        one of the three methods implemented: 'circular' or 'linear'
        [default: 'circular].
    
    '''
    s1 = np.array(image_a.shape[-2:])
    s2 = np.array(image_b.shape[-2:])
    image_a = piv_prc.normalize_intensity(image_a)
    image_b = piv_prc.normalize_intensity(image_b)
        
    if correlation_method == 'circular':
        f2a = rfft2(image_a, axes = (-2, -1))
        f2b = rfft2(image_b, axes = (-2, -1))
        corr = fftshift((
            np.sum(image_a ** 2) # im1^2
            -2*irfft2(f2a.conj() * f2b).real # -2(im1* im2)
             + irfft2(f2b.conj() * f2b).real # im2^2; this was an accident that worked :D
            ), 
            axes = (-2, -1)
        )
        
    else:
        size = s1 + s2 - 1
        fsize = 2 ** np.ceil(np.log2(size)).astype(int)
        fslice = (
            slice(0, image_a.shape[0]),
            slice((fsize[0]-s1[0])//2, (fsize[0]+s1[0])//2),
            slice((fsize[1]-s1[1])//2, (fsize[1]+s1[1])//2)
        )
        f2a = rfft2(image_a, fsize, axes=(-2, -1))
        f2b = rfft2(image_b, fsize, axes=(-2, -1))        
        corr = fftshift(
            np.sum(image_a ** 2) - 2*\
            irfft2(f2a.conj() * f2b).real +\
            irfft2(f2b.conj() * f2b).real,
            axes = (-2, -1)
        )[fslice]
    corr *= -1 # invert so gaussian peak fit can work
    #corr /= (s2[0] * s2[1]) # poor attempt at normalization
    #corr -= corr.min(axis = 0)
    return corr

This method produces better signal-to-noise and peak-to-peak ratios than standard fft-based cross correlations. It appears to be slightly more robust while having the same computational speed as OpenPIV's vectorized approach. However, as I don't currently have access to that article, I am not sure if this is the algorithm you're mentioning.

[bschumac]

Wow amazing! I will need to doublecheck this but yes your equation is very similar to what I wrote in my TIV attempt (with loops and therefore way slower...)! I haven't been looking into the greyscale differences for ages, but I am very happy to check the speed and the similarities to Kaga's approach (which is just: im1 - im2 no quadratic difference involved).

I am currently finishing my PhD therefore this is on the ToDo list for afterwards (October 21).
I have uploaded the proceedings where you can find the article of Kaga under page pp 853-857.

https://drive.google.com/file/d/1nYmsS4wc0Uaiq-WdvrKYzqA5Tv7BmsUc/view?usp=sharing

[alexlib]

@bschumac
Is the process similar to minimum squared (quadratic) differences? If so, then it can be written as the following;

def mqd(image_a, image_b, correlation_method = 'circular'):
    '''
    FFT accelarated minimum quadratic differences correlation.
    MQD, at its simplest, can be calculated by corr[x,y] = (im1[x,y] - im2[x-u, y-v])^2.
    However, looping through the interoogation windows is expensive, so FFTs
    are used to speed up the process. Using FOIL, we can break down the equation
    to im1^2  - 2(im1*im2) + im2^2, making its application much easier.
    
    Parameters
    ----------
    image_a : 3d np.ndarray, first dimension is the number of windows,
        and two last dimensions are interrogation windows of the first image

    image_b : similar

    correlation_method : string
        one of the three methods implemented: 'circular' or 'linear'
        [default: 'circular].
    
    '''
    s1 = np.array(image_a.shape[-2:])
    s2 = np.array(image_b.shape[-2:])
    image_a = piv_prc.normalize_intensity(image_a)
    image_b = piv_prc.normalize_intensity(image_b)
        
    if correlation_method == 'circular':
        f2a = rfft2(image_a, axes = (-2, -1))
        f2b = rfft2(image_b, axes = (-2, -1))
        corr = fftshift((
            np.sum(image_a ** 2) # im1^2
            -2*irfft2(f2a.conj() * f2b).real # -2(im1* im2)
             + irfft2(f2b.conj() * f2b).real # im2^2; this was an accident that worked :D
            ), 
            axes = (-2, -1)
        )
        
    else:
        size = s1 + s2 - 1
        fsize = 2 ** np.ceil(np.log2(size)).astype(int)
        fslice = (
            slice(0, image_a.shape[0]),
            slice((fsize[0]-s1[0])//2, (fsize[0]+s1[0])//2),
            slice((fsize[1]-s1[1])//2, (fsize[1]+s1[1])//2)
        )
        f2a = rfft2(image_a, fsize, axes=(-2, -1))
        f2b = rfft2(image_b, fsize, axes=(-2, -1))        
        corr = fftshift(
            np.sum(image_a ** 2) - 2*\
            irfft2(f2a.conj() * f2b).real +\
            irfft2(f2b.conj() * f2b).real,
            axes = (-2, -1)
        )[fslice]
    corr *= -1 # invert so gaussian peak fit can work
    #corr /= (s2[0] * s2[1]) # poor attempt at normalization
    #corr -= corr.min(axis = 0)
    return corr

This method produces better signal-to-noise and peak-to-peak ratios than standard fft-based cross correlations. It appears to be slightly more robust while having the same computational speed as OpenPIV's vectorized approach. However, as I don't currently have access to that article, I am not sure if this is the algorithm you're mentioning.

@ErichZimmer - have you submitted this new method to openpiv-python? I cannot find it in the code. It looks great.

[ErichZimmer]

This has not been submitted at a pull request yet (mainly because I'm still testing it on my spare time). However, once I'm done, I'll create a pull request with single pass minimum quadratic differences (MQD) and gaussian transformed phase-only correlation (GTPC) algorithms.

[ErichZimmer]

Currently, the minimum quadratic differences algorithm doesn't appear to perform as well as some commercial softwares. In most cases, it performs just as good or slightly worse than the normal FFT-based cross corrections. Furthermore, since the autocorrelation isn't really calculated for image_a (to save complexity and speed), the correlation map could not normalized by the autocorrelation peak as it doesn't exist. However, in some cases, the aformentioned minimum quadratic differences algorithm can be more robust leading to fewer outliers. In conclusion, more testing needs to be done and possible new solutions should be explored before a final decision is made.
Note, when using zero-padding, mqd and phase correlation typically out performs OpenPIV's normalized cross-correlation algorithm. The above statements were for the much faster "circular" correlations.

[ErichZimmer]

Updated the minimum quadratic differences algorithm to include normalization to [0,1].

def minimum_quadradic_differences(
    image_a, image_b, 
    correlation_method = 'circular', 
):
    '''
    FFT accelarated minimum quadratic differences correlation.
    MQD, at its simplest, can be calculated by corr = (im1 - im2)^2.
    However, looping through the interogation windows is expensive, so FFTs
    are used to speed up the process. Using FOIL, we can break down the equation
    to im1^2  - 2(im1*im2) + im2^2, making its application much easier to apply.
    Correlation is normalized to [-1,1] by division of autocorrelation peak and
    then normalized to [0,1].

    Parameters
    ----------
    image_a : 3d np.ndarray, first dimension is the number of windows,
        and two last dimensions are interrogation windows of the first image

    image_b : similar

    correlation_method : string
        one of the two methods implemented: 'circular' or 'linear'
        [default: 'circular'].
    Returns
    -------
    corr : 3d np.ndarray
        a three dimensions array for the correlation function.
    '''
    if correlation_method not in ['circular', 'linear']:
        raise ValueError(f'Correlation method not supported {correlation_method}')
    s1 = np.array(image_a.shape[-2:])
    s2 = np.array(image_b.shape[-2:])
    
    if correlation_method == 'circular':
        f2a = rfft2(image_a, axes = (-2, -1))
        f2b = rfft2(image_b, axes = (-2, -1))
        f2sum = np.sum(image_a ** 2, axis = (-2, -1))
        corr = fftshift(
            f2sum[:, np.newaxis, np.newaxis] - 2*\
            irfft2(f2a.conj() * f2b) +\
            irfft2(f2b.conj() * f2b), 
            axes = (-2, -1)
        ).real
        corr = corr / np.nanmax(corr, axis = (-2, -1))[:, np.newaxis, np.newaxis]
    else:
        size = s1 + s2 - 1
        fsize = 2 ** np.ceil(np.log2(size)).astype(int)
        fslice = (
            slice(0, image_a.shape[0]),
            slice((fsize[0]-s1[0])//2, (fsize[0]+s1[0])//2),
            slice((fsize[1]-s1[1])//2, (fsize[1]+s1[1])//2)
        )
        f2a = rfft2(image_a, fsize, axes=(-2, -1))
        f2b = rfft2(image_b, fsize, axes=(-2, -1))        
        f2sum = np.sum(image_a ** 2, axis = (-2, -1))
        corr = fftshift(
            f2sum[:, np.newaxis, np.newaxis] - 2*\
            irfft2(f2a.conj() * f2b) +\
            irfft2(f2b.conj() * f2b), 
            axes = (-2, -1)
        ).real[fslice]
        corr = corr / np.nanmax(corr, axis = (-2, -1))[:, np.newaxis, np.newaxis]
    return (-corr + 1) / 2

It appears to not be as robust as other algorithms, but may provide more precise/accurate (have to work on that term) results. As I don't have any synthetic images with exact results at hand (tried to generate some, but failed), I will have to validate my claims later.

[alexlib]

Please add some test case and comparison with the existing methods

[ErichZimmer]

@alexlib ,
Here is a simple comparison I made earlier today.

RMS and bias error equations were used from the article:
http://www.rug.nl › Chapter_2PDF
University of Groningen The flapping flight of birds Thielicke, William

[alexlib]

Looks great @ErichZimmer can we see the rest please? Somehow it is not clear to me why some displacements are so much better for mqd while other correlations are insensitive.

[ErichZimmer]

@alexlib ,
There was a bug in the script which was caused by removing an addition sign. Here is the corrected results.

A full test would be done soon as this was a proof of concept idea of an automated test script.

[alexlib]

OK, feel free to create a pull request. I'm sure that in some cases MQD can be useful.

is displacement in pixels?
Could you try a larger displacements range please?