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

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


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

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

五、SVG绘图相关

<span style="font-size:18px;">#
#===============================================================================
# OffsetPolygons (+ ancilliary functions)
#===============================================================================

FloatPoint = namedtuple('FloatPoint', 'x y')
Rect = namedtuple('FloatPoint', 'left top right bottom')

def _GetUnitNormal(pt1, pt2):
    if pt2.x == pt1.x and pt2.y == pt1.y:
        return FloatPoint(0.0, 0.0)
    dx = pt2.x - pt1.x
    dy = pt2.y - pt1.y
    f = 1.0 / math.hypot(dx, dy)
    dx = float(dx) * f
    dy = float(dy) * f
    return FloatPoint(dy, -dx)

def _BuildArc(pt, a1, a2, r, limit):
    arcFrac = abs(a2 - a1) / (2 * math.pi);
    steps = int(arcFrac * math.pi / math.acos(1 - limit / abs(r)))
    if steps < 2: steps = 2
    elif steps > 222.0 * arcFrac:
        steps = int(222.0 * arcFrac)
    
    result = []
    y = math.sin(a1)
    x = math.cos(a1)
    s = math.sin((a2-a1)/steps)
    c = math.cos((a2-a1)/steps)
    for _ in range(steps+1):
        result.append(FloatPoint(pt.x + round(x * r), pt.y + round(y * r)))
        x2 = x
        x = x * c - s * y    # cross product & dot product here ...
        y = x2 * s + y * c   # avoids repeat calls to the much slower sin() & cos()
        
    return result

def _GetBounds(pts):
    left = None
    for poly in pts:
        for pt in poly:
            left = pt.x
            top = pt.y
            right = pt.x
            bottom = pt.y
            break
        break
    
    for poly in pts:
        for pt in poly:
            if pt.x < left: left = pt.x
            if pt.x > right: right = pt.x
            if pt.y < top: top = pt.y
            if pt.y > bottom: bottom = pt.y
    if left is None: return Rect(0, 0, 0, 0)
    else: return Rect(left, top, right, bottom)

def _GetLowestPt(poly):
    # precondition: poly must not be empty
    result = poly[0]
    for pt in poly:
        if pt.y > result.y or (pt.y == result.y and pt.x < result.x):
            result = pt
    return result

def _StripDupPts(poly):
    if poly == []: return poly
    for i in range(1, len(poly)):
        if _PointsEqual(poly[i-1], poly[i]): poly.pop(i)
    i = len(poly) -1
    while i > 0 and _PointsEqual(poly[i], poly[0]):
        poly.pop(i)
        i -= 1
    return poly

def _OffsetInternal(polys, isPolygon, delta, jointype = JoinType.Square, endtype = EndType.Square, limit = 0.0): 
    
    def _DoSquare(pt):
        pt1 = Point(round(pt.x + Normals[k].x * delta), round(pt.y + Normals[k].y * delta))
        pt2 = Point(round(pt.x + Normals[j].x * delta), round(pt.y + Normals[j].y * delta))
        if (Normals[k].x*Normals[j].y-Normals[j].x*Normals[k].y) * delta >= 0:
            a1 = math.atan2(Normals[k].y, Normals[k].x)
            a2 = math.atan2(-Normals[j].y, -Normals[j].x)
            a1 = abs(a2 - a1);
            if a1 > math.pi: a1 = math.pi * 2 - a1
            dx = math.tan((math.pi - a1)/4) * abs(delta)
            
            pt1 = Point(round(pt1.x -Normals[k].y * dx), round(pt1.y + Normals[k].x * dx))
            result.append(pt1)
            pt2 = Point(round(pt2.x + Normals[j].y * dx), round(pt2.y - Normals[j].x * dx))
            result.append(pt2)
        else:
            result.append(pt1)
            result.append(pt)
            result.append(pt2)

    def _DoMiter(pt, r):
        if ((Normals[k].x* Normals[j].y - Normals[j].x * Normals[k].y) * delta >= 0):
            q = delta / r;
            result.append(Point(round(pt.x + (Normals[k].x + Normals[j].x) *q),
              round(pt.y + (Normals[k].y + Normals[j].y) *q)))
        else:
            pt1 = Point(round(pt.x + Normals[k].x * delta), \
                        round(pt.y + Normals[k].y * delta))
            pt2 = Point(round(pt.x + Normals[j].x * delta), \
                        round(pt.y + Normals[j].y * delta))
            result.append(pt1)
            result.append(pt)
            result.append(pt2)

    def _DoRound(pt, limit):
        pt1 = Point(round(pt.x + Normals[k].x * delta), \
                    round(pt.y + Normals[k].y * delta))
        pt2 = Point(round(pt.x + Normals[j].x * delta), \
                    round(pt.y + Normals[j].y * delta))
        result.append(pt1)
        if (Normals[k].x * Normals[j].y - Normals[j].x * Normals[k].y) *delta >= 0:
            if (Normals[j].x * Normals[k].x + Normals[j].y * Normals[k].y) < 0.985:
                a1 = math.atan2(Normals[k].y, Normals[k].x)
                a2 = math.atan2(Normals[j].y, Normals[j].x)
                if (delta > 0) and (a2 < a1): a2 =