diff --git a/.gitignore b/.gitignore index f01a138..7819436 100644 --- a/.gitignore +++ b/.gitignore @@ -54,5 +54,19 @@ v_*/ #cached python version for pyenv .python-version -# generated DXF files in examples -examples/*.dxf + + +# generated DXF files +*.dxf +# but include documentation DXF files +!docs/*.dxf + +# generated stl files +*.stl +# but include documentation stl files +!docs/*.stl + +# generated step files +*.step +# but include documentation step files +!docs/*.step diff --git a/README.md b/README.md index 194a997..3ae448b 100644 --- a/README.md +++ b/README.md @@ -79,6 +79,9 @@ calling out include: creating sectional views of assemblies. - `yapcad.io.step`/`yapcad.io.stl` – faceted exporters suitable for interchange with FreeCAD, slicers, and other simulation tools. +- `tools/validate_mesh.py` – CLI helper that runs `admesh`, `meshfix`, and an + optional slicer to gauge whether STL output is robust enough for CAM; see + `docs/mesh_validation.md` for usage. To build the HTML **yapCAD** documentation locally, install the documentation dependencies and run Sphinx from the project root: diff --git a/docs/mesh_validation.md b/docs/mesh_validation.md new file mode 100644 index 0000000..c60ffa3 --- /dev/null +++ b/docs/mesh_validation.md @@ -0,0 +1,66 @@ +# Mesh Validation Pipeline + +This repository now includes a small command–line helper that mirrors the +checks typically performed downstream of a CAD/boolean workflow. The goal is +to gauge whether the generated STL will survive real fabrication tools even if +the mesh is not mathematically watertight. + +## Prerequisites + +The script looks for the following utilities on `PATH`: + +| Tool | Purpose | macOS install hint | +|----------|--------------------------------------|--------------------| +| `admesh` | STL statistics (holes, volume, etc.) | `brew install admesh` | +| `meshfix`| Mesh repair via CGAL | `brew install meshfix` | + +If you prefer to stay in Python space, install the optional +[`pymeshfix`](https://pypi.org/project/PyMeshFix/) and +[`trimesh`](https://pypi.org/project/trimesh/) packages; the validation script +will automatically fall back to those when the stand-alone `meshfix` binary is +missing. + +You can pull these extras in with a single command using the optional-deps +group defined in `pyproject.toml`: + +```bash +pip install yapCAD[meshcheck] +``` +| `prusa-slicer` / `prusa-slicer-console` / `slic3r` | CLI slicer check | download from vendor; add binary to `PATH` | + +If a tool is missing the script will skip the corresponding step and include a +note in the report. + +## Usage + +1. Export an STL using the demo driver (or your own code): + + ```bash + PYTHONPATH=src python examples/solid_boolean_demo.py \ + --mode stl --operation difference --shapes box_hole \ + --output box_hole_difference + ``` + +2. Run the validator: + + ```bash + python tools/validate_mesh.py box_hole_difference.stl \ + --workdir build/mesh_checks + ``` + + This prints a summary to stdout and keeps any repaired STL/G-code files + under `build/mesh_checks`. + +3. (Optional) Request JSON output for CI: + + ```bash + python tools/validate_mesh.py box_hole_difference.stl --json + ``` + +The report includes highlights from `admesh`, the return code from `meshfix` +(with output file location), and the status of the slicer run. When the slicer +succeeds you’ll find a G-code file beside the repaired meshes. + +This pipeline gives a far more practical read on model quality than the strict +`issolidclosed` check, and mirrors the tools typically invoked before sending a +part to printers or CAM software. diff --git a/docs/solid_boolean_roadmap.md b/docs/solid_boolean_roadmap.md new file mode 100644 index 0000000..b32258b --- /dev/null +++ b/docs/solid_boolean_roadmap.md @@ -0,0 +1,24 @@ +# Boolean Engine Roadmap + +Status snapshot as of this session: + +- Native mesh boolean workflow now lives in `yapcad/boolean/native.py`, with `geom3d` re-exporting it for backward compatibility. +- Diagnostics still rely on `tools/validate_mesh.py`; no engine-specific metrics yet. +- External engine integration is underway; a `trimesh` backend is plumbed (requires an installed boolean operator such as Blender/OpenSCAD/Cork). + +Near-term tasks: + +- [x] **Code extraction** – current boolean implementation now resides in `yapcad/boolean/native.py`, referenced by `geom3d` wrappers. +- [x] **Engine selector UX** – `solid_boolean(..., engine=...)` now routes through the native engine by default, accepts `trimesh` (and optional `trimesh:backend`), and the demo exposes `--engine`; document env flags (`YAPCAD_BOOLEAN_ENGINE`, `YAPCAD_TRIMESH_BACKEND`) for benchmarking. +- [ ] **External prototype** – wrap at least one library boolean (e.g., `trimesh` or PyMeshFix) for STL solids, including geometry conversion helpers. +- [ ] **Benchmark harness** – update the demo CLI and validation scripts to iterate across registered engines, saving STL/lint outputs for comparison. + +Open questions: +- Which external kernel offers the best balance of licensing, install footprint, and mesh quality? +- Do we keep the native engine on equal footing (for offline/embedded use), or treat it as a fallback once a third-party backend is available? +- What minimal metrics should every engine report (shell count, single-facet edges, `validate_mesh` scores, render snapshots)? + +Once the refactor lands, future work will include: +- Automatic tolerance scaling per engine. +- Capturing validation output in machine-readable form (JSON) to feed regression dashboards. +- Evaluating hybrid workflows (native preprocessing, external boolean, native stitches/cleanup). diff --git a/examples/solid_boolean_demo.py b/examples/solid_boolean_demo.py new file mode 100644 index 0000000..4cbfb23 --- /dev/null +++ b/examples/solid_boolean_demo.py @@ -0,0 +1,103 @@ +"""Demonstration of yapCAD solid boolean operations. + +This script can visualise the result via OpenGL or export the mesh to STL/STEP. +""" + +import argparse +from pathlib import Path + +from yapcad.geom import point +from yapcad.geom3d import solid_boolean, translatesolid +from yapcad.geom3d_util import prism, sphere, conic, tube +from yapcad.io.stl import write_stl +from yapcad.io.step import write_step + + +def _build_solids(kind: str): + if kind == 'cube': + a = prism(2, 2, 2) + b = translatesolid(prism(2, 2, 2), point(0.75, 0.0, 0.0)) + elif kind == 'sphere': + a = sphere(2.0) + b = translatesolid(sphere(2.0), point(1.0, 0.0, 0.0)) + elif kind == 'box_hole': + a = prism(2, 2, 2) + b = conic(0.6, 0.6, 3.0, center=point(0.0, 0.0, -1.5)) + elif kind == 'tube_hole': + a = tube(outer_diameter=3.0, wall_thickness=0.5, length=4.0, + base_point=point(0.0, 0.0, -2.0)) + b = conic(0.6, 0.6, 4.2, center=point(1.2, 0.0, -2.1)) + else: + raise ValueError(f'unsupported shape kind: {kind!r}') + return a, b + + +def _render_opengl(solids, result): + from yapcad.pyglet_drawable import pygletDraw + + viewer = pygletDraw() + viewer.make_object('obj1', lighting=True, linecolor='white', + material='pearl', + position=point(-2.5, 2.5, 0)) + + viewer.make_object('obj2', lighting=True, linecolor='white', + material='obsidian', + position=point(2.5, 2.5, 0)) + + viewer.make_object('rslt', lighting=True, linecolor='white', + material='gold', + position=point(0, 0, 0)) + + viewer.draw_solid(solids[0], name='obj1') + viewer.draw_solid(solids[1], name='obj2') + viewer.draw_solid(result, name='rslt') + viewer.display() + + +def _export_stl(result, output): + path = Path(output) + if path.suffix.lower() != '.stl': + path = path.with_suffix('.stl') + write_stl(result, path) + print(f'Wrote STL to {path}') + + +def _export_step(result, output): + path = Path(output) + if path.suffix.lower() not in {'.step', '.stp'}: + path = path.with_suffix('.step') + write_step(result, path) + print(f'Wrote STEP to {path}') + + +def main(): + parser = argparse.ArgumentParser(description='Solid boolean demonstration.') + parser.add_argument('--mode', choices=['gl', 'stl', 'step'], default='gl', + help='visualise with OpenGL or export as STL/STEP') + parser.add_argument('--operation', choices=['union', 'intersection', 'difference'], + default='union', help='boolean operation to apply') + parser.add_argument('--shapes', choices=['cube', 'sphere', 'box_hole', 'tube_hole'], + default='cube', + help='primitive pair to use') + parser.add_argument('--output', default='boolean_result', + help='output basename for STL/STEP export') + parser.add_argument('--stitch', action='store_true', + help='experimentally stitch open edges in the result') + parser.add_argument('--engine', default=None, + help='boolean engine to use (default: native)') + args = parser.parse_args() + + solids = _build_solids(args.shapes) + result = solid_boolean(solids[0], solids[1], args.operation, + stitch=args.stitch, engine=args.engine) + + if args.mode == 'gl': + _render_opengl(solids, result) + elif args.mode == 'stl': + _export_stl(result, args.output) + else: + _export_step(result, args.output) + + +if __name__ == '__main__': + main() diff --git a/pyproject.toml b/pyproject.toml index e5af44e..51f1cb9 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -28,6 +28,10 @@ tests = [ "pytest", "pytest-cov" ] +meshcheck = [ + "trimesh>=4.0", + "pymeshfix>=0.16", +] [project.urls] "Homepage" = "https://github.com/rdevaul/yapCAD/" "Documentation" = "https://yapcad.readthedocs.io/en/latest/" diff --git a/src/yapcad/boolean/__init__.py b/src/yapcad/boolean/__init__.py new file mode 100644 index 0000000..d95bdb5 --- /dev/null +++ b/src/yapcad/boolean/__init__.py @@ -0,0 +1,21 @@ +from . import native as native + +__all__ = ['native'] + +try: + from . import trimesh_engine as trimesh +except Exception: # optional dependency + trimesh = None +else: + __all__.append('trimesh') + +ENGINE_REGISTRY = {'native': native} +if trimesh is not None: + ENGINE_REGISTRY['trimesh'] = trimesh + + +def get_engine(name: str): + return ENGINE_REGISTRY.get(name) + + +__all__.extend(['ENGINE_REGISTRY', 'get_engine']) diff --git a/src/yapcad/boolean/native.py b/src/yapcad/boolean/native.py new file mode 100644 index 0000000..deb2f43 --- /dev/null +++ b/src/yapcad/boolean/native.py @@ -0,0 +1,1012 @@ +"""Native boolean engine extracted from yapcad.geom3d.""" + +from __future__ import annotations + +import math + +from yapcad.geom import * +from yapcad.geom_util import * +from yapcad.xform import * +from yapcad.triangulator import triangulate_polygon +from yapcad.octtree import NTree + + +def _geom3d(): + from yapcad import geom3d as _g3 + return _g3 + +def _ensure_surface_metadata_dict(s): + if not isinstance(s, list) or len(s) < 4 or s[0] != 'surface': + raise ValueError('bad surface passed to _ensure_surface_metadata_dict') + if len(s) == 4: + s.extend([[], [], {}]) + return s[6] + if len(s) == 5: + s.extend([[], {}]) + return s[6] + if len(s) == 6: + s.append({}) + return s[6] + metadata = s[6] + if metadata is None: + metadata = {} + elif not isinstance(metadata, dict): + metadata = {'legacy_metadata': metadata} + s[6] = metadata + return metadata + + +def _triangle_bbox(tri, tol): + mins = [min(pt[i] for pt in tri) - tol for i in range(3)] + maxs = [max(pt[i] for pt in tri) + tol for i in range(3)] + return [point(mins[0], mins[1], mins[2]), + point(maxs[0], maxs[1], maxs[2])] + + +def _triangle_plane_intersection_points(tri, plane, tol): + n, d = plane + points = [] + for idx in range(3): + p0 = tri[idx] + p1 = tri[(idx + 1) % 3] + dir_vec = [p1[0] - p0[0], p1[1] - p0[1], p1[2] - p0[2]] + denom = n[0] * dir_vec[0] + n[1] * dir_vec[1] + n[2] * dir_vec[2] + if abs(denom) < tol: + continue + numer = d - (n[0] * p0[0] + n[1] * p0[1] + n[2] * p0[2]) + t = numer / denom + if t < -tol or t > 1.0 + tol: + continue + t_clamped = max(0.0, min(1.0, t)) + pt = _lerp_point(p0, p1, t_clamped) + if abs(n[0] * pt[0] + n[1] * pt[1] + n[2] * pt[2] - d) <= tol * 10.0: + points.append(pt) + return _unique_points(points, tol) if points else [] + + +def _candidate_planes_for_triangle(tri, target, tri_plane, tol): + box = _geom3d().solidbbox(target) + if box: + center = point((box[0][0] + box[1][0]) / 2.0, + (box[0][1] + box[1][1]) / 2.0, + (box[0][2] + box[1][2]) / 2.0) + extent = max(box[1][0] - box[0][0], + box[1][1] - box[0][1], + box[1][2] - box[0][2]) + else: + center = point(0, 0, 0) + extent = 0.0 + + epsilon_dot = max(extent * 1e-6, 1e-6) + + tri_bbox = _triangle_bbox(tri, tol * 2.0) + seen = set() + planes = [] + snap_candidates = [] + for surf in target[1]: + tree = surface_octree(surf) + if tree is None: + candidates = list(_iter_triangles_from_surface(surf)) + else: + elems = tree.getElements(tri_bbox) + candidates = [] + for elem in elems: + if isinstance(elem, list): + candidates.append(elem) + elif isinstance(elem, tuple): + candidates.append(elem[0]) + else: + candidates.append(elem) + if not candidates: + candidates = list(_iter_triangles_from_surface(surf)) + for cand in candidates: + plane = _triangle_plane(cand) + n, d = plane + centroid = point( + (cand[0][0] + cand[1][0] + cand[2][0]) / 3.0, + (cand[0][1] + cand[1][1] + cand[2][1]) / 3.0, + (cand[0][2] + cand[1][2] + cand[2][2]) / 3.0, + ) + vec = point(centroid[0] - center[0], + centroid[1] - center[1], + centroid[2] - center[2]) + dot_sign = n[0] * vec[0] + n[1] * vec[1] + n[2] * vec[2] + if dot_sign > epsilon_dot: + sense = -1 + elif dot_sign < -epsilon_dot: + sense = 1 + else: + center_eval = _plane_eval(plane, center) + sense = -1 if center_eval <= 0 else 1 + plane_with_sense = (n, d, sense) + key = _plane_key(plane_with_sense, tol) + if key not in seen: + seen.add(key) + planes.append(plane_with_sense) + snap_candidates.extend(_triangle_plane_intersection_points(cand, tri_plane, tol)) + unique_snap = _unique_points(snap_candidates, max(tol * 10.0, 1e-6)) if snap_candidates else [] + return planes, unique_snap + + +_DEFAULT_RAY_TOL = 1e-7 + + +def invalidate_surface_octree(s): + meta = _ensure_surface_metadata_dict(s) + meta.pop('_octree', None) + meta['_octree_dirty'] = True + + +def _bbox_overlap(box_a, box_b, tol): + return not ( + box_a[1][0] < box_b[0][0] - tol or + box_a[0][0] > box_b[1][0] + tol or + box_a[1][1] < box_b[0][1] - tol or + box_a[0][1] > box_b[1][1] + tol or + box_a[1][2] < box_b[0][2] - tol or + box_a[0][2] > box_b[1][2] + tol + ) + + +def _segment_bbox(p0, p1, pad): + return [ + point(min(p0[0], p1[0]) - pad, min(p0[1], p1[1]) - pad, min(p0[2], p1[2]) - pad), + point(max(p0[0], p1[0]) + pad, max(p0[1], p1[1]) + pad, max(p0[2], p1[2]) + pad), + ] + + +def _plane_eval(plane, p): + n, d = plane + return dot(n, p) - d + + +def _segment_plane_intersection(p1, p2, plane, tol): + n, d = plane + direction = sub(p2, p1) + denom = dot(n, direction) + if abs(denom) < tol: + return point((p1[0] + p2[0]) * 0.5, (p1[1] + p2[1]) * 0.5, (p1[2] + p2[2]) * 0.5) + t = (d - dot(n, p1)) / denom + t = max(0.0, min(1.0, t)) + return point(p1[0] + direction[0] * t, + p1[1] + direction[1] * t, + p1[2] + direction[2] * t) + + +def _dedupe_polygon(poly, tol): + if not poly: + return [] + deduped = [poly[0]] + for pt in poly[1:]: + if dist(pt, deduped[-1]) > tol: + deduped.append(pt) + if len(deduped) > 2 and dist(deduped[0], deduped[-1]) <= tol: + deduped[-1] = deduped[0] + deduped.pop() + return deduped + + +def _clip_polygon_against_plane(poly, plane, tol, keep_inside=True): + if not poly: + return [] + n, d, sense = plane + clipped = [] + prev = poly[-1] + prev_eval = _plane_eval((n, d), prev) + if sense <= 0: + prev_inside = prev_eval <= tol if keep_inside else prev_eval >= -tol + else: + prev_inside = prev_eval >= -tol if keep_inside else prev_eval <= tol + for curr in poly: + curr_eval = _plane_eval((n, d), curr) + if sense <= 0: + curr_inside = curr_eval <= tol if keep_inside else curr_eval >= -tol + else: + curr_inside = curr_eval >= -tol if keep_inside else curr_eval <= tol + if curr_inside: + if not prev_inside: + clipped.append(_segment_plane_intersection(prev, curr, (n, d), tol)) + clipped.append(point(curr)) + elif prev_inside: + clipped.append(_segment_plane_intersection(prev, curr, (n, d), tol)) + prev = curr + prev_eval = curr_eval + prev_inside = curr_inside + clipped = _dedupe_polygon(clipped, tol) + if not keep_inside and len(clipped) >= 3: + clipped = list(reversed(clipped)) + return clipped + + +def _split_polygon_by_plane(poly, plane, tol): + if not poly: + return [], [] + + n, d, sense = plane + evals = [_plane_eval((n, d), p) for p in poly] + max_eval = max(evals) + min_eval = min(evals) + + if sense <= 0: + if max_eval <= tol: + return [point(p) for p in poly], [] + if min_eval >= -tol: + return [], [[point(p) for p in poly]] + else: + if min_eval >= -tol: + return [point(p) for p in poly], [] + if max_eval <= tol: + return [], [[point(p) for p in poly]] + + inside = _clip_polygon_against_plane(poly, plane, tol, keep_inside=True) + outside = _clip_polygon_against_plane(poly, plane, tol, keep_inside=False) + outside_polys = [outside] if outside else [] + return inside, outside_polys + + +def _split_polygon_by_planes(poly, planes, tol): + inside_polys = [poly] + outside_polys = [] + for plane in planes: + next_inside = [] + for current in inside_polys: + inside, outside = _split_polygon_by_plane(current, plane, tol) + outside_polys.extend(outside) + if inside: + next_inside.append(inside) + inside_polys = next_inside + if not inside_polys: + break + return inside_polys, outside_polys + + +def _triangulate_polygon(poly, reference_normal=None): + if len(poly) < 3: + return [] + if len(poly) == 3: + tri = [point(poly[0]), point(poly[1]), point(poly[2])] + if reference_normal is not None: + v01 = sub(tri[1], tri[0]) + v02 = sub(tri[2], tri[0]) + if dot(cross(v01, v02), reference_normal) < 0: + tri[1], tri[2] = tri[2], tri[1] + return [tri] + + anchor = point(poly[0]) + triangles = [] + for i in range(1, len(poly) - 1): + tri = [anchor, point(poly[i]), point(poly[i + 1])] + if reference_normal is not None: + v01 = sub(tri[1], tri[0]) + v02 = sub(tri[2], tri[0]) + if dot(cross(v01, v02), reference_normal) < 0: + tri[1], tri[2] = tri[2], tri[1] + triangles.append(tri) + return triangles + + +def _triangle_plane(tri): + p0, n = _geom3d().tri2p0n(tri) + d = dot(n, p0) + return n, d + + +def _plane_key(plane, tol): + n, d, sense = plane + scale = 1.0 / tol + return (round(n[0] * scale), round(n[1] * scale), round(n[2] * scale), round(d * scale), sense) + + +def _sub3(a, b): + return [a[0] - b[0], a[1] - b[1], a[2] - b[2]] + + +def _dot3(a, b): + return a[0] * b[0] + a[1] * b[1] + a[2] * b[2] + + +def _cross3(a, b): + return [ + a[1] * b[2] - a[2] * b[1], + a[2] * b[0] - a[0] * b[2], + a[0] * b[1] - a[1] * b[0], + ] + + +def _mag3(v): + return math.sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]) + + +def _normalize3(v, tol): + m = _mag3(v) + if m < tol: + return [0.0, 0.0, 0.0] + return [v[0] / m, v[1] / m, v[2] / m] + + +def _lerp_point(p0, p1, t): + return point( + p0[0] + (p1[0] - p0[0]) * t, + p0[1] + (p1[1] - p0[1]) * t, + p0[2] + (p1[2] - p0[2]) * t, + ) + + + +def _point_in_triangle(pt, tri, tol): + a = [tri[0][0], tri[0][1], tri[0][2]] + b = [tri[1][0], tri[1][1], tri[1][2]] + c = [tri[2][0], tri[2][1], tri[2][2]] + p = [pt[0], pt[1], pt[2]] + + v0 = _sub3(c, a) + v1 = _sub3(b, a) + v2 = _sub3(p, a) + + dot00 = _dot3(v0, v0) + dot01 = _dot3(v0, v1) + dot02 = _dot3(v0, v2) + dot11 = _dot3(v1, v1) + dot12 = _dot3(v1, v2) + + denom = dot00 * dot11 - dot01 * dot01 + if abs(denom) < tol: + return False + inv = 1.0 / denom + u = (dot11 * dot02 - dot01 * dot12) * inv + v = (dot00 * dot12 - dot01 * dot02) * inv + return u >= -tol and v >= -tol and (u + v) <= 1.0 + tol + + +def _segment_triangle_intersection_param(p0, p1, tri, tol): + n, d = _triangle_plane(tri) + dir_vec = _sub3([p1[0], p1[1], p1[2]], [p0[0], p0[1], p0[2]]) + denom = _dot3(n, dir_vec) + if abs(denom) < tol: + return None + numer = d - (n[0] * p0[0] + n[1] * p0[1] + n[2] * p0[2]) + t = numer / denom + if t < -tol or t > 1.0 + tol: + return None + t_clamped = max(0.0, min(1.0, t)) + intersection = _lerp_point(p0, p1, t_clamped) + if not _point_in_triangle(intersection, tri, tol): + return None + return (t_clamped, intersection) + + +def _collect_segment_intersections(p0, p1, solid, tol): + results = [] + query_box = _segment_bbox(p0, p1, tol * 5) + for surf in solid[1]: + tree = surface_octree(surf) + if tree is None: + candidates = list(_iter_triangles_from_surface(surf)) + else: + elems = tree.getElements(query_box) + candidates = [] + for elem in elems: + if isinstance(elem, list): + candidates.append(elem) + elif isinstance(elem, tuple): + candidates.append(elem[0]) + else: + candidates.append(elem) + for tri in candidates: + res = _segment_triangle_intersection_param(p0, p1, tri, tol) + if res is not None: + results.append(res) + return results + + +def _unique_points(points, tol): + unique = [] + for pt in points: + candidate = point(pt) + if not any(dist(candidate, existing) <= tol for existing in unique): + unique.append(candidate) + return unique + + +def _points_to_polygon(points, reference_normal, tol): + unique = _unique_points(points, tol) + if len(unique) < 3: + return [] + + centroid = point( + sum(p[0] for p in unique) / len(unique), + sum(p[1] for p in unique) / len(unique), + sum(p[2] for p in unique) / len(unique), + ) + + basis_u = None + for pt in unique: + vec = _sub3([pt[0], pt[1], pt[2]], [centroid[0], centroid[1], centroid[2]]) + if _mag3(vec) > tol: + basis_u = _normalize3(vec, tol) + break + if basis_u is None: + return [] + + ref = _normalize3([reference_normal[0], reference_normal[1], reference_normal[2]], tol) + basis_v = _cross3(ref, basis_u) + if _mag3(basis_v) < tol: + basis_v = _cross3(basis_u, ref) + basis_v = _normalize3(basis_v, tol) + + def _angle(pt): + vec = _sub3([pt[0], pt[1], pt[2]], [centroid[0], centroid[1], centroid[2]]) + x = _dot3(vec, basis_u) + y = _dot3(vec, basis_v) + return math.atan2(y, x) + + ordered = sorted(unique, key=_angle) + compact = [ordered[0]] + for pt in ordered[1:]: + if dist(compact[-1], pt) > tol: + compact.append(pt) + if dist(compact[0], compact[-1]) <= tol: + compact[-1] = compact[0] + compact = compact[:-1] + if len(compact) < 3: + return [] + + v0 = _sub3([compact[1][0], compact[1][1], compact[1][2]], [compact[0][0], compact[0][1], compact[0][2]]) + v1 = _sub3([compact[2][0], compact[2][1], compact[2][2]], [compact[1][0], compact[1][1], compact[1][2]]) + normal = _cross3(v0, v1) + if _mag3(normal) < tol: + return [] + if _dot3(normal, [reference_normal[0], reference_normal[1], reference_normal[2]]) < 0: + compact.reverse() + + return compact + + +def _orient_triangle(tri, reference_normal, tol): + pts = [point(tri[0]), point(tri[1]), point(tri[2])] + v01 = sub(pts[1], pts[0]) + v02 = sub(pts[2], pts[0]) + normal = cross(v01, v02) + area_mag = mag(normal) + if area_mag < tol: + return None + ref = reference_normal + if ref is not None: + ref_vec = [ref[0], ref[1], ref[2]] + if dot(normal, ref_vec) < 0: + pts[1], pts[2] = pts[2], pts[1] + return pts + + +def _clip_triangle_against_solid(tri, target, tol, reference_normal): + polygon = [point(tri[0]), point(tri[1]), point(tri[2])] + centroid = point( + (tri[0][0] + tri[1][0] + tri[2][0]) / 3.0, + (tri[0][1] + tri[1][1] + tri[2][1]) / 3.0, + (tri[0][2] + tri[1][2] + tri[2][2]) / 3.0, + ) + + tri_plane = _triangle_plane(tri) + planes, snap_candidates = _candidate_planes_for_triangle(tri, target, tri_plane, tol) + + inside_polys = [] + if planes: + inside_raw, _ = _split_polygon_by_planes(polygon, planes, tol) + for poly in inside_raw: + cleaned = _dedupe_polygon([point(p) for p in poly], tol) + if len(cleaned) >= 3: + inside_polys.append(cleaned) + + if not inside_polys: + if solid_contains_point(target, centroid, tol=tol): + inside_polys = [polygon] + else: + oriented = _orient_triangle(tri, reference_normal, tol) + if oriented: + return [], [oriented], False + return [], [], False + + _, _, to_local, to_world = _geom3d().tri2p0n(tri, basis=True) + + if snap_candidates: + snap_tol = max(5e-2, tol * 1e6) + + def _snap_point(pt): + best = snap_tol + best_candidate = None + for cand in snap_candidates: + d = dist(pt, cand) + if d < best: + best = d + best_candidate = cand + return point(best_candidate) if best_candidate is not None else point(pt) + + snapped = [] + for poly in inside_polys: + snapped_poly = [_snap_point(pt) for pt in poly] + snapped_poly = _dedupe_polygon(snapped_poly, tol) + if len(snapped_poly) >= 3: + snapped.append(snapped_poly) + if snapped: + inside_polys = snapped + + def _to_local(pt): + loc = to_local.mul(pt) + return (loc[0], loc[1]) + + def _from_local(xy): + local_pt = point(xy[0], xy[1], 0.0) + world_pt = to_world.mul(local_pt) + return point(world_pt[0], world_pt[1], world_pt[2]) + + outer_loop = [_to_local(pt) for pt in polygon] + inside_loops = [[_to_local(pt) for pt in poly] for poly in inside_polys] + + inside_triangles = [] + for loop in inside_loops: + tri2d_list = triangulate_polygon(loop) + for tri2d in tri2d_list: + tri3d = [_from_local(pt) for pt in tri2d] + oriented = _orient_triangle(tri3d, reference_normal, tol) + if oriented: + inside_triangles.append(oriented) + + holes = inside_loops if inside_loops else None + outside_triangles = [] + tri2d_list = triangulate_polygon(outer_loop, holes=holes) + for tri2d in tri2d_list: + tri3d = [_from_local(pt) for pt in tri2d] + # Don't call _orient_triangle! The winding order is already preserved through: + # 1. Original triangle had correct outward normal + # 2. Clipping preserves vertex order (just removes parts) + # 3. Transform to 2D preserves order + # 4. triangulate_polygon produces CCW triangles in 2D + # 5. Transform back to 3D should maintain correct orientation + # Check for degeneracy with better quality metric: + v01 = sub(tri3d[1], tri3d[0]) + v02 = sub(tri3d[2], tri3d[0]) + cross_prod = cross(v01, v02) + area = mag(cross_prod) + if area >= tol: + # Also filter out very thin slivers by checking aspect ratio + # A degenerate sliver has tiny area relative to edge lengths + edge_lengths = [mag(v01), mag(v02), dist(tri3d[1], tri3d[2])] + max_edge = max(edge_lengths) + # For a reasonable triangle, area should be at least 1% of max_edge^2 + # This filters out slivers where one edge is extremely small + quality_threshold = 0.01 * max_edge * max_edge + if area >= max(tol, quality_threshold): + outside_triangles.append([point(tri3d[0]), point(tri3d[1]), point(tri3d[2])]) + + split = bool(inside_triangles) and bool(outside_triangles) + return inside_triangles, outside_triangles, split + + + + +def stitch_open_edges(triangles, tol): + """Experimental: attempt to close open edge loops by triangulating them. + + Parameters + ---------- + triangles : Iterable + A sequence of triangle coordinate lists ``[[x, y, z, w], ...]`` + describing a mesh with potential boundary edges. + tol : float + Tolerance used when deduplicating vertices and detecting shared + edges. Typically reuse ``_DEFAULT_RAY_TOL``. + + Returns + ------- + list + A new list of triangles with any stitched faces appended. + """ + if not triangles: + return triangles + + oriented_edges = [] + canonical_counts = {} + base_triangles = [] + for tri in triangles: + try: + n = _triangle_normal(tri) + except ValueError: + continue + pts = [point(tri[0]), point(tri[1]), point(tri[2])] + base_triangles.append(pts) + edges = ((pts[0], pts[1]), (pts[1], pts[2]), (pts[2], pts[0])) + for start, end in edges: + canonical = _canonical_edge_key(start, end) + canonical_counts[canonical] = canonical_counts.get(canonical, 0) + 1 + oriented_edges.append((start, end, n)) + + boundary = [] + for start, end, normal in oriented_edges: + if canonical_counts[_canonical_edge_key(start, end)] == 1: + boundary.append((start, end, normal)) + + if not boundary: + return base_triangles + + def edge_id(start, end): + return (_point_to_key(start), _point_to_key(end)) + + adjacency = {} + for start, end, normal in boundary: + adjacency.setdefault(_point_to_key(start), []).append((start, end, normal)) + + used = set() + loops = [] + for start, end, normal in boundary: + eid = edge_id(start, end) + if eid in used: + continue + loop = [] + normals = [] + current_start = start + current_end = end + current_normal = normal + start_key = _point_to_key(current_start) + while True: + loop.append(point(current_start)) + normals.append(current_normal) + used.add(edge_id(current_start, current_end)) + next_key = _point_to_key(current_end) + if next_key == start_key: + break + candidates = adjacency.get(next_key) + if not candidates: + loop = [] + break + next_edge = None + for candidate in candidates: + cid = edge_id(candidate[0], candidate[1]) + if cid not in used: + next_edge = candidate + break + if next_edge is None: + loop = [] + break + current_start, current_end, current_normal = next_edge + if loop and len(loop) >= 3: + loops.append((loop, normals)) + + if not loops: + return base_triangles + + stitched = list(base_triangles) + for loop_points, normals in loops: + polygon = _unique_points(loop_points, tol) + if len(polygon) < 3: + continue + avg = [0.0, 0.0, 0.0] + for n in normals: + avg[0] += n[0] + avg[1] += n[1] + avg[2] += n[2] + avg = _normalize3(avg, tol) + if _mag3(avg) < tol: + avg = normals[0] + poly = _points_to_polygon(polygon, avg, tol) + if not poly: + continue + for tri in _triangulate_polygon(poly, avg): + try: + normal = _triangle_normal(tri) + except ValueError: + continue + if _dot3(normal, [avg[0], avg[1], avg[2]]) < 0: + tri = [tri[0], tri[2], tri[1]] + try: + normal = _triangle_normal(tri) + except ValueError: + continue + stitched.append(tri) + return stitched + +def stitch_solid(sld, tol=_DEFAULT_RAY_TOL): + """Experimental helper that runs :func:`stitch_open_edges` on a solid.""" + + if not _geom3d().issolid(sld, fast=False): + raise ValueError('invalid solid passed to stitch_solid') + + triangles = list(_iter_triangles_from_solid(sld)) + stitched = stitch_open_edges(triangles, tol) + if not stitched: + return sld + + surface_result = _surface_from_triangles(stitched) + if surface_result: + return _geom3d().solid([surface_result], [], ['utility', 'stitch_open_edges']) + return sld + + + +def _boolean_fragments(source, target, tol): + outside_tris = [] + inside_tris = [] + inside_overlap = [] + + for tri in _iter_triangles_from_solid(source): + reference_normal = _triangle_normal(tri) + inside_parts, outside_parts, split = _clip_triangle_against_solid( + tri, target, tol, reference_normal + ) + if inside_parts: + inside_tris.extend(inside_parts) + if outside_parts: + outside_tris.extend(outside_parts) + if split and inside_parts: + inside_overlap.extend(inside_parts) + return outside_tris, inside_tris, inside_overlap + + + +def _ray_triangle_intersection(origin, direction, triangle, tol=_DEFAULT_RAY_TOL): + v0, v1, v2 = triangle + e1 = sub(v1, v0) + e2 = sub(v2, v0) + h = cross(direction, e2) + a = dot(e1, h) + if abs(a) < tol: + return None + f = 1.0 / a + s = sub(origin, v0) + u = f * dot(s, h) + if u < -tol or u > 1.0 + tol: + return None + q = cross(s, e1) + v = f * dot(direction, q) + if v < -tol or u + v > 1.0 + tol: + return None + t = f * dot(e2, q) + if t < -tol: + return None + w = 1.0 - u - v + if w < -tol: + return None + hit = point(origin[0] + direction[0] * t, + origin[1] + direction[1] * t, + origin[2] + direction[2] * t) + return t, hit, (u, v, w) + + +def surface_octree(s, rebuild=False): + meta = _ensure_surface_metadata_dict(s) + tree = meta.get('_octree') + if rebuild or meta.get('_octree_dirty') or tree is None: + verts = s[1] + faces = s[3] + if not faces: + meta['_octree'] = None + meta['_octree_dirty'] = False + return None + extent = 0.0 + for axis in range(3): + vals = [v[axis] for v in verts] + extent = max(extent, max(vals) - min(vals)) + mindim = max(extent / 32.0, epsilon) + tree = NTree(n=8, mindim=mindim) + for face in faces: + if len(face) != 3: + continue + tri = [verts[face[0]], verts[face[1]], verts[face[2]]] + tree.addElement(tri) + tree.updateTree() + meta['_octree'] = tree + meta['_octree_dirty'] = False + return meta['_octree'] + + +def _surface_from_triangles(triangles): + if not triangles: + return None + verts = [] + normals = [] + faces = [] + for tri in triangles: + v01 = sub(tri[1], tri[0]) + v02 = sub(tri[2], tri[0]) + if mag(cross(v01, v02)) < epsilon: + continue + n = _triangle_normal(tri) + indices = [] + for pt in tri: + verts.append(point(pt)) + normals.append([n[0], n[1], n[2], 0.0]) + indices.append(len(verts) - 1) + faces.append(indices) + surface_obj = ['surface', verts, normals, faces, [], []] + _ensure_surface_metadata_dict(surface_obj) + invalidate_surface_octree(surface_obj) + return surface_obj + + +def _iter_triangles_from_surface(surf): + verts = surf[1] + faces = surf[3] + for face in faces: + if len(face) != 3: + continue + yield [verts[face[0]], verts[face[1]], verts[face[2]]] + + +def _iter_triangles_from_solid(sld): + for surf in sld[1]: + yield from _iter_triangles_from_surface(surf) + + +def _group_hits(hits, tol): + if not hits: + return [] + hits.sort(key=lambda x: x[0]) + groups = [[hits[0]]] + for hit in hits[1:]: + if abs(hit[0] - groups[-1][-1][0]) <= tol: + groups[-1].append(hit) + else: + groups.append([hit]) + return groups + + +def _triangle_normal(tri): + _, n = _geom3d().tri2p0n(tri) + return n +def solid_contains_point(sld, p, tol=_DEFAULT_RAY_TOL): + if not _geom3d().issolid(sld, fast=False): + raise ValueError('invalid solid passed to solid_contains_point') + if not _geom3d().ispoint(p): + raise ValueError('invalid point passed to solid_contains_point') + + surfaces = sld[1] + if not surfaces: + return False + + box = _geom3d().solidbbox(sld) + if box: + expanded = [ + point(box[0][0] - tol, box[0][1] - tol, box[0][2] - tol), + point(box[1][0] + tol, box[1][1] + tol, box[1][2] + tol), + ] + if not _geom3d().isinsidebbox(expanded, p): + return False + extent = max(box[1][0] - box[0][0], + box[1][1] - box[0][1], + box[1][2] - box[0][2]) + ray_length = max(extent * 1.5, 1.0) + else: + ray_length = 1.0 + + directions = [ + vect(1, 0, 0, 0), + vect(-1, 0, 0, 0), + vect(0, 1, 0, 0), + vect(0, -1, 0, 0), + vect(0, 0, 1, 0), + vect(0, 0, -1, 0), + ] + + seen_inside = False + for direction in directions: + far_point = point(p[0] + direction[0] * ray_length, + p[1] + direction[1] * ray_length, + p[2] + direction[2] * ray_length) + query_box = _segment_bbox(p, far_point, tol * 5) + hits = [] + for surf in surfaces: + tree = surface_octree(surf) + if tree is None: + candidates = list(_iter_triangles_from_surface(surf)) + else: + elems = tree.getElements(query_box) + candidates = [] + for elem in elems: + if isinstance(elem, list): + candidates.append(elem) + elif isinstance(elem, tuple): + candidates.append(elem[0]) + else: + candidates.append(elem) + if not candidates: + candidates = list(_iter_triangles_from_surface(surf)) + for tri in candidates: + result = _ray_triangle_intersection(p, direction, tri, tol=tol) + if result is None: + continue + t_hit, hit_point, _ = result + if t_hit < -tol or t_hit > ray_length + tol: + continue + if t_hit <= tol: + return True + normal = _triangle_normal(tri) + sign = -1 if dot(normal, direction) > 0 else 1 + hits.append((t_hit, hit_point, sign)) + groups = _group_hits(hits, tol) + parity = 0 + for group in groups: + sign_sum = sum(hit[2] for hit in group) + if sign_sum == 0: + continue + parity ^= 1 + if parity == 0: + return False + seen_inside = True + return seen_inside + + +def solids_intersect(a, b, tol=_DEFAULT_RAY_TOL): + if not _geom3d().issolid(a, fast=False) or not _geom3d().issolid(b, fast=False): + raise ValueError('invalid solid passed to solids_intersect') + if not a[1] or not b[1]: + return False + + box_a = _geom3d().solidbbox(a) + box_b = _geom3d().solidbbox(b) + if box_a and box_b and not _bbox_overlap(box_a, box_b, tol): + return False + + center_a = point((box_a[0][0] + box_a[1][0]) / 2.0, + (box_a[0][1] + box_a[1][1]) / 2.0, + (box_a[0][2] + box_a[1][2]) / 2.0) + center_b = point((box_b[0][0] + box_b[1][0]) / 2.0, + (box_b[0][1] + box_b[1][1]) / 2.0, + (box_b[0][2] + box_b[1][2]) / 2.0) + if solid_contains_point(a, center_b, tol=tol) or solid_contains_point(b, center_a, tol=tol): + return True + + for tri_a in _iter_triangles_from_solid(a): + for tri_b in _iter_triangles_from_solid(b): + if _geom3d().triTriIntersect(tri_a, tri_b): + return True + return False + + +def solid_boolean(a, b, operation, tol=_DEFAULT_RAY_TOL, *, stitch=False): + if operation not in {'union', 'intersection', 'difference'}: + raise ValueError(f'unsupported solid boolean operation {operation!r}') + + outside_a, inside_a, overlap_a = _boolean_fragments(a, b, tol) + outside_b, inside_b, overlap_b = _boolean_fragments(b, a, tol) + + if operation == 'union': + result_tris = outside_a + outside_b + if not result_tris: + result_tris = inside_a if inside_a else inside_b + + # Filter out triangles in the interior of the overlap region + # For union, if a triangle's center is inside both input solids, + # it's in the interior and should not be on the surface + filtered_tris = [] + for tri in result_tris: + center = point((tri[0][0] + tri[1][0] + tri[2][0]) / 3.0, + (tri[0][1] + tri[1][1] + tri[2][1]) / 3.0, + (tri[0][2] + tri[1][2] + tri[2][2]) / 3.0) + # Use a tighter tolerance for containment check to avoid false positives + check_tol = tol * 100 + in_a = solid_contains_point(a, center, tol=check_tol) + in_b = solid_contains_point(b, center, tol=check_tol) + # Keep triangle only if it's not clearly inside both solids + if not (in_a and in_b): + filtered_tris.append(tri) + result_tris = filtered_tris + + elif operation == 'intersection': + result_tris = inside_a + inside_b + if not result_tris: + result_tris = overlap_a + overlap_b + else: # difference + reversed_inside_b = [[tri[0], tri[2], tri[1]] for tri in inside_b] + if not reversed_inside_b and overlap_b: + reversed_inside_b = [[tri[0], tri[2], tri[1]] for tri in overlap_b] + result_tris = outside_a + reversed_inside_b + + if stitch: + result_tris = stitch_open_edges(result_tris, tol) + + surface_result = _surface_from_triangles(result_tris) + if surface_result: + return _geom3d().solid([surface_result], [], ['boolean', operation]) + return _geom3d().solid([], [], ['boolean', operation]) + + + diff --git a/src/yapcad/boolean/trimesh_engine.py b/src/yapcad/boolean/trimesh_engine.py new file mode 100644 index 0000000..f32c485 --- /dev/null +++ b/src/yapcad/boolean/trimesh_engine.py @@ -0,0 +1,155 @@ +"""Trimesh-backed boolean engine for yapCAD solids. + +This engine is optional. It converts yapCAD solids to ``trimesh.Trimesh`` +instances, dispatches boolean operations via :mod:`trimesh.boolean`, and +converts the resulting mesh back into a yapCAD solid. + +Availability depends on both the ``trimesh`` package and at least one +boolean backend supported by ``trimesh`` (e.g. Blender, OpenSCAD, Cork). +""" + +from __future__ import annotations + +from typing import Iterable + +import os +import numpy as np + +try: + import trimesh +except ImportError: # pragma: no cover - optional dependency + trimesh = None # type: ignore[assignment] + +from . import native as _native + +ENGINE_NAME = "trimesh" + + +def engines_available() -> set[str]: + """Return the set of trimesh boolean backends that are operational.""" + + if trimesh is None: # pragma: no cover - optional dependency + return set() + return set(trimesh.boolean.engines_available) + + +def is_available(backend: str | None = None) -> bool: + """Check whether the engine can run (trimesh + backend present).""" + + available = engines_available() + if not available: + return False + if backend is None: + return True + return backend in available + + +def _solid_to_mesh(sld) -> "trimesh.Trimesh": + if trimesh is None: # pragma: no cover - optional dependency + raise RuntimeError("trimesh is not installed") + + triangles = list(_native._iter_triangles_from_solid(sld)) + if not triangles: + return trimesh.Trimesh(vertices=np.zeros((0, 3)), faces=np.zeros((0, 3), dtype=np.int64), process=False) + + vertex_map: dict[tuple[float, float, float], int] = {} + verts = [] + faces = [] + for tri in triangles: + face_inds = [] + for pt in tri: + key = (round(pt[0], 9), round(pt[1], 9), round(pt[2], 9)) + idx = vertex_map.get(key) + if idx is None: + idx = len(verts) + vertex_map[key] = idx + verts.append([pt[0], pt[1], pt[2]]) + face_inds.append(idx) + faces.append(face_inds) + + verts_np = np.asarray(verts, dtype=float) + faces_np = np.asarray(faces, dtype=np.int64) + mesh = trimesh.Trimesh(vertices=verts_np, faces=faces_np, process=False) + mesh.remove_unreferenced_vertices() + if not mesh.is_watertight: + mesh.merge_vertices() + mesh.remove_duplicate_faces() + mesh.remove_degenerate_faces() + mesh.process(validate=True) + try: + mesh.fix_normals() + except Exception: + pass + return mesh + + +def _mesh_to_solid(mesh: "trimesh.Trimesh", operation: str) -> list: + triangles = np.asarray(mesh.triangles) + if triangles.size == 0: + return _native._geom3d().solid([], [], ['boolean', f'{ENGINE_NAME}:{operation}']) + + tri_points = [] + for tri in triangles: + tri_points.append([ + _native.point(tri[0][0], tri[0][1], tri[0][2]), + _native.point(tri[1][0], tri[1][1], tri[1][2]), + _native.point(tri[2][0], tri[2][1], tri[2][2]), + ]) + + surface = _native._surface_from_triangles(tri_points) + if surface is None: + return _native._geom3d().solid([], [], ['boolean', f'{ENGINE_NAME}:{operation}']) + return _native._geom3d().solid([surface], [], ['boolean', f'{ENGINE_NAME}:{operation}']) + + +def solid_boolean(a, b, operation: str, tol=_native._DEFAULT_RAY_TOL, *, stitch: bool = False, + backend: str | None = None): + """Perform a boolean between ``a`` and ``b`` using trimesh.""" + + if trimesh is None: # pragma: no cover - optional dependency + raise RuntimeError("trimesh is not installed; install trimesh to enable this engine") + + available = engines_available() + if backend is not None and backend not in available: + raise RuntimeError( + f"trimesh backend '{backend}' is not available; install the appropriate binary (available: {available})" + ) + if backend is None and not available: + raise RuntimeError( + "no trimesh boolean backends are available; install Blender, OpenSCAD, Cork, or another supported engine" + ) + + mesh_a = _solid_to_mesh(a) + mesh_b = _solid_to_mesh(b) + + backup_cache = None + if backend == 'blender': + backup_cache = os.environ.get('ARCH_CACHE_LINE_SIZE') + os.environ['ARCH_CACHE_LINE_SIZE'] = '64' + + op = operation.lower() + try: + if op == 'union': + result = trimesh.boolean.union([mesh_a, mesh_b], engine=backend, check_volume=False) + elif op == 'intersection': + result = trimesh.boolean.intersection([mesh_a, mesh_b], engine=backend, check_volume=False) + elif op == 'difference': + result = trimesh.boolean.difference([mesh_a, mesh_b], engine=backend, check_volume=False) + else: + raise RuntimeError(f"unsupported boolean operation '{operation}' for trimesh engine") + except Exception as exc: # pragma: no cover - depends on external binaries + raise RuntimeError(f"trimesh boolean operation failed: {exc}") from exc + finally: + if backend == 'blender': + if backup_cache is None: + os.environ.pop('ARCH_CACHE_LINE_SIZE', None) + else: + os.environ['ARCH_CACHE_LINE_SIZE'] = backup_cache + + if result is None or result.faces.size == 0: + return _native._geom3d().solid([], [], ['boolean', f'{ENGINE_NAME}:{operation}']) + + return _mesh_to_solid(result, operation) + + +__all__ = ['ENGINE_NAME', 'is_available', 'solid_boolean', 'engines_available'] diff --git a/src/yapcad/combine.py b/src/yapcad/combine.py index 733a0be..3054d28 100644 --- a/src/yapcad/combine.py +++ b/src/yapcad/combine.py @@ -1,4 +1,5 @@ -## yapCAD boolen operation support +## yapCAD boolen operation support for 2D closed curves. For three dimensional +## boolean operations, see geom3d.py from yapcad.geom import * from yapcad.geom_util import * @@ -47,6 +48,10 @@ def _prepare_geom(self, geom_obj): else: gl = geom_obj + # If gl is a single arc/line/poly, wrap it in a list for geomlist2poly_with_holes + if isarc(gl) or isline(gl) or ispoly(gl): + gl = [gl] + try: outer, holes = geomlist2poly_with_holes(gl, self.__minang, self.__minlen, checkcont=False) except Exception: diff --git a/src/yapcad/geom3d.py b/src/yapcad/geom3d.py index 3013ed3..130520b 100644 --- a/src/yapcad/geom3d.py +++ b/src/yapcad/geom3d.py @@ -6,7 +6,21 @@ from yapcad.geom_util import * from yapcad.xform import * from functools import reduce +import os +from yapcad.octtree import NTree + from yapcad.triangulator import triangulate_polygon +import yapcad.boolean.native as _boolean_native +_DEFAULT_RAY_TOL = _boolean_native._DEFAULT_RAY_TOL +invalidate_surface_octree = _boolean_native.invalidate_surface_octree +surface_octree = _boolean_native.surface_octree +_surface_from_triangles = _boolean_native._surface_from_triangles +_iter_triangles_from_surface = _boolean_native._iter_triangles_from_surface +_iter_triangles_from_solid = _boolean_native._iter_triangles_from_solid +stitch_open_edges = _boolean_native.stitch_open_edges +stitch_solid = _boolean_native.stitch_solid +solid_contains_point = _boolean_native.solid_contains_point +solids_intersect = _boolean_native.solids_intersect """ @@ -499,7 +513,25 @@ def surfacebbox(s): if not issurface(s): raise ValueError('bad surface passed to surfacebbox') return polybbox(s[1]) - + + + + +def solid_boolean(a, b, operation, tol=_DEFAULT_RAY_TOL, *, stitch=False, engine=None): + selected_raw = engine or os.environ.get('YAPCAD_BOOLEAN_ENGINE', 'native') + backend = None + if selected_raw and ':' in selected_raw: + selected, backend = selected_raw.split(':', 1) + else: + selected = selected_raw + if selected == 'native': + return _boolean_native.solid_boolean(a, b, operation, tol=tol, stitch=stitch) + if selected == 'trimesh': + from yapcad.boolean import trimesh_engine + backend = backend or os.environ.get('YAPCAD_TRIMESH_BACKEND') + return trimesh_engine.solid_boolean(a, b, operation, tol=tol, stitch=stitch, backend=backend) + raise ValueError(f'unknown boolean engine {selected_raw!r}') + def issurface(s,fast=True): """ Check to see if ``s`` is a valid surface. @@ -511,8 +543,8 @@ def filterInds(inds,verts): return (len(list(filter(lambda x: not (isinstance(x,int) or x < 0 or x >= l), inds))) == 0) - - if not isinstance(s,list) or s[0] != 'surface' or len(s) not in (6,7): + + if not isinstance(s,list) or len(s) not in (6,7) or s[0] != 'surface': return False if fast: return True @@ -657,7 +689,7 @@ def issolid(s,fast=True): of surfaces completely bounds a volume of space without holes """ - if not isinstance(s,list) or s[0] != 'solid' or len(s) not in (4,5): + if not isinstance(s,list) or len(s) not in (4,5) or s[0] != 'solid': return False if fast: return True diff --git a/src/yapcad/geom3d_util.py b/src/yapcad/geom3d_util.py index 996c010..9a344a2 100644 --- a/src/yapcad/geom3d_util.py +++ b/src/yapcad/geom3d_util.py @@ -81,6 +81,7 @@ def makeIcoPoints(center,radius): [0,2,4],[0,4,6],[0,6,8],[0,8,10],[0,10,2] ] vertexHash = {} +_vertexHash_owner = None def addVertex(nv,nn,verts,normals): """ Utility function that takes a vertex and associated normal and a @@ -91,18 +92,29 @@ def addVertex(nv,nn,verts,normals): returns the index, and the (potentiall updated) lists """ - global vertexHash - if len(verts) == 0: + global vertexHash, _vertexHash_owner + owner_id = id(verts) + if _vertexHash_owner != owner_id or len(verts) == 0: vertexHash = {} + _vertexHash_owner = owner_id found = False - vkey = f"{nv[0]:.2f}{nv[1]:.2f}{nv[2]:.2f}" + # Normalize tiny values to avoid "-0.00" != "0.00" hash key mismatch + x = 0.0 if abs(nv[0]) < epsilon else nv[0] + y = 0.0 if abs(nv[1]) < epsilon else nv[1] + z = 0.0 if abs(nv[2]) < epsilon else nv[2] + vkey = f"{x:.2f}{y:.2f}{z:.2f}" if vkey in vertexHash: found = True inds = vertexHash[vkey] + valid_inds = [] for i in inds: + if i >= len(verts): + continue if vclose(nv,verts[i]): return i,verts,normals + valid_inds.append(i) + vertexHash[vkey] = valid_inds verts.append(nv) normals.append(nn) @@ -347,20 +359,32 @@ def conic(baser,topr,height, center=point(0,0,0),angr=10): ['procedure',call]) else: topP = add(center,point(0,0,height)) - conV = [ topP ] + baseV + # Only use perimeter vertices (baseV[1:]), skip the center point + conV = [ topP ] + baseV[1:] ll = len(conV) - conN = [[0,0,1,0]] + # Initialize all normals to a default value + conN = [[0,0,1,0] for _ in range(ll)] conF = [] - for i in range(1,ll): - p0= conV[0] - p1= conV[(i-1)%ll] - p2= conV[(i+1)%ll] + # ll = apex(1) + perimeter vertices(36) = 37 + # Perimeter vertex indices: 1 to ll-1 + num_perimeter = ll - 1 + for i in range(1, ll): + p0 = conV[0] # apex + # Wrap indices within perimeter range [1, ll-1] + prev_idx = ((i - 2) % num_perimeter) + 1 + next_idx = (i % num_perimeter) + 1 + p1 = conV[prev_idx] + p2 = conV[next_idx] - conF.append([0,i,(i+1)%ll]) - pp,n0 = tri2p0n([p0,p1,p2]) + try: + pp, n0 = tri2p0n([p0, p1, p2]) + except ValueError: + # Skip degenerate faces near the apex + continue - conN.append(n0) + conF.append([0, i, next_idx]) + conN[i] = n0 conS = surface(conV,conN,conF) @@ -384,28 +408,112 @@ def makeRevolutionSurface(contour,zStart,zEnd,steps,arcSamples=36): degStep = 360.0/arcSamples radStep = pi2/arcSamples + + # Pre-compute cos/sin values to avoid floating point errors at the seam + # Explicitly ensure that index 0 uses exact values + angle_cos = [] + angle_sin = [] + for i in range(arcSamples): + if i == 0: + angle_cos.append(1.0) + angle_sin.append(0.0) + else: + angle = i * radStep + angle_cos.append(math.cos(angle)) + angle_sin.append(math.sin(angle)) + + # Check if we need pole caps + r_start = contour(zStart) + r_end = contour(zEnd) + need_start_cap = r_start < epsilon * 10 + need_end_cap = r_end < epsilon * 10 + + # Add pole vertices if needed + start_pole_idx = None + end_pole_idx = None + if need_start_cap: + pole_point = [0.0, 0.0, zStart, 1.0] + pole_normal = [0.0, 0.0, -1.0, 0.0] # Points downward for bottom pole + start_pole_idx, sV, sN = addVertex(pole_point, pole_normal, sV, sN) + + if need_end_cap: + pole_point = [0.0, 0.0, zEnd, 1.0] + pole_normal = [0.0, 0.0, 1.0, 0.0] # Points upward for top pole + end_pole_idx, sV, sN = addVertex(pole_point, pole_normal, sV, sN) + for i in range(steps): z = i*zD+zStart r0 = contour(z) r1 = contour(z+zD) - if r0 < epsilon*10: - r0 = epsilon*10 - if r1 < epsilon*10: - r1 = epsilon*10 + + # Handle pole caps + if i == 0 and need_start_cap: + # Create triangular faces from pole to first ring + if r1 < epsilon: + r1 = epsilon + for j in range(arcSamples): + a1_idx = j + a2_idx = (j+1) % arcSamples + + pp1 = [angle_cos[a1_idx]*r1, angle_sin[a1_idx]*r1, z+zD, 1.0] + pp2 = [angle_cos[a2_idx]*r1, angle_sin[a2_idx]*r1, z+zD, 1.0] + + try: + _, n = tri2p0n([sV[start_pole_idx], pp1, pp2]) + except ValueError: + continue + + k1, sV, sN = addVertex(pp1, n, sV, sN) + k2, sV, sN = addVertex(pp2, n, sV, sN) + sF.append([start_pole_idx, k1, k2]) + continue + + if i == steps - 1 and need_end_cap: + # Create triangular faces from last ring to pole + if r0 < epsilon: + r0 = epsilon + for j in range(arcSamples): + a1_idx = j + a2_idx = (j+1) % arcSamples + + p1 = [angle_cos[a1_idx]*r0, angle_sin[a1_idx]*r0, z, 1.0] + p2 = [angle_cos[a2_idx]*r0, angle_sin[a2_idx]*r0, z, 1.0] + + try: + _, n = tri2p0n([p1, sV[end_pole_idx], p2]) + except ValueError: + continue + + k1, sV, sN = addVertex(p1, n, sV, sN) + k2, sV, sN = addVertex(p2, n, sV, sN) + sF.append([k1, end_pole_idx, k2]) + continue + + # Regular quad strips for non-pole sections + if r0 < epsilon: + r0 = epsilon + if r1 < epsilon: + r1 = epsilon + for j in range(arcSamples): - a0 = (j-1)*radStep - a1 = j*radStep - a2 = (j+1)*radStep + # Use pre-computed values with proper wrapping + a0_idx = (j-1) % arcSamples + a1_idx = j + a2_idx = (j+1) % arcSamples - p0 = [math.cos(a0)*r0,math.sin(a0)*r0,z,1.0] - p1 = [math.cos(a1)*r0,math.sin(a1)*r0,z,1.0] - p2 = [math.cos(a2)*r0,math.sin(a2)*r0,z,1.0] + p0 = [angle_cos[a0_idx]*r0, angle_sin[a0_idx]*r0, z, 1.0] + p1 = [angle_cos[a1_idx]*r0, angle_sin[a1_idx]*r0, z, 1.0] + p2 = [angle_cos[a2_idx]*r0, angle_sin[a2_idx]*r0, z, 1.0] - pp1 = [math.cos(a1)*r1,math.sin(a1)*r1,z+zD,1.0] - pp2 = [math.cos(a2)*r1,math.sin(a2)*r1,z+zD,1.0] + pp1 = [angle_cos[a1_idx]*r1, angle_sin[a1_idx]*r1, z+zD, 1.0] + pp2 = [angle_cos[a2_idx]*r1, angle_sin[a2_idx]*r1, z+zD, 1.0] + + try: + p,n = tri2p0n([p0,p2,pp1]) + except ValueError: + # Skip degenerate faces + continue - p,n = tri2p0n([p0,p2,pp1]) - k1,sV,sN = addVertex(p1,n,sV,sN) k2,sV,sN = addVertex(p2,n,sV,sN) k3,sV,sN = addVertex(pp2,n,sV,sN) @@ -575,12 +683,21 @@ def loop_area(loop): def _loft_surface(lower_loop, upper_loop, invert=False): - """Create a surface connecting two XY loops.""" + """Create a surface connecting two loops. + + Args: + lower_loop: List of points forming the lower loop (must be open, not closed) + upper_loop: List of points forming the upper loop (must be open, not closed) + invert: If True, reverse the face winding order + + Returns: + A yapCAD surface connecting the two loops with triangle strips + """ if not lower_loop or not upper_loop: raise ValueError('invalid loops passed to loft surface') - lower = [point(p) for p in lower_loop[:-1]] - upper = [point(p) for p in upper_loop[:-1]] + lower = [point(p) for p in lower_loop] + upper = [point(p) for p in upper_loop] if len(lower) != len(upper): raise ValueError('loop length mismatch in loft surface') diff --git a/tests/test_boolean_regression.py b/tests/test_boolean_regression.py new file mode 100644 index 0000000..e809974 --- /dev/null +++ b/tests/test_boolean_regression.py @@ -0,0 +1,89 @@ +#!/usr/bin/env python3 +""" +Test for Boolean operation regression found in example12.py + +This test documents a regression where Boolean union operations were crashing +due to issurface() and issolid() not properly checking list lengths before +accessing elements. + +The crash has been fixed, but there remains a logic bug in combineglist where +union operations don't correctly include geometry from both operands. +""" + +import pytest +from yapcad.geom import * +from yapcad.poly import * +from yapcad.combine import * + + +class TestBooleanRegression: + """Tests for Boolean operation bugs found in example12.py""" + + def test_boolean_no_crash_on_circles(self): + """Boolean operations on circles should not crash""" + circle1 = Circle(point(0, 0), 10) + circle2 = Circle(point(15, 0), 10) + + # This used to crash with IndexError in issolid() + nb = Boolean('union', [circle1, circle2]) + result = nb.geom + + # Should return something (even if logic is wrong) + assert isinstance(result, list) + + def test_boolean_no_crash_empty_result(self): + """Boolean operations should handle empty results without crashing""" + # Create two non-overlapping, non-contained shapes + circle1 = Circle(point(0, 0), 5) + circle2 = Circle(point(100, 0), 5) + + # Should not crash even if result handling is imperfect + nb = Boolean('union', [circle1, circle2]) + result = nb.geom + + assert isinstance(result, list) + + def test_boolean_union_includes_both_circles(self): + """Boolean union should include geometry from both operands.""" + circle1 = Circle(point(0, 0), 10) + circle2 = Circle(point(15, 0), 10) + + nb = Boolean('union', [circle1, circle2]) + result = nb.geom + + # Boolean operations convert to line segments, not arcs + # Check that result has geometry spanning both circles + assert len(result) > 0, "Result should not be empty" + + # Check that some line segments are near circle1's extent + has_circle1_geom = any( + isline(elem) and (elem[0][0] < 0 or elem[1][0] < 0) # Left of origin + for elem in result + ) + + # Check that some line segments are near circle2's extent + has_circle2_geom = any( + isline(elem) and (elem[0][0] > 20 or elem[1][0] > 20) # Right of circle2 center + for elem in result + ) + + assert has_circle1_geom, "Missing geometry from first circle" + assert has_circle2_geom, "Missing geometry from second circle" + + def test_boolean_union_bbox_spans_both(self): + """Boolean union bbox should span both operands.""" + circle1 = Circle(point(0, 0), 10) + circle2 = Circle(point(15, 0), 10) + + nb = Boolean('union', [circle1, circle2]) + result = nb.geom + + result_bbox = bbox(result) + + # Union should span from x=-10 (left of circle1) to x=25 (right of circle2) + assert result_bbox[0][0] <= -9.9, f"Left bound {result_bbox[0][0]} should be near -10" + assert result_bbox[1][0] >= 24.9, f"Right bound {result_bbox[1][0]} should be near 25" + + +if __name__ == "__main__": + pytest.main([__file__, "-v"]) diff --git a/tests/test_solid_spatial.py b/tests/test_solid_spatial.py new file mode 100644 index 0000000..17a4c3f --- /dev/null +++ b/tests/test_solid_spatial.py @@ -0,0 +1,170 @@ +import copy +import pytest + +from yapcad.geom import point +from yapcad.geom3d import ( + reversesurface, + solid_boolean, + solid_contains_point, + solids_intersect, + translatesolid, + _iter_triangles_from_solid, + tri2p0n, +) +from yapcad.geom3d_util import prism, sphere, conic, tube + + +def _offset_point(base, normal, scale): + return point(base[0] + normal[0] * scale, + base[1] + normal[1] * scale, + base[2] + normal[2] * scale) + + +def _assert_normals_outward(sld, *, label='solid', step=1e-3): + triangles = list(_iter_triangles_from_solid(sld)) + assert triangles, f"{label} produced no faces" + + max_checks = min(32, len(triangles)) + stride = max(1, len(triangles) // max_checks) + checked = 0 + + for idx, tri in enumerate(triangles): + if idx % stride != 0: + continue + try: + center, normal = tri2p0n(tri) + except ValueError: + continue + + scale = step + detected = False + flipped = False + for _ in range(6): + inside_pt = _offset_point(center, normal, -scale) + outside_pt = _offset_point(center, normal, scale) + inside = solid_contains_point(sld, inside_pt) + outside = solid_contains_point(sld, outside_pt) + if inside and not outside: + detected = True + break + if outside and not inside: + flipped = True + break + scale *= 2.0 + if flipped: + raise AssertionError( + f"normal orientation check failed for {label} at center {center[:3]}" + ) + if detected: + checked += 1 + + assert checked > 0, f"no triangles sampled for {label}" + + +def _translate(sld, delta): + return translatesolid(sld, point(delta[0], delta[1], delta[2])) + + +def _make_cavity_solid(): + outer = prism(2, 2, 2) + inner = prism(1, 1, 1) + cavity_surfaces = outer[1] + [reversesurface(s) for s in inner[1]] + return ['solid', cavity_surfaces, [], []] + + +def test_solid_contains_point_basic(): + cube = prism(2, 2, 2) + assert solid_contains_point(cube, point(0.1, 0.1, 0.1)) + assert solid_contains_point(cube, point(1.0, 0.0, 0.0)) + assert not solid_contains_point(cube, point(2.5, 0.0, 0.0)) + + +def test_solid_contains_point_cavity(): + cavity = _make_cavity_solid() + assert not solid_contains_point(cavity, point(0.1, 0.0, 0.0)) + assert solid_contains_point(cavity, point(0.9, 0.0, 0.0)) + + +@pytest.mark.slow +def test_solid_boolean_cubes(): + a = prism(2, 2, 2) + b = _translate(prism(2, 2, 2), (0.75, 0.0, 0.0)) + + union = solid_boolean(a, b, 'union') + assert solid_contains_point(union, point(-0.9, 0.0, 0.0)) + assert solid_contains_point(union, point(1.6, 0.0, 0.0)) + + intersection = solid_boolean(a, b, 'intersection') + assert solid_contains_point(intersection, point(0.6, 0.0, 0.0)) + assert not solid_contains_point(intersection, point(-0.9, 0.0, 0.0)) + + difference = solid_boolean(a, b, 'difference') + assert solid_contains_point(difference, point(-0.9, 0.0, 0.0)) + assert not solid_contains_point(difference, point(0.6, 0.0, 0.0)) + + +@pytest.mark.slow +def test_solid_boolean_spheres(): + a = sphere(2.0) + b = copy.deepcopy(a) + b = translatesolid(b, point(1.0, 0.0, 0.0)) + + union = solid_boolean(a, b, 'union') + assert solid_contains_point(union, point(-0.8, 0.0, 0.0)) + assert solid_contains_point(union, point(1.8, 0.0, 0.0)) + + intersection = solid_boolean(a, b, 'intersection') + assert solid_contains_point(intersection, point(0.5, 0.0, 0.0)) + assert not solid_contains_point(intersection, point(-1.8, 0.0, 0.0)) + + difference = solid_boolean(a, b, 'difference') + assert solid_contains_point(difference, point(-0.8, 0.0, 0.0)) + assert not solid_contains_point(difference, point(0.5, 0.0, 0.0)) + + +@pytest.mark.slow +def test_solids_intersect_detection(): + base = prism(2, 2, 2) + shifted = _translate(prism(2, 2, 2), (0.5, 0.0, 0.0)) + assert solids_intersect(base, shifted) + + far = _translate(prism(2, 2, 2), (10.0, 0.0, 0.0)) + assert not solids_intersect(base, far) + + +@pytest.mark.slow +def test_solid_boolean_sphere_normals(): + a = sphere(2.0) + b = copy.deepcopy(a) + b = translatesolid(b, point(1.0, 0.0, 0.0)) + + union = solid_boolean(a, b, 'union') + _assert_normals_outward(union, label='sphere union') + + difference = solid_boolean(a, b, 'difference') + _assert_normals_outward(difference, label='sphere difference') + + +@pytest.mark.slow +def test_solid_boolean_box_minus_cylinder(): + box = prism(2, 2, 2) + drill = conic(0.6, 0.6, 3.0, center=point(0.0, 0.0, -1.5)) + + result = solid_boolean(box, drill, 'difference') + + assert not solid_contains_point(result, point(0.0, 0.0, 0.0)) + assert solid_contains_point(result, point(0.9, 0.0, 0.0)) + _assert_normals_outward(result, label='box minus cylinder') + + +@pytest.mark.slow +def test_solid_boolean_tube_side_hole(): + shell = tube(outer_diameter=3.0, wall_thickness=0.5, length=4.0, + base_point=point(0.0, 0.0, -2.0)) + drill = conic(0.6, 0.6, 4.2, center=point(1.2, 0.0, -2.1)) + + result = solid_boolean(shell, drill, 'difference') + + assert not solid_contains_point(result, point(1.2, 0.0, 0.0)) + assert solid_contains_point(result, point(0.0, 1.3, 0.0)) + _assert_normals_outward(result, label='tube side hole') diff --git a/tools/analyze_vertices.py b/tools/analyze_vertices.py new file mode 100755 index 0000000..9e1d75c --- /dev/null +++ b/tools/analyze_vertices.py @@ -0,0 +1,231 @@ +#!/usr/bin/env python3 +"""Analyze yapCAD solids for orphaned and duplicated vertices. + +This tool operates directly on yapCAD's internal representation, +analyzing vertex usage and identifying potential geometry issues with +configurable epsilon tolerance for determining vertex proximity. + +Unlike trimesh-based diagnostics, this operates at the yapCAD level +and can detect: +- Truly orphaned vertices (not referenced by any face) +- Near-duplicate vertices (within epsilon distance) +- Vertex usage statistics per surface + +Example: + PYTHONPATH=./src python tools/analyze_vertices.py \\ + --primitive tube \\ + --params '{"outer_diameter": 3.0, "wall_thickness": 0.5, "length": 4.0}' \\ + --epsilon 1e-9 +""" + +from __future__ import annotations + +import argparse +import json +import math +from collections import defaultdict +from dataclasses import asdict, dataclass +from pathlib import Path +from typing import Any, Dict, List, Set, Tuple + +from yapcad.geom import point +from yapcad.geom3d import issolid, issurface +from yapcad.geom3d_util import prism, sphere, conic, tube + + +@dataclass +class VertexAnalysis: + """Results of vertex analysis for a yapCAD solid.""" + total_vertices: int + unique_vertices: int # After de-duplication with epsilon + orphaned_vertices: int # Not referenced by any face + duplicate_groups: int # Number of groups of duplicate vertices + max_duplicates_in_group: int # Largest duplicate group + surfaces_analyzed: int + epsilon: float + + # Per-surface statistics + surface_stats: List[Dict[str, Any]] + + def to_dict(self) -> dict: + return asdict(self) + + +def _vertex_key(v: list, epsilon: float) -> tuple: + """Create a hashable key for a vertex with epsilon precision.""" + if epsilon <= 0: + epsilon = 1e-10 + scale = 1.0 / epsilon + return ( + round(v[0] * scale), + round(v[1] * scale), + round(v[2] * scale), + ) + + +def _vertex_distance(v1: list, v2: list) -> float: + """Calculate Euclidean distance between two vertices.""" + dx = v1[0] - v2[0] + dy = v1[1] - v2[1] + dz = v1[2] - v2[2] + return math.sqrt(dx*dx + dy*dy + dz*dz) + + +def analyze_solid_vertices(solid, epsilon: float = 1e-9) -> VertexAnalysis: + """Analyze vertices in a yapCAD solid. + + Args: + solid: A yapCAD solid object + epsilon: Distance threshold for considering vertices as duplicates + + Returns: + VertexAnalysis object with detailed vertex statistics + """ + if not issolid(solid): + raise ValueError('Input must be a yapCAD solid') + + surfaces = solid[1] # Extract surfaces from solid + + # Global vertex tracking + all_vertices: List[list] = [] + vertex_to_indices: Dict[tuple, List[int]] = defaultdict(list) + referenced_indices: Set[int] = set() + + # Per-surface statistics + surface_stats = [] + + for surf_idx, surf in enumerate(surfaces): + if not issurface(surf): + continue + + surf_vertices = surf[1] + surf_faces = surf[3] + + # Track which local vertices are referenced by faces + local_referenced = set() + for face in surf_faces: + for vertex_idx in face: + local_referenced.add(vertex_idx) + + local_orphaned = len(surf_vertices) - len(local_referenced) + + # Add vertices to global list + start_idx = len(all_vertices) + for local_idx, v in enumerate(surf_vertices): + global_idx = start_idx + local_idx + all_vertices.append(v) + + # Track by epsilon-rounded key + key = _vertex_key(v, epsilon) + vertex_to_indices[key].append(global_idx) + + # Track if this vertex is referenced + if local_idx in local_referenced: + referenced_indices.add(global_idx) + + surface_stats.append({ + 'surface_index': surf_idx, + 'vertices': len(surf_vertices), + 'faces': len(surf_faces), + 'orphaned_vertices': local_orphaned, + 'referenced_vertices': len(local_referenced), + }) + + # Analyze duplicates + duplicate_groups = 0 + max_duplicates = 0 + + for key, indices in vertex_to_indices.items(): + if len(indices) > 1: + # Verify they're actually close (not just hash collisions) + verified_group = [indices[0]] + for idx in indices[1:]: + if any(_vertex_distance(all_vertices[idx], all_vertices[ref_idx]) < epsilon + for ref_idx in verified_group): + verified_group.append(idx) + + if len(verified_group) > 1: + duplicate_groups += 1 + max_duplicates = max(max_duplicates, len(verified_group)) + + total_vertices = len(all_vertices) + orphaned_vertices = total_vertices - len(referenced_indices) + unique_vertices = total_vertices - (sum(len(indices) - 1 for indices in vertex_to_indices.values() + if len(indices) > 1)) + + return VertexAnalysis( + total_vertices=total_vertices, + unique_vertices=unique_vertices, + orphaned_vertices=orphaned_vertices, + duplicate_groups=duplicate_groups, + max_duplicates_in_group=max_duplicates, + surfaces_analyzed=len(surfaces), + epsilon=epsilon, + surface_stats=surface_stats, + ) + + +def _build_primitive(name: str, params: Dict[str, Any]): + """Build a yapCAD primitive from name and parameters.""" + if name == 'sphere': + diameter = float(params.get('diameter', 2.0)) + depth = int(params.get('depth', 2)) + return sphere(diameter, depth=depth) + if name == 'prism': + length = float(params.get('length', 2.0)) + width = float(params.get('width', 2.0)) + height = float(params.get('height', 2.0)) + return prism(length, width, height) + if name == 'conic': + baser = float(params.get('base_radius', 1.0)) + topr = float(params.get('top_radius', 0.5)) + height = float(params.get('height', 2.0)) + angr = float(params.get('angular_resolution', 10.0)) + return conic(baser, topr, height, center=point(0, 0, 0), angr=angr) + if name == 'tube': + outer_diameter = float(params.get('outer_diameter', 3.0)) + wall_thickness = float(params.get('wall_thickness', 0.5)) + length = float(params.get('length', 4.0)) + include_caps = bool(params.get('include_caps', True)) + base_point = params.get('base_point') + if base_point is not None: + base_point = point(*base_point) + return tube(outer_diameter=outer_diameter, wall_thickness=wall_thickness, + length=length, base_point=base_point, include_caps=include_caps) + raise ValueError(f'unsupported primitive {name!r}') + + +def main(): + parser = argparse.ArgumentParser( + description='Analyze yapCAD solid vertices for orphans and duplicates.', + formatter_class=argparse.RawDescriptionHelpFormatter, + epilog=__doc__, + ) + parser.add_argument('--primitive', choices=['sphere', 'prism', 'conic', 'tube'], required=True) + parser.add_argument('--params', help='JSON object with primitive parameters') + parser.add_argument('--epsilon', type=float, default=1e-9, + help='Distance threshold for duplicate detection (default: 1e-9)') + parser.add_argument('--json-output', type=Path, help='Write JSON results to file') + + args = parser.parse_args() + + params: Dict[str, Any] = {} + if args.params: + try: + params = json.loads(args.params) + except json.JSONDecodeError as exc: + raise SystemExit(f'Failed to parse --params JSON: {exc}') from exc + + solid = _build_primitive(args.primitive, params) + analysis = analyze_solid_vertices(solid, epsilon=args.epsilon) + + if args.json_output: + args.json_output.parent.mkdir(parents=True, exist_ok=True) + args.json_output.write_text(json.dumps(analysis.to_dict(), indent=2)) + print(f'Results written to {args.json_output}') + else: + print(json.dumps(analysis.to_dict(), indent=2)) + + +if __name__ == '__main__': + main() diff --git a/tools/diagnose_boolean.py b/tools/diagnose_boolean.py new file mode 100755 index 0000000..dcb442d --- /dev/null +++ b/tools/diagnose_boolean.py @@ -0,0 +1,66 @@ +#!/usr/bin/env python3 +"""Diagnose boolean meshes produced by yapCAD using trimesh helpers.""" + +from __future__ import annotations + +import argparse +import json +from pathlib import Path + +from mesh_diagnostics import MeshDiagnostics, diagnose_mesh, export_boundary + +from yapcad.geom import point +from yapcad.geom3d import solid_boolean, translatesolid +from yapcad.geom3d_util import prism, sphere, conic, tube + +from yapcad.boolean.trimesh_engine import _solid_to_mesh + + +def _build_solids(kind: str): + if kind == 'cube': + a = prism(2, 2, 2) + b = translatesolid(prism(2, 2, 2), point(0.75, 0.0, 0.0)) + elif kind == 'sphere': + a = sphere(2.0) + b = translatesolid(sphere(2.0), point(1.0, 0.0, 0.0)) + elif kind == 'box_hole': + a = prism(2, 2, 2) + b = conic(0.6, 0.6, 3.0, center=point(0.0, 0.0, -1.5)) + elif kind == 'tube_hole': + a = tube(outer_diameter=3.0, wall_thickness=0.5, length=4.0, + base_point=point(0.0, 0.0, -2.0)) + b = conic(0.6, 0.6, 4.2, center=point(1.2, 0.0, -2.1)) + else: + raise ValueError(f'unsupported shape kind: {kind!r}') + return a, b + + +def main(): + parser = argparse.ArgumentParser(description='Diagnose boolean mesh quality.') + parser.add_argument('--shapes', choices=['cube', 'sphere', 'box_hole', 'tube_hole'], required=True) + parser.add_argument('--operation', choices=['union', 'intersection', 'difference'], default='union') + parser.add_argument('--engine', default=None, help='boolean engine (e.g., native, trimesh:manifold)') + parser.add_argument('--stitch', action='store_true', help='enable stitch flag for boolean call') + parser.add_argument('--boundary-output', type=Path, help='optional path to export boundary edges (PLY)') + parser.add_argument('--json', type=Path, help='write diagnostics JSON to path') + args = parser.parse_args() + + solids = _build_solids(args.shapes) + result = solid_boolean(solids[0], solids[1], args.operation, + stitch=args.stitch, engine=args.engine) + + mesh = _solid_to_mesh(result) + diagnostics = diagnose_mesh(mesh) + + if args.boundary_output: + export_boundary(mesh, args.boundary_output) + + print(json.dumps(diagnostics.to_dict(), indent=2)) + + if args.json: + args.json.parent.mkdir(parents=True, exist_ok=True) + args.json.write_text(json.dumps(diagnostics.to_dict(), indent=2)) + + +if __name__ == '__main__': + main() diff --git a/tools/diagnose_primitive.py b/tools/diagnose_primitive.py new file mode 100755 index 0000000..8e96491 --- /dev/null +++ b/tools/diagnose_primitive.py @@ -0,0 +1,82 @@ +#!/usr/bin/env python3 +"""Diagnose individual yapCAD primitives using trimesh helpers.""" + +from __future__ import annotations + +import argparse +import json +from pathlib import Path +from typing import Any, Dict + +from mesh_diagnostics import diagnose_mesh, export_boundary + +from yapcad.geom import point +from yapcad.geom3d_util import prism, sphere, conic, tube +from yapcad.boolean.trimesh_engine import _solid_to_mesh + + +def _build_primitive(name: str, params: Dict[str, Any]): + if name == 'sphere': + diameter = float(params.get('diameter', 2.0)) + depth = int(params.get('depth', 2)) + return sphere(diameter, depth=depth) + if name == 'prism': + length = float(params.get('length', 2.0)) + width = float(params.get('width', 2.0)) + height = float(params.get('height', 2.0)) + return prism(length, width, height) + if name == 'conic': + baser = float(params.get('base_radius', 1.0)) + topr = float(params.get('top_radius', 0.5)) + height = float(params.get('height', 2.0)) + angr = float(params.get('angular_resolution', 10.0)) + return conic(baser, topr, height, center=params.get('center', point(0, 0, 0)), angr=angr) + if name == 'tube': + outer_diameter = float(params.get('outer_diameter', 3.0)) + wall_thickness = float(params.get('wall_thickness', 0.5)) + length = float(params.get('length', 4.0)) + include_caps = bool(params.get('include_caps', True)) + base_point = params.get('base_point') + if base_point is not None: + base_point = point(*base_point) + return tube(outer_diameter=outer_diameter, wall_thickness=wall_thickness, + length=length, base_point=base_point, include_caps=include_caps) + raise ValueError(f'unsupported primitive {name!r}') + + +def main(): + parser = argparse.ArgumentParser(description='Diagnose individual yapCAD primitives.') + parser.add_argument('--primitive', choices=['sphere', 'prism', 'conic', 'tube'], required=True) + parser.add_argument('--params', help='JSON object with primitive parameters (e.g. {"diameter": 3.0})') + parser.add_argument('--boundary-output', type=Path, help='optional path to export boundary edges (PLY)') + parser.add_argument('--json', type=Path, help='write diagnostics JSON to path') + parser.add_argument('--export-stl', type=Path, help='optional path to export the mesh as STL') + args = parser.parse_args() + + params: Dict[str, Any] = {} + if args.params: + try: + params = json.loads(args.params) + except json.JSONDecodeError as exc: + raise SystemExit(f'Failed to parse --params JSON: {exc}') from exc + + solid = _build_primitive(args.primitive, params) + mesh = _solid_to_mesh(solid) + diagnostics = diagnose_mesh(mesh) + + if args.boundary_output: + export_boundary(mesh, args.boundary_output) + + if args.export_stl: + args.export_stl.parent.mkdir(parents=True, exist_ok=True) + mesh.export(str(args.export_stl)) + + print(json.dumps(diagnostics.to_dict(), indent=2)) + + if args.json: + args.json.parent.mkdir(parents=True, exist_ok=True) + args.json.write_text(json.dumps(diagnostics.to_dict(), indent=2)) + + +if __name__ == '__main__': + main() diff --git a/tools/find_func.py b/tools/find_func.py new file mode 100644 index 0000000..8674967 --- /dev/null +++ b/tools/find_func.py @@ -0,0 +1,42 @@ +import ast +import inspect +import sys + +def get_function_source(filename, function_name): + """Prints the source code of a specified function.""" + try: + with open(filename, "r") as source_file: + source_code = source_file.read() + except FileNotFoundError: + print(f"Error: The file '{filename}' was not found.") + return + + tree = ast.parse(source_code) + found_func = False + + for node in ast.walk(tree): + if isinstance(node, (ast.FunctionDef, ast.AsyncFunctionDef)): + if node.name == function_name: + found_func = True + + # Extract the function body using line numbers from AST + lines = source_code.splitlines() + start_line = node.lineno - 1 + end_line = node.end_lineno + + # Print the extracted lines + for line in lines[start_line:end_line]: + print(line) + return + + if not found_func: + print(f"Function '{function_name}' not found in {filename}.") + +if __name__ == "__main__": + if len(sys.argv) < 3: + print("Usage: python find_func.py ") + sys.exit(1) + + filename = sys.argv[1] + function_name = sys.argv[2] + get_function_source(filename, function_name) diff --git a/tools/mesh_diagnostics.py b/tools/mesh_diagnostics.py new file mode 100644 index 0000000..676adcf --- /dev/null +++ b/tools/mesh_diagnostics.py @@ -0,0 +1,107 @@ +"""Shared mesh-diagnostics helpers for yapCAD solids.""" + +from __future__ import annotations + +from dataclasses import asdict, dataclass +from pathlib import Path +from typing import Optional + +try: + import numpy as np +except ImportError as exc: # pragma: no cover - optional dependency + raise SystemExit("numpy is required for mesh diagnostics (pip install numpy)") from exc + +try: + import trimesh +except ImportError as exc: # pragma: no cover - optional dependency + raise SystemExit("trimesh is required for mesh diagnostics (pip install trimesh[easy])") from exc + + +@dataclass +class MeshDiagnostics: + triangles: int + vertices: int + is_watertight: bool + is_volume: bool + euler_number: int + boundary_edges: int + self_intersections: int + min_edge_length: Optional[float] + max_edge_length: Optional[float] + orphan_vertices: int + + def to_dict(self) -> dict: + return asdict(self) + + +def _boundary_mask(mesh: "trimesh.Trimesh") -> Optional[np.ndarray]: + try: + unique = mesh.edges_unique + inverse = mesh.edges_unique_inverse + if unique is None or inverse is None: + return None + counts = np.bincount(inverse, minlength=len(unique)) + if counts.size == 0: + return np.zeros(0, dtype=bool) + return counts == 1 + except Exception: + return None + + +def diagnose_mesh(mesh: "trimesh.Trimesh") -> MeshDiagnostics: + boundary_mask = _boundary_mask(mesh) + boundary_count = int(boundary_mask.sum()) if boundary_mask is not None else 0 + + self_intersections = 0 + try: + faces = mesh.self_intersecting_faces + if faces is not None: + self_intersections = int(len(faces)) + except Exception: + pass + + min_len = max_len = None + try: + lengths = mesh.edges_unique_length + if lengths is not None and len(lengths) > 0: + min_len = float(np.min(lengths)) + max_len = float(np.max(lengths)) + except Exception: + pass + + orphan_vertices = 0 + try: + vf = mesh.vertex_faces + if vf is not None: + orphan_vertices = int((vf < 0).sum()) + except Exception: + orphan_vertices = 0 + + is_volume = bool(getattr(mesh, "is_volume", mesh.is_watertight)) + + return MeshDiagnostics( + triangles=int(len(mesh.faces)) if mesh.faces is not None else 0, + vertices=int(len(mesh.vertices)) if mesh.vertices is not None else 0, + is_watertight=bool(mesh.is_watertight), + is_volume=is_volume, + euler_number=int(getattr(mesh, "euler_number", 0)), + boundary_edges=boundary_count, + self_intersections=int(self_intersections), + min_edge_length=min_len, + max_edge_length=max_len, + orphan_vertices=int(orphan_vertices), + ) + + +def export_boundary(mesh: "trimesh.Trimesh", path: Path) -> Optional[Path]: + mask = _boundary_mask(mesh) + if mask is None or not mask.any(): + return None + segments = mesh.edges_unique[mask] + if segments.size == 0: + return None + line = trimesh.load_path(mesh.vertices[segments]) + path.parent.mkdir(parents=True, exist_ok=True) + line.export(str(path)) + return path + diff --git a/tools/validate_mesh.py b/tools/validate_mesh.py new file mode 100644 index 0000000..34a3f1e --- /dev/null +++ b/tools/validate_mesh.py @@ -0,0 +1,353 @@ +#!/usr/bin/env python3 +"""Run CAM-oriented linting on triangular meshes. + +This script chains together three pragmatic checks that mirror the +downstream tools typically used in fabrication workflows: + +1. ``admesh`` for statistics about the STL (holes, connected components, + facet orientation, volume, bounding box, etc.). +2. ``meshfix`` for geometry healing (fills small holes, removes dangling + triangles, produces a cleaned STL for further work). +3. A slicer front-end such as ``prusa-slicer``/``prusa-slicer-console``/ + ``slic3r``/``curaengine`` to prove the model slices successfully. + +Each step is optional – if a tool is not installed the script will emit a +warning and carry on. The goal is to provide a quick report that correlates +better with real-world CAM robustness than the strict, topological +``issolidclosed`` test. + +Example +======= + +.. code-block:: bash + + PYTHONPATH=src python examples/solid_boolean_demo.py \ + --mode stl --operation difference --shapes box_hole \ + --output box_hole_difference + + python tools/validate_mesh.py box_hole_difference.stl \ + --workdir build/mesh_checks + +The command above keeps all intermediate files (repaired STL, generated +G-code, log snippets) under ``build/mesh_checks`` so they can be inspected +or archived alongside the STEP export. +""" + +from __future__ import annotations + +import argparse +import json +import os +import shutil +import subprocess +import sys +import tempfile +from pathlib import Path +from typing import Iterable, Optional + + +def _run_pymeshfix(stl: Path, output: Path) -> Optional[dict[str, object]]: + """Attempt to repair ``stl`` using the PyMeshFix Python bindings. + + Returns a summary dictionary if successful, otherwise ``None`` when the + dependency is missing or the mesh could not be loaded. + """ + + try: + import pymeshfix # type: ignore[import-untyped] + import trimesh # type: ignore[import-untyped] + except ImportError: + return None + + try: + mesh = trimesh.load(str(stl), force='mesh') + except Exception as exc: # pragma: no cover - depends on trimesh internals + return { + "returncode": 1, + "stdout": "", + "stderr": f"Failed to load STL with trimesh: {exc}", + } + + if not isinstance(mesh, trimesh.Trimesh) or mesh.vertices.size == 0: + return { + "returncode": 1, + "stdout": "", + "stderr": "Input did not contain a valid triangular mesh", + } + + fixer = pymeshfix.MeshFix(mesh.vertices, mesh.faces) + fixer.repair(verbose=False) + try: + fixer.write_stl(str(output)) + except AttributeError: + # Older PyMeshFix releases expose the repaired data via ``v``/``f``. + repaired = trimesh.Trimesh(fixer.v, fixer.f, process=False) + repaired.export(str(output)) + return { + "returncode": 0, + "stdout": "PyMeshFix repair completed", + "stderr": "", + "output": str(output), + } + + +ToolName = str + + +def which(names: Iterable[ToolName]) -> Optional[Path]: + """Return the first executable in *names* that exists on PATH.""" + + for name in names: + path = shutil.which(name) + if path: + return Path(path) + return None + + +def run(cmd: Iterable[str], *, cwd: Optional[Path] = None) -> subprocess.CompletedProcess: + """Execute *cmd*, capturing stdout/stderr for later inspection.""" + + result = subprocess.run( + list(cmd), + cwd=str(cwd) if cwd else None, + check=False, + text=True, + encoding="utf-8", + errors="replace", + capture_output=True, + ) + return result + + +def parse_admesh_report(output: str) -> dict[str, str]: + """Extract a few useful numbers from ``admesh`` output.""" + + summary: dict[str, str] = {} + for line in output.splitlines(): + line = line.strip() + if not line: + continue + for prefix in ( + "Number of facets:", + "Volume:", + "Number of shells:", + "Number of pieces of geometry:", + "Number of edges (with 2 facets):", + "Number of edges (with 1 facet):", + ): + if line.startswith(prefix): + summary[prefix.rstrip(':')] = line.split(':', 1)[1].strip() + return summary + + +def slicer_command(output_dir: Path, stl: Path) -> Optional[list[str]]: + """Return a slicer command if a supported engine is installed.""" + + slicers: list[tuple[list[str], list[str]]] = [ + (["prusa-slicer-console", "prusa-slicer", "PrusaSlicer"], ["--help"]), + (["slic3r"], ["--help"]), + (["curaengine"], ["--help"]), + ] + + for names, probe_args in slicers: + binary = which(names) + if not binary: + continue + + # Make sure the executable is functional by running the probe command. + probe = run([str(binary), *probe_args]) + if probe.returncode != 0: + continue + + gcode = output_dir / (stl.stem + ".gcode") + + if binary.name.startswith("curaengine"): + profile = os.environ.get("CURA_PROFILE") + if not profile: + continue + return [str(binary), "slice", "-v", "-o", str(gcode), "-j", profile, "-l", str(stl)] + + if binary.name.startswith("slic3r"): + return [ + str(binary), + str(stl), + "--output", + str(gcode), + "--loglevel", + "0", + ] + + # Default for PrusaSlicer variants. + return [ + str(binary), + "--export-gcode", + str(stl), + "--output", + str(gcode), + ] + + return None + + +def main(argv: list[str] | None = None) -> int: + parser = argparse.ArgumentParser( + description="Run mesh linting/healing/slicing checks in one pass.", + formatter_class=argparse.ArgumentDefaultsHelpFormatter, + ) + parser.add_argument("stl", type=Path, help="Path to the STL file to validate") + parser.add_argument( + "--workdir", + type=Path, + help="Directory to store intermediate files (default: temporary directory).", + ) + parser.add_argument( + "--json", + action="store_true", + help="Emit a machine-readable JSON summary instead of human text.", + ) + + args = parser.parse_args(argv) + stl_path = args.stl + + if not stl_path.exists(): + parser.error(f"STL file not found: {stl_path}") + + if not stl_path.suffix.lower() == ".stl": + parser.error("The input file must be an STL") + + if args.workdir: + workdir = args.workdir + workdir.mkdir(parents=True, exist_ok=True) + cleanup = False + else: + workdir = Path(tempfile.mkdtemp(prefix="yapcad-meshcheck-")) + cleanup = True + + summary: dict[str, object] = { + "input": str(stl_path.resolve()), + "workdir": str(workdir.resolve()), + "admesh": None, + "meshfix": None, + "slicer": None, + "notes": [], + } + + # 1. admesh ------------------------------------------ + admesh_bin = which(["admesh"]) + if admesh_bin: + repaired_stl = workdir / f"{stl_path.stem}-admesh.stl" + cmd = [str(admesh_bin), str(stl_path)] + result = run(cmd) + admesh_info = { + "returncode": result.returncode, + "stdout": result.stdout, + "stderr": result.stderr, + } + admesh_info.update(parse_admesh_report(result.stdout)) + summary["admesh"] = admesh_info + else: + summary["notes"].append("admesh not found on PATH; skipping mesh statistics") + + # 2. meshfix ----------------------------------------- + meshfix_bin = which(["meshfix"]) + if meshfix_bin: + meshfix_out = workdir / f"{stl_path.stem}-meshfix.stl" + cmd = [str(meshfix_bin), str(stl_path), "-o", str(meshfix_out)] + result = run(cmd) + summary["meshfix"] = { + "returncode": result.returncode, + "stdout": result.stdout, + "stderr": result.stderr, + "output": str(meshfix_out), + } + if result.returncode != 0: + summary["notes"].append("meshfix reported an error; see stderr for details") + else: + meshfix_out = workdir / f"{stl_path.stem}-pymeshfix.stl" + pymeshfix_info = _run_pymeshfix(stl_path, meshfix_out) + if pymeshfix_info is not None: + summary["meshfix"] = pymeshfix_info + if pymeshfix_info.get("returncode") != 0: + summary["notes"].append("PyMeshFix reported an issue; see stderr for details") + else: + summary["notes"].append("meshfix not found and PyMeshFix unavailable; skipping repair step") + + # 3. slicer ------------------------------------------ + slicer_cmd = slicer_command(workdir, stl_path) + if slicer_cmd: + result = run(slicer_cmd) + summary["slicer"] = { + "command": slicer_cmd, + "returncode": result.returncode, + "stdout": result.stdout, + "stderr": result.stderr, + } + if result.returncode != 0: + summary["notes"].append("Slicer failed – inspect stdout/stderr for diagnostics") + else: + summary["notes"].append( + "No supported slicer (prusa-slicer, slic3r, curaengine) found; skipping G-code export" + ) + + if args.json: + print(json.dumps(summary, indent=2)) + else: + lines = [f"Mesh validation report for {stl_path}"] + lines.append("Output directory: " + str(workdir)) + lines.append("") + + if summary["admesh"]: + lines.append("admesh summary:") + admesh_info = summary["admesh"] # type: ignore[assignment] + for key in ( + "Number of facets", + "Number of facets connected to:", + "Number of edges (with 1 facet)", + "Number of edges (with 2 facets)", + "Number of shells", + "Volume", + ): + value = admesh_info.get(key) + if value is not None: + lines.append(f" {key}: {value}") + lines.append("") + else: + lines.append("admesh: skipped") + lines.append("") + + if summary["meshfix"]: + meshfix_info = summary["meshfix"] # type: ignore[assignment] + lines.append("meshfix: return code %s" % meshfix_info["returncode"]) + lines.append(" output: %s" % meshfix_info["output"]) + lines.append("") + else: + lines.append("meshfix: skipped") + lines.append("") + + if summary["slicer"]: + slicer_info = summary["slicer"] # type: ignore[assignment] + cmd_str = " ".join(str(x) for x in slicer_info.get("command", [])) + lines.append("slicer command: " + cmd_str) + lines.append(" return code: %s" % slicer_info["returncode"]) + lines.append("") + else: + lines.append("slicer: skipped") + lines.append("") + + if summary["notes"]: + lines.append("Notes:") + for note in summary["notes"]: + lines.append(" - " + note) + else: + lines.append("No warnings.") + + print("\n".join(lines)) + + if cleanup: + shutil.rmtree(workdir, ignore_errors=True) + + return 0 + + +if __name__ == "__main__": + sys.exit(main())