delaunay_triangulation.cpp

#include "delaunay_triangulation.hpp"
#include <algorithm>

bool point_within_triangle(const glm::vec2 &p, const glm::vec2 &a, const glm::vec2 &b, const glm::vec2 &c) {
    auto sign = [](const glm::vec2 &p1, const glm::vec2 &p2, const glm::vec2 &p3) {
        return (p1.x - p3.x) * (p2.y - p3.y) - (p2.x - p3.x) * (p1.y - p3.y);
    };

    float d1 = sign(p, a, b);
    float d2 = sign(p, b, c);
    float d3 = sign(p, c, a);

    bool has_neg = (d1 < 0) || (d2 < 0) || (d3 < 0);
    bool has_pos = (d1 > 0) || (d2 > 0) || (d3 > 0);

    return !(has_neg && has_pos);
}

std::pair<float, glm::vec2> normalize_points(std::vector<glm::vec2> &points) {
    if (points.empty())
        return {1.0f, glm::vec2(0, 0)};

    float xmin = points[0].x, xmax = points[0].x;
    float ymin = points[0].y, ymax = points[0].y;

    for (const auto &p : points) {
        xmin = std::min(xmin, p.x);
        xmax = std::max(xmax, p.x);
        ymin = std::min(ymin, p.y);
        ymax = std::max(ymax, p.y);
    }

    float width = xmax - xmin;
    float height = ymax - ymin;
    float scale = std::max(width, height);

    if (scale > 0) {
        for (auto &p : points) {
            p.x = (p.x - xmin) / scale;
            p.y = (p.y - ymin) / scale;
        }
    }

    return {scale, glm::vec2(xmin, ymin)};
}

void sort_points_by_spatial_proximity(std::vector<glm::vec2> &points) {
    int n = points.size();
    int n_bin_rows = static_cast<int>(std::ceil(std::pow(n, 0.25)));

    std::vector<std::pair<int, glm::vec2>> indexed;
    indexed.reserve(n);

    for (const auto &p : points) {
        int row = static_cast<int>(p.y * n_bin_rows * 0.999f);
        int col = static_cast<int>(p.x * n_bin_rows * 0.999f);
        int bin = (row % 2) ? (row + 1) * n_bin_rows - col : row * n_bin_rows + col + 1;
        indexed.emplace_back(bin, p);
    }

    std::sort(indexed.begin(), indexed.end(), [](const auto &a, const auto &b) { return a.first < b.first; });

    for (size_t i = 0; i < points.size(); ++i) {
        points[i] = indexed[i].second;
    }
}

int find_containing_triangle(const glm::vec2 &p, int start_idx, const std::vector<glm::vec2> &points,
                             const std::vector<TriangleAdj> &triangles) {
    int j = start_idx;

    while (true) {
        const auto &tri = triangles[j];
        const glm::vec2 &v0 = points[tri.vertices[0]];
        const glm::vec2 &v1 = points[tri.vertices[1]];
        const glm::vec2 &v2 = points[tri.vertices[2]];

        if (point_within_triangle(p, v0, v1, v2)) {
            return j;
        }

        // move towards the point
        glm::vec2 ab = v1 - v0;
        glm::vec2 bc = v2 - v1;
        glm::vec2 ca = v0 - v2;
        glm::vec2 ap = p - v0;
        glm::vec2 bp = p - v1;
        glm::vec2 cp = p - v2;

        float s1 = ap.x * ab.y - ap.y * ab.x;
        float s2 = bp.x * bc.y - bp.y * bc.x;
        float s3 = cp.x * ca.y - cp.y * ca.x;

        if (s1 > 0 && s1 >= s2 && s1 >= s3)
            j = tri.adjacencies[0];
        else if (s2 > 0 && s2 >= s1 && s2 >= s3)
            j = tri.adjacencies[1];
        else if (s3 > 0 && s3 >= s1 && s3 >= s2)
            j = tri.adjacencies[2];
    }
}

bool in_circumcircle(const glm::vec2 &p, const glm::vec2 &v1, const glm::vec2 &v2, const glm::vec2 &v3) {
    float cosa = (v1.x - v3.x) * (v2.x - v3.x) + (v1.y - v3.y) * (v2.y - v3.y);
    float cosb = (v2.x - p.x) * (v1.x - p.x) + (v2.y - p.y) * (v1.y - p.y);
    float sina = (v1.x - v3.x) * (v2.y - v3.y) - (v1.y - v3.y) * (v2.x - v3.x);
    float sinb = (v2.x - p.x) * (v1.y - p.y) - (v2.y - p.y) * (v1.x - p.x);

    return (cosa < 0 && cosb < 0) || (-cosa * sinb > cosb * sina);
}

DelaunayResult delaunay_triangulate(std::vector<glm::vec2> points) {
    if (points.size() < 3) {
        return DelaunayResult{points, {}};
    }

    int num_original_points = points.size();

    // normalize and sort points
    auto [scale, offset] = normalize_points(points);
    sort_points_by_spatial_proximity(points);

    // add super-triangle
    points.push_back(glm::vec2(-100, -100));
    points.push_back(glm::vec2(100, -100));
    points.push_back(glm::vec2(0, 100));

    // initialize with super-triangle
    std::vector<TriangleAdj> triangles;
    TriangleAdj super_tri;
    super_tri.vertices = {num_original_points, num_original_points + 1, num_original_points + 2};
    triangles.push_back(super_tri);

    // insert points one by one
    for (int i = 0; i < num_original_points; ++i) {
        int containing_tri = find_containing_triangle(points[i], triangles.size() - 1, points, triangles);

        // split containing triangle into three
        TriangleAdj &old_tri = triangles[containing_tri];
        int v0 = old_tri.vertices[0];
        int v1 = old_tri.vertices[1];
        int v2 = old_tri.vertices[2];
        int adj0 = old_tri.adjacencies[0];
        int adj1 = old_tri.adjacencies[1];
        int adj2 = old_tri.adjacencies[2];

        // create two new triangles
        TriangleAdj tri1, tri2;
        tri1.vertices = {i, v1, v2};
        tri2.vertices = {i, v2, v0};

        int idx1 = triangles.size();
        int idx2 = triangles.size() + 1;

        // update adjacencies
        if (adj0 >= 0) {
            for (int m = 0; m < 3; ++m) {
                if (triangles[adj0].adjacencies[m] == containing_tri) {
                    triangles[adj0].adjacencies[m] = containing_tri;
                    break;
                }
            }
        }
        if (adj1 >= 0) {
            for (int m = 0; m < 3; ++m) {
                if (triangles[adj1].adjacencies[m] == containing_tri) {
                    triangles[adj1].adjacencies[m] = idx1;
                    break;
                }
            }
        }
        if (adj2 >= 0) {
            for (int m = 0; m < 3; ++m) {
                if (triangles[adj2].adjacencies[m] == containing_tri) {
                    triangles[adj2].adjacencies[m] = idx2;
                    break;
                }
            }
        }

        tri1.adjacencies = {containing_tri, adj1, idx2};
        tri2.adjacencies = {idx1, adj2, containing_tri};

        old_tri.vertices = {i, v0, v1};
        old_tri.adjacencies = {idx2, adj0, idx1};

        triangles.push_back(tri1);
        triangles.push_back(tri2);

        // lawson flip algorithm
        std::vector<int> stack;
        if (old_tri.adjacencies[1] >= 0)
            stack.push_back(containing_tri);
        if (tri1.adjacencies[1] >= 0)
            stack.push_back(idx1);
        if (tri2.adjacencies[1] >= 0)
            stack.push_back(idx2);

        while (!stack.empty()) {
            int L = stack.back();
            stack.pop_back();

            TriangleAdj &tri_L = triangles[L];
            int R = tri_L.adjacencies[1];
            if (R < 0)
                continue;

            // find opposite vertex
            int opp_vert = -1, opp_vert_id = -1;
            for (int k = 0; k < 3; ++k) {
                int v = triangles[R].vertices[k];
                if (v != tri_L.vertices[1] && v != tri_L.vertices[2]) {
                    opp_vert = v;
                    opp_vert_id = k;
                    break;
                }
            }

            if (opp_vert < 0)
                continue;

            // check circumcircle condition
            if (in_circumcircle(points[i], points[tri_L.vertices[2]], points[tri_L.vertices[1]], points[opp_vert])) {
                // perform edge flip
                TriangleAdj &tri_R = triangles[R];
                int C = tri_L.adjacencies[2];
                int A = tri_R.adjacencies[(opp_vert_id + 2) % 3];

                // update adjacencies of neighbors
                if (A >= 0) {
                    for (int m = 0; m < 3; ++m) {
                        if (triangles[A].adjacencies[m] == R) {
                            triangles[A].adjacencies[m] = L;
                            break;
                        }
                    }
                }
                if (C >= 0) {
                    for (int m = 0; m < 3; ++m) {
                        if (triangles[C].adjacencies[m] == L) {
                            triangles[C].adjacencies[m] = R;
                            break;
                        }
                    }
                }

                // update triangle r
                tri_R.vertices[(opp_vert_id + 2) % 3] = i;
                for (int m = 0; m < 3; ++m) {
                    if (tri_R.adjacencies[m] == L) {
                        tri_R.adjacencies[m] = C;
                        break;
                    }
                }
                for (int m = 0; m < 3; ++m) {
                    if (tri_R.adjacencies[m] == A) {
                        tri_R.adjacencies[m] = L;
                        break;
                    }
                }
                while (tri_R.vertices[0] != i) {
                    std::rotate(tri_R.vertices.begin(), tri_R.vertices.begin() + 1, tri_R.vertices.end());
                    std::rotate(tri_R.adjacencies.begin(), tri_R.adjacencies.begin() + 1, tri_R.adjacencies.end());
                }

                // update triangle l
                tri_L.vertices[2] = opp_vert;
                for (int m = 0; m < 3; ++m) {
                    if (tri_L.adjacencies[m] == C) {
                        tri_L.adjacencies[m] = R;
                        break;
                    }
                }
                for (int m = 0; m < 3; ++m) {
                    if (tri_L.adjacencies[m] == R) {
                        tri_L.adjacencies[m] = A;
                        break;
                    }
                }

                // add to stack
                if (tri_L.adjacencies[1] >= 0)
                    stack.push_back(L);
                if (tri_R.adjacencies[1] >= 0)
                    stack.push_back(R);
            }
        }
    }

    // remove triangles containing super-triangle vertices
    std::vector<TriangleAdj> final_triangles;
    for (const auto &tri : triangles) {
        bool uses_super_vertex = false;
        for (int v : tri.vertices) {
            if (v >= num_original_points) {
                uses_super_vertex = true;
                break;
            }
        }
        if (!uses_super_vertex) {
            final_triangles.push_back(tri);
        }
    }

    // denormalize points
    points.resize(num_original_points);
    for (auto &p : points) {
        p.x = p.x * scale + offset.x;
        p.y = p.y * scale + offset.y;
    }

    return DelaunayResult{points, final_triangles};
}