--- a/python/moving.py	Tue Jul 16 01:36:50 2013 -0400
+++ b/python/moving.py	Tue Jul 16 17:00:17 2013 -0400
@@ -7,7 +7,13 @@
 from math import sqrt
 from numpy import median
 
-#from shapely.geometry import Polygon
+try:
+    from shapely.geometry import Polygon, Point as shapelyPoint
+    from shapely.prepared import prep
+    shapelyAvailable = True
+except ImportError:
+    print('Shapely library could not be loaded')
+    shapelyAvailable = False
 
 __metaclass__ = type
 
@@ -197,31 +203,35 @@
     def asint(self):
         return Point(int(self.x), int(self.y))
 
+    if shapelyAvailable:
+        def asShapely(self):
+            return shapelyPoint(self.x, self.y)
+
     def project(self, homography):
         from numpy.core.multiarray import array
         projected = cvutils.projectArray(homography, array([[self.x], [self.y]]))
         return Point(projected[0], projected[1])
 
-    def inPolygon(self, poly):
-        '''Returns if the point x, y is inside the polygon.
-        The polygon is defined by the ordered list of points in poly
-        
+    def inPolygonNoShapely(self, polygon):
+        '''Indicates if the point x, y is inside the polygon
+        (array of Nx2 coordinates of the polygon vertices)
+
         taken from http://www.ariel.com.au/a/python-point-int-poly.html
 
-        Use points_inside_poly from matplotlib.nxutils'''
+        Use Polygon.contains if Shapely is installed'''
 
-        n = len(poly);
+        n = polygon.shape[0];
         counter = 0;
 
-        p1 = poly[0];
+        p1 = polygon[0,:];
         for i in range(n+1):
-            p2 = poly[i % n];
-            if self.y > min(p1.y,p2.y):
-                if self.y <= max(p1.y,p2.y):
-                    if self.x <= max(p1.x,p2.x):
-                        if p1.y != p2.y:
-                            xinters = (self.y-p1.y)*(p2.x-p1.x)/(p2.y-p1.y)+p1.x;
-                        if p1.x == p2.x or self.x <= xinters:
+            p2 = polygon[i % n,:];
+            if self.y > min(p1[1],p2[1]):
+                if self.y <= max(p1[1],p2[1]):
+                    if self.x <= max(p1[0],p2[0]):
+                        if p1[1] != p2[1]:
+                            xinters = (self.y-p1[1])*(p2[0]-p1[0])/(p2[1]-p1[1])+p1[0];
+                        if p1[0] == p2[0] or self.x <= xinters:
                             counter+=1;
             p1=p2
         return (counter%2 == 1);
@@ -245,6 +255,13 @@
         from matplotlib.pyplot import scatter
         scatter([p.x for p in points],[p.y for p in points], **kwargs)
 
+if shapelyAvailable:
+    def pointsInPolygon(points, polygon):
+        '''Optimized tests of a series of points within (Shapely) polygon '''
+        prepared_polygon = prep(polygon)
+        return filter(prepared_polygon.contains, points)
+
+
 class NormAngle(object):
     '''Alternate encoding of a point, by its norm and orientation'''
 
@@ -525,19 +542,24 @@
         else:
             return None
     
-    def getTrajectoryInPolygon(self, polygon):
-        '''Returns the set of points inside the polygon
+    def getTrajectoryInPolygonNoShapely(self, polygon):
+        '''Returns the trajectory built with the set of points inside the polygon
         (array of Nx2 coordinates of the polygon vertices)'''
-        import matplotlib.nxutils as nx
         traj = Trajectory()
-        result = nx.points_inside_poly(self.asArray().T, polygon)
-        for i in xrange(self.length()):
-            if result[i]:
-                traj.addPositionXY(self.positions[0][i], self.positions[1][i])
+        for p in self:
+            if p.inPolygonNoShapely(polygon):
+                traj.addPosition(p)
         return traj
 
-    # version 2: use shapely polygon contains
-
+    if shapelyAvailable:
+        def getTrajectoryInPolygon(self, polygon):
+            '''Returns the trajectory built with the set of points inside the (shapely) polygon'''
+            traj = Trajectory()
+            points = [p.asShapely() for p in self]
+            for p in pointsInPolygon(points, polygon):
+                    traj.addPositionXY(p.x, p.y)
+            return traj
+    
     @staticmethod
     def lcss(t1, t2, lcss):
         return lcss.compute(t1, t2)