-
Notifications
You must be signed in to change notification settings - Fork 0
/
get_triangulation.py
97 lines (80 loc) · 2.81 KB
/
get_triangulation.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
"""
Currently creates random vertices and shows a possible triangulation for them.
"""
import numpy as np
import matplotlib.pyplot as plt
from random import random, seed
from scipy.spatial import Delaunay
from itertools import combinations
def area(t, points):
# Calculate area of triangle from its vertices
p0 = points[t[0]]
p1 = points[t[1]]
p2 = points[t[2]]
A = abs(p0[0]*(p1[1]-p2[1]) + p1[0]*(p2[1]-p0[1]) + p2[0]*(p0[1]-p1[1]))/2
return A
def side_lengths(t, points):
sl = []
for ind in range(len(t)):
p1 = points[t[ind]]
p2 = points[t[(ind+1) % len(t)]]
s = np.sqrt((p1[0]-p2[0])**2+(p1[1]-p2[1])**2)
sl.append(s)
return sl
def check_for_small_tri(tri, points, bad_const):
bad_tri = []
for t in tri:
A = area(t, points)
s = side_lengths(t, points)
check = A/max(s)**2
if check < bad_const:
bad_tri.append(t)
return bad_tri
def get_good_triangles(n, rand_seed=None, bad_const=0.03):
"""
Make sure the ratio of the triangles area to the length of its longest side
squared is greater than some constant.
"""
if rand_seed is not None:
seed(rand_seed)
x = [random() for i in range(n)]
y = [random() for i in range(n)]
points = np.array(list(zip(x, y)))
tri = Delaunay(points).simplices
bad_tri = check_for_small_tri(tri, points, bad_const)
check_removal = set()
for t in bad_tri:
check_removal |= set(t)
# See if removing a single point from the bad triangles will make it so there are none.
worked = len(bad_tri) == 0
for r in check_removal:
check_points = np.concatenate((points[:r], points[r+1:]))
check_tri = Delaunay(check_points).simplices
if len(check_for_small_tri(check_tri, check_points, bad_const)) == 0:
points = check_points
tri = check_tri
worked = True
break
for rem_num in range(2, len(check_removal)):
if not worked:
for r in combinations(check_removal, 2):
# check_points = np.concatenate((points[:r], points[r+1:]))
check_points = np.delete(points, r, axis=0)
check_tri = Delaunay(check_points).simplices
if len(check_for_small_tri(check_tri, check_points, bad_const)) == 0:
points = check_points
tri = check_tri
worked = True
break
if worked:
break
if not worked:
print('Wasnt able to fix')
return points, tri
if __name__ == '__main__':
points, tri = get_good_triangles(30, rand_seed=154)
plt.triplot(points[:, 0], points[:, 1], tri)
plt.plot(points[:, 0], points[:, 1], 'o')
for i, p in enumerate(points):
plt.annotate(i, p)
plt.show()