"""One sealed, single-piece surface around a model, whatever state its triangles are in."""

import numpy as np
from scipy import ndimage
from skimage import measure

PAD = 3


def solid_grid(vertices, faces, resolution=256):
    """The model as filled voxels: (solid, voxel size, grid origin)."""
    lo, hi = vertices.min(0), vertices.max(0)
    h = (hi - lo).max() / resolution
    shape = np.ceil((hi - lo) / h).astype(int) + 2 * PAD + 1

    # The wall: sample every triangle at half a voxel so none is skipped.
    a, b, c = vertices[faces[:, 0]], vertices[faces[:, 1]], vertices[faces[:, 2]]
    longest = np.max(np.stack([np.linalg.norm(b - a, axis=1), np.linalg.norm(c - a, axis=1), np.linalg.norm(c - b, axis=1)]), axis=0)
    n = np.maximum(1, np.ceil(longest / (0.5 * h))).astype(int)
    wall = np.zeros(shape, dtype=bool)
    for steps in np.unique(n):
        sel = n == steps
        i, j = np.meshgrid(np.arange(steps + 1), np.arange(steps + 1), indexing="ij")
        keep = i + j <= steps
        wi, wj = (i[keep] / steps)[None, :, None], (j[keep] / steps)[None, :, None]
        pts = a[sel][:, None] + (b[sel] - a[sel])[:, None] * wi + (c[sel] - a[sel])[:, None] * wj
        idx = np.floor((pts.reshape(-1, 3) - lo) / h).astype(int) + PAD
        wall[idx[:, 0], idx[:, 1], idx[:, 2]] = True

    # Thicken, flood the outside from a corner, take the voxel back.
    thick = ndimage.binary_dilation(wall)
    labels, _ = ndimage.label(~thick)
    outside = labels == labels[0, 0, 0]
    solid = ~outside
    solid &= ~(ndimage.binary_dilation(outside) & ~wall)
    return solid, h, lo


def ball(radius):
    r = int(radius)
    x, y, z = np.mgrid[-r : r + 1, -r : r + 1, -r : r + 1]
    return x * x + y * y + z * z <= radius * radius + 0.5


def thicken(solid, radius, where=None):
    """Parts thinner than a ball of `radius` voxels grown to that ball's width.

    A remesher cannot put faces on both sides of a sheet thinner than its
    edges: the two sides snap together and leave holes and loose shards. A
    spoiler, a handle or a blade thinner than about one face is widened to one
    face here, which changes the look by less than a face does.

    Sharp outer edges and corners look thin to this test too, and would grow
    a bead; `where` (a mask of the grid) limits the widening to the places
    that need it.
    """
    if radius < 1:
        return solid
    structure = ball(radius)
    thin = solid & ~ndimage.binary_opening(solid, structure=structure)
    if where is not None:
        thin &= where
    if not thin.any():
        return solid
    return solid | ndimage.binary_dilation(thin, structure=structure)


def near(shape, h, lo, points, reach):
    """A mask of the grid cells within `reach` (model units) of any of `points`."""
    mask = np.zeros(shape, dtype=bool)
    idx = np.clip(np.floor((np.asarray(points) - lo) / h).astype(int) + PAD, 0, np.array(shape) - 1)
    mask[idx[:, 0], idx[:, 1], idx[:, 2]] = True
    return ndimage.distance_transform_edt(~mask) * h <= reach


def grid_surface(solid, h, lo):
    """The sealed surface of a voxel solid: (vertices, triangles)."""
    # A smooth field, so marching cubes gives a surface without staircase.
    field = ndimage.gaussian_filter(solid.astype(np.float32), sigma=1.0)
    verts, tris, _, _ = measure.marching_cubes(field, level=0.5, allow_degenerate=False)
    verts = (verts - PAD) * h + lo
    # Marching cubes winds faces toward the low side of the field, which is the
    # outside here; turn them so normals point out, and check by volume.
    tris = tris[:, [0, 2, 1]]
    a, b, c = verts[tris[:, 0]], verts[tris[:, 1]], verts[tris[:, 2]]
    if np.einsum("ij,ij->i", a, np.cross(b, c)).sum() < 0:
        tris = tris[:, [0, 2, 1]]
    return verts.astype(np.float32), tris.astype(np.uint32)


def solid_surface(vertices, faces, resolution=256, min_thickness=0.0):
    """Sealed surface around the model; parts thinner than `min_thickness` (model units) are widened to it."""
    solid, h, lo = solid_grid(vertices, faces, resolution)
    # Anything under two voxels would be blurred away by the smoothing: that is the floor.
    radius = max(1.0, 0.5 * min_thickness / h)
    solid = thicken(solid, radius)
    v, f = grid_surface(solid, h, lo)
    return v, f, h