base.cpp

#include <ultimaille/all.h>
#include <string>
#define FOR(i, n) for(int i = 0; i < static_cast<int>(n); i++)

using namespace UM;

vec<3> circumcenter(const vec<3>& A,
                    const vec<3>& B,
                    const vec<3>& C)
{ // CHATGPT
    const vec<3> u = B - A;
    const vec<3> v = C - A;
    const vec<3> n = cross(u, v);
    const double n2 = n.norm2();
    assert(n2 > 0.0);

    return A
         + ((cross(v, n) * u.norm2())
         +  (cross(n, u) * v.norm2()))
         / (2.0 * n2);
}

double triangle_area(const vec<3>& A,
                     const vec<3>& B,
                     const vec<3>& C)
{
    return cross((B - A), (C - A)).norm() / 2;
}

void compute_barycentre(Triangles& m, PointAttribute<vec4>& barycentres) {
    for (auto A : m.iter_vertices()) {
        for (auto he : A.iter_halfedges()) {
            auto f = he.facet();
            assert(f.active());
            auto B = he.to();
            auto C = he.next().to();
            auto O = circumcenter(A, B, C);
            auto E = (A.pos() + B.pos()) / 2;
            auto F = (A.pos() + C.pos()) / 2;
            auto area = triangle_area(A, O, E) + triangle_area(A, F, O);
            auto M = (A.pos() + F + O + E) / 4;
            barycentres[A] += M.xyz0() * area + vec4(0., 0., 0., area);
        }
    }
}


std::tuple<double, bool> compute_angle(Surface::Halfedge he) {
    auto ho = he.opposite();
    auto bot = he.from();
    auto top = he.to();
    auto left = he.next().to();
    auto right = ho.next().to();
    vec3 u = left.pos() - top.pos();
    vec3 v = left.pos() - bot.pos();
    double angle = atan2(cross(u, v).norm(), u*v);
    v = right.pos() - top.pos();
    double angle2 = atan2(cross(u, v).norm(), u*v);
    u = left.pos() - bot.pos();
    v = right.pos() - bot.pos();
    double angle3 = atan2(cross(u, v).norm(), u*v);
    return {angle, angle > angle2 && angle > angle3};
}

bool flip_largest_edge(Triangles& m, FacetAttribute<int>& flipped) {
    Surface::Halfedge largest = {m, -1};
    double largest_angle = 0;
    for (auto he : m.iter_halfedges()) {
        if (!he.active()) continue;
        if (!he.opposite().active()) continue;        
        auto [angle, ok] = compute_angle(he);
        if (!ok) continue;
        if (angle >= largest_angle) {
            largest = he;
            largest_angle = angle;
        }
    }
    if (largest != -1) {
        auto he = largest;
        auto opposite = he.opposite();
        auto vbas = he.from();
        auto vhaut = he.to();
        auto vgauche = he.next().to();
        auto vdroite = opposite.next().to();
        m.conn->active[he.facet()] = false;
        m.conn->active[opposite.facet()] = false;
        flipped[m.conn->create_facet({vbas, vdroite, vgauche})] = 1;
        flipped[m.conn->create_facet({vdroite, vhaut, vgauche})] = 1;
        std::cout << std::to_string(largest) << "\n";
        return true;
    }
    return false;
}


int main(int argc, char** argv) {
    UM::Triangles m;
    SurfaceAttributes attributes = read_by_extension("../after_distorted.geogram", m);

    FacetAttribute<int> flipped(m);
    PointAttribute<vec4> barycentres(m);
    m.connect();
    FOR(i, 100) {
        while (flip_largest_edge(m, flipped)) {}
        compute_barycentre(m, barycentres);
        for (auto v : m.iter_vertices()) {
            if (v.on_boundary()) continue;
            m.points[v] = barycentres[v].xyz() / barycentres[v][3];
        }
        m.compact();
        std::string s = ".geogram";
        s.insert(0, i, '1');
        write_by_extension("../out/" + s, m);
    }
    PointAttribute<double> distances(m);
    // for (auto v : m.iter_vertices()) {
    //     distances[v] = (barycentres[v].xyz()/barycentres[v][3] - v.pos()).norm();
    // }

    // m.compact();
    // write_by_extension("after_flipped.geogram", m, {{"distances", distances}});
    return 0;
}