diff --git a/nodes/analysis/__init__.py b/nodes/analysis/__init__.py index d2ef3334..c18da7a7 100644 --- a/nodes/analysis/__init__.py +++ b/nodes/analysis/__init__.py @@ -23,6 +23,8 @@ from .ultimate_mesh_inspection import NODE_CLASS_MAPPINGS as ULTIMATE_INSPECT_MA from .ultimate_mesh_inspection import NODE_DISPLAY_NAME_MAPPINGS as ULTIMATE_INSPECT_DISPLAY from .compute_curvature import NODE_CLASS_MAPPINGS as COMPUTE_CURVATURE_MAPPINGS from .compute_curvature import NODE_DISPLAY_NAME_MAPPINGS as COMPUTE_CURVATURE_DISPLAY +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 # Combine all mappings NODE_CLASS_MAPPINGS = { @@ -36,6 +38,7 @@ NODE_CLASS_MAPPINGS = { **EXTRACT_LARGEST_MAPPINGS, **ULTIMATE_INSPECT_MAPPINGS, **COMPUTE_CURVATURE_MAPPINGS, + **FILLET_FIELD_MAPPINGS, } NODE_DISPLAY_NAME_MAPPINGS = { @@ -49,6 +52,7 @@ NODE_DISPLAY_NAME_MAPPINGS = { **EXTRACT_LARGEST_DISPLAY, **ULTIMATE_INSPECT_DISPLAY, **COMPUTE_CURVATURE_DISPLAY, + **FILLET_FIELD_DISPLAY, } __all__ = ['NODE_CLASS_MAPPINGS', 'NODE_DISPLAY_NAME_MAPPINGS'] diff --git a/nodes/analysis/compute_fillet_field.py b/nodes/analysis/compute_fillet_field.py new file mode 100644 index 00000000..9e3de176 --- /dev/null +++ b/nodes/analysis/compute_fillet_field.py @@ -0,0 +1,204 @@ +# SPDX-License-Identifier: GPL-3.0-or-later +# Copyright (C) 2025 ComfyUI-GeometryPack Contributors + +"""Compute a noise-robust "fillet field": where a surface is genuinely, smoothly +curved (rolling-ball fillets, cylinders, blends) as opposed to flat-but-noisy or +sharp-edged. + +The discriminator a fillet needs is NOT raw per-edge curvature -- noise lives at +that scale too, so a noisy flat and a noisy fillet look identical edge-to-edge. +The fix is SCALE: estimate curvature with a least-squares quadric fit over a +k-ring neighborhood (libigl `principal_curvature`). The fit averages out +high-frequency noise (random -> cancels) while the fillet's steady, coherent bend +survives (consistent -> accumulates). So: + + * noisy flat -> fitted curvature ~ 0 -> fillet_field ~ 0 (let it smooth) + * noisy fillet -> fitted curvature ~ 1/r > 0 -> fillet_field high (protect) + * sharp edge -> handled by crease angle elsewhere + +The tightest principal curvature |k|max is the rolling-ball curvature: its +reciprocal `fillet_radius = 1/|k|max` IS the local rolling-ball radius (in model +units). Multiplying by the bounding-box diagonal makes it dimensionless / +scale-free (the "global" threshold value): `fillet_field = |k|max * bbox_diag`, +so the same part at 1 mm or 1 km reads identically. + +Fields attached (as vertex_attributes, viewable / thresholdable downstream): + fillet_field = |k|max * bbox_diag (dimensionless, primary; threshold this) + fillet_radius = 1/|k|max (model units; the local rolling-ball radius) + curvedness = sqrt((k1^2+k2^2)/2) (raw magnitude, 1/length) + is_fillet = 1 where fillet_radius <= fillet_radius_frac * bbox_diag else 0 + (the ready-to-use "protect these faces" mask) +""" + +import logging + +import numpy as np +import trimesh as trimesh_module +from comfy_api.latest import io + +log = logging.getLogger("geometrypack") + + +class ComputeFilletFieldNode(io.ComfyNode): + """Noise-robust detector of genuinely-curved (fillet/blend) regions.""" + + @classmethod + def define_schema(cls): + return io.Schema( + node_id="GeomPackComputeFilletField", + display_name="Compute Fillet Field", + category="geompack/analysis", + description=( + "Find where a surface is genuinely smoothly curved (rolling-ball fillets, " + "blends, cylinders) vs flat-but-noisy -- robust to noise ON the fillet.\n" + "\n" + "Uses a least-squares quadric fit (libigl principal_curvature) over a k-ring " + "set by 'radius': the fit averages out high-frequency noise while the fillet's " + "steady bend survives. The tightest principal curvature is the rolling-ball " + "curvature, so fillet_radius = 1/|k|max is the local blend radius.\n" + "\n" + "Outputs scalar vertex fields: fillet_field (dimensionless = |k|max*bbox_diag, " + "the scale-free value to threshold), fillet_radius (model units), curvedness, " + "and is_fillet (a 0/1 protect mask for curves up to fillet_radius_frac of the " + "object). Feed the mask/field to Preview Mesh Boundaries to eyeball it, or use " + "it to exempt fillet regions from guided-normal sharpening.\n" + "\n" + "radius: bigger = more noise-robust but blurs tighter fillets (needs noise_scale " + "< radius < fillet_width to work). fillet_radius_frac: protect curves whose " + "rolling-ball radius is up to this fraction of the bounding box." + ), + is_output_node=True, + inputs=[ + io.Custom("TRIMESH").Input("trimesh"), + io.Int.Input("radius", default=6, min=1, max=16, step=1, tooltip=( + "Quadric-fit neighborhood (k-ring, in avg edge lengths). This is the SCALE " + "knob: it must be larger than the noise wavelength (so noise averages out) " + "but smaller than the fillet width (so it isn't blurred across). ~6 default; " + "raise for noisy scans, lower for crisp CAD with tight fillets.")), + io.Float.Input("fillet_radius_frac", default=0.1, min=0.0001, max=2.0, step=0.01, tooltip=( + "is_fillet mask cutoff: flag a vertex as a fillet/curve to protect when its " + "local rolling-ball radius is <= this fraction of the bounding-box diagonal. " + "0.1 = protect anything curving tighter than 10% of the object. Scale-free.")), + io.Int.Input("smoothing_iterations", default=2, min=0, max=50, step=1, tooltip=( + "Extra cotangent-Laplacian diffusion of the field after fitting. A couple of " + "passes clean up residual speckle on the mask without re-fitting. 0 = off.")), + io.Float.Input("clamp_percentile", default=1.0, min=0.0, max=20.0, step=0.5, tooltip=( + "Clip the curvature fields to their [p, 100-p] percentile so sliver triangles " + "don't blow out the range / the bbox-normalization. 0 = off.")), + io.Combo.Input("preclean", options=["true", "false"], default="true", tooltip=( + "Merge duplicate vertices, drop degenerate faces, fix normals before fitting " + "(recommended -- degenerate triangles are the main quadric-fit failure mode).")), + ], + outputs=[ + io.Custom("TRIMESH").Output(display_name="mesh_with_fillet_field"), + io.String.Output(display_name="info"), + ], + ) + + @classmethod + def execute(cls, trimesh, radius=6, fillet_radius_frac=0.1, + smoothing_iterations=2, clamp_percentile=1.0, preclean="true"): + import igl + import scipy.sparse as sp + + 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-face cleanup skipped: %s", e) + try: + trimesh_module.repair.fix_normals(mesh) + except Exception as e: + log.debug("fix_normals skipped: %s", e) + + V = np.ascontiguousarray(mesh.vertices, dtype=np.float64) + F = np.ascontiguousarray(mesh.faces, dtype=np.int64) + n = len(V) + + # Scale reference: bounding-box diagonal makes the field dimensionless. + bbox = V.max(axis=0) - V.min(axis=0) + diag = float(np.linalg.norm(bbox)) + if diag <= 0: + diag = 1.0 + log.info("Compute Fillet Field: %d verts %d faces, radius=%d, bbox_diag=%.6g", + n, len(F), radius, diag) + + # Noise-robust principal curvatures via least-squares quadric fit. + out = igl.principal_curvature(V, F, int(max(1, radius))) + k1 = np.asarray(out[2], dtype=np.float64) # max principal + k2 = np.asarray(out[3], dtype=np.float64) # min principal + + curvedness = np.sqrt(0.5 * (k1 * k1 + k2 * k2)) + kmax_abs = np.maximum(np.abs(k1), np.abs(k2)) # rolling-ball curvature + + fields = { + "curvedness": curvedness, + "fillet_field": kmax_abs * diag, # dimensionless, primary + } + + # Optional Laplacian diffusion (denoise the field, not the geometry). + if smoothing_iterations > 0: + M = igl.massmatrix(V, F, igl.MASSMATRIX_TYPE_VORONOI) + Mdiag = np.asarray(M.diagonal(), dtype=np.float64) + Mdiag[Mdiag <= 0] = 1e-12 + Minv = sp.diags(1.0 / Mdiag) + L = igl.cotmatrix(V, F) + lam = 0.5 + for name in list(fields.keys()): + f = np.asarray(fields[name], dtype=np.float64).copy() + for _ in range(int(smoothing_iterations)): + f = f + lam * np.asarray(Minv @ (L @ f)) + fields[name] = f + + # Percentile clamp so slivers don't dominate. + if clamp_percentile and clamp_percentile > 0: + p = float(clamp_percentile) + for name in list(fields.keys()): + lo, hi = np.percentile(fields[name], [p, 100.0 - p]) + fields[name] = np.clip(fields[name], lo, hi) + + # Derived radius + protect mask, from the (possibly smoothed) fillet_field. + ff = np.asarray(fields["fillet_field"], dtype=np.float64) + fillet_radius = diag / np.maximum(ff, 1e-12) # = 1/|k|max in model units + is_fillet = (fillet_radius <= float(fillet_radius_frac) * diag).astype(np.float32) + fields["fillet_radius"] = fillet_radius + fields["is_fillet"] = is_fillet + + for name, f in fields.items(): + mesh.vertex_attributes[name] = np.ascontiguousarray(f, dtype=np.float32) + + mesh.metadata = (mesh.metadata.copy() if mesh.metadata else {}) + mesh.metadata["fillet_field"] = { + "radius": int(radius), "bbox_diag": diag, + "fillet_radius_frac": float(fillet_radius_frac), + "smoothing_iterations": int(smoothing_iterations), + } + + n_fillet = int(is_fillet.sum()) + pct = 100.0 * n_fillet / max(1, n) + cutoff_r = float(fillet_radius_frac) * diag + lines = [ + "Compute Fillet Field (quadric fit, radius=%d)" % int(radius), + f"verts={n:,} faces={len(F):,} | bbox_diag={diag:.6g}", + f"protect cutoff: rolling-ball radius <= {cutoff_r:.6g} " + f"({fillet_radius_frac:g} x bbox) -> fillet_field >= {1.0/float(fillet_radius_frac):.3g}", + f"flagged as fillet/curve: {n_fillet:,} verts ({pct:.1f}%)", + "", + f"fillet_field min={float(ff.min()):.4g} max={float(ff.max()):.4g} mean={float(ff.mean()):.4g}", + f"fillet_radius min={float(fillet_radius.min()):.4g} max={float(fillet_radius.max()):.4g} (model units)", + f"curvedness min={float(fields['curvedness'].min()):.4g} max={float(fields['curvedness'].max()):.4g} (1/length)", + ] + info = "\n".join(lines) + log.info("Compute Fillet Field done: %s", " | ".join(lines[1:4])) + + return io.NodeOutput(mesh, info, ui={"text": [info]}) + + +NODE_CLASS_MAPPINGS = {"GeomPackComputeFilletField": ComputeFilletFieldNode} +NODE_DISPLAY_NAME_MAPPINGS = {"GeomPackComputeFilletField": "Compute Fillet Field"} diff --git a/nodes/analysis/ultimate_mesh_inspection.py b/nodes/analysis/ultimate_mesh_inspection.py index c81c90ad..99c4a26f 100644 --- a/nodes/analysis/ultimate_mesh_inspection.py +++ b/nodes/analysis/ultimate_mesh_inspection.py @@ -29,12 +29,24 @@ metadata['file_path']) without processing for those Tier-0 counts. import json import logging +import os +import uuid import numpy as np import trimesh as trimesh_module from comfy_api.latest import io +# vertex (point) fields vs face (cell) fields -- used to qualify field names for the +# vtk.js viewer, which prefixes point data with "point:" and cell data with "cell:". +_VERTEX_FIELDS = {"boundary_vertex", "nonmanifold_vertex", "valence"} + + +def _qualify(field): + if not field: + return None + return ("point:" if field in _VERTEX_FIELDS else "cell:") + field + log = logging.getLogger("geometrypack") RAD2DEG = 180.0 / np.pi @@ -66,6 +78,31 @@ def _nonmanifold_and_boundary_edges(mesh): return nm_edges, bnd_edges +def _winding_inconsistent_faces(mesh): + """Per-face mask of faces touching an inconsistently-wound manifold edge. + + For two faces sharing an undirected edge {a,b}, consistent orientation means one + traverses a->b and the other b->a (opposite directed edges). If BOTH traverse it + the same way, the winding flips across that edge -> both faces flagged. + """ + from trimesh.grouping import group_rows + try: + F = len(mesh.faces) + es = mesh.edges_sorted # (3F,2) undirected + ed = mesh.edges # (3F,2) directed (face order) + ef = mesh.edges_face # (3F,) face index per directed edge + pairs = group_rows(es, require_count=2) # (k,2) row-index pairs of manifold edges + mask = np.zeros(F, dtype=bool) + if len(pairs): + same_dir = np.all(ed[pairs[:, 0]] == ed[pairs[:, 1]], axis=1) # True = flipped + bad = pairs[same_dir] + if len(bad): + mask[np.unique(ef[bad.reshape(-1)])] = True + return mask + except Exception: + return None + + def _boundary_loops(bnd_edges): """Count boundary loops (connected components of the boundary-edge graph) via union-find.""" if len(bnd_edges) == 0: @@ -151,8 +188,9 @@ class UltimateMeshInspection(io.ComfyNode): "duplicate/degenerate counts reflect the file, not trimesh's auto-clean."), io.Boolean.Input("check_nonmanifold_vertices", default=True, optional=True, tooltip="Detect bowtie / non-manifold vertices (per-vertex; skipped above ~400k faces)."), - io.Boolean.Input("check_self_intersection", default=False, optional=True, - tooltip="Detect self-intersecting face pairs (expensive; needs pymeshlab)."), + io.Boolean.Input("check_self_intersection", default=True, optional=True, + tooltip="Detect self-intersecting faces (uses pymeshlab; can be slow on huge " + "meshes -- turn off for speed). Lights up the offending faces."), ], outputs=[ io.Custom("TRIMESH").Output(display_name="mesh"), @@ -182,8 +220,8 @@ class UltimateMeshInspection(io.ComfyNode): S = {} # structured stats rows = [] # (group, label, value, field_or_None) - def row(group, label, value, field=None): - rows.append((group, label, value, field)) + def row(group, label, value, field=None, focus=None): + rows.append((group, label, value, field, focus)) # ---- Tier 0: numerics / sanity ---- nonfinite = int(np.count_nonzero(~np.isfinite(verts))) @@ -204,7 +242,11 @@ class UltimateMeshInspection(io.ComfyNode): if F == 0: report = cls._format(rows, raw_note) - return io.NodeOutput(mesh, report, json.dumps(S)) + ui_rows = [{"group": g, "label": l, "value": str(v), "field": _qualify(f), "focus": fo} + for (g, l, v, f, fo) in rows] + return io.NodeOutput(mesh, report, json.dumps(S), + ui={"mesh_file": [], "report": [ui_rows], + "vertex_count": [V], "face_count": [0], "is_watertight": [False]}) faces = np.asarray(mesh.faces) @@ -240,12 +282,20 @@ class UltimateMeshInspection(io.ComfyNode): is_sliver = (face_min_ang < 10.0) & ~is_degenerate is_needle = (aspect > 20.0) & ~is_degenerate is_cap = (face_max_ang > 150.0) & ~is_degenerate + fcent = tv.mean(1) # (F,3) face centroids (for camera focus on a worst face) + worst_min_face = int(np.argmin(face_min_ang)) worst_q_face = int(np.argmin(Q)) + worst_aspect_face = int(np.argmax(aspect)) + worst_area_face = int(np.argmin(areas)) + is_worst_min_angle = np.zeros(F, np.float32); is_worst_min_angle[worst_min_face] = 1.0 is_worst_quality = np.zeros(F, np.float32); is_worst_quality[worst_q_face] = 1.0 + is_worst_aspect = np.zeros(F, np.float32); is_worst_aspect[worst_aspect_face] = 1.0 + is_smallest_area = np.zeros(F, np.float32); is_smallest_area[worst_area_face] = 1.0 amn, amx, ame, ast = _stats(face_min_ang) qmn, qmx, qme, qst = _stats(Q) armn, armx, arme, arst = _stats(aspect) + aa_mn, aa_mx, aa_me, aa_st = _stats(areas) pct_lt30 = float(100.0 * np.count_nonzero(face_min_ang < 30.0) / F) pct_q_lt01 = float(100.0 * np.count_nonzero(Q < 0.1) / F) S["element_quality"] = { @@ -255,12 +305,20 @@ class UltimateMeshInspection(io.ComfyNode): "degenerate": int(is_degenerate.sum()), "sliver": int(is_sliver.sum()), "needle": int(is_needle.sum()), "cap": int(is_cap.sum()), } - row("Element quality", "min angle (deg)", f"worst {amn:.2f} | mean {ame:.2f} +/- {ast:.2f}", "min_angle_deg") + row("Element quality", "min angle - worst", f"{amn:.2f} deg (face {worst_min_face})", + "is_worst_min_angle", fcent[worst_min_face].tolist()) + row("Element quality", "min angle - mean", f"{ame:.2f} +/- {ast:.2f} deg", "min_angle_deg") row("Element quality", "% faces min-angle < 30", f"{pct_lt30:.2f}%", "min_angle_deg") - row("Element quality", "triangle Q (1=ideal)", f"worst {qmn:.3f} | mean {qme:.3f} +/- {qst:.3f}", "tri_quality") + row("Element quality", "triangle Q - worst", f"{qmn:.3f} (face {worst_q_face})", + "is_worst_quality", fcent[worst_q_face].tolist()) + row("Element quality", "triangle Q - mean", f"{qme:.3f} +/- {qst:.3f} (1=ideal)", "tri_quality") row("Element quality", "% faces Q < 0.1", f"{pct_q_lt01:.2f}%", "tri_quality") - row("Element quality", "aspect ratio", f"worst {armx:.2f} | mean {arme:.2f} +/- {arst:.2f}", "aspect_ratio") - row("Element quality", "worst triangle (by Q)", f"face {worst_q_face} (Q={Q[worst_q_face]:.3g})", "is_worst_quality") + row("Element quality", "aspect ratio - worst", f"{armx:.2f} (face {worst_aspect_face})", + "is_worst_aspect", fcent[worst_aspect_face].tolist()) + row("Element quality", "aspect ratio - mean", f"{arme:.2f} +/- {arst:.2f}", "aspect_ratio") + row("Element quality", "area - smallest", f"{aa_mn:.3g} (face {worst_area_face})", + "is_smallest_area", fcent[worst_area_face].tolist()) + row("Element quality", "area - mean", f"{aa_me:.3g} +/- {aa_st:.3g}", "face_area") row("Degenerate", "zero-area", int(is_degenerate.sum()), "is_degenerate") row("Degenerate", "slivers (min-ang<10)", int(is_sliver.sum()), "is_sliver") row("Degenerate", "needles (aspect>20)", int(is_needle.sum()), "is_needle") @@ -285,6 +343,7 @@ class UltimateMeshInspection(io.ComfyNode): winding_ok = bool(mesh.is_winding_consistent) except Exception: winding_ok = None + wind_mask = _winding_inconsistent_faces(mesh) try: watertight = bool(mesh.is_watertight) except Exception: @@ -312,7 +371,13 @@ class UltimateMeshInspection(io.ComfyNode): "boundary_edges": n_boundary, "boundary_loops": n_loops, "nonmanifold_edges": n_nm_edges, "components": n_comp} row("Topology", "watertight", watertight, "touches_boundary") - row("Topology", "winding consistent", winding_ok) + if wind_mask is not None and wind_mask.any(): + nwf = int(wind_mask.sum()) + row("Topology", "winding consistent", f"{winding_ok} ({nwf} flipped faces)", + "winding_inconsistent", fcent[int(np.argmax(wind_mask))].tolist()) + else: + row("Topology", "winding consistent", winding_ok, + "winding_inconsistent" if wind_mask is not None else None) row("Topology", "boundary edges", n_boundary, "touches_boundary") row("Topology", "boundary loops (holes)", n_loops, "touches_boundary") if loop_len_note: @@ -344,7 +409,12 @@ class UltimateMeshInspection(io.ComfyNode): irregular = int(np.count_nonzero((valence != 6) & interior)) else: pct_v6, irregular = 0.0, 0 - S["regularity"] = {"pct_interior_valence6": pct_v6, "irregular_interior_vertices": irregular} + max_val_vert = int(np.argmax(valence)) if V else 0 + max_val = int(valence.max()) if V else 0 + S["regularity"] = {"pct_interior_valence6": pct_v6, "irregular_interior_vertices": irregular, + "max_valence": max_val} + row("Regularity", "max vertex valence", f"{max_val} (vertex {max_val_vert})", + "valence", verts[max_val_vert].tolist() if V else None) row("Regularity", "% interior verts valence 6", f"{pct_v6:.1f}%", "valence") row("Regularity", "irregular interior verts", irregular, "valence") @@ -387,11 +457,17 @@ class UltimateMeshInspection(io.ComfyNode): row("Topology", "non-manifold (bowtie) verts", int(nmv.sum()), "nonmanifold_vertex") S["topology"]["nonmanifold_vertices"] = int(nmv.sum()) - # ---- Tier 4: self-intersection (gated) ---- + # ---- Tier 4: self-intersection (gated; default on) ---- + si_mask = None if check_self_intersection: - n_si = cls._self_intersections(mesh) - row("Geometry", "self-intersecting face pairs", n_si if n_si is not None else "n/a (need pymeshlab)") - S["geometry"] = {"self_intersecting_pairs": n_si} + si_mask = cls._self_intersections(mesh) # (F,) bool or None + if si_mask is None: + row("Geometry", "self-intersecting faces", "n/a (pymeshlab unavailable)") + else: + n_si = int(si_mask.sum()) + focus = fcent[int(np.argmax(si_mask))].tolist() if n_si else None + row("Geometry", "self-intersecting faces", n_si, "is_self_intersecting", focus) + S["geometry"] = {"self_intersecting_faces": n_si} # ---- bake highlight fields onto the mesh ---- try: @@ -424,6 +500,14 @@ class UltimateMeshInspection(io.ComfyNode): mesh.face_attributes["is_needle"] = is_needle.astype(np.float32) mesh.face_attributes["is_cap"] = is_cap.astype(np.float32) mesh.face_attributes["is_worst_quality"] = is_worst_quality + mesh.face_attributes["is_worst_min_angle"] = is_worst_min_angle + mesh.face_attributes["is_worst_aspect"] = is_worst_aspect + mesh.face_attributes["face_area"] = areas.astype(np.float32) + mesh.face_attributes["is_smallest_area"] = is_smallest_area + mesh.face_attributes["is_self_intersecting"] = ( + si_mask.astype(np.float32) if si_mask is not None else np.zeros(F, np.float32)) + mesh.face_attributes["winding_inconsistent"] = ( + wind_mask.astype(np.float32) if wind_mask is not None else np.zeros(F, np.float32)) mesh.face_attributes["touches_boundary"] = boundary_face mesh.face_attributes["touches_nonmanifold_edge"] = nm_face mesh.face_attributes["is_inverted"] = is_inverted @@ -434,19 +518,43 @@ class UltimateMeshInspection(io.ComfyNode): log.warning("Could not bake all fields: %s", e) report = cls._format(rows, raw_note) + + # export a VTP (with all baked fields) for the embedded vtk.js viewer + mesh_file = None + try: + import folder_paths + from ..visualization._vtp_export import export_mesh_with_scalars_vtp + out_dir = folder_paths.get_output_directory() + mesh_file = f"ultimate_inspect_{uuid.uuid4().hex[:8]}.vtp" + export_mesh_with_scalars_vtp(mesh, os.path.join(out_dir, mesh_file)) + except Exception as e: + log.warning("VTP export for viewer failed: %s", e) + + ui_rows = [{"group": g, "label": l, "value": str(v), "field": _qualify(f), "focus": fo} + for (g, l, v, f, fo) in rows] + ui_data = { + "mesh_file": [mesh_file] if mesh_file else [], + "report": [ui_rows], + "vertex_count": [V], + "face_count": [F], + "is_watertight": [bool(watertight)], + } log.info("[UltimateMeshInspection] V=%d F=%d watertight=%s genus=%s", V, F, watertight, S["topology"].get("genus")) - return io.NodeOutput(mesh, report, json.dumps(S)) + return io.NodeOutput(mesh, report, json.dumps(S), ui=ui_data) @staticmethod def _self_intersections(mesh): + """Return a per-face boolean mask of self-intersecting faces (pymeshlab), or None.""" try: import pymeshlab as ml ms = ml.MeshSet() ms.add_mesh(ml.Mesh(np.asarray(mesh.vertices), np.asarray(mesh.faces))) - ms.compute_selection_by_self_intersections_per_face() - m = ms.current_mesh() - return int(np.count_nonzero(m.face_selection_array())) + try: + ms.compute_selection_by_self_intersections_per_face() + except Exception: + ms.apply_filter("compute_selection_by_self_intersections_per_face") # older API name + return np.asarray(ms.current_mesh().face_selection_array(), dtype=bool) except Exception as e: log.info("self-intersection check unavailable: %s", e) return None @@ -455,7 +563,7 @@ class UltimateMeshInspection(io.ComfyNode): def _format(rows, note): lines = [f"=== Ultimate Mesh Inspection{note} ===", ""] cur = None - for group, label, value, field in rows: + for group, label, value, field, _focus in rows: if group != cur: lines.append(f"[{group}]") cur = group diff --git a/nodes/smoothing/sharpen_backends/guided_normal_gpu.py b/nodes/smoothing/sharpen_backends/guided_normal_gpu.py index 6927937d..fd297768 100644 --- a/nodes/smoothing/sharpen_backends/guided_normal_gpu.py +++ b/nodes/smoothing/sharpen_backends/guided_normal_gpu.py @@ -86,6 +86,76 @@ def _pad_inner_edges(neighbors, adj_pairs, face_to_adj, m): return IE, IEM +def _build_topology_vectorized(F, adj, n, m, rings): + """Vectorized equivalent of the per-face Python topology build. + + Produces the SAME padded arrays as _pad_patches + _pad_inner_edges, but via + sparse boolean reachability instead of O(m*K*deg) interpreted set loops -- + which hang on large meshes at rings>=2. + + Patch(f) = faces within `rings` shared-vertex hops of f (incl. self), exactly + the BFS the Python version grows. Inner edges of a patch = adjacency edges + whose BOTH faces lie in the patch. Index order within a row may differ from + the Python version, but every downstream use is an order-independent reduce + (sum / amax / argmin over sorted indices), so results are unchanged. + + Returns (P[m,K], M[m,K] bool, K, IE[m,Emax], IEM[m,Emax] bool). + """ + import scipy.sparse as sp + + rings = max(1, int(rings)) + # face<->vertex incidence; FF = face-face "shares >=1 vertex" (incl. self) + rows = np.repeat(np.arange(m, dtype=np.int64), 3) + A = sp.csr_matrix((np.ones(3 * m, dtype=np.int8), (rows, F.ravel())), shape=(m, n)) + FF = (A @ A.T).astype(bool).tocsr() + S = FF + for _ in range(rings - 1): + S = (S @ FF).astype(bool) + S = S.tocsr() + S.sort_indices() + + # --- padded patches from S rows (column indices are sorted per row) --- + cnts = np.diff(S.indptr).astype(np.int64) + K = int(cnts.max()) if m else 1 + K = max(1, K) + P = np.zeros((m, K), dtype=np.int64) + M = np.zeros((m, K), dtype=bool) + rep = np.repeat(np.arange(m, dtype=np.int64), cnts) + pos = np.arange(S.nnz, dtype=np.int64) - np.repeat(S.indptr[:-1].astype(np.int64), cnts) + P[rep, pos] = S.indices + M[rep, pos] = True + + # --- inner edges: edge ai is inner to face i iff i in patch(fa)&patch(fb) --- + # (S symmetric -> fa in patch(i) <=> i in patch(fa)); chunk over edges to cap memory. + Scsc = S.tocsc() + fa = np.ascontiguousarray(adj[:, 0]) + fb = np.ascontiguousarray(adj[:, 1]) + E = len(adj) + step = max(1, 2_000_000 // max(1, K)) + ir_parts, ic_parts = [], [] + for s in range(0, E, step): + e = min(E, s + step) + I = Scsc[:, fa[s:e]].multiply(Scsc[:, fb[s:e]]).tocoo() + ir_parts.append(I.row.astype(np.int64)) + ic_parts.append(I.col.astype(np.int64) + s) # local col -> global edge id + ir = np.concatenate(ir_parts) if ir_parts else np.zeros(0, np.int64) + ic = np.concatenate(ic_parts) if ic_parts else np.zeros(0, np.int64) + + ecnts = np.bincount(ir, minlength=m).astype(np.int64) + Emax = int(ecnts.max()) if len(ecnts) else 1 + Emax = max(1, Emax) + IE = np.zeros((m, Emax), dtype=np.int64) + IEM = np.zeros((m, Emax), dtype=bool) + order = np.argsort(ir, kind="stable") + ir, ic = ir[order], ic[order] + starts = np.zeros(m + 1, dtype=np.int64) + np.cumsum(ecnts, out=starts[1:]) + pos2 = np.arange(len(ir), dtype=np.int64) - np.repeat(starts[:-1], ecnts) + IE[ir, pos2] = ic + IEM[ir, pos2] = True + return P, M, K, IE, IEM + + def _guided_normal_gpu(mesh, normal_iterations, vertex_iterations, sigma_s, sigma_r, neighborhood_rings=1): import torch @@ -102,19 +172,27 @@ def _guided_normal_gpu(mesh, normal_iterations, vertex_iterations, sigma_s, sigm n = len(V0) # --- one-time topology (CPU), constant across iterations --- + # Vectorized sparse build (bit-equivalent sets); the old pure-Python per-face + # set loops hang on large meshes at rings>=2. Fall back to Python on failure. t_topo = time.perf_counter() - vtf = _build_vertex_to_faces(n, Ff) - neighbors = _build_vertex_based_face_neighbors(Ff, vtf, include_central=True, rings=neighborhood_rings) - face_to_adj = [[] for _ in range(m)] - for ai in range(len(adj)): - fa, fb = adj[ai] - face_to_adj[fa].append(ai) - face_to_adj[fb].append(ai) - P_np, M_np, K = _pad_patches(neighbors, m) - IE_np, IEM_np = _pad_inner_edges(neighbors, adj, face_to_adj, m) - log.info("[guided_normal_gpu] topology build (CPU): %.3fs (m=%d faces, n=%d verts, " + try: + P_np, M_np, K, IE_np, IEM_np = _build_topology_vectorized(Ff, adj, n, m, neighborhood_rings) + topo_mode = "vectorized" + except Exception as ex: + log.warning("[guided_normal_gpu] vectorized topology failed (%s); using Python fallback", ex) + vtf = _build_vertex_to_faces(n, Ff) + neighbors = _build_vertex_based_face_neighbors(Ff, vtf, include_central=True, rings=neighborhood_rings) + face_to_adj = [[] for _ in range(m)] + for ai in range(len(adj)): + fa, fb = adj[ai] + face_to_adj[fa].append(ai) + face_to_adj[fb].append(ai) + P_np, M_np, K = _pad_patches(neighbors, m) + IE_np, IEM_np = _pad_inner_edges(neighbors, adj, face_to_adj, m) + topo_mode = "python-fallback" + log.info("[guided_normal_gpu] topology build (%s): %.3fs (m=%d faces, n=%d verts, " "K=%d, Emax=%d, rings=%d)", - time.perf_counter() - t_topo, m, n, K, IE_np.shape[1], neighborhood_rings) + topo_mode, time.perf_counter() - t_topo, m, n, K, IE_np.shape[1], neighborhood_rings) dev = _torch_device() eps = 1e-12 diff --git a/nodes/visualization/_vtp_export.py b/nodes/visualization/_vtp_export.py index 62d2e774..1628c0e1 100644 --- a/nodes/visualization/_vtp_export.py +++ b/nodes/visualization/_vtp_export.py @@ -22,6 +22,109 @@ log = logging.getLogger("geometrypack") def export_mesh_with_scalars_vtp(trimesh: trimesh_module.Trimesh, filepath: str): + """Export trimesh to a .vtp with scalar attributes, BINARY + uncompressed. + + Binary (base64) is ~3-4x smaller than ASCII and is parsed by vtk.js as a typed + array instead of parseFloat-ing millions of text numbers on the main thread -- + that text parse is what froze the browser on large meshes. It is LOSSLESS (full + float32 / int connectivity). MUST be uncompressed: the bundled vtk.js has no + zlib, so it can read ascii/binary VTP but not zlib-compressed VTP. + + Falls back to the ASCII writer if the `vtk` module is unavailable or anything + goes wrong, so exports never hard-fail. + """ + try: + _export_binary_vtp(trimesh, filepath) + except Exception as e: + log.warning("Binary VTP export failed (%s); falling back to ASCII writer", e) + _export_ascii_vtp(trimesh, filepath) + + +def _vtk_ready(name, arr): + """Contiguous, dtype-trimmed copy for VTK: 0/1 masks -> uint8, else float32.""" + is_mask = (name.startswith("is_") or name.startswith("touches_") + or name in ("boundary_vertex", "nonmanifold_vertex", "winding_inconsistent")) + dt = np.uint8 if is_mask else np.float32 + return np.ascontiguousarray(np.asarray(arr).astype(dt)) + + +def _export_binary_vtp(trimesh: trimesh_module.Trimesh, filepath: str): + """Write a binary, uncompressed .vtp via the vtk XML writer (lossless).""" + import vtk + from vtk.util import numpy_support as nps + + is_pc = is_point_cloud(trimesh) + verts = np.ascontiguousarray(np.asarray(trimesh.vertices, dtype=np.float32)) + n_verts = len(verts) + + poly = vtk.vtkPolyData() + pts = vtk.vtkPoints() + pts.SetData(nps.numpy_to_vtk(verts, deep=1)) + poly.SetPoints(pts) + + if is_pc: + offsets = np.arange(0, n_verts + 1, dtype=np.int64) + conn = np.arange(n_verts, dtype=np.int64) + ca = vtk.vtkCellArray() + ca.SetData(nps.numpy_to_vtkIdTypeArray(offsets, deep=1), + nps.numpy_to_vtkIdTypeArray(conn, deep=1)) + if n_verts < 2_000_000_000 and hasattr(ca, "ConvertTo32BitStorage"): + ca.ConvertTo32BitStorage() # int32 connectivity ~halves the cell bytes + poly.SetVerts(ca) + else: + faces = np.asarray(trimesh.faces) + n_faces = len(faces) + conn = np.ascontiguousarray(faces.reshape(-1).astype(np.int64)) + offsets = np.arange(0, 3 * (n_faces + 1), 3, dtype=np.int64) + ca = vtk.vtkCellArray() + ca.SetData(nps.numpy_to_vtkIdTypeArray(offsets, deep=1), + nps.numpy_to_vtkIdTypeArray(conn, deep=1)) + if n_verts < 2_000_000_000 and hasattr(ca, "ConvertTo32BitStorage"): + ca.ConvertTo32BitStorage() # int32 connectivity ~halves the cell bytes + poly.SetPolys(ca) + + # PointData: vertex normals (built-in) + vertex_attributes + point_data = poly.GetPointData() + if not is_pc and hasattr(trimesh, "vertex_normals"): + try: + nrm = np.asarray(trimesh.vertex_normals, dtype=np.float32) + if nrm.shape == (n_verts, 3): + a = nps.numpy_to_vtk(np.ascontiguousarray(nrm), deep=1) + a.SetName("normals") + point_data.AddArray(a) + except Exception as e: + log.warning("Could not export normals: %s", e) + if getattr(trimesh, "vertex_attributes", None): + for name, vals in trimesh.vertex_attributes.items(): + v = np.asarray(vals) + if v.ndim > 1 and v.shape[1] > 4: + continue + a = nps.numpy_to_vtk(_vtk_ready(name, v), deep=1) + a.SetName(name) + point_data.AddArray(a) + + # CellData: face_attributes (meshes only) + if not is_pc and getattr(trimesh, "face_attributes", None): + cell_data = poly.GetCellData() + for name, vals in trimesh.face_attributes.items(): + v = np.asarray(vals) + if v.ndim > 1 and v.shape[1] > 4: + continue + a = nps.numpy_to_vtk(_vtk_ready(name, v), deep=1) + a.SetName(name) + cell_data.AddArray(a) + + writer = vtk.vtkXMLPolyDataWriter() + writer.SetFileName(filepath) + writer.SetInputData(poly) + writer.SetDataModeToBinary() # base64 binary DataArrays (vtk.js reads these) + writer.SetCompressorTypeToNone() # NO zlib -- vtk.js cannot inflate + writer.Write() + log.info("Exported binary VTP: %s (%d verts, %d %s)", filepath, n_verts, + n_verts if is_pc else len(trimesh.faces), "points" if is_pc else "faces") + + +def _export_ascii_vtp(trimesh: trimesh_module.Trimesh, filepath: str): """ Export trimesh to VTK PolyData XML format (.vtp) with scalar attributes. diff --git a/nodes/visualization/preview_mesh_boundaries.py b/nodes/visualization/preview_mesh_boundaries.py index f17468ec..9d1ddbdc 100644 --- a/nodes/visualization/preview_mesh_boundaries.py +++ b/nodes/visualization/preview_mesh_boundaries.py @@ -68,6 +68,64 @@ def _edge_values(mesh, field_name, reduction): return np.asarray(adj_edges, dtype=np.int64), np.asarray(vals, dtype=np.float64), used, avail +def _edge_clusters(edges, min_edges): + """Group feature edges into connected clusters; keep the big ones, color each. + + Builds the undirected feature-edge graph, finds connected components, and + keeps every component with MORE than `min_edges` edges. No cycle / genus / + closure check -- chains, loops and branching webs are all kept as long as + they're large enough. Each kept cluster gets a distinct 1-based id (for + per-cluster coloring). + + edges: int array [M,2] of vertex-index pairs (undirected). + Returns (cluster_edges[K,2], cluster_id[K] 1-based, stats dict). + """ + import numpy as np + from collections import defaultdict + + adj = defaultdict(set) + seen = set() + for a, b in edges: + a, b = int(a), int(b) + if a == b: + continue + key = (a, b) if a < b else (b, a) + if key in seen: + continue + seen.add(key) + adj[a].add(b) + adj[b].add(a) + + visited = set() + out_edges, out_id = [], [] + n_clusters = n_dropped = 0 + + for start in list(adj.keys()): + if start in visited: + continue + comp, stack = [], [start] + visited.add(start) + while stack: + v = stack.pop() + comp.append(v) + for w in adj[v]: + if w not in visited: + visited.add(w) + stack.append(w) + comp_edges = [(v, w) for v in comp for w in adj[v] if v < w] + if len(comp_edges) <= min_edges: + n_dropped += 1 + continue + n_clusters += 1 + for e in comp_edges: + out_edges.append(e) + out_id.append(n_clusters) + + ce = np.asarray(out_edges, dtype=np.int64).reshape(-1, 2) + cid = np.asarray(out_id, dtype=np.int64) + return ce, cid, {"clusters": n_clusters, "dropped": n_dropped} + + class PreviewMeshBoundaries(io.ComfyNode): """Threshold mesh edges by an adjacent-face metric (dihedral, etc.) and preview in VTK.js.""" @@ -87,6 +145,15 @@ class PreviewMeshBoundaries(io.ComfyNode): io.Combo.Input("reduction", options=["auto", "angle", "abs_diff", "l2"], default="auto", tooltip="How to combine the two faces' values into one per-edge number. " "angle = degrees between vectors (dihedral for normals)."), + io.Combo.Input("mode", options=["edges", "loops"], default="edges", + tooltip="edges = show every passing edge. loops = group passing edges into " + "connected clusters, drop small ones (<= min_edges), and color each " + "remaining cluster a distinct color. No cycle/closure check -- chains, " + "loops and branching webs are all kept. Coincident vertices are merged " + "first so clusters connect across split hard edges."), + io.Int.Input("min_edges", default=10, min=0, max=100000, step=1, + tooltip="loops mode only: keep clusters with MORE than this many edges. " + "Smaller clusters are discarded as noise."), io.Float.Input("threshold", default=30.0, min=0.0, max=100000.0, step=1.0, tooltip="Edges pass when their value meets the threshold (dihedral in degrees)."), io.Combo.Input("comparison", options=[">=", "<="], default=">=", @@ -99,22 +166,49 @@ class PreviewMeshBoundaries(io.ComfyNode): @classmethod def execute(cls, mesh, face_field="face_normals", reduction="auto", - threshold=30.0, comparison=">=", show_surface=True): + threshold=30.0, comparison=">=", mode="edges", min_edges=10, + show_surface=True): import numpy as np import pyvista as pv import folder_paths - adj_edges, vals, used_field, avail = _edge_values(mesh, face_field, reduction) + # loops mode needs coincident vertices merged so feature edges share + # indices and can connect into closed cycles (also fixes missing + # face_adjacency across split hard edges). Work on a copy. + work = mesh + if mode == "loops": + try: + work = mesh.copy() + work.merge_vertices() + except Exception as e: + log.warning("[PreviewMeshBoundaries] merge_vertices failed: %s", e) + work = mesh + + adj_edges, vals, used_field, avail = _edge_values(work, face_field, reduction) if comparison == "<=": passing = vals <= threshold else: passing = vals >= threshold edges = adj_edges[passing] if len(adj_edges) else np.zeros((0, 2), np.int64) ev = vals[passing] if len(vals) else np.zeros(0) + + cluster_id_cell = None + cluster_stats = None + if mode == "loops": + # value lookup keyed by the undirected edge, to keep edge_value after filtering + val_map = {} + for (a, b), v in zip(edges, ev): + a, b = int(a), int(b) + val_map[(a, b) if a < b else (b, a)] = float(v) + edges, cluster_ids, cluster_stats = _edge_clusters(edges, int(min_edges)) + ev = np.array([val_map.get((int(a), int(b)) if a < b else (int(b), int(a)), 0.0) + for a, b in edges], dtype=np.float64) + cluster_id_cell = cluster_ids.astype(np.float32) + K = int(len(edges)) - points = np.asarray(mesh.vertices, dtype=np.float64) - faces = np.asarray(getattr(mesh, "faces", np.zeros((0, 3), int)), dtype=np.int64) + points = np.asarray(work.vertices, dtype=np.float64) + faces = np.asarray(getattr(work, "faces", np.zeros((0, 3), int)), dtype=np.int64) F = int(len(faces)) if show_surface else 0 # Build explicitly: pv.PolyData(points) alone would add one VERT cell per @@ -135,6 +229,12 @@ class PreviewMeshBoundaries(io.ComfyNode): boundary.append(np.zeros(F)); edge_value.append(np.zeros(F)) combined.cell_data["boundary"] = np.concatenate(boundary).astype(np.float32) combined.cell_data["edge_value"] = np.concatenate(edge_value).astype(np.float32) + if cluster_id_cell is not None: + cid = [cluster_id_cell] if K else [] + if F: + cid.append(np.zeros(F, np.float32)) + if cid: + combined.cell_data["cluster_id"] = np.concatenate(cid).astype(np.float32) tmp = folder_paths.get_temp_directory() os.makedirs(tmp, exist_ok=True) @@ -142,14 +242,28 @@ class PreviewMeshBoundaries(io.ComfyNode): # ASCII VTP -> reliably parsed by the VTK.js XMLPolyDataReader (incl. Lines). combined.save(os.path.join(tmp, filename), binary=False) - summary = (f"boundaries: {K} edge(s) {comparison} {threshold:g} via " - f"reduce={reduction} on '{used_field}'. " - f"available face fields: {avail or ['(none)']}") + fields = ["boundary", "edge_value"] + if cluster_id_cell is not None: + fields.append("cluster_id") + + if mode == "loops": + s = cluster_stats or {"clusters": 0, "dropped": 0} + summary = (f"loops: {s['clusters']} cluster(s) kept ({K} edge(s), each >{int(min_edges)} " + f"edges) at {comparison} {threshold:g} via reduce={reduction} on " + f"'{used_field}'. dropped {s['dropped']} small cluster(s). " + f"color by 'cluster_id' to tell them apart. " + f"available face fields: {avail or ['(none)']}") + if s["clusters"] == 0: + summary += (" -- no clusters this big: lower min_edges or the threshold.") + else: + summary = (f"boundaries: {K} edge(s) {comparison} {threshold:g} via " + f"reduce={reduction} on '{used_field}'. " + f"available face fields: {avail or ['(none)']}") log.info("[PreviewMeshBoundaries] %s", summary) return io.NodeOutput(ui={ "mesh_file": [filename], "boundary_edges": [K], - "field_names": [["boundary", "edge_value"]], + "field_names": [fields], "summary": [summary], })