[从头学数学] 第286节 [计算几何] 多边形的布尔运算（上）

剧情提要：
阿伟看到了一本比较有趣的书，是关于《计算几何》的，2008年由北清派出版。很好奇
它里面讲了些什么，就来看看啦。


正剧开始：
星历2016年10月21日 14:37:23, 银河系厄尔斯星球中华帝国江南行省。

[工程师阿伟]正在和[机器小伟]一起研究[计算几何]]。

一、生成随机多边形

<span style="font-size:18px;">#import mathfrom collections import namedtuplefrom decimal import Decimal, getcontextimport refrom random import randint# from subprocess import callimport osgetcontext().prec = 8horizontal = Decimal('-Infinity')Point = namedtuple('Point', 'x y')DoublePoint = namedtuple('DoublePoint', 'x y')#生成随机多边形def RandomPoly(maxWidth, maxHeight, vertCnt):    result = []    for _ in range(vertCnt):        result.append(Point(randint(0, maxWidth), randint(0, maxHeight)))    return resultdef tmp():    scale = 1;    a = RandomPoly(640 * scale, 480 * scale, 10);    print('a = ', a);    b = RandomPoly(640 * scale, 480 * scale, 16);    print('b = ', b);if __name__ == '__main__':    tmp();>>> a =  [Point(x=376, y=223), Point(x=471, y=386), Point(x=110, y=381), Point(x=206, y=85), Point(x=190, y=41), Point(x=547, y=249), Point(x=501, y=75), Point(x=501, y=183), Point(x=497, y=240), Point(x=528, y=83)]b =  [Point(x=432, y=428), Point(x=558, y=145), Point(x=64, y=52), Point(x=381, y=467), Point(x=460, y=321), Point(x=614, y=318), Point(x=173, y=80), Point(x=301, y=399), Point(x=431, y=129), Point(x=82, y=79), Point(x=449, y=71), Point(x=300, y=29), Point(x=558, y=225), Point(x=5, y=87), Point(x=613, y=327), Point(x=342, y=353)]>>> #</span>

二、有关数据的文件操作

<span style="font-size:18px;">##读取/保存多边形数据到文件def LoadFile1(lines):    # File type 1: first line is total polygons count and subsequent lines     # contain the polygon vertex count followed by its coords     try:        polygons = []        poly = []        for l in lines:            vals = re.split(' |, |,', l.strip())            if len(vals) < 2:                  if (len(poly)  > 2):                    polygons.append(poly)                poly = []            else:                 poly.append(Point(int(vals[0]), int(vals[1])))        if (len(poly)  > 2):            polygons.append(poly)        return polygons    except:        return Nonedef LoadFile2(lines):    # File type 2: vertex coords on consecutive lines for each polygon     # where each polygon is separated by an empty line     try:        polygons = []        poly = []        for l in lines:            l = l.strip()            if (l == ''):                 if (len(poly)  > 2):                    polygons.append(poly)                poly = []            else:                 vals = re.split(' |, |,', l)                poly.append(Point(int(vals[0]), int(vals[1])))        if (len(poly)  > 2):            polygons.append(poly)        return polygons    except:        return Nonedef LoadFile(filename):    try:        f = open(filename, 'r')        try:            lines = f.readlines()        finally:            f.close()        # pick file type from format of first line ...        if len(lines) == 0: return []        elif not ',' in lines[0]: return LoadFile1(lines)        else: return LoadFile2(lines)    except:        return None    def SaveToFile(filename, polys, scale = 1.0):    invScale = 1.0 / scale    try:        f = open(filename, 'w')        try:            if invScale == 1:                for poly in polys:                    for pt in poly:                        f.write("{0}, {1}\n".format(pt.x, pt.y))                    f.write("\n")            else:                for poly in polys:                    for pt in poly:                        f.write("{0:.4f}, {1:.4f}\n".format(pt.x * invScale, pt.y * invScale))                    f.write("\n")        finally:            f.close()    except:        return#测试多边形数据生成def tmp_0():    scale = 1;    a = RandomPoly(640 * scale, 480 * scale, 10);    print('a = ', a);    b = RandomPoly(640 * scale, 480 * scale, 16);    print('b = ', b);#测试多边形数据写入文件def tmp_1():    scaleExp = 0    scale = math.pow(10, scaleExp)    invScale = 1.0 / scale       subj, clip = [], []    subj.append(RandomPoly(640 * scale, 480 * scale, 10))    clip.append(RandomPoly(640 * scale, 480 * scale, 16))    SaveToFile('./subj2.txt', subj, scale)    SaveToFile('./clip2.txt', clip, scale)    print('写入文件完毕。');>>> subj =  [[Point(x=199, y=133), Point(x=577, y=464), Point(x=468, y=386), Point(x=130, y=128), Point(x=336, y=140), Point(x=324, y=418), Point(x=91, y=223), Point(x=67, y=455), Point(x=168, y=458), Point(x=639, y=227)]]clip =  [[Point(x=233, y=423), Point(x=193, y=76), Point(x=369, y=373), Point(x=538, y=416), Point(x=155, y=301), Point(x=215, y=323), Point(x=82, y=473), Point(x=136, y=405), Point(x=466, y=145), Point(x=250, y=336), Point(x=480, y=218), Point(x=135, y=320), Point(x=492, y=135), Point(x=40, y=361), Point(x=467, y=323), Point(x=347, y=26)]]>>> #获取多边形坐标点数据def tmp():    subj = LoadFile('./subj2.txt')    clip = LoadFile('./clip2.txt')    vert_subj = [];    vert_clip = [];    for i in range(len(subj[0])):        vert_subj.append([subj[0][i].x, subj[0][i].y]);    for i in range(len(clip[0])):        vert_clip.append([clip[0][i].x, clip[0][i].y]);    print('vert_subj = ', vert_subj);    print('vert_clip = ', vert_clip);>>> vert_subj =  [[199, 133], [577, 464], [468, 386], [130, 128], [336, 140], [324, 418], [91, 223], [67, 455], [168, 458], [639, 227]]vert_clip =  [[233, 423], [193, 76], [369, 373], [538, 416], [155, 301], [215, 323], [82, 473], [136, 405], [466, 145], [250, 336], [480, 218], [135, 320], [492, 135], [40, 361], [467, 323], [347, 26]]>>> #</span>

<span style="font-size:18px;">//if (1) {var r = 20;              config.setSector(8,15,7,2);                config.graphPaper2D(0, 0, r);              config.axis2D(0, 0, 450, 1.2);   //坐标轴设定                var scaleX = 2*r, scaleY = 2*r;                  var spaceX = 50, spaceY = 50;                   var xS = 0, xE = 640;                  var yS = 0, yE = 480;                  config.axisSpacing(xS, xE, spaceX, scaleX, 'X');                    config.axisSpacing(yS, yE, spaceY, scaleY, 'Y');                                        var transform = new Transform();    var lable = [];for (var i = 0; i < 100; i++) {lable.push(i.toFixed(0));}var colorArray = ['red', 'orange', 'yellow', 'green', 'cyan', 'blue', 'purple'];var color = 0;var seg = [];//主多边形seg = transform.scale($vert_subj, scaleX/spaceX, scaleY/spaceY);shape.pointDraw([].concat(seg), 'red', 1, 1, lable);//从多边形seg = transform.scale($vert_clip, scaleX/spaceX, scaleY/spaceY);shape.pointDraw([].concat(seg), 'blue', 1, 1, lable);plot.fillText('图：主多边形和从多边形的顶点', 10, 60, 300);}//</span>

三、外围方法与属性

<span style="font-size:18px;">#class ClipType: (Intersection, Union, Difference, Xor) = range(4)class PolyType:    (Subject, Clip) = range(2)class PolyFillType: (EvenOdd, NonZero, Positive, Negative) = range(4)class JoinType: (Square, Round, Miter) = range(3)class EndType: (Closed, Butt, Square, Round) = range(4)class EdgeSide: (Left, Right) = range(2)class Protects: (Neither, Left, Right, Both) = range(4)class Direction: (LeftToRight, RightToLeft) = range(2)class LocalMinima(object):    leftBound = rightBound = nextLm = None    #根据坐标的y值确定它的等高扫描线的左右范围    def __init__(self, y, leftBound, rightBound):        self.y = y        self.leftBound = leftBound        self.rightBound = rightBoundclass Scanbeam(object):    __slots__ = ('y','nextSb')    def __init__(self, y, nextSb = None):        self.y = y        self.nextSb = nextSb    def __repr__(self):        s = 'None'        if self.nextSb is not None: s = '<obj>'        return "(y:%i, nextSb:%s)" % (self.y, s)#交点class IntersectNode(object):    __slots__ = ('e1','e2','pt','nextIn')    def __init__(self, e1, e2, pt):        self.e1 = e1        self.e2 = e2        self.pt = pt        self.nextIn = Noneclass OutPt(object):    __slots__ = ('idx','pt','prevOp','nextOp')    def __init__(self, idx, pt):        #从最外圈数进去第几个环，从-1开始计数        self.idx = idx        self.pt = pt        self.prevOp = None        self.nextOp = Noneclass OutRec(object):    __slots__ = ('idx','bottomPt','isHole','FirstLeft', 'pts','PolyNode')    def __init__(self, idx):        self.idx = idx        self.bottomPt = None        self.isHole = False        self.FirstLeft = None        self.pts = None        self.PolyNode = Noneclass JoinRec(object):    __slots__ = ('pt1a','pt1b','poly1Idx','pt2a', 'pt2b','poly2Idx')class HorzJoin(object):    edge = None    savedIdx = 0    prevHj = None    nextHj = None    def __init__(self, edge, idx):        self.edge = edge        self.savedIdx = idx        #===============================================================================# Unit global functions ...#===============================================================================def IntsToPoints(ints):    result = []    for i in range(0, len(ints), 2):        result.append(Point(ints[i], ints[i+1]))    return resultdef Area(polygon):    # see http://www.mathopenref.com/coordpolygonarea2.html    highI = len(polygon) - 1    A = (polygon[highI].x + polygon[0].x) * (polygon[0].y - polygon[highI].y)    for i in range(highI):        A += (polygon[i].x + polygon[i+1].x) * (polygon[i+1].y - polygon[i].y)    return float(A) / 2def Orientation(polygon):    return Area(polygon) > 0.0#===============================================================================# Edge class #===============================================================================class Edge(object):    def __init__(self):#界定范围的六个量        #边是从y轴从小到大划线, top代表边的y值较小的点        self.xBot, self.yBot, self.xCurr, self.yCurr, = 0, 0, 0, 0        self.xTop, self.yTop = 0, 0#斜率相关三个量        self.dx, self.deltaX , self.deltaY = Decimal(0), Decimal(0), Decimal(0)#类型两个量        self.polyType = PolyType.Subject         self.side = EdgeSide.Left#穿插回转指示量三个        self.windDelta, self.windCnt, self.windCnt2 = 0, 0, 0 #        self.outIdx = -1#相关边和线段七个量        self.nextE, self.prevE, self.nextInLML = None, None, None        self.prevInAEL, self.nextInAEL, self.prevInSEL, self.nextInSEL = None, None, None, None            def __repr__(self):        return "(%i,%i . %i,%i {dx:%0.2f} %i {%x})" % \            (self.xBot, self.yBot, self.xTop, self.yTop, self.dx, self.outIdx, id(self))    def __str__(self):        return "[[%i, %i], [%i, %i]]" % (self.xBot, self.yBot, self.xTop, self.yTop)        #===============================================================================# ClipperBase class (+ data structs & ancilliary functions)#===============================================================================#是否相同点def _PointsEqual(pt1, pt2):    return (pt1.x == pt2.x) and (pt1.y == pt2.y)#是否相同斜率def _SlopesEqual(pt1, pt2, pt3, pt4 = None):    if pt4 is None:        return (pt1.y-pt2.y)*(pt2.x-pt3.x) == (pt1.x-pt2.x)*(pt2.y-pt3.y)    else:        return (pt1.y-pt2.y)*(pt3.x-pt4.x) == (pt1.x-pt2.x)*(pt3.y-pt4.y)def _SlopesEqual2(e1, e2):    return e1.deltaY * e2.deltaX == e1.deltaX * e2.deltaY#边的斜率def _SetDx(e):    #两端点的dx    e.deltaX = Decimal(e.xTop - e.xBot)    #两端点的dy    e.deltaY = Decimal(e.yTop - e.yBot)    if e.deltaY == 0: e.dx = horizontal    else: e.dx = e.deltaX/e.deltaYdef _SwapSides(e1, e2):    side    = e1.side    e1.side = e2.side    e2.side = sidedef _SwapPolyIndexes(e1, e2):    idx       = e1.outIdx    e1.outIdx = e2.outIdx    e2.outIdx = idxdef _InitEdge(e, eNext, ePrev, pt, polyType):    e.nextE = eNext    e.prevE = ePrev    e.xCurr = pt.x    e.yCurr = pt.y    #边是从y轴从小到大划线, top代表边的y值较小的点    if e.yCurr >= e.nextE.yCurr:        e.xBot = e.xCurr        e.yBot = e.yCurr        e.xTop = e.nextE.xCurr        e.yTop = e.nextE.yCurr        e.windDelta = 1    else:        e.xTop = e.xCurr        e.yTop = e.yCurr        e.xBot = e.nextE.xCurr        e.yBot = e.nextE.yCurr        e.windDelta = -1    _SetDx(e)    e.outIdx = -1    e.PolyType = polyTypedef _SwapX(e):    e.xCurr = e.xTop    e.xTop = e.xBot    e.xBot = e.xCurr#===============================================================================# PolyNode & PolyTree classes (+ ancilliary functions)#===============================================================================class PolyNode(object):    """Node of PolyTree"""        def __init__(self):        self.Contour = []        self.Childs = []        self.Parent = None        self.Index = 0        self.ChildCount = 0        def IsHole(self):        result = True        while (self.Parent is not None):            result = not result            self.Parent = self.Parent.Parent        return result        def GetNext(self):        if (self.ChildCount > 0):            return self.Childs[0]        else:            return self._GetNextSiblingUp()        def _AddChild(self, node):        self.Childs.append(node)        node.Index = self.ChildCount        node.Parent = self        self.ChildCount += 1        def _GetNextSiblingUp(self):        if (self.Parent is None):            return None        elif (self.Index == self.Parent.ChildCount - 1):            return self.Parent._GetNextSiblingUp()        else:            return self.Parent.Childs[self.Index +1]class PolyTree(PolyNode):    """Container for PolyNodes"""    def __init__(self):        PolyNode.__init__(self)        self._AllNodes = []            def Clear(self):        self._AllNodes = []        self.Childs = []        self.ChildCount = 0        def GetFirst(self):        if (self.ChildCount > 0):            return self.Childs[0]        else:            return None        def Total(self):        return len(self._AllNodes)def _AddPolyNodeToPolygons(polynode, polygons):    """Internal function for PolyTreeToPolygons()"""    if (len(polynode.Contour) > 0):        polygons.append(polynode.Contour)    for i in range(polynode.ChildCount):        _AddPolyNodeToPolygons(polynode.Childs[i], polygons)def PolyTreeToPolygons(polyTree):    result = []    _AddPolyNodeToPolygons(polyTree, result)    return resultclass ClipperBase(object):    def __init__(self):        self._EdgeList      = []       # 2D array        self._LocalMinList  = None     # single-linked list of LocalMinima        self._CurrentLocMin = None            def _InsertLocalMinima(self, lm):        if self._LocalMinList is None:            self._LocalMinList = lm        elif lm.y >= self._LocalMinList.y:            lm.nextLm = self._LocalMinList            self._LocalMinList = lm        else:            tmp = self._LocalMinList            while tmp.nextLm is not None and lm.y < tmp.nextLm.y:                    tmp = tmp.nextLm            lm.nextLm = tmp.nextLm            tmp.nextLm = lm    def _AddBoundsToLML(self, e):        e.nextInLML = None        e = e.nextE        while True:            if e.dx == horizontal:                if (e.nextE.yTop < e.yTop) and (e.nextE.xBot > e.prevE.xBot): break                if (e.xTop != e.prevE.xBot): _SwapX(e)                e.nextInLML = e.prevE            elif e.yBot == e.prevE.yBot: break            else: e.nextInLML = e.prevE            e = e.nextE        if e.dx == horizontal:            if (e.xBot != e.prevE.xBot): _SwapX(e)            lm = LocalMinima(e.prevE.yBot, e.prevE, e)        elif (e.dx < e.prevE.dx):            lm = LocalMinima(e.prevE.yBot, e.prevE, e)        else:            lm = LocalMinima(e.prevE.yBot, e, e.prevE)        lm.leftBound.side = EdgeSide.Left        lm.rightBound.side = EdgeSide.Right        self._InsertLocalMinima(lm)        while True:            if e.nextE.yTop == e.yTop and e.nextE.dx != horizontal: break            e.nextInLML = e.nextE            e = e.nextE            if e.dx == horizontal and e.xBot != e.prevE.xTop: _SwapX(e)        return e.nextE    def _Reset(self):        lm = self._LocalMinList        if lm is not None: self._CurrentLocMin = lm        while lm is not None:            e = lm.leftBound            while e is not None:                e.xCurr    = e.xBot                e.yCurr    = e.yBot                e.side     = EdgeSide.Left                e.outIdx = -1                e = e.nextInLML            e = lm.rightBound            while e is not None:                e.xCurr    = e.xBot                e.yCurr    = e.yBot                e.side     = EdgeSide.Right                e.outIdx = -1                e = e.nextInLML            lm = lm.nextLm                def AddPolygon(self, polygon, polyType):        ln = len(polygon)        if ln < 3: return False        pg = polygon[:]        j = 0        # remove duplicate points and co-linear points        for i in range(1, len(polygon)):            if _PointsEqual(pg[j], polygon[i]):                 continue            elif (j > 0) and _SlopesEqual(pg[j-1], pg[j], polygon[i]):                if _PointsEqual(pg[j-1], polygon[i]): j -= 1            else: j += 1            pg[j] = polygon[i]        if (j < 2): return False        # remove duplicate points and co-linear edges at the loop around        # of the start and end coordinates ...        ln = j +1        while (ln > 2):            if _PointsEqual(pg[j], pg[0]): j -= 1            elif _PointsEqual(pg[0], pg[1]) or _SlopesEqual(pg[j], pg[0], pg[1]):                pg[0] = pg[j]                j -= 1            elif _SlopesEqual(pg[j-1], pg[j], pg[0]): j -= 1            elif _SlopesEqual(pg[0], pg[1], pg[2]):                for i in range(2, j +1): pg[i-1] = pg[i]                j -= 1            else: break            ln -= 1        if ln < 3: return False        edges = []        for i in range(ln):            edges.append(Edge())        edges[0].xCurr = pg[0].x        edges[0].yCurr = pg[0].y        _InitEdge(edges[ln-1], edges[0], edges[ln-2], pg[ln-1], polyType)        for i in range(ln-2, 0, -1):            _InitEdge(edges[i], edges[i+1], edges[i-1], pg[i], polyType)        _InitEdge(edges[0], edges[1], edges[ln-1], pg[0], polyType)        e = edges[0]        eHighest = e        while True:            e.xCurr = e.xBot            e.yCurr = e.yBot            if e.yTop < eHighest.yTop: eHighest = e            e = e.nextE            if e == edges[0]: break        # make sure eHighest is positioned so the following loop works safely ...        if eHighest.windDelta > 0: eHighest = eHighest.nextE        if eHighest.dx == horizontal: eHighest = eHighest.nextE        # finally insert each local minima ...        e = eHighest        while True:            e = self._AddBoundsToLML(e)            if e == eHighest: break        self._EdgeList.append(edges)    def AddPolygons(self, polygons, polyType):        result = False        for p in polygons:            if self.AddPolygon(p, polyType): result = True        return result    def Clear(self):        self._EdgeList = []        self._LocalMinList    = None        self._CurrentLocMin = None    def _PopLocalMinima(self):        if self._CurrentLocMin is not None:            self._CurrentLocMin = self._CurrentLocMin.nextLmclass Clipper(ClipperBase):    def __init__(self):        ClipperBase.__init__(self)#测试生成边集def tmp():    vert_subj =  [[199, 133], [577, 464], [468, 386], [130, 128], [336, 140], [324, 418], [91, 223], [67, 455], [168, 458], [639, 227]]    vert_clip =  [[233, 423], [193, 76], [369, 373], [538, 416], [155, 301], [215, 323], [82, 473], [136, 405], [466, 145], [250, 336], [480, 218], [135, 320], [492, 135], [40, 361], [467, 323], [347, 26]]    c = Clipper()    pft = PolyFillType.EvenOdd    subj, clip = [], [];    for i in range(len(vert_subj)):        subj.append(Point(vert_subj[i][0], vert_subj[i][1]));    for i in range(len(vert_clip)):        clip.append(Point(vert_clip[i][0], vert_clip[i][1]));        c.AddPolygons([subj], PolyType.Subject)    c.AddPolygons([clip], PolyType.Clip);    for i in range(len(c._EdgeList)):        s = '[';        print(len(c._EdgeList[i]));        for j in range(len(c._EdgeList[i])):            s += str(c._EdgeList[i][j]);            s += ',';        s += ']';        print(s);>>> 10[[[577, 464], [199, 133]],[[577, 464], [468, 386]],[[468, 386], [130, 128]],[[336, 140], [130, 128]],[[324, 418], [336, 140]],[[324, 418], [91, 223]],[[67, 455], [91, 223]],[[168, 458], [67, 455]],[[168, 458], [639, 227]],[[639, 227], [199, 133]],]16[[[233, 423], [193, 76]],[[369, 373], [193, 76]],[[538, 416], [369, 373]],[[538, 416], [155, 301]],[[215, 323], [155, 301]],[[82, 473], [215, 323]],[[82, 473], [136, 405]],[[136, 405], [466, 145]],[[250, 336], [466, 145]],[[250, 336], [480, 218]],[[135, 320], [480, 218]],[[135, 320], [492, 135]],[[40, 361], [492, 135]],[[40, 361], [467, 323]],[[467, 323], [347, 26]],[[233, 423], [347, 26]],]>>>#</span>

<span style="font-size:18px;">//if (1) {var r = 20;              config.setSector(8,15,7,2);                config.graphPaper2D(0, 0, r);              config.axis2D(0, 0, 450, 1.2);   //坐标轴设定                var scaleX = 2*r, scaleY = 2*r;                  var spaceX = 50, spaceY = 50;                   var xS = 0, xE = 640;                  var yS = 0, yE = 480;                  config.axisSpacing(xS, xE, spaceX, scaleX, 'X');                    config.axisSpacing(yS, yE, spaceY, scaleY, 'Y');                                        var transform = new Transform();    var lable = [];for (var i = 0; i < 100; i++) {lable.push(i.toFixed(0));}var colorArray = ['red', 'orange', 'yellow', 'green', 'cyan', 'blue', 'purple'];var color = 0;var seg = [];var subjEdges = $subj_edge.length;var clipEdges = $clip_edge.length;for (var i = 0; i < subjEdges; i++) {seg = transform.scale($subj_edge[i], scaleX/spaceX, scaleY/spaceY);shape.multiLineDraw([].concat(seg), 'red');}for (var i = 0; i < clipEdges; i++) {seg = transform.scale($clip_edge[i], scaleX/spaceX, scaleY/spaceY);shape.multiLineDraw([].concat(seg), 'blue');}//主多边形seg = transform.scale($vert_subj, scaleX/spaceX, scaleY/spaceY);shape.pointDraw([].concat(seg), 'red', 1, 1, lable);//从多边形seg = transform.scale($vert_clip, scaleX/spaceX, scaleY/spaceY);shape.pointDraw([].concat(seg), 'blue', 1, 1, lable);plot.fillText('图：主多边形和从多边形的边', 10, 60, 300);}//</span>

本节到此结束，欲知后事如何，请看下回分解。