-
Notifications
You must be signed in to change notification settings - Fork 4
/
isoline.py
252 lines (219 loc) · 10.6 KB
/
isoline.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
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
from __future__ import annotations
from dataclasses import dataclass
import numpy as np
from .cell import Cell, build_tree
from .point import Func, Point, ValuedPoint, binary_search_zero
def plot_isoline(
fn: Func,
pmin: Point,
pmax: Point,
min_depth: int = 5,
max_quads: int = 10000,
tol: np.ndarray | None = None,
) -> list[list[Point]]:
"""Get the curve representing fn([x,y])=0 on pmin[0] ≤ x ≤ pmax[0] ∩ pmin[1] ≤ y ≤ pmax[1]
Returns as a list of curves, where each curve is a list of points"""
pmin = np.asarray(pmin)
pmax = np.asarray(pmax)
if tol is None:
tol = (pmax - pmin) / 1000
else:
tol = np.asarray(tol)
quadtree = build_tree(2, fn, pmin, pmax, min_depth, max_quads, tol)
triangles = Triangulator(quadtree, fn, tol).triangulate()
return CurveTracer(triangles, fn, tol).trace()
@dataclass
class Triangle:
vertices: list[ValuedPoint]
""" The order of triangle "next" is such that, when walking along the isoline in the direction of next,
you keep positive function values on your right and negative function values on your left."""
next: Triangle | None = None
next_bisect_point: ValuedPoint | None = None
prev: Triangle | None = None
visited: bool = False
def four_triangles(
a: ValuedPoint, b: ValuedPoint, c: ValuedPoint, d: ValuedPoint, center: ValuedPoint
) -> tuple[Triangle, Triangle, Triangle, Triangle]:
"""a,b,c,d should be clockwise oriented, with center on the inside of that quad"""
return (
Triangle([a, b, center]),
Triangle([b, c, center]),
Triangle([c, d, center]),
Triangle([d, a, center]),
)
class Triangulator:
"""While triangulating, also compute the isolines.
Divides each quad into 8 triangles from the quad's center. This simplifies
adjacencies between triangles for the general case of multiresolution quadtrees.
Based on Manson, Josiah, and Scott Schaefer. "Isosurfaces
over simplicial partitions of multiresolution grids." Computer Graphics Forum.
Vol. 29. No. 2. Oxford, UK: Blackwell Publishing Ltd, 2010.
(https://people.engr.tamu.edu/schaefer/research/iso_simplicial.pdf), but this
does not currently implement placing dual vertices based on the gradient.
"""
def __init__(self, root: Cell, fn: Func, tol: np.ndarray) -> None:
self.triangles: list[Triangle] = []
self.hanging_next: dict[bytes, Triangle] = {}
self.root = root
self.fn = fn
self.tol = tol
def triangulate(self) -> list[Triangle]:
self.triangulate_inside(self.root)
return self.triangles
def triangulate_inside(self, quad: Cell) -> None:
if quad.children:
for child in quad.children:
self.triangulate_inside(child)
self.triangulate_crossing_row(quad.children[0], quad.children[1])
self.triangulate_crossing_row(quad.children[2], quad.children[3])
self.triangulate_crossing_col(quad.children[0], quad.children[2])
self.triangulate_crossing_col(quad.children[1], quad.children[3])
def triangulate_crossing_row(self, a: Cell, b: Cell) -> None:
"""Quad b should be to the right (greater x values) than quad a"""
if a.children and b.children:
self.triangulate_crossing_row(a.children[1], b.children[0])
self.triangulate_crossing_row(a.children[3], b.children[2])
elif a.children:
self.triangulate_crossing_row(a.children[1], b)
self.triangulate_crossing_row(a.children[3], b)
elif b.children:
self.triangulate_crossing_row(a, b.children[0])
self.triangulate_crossing_row(a, b.children[2])
else:
# a and b are minimal 2-cells
face_dual_a = self.get_face_dual(a)
face_dual_b = self.get_face_dual(b)
# Add the four triangles from the centers of a and b to the shared edge between them
if a.depth < b.depth:
# b is smaller
edge_dual = self.get_edge_dual(b.vertices[2], b.vertices[0])
triangles = four_triangles(b.vertices[2], face_dual_b, b.vertices[0], face_dual_a, edge_dual)
else:
edge_dual = self.get_edge_dual(a.vertices[3], a.vertices[1])
triangles = four_triangles(a.vertices[3], face_dual_b, a.vertices[1], face_dual_a, edge_dual)
self.add_four_triangles(triangles)
def triangulate_crossing_col(self, a: Cell, b: Cell) -> None:
"""Mostly a copy-paste of triangulate_crossing_row. For n-dimensions, want to pass a
dir index into a shared triangulate_crossing_dir function instead"""
if a.children and b.children:
self.triangulate_crossing_col(a.children[2], b.children[0])
self.triangulate_crossing_col(a.children[3], b.children[1])
elif a.children:
self.triangulate_crossing_col(a.children[2], b)
self.triangulate_crossing_col(a.children[3], b)
elif b.children:
self.triangulate_crossing_col(a, b.children[0])
self.triangulate_crossing_col(a, b.children[1])
else:
# a and b are minimal 2-cells
face_dual_a = self.get_face_dual(a)
face_dual_b = self.get_face_dual(b)
# Add the four triangles from the centers of a and b to the shared edge between them
if a.depth < b.depth:
# b is smaller
edge_dual = self.get_edge_dual(b.vertices[0], b.vertices[1])
triangles = four_triangles(b.vertices[0], face_dual_b, b.vertices[1], face_dual_a, edge_dual)
else:
edge_dual = self.get_edge_dual(a.vertices[2], a.vertices[3])
triangles = four_triangles(a.vertices[2], face_dual_b, a.vertices[3], face_dual_a, edge_dual)
self.add_four_triangles(triangles)
def add_four_triangles(self, triangles: tuple[Triangle, Triangle, Triangle, Triangle]) -> None:
for i in range(4):
self.next_sandwich_triangles(triangles[i], triangles[(i + 1) % 4], triangles[(i + 2) % 4])
self.triangles.extend(triangles)
def set_next(self, tri1: Triangle, tri2: Triangle, vpos: ValuedPoint, vneg: ValuedPoint) -> None:
if not vpos.val > 0 >= vneg.val:
return
intersection, is_zero = binary_search_zero(vpos, vneg, self.fn, self.tol)
if not is_zero:
return
tri1.next_bisect_point = intersection
tri1.next = tri2
tri2.prev = tri1
def next_sandwich_triangles(self, a: Triangle, b: Triangle, c: Triangle) -> None:
"""Find the "next" triangle for the triangle b. See Triangle for a description of the curve orientation.
We assume the triangles are oriented such that they share common vertices center←a[2]≡b[2]≡c[2]
and x←a[1]≡b[0], y←b[1]≡c[0]"""
center = b.vertices[2]
x = b.vertices[0]
y = b.vertices[1]
# Simple connections: inside the same four triangles
# (Group 0 with negatives)
if center.val > 0 >= y.val:
self.set_next(b, c, center, y)
# (Group 0 with negatives)
if x.val > 0 >= center.val:
self.set_next(b, a, x, center)
# More difficult connections: complete a hanging connection
# or wait for another triangle to complete this
# We index using (double) the midpoint of the hanging edge
id = (x.pos + y.pos).data.tobytes()
# (Group 0 with negatives)
if y.val > 0 >= x.val:
if id in self.hanging_next:
self.set_next(b, self.hanging_next[id], y, x)
del self.hanging_next[id]
else:
self.hanging_next[id] = b
elif y.val <= 0 < x.val:
if id in self.hanging_next:
self.set_next(self.hanging_next[id], b, x, y)
del self.hanging_next[id]
else:
self.hanging_next[id] = b
def get_edge_dual(self, p1: ValuedPoint, p2: ValuedPoint) -> ValuedPoint:
"""Returns the dual point on an edge p1--p2"""
if (p1.val > 0) != (p2.val > 0):
# The edge crosses the isoline, so take the midpoint
return ValuedPoint.midpoint(p1, p2, self.fn)
dt = 0.01
# We intersect the planes with normals <∇f(p1), -1> and <∇f(p2), -1>
# move slightly from p1 to p2. df = ∆f, so ∆f/∆t = 100*df1 near p1
df1 = self.fn(p1.pos * (1 - dt) + p2.pos * dt)
# move slightly from p2 to p1. df = ∆f, so ∆f/∆t = -100*df2 near p2
df2 = self.fn(p1.pos * dt + p2.pos * (1 - dt))
# (Group 0 with negatives)
if (df1 > 0) == (df2 > 0):
# The function either increases → ← or ← →, so a lerp would shoot out of bounds
# Take the midpoint
return ValuedPoint.midpoint(p1, p2, self.fn)
else:
# Increases → 0 → or ← 0 ←
v1 = ValuedPoint(p1.pos, df1)
v2 = ValuedPoint(p2.pos, df2)
return ValuedPoint.intersectZero(v1, v2, self.fn)
def get_face_dual(self, quad: Cell) -> ValuedPoint:
# TODO: proper face dual
return ValuedPoint.midpoint(quad.vertices[0], quad.vertices[-1], self.fn)
class CurveTracer:
active_curve: list[ValuedPoint]
def __init__(self, triangles: list[Triangle], fn: Func, tol: np.ndarray) -> None:
self.triangles = triangles
self.fn = fn
self.tol = tol
def trace(self) -> list[list[Point]]:
curves: list[list[ValuedPoint]] = []
for triangle in self.triangles:
if not triangle.visited and triangle.next is not None:
self.active_curve = []
self.march_triangle(triangle)
# triangle.next is not None, so there should be at least one segment
curves.append(self.active_curve)
return [[v.pos for v in curve] for curve in curves]
def march_triangle(self, triangle: Triangle) -> None:
start_triangle = triangle
closed_loop = False
# Iterate backwards to the start of a connected curve
while triangle.prev is not None:
triangle = triangle.prev
if triangle is start_triangle:
closed_loop = True
break
while triangle is not None and not triangle.visited:
if triangle.next_bisect_point is not None:
self.active_curve.append(triangle.next_bisect_point)
triangle.visited = True
triangle = triangle.next
if closed_loop:
# close back the loop
self.active_curve.append(self.active_curve[0])