|
| 1 | +#!/usr/bin/env python3 |
| 2 | +"""Independent STL topology/bounds check for tq-threads release proofs.""" |
| 3 | + |
| 4 | +from __future__ import annotations |
| 5 | + |
| 6 | +import argparse |
| 7 | +import math |
| 8 | +import struct |
| 9 | +import sys |
| 10 | +from collections import Counter |
| 11 | +from pathlib import Path |
| 12 | + |
| 13 | + |
| 14 | +def _read_stl(path: Path) -> list[tuple[tuple[float, float, float], ...]]: |
| 15 | + data = path.read_bytes() |
| 16 | + tris: list[tuple[tuple[float, float, float], ...]] = [] |
| 17 | + |
| 18 | + if len(data) >= 84: |
| 19 | + n = struct.unpack_from("<I", data, 80)[0] |
| 20 | + if 84 + 50 * n == len(data): |
| 21 | + off = 84 |
| 22 | + for _ in range(n): |
| 23 | + vals = struct.unpack_from("<12fH", data, off) |
| 24 | + tris.append((vals[3:6], vals[6:9], vals[9:12])) |
| 25 | + off += 50 |
| 26 | + return tris |
| 27 | + |
| 28 | + verts: list[tuple[float, float, float]] = [] |
| 29 | + for line in data.decode("utf-8", errors="ignore").splitlines(): |
| 30 | + parts = line.strip().split() |
| 31 | + if len(parts) == 4 and parts[0] == "vertex": |
| 32 | + verts.append(tuple(float(x) for x in parts[1:4])) |
| 33 | + if len(verts) % 3: |
| 34 | + raise ValueError(f"{path}: ASCII STL has an incomplete triangle") |
| 35 | + for i in range(0, len(verts), 3): |
| 36 | + tris.append((verts[i], verts[i + 1], verts[i + 2])) |
| 37 | + return tris |
| 38 | + |
| 39 | + |
| 40 | +def _key(v: tuple[float, float, float], scale: float) -> tuple[int, int, int]: |
| 41 | + return tuple(round(c / scale) for c in v) |
| 42 | + |
| 43 | + |
| 44 | +def analyze(path: Path, tol: float) -> dict[str, float | int | bool]: |
| 45 | + tris = _read_stl(path) |
| 46 | + if not tris: |
| 47 | + raise ValueError(f"{path}: no triangles found") |
| 48 | + |
| 49 | + edges: Counter[tuple[tuple[int, int, int], tuple[int, int, int]]] = Counter() |
| 50 | + xs: list[float] = [] |
| 51 | + ys: list[float] = [] |
| 52 | + zs: list[float] = [] |
| 53 | + |
| 54 | + for tri in tris: |
| 55 | + keys = [_key(v, tol) for v in tri] |
| 56 | + for a, b in ((0, 1), (1, 2), (2, 0)): |
| 57 | + edge = tuple(sorted((keys[a], keys[b]))) |
| 58 | + edges[edge] += 1 |
| 59 | + for x, y, z in tri: |
| 60 | + xs.append(x) |
| 61 | + ys.append(y) |
| 62 | + zs.append(z) |
| 63 | + |
| 64 | + bad_edges = sum(1 for count in edges.values() if count != 2) |
| 65 | + max_radius = max(math.hypot(x, y) for x, y in zip(xs, ys)) |
| 66 | + min_radius = min(math.hypot(x, y) for x, y in zip(xs, ys) if math.hypot(x, y) > tol) |
| 67 | + |
| 68 | + return { |
| 69 | + "triangles": len(tris), |
| 70 | + "unique_edges": len(edges), |
| 71 | + "bad_edges": bad_edges, |
| 72 | + "manifold": bad_edges == 0, |
| 73 | + "diameter": 2 * max_radius, |
| 74 | + "min_diameter": 2 * min_radius, |
| 75 | + "height": max(zs) - min(zs), |
| 76 | + "z_min": min(zs), |
| 77 | + "z_max": max(zs), |
| 78 | + } |
| 79 | + |
| 80 | + |
| 81 | +def main() -> int: |
| 82 | + parser = argparse.ArgumentParser() |
| 83 | + parser.add_argument("stl", type=Path) |
| 84 | + parser.add_argument("--tol", type=float, default=1e-5) |
| 85 | + parser.add_argument("--expect-manifold", action="store_true") |
| 86 | + parser.add_argument("--max-diameter", type=float) |
| 87 | + parser.add_argument("--min-diameter", type=float) |
| 88 | + parser.add_argument("--expect-height", type=float) |
| 89 | + parser.add_argument("--height-tol", type=float, default=0.08) |
| 90 | + args = parser.parse_args() |
| 91 | + |
| 92 | + result = analyze(args.stl, args.tol) |
| 93 | + for key, value in result.items(): |
| 94 | + print(f"{key}: {value}") |
| 95 | + |
| 96 | + failures: list[str] = [] |
| 97 | + if args.expect_manifold and not result["manifold"]: |
| 98 | + failures.append(f"expected a closed 2-manifold, found {result['bad_edges']} bad edges") |
| 99 | + if args.max_diameter is not None and result["diameter"] > args.max_diameter: |
| 100 | + failures.append(f"diameter {result['diameter']:.6f} > {args.max_diameter:.6f}") |
| 101 | + if args.min_diameter is not None and result["diameter"] < args.min_diameter: |
| 102 | + failures.append(f"diameter {result['diameter']:.6f} < {args.min_diameter:.6f}") |
| 103 | + if args.expect_height is not None: |
| 104 | + delta = abs(result["height"] - args.expect_height) |
| 105 | + if delta > args.height_tol: |
| 106 | + failures.append(f"height {result['height']:.6f} differs from {args.expect_height:.6f} by {delta:.6f}") |
| 107 | + |
| 108 | + if failures: |
| 109 | + for failure in failures: |
| 110 | + print(f"FAIL: {failure}", file=sys.stderr) |
| 111 | + return 1 |
| 112 | + print("PASS") |
| 113 | + return 0 |
| 114 | + |
| 115 | + |
| 116 | +if __name__ == "__main__": |
| 117 | + raise SystemExit(main()) |
0 commit comments