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

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

Answer 1

[Balavignesh]

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

Kindly have a look at the following documentation links to have more information on:

• 'meshgrid': https://www.mathworks.com/help/matlab/ref/meshgrid.html
• 'conv2': https://www.mathworks.com/help/matlab/ref/conv2.html

Hope this helps!

Balavignesh

[Donovan Allum]

Hi Balavignesh,

I think our approaches are quite similar except

1. you make use of conv2, and
2. You only look at a single central point.

My goal is to loop over all possible centres (this is the definition of autocorrelation functions).

I could easily adapt your code by adding the outer loops and perhaps your use of conv2 is faster.

I think my 'myfilt' array is the same as your 'shell_mask' array.

What would be ideal would be to vectorize the loop over bins. But its possible my memory requirements will skyrocket.