veclavtalica/mjestecko
1
1
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

push .static contents

Browse Source
This commit is contained in:
veclav talica
2024-02-10 22:38:35 +05:00
parent bb479c9371
commit 760f947c00
8 changed files with 482 additions and 1 deletions
Download Patch File Download Diff File Expand all files Collapse all files

286
articles/incremental-delaunay/.static/incremental-delaunay.zig.txt Normal file
View File

@ -0,0 +1,286 @@
//! Based on: https://www.cs.umd.edu/class/spring2020/cmsc754/Lects/lect13-delaun-alg.pdf
//! Optimizations involved:
//! - Cached neighbors for traversal.
//! - Minimal memory footprint.
//! - Cached circumferences.
//! - No circumference calculations for new subdivisions, - circumferences of neighbors are used instead.
//! - Lazy circumference calculation, as some places might not be neighboring new subdivisions.
//! - Extensive use of vectorization.
//! - Care given to linear access of memory.
// todo: This method allows zero area triangles, we need to eliminate them.
// Points that lie on edges can be detected in pointRelation function by == 0 comparison.
const std = @import("std");
// Could be redefined as pleased, but i consider these to be most sensical for given implementation.
pub const VertexComponent = f32;
pub const Vertex = @Vector(2, VertexComponent);
pub const Index = u15;
pub const Area = GenericArea(VertexComponent);
pub const Builder = struct {
triangles: std.ArrayList(Triangle),
vertices: std.ArrayList(Vertex),
allocator: std.mem.Allocator,
// todo: init with expected amount of points to preallocate beforehand.
pub fn init(allocator: std.mem.Allocator, area: Area) !@This() {
var triangles = try std.ArrayList(Triangle).initCapacity(allocator, 2);
errdefer triangles.deinit();
var vertices = try std.ArrayList(Vertex).initCapacity(allocator, 4);
errdefer vertices.deinit();
try vertices.ensureUnusedCapacity(4);
try triangles.ensureUnusedCapacity(2);
for (area.corners()) |corner|
vertices.append(corner) catch unreachable;
triangles.append(Triangle{
.points = [3]Index{ 0, 2, 1 },
.neighbors = [3]?Index{ null, 1, null },
}) catch unreachable;
triangles.append(Triangle{
.points = [3]Index{ 3, 1, 2 },
.neighbors = [3]?Index{ null, 0, null },
}) catch unreachable;
return .{
.triangles = triangles,
.vertices = vertices,
.allocator = allocator,
};
}
pub fn insertAtRandom(self: *@This(), point: Vertex, generator: std.rand.Random) !void {
// Find a triangle the point lies starting from some random triangle.
var abc_index: Index = @intCast(generator.int(Index) % self.triangles.items.len);
var abc = &self.triangles.items[abc_index];
var relation = abc.pointRelation(self.vertices, point);
while (relation != .contained) {
abc_index = abc.neighbors[@intCast(@intFromEnum(relation))].?;
abc = &self.triangles.items[abc_index];
relation = abc.pointRelation(self.vertices, point);
}
// Allocate two new triangles, as well as new vertex.
const new_vertex_index: Index = @intCast(self.vertices.items.len);
try self.vertices.append(point);
const pbc_index: Index = @intCast(self.triangles.items.len);
const apc_index: Index = @intCast(self.triangles.items.len + 1);
try self.triangles.ensureUnusedCapacity(2);
// Divide the abc triangle into three.
abc = &self.triangles.items[abc_index];
// Insert pbc.
self.triangles.append(Triangle{
.points = [3]Index{ new_vertex_index, abc.points[1], abc.points[2] },
.neighbors = [3]?Index{ abc_index, abc.neighbors[1], apc_index },
}) catch unreachable;
// Insert apc.
self.triangles.append(Triangle{
.points = [3]Index{ abc.points[0], new_vertex_index, abc.points[2] },
.neighbors = [3]?Index{ abc_index, pbc_index, abc.neighbors[2] },
}) catch unreachable;
// Update neighbors to be aware of new triangles.
inline for (abc.neighbors[1..], [2]Index{ pbc_index, apc_index }) |n, e|
if (n) |i| {
const p = &self.triangles.items[i];
p.neighbors[p.neighborPosition(abc_index)] = e;
};
// Existing abc is reused.
abc.points[2] = new_vertex_index;
abc.neighbors[1] = pbc_index;
abc.neighbors[2] = apc_index;
abc.circumference = null;
// Recursively adjust edges of triangles so that circumferences are only encasing 3 points at a time.
// todo: Try inlining initial calls via @call(.always_inline, ...).
self.trySwapping(abc_index, 0);
self.trySwapping(pbc_index, 1);
self.trySwapping(apc_index, 2);
}
fn trySwapping(self: @This(), triangle_index: Index, edge: u2) void {
// First find opposite to edge point that lies in neighbor.
const triangle = &self.triangles.items[triangle_index];
const neighbor_index = triangle.neighbors[edge];
if (neighbor_index == null)
return;
const neighbor = &self.triangles.items[neighbor_index.?];
if (neighbor.circumference == null)
neighbor.circumference = Triangle.Circumference.init(neighbor.*, self.vertices);
// Position of neighbor's point opposide to shared with triangle edge.
const point_order = neighbor.nextAfter(triangle.points[edge]);
const point_index = neighbor.points[point_order];
const prev_edge = if (edge == 0) 2 else edge - 1;
if (neighbor.doesFailIncircleTest(self.vertices.items[triangle.points[prev_edge]])) {
// Incircle test failed, swap edges of two triangles and then try swapping newly swapped ones.
const next_edge = (edge + 1) % 3;
const next_point_order = (point_order + 1) % 3;
const prev_point_order = if (point_order == 0) 2 else point_order - 1;
// Update neighbors of triangles in which edge was swapped.
if (triangle.neighbors[next_edge]) |i| {
const n = &self.triangles.items[i];
n.neighbors[n.neighborPosition(triangle_index)] = neighbor_index.?;
}
if (neighbor.neighbors[prev_point_order]) |i| {
const n = &self.triangles.items[i];
n.neighbors[n.neighborPosition(neighbor_index.?)] = triangle_index;
}
const neighbor_prev_point_order_neighbor_index_cache = neighbor.neighbors[prev_point_order];
neighbor.points[prev_point_order] = triangle.points[prev_edge];
neighbor.neighbors[next_point_order] = triangle.neighbors[next_edge];
neighbor.neighbors[prev_point_order] = triangle_index;
neighbor.circumference = null;
triangle.points[next_edge] = point_index;
triangle.neighbors[next_edge] = neighbor_index.?;
triangle.neighbors[edge] = neighbor_prev_point_order_neighbor_index_cache;
triangle.circumference = null;
self.trySwapping(triangle_index, edge);
self.trySwapping(neighbor_index.?, point_order);
}
}
};
const Triangle = struct {
// References to vertices it's composed of, named abc, in CCW orientation.
points: [3]Index,
// References to triangles that are on other side of any edge, if any.
// Order is: ab, bc, ca
neighbors: [3]?Index,
// Lazily calculated and cached for incircle tests.
circumference: ?Circumference = null,
pub const Circumference = struct {
center: Vertex,
radius_squared: VertexComponent, // todo: Way to get a type capable of holding squared values.
pub fn init(triangle: Triangle, vertices: std.ArrayList(Vertex)) @This() {
const a = vertices.items[triangle.points[0]];
const b = vertices.items[triangle.points[1]];
const c = vertices.items[triangle.points[2]];
const ab: Vertex = @splat(magnitudeSquared(a));
const cd: Vertex = @splat(magnitudeSquared(b));
const ef: Vertex = @splat(magnitudeSquared(c));
const cmb = @shuffle(VertexComponent, c - b, undefined, [2]i32{ 1, 0 });
const amc = @shuffle(VertexComponent, a - c, undefined, [2]i32{ 1, 0 });
const bma = @shuffle(VertexComponent, b - a, undefined, [2]i32{ 1, 0 });
const center = ((ab * cmb + cd * amc + ef * bma) / (a * cmb + b * amc + c * bma)) / @as(Vertex, @splat(2));
return .{
.center = center,
.radius_squared = magnitudeSquared(a - center),
};
}
};
// todo: Try perpendicular dot product approach.
pub fn pointRelation(self: @This(), vertices: std.ArrayList(Vertex), point: Vertex) enum(u2) {
outside_ab = 0,
outside_bc = 1,
outside_ca = 2,
contained = 3,
} {
const a = vertices.items[self.points[0]];
const b = vertices.items[self.points[1]];
const c = vertices.items[self.points[2]];
// https://stackoverflow.com/questions/1560492/how-to-tell-whether-a-point-is-to-the-right-or-left-side-of-a-line
const p = point;
// Calculate cross products for all edges at once.
const q = @Vector(12, VertexComponent){ b[0], b[1], c[0], c[1], a[0], a[1], p[1], p[0], p[1], p[0], p[1], p[0] };
const w = @Vector(12, VertexComponent){ a[0], a[1], b[0], b[1], c[0], c[1], a[1], a[0], b[1], b[0], c[1], c[0] };
const e = q - w;
const r = @shuffle(VertexComponent, e, undefined, [6]i32{ 0, 1, 2, 3, 4, 5 });
const t = @shuffle(VertexComponent, e, undefined, [6]i32{ 6, 7, 8, 9, 10, 11 });
const y = r * t;
const u = @shuffle(VertexComponent, y, undefined, [3]i32{ 0, 2, 4 });
const i = @shuffle(VertexComponent, y, undefined, [3]i32{ 1, 3, 5 });
const o = (u - i) > @Vector(3, VertexComponent){ 0, 0, 0 };
// const o = (u - i) <= @Vector(3, VertexComponent){ 0, 0, 0 };
// if (@reduce(.And, o))
// return .contained
// else if (!o[0])
// return .outside_ab
// else if (!o[1])
// return .outside_bc
// else
// return .outside_ca;
const mask = @as(u3, @intFromBool(o[0])) << 2 | @as(u3, @intFromBool(o[1])) << 1 | @as(u3, @intFromBool(o[2]));
return @enumFromInt(@clz(mask));
}
pub inline fn doesFailIncircleTest(self: @This(), point: Vertex) bool {
return magnitudeSquared(self.circumference.?.center - point) < self.circumference.?.radius_squared;
}
// todo: Shouldn't be here.
pub inline fn magnitudeSquared(p: Vertex) VertexComponent {
return @reduce(.Add, p * p);
}
// Finds which point comes after given one, by index, CCW.
// Used to translate point names when traveling between neighbors.
pub inline fn nextAfter(self: @This(), point_index: Index) u2 {
inline for (self.points, 0..) |p, i|
if (point_index == p)
return @intCast((i + 1) % 3);
unreachable;
}
pub inline fn neighborPosition(self: @This(), triangle_index: Index) usize {
inline for (self.neighbors, 0..) |n, i|
if (triangle_index == n)
return i;
unreachable;
}
};
pub fn GenericArea(comptime T: type) type {
return struct {
// note: Upper-left origin is assumed, if second point lies left or up of first it willn't work.
xyxy: @Vector(4, T),
/// Order: Upperleft, upperright, bottomleft, bottomright.
pub fn corners(self: @This()) [4]@Vector(2, T) {
return [4]@Vector(2, T){
@Vector(2, T){ self.xyxy[0], self.xyxy[1] },
@Vector(2, T){ self.xyxy[2], self.xyxy[1] },
@Vector(2, T){ self.xyxy[0], self.xyxy[3] },
@Vector(2, T){ self.xyxy[2], self.xyxy[3] },
};
}
};
}

BIN
articles/incremental-delaunay/.static/web.png (Stored with Git LFS) Normal file
View File

Binary file not shown.