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])
|