In [1]:

%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import pandas as pd
import os
import scipy.optimize as so
os.chdir('/Users/q6600sl/Documents/Work/vial_detector')
from vial_detector import *

1, Load the video¶

The detection is done using trackpy package, which contains its own video reading method based on PyAV. The problem is that you will be stucked with standard framerate (25fps).

Therfore, vial_detector currently reads video independently from trackpy and pyAV

In [2]:

D = detector('./examples/Test_data.h264')

In [3]:

print 'Video Dimension:', D.img_dim
print 'Video Framerate:', D.frame_rate

Video Dimension: (1280, 960)
Video Framerate: 25

This is a very noisy video, not something you should have done in the first place. It is a wise idea to get the photography part as clean and consistent in the first place, which will save a lot of effort later on. How well will the package perform on a low quality video?

In [4]:

#Show a frame of the video
plt.imshow(D.img_stack[0])

Out[4]:

<matplotlib.image.AxesImage at 0x114d7d110>

2, Set the region of interest (ROI)¶

Of course it is easier to do this in a GUI. See the GUI example (in the /example folder)

In [5]:

D.show_ROI((500, 740), (370, 540))

In [6]:

D.set_ROI((500, 740), (370, 540))

3, Choose a set of early frames to caluclate the background image¶

In [7]:

D.show_30frames()

In [8]:

D.set_backgrd((0,5))

In [9]:

plt.imshow(D.bkgrd_image, cmap=cm.Greys_r)

Out[9]:

<matplotlib.image.AxesImage at 0x12044a490>

If a particle has not moved relative to its original position in the background image, either the particle is not a fly (such as background feature), or indeed a fly that has not moved. These particle will have low intensity in the background substracted image but high intensity in the original image. On the other hand, a moving fly will present itself as a dark particle in both the orginial image and the background substracted image. Therefore, the flies that have moved from the initial position can be distinguished from backgound features and non-moving flies.

In [10]:

fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(13, 6))
ax1.imshow(D.crtd_image[100], cmap=cm.Greys_r)
ax2.imshow(D.img_stack[100], cmap=cm.Greys_r)
ax1.set_title('Background substracted image')
ax2.set_title('Original image')

Out[10]:

<matplotlib.text.Text at 0x115b10bd0>

4, Detection¶

In [11]:

df = D.bckgrd_pdetect(thld     = 125, 
                      diameter = 5, 
                      maxsize  = 7,
                      minmass  = 300,
                      invert   = True)

Analyzing frame 149

In [12]:

df[1].head() #result is a list of pandas.DataFrame, one for each frame of video

Out[12]:

	x	y	mass	size	ecc	signal	ep	ST_dsty	True_particle	Timestamp	frame
0	111.726225	2.250720	347	1.303951	0.070628	113.494304	0.039308	0	False	0.04	1
1	218.016620	2.304709	361	1.392501	0.275820	113.494304	0.037043	6	False	0.04	1
2	224.869779	2.965602	407	1.381699	0.155834	113.494304	0.037293	20	False	0.04	1
3	178.995485	3.916479	443	1.408616	0.207425	113.494304	0.036683	8	False	0.04	1
4	192.909385	3.987055	309	1.395787	0.060083	113.494304	0.036968	1	False	0.04	1

use .diagnostic_plot method to see how good the tracking works and adjust detector parameter accordingally. Red dots are moving flies and cross dots are all the detected particles. incread thld value decrease the number of Red dots (more stringent). diameter, maxsize control size of particles. minimass determines the minimal density cutoff, and these 3 pareameter control the stringency of the particle detection process. invert = True means to detect dark dots, invert = False on the other hand means to detect bright dots.

In [13]:

fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(13, 6))
D.diagnostic_plot(10, ax=ax1)
D.diagnostic_plot(140, ax=ax2)
ax1.set_title('Frame #11')
ax2.set_title('Frame #141')

Out[13]:

<matplotlib.text.Text at 0x120e84110>

5, Plotting the mean vertical distance and the no. of flies detected for each frame¶

In [14]:

big_df = pd.concat(df)
big_df = big_df[big_df.True_particle]
grp_by = big_df.groupby('Timestamp')

In [15]:

fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(13, 4))
grp_by.y.mean().plot(ax=ax1)
grp_by.y.count().plot(ax=ax2)
ax1.set_ylabel('Mean vertical position')
ax2.set_ylabel('Count of flies that have moved')

Out[15]:

<matplotlib.text.Text at 0x11ed2b4d0>

6, Fit a simple model¶

A simple model: assuming there are total nuber of $n$ flies, and the probobably of moving away from its initial postion is $r$ per second, at $t$ time, the number of flies that have moved is $f(t)$

$f(t) = n(1-(1-r)^t)$

Ideally the model should be fit with MLE, but I will show a simple least mean square fit here to demonstrate the idea

In [16]:

F = lambda t, r, n: n*(1-(1-r)**t)
count_df = grp_by.y.count()
p, vcov  = so.curve_fit(F,
                        count_df.index.values,
                        count_df, p0=(0.25, 25.))
print p
print vcov

[  0.84246412  16.23739641]
[[ 0.00098015 -0.0042894 ]
 [-0.0042894   0.07968434]]

In [17]:

_x = np.linspace(0, 6)
plt.plot(_x, F(_x, *p), label='Model Fit')
grp_by.y.count().plot(ax=plt.gca(), label='Observed Data')
plt.ylabel('Count of flies that have moved')
plt.legend(loc=4)

Out[17]:

<matplotlib.legend.Legend at 0x11ed6e790>

Therefore, the flies have a probablity of 0.84 to move every seconds.

7, Tracking the particles¶

Tracking method thats at least two parameters: search_range and length_cutoff. 1, Controls the max. distance where two particles will be joined in a track. 2, Minimal track length in terms of number of frames. The resultant is a DataFrame containing the tracks.

In [18]:

tracks = D.particle_tracking(5, 25, memory=5)

Frame 149: 20 trajectories present

In [19]:

ax = plt.subplot(111)
D.diagnostic_plot(-1, ax=ax)
tp.plot_traj(tracks, ax=ax)
ax.set_xlim(0, D.img_stack[-1].shape[-1])
ax.set_ylim(0, D.img_stack[-1].shape[0])
ax.invert_yaxis()
plt.title('Tracks of the particles')

Out[19]:

<matplotlib.text.Text at 0x11fec5110>

In [20]:

tracks.head()

Out[20]:

	x	y	mass	size	ecc	signal	ep	ST_dsty	True_particle	Timestamp	frame	particle
frame
6	9.859101	155.828637	1313	1.345212	0.137353	121.245556	0.040370	152	True	0.24	6	6
7	10.002297	155.164625	1306	1.360119	0.126564	122.608873	0.038575	168	True	0.28	7	6
8	9.061295	155.018595	1452	1.367671	0.132276	124.415071	0.038000	173	True	0.32	8	6
9	9.862577	155.101840	1630	1.364313	0.153019	125.492175	0.037421	167	True	0.36	9	6
10	9.060841	154.946350	1808	1.367488	0.115033	126.813497	0.036893	192	True	0.40	10	6

1, Load the video¶

2, Set the region of interest (ROI)¶

3, Choose a set of early frames to caluclate the background image¶

4, Detection¶

5, Plotting the mean vertical distance and the no. of flies detected for each frame¶

6, Fit a simple model¶

7, Tracking the particles¶