From b49f4e08b5461c46a67b788a3d18c0dd01586ac2 Mon Sep 17 00:00:00 2001 From: PozzettiAndrea Date: Wed, 17 Jun 2026 09:03:48 +0200 Subject: [PATCH] Add Region Growing node + Decrease-Gaussian sharpen backend Region Growing (analysis): seed faces flattest-first and grow regions by a tunable similarity metric -- criterion dropdown (normals / normals+area / curvature+normals / curvature) with per-term weights and a threshold. Greedy flood-fill against the region's running area-weighted representative; small regions merged into their strongest neighbour. Outputs per-face region_id. On 101 (CAD): by normals it recovers the planar structure (top region ~79k faces, ~70% of faces in large regions) in ~6s. Decrease Gaussian (sharpen backend): gradient flow on the Gaussian-curvature energy sum K_i^2 (K_i = vertex angle defect) + Tikhonov anchor, so |K| is driven toward zero -> locally developable (planes/cylinders/cones). Exact analytic angle-defect gradient (verified vs finite differences to 1e-10); every step passes through the foldless signed-area barrier. Verified: monotone decrease of mean|K|, total K conserved (Gauss-Bonnet, 2*pi*chi), watertight/foldless/volume-stable. Wired into Sharpen Mesh. (Includes segment_patches.py, previously-untracked, so the analysis __init__ imports resolve on a clean checkout.) Co-Authored-By: Claude Opus 4.8 --- nodes/analysis/__init__.py | 8 + nodes/analysis/region_growing.py | 250 ++++++++++++++++++ nodes/analysis/segment_patches.py | 192 ++++++++++++++ nodes/smoothing/sharpen.py | 19 ++ nodes/smoothing/sharpen_backends/__init__.py | 3 + .../sharpen_backends/decrease_gaussian.py | 184 +++++++++++++ 6 files changed, 656 insertions(+) create mode 100644 nodes/analysis/region_growing.py create mode 100644 nodes/analysis/segment_patches.py create mode 100644 nodes/smoothing/sharpen_backends/decrease_gaussian.py diff --git a/nodes/analysis/__init__.py b/nodes/analysis/__init__.py index abf330b8..2a025f51 100644 --- a/nodes/analysis/__init__.py +++ b/nodes/analysis/__init__.py @@ -27,6 +27,10 @@ from .compute_fillet_field import NODE_CLASS_MAPPINGS as FILLET_FIELD_MAPPINGS from .compute_fillet_field import NODE_DISPLAY_NAME_MAPPINGS as FILLET_FIELD_DISPLAY from .segment_by_curvature import NODE_CLASS_MAPPINGS as SEGMENT_CURVATURE_MAPPINGS from .segment_by_curvature import NODE_DISPLAY_NAME_MAPPINGS as SEGMENT_CURVATURE_DISPLAY +from .segment_patches import NODE_CLASS_MAPPINGS as SEGMENT_PATCHES_MAPPINGS +from .segment_patches import NODE_DISPLAY_NAME_MAPPINGS as SEGMENT_PATCHES_DISPLAY +from .region_growing import NODE_CLASS_MAPPINGS as REGION_GROWING_MAPPINGS +from .region_growing import NODE_DISPLAY_NAME_MAPPINGS as REGION_GROWING_DISPLAY # Combine all mappings NODE_CLASS_MAPPINGS = { @@ -42,6 +46,8 @@ NODE_CLASS_MAPPINGS = { **COMPUTE_CURVATURE_MAPPINGS, **FILLET_FIELD_MAPPINGS, **SEGMENT_CURVATURE_MAPPINGS, + **SEGMENT_PATCHES_MAPPINGS, + **REGION_GROWING_MAPPINGS, } NODE_DISPLAY_NAME_MAPPINGS = { @@ -57,6 +63,8 @@ NODE_DISPLAY_NAME_MAPPINGS = { **COMPUTE_CURVATURE_DISPLAY, **FILLET_FIELD_DISPLAY, **SEGMENT_CURVATURE_DISPLAY, + **SEGMENT_PATCHES_DISPLAY, + **REGION_GROWING_DISPLAY, } __all__ = ['NODE_CLASS_MAPPINGS', 'NODE_DISPLAY_NAME_MAPPINGS'] diff --git a/nodes/analysis/region_growing.py b/nodes/analysis/region_growing.py new file mode 100644 index 00000000..c31726f8 --- /dev/null +++ b/nodes/analysis/region_growing.py @@ -0,0 +1,250 @@ +# SPDX-License-Identifier: GPL-3.0-or-later +# Copyright (C) 2025 ComfyUI-GeometryPack Contributors + +"""Region growing segmentation: seed faces and grow regions by SIMILARITY. + +Unlike segment_patches (which cuts at sharp dihedral edges), this clusters faces +by a tunable similarity metric -- a face joins a region while it stays close to +the region's running representative in: + - normal direction (group co-planar / smoothly-varying faces) + - + face area weighting (large faces dominate the representative -> robust + region normals, the VSA-style area-weighted L2,1 metric) + - + curvature (also require similar mean curvature -> separate a + flat from a fillet even when their normals agree at + the seam) +Seeds are taken flattest-first (lowest curvature), so big planar regions form +first and creases/fillets fall out as their own regions. Each component term has +its own weight. Outputs per-face `region_id` (1..R) for colouring / Split By Field +/ per-region primitive fitting. +""" + +import logging +from collections import deque + +import numpy as np +import trimesh as trimesh_module +from comfy_api.latest import io + +log = logging.getLogger("geometrypack") + +_CRITERIA = { + "normals": (True, False, False), + "normals + area": (True, False, True), + "curvature + normals": (True, True, False), + "curvature": (False, True, False), +} + + +def _face_mean_curvature(V, F): + """Per-face |mean curvature|: vertex mean-curvature (cotan Laplacian) averaged + to faces. Robust, cheap, and only used for seed ordering + the curvature term.""" + import igl + L = igl.cotmatrix(V, F) + M = igl.massmatrix(V, F, igl.MASSMATRIX_TYPE_BARYCENTRIC) + Minv = 1.0 / np.clip(np.asarray(M.diagonal()), 1e-12, None) + Hv = 0.5 * np.linalg.norm(Minv[:, None] * (L @ V), axis=1) # per-vertex |H| + return Hv[F].mean(axis=1) # per-face |H| + + +def _build_face_neighbors_csr(adj_pairs, m): + """CSR (indptr, indices) face-adjacency from the shared-edge face pairs.""" + a = adj_pairs[:, 0] + b = adj_pairs[:, 1] + src = np.concatenate([a, b]) + dst = np.concatenate([b, a]) + order = np.argsort(src, kind="stable") + src_s = src[order] + dst_s = dst[order] + indptr = np.zeros(m + 1, dtype=np.int64) + np.add.at(indptr, src_s + 1, 1) + np.cumsum(indptr, out=indptr) + return indptr, dst_s.astype(np.int64) + + +def _region_grow(N, areas, Hface, indptr, indices, use_normal, use_curv, use_area, + w_normal, w_curv, threshold): + """Greedy flood-fill region growing. A face joins a region while its weighted + dissimilarity to the region's running (area-weighted) representative stays + below `threshold`. Seeds flattest-first. Returns per-face region id (1..R).""" + m = len(N) + cscale = float(Hface.std()) + 1e-9 + order = np.argsort(Hface, kind="stable") # flattest faces seed first + labels = np.full(m, -1, dtype=np.int64) + rid = 0 + for seed in order: + if labels[seed] >= 0: + continue + rid += 1 + labels[seed] = rid + w0 = float(areas[seed]) if use_area else 1.0 + n_acc = N[seed] * w0 # running (weighted) normal sum + c_acc = float(Hface[seed]) * w0 # running weighted curvature sum + w_sum = w0 + q = deque((seed,)) + while q: + g = q.popleft() + n_rep = n_acc / (np.linalg.norm(n_acc) + 1e-12) + c_rep = c_acc / w_sum + for idx in range(indptr[g], indptr[g + 1]): + h = indices[idx] + if labels[h] >= 0: + continue + d = 0.0 + if use_normal: + d += w_normal * (1.0 - float(np.dot(N[h], n_rep))) + if use_curv: + d += w_curv * abs(float(Hface[h]) - c_rep) / cscale + if d <= threshold: + labels[h] = rid + wh = float(areas[h]) if use_area else 1.0 + n_acc = n_acc + N[h] * wh + c_acc += float(Hface[h]) * wh + w_sum += wh + q.append(h) + return labels + + +def _merge_small(labels, indptr, indices, areas, min_faces): + """Merge regions smaller than min_faces into the adjacent region they share the + most faces with. Iterates to a fixed point; relabels 1..R contiguously.""" + if min_faces <= 1: + uniq, inv = np.unique(labels, return_inverse=True) + return (inv + 1).astype(np.int64) + labels = labels.copy() + for _ in range(1000): + uniq, inv = np.unique(labels, return_inverse=True) + labels = inv.astype(np.int64) + counts = np.bincount(labels) + small = np.where(counts < min_faces)[0] + if len(small) == 0: + break + small_set = set(small.tolist()) + changed = False + for r in small: + members = np.where(labels == r)[0] + # tally neighbour-region face counts + nbr = {} + for g in members: + for idx in range(indptr[g], indptr[g + 1]): + rl = labels[indices[idx]] + if rl != r: + nbr[rl] = nbr.get(rl, 0) + 1 + if nbr: + best = max(nbr, key=nbr.get) + labels[members] = best + changed = True + if not changed: + break + uniq, inv = np.unique(labels, return_inverse=True) + return (inv + 1).astype(np.int64) + + +class RegionGrowingNode(io.ComfyNode): + """Seed + grow regions, clustering faces by a tunable similarity metric.""" + + @classmethod + def define_schema(cls): + return io.Schema( + node_id="GeomPackRegionGrowing", + display_name="Region Growing", + category="geompack/analysis", + description=( + "Segment a mesh by SEEDING faces (flattest first) and GROWING regions: a " + "face joins while it stays similar to the region's running representative. " + "Choose what 'similar' means in `criterion` (normals / normals+area / " + "curvature+normals / curvature) and tune each component's weight + the " + "`threshold` (lower = more, smaller regions). Outputs per-face region_id " + "(1..R) -- colour it, Split By Field, or fit a primitive per region." + ), + is_output_node=True, + inputs=[ + io.Custom("TRIMESH").Input("trimesh"), + io.Combo.Input("criterion", + options=list(_CRITERIA.keys()), default="curvature + normals", + tooltip=( + "What makes two faces belong to the same region:\n" + " normals = co-planar / smoothly-varying normals only;\n" + " normals + area = same, but big faces dominate the region normal " + "(area-weighted, robust on uneven tessellation);\n" + " curvature + normals = also require similar mean curvature (separates a " + "flat from a fillet that meet tangentially);\n" + " curvature = group by curvature magnitude alone.")), + io.Float.Input("threshold", default=0.05, min=0.0005, max=2.0, step=0.005, display_mode="number", tooltip=( + "Similarity tolerance to ADD a face to the growing region. LOWER = stricter " + "= more, smaller regions; HIGHER = looser = fewer, larger regions. The " + "normal term is 1-cos(angle) (0.015~10deg, 0.06~20deg, 0.13~30deg); the " + "curvature term is in units of the mesh's curvature std. Default 0.05.")), + io.Float.Input("normal_weight", default=1.0, min=0.0, max=20.0, step=0.1, tooltip=( + "Weight of the NORMAL-difference term (1-cos angle between the face normal " + "and the region's representative normal). Higher = normals matter more. " + "Ignored unless the criterion includes normals.")), + io.Float.Input("curvature_weight", default=1.0, min=0.0, max=20.0, step=0.1, tooltip=( + "Weight of the CURVATURE-difference term (|H_face - H_region| / curvature " + "std). Higher = curvature matters more, so a fillet won't merge into the " + "abutting flat. Ignored unless the criterion includes curvature.")), + io.Int.Input("min_region_size", default=10, min=1, max=1000000, step=1, tooltip=( + "Regions with fewer faces than this are merged into the neighbour they " + "share the most boundary with (de-speckle). 1 = keep every region.")), + ], + outputs=[ + io.Custom("TRIMESH").Output(display_name="segmented_mesh"), + io.String.Output(display_name="info"), + ], + ) + + @classmethod + def execute(cls, trimesh, criterion="curvature + normals", threshold=0.05, + normal_weight=1.0, curvature_weight=1.0, min_region_size=10): + mesh = trimesh.copy() + try: + mesh.merge_vertices() + mesh.update_faces(mesh.nondegenerate_faces()) + mesh.remove_unreferenced_vertices() + except Exception as e: + log.debug("preclean skipped: %s", e) + + V = np.ascontiguousarray(mesh.vertices, dtype=np.float64) + F = np.ascontiguousarray(mesh.faces, dtype=np.int64) + m = len(F) + use_normal, use_curv, use_area = _CRITERIA[criterion] + + N = np.asarray(mesh.face_normals, dtype=np.float64) + areas = np.asarray(mesh.area_faces, dtype=np.float64) + Hface = _face_mean_curvature(V, F) # also used for seed ordering + + adj = np.asarray(mesh.face_adjacency) + if len(adj) == 0: + raise ValueError("Mesh has no face adjacency (disconnected or degenerate).") + indptr, indices = _build_face_neighbors_csr(adj, m) + + labels = _region_grow(N, areas, Hface, indptr, indices, + use_normal, use_curv, use_area, + float(normal_weight), float(curvature_weight), float(threshold)) + labels = _merge_small(labels, indptr, indices, areas, int(min_region_size)) + n_reg = int(labels.max()) + + mesh.face_attributes["region_id"] = labels.astype(np.int64) + # also a per-vertex copy (majority face label) for vertex-field consumers + vlab = np.zeros(len(V), dtype=np.float32) + cnt = np.zeros(len(V)) + for k in range(3): + np.add.at(vlab, F[:, k], labels) + np.add.at(cnt, F[:, k], 1.0) + mesh.vertex_attributes["region_id"] = (vlab / np.maximum(cnt, 1.0)).astype(np.float32) + + sizes = np.bincount(labels)[1:] + info = ( + f"Region Growing:\n\n" + f"Faces: {m:,} | regions: {n_reg:,}\n" + f"criterion: {criterion} (n_w={normal_weight}, c_w={curvature_weight})\n" + f"threshold: {threshold} | min_region_size: {min_region_size}\n\n" + f"region face-count: min {int(sizes.min())}, median {int(np.median(sizes))}, " + f"max {int(sizes.max())}\n" + f"Field: face_attributes['region_id'] (1..{n_reg})" + ) + log.info("Region Growing: %d regions (%s, thr=%.3f)", n_reg, criterion, threshold) + return io.NodeOutput(mesh, info, ui={"text": [info]}) + + +NODE_CLASS_MAPPINGS = {"GeomPackRegionGrowing": RegionGrowingNode} +NODE_DISPLAY_NAME_MAPPINGS = {"GeomPackRegionGrowing": "Region Growing"} diff --git a/nodes/analysis/segment_patches.py b/nodes/analysis/segment_patches.py new file mode 100644 index 00000000..3e531bd3 --- /dev/null +++ b/nodes/analysis/segment_patches.py @@ -0,0 +1,192 @@ +# SPDX-License-Identifier: GPL-3.0-or-later +# Copyright (C) 2025 ComfyUI-GeometryPack Contributors + +"""Segment a mesh into CADable patches and OUTPUT them (mesh + per-face region_id). + +Pipeline: + 1. feature edges = dihedral angle >= crease_angle_deg (optionally broadened over + `wide_rings` so chamfers / fine fillets read as one turn, not many tiny ones). + 2. flood-fill faces into regions, cutting at those feature edges. + 3. CLEAN UP: iteratively merge each region smaller than min_patch_faces (or + min_patch_area_pct of the surface) into the neighbour it shares the most + boundary with. This is what turns a speckled segmentation into a usable set + of patches. + +Outputs the mesh with face_attributes['region_id'] (1..R) so you can colour it in +Preview Mesh (fields), split it with Split By Field, or fit primitives per patch. +""" + +import logging + +import numpy as np +import trimesh as trimesh_module +from comfy_api.latest import io + +from ..visualization.preview_mesh_boundaries import _edge_values, _face_regions + +log = logging.getLogger("geometrypack") + + +def _merge_small_regions(region_id, adj_pairs, wall_mask, face_areas, min_faces, min_area): + """Merge regions below the size thresholds into their strongest neighbour + (most shared boundary edges, ties broken by neighbour size). Iterates to a + fixed point. Returns cleaned, contiguously-relabelled region_id (1..R).""" + from collections import defaultdict + + labels = np.asarray(region_id, dtype=np.int64).copy() + adj_pairs = np.asarray(adj_pairs).reshape(-1, 2) + wall = np.asarray(wall_mask, dtype=bool) + walls = adj_pairs[wall] if len(adj_pairs) else np.zeros((0, 2), np.int64) + areas = np.asarray(face_areas, dtype=np.float64) + + def relabel(lbl): + uniq, inv = np.unique(lbl, return_inverse=True) + return (inv + 1).astype(np.int64) + + for _ in range(10000): # converges fast (region count strictly drops per pass) + labels = relabel(labels) + n = int(labels.max()) if len(labels) else 0 + fcount = np.bincount(labels, minlength=n + 1) + acount = np.bincount(labels, weights=areas, minlength=n + 1) + + small = [r for r in range(1, n + 1) + if fcount[r] > 0 and (fcount[r] < min_faces or acount[r] < min_area)] + if not small: + break + + # region-adjacency strength from the wall edges (boundary between regions) + rl = labels[walls[:, 0]] if len(walls) else np.zeros(0, np.int64) + rr = labels[walls[:, 1]] if len(walls) else np.zeros(0, np.int64) + m = rl != rr + nbr = defaultdict(lambda: defaultdict(int)) + for a, b in zip(rl[m].tolist(), rr[m].tolist()): + nbr[a][b] += 1 + nbr[b][a] += 1 + + # union-find merge of every small region into its strongest neighbour + parent = list(range(n + 1)) + + def find(x): + r = x + while parent[r] != r: + r = parent[r] + while parent[x] != r: + parent[x], x = r, parent[x] + return r + + merged = False + for r in sorted(small, key=lambda r: fcount[r]): + if not nbr[r]: + continue # isolated region with no neighbour: leave it + target = max(nbr[r].items(), key=lambda kv: (kv[1], fcount[kv[0]]))[0] + ra, rt = find(r), find(target) + if ra != rt: + parent[ra] = rt + merged = True + if not merged: + break + labels = np.array([find(int(l)) for l in labels], dtype=np.int64) + + return relabel(labels) + + +class SegmentPatchesNode(io.ComfyNode): + """Flood-fill a mesh into CADable patches, merge speckle, output region_id.""" + + @classmethod + def define_schema(cls): + return io.Schema( + node_id="GeomPackSegmentPatches", + display_name="Segment Mesh into Patches", + category="geompack/analysis", + description=( + "Segment a mesh into CADable surface patches and output them as a per-face " + "region_id field (1..R). Feature edges = dihedral >= crease_angle_deg cut the " + "mesh; faces flood-fill into patches; then small/speckle patches are merged into " + "their strongest neighbour so you get a clean set instead of fragments.\n" + "\n" + "wide_rings broadens the normals first (chamfers / fine fillets become one patch " + "border instead of many tiny ones). min_patch_faces / min_patch_area_pct control " + "how aggressively speckle is absorbed.\n" + "\n" + "Color the output in Preview Mesh (fields -> region_id), split it per-patch with " + "Split By Field, or feed it to primitive fitting." + ), + is_output_node=True, + inputs=[ + io.Custom("TRIMESH").Input("trimesh"), + io.Float.Input("crease_angle_deg", default=30.0, min=0.0, max=180.0, step=1.0, + tooltip="Dihedral angle (deg) above which an edge cuts patches apart."), + io.Int.Input("wide_rings", default=0, min=0, max=6, step=1, + tooltip="Broaden each face's normal over this many rings before measuring the " + "dihedral, so chamfers / fine fillets read as a single border. 0 = off, " + "2-4 catches wider blends (also rounds genuinely-sharp single edges)."), + io.Int.Input("min_patch_faces", default=20, min=0, max=1000000, step=1, + tooltip="Merge any patch with fewer than this many faces into a neighbour. " + "Raise to absorb more speckle."), + io.Float.Input("min_patch_area_pct", default=0.05, min=0.0, max=50.0, step=0.01, + tooltip="Also merge any patch whose area is below this percent of the total " + "surface area. Catches thin/sliver patches that have many faces."), + io.Combo.Input("preclean", options=["true", "false"], default="true", + tooltip="Merge duplicate vertices + drop degenerate faces first (recommended; " + "needed for correct face adjacency)."), + ], + outputs=[ + io.Custom("TRIMESH").Output(display_name="mesh_with_patches"), + io.String.Output(display_name="info"), + ], + ) + + @classmethod + def execute(cls, trimesh, crease_angle_deg=30.0, wide_rings=0, + min_patch_faces=20, min_patch_area_pct=0.05, preclean="true"): + mesh = trimesh.copy() + if preclean == "true": + try: + mesh.merge_vertices() + except Exception as e: + log.debug("merge_vertices skipped: %s", e) + try: + mesh.update_faces(mesh.nondegenerate_faces()) + mesh.remove_unreferenced_vertices() + except Exception as e: + log.debug("degenerate cleanup skipped: %s", e) + + n_faces = int(len(mesh.faces)) + adj_pairs = np.asarray(getattr(mesh, "face_adjacency", np.zeros((0, 2), int))) + _, vals, used, _ = _edge_values(mesh, "face_normals", "angle", wide_dihedral=int(wide_rings)) + passing = vals >= float(crease_angle_deg) + + region_id, n_raw = _face_regions(n_faces, adj_pairs, passing) + + areas = np.asarray(mesh.area_faces, dtype=np.float64) + total_area = float(areas.sum()) or 1.0 + min_area = (float(min_patch_area_pct) / 100.0) * total_area + region_id = _merge_small_regions(region_id, adj_pairs, passing, areas, + int(min_patch_faces), min_area) + n_final = int(region_id.max()) if len(region_id) else 0 + + mesh.face_attributes["region_id"] = region_id.astype(np.float32) + mesh.metadata = (mesh.metadata.copy() if mesh.metadata else {}) + mesh.metadata["patches"] = { + "crease_angle_deg": float(crease_angle_deg), "wide_rings": int(wide_rings), + "n_patches": n_final, "n_raw": int(n_raw), + } + + # patch-size summary (top few) + counts = np.bincount(region_id.astype(np.int64)) + sizes = sorted(counts[1:].tolist(), reverse=True) + biggest = ", ".join(str(s) for s in sizes[:8]) + info = ( + f"Segment Mesh into Patches\n" + f"faces={n_faces:,} feature edges={int(passing.sum()):,} ({used})\n" + f"patches: {n_raw:,} raw -> {n_final:,} after merging speckle\n" + f" (min {int(min_patch_faces)} faces / {min_patch_area_pct:g}% area)\n" + f"largest patches (faces): {biggest}" + ) + log.info("[SegmentPatches] %s", info.replace("\n", " | ")) + return io.NodeOutput(mesh, info, ui={"text": [info]}) + + +NODE_CLASS_MAPPINGS = {"GeomPackSegmentPatches": SegmentPatchesNode} +NODE_DISPLAY_NAME_MAPPINGS = {"GeomPackSegmentPatches": "Segment Mesh into Patches"} diff --git a/nodes/smoothing/sharpen.py b/nodes/smoothing/sharpen.py index a9b15c38..6e4fde12 100644 --- a/nodes/smoothing/sharpen.py +++ b/nodes/smoothing/sharpen.py @@ -55,6 +55,7 @@ class SharpenMeshNode(io.ComfyNode): "fast_effective": "GeomPackSharpen_FastEffective", "non_iterative": "GeomPackSharpen_NonIterative", "curvature_guided": "GeomPackSharpen_CurvatureGuided", + "decrease_gaussian": "GeomPackSharpen_DecreaseGaussian", } @classmethod @@ -237,6 +238,24 @@ class SharpenMeshNode(io.ComfyNode): "more feature-preserving." )), ]), + io.DynamicCombo.Option("decrease_gaussian", [ + io.Int.Input("iterations", default=20, min=1, max=2000, step=1, tooltip=( + "Gradient-descent steps on the Gaussian-curvature energy sum K_i^2 " + "(K_i = vertex angle defect). More steps = lower |Gaussian curvature| " + "= more DEVELOPABLE (planes/cylinders/cones, K->0). Gauss-Bonnet keeps " + "the topological total fixed, so the flow flattens the bulk and the " + "anchor keeps the overall shape.")), + io.Float.Input("strength", default=0.05, min=0.001, max=1.0, step=0.001, display_mode="number", tooltip=( + "Per-step size as a fraction of average edge length (largest vertex " + "move ~ strength*edge). The foldless barrier prevents triangle " + "inversion regardless. Default 0.05.")), + io.Float.Input("anchor_weight", default=0.5, min=0.0, max=50.0, step=0.01, display_mode="number", tooltip=( + "Tikhonov lambda: how strongly vertices stay at their ORIGINAL " + "positions while K is reduced. LOWER = flatten harder (bigger shape " + "change toward developable); HIGHER = stay close to input. Default 0.5.")), + io.Combo.Input("use_gpu", options=["true", "false"], default="true", tooltip=( + "Run on CUDA (recommended). false = CPU torch.")), + ]), io.DynamicCombo.Option("curvature_guided", [ io.Combo.Input("regularizer", options=["tv", "bilateral"], default="tv", tooltip=( "tv = TOTAL-VARIATION on the curvature field (Chambolle-Pock) -> " diff --git a/nodes/smoothing/sharpen_backends/__init__.py b/nodes/smoothing/sharpen_backends/__init__.py index 7c486014..ca565e08 100644 --- a/nodes/smoothing/sharpen_backends/__init__.py +++ b/nodes/smoothing/sharpen_backends/__init__.py @@ -13,6 +13,7 @@ from .guided_normal import NODE_CLASS_MAPPINGS as GN_MAPS, NODE_DISPLAY_NAME_MAP from .fast_effective import NODE_CLASS_MAPPINGS as FE_MAPS, NODE_DISPLAY_NAME_MAPPINGS as FE_DISP from .non_iterative import NODE_CLASS_MAPPINGS as NI_MAPS, NODE_DISPLAY_NAME_MAPPINGS as NI_DISP from .curvature_guided import NODE_CLASS_MAPPINGS as CG_MAPS, NODE_DISPLAY_NAME_MAPPINGS as CG_DISP +from .decrease_gaussian import NODE_CLASS_MAPPINGS as DG_MAPS, NODE_DISPLAY_NAME_MAPPINGS as DG_DISP NODE_CLASS_MAPPINGS = {} NODE_DISPLAY_NAME_MAPPINGS = {} @@ -25,6 +26,7 @@ NODE_CLASS_MAPPINGS.update(GN_MAPS) NODE_CLASS_MAPPINGS.update(FE_MAPS) NODE_CLASS_MAPPINGS.update(NI_MAPS) NODE_CLASS_MAPPINGS.update(CG_MAPS) +NODE_CLASS_MAPPINGS.update(DG_MAPS) NODE_DISPLAY_NAME_MAPPINGS.update(TS_DISP) NODE_DISPLAY_NAME_MAPPINGS.update(UM_DISP) @@ -34,5 +36,6 @@ NODE_DISPLAY_NAME_MAPPINGS.update(GN_DISP) NODE_DISPLAY_NAME_MAPPINGS.update(FE_DISP) NODE_DISPLAY_NAME_MAPPINGS.update(NI_DISP) NODE_DISPLAY_NAME_MAPPINGS.update(CG_DISP) +NODE_DISPLAY_NAME_MAPPINGS.update(DG_DISP) __all__ = ['NODE_CLASS_MAPPINGS', 'NODE_DISPLAY_NAME_MAPPINGS'] diff --git a/nodes/smoothing/sharpen_backends/decrease_gaussian.py b/nodes/smoothing/sharpen_backends/decrease_gaussian.py new file mode 100644 index 00000000..3d88f9f0 --- /dev/null +++ b/nodes/smoothing/sharpen_backends/decrease_gaussian.py @@ -0,0 +1,184 @@ +# SPDX-License-Identifier: GPL-3.0-or-later +# Copyright (C) 2025 ComfyUI-GeometryPack Contributors + +"""Sharpen backend that literally DECREASES Gaussian curvature (toward developable). + +Discrete Gaussian curvature at a vertex is the angle defect K_i = 2*pi - sum of the +incident triangle angles. This backend does gradient-descent on the Gaussian-curvature +energy + + E = sum_i K_i^2 + lambda * ||x - x0||^2 + +so the pointwise |K| is driven toward zero -> the surface becomes locally DEVELOPABLE +(planes / cylinders / cones, K=0), which is what most CAD surfaces are. (Gauss-Bonnet +fixes the TOTAL sum_i K_i = 2*pi*chi, so the flow spreads that topological remainder +thinly and flattens the bulk; the anchor lambda keeps it near the input so the shape +isn't destroyed.) The exact analytic angle-defect gradient is used (verified vs finite +differences), and every step is passed through the same signed-area barrier as the +foldless guided_normal update, so no triangle inverts. + +Related: developable surface flows (Stein, Grinspun, Crane 2018). Here we use the simple +L2 angle-defect energy (decrease |K|); an L1 variant would instead CONCENTRATE K onto +sparse seams.""" + +import logging + +import numpy as np +import trimesh as trimesh_module +from comfy_api.latest import io + +log = logging.getLogger("geometrypack") + +_TWO_PI = 2.0 * np.pi + + +def _boundary_vertices(F, n): + """Indices of vertices on a boundary edge (edge used by exactly one face).""" + E = np.sort(np.vstack([F[:, [0, 1]], F[:, [1, 2]], F[:, [2, 0]]]), axis=1) + uniq, cnt = np.unique(E, axis=0, return_counts=True) + be = uniq[cnt == 1] + bv = np.zeros(n, dtype=bool) + if len(be): + bv[be.ravel()] = True + return bv + + +def _decrease_gaussian_gpu(mesh, iterations, strength, anchor_weight, use_gpu): + import torch + from .guided_normal_gpu import _anti_flip_step_torch, _torch_device + + V0 = np.asarray(mesh.vertices, dtype=np.float64) + Ff = np.asarray(mesh.faces, dtype=np.int64) + if len(V0) == 0 or len(Ff) == 0: + return None, "Empty mesh (no vertices or faces).", "cpu" + + n = len(V0) + dev = _torch_device() if use_gpu else __import__("torch").device("cpu") + eps = 1e-12 + bv = _boundary_vertices(Ff, n) + + V = torch.as_tensor(V0, dtype=torch.float32, device=dev) + X0 = V.clone() + F = torch.as_tensor(Ff, dtype=torch.long, device=dev) + i0, i1, i2 = F[:, 0], F[:, 1], F[:, 2] + bmask = torch.as_tensor(bv, device=dev) + mean_edge = float(np.mean(np.linalg.norm(V0[Ff[:, 1]] - V0[Ff[:, 0]], axis=1))) + lam = float(anchor_weight) + + def ang_grad(p0, p1, p2): + """Angle at p0 between edges p0->p1, p0->p2, and its gradient w.r.t. p0,p1,p2.""" + e1 = p1 - p0 + e2 = p2 - p0 + N = torch.cross(e1, e2, dim=1) + N = N / (N.norm(dim=1, keepdim=True) + eps) + u = e1 / (e1.norm(dim=1, keepdim=True) + eps) + w = e2 / (e2.norm(dim=1, keepdim=True) + eps) + ang = torch.arccos(torch.clamp((u * w).sum(1), -1.0 + 1e-7, 1.0 - 1e-7)) + g1 = torch.cross(N, p0 - p1, dim=1) / ((p0 - p1).pow(2).sum(1, keepdim=True) + eps) + g2 = torch.cross(N, p2 - p0, dim=1) / ((p2 - p0).pow(2).sum(1, keepdim=True) + eps) + g0 = -(g1 + g2) + return ang, g0, g1, g2 + + for _ in range(int(iterations)): + a, b, c = V[i0], V[i1], V[i2] + al, gaa, gab, gac = ang_grad(a, b, c) # angle at vertex i0 + be, gbb, gbc, gba = ang_grad(b, c, a) # angle at vertex i1 (p0=b,p1=c,p2=a) + ga, gcc, gca, gcb = ang_grad(c, a, b) # angle at vertex i2 + + K = torch.full((n,), _TWO_PI, device=dev) + K.index_add_(0, i0, -al) + K.index_add_(0, i1, -be) + K.index_add_(0, i2, -ga) + K[bmask] = 0.0 # boundary defect isn't curvature + Ka, Kb, Kc = K[i0], K[i1], K[i2] + + # dE/dx = sum_v 2 K_v * grad K_v, grad K_v = -sum_{angle at v} grad(angle) + grad = torch.zeros((n, 3), device=dev) + ca = (2.0 * Ka).unsqueeze(1) + grad.index_add_(0, i0, ca * (-gaa)); grad.index_add_(0, i1, ca * (-gab)); grad.index_add_(0, i2, ca * (-gac)) + cb = (2.0 * Kb).unsqueeze(1) + grad.index_add_(0, i1, cb * (-gbb)); grad.index_add_(0, i2, cb * (-gbc)); grad.index_add_(0, i0, cb * (-gba)) + cc = (2.0 * Kc).unsqueeze(1) + grad.index_add_(0, i2, cc * (-gcc)); grad.index_add_(0, i0, cc * (-gca)); grad.index_add_(0, i1, cc * (-gcb)) + + grad = grad + 2.0 * lam * (V - X0) # Tikhonov anchor (keeps shape / gauge) + grad[bmask] = 0.0 + + # descent step: cap the LARGEST vertex move to strength*edge (conservative + # but guarantees monotone descent on this stiff non-convex energy; median- + # based stepping overshoots and drives |K| back up). Barrier => no fold. + gmax = grad.norm(dim=1).max() + eta = strength * mean_edge / (gmax + eps) + move = -eta * grad + V = V + _anti_flip_step_torch(V, F, move) # foldless barrier + + out = trimesh_module.Trimesh(vertices=V.detach().cpu().numpy().astype(np.float64), + faces=Ff.astype(np.int32), process=False) + return out, "", str(dev) + + +class SharpenDecreaseGaussianNode(io.ComfyNode): + """Decrease Gaussian curvature (developable flow) sharpening backend.""" + + @classmethod + def define_schema(cls): + return io.Schema( + node_id="GeomPackSharpen_DecreaseGaussian", + display_name="Sharpen Decrease Gaussian (backend)", + category="geompack/smoothing", + is_dev_only=True, + is_output_node=True, + inputs=[ + io.Custom("TRIMESH").Input("trimesh"), + io.Int.Input("iterations", default=20, min=1, max=2000, step=1, tooltip=( + "Gradient-descent steps on the Gaussian-curvature energy sum K_i^2 " + "(K_i = vertex angle defect). More steps = lower |Gaussian curvature| = " + "more developable (flatter / more cylinder-and-cone-like).")), + io.Float.Input("strength", default=0.05, min=0.001, max=1.0, step=0.001, display_mode="number", tooltip=( + "Per-step size, as a fraction of the average edge length (the largest " + "vertex move each step is ~strength * edge). Higher = faster but coarser; " + "the foldless barrier still prevents triangle inversion. Default 0.05.")), + io.Float.Input("anchor_weight", default=0.5, min=0.0, max=50.0, step=0.01, display_mode="number", tooltip=( + "Tikhonov lambda -- how strongly vertices are held to their ORIGINAL " + "positions while Gaussian curvature is reduced. LOWER = flatten harder " + "(bigger shape change); HIGHER = stay close to input. Default 0.5.")), + io.Combo.Input("use_gpu", options=["true", "false"], default="true", tooltip=( + "Run on CUDA (recommended). false = CPU torch.")), + ], + outputs=[ + io.Custom("TRIMESH").Output(display_name="sharpened_mesh"), + io.String.Output(display_name="info"), + ], + ) + + @classmethod + def execute(cls, trimesh, iterations=20, strength=0.05, anchor_weight=0.5, use_gpu="true"): + import time + iv, ifc = len(trimesh.vertices), len(trimesh.faces) + log.info("Backend: decrease_gaussian | %d verts %d faces | iters=%d strength=%.3f anchor=%.3f gpu=%s", + iv, ifc, iterations, strength, anchor_weight, use_gpu) + t0 = time.perf_counter() + sharpened, error, dev = _decrease_gaussian_gpu( + trimesh, iterations, strength, anchor_weight, use_gpu == "true") + elapsed = time.perf_counter() - t0 + if sharpened is None: + raise ValueError(f"Sharpening failed (decrease_gaussian): {error}") + + if hasattr(trimesh, "metadata") and trimesh.metadata: + sharpened.metadata = trimesh.metadata.copy() + sharpened.metadata["sharpening"] = { + "algorithm": "decrease_gaussian", "device": str(dev), + "iterations": iterations, "strength": strength, "anchor_weight": anchor_weight, + } + disp = np.linalg.norm(np.asarray(sharpened.vertices) - np.asarray(trimesh.vertices), axis=1) + info = (f"Sharpen Mesh Results (decrease_gaussian, device={dev}):\n\n" + f"Iterations: {iterations}\nStrength: {strength}\nAnchor weight: {anchor_weight}\n" + f"Time: {elapsed:.2f}s\n\nVertices: {iv:,} (unchanged)\nFaces: {ifc:,} (unchanged)\n\n" + f"Displacement:\n Average: {float(np.mean(disp)):.6f}\n Maximum: {float(np.max(disp)):.6f}\n") + if len(sharpened.vertices) == len(trimesh.vertices): + sharpened.vertex_attributes["sharpen_displacement_magnitude"] = disp.astype(np.float32) + return io.NodeOutput(sharpened, info, ui={"text": [info]}) + + +NODE_CLASS_MAPPINGS = {"GeomPackSharpen_DecreaseGaussian": SharpenDecreaseGaussianNode} +NODE_DISPLAY_NAME_MAPPINGS = {"GeomPackSharpen_DecreaseGaussian": "Sharpen Decrease Gaussian (backend)"}