Understanding correlation in image processing using examples

Question

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.

Answer 1

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]:

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

def simple_convolve(img, mask):
    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):
            for curr_mask_row in range(0, mask.shape[0]):
                for curr_mask_col in range(0, mask.shape[1]):
                    output[curr_row, curr_col] += img[curr_row + curr_mask_row, curr_col + curr_mask_col] * mask[curr_mask_row, curr_mask_col]
    return output
​
def better_convolve(img, mask):
    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])

    mask_mean  = mask.mean()
    mask_sigma = (((mask - mask_mean)**2).sum())**0.5
    mask = mask - mask_mean
    mask = mask/mask_sigma

    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
            for curr_mask_row in range(0, mask.shape[0]):
                for curr_mask_col in range(0, mask.shape[1]):
                    ctr = ctr + 1.0
                    curr_img_mean += img[curr_row + curr_mask_row, curr_col + curr_mask_col]
            curr_img_mean = curr_img_mean / ctr

            # Compute the std value of the sub-image of img
            ctr = 0.0
            for curr_mask_row in range(0, mask.shape[0]):
                for curr_mask_col in range(0, mask.shape[1]):
                    ctr = ctr + 1.0
                    curr_img_sigma += (img[curr_row + curr_mask_row, curr_col + curr_mask_col] - curr_img_mean)**2.0
            curr_img_sigma = curr_img_sigma**0.5

            for curr_mask_row in range(0, mask.shape[0]):
                for curr_mask_col in range(0, mask.shape[1]):         
                    x = (img[curr_row + curr_mask_row, curr_col + curr_mask_col] - curr_img_mean)
                    if curr_img_sigma == 0.0:
                        x = 0.0
                    else:
                        x = x/curr_img_sigma

                    output[curr_row, curr_col] += x * mask[curr_mask_row, curr_mask_col]
    return output
In [5]:

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

better_convolve(img, mask)
Out[6]:
array([[ 0.        ,  0.        ,  0.        ],
       [-1.        , -0.57735027,  0.        ],
       [ 1.        ,  0.57735027,  0.        ]])

Answer 2

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

Answer 3

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     | ?     | ? | ? |
+-------+-------+---+---+
| ?     | ?     | ? | ? |
+-------+-------+---+---+
| ?     | ?     | ? | ? |
+-------+-------+---+---+
| ?     | ?     | ? | ? |
+-------+-------+---+---+

Answer 4

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.