【发布时间】:2015-11-10 15:09:03
【问题描述】:
我正在使用 Triangle 模块生成受约束的 Delaunay 三角剖分(可在 http://www.lfd.uci.edu/~gohlke/pythonlibs/ 获得)。在某些情况下,函数 triangle.triangulate 崩溃,并且 Windows 显示“Python.exe 已停止响应”。我尝试过使用 try/except 结构,如下例所示,它会不一致地崩溃。我假设三角测量过程的一部分是随机的(参见下面的“次要问题”),但我并不完全确定。
这种不一致相当令人担忧。我在 Windows 上。
from shapely.geometry import Polygon, MultiPolygon, LineString
import numpy as np
import triangle
import matplotlib.pyplot as plt
import matplotlib.tri as tri
import random
def meshXOR(polyRef, shape, otherVerts, otherSegs, otherHole):
verts = []
verts3 = []
segs = []
outerLength = len(polyRef[shape[0]])
for i in range(outerLength-1):#-1 because xorpoly duplicates first point into last spot
verts.append((polyRef[shape[0]][i][0],polyRef[shape[0]][i][1])) #append the point to the verts array which will be fed into the delaunay triangulator
if i == outerLength - 2:
segs.append([i,0])
else:
segs.append([i,i+1])
h = []
for cInd in shape[1]:
shift = len(verts)
innerLength = len(polyRef[cInd])
for i in range(innerLength-1):
verts.append((polyRef[cInd][i][0],polyRef[cInd][i][1]))
if i == innerLength - 2:
segs.append([i+shift,shift])
else:
segs.append([i+shift,i+1+shift])
h += list(Polygon(polyRef[cInd]).representative_point().coords)
print 'verts are: ', verts
#output: verts are: [(0.0, 5.0), (0.0, 10.0), (10.0, 10.0), (10.0, 0.0), (0.0, 0.0), (0.0, 5.0), (7.0, 3.0), (7.0, 7.0)]
print 'segs are: ', segs
#output: segs are: [[0, 1], [1, 2], [2, 3], [3, 4], [4, 0], [5, 6], [6, 7], [7, 5]]
print 'holes are: ', h
#output: holes are: [(5.25, 6.0)]
print 'verts: ', verts == otherVerts
print 'segs: ', segs == otherSegs
print 'hole: ', h == otherHole
return verts, segs, h
pA = Polygon([[0.0,0.0],[10.0,0.0],[10.0,10.0],[0.0,10.0]])
pB = Polygon([[0.0,5.0],[7.0,3.0],[7.0,7.0]])
xorPoly = pA.symmetric_difference(pB)
if xorPoly.geom_type == 'Polygon': xorPoly = MultiPolygon([xorPoly])
otherVerts = [(0.0, 5.0), (0.0, 10.0), (10.0, 10.0), (10.0, 0.0), (0.0, 0.0), (0.0, 5.0), (7.0, 3.0), (7.0, 7.0)]
otherSegs = [[0, 1], [1, 2], [2, 3], [3, 4], [4, 0], [5, 6], [6, 7], [7, 5]]
otherHole = [(5.25,6.0)]
xorPolys = []
shapes = []
for poly in xorPoly:
shapes.append([len(xorPolys), [], len(shapes)])
xorPolys.append(list(poly.exterior.coords))
for ip in poly.interiors:
shapes[-1][1].append(len(xorPolys))
xorPolys.append(list(ip.coords))
try:
verts, segs, holes = meshXOR(xorPolys, shapes[0], otherVerts, otherSegs, otherHole) # i even tried placing it here
except:
print 'failed'
if len(holes)>0:
A = dict(vertices = np.asarray(verts), segments = np.asarray(segs), holes = holes)
else:
A = dict(vertices = np.asarray(verts), segments = np.asarray(segs))
print 'about to tri'
B = triangle.triangulate(A, opts = 'pi') #this is the step that try/except doesn't work on
print 'completed tri'
try:
B_t = B["triangles"].tolist()
except:
print 'no trianlges'
if B_t != []:
cols = []
import random
for tri in B_t:
cols.append(random.random())
plt.figure()
plt.gca().set_aspect('equal')
xy = np.asarray(verts)
plt.tripcolor(xy[:,0], xy[:,1], B_t, facecolors=np.array(cols))
#for tri in B_t:
#print 'tri is: ', [verts[t] for t in tri]
plt.show()
else:
print 'no triangles'
我的问题:有没有办法做类似“尝试/排除”结构的事情来捕捉这个错误?或者,我对三角形模块做错了吗?
编辑:
解决方案:Gamrix 的 回复中的引用引出了解决方案。如果点之间的差异太小(欧几里得距离),则三角函数会崩溃。删除间隔小于 1e-12 的点解决了这个问题。
【问题讨论】:
标签: python crash triangulation delaunay