# Understanding correlation in image processing using examples

I have been trying to understand the intuition behind the correlation. I understand from this guide that correlation is a way to detect a known waveform in a noisy background. I have seen from this SO answer that correlation can be used for template matching. So I decided to try it and hence the reason for my confusion.

Consider a 2D array,

+---+---+---+---+
| 1 | 1 | 1 | 1 |
+---+---+---+---+
| 1 | 1 | 1 | 1 |
+---+---+---+---+
| 0 | 0 | 1 | 1 |
+---+---+---+---+
| 1 | 1 | 1 | 1 |
+---+---+---+---+

I did correlation without padding with the below array,

+---+---+
| 0 | 0 |
+---+---+
| 1 | 1 |
+---+---+

The output I got is,

+---+---+---+
| 2 | 2 | 2 |
+---+---+---+
| 0 | 1 | 2 |
+---+---+---+
| 2 | 2 | 2 |
+---+---+---+

I was expecting a maximum value where the smaller array was present at the big array. But I got multiple peak values. Can you help me understand why this happened and how to overcome it?

Thank you.

• You picked a tough example. Short answer: Change your "0"s to another value, e.g., 2, and it should work much better. What's really happening: your signals are not zero mean, correlation requires to center the signals (i.e., subtract means). Nov 19 '19 at 12:02
• @Florian Correlation does not require the subtraction of the mean, but covariance does Nov 19 '19 at 13:42
• What do you mean by "I was expecting a maximum value where the smaller array was present at the big array"? Nov 19 '19 at 13:42
• @Engineer: Pearson correlation is using centering. Not to be confused with the correlation matrix, which is not centered. I was referring to correlation in the sense of a correlation coefficient (as a measure of "similarity"). Nov 19 '19 at 16:53
• The question seems to be about cross-correlation, based on the cited reference. @tpk there are other definitions of correlation (like Pearson's mentioned above) so it is good to specify. Nov 19 '19 at 18:11

As it was already posted multiple times: The problem comes from an inaccurate definition of correlation in your application.

The Pearson correlation coefficient does require the data to be

1. centered, ie the mean must be subtracted
2. normalized, ie the data must be divided by the standard deviation

This centering and normalization must be done for the mask as well for each sub-matrix of your larger matrix.

In your example, you would end up with a correlation matrix as:

$$\left( \begin{matrix} 0.0&0.0&0.0\\ -1.0&-0.577&0.0\\ 1.0&0.577&0.0 \end{matrix} \right)$$

Lets have a look at the matrix entries:

• 0.0 They are there, since there are sub-matrices with all ones in them. The mean of such a matrix is 1 and therefore it is centered to the corresponding matrix with only zeros. Hence the correlation is zero.

• 1.0 is the correlation that you would like to see. The corresponding sub-matrix fits the mask perfectly, hence full correlation, i.e. 1

• -1.0 is the correlation of the sub-matrix that fits the mask, but the sub-matrix entries are switched (i.e. the sub-matrix is 0 where the mask is 1 and vice versa). This is a perfect anti-correlation, i.e. -1

• The 0.577 cases: They occur where the sub-matrix deviates from the mask by one number (i.e. there is a 1 where the mask has a 0 and vice versa). It also happens to be a symmetric case in your example, therefore there is a pair that correlates positively with one switch and one that correlates negatively with a corresponding switch.

Here is a pretty ugly brute-force code to illustrate it. It is neither optimized in any way, nor tested for any other case than yours. I hope it is somewhat correct ;)

import numpy as np
In [2]:

img = np.ones([4,4])
img[2,0:2] = 0.0
print(img)
[[1. 1. 1. 1.]
[1. 1. 1. 1.]
[0. 0. 1. 1.]
[1. 1. 1. 1.]]
In [3]:

[[0. 0.]
[1. 1.]]
In [4]:

max_row    = img.shape[0] - mask.shape[0] + 1
max_col    = img.shape[1] - mask.shape[1] + 1

output = np.zeros([max_row, max_col])

for curr_row in range(0, max_row):
for curr_col in range(0, max_col):
return output
​
max_row    = img.shape[0] - mask.shape[0] + 1
max_col    = img.shape[1] - mask.shape[1] + 1

output = np.zeros([max_row, max_col])

curr_img_mean  = 0.0
curr_img_sigma = 0.0

for curr_row in range(0, max_row):
for curr_col in range(0, max_col):

curr_img_mean  = 0.0
curr_img_sigma = 0.0

# Compute the mean value of the sub-image of img
ctr = 0.0
ctr = ctr + 1.0
curr_img_mean = curr_img_mean / ctr

# Compute the std value of the sub-image of img
ctr = 0.0
ctr = ctr + 1.0
curr_img_sigma = curr_img_sigma**0.5

if curr_img_sigma == 0.0:
x = 0.0
else:
x = x/curr_img_sigma

return output
In [5]:

Out[5]:
array([[2., 2., 2.],
[0., 1., 2.],
[2., 2., 2.]])
In [6]:

Out[6]:
array([[ 0.        ,  0.        ,  0.        ],
[-1.        , -0.57735027,  0.        ],
[ 1.        ,  0.57735027,  0.        ]])
• @Florain's answer is equally good and, it is crisp and short. Accepting this answer as it is more descriptive.
– tpk
Nov 20 '19 at 12:49
• No complaints, this one is definitely more verbose :) Nov 20 '19 at 14:03

You picked a tough example.

Short answer: Change your "0"s to another value, e.g., 2, and it should work much better.

What's really happening: your signals are not zero mean, correlation requires to center the signals (i.e., subtract means). Example, since it's easer to understand in 1-D: say you want to find the pattern p=[0,2,2,0] in the sequence s=[2,1,1,1,0,1,1,0,1,1,1,2]. You'd expect to find it in the middle (note that scaling of the pattern is ignored in correlations). Now compare:

conv(s,p,'same')
>> 6   4   4   2   2   4   2   2   4   4   6   4

with

conv(s-mean(s),p-mean(p),'same')
>> 1  -1   1  -1  -1   2  -1  -1   1  -1   1   1

There's your peak, right in the middle :)

• Convolution is not always the same thing as correlation Nov 19 '19 at 18:01
• That's true, thanks for noticing it. I was trying to simplify here to convey the basic concept. Nov 20 '19 at 14:01

This might be hard with trying to illustrate this way but here I go. Take your original image and put the filter in the upper left corner like so:

+-------+-------+---+---+
| 1 [0] | 1 [0] | 1 | 1 |
+-------+-------+---+---+
| 1 [1] | 1 [1] | 1 | 1 |
+-------+-------+---+---+
| 0     | 0     | 1 | 1 |
+-------+-------+---+---+
| 1     | 1     | 1 | 1 |
+-------+-------+---+---+

The bracketed values are the values of the filter. Now compute the correlation: $$1*0 + 1*0 + 1*1 + 1*1 = 2$$. So we assign the value $$2$$ to the corresponding location in the result image:

+---+---+---+
| 2 | ? | ? |
+---+---+---+
| ? | ? | ? |
+---+---+---+
| ? | ? | ? |
+---+---+---+

Now we shift the filter, lets shift one pixel to the right.

+---+-------+-------+---+
| 1 | 1 [0] | 1 [0] | 1 |
+---+-------+-------+---+
| 1 | 1 [1] | 1 [1] | 1 |
+---+-------+-------+---+
| 0 | 0     | 1     | 1 |
+---+-------+-------+---+
| 1 | 1     | 1     | 1 |
+---+-------+-------+---+

Again, we compute the correlation: $$1*0 + 1*0 + 1*1 + 1*1 = 2$$. And again assign the pixel in the correlation image:

+---+---+---+
| 2 | 2 | ? |
+---+---+---+
| ? | ? | ? |
+---+---+---+
| ? | ? | ? |
+---+---+---+

And you continue this until you have shifted the filter all the way through the original image. I hope this helps!

EDIT:

Since another answer has disputed mine, I will expand. In image processing, correlation is defined as: $$Y(r, c) = \sum_{r'=-h}^{h} \sum_{c'=-w}^{w} I(r+r', c+c') F(r', c')$$, where the image $$I$$ has dimensions $$(2h+1) \times (2w+1)$$ and the filter $$F$$ is assumed to be using some sort of padding (https://www2.cs.duke.edu/courses/fall15/compsci527/notes/convolution-filtering.pdf). Intuitively, the equation is doing what I already outlined in my answer: put the filter over the image, calculate the sum of product of pixels, shift the filter through the image. In fact, MATLAB is a popular tool to perform image processing tasks and they implement a function called imfilter which you can try out to check your answers! They provide an easy to understand tutorial on correlation (https://www.mathworks.com/help/images/what-is-image-filtering-in-the-spatial-domain.html).

image = [1 1 1 1; 1 1 1 1; 0 0 1 1; 1 1 1 1];
filt = [0 0; 1 1];
imfilter(image, filt)

MATLAB actually computes additional values along the boundaries like this:

[0]     [0]
+-------+-------+---+---+
| 1 [1] | 1 [1] | 1 | 1 |
+-------+-------+---+---+
| 1     | 1     | 1 | 1 |
+-------+-------+---+---+
| 0     | 0     | 1 | 1 |
+-------+-------+---+---+
| 1     | 1     | 1 | 1 |
+-------+-------+---+---+

And by default it uses zero padding to do the same calculation as I have shown in my answer: $$0*0 + 0*0 + 1*1 + 1*1 = 2$$. And then they fill the corresponding result pixel. The result is that the output has the same size as the input image.

+-------+-------+---+---+
| 2     | ?     | ? | ? |
+-------+-------+---+---+
| ?     | ?     | ? | ? |
+-------+-------+---+---+
| ?     | ?     | ? | ? |
+-------+-------+---+---+
| ?     | ?     | ? | ? |
+-------+-------+---+---+
• Great to see answers from new comers to our community. +1.
– Royi
Nov 19 '19 at 21:20
• +1 for sharing how correlation in image processing works. This answer has made it more clear. My question was to understand how correlation is used for template matching and why it failed in this particular example in the question.
– tpk
Nov 20 '19 at 12:37

That is because correlation (and convolution) are not meant to "match" exactly a given pattern. They are multiplicative operators in their nature so they are strongly related to signal amplitude if you multiply the reference signal by N, the output gets twice bigger.

for instance your operator will return a peak twice higher when encountering this piece of signal:

+---+---+
| 0 | 0 |
+---+---+
| 2 | 2 |
+---+---+

than when encountering this one:

+---+---+
| 0 | 0 |
+---+---+
| 1 | 1 |
+---+---+

Moreover, with your operator the result when encountering this signal:

+---+---+
|999|999|
+---+---+
| 1 | 1 |
+---+---+

will be the same as above because of the zeros.

In fact your operator acts like a low-pass horizontal filter in this case. As stated by @M529 the Pearson Correlation Coefficient better suits your application.

Another interesting operator would be the Mean Square Error or any other error function. Its output is 0 when the signal matches perfectly, and >0 when the error gets bigger.

• Great to see answers from new comers to our community. +1.
– Royi
Nov 19 '19 at 21:20