“判断一个点是否在一个多边形里”,一开始以为是个挺难的问题,但Google了一下之后发现其实蛮简单,所用到的算法叫做“Ray-casting Algorithm”,中文应该叫“光线投射算法”,这是维基百科的描述:[维基百科]
简单地说可以这么判断:从这个点引出一根“射线”,与多边形的任意若干条边相交,累计相交的边的数目,如果是奇数,那么点就在多边形内,否则点就在多边形外。
from collections import namedtuple from pprint import pprint as pp import sys Pt = namedtuple('Pt', 'x, y') # Point Edge = namedtuple('Edge', 'a, b') # Polygon edge from a to b Poly = namedtuple('Poly', 'name, edges') # Polygon _eps = 0.00001 _huge = sys.float_info.max _tiny = sys.float_info.min def rayintersectseg(p, edge): ''' takes a point p=Pt() and an edge of two endpoints a,b=Pt() of a line segment returns boolean ''' a,b = edge if a.y > b.y: a,b = b,a if p.y == a.y or p.y == b.y: p = Pt(p.x, p.y + _eps) intersect = False if (p.y > b.y or p.y < a.y) or ( p.x > max(a.x, b.x)): return False if p.x < min(a.x, b.x): intersect = True else: if abs(a.x - b.x) > _tiny: m_red = (b.y - a.y) / float(b.x - a.x) else: m_red = _huge if abs(a.x - p.x) > _tiny: m_blue = (p.y - a.y) / float(p.x - a.x) else: m_blue = _huge intersect = m_blue >= m_red return intersect def _odd(x): return x%2 == 1 def ispointinside(p, poly): ln = len(poly) return _odd(sum(rayintersectseg(p, edge) for edge in poly.edges )) def polypp(poly): print "\n Polygon(name='%s', edges=(" % poly.name print ' ', ',\n '.join(str(e) for e in poly.edges) + '\n ))' if __name__ == '__main__': polys = [ Poly(name='square', edges=( Edge(a=Pt(x=0, y=0), b=Pt(x=10, y=0)), Edge(a=Pt(x=10, y=0), b=Pt(x=10, y=10)), Edge(a=Pt(x=10, y=10), b=Pt(x=0, y=10)), Edge(a=Pt(x=0, y=10), b=Pt(x=0, y=0)) )), Poly(name='square_hole', edges=( Edge(a=Pt(x=0, y=0), b=Pt(x=10, y=0)), Edge(a=Pt(x=10, y=0), b=Pt(x=10, y=10)), Edge(a=Pt(x=10, y=10), b=Pt(x=0, y=10)), Edge(a=Pt(x=0, y=10), b=Pt(x=0, y=0)), Edge(a=Pt(x=2.5, y=2.5), b=Pt(x=7.5, y=2.5)), Edge(a=Pt(x=7.5, y=2.5), b=Pt(x=7.5, y=7.5)), Edge(a=Pt(x=7.5, y=7.5), b=Pt(x=2.5, y=7.5)), Edge(a=Pt(x=2.5, y=7.5), b=Pt(x=2.5, y=2.5)) )), Poly(name='strange', edges=( Edge(a=Pt(x=0, y=0), b=Pt(x=2.5, y=2.5)), Edge(a=Pt(x=2.5, y=2.5), b=Pt(x=0, y=10)), Edge(a=Pt(x=0, y=10), b=Pt(x=2.5, y=7.5)), Edge(a=Pt(x=2.5, y=7.5), b=Pt(x=7.5, y=7.5)), Edge(a=Pt(x=7.5, y=7.5), b=Pt(x=10, y=10)), Edge(a=Pt(x=10, y=10), b=Pt(x=10, y=0)), Edge(a=Pt(x=10, y=0), b=Pt(x=2.5, y=2.5)) )), Poly(name='exagon', edges=( Edge(a=Pt(x=3, y=0), b=Pt(x=7, y=0)), Edge(a=Pt(x=7, y=0), b=Pt(x=10, y=5)), Edge(a=Pt(x=10, y=5), b=Pt(x=7, y=10)), Edge(a=Pt(x=7, y=10), b=Pt(x=3, y=10)), Edge(a=Pt(x=3, y=10), b=Pt(x=0, y=5)), Edge(a=Pt(x=0, y=5), b=Pt(x=3, y=0)) )), ] testpoints = (Pt(x=5, y=5), Pt(x=5, y=8), Pt(x=-10, y=5), Pt(x=0, y=5), Pt(x=10, y=5), Pt(x=8, y=5), Pt(x=10, y=10)) print "\n TESTING WHETHER POINTS ARE WITHIN POLYGONS" for poly in polys: polypp(poly) print ' ', '\t'.join("%s: %s" % (p, ispointinside(p, poly)) for p in testpoints[:3]) print ' ', '\t'.join("%s: %s" % (p, ispointinside(p, poly)) for p in testpoints[3:6]) print ' ', '\t'.join("%s: %s" % (p, ispointinside(p, poly)) for p in testpoints[6:])
public static class RayCastingAlgorithm { public static bool IsWithin(Point pt, IList<Point> polygon, bool noneZeroMode) { int ptNum = polygon.Count(); if (ptNum < 3) { return false; } int j = ptNum - 1; bool oddNodes = false; int zeroState = 0; for (int k = 0; k < ptNum; k++) { Point ptK = polygon[k]; Point ptJ = polygon[j]; if (((ptK.Y > pt.Y) != (ptJ.Y > pt.Y)) && (pt.X < (ptJ.X - ptK.X) * (pt.Y - ptK.Y) / (ptJ.Y - ptK.Y) + ptK.X)) { oddNodes = !oddNodes; if (ptK.Y > ptJ.Y) { zeroState++; } else { zeroState--; } } j = k; } return noneZeroMode?zeroState!=0:oddNodes; } }