您的位置:首页 > 其它

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

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

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

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





一、生成随机多边形

<span style="font-size:18px;">#
import math
from collections import namedtuple
from decimal import Decimal, getcontext
import re
from random import randint
# from subprocess import call
import os

getcontext().prec = 8
horizontal = 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 result

def 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 None

def 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 None

def 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 = rightBound

class 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 = None

class OutPt(object):
__slots__ = ('idx','pt','prevOp','nextOp')
def __init__(self, idx, pt):
#从最外圈数进去第几个环,从-1开始计数
self.idx = idx
self.pt = pt
self.prevOp = None
self.nextOp = None

class 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 = None

class 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 result

def 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) / 2

def 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.deltaY

def _SwapSides(e1, e2):
side = e1.side
e1.side = e2.side
e2.side = side

def _SwapPolyIndexes(e1, e2):
idx = e1.outIdx
e1.outIdx = e2.outIdx
e2.outIdx = idx

def _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 = polyType

def _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 result

class 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.nextLm

class 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>




本节到此结束,欲知后事如何,请看下回分解。
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签: