How to calculate cross-correlation on a 2D domain

Specific Domains Statistics
1

Hi everyone,

I’m studying a scalar field u(x,t) on a 2D box LxL, with periodic boundary conditions.
What I would like to do is to calculate the cross-correlation between u(x,t1) and u(x+r,t1) at varying of r, where r is the euclidean distance.

After discretized the domain (then now it could be see as a lattice), my problem is about how to reach, for each node of the lattice, the points on the circle of radius r.

I have implemented a simple way to do it, but it doesn’t work with the periodic boundary conditions:

r=4 #for instance
dd=1:N
kk=[Int64[] for i in dd]

for i in dd
    
    
    for j in dd
        if norm(A[:,i]-A[:,j], 2)==r
            append!(kk[i],[j])
        end
    end
end

Where A is the 2xN matrix containing the coordinates of each node of the lattice. Here the kk vector contains for each node his “neighbors” at the euclidean distance r, but it not consider the bc.

I’m sure that exist a better way to solve this probem, every help is welcome, thanks!

2

Hi @Frazze ,

are you trying to do an Empirical Variogram?
If so, I think GitHub - JuliaEarth/GeoStats.jl: An extensible framework for geospatial data science and geostatistical modeling fully written in Julia has an implementation.

I think you probably want some tolerance on the exact distance if you’ve discretized your data (i.e. sth like abs(norm(A[:,i]-A[:,j], 2) - r) < ϵ

For the periodic boundary conditions, maybe GitHub - Vexatos/CircularArrays.jl: Multi-dimensional arrays with fixed size and circular indexing. could be useful.

1 Like
3

Ah, I think I misunderstood. You’re only looking for 1-dimensional distance? If so, disregard my previous comment. But CircularArrays.jl might still be useful. Or just use the mod1() function (which is part of base julia).

4

Thank @trahflow for the answer,

I only want to know what points are at distance r from a fixed point (x1,y1), and my troubles are at the boundary of my domain, here the norm can’t see all the points at the correct distance.
Now I will take a look at CircularArrays.jl, but how do you think it could be used?

Thanks a lot for the help

5

@Frazze do you need to preserve the flat coordinates with spherical boundary conditions or you are after an actual sphere geometry? Notice that you can perform the variogram calculations linked above on spheres just fine. Try to define your domain as domain = triangulate(Sphere((0.,0.), 1.0)) in GeoStats.jl and see if that addresses the problem.

We are still polishing the features of the project related to manifolds and non-Euclidean geometry, but maybe you can cook something up already. Try to read the GeoStats.jl and Meshes.jl docs to see if you find something useful.

6

Sorry I went way too quick over your question and thus misunderstood :wink:

is the follownig what you’re looking for?

First define cyclic difference for each coordinate:

function cyclic_coord_diff(a, b)
     d = abs(a - b)
     if d > L/2
         d -= L
     end
     d
 end

then define a distance function like:

distance(a, b) = norm(cyclic_coord_diff.(a, b), 2)

(note the broadcasting of cyclic_coord_diff)

7

The coordinates are only a “trick” for me to try to solve the problem. I have a 2D domain where some ODEs are calculated. Then I want to calculate the cross correlation of the solution between all the points of the box, but I wish to maintain the regular square-lattice discretization.

8

Thank you so much for the help!
I think that this should be the correct way to solve the problem.

9

From experience I would recommend the following approach:

  1. take spatial 2d fourier transform of each field
  2. multiply them for each point in k-space
  3. perform a 1d inverse Hankel transform w.r.t. the norm of k

This elegantly solves the problem of taking a polar average over a cartesian grid.

4 Likes
10

I have decided to implement this approach because the issue seems to pop up and the details (esp normalization and scaling) tend to be quite tricky. You can find the prototype here

https://github.com/gwater/PairCorrelation.jl/blob/be746b8920c637cae5c652bad16b7d8f88c6a679/src/PairCorrelation.jl#L47

11

Thank you so much for the advices! I will take a look at PairCorrelation.jl

12

I’m sorry @josuagrw, I take this post to do another question, still related to the topic.

What if I would like to calculate the power spectra (then the ft of spatial correlation) of mi scalar field at one fixed time t just having the solution of the field u(x,y,t)?

There is a way to do it without calculate before the correlation?

Thank you so much, I’m really new with this stuff…

13

First off: With periodic boundary conditions, the power spectrum is discrete because only specific wavelengths fulfill that boundary condition.

Once you realize that, it actually makes everything easier: The spectral density (power spectrum) is just the square absolute value of the fourier transform, so

ks = rfftkspace2(…)
power_spectrum = abs2.(rfft(data))

You can then decide how you want to plot those discrete data points (for example pick a specific axis or collapse everything onto the norm of k).

1 Like
14

Thank you so much @josuagrw !

I think I understand almost everything. ks contains the kpoints in the 2d kspace?

How can I do If I would like to plot the spectrum collapsing onto the norm of k?

15

correct

there are different ways to do it. The naive way is

using Plots

scatter(
    reshape(norm.(ks), :),
    reshape(power_spectrum, :)
)

which typically looks quite messy (depending on how isotropic your data really is). If you want / need to clean that up more, it really depends on what kind of data you are analysing and how you want to present it.

16

Ok, being k periodic I have truncate the power spectrum after the first period.

Doing it the plot shows:

PowerSpectrum

And this seems to be consistent with what I’m studying.

Thank you for all the advices!

1 Like
17

Hi @Frazze sorry I forgot to tell you about fftshift(): Basically the modes are not ordered sequentially in memory, so plot(power_spectrum) will give you a wrong picture.

This is why it is important to use fftshift or fftfreq (or their equivalentes rfftshift and rfftfreq)

18

Ok thank you,
Just another thing that it’s not so clear for me: What about the array dimension of power_spectrum?
I can’t understand why, starting from a data array of dim= 16384 (a domain 128x128 reshaped) with rfft I reach a vector of dim=8193, in particular this vector dim is different than the k_norm array dim, then I don’t know how to represent it…

19

Basically rfft halves the memory because the Fourier transform of a real input has a certain symmetry. Look here to see how you can calculate the modes in the 2D case

https://github.com/gwater/PairCorrelation.jl/blob/be746b8920c637cae5c652bad16b7d8f88c6a679/src/PairCorrelation.jl#L10

1 Like