"""Printability check on the exported STLs (already in print orientation, bed at z = 0).

Reports downward-facing surfaces that are not on the bed, split into
  - bridges     (horizontal ceilings, |n.z| > 0.98): fine if the span is short
  - overhangs   (steeper than 45 deg from vertical, i.e. n.z < -0.707)
and estimates the longest unsupported span of each bridge island.
"""
import sys
import struct
import numpy as np
from scipy.sparse.csgraph import connected_components
from scipy.sparse import coo_matrix


def load_stl(path):
    with open(path, "rb") as f:
        f.read(80)
        n = struct.unpack("<I", f.read(4))[0]
        dt = np.dtype([("n", "<f4", 3), ("v", "<f4", (3, 3)), ("a", "<u2")])
        return np.fromfile(f, dtype=dt, count=n)["v"].astype(float)


def islands(T, mask):
    idx = np.nonzero(mask)[0]
    if len(idx) == 0:
        return []
    V = np.round(T[idx].reshape(-1, 3), 4)
    _, inv = np.unique(V, axis=0, return_inverse=True)
    inv = inv.reshape(-1, 3)
    rows = np.repeat(np.arange(len(idx)), 3)
    g = coo_matrix((np.ones(len(rows)), (rows, inv.ravel())), shape=(len(idx), inv.max() + 1))
    adj = g @ g.T
    n, lab = connected_components(adj, directed=False)
    return [idx[lab == k] for k in range(n)]


def report(path, max_bridge=12.0):
    T = load_stl(path)
    zmin = T[..., 2].min()
    n = np.cross(T[:, 1] - T[:, 0], T[:, 2] - T[:, 0])
    area = 0.5 * np.linalg.norm(n, axis=1)
    n = n / (np.linalg.norm(n, axis=1, keepdims=True) + 1e-12)
    zc = T[..., 2].mean(1)
    on_bed = zc < zmin + 0.05
    down = (n[:, 2] < -0.707) & ~on_bed
    bridge = down & (n[:, 2] < -0.98)
    steep = down & ~bridge
    print(f"\n{path.split('/')[-1]}: {len(T)} tris, height {np.ptp(T[..., 2]):.1f} mm")
    print(f"  overhang >45deg (non-horizontal): {area[steep].sum():7.1f} mm2")
    print(f"  horizontal ceilings (bridges):    {area[bridge].sum():7.1f} mm2")
    bad = []
    for isl in islands(T, bridge):
        P = T[isl].reshape(-1, 2 if False else 3)[:, :2]
        z = T[isl][..., 2].mean()
        a = area[isl].sum()
        # span = smaller side of the oriented bounding box (via PCA)
        c = P - P.mean(0)
        if len(c) > 2:
            _, _, vt = np.linalg.svd(c, full_matrices=False)
            proj = c @ vt.T
            ext = np.ptp(proj, axis=0)
            span = ext.min() if len(ext) > 1 else 0
        else:
            span = 0
        if a > 0.5:
            bad.append((span, a, z, P.mean(0)))
    bad.sort(key=lambda b: -b[0])
    for span, a, z, c in bad[:12]:
        flag = "  <-- long" if span > max_bridge else ""
        print(f"    bridge at z={z:5.2f}  area {a:6.1f} mm2  span~{span:5.1f} mm  near ({c[0]:.1f},{c[1]:.1f}){flag}")
    for isl in islands(T, steep):
        a = area[isl].sum()
        if a > 2.0:
            P = T[isl].reshape(-1, 3)
            print(f"    steep overhang area {a:6.1f} mm2 z[{P[:, 2].min():.1f},{P[:, 2].max():.1f}] near ({P[:, 0].mean():.1f},{P[:, 1].mean():.1f})")


def layer_masks(T, layer=0.2, res=0.1):
    """Slice the mesh at every layer centre and fill each slice (even-odd scanlines) on a res grid."""
    lo = T.reshape(-1, 3).min(0)
    hi = T.reshape(-1, 3).max(0)
    nx, ny = int((hi[0] - lo[0]) / res) + 3, int((hi[1] - lo[1]) / res) + 3
    tzmin, tzmax = T[..., 2].min(1), T[..., 2].max(1)
    for z in np.arange(lo[2] + layer / 2, hi[2], layer):
        tri = T[(tzmin < z) & (tzmax > z)]
        d = tri[..., 2] - z
        ends = []
        for a, b in ((0, 1), (1, 2), (2, 0)):
            cross = d[:, a] * d[:, b] < 0
            t = d[:, a] / np.where(cross, d[:, a] - d[:, b], 1.0)
            ends.append((cross, tri[:, a, :2] + (tri[:, b, :2] - tri[:, a, :2]) * t[:, None]))
        segs = []
        for (c1, p1), (c2, p2) in ((ends[0], ends[1]), (ends[1], ends[2]), (ends[2], ends[0])):
            m = c1 & c2
            segs.append(np.stack([p1[m], p2[m]], 1))
        S = np.concatenate(segs) if segs else np.zeros((0, 2, 2))
        acc = np.zeros((ny, nx + 1), np.int32)
        if len(S):
            rows = lo[1] + (np.arange(ny) + 0.5) * res
            y1, y2 = S[:, 0, 1][:, None], S[:, 1, 1][:, None]
            hit = (y1 <= rows) != (y2 <= rows)                      # half-open rule
            si, ri = np.nonzero(hit)
            x1, x2 = S[si, 0, 0], S[si, 1, 0]
            t = (rows[ri] - S[si, 0, 1]) / (S[si, 1, 1] - S[si, 0, 1])
            x = x1 + (x2 - x1) * t
            ci = np.clip(np.ceil((x - lo[0]) / res - 0.5).astype(int), 0, nx)
            np.add.at(acc, (ri, ci), 1)
        yield z, (np.cumsum(acc, axis=1)[:, :nx] % 2) == 1, lo


def floating_islands(path, min_area=0.3):
    """Regions of a layer that touch nothing in the layer below: they would be printed in mid-air."""
    from scipy.ndimage import label
    T = load_stl(path)
    prev, found = None, []
    for z, mask, lo in layer_masks(T):
        if prev is not None:
            lab, n = label(mask, structure=np.ones((3, 3)))
            for k in range(1, n + 1):
                comp = lab == k
                if not (comp & prev).any():
                    a = comp.sum() * 0.01
                    if a >= min_area:
                        r, c = np.nonzero(comp)
                        found.append((z, a, lo[0] + c.mean() * 0.1, lo[1] + r.mean() * 0.1))
        prev = mask
    name = path.split("/")[-1]
    if found:
        for z, a, x, y in found:
            print(f"  {name}: FLOATING island at z={z:5.2f}  area {a:5.1f} mm2  near ({x:.1f},{y:.1f})")
    else:
        print(f"  {name}: no floating islands")
    return found


if __name__ == "__main__":
    for p in sys.argv[1:]:
        report(p)
    print("\nfloating islands (layer 0.2 mm, grid 0.1 mm):")
    for p in sys.argv[1:]:
        floating_islands(p)