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
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.
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.
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.
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_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);
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,:);
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;
% 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));
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
What are you actually working on?
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:
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()
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.