Mercurial Hosting > traffic-intelligence

changeset 372:349eb1e09f45

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Cleaned the methods/functions indicating if a point is in a polygon In general, shapely should be used, especially for lots of points: from shapely.geometry import Polygon, Point poly = Polygon(array([[0,0],[0,1],[1,1],[1,0]])) p = Point(0.5,0.5) poly.contains(p) -> returns True poly.contains(Point(-1,-1)) -> returns False You can convert a moving.Point to a shapely point: p = moving.Point(1,2) p.asShapely() returns the equivalent shapely point If you have several points to test, use moving.pointsInPolygon(points, polygon) where points are moving.Point and polygon is a shapely polygon.
author Nicolas Saunier <nicolas.saunier@polymtl.ca>
date Tue, 16 Jul 2013 17:00:17 -0400
parents 924e38c9f70e
children d0b86ed50f32
files python/moving.py python/tests/moving.txt python/utils.py
diffstat 3 files changed, 89 insertions(+), 56 deletions(-) [+]
line wrap: on
line diff
--- 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)
--- a/python/tests/moving.txt	Tue Jul 16 01:36:50 2013 -0400
+++ b/python/tests/moving.txt	Tue Jul 16 17:00:17 2013 -0400
@@ -49,9 +49,9 @@
 >>> Point.distanceNorm2(Point(3,4),Point(1,7))
 3.605551275463989
 
->>> Point(3,2).inPolygon([Point(0,0),Point(1,0),Point(1,1),Point(0,1)])
+>>> Point(3,2).inPolygonNoShapely(np.array([[0,0],[1,0],[1,1],[0,1]]))
 False
->>> Point(3,2).inPolygon([Point(0,0),Point(4,0),Point(4,3),Point(0,3)])
+>>> Point(3,2).inPolygonNoShapely(np.array([[0,0],[4,0],[4,3],[0,3]]))
 True
 
 >>> predictPositionNoLimit(10, Point(0,0), Point(1,1)) # doctest:+ELLIPSIS
@@ -70,9 +70,9 @@
 True
 >>> t1[1]
 (1.500000,3.500000)
->>> t1.getTrajectoryInPolygon(np.array([[0,0],[4,0],[4,3],[0,3]]))
+>>> t1.getTrajectoryInPolygonNoShapely(np.array([[0,0],[4,0],[4,3],[0,3]]))
 (0.500000,0.500000)
->>> t1.getTrajectoryInPolygon(np.array([[10,10],[14,10],[14,13],[10,13]])).length()
+>>> t1.getTrajectoryInPolygonNoShapely(np.array([[10,10],[14,10],[14,13],[10,13]])).length()
 0
 
 >>> from utils import LCSS
--- a/python/utils.py	Tue Jul 16 01:36:50 2013 -0400
+++ b/python/utils.py	Tue Jul 16 17:00:17 2013 -0400
@@ -218,13 +218,14 @@
         self.lengthFunc = lengthFunc
         self.alignmentShift = 0
 
-    def similarities(self, l1, l2):
+    def similarities(self, l1, l2, shift=0):
         from numpy import zeros, int as npint
         n1 = len(l1)
         n2 = len(l2)
         self.similarityTable = zeros((n1+1,n2+1), dtype = npint)
         for i in xrange(1,n1+1):
-            for j in xrange(max(1,i-self.delta),min(n2+1,i+self.delta+1)):
+            for j in xrange(max(1,i-shift-self.delta),min(n2+1,i-shift+self.delta+1)):
+                #print max(1,i-shift-self.delta),min(n2+1,i-shift+self.delta+1)
                 if self.similarityFunc(l1[i-1], l2[j-1]):
                     self.similarityTable[i,j] = self.similarityTable[i-1,j-1]+1
                 else:
@@ -249,51 +250,61 @@
 
         eg distance(p1, p2) < epsilon
         '''
-        from numpy import argmax
+        #from numpy import argmax
         self.similarities(l1, l2)
-        self.similarityTable = self.similarityTable[:, :min(len(l2), len(l1)+self.delta)+1]
+        #self.similarityTable = self.similarityTable[:, :min(len(l2), len(l1)+self.delta)+1]
         if computeSubSequence:
             self.subSequenceIndices = self.subSequence(len(l1), len(l2))
-        return self.similarityTable[-1,-1]
+        return self.similarityTable.max()
 
     def _compute(self, _l1, _l2, computeSubSequence = False):
         '''returns the best matching if using a finite delta by shiftinig the series alignments'''
+        if len(_l2) < len(_l1): # l1 is the shortest
+            l1 = _l2
+            l2 = _l1
+            revertIndices = True
+        else:
+            l1 = _l1
+            l2 = _l2
+            revertIndices = False
+        n1 = len(l1)
+        n2 = len(l2)
+
         if self.aligned:
-            from numpy import argmax
-            if len(_l2) < len(_l1): # l1 is the shortest
-                l1 = _l2
-                l2 = _l1
-                revertIndices = True
-            else:
-                l1 = _l1
-                l2 = _l2
-                revertIndices = False
-            n1 = len(l1)
-            n2 = len(l2)
             # for i in xrange(min(delta,n1), max(n1+n2-delta, n2+1)): # i is the alignment of the end of l1 in l2
             #     print l1[min(-i-1,n1):] # min(n1+n2-i,n1)
             #     print l2[max(0,i-n1):]
             #     print LCSS(l1[min(-i-1,n1):], l2[max(0,i-n1):], similarityFunc, delta)
             lcssValues = {}
             similarityTables = {}
-            for i in xrange(min(self.delta,n1), max(n1+n2-self.delta, n2+1)):
+            #for i in xrange(min(self.delta,n1), max(n1+n2-self.delta, n2+1)):
+            for i in xrange(-max(n1,n2)-self.delta, +max(n1,n2)+self.delta):
                 #print l1[min(-i-1,n1):] # min(n1+n2-i,n1)
                 #print l2[max(0,i-n1):]
-                lcssValues[i] = self.computeLCSS(l1[min(-i-1,n1):], l2[max(0,i-n1):])
+                #lcssValues[i] = self.computeLCSS(l1[min(-i-1,n1):], l2[max(0,i-n1):])
                 #print i, lcssValues[i]
+                lcssValues[i] = self.similarities(l1, l2, i)
                 similarityTables[i] = self.similarityTable
-            imax = argMaxDict(lcssValues)
-            self.similarityTable = similarityTables[imax]
-            self.subSequenceIndices = self.subSequence(self.similarityTable.shape[0]-1, self.similarityTable.shape[1]-1) 
+                #print i
+                print self.similarityTable
+            self.alignmentShift = argMaxDict(lcssValues)
+            self.similarityTable = similarityTables[self.alignmentShift]
+            # do the subsequence computation here, once similarityTable is set
+            #self.subSequenceIndices = self.subSequence(self.similarityTable.shape[0]-1, self.similarityTable.shape[1]-1) 
+            #lcss = lcssValues[imax]
+        else:
+            self.alignmentShift = 0
+            self.similarities(l1, l2)
+        self.similarityTable = self.similarityTable[:, :min(len(l2), len(l1)+self.delta)+1]
+        if computeSubSequence:
+            self.subSequenceIndices = self.subSequence(self.similarityTable.shape[0]-1, self.similarityTable.shape[1]-1)
             if revertIndices:
-                self.subSequenceIndices = [(j+imax-n1,i) for i,j in self.subSequenceIndices]
-                self.alignmentShift = imax-n1
+                self.subSequenceIndices = [(j+self.alignmentShift,i) for i,j in self.subSequenceIndices]
+            #self.alignmentShift = imax-n1
             else:
-                self.subSequenceIndices = [(i+n1-imax,j) for i,j in self.subSequenceIndices]
-                self.alignmentShift = n1-imax
-            return lcssValues[imax]
-        else:
-            return self.computeLCSS(_l1, _l2, computeSubSequence)
+                self.subSequenceIndices = [(i+self.alignmentShift,j) for i,j in self.subSequenceIndices]
+            #self.alignmentShift = n1-imax
+        return self.similarityTable[-1,-1]
 
     def compute(self, l1, l2, computeSubSequence = False):
         '''get methods are to be shadowed in child classes '''