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
# Load a standard photo
im = IM.open('/Users/YRZ/Downloads/1.jpg')
# Convert it to gray-scale image
im = im.convert("L")
plt.imshow(im, cmap = plt.cm.gray_r)
plt.show()
# Check image size
print(im.size)
# rotate 0 degrees
rot = im.rotate(0)
plt.imshow(rot, cmap = plt.cm.gray_r)
plt.show()
# 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()
# 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()
# 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()
# 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))
# Set threshold
THRESHOLD = 50
# 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)
# 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()
# 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.
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))
# 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()