What is the fundamental difference between linear and circular method of find correlations?

[Ivan Nepomnyashchikh]

Hello, dear OpenPIV community.

I have a rather simple python code which boils down to the following sequence:

1. find correlations with pyprocess.fft_correlate_images()
2. convert them to displacements with pyprocess.correlation_to_displacement()
3. find SNR with pyprocess.sig2noise_ratio()
4. validate SNR map with validation.sig2noise_val()
5. replace outliers with filters.replace_outliers()

I am shooting at about 92-95% of the valid vectors based on SNR. To achieve that, the maximum SNR threshold I could use was 1.003. Which is rather low.

When I changed the default settings of pyprocess.fft_correlate_images() from "correlation_method="circular" to

"correlation_method="linear", I managed to increase my SNR threshold to 2 which would consistently result in more than 95% of the valid vectors for every PIV pair in my dataset.

That being said, the velocity fields did not change much visually.

Obviously, the SNR threshold of 2 is a dream. But I do not understand why I got it by changing the correlation method from circular to linear. And I want to understand what I am doing.

I looked up the source code and respective numpy libraries. The algorithm is understandable to me in general. But I don't understand the mathematics behind it.

I am wondering if somebody can explain to me, please, what is the fundamental (i.e., mathematical or physical) difference between the linear and circular correlation methods, why the circular correlation method was picked to be the default one and why would I get much better SNR with the linear method.

Thank you in advance.

Ivan

[Alex Liberzon]

Hi Ivan 

I will split the answer into several replies. Please convert some of it to useful and practical examples we can add to our documentation or openpiv-python-examples

Answer no. 1: signal-to-noise ratio values: 

There is no single value that can be used to threshold invalid vectors. This is because different correlation methods provide different values of correlation maps (the result of cross-correlation as explained in the following answer) and different definitions: peak to peak or peak to mean. Therefore, changing the correlation method or the S/N ratio definition will create a new distribution of values for a given image pair or set. I suggest always running the small subset of images, creating a histogram of S2N values, and then choosing a threshold that separates less and more suspicious vectors. BTW, these vectors LOOK GOOD most of the time, so the high or low S2N is not yet a definition of an invalid vector. 

I recommend the book by Raffel. It contains everything one needs to know about PIV, at least in the context of OpenPIV. 

This quick source that I found is reasonably correct https://arxiv.org/pdf/1405.3023

[Alex Liberzon]

Answer no. 2

Circular vs Linear correlation. 

Cross-correlation basic definition we use is called 'direct cross-correlation,' and it is defined in the space domain. We take two finite signals, A and B, each of length N (let's say a single row of intensity values in the image) and shift B over A from left to right such that the first B touches that A from the left 1 pixel, so the shift is -N and then shifts and overlaps completely (zero shift) and to the right to the last pixel touches A - +N. At each shift, we multiply the intensity of overlapped regions between part of A and part of B and get a value of cross-correlation. The result is a signal of 2N-1 length. if we do it A with A we get correlation and obviously the highest value will be at the shift = 0, i.e. complete overlap. This animation is great, thanks to Zhigunov

https://commons.wikimedia.org/wiki/File:Cross_Correlation_Animation.gif

You can find this algorithm in OpenPIV here: 

    # elif correlation_method == 'direct':
    #     return convolve2d(normalize_intensity(window_a),
    #                       normalize_intensity(window_b[::-1, ::-1, ::-1]), 'full')
    else:
        raise ValueError("method is not implemented")

As you see, we eventually dropped it out because it is slower than the FFT-based ones. Its definition is rather clear in Scipy correlate

https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.correlate.html

Maybe we need to learn from Scipy how to choose automatically which method is faster. But the general agreement that FFT-based ones are faster. 

When the PIV community moved to FFT-based correlations, it becomes clear that the FFT based cross-correlation is a sort of "periodic" or "circular" assumption, i.e. the finite signal of length N is presumably periodic and will repeat itself again and again. Sort of imagine one interrogation window repeats itself over 32 pixels again and again. In this case the circular method, i.e. using the original signal of length N and creating a result of 2N-1 will give the correct result. The proof of this is out of the scope of this answer, I recommend the book, but you can find it here also: 

http://www.learnpiv.org/fft/

So, to adjust the FFT-based correlation algorithm to the finite length non-periodic interrogation windows, the signal processing community suggests to pad zeros of at least length N on each side of the window - this means that the basic "period" is much longer and the inner (middle) part is our signal and then the FFT-based will correctly represent the shift close to the center of the window, less than 1/2 of window  size and will tame down the high frequencies coming from the edges of the window (due to non-periodicity of these finite windows). So linear one is the same formula as the circular one, but with zero-padded windows

Erich Zimmer asked to add both algorithms to the OpenPIV - he claimed that in some cases circular is faster (smaller windows) and if the shift is not too large, then maybe circular can be even more accurate. In the original OpenPIV (in Matlab and Python) we had only zero-padded version - as Raffel's book recommends. 

You can find the discussion here: 

https://github.com/OpenPIV/openpiv-python/pull/208

If you could - please convert this answer to a good tutorial and we'll add it to the documentation and examples. 

Thanks

Alex

On Saturday, July 6, 2024 at 2:20:22 AM UTC+3 nepo...@oregonstate.edu wrote:

[Alex Liberzon]

Actually ChatGPT4o gave me a reasonable answer and even added a small OpenPIV example :) 

Zero-padding is a common technique used before performing FFT-based cross-correlation to improve the accuracy and interpretability of the results. Here's an explanation of why zero-padding is important:

### Purpose of Zero-Padding

1. **Increase Resolution in the Frequency Domain**:
   - **Interpolation**: Zero-padding the signal before taking the FFT increases the number of points in the FFT result, effectively interpolating the frequency spectrum. This helps achieve finer resolution in the frequency domain, making it easier to identify the peak correlation more accurately.
   - **More Accurate Peak Detection**: With higher resolution, the position of the peak in the cross-correlation function can be determined more precisely, which is particularly important for accurate time delay estimation.

2. **Avoid Circular Convolution Artifacts**:
   - **Linear Convolution**: In FFT-based cross-correlation, the use of FFT inherently assumes circular convolution. Zero-padding helps to mitigate the circular convolution artifacts, making the result more representative of the true linear convolution. This is because zero-padding effectively separates the instances of the signals, avoiding overlap that can cause spurious correlations.
   - **Signal Separation**: By padding with zeros, the effective length of the signals is increased, ensuring that the end of one signal does not interfere with the start of the other, which can distort the cross-correlation result.

3. **Handling Boundary Effects**:
   - **Reducing Edge Effects**: Zero-padding helps to reduce edge effects by adding extra space around the signal. This ensures that the FFT does not wrap around the edges of the signal, which can introduce artifacts and distortions in the cross-correlation computation.

### Demonstration with OpenPIV-Python

Let's demonstrate the effect of zero-padding in FFT-based cross-correlation using OpenPIV-Python. We'll perform cross-correlation on a pair of signals with and without zero-padding.

```python
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import correlate
import openpiv

# Create example signals
N = 100  # Number of samples
t = np.linspace(0, 1, N, endpoint=False)
f = np.sin(2 * np.pi * t)  # Sine wave
g = np.sin(2 * np.pi * t + np.pi / 4)  # Phase-shifted sine wave

# Compute linear cross-correlation without zero-padding
linear_corr = correlate(f, g, mode='full')

# Zero-pad the signals
padding_length = N
f_padded = np.pad(f, (0, padding_length), 'constant')
g_padded = np.pad(g, (0, padding_length), 'constant')

# Compute FFT-based cross-correlation with zero-padding
f_fft = np.fft.fft(f_padded)
g_fft = np.fft.fft(g_padded)
circular_corr = np.fft.ifft(f_fft * np.conj(g_fft)).real

# Generate lag arrays
lags_linear = np.arange(-N + 1, N)
lags_circular = np.arange(-N + 1 + padding_length, padding_length)

# Plot the results
plt.figure(figsize=(12, 8))

plt.subplot(2, 1, 1)
plt.plot(lags_linear, linear_corr)
plt.title('Linear Cross-Correlation without Zero-Padding')
plt.xlabel('Lag')
plt.ylabel('Correlation')
plt.grid()

plt.subplot(2, 1, 2)
plt.plot(lags_circular, circular_corr)
plt.title('Circular Cross-Correlation with Zero-Padding')
plt.xlabel('Lag')
plt.ylabel('Correlation')
plt.grid()

plt.tight_layout()
plt.show()
```

### Observations
- **Without Zero-Padding**: The linear cross-correlation without zero-padding might have lower resolution and could suffer from edge effects, making it harder to accurately determine the peak correlation.
- **With Zero-Padding**: The FFT-based cross-correlation with zero-padding provides higher resolution and mitigates circular convolution artifacts, leading to more accurate and interpretable results.

### Conclusion
Zero-padding before FFT-based cross-correlation is crucial for achieving higher resolution in the frequency domain, avoiding artifacts from circular convolution, and reducing boundary effects. This practice ensures more accurate and reliable cross-correlation results, which is particularly important in applications requiring precise time delay estimation and signal analysis.

[Ivan Nepomnyashchikh]

Thank you very much, Dr. Liberzon!

This explains a lot! I will think more about all this. I can try and create an example as you are suggesting, but later. Right now I am learning how to create Jupiter notebooks and make examples in them. I have finally (after a year of on and off work) completed the VortexFitting function for the PIVPY package and now I am trying to make a Jupiter Notebook example for it. Once I am done with that, I will make a pull request for my Gamma1 and Gamma2 implementations for PIVPY. And after that I will do this comparison between linear and circular correlations.

Ivan

[Ivan Nepomnyashchikh]

Hello, Dr. Liberzon.

I converted this discussion into a Jupiter Notebook. It is ready for a pull request. But I have already submitted another pull-request for another notebook which has been waiting your approval. I am wondering if you can accept that pull request, please, and I will submit this one. Also, the pull request for the median filter in windef.py is still waiting your approval as well.

Thank you.

Ivan