From ae699565bb8cffa3e0f86b2db14a9f7a72e8c34e Mon Sep 17 00:00:00 2001 From: Stephen Gutekanst Date: Fri, 16 Sep 2022 10:17:37 -0700 Subject: [PATCH] trimesh2d: fix correctness issues, algo correctly implemented Signed-off-by: Stephen Gutekanst --- libs/trimesh2d/README.md | 2 +- libs/trimesh2d/src/main.zig | 148 +++++++++++++++++++++++------------- 2 files changed, 97 insertions(+), 53 deletions(-) diff --git a/libs/trimesh2d/README.md b/libs/trimesh2d/README.md index 8964d524..87183684 100644 --- a/libs/trimesh2d/README.md +++ b/libs/trimesh2d/README.md @@ -41,7 +41,7 @@ pub fn main() { var polygon = std.ArrayListUnmanaged(f32){}; // append your polygon vertices: - // try polygon.append(1.0); + // try polygon.append(allocator, 1.0); var out_triangles = std.ArrayListUnmanaged(u32){}; var processor = trimesh2d.Processor(f32){}; diff --git a/libs/trimesh2d/src/main.zig b/libs/trimesh2d/src/main.zig index c5ea4031..3732f28a 100644 --- a/libs/trimesh2d/src/main.zig +++ b/libs/trimesh2d/src/main.zig @@ -4,6 +4,8 @@ const std = @import("std"); /// (call reset between process calls.) The type T denotes e.g. f16, f32, or f64 vertices. pub fn Processor(comptime T: type) type { return struct { + const Vec2 = @Vector(2, T); + // Doubly linked list for fast polygon inspection prev: std.ArrayListUnmanaged(u32) = .{}, next: std.ArrayListUnmanaged(u32) = .{}, @@ -16,14 +18,14 @@ pub fn Processor(comptime T: type) type { /// Resets the processor, clearing the internal buffers and preparing it for processing a /// new polygon. - pub fn reset(self: *Processor) void { + pub fn reset(self: *@This()) void { self.prev.clearRetainingCapacity(); self.next.clearRetainingCapacity(); self.ears.clearRetainingCapacity(); self.is_ear.clearRetainingCapacity(); } - pub fn deinit(self: *Processor, allocator: std.mem.Allocator) void { + pub fn deinit(self: *@This(), allocator: std.mem.Allocator) void { self.prev.deinit(allocator); self.next.deinit(allocator); self.ears.deinit(allocator); @@ -32,96 +34,118 @@ pub fn Processor(comptime T: type) type { /// Processes a simple polygon (no holes) into triangles in linear time, writing the /// triangles to out_triangles (indices into polygon vertices list.) + /// + /// The polygons must be sorted in counter-clockwise order. pub fn process( - self: *Processor, + self: *@This(), allocator: std.mem.Allocator, - polygon: std.ArrayListUnmanaged(T), + // TODO(trimesh2d): make this a slice? + polygon: std.ArrayListUnmanaged(Vec2), out_triangles: *std.ArrayListUnmanaged(u32), ) error{OutOfMemory}!void { - if (polygon.len < 3) { + if (polygon.items.len < 3) { return; } // Ensure our doubly linked list and ears list are large enough. - const size = polygon.len; - try self.prev.ensureTotalCapacity(allocator, size); - try self.next.ensureTotalCapacity(allocator, size); + const size = polygon.items.len; + try self.prev.resize(allocator, size); + try self.next.resize(allocator, size); try self.ears.ensureTotalCapacity(allocator, size); try self.is_ear.resize(allocator, size); + for (self.is_ear.items) |*v| v.* = false; // Fill prev list with prior-index values, e.g.: // [4, 0, 1, 2, 3] - for (self.prev.items) |_, i| self.prev.items[i] = if (i == 0) size - 1 else i - 1; + for (self.prev.items) |_, i| self.prev.items[i] = @intCast(u32, if (i == 0) size - 1 else i - 1); // Fill next list with next-index values, e.g.: // [1, 2, 3, 4, 0] - for (self.prev.items) |_, i| self.prev.items[i] = if (i == self.prev.items.len - 1) size - 1 else i + 1; + for (self.next.items) |_, i| self.next.items[i] = @intCast(u32, if (i == size - 1) 0 else i + 1); // Detect all safe ears in O(n). // This amounts to finding all convex vertices but the endpoints of the constrained edge var curr: u32 = 1; - while (cur < size - 1) : (curr += 1) { - // NOTE: the polygon may contain dangling edges, so !arrayListElementsEqual(prev, next) - // avoids need to even do the more expensive ear test for them below. - if (!arrayListElementsEqual(prev, next) and orient2d( - // TODO(trimesh2d): all this code would be simpler if we had a poly index helper - // which returned a @Vector2(2, T) - @Vector(2, T){ poly[self.prev.items[curr]], poly[self.prev.items[curr] + 1] }, - @Vector(2, T){ poly[curr], poly[curr + 1] }, - @Vector(2, T){ poly[self.next.items[curr]], poly[self.next.items[curr] + 1] }, - )) { - try self.ears.append(curr); + while (curr < size - 1) : (curr += 1) { + // NOTE: the polygon may contain dangling edges + if (orient2d( + polygon.items[self.prev.items[curr]], + polygon.items[curr], + polygon.items[self.next.items[curr]], + ) > 0) { + try self.ears.append(allocator, curr); self.is_ear.items[curr] = true; } } // Progressively delete all ears, updating the data structure - const length = size; - while (true) { - const curr = self.ears.pop(); + while (self.ears.items.len > 0) { + curr = self.ears.pop(); // make new tri - try out_triangles.append(self.prev.items[curr]); - try out_triangles.append(curr); - try out_triangles.append(self.next.items[curr]); + try out_triangles.append(allocator, self.prev.items[curr]); + try out_triangles.append(allocator, curr); + try out_triangles.append(allocator, self.next.items[curr]); // exclude curr from the polygon, connecting prev and next self.next.items[self.prev.items[curr]] = self.next.items[curr]; self.prev.items[self.next.items[curr]] = self.prev.items[curr]; - length -= 1; - if (length < 3) return; // last triangle - // check if prev and next have become new ears if (!self.is_ear.items[self.prev.items[curr]] and self.prev.items[curr] != 0) { if (self.prev.items[self.prev.items[curr]] != self.next.items[curr] and orient2d( - @Vector(2, T){ poly[self.prev.items[self.prev.items[curr]]], poly[(self.prev.items[self.prev.items[curr]]) + 1] }, - @Vector(2, T){ poly[self.prev.items[curr]], poly[self.prev.items[curr] + 1] }, - @Vector(2, T){ poly[self.next.items[curr]], poly[self.next.items[curr] + 1] }, + polygon.items[self.prev.items[self.prev.items[curr]]], + polygon.items[self.prev.items[curr]], + polygon.items[self.next.items[curr]], ) > 0) { - try self.ears.append(self.prev.items[curr]); + try self.ears.append(allocator, self.prev.items[curr]); self.is_ear.items[self.prev.items[curr]] = true; } } if (!self.is_ear.items[self.next.items[curr]] and self.next.items[curr] < size - 1) { if (self.next.items[self.next.items[curr]] != self.prev.items[curr] and orient2d( - @Vector(2, T){ poly[self.prev.items[curr]], poly[(self.prev.items[curr]) + 1] }, - @Vector(2, T){ poly[self.next.items[curr]], poly[self.next.items[curr] + 1] }, - @Vector(2, T){ poly[self.next.items[self.next.items[curr]]], poly[self.next.items[self.next.items[curr]] + 1] }, + polygon.items[self.prev.items[curr]], + polygon.items[self.next.items[curr]], + polygon.items[self.next.items[self.next.items[curr]]], ) > 0) { - try self.ears.append(self.next.items[curr]); + try self.ears.append(allocator, self.next.items[curr]); self.is_ear.items[self.next.items[curr]] = true; } } } } + pub fn sort(allocator: std.mem.Allocator, polygon: std.ArrayListUnmanaged(Vec2)) !std.ArrayListUnmanaged(Vec2) { + var max_dist: f32 = 0; + var extrema_start: usize = undefined; + var extrema_end: usize = undefined; + var i: usize = 0; + while (i < polygon.items.len) : (i += 1) { + var next_index = (i + 1) % polygon.items.len; + var p0 = polygon.items[i]; + var p1 = polygon.items[next_index]; + var dist = std.math.hypot(T, p1[0] - p0[0], p1[1] - p0[1]); + if (dist > max_dist) { + max_dist = dist; + extrema_start = i; + extrema_end = next_index; + } + } + + var sorted = std.ArrayListUnmanaged(Vec2){}; + i = extrema_end; + while (i < polygon.items.len) : (i += 1) try sorted.append(allocator, polygon.items[i]); + i = 0; + while (i <= extrema_start) : (i += 1) try sorted.append(allocator, polygon.items[i]); + return sorted; + } + /// Inexact geometric predicate. /// Basically Shewchuk's orient2dfast() fn orient2d( - pa: @Vector(2, T), - pb: @Vector(2, T), - pc: @Vector(2, T), + pa: Vec2, + pb: Vec2, + pc: Vec2, ) T { const acx = pa[0] - pc[0]; const bcx = pb[0] - pc[0]; @@ -129,20 +153,40 @@ pub fn Processor(comptime T: type) type { const bcy = pb[1] - pc[1]; return acx * bcy - acy * bcx; } - - fn arrayListElementsEqual( - a: std.ArrayListUnmanaged(u32), - b: std.ArrayListUnmanaged(u32), - ) bool { - if (a.len != b.len) return false; - for (a.items) |aa, i| { - if (b.items[i] != aa) return false; - } - return true; - } }; } +test "simple" { + const allocator = std.testing.allocator; + const Vec2 = @Vector(2, f32); + + var polygon = std.ArrayListUnmanaged(Vec2){}; + defer polygon.deinit(allocator); + // CCW + try polygon.append(allocator, Vec2{ 0.0, 0.0 }); // bottom-left + try polygon.append(allocator, Vec2{ 1.0, 0.0 }); // bottom-right + try polygon.append(allocator, Vec2{ 1.0, 1.0 }); // top-right + try polygon.append(allocator, Vec2{ 0.0, 1.0 }); // top-left + + var out_triangles = std.ArrayListUnmanaged(u32){}; + defer out_triangles.deinit(allocator); + var processor = Processor(f32){}; + defer processor.deinit(allocator); + + // Process a polygon. + try processor.process(allocator, polygon, &out_triangles); + + // out_triangles has indices into polygon.items of our triangle vertices. + // If desired, call .reset() and call .process() again! Internal buffers will be reused. + try std.testing.expectEqual(@as(usize, 6), out_triangles.items.len); + try std.testing.expectEqual(@as(u32, 1), out_triangles.items[0]); // bottom-right + try std.testing.expectEqual(@as(u32, 2), out_triangles.items[1]); // top-right + try std.testing.expectEqual(@as(u32, 3), out_triangles.items[2]); // top-left + try std.testing.expectEqual(@as(u32, 0), out_triangles.items[3]); // bottom-left + try std.testing.expectEqual(@as(u32, 1), out_triangles.items[4]); // bottom-right + try std.testing.expectEqual(@as(u32, 3), out_triangles.items[5]); // top-left +} + test { std.testing.refAllDeclsRecursive(@This()); }