"""Convex hull straight from an FCP descriptor: no decode grid, no mesh, no qhull.

The hull of the decoded solid is the hull of its solid rings (holes are inside; a filled section has the hull of its
outline), with each part's end rings held to the part's outer edges (`_parts.part_range`). Per level (a layer centre or
a part edge): the 2D hull of its ring samples and its support point in M directions. Per direction: the upper concave
envelope of (z, support) over the levels; a hull vertex whose normal has that azimuth lies on it. Neighbouring
directions' chains are joined by z and the end levels capped: a closed triangle surface, convex up to the M directions.
"""
import importlib

import numba
import numpy as np

D = importlib.import_module("fcp.decode")


@numba.njit(cache=True)
def _hull2d(x, y, idx, buf):
    """Andrew's monotone chain over the points idx; CCW hull indices into buf, returns their count."""
    o = idx[np.argsort(y[idx], kind="mergesort")]
    o = o[np.argsort(x[o], kind="mergesort")]
    n = len(o)
    k = 0
    for i in range(n):
        p = o[i]
        while k >= 2 and (x[buf[k-1]]-x[buf[k-2]])*(y[p]-y[buf[k-2]]) - (y[buf[k-1]]-y[buf[k-2]])*(x[p]-x[buf[k-2]]) <= 0:
            k -= 1
        buf[k] = p
        k += 1
    t = k + 1
    for i in range(n - 2, -1, -1):
        p = o[i]
        while k >= t and (x[buf[k-1]]-x[buf[k-2]])*(y[p]-y[buf[k-2]]) - (y[buf[k-1]]-y[buf[k-2]])*(x[p]-x[buf[k-2]]) <= 0:
            k -= 1
        buf[k] = p
        k += 1
    return k - 1 if n > 1 else n


@numba.njit(cache=True)
def _loft(x, y, start, zl, M):
    L = len(zl)
    ca = np.empty(M)
    sa = np.empty(M)
    for m in range(M):
        ca[m] = np.cos(2 * np.pi * m / M)
        sa[m] = np.sin(2 * np.pi * m / M)
    S = np.empty((L, M))
    arg = np.empty((L, M), np.int64)
    buf = np.empty(int((start[1:] - start[:-1]).max()) + 2, np.int64)
    for l in range(L):  # support points: a rotating pointer around the level's CCW hull
        k = _hull2d(x, y, np.arange(start[l], start[l + 1]), buf)
        j = 0
        bv = -np.inf
        for jj in range(k):
            v = x[buf[jj]] * ca[0] + y[buf[jj]] * sa[0]
            if v > bv:
                bv = v
                j = jj
        for m in range(M):
            v = x[buf[j]] * ca[m] + y[buf[j]] * sa[m]
            while True:
                jn = j + 1 if j + 1 < k else 0
                vn = x[buf[jn]] * ca[m] + y[buf[jn]] * sa[m]
                if vn > v:
                    j = jn
                    v = vn
                else:
                    break
            S[l, m] = v
            arg[l, m] = buf[j]
    chain = np.empty((M, L), np.int64)  # per direction: the levels on the upper envelope of (z, support)
    clen = np.zeros(M, np.int64)
    for m in range(M):
        k = 0
        for l in range(L):
            while k >= 2:
                a, c = chain[m, k - 2], chain[m, k - 1]
                if (zl[c] - zl[a]) * (S[l, m] - S[a, m]) - (S[c, m] - S[a, m]) * (zl[l] - zl[a]) >= 0:
                    k -= 1
                else:
                    break
            chain[m, k] = l
            k += 1
        clen[m] = k
    vid = -np.ones(len(x), np.int64)
    nv = 0
    for m in range(M):
        for i in range(clen[m]):
            s = arg[chain[m, i], m]
            if vid[s] < 0:
                vid[s] = nv
                nv += 1
    V = np.empty((nv, 3))
    for l in range(L):
        for s in range(start[l], start[l + 1]):
            if vid[s] >= 0:
                V[vid[s], 0] = x[s]
                V[vid[s], 1] = y[s]
                V[vid[s], 2] = zl[l]
    T = np.empty((2 * M * L + 2 * M, 3), np.int64)
    nt = 0
    for m in range(M):  # bands: neighbouring chains merged by z
        m2 = m + 1 if m + 1 < M else 0
        i = 0
        j = 0
        na = clen[m]
        nb = clen[m2]
        while i < na - 1 or j < nb - 1:
            a = vid[arg[chain[m, i], m]]
            b = vid[arg[chain[m2, j], m2]]
            if j == nb - 1 or (i < na - 1 and zl[chain[m, i + 1]] <= zl[chain[m2, j + 1]]):
                c = vid[arg[chain[m, i + 1], m]]
                i += 1
            else:
                c = vid[arg[chain[m2, j + 1], m2]]
                j += 1
            if a != b and b != c and a != c:
                T[nt, 0] = a
                T[nt, 1] = b
                T[nt, 2] = c
                nt += 1
    for e in range(2):  # caps: fans over the end levels' support polygons
        l = 0 if e == 0 else L - 1
        a = vid[arg[l, 0]]
        for m in range(1, M - 1):
            b = vid[arg[l, m]]
            c = vid[arg[l, m + 1]]
            if a != b and b != c and a != c:
                T[nt, 0] = a
                T[nt, 1] = c if e == 0 else b
                T[nt, 2] = b if e == 0 else c
                nt += 1
    return V, T[:nt]


def levels(d, n_t=128):
    """The solid rings as samples grouped by level (ascending z): x, y, start offsets, level heights (canonical)."""
    c = D.grid(d.layers)
    old = D._n_t
    D._n_t = lambda c: n_t
    try:
        r = D._axis_rings(d.parts, d.coef(), c, d.layers, np, interp=d.interp)
    finally:
        D._n_t = old
    h, s = r.h, r.sign > 0
    first = np.r_[True, h[1:] != h[:-1] + 1]  # rings come part by part, heights consecutive within a part
    last = np.r_[h[1:] != h[:-1] + 1, True]
    ds = (c[1] - c[0]) / 2
    sel = np.r_[s, s[first], s[last]]
    XY = np.concatenate([r.XY, r.XY[first], r.XY[last]])[sel]
    z = np.concatenate([c[h], c[h[first]] - ds, c[h[last]] + ds])[sel].astype(np.float64)
    o = np.argsort(z, kind="stable")
    XY, z = XY[o], z[o]
    zl, cnt = np.unique(z, return_counts=True)
    start = np.r_[0, np.cumsum(cnt * XY.shape[2])]
    return XY[:, 0].ravel().astype(np.float64), XY[:, 1].ravel().astype(np.float64), start, zl


def convex_hull(d, M=64, n_t=128):
    """World vertices (V, 3) and triangles (T, 3), outward (counter-clockwise seen from outside)."""
    x, y, start, zl = levels(d, n_t)
    V, T = _loft(x, y, start, zl, M)
    if np.linalg.det(d.pose.rotation) * np.prod(d.pose.scale) < 0:
        T = T[:, ::-1]
    return d.pose.to_world(V), T