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