Step Linking
Back to 2014.
This python script reads in the lists of particles for each image generated by the python image processing step and outputs a list of steps in "steps.txt" and also images of the particle associations, in pixel units, in each step between figures in "match*.png". In version 1.1, the match.png has color coding to show which steps are marked as outliers, to enable checking of whether they are misindentified steps or actual steps. Also the associations are labeled so you can match them to the steps.txt list.
Version 1.1 (Some improvements and bug fixes pertaining to particular cases of user input.)
# -*- coding: utf-8 -*-
import numpy as np
import pylab as plt
from glob import glob
#Adjustable parameters:
image_width=1600
image_height=1200
min_pix = 40 #minimum number of pixels a particle must have
max_dist = 60 #particle association radius (pixel units)
max_pix_diff = 10 #maximum difference in pixel count in particle association
max_num_steps = 120 #maximum number of steps to find (b/f removing outliers)
outlier_threshold = 35 #maximum dx or dy allowed for step consult ...
# xy_step_histogram.png & variance calculation. ...
# For ~100 steps a 3.5 sigma cut seems reasonable.
#Read in the particle locations and pixel count for all images
filelist = glob('*.txt')
particle_lists=[] #list of lists for all images
for filename in filelist:
particle_lists.append(np.loadtxt(filename))
dx = []
dy = []
nparticles=0
for i in range(len(particle_lists)-1):
matches = []
print "Processing steps between %s and %s." % (filelist[i], filelist[i+1])
for p1 in particle_lists[i]:
displacements = []
diffs = []
p2s = []
if p1[2] < min_pix: #Skip if pixels in p1 are too few
continue
for p2 in particle_lists[i+1]:
if p2[2] <= min_pix: #Skip if pixels in p2 are too few
continue
dist = np.sqrt((p1[0]-p2[0])**2 + (p1[1]-p2[1])**2)
if dist <= max_dist: # Are p1 and p2 associated?
displacements.append([(p1[0]-p2[0]),(p1[1]-p2[1])])
diffs.append(abs(p1[2]-p2[2]))
p2s.append(p2)
#Veto p1 if it has no associations or if it has more than one
if len(displacements)==0 or len(displacements) > 1:
#print displacements
continue
if diffs[0] > max_pix_diff:
continue
#If the particle 2 has already been matched to another p1 veto both p1s
double_match=False
for k in range(len(matches)):
if (p2s[0][0] == matches[k][0]) and (p2s[0][1] == matches[k][1]):
matches.pop(k)
double_match=True
break
if double_match: continue
#Otherwise we've found a good association! Record the displacements.
matches.append(p2s[0])
dy.append(displacements[0][0])
dx.append(displacements[0][1])
#Plot the association in a scatter plot "matches*.png"
dist = np.sqrt(displacements[0][0]**2 + displacements[0][1]**2)
if dist < outlier_threshold:
nparticles+=1
plt.scatter([p1[1], (p1[1]-dx[-1])], [image_height- p1[0],image_height-(p1[0]-dy[-1])], color='blue')
plt.plot([p1[1], (p1[1]-dx[-1])], [image_height- p1[0],image_height-(p1[0]-dy[-1])], color='blue')
plt.text(p1[1]+5, image_height- p1[0]+5, "%d" % nparticles)
else: #Plot outliers but don't number them b/c we want to reference the output list from the figure
plt.scatter([p1[1], (p1[1]-dx[-1])], [image_height- p1[0],image_height-(p1[0]-dy[-1])], color='red')
plt.plot([p1[1], (p1[1]-dx[-1])], [image_height- p1[0],image_height-(p1[0]-dy[-1])], color='red')
plt.xlim([0,image_width])
plt.ylim([0,image_height])
if len(dx)>max_num_steps: #Stop after we have 100
break
plt.savefig('matches_%s_%s.png' % (filelist[i].split('.')[0],filelist[i+1].split('.')[0]))
plt.clf()
if len(dx)>max_num_steps: #Stop after we have 100
break
#Convert from list type to numpy array and calibrate
dx = np.array(dx)
dy = np.array(dy)
print "Standard deviation in dx:", dx.std()
print "Standard deviation in dy:", dy.std()
#Plot histograms
plt.hist(dx, bins=15, range=[-max_dist,max_dist], histtype='step', color='blue', label=r'$\Delta$x')
plt.hist(dy, bins=15, range=[-max_dist,max_dist], histtype='step', color='red', label=r'$\Delta$y')
plt.xlabel('Distance (pixels)')
plt.legend()
plt.savefig('xy_step_histogram.png')
plt.clf()
#Remove outlier steps: require dx AND dy be less than outlier threshold
inds = np.where((abs(dx)<outlier_threshold)*(abs(dy)<outlier_threshold))
dx = dx[inds]
dy = dy[inds]
print "%d steps remain after outlier cut." % len(dx)
print "New standard deviation in dx:", dx.std()
print "New standard deviation in dy:", dy.std()
#Write out steps:
np.savetxt('steps.txt', np.array([dx,dy]).transpose())
Version 1.0 (Original version)
# -*- coding: utf-8 -*-
import numpy as np
import pylab as plt
from glob import glob
#Adjustable parameters:
image_width=1600
image_height=1200
min_pix = 40 #minimum number of pixels a particle must have
max_dist = 60 #particle association radius (pixel units)
max_pix_diff = 10 #maximum difference in pixel count in particle association
max_num_steps = 120 #maximum number of steps to find (b/f removing outliers)
outlier_threshold = 35 #maximum dx or dy allowed for step consult ...
# xy_step_histogram.png & variance calculation. ...
# For ~100 steps a 3.5 sigma cut seems reasonable.
#Read in the particle locations and pixel count for all images
filelist = glob('*.txt')
particle_lists=[] #list of lists for all images
for filename in filelist:
particle_lists.append(np.loadtxt(filename))
dx = []
dy = []
for i in range(len(particle_lists)-1):
print "Processing steps between %s and %s." % (filelist[i], filelist[i+1])
for p1 in particle_lists[i]:
displacements = []
if p1[2] < min_pix: #Skip if pixels in p1 are too few
continue
for p2 in particle_lists[i+1]:
if p2[2] <= min_pix: #Skip if pixels in p2 are too few
continue
if abs(p1[2]-p2[2])>max_pix_diff:
continue
dist = np.sqrt((p1[0]-p2[0])**2 + (p1[1]-p2[1])**2)
if dist <= max_dist: # Are p1 and p2 associated?
displacements.append([(p1[0]-p2[0]),(p1[1]-p2[1])])
#Veto p1 if it has no associations or if it has more than one
if len(displacements)==0 or len(displacements) > 1:
#print displacements
continue
#Otherwise we've found a good association! Record the displacements.
dy.append(displacements[0][0])
dx.append(displacements[0][1])
#Plot the association in a scatter plot "matches*.png"
plt.scatter([p1[1], (p1[1]-dx[-1])], [image_height- p1[0],image_height-(p1[0]-dy[-1])])
plt.xlim([0,image_width])
plt.ylim([0,image_height])
if len(dx)>max_num_steps: #Stop after we have 100
break
if len(dx)>max_num_steps: #Stop after we have 100
break
plt.savefig('matches_%s_%s.png' % (filelist[i].split('.')[0],filelist[i+1].split('.')[0]))
plt.clf()
plt.clf()
#Convert from list type to numpy array and calibrate
dx = np.array(dx)
dy = np.array(dy)
print "Standard deviation in dx:", dx.std()
print "Standard deviation in dy:", dy.std()
#Plot histograms
plt.hist(dx, bins=15, range=[-max_dist,max_dist], histtype='step', color='blue', label=r'$\Delta$x')
plt.hist(dy, bins=15, range=[-max_dist,max_dist], histtype='step', color='red', label=r'$\Delta$y')
plt.xlabel('Distance (pixels)')
plt.legend()
plt.savefig('xy_step_histogram.png')
plt.clf()
#Remove outlier steps: require dx AND dy be less than outlier threshold
inds = np.where((abs(dx)<outlier_threshold)*(abs(dy)<outlier_threshold))
dx = dx[inds]
dy = dy[inds]
print "%d steps remain after outlier cut." % len(dx)
#Write out steps:
np.savetxt('steps.txt', np.array([dx,dy]).transpose())