- 'meshgrid': https://www.mathworks.com/help/matlab/ref/meshgrid.html
- 'conv2': https://www.mathworks.com/help/matlab/ref/conv2.html
How to calculate a 2D spatial correlation function for scalar data on a 2D grid?
31 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
I have a 2D velocity field
. I want to take
and calculate
where the angled brackets denote an average over all grid points in a concentric shell centred at
with radius
. I currently have some code that does this but since I'm working with a 512x512 grid it takes ~8 hours to run. This isn't surprising because the number of operations scales like
, where N is the number of grid points in each direction. But the code I'm using uses loops and really I'm wondering if it can't be made more efficient with the use of some combination of xcorr2 and conv2 (I am aware that xcorr2 is basically just a specific application of conv2). I suspect it is possible but I do not understand exactly what these functions do and how they might speed up my problem.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1112050/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1112055/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1112060/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1112065/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1112070/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1112075/image.png)
At the bottom is the code I currently have for a very small w (Note that data is NOT periodic)
I am aware that currently I think my code is "over counting" because wi*wj and wj*wi are both included so ideas to make this more efficient would also be helpful.
clear; close all
x = [1 1 1; 2 2 2; 3 3 3];
x = x(:);
y = [1 2 3; 1 2 3; 1 2 3];
y = y(:);
w = [0.8147 0.9134 0.2785; 0.9058 0.6324 0.5469; 0.1270 0.0975 0.9575];
sz = size(x);
Nx = sz(1);
Ny = sz(2);
w = w(:);
Nbins=3;
maxr=max(x);
dr = maxr/Nbins;
radfield = zeros(1,Nbins);
count = radfield;
kk=0;
for ii=1:Nx
for jj=1:Ny
%Define central point
kk = kk+1;
vec = [1:(kk-1) (kk+1):Nx*Ny];
xij = x(kk);
yij = y(kk);
%Create a vector which gives the distance between (xij,yij)
%and every other grid point, store as a vector.
rmag=sqrt((x(vec)-xij).^2+(y(vec)-yij).^2);
%Store velocity at (xij,yij)
velij = w(kk);
%Calculate the product with velij with every other grid point
velvel = velij.*w(vec);
%Average over concentric shells of thickness dr.
for ll = 1:Nbins
myfilt=(rmag<ll*dr&rmag>(ll-1)*dr); % the filter for the shell
radfield(ll)=radfield(ll)+sum(velvel.*(myfilt)); % the mean over the shell
count(ll)=count(ll)+sum(myfilt);
end
end
end
r = (1:Nbins)*dr-0.5*dr;
radfield
r
0 Commenti
Risposte (1)
Balavignesh
il 14 Dic 2023
Hi Donavan,
As per my understanding, you have a 2-dimensional velocity field and your goal is to calculate the average velocity over concentric shells centered at each grid point in a 2D velocity field. The computational challenge arises from the size of the grid (512 X 512), which results in a large number of operations, especially when using nested loops.
Instead of looping through each grid point to create a shell mask, I would suggest you create a single "kernel" matrix that represents the shell. This kernel will have same dimensions as your velocity field, with values of 1 inside the shell and 0 outside.
You could use 'conv2' function to apply this kernel acorss the entrie velocity field in a single operation. This operation sums the products of the velocity field and the kernel at each point, yielding a matrix of summed velocities. You could normalize the shell area by dividing the summed velocities with the area of the shell.
The following example code may help you understand this:
% Define the grid size
N = 512;
w = rand(N); % Example 2D velocity field
% Define the number of bins and the maximum radius
Nbins = 10; % Example number of bins
maxr = N / 2; % Example maximum radius
dr = maxr / Nbins;
% Preallocate the radial field array
radfield = zeros(1, Nbins);
% Create a grid of distances from the center
[X, Y] = meshgrid(1:N, 1:N);
Xc = N / 2;
Yc = N / 2;
r = sqrt((X - Xc).^2 + (Y - Yc).^2);
% Calculate the radial field for each bin
for bin = 1:Nbins
% Create the kernel for the shell
shell_mask = (r >= (bin - 1) * dr) & (r < bin * dr);
% Sum the velocity field within the shell using conv2
shell_sum = conv2(w, double(shell_mask), 'same');
% Count the number of points in the shell
shell_count = conv2(ones(size(w)), double(shell_mask), 'same');
% Calculate the average velocity within the shell
radfield(bin) = sum(shell_sum(shell_mask)) / sum(shell_count(shell_mask));
end
disp(radfield)
Kindly have a look at the following documentation links to have more information on:
Hope this helps!
Balavignesh
Vedere anche
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!