Ghost imaging¶
Images are saved into compressed numpy arrays, they can be load with the np.load() function.
In [139]:
import numpy as np
import matplotlib.pyplot as plt
npzdir = '/home/valerio/1slit_thermal_npz'
npzfile = 'luce_0000.npz'
imagearray = np.load(npzdir+'/'+npzfile)
Remember that a numpy npz files is a dictionary-like container: it means that the array(s) inside it are accessible through a key. The key list is shown with the files method.
In [140]:
print('The keys in the ',npzfile,'file are', imagearray.files)
Then the image is accessible using the 'arr_0' key:
In [141]:
image = imagearray['arr_0']
print('The image contained in the',npzfile,' has a size of :',image.shape)
Let's write a small function to open an image from a npz file for later use
In [142]:
def load_image(f = 'luce_0000.npz'):
return np.load(f)['arr_0']
Plot one image and its histogram (use a function for later use)
In [143]:
image = load_image('luce_0000.npz')
def plotimage(image) :
fig, ax = plt.subplots(1,2, figsize=(18,5))
# left plot (image)
ax[0].set(title=f'Image - size = {image.shape[0]}x{image.shape[1]} pixels', xlabel = 'horiz', ylabel='vert')
h = ax[0].imshow(image, cmap='gray')
plt.colorbar(h, ax=ax[0])
# right plot (histogram)
ax[1].set(title='Image histogram', xlabel = 'pixel value', ylabel='nr. of pixels')
ax[1].hist(image.flatten(), bins=20, color='gray', log=True)
return fig, ax
plotimage(image)
plt.show()
Dark images¶
Some dark images are included in the data set. In order to be able to detect them we have to fix some threshold and define what "dark" is. Let's sum the pixel values for a well enlighted region (e.g. the right side) and check their distribution:
In [145]:
light_list = [] # prepare a list for later cross-check
# check the first 200 images
for i in range(200):
# open the image
im = load_image(f'{npzdir}/luce_{i:04}.npz')
# sum the pixel values in a given region and store the results in a list
light = im[200:500, 500:800].sum()
light_list.append(light)
plt.hist(light_list, 50)
darkthr = 200000 # <-- threshold to detect dark images
plt.axvline(darkthr, color='r', linestyle='--')
plt.annotate('dark', xy=(darkthr, 17), xytext=(darkthr-50000,17),ha='left', va='center', arrowprops=dict(arrowstyle="<-",
connectionstyle="arc3", color='red'), fontsize=14)
plt.annotate('not dark', xy=(darkthr, 17), xytext=(darkthr+70000,17),ha='right', va='center', arrowprops=dict(arrowstyle="<-",
connectionstyle="arc3", color='red'), fontsize=14)
plt.show()
Let's check which images have a sum below 200000
In [146]:
print('List of the (probably) dark images:')
for i,l in enumerate(light_list):
if l < 200000:
print('Image nr. ',i,' -> total light ',l)
Let's check for example image nr. 41 and 133 (we can use the previously defined function to do the plot):
In [147]:
image = load_image(f'{npzdir}/luce_0041.npz')
plotimage(image)
plt.show()
image = load_image(f'{npzdir}/luce_0133.npz')
plotimage(image)
plt.show()
Dark images has to be excluded from the analysis.
A "dark run" is also available in single image, each pixel contains the average value
In [148]:
image = load_image(f'{npzdir}/luce_0040.npz')
average_dark = np.loadtxt('media_buio0.txt')
plotimage(average_dark)
plt.show()
The bucket region¶
Sum of 50 images to have a better idea on where the 2 light beams are located
In [149]:
from matplotlib.patches import Rectangle
total = np.zeros(image.shape)
for i in range(50):
# open the image
im = load_image(f'{npzdir}/luce_{i:04}.npz')
total = total + im
bx1, bx2, by1, by2 = 80, 140, 260, 530
bucket = image[by1:by2, bx1:bx2].sum()
fig, ax = plt.subplots(1,3, figsize=(20,5))
# projection along the 1st axis (rows). All columns will be squeezed into a single row (1 array)
projection1 = total.sum(axis=0)
# projection along the 2nd axis (columns). All rows will be squeezed into a single columns
# It's better to perform this sum excluding the much more enlighted region of the image
projection2 = total[:,:200].sum(axis=1)
ax[0].plot(projection1, '.')
ax[0].axvline(bx1, color='red', linestyle='--')
ax[0].axvline(bx2, color='red', linestyle='--')
ax[1].plot(projection2, '.')
ax[1].axvline(by1, color='red', linestyle='--')
ax[1].axvline(by2, color='red', linestyle='--')
ax[2].imshow(total, cmap='gray')
ax[2].add_patch(Rectangle((bx1,by1), width = bx2-bx1, height = by2-by1, ec='r', fc='none'))
plt.show()
Main analysys¶
Loop on all the images:
- do not consider the dark images and store them separately
- use different (rectangular shaped) buckets with different size. From a single pixel to the whole bucket region
In this version of the analysis there is no background subtraction
In [150]:
# file location
npzdir = '/home/valerio/1slit_thermal_npz'
npzfile = 'luce_0000.npz'
# Function do small N to tune the analysis, then N = 10000 for the final version
def doanalysis(N):
# prepare the empy images to store everything (open the first one to properly set the shape)
im0 = load_image(f'{npzdir}/luce_0000.npz')
dark, total, total_corr1, total_corr2, total_corr3 = np.zeros((5, im0.shape[0], im0.shape[1]))
# counters
Ndark = buck1_sum = buck2_sum = buck3_sum = 0
# list of the index of the dark images (for later usage)
dark_list = []
# list of pixel values (to see the light distribution for some pixels)
some_pixel_values = []
# main loop
for i in range(N):
im = load_image(f'{npzdir}/luce_{i:04}.npz')
# store some pixel values (2 in the bucket, 2 in the dark region, 2 in the main beam region)
py1, px1 = 399, 109 # bucket 1
py2, px2 = 327, 118 # bucket 2
py3, px3 = 63, 165 # dark 1
py4, px4 = 710, 178 # dark 2
py5, px5 = 220, 575 # light 1
py6, px6 = 499, 804 # light 2
some_pixel_values.append([ im[py1,px1], im[py2,px2], im[py3,px3], im[py4,px4], im[py5,px5], im[py6,px6] ])
# dark check
is_dark = False
if im[200:500, 500:800].sum() < 200000 :
dark = dark + im
Ndark += 1 # counter of the dark images
is_dark = True
dark_list.append(i)
# buckets
buck1 = im[400,110] # single pixel value
buck2 = im[395:405, 105:115].sum() # sum of 10x10 pixels
buck3 = im[by1:by2, bx1:bx2].sum() # entire bucket region
if not is_dark:
buck1_sum += buck1
buck2_sum += buck2
buck3_sum += buck3
# normal image
total = total + im
# correlation images (without background subtraction)
total_corr1 = total_corr1 + im * buck1
total_corr2 = total_corr2 + im * buck2
total_corr3 = total_corr3 + im * buck3
if i%100 == 0 :
# saving results to a file (compute the average) every 100 images
# it can be useful to observe how the final results is growing as
# the statistic increases.
# Warning -> a lot of data! Increase 100 to 200 or more to save space
# or put to 10000 to save just 1 final file
print(f'{i} processed images ({Ndark} dark)')
if Ndark>0:
np.savez_compressed(f'results/results{i}.npz',
normal=total / (i-Ndark),
dark=dark/Ndark,
dark_list = np.array(dark_list),
ghost1=total_corr1/buck1_sum,
ghost2=total_corr2/buck2_sum,
ghost3=total_corr3/buck3_sum,
counts = np.array([i,Ndark]),
some_pixel_values = np.array(some_pixel_values)
)
else:
print('No dark images. Not saving results.')
doanalysis(201) # put 10001 to process everything
Results¶
In [152]:
f = np.load('resultsKEEP/results10000.npz')
normal, ghost1, ghost2, ghost3, dark, Ntot, Ndark= f['normal'], f['ghost1'], f['ghost2'], f['ghost3'], f['dark'], f['counts'][0], f['counts'][1]
fig, ax = plt.subplots(2,3, figsize=(24,10))
ax[0][0].set_title(f'normal {Ntot-Ndark} images')
h = ax[0][0].imshow(normal, cmap='gray')
plt.colorbar(h, ax = ax[0][0])
ax[0][1].set_title(f'dark ({Ndark} images)')
h = ax[0][1].imshow(dark, cmap='gray')
plt.colorbar(h, ax = ax[0][1])
ax[0][2].set_title('normal-dark')
h = ax[0][2].imshow(normal-dark, cmap='gray', vmin=0, vmax = normal.max())
plt.colorbar(h, ax = ax[0][2])
ax[1][0].set_title('ghost1/normal (bucket 1 pixel)')
h = ax[1][0].imshow( ghost1/normal, cmap='gray')
plt.colorbar(h, ax = ax[1][0])
ax[1][1].set_title('ghost2/normal (bucket 10x10 pixels)')
h = ax[1][1].imshow( ghost2 / normal, cmap='gray')
plt.colorbar(h, ax = ax[1][1])
ax[1][2].set_title('ghost3 (big bucket)')
h = ax[1][2].imshow( ghost3/normal, cmap='gray')
plt.colorbar(h, ax = ax[1][2])
plt.show()
Distribution of the light value on the 6 chosen pixels¶
In [153]:
f = np.load('resultsKEEP/results10000.npz')
a = f['some_pixel_values']
fig, ax = plt.subplots(2,3, figsize=(15,7))
titles = ['bucket 1', 'bucket 2', 'dark 1', 'dark 2', 'light 1', 'light 2']
for i in range(6):
ir, ic = i%2, int(i/2)
ax[ir][ic].hist(a[:,i],40, log=True, range=(0,40), color='gray')
ax[ir][ic].set_title(f'{titles[i]}')
plt.show()
Saving some pngs just for creating a gif (using imagemagick)¶
In [ ]:
## Before: create a png directory
def savepng(N) :
filename = f'resultsKEEP/results{N}.npz'
f = np.load(filename)
normal, ghost2 = f['normal'], f['ghost2']
fig, ax = plt.subplots(1,1, figsize=(8,5))
ax.set_title(f'Nr of images = {N}')
h = ax.imshow( ghost2 / normal, cmap='gray')
plt.colorbar(h, ax = ax)
fig.savefig(f'png/png{N}.png')
plt.close(fig)
for i in range(100,10100,100) :
savepng(i)
After creating the pngs the animation can be generated using:
convert -delay 15 -loop 0 $(awk 'BEGIN{for(i=100;i<10100;i+=100)printf("png/png"i".png ")}') anim.gif
Resulting in:
