In this Notebook we analyse movies taken from a 2D particle simulator used for educational purposes. In this device, small magnetic and repulsive discs move on a flat surface due to air flowing from the sides as well as from underneath to reduce friction (see figure below). The particles are confined by repulsive, magnetic bars, located on the four sides of the square surface.
shift+return.() brackets and press shift+tab-tab.In the first part we explore the radial distribution function (RDF), $g(r)$.
Flow of program:
ffmpeg).h5 format)More here info about the particle tracker here.
In [53]:
# load modules required for the analysis
from __future__ import division, unicode_literals, print_function # for compatibility with Python 2 and 3
import matplotlib as mpl
import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np, pandas as pd
from pandas import DataFrame, Series
import os, pims, scipy, trackpy as tp
import base64
from scipy.spatial import distance
from IPython.display import HTML
plt.rcParams.update({'font.size': 16, 'figure.figsize': [8.0, 6.0]})
# this function is used to visualize videos in the notebook
def video(file, mimetype='mp4'):
""" Show given video file """
video_encoded = base64.b64encode( open(file, "r+b").read() )
return HTML(data='''<video alt="test" controls>
<source src="data:video/mp4;base64,{0}" type="video/mp4" />
</video>'''.format(video_encoded.decode('ascii')))
The following videos of particle simulations were recorded using a smart phone and these will be the basis of the following analysis.
Here we create an array containing the filenames of the movies and two parameters that will be used for the image analysis:
In [54]:
movies = [
{'file': 'movies/N20-light.mp4', 'percentile': 60, 'diameter': 53}, # dictionary for movie 0
{'file': 'movies/N40-light.mp4', 'percentile': 60, 'diameter': 53}, # dictionary for movie 1
{'file': 'movies/N55-light.mp4', 'percentile': 60, 'diameter': 53} # dictionary for movie 2
] # more movies can be added as needed...
In [55]:
movie = movies[1] # select movie file here by index, starting from 0
moviefile = movie['file'] # select key 'file' of the dictionary
Here we visualize the movie.
In [56]:
video(moviefile)
Out[56]:
In [57]:
imgdir=os.path.splitext(moviefile)[0]
if not os.path.exists(imgdir):
os.makedirs(imgdir)
!ffmpeg -i $moviefile -f image2 -vcodec mjpeg $imgdir/img-%03d.jpg -v 0
In [58]:
frames = pims.ImageSequence(imgdir+'/img*.jpg', as_grey=True) # load all frames
print('read',len(frames),'frames.')
f = tp.locate(frames[0], diameter=movie['diameter'], invert=True, percentile=movie['percentile'])
print(f.tail())
# this is how we may separate light and heavy particles
heavy = f[f['mass']>50000] # heavy particles
light = f[f['mass']<50000] # light particles
tp.annotate(f, frames[0])
plt.show()
In [59]:
trjfile=os.path.splitext(moviefile)[0]+'.h5'
if os.path.isfile(trjfile):
print('opening existing trajectory file', trjfile)
else:
with tp.PandasHDFStore(trjfile) as s:
cnt=0
for image in frames:
cnt=cnt+1
print( 'frame %d/%d.' % (cnt, len(frames)), end=' ')
features = tp.locate(
image, diameter=movie['diameter'], percentile=movie['percentile'], invert=True)
print('number of particles =',len(features))
s.put(features[['x','y','mass','frame']])
In this section we calculate all distances between all particles for each frame. These are then binned into a histogram to give the probability of observing a particular separation.
At the same time we sample the distribution for ideal particles by simply generating random positions and perform the same analysis as for the "real" particles.
In [60]:
dist = np.ndarray(shape=(0,0)) # movie distribution
udist = np.ndarray(shape=(0,0)) # uniform distribution
with tp.PandasHDFStore(trjfile) as s:
# measure all distances between discs
for frame in s:
dist = np.append( dist, distance.pdist( frame[['x','y']] ))
data = s.dump() # keep full trajectory in `data`
# find box corners
xmin,xmax = min(data.x), max(data.x)
ymin,ymax = min(data.y), max(data.y)
# measure distances for uniform, random distribution (no correlations)
# (number of points/particles is unimportant.)
x = np.random.randint(xmin, xmax+1, 4000)
y = np.random.randint(ymin, ymax+1, 4000)
p = np.array([x, y]).T
udist = np.append( udist, distance.pdist(p) )
hist = plt.hist(dist, bins=150, normed=True, range=[0,700], histtype='step', color='black', label='real')
uhist = plt.hist(udist, bins=150, normed=True, range=[0,700], histtype='step', color='red', label='ideal')
plt.legend(loc=0,frameon=False)
plt.xlabel('distance (pixels)')
plt.ylabel('probability')
We have now calculated the distance distribution, hist, from the simulated particles from the movie, as well as for a uniform distribution of $N$ particles, uhist. The radial distribution function is simply the ratio between the two.
This means that if the particles were behaving ideally (which they don't), $g(r)$ would be unity for all separations, $r$. After plotting, the final rdf is saved to disk.
In [61]:
r = hist[1][:len(hist[0])]
g = hist[0] / uhist[0]
plt.plot( r, g, 'k-')
plt.xlabel('$r$ (pixels)')
plt.ylabel('$g(r)$')
plt.title('Radial distribution function (RDF)')
rdffile=os.path.splitext(moviefile)[0]+'.rdf.dat'
np.savetxt(rdffile, np.array([r,g]).T, header='rdf from '+moviefile)
In [62]:
for d in movies:
name = os.path.splitext(d['file'])[0]
rdffile = name+'.rdf.dat'
if os.path.isfile(rdffile):
r, g = np.loadtxt(rdffile, unpack=True )
plt.plot(r, g, '-', label=os.path.basename(name), lw=2)
plt.legend(loc=0,frameon=False)
plt.xlabel('$r$')
plt.ylabel('$g(r)$')
plt.xlim([40,300])
plt.ylim([0,4])
Out[62]:
np.log() as in pmf=-np.log(g)).
In [63]:
t = tp.link_df(data, search_range=20, memory=3)
tp.plot_traj(t)
Out[63]:
In [64]:
em = tp.emsd(t, mpp=20./600., fps=24) # mpp=cm per pixel
In [65]:
plt.plot(em.index, em, 'o')
plt.xscale('log')
plt.yscale('log')
plt.ylabel(r'$\langle \Delta r^2 \rangle$ [cm$^2$]')
plt.xlabel('lag time $t$ [s]')
Out[65]:
Mikael Lund, Department of Theretical Chemistry
Lund University, POB 124, SE-22100 Lund, Sweden
http://www.teokem.lu.se