Python - Normalized cross-correlation to measure similarites in 2 images I'm trying to measure per-pixel similarities in two images (same array shape and type) using Python.

In many scientific papers (like this one), normalized cross-correlation is used. Here's an image from the ict paper showing the wanted result:



(b) and (c) are the 2 input images, and (d) is the per-pixel confidence.

I only used OpenCV before to do template matching with normalized cross correlation using cv2.matchTemplate function, but in this case it seems to be a really different use of cross correlation.

Do you know if I can approch this result using Python and image processing libraries (numpy, openCV, sciPy etc...), and the logic behind this use of cross correlation ? 

Thanks.
I guess you can compute for each pixel the correlation coefficient between patches centered on this pixel in the two images of interest. Here is an example where I downloaded the figure attached here and tried to compute the correlation in such a way. The output looks different from the one of the article, but it was to be expected since the resolution is very different.

from skimage import io, feature
from scipy import ndimage
import numpy as np

def correlation_coefficient(patch1, patch2):
    product = np.mean((patch1 - patch1.mean()) * (patch2 - patch2.mean()))
    stds = patch1.std() * patch2.std()
    if stds == 0:
        return 0
    else:
        product /= stds
        return product

im = io.imread('faces.jpg', as_grey=True)

im1 = im[16:263, 4:146]
sh_row, sh_col = im1.shape
im2 = im[16:263, 155:155+sh_col]

# Registration of the two images
translation = feature.register_translation(im1, im2, upsample_factor=10)[0]
im2_register = ndimage.shift(im2, translation)

d = 1

correlation = np.zeros_like(im1)

for i in range(d, sh_row - (d + 1)):
    for j in range(d, sh_col - (d + 1)):
        correlation[i, j] = correlation_coefficient(im1[i - d: i + d + 1,
                                                        j - d: j + d + 1],
                                                    im2[i - d: i + d + 1,
                                                        j - d: j + d + 1])

io.imshow(correlation, cmap='gray')
io.show()




The code above is a naive and slow implementation of the correlation, as the two for loops are very slow. For faster execution, you could for example port the script to Cython.

In the article, I think the idea is to measure whether face expressions look similar or not. If a pixel has a large correlation index between two images, it means that the region of the face where this pixel is located does not change much between the images. 

If you use this method on good-resolution images, you should increase the patch size for more accurate results (d=2 or 3).