"""Slice Wrap — algorithmes (numpy seul, aucune dépendance externe).

Enveloppe fermée d'un maillage troué / à doubles parois :
  1. échantillonnage de la surface → voxels occupés (pas = pitch)
  2. fermeture morphologique (dilatation puis érosion, voisinage 6) → bouche les petits trous
  3. remplissage : 3D global + tranche par tranche selon X, Y, Z (union). Un trou de la coque ne fait
     « fuir » que les tranches qu'il traverse, alors qu'un remplissage 3D seul fuit entièrement.
  4. ouverture (retire les filaments)
  5. surface en cubes → Voxel Remesh de Blender (le volume est déjà plein : pas de fuite) → plus grande pièce
     → lissage de Taubin

Les calculs longs sont des **générateurs** : ils cèdent la main régulièrement en produisant
(fraction 0–1, texte ou None) et renvoient leur résultat par `return`. L'opérateur modal les fait avancer par
tranches (Blender reste affiché, barre de progression, Échap pour annuler) ; `run` les exécute d'une traite.
"""
import numpy as np
import bpy
from mathutils import Matrix

# ----------------------------------------------------------------------------- progression


def run(gen):
    """Exécute un générateur de calcul jusqu'au bout et renvoie son résultat."""
    while True:
        try:
            next(gen)
        except StopIteration as e:
            return e.value


def sub(gen, a, b):
    """Relaie la progression de `gen` dans l'intervalle [a, b] du parent ; renvoie le résultat de `gen`.
    Usage : `res = yield from sub(enfant(...), 0.2, 0.5)`."""
    try:
        while True:
            try:
                f, text = next(gen)
            except StopIteration as e:
                return e.value
            yield a + (b - a) * min(max(f, 0.0), 1.0), text
    finally:
        gen.close()        # annulation : le nettoyage de l'enfant s'exécute tout de suite


# ----------------------------------------------------------------------------- voxels


def _shift(a, axis, step):
    """Décale un tableau booléen de `step` (±1) le long de `axis`, sans rebouclage (bord = False)."""
    out = np.zeros_like(a)
    src = [slice(None)] * 3
    dst = [slice(None)] * 3
    if step > 0:
        src[axis], dst[axis] = slice(0, -step), slice(step, None)
    else:
        src[axis], dst[axis] = slice(-step, None), slice(0, step)
    out[tuple(dst)] = a[tuple(src)]
    return out


def dilate(a, iterations=1, axes=(0, 1, 2)):
    for _ in range(iterations):
        b = a.copy()
        for ax in axes:
            b |= _shift(a, ax, 1)
            b |= _shift(a, ax, -1)
        a = b
    return a


def erode(a, iterations=1, axes=(0, 1, 2)):
    return ~dilate(~a, iterations, axes)


def _outside(walls, axes):
    """Vide relié au bord de la grille, propagé seulement le long de `axes` (2 axes = tranche par tranche)."""
    seed = np.zeros_like(walls)
    for ax in axes:
        idx = [slice(None)] * 3
        idx[ax] = 0
        seed[tuple(idx)] = True
        idx[ax] = -1
        seed[tuple(idx)] = True
    free = ~walls
    out = seed & free
    est = max(walls.shape) * 0.5          # nombre de passes probable (estimation pour la progression)
    i = 0
    while True:
        grown = dilate(out, 1, axes) & free
        if np.array_equal(grown, out):
            return out
        out = grown
        i += 1
        yield min(i / est, 0.95), None


def fill(walls, slice_axes=(0, 1, 2)):
    """Volume plein : 3D global ∪ remplissage 2D tranche par tranche pour chaque axe demandé (générateur)."""
    passes = [(0, 1, 2)] + [tuple(a for a in (0, 1, 2) if a != ax) for ax in slice_axes]
    solid = None
    for k, axes in enumerate(passes):
        out = yield from sub(_outside(walls, axes), k / len(passes), (k + 1) / len(passes))
        solid = ~out if solid is None else solid | ~out
    return solid


# ----------------------------------------------------------------------------- maillage ↔ voxels


def mesh_world_arrays(obj, depsgraph):
    """Sommets (monde) et triangles du maillage évalué."""
    ev = obj.evaluated_get(depsgraph)
    me = ev.to_mesh()
    me.calc_loop_triangles()
    nv, nt = len(me.vertices), len(me.loop_triangles)
    co = np.empty(nv * 3, np.float64)
    me.vertices.foreach_get("co", co)
    tri = np.empty(nt * 3, np.int64)
    me.loop_triangles.foreach_get("vertices", tri)
    ev.to_mesh_clear()
    M = np.array(obj.matrix_world)
    co = co.reshape(-1, 3) @ M[:3, :3].T + M[:3, 3]
    return co, tri.reshape(-1, 3)


def sample_surface(co, tri, pitch, max_points=40_000_000):
    """Sommets + points répartis sur les triangles plus grands qu'un demi-voxel, pour que la coque voxel soit
    continue. Sur un maillage déjà dense (scan, IA), seuls les sommets servent : rien n'est dupliqué."""
    e0 = np.linalg.norm(co[tri[:, 1]] - co[tri[:, 0]], axis=1)
    e1 = np.linalg.norm(co[tri[:, 2]] - co[tri[:, 1]], axis=1)
    e2 = np.linalg.norm(co[tri[:, 0]] - co[tri[:, 2]], axis=1)
    k = np.clip(np.ceil(np.maximum(e0, np.maximum(e1, e2)) / (0.5 * pitch)).astype(np.int64), 1, 64)
    del e0, e1, e2
    big = k >= 2                                  # k = 1 : les 3 points seraient les sommets eux-mêmes
    n = ((k[big] + 1) * (k[big] + 2) // 2).sum() + len(co)
    if n > max_points:
        raise RuntimeError(f"{n:,} points d'échantillonnage : augmentez la taille de voxel")
    pts = [co]
    tb, kb = tri[big], k[big]
    for kk in np.unique(kb):
        t = tb[kb == kk]
        a, b, c = co[t[:, 0]], co[t[:, 1]], co[t[:, 2]]
        i, j = np.meshgrid(np.arange(kk + 1), np.arange(kk + 1), indexing="ij")
        m = (i + j) <= kk
        u, v = (i[m] / kk)[None, :, None], (j[m] / kk)[None, :, None]
        pts.append((a[:, None] + u * (b - a)[:, None] + v * (c - a)[:, None]).reshape(-1, 3))
    return np.vstack(pts)


def voxelize(points, pitch, pad):
    lo = points.min(0) - pad * pitch
    idx = np.floor((points - lo) / pitch).astype(np.int64)
    shape = tuple(idx.max(0) + pad + 1)
    occ = np.zeros(shape, bool)
    occ[idx[:, 0], idx[:, 1], idx[:, 2]] = True
    return occ, lo


def blocky_mesh(solid, lo, pitch, name):
    """Surface en cubes du volume plein (faces entre voxel plein et vide), sommets partagés."""
    s = np.pad(solid, 1)
    lo = lo - pitch
    quads = []
    corner = {  # pour chaque axe : 4 coins de la face (dans le plan des 2 autres axes)
        0: [(0, 0, 0), (0, 1, 0), (0, 1, 1), (0, 0, 1)],
        1: [(0, 0, 0), (0, 0, 1), (1, 0, 1), (1, 0, 0)],
        2: [(0, 0, 0), (1, 0, 0), (1, 1, 0), (0, 1, 0)],
    }
    for ax in (0, 1, 2):
        a = s
        b = np.roll(s, -1, axis=ax)
        for flip, mask in ((False, a & ~b), (True, ~a & b)):
            p = np.argwhere(mask)
            if not len(p):
                continue
            p[:, ax] += 1                                   # la face est sur le plan i+1
            cs = corner[ax][::-1] if flip else corner[ax]
            quads.append(np.stack([p + np.array(c) for c in cs], axis=1))
    q = np.concatenate(quads)                               # (F, 4, 3) coins en indices de grille
    flat = q.reshape(-1, 3)
    uniq, inv = np.unique(flat, axis=0, return_inverse=True)
    verts = uniq * pitch + lo
    faces = inv.reshape(-1, 4)
    me = bpy.data.meshes.new(name)
    me.from_pydata(verts.astype(np.float32), [], faces.astype(np.int32))
    me.validate()
    return me


# ----------------------------------------------------------------------------- lissage / nettoyage


def _edges(me):
    e = np.empty(len(me.edges) * 2, np.int64)
    me.edges.foreach_get("vertices", e)
    return e.reshape(-1, 2)


def _neighbor_mean(co, e):
    acc = np.zeros_like(co)
    cnt = np.zeros(len(co))
    np.add.at(acc, e[:, 0], co[e[:, 1]])
    np.add.at(acc, e[:, 1], co[e[:, 0]])
    np.add.at(cnt, e[:, 0], 1)
    np.add.at(cnt, e[:, 1], 1)
    return acc / np.maximum(cnt, 1)[:, None], cnt


def taubin(me, iterations=20, lamb=0.5, mu=-0.53):
    """Lissage de Taubin : efface les marches et le moiré sans faire rétrécir la forme (générateur)."""
    if iterations <= 0:
        return
    co = np.empty(len(me.vertices) * 3)
    me.vertices.foreach_get("co", co)
    co = co.reshape(-1, 3)
    e = _edges(me)
    for i in range(iterations):
        for f in (lamb, mu):
            m, _ = _neighbor_mean(co, e)
            co = co + f * (m - co)
        yield (i + 1) / iterations, None
    me.vertices.foreach_set("co", co.ravel())
    me.update()


def keep_largest_island(me):
    """Supprime les îlots détachés (garde la plus grande pièce du maillage) (générateur)."""
    import bmesh
    bm = bmesh.new()
    try:
        bm.from_mesh(me)
        nv = max(len(bm.verts), 1)
        seen, best = set(), []
        last = 0
        for v in bm.verts:
            if v in seen:
                continue
            comp, stack = [], [v]
            seen.add(v)
            while stack:
                x = stack.pop()
                comp.append(x)
                for ed in x.link_edges:
                    w = ed.other_vert(x)
                    if w not in seen:
                        seen.add(w)
                        stack.append(w)
            if len(comp) > len(best):
                best = comp
            if len(seen) - last >= 50_000:
                last = len(seen)
                yield 0.9 * len(seen) / nv, None
        keep = set(best)
        bmesh.ops.delete(bm, geom=[v for v in bm.verts if v not in keep], context="VERTS")
        bm.to_mesh(me)
    finally:
        bm.free()


def despike(me, threshold, iterations=3, ratio=3.0, strength=1.0, weights=None):
    """Filtre anti-pointes (générateur ; le maillage n'est écrit qu'à la fin). Un sommet est une pointe si
      - son écart à la moyenne de ses voisins dépasse `threshold`, ET
      - cet écart est `ratio` fois plus grand que l'écart moyen de ses voisins (une pointe est isolée ;
        dans une zone simplement courbe — oreille, paupière — tous les voisins dévient pareil).
    La pointe est ramenée vers la moyenne de ses voisins de `strength` × `weights` (poids 0–1 par sommet,
    None = 1 partout) : le reste du relief n'est pas lissé. Renvoie le nombre de corrections."""
    co = np.empty(len(me.vertices) * 3)
    me.vertices.foreach_get("co", co)
    co = co.reshape(-1, 3)
    e = _edges(me)
    w = np.full(len(co), float(strength)) if weights is None else np.clip(weights, 0.0, 1.0) * strength
    active = w > 0
    total = 0
    for it in range(iterations):
        m, cnt = _neighbor_mean(co, e)
        d = np.linalg.norm(co - m, axis=1)
        acc = np.zeros(len(co))
        np.add.at(acc, e[:, 0], d[e[:, 1]])
        np.add.at(acc, e[:, 1], d[e[:, 0]])
        dn = acc / np.maximum(cnt, 1)
        bad = active & (d > threshold) & (d > ratio * dn) & (cnt > 0)
        if not bad.any():
            break
        co[bad] += (m[bad] - co[bad]) * w[bad, None]
        total += int(bad.sum())
        yield (it + 1) / iterations, None
    me.vertices.foreach_set("co", co.ravel())
    me.update()
    return total


def group_weights(me, group_indices):
    """Poids max par sommet sur les groupes de sommets donnés (indices) ; None si aucun groupe (générateur)."""
    if not group_indices:
        return None
    gi = set(group_indices)
    nv = max(len(me.vertices), 1)
    w = np.zeros(len(me.vertices))
    for v in me.vertices:
        best = 0.0
        for g in v.groups:
            if g.group in gi and g.weight > best:
                best = g.weight
        w[v.index] = best
        if v.index % 50_000 == 49_999:
            yield v.index / nv, None
    return w


def fix_missed_rays(me, before, inflate, low=0.3, deep=2.0, max_iter=500):
    """Répare les ratés de projection, connus exactement en comparant chaque sommet avant / après le Shrinkwrap :
      - raté « vide »  : le sommet n'a presque pas bougé (< low × gonflement) → aiguille vers l'extérieur ;
      - raté « profond » : il a bougé plus de `deep` fois le déplacement moyen de ses voisins → creux (paroi interne).
    Les sommets ratés sont reconstruits par diffusion depuis leurs voisins réussis. Renvoie le nombre réparé
    (générateur ; le maillage n'est écrit qu'à la fin)."""
    co = np.empty(len(me.vertices) * 3)
    me.vertices.foreach_get("co", co)
    co = co.reshape(-1, 3)
    disp = np.linalg.norm(co - before, axis=1)
    e = _edges(me)
    acc = np.zeros(len(co))
    cnt = np.zeros(len(co))
    np.add.at(acc, e[:, 0], disp[e[:, 1]])
    np.add.at(acc, e[:, 1], disp[e[:, 0]])
    np.add.at(cnt, e[:, 0], 1)
    np.add.at(cnt, e[:, 1], 1)
    disp_n = acc / np.maximum(cnt, 1)
    bad = (disp < low * inflate) | (disp > deep * np.maximum(disp_n, low * inflate))
    nbad = int(bad.sum())
    known = ~bad
    for _ in range(max_iter):
        if known.all():
            break
        s_acc = np.zeros_like(co)
        s_cnt = np.zeros(len(co))
        kf = known.astype(float)
        np.add.at(s_acc, e[:, 0], co[e[:, 1]] * kf[e[:, 1], None])
        np.add.at(s_acc, e[:, 1], co[e[:, 0]] * kf[e[:, 0], None])
        np.add.at(s_cnt, e[:, 0], kf[e[:, 1]])
        np.add.at(s_cnt, e[:, 1], kf[e[:, 0]])
        fill = ~known & (s_cnt > 0)
        co[fill] = s_acc[fill] / s_cnt[fill, None]
        known |= fill
        yield 0.95 * float(known.mean()), None
    for _ in range(3):   # adoucit la reprise sur les zones reconstruites seulement
        m, _c = _neighbor_mean(co, e)
        co[bad] = m[bad]
    me.vertices.foreach_set("co", co.ravel())
    me.update()
    return nbad


def mean_edge_length(me):
    co = np.empty(len(me.vertices) * 3)
    me.vertices.foreach_get("co", co)
    co = co.reshape(-1, 3)
    e = _edges(me)
    return float(np.linalg.norm(co[e[:, 0]] - co[e[:, 1]], axis=1).mean())


# ----------------------------------------------------------------------------- enveloppe complète


def unique_object_name(base):
    """`base` s'il est libre, sinon base.001, base.002… (comme Blender) : on n'écrase jamais un objet existant."""
    if base not in bpy.data.objects:
        return base
    n = 1
    while f"{base}.{n:03d}" in bpy.data.objects:
        n += 1
    return f"{base}.{n:03d}"


def build_envelope(src, pitch, close_iter=3, slice_axes=(0, 1, 2), open_iter=1, smooth_iter=20,
                   largest=True, name=None, report=print, collection=None, replace=None):
    """Crée l'objet enveloppe de `src` sous un nom libre dérivé de `name` (numérotation .001… si besoin), rangé dans
    `collection` (ou la fonction qui la renvoie ; None = collections de la source). `replace` = objet existant à recalculer sur place
    (Regenerate) : seul cas où un maillage existant est remplacé. Renvoie l'objet (générateur).
    La scène n'est modifiée qu'à la toute fin : une annulation en cours de route ne laisse aucun demi-résultat."""
    import time
    t0 = time.time()
    yield 0.0, "Lecture du maillage source…"
    dg = bpy.context.evaluated_depsgraph_get()
    co, tri = mesh_world_arrays(src, dg)
    yield 0.03, "Échantillonnage de la surface…"
    pts = sample_surface(co, tri, pitch)
    del co, tri
    pad = close_iter + 3
    yield 0.06, "Voxelisation…"
    occ, lo = voxelize(pts, pitch, pad)
    del pts
    report(f"Grille {occ.shape} — {int(occ.sum()):,} voxels de surface")
    yield 0.08, "Fermeture des trous…"
    walls = erode(dilate(occ, close_iter), close_iter) if close_iter > 0 else occ
    del occ
    yield 0.1, "Remplissage tranche par tranche…"
    solid = yield from sub(fill(walls, slice_axes), 0.1, 0.6)
    del walls
    if open_iter > 0:
        yield 0.6, "Retrait des filaments…"
        solid = dilate(erode(solid, open_iter), open_iter)
    report(f"Volume plein : {int(solid.sum()):,} voxels ({time.time() - t0:.1f} s)")
    if replace is not None:
        name = replace.name
    else:
        name = unique_object_name(name or f"{src.name}_slicewrapped")
    yield 0.62, "Surface du volume…"
    me = blocky_mesh(solid, lo, pitch, name)
    del solid
    new_me = tmp = None
    try:
        # Voxel Remesh de Blender sur un volume DÉJÀ plein et fermé → aucune fuite possible.
        # Objet temporaire : le Remesh est évalué hors de l'objet final.
        yield 0.66, "Voxel Remesh (Blender)…"
        tmp = bpy.data.objects.new(name + "_sw_tmp", me)
        bpy.context.scene.collection.objects.link(tmp)
        mod = tmp.modifiers.new("SliceWrap_Remesh", "REMESH")
        mod.mode = "VOXEL"
        mod.voxel_size = pitch
        mod.adaptivity = 0.0
        dg = bpy.context.evaluated_depsgraph_get()
        new_me = bpy.data.meshes.new_from_object(tmp.evaluated_get(dg), depsgraph=dg)
        bpy.data.objects.remove(tmp)
        tmp = None
        bpy.data.meshes.remove(me)
        me = None
        if largest:
            yield 0.75, "Plus grande pièce…"
            yield from sub(keep_largest_island(new_me), 0.75, 0.85)
        yield 0.85, "Lissage…"
        yield from sub(taubin(new_me, smooth_iter), 0.85, 1.0)
        new_me.polygons.foreach_set("use_smooth", [True] * len(new_me.polygons))
        new_me.update()
    except BaseException:          # erreur ou annulation : on retire tout ce qui a été créé
        if tmp is not None:
            bpy.data.objects.remove(tmp)
        for m in (me, new_me):
            if m is not None:
                bpy.data.meshes.remove(m)
        raise
    if replace is not None:
        obj = replace
        old_me = obj.data
        obj.data = new_me
        if old_me.users == 0:
            bpy.data.meshes.remove(old_me)
    else:
        obj = bpy.data.objects.new(name, new_me)
        if callable(collection):          # collection créée seulement à la fin (rien de laissé si annulation)
            collection = collection()
        for c in ([collection] if collection is not None else src.users_collection):
            c.objects.link(obj)
    new_me.name = name
    obj.matrix_world = Matrix.Identity(4)
    report(f"Enveloppe {obj.name} : {len(obj.data.polygons):,} faces ({time.time() - t0:.1f} s)")
    return obj


# ----------------------------------------------------------------------------- analyse de forme

SHAPES = ("Flat", "Bulge", "Hollow", "Ridge", "Fold", "Saddle")
SHAPE_COLORS = np.array([            # couleurs d'affichage des tags (linéaire)
    (0.55, 0.55, 0.55),   # Flat    gris
    (0.90, 0.25, 0.15),   # Bulge   rouge   (dôme convexe)
    (0.15, 0.35, 0.90),   # Hollow  bleu    (cuvette concave)
    (0.95, 0.75, 0.10),   # Ridge   jaune   (crête convexe, une direction)
    (0.10, 0.75, 0.85),   # Fold    cyan    (pli / vallée concave, une direction)
    (0.70, 0.25, 0.85),   # Saddle  violet  (selle)
])


def _diffuse(vals, e, cnt, iterations):
    """Moyenne sur le voisinage (k anneaux) d'un champ par sommet : (V, ...) → même forme (générateur)."""
    flat = vals.reshape(len(vals), -1)
    for i in range(iterations):
        acc = flat.copy()
        np.add.at(acc, e[:, 0], flat[e[:, 1]])
        np.add.at(acc, e[:, 1], flat[e[:, 0]])
        flat = acc / (cnt + 1.0)[:, None]
        yield (i + 1) / max(iterations, 1), None
    return flat.reshape(vals.shape)


def analyze_shape(co, normals, e, rings=4, flat_threshold=0.15, aniso_threshold=0.55):
    """Analyse locale de la forme sur un voisinage de `rings` anneaux.
    Tenseur de covariance des normales T = moyenne(n nᵀ) ; valeurs propres λ1 ≥ λ2 ≥ λ3 :
      - variation = λ2 + λ3 (0 = plat)              → `detail` 0–1 (normalisé au 95e centile)
      - anisotropie = (λ2 − λ3) / (λ2 + λ3)          → ~1 : la surface tourne dans UNE direction (pli, crête, doigt)
      - direction du pli = vecteur propre de λ3      → direction le long de laquelle la surface ne tourne pas
    Courbure moyenne signée H ≈ −(moyenne des voisins − p)·n / arête² (> 0 convexe, < 0 concave).
    Tags : Flat, Bulge (dôme), Hollow (cuvette), Ridge (crête), Fold (pli), Saddle (selle).
    Renvoie dict(detail, curvature, anisotropy, shape, flow) (générateur)."""
    n = normals / np.maximum(np.linalg.norm(normals, axis=1, keepdims=True), 1e-12)
    cnt = np.zeros(len(co))
    np.add.at(cnt, e[:, 0], 1)
    np.add.at(cnt, e[:, 1], 1)
    T = n[:, :, None] * n[:, None, :]
    T = yield from sub(_diffuse(T, e, cnt, rings), 0.0, 0.7)
    w, v = np.linalg.eigh(T)                       # valeurs propres croissantes : w[:,0] = λ3
    lam3, lam2 = np.clip(w[:, 0], 0, None), np.clip(w[:, 1], 0, None)
    var = lam2 + lam3
    aniso = (lam2 - lam3) / np.maximum(var, 1e-12)
    flow = v[:, :, 0]                              # vecteur propre de λ3
    m, _ = _neighbor_mean(co, e)
    el = np.linalg.norm(co[e[:, 0]] - co[e[:, 1]], axis=1).mean()
    H = -np.einsum("ij,ij->i", m - co, n) / max(el * el, 1e-12)
    H = yield from sub(_diffuse(H, e, cnt, rings), 0.75, 1.0)
    scale = max(np.percentile(var, 95), 1e-12)
    detail = np.clip(var / scale, 0, 1)
    hs = max(np.percentile(np.abs(H), 95), 1e-12)
    curv = np.clip(H / hs, -1, 1)
    shape = np.zeros(len(co), np.int32)
    active = detail >= flat_threshold
    iso = aniso < aniso_threshold
    shape[active & iso & (curv > 0.15)] = 1
    shape[active & iso & (curv < -0.15)] = 2
    shape[active & ~iso & (curv >= 0)] = 3
    shape[active & ~iso & (curv < 0)] = 4
    shape[active & iso & (np.abs(curv) <= 0.15)] = 5
    return {"detail": detail, "curvature": curv, "anisotropy": aniso, "shape": shape, "flow": flow}


def mesh_arrays(me):
    co = np.empty(len(me.vertices) * 3)
    me.vertices.foreach_get("co", co)
    nr = np.empty(len(me.vertices) * 3)
    me.vertex_normals.foreach_get("vector", nr)
    return co.reshape(-1, 3), nr.reshape(-1, 3), _edges(me)


def write_analysis(obj, res, count):
    """Écrit l'analyse sur les `count` premiers sommets du maillage de base de `obj` (les sommets d'origine
    viennent en premier dans un maillage subdivisé) : attributs, couleurs d'affichage, groupes de sommets.
    Les couleurs portent le suffixe _Color : groupes de sommets et attributs partagent le même espace de noms."""
    me = obj.data
    nb = len(me.vertices)
    k = min(nb, count)

    def put(name, dtype, vals):
        if name in me.attributes:
            me.attributes.remove(me.attributes[name])
        a = me.attributes.new(name, dtype, "POINT")
        full = np.zeros((nb,) + vals.shape[1:], vals.dtype)
        full[:k] = vals[:k]
        a.data.foreach_set("value" if dtype != "FLOAT_VECTOR" else "vector", full.ravel())

    put("sw_detail", "FLOAT", res["detail"].astype(np.float32))
    put("sw_curvature", "FLOAT", res["curvature"].astype(np.float32))
    put("sw_anisotropy", "FLOAT", res["anisotropy"].astype(np.float32))
    put("sw_shape", "INT", res["shape"].astype(np.int32))
    put("sw_flow", "FLOAT_VECTOR", res["flow"].astype(np.float32))

    def color(name, rgb):
        if name in me.color_attributes:
            me.color_attributes.remove(me.color_attributes[name])
        ca = me.color_attributes.new(name, "FLOAT_COLOR", "POINT")
        full = np.ones((nb, 4), np.float32)
        full[:k, :3] = rgb[:k]
        ca.data.foreach_set("color", full.ravel())
        return ca

    d = res["detail"][:, None]
    ramp = np.where(d < 0.5, (1 - 2 * d) * np.array([0.05, 0.15, 0.6]) + 2 * d * np.array([0.1, 0.8, 0.2]),
                    (2 - 2 * d) * np.array([0.1, 0.8, 0.2]) + (2 * d - 1) * np.array([0.9, 0.1, 0.05]))
    color("SW_Detail_Color", ramp.astype(np.float32))
    color("SW_Shape_Color", SHAPE_COLORS[res["shape"]].astype(np.float32))
    # groupes de sommets : densité + un groupe par tag (utilisables comme filtres Despike, masques, etc.)
    groups = {"SW_Detail": res["detail"]}
    for i, nm in enumerate(SHAPES[1:], start=1):
        groups[f"SW_{nm}"] = (res["shape"] == i).astype(float)
    for nm, w in groups.items():
        vg = obj.vertex_groups.get(nm) or obj.vertex_groups.new(name=nm)
        vg.remove(list(range(nb)))
        w = np.asarray(w[:k], float)
        for val in np.unique(np.round(w, 2)):
            if val <= 0:
                continue
            idx = np.nonzero(np.round(w, 2) == val)[0].tolist()
            vg.add(idx, float(val), "REPLACE")
    me.update()
    return {nm: int((res["shape"][:k] == i).sum()) for i, nm in enumerate(SHAPES)}