Multi-Photon Correlation Statistics

In [1]:
from PIL import Image as IM
import numpy as np
import random
import matplotlib.pyplot as plt

%matplotlib inline
plt.rcParams['figure.figsize'] = (10.0, 8.0) # set default size of plots
plt.rcParams['image.interpolation'] = 'nearest'
plt.rcParams['image.cmap'] = 'gray'

%load_ext autoreload
%autoreload 2

Tuning Hyperparameters

In [2]:
# Load a standard photo
im = IM.open('/Users/YRZ/Downloads/1.jpg')
In [3]:
# Convert it to gray-scale image
im = im.convert("L")
plt.imshow(im, cmap = plt.cm.gray_r)
plt.show()
In [4]:
# Check image size
print(im.size)

(1024, 1024)
In [5]:
# rotate 0 degrees
rot = im.rotate(0)
plt.imshow(rot, cmap = plt.cm.gray_r)
plt.show()
In [6]:
# mask(left, upper, right, lower) CAUTION: Delta's should be ingeters!
mask = (445, 384, 625, 520)
region = im.crop(mask)
plt.imshow(region, cmap = plt.cm.gray_r)
plt.show()

# Detailed check: 
# Adjust the mask until detailed sub-images become perfect.

# Num of Rows and Cols
NUM_R, NUM_C = 4, 5
# Delta length of Height and Width
DELTA_H = (mask[3] - mask[1]) / NUM_R
DELTA_W = (mask[2] - mask[0]) / NUM_C
plt.title("first row")
plt.imshow(region.crop((0,0,DELTA_W*NUM_C,DELTA_H)), cmap = plt.cm.gray_r)
plt.show()
plt.title("first col")
plt.imshow(region.crop((0,0,DELTA_W,DELTA_H*NUM_R)), cmap = plt.cm.gray_r)
plt.show()
plt.title("last row")
plt.imshow(region.crop((0,DELTA_H*(NUM_R-1),DELTA_W*NUM_C,DELTA_H*NUM_R)), cmap = plt.cm.gray_r)
plt.show()
plt.title("last col")
plt.imshow(region.crop((DELTA_W*(NUM_C-1),0,DELTA_W*NUM_C,DELTA_H*NUM_R)), cmap = plt.cm.gray_r)
plt.show()
In [7]:
# Filters
from PIL import ImageFilter

def filter(img):
#     fimg = img.filter(ImageFilter.EDGE_ENHANCE)
    fimg = img.filter(ImageFilter.MedianFilter(3))
#     fimg = fimg.filter(ImageFilter.SHARPEN)
    return fimg

imfilter = filter(region)
plt.imshow(imfilter, cmap = plt.cm.gray_r)
plt.show()

See a sample

In [8]:
# Load a sample
smp = IM.open('/Users/YRZ/Downloads/photos/minibatch/DATA_FASTEST762.jpg')

# Check Size
print(smp.size)
smp = smp.convert('L')
smp = smp.crop(mask)

# Check Original Image in Target Region
plt.title('original image')
plt.imshow(smp, cmap = plt.cm.gray)
plt.show()

# Check Filtered Image in Target Region
smp = filter(smp)
plt.title('after filter')
plt.imshow(smp, cmap = plt.cm.gray_r)
plt.show()

(1024, 1024)
In [9]:
# Segment into sub-images
def secs(img, sub_imgs):
    for i in range(NUM_R):
        for j in range(NUM_C):
            sub_imgs.append(img.crop((DELTA_W * j, DELTA_H * i, DELTA_W * (j + 1), DELTA_H * (i + 1))))
    
sub_smps = []
secs(smp, sub_smps)

# Convert to vectors
sub_smps_datas = []
for s in sub_smps:
    sub_smps_datas.append(np.array(s))
#     sub_smps_datas.append(np.matrix(s.getdata(), dtype='float') / 255.0)
#     plt.imshow(s, cmap = plt.cm.gray)
#     plt.show()
#     print(s.size)

# Find a section to see the best threshold
plt.imshow(sub_smps[1], cmap = plt.cm.gray_r )
plt.show()
tmp = np.reshape(sub_smps_datas[1], (int(DELTA_W), int(DELTA_H)))
print('Thus we think the best threshold is about', np.max(tmp))

Thus we think the best threshold is about 57
In [10]:
# Set threshold
THRESHOLD = 50
In [11]:
# Count Correlations!

# Here, corr is a vector encoding photon distribution
# Order: form the upper left to the bottom right.
corr = np.zeros(NUM_R*NUM_C)
ind = 0
for s in sub_smps_datas:
    if np.max(s) >= THRESHOLD:
        corr[ind] = 1
    ind += 1

print(corr)

[ 0.  1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.
  1.  0.]
In [12]:
# Plot Correlation

import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D

mpl.rcParams['font.size'] = 10

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

dx = np.ones(10)
dy = np.ones(10)

# Question: Should we always take the case of two photons hitting the same place into account?
# Answer1: Yes.
zvals = []
for x in range(1, NUM_R * NUM_C + 1):
    xs = [x - 1]* (NUM_R * NUM_C)
    ys = range(1, NUM_R * NUM_C + 1)
    dz = corr[x - 1] * corr[range(NUM_R * NUM_C)]
    zvals.append(dz)
    zs = np.zeros(NUM_R * NUM_C)
    ax.bar3d(xs, ys, zs, dx = 0.80*np.ones(NUM_R * NUM_C), 
             dy = 0.80*np.ones(NUM_R * NUM_C), dz = dz, color='#00ceaa', alpha=0.8)

ax.set_xlabel('Position1')
ax.set_ylabel('Position2')
ax.set_zlabel('Frequency')

plt.show()

# cmap = plt.cm.coolwarm
cmap = mpl.colors.LinearSegmentedColormap.from_list('my_colormap',['#003333','#FFFF33','#99FFFF'], 256)

img_show = plt.imshow(zvals,interpolation='nearest',
                    cmap = cmap,
                    origin='lower')
plt.colorbar(img_show ,cmap = cmap)
plt.show()
In [13]:
# Question: Should we always take the case of two photons hitting the same place into account?
# Answer2: No. We only take account of overlapping when there is only one hitten place.

mpl.rcParams['font.size'] = 10

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

dx = np.ones(10)
dy = np.ones(10)

zvals = []
for x in range(1, NUM_R * NUM_C + 1):
    xs = [x - 1]* (NUM_R * NUM_C)
    ys = range(1, NUM_R * NUM_C + 1)
    dz = corr[x - 1] * corr[range(NUM_R * NUM_C)]
    if np.sum(corr) > 1: 
        dz[x - 1] = 0
    zvals.append(dz)
    zs = np.zeros(NUM_R * NUM_C)
    ax.bar3d(xs, ys, zs, dx = 0.80*np.ones(NUM_R * NUM_C), 
             dy = 0.80*np.ones(NUM_R * NUM_C), dz = dz, color='#00ceaa', alpha=0.8)

ax.set_xlabel('Position1')
ax.set_ylabel('Position2')
ax.set_zlabel('Frequency')

plt.show()

img_show = plt.imshow(zvals,interpolation='nearest',
                    cmap = cmap,
                    origin='lower')
plt.colorbar(img_show ,cmap = cmap)
plt.show()
# I prefer Answer2.

Pipeline

In [14]:
import os

# set file pattern and path
FlagStr= 'jpg'
# FindPath = '/Users/YRZ/Downloads/photos/minibatch'
FindPath = '/Users/YRZ/Downloads/photos'

# Piple function to get correlation vector
def get_corr(img_name):
    img = IM.open(img_name)
    img = img.convert('L')
    img = img.crop(mask)
    img = filter(img)
    sub_imgs = []
    secs(img, sub_imgs)
    sub_imgs_datas = []
    for s in sub_imgs:
        sub_imgs_datas.append(np.array(s))
    corr = np.zeros(NUM_R*NUM_C)
    ind = 0
    for s in sub_imgs_datas:
        if np.max(s) >= THRESHOLD:
            corr[ind] = 1
        ind += 1
    return corr

# Visit all files in the directory 
container = []
FileNames = os.listdir(FindPath)
for fn in FileNames:
    if (len(FlagStr)>0) and (FlagStr in fn):
        fullfilename = os.path.join(FindPath,fn)
        container.append(get_corr(fullfilename))
In [15]:
# Plot the final result!

mpl.rcParams['font.size'] = 10

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

dx = np.ones(10)
dy = np.ones(10)

zvals = []
for x in range(1, NUM_R * NUM_C + 1):
    xs = [x - 1]* (NUM_R * NUM_C)
    ys = range(1, NUM_R * NUM_C + 1)
    dz = np.zeros(NUM_R * NUM_C)
    for corr_count in container:
        dz_tmp = corr_count[x - 1] * corr_count[range(NUM_R * NUM_C)]
        if np.sum(corr) > 1: 
            dz_tmp[x - 1] = 0
        dz += dz_tmp
    zs = np.zeros(NUM_R * NUM_C)
    zvals.append(dz)
    ax.bar3d(xs, ys, zs, dx = 0.80*np.ones(NUM_R * NUM_C), 
             dy = 0.80*np.ones(NUM_R * NUM_C), dz = dz, color='#00ceaa', alpha=0.8)

ax.set_xlabel('Position1')
ax.set_ylabel('Position2')
ax.set_zlabel('Frequency')

plt.show()

img_show = plt.imshow(zvals,interpolation='nearest',
                    cmap = cmap,
                    origin='lower')
plt.colorbar(img_show, cmap = cmap)
plt.show()
In [ ]: