/*
 * decimate.c -- quadric-error edge-collapse simplification for ez3d.
 *
 * Classic Garland & Heckbert QEM in the practical shape popularized by
 * Forstmann's fast-quadric: every vertex carries a symmetric 4x4 quadric
 * summed from its (area-weighted) incident face planes; an edge collapse
 * costs v^T (Q_u + Q_v) v evaluated at the best contraction point; a
 * min-heap always collapses the globally cheapest edge until the cheapest
 * remaining cost exceeds the budget. Heap entries are validated lazily
 * (both endpoints alive, edge still present) and their error is recomputed
 * on pop, because earlier collapses grow the quadrics -- and since costs
 * only ever grow, an entry that re-pops above budget is discarded for
 * good.
 *
 * Manifoldness is preserved by refusing collapses that violate the edge
 * link condition: collapsing (u,v) is legal only when the neighbor sets of
 * u and v intersect exactly in the apex vertices of the triangles sharing
 * the edge. Additionally an interior edge between two boundary (rim)
 * vertices is never collapsed, which keeps the print outline a single
 * closed loop.
 *
 * The error metric is area-weighted (standard: flat regions collapse
 * first), so the budget in mm is a nominal scale for the level ladder in
 * main.c, not a hard geometric tolerance.
 */
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "decimate.h"
#include "triangle.h"

#define DC_NONE UINT32_MAX

static void *dc_alloc(size_t n, size_t sz)
{
    void *p = calloc(n, sz);
    if (!p) {
        fprintf(stderr, "ez3d: decimate: out of memory\n");
        exit(2);
    }
    return p;
}

/* ---- symmetric 4x4 quadric ------------------------------------------ */

typedef struct {
    double xx, xy, xz, xw;
    double yy, yz, yw;
    double zz, zw;
    double ww;
} quadsym_t;

/*
 * Add one 4x4 matrix to another.
 */
static void quad_add(quadsym_t *d, const quadsym_t *s)
{
    d->xx += s->xx; d->xy += s->xy; d->xz += s->xz; d->xw += s->xw;
    d->yy += s->yy; d->yz += s->yz; d->yw += s->yw;
    d->zz += s->zz; d->zw += s->zw;
    d->ww += s->ww;
}

/*
 * What does this do?
 */
static double quad_eval(const quadsym_t *q, const vertex_t *v)
{
    const double x = v->x;
    const double y = v->y;
    const double z = v->z;

    return q->xx * x * x + 2 * q->xy * x * y + 2 * q->xz * x * z
           + 2 * q->xw * x
           + q->yy * y * y + 2 * q->yz * y * z + 2 * q->yw * y
           + q->zz * z * z + 2 * q->zw * z
           + q->ww;
}

/* Contraction point of minimum quadric error: solve the 3x3 normal system
 * of the summed quadric; fall back if it is singular. */
static int quad_optimum(const quadsym_t *q, vertex_t *out)
{
    double a11 = q->xx, a12 = q->xy, a13 = q->xz;
    double a22 = q->yy, a23 = q->yz, a33 = q->zz;
    double b1 = -q->xw, b2 = -q->yw, b3 = -q->zw;
    double det, px, py, pz;

    det = a11 * (a22 * a33 - a23 * a23)
          - a12 * (a12 * a33 - a23 * a13)
          + a13 * (a12 * a23 - a22 * a13);
    if (fabs(det) < 1e-12 * (fabs(a11 * a22 * a33) + 1.0))
        return 0;

    /*
     * Do the calculation as double, for best accuracy.
     */
    px = (b1 * (a22 * a33 - a23 * a23) - a12 * (b2 * a33 - a23 * b3)
          + a13 * (b2 * a23 - a22 * b3)) / det;
    py = (a11 * (b2 * a33 - a23 * b3) - b1 * (a12 * a33 - a13 * a23)
          + a13 * (a12 * b3 - b2 * a13)) / det;
    pz = (a11 * (a22 * b3 - b2 * a23) - a12 * (a12 * b3 - b2 * a13)
          + b1 * (a12 * a23 - a22 * a13)) / det;
    if (!isfinite(px) || !isfinite(py) || !isfinite(pz))
        return 0;

    /*
     * Final vertex must be float because STL can only output this.
     */
    out->x = (float)px;
    out->y = (float)py;
    out->z = (float)pz;
    return 1;
}

/*
 * Best contraction position and its quadric cost for edge (u,v). The
 * solved optimum is only trusted when it stays inside the endpoints'
 * bounding box: a near-singular summed quadric has a free direction in
 * which the "optimal" point runs arbitrarily far away at zero cost.
 * `on_segment` (crease edges) restricts candidates to points on the edge
 * itself, keeping protected vertices on their feature line.
 */
static double edge_cost(const quadsym_t *qu, const quadsym_t *qv,
                        const vertex_t *pu, const vertex_t *pv,
                        int on_segment, vertex_t *out)
{
    quadsym_t q = *qu;
    vertex_t mid;
    vertex_t p;
    double best = HUGE_VAL;

    quad_add(&q, qv);
    vertex_mid(pu, pv, &mid);
    for (int c = 0; c < 3; c++) {
        const vertex_t *pt = (c == 0) ? pu : (c == 1) ? pv : &mid;
        double e = quad_eval(&q, pt);

        if (e < best) {
            best = e;
            *out = *pt;
        }
    }
    if (!on_segment && quad_optimum(&q, &p)) {
        if (vertex_cube_is_inside(pu, pv, &p)) {
            double e = quad_eval(&q, &p);

            if (e < best) {
                best = e;
                *out = p;
            }
        }
    }
    return best;
}

/* ---- min-heap of candidate collapses --------------------------------- */

typedef struct {
    double err;
    uint32_t u, v;
} hent;

typedef struct {
    hent *e;
    size_t n, cap;
} heap;

static void heap_push(heap *h, double err, uint32_t u, uint32_t v)
{
    if (h->n == h->cap) {
        h->cap = h->cap ? h->cap * 2 : 1 << 16;
        h->e = realloc(h->e, h->cap * sizeof *h->e);
        if (!h->e) {
            fprintf(stderr, "ez3d: decimate: out of memory\n");
            exit(2);
        }
    }
    size_t i = h->n++;

    h->e[i].err = err;
    h->e[i].u = u;
    h->e[i].v = v;
    while (i > 0) {
        size_t up = (i - 1) / 2;

        if (h->e[up].err <= h->e[i].err)
            break;
        hent tmp = h->e[up];
        h->e[up] = h->e[i];
        h->e[i] = tmp;
        i = up;
    }
}

static hent heap_pop(heap *h)
{
    hent top = h->e[0];
    size_t i = 0;

    h->e[0] = h->e[--h->n];
    for (;;) {
        size_t l = 2 * i + 1, r = l + 1, s = i;

        if (l < h->n && h->e[l].err < h->e[s].err)
            s = l;
        if (r < h->n && h->e[r].err < h->e[s].err)
            s = r;
        if (s == i)
            break;
        hent tmp = h->e[s];
        h->e[s] = h->e[i];
        h->e[i] = tmp;
        i = s;
    }
    return top;
}

/* ---- per-vertex incident triangles ----------------------------------- */

typedef struct {
    uint32_t *t;
    uint32_t n, cap;
} vtrilist;

static void vt_push(vtrilist *l, uint32_t tri)
{
    if (l->n == l->cap) {
        l->cap = l->cap ? l->cap * 2 : 8;
        l->t = realloc(l->t, l->cap * sizeof *l->t);
        if (!l->t) {
            fprintf(stderr, "ez3d: decimate: out of memory\n");
            exit(2);
        }
    }
    l->t[l->n++] = tri;
}

/*
 * Does triangle t contain vertex x?
 * return 0, 1, 2 as the vertex index within the triangle if found.
 * return -1 if it fails, and not found.
 */
static int tri_has(const trimesh_t *m, uint32_t t, uint32_t x)
{
    for (int k = 0; k < 3; k++)
        if (m->tri[t].v[k] == x)
            return k;
    return -1;
}

/* ---- live-position index ---------------------------------------------
 * Chains of vertex ids over exact float coordinates. Two live vertices
 * must never share a position: a slicer (and our verifier) merge them
 * into one point, gluing their triangle fans into a non-manifold. The
 * decimator keeps this index current as vertices move. */

typedef struct {
    uint32_t *bucket;
    uint32_t *next;    // per-vertex chain link
    size_t nb;         // power of two
} poshash_t;

static uint64_t pos_hash(vertex_t *v)
{
    float key[3] = { v->x, v->y, v->z };
    const unsigned char *b = (const void *)key;
    uint64_t h = 1469598103934665603ull;

    for (size_t i = 0; i < 12; i++) {
        h ^= b[i];
        h *= 1099511628211ull;
    }
    return h;
}

static void pos_init(poshash_t *ph, size_t nv)
{
    ph->nb = 1 << 12;
    while (ph->nb < nv * 2)
        ph->nb *= 2;
    ph->bucket = dc_alloc(ph->nb, 4);
    for (size_t i = 0; i < ph->nb; i++)
        ph->bucket[i] = DC_NONE;
    ph->next = dc_alloc(nv, 4);
}

static uint32_t pos_lookup(const poshash_t *ph, const trimesh_t *m, vertex_t *v)
{
    size_t j = pos_hash(v) & (ph->nb - 1);

    for (uint32_t i = ph->bucket[j]; i != DC_NONE; i = ph->next[i]) {
        if (m->v[i].x == v->x && m->v[i].y == v->y && m->v[i].z == v->z)
            return i;
    }
    return DC_NONE;
}

static void pos_remove(poshash_t *ph, const trimesh_t *m, uint32_t vi)
{
    vertex_t *v = &m->v[vi];
    size_t j = pos_hash(v) & (ph->nb - 1);
    uint32_t *p = &ph->bucket[j];

    while (*p != DC_NONE) {
        if (*p == vi) {
            *p = ph->next[vi];
            return;
        }
        p = &ph->next[*p];
    }
}

static void pos_insert(poshash_t *ph, const trimesh_t *m, uint32_t vi)
{
    size_t j = pos_hash(&m->v[vi]) & (ph->nb - 1);

    ph->next[vi] = ph->bucket[j];
    ph->bucket[j] = vi;
}

/* ---- edge census (initial; keys are stable ids) ---------------------- */

typedef struct {
    uint32_t *a, *b, *cnt;     // a < b; cnt == DC_NONE means empty slot
    uint32_t *t1, *t2;         // adjacent triangles (t2 unused if cnt<2)
    size_t cap, count;
} edgeset_t;

static uint64_t edge_hash(uint32_t a, uint32_t b)
{
    uint64_t h = 1469598103934665603ull;

    h ^= a;
    h *= 1099511628211ull;
    h ^= b;
    h *= 1099511628211ull;
    return h;
}

static void es_init(edgeset_t *es, size_t expected)
{
    size_t cap = 1 << 12;

    while (cap < expected * 2)
        cap *= 2;
    es->cap = cap;
    es->count = 0;
    es->a = dc_alloc(cap, 4);
    es->b = dc_alloc(cap, 4);
    es->cnt = dc_alloc(cap, 4);
    es->t1 = dc_alloc(cap, 4);
    es->t2 = dc_alloc(cap, 4);
    memset(es->cnt, 0xFF, cap * 4);
}

static void es_grow(edgeset_t *es)
{
    edgeset_t old = *es;

    es->cap *= 2;
    es->a = dc_alloc(es->cap, 4);
    es->b = dc_alloc(es->cap, 4);
    es->cnt = dc_alloc(es->cap, 4);
    es->t1 = dc_alloc(es->cap, 4);
    es->t2 = dc_alloc(es->cap, 4);
    memset(es->cnt, 0xFF, es->cap * 4);
    for (size_t i = 0; i < old.cap; i++)
        if (old.cnt[i] != DC_NONE) {
            size_t j = edge_hash(old.a[i], old.b[i]) & (es->cap - 1);

            while (es->cnt[j] != DC_NONE)
                j = (j + 1) & (es->cap - 1);
            es->a[j] = old.a[i];
            es->b[j] = old.b[i];
            es->cnt[j] = old.cnt[i];
            es->t1[j] = old.t1[i];
            es->t2[j] = old.t2[i];
        }
    free(old.a);
    free(old.b);
    free(old.cnt);
    free(old.t1);
    free(old.t2);
}

/*
 * Count one more triangle on undirected edge (u,v); a<b on insertion.
 */
static void es_add(edgeset_t *es, uint32_t u, uint32_t v, uint32_t tri)
{
    if (u > v) {
        uint32_t tmp = u;
        u = v;
        v = tmp;
    }
    if ((es->count + 1) * 10 >= es->cap * 7)
        es_grow(es);
    size_t j = edge_hash(u, v) & (es->cap - 1);

    while (es->cnt[j] != DC_NONE) {
        if (es->a[j] == u && es->b[j] == v) {
            if (es->cnt[j] == 1)
                es->t2[j] = tri;
            es->cnt[j]++;
            return;
        }
        j = (j + 1) & (es->cap - 1);
    }
    es->a[j] = u;
    es->b[j] = v;
    es->t1[j] = tri;
    es->cnt[j] = 1;
    es->count++;
}

/* Unit normal of triangle t from its current vertex positions. */
static void face_normal(const trimesh_t *m, uint32_t t, double out[3])
{
    const triangle_t *tp = &m->tri[t];
    const vertex_t *a = &m->v[tp->v[0]];
    const vertex_t *b = &m->v[tp->v[1]];
    const vertex_t *c = &m->v[tp->v[2]];
    dvec_t norm;
    tri_normal_vec(a, b, c, &norm);
    double len = tri_normal_len(&norm);

    if (len > 0.0) {
        out[0] = norm.x / len;
        out[1] = norm.y / len;
        out[2] = norm.z / len;
    } else {
        out[0] = out[1] = out[2] = 0.0;
    }
}

/* Edges whose adjacent faces meet sharper than this are creases (feature
 * lines): the cap-to-wall edge is 90 degrees, terrain-to-wall typically
 * 60-90, smooth terrain interior well under 20. */
#define CREASE_DOT 0.766 /* cos(40 degrees) */

/* Crease (rim) edges longer than this are never collapsed. Keeping a rim
 * vertex at least every ~2 units bounds how far a triangle anchored on
 * the rim can chord across a concave corner such as a narrow cut-out;
 * without the cap, an entire rim run can merge down to its two end
 * corners and the attached terrain/cap triangles then bridge the mouth. */
#define MAX_CREASE_COLLAPSE 2.0

/* Upper limit on pending heap entries (~1 GB at 16 bytes each). The
 * neighbor dedup below keeps the heap proportional to real candidates,
 * so this is a safety net, not a tuning knob: when it trips, further
 * candidates are dropped and the mesh keeps whatever detail it has. */
#define MAX_HEAP_ENTRIES (64ull << 20)

/* Refuse collapses into a vertex whose incident-triangle fan is larger
 * than this. Every pop scans the endpoint fans (link condition, fold
 * preview, re-push), so unbounded fans make the collapse loop quadratic
 * in time on big flat caps; the cap bounds per-collapse work. Flat bulk
 * is much better handled by --simplify before decimation anyway. */
#define MAX_FAN_TRIS 128

size_t ez3d_decimate(trimesh_t *m, double err_mm)
{
    const size_t nv0 = m->nv, ntri0 = m->ntri;
    quadsym_t *Q;
    dvec_t *fn;           /* ntri0 unit face normals */
    vtrilist *vt;
    uint8_t *vdead, *tdead, *crease;
    uint32_t *stamp;
    uint32_t *tmark;        /* per-collapse dedup of vt[] duplicates */
    poshash_t ph;
    uint32_t gen = 0;
    heap h = { NULL, 0, 0 };
    edgeset_t es;
    size_t ntri_alive = 0;
    double budget = err_mm * err_mm;
    int nonmanifold = 0;
    int warned_heap_full = 0;

    if (nv0 == 0 || ntri0 == 0)
        return ntri0;

    Q = dc_alloc(nv0, sizeof(quadsym_t));
    vt = dc_alloc(nv0, sizeof(vtrilist));
    vdead = dc_alloc(nv0, 1);
    tdead = dc_alloc(ntri0, 1);
    crease = dc_alloc(nv0, 1);
    stamp = dc_alloc(nv0, sizeof *stamp);
    tmark = dc_alloc(ntri0, sizeof *tmark);
    fn = dc_alloc(ntri0, sizeof *fn);
    pos_init(&ph, nv0);

    /* Face planes -> vertex quadrics + adjacency. Zero-area input
     * triangles (there should be none) are dropped up front. */
    for (uint32_t t = 0; t < ntri0; t++) {

        const triangle_t *tp = &m->tri[t];
        const vertex_t *a = &m->v[tp->v[0]];
        const vertex_t *b = &m->v[tp->v[1]];
        const vertex_t *c = &m->v[tp->v[2]];
        dvec_t norm;
        tri_normal_vec(a, b, c, &norm);

        double len = tri_normal_len(&norm);
        double px, py, pz, d;
        quadsym_t q;

        if (len <= 0.0) {
            tdead[t] = 1;
            continue;
        }
        norm.x /= len;
        norm.y /= len;
        norm.z /= len;
        fn[t].x = norm.x;
        fn[t].y = norm.y;
        fn[t].z = norm.z;
        d = -(norm.x * a->x + norm.y * a->y + norm.z * a->z); // Dot-product
        /*
         * area-weighted plane coefficients
         */
        px = 0.5 * len * norm.x;
        py = 0.5 * len * norm.y;
        pz = 0.5 * len * norm.z;
        q.xx = px * norm.x; q.xy = px * norm.y; q.xz = px * norm.z; q.xw = px * d;
        q.yy = py * norm.y; q.yz = py * norm.z; q.yw = py * d;
        q.zz = pz * norm.z; q.zw = pz * d;
        q.ww = 0.5 * len * d * d;
        for (int kk = 0; kk < 3; kk++) {
            quad_add(&Q[m->tri[t].v[kk]], &q);
            vt_push(&vt[m->tri[t].v[kk]], t);
        }
        ntri_alive++;
    }

    /*
     * Index every live position.
     * The mesher deduped vertices exactly, so initial positions are unique.
     */
    for (uint32_t vi = 0; vi < nv0; vi++)
        pos_insert(&ph, m, vi);

    /*
     * Edge census: count triangles per undirected edge and mark crease
     * vertices -- endpoints of edges whose faces meet at a sharp dihedral
     * angle. In a closed heightmap solid the creases ARE the print rim,
     * hole outlines and wall tops; protecting them keeps cut-out regions
     * from being bridged or slumped by aggressive collapses.
     */
    es_init(&es, ntri0 * 3);
    for (uint32_t t = 0; t < ntri0; t++) {
        if (tdead[t])
            continue;
        for (int kk = 0; kk < 3; kk++) {
            uint32_t u = m->tri[t].v[kk];
            uint32_t v = m->tri[t].v[(kk + 1) % 3];

            if (u != v)
                es_add(&es, u, v, t);
        }
    }
    for (size_t j = 0; j < es.cap; j++) {
        if (es.cnt[j] == DC_NONE)
            continue;
        if (es.cnt[j] > 2) {
            nonmanifold = 1;
        } else if (es.cnt[j] == 2) {
            const dvec_t *n1 = &fn[es.t1[j]];
            const dvec_t *n2 = &fn[es.t2[j]];
            double dot = tri_dot_prod(n1, n2);

            if (dot < CREASE_DOT) {
                crease[es.a[j]] = 1;
                crease[es.b[j]] = 1;
            }
        }
    }
    if (nonmanifold) {
        /* Not produced by our mesher; refuse rather than corrupt. */
        fprintf(stderr, "ez3d: decimate: mesh is not edge-manifold, "
                        "leaving it unchanged\n");
        goto out;
    }

    /* Seed the heap with every edge already inside the budget. Costs only
     * grow as collapses accumulate, so anything above budget now stays
     * above budget forever. */
    for (size_t j = 0; j < es.cap; j++) {
        if (es.cnt[j] == DC_NONE)
            continue;
        uint32_t u = es.a[j], v = es.b[j];
        vertex_t unused;
        double e = edge_cost(&Q[u], &Q[v], &m->v[u], &m->v[v], 0, &unused);

        if (e <= budget)
            heap_push(&h, e, u, v);
    }

    /* Collapse loop. */
    while (h.n > 0) {
        hent e = heap_pop(&h);
        uint32_t u = e.u, v = e.v;
        vertex_t candidate;
        double cost;
        uint32_t apex0 = DC_NONE, apex1 = DC_NONE;
        int nshared = 0, legal = 1;
        uint32_t w;

        if (vdead[u] || vdead[v])
            continue;
        if (vt[u].n > MAX_FAN_TRIS || vt[v].n > MAX_FAN_TRIS)
            continue; /* mega-fan: leave it, see MAX_FAN_TRIS */
        cost = edge_cost(&Q[u], &Q[v], &m->v[u], &m->v[v], 0, &candidate);
        if (cost > budget)
            continue; // permanently: costs never decrease

        /*
         * Which triangles currently share the edge (u,v)? Compacting the
         * dead entries out of vt[u] as we go keeps scans short; tmark
         * dedups triangles that entered vt[u] twice via earlier collapses.
         */
        uint32_t g_shared = ++gen;
        uint32_t sh_t0 = 0, sh_t1 = 0;

        w = 0;
        for (uint32_t i = 0; i < vt[u].n; i++) {
            uint32_t t = vt[u].t[i];

            if (tdead[t] || tmark[t] == g_shared)
                continue;
            tmark[t] = g_shared;
            vt[u].t[w++] = t;
            int ku = tri_has(m, t, u);
            int kv = tri_has(m, t, v);
            if (ku >= 0 && kv >= 0) { // Yes, it has the edge (u, v)
                uint32_t apex = m->tri[t].v[3 - ku - kv];

                if (nshared == 0) {
                    apex0 = apex;
                    sh_t0 = t;
                } else {
                    apex1 = apex;
                    sh_t1 = t;
                }
                nshared++;
            }
        }
        vt[u].n = w;
        if (nshared == 0)
            continue; // edge vanished in an earlier collapse
        if (nshared > 2)
            continue; // degenerate configuration; leave it alone

        /* Is this still a crease edge? Recompute from current geometry,
         * since collapses may have flattened or sharpened the faces. */
        int is_crease = 0;

        if (nshared == 2) {
            double n1[3], n2[3];

            face_normal(m, sh_t0, n1);
            face_normal(m, sh_t1, n2);
            if (n1[0] * n2[0] + n1[1] * n2[1] + n1[2] * n2[2] < CREASE_DOT)
                is_crease = 1;
        }
        /* Crease vertices may only merge along a crease edge, and the
         * merged point stays on the crease line. This keeps rims, wall
         * tops and hole outlines intact AND freezes the first ring of
         * terrain triangles around each rim: without that, triangles
         * anchored at a rim corner can chord across a narrow cut-out
         * (quadrics are blind to empty space -- holes contribute no
         * planes to the error metric). Long rim runs are also kept dense
         * enough that no attached triangle can bridge a concave mouth. */
        if ((crease[u] || crease[v]) && !is_crease)
            continue;
        if (!is_crease) {
            /* First-ring freeze: an interior vertex adjacent to any
             * crease vertex never merges. Otherwise interior chains can
             * walk all the way around a narrow cut-out and drag a
             * rim-anchored triangle's far vertex across the mouth (the
             * flat bottom cap is especially prone to this: its interior
             * collapses are quadric-error free). */
            int near = 0;

            for (uint32_t i = 0; i < vt[u].n && !near; i++) {
                uint32_t t = vt[u].t[i];

                if (tdead[t])
                    continue;
                for (int kk = 0; kk < 3; kk++) {
                    uint32_t nb = m->tri[t].v[kk];

                    if (nb != u && nb != v && crease[nb]) {
                        near = 1;
                        break;
                    }
                }
            }
            for (uint32_t i = 0; i < vt[v].n && !near; i++) {
                uint32_t t = vt[v].t[i];

                if (tdead[t])
                    continue;
                for (int kk = 0; kk < 3; kk++) {
                    uint32_t nb = m->tri[t].v[kk];

                    if (nb != u && nb != v && crease[nb]) {
                        near = 1;
                        break;
                    }
                }
            }
            if (near)
                continue;
        }
        if (is_crease) {
            dvec_t delta;
            vertex_diff(&m->v[u], &m->v[v], &delta);
            if (delta.x * delta.x + delta.y * delta.y + delta.z * delta.z
                > MAX_CREASE_COLLAPSE * MAX_CREASE_COLLAPSE)
                continue;
        }
        cost = edge_cost(&Q[u], &Q[v], &m->v[u], &m->v[v],
                         is_crease, &candidate);
        if (cost > budget)
            continue; /* permanently: costs never decrease */

        /* Link condition: N(u) ^ N(v) must be exactly the apex set. */
        gen++;
        for (uint32_t i = 0; i < vt[u].n; i++) {
            uint32_t t = vt[u].t[i];

            for (int kk = 0; kk < 3; kk++) {
                uint32_t nb = m->tri[t].v[kk];

                if (nb != u && nb != v)
                    stamp[nb] = gen;
            }
        }
        for (uint32_t i = 0; i < vt[v].n && legal; i++) {
            uint32_t t = vt[v].t[i];

            if (tdead[t])
                continue;
            for (int kk = 0; kk < 3 && legal; kk++) {
                uint32_t nb = m->tri[t].v[kk];

                if (nb == u || nb == v)
                    continue;
                if (stamp[nb] == gen && nb != apex0 && nb != apex1)
                    legal = 0; // shared neighbor beyond the apexes
            }
        }
        if (!legal)
            continue;

        /* Degenerate/fold preview: every triangle that survives the
         * collapse (rewritten with v->u at p) must keep a sane area and
         * stay within ~87 degrees of its ORIGINAL orientation. The
         * absolute bound matters: comparing against the pre-collapse
         * orientation only would let a chain of individually-legal
         * collapses walk a terrain triangle past vertical -- a fold.
         * Folded facets face into the model and render as dark triangles
         * in shader-based viewers. */
        for (uint32_t pass = 0; pass < 2 && legal; pass++) {
            vtrilist *l = pass ? &vt[v] : &vt[u];
            uint32_t other = pass ? u : v;

            for (uint32_t i = 0; i < l->n && legal; i++) {
                uint32_t t = l->t[i];

                if (tdead[t])
                    continue;
                if (tri_has(m, t, other) >= 0)
                    continue; /* spans (u,v): dies with the edge */

                vertex_t *pv[3];
                for (int kk = 0; kk < 3; kk++) {
                    uint32_t x = m->tri[t].v[kk];

                    pv[kk] = (x == u || x == v) ? &candidate : &m->v[x];
                }

                dvec_t nn;
                tri_normal_vec(pv[0], pv[1], pv[2], &nn);

                double fnz = fn[t].z, len = tri_normal_len(&nn);
                if (len < 1e-6 || (nn.x * fn[t].x + nn.y * fn[t].y
                                   + nn.z * fn[t].z) / len < 0.05)
                    legal = 0;
                /* Heightfield surfaces never change facing: terrain
                 * stays upward-facing, the cap stays downward-facing
                 * (walls start exactly vertical, fn.z == 0, and only get
                 * the 87-degree rule above). */
                if (legal && fnz > 0.0 && nn.z < 0.0)
                    legal = 0;
                if (legal && fnz < 0.0 && nn.z > 0.0)
                    legal = 0;
            }
        }
        if (!legal)
            continue;

        /* Never land on another live vertex's exact position: a slicer
         * (and our verifier) merge coincident vertices, which would glue
         * two triangle fans into a non-manifold point. Landing on u or v
         * themselves is fine (that is an endpoint collapse). */
        {
            uint32_t hit = pos_lookup(&ph, m, &candidate);

            if (hit != DC_NONE && hit != u && hit != v)
                continue;
        }

        /* Collapse v into u at the optimal point. */
        uint32_t g_rw = ++gen;

        pos_remove(&ph, m, u);
        pos_remove(&ph, m, v);
        m->v[u] = candidate;
        pos_insert(&ph, m, u);
        quad_add(&Q[u], &Q[v]);
        crease[u] |= crease[v];
        vdead[v] = 1;
        for (uint32_t i = 0; i < vt[v].n; i++) {
            uint32_t t = vt[v].t[i];

            if (tdead[t] || tmark[t] == g_rw)
                continue;
            tmark[t] = g_rw;
            if (tri_has(m, t, u) >= 0) {
                tdead[t] = 1; /* the triangle spanning (u,v) */
                ntri_alive--;
                continue;
            }
            for (int kk = 0; kk < 3; kk++)
                if (m->tri[t].v[kk] == v)
                    m->tri[t].v[kk] = u;
            vt_push(&vt[u], t);
        }
        /* New candidate edges around the merged vertex. Each distinct
         * neighbour is pushed at most once: without the stamp dedup a
         * vertex whose fan has grown large re-pushes the same edge once
         * per shared triangle, and the heap grows quadratically with the
         * number of collapses -- large flat inputs then exhaust memory
         * long before the budget is reached. */
        uint32_t g_push = ++gen;
        int heap_full = 0;

        for (uint32_t i = 0; i < vt[u].n; i++) {
            uint32_t t = vt[u].t[i];

            if (tdead[t])
                continue;
            for (int kk = 0; kk < 3; kk++) {
                uint32_t nb = m->tri[t].v[kk];
                vertex_t unused;
                double ec;

                if (nb == u || vdead[nb])
                    continue;
                if (stamp[nb] == g_push)
                    continue; /* already queued in this round */
                stamp[nb] = g_push;
                if (h.n >= MAX_HEAP_ENTRIES) {
                    heap_full = 1;
                    continue;
                }
                ec = edge_cost(&Q[u], &Q[nb], &m->v[u], &m->v[nb], 0, &unused);
                if (ec <= budget)
                    heap_push(&h, ec, u, nb);
            }
        }
        if (heap_full && !warned_heap_full) {
            fprintf(stderr, "ez3d: decimate: candidate heap full (%llu "
                            "entries), keeping current detail\n",
                    (unsigned long long)MAX_HEAP_ENTRIES);
            warned_heap_full = 1;
        }
    }

out:
    /* Compact: rebuild the trimesh with only live vertices/faces. */
    if (ntri_alive > 0) {
        vertex_t *nv = dc_alloc(nv0 + 1, sizeof *nv);
        uint32_t *remap = dc_alloc(nv0, 4);
        uint32_t nvn = 0;
        trimesh_t fresh;
        size_t kept = 0;

        for (size_t vi = 0; vi < nv0; vi++)
            remap[vi] = vdead[vi] ? DC_NONE : nvn++;
        for (size_t vi = 0; vi < nv0; vi++)
            if (!vdead[vi]) {
                nv[remap[vi]] = m->v[vi];
            }
        tm_init(&fresh);
        /* Survivors keep their own ids even when two of them landed on
         * identical coordinates: deduping them would glue two triangle
         * fans into one point and create over-shared edges. No live
         * triangle can have two coincident corners -- the area preview
         * rejects those collapses -- so appending is safe. */
        for (size_t vi = 0; vi < nvn; vi++)
            tm_append(&fresh, &nv[vi]);

       for (uint32_t t = 0; t < ntri0; t++) {
            const triangle_t *tp = &m->tri[t];
            triangle_t t2;

            if (tdead[t])
                continue;
            t2.v[0] = remap[tp->v[0]];
            t2.v[1] = remap[tp->v[1]];
            t2.v[2] = remap[tp->v[2]];
            if (t2.v[0] == t2.v[1] || t2.v[1] == t2.v[2] || t2.v[2] == t2.v[0])
                continue; /* corners collapsed onto one point */
            tm_tri(&fresh, &t2);
            kept++;
        }
        tm_free(m);
        *m = fresh;
        ntri_alive = kept;
        free(nv);
        free(remap);
    }

    free(h.e);
    free(Q);
    free(fn);
    for (size_t vi = 0; vi < nv0; vi++)
        free(vt[vi].t);
    free(vt);
    free(vdead);
    free(tdead);
    free(crease);
    free(stamp);
    free(tmark);
    free(ph.bucket);
    free(ph.next);
    free(es.a);
    free(es.b);
    free(es.cnt);
    free(es.t1);
    free(es.t2);
    return ntri_alive;
}