In [1]:
import trackpy as tp
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import math
import seaborn
%pylab inline
In [2]:
traj_file = '/home/viva/Data/testing_code/filtered_data_tracer+janus_3%_H2O2_13(green)2016-06-14_pickled.pkl'
In [3]:
t1 = pd.read_pickle(traj_file)
In [4]:
tp.subpx_bias(t1)
## Note that there is subpixel bias so this might not be the best tracking.
Out[4]:
In [5]:
## concatenate a new numerical column to a matrix
def put_z_position_in_matrix(mat2D, z=0):
z_position = np.zeros(len(mat2D)) + z
z_position = np.matrix(z_position)
mat3D = np.concatenate((mat2D.T, z_position))
mat3D = mat3D.T
return mat3D
## Check to see if dataframe has z column; otherwise assume z=0.
def get_3D_matrix_from_dataframe(df, xlabel='x',ylabel='y',zlabel='z'):
try:
matrix = np.mat(df[[xlabel,ylabel,zlabel]])
except KeyError:
matrix = np.mat(df[[xlabel,ylabel]])
matrix = put_z_position_in_matrix(matrix,0)
return matrix
## The variable A_3D will be a matrix consisting of
## all coordinates in frame i
## whose particle is also tracked in frame f.
## The variable B_3D will be a matrix consisting of
## all coordinates in frame i
## whose particle is also tracked in frame f.
## This function currently assumes the particles tracked in the image frame
## are all at the same z.
def matrices_from_dataframe(t1, framei, framef=None, z=0):
# set default for final frame
if framef == None:
framef = framei+1
# an inner merge will drop any rows for
# particles that are not in both frames
AB = pd.merge(t1[t1['frame'] == framei],
t1[t1['frame'] == framef],
how='inner',
on='particle',
suffixes=('_i','_f'))
# Pull out the coordinates and convert to matrices.
# If z positions are not available, they are set to zero.
A_3D = get_3D_matrix_from_dataframe(AB, xlabel='x_i',ylabel='y_i',zlabel='z_i')
B_3D = get_3D_matrix_from_dataframe(AB, xlabel='x_f',ylabel='y_f',zlabel='z_f')
assert len(A_3D) == len(B_3D)
return A_3D, B_3D
## Given a matrix B which
## has experienced rotation R and translation t,
## undo that transformation.
def rotational_drift_subtraction(B, R, t):
n = len(B)
drift_subtracted = R.T * (B.T - np.tile(t,(1,n)))
drift_subtracted = drift_subtracted.T
return drift_subtracted
## This function is copied from http://nghiaho.com/uploads/code/rigid_transform_3D.py_
# Input: expects Nx3 matrix of points
# Returns R,t
# R = 3x3 rotation matrix
# t = 3x1 column vector
def rigid_transform_3D(A, B):
assert len(A) == len(B)
N = A.shape[0]; # total points
centroid_A = np.mean(A, axis=0)
centroid_B = np.mean(B, axis=0)
# centre the points
AA = A - np.tile(centroid_A, (N, 1))
BB = B - np.tile(centroid_B, (N, 1))
# dot is matrix multiplication for array
H = np.transpose(AA) * BB
U, S, Vt = np.linalg.svd(H)
R = Vt.T * U.T
# special reflection case
if np.linalg.det(R) < 0:
print "Reflection detected"
Vt[2,:] *= -1
R = Vt.T * U.T
t = -R*centroid_A.T + centroid_B.T
#print t
return R, t
# Calculate the axis and angle of rotation for a given rotation matrix R
def axis_angle(R):
h = R[2,1]
f = R[1,2]
c = R[0,2]
g = R[2,0]
d = R[1,0]
b = R[0,1]
# axis of rotation
axis = [h-f, c-g, d-b]
# angle of rotation, in radians
angle = np.arccos((np.trace(R) - 1)/2)
## different way to calculate angle
# axis_length = np.linalg.norm(axis)
# angle = np.arcsin(axis_length/2)
return np.mat(axis), angle # in radians
def print_head(matrix, max_printable_length = 10):
if len(matrix)>max_printable_length:
print matrix[0:max_printable_length]
print "..."
else:
print matrix
In [6]:
t1.head()
Out[6]:
In [7]:
t1.frame.min()
Out[7]:
In [8]:
t1.frame.max()
Out[8]:
In [9]:
len(t1.frame.unique())
Out[9]:
In [10]:
## How many frames to drift-subtract?
max_frame = t1.frame.max()
max_frame = 200 # for debugging
# Set this to false if you don't care about calculating tm
# and just want to plot the drift.
# Drift subtraction is slow; there's a loop inside a loop.
do_drift_subtraction = True;
## Initialize loop
if do_drift_subtraction:
tm=0 # overwrite tm
del tm # erase tm
prev_frame = None;
R_list = []
t_list = []
x_drifts = []
y_drifts = []
z_drifts = []
axis_list = []
angle_list = []
frame_list = []
verbose = False
labelx = 'x'
labely = 'y'
labelz = 'z'
transformed_str = '_drift_subtracted'
labelx2 = labelx + transformed_str
labely2 = labely + transformed_str
labelz2 = labelz + transformed_str
labelnote = 'relative_to_frame'
for current_frame in sort(t1.frame.unique()):
if current_frame > max_frame:
break;
if verbose:
print "Frame ", current_frame
if prev_frame is None:
relative_to = current_frame;
prev_frame = current_frame;
continue; # skip first frame
assert prev_frame is not None
# A is a shorthand for the previous frame.
# B is a shorthand for the current frame.
# Get raw coordinates from current frame and previous frame
A_3D, B_3D = matrices_from_dataframe(t1, prev_frame, current_frame)
# Figure out the transformation that occured between frames
ret_R, ret_t = rigid_transform_3D(A_3D, B_3D);
# Save a copy of the transformation
R_list.append(ret_R)
t_list.append(ret_t)
x_drifts.append(np.array(ret_t)[0][0])
y_drifts.append(np.array(ret_t)[1][0])
z_drifts.append(np.array(ret_t)[2][0])
current_axis,current_angle = axis_angle(ret_R)
axis_list.append(current_axis)
angle_list.append(current_angle)
frame_list.append(current_frame)
if do_drift_subtraction:
## Do the rotational drift subtraction.
## I need to do this with all particles in current frame,
## not just the ones that also appear in previous frame.
B_dataframe = t1[t1['frame'] == current_frame].copy()
B = get_3D_matrix_from_dataframe(B_dataframe)
for R,t in zip(reversed(R_list),reversed(t_list)):
if verbose:
print "undoing transformation"
print R
B = rotational_drift_subtraction(B, R, t)
# This is rather brute force,
# but I wanted to make sure I'm correct first.
# The better thing to do is probably to calculate
# the total transformation before transforming the coordinates.
## Record the drift-subtracted coordinates
# (i.e. Put the transformed data in the dataframe)
x_sub_data = np.array(B[:,0]).T[0]
y_sub_data = np.array(B[:,1]).T[0]
z_sub_data = np.array(B[:,2]).T[0]
B_dataframe[labelx2]=x_sub_data
B_dataframe[labely2]=y_sub_data
if not np.array_equal(z_sub_data, np.zeros_like(z_sub_data)):
## Not tested with a z column
B_dataframe[labelz2]=z_sub_data
num_new_cols = 4
else:
## no z data
num_new_cols = 3
B_dataframe[labelnote] = relative_to
try:
tm = pd.concat([tm, B_dataframe])
except NameError:
# Initialize tm
tm = B_dataframe.copy()
prev_frame = current_frame;
# end loop
## Rename some columns in tm
if do_drift_subtraction:
# Put the new columns up front
cols = tm.columns.tolist()
cols = cols[-num_new_cols:]+cols[:-num_new_cols]
tm = tm.reindex(columns=cols);
## Rename raw columns
tm = tm.rename(index=str,
columns={labelx: labelx + "_raw",
labely: labely + "_raw"});
tm = tm.rename(index=str,
columns={labelx2: labelx,
labely2: labely});
if num_new_cols == 4:
## Not tested with a z column
tm = tm.rename(index=str,
columns={labelz: labelz + "_raw"});
tm = tm.rename(index=str,
columns={labelz2: labelz});
In [11]:
tm.tail()
Out[11]:
In [12]:
len(t1.frame[0:max_frame])
Out[12]:
In [13]:
len(frame_list)
Out[13]:
In [14]:
plt.plot(frame_list,angle_list)
plt.title(traj_file + '\nAngular drift\n')
plt.xlabel('Frame')
plt.ylabel('Angular drift [radians]')
Out[14]:
In [15]:
plt.plot(frame_list,x_drifts, label="x")
plt.plot(frame_list,y_drifts, label="y")
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.title(traj_file + '\nTranslational drift\n')
plt.xlabel('Frame')
plt.ylabel('Translational drift [pixels]')
Out[15]:
In [16]:
pylab.axis('equal')
plt.title(traj_file + '\nRaw trajectories\n')
ax = tp.plot_traj(t1[t1['frame'] < max_frame], legend=False)
In [17]:
pylab.axis('equal')
plt.title(traj_file + '\nTrajectories with rotational drift subtraction\n')
ax = tp.plot_traj(tm, legend=False)
In [ ]: