-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathdelaunator_2d.py
100 lines (93 loc) · 3.45 KB
/
delaunator_2d.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
98
99
100
import taichi as ti
import taichi.math as tm
import argparse
def parse_args():
parser = argparse.ArgumentParser()
parser.add_argument('--device', type=str, default='gpu')
parser.add_argument('--case', type=int, default=0)
parser.add_argument('--dir', type=str, default='output')
parser.add_argument('--frame', type=int, default=100)
parser.add_argument('--iter', type=int, default=40)
return parser.parse_args()
cmd_args = parse_args()
ti.init(arch=ti.cuda if cmd_args.device == 'gpu' else ti.cpu, short_circuit_operators=True)
Nmax = 100000
N_neighbor = 800
chain_pre = ti.field(dtype=int, shape=(Nmax, N_neighbor))
chain_nxt = ti.field(dtype=int, shape=(Nmax, N_neighbor))
node_queue = ti.field(dtype=int, shape=(Nmax, N_neighbor * 3))
@ti.func
def angle_of_vec2(a: tm.vec2, b: tm.vec2) -> float:
res = 0.
a_len = tm.length(a)
b_len = tm.length(b)
if a_len > 0. and b_len > 0.:
res = tm.acos(tm.clamp(tm.dot(a, b) / (a_len * b_len), -1., 1.))
return res
@ti.func
def need_flip(cur_id: int, u: int, x: ti.template()) -> ti.i8:
res = 0
if angle_of_vec2(-x[cur_id, chain_pre[cur_id, u]], x[cur_id, u] - x[cur_id, chain_pre[cur_id, u]])\
+ angle_of_vec2(-x[cur_id, chain_nxt[cur_id, u]], x[cur_id, u] - x[cur_id, chain_nxt[cur_id, u]]) > tm.pi:
res = 1
return res
@ti.func
def get_local_mesh(cur_id: int, n_neighbor: ti.template(), x: ti.template(), ring_ids: ti.template()):
for i in range(n_neighbor[cur_id]):
ring_ids[cur_id, i] = i
# Sort by polar angle
while True:
has_swaped = 0
for i in range(n_neighbor[cur_id] - 1):
u = ring_ids[cur_id, i]
u_nxt = ring_ids[cur_id, i + 1]
angle = tm.atan2(x[cur_id, u][1], x[cur_id, u][0])
angle_nxt = tm.atan2(x[cur_id, u_nxt][1], x[cur_id, u_nxt][0])
if angle > angle_nxt or (angle == angle_nxt and tm.length(x[cur_id, u]) > tm.length(x[cur_id, u_nxt])):
tmp = ring_ids[cur_id, i + 1]
ring_ids[cur_id, i + 1] = ring_ids[cur_id, i]
ring_ids[cur_id, i] = tmp
has_swaped = 1
if has_swaped == 0:
break
# Build chain
for i in range(n_neighbor[cur_id] - 1):
chain_nxt[cur_id, ring_ids[cur_id, i]] = ring_ids[cur_id, i + 1]
chain_pre[cur_id, ring_ids[cur_id, i + 1]] = ring_ids[cur_id, i]
chain_nxt[cur_id, ring_ids[cur_id, n_neighbor[cur_id] - 1]] = ring_ids[cur_id, 0]
chain_pre[cur_id, ring_ids[cur_id, 0]] = ring_ids[cur_id, n_neighbor[cur_id] - 1]
# BFS for edge flip
q_st = 0
q_en = 0
for i in range(n_neighbor[cur_id]):
if need_flip(cur_id, ring_ids[cur_id, i], x):
node_queue[cur_id, q_en] = ring_ids[cur_id, i]
q_en += 1
while q_st < q_en:
u = node_queue[cur_id, q_st]
q_st += 1
if chain_nxt[cur_id, u] == -1 or not need_flip(cur_id, u, x):
continue
chain_nxt[cur_id, chain_pre[cur_id, u]] = chain_nxt[cur_id, u]
chain_pre[cur_id, chain_nxt[cur_id, u]] = chain_pre[cur_id, u]
if need_flip(cur_id, chain_nxt[cur_id, u], x):
node_queue[cur_id, q_en] = chain_nxt[cur_id, u]
q_en += 1
if need_flip(cur_id, chain_pre[cur_id, u], x):
node_queue[cur_id, q_en] = chain_pre[cur_id, u]
q_en += 1
chain_nxt[cur_id, u] = chain_pre[cur_id, u] = -1
# Get final ring and stores in ring_ids
u = -1
for _ in range(1):
for i in range(n_neighbor[cur_id]):
if chain_nxt[cur_id, i] != -1:
u = i
break
n_neighbor[cur_id] = 0
while True:
ring_ids[cur_id, n_neighbor[cur_id]] = u
n_neighbor[cur_id] += 1
u = chain_nxt[cur_id, u]
if u == ring_ids[cur_id, 0]:
break