Cross-correlation using the FFT

1,880 views
Skip to first unread message

Theo Käufer

unread,
Aug 21, 2019, 10:46:42 PM8/21/19
to PIVlab

I have spent some time analysing the code of the FFTmulti function and it is really impressive. But I got one question concerning the actual cross-correlation in the frequency domain.

 

The cross-correlation is done in the frequency domain by this line of code:

 

result_conv = fftshift(fftshift(real(ifft2(conj(fft2(image1_cut)).*fft2(image2_cut))), 1), 2);

 

This line of code is based on the cross-correlation theorem which states that the cross correlation of two functions in the time or spatial domain can be obtained by multiplying the functions in the frequency domain.

For discrete signals this line of code delivers the circular cross-correlation of the signals. The circular cross correlation is in general not the same as the linear cross-correlation which is normally use to determine the particle displacement.

But the linear cross-correlation can be calculated using the equation of the circular cross-correlation if the signals are zero-padded to the size S=size(signal1)+size(signal2)-1. This is used in Fast Correlation algorithms.

 

I went through the code several times and so far, I could not see any zero-padding of the image1_cut or image2_cut. Now I wonder if PIVlab is using linear cross-correlation to determine the particle displacement or if I got something wrong?

If some more information are needed i will be glad to provide them but I want to keep it simple at the first step.


Thank you in advance


William Thielicke

unread,
Aug 22, 2019, 9:14:33 AM8/22/19
to PIVlab
Hi Theo, thanks a lot for your comment. I must admit that I wrote this line of code about 10 years ago, and I don't remember why I did it like this... There is no zero-padding, you are right. I looked up this subject on the web and found this:
https://dsp.stackexchange.com/questions/736/how-do-i-implement-cross-correlation-to-prove-two-audio-files-are-similar

I guess you already had a look at it. However, if I follow the "instructions" given there, and zero pad the sub-images before doing the correlation, then the result is pretty bad... Have you tried it yourself?

Have you checked my chapter on PIVLab here?
http://www.rug.nl/research/portal/files/14094707/Chapter_2.pdf
On one page (I can't find it now, because the internet in my train is super slow), I am showing pictures comparing the cross-correlation matrix of DFT vs DCC. You can see that the correlation peak flips around if the displacement is too large. Does that add any valuable info to this discussion...? Could you check this and also modify the correlation so that it gives a good 'non-circular' correlation?

Just to make sure that nobody gets scared: All the algorithms were extensively tested and perform at least as good as commercial piv software.

Theo Käufer

unread,
Aug 22, 2019, 12:55:46 PM8/22/19
to PIVlab
Hi William, thank you for the fast response and the helpful links.

I modified the code so that the cross-correlation multiplication is performed with zero-padded images.

image1_cut=padarray(image1_cut,[interrogationarea-1 interrogationarea-1],0,'post'); %zero-pad the interrogation areas to size S= 2*interrogationarea-1
image2_cut=padarray(image2_cut,[interrogationarea-1 interrogationarea-1],0,'post'); %zero-pad the interrogation areas to size S= 2*interrogationarea-1

result_conv = fftshift(fftshift(real(ifft2(conj(fft2(image1_cut)).*fft2(image2_cut))), 1), 2); % do the  cross-correlation in the frequency domain
result_conv =result_conv((interrogationarea/2)+1:(3*interrogationarea/2),(interrogationarea/2)+1:(3*interrogationarea/2),:);

I cut the result to the size of an interrogation window (eg 64x64). I had to do this to use the rest of the code in the FFTmulti function since the original
output had the size of an interrogation window. This should make no difference for displacements smaller than interrogation window/2. I
modified the SubPixOfffset to 0 when I cut the results to obtain the matrix entries I needed. Alternative:
result_conv =result_conv((interrogationarea/2):(3*interrogationarea/2)-1,(interrogationarea/2):(3*interrogationarea/2)-1,:);  with SubPixOff= 1 gives the same result.
In fact, the zero padded multiplication in the frequency domain with cut results should give the same result as
convr=conv2(image1_cut,rot90(image2_cut,2),'same'),(correlation is equal to convolution with reversed signal) I tried it on some 64x64 images and it does.

In the next step I plotted the vector field of the displacement obtained by the cross-correlation with and without zero padding with the exact same settings. The results look similar but not the same. I attached them to this post.
I also plotted a surface plot of the visualize the difference of theabsolute value of the displacement depending on the image position. I also attached this plot to the post. The maximum displacement is about 6 px/frame and the z-scaling is in px/frame just for the context. The flow analysed is a jet coming entering from the right side of the image. If you want I can send you the image I used for evaluation.

You wrote: On one page (I can't find it now, because the internet in my train is super slow), I am showing pictures comparing the cross-correlation matrix of DFT vs DCC. You can see that the correlation peak flips around if the displacement is too large. Does that add any valuable info to this discussion...? Could you check this and also modify the correlation so that it gives a good 'non-circular' correlation?

I checked this and I think this should make no difference as long as your displacement is smaller then interrogation window/2.

For my investigations i assumed that i can use the rest of the code of the function without further modifying apart from the zerpadding and cutting the result and the zero-padding offset. So far i have not noticed any reason why it should not work.
Please correct me when I am wrong with that assumption or i made a mistake with the zero-padding.

Best regards

Theo

vectorfield_no_zeropad.png
vectorfield_zeropad.png
Deviation_betwenn_corr_methods.png

William Thielicke

unread,
Aug 22, 2019, 1:19:57 PM8/22/19
to PIVlab
Hi, I think we need to compare the correlation matrices and not the vector map in order to learn the (potential) difference of what you are suggesting. Because there are many things going on in piv_fftmulti that affect the final output. For example, the same correlation function is performed multiple times in piv_fftmulti, and you would need to modify the code in every instance.
I think your zero-padding is not as intended (set a breakpoint in piv_fftmulti and display image1_cut).
Zero-padding should be like this in my opinion:
padarray(image1_cut, [interrogationarea/2 interrogationarea/2],0,'both')

I personally have the feeling that there will not be a difference as long as the displacement is smaller than interrogationarea/2.
And that will always be the case when you use more than one pass.

Theo Käufer

unread,
Aug 22, 2019, 2:59:22 PM8/22/19
to PIV...@googlegroups.com
Hi, again thank you for the fast respsonse. I forgot to mention that I did just one single pass for testing in order to keep things as simple as possible.

Referring to this page of the matlab documentation and based on my (little) personal experience I think that the zero padding has to be done at the end of the signal. I achieve the right results this way.

For further simplification I copied result_conv = fftshift(fftshift(real(ifft2(conj(fft2(A)).*fft2(B))), 1), 2);

and used it to calculate the cross correlation of 2 synthetic particle images A and B with the size of 32x32 px. I compared the resulting correlation matrix with

the result of the zero-padded version of this multiplication result_conv = fftshift(fftshift(real(ifft2(conj(fft2(A_pad)).*fft2(B_pad))), 1), 2)

and the cross-correlation of B and A in the spatial domain. For comparison I used imagesc(*correlation map*).

 

The visualization of the correlation matrix shows small differences between the multiplication used in FFTmulti and the zero-padded version of it with cutted correlation matrix

The code I used is below. I attached the Images I used so you can try yourself if you want to.


%read example images
A=im2double(imread('Image1.png'));
B=im2double(imread('Image2.png'));

%A and B have the same size

%pad the A and B to size 2*A-1 and 2*B-1
A_pad=padarray(A,[size(A,1)-1 size(A,2)-1],0,'post');
B_pad=padarray(B,[size(B,1)-1 size(B,2)-1],0,'post');

%%%%%%%%%%%%%%% Cross-correlation using FFT with padded images %%%%%%%%%%%%%
result_conv_zeropad = fftshift(fftshift(real(ifft2(conj(fft2(A_pad)).*fft2(B_pad))),1),2);

%%%%%%%%%%%%%%% Cross-correlation using the defaul FFTmulti Cross-correlating function%%%%%%%%%%%%%
result_conv_no_zeropad = fftshift(fftshift(real(ifft2(conj(fft2(A)).*fft2(B))),1),2);

%%%%% Cross correlation in the spatial domain using the MATLAB xcorr2 function
result_xcorr2=xcorr2(B,A);

%%%%%%%% plots of the correlation maps %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1)
title('Correlation maps based correlation matrices')
subplot(1,3,1)
%%%%%%%%%% with zero padding%%%%%%%%%%%%%%%%%%%%%
imagesc(result_conv_zeropad)
title('fft cross-correlation with zeropadding')

daspect([1 1 1 ])

subplot(1,3,2)
%%%%%%%%%% without zero padding%%%%%%%%%%%%%%%%%%%%%
imagesc(result_conv_no_zeropad)
title('fft cross-correlation without zeropadding')
daspect([1 1 1 ])

subplot(1,3,3)
%%%%%%%%%%% in the spatial domain %%%%%%%%%%%%%%%%%
imagesc(result_xcorr2)
title('spatial cross-correlation')
daspect([1 1 1 ])


%%%%%%%%%%%%%%%% figures with cutted correlation matrices%%%%%%%%%%%%%%%%

figure(2)
title('Correlation maps based on cutted correlation matrices')
subplot(2,2,1)
%%%%%%%%%% with zero padding%%%%%%%%%%%%%%%%%%%%%
imagesc(result_conv_zeropad(size(A,1)/2+1:3*size(A,1)/2,size(A,2)/2+1:3*size(A,2)/2))
title('fft cross-correlation with zeropadding')

daspect([1 1 1 ])

subplot(2,2,2)
%%%%%%%%%% without zero padding%%%%%%%%%%%%%%%%%%%%%
imagesc(result_conv_no_zeropad)
title('fft cross-correlation without zeropadding')
daspect([1 1 1 ])

subplot(2,2,3)
%%%%%%%%%%% in the spatial domain %%%%%%%%%%%%%%%%%
imagesc(result_xcorr2(size(A,1)/2+1:3*size(A,1)/2,size(A,2)/2+1:3*size(A,2)/2))
title('spatial cross-correlation')
daspect([1 1 1 ])

%%%%%%%%%%% with convolution and flipped signal using%%%%%%%%%%%%%%%%%
result_conv2=conv2(B,rot90(A,2),'same');
subplot(2,2,4)
imagesc(result_conv2)
title('spatial cross-correlation by convolution with a flipped signal')
daspect([1 1 1 ])

Concerning the size of the correlation matrix I agree with you that it will not make any difference as long as the displacement is small enough. But in my opinion, this has no impact on the difference between linear and circular cross-correlation.


Edit: I forgot to attached the plotted correlation maps.

Image1.png
Image2.png
correlation_matrices.png
cutted_correlation_matrices.png

William Thielicke

unread,
Aug 23, 2019, 5:21:50 AM8/23/19
to PIVlab
Dear Theo, thanks a lot for your input. I will try this soon. I will check the influence on the correlation matrix for different displacements, and then I will perform some tests on the accuracy of the modified correlation with lots of synthetic images. That may take some days, but I will post the progress here!

William Thielicke

unread,
Aug 27, 2019, 1:54:01 PM8/27/19
to PIVlab
Hi Theo, I am just about to start analyzing the effect of the modification, but this is one of the first results, and I wonder if I did something wrong....:
(PIVlab is set to run 1 pass only)

Original Code:
image1_cut = image1_roi(ss1);
image2_cut = image2_roi(ss1);

%do fft2:
result_conv = fftshift(fftshift(real(ifft2(conj(fft2(image1_cut)).*fft2(image2_cut))), 1), 2);

Padding code:
image1_cut = image1_roi(ss1);
image2_cut = image2_roi(ss1);

% padding (faster than als padarray):
image1_cut=[image1_cut zeros(interrogationarea,interrogationarea-1,size(image1_cut,3)); zeros(interrogationarea-1,2*interrogationarea-1,size(image1_cut,3))];
image2_cut=[image2_cut zeros(interrogationarea,interrogationarea-1,size(image2_cut,3)); zeros(interrogationarea-1,2*interrogationarea-1,size(image2_cut,3))];

%do fft2:
result_conv = fftshift(fftshift(real(ifft2(conj(fft2(image1_cut)).*fft2(image2_cut))), 1), 2);

%cropping of correlation matrix:
result_conv =result_conv((interrogationarea/2):(3*interrogationarea/2)-1,(interrogationarea/2):(3*interrogationarea/2)-1,:); 

I am using the PIVlab_Oseen_z_01/02 example images, and set the interrogationa area to 8 px and step to 8 px. The results are attached. It seems like the padded code gives less valid vectors. Any idea what is going on? Something wrong in the code? The problem is not happening when e.g. interrogation area is set to 64 and step to 8.

I suspect that something in these lines might go wrong, maybe you are willing to assist to find the problem:

minres = permute(repmat(squeeze(min(min(result_conv))), [1, size(result_conv, 1), size(result_conv, 2)]), [2 3 1]);
deltares = permute(repmat(squeeze(max(max(result_conv))-min(min(result_conv))),[ 1, size(result_conv, 1), size(result_conv, 2)]), [2 3 1]);
result_conv = ((result_conv-minres)./deltares)*255;

or

% we need only one peak from each couple pictures
[z1, zi] = sort(z);
dz1 = [z1(1); diff(z1)];
i0 = find(dz1~=0);
x1 = x(zi(i0));
y1 = y(zi(i0));
z1 = z(zi(i0));

or...?

padding_002.jpg
no_padding_001.jpg

William Thielicke

unread,
Aug 27, 2019, 2:45:21 PM8/27/19
to PIVlab
I read this (Fundamentals of digital particle
image velocimetry
J Westerweel 1997), it sounds like they think zero-padding is not necessary:

One should be aware of the fact that the DFT applies
to periodic signals. The correlation domain for (43) ranges
from −N + 1 to N for each component. Hence, the DFT
should be carried out on a 2N × 2N domain. This can
be accomplished by padding I 0 and I 00 with zeros. When
the DFT is carried out over a smaller domain, then the
part of the correlation that is not resolved is folded back
onto the correlation. This is comparable to ‘aliasing’ in the
frequency domain.
When the displacement complies with the PIV design
rule in (33), i.e. |sD| < 1
4DI , then the correlation is zero for
|r|, |s| > 1
4N. In that case, zero padding is not required,
and the DFT can be carried out on an N × N domain.
However, one should be cautious when evaluating multiple-
exposure recordings.

Theo Käufer

unread,
Aug 27, 2019, 5:38:14 PM8/27/19
to PIV...@googlegroups.com

Hi William, thank you for the citation of the article.

 

I further searched for this topic again in this book: Particle image velocimetry: A practical guide (Raffel et al. 2018)

I found some information in section 5.3.1.1

 

They state that zero-padding can even deteriorate the correlation signal.

 

 

"One of these methods, zero padding, which entails extending the sample size to
four times the original size by filling in zeroes, will perform poorly because the data
(i.e. image sample) generally consists of a nonzero (noisy) background on which the
signal (i.e. particle images) is overlaid. The edge discontinuity brought about in the
zero padding process contaminates the spectra of the data with high frequency noise
which in turn deteriorates the cross-correlation signal."

 

The refer to the method you used in the FFTmulti function as circular correlation.

 

"If all of the above points are properly handled, an FFT-based interrogation algorithm
as shown in Fig. 5.10 will reliably provide the necessary correlation data from
which the displacement data can be retrieved. For the reasons given above, this
implementation of the cross-correlation function is sometimes referred to as circular
cross-correlation compared to the linear cross-correlation of Eq. (5.1)."

 

That answers my original question if PIVlab uses linear correlation to determine the displacement.

 

If you still want to further investigate the difference in the results between circular and linear correlation let me know. I would be glad to help you!

 

Otherwise, thank you for the discussion and the great tool you developed. It helped me a lot!

 

Best regards

 

Theo

William Thielicke

unread,
Aug 28, 2019, 2:08:31 AM8/28/19
to PIVlab
Dear Theo, ok great, so I was lucky and chose the best approach when I was programming this. Most likely, I had the information from the same book, but never thought about this in detail. It was good to check and confirm the algorithm again!

What are you actually working on?

thielicke...@gmail.com

unread,
Aug 28, 2019, 3:50:06 AM8/28/19
to PIVlab
If the edge between the nonzero background and the padded zeros is the problem, I will try to smooth this edge... Just because I want to know the difference. I'll report here.

Theo Käufer

unread,
Aug 28, 2019, 6:08:48 AM8/28/19
to PIVlab
Hello William, I am writing a little report about the function of PIV evaluation (Cross-correlation, image deformation, multigrid) and implementing PIVlab functions into Python code since I want to switch from MATLAB to Python and the available OpenPIV is leacking some advanced techniques. (OpenPIV is doing the zero padding)

Therefore, this discussion helped me a lot. I am courios about the results of your smoothed zeropadding.

thielicke...@gmail.com

unread,
Aug 29, 2019, 5:00:51 AM8/29/19
to PIVlab
Here are the results of the smoothing (or maybe better: 'border-smearing') of the zero-padded images. It's becoming a bit better, but not as good as the non-zero-padded data.

Code:
image1_cut = image1_roi(ss1);
image2_cut = image2_roi(ss1);

% padding (faster than als padarray):
image1_cut=[image1_cut zeros(interrogationarea,interrogationarea-1,size(image1_cut,3)); zeros(interrogationarea-1,2*interrogationarea-1,size(image1_cut,3))];
image2_cut=[image2_cut zeros(interrogationarea,interrogationarea-1,size(image2_cut,3)); zeros(interrogationarea-1,2*interrogationarea-1,size(image2_cut,3))];

% border smoothing to remove High-Frequency noise....
blank=floor(interrogationarea/2)
image1_cut(1:interrogationarea+blank,interrogationarea+1:interrogationarea+blank,:)=nan;
image1_cut(interrogationarea+1:interrogationarea+blank,1:interrogationarea+blank,:)=nan;
image2_cut(1:interrogationarea+blank,interrogationarea+1:interrogationarea+blank,:)=nan;
image2_cut(interrogationarea+1:interrogationarea+blank,1:interrogationarea+blank,:)=nan;
for numsi=1:size(image1_cut,3)
   image1_cut(:,:,numsi)=inpaint_nans(image1_cut(:,:,numsi),4);
   image2_cut(:,:,numsi)=inpaint_nans(image2_cut(:,:,numsi),4);
end
% do fft2:
result_conv = fftshift(fftshift(real(ifft2(conj(fft2(image1_cut)).*fft2(image2_cut))), 1), 2);
% cropping of correlation matrix:
00_original.jpg
00_original_corrmat.jpg
01_padding.jpg
01_padding_corrmat.jpg
02_padding+edgesmearing.jpg
02_padding+edgesmearing_corrmat.jpg

alex.l...@gmail.com

unread,
Aug 29, 2019, 6:30:27 AM8/29/19
to PIVlab
Thanks, William for mentioning this conversation. 

Theo - just to make sure that you get the right answers - I'd like to invite you to the openpiv-users forum https://groups.google.com/forum/#!forum/openpiv-users
and OpenPIV does have multi-pass, it's just written in a pre-compiled Cython for the sake of speed. 

I learned quite a lot from this discussion, both about the way PIVLab was written and about the small details of the FFT-based correlation. I suggest to discuss the two things separately:

1. linear vs circular correlation using FFT in a general sense
2. FFT-based correlation in OpenPIV and PIVLab (and MATPIV and other open source projects - don't know about TSI, Dantec, Lavision code) and the small tricks used to get the best (! - we have to work on this definition) vector fields. 

My take on 1 - of course if we take a window of 32 pixels and perform an FFT and IFFT on it - we get a circular correlation in the sense that we can replicate this window and it will give the same frequencies so it's periodic. "to kill" the periodicity - we pad with zeros and then we kind of get an approximation to linear correlation. I am not sure what has to be the length of zero padding to get to the same(!) as just as linear shift and multiplication loop, but I guess double size helps. I teach experimental course and there I explain that zero-padding is only improving resolution in proximity of the spectral peaks because it kills the effect of sharp edges on the FFT. Sharp edges are understood in this periodic sense as connection of the last pixel with the first pixel (it's a period of 32 pixels) and this introduces a lot of leakage around the peaks of significant energy. Some energy is lost to the "wrong" frequencies. 

Regarding the practical use of this for PIV - the background and zero-padding do create some additional problem and therefore in some codes the background is subtracted first - is it so in PIVLab? I remember it's so in MATPIV and so in OpenPIV in Python and Matlab (and if I'm not wrong also in C++) 

I attach below the script from OpenPIV to demonstrate it: 

a) subtract the mean - so the background is now almost zero (sometimes it's even negative ...) 
b) pad with zeros, using the FFT of a double - interrogation window - size (shape + shape is 64 pixels for 32 pixel windows)
c) FFT and IFFT and pick up the central region. 
see the code below


Now to the point 2 - philosophically speaking - none of the correlation methods or other pattern matching methods using particles can claim that it is the only right way to estimate displacement between two images. It is kind of obvious that all are just approximations to the mathematical model that says that two images have greatest similarity at the point of the most probable displacement. We could also introduce rotation and shear and intensity adjustments, Gaussian filtering before and after, intensity capping and whatever else you can invent to "improve" the situation. I believe that probably some adjustments are better for small displacements, some are better for sharp images and some are better for high dynamic range cases. PIV is an experimental method and it will remain such in the sense that we have to look at the final output and be able to estimate its uncertainty versus other measurement system or some solid physics. Of course we do not want some simple bug or a biased estimator to destroy all the vectors, but as it was discussed above - probably one gets somewhat better result at smaller vectors and a bit worse for large ones and then a mutlipass takes another take and somehow improves the situation iteratively. I guess it could be shown also that for some choice of parameters, also multipass is worse than a singlepass :) 

So these are just my humble $0.02 on the topic. 
Alex

    if corr_method == 'fft':
       window_b = np.conj(window_b[::-1, ::-1])
       if nfftx is None:
           nfftx = nextpower2(window_b.shape[0] + window_a.shape[0])  
       if nffty is None:
           nffty = nextpower2(window_b.shape[1] + window_a.shape[1])
       
        f2a = rfft2(normalize_intensity(window_a), s=(nfftx, nffty))
       f2b = rfft2(normalize_intensity(window_b), s=(nfftx, nffty))
       corr = irfft2(f2a * f2b).real
       corr = corr[:window_a.shape[0] + window_b.shape[0],
                    :window_b.shape[1] + window_a.shape[1]]
       return corr
   elif corr_method == 'direct':
       return convolve2d(normalize_intensity(window_a),
       normalize_intensity(window_b[::-1, ::-1]), 'full')
   else:
       raise ValueError('method is not implemented')


def normalize_intensity(window):
   """Normalize interrogation window by removing the mean value.
   Parameters
   ----------
   window :  2d np.ndarray
       the interrogation window array
   Returns
   -------
   window :  2d np.ndarray
       the interrogation window array, with mean value equal to zero.
   """
   return window - window.mean()


thielicke...@gmail.com

unread,
Aug 29, 2019, 7:17:06 AM8/29/19
to PIVlab
Hi Alex, thanks for your input and explanations! Oh yes... subtracting the background makes a lot of sense :-D I didn't think of this simple solution somehow. I will try this too (just because it's fun to play around with the algorithms again), maybe a combination with the 'edge smearing' might even help...

In my opinion, the "best correlation algorithm" is different for the first and for subsequent passes: The first pass should be as robust as possible, and it doesn't even matter if it is wrong by 0.5 pixels. Later passes need to be less robust, but they should be as precise as possible.

So I think for the first pass of a multi-pass analysis, we can use some really hardcore image filtering. Recently, I implemented "5x repeated correlation", which makes an analysis of bad image data more robust. Maybe something like "20x repeated correlation" would make sense for the first pass of a multipass analysis...

Cheers,
William

alex.l...@gmail.com

unread,
Aug 29, 2019, 1:39:51 PM8/29/19
to PIVlab
Agree. 

5x repeated correlation is nice - maybe it's even better to shift by non-integer shifts not to fall into stronger peak locking. It is very similar to what is used in microPIV called ensemble correlation. If the recording is fast enough, the flow velocity does not change for 4-5 frames. So they do correlation on 4-5 frames and then ensemble averaged it - in the correlation space - to get a better peak. 

Nice idea, we probably could implement it in Python. 

Theo Käufer

unread,
Sep 1, 2019, 10:55:42 AM9/1/19
to PIVlab

Hello Alex and William,

 

@William thank you for sharing your results its interesting to see that the smoothing leads to better results.

 

@Alex thank you for invitation to the OpenPIV forum and mentioning the multi-pass function implemented in OpenPIV. I wasn't aware about it.

 

1. In my opinion the signals should be zero-padded to the size of M+N-1 where M is the size of the first window and N is the size of the second window since the size of the direct linear cross-correlation is M+N-1. This is just enough zero-padding to avoid the impact of the periodicity and also the size of the invers Fourier transform of the elementwise multiplication int the frequency domain of the zero-padded windows is M+N-1. This is used in Fast Convolution/Cross-Correlation algorithms.

I think you can in theory even detect a displacement that is bigger than window size/2 it the correlation peak is clear enough using this since you receive the full correlation matrix.

However, following the rule of thumb that the maximum displacement should be less than 1/4 of the window size we can crop the resulting correlation matrix to the window size (or even below?).

 

In chapter 5.3.1.1  Raffel, Markus, et al. Particle image velocimetry: a practical guide. Springer, 2018. are mentioning some conditions under which circular cross-correlation is providing sufficient displacement information.

 

They mention a maximum displacement which is about 1/4 of the window size and they suggest implementing a weighting function to avoid a bias effect because there is less correlation data available at the borders of the signal. However, I think things will change if multipass is used.

 

2. Thank you for mentioning that. I agree that there is no perfect solution for every situation. I also want to stress that further "improvement" often comes with more computational efforts that should not be underestimated!

 

Let me summarize the difference in the approaches of PIVlab and OpenPIV to calculate the Cross correlation:

 

PIVlab is doing no zero-padding and therefore does circular cross-correlation.

 

OpenPIV is doing zero-padding and therefore linear cross correlation but does background subtraction to avoid the errors caused by the high frequency noise.

 

I am really interested in a comparison between PIVlab and OpenPIV under the similar conditions and I will most likely perform one by my own, but it will probably take me some time. I will let you know about the results!

 

Concerning the multipass I completely agree with you.


Best regards,

Theo

William Thielicke

unread,
Sep 4, 2019, 5:18:05 AM9/4/19
to PIVlab
It is not over yet... ;-D
I generated some very bad PIV images with a lot of particle pair loss and noise. Then, I quickly compared the algorithms (this time zero-padding WITH subtracting the mean before the correlation), and the zero-padded algorithm seems to be more robust (see images). However, it is also by a factor of 8.9 slower (mainly because of the larger size of the matrices for the correlation).
In the next days, I will do more comparisons on robustness and precision, then I'll decide if I accept the almost 10x slower speed.
padding_001.jpg
no_padding_002.jpg

William Thielicke

unread,
Sep 4, 2019, 10:49:49 AM9/4/19
to PIV...@googlegroups.com
I found some old accuracy test scripts, so the comparison went quicker than I thought...
These are the results for a single pass analysis with 24 px interrogation area on artificial images with 3px diameter particles. Every image shows the results for padded vs. non-padded correlation.

Most likely, I will add a slider to PIVlab, where the user can choose between speed and quality of the analysis. The slider will enable / disable different parameters that affect speed and quality, e.g. spline vs. linear window deformation; circular vs. linear correlation; repeated correlation; and so on...

noise_rms.jpg
noise_wrong.jpg
noise_bias.jpg
shift_bias.jpg
shift_rms.jpg
loss_bias.jpg
loss_rms.jpg

Theo Käufer

unread,
Sep 5, 2019, 2:18:03 PM9/5/19
to PIVlab
The results are very intresting. Have you performed a multipass interrogation with the different correlation methods and compared the vectorfields?
Maybe the accuary missing in the first pass can get compensated by the following passes. I think you can even perform some more passes with the non-zeropadding method and still be below the efforts off the zeropadding version.
What do you think?

William Thielicke

unread,
Sep 6, 2019, 5:39:34 AM9/6/19
to PIV...@googlegroups.com
Hi, here are the results of a more detailed test:
Linear correlation is only 2.5x slower (the above results were valid for linear correlation with repeated correlation enabled). And linear correlation is better (less bias error, less rms error, more robust). I will check if it helps to enable this only in the final pass.

Edit: There is hardly any difference when linear correlation and / or repeated correlation is only enabled in the final pass. So the latest PIVlab version does this only in the final pass. That saves 50% of the computing time.
The reson is that spurious vectors in earlier passes are filtered and interpolated - so a less robust correlation algorithm is not really a big problem. The following algorithms will start searching at the correct position.

me15...@smail.iitm.ac.in

unread,
Sep 6, 2019, 6:13:12 AM9/6/19
to PIVlab
Hello William, 

Is the correlation updates will be available in the update version of PIV Lab? 

When we are suppose to get the PIV lab 2.1.

Thanks!

Regards

Aaditya

William Thielicke

unread,
Sep 6, 2019, 7:36:54 AM9/6/19
to PIVlab
Hi, yes this will be included in the coming release. There are quite a number of changes, and I am currently having some small problems with saving and loading PIVLab settings files. I hope I can finish everything this weekend.

William Thielicke

unread,
Sep 6, 2019, 12:51:56 PM9/6/19
to PIVlab
Reply all
Reply to author
Forward
0 new messages