Step Linking

# -*- coding: utf-8 -*-
import numpy as np
import pylab as plt
from glob import glob

#Adjustable parameters:
min_pix = 40 #minimum number of pixels a particle must have
max_dist = 50 #particle association radius (pixel units)
max_pix_diff = 10 #maximum difference in pixels in particle association
max_num_steps = 120 #maximum number of steps to find (b/f removing outliers)
outlier_threshold = 40 #maximum dx or dy allowed for step consult xy_step_histogram.png

#read in  the particle locations for all images
filelist = glob('*.txt')
particle_lists=[] #list of 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])
        plt.scatter([p1[1], (p1[1]-dx[-1])], [1200- p1[0],1200-(p1[0]-dy[-1])])
        plt.xlim([0,1600])
        plt.ylim([0,1200])
        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=10, range=[-max_dist,max_dist], histtype='step', color='blue', label=r'$\Delta$x')
plt.hist(dy, bins=10, range=[-max_dist,max_dist], histtype='step', color='red', label=r'$\Delta$y')
plt.xlabel(r'Distance ($\mu$m)')
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())