Mercurial Hosting > traffic-intelligence
changeset 372:349eb1e09f45
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 '''