1. Olivier Debeir
  2. ivctrack

Commits

Olivier Debeir  committed 3fcaf5b

add speed measurement

  • Participants
  • Parent commits 5b77bdd
  • Branches master
  • Tags v0.1.1

Comments (0)

Files changed (5)

File docs/example_measure.rst

View file
+Trajectory analysis
+-----------------------------
+
+Cell speed analysis : average speed, MRDO and hull speed :
+
+.. ipython::
+
+    In [28]: import matplotlib.pyplot as plt
+
+    In [32]: from ivctrack.hdf5_read import get_hdf5_data
+
+    In [32]: from ivctrack.measurement import speed_feature_extraction
+
+    In [9]: hdf5_filename = '../test/data/test_rev.hdf5'
+
+    In [29]: feat,data = speed_feature_extraction(hdf5_filename)
+
+    In [30]: print feat
+
+    In [31]: plt.scatter(data[:,1],data[:,3])
+
+    @savefig plot_measure1.png width=4in
+    In [34]: plt.draw()
+
+    In [34]: plt.figure()
+
+    In [31]: plt.hist(data[:,2:4])
+
+    @savefig plot_measure2.png width=4in
+    In [34]: plt.draw()

File docs/examples.rst

View file
     example_experiment.rst
     example_read_hdf5.rst
     example_plot.rst
+    example_measure.rst
 
 
 

File ivctrack/__init__.py

View file
 import cellmodel
 import helpers
 import meanshift
-import hdf5_read
+import hdf5_read
+import measurement

File ivctrack/measurement.py

View file
+# -*- coding: utf-8 -*-
+'''Trajectory feature extraction
+'''
+__author__ = 'Copyright (C) 2012, Olivier Debeir <odebeir@ulb.ac.be>'
+__license__ ="""
+pyrankfilter is a python module that implements 2D numpy arrays rank filters, the filter core is C-code
+compiled on the fly (with an ad-hoc kernel).
+
+Copyright (C) 2012  Olivier Debeir
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+"""
+import numpy as npy
+from hdf5_read import get_hdf5_data
+from quickhull2d import qhull
+from scipy.spatial.distance import pdist,squareform
+
+def inst_speed(xy):
+    s = npy.sqrt(npy.sum(npy.diff(xy,axis=0)**2,axis=1))
+    return s
+
+def path_length(xy):
+    return npy.sum(inst_speed(xy))
+
+def avg_speed(xy):
+    return path_length(xy)/xy.shape[0]
+
+def mrdo_speed(xy):
+    xy0 = xy[0,:]
+    d = npy.sqrt(npy.sum((xy-xy0)**2,axis=1))
+    d_max = npy.max(d)
+    return d_max/xy.shape[0]
+
+def hull_speed(xy):
+    hull_xy = qhull(xy)
+    d = pdist(hull_xy, 'euclidean')
+    dmax_idx = npy.argmax(d)
+    return d[dmax_idx]/xy.shape[0]
+
+def plot_hull_speed(xy):
+    import matplotlib.pyplot as plt
+
+    hull_xy = qhull(xy)
+    d = pdist(hull_xy, 'euclidean')
+    dmax_idx = npy.argmax(d)
+
+    r = npy.zeros_like(d)
+    r[dmax_idx] = 1
+    sd = squareform(r)
+
+    ext_idx = npy.nonzero(sd)[1]
+    p0 = hull_xy[ext_idx[0],:]
+    p1 = hull_xy[ext_idx[1],:]
+
+    fig = plt.figure()
+    ax = fig.add_subplot(111, aspect='equal')
+    plt.scatter(xy[:,0],xy[:,1])
+    plt.plot(hull_xy[:,0],hull_xy[:,1])
+    plt.plot([p0[0],p1[0]],[p0[1],p1[1]])
+    plt.show()
+    return 0
+
+def speed_feature_extraction(hdf5_filename):
+
+    features,data = get_hdf5_data(hdf5_filename,fields=['center'])
+    measures = []
+    feat_name = ['path_length','avg_speed','mrdo','hull_speed']
+    for d in data:
+        xy = d['center']
+        pl = path_length(xy)
+        avg = avg_speed(xy)
+        mrdo = mrdo_speed(xy)
+        h = hull_speed(xy)
+        measures.append((pl,avg,mrdo,h))
+    measures = npy.asarray(measures)
+
+    return (feat_name,measures)
+
+if __name__=='__main__':
+
+    import matplotlib.pyplot as plt
+
+    hdf5_filename = '../test/data/test_rev.hdf5'
+
+    feat,data = speed_feature_extraction(hdf5_filename)
+    print feat
+    print data
+
+    plt.scatter(data[:,1],data[:,3])
+    plt.show()

File ivctrack/quickhull2d.py

View file
+# Copyright (c) 2012 the authors listed at the following URL, and/or
+# the authors of referenced articles or incorporated external code:
+# http://en.literateprograms.org/Quickhull_(Python,_arrays)?action=history&offset=20091103134026
+# 
+# Permission is hereby granted, free of charge, to any person obtaining
+# a copy of this software and associated documentation files (the
+# "Software"), to deal in the Software without restriction, including
+# without limitation the rights to use, copy, modify, merge, publish,
+# distribute, sublicense, and/or sell copies of the Software, and to
+# permit persons to whom the Software is furnished to do so, subject to
+# the following conditions:
+# 
+# The above copyright notice and this permission notice shall be
+# included in all copies or substantial portions of the Software.
+# 
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+# IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
+# CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
+# TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
+# SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+# 
+# Retrieved from: http://en.literateprograms.org/Quickhull_(Python,_arrays)?oldid=16555
+
+from numpy import *
+
+link = lambda a,b: concatenate((a,b[1:]))
+edge = lambda a,b: concatenate(([a],[b]))
+def qhull(sample):
+    def dome(sample,base): 
+        h, t = base
+        dists = dot(sample-h, dot(((0,-1),(1,0)),(t-h)))
+        outer = repeat(sample, dists>0, 0)
+
+        if len(outer):
+            pivot = sample[argmax(dists)]
+            return link(dome(outer, edge(h, pivot)),
+                        dome(outer, edge(pivot, t)))
+        else:
+            return base
+    if len(sample) > 2:
+        axis = sample[:,0]
+        base = take(sample, [argmin(axis), argmax(axis)], 0)
+        return link(dome(sample, base),
+                    dome(sample, base[::-1]))
+    else:
+        return sample
+
+if __name__ == "__main__":
+    #sample = 10*array([(x,y) for x in arange(10) for y in arange(10)])
+    sample = 100*random.random((32,2))
+    hull = qhull(sample)
+	
+    print "%!"
+    print "100 500 translate 2 2 scale 0 0 moveto"
+    print "/tick {moveto 0 2 rlineto 0 -4 rlineto 0 2 rlineto"
+    print "              2 0 rlineto -4 0 rlineto 2 0 rlineto} def"
+    for (x,y) in sample:
+        print x, y, "tick"
+    print "stroke"
+    print hull[0,0], hull[0,1], "moveto"
+    for (x,y) in hull[1:]:
+        print x, y, "lineto"
+    print "closepath stroke showpage"