using System;
using System.Collections.Generic;
using UnityEngine;
#if UNITY_EDITOR
using UnityEditor;
using UnityEditor.SceneManagement;
#endif

/// <summary>
/// Editor tool that decomposes an environment mesh (typically an open-shell / non-watertight
/// furniture mesh) into a small set of CONVEX MeshColliders via surface clustering, so
/// <see cref="ConcealmentMeter"/> can treat the furniture as "embeddable" through its
/// convex-collider ClosestPoint path WITHOUT a single convex hull filling the open compartments.
///
/// Approach (no voxels, no native plugin — "approximate convex decomposition" on the surface):
///  1. Weld coincident vertices so connectivity is detected across split-normal/UV seams.
///  2. Grow connected triangle regions into SLABS: a region keeps absorbing edge-adjacent
///     triangles while it stays thinner than <see cref="maxSlabThickness"/> along its flattest
///     axis (PCA minor axis). This spans both faces of one board (they are within a board's
///     thickness) but refuses to jump the gap to the next board (that would thicken the region).
///  3. Convex-hull each slab into its own convex collider. Boards become thin convex slabs;
///     the empty space between boards is covered by no collider, so it stays hollow.
///
/// The colliders are real child GameObjects referencing SAVED hull-mesh assets, so they survive
/// domain reload / play-mode / scene save (unlike an in-memory-only bake).
///
/// Runtime-inert: the component exists only so <see cref="ConcealmentMeter"/> can detect it in a
/// collider's parent and skip the furniture's own concave MeshCollider (exactly like a baked
/// <c>ConcealmentVolume</c> does) — the generated convex slabs are the authoritative embedding shape.
/// </summary>
[DisallowMultipleComponent]
public class ConcealmentColliderBaker : MonoBehaviour
{
    [Tooltip("Optional. If empty, the shared mesh of a MeshFilter on this object (or a child) is used.")]
    [SerializeField] private Mesh sourceMesh;

    [Header("Clustering")]
    [Tooltip("Vertices closer than this (local units) are welded, so shared edges are found across split-normal/UV seams. ~ a hair below the smallest real feature.")]
    [Min(1e-4f)] [SerializeField] private float weldEpsilon = 1e-3f;

    [Tooltip("A connected region keeps growing while it stays thinner than this (local units) along its flattest axis. ~ board thickness. Larger => fewer, chunkier colliders (risk of bridging compartments); smaller => more, thinner slabs.")]
    [Min(1e-3f)] [SerializeField] private float maxSlabThickness = 0.12f;

    [Tooltip("Flat/thin clusters are inflated to at least this thickness (local units) so the convex collider has usable volume for ClosestPoint tests.")]
    [Min(1e-3f)] [SerializeField] private float minThickness = 0.04f;

    [Tooltip("A slab is recursively split until its convex-hull footprint is at least this fraction covered by real surface (0..1). Guards against one hull bridging open gaps between coplanar members (e.g. a table side = 2 legs + top rail merge into a flat panel whose hull fills the empty gap). Lower => tolerate more concavity/bridging; higher => split into more, tighter surface-hugging pieces. Solid boards (~fully covered) are never split.")]
    [Range(0f, 1f)] [SerializeField] private float hullFillThreshold = 0.6f;

    [Tooltip("Clusters with fewer triangles than this are discarded as noise.")]
    [Min(1)] [SerializeField] private int minClusterTriangles = 1;

    [Tooltip("Hard cap on generated colliders. The meter's broadphase is capped at 64 colliders TOTAL across all environment in range, so keep this modest. Largest clusters win.")]
    [Min(1)] [SerializeField] private int maxColliders = 24;

    [Header("Output")]
    [Tooltip("If true, generated colliders inherit the source object's layer (so they land on the meter's environment mask). Otherwise 'Collider Layer' is used.")]
    [SerializeField] private bool useSourceLayer = true;
    [SerializeField] private int colliderLayer;

    [SerializeField, HideInInspector] private List<string> generatedAssetPaths = new List<string>();

    private const string ContainerName = "ConcealmentColliders (generated)";
    private const string AssetFolder = "Assets/Generated/ConcealmentColliders";
    private const int MaxRefineDepth = 8;

#if UNITY_EDITOR
    [ContextMenu("Bake Convex Colliders")]
    public void Bake()
    {
        Mesh mesh = sourceMesh;
        Transform meshSpace = transform;
        if (!mesh)
        {
            MeshFilter mf = GetComponentInChildren<MeshFilter>();
            if (mf && mf.sharedMesh) { mesh = mf.sharedMesh; meshSpace = mf.transform; }
        }
        if (!mesh)
        {
            Debug.LogError($"[ConcealmentColliderBaker] '{name}': no source mesh (assign one or add a MeshFilter).", this);
            return;
        }

        Clear(); // remove previous container + assets first

        string baseName = ResolveAssetBaseName();

        Vector3[] rawVerts = mesh.vertices;
        int[] rawTris = mesh.triangles;

        // 1. Weld coincident vertices.
        WeldVertices(rawVerts, weldEpsilon, out List<Vector3> wv, out int[] map);

        // Welded triangles (drop degenerate).
        var triA = new List<int>(); var triB = new List<int>(); var triC = new List<int>();
        for (int i = 0; i + 2 < rawTris.Length; i += 3)
        {
            int a = map[rawTris[i]], b = map[rawTris[i + 1]], c = map[rawTris[i + 2]];
            if (a == b || b == c || a == c) continue;
            triA.Add(a); triB.Add(b); triC.Add(c);
        }
        int triCount = triA.Count;
        if (triCount == 0)
        {
            Debug.LogWarning($"[ConcealmentColliderBaker] '{name}': mesh has no valid triangles after welding.", this);
            return;
        }

        // 2. Edge adjacency (welded edges shared by triangles).
        var edgeToTri = new Dictionary<long, int>(triCount * 3);
        var neighbors = new List<int>[triCount];
        for (int t = 0; t < triCount; t++) neighbors[t] = new List<int>();
        for (int t = 0; t < triCount; t++)
        {
            LinkEdge(edgeToTri, neighbors, triA[t], triB[t], t);
            LinkEdge(edgeToTri, neighbors, triB[t], triC[t], t);
            LinkEdge(edgeToTri, neighbors, triC[t], triA[t], t);
        }

        // 3. Grow connected triangles into thickness-bounded slabs.
        //    Disconnected boards never link (separate components); connected assemblies split
        //    wherever a region would exceed maxSlabThickness.
        var clusters = new List<List<int>>();
        bool[] assigned = new bool[triCount];
        var queue = new Queue<int>();
        var vset = new HashSet<int>();
        var vpos = new List<Vector3>();
        var added = new List<int>();

        for (int s = 0; s < triCount; s++)
        {
            if (assigned[s]) continue;

            var cluster = new List<int> { s };
            assigned[s] = true;
            vset.Clear(); vpos.Clear();
            AddTriVerts(vset, vpos, wv, triA[s], triB[s], triC[s], null);

            queue.Clear();
            for (int k = 0; k < neighbors[s].Count; k++) queue.Enqueue(neighbors[s][k]);

            while (queue.Count > 0)
            {
                int t = queue.Dequeue();
                if (assigned[t]) continue;

                int before = vpos.Count;
                AddTriVerts(vset, vpos, wv, triA[t], triB[t], triC[t], added);

                float th = SlabThickness(vpos, out _);
                if (th <= maxSlabThickness)
                {
                    assigned[t] = true;
                    cluster.Add(t);
                    for (int k = 0; k < neighbors[t].Count; k++)
                        if (!assigned[neighbors[t][k]]) queue.Enqueue(neighbors[t][k]);
                }
                else
                {
                    // Revert the tentatively added verts.
                    for (int k = 0; k < added.Count; k++) vset.Remove(added[k]);
                    if (vpos.Count > before) vpos.RemoveRange(before, vpos.Count - before);
                }
            }

            clusters.Add(cluster);
        }

        // 3b. Refine: split slabs whose convex hull would bridge open concavities. A table side is
        //     2 legs + a top rail — all coplanar and thin, so they merge into ONE slab whose hull
        //     fills the empty gap between the legs. Each slab is recursively halved along its larger
        //     in-plane axis until its hull footprint is mostly covered by real surface. Solid boards
        //     (the shelf) are already ~fully covered => returned untouched.
        var refined = new List<List<int>>(clusters.Count);
        for (int c = 0; c < clusters.Count; c++)
            RefineCluster(clusters[c], wv, triA, triB, triC, 0, refined);
        clusters = refined;

        // 4. Filter noise + cap count (keep largest slabs).
        clusters.RemoveAll(c => c.Count < minClusterTriangles);
        clusters.Sort((x, y) => y.Count.CompareTo(x.Count));
        int dropped = 0;
        if (clusters.Count > maxColliders)
        {
            dropped = clusters.Count - maxColliders;
            clusters.RemoveRange(maxColliders, dropped);
        }

        // 5. Build a convex collider per slab.
        EnsureFolder();
        Transform container = new GameObject(ContainerName).transform;
        container.SetParent(meshSpace, false);
        container.localPosition = Vector3.zero;
        container.localRotation = Quaternion.identity;
        container.localScale = Vector3.one;
        int outLayer = useSourceLayer ? gameObject.layer : colliderLayer;

        int made = 0;
        for (int ci = 0; ci < clusters.Count; ci++)
        {
            // Unique verts of this cluster.
            var pset = new HashSet<int>();
            var pts = new List<Vector3>();
            List<int> cl = clusters[ci];
            for (int j = 0; j < cl.Count; j++)
            {
                int t = cl[j];
                if (pset.Add(triA[t])) pts.Add(wv[triA[t]]);
                if (pset.Add(triB[t])) pts.Add(wv[triB[t]]);
                if (pset.Add(triC[t])) pts.Add(wv[triC[t]]);
            }

            float thick = SlabThickness(pts, out Vector3 normal);
            List<Vector3> hullPts;
            if (thick < minThickness)
            {
                // Inflate a (near-)flat slab into a thin prism so the hull has volume.
                float half = minThickness * 0.5f;
                hullPts = new List<Vector3>(pts.Count * 2);
                for (int k = 0; k < pts.Count; k++)
                {
                    hullPts.Add(pts[k] + normal * half);
                    hullPts.Add(pts[k] - normal * half);
                }
            }
            else hullPts = pts;

            if (!ComputeConvexHull(hullPts.ToArray(), out List<Vector3> hv, out List<int> ht) || hv.Count < 4)
                continue;

            Mesh cm = new Mesh { name = $"{baseName}_{made}" };
            cm.SetVertices(hv);
            cm.SetTriangles(ht, 0);
            cm.RecalculateBounds();

            string path = $"{AssetFolder}/{baseName}_{made}.asset";
            AssetDatabase.CreateAsset(cm, path);
            generatedAssetPaths.Add(path);

            GameObject go = new GameObject($"hull_{made}");
            go.transform.SetParent(container, false);
            go.layer = outLayer;
            MeshCollider col = go.AddComponent<MeshCollider>();
            col.sharedMesh = cm;
            col.convex = true;
            col.isTrigger = false;

            made++;
        }

        AssetDatabase.SaveAssets();
        EditorUtility.SetDirty(this);
        UnityEditor.SceneManagement.EditorSceneManager.MarkSceneDirty(gameObject.scene);

        Debug.Log($"[ConcealmentColliderBaker] '{name}': baked {made} convex collider(s)" +
                  (dropped > 0 ? $" (+{dropped} smaller slab(s) dropped by Max Colliders={maxColliders})" : "") +
                  ". Save the scene to persist.", this);
    }

    [ContextMenu("Clear Baked Colliders")]
    public void Clear()
    {
        for (int i = 0; i < generatedAssetPaths.Count; i++)
            if (!string.IsNullOrEmpty(generatedAssetPaths[i]))
                AssetDatabase.DeleteAsset(generatedAssetPaths[i]);
        generatedAssetPaths.Clear();

        var all = GetComponentsInChildren<Transform>(true);
        for (int i = 0; i < all.Length; i++)
            if (all[i] && all[i] != transform && all[i].name == ContainerName)
                DestroyImmediate(all[i].gameObject);

        EditorUtility.SetDirty(this);
        if (gameObject.scene.IsValid())
            UnityEditor.SceneManagement.EditorSceneManager.MarkSceneDirty(gameObject.scene);
    }

    private static void EnsureFolder()
    {
        if (!AssetDatabase.IsValidFolder("Assets/Generated"))
            AssetDatabase.CreateFolder("Assets", "Generated");
        if (!AssetDatabase.IsValidFolder(AssetFolder))
            AssetDatabase.CreateFolder("Assets/Generated", "ConcealmentColliders");
    }

    /// <summary>
    /// Readable, deterministic base name for the generated collider-mesh assets (replaces the old GUID):
    ///  - Prefab stage: "{assetPath without 'Assets/' and '.prefab'}/{path below the prefab root}"
    ///  - Scene:        "{sceneName}/{full hierarchy path from the scene root}"
    /// Sanitized to a filesystem-safe token; on collision with a DIFFERENT baker's assets, "_2", "_3", … is appended.
    /// </summary>
    private string ResolveAssetBaseName()
    {
        string raw;
        PrefabStage stage = PrefabStageUtility.GetPrefabStage(gameObject);
        if (stage != null)
        {
            string ap = stage.assetPath;
            if (ap.StartsWith("Assets/")) ap = ap.Substring("Assets/".Length);
            if (ap.EndsWith(".prefab")) ap = ap.Substring(0, ap.Length - ".prefab".Length);
            // Path below the prefab root — the root's name already lives in the asset path.
            string hier = GetHierarchyPath(transform, stage.prefabContentsRoot.transform);
            raw = hier.Length == 0 ? ap : ap + "/" + hier;
        }
        else
        {
            string scene = string.IsNullOrEmpty(gameObject.scene.name) ? "Untitled" : gameObject.scene.name;
            raw = scene + "/" + GetHierarchyPath(transform, null);
        }

        string baseName = SanitizeFileName(raw);
        if (baseName.Length == 0) baseName = "Bake";

        // Clear() already removed THIS baker's assets, so a hit here means another baker resolved to the
        // same name -> disambiguate. The per-collider "_0" file is the presence marker of a base name.
        string candidate = baseName;
        for (int dup = 2; System.IO.File.Exists($"{AssetFolder}/{candidate}_0.asset"); dup++)
            candidate = $"{baseName}_{dup}";
        return candidate;
    }

    /// <summary>Slash-joined transform names from <paramref name="t"/> up to (but excluding) <paramref name="stopAt"/> (null = up to the root).</summary>
    private static string GetHierarchyPath(Transform t, Transform stopAt)
    {
        var names = new List<string>();
        for (Transform cur = t; cur && cur != stopAt; cur = cur.parent)
            names.Add(cur.name);
        names.Reverse();
        return string.Join("/", names);
    }

    /// <summary>Collapses every run of non-[A-Za-z0-9-] characters to a single '_' and trims leading/trailing '_'.</summary>
    private static string SanitizeFileName(string s)
    {
        var sb = new System.Text.StringBuilder(s.Length);
        bool pendingSep = false;
        for (int i = 0; i < s.Length; i++)
        {
            char c = s[i];
            bool keep = (c >= 'a' && c <= 'z') || (c >= 'A' && c <= 'Z') || (c >= '0' && c <= '9') || c == '-';
            if (keep)
            {
                if (pendingSep && sb.Length > 0) sb.Append('_');
                pendingSep = false;
                sb.Append(c);
            }
            else pendingSep = true;
        }
        return sb.ToString();
    }
#endif

    // ---- Welding & clustering helpers (runtime-safe: no editor API) ----

    private static void WeldVertices(Vector3[] src, float eps, out List<Vector3> welded, out int[] map)
    {
        welded = new List<Vector3>(src.Length);
        map = new int[src.Length];
        float inv = 1f / Mathf.Max(eps, 1e-6f);
        var lut = new Dictionary<long, int>(src.Length);
        for (int i = 0; i < src.Length; i++)
        {
            Vector3 p = src[i];
            int qx = Mathf.RoundToInt(p.x * inv);
            int qy = Mathf.RoundToInt(p.y * inv);
            int qz = Mathf.RoundToInt(p.z * inv);
            long key = ((long)(qx & 0x1FFFFF) << 42) | ((long)(qy & 0x1FFFFF) << 21) | (long)(qz & 0x1FFFFF);
            if (lut.TryGetValue(key, out int idx)) { map[i] = idx; }
            else { idx = welded.Count; lut[key] = idx; welded.Add(p); map[i] = idx; }
        }
    }

    private static void LinkEdge(Dictionary<long, int> edgeToTri, List<int>[] neighbors, int u, int v, int tri)
    {
        long key = EdgeKey(u, v);
        if (edgeToTri.TryGetValue(key, out int other))
        {
            neighbors[tri].Add(other);
            neighbors[other].Add(tri);
            edgeToTri[key] = tri; // chain non-manifold edges (>2 tris) so all get linked
        }
        else edgeToTri[key] = tri;
    }

    private static void AddTriVerts(HashSet<int> set, List<Vector3> pos, List<Vector3> wv, int a, int b, int c, List<int> addedOut)
    {
        addedOut?.Clear();
        if (set.Add(a)) { pos.Add(wv[a]); addedOut?.Add(a); }
        if (set.Add(b)) { pos.Add(wv[b]); addedOut?.Add(b); }
        if (set.Add(c)) { pos.Add(wv[c]); addedOut?.Add(c); }
    }

    /// <summary>Extent of the point cloud along its PCA minor axis (= slab thickness), plus that axis.</summary>
    private static float SlabThickness(List<Vector3> pts, out Vector3 normal)
    {
        normal = Vector3.up;
        int n = pts.Count;
        if (n < 2) return 0f;

        Vector3 mean = Vector3.zero;
        for (int i = 0; i < n; i++) mean += pts[i];
        mean /= n;

        double c00 = 0, c01 = 0, c02 = 0, c11 = 0, c12 = 0, c22 = 0;
        for (int i = 0; i < n; i++)
        {
            double x = pts[i].x - mean.x, y = pts[i].y - mean.y, z = pts[i].z - mean.z;
            c00 += x * x; c01 += x * y; c02 += x * z;
            c11 += y * y; c12 += y * z; c22 += z * z;
        }

        var a = new double[3, 3] { { c00, c01, c02 }, { c01, c11, c12 }, { c02, c12, c22 } };
        Jacobi3(a, out double[,] v, out double[] d);

        int mi = 0;
        if (d[1] < d[mi]) mi = 1;
        if (d[2] < d[mi]) mi = 2;
        normal = new Vector3((float)v[0, mi], (float)v[1, mi], (float)v[2, mi]);
        if (normal.sqrMagnitude < 1e-12f) { normal = Vector3.up; return 0f; }
        normal.Normalize();

        float mn = float.MaxValue, mx = float.MinValue;
        for (int i = 0; i < n; i++)
        {
            float dp = Vector3.Dot(pts[i], normal);
            if (dp < mn) mn = dp;
            if (dp > mx) mx = dp;
        }
        return mx - mn;
    }

    // ---- Concavity refinement (split slabs whose hull would bridge open gaps) ----

    /// <summary>
    /// Recursively halves a slab cluster along its larger in-plane axis until its hull footprint is
    /// at least <see cref="hullFillThreshold"/> covered by real surface (or it gets too small / too deep).
    /// Prevents one convex hull from bridging the empty space between disjoint coplanar members (table legs).
    /// </summary>
    private void RefineCluster(List<int> cluster, List<Vector3> wv, List<int> triA, List<int> triB, List<int> triC, int depth, List<List<int>> outClusters)
    {
        if (cluster.Count <= 2 || depth >= MaxRefineDepth)
        {
            outClusters.Add(cluster);
            return;
        }

        float fill = ClusterFill(cluster, wv, triA, triB, triC, out Vector3 axis);
        if (fill >= hullFillThreshold)
        {
            outClusters.Add(cluster);
            return;
        }

        // Halve at the centroid median along the larger in-plane axis. Recursion + the fill test
        // converge on surface-hugging pieces; a solid board is >= threshold at once and never splits.
        int m = cluster.Count;
        var proj = new float[m];
        for (int j = 0; j < m; j++)
            proj[j] = Vector3.Dot((wv[triA[cluster[j]]] + wv[triB[cluster[j]]] + wv[triC[cluster[j]]]) * (1f / 3f), axis);
        var sorted = (float[])proj.Clone();
        Array.Sort(sorted);
        float median = sorted[m / 2];

        var left = new List<int>();
        var right = new List<int>();
        for (int j = 0; j < m; j++)
        {
            if (proj[j] < median) left.Add(cluster[j]); else right.Add(cluster[j]);
        }
        if (left.Count == 0 || right.Count == 0)
        {
            outClusters.Add(cluster); // all centroids coincide on the axis — can't halve, keep as-is
            return;
        }

        RefineCluster(left, wv, triA, triB, triC, depth + 1, outClusters);
        RefineCluster(right, wv, triA, triB, triC, depth + 1, outClusters);
    }

    /// <summary>
    /// Fraction (0..1) of the cluster's in-plane bounding footprint actually covered by surface, via a
    /// coverage grid on the PCA plane (union of the projected triangles — single-/double-sided agnostic).
    /// Low fill => the convex hull would enclose empty space. Outputs the larger in-plane axis (split dir).
    /// </summary>
    private static float ClusterFill(List<int> cluster, List<Vector3> wv, List<int> triA, List<int> triB, List<int> triC, out Vector3 splitAxis)
    {
        var pts = new List<Vector3>(cluster.Count * 3);
        var seen = new HashSet<int>();
        for (int j = 0; j < cluster.Count; j++)
        {
            int t = cluster[j];
            if (seen.Add(triA[t])) pts.Add(wv[triA[t]]);
            if (seen.Add(triB[t])) pts.Add(wv[triB[t]]);
            if (seen.Add(triC[t])) pts.Add(wv[triC[t]]);
        }

        PcaFrame(pts, out Vector3 u, out Vector3 v, out _);

        float umin = float.MaxValue, umax = float.MinValue, vmin = float.MaxValue, vmax = float.MinValue;
        for (int i = 0; i < pts.Count; i++)
        {
            float pu = Vector3.Dot(pts[i], u), pv = Vector3.Dot(pts[i], v);
            if (pu < umin) umin = pu; if (pu > umax) umax = pu;
            if (pv < vmin) vmin = pv; if (pv > vmax) vmax = pv;
        }
        float eu = Mathf.Max(umax - umin, 1e-5f), ev = Mathf.Max(vmax - vmin, 1e-5f);
        splitAxis = eu >= ev ? u : v;

        float cell = Mathf.Max(Mathf.Max(eu, ev) / 32f, 1e-4f);
        int NU = Mathf.Clamp(Mathf.CeilToInt(eu / cell), 1, 40);
        int NV = Mathf.Clamp(Mathf.CeilToInt(ev / cell), 1, 40);
        var cov = new bool[NU, NV];
        for (int j = 0; j < cluster.Count; j++)
        {
            int t = cluster[j];
            RasterTri(cov, NU, NV,
                ToGrid(wv[triA[t]], u, v, umin, vmin, eu, ev, NU, NV),
                ToGrid(wv[triB[t]], u, v, umin, vmin, eu, ev, NU, NV),
                ToGrid(wv[triC[t]], u, v, umin, vmin, eu, ev, NU, NV));
        }

        int covered = 0;
        for (int a = 0; a < NU; a++) for (int b = 0; b < NV; b++) if (cov[a, b]) covered++;
        return (float)covered / (NU * NV);
    }

    private static Vector2 ToGrid(Vector3 p, Vector3 u, Vector3 v, float umin, float vmin, float eu, float ev, int NU, int NV)
        => new Vector2((Vector3.Dot(p, u) - umin) / eu * NU, (Vector3.Dot(p, v) - vmin) / ev * NV);

    /// <summary>Marks every grid cell whose center falls inside the 2D triangle (edge-on triangles paint their thin bbox strip).</summary>
    private static void RasterTri(bool[,] cov, int NU, int NV, Vector2 a, Vector2 b, Vector2 c)
    {
        float minx = Mathf.Min(a.x, Mathf.Min(b.x, c.x)), maxx = Mathf.Max(a.x, Mathf.Max(b.x, c.x));
        float miny = Mathf.Min(a.y, Mathf.Min(b.y, c.y)), maxy = Mathf.Max(a.y, Mathf.Max(b.y, c.y));
        int x0 = Mathf.Clamp(Mathf.FloorToInt(minx), 0, NU - 1), x1 = Mathf.Clamp(Mathf.CeilToInt(maxx), 0, NU - 1);
        int y0 = Mathf.Clamp(Mathf.FloorToInt(miny), 0, NV - 1), y1 = Mathf.Clamp(Mathf.CeilToInt(maxy), 0, NV - 1);

        float d = (b.y - c.y) * (a.x - c.x) + (c.x - b.x) * (a.y - c.y);
        if (Mathf.Abs(d) < 1e-9f)
        {
            for (int x = x0; x <= x1; x++) for (int y = y0; y <= y1; y++) cov[x, y] = true;
            return;
        }
        for (int x = x0; x <= x1; x++)
            for (int y = y0; y <= y1; y++)
            {
                float px = x + 0.5f, py = y + 0.5f;
                float w0 = ((b.y - c.y) * (px - c.x) + (c.x - b.x) * (py - c.y)) / d;
                float w1 = ((c.y - a.y) * (px - c.x) + (a.x - c.x) * (py - c.y)) / d;
                float w2 = 1f - w0 - w1;
                if (w0 >= -1e-3f && w1 >= -1e-3f && w2 >= -1e-3f) cov[x, y] = true;
            }
    }

    /// <summary>PCA orthonormal frame of a point cloud: u = major axis, v = mid, n = minor (slab normal). Eigenvectors sorted by variance desc.</summary>
    private static void PcaFrame(List<Vector3> pts, out Vector3 u, out Vector3 v, out Vector3 n)
    {
        u = Vector3.right; v = Vector3.up; n = Vector3.forward;
        int cnt = pts.Count;
        if (cnt < 2) return;

        Vector3 mean = Vector3.zero;
        for (int i = 0; i < cnt; i++) mean += pts[i];
        mean /= cnt;

        double c00 = 0, c01 = 0, c02 = 0, c11 = 0, c12 = 0, c22 = 0;
        for (int i = 0; i < cnt; i++)
        {
            double x = pts[i].x - mean.x, y = pts[i].y - mean.y, z = pts[i].z - mean.z;
            c00 += x * x; c01 += x * y; c02 += x * z;
            c11 += y * y; c12 += y * z; c22 += z * z;
        }

        var mat = new double[3, 3] { { c00, c01, c02 }, { c01, c11, c12 }, { c02, c12, c22 } };
        Jacobi3(mat, out double[,] vec, out double[] d);

        int hi = 0, mid = 1, lo = 2;
        if (d[hi] < d[mid]) { int t = hi; hi = mid; mid = t; }
        if (d[mid] < d[lo]) { int t = mid; mid = lo; lo = t; }
        if (d[hi] < d[mid]) { int t = hi; hi = mid; mid = t; }

        u = AxisOr(vec, hi, Vector3.right);
        v = AxisOr(vec, mid, Vector3.up);
        n = AxisOr(vec, lo, Vector3.forward);
    }

    private static Vector3 AxisOr(double[,] vec, int col, Vector3 fallback)
    {
        Vector3 a = new Vector3((float)vec[0, col], (float)vec[1, col], (float)vec[2, col]);
        return a.sqrMagnitude < 1e-12f ? fallback : a.normalized;
    }

    /// <summary>Cyclic Jacobi eigendecomposition of a symmetric 3x3 matrix (eigenvectors in columns of v).</summary>
    private static void Jacobi3(double[,] a, out double[,] v, out double[] d)
    {
        v = new double[3, 3] { { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } };
        for (int sweep = 0; sweep < 24; sweep++)
        {
            double off = Math.Abs(a[0, 1]) + Math.Abs(a[0, 2]) + Math.Abs(a[1, 2]);
            if (off < 1e-18) break;
            JacobiRotate(a, v, 0, 1);
            JacobiRotate(a, v, 0, 2);
            JacobiRotate(a, v, 1, 2);
        }
        d = new double[] { a[0, 0], a[1, 1], a[2, 2] };
    }

    private static void JacobiRotate(double[,] a, double[,] v, int p, int q)
    {
        double apq = a[p, q];
        if (Math.Abs(apq) < 1e-20) return;

        double app = a[p, p], aqq = a[q, q];
        double phi = 0.5 * Math.Atan2(2.0 * apq, aqq - app);
        double c = Math.Cos(phi), s = Math.Sin(phi);

        a[p, p] = c * c * app - 2.0 * s * c * apq + s * s * aqq;
        a[q, q] = s * s * app + 2.0 * s * c * apq + c * c * aqq;
        a[p, q] = 0; a[q, p] = 0;

        int r = 3 - p - q;
        double arp = a[r, p], arq = a[r, q];
        double nrp = c * arp - s * arq;
        double nrq = s * arp + c * arq;
        a[r, p] = nrp; a[p, r] = nrp;
        a[r, q] = nrq; a[q, r] = nrq;

        for (int k = 0; k < 3; k++)
        {
            double vkp = v[k, p], vkq = v[k, q];
            v[k, p] = c * vkp - s * vkq;
            v[k, q] = s * vkp + c * vkq;
        }
    }

    // ---- Convex hull (incremental; outputs a triangulated hull mesh) ----

    private struct HFace { public int a, b, c; public Vector3 n; }

    private static bool ComputeConvexHull(Vector3[] pts, out List<Vector3> hullVertsOut, out List<int> hullTrisOut)
    {
        hullVertsOut = new List<Vector3>();
        hullTrisOut = new List<int>();

        int n = pts.Length;
        if (n < 4) return false;

        Vector3 bmn = pts[0], bmx = pts[0];
        for (int i = 1; i < n; i++) { bmn = Vector3.Min(bmn, pts[i]); bmx = Vector3.Max(bmx, pts[i]); }
        float diag = (bmx - bmn).magnitude;
        if (diag < 1e-6f) return false;
        float eps = diag * 1e-6f;

        int minX = 0, maxX = 0;
        for (int i = 1; i < n; i++) { if (pts[i].x < pts[minX].x) minX = i; if (pts[i].x > pts[maxX].x) maxX = i; }
        int i0 = minX, i1 = maxX;
        if (i0 == i1) return false;

        Vector3 lineDir = pts[i1] - pts[i0];
        float lineLen = lineDir.magnitude;
        if (lineLen < eps) return false;
        lineDir /= lineLen;

        int i2 = -1; float best = eps;
        for (int i = 0; i < n; i++)
        {
            float dd = Vector3.Cross(pts[i] - pts[i0], lineDir).magnitude;
            if (dd > best) { best = dd; i2 = i; }
        }
        if (i2 < 0) return false;

        Vector3 nrm0 = Vector3.Cross(pts[i1] - pts[i0], pts[i2] - pts[i0]);
        float nl = nrm0.magnitude;
        if (nl < eps) return false;
        nrm0 /= nl;

        int i3 = -1; best = eps;
        for (int i = 0; i < n; i++)
        {
            float dd = Mathf.Abs(Vector3.Dot(pts[i] - pts[i0], nrm0));
            if (dd > best) { best = dd; i3 = i; }
        }
        if (i3 < 0) return false;

        Vector3 interior = (pts[i0] + pts[i1] + pts[i2] + pts[i3]) * 0.25f;

        var faces = new List<HFace>(4 + n);
        AddHullFace(faces, pts, interior, i0, i1, i2);
        AddHullFace(faces, pts, interior, i0, i1, i3);
        AddHullFace(faces, pts, interior, i0, i2, i3);
        AddHullFace(faces, pts, interior, i1, i2, i3);

        var visible = new List<int>();
        var edgeCount = new Dictionary<long, int>();
        var horizonA = new List<int>();
        var horizonB = new List<int>();

        for (int pi = 0; pi < n; pi++)
        {
            if (pi == i0 || pi == i1 || pi == i2 || pi == i3) continue;
            Vector3 p = pts[pi];

            visible.Clear();
            for (int f = 0; f < faces.Count; f++)
            {
                HFace hf = faces[f];
                if (Vector3.Dot(hf.n, p - pts[hf.a]) > eps) visible.Add(f);
            }
            if (visible.Count == 0) continue;

            edgeCount.Clear();
            for (int vi = 0; vi < visible.Count; vi++)
            {
                HFace hf = faces[visible[vi]];
                AccEdge(edgeCount, hf.a, hf.b);
                AccEdge(edgeCount, hf.b, hf.c);
                AccEdge(edgeCount, hf.c, hf.a);
            }
            horizonA.Clear(); horizonB.Clear();
            for (int vi = 0; vi < visible.Count; vi++)
            {
                HFace hf = faces[visible[vi]];
                if (EdgeIsHorizon(edgeCount, hf.a, hf.b)) { horizonA.Add(hf.a); horizonB.Add(hf.b); }
                if (EdgeIsHorizon(edgeCount, hf.b, hf.c)) { horizonA.Add(hf.b); horizonB.Add(hf.c); }
                if (EdgeIsHorizon(edgeCount, hf.c, hf.a)) { horizonA.Add(hf.c); horizonB.Add(hf.a); }
            }

            visible.Sort();
            for (int vi = visible.Count - 1; vi >= 0; vi--) faces.RemoveAt(visible[vi]);

            for (int h = 0; h < horizonA.Count; h++)
                AddHullFace(faces, pts, interior, horizonA[h], horizonB[h], pi);
        }

        if (faces.Count < 4) return false;

        var remap = new Dictionary<int, int>();
        for (int f = 0; f < faces.Count; f++)
        {
            HFace hf = faces[f];
            hullTrisOut.Add(RemapVert(remap, hullVertsOut, pts, hf.a));
            hullTrisOut.Add(RemapVert(remap, hullVertsOut, pts, hf.b));
            hullTrisOut.Add(RemapVert(remap, hullVertsOut, pts, hf.c));
        }
        return true;
    }

    private static void AddHullFace(List<HFace> faces, Vector3[] pts, Vector3 interior, int a, int b, int c)
    {
        Vector3 nn = Vector3.Cross(pts[b] - pts[a], pts[c] - pts[a]);
        float l = nn.magnitude;
        if (l < 1e-20f) return;
        nn /= l;
        if (Vector3.Dot(nn, pts[a] - interior) < 0f) nn = -nn;
        faces.Add(new HFace { a = a, b = b, c = c, n = nn });
    }

    private static long EdgeKey(int u, int v)
    {
        int lo = u < v ? u : v, hi = u < v ? v : u;
        return ((long)lo << 32) | (uint)hi;
    }

    private static void AccEdge(Dictionary<long, int> ec, int u, int v)
    {
        long k = EdgeKey(u, v);
        ec.TryGetValue(k, out int c);
        ec[k] = c + 1;
    }

    private static bool EdgeIsHorizon(Dictionary<long, int> ec, int u, int v)
    {
        ec.TryGetValue(EdgeKey(u, v), out int c);
        return c == 1;
    }

    private static int RemapVert(Dictionary<int, int> remap, List<Vector3> outV, Vector3[] pts, int idx)
    {
        if (remap.TryGetValue(idx, out int r)) return r;
        r = outV.Count; remap[idx] = r; outV.Add(pts[idx]);
        return r;
    }
}