How to calculate a 2D spatial correlation function for scalar data on a 2D grid?

Question

Donovan Allum 2022-8-30

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1790985-how-to-calculate-a-2d-spatial-correlation-function-for-scalar-data-on-a-2d-grid

评论： Donovan Allum 2023-12-16

在 MATLAB Online 中打开

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.

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
radfield = 1×3
         0    5.6914    7.0944
r
r = 1×3
    0.5000    1.5000    2.5000

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Balavignesh 2023-12-14

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1790985-how-to-calculate-a-2d-spatial-correlation-function-for-scalar-data-on-a-2d-grid#answer_1371309

在 MATLAB Online 中打开

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)
    0.4973    0.5001    0.4998    0.4999    0.5003    0.4998    0.4996    0.4997    0.4997    0.4998