mjestecko/articles/incremental-delaunay/.static/incremental-delaunay.zig.txt

//! 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] },
            };
        }
    };
}