583 lines
24 KiB
Python
583 lines
24 KiB
Python
"""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)}
|