# How to calculate cross-correlation on a 2D domain

**URL:** <https://discourse.julialang.org/t/how-to-calculate-cross-correlation-on-a-2d-domain/75161>\
**Category:** Statistics\
**Tags:** question\
**Created:** [January 25, 2022, 9:47am UTC](https://discourse.julialang.org/t/how-to-calculate-cross-correlation-on-a-2d-domain/75161 "2022-01-25T09:47:05Z")\
**Posts on this page:** 19\
**Page:** 1

<div class="post-metadata">

**Author:** ![Frazze](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/frazze/32/25142_2.png) [@Frazze](https://discourse.julialang.org/u/Frazze)\
**Post date:** [January 25, 2022, 9:47am UTC](https://discourse.julialang.org/t/how-to-calculate-cross-correlation-on-a-2d-domain/75161/1 "2022-01-25T09:47:05Z")

</div>

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:

```julia
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!

---

<div class="post-metadata">

**Author:** ![trahflow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/trahflow/32/30585_2.png) [@trahflow](https://discourse.julialang.org/u/trahflow)\
**Post date:** [January 25, 2022, 10:24am UTC](https://discourse.julialang.org/t/how-to-calculate-cross-correlation-on-a-2d-domain/75161/2 "2022-01-25T10:24:26Z")

</div>

Hi @Frazze ,

are you trying to do an [Empirical Variogram](https://en.wikipedia.org/wiki/Variogram#Empirical_variogram)?  
If so, I think [GitHub - JuliaEarth/GeoStats.jl: An extensible framework for geospatial data science and geostatistical modeling fully written in Julia](https://github.com/JuliaEarth/GeoStats.jl) has an implementation.

> [@Frazze](#):
>
> ` if norm(A[:,i]-A[:,j], 2)==r`

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.](https://github.com/Vexatos/CircularArrays.jl) could be useful.

---

<div class="post-metadata">

**Author:** ![trahflow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/trahflow/32/30585_2.png) [@trahflow](https://discourse.julialang.org/u/trahflow)\
**Post date:** [January 25, 2022, 10:30am UTC](https://discourse.julialang.org/t/how-to-calculate-cross-correlation-on-a-2d-domain/75161/3 "2022-01-25T10:30:19Z")

</div>

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).

---

<div class="post-metadata">

**Author:** ![Frazze](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/frazze/32/25142_2.png) [@Frazze](https://discourse.julialang.org/u/Frazze)\
**Post date:** [January 25, 2022, 10:48am UTC](https://discourse.julialang.org/t/how-to-calculate-cross-correlation-on-a-2d-domain/75161/4 "2022-01-25T10:48:01Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [January 25, 2022, 10:56am UTC](https://discourse.julialang.org/t/how-to-calculate-cross-correlation-on-a-2d-domain/75161/5 "2022-01-25T10:56:51Z")

</div>

@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.

---

<div class="post-metadata">

**Author:** ![trahflow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/trahflow/32/30585_2.png) [@trahflow](https://discourse.julialang.org/u/trahflow)\
**Post date:** [January 25, 2022, 11:13am UTC](https://discourse.julialang.org/t/how-to-calculate-cross-correlation-on-a-2d-domain/75161/6 "2022-01-25T11:13:02Z")

</div>

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

is the follownig what you’re looking for?

First define cyclic difference for each coordinate:

```julia
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:

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

```

(note the broadcasting of `cyclic_coord_diff`)

---

<div class="post-metadata">

**Author:** ![Frazze](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/frazze/32/25142_2.png) [@Frazze](https://discourse.julialang.org/u/Frazze)\
**Post date:** [January 25, 2022, 11:17am UTC](https://discourse.julialang.org/t/how-to-calculate-cross-correlation-on-a-2d-domain/75161/7 "2022-01-25T11:17:52Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Frazze](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/frazze/32/25142_2.png) [@Frazze](https://discourse.julialang.org/u/Frazze)\
**Post date:** [January 25, 2022, 12:07pm UTC](https://discourse.julialang.org/t/how-to-calculate-cross-correlation-on-a-2d-domain/75161/8 "2022-01-25T12:07:49Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [January 25, 2022, 1:23pm UTC](https://discourse.julialang.org/t/how-to-calculate-cross-correlation-on-a-2d-domain/75161/9 "2022-01-25T13:23:42Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [January 27, 2022, 12:16pm UTC](https://discourse.julialang.org/t/how-to-calculate-cross-correlation-on-a-2d-domain/75161/10 "2022-01-27T12:16:27Z")

</div>

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](https://github.com/gwater/PairCorrelation.jl/blob/be746b8920c637cae5c652bad16b7d8f88c6a679/src/PairCorrelation.jl#L47)

---

<div class="post-metadata">

**Author:** ![Frazze](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/frazze/32/25142_2.png) [@Frazze](https://discourse.julialang.org/u/Frazze)\
**Post date:** [January 27, 2022, 1:05pm UTC](https://discourse.julialang.org/t/how-to-calculate-cross-correlation-on-a-2d-domain/75161/11 "2022-01-27T13:05:23Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Frazze](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/frazze/32/25142_2.png) [@Frazze](https://discourse.julialang.org/u/Frazze)\
**Post date:** [February 3, 2022, 9:51am UTC](https://discourse.julialang.org/t/how-to-calculate-cross-correlation-on-a-2d-domain/75161/12 "2022-02-03T09:51:19Z")

</div>

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…

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [February 3, 2022, 11:36am UTC](https://discourse.julialang.org/t/how-to-calculate-cross-correlation-on-a-2d-domain/75161/13 "2022-02-03T11:36:44Z")

</div>

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

```julia
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).

---

<div class="post-metadata">

**Author:** ![Frazze](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/frazze/32/25142_2.png) [@Frazze](https://discourse.julialang.org/u/Frazze)\
**Post date:** [February 4, 2022, 8:50am UTC](https://discourse.julialang.org/t/how-to-calculate-cross-correlation-on-a-2d-domain/75161/14 "2022-02-04T08:50:03Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [February 4, 2022, 9:54am UTC](https://discourse.julialang.org/t/how-to-calculate-cross-correlation-on-a-2d-domain/75161/15 "2022-02-04T09:54:12Z")

</div>

> [@Frazze](#):
>
> I think I understand almost everything. ks contains the kpoints in the 2d kspace?

correct

> [@Frazze](#):
>
> How can I do If I would like to plot the spectrum collapsing onto the norm of k?

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

```julia
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.

---

<div class="post-metadata">

**Author:** ![Frazze](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/frazze/32/25142_2.png) [@Frazze](https://discourse.julialang.org/u/Frazze)\
**Post date:** [February 4, 2022, 10:59am UTC](https://discourse.julialang.org/t/how-to-calculate-cross-correlation-on-a-2d-domain/75161/16 "2022-02-04T10:59:56Z")

</div>

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

Doing it the plot shows:

![PowerSpectrum](https://global.discourse-cdn.com/julialang/original/3X/9/1/91c548d5e39f0e7c322fff77df98148c3dd4a7e1.png)

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

Thank you for all the advices!

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [February 4, 2022, 11:30am UTC](https://discourse.julialang.org/t/how-to-calculate-cross-correlation-on-a-2d-domain/75161/17 "2022-02-04T11:30:01Z")

</div>

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](https://juliamath.github.io/AbstractFFTs.jl/stable/api/#AbstractFFTs.fftshift) _or_ [fftfreq](https://juliamath.github.io/AbstractFFTs.jl/stable/api/#AbstractFFTs.fftfreq) (or their equivalentes rfftshift and rfftfreq)

---

<div class="post-metadata">

**Author:** ![Frazze](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/frazze/32/25142_2.png) [@Frazze](https://discourse.julialang.org/u/Frazze)\
**Post date:** [February 6, 2022, 1:08pm UTC](https://discourse.julialang.org/t/how-to-calculate-cross-correlation-on-a-2d-domain/75161/18 "2022-02-06T13:08:20Z")

</div>

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…

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [February 6, 2022, 8:23pm UTC](https://discourse.julialang.org/t/how-to-calculate-cross-correlation-on-a-2d-domain/75161/19 "2022-02-06T20:23:50Z")

</div>

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](https://github.com/gwater/PairCorrelation.jl/blob/be746b8920c637cae5c652bad16b7d8f88c6a679/src/PairCorrelation.jl#L10)