In [1]:

    

# Before running this notebook, you need to specify the number of engines (cpus) to use.
# for example, I ran the following at the terminal:
#  ipcluster start -n 7

import matplotlib as mpl
import matplotlib.pyplot as plt

import numpy as np
import pandas as pd

import pims
import trackpy as tp
import os
import time

%pylab inline


    

    




Populating the interactive namespace from numpy and matplotlib

In [2]:

    

from __future__ import division  # this makes mathematical division work better

myhome = '/home/viva'
datafolder = 'group/viva/Data'
data_date = '2016-01-26'
movienumber = '02'
moviename = '02_60x1.0_Hematite_swimmers_with_Lumencor_violet_100%'
filename = os.path.join(myhome, datafolder, data_date, moviename + '.avi')

darkcount_filename = os.path.join(myhome, datafolder, data_date, 'darkcount.avi')

filename


    

    
Out[2]:





'/home/viva/group/viva/Data/2016-01-26/02_60x1.0_Hematite_swimmers_with_Lumencor_violet_100%.avi'

In [8]:

    

frames = pims.Video(filename, as_grey=True)
#darkframes = pims.Video(darkcount_filename, as_grey=True)

# MacGruberscope 60x1.0
scaling = 150.0/1182.0  # microns per pixel, measured 2016-02-02
fps = 12.007 
frametime = 1000/fps    # milliseconds

bg_flag = False;   # the movie has not been backgrounded yet


    

In [4]:

    

import datetime
today = datetime.date.today().isoformat()

myanalysisfolder = 'group/viva/Analysis'
thismovieanalysisfolder = os.path.join(myhome, 
                                 myanalysisfolder, 
                                 today,
                                 'data_taken_' + data_date + ',_movie' + moviename)

thismovieanalysisfolder


    

    
Out[4]:





'/home/viva/group/viva/Analysis/2016-03-01/data_taken_2016-01-26,_movie02_60x1.0_Hematite_swimmers_with_Lumencor_violet_100%'

In [5]:

    

import holopy as hp
hp.show(rawframes[100])
plt.title('A frame from the movie\n' + moviename +'\n')


    

    
Out[5]:





<matplotlib.text.Text at 0x7f05689fd350>






    

In [6]:

    

rawframes[0].shape


    

    
Out[6]:





(1024, 1280)

In [12]:

    

# find bright spots in a frame.
# featuresize must be odd.
# read up on this in the trackpy literature.
i = 50
featuresize = 13
plt.title('frame ' + str(i) + ', seeking spots of diameter ' 
          + str(featuresize) + ' pixels\n' + moviename +'\n')
f1 = tp.locate(frames[i], diameter=featuresize, invert=True, minmass=400)
tp.annotate(f1, frames[i])
tp.subpx_bias(f1)
#f1


    

    

















    
Out[12]:





array([[<matplotlib.axes.AxesSubplot object at 0x7f0566d87490>,
        <matplotlib.axes.AxesSubplot object at 0x7f056692b410>]], dtype=object)






    

In [13]:

    

# Following the recommendations of the trackpy walkthrough, I am using ipyparallel.

from ipyparallel import Client

client = Client()
view = client.load_balanced_view()

client.ids


    

    
Out[13]:





[0, 1, 2, 3, 4, 5, 6]

In [18]:

    

curried_locate = lambda image: tp.locate(image, featuresize, invert=True, minmass=400)


    

In [19]:

    

%%px

import trackpy as tp

# The above will import trackpy again, but now with parallel processing.
# https://ipython.org/ipython-doc/3/parallel/magics.html


    

In [20]:

    

amr = view.map_async(curried_locate, frames[:100]) # take out the 100 if you want all frames
amr.wait_interactive()


    

    




 100/100 tasks finished after    0 s
done

In [21]:

    

time1 = time.time()
serial_result = list(map(curried_locate, frames[:100]))  # take out the 100 if you want all frames
elapsed = time.time() - time1
f_coords = pd.concat(serial_result, axis=0, join='outer')

print 'Found all coordinates in ' + str(elapsed/60.0) + ' minutes.'


    

    




Found all coordinates in 2.24124991496 minutes.

In [30]:

    

serial_result[0].head()


    

    
Out[30]:






  

    

      

      
x
      
y
      
mass
      
size
      
ecc
      
signal
      
raw_mass
      
ep
      
frame
    
  
  

    

      
0
      
 661.857917
      
 14.143296
      
 3756.530743
      
 3.838871
      
 0.010585
      
 47.453238
      
 17367
      
 0.030774
      
 0
    
    

      
1
      
 933.619733
      
 18.730041
      
 4014.488751
      
 3.683934
      
 0.009914
      
 58.764766
      
 18108
      
 0.026133
      
 0
    
    

      
2
      
 281.656695
      
 25.215503
      
 3883.164674
      
 3.671493
      
 0.058993
      
 59.316547
      
 17477
      
 0.029984
      
 0
    
    

      
3
      
 717.012658
      
 39.087101
      
 3731.148778
      
 3.730961
      
 0.023879
      
 54.626402
      
 18237
      
 0.025465
      
 0
    
    

      
4
      
 383.427110
      
 43.957255
      
 3010.797591
      
 3.589836
      
 0.080933
      
 51.591602
      
 16473
      
 0.039165
      
 0
    
  

5 rows × 9 columns

In [22]:

    

# tell me how many frames are in the movie
nframes = f_coords['frame'].max() - f_coords['frame'].min() + 1
nframes
# This will error out if your movie wasn't a pims movie 
# because tp.locate only puts the frame column on a pims frame.


    

    
Out[22]:





100

In [23]:

    

# We have just built a list of coordinates called f_coords where we have seen particles. '
# Now we want to link these together from one frame to the next 
# so we can identify the trajectory for each particle.

# Documentation: http://soft-matter.github.io/trackpy/generated/trackpy.link_df.html

t = tp.link_df(features=f_coords, search_range=10, memory=3)
#t = pd.read_pickle('t.pkl')

# search_range gives the maximum distance features can move between frames. 
#              I think it's measured in pixels.
# memory gives the maximum number of frames during which a feature can vanish, 
#        then reappear nearby, and still be considered the same particle.
# This will run faster if the numba package is available.


    

    




Frame 99: 65 trajectories present

In [24]:

    

trajectory_plot = tp.plot_traj(t, superimpose = frames[0], label=False)


    

In [25]:

    

# only keep trajectories that last at least this many frames
t1 = tp.filter_stubs(t, 50)
# Compare the number of particles in the unfiltered and filtered data.
print 'Before:', t['particle'].nunique()
print 'After:', t1['particle'].nunique()


    

    




Before: 109
After: 66

In [26]:

    

tracks = t1['particle'].astype(int).unique()

print size(tracks)


    

    




66

In [27]:

    

trajectory_plot = tp.plot_traj(t1, superimpose = frames[0], label=False)


    

In [28]:

    

try:
    axes().set_aspect('equal', 'datalim') 
except:
    pass
trajectory_plot = tp.plot_traj(t1, mpp=scaling)

#savefig()


    

In [29]:

    

d = tp.compute_drift(t1, smoothing=15)

#plt.figure()
d.plot(grid=False)
plt.title('Drift in ' + moviename + '\n')


    

    
Out[29]:





<matplotlib.text.Text at 0x7f0565d4dd10>






    

In [ ]:

    




    

In [ ]: