Phase correlation for added robustness

[ErichZimmer]

Is your feature request related to a problem? Please describe. Sometimes standard FFT correlations and normalized FFT correlations can made be more robust (ignoring zero padding) for lower quality homemade/non-commercial DPIV systems and in general noisy environments with non-uniform illumination. The use of filters mostly mitigate this issue, however a new cross correlation algorithm is proposed that is more resilient against said abnormalities and could produce better results without the additional computational costs of filters (still, pre-processing is recommended).

Describe the solution you'd like To increase the robustness of the cross correlations, phase filtering is applied (phase correlation). In the simplest form, phase correlation can be computed by calculating the inverse Fourier transform of the cross power spectrum. This requires a small modification to the existing OpenPIV algorithm which can be seen below.

def phase_correlation(aa, bb):
    aa = rfft2(aa).conj()
    bb = rfft2(bb)
    ab = aa * bb
    R = ab / np.absolute(ab) # calculation of cross power spectrum
    corr = fftshift(
        irfft2(
            R
        ).real, axes=(-2, -1)
    )
    return corr

However, discrepancies are caused by near zero values so filtering or windowing could be used to improve the robustness further. Furthermore, convolving a gaussian filter with the cross power spectrum seems to yield better results when there is a noisy background, suggesting that filtering in the Fourier domain is indeed useful.

Describe alternatives you've considered Don't do anything as OpenPIV is already good and robust.

Additional context I ran out of time so I'll have to finish this feature request later :(

Regards, Zimmer

[alexlib]

@ErichZimmer I favor accepting new correlation methods into openpiv-python. Thanks for your suggestions. Let's just make sure we can explain the user how to use one of those and when. Already now we have sometimes questions on what do with the linear vs circular or single pass vs multi pass or window deformation by @TKaeufer. In any case, the phase correlation was published to do great things for some cases and we have to show this - a simple decision tree for the user which correlation to choose and what are each advantage/disadvantage.

[ErichZimmer]

@alexlib Here is what I got so far for the phase correlation function.

def phase_correlation(
    image_a, image_b,
    correlation_method = 'circular',
    algorithm = 'spof',
    sigma = 1
):
    '''
    Phase filtering to produce a phase-only correlation. Two methods
    are implemented here: Gaussian filtered phase correlation and
    "symmetric" phase correlation. Which one is better, who knows.
    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].
    algorithm : string
        one of three filters implemented: 'phat', 'gphat' or 'spof'
    sigma: float
        used in gphat to widen the peak for better subpixel estimation
    '''
    if algorithm not in ['spof', 'phat', 'gphat']:
        raise ValueError('Algorithm not supported.')
    s1 = np.array(image_a.shape[-2:])
    s2 = np.array(image_b.shape[-2:])
    if correlation_method == 'linear':
        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)).conj()
        f2b = rfft2(image_b, fsize, axes=(-2, -1))
        r = f2a * f2b
        if algorithm == 'spof':
            r /= (np.sqrt(np.absolute(f2a) * np.absolute(f2b)) + 1e-10)
            corr = fftshift(irfft2(r).real, axes=(-2, -1))[fslice] / (s2[0] * s2[1])
            corr = np.clip(corr, 0, 1)
        elif algorithm == 'phat':
            r /= np.abs(r + 1e-10)
            corr = fftshift(irfft2(r).real, axes=(-2, -1))[fslice]
            corr = np.clip(corr, 0, 1)
        else:
            r /= np.abs(r + 1e-10)
            w1 = get_window(('gaussian', sigma), r[0,:,:].shape[0])
            w2 = get_window(('gaussian', sigma), r[0,:,:].shape[1]*2 -1)
            window = np.outer(w1, w2)
            corr = irfft2(r * rfft2(window), axes=(-2, -1))[fslice]
        return corr
    elif correlation_method == 'circular':
        f2a = rfft2(image_a).conj()
        f2b = rfft2(image_b)
        r = f2a * f2b
        if algorithm == 'spof':
            r /= (np.sqrt(np.absolute(f2a) * np.absolute(f2b)) + 1e-10)
            corr = fftshift(irfft2(r).real, axes=(-2, -1)) / (s2[0] * s2[1])
            corr = np.clip(corr, 0, 1)
        elif algorithm == 'phat':
            r /= np.abs(r + 1e-10)
            corr = fftshift(irfft2(r).real, axes=(-2, -1))
            corr = np.clip(corr, 0, 1)
        else:
            w1 = get_window(('gaussian', sigma), s2[0])
            w2 = get_window(('gaussian', sigma), s2[1])
            window = np.outer(w1, w2)
            r /= np.absolute(r+1e-10)
            corr = irfft2(r * rfft2(window), axes = (-2, -1))
        return corr
    else:
        print("Correlation method is not implemented.")

Sorry if there is some inconsistencies as I used a phone on this comment.

[alexlib]

What is spof and phat and what are the papers/books that describe these methods. What are the differences? Advantages/disadvantages? We need at least a couple of test cases to show the effect of these algorithms. Do you have some results? Even a separate repo is fine.

[alexlib]

And gphat of course :) @ErichZimmer thanks

[ErichZimmer]

I think we should stay with "symmetric" phase-only filtering (SPOF) as it is more robust than the standard phase-only filtering and its gaussian transformed version (mainly because of very poor implementation). It also produces a wider peak than the standard phase-only filter (at least during my little "play session") which seems to accept peak fitting functions pretty well. I'll be making a repository soon to demonstrate the robustness and when to/not to use phase correlation. References as enquired: Symmetric phase-only filtering: 10.1088/0957-0233/16/3/001 Gaussian transformed phase-only correlation: 10.1007/s00348-008-0492-6 Note, my implementation of "gphat" is not correct, which lead to my current findings. Robust phase correlation: 10.1088/0957-0233/20/5/055401 (An interesting derivative of phase correlation)

[ErichZimmer]

I think I managed to get everything working, and for an idealized image pair, phase correlation is slightly more accurate at times, but the original FFT cross-correlation seems to be slightly more stable. I still messed something up, so please ignore this comment :(

[alexlib]

Let's do a pull request with these tests included. We need it for the future extensions as well

[ErichZimmer]

Here is the corrected results. I threw JPIV in the mix for a random comparison. There seems to be a slight misalignment in the x-axis, but I'll have to fix that later.

For the pull requests, I have the new correlation functions split off in a different branch so there would be two pull requests, one for the correlation functions and one for the vectorized functions. I just need a stronger internet so I can push them to my fork, which might take some time. :(

[alexlib]

Please when it's ready, create a pull request to here https://github.com/OpenPIV/openpiv-python/pulls

[ErichZimmer]

I finished the testing. To solve the offset axis issue, I created 64 synthetic images and made image pairs by shifting the original image by a known distance. The displacement shift increases 0.0625 pixels per image pair. To get rid of any edge issues, all vectors two units from the edges were removed. All images were analyzed with multipass window deformation and windowing of 64x64>32x32>24x24 and 50% overlap. It's nice to see the errors go to zero or near zero at integer displacements, which verifies that I finally used the correct equations for bias and RMS errors.

[ErichZimmer]

With some more testing with "symmetric phase correlation", it performs about the same as the standard FFT cross correlation. I haven't done any extensive testing yet, but here are the preliminary results. For comparison, Fluere was added for additional comparisons and validations.

OpenPIV script settings: preprocessing: None windowing: None passes: 3 interrogation windows: 64>32>24 overlap: 50% deformation method: second image deformation algorithm: B-Splines deformation degree: 3 pass validation: local median with threshold of 3, none on last pass smoothing: 1.6 on first pass, 0.8 on second pass, none on last pass

Fluere: preprocessing: None windowing: None passes: 3 interrogation windows: 64>32>24 overlap: 50% deformation method: second image (Only option) deformation algorithm: sinc (7x7 kernel) deformation degree: N/A pass validation: local median with threshold of 3, none on last pass

[alexlib]

Let's add phase correlation as an additional option.

[ErichZimmer]

Should have a pull request for new filters and phase correlation. Here is the results for normalized phase correlation.

[alexlib]

Yes please. We need some tutorial for the user on how to choose one of the correlation types or what are pros and cons of each.

[ErichZimmer]

The "symmetric" phase correlation didn't show too much of a benefit over normalized correlations, and users might get confused unless we find where phase filtering shines. Until then, I don't plan on implementing it (unless I get the more advanced phase filtering functions to work).

[ErichZimmer]

Closing this issue for now as I am not persuing phase correlation methods.