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

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


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

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

<span style="font-size:18px;">#
#===============================================================================
# Clipper class (+ data structs & ancilliary functions)
#===============================================================================
def _IntersectPoint(edge1, edge2):
    if _SlopesEqual2(edge1, edge2):
        if (edge2.ybot > edge1.ybot): y = edge2.ybot 
        else: y = edge1.ybot
        return Point(0, y), False
    if edge1.dx == 0:
        x = edge1.xBot
        if edge2.dx == horizontal:
            y = edge2.yBot
        else:
            b2 = edge2.yBot - Decimal(edge2.xBot)/edge2.dx
            y = round(Decimal(x)/edge2.dx + b2)
    elif edge2.dx == 0:
        x = edge2.xBot
        if edge1.dx == horizontal:
            y = edge1.yBot
        else:
            b1 = edge1.yBot - Decimal(edge1.xBot)/edge1.dx
            y = round(Decimal(x)/edge1.dx + b1)
    else:
        b1 = edge1.xBot - edge1.yBot * edge1.dx
        b2 = edge2.xBot - edge2.yBot * edge2.dx
        m    = Decimal(b2-b1)/(edge1.dx - edge2.dx)
        y    = round(m)
        if math.fabs(edge1.dx) < math.fabs(edge2.dx):
            x = round(edge1.dx * m + b1)
        else:
            x = round(edge2.dx * m + b2)
    if (y < edge1.yTop) or (y < edge2.yTop):
        if (edge1.yTop > edge2.yTop):
            return Point(edge1.xTop,edge1.yTop), _TopX(edge2, edge1.yTop) < edge1.xTop
        else:
            return Point(edge2.xTop,edge2.yTop), _TopX(edge1, edge2.yTop) > edge2.xTop
    else:
        return Point(x,y), True

def _TopX(e, currentY):
    if currentY == e.yTop: return e.xTop
    elif e.xTop == e.xBot: return e.xBot
    else: return e.xBot + round(e.dx * Decimal(currentY - e.yBot))

def _E2InsertsBeforeE1(e1,e2):
    if (e2.xCurr == e1.xCurr): 
        if (e2.yTop > e1.yTop):
            return e2.xTop < _TopX(e1, e2.yTop) 
        return e1.xTop > _TopX(e2, e1.yTop) 
    else: 
        return e2.xCurr < e1.xCurr

def _IsMinima(e):
    return e is not None and e.prevE.nextInLML != e and e.nextE.nextInLML != e

def _IsMaxima(e, y):
    return e is not None and e.yTop == y and e.nextInLML is None

def _IsIntermediate(e, y):
    return e.yTop == y and e.nextInLML is not None

def _GetMaximaPair(e):
    if not _IsMaxima(e.nextE, e.yTop) or e.nextE.xTop != e.xTop:
        return e.prevE
    else:
        return e.nextE

def _GetnextInAEL(e, direction):
    if direction == Direction.LeftToRight: return e.nextInAEL
    else: return e.prevInAEL

def _ProtectLeft(val):
    if val: return Protects.Both
    else: return Protects.Right

def _ProtectRight(val):
    if val: return Protects.Both
    else: return Protects.Left

def _GetDx(pt1, pt2):
    if (pt1.y == pt2.y): return horizontal
    else: return Decimal(pt2.x - pt1.x)/(pt2.y - pt1.y)

def _Param1RightOfParam2(outRec1, outRec2):
    while outRec1 is not None:
        outRec1 = outRec1.FirstLeft
        if outRec1 == outRec2: return True
    return False

def _FirstParamIsbottomPt(btmPt1, btmPt2):
    p = btmPt1.prevOp
    while _PointsEqual(p.pt, btmPt1.pt) and (p != btmPt1): p = p.prevOp
    dx1p = abs(_GetDx(btmPt1.pt, p.pt))
    p = btmPt1.nextOp
    while _PointsEqual(p.pt, btmPt1.pt) and (p != btmPt1): p = p.nextOp
    dx1n = abs(_GetDx(btmPt1.pt, p.pt))

    p = btmPt2.prevOp
    while _PointsEqual(p.pt, btmPt2.pt) and (p != btmPt2): p = p.prevOp
    dx2p = abs(_GetDx(btmPt2.pt, p.pt))
    p = btmPt2.nextOp
    while _PointsEqual(p.pt, btmPt2.pt) and (p != btmPt2): p = p.nextOp
    dx2n = abs(_GetDx(btmPt2.pt, p.pt))
    return (dx1p >= dx2p and dx1p >= dx2n) or (dx1n >= dx2p and dx1n >= dx2n)

def _GetBottomPt(pp):
    dups = None
    p = pp.nextOp
    while p != pp:
        if p.pt.y > pp.pt.y:
            pp = p
            dups = None
        elif p.pt.y == pp.pt.y and p.pt.x <= pp.pt.x:
            if p.pt.x < pp.pt.x:
                dups = None
                pp = p
            else:
                if p.nextOp != pp and p.prevOp != pp: dups = p
        p = p.nextOp
    if dups is not None:
        while dups != p:
            if not _FirstParamIsbottomPt(p, dups): pp = dups
            dups = dups.nextOp
            while not _PointsEqual(dups.pt, pp.pt): dups = dups.nextOp
    return pp

def _GetLowermostRec(outRec1, outRec2):
    if (outRec1.bottomPt is None): 
        outPt1 = _GetBottomPt(outRec1.pts)
    else: outPt1 = outRec1.bottomPt
    if (outRec2.bottomPt is None): 
        outPt2 = _GetBottomPt(outRec2.pts)
    else: outPt2 = outRec2.bottomPt
    if (outPt1.pt.y > outPt2.pt.y): return outRec1
    elif (outPt1.pt.y < outPt2.pt.y): return outRec2
    elif (outPt1.pt.x < outPt2.pt.x): return outRec1
    elif (outPt1.pt.x > outPt2.pt.x): return outRec2
    elif (outPt1.nextOp == outPt1): return outRec2
    elif (outPt2.nextOp == outPt2): return outRec1
    elif _FirstParamIsbottomPt(outPt1, outPt2): return outRec1
    else: return outRec2

def _SetHoleState(e, outRec, polyOutList):
    isHole = False
    e2 = e.prevInAEL
    while e2 is not None:
        if e2.outIdx >= 0:
            isHole = not isHole
            if outRec.FirstLeft is None:
                outRec.FirstLeft = polyOutList[e2.outIdx]
        e2 = e2.prevInAEL
    outRec.isHole = isHole

def _PointCount(pts):
    if pts is None: return 0
    p = pts
    result = 0
    while True:
        result += 1
        p = p.nextOp
        if p == pts: break
    return result

def _PointIsVertex(pt, outPts):
    op = outPts
    while True:
        if _PointsEqual(op.pt, pt): return True
        op = op.nextOp
        if op == outPts: break
    return False
               
def _ReversePolyPtLinks(pp):
    if pp is None: return
    pp1 = pp
    while True:
        pp2 = pp1.nextOp
        pp1.nextOp = pp1.prevOp
        pp1.prevOp = pp2;
        pp1 = pp2
        if pp1 == pp: break

def _FixupOutPolygon(outRec):
    lastOK = None
    outRec.bottomPt = None
    pp = outRec.pts
    while True:
        if pp.prevOp == pp or pp.nextOp == pp.prevOp:
            outRec.pts = None
            return
        if _PointsEqual(pp.pt, pp.nextOp.pt) or \
                _SlopesEqual(pp.prevOp.pt, pp.pt, pp.nextOp.pt):
            lastOK = None
            pp.prevOp.nextOp = pp.nextOp
            pp.nextOp.prevOp = pp.prevOp
            pp = pp.prevOp
        elif pp == lastOK: break
        else:
            if lastOK is None: lastOK = pp
            pp = pp.nextOp
    outRec.pts = pp

def _FixHoleLinkage(outRec):
    if outRec.FirstLeft is None or \
        (outRec.isHole != outRec.FirstLeft.isHole and \
            outRec.FirstLeft.pts is not None): return
    orfl = outRec.FirstLeft
    while orfl is not None and \
            (orfl.isHole == outRec.isHole or orfl.pts is None):
        orfl = orfl.FirstLeft
    outRec.FirstLeft = orfl
    
def _GetOverlapSegment(pt1a, pt1b, pt2a, pt2b):
    # precondition: segments are co-linear
    if abs(pt1a.x - pt1b.x) > abs(pt1a.y - pt1b.y):
        if pt1a.x > pt1b.x: tmp = pt1a; pt1a = pt1b; pt1b = tmp
        if pt2a.x > pt2b.x: tmp = pt2a; pt2a = pt2b; pt2b = tmp
        if (pt1a.x > pt2a.x): pt1 = pt1a
        else: pt1 = pt2a
        if (pt1b.x < pt2b.x): pt2 = pt1b
        else: pt2 = pt2b
        return pt1, pt2, pt1.x < pt2.x
    else:
        if pt1a.y < pt1b.y: tmp = pt1a; pt1a = pt1b; pt1b = tmp 
        if pt2a.y < pt2b.y: tmp = pt2a; pt2a = pt2b; pt2b = tmp
        if (pt1a.y < pt2a.y): pt1 = pt1a 
        else: pt1 = pt2a
        if (pt1b.y > pt2b.y): pt2 = pt1b 
        else: pt2 = pt2b
        return pt1, pt2, pt1.y > pt2.y

    
def _FindSegment(outPt, pt1, pt2):
    if outPt is None: return outPt, pt1, pt2, False
    pt1a = pt1; pt2a = pt2
    outPt2 = outPt
    while True:
        if _SlopesEqual(pt1a, pt2a, outPt.pt, outPt.prevOp.pt) and _SlopesEqual(pt1a, pt2a, outPt.pt):
            pt1, pt2, overlap = _GetOverlapSegment(pt1a, pt2a, outPt.pt, outPt.prevOp.pt)
            if overlap: return outPt, pt1, pt2, True
        outPt = outPt.nextOp
        if outPt == outPt2: return outPt, pt1, pt2, False

def _Pt3IsBetweenPt1AndPt2(pt1, pt2, pt3):
    if _PointsEqual(pt1, pt3) or _PointsEqual(pt2, pt3): return True
    elif pt1.x != pt2.x: return (pt1.x < pt3.x) == (pt3.x < pt2.x)
    else: return (pt1.y < pt3.y) == (pt3.y < pt2.y)

def _InsertPolyPtBetween(outPt1, outPt2, pt):
    if outPt1 == outPt2: raise Exception("JoinError")
    result = OutPt(outPt1.idx, pt)
    if outPt2 == outPt1.nextOp:
        outPt1.nextOp = result
        outPt2.prevOp = result
        result.nextOp = outPt2
        result.prevOp = outPt1
    else:
        outPt2.nextOp = result
        outPt1.prevOp = result
        result.nextOp = outPt1
        result.prevOp = outPt2
    return result

def _PointOnLineSegment(pt, linePt1, linePt2):
    return ((pt.x == linePt1.x) and (pt.y == linePt1.y)) or \
        ((pt.x == linePt2.x) and (pt.y == linePt2.y)) or \
        (((pt.x > linePt1.x) == (pt.x < linePt2.x)) and \
        ((pt.y > linePt1.y) == (pt.y < linePt2.y)) and \
        ((pt.x - linePt1.x) * (linePt2.y - linePt1.y) == \
        (linePt2.x - linePt1.x) * (pt.y - linePt1.y)))

def _PointOnPolygon(pt, pp):
    pp2 = pp;
    while True:
        if (_PointOnLineSegment(pt, pp2.pt, pp2.nextOp.pt)):
            return True
        pp2 = pp2.nextOp
        if (pp2 == pp): return False

def _PointInPolygon(pt, outPt): 
    result = False
    outPt2 = outPt
    while True:
        if ((((outPt2.pt.y <= pt.y) and (pt.y < outPt2.prevOp.pt.y)) or \
            ((outPt2.prevOp.pt.y <= pt.y) and (pt.y < outPt2.pt.y))) and \
            (pt.x < (outPt2.prevOp.pt.x - outPt2.pt.x) * (pt.y - outPt2.pt.y) / \
            (outPt2.prevOp.pt.y - outPt2.pt.y) + outPt2.pt.x)): result = not result
        outPt2 = outPt2.nextOp
        if outPt2 == outPt: break

def _Poly2ContainsPoly1(outPt1, outPt2):
    pt = outPt1
    if (_PointOnPolygon(pt.pt, outPt2)):
        pt = pt.nextOp
        while (pt != outPt1 and _PointOnPolygon(pt.pt, outPt2)):
            pt = pt.nextOp
        if (pt == outPt1): return True
    return _PointInPolygon(pt.pt, outPt2)    
    
def _EdgesAdjacent(inode):
    return (inode.e1.nextInSEL == inode.e2) or \
        (inode.e1.prevInSEL == inode.e2)

def _UpdateOutPtIdxs(outrec):
    op = outrec.pts
    while True:
        op.idx = outrec.idx
        op = op.prevOp
        if (op == outrec.pts): break

class Clipper(ClipperBase):
    #算上__init__，共有49个方法
        #__init__
        #_Reset
        #Clear
        #_InsertScanbeam
        #_PopScanbeam
        #_SetWindingCount
        #_IsEvenOddFillType
        #_IsEvenOddAltFillType
        #_IsContributing
        #_AddEdgeToSEL
        #_CopyAELToSEL
        #_InsertEdgeIntoAEL
        #_InsertLocalMinimaIntoAEL
        #_SwapPositionsInAEL
        #_SwapPositionsInSEL
        #_IsTopHorz
        #_ProcessHorizontal
        #_ProcessHorizontals
        #_AddJoin
        #_FixupJoinRecs
        #_AddHorzJoin
        #_InsertIntersectNode
        #_ProcessIntersections
        #_BuildIntersectList
        #_ProcessIntersectList
        #_DeleteFromAEL
        #_DeleteFromSEL
        #_IntersectEdges
        #_DoMaxima
        #_UpdateEdgeIntoAEL
        #_AddLocalMinPoly
        #_AddLocalMaxPoly
        #_CreateOutRec
        #_AddOutPt
        #_AppendPolygon
        #_FixupIntersectionOrder
        #_ProcessEdgesAtTopOfScanbeam
        #_Area
        #_JoinPoints
        #_FixupFirstLefts1
        #_FixupFirstLefts2
        #_GetOutRec
        #_JoinCommonEdges
        #_DoSimplePolygons
        #_ExecuteInternal
        #Execute
        #Execute2
        #_BuildResult
        #_BuildResult2
    

    def __init__(self):
        ClipperBase.__init__(self)

        self.ReverseOutput     = False
        self.ForceSimple       = False
        
        self._PolyOutList = []        
        self._ClipType         = ClipType.Intersection
        self._Scanbeam         = None
        self._ActiveEdges      = None
        self._SortedEdges      = None
        self._IntersectNodes   = None
        self._ClipFillType     = PolyFillType.EvenOdd
        self._SubjFillType     = PolyFillType.EvenOdd
        self._ExecuteLocked    = False
        self._UsingPolyTree    = False
        self._JoinList         = None
        self._HorzJoins        = None

    #以下4个方法
    def _Reset(self):
        ClipperBase._Reset(self)
        self._Scanbeam = None
        self._PolyOutList = []
        lm = self._LocalMinList
        while lm is not None:
            self._InsertScanbeam(lm.y)
            lm = lm.nextLm

    def Clear(self):
        self._PolyOutList = []
        ClipperBase.Clear(self)

    def _InsertScanbeam(self, y):
        if self._Scanbeam is None:
            self._Scanbeam = Scanbeam(y)
        elif y > self._Scanbeam.y:
            self._Scanbeam = Scanbeam(y, self._Scanbeam)
        else:
            sb = self._Scanbeam
            while sb.nextSb is not None and y <= sb.nextSb.y:
                sb = sb.nextSb
            if y == sb.y: return
            newSb = Scanbeam(y, sb.nextSb)
            sb.nextSb = newSb

    def _PopScanbeam(self):
        result = self._Scanbeam.y
        self._Scanbeam = self._Scanbeam.nextSb
        return result

    #以下3个方法, 有关属性
    def _SetWindingCount(self, edge):
        e = edge.prevInAEL
        while e is not None and e.PolyType != edge.PolyType:
            e = e.prevInAEL
        if e is None:
            edge.windCnt = edge.windDelta
            edge.windCnt2 = 0
            e = self._ActiveEdges
        elif self._IsEvenOddFillType(edge):
            edge.windCnt = 1
            edge.windCnt2 = e.windCnt2
            e = e.nextInAEL
        else:
            if e.windCnt * e.windDelta < 0:
                if (abs(e.windCnt) > 1):
                    if (e.windDelta * edge.windDelta < 0): edge.windCnt = e.windCnt
                    else: edge.windCnt = e.windCnt + edge.windDelta
                else:
                    edge.windCnt = e.windCnt + e.windDelta + edge.windDelta
            elif (abs(e.windCnt) > 1) and (e.windDelta * edge.windDelta < 0):
                edge.windCnt = e.windCnt
            elif e.windCnt + edge.windDelta == 0:
                edge.windCnt = e.windCnt
            else:
                edge.windCnt = e.windCnt + edge.windDelta
            edge.windCnt2 = e.windCnt2
            e = e.nextInAEL
        # update windCnt2 ...
        if self._IsEvenOddAltFillType(edge):
            while (e != edge):
                if edge.windCnt2 == 0: edge.windCnt2 = 1
                else: edge.windCnt2 = 0
                e = e.nextInAEL
        else:
            while (e != edge):
                edge.windCnt2 += e.windDelta
                e = e.nextInAEL

    def _IsEvenOddFillType(self, edge):
        if edge.PolyType == PolyType.Subject:
            return self._SubjFillType == PolyFillType.EvenOdd
        else:
            return self._ClipFillType == PolyFillType.EvenOdd

    def _IsEvenOddAltFillType(self, edge):
        if edge.PolyType == PolyType.Subject:
            return self._ClipFillType == PolyFillType.EvenOdd
        else:
            return self._SubjFillType == PolyFillType.EvenOdd

    #交、并、差、异或的每种运算所保留下来的点，在这里得到区别对待
    def _IsContributing(self, edge):
        if edge.PolyType == PolyType.Subject:
            pft = self._SubjFillType
            pft2 = self._ClipFillType
        else:
            pft = self._ClipFillType
            pft2 = self._SubjFillType
        if pft &#