# Cross-correlation function in Julia vs. Python

**URL:** <https://discourse.julialang.org/t/cross-correlation-function-in-julia-vs-python/7573>\
**Category:** Statistics\
**Tags:** question, statistics\
**Created:** [December 6, 2017, 5:02pm UTC](https://discourse.julialang.org/t/cross-correlation-function-in-julia-vs-python/7573 "2017-12-06T17:02:09Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![Duls91](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/duls91/32/2858_2.png) [@Duls91](https://discourse.julialang.org/u/Duls91)\
**Post date:** [December 6, 2017, 5:02pm UTC](https://discourse.julialang.org/t/cross-correlation-function-in-julia-vs-python/7573/1 "2017-12-06T17:02:09Z")

</div>

Hi there,

I am new to Julia. I want to perform a cross-correlation with two 2d arrays (both 5X5). Using crosscor() (StatsBase.jl) returns a 9x5x5 array. What I would like is the cross-correlation matrix of dims M+N-1, such as the output of the scipy.signal.correlate2d() from Python, giving a 9x9 Matrix.  
Here is my code:

```julia
a = ones(5,5)
b = ones(5,5)
crosscor(a,b, demean=false)
Output:
9×5×5 Array{Float64,3}:

```

What kind of cross-correlation is implemented here? I am having a hard time grasping this.  
If I use the 2D convolution with a reversed array, I do get what I want:  
`conv2(a,b[end:-1:1,end:-1:1])`  
Output:  
9×9 Array{Float64,2}

but the implementation is FFT-based. I would like a direct convolution, but I cannot find one in Julia.

hope you can help me clarify this! Thanks!

---

<div class="post-metadata">

**Author:** ![ExpandingMan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/expandingman/32/866_2.png) [@ExpandingMan](https://discourse.julialang.org/u/ExpandingMan)\
**Post date:** [December 6, 2017, 7:32pm UTC](https://discourse.julialang.org/t/cross-correlation-function-in-julia-vs-python/7573/2 "2017-12-06T19:32:18Z")

</div>

See the documentation by doing `?crosscor` in the REPL.

I agree that the documentation in this case is not as clear as it could be. The dimension of length 9 is showing you expectation values for different “lags”, i.e. I think what you are seeing here in the vector case would be.  
 \langle (x\_{i}-\langle x\rangle)(x\_{i+j} - \langle x\rangle)\rangle  
for a range of j (given in the documentation), but it would be really nice if the exact expression were written in the documentation. The remaining two indices are showing you these values for each pair of columns.

---

<div class="post-metadata">

**Author:** ![GunnarFarneback](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gunnarfarneback/32/1827_2.png) [@GunnarFarneback](https://discourse.julialang.org/u/GunnarFarneback)\
**Post date:** [December 6, 2017, 11:10pm UTC](https://discourse.julialang.org/t/cross-correlation-function-in-julia-vs-python/7573/3 "2017-12-06T23:10:01Z")

</div>

This can be solved by `imfilter` and `padarray` from the `ImageFiltering` package but if your arrays are small you will likely get perfectly decent performance by implementing the function yourself, just writing the nested loops involved.

---

<div class="post-metadata">

**Author:** ![Duls91](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/duls91/32/2858_2.png) [@Duls91](https://discourse.julialang.org/u/Duls91)\
**Post date:** [December 7, 2017, 3:30pm UTC](https://discourse.julialang.org/t/cross-correlation-function-in-julia-vs-python/7573/4 "2017-12-07T15:30:28Z")

</div>

Thank you for the answers so far. I see that the concept of cross-correlation here is quite different then I need it to be. I would like to see whether one image a is part of the other image b. So I would slide one image over the other and see for which displacement I get the highest cross-correlation value.

Python gives the following output:  
a = np.ones((3,3))  
b = np.ones((3,3))  
c = ss.correlate2d(a, b)

```julia
array([[1., 2., 3., 2., 1.], 
          [2., 4., 6., 4., 2.],
          [3., 6., 9., 6., 3.], 
          [2., 4., 6., 4., 2.], 
          [1., 2., 3., 2., 1.]])

```

So the peak in my correlation matrix is 9, as the center of it, meaning that without displacement, the two images are the most “similar”.

Julia on the other hand gives me the following:

```julia
a = ones(3,3)
b = ones(3,3)
crosscor(a,b, demean=false)*3

[:, :, 1] = 1.0 1.0 1.0 
              2.0 2.0 2.0 
              3.0 3.0 3.0 
              2.0 2.0 2.0 
              1.0 1.0 1.0 
[:, :, 2] = 1.0 1.0 1.0 
              2.0 2.0 2.0 
              3.0 3.0 3.0 
              2.0 2.0 2.0 
              1.0 1.0 1.0

```

And so on…So it does not tell me anything.  
I might implement a function resembling the Python output. I would just like to know the difference in the Python CC and the Julia one.

---

<div class="post-metadata">

**Author:** ![ExpandingMan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/expandingman/32/866_2.png) [@ExpandingMan](https://discourse.julialang.org/u/ExpandingMan)\
**Post date:** [December 7, 2017, 3:46pm UTC](https://discourse.julialang.org/t/cross-correlation-function-in-julia-vs-python/7573/5 "2017-12-07T15:46:46Z")

</div>

I’m definitely not an image processing person by any means, but it might be worth thinking about whether this is really an effective measure for similarity of images. I’d be doubtful about how far you could take the concept of “lag” or displacement here especially since there is no sense in which x\_{i+1} follows x\_{i} in a sequence if they lie on the boundary of a row of pixels.

Just thinking out loud here, I would think that one thing you could do would be to subtract the images from each other and perform some sort of test to see if the difference is “random noise”. Of course this is easier said than done, and perhaps it would be a good place for some sort of dimensionality reduction.

Having said all that. I would think that the conceptually simplest way to do what you are trying to do would be to just vectorize the image. Then you can just use `crosscor` and know exactly what you are looking at.

You also might want to check out the various packages [here](https://github.com/JuliaImages).

---

<div class="post-metadata">

**Author:** ![Duls91](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/duls91/32/2858_2.png) [@Duls91](https://discourse.julialang.org/u/Duls91)\
**Post date:** [December 7, 2017, 4:06pm UTC](https://discourse.julialang.org/t/cross-correlation-function-in-julia-vs-python/7573/6 "2017-12-07T16:06:26Z")

</div>

This cross-correlation principle is commonly used for image similarity (like in particle image velocimetry - the thing I want to implement in Julia finally). But I could certainly look also into other “similarity”, i.e. “disimilarity” measures like the subtracting that you proposed.  
Also, I will look more into the vectorization and then try it out.

Thanks for the ideas and the link!

---

<div class="post-metadata">

**Author:** ![ExpandingMan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/expandingman/32/866_2.png) [@ExpandingMan](https://discourse.julialang.org/u/ExpandingMan)\
**Post date:** [December 7, 2017, 4:10pm UTC](https://discourse.julialang.org/t/cross-correlation-function-in-julia-vs-python/7573/7 "2017-12-07T16:10:22Z")

</div>

> [@Duls91](#):
>
> This cross-correlation principle is commonly used for image similarity

Like I said, I’m a layman, so I defer to your judgment. The reason I questioned it is because it is far from obvious to me what the expectations about what the cross-correlation would be in the general case.

---

<div class="post-metadata">

**Author:** ![GunnarFarneback](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gunnarfarneback/32/1827_2.png) [@GunnarFarneback](https://discourse.julialang.org/u/GunnarFarneback)\
**Post date:** [December 7, 2017, 4:25pm UTC](https://discourse.julialang.org/t/cross-correlation-function-in-julia-vs-python/7573/8 "2017-12-07T16:25:54Z")

</div>

> [@ExpandingMan](#):
>
> I’m definitely not an image processing person by any means, but it might be worth thinking about whether this is really an effective measure for similarity of images.

It’s a standard image processing technique, except you need to do some normalization to avoid getting the idea that a white area in the image is most similar to your template. E.g. like this:

```julia
using Images

function correlate_image(img, template)
    z = imfilter(Float64, img, centered(template), Algorithm.FIR())
    q = sqrt.(imfilter(Float64, img.^2, centered(ones(size(template))),
                       Algorithm.FIR()))
    m = sqrt(sum(template.^2))
    corr = z ./ (m * max.(eps(), q))
    return corr
end

```

For illustration:

```julia
julia> x=rand(1:9, 8, 8)
8×8 Array{Int64,2}:
 7 2 2 9 7 2 4 7
 3 7 9 5 3 2 9 5
 1 9 1 4 1 4 2 8
 7 6 3 1 8 8 6 7
 4 9 4 4 1 6 6 2
 4 8 7 1 8 2 1 8
 5 4 2 5 5 3 4 1
 3 8 3 8 4 5 7 8

julia> y = x[2:4, 5:7]
3×3 Array{Int64,2}:
 3 2 9
 1 4 2
 8 8 6

julia> c = correlate_image(x,y)
8×8 Array{Float64,2}:
 0.661888 0.737724 0.908773 0.822907 0.602351 0.699639 0.862974 0.875247
 0.573246 0.637993 0.802854 0.615533 0.602121 0.626725 0.815854 0.789346
 0.884582 0.892484 0.633597 0.678266 0.740834 1.0 0.81138 0.83531 
 0.827208 0.685541 0.795446 0.605092 0.714524 0.644814 0.824203 0.741265
 0.772589 0.803458 0.695693 0.8739 0.736878 0.742672 0.715778 0.796131
 0.870059 0.706784 0.678821 0.599637 0.903347 0.756111 0.604329 0.58481 
 0.847023 0.830675 0.713355 0.892543 0.76089 0.728849 0.960917 0.8947  
 0.74215 0.778543 0.781013 0.849 0.773851 0.82505 0.754464 0.763592

julia> c[3,6]
1.0

```

Notes:

- You can substitute ImageFiltering for Images if you don’t want the whole package.
- You rather want to force spatial filtering rather than FFT in case you have a sufficiently large area of zeros in your image, or you have to figure out exactly how to deal with zeros being corrupted into epsilons of varying signs.
- Never mind the variable names. That code was never meant to go public and I won’t make the effort to come up with better ones.

---

<div class="post-metadata">

**Author:** ![eljungsk](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eljungsk/32/2865_2.png) [@eljungsk](https://discourse.julialang.org/u/eljungsk)\
**Post date:** [December 7, 2017, 8:49pm UTC](https://discourse.julialang.org/t/cross-correlation-function-in-julia-vs-python/7573/9 "2017-12-07T20:49:09Z")

</div>

If your goal is to implement PIV, I guess that the two input arrays are of the same size and that you want to perform the correlation a large number of times for similar inputs.

I implemented a simple multipass PIV code for fun a year ago, and found that the performance could be significantly increased by writing a custom correlation function using FFT. Something along the lines of this might work for you (no guarantees of correctness 🙂):

```julia
# The size of your interrogation area (assumed here to be square)
interrogationArea = 32

# Dummy data
a = rand(interrogationArea, interrogationArea)
b = rand(interrogationArea, interrogationArea)

# Plan the forward and inverse transforms for better performance
p = plan_fft(rand(Float64, interrogationArea, interrogationArea), flags = FFTW.MEASURE)
ip = plan_ifft(rand(Complex128, interrogationArea, interrogationArea), flags = FFTW.MEASURE)

function correlate_fft(p, ip, a, b)
       corr = fftshift(real(ip*(conj(p*a).*(p*b))))
end

```

A quick benchmark show a speedup of roughly 50x compared to the ‘’‘correlate\_images’‘’ function proposed above, but this solution is of course much less general.

---

<div class="post-metadata">

**Author:** ![Duls91](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/duls91/32/2858_2.png) [@Duls91](https://discourse.julialang.org/u/Duls91)\
**Post date:** [December 8, 2017, 12:59pm UTC](https://discourse.julialang.org/t/cross-correlation-function-in-julia-vs-python/7573/10 "2017-12-08T12:59:02Z")

</div>

Thanks for the code and your reply!

> [@eljungsk](#):
>
> I guess that the two input arrays are of the same size and that you want to perform the correlation a large number of times for similar inputs.

Yes, it will be a lot of calculations. I’ll be going for some GPU computing for this.

> [@eljungsk](#):
>
> A quick benchmark show a speedup of roughly 50x compared to the ‘’‘correlate\_images’’’ function proposed above, but this solution is of course much less general.

I want to implement both direct and more faster methods, but would like to start with the more general one and test my PIV performance with all the different nuances.

---

<div class="post-metadata">

**Author:** ![Duls91](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/duls91/32/2858_2.png) [@Duls91](https://discourse.julialang.org/u/Duls91)\
**Post date:** [December 8, 2017, 1:03pm UTC](https://discourse.julialang.org/t/cross-correlation-function-in-julia-vs-python/7573/11 "2017-12-08T13:03:13Z")

</div>

Thank you very much for the code contribution!

It is very close to what I want to have. The resulting matrix should be 10x10 (8 + 3 - 1). But I will play around with the code and do the adjustments.

---

<div class="post-metadata">

**Author:** ![jianghaizhu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jianghaizhu/32/1866_2.png) [@jianghaizhu](https://discourse.julialang.org/u/jianghaizhu)\
**Post date:** [February 26, 2018, 4:17pm UTC](https://discourse.julialang.org/t/cross-correlation-function-in-julia-vs-python/7573/12 "2018-02-26T16:17:29Z")

</div>

I am trying to follow your approach to calculate cross-correlations using FFT. In my case, I have millions of images to compare. Thus, the function you provided allocates way too much memories and GC would be very busy. In your code, every step will allocate a temporary array. I tried to eliminate the allocations as much as possible.

```nohighlight
function corr_fft(P, IP, A, B)
	cc = fftshift(real.(IP*(conj(P*A) .* (P*B))))
end

fftshift!(y, x) = circshift!(y, x, div.([size(x)...],2))

function corr_fft1!(cc, P, IP, A, B, XC)
	FA = P*A
	conj!(FA)
	FB = P*B
	FA .= FA .* FB
	A_mul_B!(FB, IP, FA)
	XC .= real.(FB)
	fftshift!(cc, XC)
end

```

In this version, `cc` and `XC` are pre-allocated real arrays outside of the function. But I still have to allocate `FA` and `FB`, two complex arrays, for each call. Unfortunately, if one wants to use `A_mul_B!(Y, plan, X)`, the input image `X` has to be a complex array as well. But if the number of calls to this function is greater than the number of images, it might be worth to convert all images to complex arrays first and use the following version.

```nohighlight
function corr_fft2!(cc, P, IP, AC, BC, FA, FB, XC)
	A_mul_B!(FA, P, AC)
	conj!(FA)
	A_mul_B!(FB, P, BC)
	FA .= FA .* FB
	A_mul_B!(FB, IP, FA)
	XC .= real.(FB)
	fftshift!(cc, XC)
end

```

Here, `cc` and `XC` are pre-allocated real arrays, `AC` and `BC` are converted complex images, and `FA` and `FB` are pre-allocated complex arrays. For `100 x 100` pixel images, I benchmarked the three versions in my MacBook pro. `AC` and `BC` are complex arrays converted from `A` and `B`.

```nohighlight
julia> @btime corr_fft($P, $IP, $A, $B)
  204.048 μs (22 allocations: 625.95 KiB)

julia> @btime corr_fft1!($cc, $P, $IP, $A, $B, $XC)
  193.094 μs (11 allocations: 313.03 KiB)

julia> @btime corr_fft2!($cc, $P, $IP, $AC, $BC, $FA, $FB, $XC)
  177.746 μs (3 allocations: 224 bytes)

```

The last version is much faster than the first two when loops through millions of images. Now, if only there is a way to directly do `A_mul_B!(Y, plan, X)` without requiring converting all my images to complex arrays, that would be super! I wonder if there is any method in FFTW package to do that. I didn’t find any in the manual.

---

<div class="post-metadata">

**Author:** ![jwu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jwu/32/3775_2.png) [@jwu](https://discourse.julialang.org/u/jwu)\
**Post date:** [March 21, 2018, 7:29pm UTC](https://discourse.julialang.org/t/cross-correlation-function-in-julia-vs-python/7573/13 "2018-03-21T19:29:49Z")

</div>

FYI, our lab have an optimized 2D implementation. It was already used to compare millions of images.  
[https://github.com/seung-lab/Alembic/blob/master/src/math/convolve\_inplace.jl#L520](https://github.com/seung-lab/Alembic/blob/master/src/math/convolve_inplace.jl#L520)