"""Prototype (bench only, the package is untouched): one scale for everything.

The encoder works at one scale eps, the finest detail the input supports: the point spacing s for a clean cloud, C sigma
when the noise sigma is larger (eps = max(s, C sigma), sigma given: the noise is a parameter). With f = eps / s:
  * surface: the surf slicer's robust-plane inlier distance (TAU), RIMLS support (H_SUPPORT) and band (BAND) times f;
  * small loops: the minimum loop area (MIN_LOOP spacings^2) times f^2;
  * parts (cuts="pointtest"): a track is cut only where one height fit fails to explain the points (pointfit.py, its
    margin of 2 spacings; scaling it by f was measured and changes little);
  * thin walls (thin=True, off by default): thinwall.py at the spacing scale (helps the clean loaf pan, costs a little
    under noise; scaling it by f deletes geometry).
For a clean cloud f = 1: the package's slicer with the point test.
Measured (report/eps_components.txt, eps_variants.txt, eps_minloop.txt): under per-point noise 5-30x fewer parts than
the package at about the noise level of error; C = 4; the minimum loop as eps^2 beats eps and a fixed one.
estimate_sigma (a local quadric fit) reads thick shapes' noise but takes thin walls for noise (report/sigma_estimators.txt).

    from eps import encode_eps
    d = encode_eps(points, up=(0, 0, 1), sigma=0.006 * diag)   # the noise, world units (0: clean)
"""
import sys

import numpy as np
from scipy.spatial import cKDTree

sys.path.insert(0, "/home/joao/code/fcp_devel/api")
sys.path.insert(0, __file__.rsplit("/", 1)[0])
import fcp
from fcp import _slice, _surf
from fcp._slice import PointSlicer

C = 4.0  # eps = max(spacing, C sigma)
K_SIGMA = 64  # neighbours of the plane fits of estimate_sigma
LAST = {}


def estimate_sigma(P, k=24, n_query=3000, seed=0):
    """The noise sigma (world units): the median RMS residual of a local quadratic height fit over k neighbours (degrees
    of freedom corrected). Reads thick shapes' noise; a thin wall's two faces read as noise (no auto mode yet)."""
    P = np.asarray(P, float)
    q = np.random.default_rng(seed).choice(len(P), min(n_query, len(P)), replace=False)
    _, nb = cKDTree(P).query(P[q], k)
    out = []
    for ids in nb:
        Q = P[ids] - P[ids].mean(0)
        _, v = np.linalg.eigh(Q.T @ Q)
        L = Q @ v
        x, y, z = L[:, 2], L[:, 1], L[:, 0]
        A = np.c_[np.ones(k), x, y, x * x, x * y, y * y]
        c, *_ = np.linalg.lstsq(A, z, rcond=None)
        r = z - A @ c
        out.append(np.sqrt(r @ r / (k - 6)))
    return float(np.median(out))


def encode_eps(P, up=None, sigma=0.0, thin=False, cuts="pointtest", device="cuda", **kw):
    P = np.asarray(P, float)
    if sigma == "auto":
        sigma = estimate_sigma(P)
    s = float(PointSlicer(P, up=up, device=device).spacing)
    f = max(1.0, C * float(sigma) / s)
    LAST.update(sigma=float(sigma), spacing=s, f=f)
    saved_surf = {k: getattr(_surf, k) for k in ("TAU", "H_SUPPORT", "BAND")}
    saved_loop = _slice.MIN_LOOP
    try:
        for k, v in saved_surf.items():
            setattr(_surf, k, v * f)
        _slice.MIN_LOOP = saved_loop * f * f
        if thin:
            import thinwall as tw

            return tw.encode_thinwall(P, up=up, device=device, cuts=cuts, **kw)
        if cuts == "pointtest":
            import pointfit

            return pointfit.encode_pointtest(P, up=up, device=device, **kw)
        return fcp.encode(P, up=up, device=device, **kw)
    finally:
        for k, v in saved_surf.items():
            setattr(_surf, k, v)
        _slice.MIN_LOOP = saved_loop