"""Prototype (bench only, the package is untouched): the GWN slicer.

Inside / outside from one global field instead of a per-layer flood: the mollified generalized winding number of the
oriented, area-weighted points (Barill et al. 2018; the iso-value of DWG, Liu et al. 2025; the orientation objective
of GCNO, Xu et al. 2023 / Gotsman & Hormann 2024, here on patches):

  w(x) = 1/(4 pi) sum_i a_i n_i . (p_i - x) / (|p_i - x|^2 + eps_i^2)^(3/2)

  a_i     each point's area (pi r_16^2 / 15, as the package's pose weights): the field does not see the density;
  eps_i   the mollifier: max(C_EPS * local spacing_i, C_SIGMA * sigma) -- the one scale knob. Small: the field follows
          the points (a flat sheet's 0.5 crossing is exactly on the sheet for every eps); large: smooth, no foam,
          detail under eps is lost;
  n_i     surf's oriented directions (`_surf.Surf`), then patch flips that make w two-valued (0 outside, 1 inside) on
          probes off the surface (double-well energy w^2 (1 - w)^2), then per-point flips by the jump of w across each
          point (its own term left out);
  tau     the iso-value: the median of w at the points (DWG), so a global area bias does not move the contour.

Per layer the raster of tau - w (positive outside) goes through surf's own marching squares and loop chaining; loops
under MIN_LOOP eps^2 go. The rest of the encoder is the package's.

    from gwn import encode_gwn
    d = encode_gwn(points, up=(0, 0, 1), c_eps=1.0, sigma=0.0)
"""
import sys
import time

import numpy as np

sys.path.insert(0, "/home/joao/code/fcp_devel/api")
import cupy as cp  # noqa: E402
import fcp  # noqa: E402
from fcp import _slice, _surf  # noqa: E402
from fcp._rings import host  # noqa: E402
from fcp._slice import PointSlicer  # noqa: E402

_encode = sys.modules["fcp.encode"]

C_EPS = 0.5  # eps_i = C_EPS x local spacing (clean)
C_SIGMA = 3.0  # ... at least C_SIGMA x sigma (sigma given, world units)
PIX = 0.5  # raster pixel, spacings (surf's)
MIN_LOOP = 2.0  # loops under MIN_LOOP x max(spacing, eps_med)^2 go (the package's 2 spacings^2 when eps <= spacing)
PROBE = 1.5  # probes at PROBE eps_i off each point, both sides
N_PROBE = 6000  # points probed for the orientation (both sides), and as many empty-space probes
MIN_PATCH = 30  # patches with fewer points are not flipped as a whole (the per-point pass decides them)
MAX_PATCHES = 400  # largest patches that are variables of the patch flips
POINT_ROUNDS = 0  # per-point flip rounds (measured: they flip thin walls; off)
LOG = {}
ISO = "local"  # "global": tau = the median of w at the points (DWG); "local": tau(x) = Shepard average of w at the points
TAU_MIN = 0.05  # the iso-value of an indicator lies in (0, 1): clamped to [TAU_MIN, 1 - TAU_MIN]
SHEP_POW = 4.0  # Shepard weights a / (d^2 + eps^2)^(SHEP_POW / 2)
SHEP_EPS = 1.0  # ... with eps x SHEP_EPS
THIN = 0  # 1: thin sheets (points whose neighbours' normals cancel) also count as inside within H_THIN eps of the
#           points (a band from the normalized density; the unsigned "metric" part beside the winding "phase")
THIN_BY = "contrast"  # "contrast": |w(p - d n) - w(p + d n)| (d = PROBE eps) below CT_LO: thin, above CT_HI: not
#                      (the field cannot tell the two sides apart there, whatever the normals' signs); "normals": below
CT_LO, CT_HI = 0.15, 0.35
CR_LO, CR_HI = 0.3, 0.6  # |mean normal| of the neighbours (weights a / (d^2 + eps^2)^4): thin below CR_LO, not above CR_HI
H_THIN = 1.0  # band half-thickness of a thin sheet, eps
THIN_EPS = 1.0  # density / thinness kernels: eps x THIN_EPS
SPLAT = 1  # 1: each point's dipole spread over a disc in its tangent plane (its centre and a hexagon of SPLAT_R local
#            spacings, area a / 7 each) with the small mollifier SPLAT_EPS local spacings: wide along the surface (the gaps
#            between samples close), sharp across it (walls of about a spacing stay apart)
SPLAT_R = 1.25
SPLAT_EPS = 0.25
CLOSE = 0.0  # > 0: each layer's inside region is closed (dilated, then eroded) by CLOSE spacings; the pixels the
#             closing adds (gaps narrower than 2 CLOSE spacings: a thin wall opened by a sampling gap) become just inside
CLOSE3 = 0.0  # > 0: the same closing in 3D (all layers on one grid, the layer step as the third spacing), CLOSE3 spacings
CLOSE_IN = 1e-3  # their field value (inside, small: the contour stays at the pixel's edge)
NORMALS = "surf"  # "surf": surf's directions (then the patch flips); "grad": -grad w at the points after them (DWG)

_SRC = r"""
extern "C" __global__ void gwn(const float* X, const float* V, const float* E2, const int m,
                               const float* Q, const long long q, float* out) {
    // w at each query: the mollified dipole sum over all points; points in shared-memory tiles
    __shared__ float sx[256], sy[256], sz[256], vx[256], vy[256], vz[256], se[256];
    long long t = blockDim.x * (long long)blockIdx.x + threadIdx.x;
    float qx = 0.f, qy = 0.f, qz = 0.f;
    if (t < q) { qx = Q[3 * t]; qy = Q[3 * t + 1]; qz = Q[3 * t + 2]; }
    double acc = 0.0;
    for (int base = 0; base < m; base += 256) {
        int k = base + threadIdx.x;
        if (k < m) {
            sx[threadIdx.x] = X[3 * k]; sy[threadIdx.x] = X[3 * k + 1]; sz[threadIdx.x] = X[3 * k + 2];
            vx[threadIdx.x] = V[3 * k]; vy[threadIdx.x] = V[3 * k + 1]; vz[threadIdx.x] = V[3 * k + 2];
            se[threadIdx.x] = E2[k];
        }
        __syncthreads();
        int n = min(256, m - base);
        float a = 0.f;
        for (int j = 0; j < n; ++j) {
            float dx = sx[j] - qx, dy = sy[j] - qy, dz = sz[j] - qz;
            float r2 = dx * dx + dy * dy + dz * dz + se[j];
            float inv = rsqrtf(r2);
            a += (vx[j] * dx + vy[j] * dy + vz[j] * dz) * inv * inv * inv;
        }
        acc += a;
        __syncthreads();
    }
    if (t < q) out[t] = (float)(acc * 0.07957747154594767);  // 1 / (4 pi)
}
"""
_K = {}


def _kernel():
    if "gwn" not in _K:
        _K["gwn"] = cp.RawKernel(_SRC, "gwn")
    return _K["gwn"]


def wind(X, V, E2, Q, chunk=1 << 20):
    """w at the queries Q (q, 3) of the points X with dipoles V = a n and squared mollifiers E2 (CuPy arrays)."""
    Q = cp.ascontiguousarray(cp.asarray(Q, cp.float32))
    out = cp.empty(len(Q), cp.float32)
    for a in range(0, len(Q), chunk):
        b = min(len(Q), a + chunk)
        sub = Q[a:b]
        _kernel()(((b - a + 255) // 256,), (256,), (X, V, E2, np.int32(len(X)), sub, np.int64(b - a), out[a:b]))
    return out

_SRC_G = r"""
extern "C" __global__ void gwn_grad(const float* X, const float* V, const float* E2, const int m,
                                    const float* Q, const long long q, float* out) {
    // w and grad w at each query (out: 4 per query)
    __shared__ float sx[256], sy[256], sz[256], vx[256], vy[256], vz[256], se[256];
    long long t = blockDim.x * (long long)blockIdx.x + threadIdx.x;
    float qx = 0.f, qy = 0.f, qz = 0.f;
    if (t < q) { qx = Q[3 * t]; qy = Q[3 * t + 1]; qz = Q[3 * t + 2]; }
    double w = 0.0, gx = 0.0, gy = 0.0, gz = 0.0;
    for (int base = 0; base < m; base += 256) {
        int k = base + threadIdx.x;
        if (k < m) {
            sx[threadIdx.x] = X[3 * k]; sy[threadIdx.x] = X[3 * k + 1]; sz[threadIdx.x] = X[3 * k + 2];
            vx[threadIdx.x] = V[3 * k]; vy[threadIdx.x] = V[3 * k + 1]; vz[threadIdx.x] = V[3 * k + 2];
            se[threadIdx.x] = E2[k];
        }
        __syncthreads();
        int n = min(256, m - base);
        float a = 0.f, bx = 0.f, by = 0.f, bz = 0.f;
        for (int j = 0; j < n; ++j) {
            float dx = sx[j] - qx, dy = sy[j] - qy, dz = sz[j] - qz;
            float r2 = dx * dx + dy * dy + dz * dz + se[j];
            float inv = rsqrtf(r2), i3 = inv * inv * inv, i5 = i3 * inv * inv;
            float vd = vx[j] * dx + vy[j] * dy + vz[j] * dz;
            a += vd * i3;
            bx += -vx[j] * i3 + 3.f * vd * dx * i5;
            by += -vy[j] * i3 + 3.f * vd * dy * i5;
            bz += -vz[j] * i3 + 3.f * vd * dz * i5;
        }
        w += a; gx += bx; gy += by; gz += bz;
        __syncthreads();
    }
    if (t < q) {
        const double c = 0.07957747154594767;
        out[4 * t] = (float)(w * c); out[4 * t + 1] = (float)(gx * c); out[4 * t + 2] = (float)(gy * c);
        out[4 * t + 3] = (float)(gz * c);
    }
}
"""


def wind_grad(X, V, E2, Q, chunk=1 << 20):
    """(w, grad w) at the queries Q."""
    if "gwn_grad" not in _K:
        _K["gwn_grad"] = cp.RawKernel(_SRC_G, "gwn_grad")
    Q = cp.ascontiguousarray(cp.asarray(Q, cp.float32))
    out = cp.empty((len(Q), 4), cp.float32)
    for a in range(0, len(Q), chunk):
        b = min(len(Q), a + chunk)
        _K["gwn_grad"](((b - a + 255) // 256,), (256,), (X, V, E2, np.int32(len(X)), Q[a:b], np.int64(b - a), out[a:b]))
    return out[:, 0], out[:, 1:]

_SRC_T = r"""
extern "C" __global__ void shepard(const float* X, const float* A, const float* E2, const float* Wv, const int m,
                                   const float* Q, const long long q, const float pw, float* out) {
    // the local iso-value: sum g w / sum g with g = a / (|p - x|^2 + e^2)^(pw / 2)
    __shared__ float sx[256], sy[256], sz[256], sa[256], se[256], sw[256];
    long long t = blockDim.x * (long long)blockIdx.x + threadIdx.x;
    float qx = 0.f, qy = 0.f, qz = 0.f;
    if (t < q) { qx = Q[3 * t]; qy = Q[3 * t + 1]; qz = Q[3 * t + 2]; }
    double num = 0.0, den = 0.0;
    for (int base = 0; base < m; base += 256) {
        int k = base + threadIdx.x;
        if (k < m) {
            sx[threadIdx.x] = X[3 * k]; sy[threadIdx.x] = X[3 * k + 1]; sz[threadIdx.x] = X[3 * k + 2];
            sa[threadIdx.x] = A[k]; se[threadIdx.x] = E2[k]; sw[threadIdx.x] = Wv[k];
        }
        __syncthreads();
        int n = min(256, m - base);
        float a = 0.f, b = 0.f;
        for (int j = 0; j < n; ++j) {
            float dx = sx[j] - qx, dy = sy[j] - qy, dz = sz[j] - qz;
            float g = sa[j] * powf(dx * dx + dy * dy + dz * dz + se[j], -0.5f * pw);
            a += g * sw[j]; b += g;
        }
        num += a; den += b;
        __syncthreads();
    }
    if (t < q) out[t] = (float)(num / den);
}
"""


def shepard(X, A, E2, Wv, Q, pw=4.0, chunk=1 << 20):
    """The Shepard average of the point values Wv at the queries Q (weights a / (d^2 + e^2)^(pw/2))."""
    if "shepard" not in _K:
        _K["shepard"] = cp.RawKernel(_SRC_T, "shepard")
    Q = cp.ascontiguousarray(cp.asarray(Q, cp.float32))
    out = cp.empty(len(Q), cp.float32)
    for a in range(0, len(Q), chunk):
        b = min(len(Q), a + chunk)
        _K["shepard"](((b - a + 255) // 256,), (256,), (X, A, E2, Wv, np.int32(len(X)), Q[a:b], np.int64(b - a),
                                                       np.float32(pw), out[a:b]))
    return out

_SRC_D = r"""
extern "C" __global__ void density(const float* X, const float* A, const float* E2, const int m,
                                   const float* Q, const long long q, float* out) {
    // rho(x) = sum a / (2 pi e^2) exp(-|p - x|^2 / (2 e^2)): 1 on a flat sheet, falls to exp(-z^2/2e^2) off it
    __shared__ float sx[256], sy[256], sz[256], sa[256], se[256];
    long long t = blockDim.x * (long long)blockIdx.x + threadIdx.x;
    float qx = 0.f, qy = 0.f, qz = 0.f;
    if (t < q) { qx = Q[3 * t]; qy = Q[3 * t + 1]; qz = Q[3 * t + 2]; }
    double acc = 0.0;
    for (int base = 0; base < m; base += 256) {
        int k = base + threadIdx.x;
        if (k < m) {
            sx[threadIdx.x] = X[3 * k]; sy[threadIdx.x] = X[3 * k + 1]; sz[threadIdx.x] = X[3 * k + 2];
            sa[threadIdx.x] = A[k]; se[threadIdx.x] = E2[k];
        }
        __syncthreads();
        int n = min(256, m - base);
        float a = 0.f;
        for (int j = 0; j < n; ++j) {
            float dx = sx[j] - qx, dy = sy[j] - qy, dz = sz[j] - qz;
            float d2 = dx * dx + dy * dy + dz * dz;
            if (d2 < 25.f * se[j]) a += sa[j] / (6.283185307f * se[j]) * __expf(-0.5f * d2 / se[j]);
        }
        acc += a;
        __syncthreads();
    }
    if (t < q) out[t] = (float)acc;
}
"""


def density(X, A, E2, Q, chunk=1 << 20):
    """The normalized point density at the queries (1 on a flat sheet sampled at the areas A)."""
    if "density" not in _K:
        _K["density"] = cp.RawKernel(_SRC_D, "density")
    Q = cp.ascontiguousarray(cp.asarray(Q, cp.float32))
    out = cp.empty(len(Q), cp.float32)
    for a in range(0, len(Q), chunk):
        b = min(len(Q), a + chunk)
        _K["density"](((b - a + 255) // 256,), (256,), (X, A, E2, np.int32(len(X)), Q[a:b], np.int64(b - a), out[a:b]))
    return out


def _self_term(V, E2, delta, n_sign):
    """One point's own contribution at p + n_sign delta n (its dipole along n): -n_sign a delta / (delta^2+e^2)^1.5 / 4pi."""
    a = cp.sqrt((V * V).sum(1))
    return -n_sign * a * delta / (delta * delta + E2) ** 1.5 / (4 * np.pi)


def _patches(X, N, idx):
    """Connected patches of the k-NN graph along one consistently oriented sheet (host labels)."""
    from scipy.sparse import coo_matrix
    from scipy.sparse.csgraph import connected_components

    m, k = idx.shape
    j = idx[:, 1:]
    i = cp.repeat(cp.arange(m)[:, None], k - 1, 1)
    d = X[j] - X[i][:, :]
    dn = cp.sqrt((d * d).sum(2)) + 1e-30
    ok = ((N[i] * N[j]).sum(2) > _slice.SMOOTH)
    ok &= (cp.abs((d * N[i]).sum(2)) / dn < _slice.TANGENT) & (cp.abs((d * N[j]).sum(2)) / dn < _slice.TANGENT)
    a, b = host(i[ok]), host(j[ok])
    G = coo_matrix((np.ones(len(a)), (a, b)), shape=(m, m))
    return connected_components(G, directed=False)[1]


def _double_well(w):
    return (w * w) * (1 - w) * (1 - w)


class GWN:
    """The oriented points of one cloud (surf's metric canonical frame), their mollified winding field, loops()."""

    def __init__(self, ps: PointSlicer, c_eps=C_EPS, sigma=0.0, orient=True, seed=0, point_rounds=None):
        S = ps._surf
        X, N, idx = S.X, S.N.copy(), S.idx
        s = float(ps.spacing)
        m, k = idx.shape
        t0 = time.perf_counter()
        # local area and spacing from the k-th neighbour (the package's pose weights)
        dk = cp.sqrt(((X[idx[:, -1]] - X) ** 2).sum(1))
        area = np.pi * dk * dk / (k - 1)
        loc = cp.sqrt(area)
        eps = cp.maximum(c_eps * loc, cp.float32(C_SIGMA * sigma))
        self.X, self.s, self.eps = X, s, eps
        self.E2 = cp.ascontiguousarray((eps * eps).astype(cp.float32))
        self.area = area.astype(cp.float32)
        self.eps_med = float(cp.median(eps))
        rng = np.random.default_rng(seed)
        LOG.clear()
        LOG.update(spacing=s, eps_med=self.eps_med, eps_ratio=self.eps_med / s)
        self.point_rounds = POINT_ROUNDS if point_rounds is None else point_rounds
        if orient:
            N = self._orient(X, N, idx, rng)
        self.V = cp.ascontiguousarray((N * self.area[:, None]).astype(cp.float32))
        if NORMALS == "grad":
            _, g = wind_grad(X, self.V, self.E2, X)
            gn = cp.sqrt((g * g).sum(1, keepdims=True))
            N = cp.where(gn > 0, -g / cp.maximum(gn, 1e-30), N).astype(cp.float32)
            self.V = cp.ascontiguousarray((N * self.area[:, None]).astype(cp.float32))
        self.N = N
        self.FX, self.FV, self.FE2 = X, self.V, self.E2  # the field's dipoles (the splats when SPLAT)
        if SPLAT:
            # tangent basis per point
            a = cp.where(cp.abs(N[:, 0:1]) < 0.9, cp.asarray([[1, 0, 0]], cp.float32), cp.asarray([[0, 1, 0]], cp.float32))
            t1 = cp.cross(N, a); t1 /= cp.sqrt((t1 * t1).sum(1, keepdims=True))
            t2 = cp.cross(N, t1)
            R = (SPLAT_R * loc)[:, None]
            offs = [cp.zeros_like(X)] + [R * (np.cos(k * np.pi / 3) * t1 + np.sin(k * np.pi / 3) * t2) for k in range(6)]
            self.FX = cp.ascontiguousarray(cp.concatenate([X + o for o in offs]).astype(cp.float32))
            self.FV = cp.ascontiguousarray(cp.concatenate([self.V / 7] * 7).astype(cp.float32))
            e2 = (cp.maximum(SPLAT_EPS * loc, cp.float32(C_SIGMA * sigma)) ** 2).astype(cp.float32)
            self.FE2 = cp.ascontiguousarray(cp.concatenate([e2] * 7))
        wp = wind(self.FX, self.FV, self.FE2, X)
        self.wp = cp.ascontiguousarray(wp.astype(cp.float32))
        self.SE2 = cp.ascontiguousarray((self.E2 * SHEP_EPS * SHEP_EPS).astype(cp.float32))
        if THIN:
            self.TE2 = cp.ascontiguousarray((self.E2 * THIN_EPS * THIN_EPS).astype(cp.float32))
            if THIN_BY == "contrast":
                d = (PROBE * self.eps)[:, None]
                ct = cp.abs(wind(X, self.V, self.E2, X - d * N) - wind(X, self.V, self.E2, X + d * N))
                th = cp.clip((CT_HI - ct) / (CT_HI - CT_LO), 0, 1)
            else:
                nav = cp.stack([shepard(X, self.area, self.TE2, cp.ascontiguousarray(N[:, a]), X, 8.0) for a in range(3)], 1)
                cr = cp.sqrt((nav * nav).sum(1))
                th = cp.clip((CR_HI - cr) / (CR_HI - CR_LO), 0, 1)
            self.thin = cp.ascontiguousarray(th.astype(cp.float32))
            LOG.update(thin_share=float((self.thin > 0.5).mean()))
        self.tau = float(cp.median(wp))
        LOG.update(tau=self.tau, w_p10=float(cp.percentile(wp, 10)), w_p90=float(cp.percentile(wp, 90)),
                   t_init=time.perf_counter() - t0)

    # ---------------------------------------------------------------- orientation
    def _probes(self, X, N, rng):
        m = len(X)
        sel = cp.asarray(rng.choice(m, min(m, N_PROBE), replace=False))
        d = (PROBE * self.eps[sel])[:, None]
        inner, outer = X[sel] - d * N[sel], X[sel] + d * N[sel]
        # empty space: uniform in the padded box, farther than 3 eps_med from every point
        lo, hi = host(X.min(0)), host(X.max(0))
        pad = 0.1 * (hi - lo) + 3 * self.eps_med
        U = rng.uniform(lo - pad, hi + pad, (4 * N_PROBE, 3)).astype(np.float32)
        from scipy.spatial import cKDTree

        dist = cKDTree(host(X)).query(U)[0]
        U = U[dist > 3 * self.eps_med][:N_PROBE]
        out_box = ((U < lo) | (U > hi)).any(1)  # beyond the cloud's box: outside for certain
        return sel, inner, outer, cp.asarray(U), cp.asarray(out_box)

    def _orient(self, X, N, idx, rng):
        t0 = time.perf_counter()
        lab = _patches(X, N, idx)
        n_patch = int(lab.max()) + 1
        size = np.bincount(lab, minlength=n_patch)
        big = np.argsort(-size)[:MAX_PATCHES]
        big = big[size[big] >= MIN_PATCH]
        sel, inner, outer, U, out_box = self._probes(X, N, rng)
        Q = cp.concatenate([inner, outer, U])
        V = (N * self.area[:, None]).astype(cp.float32)
        lab_d = cp.asarray(lab)
        # columns: the field of each variable patch; the rest (small patches, fixed) in one column
        cols = []
        var = np.zeros(n_patch, bool)
        var[big] = True
        for p in big:
            mask = lab_d == int(p)
            cols.append(wind(cp.ascontiguousarray(X[mask]), cp.ascontiguousarray(V[mask]),
                             cp.ascontiguousarray(self.E2[mask]), Q))
        rest = ~cp.asarray(var)[lab_d]
        fixed = wind(cp.ascontiguousarray(X[rest]), cp.ascontiguousarray(V[rest]), cp.ascontiguousarray(self.E2[rest]), Q) \
            if bool(rest.any()) else cp.zeros(len(Q), cp.float32)
        W = cp.stack(cols, 1) if cols else cp.zeros((len(Q), 0), cp.float32)
        nq = len(sel)
        far = cp.zeros(len(Q), bool)
        far[2 * nq:] = out_box  # these want w = 0 (w^2), the rest the double well

        def energy(v):  # per probe (columns of v: candidate fields)
            dw = _double_well(v)
            return cp.where(far[:, None] if v.ndim == 2 else far, v * v, dw)
        # weights: probes beside the points want inner 1 / outer 0 (the double well lets each choose), empty space too
        sgn = cp.ones(W.shape[1], cp.float32)
        w = fixed + W @ sgn
        E0 = float(energy(w).sum())
        flips = 0
        for _ in range(4 * len(big) + 10):
            trial = w[:, None] - 2 * W * sgn[None, :]
            dE = energy(trial).sum(0) - energy(w).sum()
            j = int(cp.argmin(dE))
            if len(big) == 0 or float(dE[j]) >= -1e-3 * max(1.0, E0 / len(Q)):
                break
            w = w - 2 * W[:, j] * sgn[j]
            sgn[j] = -sgn[j]
            flips += 1
        E1 = float(energy(w).sum())
        mult = cp.ones(n_patch, cp.float32)
        mult[cp.asarray(big)] = sgn
        N = N * mult[lab_d][:, None]
        # all inverted? the jump of w across the points says which side is inside
        V = cp.ascontiguousarray((N * self.area[:, None]).astype(cp.float32))
        d = (PROBE * self.eps)[:, None]
        jump = wind(X, V, self.E2, X - d * N) - wind(X, V, self.E2, X + d * N)
        gflip = float(cp.median(jump)) < 0
        if gflip:
            N = -N
        # per point: the jump of w across the point (its own term out) decides its sign
        pflips = []
        for _ in range(self.point_rounds):
            V = cp.ascontiguousarray((N * self.area[:, None]).astype(cp.float32))
            d = (PROBE * self.eps)[:, None]
            wi = wind(X, V, self.E2, X - d * N) - _self_term(V, self.E2, d[:, 0], -1)
            wo = wind(X, V, self.E2, X + d * N) - _self_term(V, self.E2, d[:, 0], +1)
            bad = (wo - wi) > 0  # more inside on the outer side: flip
            nb = int(bad.sum())
            pflips.append(nb)
            if nb == 0:
                break
            N = cp.where(bad[:, None], -N, N)
        LOG.update(global_flip=int(gflip), patches=n_patch, patch_vars=len(big), patch_flips=flips, E0=E0, E1=E1, point_flips=pflips,
                   t_orient=time.perf_counter() - t0, probes=int(len(Q)), probe_inner_outer=nq)
        return N

    # ---------------------------------------------------------------- slicing
    def loops(self, heights_m):
        """Per metric height, the closed loops (metric plane coordinates, float64) of the tau level of w."""
        B = _surf.backend("cuda")
        _, _, chain = _surf._host_loops()
        F = np.float32
        s, X = self.s, self.X
        pix = PIX * s
        reach = 3 * self.eps_med + s  # points beyond this from the plane put no crossing on it
        hm = np.asarray(heights_m, np.float64)
        nL = len(hm)
        o = cp.argsort(X[:, 2])
        Xz = cp.ascontiguousarray(X[o])
        z = cp.ascontiguousarray(Xz[:, 2])
        a = cp.searchsorted(z, cp.asarray(hm - reach, F)).astype(cp.int64)
        b = cp.searchsorted(z, cp.asarray(hm + reach, F), side="right").astype(cp.int64)
        box = host(B.slab_boxes(Xz, a, b)).astype(np.float64).reshape(nL, 4)
        live = host(b - a) >= 3
        pad = 3 * self.eps_med + 2 * pix
        lo = box[:, :2] - pad
        HW = np.zeros((nL, 2), np.int64)
        HW[live] = np.ceil((box[live, 2:] + pad - lo[live]) / pix).astype(np.int64) + 1
        off = np.concatenate([[0], np.cumsum(HW[:, 0] * HW[:, 1])])
        coff = np.concatenate([[0], np.cumsum(np.where(live, (HW[:, 0] + 1) * (HW[:, 1] + 1), 0))])
        out = [[] for _ in range(nL)]
        total, ctot = int(off[-1]), int(coff[-1])
        if total == 0:
            return out
        # pixel centres, all layers
        Ls = np.repeat(np.arange(nL), HW[:, 0] * HW[:, 1])
        loc = np.arange(total) - off[Ls]
        Wd = HW[Ls, 1]
        i, j = loc // Wd, loc % Wd
        Qp = np.column_stack([lo[Ls, 0] + (i + 0.5) * pix, lo[Ls, 1] + (j + 0.5) * pix, hm[Ls]]).astype(F)
        t0 = time.perf_counter()
        Qd = cp.asarray(Qp)
        w = wind(self.FX, self.FV, self.FE2, Qd)
        if ISO == "local":
            fld = (cp.clip(shepard(X, self.area, self.SE2, self.wp, Qd, SHEP_POW), TAU_MIN, 1 - TAU_MIN) - w).astype(F)
        else:
            fld = (F(min(max(self.tau, TAU_MIN), 1 - TAU_MIN)) - w).astype(F)
        if THIN and float(self.thin.max()) > 0:
            beta = shepard(X, self.area, self.TE2, self.thin, Qd, 8.0)
            rho = density(X, self.area, self.TE2, Qd)
            band = beta * (rho - F(np.exp(-0.5 * H_THIN * H_THIN)))
            fld = cp.minimum(fld, -band).astype(F)
        if CLOSE > 0:
            from cupyx.scipy.ndimage import distance_transform_edt as edt

            rp = CLOSE * s / pix  # pixels
            added = 0
            for L in np.nonzero(live)[0]:
                H, W = int(HW[L, 0]), int(HW[L, 1])
                Fl = fld[off[L]:off[L + 1]].reshape(H, W)
                inside = Fl < 0
                if not bool(inside.any()):
                    continue
                dil = edt(~inside) <= rp
                closed = edt(dil) > rp
                gap = closed & ~inside
                Fl[gap] = F(-CLOSE_IN)
                added += int(gap.sum())
            LOG["close_px"] = added
        if CLOSE3 > 0 and live.sum() > 1:
            from cupyx.scipy.ndimage import distance_transform_edt as edt

            idx = np.nonzero(live)[0]
            g0 = lo[idx].min(0)
            hi = (lo[idx] + HW[idx] * pix).max(0)
            GH, GW = (np.ceil((hi - g0) / pix).astype(int) + 1)
            Lmin, Lmax = idx.min(), idx.max()
            nz = Lmax - Lmin + 1
            vol = cp.ones((nz, int(GH), int(GW)), bool)  # outside
            for L in idx:
                H, W = int(HW[L, 0]), int(HW[L, 1])
                i0, j0 = (np.round((lo[L] - g0) / pix)).astype(int)
                vol[L - Lmin, i0:i0 + H, j0:j0 + W] = fld[off[L]:off[L + 1]].reshape(H, W) >= 0
            dz = float(np.median(np.diff(hm))) if nL > 1 else pix
            r = CLOSE3 * s
            samp = (dz, pix, pix)
            dil = edt(vol, sampling=samp) <= r  # within r of the inside
            closed = edt(dil, sampling=samp) > r
            gap3 = closed & vol
            added = 0
            for L in idx:
                H, W = int(HW[L, 0]), int(HW[L, 1])
                i0, j0 = (np.round((lo[L] - g0) / pix)).astype(int)
                gp = gap3[L - Lmin, i0:i0 + H, j0:j0 + W]
                Fl = fld[off[L]:off[L + 1]].reshape(H, W)
                Fl[gp] = F(-CLOSE_IN)
                added += int(gp.sum())
            LOG["close3_px"] = added
        LOG["t_field"] = time.perf_counter() - t0
        LOG["pixels"] = total
        seg = B.segments(ctot, nL, cp.asarray(coff), cp.asarray(HW.astype(np.int32)), cp.asarray(off), fld, F(1.0))
        if seg is None:
            return out
        layer, eid, pt = (host(v) for v in seg)
        pt = pt.astype(np.float64)
        bounds = np.searchsorted(layer, np.arange(nL + 1))
        min_area = MIN_LOOP * max(self.s, self.eps_med) ** 2
        for L in range(nL):
            u, v = bounds[L], bounds[L + 1]
            if v - u < 2:
                continue
            loff, P, area = chain(eid[u:v], pt[u:v])
            keep = (area * pix * pix >= min_area) & (np.diff(loff) >= 3)
            if keep.any():
                Qd = lo[L] + (P - 0.5) * pix
                out[L] += [Qd[loff[c]:loff[c + 1]] for c in np.nonzero(keep)[0]]
        return out


def slicer_gwn(P, up=None, **kw):
    """A PointSlicer (surf: pose, spacing, oriented directions) whose loops come from the GWN field."""
    ps = PointSlicer(np.asarray(P, float), up=up, device="cuda", slicer="surf")
    g = GWN(ps, **kw)
    sc = ps.pose.scale.astype(np.float64)

    def loops(heights):
        L = g.loops(np.asarray(heights, np.float64) * sc[2])
        return [[(Q / sc[:2]).astype(np.float32) for Q in ls] for ls in L]

    ps.loops = loops
    ps.gwn = g
    return ps


FIT_K = 0.0  # > 0: the cut tolerance follows max(spacing, FIT_K sigma) (fit_tol "auto" reads the slicer's spacing)


def encode_gwn(P, up=None, cuts=None, layers=128, c_eps=C_EPS, sigma=0.0, orient="auto", point_rounds=None, **kw):
    """fcp.encode of the points P with the GWN slicer (cuts="pointtest": the point test of ../fcp_2710_tests)."""
    o = _encode._slicer

    def patched(geometry, up_, normals, device, slicer="auto"):
        ori = (float(sigma) == 0.0) if orient == "auto" else orient  # the patch flips assume clean normals
        ps = slicer_gwn(geometry, up=up_, c_eps=c_eps, sigma=sigma, orient=ori, point_rounds=point_rounds)
        if FIT_K > 0:
            ps.spacing = max(float(ps.spacing), FIT_K * float(sigma))
        return ps, {"input": "points", "points": int(len(geometry)), "spacing": ps.spacing, "normals_given": False,
                    "slicer": "gwn", "unsure": ps.unsure}

    _encode._slicer = patched
    try:
        if cuts == "pointtest":
            sys.path.insert(0, "/home/joao/code/data/meshes/fcp_2710_tests/scripts")
            import pointfit

            return pointfit.encode_pointtest(P, up=up, layers=layers, **kw)
        return fcp.encode(np.asarray(P, float), up=up, layers=layers, device="cuda", **kw)
    finally:
        _encode._slicer = o


def estimate_sigma(ps: PointSlicer, c_est=2.0, rounds=6, k_med=16, normals="grad", return_normals=False):
    """The noise sigma (metric units) from the winding field itself: each point's offset from the mean surface is read
    off its w through the flat-sheet profile w(d) = 1/2 + d / (2 sqrt(d^2 + eps^2)) (inverted), against the median w
    of its k neighbours (the slow part: curvature, thin walls); sigma = 1.4826 x the median |offset|. eps starts at
    c_est spacings and is raised to 4 sigma while sigma grows (the profile saturates when sigma > eps).
    normals="grad": each round the directions become -grad w at the points (DWG's update; smooth at eps, consistent
    with the field), "surf": surf's directions throughout."""
    S = ps._surf
    X, N, idx = S.X, S.N, S.idx[:, :k_med]
    m, k = S.idx.shape
    dk = cp.sqrt(((X[S.idx[:, -1]] - X) ** 2).sum(1))
    area = (np.pi * dk * dk / (k - 1)).astype(cp.float32)
    loc = cp.sqrt(area)
    sig, hist = 0.0, []
    for _ in range(rounds):
        eps = cp.maximum(c_est * loc, cp.float32(4 * sig)).astype(cp.float32)
        E2 = cp.ascontiguousarray(eps * eps)
        V = cp.ascontiguousarray((N * area[:, None]).astype(cp.float32))
        if normals == "grad":
            _, g = wind_grad(X, V, E2, X)
            gn = cp.sqrt((g * g).sum(1, keepdims=True))
            N = cp.where(gn > 0, -g / cp.maximum(gn, 1e-30), N).astype(cp.float32)
            V = cp.ascontiguousarray((N * area[:, None]).astype(cp.float32))
        w = wind(X, V, E2, X)
        r = w - cp.median(w[idx], axis=1)
        r = cp.clip(r, -0.45, 0.45)
        d = eps * 2 * r / cp.sqrt(1 - 4 * r * r)
        new = float(1.4826 * cp.median(cp.abs(d)))
        hist.append(new)
        if new <= 1.02 * sig:
            sig = max(sig, new)
            break
        sig = new
    LOG["sigma_hist"] = hist
    return (sig, N) if return_normals else sig