Computes brightness covariance of two images and plots corresponding 2D histogram

pixel-covariance.py

#!/usr/bin/env python
#-*- coding: utf-8 -*-

"""

Program accepts two parameters: image1 image2

Both files have to have the same dimension. Any format accepted by scipy is possible. 
They may be grayscale or RGB, in the latter case the R+G+B value is taken.

By default, the images are slightly blurred to reduce noise (see user settings section below) 

Outputs:
    1) The calculated covariance of all pixels, is appended into 'covariance.txt' in the format: 
            <image1>   <image2>   <covariance>
        For same images, covariance should be 1, for random noise, it should be close to 0. For negatives, 
        covariance is -1.
        You may run this script multiple times and the results gather in this file for further 
        processing

    2) 2D covariance map is saved into the 'corr<image1><image2>.png'
        Below you may change the number of bins of this 2d histogram.

    3) Additionally, total number of "counts" is saved. For this, we assume that both images were normalized; 
        that is, if  one image contains only black and white, the white pixels are assumed to represent one 
        count each. If there are 30 levels, the brightest pixel is assumed to represent 30 counts, etc.

(c) 2017, Filip Dominec, MIT Licensed
"""


## Import common moduli
import matplotlib, sys, os, time
import matplotlib.pyplot as plt
plt.figure(figsize=(8,8))
import numpy as np
from scipy.constants import c, hbar, pi


## =========  User settings: ==========
numbin  = 32        

#blurkernel = [1]
blurkernel = [.25,.5,.25]
#blurkernel = [1,3,5,3,1]
#blurkernel = np.exp(-np.linspace(-2.,2.,11, dtype=np.float32)**2)

## =========  /User settings ==========




## Load images (of the same size!)
from scipy import misc
im1 = misc.imread(sys.argv[1])
im2 = misc.imread(sys.argv[2])
name1, name2 = (os.path.split(sys.argv[1])[1], os.path.split(sys.argv[2])[1]) ## convenient names

## Mage sure images are grayscale 
try: im1 = im1.sum(axis=2)
except: pass
try: im2 = im2.sum(axis=2) 
except: pass


## Un-normalization of brightness levels 
#print("numbers of distinct brightness (R+G+B) levels (img1,img2):", len(np.unique(im1)), len(np.unique(im2)))
norm1 = np.max(im1) / (len(np.unique(im1))-1) ## assume normalized image: the more brightness levels, the greater weight 
norm2 = np.max(im2) / (len(np.unique(im2))-1) ## note that this is   unreliable if outlier-bright points present
#print('norms', norm1, norm2)
#print('sums', np.sum(im1), np.sum(im2))

#print("# of unique values:", len(np.unique(im1)), len(np.unique(im2)))
with open('counts.txt', 'a') as stats:
    stats.write("Total detections for %s: %g" % (name1, np.sum(im1)/norm1)) 
    stats.write("Total detections for %s: %g" % (name2, np.sum(im2)/norm2))

## 
from scipy.signal import sepfir2d
im1c = sepfir2d(im1, blurkernel, blurkernel)
im2c = sepfir2d(im2, blurkernel, blurkernel)

maxrange = max(np.max(im1c), np.max(im2c))
corr = np.zeros([numbin+1,numbin+1])

for x in range(im1c.shape[0]):
    for y in range(im1c.shape[1]):
        corr[int(im1c[x,y]/maxrange*numbin), int(im2c[x,y]/maxrange*numbin)] += 1

avg1,avg2 = np.sum(im1c)/im1c.size, np.sum(im2c)/im2c.size
var1, var2 = np.sum((im1c-avg1)**2), np.sum((im2c-avg2)**2)
cov       = (np.sum((im1c-avg1)*(im2c-avg2)) - avg1*avg2) / (var1*var2)**.5
with open('covariance.txt', 'a') as covfile:
    covfile.write('{}\t{}\t{}\n'.format(name1, name2, cov))

#plt.imshow(im1); #plt.show()
corr[0,0] = 0
plt.imshow(np.log10(corr+1)); #plt.show()
plt.ylim((-0.6,numbin+.6)); plt.yscale('linear')
plt.xlim((-0.6,numbin+.6)); plt.xscale('linear')
plt.ylabel(sys.argv[1])
plt.xlabel(sys.argv[2])
#plt.colorbar()

plt.savefig("corr%s%s.png" % (name1,name2), bbox_inches='tight')
#plt.imshow(np.log10(im1c+1)); plt.show()
#print(np.max(corr))

An example of two correlated luminescence images from our electron microscope:

The result is quite a good, but not full, covariance:
uv-band.jpg blue-band.jpg 0.891632371754