/*
 * verify.c -- STL sanity checker for ez3d (binary STL only).
 *
 * Re-reads the file we just wrote and checks the invariants a slicer cares
 * about:
 *   - the file really holds the facet count its header claims
 *   - every edge is shared by exactly two facets (edge-manifold)
 *   - degenerate facets (repeated vertices) are counted
 *
 * The manifold rule is relaxed in two ways for topology a heightfield
 * mesher legitimately produces: an edge shared by only one facet passes
 * if some vertex lies on it within `tol` world units (T-junction), and an
 * edge shared by four facets traversed two in each direction passes as a
 * pinch, where the solid narrows to a line but stays closed and oriented.
 *
 * Vertices live in an append-only array, so vertex ids are stable for the
 * lifetime of the check; both hash indexes (exact-identity and spatial)
 * are separate bucket arrays that re-link those stable ids when they grow.
 */
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "verify.h"

#define NONE UINT32_MAX
#define FNV_BASIS 1469598103934665603ull

static uint64_t hash_bytes(const void *p, size_t n, uint64_t h)
{
    const unsigned char *b = p;
    for (size_t i = 0; i < n; i++) {
        h ^= b[i];
        h *= 1099511628211ull; /* FNV-1a */
    }
    return h;
}

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

/* ---- vertex space: append-only vertices + two bucket-chain indexes --- */

typedef struct {
    float x, y, z;
    uint32_t next_id;      /* chain in the exact-identity index */
    uint32_t next_xy;      /* chain in the spatial (integer cell) index */
} vtx;

typedef struct {
    vtx *v;
    size_t count, cap;     /* append-only: ids never change */
    uint32_t *idb;         /* exact (x,y,z) buckets */
    size_t nidb;           /* power of two */
    uint32_t *xyb;         /* integer-cell buckets */
    size_t nxyb;           /* power of two */
} vspace;

static uint64_t hash_xyz(const vtx *v)
{
    return hash_bytes(&v->x, 12, FNV_BASIS);
}

static uint64_t hash_cell(int32_t cx, int32_t cy)
{
    return hash_bytes(&(struct { int32_t a, b; }){ cx, cy }, 8, FNV_BASIS);
}

static void vspace_init(vspace *s)
{
    s->cap = 1 << 10;
    s->count = 0;
    s->v = xcalloc(s->cap, sizeof *s->v);
    s->nidb = 1 << 10;
    s->idb = xcalloc(s->nidb, sizeof *s->idb);
    for (size_t i = 0; i < s->nidb; i++)
        s->idb[i] = NONE;
    s->nxyb = 1 << 10;
    s->xyb = xcalloc(s->nxyb, sizeof *s->xyb);
    for (size_t i = 0; i < s->nxyb; i++)
        s->xyb[i] = NONE;
}

static void idb_grow(vspace *s)
{
    free(s->idb);
    s->nidb *= 2;
    s->idb = xcalloc(s->nidb, sizeof *s->idb);
    for (size_t i = 0; i < s->nidb; i++)
        s->idb[i] = NONE;
    for (size_t i = 0; i < s->count; i++) {
        size_t j = hash_xyz(&s->v[i]) & (s->nidb - 1);
        s->v[i].next_id = s->idb[j];
        s->idb[j] = (uint32_t)i;
    }
}

static void xyb_grow(vspace *s)
{
    free(s->xyb);
    s->nxyb *= 2;
    s->xyb = xcalloc(s->nxyb, sizeof *s->xyb);
    for (size_t i = 0; i < s->nxyb; i++)
        s->xyb[i] = NONE;
    for (size_t i = 0; i < s->count; i++) {
        int32_t cx = (int32_t)floorf(s->v[i].x);
        int32_t cy = (int32_t)floorf(s->v[i].y);
        size_t j = hash_cell(cx, cy) & (s->nxyb - 1);
        s->v[i].next_xy = s->xyb[j];
        s->xyb[j] = (uint32_t)i;
    }
}

static void vspace_push_spatial(vspace *s, uint32_t id)
{
    int32_t cx = (int32_t)floorf(s->v[id].x);
    int32_t cy = (int32_t)floorf(s->v[id].y);
    size_t j = hash_cell(cx, cy) & (s->nxyb - 1);

    s->v[id].next_xy = s->xyb[j];
    s->xyb[j] = id;
}

/* Return the id of the vertex equal to (x,y,z), appending it if new. */
static uint32_t vspace_get(vspace *s, float x, float y, float z)
{
    vtx probe = { x, y, z, NONE, NONE };
    size_t j = hash_xyz(&probe) & (s->nidb - 1);

    for (uint32_t i = s->idb[j]; i != NONE; i = s->v[i].next_id)
        if (s->v[i].x == x && s->v[i].y == y && s->v[i].z == z)
            return i;

    if (s->count == s->cap) {
        s->cap *= 2;
        s->v = realloc(s->v, s->cap * sizeof *s->v);
        if (!s->v) {
            fprintf(stderr, "verify: out of memory\n");
            exit(2);
        }
    }
    if ((s->count + 1) * 10 >= s->nidb * 7)
        idb_grow(s);
    if ((s->count + 1) * 10 >= s->nxyb * 7)
        xyb_grow(s);

    j = hash_xyz(&probe) & (s->nidb - 1);
    probe.next_id = s->idb[j];
    s->idb[j] = (uint32_t)s->count;
    s->v[s->count] = probe;
    vspace_push_spatial(s, (uint32_t)s->count);
    return (uint32_t)s->count++;
}

/* ---- edge table: open addressing, key = (a,b) vertex ids, a < b ------ */
/* Keys are stable vertex ids, so rehashing this table is safe. Besides  */
/* the shared count, ab/ba record how often each directed traversal of   */
/* the edge occurs: a consistently oriented (outward) mesh traverses a   */
/* shared edge once in each direction.                                   */

typedef struct {
    uint32_t a, b;
    uint32_t ab, ba;
    uint8_t used;
} edge_ent;

typedef struct {
    edge_ent *e;
    size_t cap;            /* power of two */
    size_t count;
} etab;

static void etab_init(etab *t)
{
    t->cap = 1 << 12;
    t->count = 0;
    t->e = xcalloc(t->cap, sizeof *t->e);
}

static void etab_grow(etab *t)
{
    edge_ent *old = t->e;
    size_t oldcap = t->cap;

    t->cap *= 2;
    t->e = xcalloc(t->cap, sizeof *t->e);
    for (size_t i = 0; i < oldcap; i++)
        if (old[i].used) {
            size_t j = hash_bytes(&old[i].a, 8, FNV_BASIS) & (t->cap - 1);
            while (t->e[j].used)
                j = (j + 1) & (t->cap - 1);
            t->e[j] = old[i];
        }
    free(old);
}

/* Record the triangle edge x->y. */
static void etab_add(etab *t, uint32_t x, uint32_t y)
{
    int fwd = x < y;
    uint64_t h;
    uint32_t a = fwd ? x : y;
    uint32_t b = fwd ? y : x;

    if ((t->count + 1) * 10 >= t->cap * 7)
        etab_grow(t);
    h = hash_bytes(&(struct { uint32_t a, b; }){ a, b }, 8, FNV_BASIS);
    size_t j = (size_t)h & (t->cap - 1);
    while (t->e[j].used) {
        if (t->e[j].a == a && t->e[j].b == b) {
            if (fwd)
                t->e[j].ab++;
            else
                t->e[j].ba++;
            return;
        }
        j = (j + 1) & (t->cap - 1);
    }
    t->e[j].used = 1;
    t->e[j].a = a;
    t->e[j].b = b;
    t->e[j].ab = fwd ? 1u : 0u;
    t->e[j].ba = fwd ? 0u : 1u;
    t->count++;
}

/* Distance from point p to segment ab. */
static double dist_point_seg(const float p[3], const float a[3],
                             const float b[3])
{
    double ux = b[0] - a[0], uy = b[1] - a[1], uz = b[2] - a[2];
    double vx = p[0] - a[0], vy = p[1] - a[1], vz = p[2] - a[2];
    double uu = ux * ux + uy * uy + uz * uz;
    double t = uu > 0.0 ? (vx * ux + vy * uy + vz * uz) / uu : 0.0;

    if (t < 0.0)
        t = 0.0;
    if (t > 1.0)
        t = 1.0;
    double dx = vx - t * ux, dy = vy - t * uy, dz = vz - t * uz;
    return sqrt(dx * dx + dy * dy + dz * dz);
}

/* Does some vertex lie on edge ab within tol (T-junction closure)? */
static int edge_closed_by_vertex(const vspace *s, uint32_t ai, uint32_t bi,
                                 double tol)
{
    const vtx *a = &s->v[ai], *b = &s->v[bi];
    int32_t x0 = (int32_t)floor(fmin(a->x, b->x));
    int32_t x1 = (int32_t)floor(fmax(a->x, b->x));
    int32_t y0 = (int32_t)floor(fmin(a->y, b->y));
    int32_t y1 = (int32_t)floor(fmax(a->y, b->y));

    for (int32_t cy = y0; cy <= y1; cy++)
        for (int32_t cx = x0; cx <= x1; cx++) {
            size_t j = hash_cell(cx, cy) & (s->nxyb - 1);
            for (uint32_t vi = s->xyb[j]; vi != NONE;
                 vi = s->v[vi].next_xy) {
                if (vi == ai || vi == bi)
                    continue;
                int32_t vx = (int32_t)floorf(s->v[vi].x);
                int32_t vy = (int32_t)floorf(s->v[vi].y);
                if (vx != cx || vy != cy)
                    continue; /* hash collision from another cell */
                float p[3] = { s->v[vi].x, s->v[vi].y, s->v[vi].z };
                float pa[3] = { a->x, a->y, a->z };
                float pb[3] = { b->x, b->y, b->z };
                if (dist_point_seg(p, pa, pb) <= tol)
                    return 1;
            }
        }
    return 0;
}

int ez3d_verify_stl(const char *path, double tol)
{
    FILE *f = fopen(path, "rb");
    if (!f) {
        fprintf(stderr, "verify: cannot open '%s'\n", path);
        return 1;
    }

    unsigned char hdr[84];
    if (fread(hdr, 1, 84, f) != 84) {
        fprintf(stderr, "verify: '%s' too small for a binary STL\n", path);
        fclose(f);
        return 1;
    }
    if (memcmp(hdr, "solid", 5) == 0) {
        fprintf(stderr, "verify: '%s' looks like ASCII STL; verification "
                        "supports binary output only\n",
                path);
        fclose(f);
        return 1;
    }
    uint32_t nfacets;
    memcpy(&nfacets, hdr + 80, 4);

    if (fseek(f, 0, SEEK_END) != 0) {
        fclose(f);
        return 1;
    }
    long fsize = ftell(f);
    if (fsize != 84L + 50L * nfacets) {
        fprintf(stderr, "verify: file size %ld does not match header count "
                        "%u (expected %ld)\n",
                fsize, nfacets, 84L + 50L * nfacets);
        fclose(f);
        return 1;
    }
    if (fseek(f, 84, SEEK_SET) != 0) {
        fclose(f);
        return 1;
    }

    vspace vs;
    etab et;
    vspace_init(&vs);
    etab_init(&et);

    float bbmin[3] = { INFINITY, INFINITY, INFINITY };
    float bbmax[3] = { -INFINITY, -INFINITY, -INFINITY };
    long degenerate = 0;
    double volume = 0.0;

    for (uint32_t i = 0; i < nfacets; i++) {
        unsigned char rec[50];
        if (fread(rec, 1, 50, f) != 50) {
            fprintf(stderr, "verify: truncated at facet %u\n", i);
            fclose(f);
            return 1;
        }
        float v[9];
        memcpy(v, rec + 12, 36); /* skip the stored normal */

        uint32_t id[3];
        for (int k = 0; k < 3; k++) {
            id[k] = vspace_get(&vs, v[3 * k], v[3 * k + 1], v[3 * k + 2]);
            for (int c = 0; c < 3; c++) {
                float val = v[3 * k + c];
                if (val < bbmin[c])
                    bbmin[c] = val;
                if (val > bbmax[c])
                    bbmax[c] = val;
            }
        }
        if (id[0] == id[1] || id[1] == id[2] || id[2] == id[0])
            degenerate++;
        else {
            /* zero-area facets: three distinct but collinear vertices */
            double ux = v[3] - v[0], uy = v[4] - v[1], uz = v[5] - v[2];
            double wx = v[6] - v[0], wy = v[7] - v[1], wz = v[8] - v[2];
            double nx = uy * wz - uz * wy;
            double ny = uz * wx - ux * wz;
            double nz2 = ux * wy - uy * wx;

            if (sqrt(nx * nx + ny * ny + nz2 * nz2) < 1e-9)
                degenerate++;
        }
        /* signed-volume contribution (positive when wound outward) */
        volume += (v[0] * (v[4] * v[8] - v[5] * v[7])
                   + v[1] * (v[5] * v[6] - v[3] * v[8])
                   + v[2] * (v[3] * v[7] - v[4] * v[6])) / 6.0;

        etab_add(&et, id[0], id[1]);
        etab_add(&et, id[1], id[2]);
        etab_add(&et, id[2], id[0]);
    }
    fclose(f);

    long exact2 = 0, over2 = 0, singles = 0, tclosed = 0, badorient = 0;
    long pinch = 0;
    for (size_t j = 0; j < et.cap; j++) {
        if (!et.e[j].used)
            continue;
        uint32_t total = et.e[j].ab + et.e[j].ba;

        if (total == 2) {
            exact2++;
            if (et.e[j].ab != 1 || et.e[j].ba != 1)
                badorient++; /* both traversals same way: flipped faces */
        } else if (total == 4 && et.e[j].ab == 2 && et.e[j].ba == 2) {
            /*
             * Pinch edge: the solid narrows to a line (e.g. a cell where
             * only two opposite quarter triangles survive the cut). Four
             * facets meet along the edge, two traversing each way, so the
             * surface is still closed and consistently oriented -- it is
             * just not edge-manifold here. Like a T-junction, this is
             * benign for slicers, so count it but do not fail.
             */
            pinch++;
        } else if (total > 2) {
            over2++;
        } else {
            singles++;
            if (edge_closed_by_vertex(&vs, et.e[j].a, et.e[j].b, tol))
                tclosed++;
        }
    }

    long unmatched = singles - tclosed;
    printf("verify: %u facets, %zu vertices, %ld degenerate\n",
           nfacets, vs.count, degenerate);
    printf("verify: bbox x [%.6g, %.6g] y [%.6g, %.6g] z [%.6g, %.6g]\n",
           bbmin[0], bbmax[0], bbmin[1], bbmax[1], bbmin[2], bbmax[2]);
    printf("verify: %ld edges shared by 2 facets, %ld T-junction edges "
           "(closed within %g), %ld pinch edges, %ld over-shared, %ld unmatched\n",
           exact2, tclosed, tol, pinch, over2, unmatched);
    printf("verify: orientation: %ld inconsistently wound shared edges, "
           "signed volume %.6g\n",
           badorient, volume);

    int pass = (nfacets > 0) && over2 == 0 && unmatched == 0 &&
               degenerate == 0 && badorient == 0 && volume > 0.0;
    printf("verify: %s\n", pass ? "PASS" : "FAIL");

    free(vs.v);
    free(vs.idb);
    free(vs.xyb);
    free(et.e);
    return pass ? 0 : 1;
}