From e1a9d9bd338a189da8ab6619a53dfbbd831c4df5 Mon Sep 17 00:00:00 2001 From: Emmanuel FARHI Date: Sun, 18 Oct 2026 20:08:58 +0900 Subject: [PATCH 1/8] mcdisplay:cad: rewrite with more output formats --- tools/Python/mcdisplay/cad/mcdisplay.py | 1078 ++++++++++++++++++++--- 1 file changed, 936 insertions(+), 142 deletions(-) diff --git a/tools/Python/mcdisplay/cad/mcdisplay.py b/tools/Python/mcdisplay/cad/mcdisplay.py index 5063949900..55f8a9b4ba 100755 --- a/tools/Python/mcdisplay/cad/mcdisplay.py +++ b/tools/Python/mcdisplay/cad/mcdisplay.py @@ -1,159 +1,951 @@ #!/usr/bin/env python3 # -*- coding: utf-8 -*- ''' -Converts MCDISPLAY geometry information, read from an instrument binary, to STEP format. -Uses the CadQuery package, https://cadquery.readthedocs.io, and its interface to -OpenCASCADE to handle some operations, most importantly, the saving to STEP file. +Converts MCDISPLAY geometry information, read from an instrument binary, +to common CAD / 3D formats WITHOUT CadQuery / OpenCASCADE. + +Supported output formats: glTF (.gltf), GLB, OBJ (+MTL), STL, PLY, VRML97 (.wrl), +X3D (also as legacy 'xml'), VTK legacy (.vtk), STEP (AP214 faceted BREP, pure +Python), and VTKJS (.vtkjs zip of the vtk.js JSON scene format). +No CadQuery/OpenCASCADE is required for any format. + +Each instrument component is assigned a color from a cycling palette, so +adjacent components are visually distinguishable in the output. +Color-capable formats: glTF, GLB, OBJ (via companion .mtl), PLY (vertex colors), +VRML97, X3D. STL and legacy VTK carry no color information. + +Only dependency: numpy (already required by mccodelib.instrgeom). + +All primitives are tessellated in pure Python/numpy: + - lines / dashed lines / multilines -> polylines (wireframe) + - circles / new circles -> polylines + - rectangles -> quad (mesh + outline) + - boxes -> 12-triangle cube + edges + - spheres, cones, cylinders, discs, annuli -> triangulated meshes + - polyhedron -> triangulated faces (fan) +Every draw call is placed in global coordinates via the component m4 matrix. ''' -__author__ = "Gregory Tucker" -__date__ = "2022-11-11" +__author__ = "Gregory Tucker (CadQuery version); dependency-free version" +__date__ = "2026-10-17" import sys +import math +import struct from pathlib import Path import subprocess sys.path.append(str(Path(__file__).resolve().parent.parent.parent)) -from mccodelib.instrgeom import Vector3d, DrawLine, DrawMultiline, DrawCircle +from mccodelib.instrgeom import (Vector3d, DrawMagnify, DrawLine, DrawDashedLine, + DrawMultiline, DrawRectangle, DrawBox, DrawCircle, + DrawNewCircle, DrawSphere, DrawCone, DrawCylinder, + DrawDisc, DrawAnnulus, DrawPolygon, DrawPolyhedron) from mccodelib import mccode_config +import numpy as np -def vector_to_tuple(v: Vector3d): - return v.x, v.y, v.z - - -def transform_shape(global_transformation_matrix, local_shape): - """Perform the coordinate transformation for a shape - - >>> placed = transform_shape(global_transformation_matrix, local_shape) - - Including handling of various possible input shape types and transformation types - """ - # Import types which the global transformation matrix could be: - from cadquery.occ_impl.geom import Matrix - from OCP.gp import gp_Trsf, gp_GTrsf # plus more transformatation types? - # Import types which the local shape could be: - from cadquery.occ_impl.shapes import Compound, Shape - from cadquery.cq import Workplane - from OCP.TopoDS import TopoDS_Compound, TopoDS_Shape - - # Get a handle to the underlying OpenCascade shape object, depending on what the input type is - global_shape = None - if isinstance(local_shape, Workplane): - global_shape = local_shape.toOCC() - elif isinstance(local_shape, (Shape, Compound)): - global_shape = local_shape.wrapped - elif insinstance(local_shape, (TopoDS_Shape, TopoDS_Compound)): - # We are already a suitable type for local to global transformation - global_shape = local_shape - else: - print(f"Unknown shape type {type(local_shape)}, attempting the 'wrapped' property, expect errors!") - global_shape = local_shape.wrapped - - transformation = None - if isinstance(global_transformation_matrix, Matrix): - transformation = global_transformation_matrix.wrapped - elif isinstance(global_transformation_matrix, (gp_Trsf, gp_GTrsf)): - # Already a suitable transformation type - transformation = global_transformation_matrix - else: - gtmt = type(global_transformation_matrix) - print(f"Unknown transformation implementation type {gtmt}, attemping the 'wrapped' property, expect errors!") - transformation = global_transformation_matrix.wrapped - - # Now decide which transformation algorithm can be used based on what the input transformation represents - if isinstance(transformation, gp_Trsf): - # This is a unitary transformation and can only reposition/reorient the shape: - from OCP.BRepBuilderAPI import BRepBuilderAPI_Transform as Trans - global_shape = Trans(global_shape, transformation).Shape() - elif isinstance(transformation, gp_GTrsf): - # This is a general transformation which can additionally rescale/stretch/skew the shape: - from OCP.BRepBuilderAPI import BRepBuilderAPI_GTransform as GTrans - global_shape = GTrans(global_shape, transformation).Shape() - else: - print(f"Unknown transformation type {type(matrix)}. No transformation applied!") - - # (Re)wrap the OpenCascade shape back to a cadquery.Shape - shape = Shape(global_shape) - return shape - - -def circle_to_brep(comp: DrawCircle): - from cadquery import Workplane - plane = comp.plane.upper() - center = vector_to_tuple(comp.center) - radius = comp.radius # / meter - w = Workplane(plane, origin=center).circle(radius) - return w - - -def multiline_to_brep(comp: DrawMultiline): - from cadquery import Workplane - vs = [vector_to_tuple(x) for x in comp.points] - - bad = [vs[i] == vs[i+1] for i in range(len(vs)-1)] - while any(bad): - first = [idx for idx, b in enumerate(bad) if b][0] - vs = vs[:first] + vs[first+1:] - bad = [vs[i] == vs[i+1] for i in range(len(vs)-1)] - - w = None - if len(vs) > 1: - w = Workplane('XY').polyline(vs) - return w - - -def drawcalls_to_shape(drawcalls): - from cadquery import Assembly - a = Assembly() - for d in drawcalls: - if isinstance(d, DrawCircle): - brep = circle_to_brep(d) - elif isinstance(d, (DrawLine, DrawMultiline)): - brep = multiline_to_brep(d) - else: - print(f"Drawing of {type(d)} not supported; please extend {__file__}!") - print('\n'.join([f'\t{x}' for x in dir(d) if x[0] != '_'])) - brep = None - if brep is not None: - a = a.add(brep) - a = a.toCompound() - return a +# --------------------------------------------------------------------------- +# Component color palette (cycled per component) +# --------------------------------------------------------------------------- + +# RGB tuples in [0, 1]; chosen to be reasonably distinguishable in sequence. +COLOR_PALETTE = [ + (0.90, 0.10, 0.10), # red + (0.10, 0.55, 0.10), # green + (0.15, 0.35, 0.95), # blue + (1.00, 0.70, 0.00), # orange + (0.60, 0.20, 0.80), # purple + (0.00, 0.80, 0.80), # cyan + (1.00, 0.45, 0.70), # pink + (0.45, 0.35, 0.20), # brown + (0.50, 0.50, 0.50), # gray + (0.20, 0.90, 0.40), # lime +] + + +def color_for(index): + ''' cycling palette color for component number `index` ''' + return COLOR_PALETTE[index % len(COLOR_PALETTE)] + + +def _rgb01(c): + return ' '.join(f'{max(0.0, min(1.0, x)):.4g}' for x in c) + + +def _rgb255(c): + return ' '.join(str(int(round(max(0.0, min(1.0, x)) * 255))) for x in c) + + +# --------------------------------------------------------------------------- +# Small geometry containers +# --------------------------------------------------------------------------- + +class Mesh: + ''' triangle soup: verts (N,3) float64, tris (M,3) int ''' + def __init__(self, verts=None, tris=None): + self.verts = np.asarray(verts, dtype=np.float64).reshape((-1, 3)) + self.tris = np.asarray(tris, dtype=np.int64).reshape((-1, 3)) + if self.tris.size and self.tris.max() >= len(self.verts): + raise ValueError("Mesh: triangle index out of range") + + def empty(self): + return len(self.tris) == 0 + + def transformed(self, m4): + return Mesh(transform_points(m4, self.verts), self.tris) + + def merge(self, other): + if other is None or other.empty(): + return self + if self.empty(): + return other + tris = np.vstack((self.tris, other.tris + len(self.verts))) + return Mesh(np.vstack((self.verts, other.verts)), tris) + +class Part: + ''' one named component: an optional surface mesh plus wireframe polylines ''' + def __init__(self, name, color=(0.8, 0.8, 0.8)): + self.name = name + self.color = color # (r, g, b) in [0, 1] + self.mesh = None + self.polylines = [] # list of (N,3) arrays -def component_to_brep(comp): - from cadquery import Matrix - orient = Matrix([[comp.m4[i+j*4] for i in range(4)] for j in range(3)]) - # A more-complex system would check if the draw calls of the component are, e.g., - # - three perpendicular circles sharing center and radius -> a sphere - # - two circles connected by lines -> a cyllinder - # - any number of lines forming closed faces - # - any number of closed faces forming closed shells - return transform_shape(orient, drawcalls_to_shape(comp.drawcalls)) + def add_mesh(self, mesh): + self.mesh = mesh if self.mesh is None else self.mesh.merge(mesh) + def add_polyline(self, pts): + pts = np.asarray(pts, dtype=np.float64).reshape((-1, 3)) + if len(pts) >= 2: + self.polylines.append(pts) -def instrument_to_assembly(inst): - from cadquery import Assembly - a = Assembly() - for component in inst.components: - a= a.add(component_to_brep(component), name = component.name) - return a + def transformed(self, m4): + p = Part(self.name, self.color) + if self.mesh is not None: + p.mesh = self.mesh.transformed(m4) + p.polylines = [transform_points(m4, pl) for pl in self.polylines] + return p + def is_empty(self): + return (self.mesh is None or self.mesh.empty()) and not self.polylines + + +def transform_points(m4, pts): + ''' apply a row-major 4x4 (3x4 usable) matrix to an (N,3) point array ''' + pts = np.asarray(pts, dtype=np.float64).reshape((-1, 3)) + r = np.array(m4, dtype=np.float64).reshape(4, 4)[:3, :3] + t = np.array(m4, dtype=np.float64).reshape(4, 4)[:3, 3] + return pts @ r.T + t + + +# --------------------------------------------------------------------------- +# Tessellation of the mccodelib draw calls (local coordinates) +# --------------------------------------------------------------------------- + +def _basis_from_axis(axis): + ''' orthonormal (w, u, v) with w parallel to axis ''' + w = np.asarray(axis, dtype=np.float64).reshape(3) + n = np.linalg.norm(w) + if n < 1e-12: + w = np.array([0.0, 0.0, 1.0]) + else: + w = w / n + a = np.array([1.0, 0.0, 0.0]) if abs(w[0]) < 0.9 else np.array([0.0, 1.0, 0.0]) + u = np.cross(w, a) + u /= np.linalg.norm(u) + v = np.cross(w, u) + return w, u, v + + +def _circle_pts(center, radius, axis, segments=64, closed=True): + ''' points on a circle of given center, radius and axis (default z) ''' + c = np.asarray(center, dtype=np.float64).reshape(3) + if axis is not None: + _, u, v = _basis_from_axis(axis) + else: + u, v = np.array([1.0, 0.0, 0.0]), np.array([0.0, 1.0, 0.0]) + theta = np.linspace(0.0, 2.0 * np.pi, segments, endpoint=not closed) + pts = c + radius * (np.outer(np.cos(theta), u) + np.outer(np.sin(theta), v)) + if closed: # close the loop + pts = np.vstack((pts, pts[0])) + return pts + + +def _fan(center, ring_pts): + ''' triangulate a (closed, non-repeated) ring into a fan from center ''' + c = np.asarray(center, dtype=np.float64).reshape((1, 3)) + ring = np.asarray(ring_pts, dtype=np.float64).reshape((-1, 3)) + n = len(ring) + verts = np.vstack((c, ring)) + tris = np.array([[0, 1 + i, 1 + (i + 1) % n] for i in range(n)], dtype=np.int64) + return Mesh(verts, tris) + + +def _strip(ring_a, ring_b, flip=False): + ''' triangulate the side wall between two same-length rings ''' + a = np.asarray(ring_a, dtype=np.float64).reshape((-1, 3)) + b = np.asarray(ring_b, dtype=np.float64).reshape((-1, 3)) + n = len(a) + verts = np.vstack((a, b)) + off = n + tris = [] + for i in range(n): + j = (i + 1) % n + if flip: + tris.append([i, j, off + j]) + tris.append([i, off + j, off + i]) + else: + tris.append([i, off + j, j]) + tris.append([i, off + i, off + j]) + return Mesh(verts, np.array(tris, dtype=np.int64)) + + +def _v3(p): + return (p.x, p.y, p.z) + + +def tessellate(d, segments=64): + ''' convert one draw call into (mesh_or_None, [polylines]) in local coords ''' + if isinstance(d, DrawMagnify): + return None, [] + if isinstance(d, DrawPolygon): + return None, [] # awaiting C code upstream; nothing to draw + if isinstance(d, (DrawLine,)): + p1, p2 = _v3(d.point_1), _v3(d.point_2) + return None, [np.array([p1, p2])] + if isinstance(d, DrawDashedLine): + p1, p2 = np.asarray(_v3(d.point_1)), np.asarray(_v3(d.point_2)) + n = max(2, int(np.linalg.norm(p2 - p1) / 0.05) + 1) # 5 cm dashes + t = np.linspace(0, 1, n) + seg = [p1 + (p2 - p1) * ti for ti in t] + lines = [np.array([seg[i], seg[i + 1]]) for i in range(0, len(seg) - 1, 2)] + return None, lines + if isinstance(d, DrawMultiline): + return None, [np.array([_v3(p) for p in d.points])] + if isinstance(d, DrawRectangle): + plane = d.plane.lower() + c = _v3(d.center) + w, h = d.width / 2.0, d.height / 2.0 + if plane == 'xy': + quad = [(c[0]-w, c[1]-h, c[2]), (c[0]+w, c[1]-h, c[2]), + (c[0]+w, c[1]+h, c[2]), (c[0]-w, c[1]+h, c[2])] + elif plane == 'yz': + quad = [(c[0], c[1]-w, c[2]-h), (c[0], c[1]-w, c[2]+h), + (c[0], c[1]+w, c[2]+h), (c[0], c[1]+w, c[2]-h)] + else: # xz + quad = [(c[0]-w, c[1], c[2]-h), (c[0]+w, c[1], c[2]-h), + (c[0]+w, c[1], c[2]+h), (c[0]-w, c[1], c[2]+h)] + quad = np.array(quad) + closed = np.vstack((quad, quad[0])) + mesh = Mesh(quad, [[0, 1, 2], [0, 2, 3]]) + return mesh, [closed] + if isinstance(d, DrawBox): + c = _v3(d.center) + x, y, z = d.xwidth / 2.0, d.yheight / 2.0, d.zlength / 2.0 + if getattr(d, 'thickness', 0) and d.thickness > 0: + t = d.thickness / 2.0 + # hollow box: keep the outline of the two faces (simple wireframe) + lines = [] + for s in (+1, -1): + r = [(c[0]-x, c[1]-y, c[2]+s*z), (c[0]+x, c[1]-y, c[2]+s*z), + (c[0]+x, c[1]+y, c[2]+s*z), (c[0]-x, c[1]+y, c[2]+s*z)] + lines.append(np.vstack((r, r[0]))) + return None, lines + V = np.array([(sx*x+c[0], sy*y+c[1], sz*z+c[2]) + for sz in (-1, 1) for sy in (-1, 1) for sx in (-1, 1)]) + # faces: bottom(z-), top(z+), sides + F = [[0, 1, 3, 2], [4, 6, 7, 5], [0, 4, 5, 1], [2, 3, 7, 6], [0, 2, 6, 4], [1, 5, 7, 3]] + tris, lines = [], [] + for f in F: + tris += [[f[0], f[1], f[2]], [f[0], f[2], f[3]]] + lines.append(np.vstack((V[f], V[f[0]]))) + return Mesh(V, tris), lines + if isinstance(d, DrawCircle): + axis = {'xy': (0, 0, 1), 'xz': (0, 1, 0), 'yz': (1, 0, 0), + 'yx': (0, 0, 1), 'zx': (0, 1, 0), 'zy': (1, 0, 0)}.get(d.plane.lower(), (0, 0, 1)) + return None, [_circle_pts(_v3(d.center), d.radius, axis, segments)] + if isinstance(d, DrawNewCircle): + return None, [_circle_pts(_v3(d.center), d.radius, _v3(d.axis_vector), segments)] + if isinstance(d, DrawSphere): + c = np.asarray(_v3(d.center)) + r = d.radius + nu, nv = segments, max(8, segments // 2) + u = np.linspace(0, 2 * np.pi, nu, endpoint=False) + v = np.linspace(0, np.pi, nv) + uu, vv = np.meshgrid(u, v) + pts = np.column_stack((c[0] + r * np.cos(uu) * np.sin(vv), + c[1] + r * np.sin(uu) * np.sin(vv), + c[2] + r * np.cos(vv))).reshape(-1, 3) + tris = [] + for i in range(nv - 1): + for j in range(nu): + a = i * nu + j + b = i * nu + (j + 1) % nu + cc = (i + 1) * nu + j + dd = (i + 1) * nu + (j + 1) % nu + tris += [[a, b, cc], [b, dd, cc]] + return Mesh(pts, tris), [] + if isinstance(d, DrawCone): + c = np.asarray(_v3(d.center)) + w, u, v = _basis_from_axis(_v3(d.axis_vector)) + th = np.linspace(0, 2*np.pi, segments, endpoint=False) + ring = c + d.radius * (np.outer(np.cos(th), u) + np.outer(np.sin(th), v)) + apex = (c + w * d.height).reshape((1, 3)) + mesh = _strip(ring, np.repeat(apex, len(ring), axis=0), flip=True) + closed = np.vstack((ring, ring[0])) + return mesh, [np.vstack((closed, apex[0], closed[0]))] + if isinstance(d, DrawCylinder): + c = np.asarray(_v3(d.center)) + w, u, v = _basis_from_axis(_v3(d.axis_vector)) + r, h = d.radius, d.height + th = np.linspace(0, 2 * np.pi, segments, endpoint=False) + ring = c + r * (np.outer(np.cos(th), u) + np.outer(np.sin(th), v)) + bot = ring - w * (h / 2.0) + top = ring + w * (h / 2.0) + mesh = _strip(bot, top) + mesh = mesh.merge(_fan(c - w * (h / 2.0), bot)) + mesh = mesh.merge(_fan(c + w * (h / 2.0), top)) + t = getattr(d, 'thickness', 0.0) or 0.0 + lines = [] + if 0 < t < r: # hollow tube: add inner wall lines + inring = c + (r - t) * (np.outer(np.cos(th), u) + np.outer(np.sin(th), v)) + for rr in (bot, top, inring - w * (h / 2.0), inring + w * (h / 2.0)): + lines.append(np.vstack((rr, rr[0]))) + else: + lines = [np.vstack((bot, bot[0])), np.vstack((top, top[0]))] + return mesh, lines + if isinstance(d, DrawDisc): + th = np.linspace(0, 2 * np.pi, segments, endpoint=False) + c = np.asarray(_v3(d.center)) + # build in the plane defined by axis_vector + _, u, v = _basis_from_axis(_v3(d.axis_vector)) + ring = c + d.radius * (np.outer(np.cos(th), u) + np.outer(np.sin(th), v)) + mesh = _fan(c, ring) + return mesh, [np.vstack((ring, ring[0]))] + if isinstance(d, DrawAnnulus): + c = np.asarray(_v3(d.center)) + _, u, v = _basis_from_axis(_v3(d.axis_vector)) + th = np.linspace(0, 2 * np.pi, segments, endpoint=False) + outer = c + d.outer_radius * (np.outer(np.cos(th), u) + np.outer(np.sin(th), v)) + inner = c + d.inner_radius * (np.outer(np.cos(th), u) + np.outer(np.sin(th), v)) + mesh = _strip(outer, inner, flip=True) + return mesh, [np.vstack((outer, outer[0])), np.vstack((inner, inner[0]))] + if isinstance(d, DrawPolyhedron): + return _polyhedron(d.faces_vertices) + print(f"Drawing of {type(d)} not supported; please extend {__file__}!") + return None, [] + + +def _polyhedron(faces_vertices): + ''' best-effort polyhedron: accepts [ [ (x,y,z)... ] per face ] or + [ [(x,y,z)...], ... ] style structures produced by instrparser ''' + try: + faces = faces_vertices + # accept (vertices, faces-as-index-lists) tuple + if isinstance(faces, tuple) and len(faces) == 2: + verts, idxs = faces + verts = np.asarray(verts, dtype=np.float64).reshape((-1, 3)) + faces = [verts[np.asarray(ix).astype(int)] for ix in idxs] + meshes, lines = [], [] + for f in faces: + f = np.asarray(f, dtype=np.float64).reshape((-1, 3)) + if len(f) < 3: + if len(f) == 2: + lines.append(f) + continue + tris = np.array([[0, i, i + 1] for i in range(1, len(f) - 1)], dtype=np.int64) + meshes.append(Mesh(f, tris)) + lines.append(np.vstack((f, f[0]))) + mesh = None + for m in meshes: + mesh = m if mesh is None else mesh.merge(m) + return mesh, lines + except Exception as e: + print(f"Polyhedron interpretation failed ({e}); skipped") + return None, [] + + +# --------------------------------------------------------------------------- +# instrument -> parts +# --------------------------------------------------------------------------- + +def component_to_part(comp, color): + part = Part(comp.name, color) + for d in comp.drawcalls: + mesh, lines = tessellate(d) + if mesh is not None: + part.add_mesh(mesh) + for l in lines: + part.add_polyline(l) + return part.transformed(comp.m4) + + +def instrument_to_parts(inst): + ''' convert components to parts, assigning a cycling color per component ''' + parts = [] + ci = 0 # color index over all components (so colors differ even if some are empty) + for comp in inst.components: + part = component_to_part(comp, color_for(ci)) + ci += 1 + if not part.is_empty(): + parts.append(part) + return parts + + +# --------------------------------------------------------------------------- +# Exporters (all pure Python) +# --------------------------------------------------------------------------- + +def _fmt(x): + return f'{x:.6g}' + + +def export_obj(parts, path): + ''' OBJ with a companion .mtl carrying one material (color) per component ''' + mtl_path = str(path)[:-len('.obj')] + '.mtl' if str(path).endswith('.obj') \ + else str(path) + '.mtl' + + def mtl_name(i, p): + # material names must not contain whitespace or '#' + return 'comp_' + ''.join(ch if ch.isalnum() else '_' for ch in p.name) + f'_{i}' + + with open(mtl_path, 'w') as mf: + mf.write('# McCode mcdisplay export (MTL)\n') + for i, p in enumerate(parts): + r, g, b = p.color + mf.write(f'newmtl {mtl_name(i, p)}\n') + mf.write(f'Kd {r:.4g} {g:.4g} {b:.4g}\nKa 0 0 0\nKs 0.2 0.2 0.2\nNs 10\n\n') + with open(path, 'w') as f: + f.write('# McCode mcdisplay export (OBJ)\n') + f.write(f'mtllib {Path(mtl_path).name}\n') + # OBJ vertex indices are global across the whole file: one running offset + voff = 1 + for i, p in enumerate(parts): + f.write(f'o {p.name}\ng {p.name}\n') + for pl in p.polylines: + f.write(''.join(f'v {_fmt(x)} {_fmt(y)} {_fmt(z)}\n' for x, y, z in pl)) + f.write('l ' + ' '.join(str(voff + i2) for i2 in range(len(pl))) + '\n') + voff += len(pl) + if p.mesh is not None and not p.mesh.empty(): + f.write(f'usemtl {mtl_name(i, p)}\n') + f.write(''.join(f'v {_fmt(x)} {_fmt(y)} {_fmt(z)}\n' for x, y, z in p.mesh.verts)) + for a, b, c in p.mesh.tris: + f.write(f'f {voff+a} {voff+b} {voff+c}\n') + voff += len(p.mesh.verts) + + +def export_stl(parts, path): + tris = [] + for p in parts: + if p.mesh is not None and not p.mesh.empty(): + tris.append(p.mesh.verts[p.mesh.tris]) + if tris: + tris = np.vstack(tris) + else: + tris = np.zeros((0, 3, 3)) + with open(path, 'wb') as f: + f.write(b'mcdisplay STL export'.ljust(80, b'\0')) + f.write(struct.pack(' 0 else n + f.write(struct.pack('<3f', *n)) + for v in t: + f.write(struct.pack('<3f', *v)) + f.write(struct.pack('\n') + f.write('\n') + f.write('\n') + for p in parts: + nm = p.name.replace('"', '') + mat = (f'') + if p.polylines: + pts, idx, off = [], [], 0 + for pl in p.polylines: + for v in pl: + pts.append('%.6g %.6g %.6g' % tuple(v)) + idx.extend([off + i for i in range(len(pl))] + [-1]) + off += len(pl) + f.write(f'{mat}') + f.write(f'' + f'') + f.write(f'\n') + if p.mesh is not None and not p.mesh.empty(): + m = p.mesh + pts = ' '.join('%.6g %.6g %.6g' % tuple(v) for v in m.verts) + idx = ' '.join(map(str, sum(([int(a), int(b), int(c), -1] for a, b, c in m.tris), []))) + f.write(f'{mat}') + f.write(f'' + f'') + f.write(f'\n') + f.write('\n\n') + + +def export_vtk(parts, path): + pts, lines, faces = [], [], [] + for p in parts: + for pl in p.polylines: + base = len(pts) + pts.extend(map(tuple, pl)) + lines.append(len(pl)) + lines.extend(range(base, base + len(pl))) + if p.mesh is not None and not p.mesh.empty(): + base = len(pts) + pts.extend(map(tuple, p.mesh.verts)) + for t in p.mesh.tris: + faces.append(3) + faces.extend(base + int(i) for i in t) + with open(path, 'w') as f: + f.write('# vtk DataFile Version 3.0\nMcCode mcdisplay export\nASCII\nDATASET POLYDATA\n') + f.write(f'POINTS {len(pts)} double\n') + for v in pts: + f.write(f'{_fmt(v[0])} {_fmt(v[1])} {_fmt(v[2])}\n') + nlpoly = sum(len(pp.polylines) for pp in parts) + f.write(f'LINES {nlpoly} {len(lines)}\n') + f.write(' '.join(map(str, lines)) + '\n') + nface = sum(1 for x in faces if x == 3) + f.write(f'POLYS {nface} {len(faces)}\n') + f.write(' '.join(map(str, faces)) + '\n') + # per-point component color scalars (COLORSCALAR names are auto-loaded by many VTK viewers) + f.write(f'POINT_DATA {len(pts)}\n') + f.write('COLORSCALARS component_color 3\n') + for p in parts: + npts = sum(len(pl) for pl in p.polylines) + if p.mesh is not None and not p.mesh.empty(): + npts += len(p.mesh.verts) + for _ in range(npts): + f.write(' '.join(f'{c:.4g}' for c in p.color) + '\n') + + +def _tube_mesh(p0, p1, radius, segments=6): + ''' open cylinder (no caps) between two points, as a triangle mesh; + viewers see the inside via doubleSided materials ''' + p0 = np.asarray(p0, dtype=np.float64).reshape(3) + p1 = np.asarray(p1, dtype=np.float64).reshape(3) + w, u, v = _basis_from_axis(p1 - p0) + th = np.linspace(0, 2 * np.pi, segments, endpoint=False) + ring = radius * (np.outer(np.cos(th), u) + np.outer(np.sin(th), v)) + return _strip(p0 + ring, p1 + ring) + + +def _gltf_buffers(parts): + ''' build accessor+buffers for all parts; returns glTF dict + binary blob ''' + bin_chunks = [] + accessors = [] + buffviews = [] + materials = [] + + def add(arr, dtype, target): + a = np.ascontiguousarray(arr) + off = sum(len(b) for b in bin_chunks) + # pad to 4 bytes + pad = (4 - (off % 4)) % 4 + if pad: + bin_chunks.append(b'\0' * pad) + off += pad + data = a.astype(dtype).tobytes() + bin_chunks.append(data) + buffviews.append({'buffer': 0, 'byteOffset': off, 'byteLength': len(data), + 'target': target}) + comp = {np.dtype('float32'): 5126, np.dtype('uint32'): 5125}[np.dtype(dtype)] + if a.ndim == 2 and a.shape[1] == 3: + atype, amin, amax = 'VEC3', a.min(axis=0), a.max(axis=0) + else: + atype, amin, amax = 'SCALAR', np.array([a.min()]), np.array([a.max()]) + acc = {'bufferView': len(buffviews) - 1, 'componentType': comp, + 'count': int(len(a)), 'type': atype, + 'min': [float(x) for x in amin], 'max': [float(x) for x in amax]} + accessors.append(acc) + return len(accessors) - 1 + + def add_material(color): + materials.append({'pbrMetallicRoughness': { + 'baseColorFactor': [float(c) for c in color] + [1.0], + 'metallicFactor': 0.1, 'roughnessFactor': 0.7}, + 'doubleSided': True}) # wireframe tubes have no caps + return len(materials) - 1 + + nodes, meshes = [], [] + # Wireframe lines are exported as thin tubes (triangle meshes) instead of + # glTF line primitives (mode 1): Blender's glTF-IO importer crashes on line + # primitives with a numpy concatenation error, and web viewers (three.js, + # f3d, castle-model-viewer) support them inconsistently. Tubes load as + # plain triangles everywhere. Radius scales with the scene size. + all_pts = [] + for p in parts: + if p.mesh is not None and not p.mesh.empty(): + all_pts.append(p.mesh.verts) + all_pts.extend(p.polylines) + diag = float(np.linalg.norm(np.vstack(all_pts).max(axis=0) - np.vstack(all_pts).min(axis=0))) \ + if all_pts else 1.0 + tube_r = max(diag * 0.00002, 1e-6) + for i, p in enumerate(parts): + mat_idx = add_material(p.color) + mesh = p.mesh + if p.polylines: + tubes = None + for pl in p.polylines: + # merge consecutive collinear segments into one tube + runs = [[pl[0], pl[1]]] + for j in range(1, len(pl) - 1): + a, b = pl[j], pl[j + 1] + prev = runs[-1] + d_prev = prev[1] - prev[0] + d_cur = b - a + n_p, n_c = np.linalg.norm(d_prev), np.linalg.norm(d_cur) + if n_p > 0 and n_c > 0 and abs(np.dot(d_prev / n_p, d_cur / n_c) - 1.0) < 1e-9: + prev[1] = b # collinear continuation: extend the previous run + else: + runs.append([a, b]) + for a, b in runs: + t = _tube_mesh(a, b, tube_r) + tubes = t if tubes is None else tubes.merge(t) + if tubes is not None: + mesh = tubes if mesh is None else mesh.merge(tubes) + if mesh is not None and not mesh.empty(): + pv = add(mesh.verts, 'float32', 34962) + pi = add(mesh.tris.reshape(-1), 'uint32', 34963) + meshes.append({'primitives': [{'attributes': {'POSITION': pv}, 'indices': pi, + 'mode': 4, 'material': mat_idx}]}) + nodes.append({'mesh': len(meshes) - 1, 'name': p.name}) + blob = b''.join(bin_chunks) + gltf = {'asset': {'version': '2.0', 'generator': 'mcdisplay (dependency-free)'}, + 'scene': 0, + 'scenes': [{'nodes': list(range(len(nodes)))}], + 'nodes': nodes, 'meshes': meshes, 'materials': materials, + 'buffers': [{'byteLength': len(blob)}], + 'bufferViews': buffviews, 'accessors': accessors} + return gltf, blob + + +def export_gltf(parts, path): + gltf, blob = _gltf_buffers(parts) + gltf['buffers'][0]['uri'] = Path(path).name + '.bin' + with open(path, 'w') as f: + import json + json.dump(gltf, f, separators=(',', ':')) + with open(str(path) + '.bin', 'wb') as f: + f.write(blob) + + +def export_glb(parts, path): + import json + gltf, blob = _gltf_buffers(parts) + gltf['buffers'][0]['bufferView'] = len(gltf['bufferViews']) + gltf['bufferViews'].append({'buffer': 0, 'byteOffset': 0, 'byteLength': len(blob)}) + js = json.dumps(gltf, separators=(',', ':')).encode() + while len(js) % 4: + js += b' ' + while len(blob) % 4: + blob += b'\0' + with open(path, 'wb') as f: + f.write(b'glTF' + struct.pack(' faceted BREP: OPEN_SHELL of ADVANCED_FACEs, + each triangle a PLANE bounded by a POLY_LOOP of 3 VERTEX_POINTs, + wrapped in a SHELL_BASED_SURFACE_MODEL + (MANIFOLD_SURFACE_SHAPE_REPRESENTATION) + - wireframe-only components -> GEOMETRICALLY_BOUNDED_WIREFRAME_SHAPE_ + REPRESENTATION of a GEOMETRIC_CURVE_SET of POLYLINEs + Component colors are attached as STYLED_ITEMs (COLOUR_RGB) where the + reader supports it (e.g. FreeCAD). Coordinates are in metres. ''' + entities = [] + nid = 0 + + def E(text): + nonlocal nid + nid += 1 + entities.append(f'#{nid}={text};') + return nid + + def pt(p): + return E(f"CARTESIAN_POINT('',({_fmt(p[0])},{_fmt(p[1])},{_fmt(p[2])}))") + + # --- global contexts and units --- + appctx = E("APPLICATION_CONTEXT('core data for automotive mechanical design processes')") + E(f"APPLICATION_PROTOCOL_DEFINITION('international standard','automotive_design',2010,{appctx})") + lu = E('SI_UNIT($,.METRE.)') + pau = E('SI_UNIT($,.RADIAN.)') + sau = E('SI_UNIT($,.STERADIAN.)') + unc = E(f"UNCERTAINTY_MEASURE_WITH_UNIT(LENGTH_MEASURE(1.E-07),{lu}," + f"'distance_accuracy_value','')") + # one complex instance for the representation context + nid += 1 + ctx = nid + entities.append(f"#{ctx}=(GEOMETRIC_REPRESENTATION_CONTEXT(3)" + f"GLOBAL_UNCERTAINTY_ASSIGNED_CONTEXT(({unc}))" + f"GLOBAL_UNIT_ASSIGNED_CONTEXT(({lu},{pau},{sau}))" + f"REPRESENTATION_CONTEXT('',''));" ) + + # --- one product per component --- + for p in parts: + nm = p.name.replace("'", "''").replace(';', '').replace('#', '') + pctx = E(f"PRODUCT_CONTEXT('',{appctx},'mechanical')") + prod = E(f"PRODUCT('{nm}','{nm}','',({pctx}))") + form = E(f"PRODUCT_DEFINITION_FORMATION_WITH_SPECIFIED_SOURCE('','',{prod},.NOT_KNOWN.)") + pdctx = E(f"PRODUCT_DEFINITION_CONTEXT('part definition',{appctx},'design')") + pdef = E(f"PRODUCT_DEFINITION('design','',{form},{pdctx})") + pds = E(f"PRODUCT_DEFINITION_SHAPE('','',{pdef})") + + styled = [] # styled item ids for this component + rep = None + if p.mesh is not None and not p.mesh.empty(): + faces = [] + for tri in p.mesh.verts[p.mesh.tris]: + a, b, c = tri + n = np.cross(b - a, c - a) + ln = np.linalg.norm(n) + if ln < 1e-12: + continue # skip degenerate triangles + n = n / ln + pa, pb, pc = pt(a), pt(b), pt(c) + va = E(f'VERTEX_POINT(\'\',{pa})') + vb = E(f'VERTEX_POINT(\'\',{pb})') + vc = E(f'VERTEX_POINT(\'\',{pc})') + loop = E(f'POLY_LOOP(\'\',({va},{vb},{vc}),.T.)') + fb = E(f'FACE_OUTER_BOUND(\'\',{loop},.T.)') + dn = E(f"DIRECTION('',({_fmt(n[0])},{_fmt(n[1])},{_fmt(n[2])}))") + ref = np.array([1.0, 0.0, 0.0]) if abs(n[0]) < 0.9 else np.array([0.0, 1.0, 0.0]) + dr = np.cross(n, ref) + dr /= np.linalg.norm(dr) + dd = E(f"DIRECTION('',({_fmt(dr[0])},{_fmt(dr[1])},{_fmt(dr[2])}))") + ax = E(f'AXIS2_PLACEMENT_3D(\'\',{pa},{dn},{dd})') + pl = E(f'PLANE(\'\',{ax})') + faces.append(E(f'ADVANCED_FACE(\'\',({fb}),{pl},.T.)')) + if faces: + shell = E(f"OPEN_SHELL('',({','.join(map(str, faces))}))") + sbsm = E(f'SHELL_BASED_SURFACE_MODEL(\'\',({shell}))') + rep = E(f'MANIFOLD_SURFACE_SHAPE_REPRESENTATION(\'\',({sbsm}),{ctx})') + # attach the component color to every face + col = E(f"COLOUR_RGB('',{p.color[0]:.4g},{p.color[1]:.4g},{p.color[2]:.4g})") + fasc = E(f'FILL_AREA_STYLE_COLOUR(\'\',{col})') + fas = E(f'FILL_AREA_STYLE(\'\',({fasc}))') + ssf = E(f'SURFACE_STYLE_FILL_AREA({fas})') + sss = E(f"SURFACE_SIDE_STYLE('',({ssf}))") + ssu = E(f'SURFACE_STYLE_USAGE(.BOTH.,{sss})') + psa = E(f'PRESENTATION_STYLE_ASSIGNMENT(({ssu}))') + for f_id in faces: + styled.append(E(f"STYLED_ITEM('',({psa}),{f_id})")) + if rep is None and p.polylines: + # wireframe-only component + pols = [] + for pl in p.polylines: + ids = ','.join(str(pt(v)) for v in pl) + pols.append(E(f"POLYLINE('',({ids}))")) + gcs = E(f"GEOMETRIC_CURVE_SET('',({','.join(map(str, pols))}))") + rep = E(f'GEOMETRICALLY_BOUNDED_WIREFRAME_SHAPE_REPRESENTATION(\'\',({gcs}),{ctx})') + if rep is None: + continue # nothing representable in STEP + E(f'SHAPE_DEFINITION_REPRESENTATION({pds},{rep})') + if styled: + E(f"MECHANICAL_DESIGN_GEOMETRIC_PRESENTATION_REPRESENTATION(''," + f"({','.join(map(str, styled))}),{ctx})") + + with open(path, 'w') as f: + f.write('ISO-10303-21;\n') + f.write('HEADER;\n') + f.write("FILE_DESCRIPTION(('MCDISPLAY geometry export'),'2;1');\n") + f.write(f"FILE_NAME('{Path(path).name}','{__date__}',('mcdisplay')," + f"('McCode'),'mcdisplay dependency-free','');\n") + f.write("FILE_SCHEMA(('AUTOMOTIVE_DESIGN { 1 0 10303 214 1 1 1 1 }'));\n") + f.write('ENDSEC;\nDATA;\n') + f.write('\n'.join(entities) + '\n') + f.write('ENDSEC;\nEND-ISO-10303-21;\n') + + +def export_vtkjs(parts, path): + ''' Pure-Python VTKJS export -- no CadQuery/VTK needed. + + Writes a .vtkjs (zip) in the VTK JSON scene format used by the vtk.js + HttpSceneLoader / ParaView Glance viewers: + index.json scene description (one entry per component) + /index.json vtkPolyData dataset for component n, with binary + arrays base64-encoded (float32 points, uint32 + connectivity), plus a per-component RGB point + data array for the cycling colors. ''' + import base64 + import json + import zipfile + + def b64(arr, dtype): + return base64.b64encode(np.ascontiguousarray(arr, dtype=dtype).tobytes()).decode() + + def ranges(pts): + out = [] + for c, nm in ((0, 'X'), (1, 'Y'), (2, 'Z')): + col = pts[:, c] + out.append({'min': float(col.min()), 'max': float(col.max()), + 'component': c, 'name': nm}) + return out + + scene = [] + datasets = [] + for i, p in enumerate(parts): + pts, lines, polys = [], [], [] + for pl in p.polylines: + base = len(pts) + pts.extend(map(tuple, pl)) + lines.append(len(pl)) + lines.extend(range(base, base + len(pl))) + if p.mesh is not None and not p.mesh.empty(): + m = p.mesh + base = len(pts) + pts.extend(map(tuple, m.verts)) + for t in m.tris: + polys += [3] + [base + int(x) for x in t] + if not pts: + continue + pts = np.asarray(pts, dtype=np.float64).reshape((-1, 3)) + ds = { + 'metadata': {'name': p.name}, + 'vtkClass': 'vtkPolyData', + 'points': {'vtkClass': 'vtkPoints', 'name': '_points', + 'numberOfComponents': 3, 'size': int(len(pts) * 3), + 'dataType': 'Float32Array', 'buffer': b64(pts, ' Date: Sun, 18 Oct 2026 22:40:40 +0900 Subject: [PATCH 2/8] mcdisplay-cad: better implementation. faster. --- tools/Python/mcdisplay/cad/mcdisplay.py | 51 +++++++++++++++---------- 1 file changed, 31 insertions(+), 20 deletions(-) diff --git a/tools/Python/mcdisplay/cad/mcdisplay.py b/tools/Python/mcdisplay/cad/mcdisplay.py index 55f8a9b4ba..b92a4008f8 100755 --- a/tools/Python/mcdisplay/cad/mcdisplay.py +++ b/tools/Python/mcdisplay/cad/mcdisplay.py @@ -2,12 +2,11 @@ # -*- coding: utf-8 -*- ''' Converts MCDISPLAY geometry information, read from an instrument binary, -to common CAD / 3D formats WITHOUT CadQuery / OpenCASCADE. +to common CAD / 3D formats. Supported output formats: glTF (.gltf), GLB, OBJ (+MTL), STL, PLY, VRML97 (.wrl), X3D (also as legacy 'xml'), VTK legacy (.vtk), STEP (AP214 faceted BREP, pure Python), and VTKJS (.vtkjs zip of the vtk.js JSON scene format). -No CadQuery/OpenCASCADE is required for any format. Each instrument component is assigned a color from a cycling palette, so adjacent components are visually distinguishable in the output. @@ -25,7 +24,7 @@ - polyhedron -> triangulated faces (fan) Every draw call is placed in global coordinates via the component m4 matrix. ''' -__author__ = "Gregory Tucker (CadQuery version); dependency-free version" +__author__ = "Gregory Tucker (CadQuery version); Emmanuel Farhi (drop CadQuery depends)" __date__ = "2026-10-17" import sys @@ -369,10 +368,10 @@ def _polyhedron(faces_vertices): # instrument -> parts # --------------------------------------------------------------------------- -def component_to_part(comp, color): +def component_to_part(comp, color, segments=64): part = Part(comp.name, color) for d in comp.drawcalls: - mesh, lines = tessellate(d) + mesh, lines = tessellate(d, segments) if mesh is not None: part.add_mesh(mesh) for l in lines: @@ -380,12 +379,12 @@ def component_to_part(comp, color): return part.transformed(comp.m4) -def instrument_to_parts(inst): +def instrument_to_parts(inst, segments=64): ''' convert components to parts, assigning a cycling color per component ''' parts = [] ci = 0 # color index over all components (so colors differ even if some are empty) for comp in inst.components: - part = component_to_part(comp, color_for(ci)) + part = component_to_part(comp, color_for(ci), segments) ci += 1 if not part.is_empty(): parts.append(part) @@ -595,8 +594,8 @@ def export_vtk(parts, path): f.write(' '.join(f'{c:.4g}' for c in p.color) + '\n') -def _tube_mesh(p0, p1, radius, segments=6): - ''' open cylinder (no caps) between two points, as a triangle mesh; +def _tube_mesh(p0, p1, radius, segments=4): + ''' open 4-sided cylinder (no caps) between two points, as a triangle mesh; viewers see the inside via doubleSided materials ''' p0 = np.asarray(p0, dtype=np.float64).reshape(3) p1 = np.asarray(p1, dtype=np.float64).reshape(3) @@ -656,9 +655,12 @@ def add_material(color): all_pts.extend(p.polylines) diag = float(np.linalg.norm(np.vstack(all_pts).max(axis=0) - np.vstack(all_pts).min(axis=0))) \ if all_pts else 1.0 - tube_r = max(diag * 0.00002, 1e-6) - for i, p in enumerate(parts): - mat_idx = add_material(p.color) + tube_r = max(diag * 2e-5, 1e-6) + # Merge all components sharing the same color into ONE mesh: web viewers + # (castle-model-viewer, three.js) spend memory and time per node/mesh/draw + # call, so a few large meshes load much faster than hundreds of small ones. + by_color = {} # color tuple -> [merged_mesh, [component names]] + for p in parts: mesh = p.mesh if p.polylines: tubes = None @@ -681,11 +683,18 @@ def add_material(color): if tubes is not None: mesh = tubes if mesh is None else mesh.merge(tubes) if mesh is not None and not mesh.empty(): - pv = add(mesh.verts, 'float32', 34962) - pi = add(mesh.tris.reshape(-1), 'uint32', 34963) - meshes.append({'primitives': [{'attributes': {'POSITION': pv}, 'indices': pi, - 'mode': 4, 'material': mat_idx}]}) - nodes.append({'mesh': len(meshes) - 1, 'name': p.name}) + entry = by_color.setdefault(tuple(p.color), [None, []]) + entry[0] = mesh if entry[0] is None else entry[0].merge(mesh) + entry[1].append(p.name) + nodes, meshes = [], [] + for gi, (color, (mesh, names)) in enumerate(by_color.items()): + mat_idx = add_material(color) + pv = add(mesh.verts, 'float32', 34962) + pi = add(mesh.tris.reshape(-1), 'uint32', 34963) + meshes.append({'primitives': [{'attributes': {'POSITION': pv}, 'indices': pi, + 'mode': 4, 'material': mat_idx}]}) + nm = names[0] if len(names) == 1 else f'color_group_{gi} ({len(names)} components)' + nodes.append({'mesh': len(meshes) - 1, 'name': nm}) blob = b''.join(bin_chunks) gltf = {'asset': {'version': '2.0', 'generator': 'mcdisplay (dependency-free)'}, 'scene': 0, @@ -946,6 +955,7 @@ def main(instr=None, dirname=None, **kwds): basicConfig(level=INFO) fmt = kwds.get('format', 'stl') + segments = int(kwds.pop('segments', 64) or 64) # tessellation resolution # output directory if dirname is None: @@ -962,7 +972,7 @@ def main(instr=None, dirname=None, **kwds): root.mkdir(parents=True, exist_ok=True) # Convert the instrument into tessellated parts, then save - parts = instrument_to_parts(instrument) + parts = instrument_to_parts(instrument, segments) output = str(root.joinpath(f"{Path(instr).stem}{EXPORTERS[fmt][1]}")) EXPORTERS[fmt][0](parts, output) print(f"Wrote {output} ({len(parts)} components)") @@ -977,7 +987,7 @@ def output_format(astring): if f in ('step', 'stl', 'xml', 'vrml', 'gltf', 'vtkjs', 'obj', 'ply', 'glb', 'x3d', 'vtk'): return f print(f"Unknown file format {astring}. Choose between STEP, STL, VRML, GLTF, GLB, " - f"OBJ, PLY, X3D, VTK, VTKJS, and XML (X3D)") + f"OBJ, PLY, XML/X3D, VTK, VTKJS") return 'stl' @@ -988,7 +998,8 @@ def output_format(astring): parser.add_argument('--dirname', help='output directory name override') parser.add_argument('-n', '--ncount', dest='n', type=float, default=0, help='Number of particles to simulate') - parser.add_argument('-f', '--format', dest='format', type=output_format, default='stl', help='Output CAD file format (STEP, STL, VRML, GLTF, GLB, OBJ, PLY, X3D, VTK, and XML/X3D)') + parser.add_argument('-s', '--segments', dest='segments', type=int, default=16, help='Tessellation resolution (circle segments); lower = smaller/lighter files, e.g. 32 or 16. Default: 16') + parser.add_argument('-f', '--format', dest='format', type=output_format, default='glb', help='Output CAD file format (STEP, STL, VRML, GLTF, GLB, OBJ, PLY, XML/X3D, VTK and VTKJS). Default: GLB') args, unknown = parser.parse_known_args() # if --inspect --first or --last are given after instr, the remaining args become "unknown", From a712abcb109e12cc9cf15edec8859d5692a58431 Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Sun, 18 Oct 2026 17:09:30 +0200 Subject: [PATCH 3/8] mcdisplay-cad: drop unused cadquery test-import from wrapper Co-Authored-By: Claude Opus 5.5 --- tools/Python/mcdisplay/cad/mcdisplay.in | 7 ------- 1 file changed, 7 deletions(-) diff --git a/tools/Python/mcdisplay/cad/mcdisplay.in b/tools/Python/mcdisplay/cad/mcdisplay.in index 01d2e885ed..cbe0f6eecb 100755 --- a/tools/Python/mcdisplay/cad/mcdisplay.in +++ b/tools/Python/mcdisplay/cad/mcdisplay.in @@ -17,13 +17,6 @@ canrun() { if ! [ -x "${UTILDIR}/${TOOL}.py" ]; then exit 127; fi - - modules="cadquery" - cmd="" - for name in ${modules}; do - cmd="${cmd}import ${name}; " - done - python3 -c "${cmd}" } if ( canrun ); then From 79a21b2e751a960a5d87bf0e684570e6db53ac22 Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Sun, 18 Oct 2026 18:04:15 +0200 Subject: [PATCH 4/8] mcdisplay-cad: use --trace=2, parse polyhedron JSON, valid STEP output, add --tubes - run mcrun with --trace=2 so solid primitives (box/cylinder/cone/sphere/polyhedron) are emitted - parse the MCDISPLAY polyhedron JSON (vertices + faces) instead of failing on it - STEP: write #-references, CARTESIAN_POINT POLY_LOOPs in FACE_SURFACEs, complex SI units, REALs with decimal point, 7-field FILE_NAME (now reads cleanly in OpenCASCADE/FreeCAD) - --tubes: draw lines as cylinders with corner spheres, as in the mcdisplay-plopp 3D view Co-Authored-By: Claude Opus 5.5 --- tools/Python/mcdisplay/cad/mcdisplay.py | 128 ++++++++++++++++-------- 1 file changed, 87 insertions(+), 41 deletions(-) diff --git a/tools/Python/mcdisplay/cad/mcdisplay.py b/tools/Python/mcdisplay/cad/mcdisplay.py index b92a4008f8..c44f5253ea 100755 --- a/tools/Python/mcdisplay/cad/mcdisplay.py +++ b/tools/Python/mcdisplay/cad/mcdisplay.py @@ -28,6 +28,7 @@ __date__ = "2026-10-17" import sys +import json import math import struct from pathlib import Path @@ -335,29 +336,19 @@ def tessellate(d, segments=64): return None, [] -def _polyhedron(faces_vertices): - ''' best-effort polyhedron: accepts [ [ (x,y,z)... ] per face ] or - [ [(x,y,z)...], ... ] style structures produced by instrparser ''' +def _polyhedron(json_str): + ''' polyhedron from the --trace=2 JSON: {"vertices": [[x,y,z],...], "faces": [{"face": [i,j,k,...]},...]} ''' try: - faces = faces_vertices - # accept (vertices, faces-as-index-lists) tuple - if isinstance(faces, tuple) and len(faces) == 2: - verts, idxs = faces - verts = np.asarray(verts, dtype=np.float64).reshape((-1, 3)) - faces = [verts[np.asarray(ix).astype(int)] for ix in idxs] - meshes, lines = [], [] - for f in faces: - f = np.asarray(f, dtype=np.float64).reshape((-1, 3)) + data = json.loads(json_str) + V = np.asarray(data['vertices'], dtype=np.float64).reshape((-1, 3)) + mesh, lines = None, [] + for f in data['faces']: + f = V[np.asarray(f['face'], dtype=np.int64)] if len(f) < 3: - if len(f) == 2: - lines.append(f) continue - tris = np.array([[0, i, i + 1] for i in range(1, len(f) - 1)], dtype=np.int64) - meshes.append(Mesh(f, tris)) - lines.append(np.vstack((f, f[0]))) - mesh = None - for m in meshes: + m = Mesh(f, [[0, i, i + 1] for i in range(1, len(f) - 1)]) mesh = m if mesh is None else mesh.merge(m) + lines.append(np.vstack((f, f[0]))) return mesh, lines except Exception as e: print(f"Polyhedron interpretation failed ({e}); skipped") @@ -605,6 +596,55 @@ def _tube_mesh(p0, p1, radius, segments=4): return _strip(p0 + ring, p1 + ring) +_G = (1 + 5 ** 0.5) / 2 # icosahedron: the sphere at each polyline corner +ICO_V = np.array([(-1, _G, 0), (1, _G, 0), (-1, -_G, 0), (1, -_G, 0), (0, -1, _G), (0, 1, _G), + (0, -1, -_G), (0, 1, -_G), (_G, 0, -1), (_G, 0, 1), (-_G, 0, -1), (-_G, 0, 1)]) / np.hypot(1, _G) +ICO_F = np.array([(0, 11, 5), (0, 5, 1), (0, 1, 7), (0, 7, 10), (0, 10, 11), (1, 5, 9), (5, 11, 4), + (11, 10, 2), (10, 7, 6), (7, 1, 8), (3, 9, 4), (3, 4, 2), (3, 2, 6), (3, 6, 8), + (3, 8, 9), (4, 9, 5), (2, 4, 11), (6, 2, 10), (8, 6, 7), (9, 8, 1)]) + + +def _tubes(xyz, r, k=6): + ''' NaN-separated polylines -> Mesh: a k-sided open cylinder per segment, a sphere per corner + (same construction as the mcdisplay-plopp 3D view) ''' + a, b = xyz[:-1], xyz[1:] + keep = np.isfinite(a).all(1) & np.isfinite(b).all(1) & (np.linalg.norm(b - a, axis=1) > 0) + a, b = a[keep], b[keep] + d = (b - a) / np.linalg.norm(b - a, axis=1, keepdims=True) + u = np.cross(d, np.where(abs(d[:, :1]) < 0.9, [1., 0, 0], [0., 1, 0])) # helper axis not parallel to d + u /= np.linalg.norm(u, axis=1, keepdims=True) + v = np.cross(d, u) + th = 2 * np.pi * np.arange(k) / k + ring = r * (np.cos(th)[None, :, None] * u[:, None] + np.sin(th)[None, :, None] * v[:, None]) # (n, k, 3) + cyl = np.concatenate([a[:, None] + ring, b[:, None] + ring], axis=1).reshape(-1, 3) + j, j1 = np.arange(k), (np.arange(k) + 1) % k + side = np.concatenate([np.stack([j, j1, k + j], 1), np.stack([j1, k + j1, k + j], 1)]) + cyl_f = (side[None] + 2 * k * np.arange(len(a))[:, None, None]).reshape(-1, 3) + + d1, d2 = xyz[1:-1] - xyz[:-2], xyz[2:] - xyz[1:-1] # joints only where the polyline turns > ~20 deg; + with np.errstate(invalid='ignore', divide='ignore'): # NaN breaks, ends and smooth arcs get none + cos = (d1 * d2).sum(1) / (np.linalg.norm(d1, axis=1) * np.linalg.norm(d2, axis=1)) + pts = xyz[1:-1][cos < 0.94] + sph = (pts[:, None] + r * ICO_V[None]).reshape(-1, 3) + sph_f = (ICO_F[None] + len(cyl) + 12 * np.arange(len(pts))[:, None, None]).reshape(-1, 3) + return Mesh(np.concatenate([cyl, sph]), np.concatenate([cyl_f, sph_f])) + + +def lines_to_tubes(parts, thickness=1e-4): + ''' replace every part's polylines by solid tubes; radius = thickness x scene bounding-box diagonal ''' + # ponytail: fixed thickness as in mcdisplay-plopp; expose it on the CLI if 1e-4 doesn't suit + pts = [p.mesh.verts for p in parts if p.mesh is not None] + [pl for p in parts for pl in p.polylines] + if not pts: + return + pts = np.vstack(pts) + r = thickness * np.linalg.norm(pts.max(0) - pts.min(0)) + nan = np.full((1, 3), np.nan) + for p in parts: + if p.polylines: + p.add_mesh(_tubes(np.vstack([np.vstack((pl, nan)) for pl in p.polylines]), r)) + p.polylines = [] + + def _gltf_buffers(parts): ''' build accessor+buffers for all parts; returns glTF dict + binary blob ''' bin_chunks = [] @@ -709,14 +749,12 @@ def export_gltf(parts, path): gltf, blob = _gltf_buffers(parts) gltf['buffers'][0]['uri'] = Path(path).name + '.bin' with open(path, 'w') as f: - import json json.dump(gltf, f, separators=(',', ':')) with open(str(path) + '.bin', 'wb') as f: f.write(blob) def export_glb(parts, path): - import json gltf, blob = _gltf_buffers(parts) gltf['buffers'][0]['bufferView'] = len(gltf['bufferViews']) gltf['bufferViews'].append({'buffer': 0, 'byteOffset': 0, 'byteLength': len(blob)}) @@ -735,8 +773,8 @@ def export_step(parts, path): ''' Pure-Python STEP (AP214-ish) writer -- no CadQuery/OpenCASCADE needed. Each component becomes one PRODUCT: - - components with a mesh -> faceted BREP: OPEN_SHELL of ADVANCED_FACEs, - each triangle a PLANE bounded by a POLY_LOOP of 3 VERTEX_POINTs, + - components with a mesh -> faceted BREP: OPEN_SHELL of FACE_SURFACEs, + each triangle a PLANE bounded by a POLY_LOOP of 3 CARTESIAN_POINTs, wrapped in a SHELL_BASED_SURFACE_MODEL (MANIFOLD_SURFACE_SHAPE_REPRESENTATION) - wireframe-only components -> GEOMETRICALLY_BOUNDED_WIREFRAME_SHAPE_ @@ -750,23 +788,31 @@ def E(text): nonlocal nid nid += 1 entities.append(f'#{nid}={text};') - return nid + return f'#{nid}' + + def r(x): + # STEP REAL needs a decimal point: 0 -> 0., 1E-05 -> 1.E-05 + m, e, x = f'{x:.6G}'.partition('E') + return (m if '.' in m else m + '.') + e + x + + def r3(p): + return f'({r(p[0])},{r(p[1])},{r(p[2])})' def pt(p): - return E(f"CARTESIAN_POINT('',({_fmt(p[0])},{_fmt(p[1])},{_fmt(p[2])}))") + return E(f"CARTESIAN_POINT('',{r3(p)})") # --- global contexts and units --- appctx = E("APPLICATION_CONTEXT('core data for automotive mechanical design processes')") E(f"APPLICATION_PROTOCOL_DEFINITION('international standard','automotive_design',2010,{appctx})") - lu = E('SI_UNIT($,.METRE.)') - pau = E('SI_UNIT($,.RADIAN.)') - sau = E('SI_UNIT($,.STERADIAN.)') + lu = E('(LENGTH_UNIT()NAMED_UNIT(*)SI_UNIT($,.METRE.))') + pau = E('(NAMED_UNIT(*)PLANE_ANGLE_UNIT()SI_UNIT($,.RADIAN.))') + sau = E('(NAMED_UNIT(*)SI_UNIT($,.STERADIAN.)SOLID_ANGLE_UNIT())') unc = E(f"UNCERTAINTY_MEASURE_WITH_UNIT(LENGTH_MEASURE(1.E-07),{lu}," f"'distance_accuracy_value','')") # one complex instance for the representation context nid += 1 - ctx = nid - entities.append(f"#{ctx}=(GEOMETRIC_REPRESENTATION_CONTEXT(3)" + ctx = f'#{nid}' + entities.append(f"{ctx}=(GEOMETRIC_REPRESENTATION_CONTEXT(3)" f"GLOBAL_UNCERTAINTY_ASSIGNED_CONTEXT(({unc}))" f"GLOBAL_UNIT_ASSIGNED_CONTEXT(({lu},{pau},{sau}))" f"REPRESENTATION_CONTEXT('',''));" ) @@ -793,25 +839,22 @@ def pt(p): continue # skip degenerate triangles n = n / ln pa, pb, pc = pt(a), pt(b), pt(c) - va = E(f'VERTEX_POINT(\'\',{pa})') - vb = E(f'VERTEX_POINT(\'\',{pb})') - vc = E(f'VERTEX_POINT(\'\',{pc})') - loop = E(f'POLY_LOOP(\'\',({va},{vb},{vc}),.T.)') + loop = E(f"POLY_LOOP('',({pa},{pb},{pc}))") fb = E(f'FACE_OUTER_BOUND(\'\',{loop},.T.)') - dn = E(f"DIRECTION('',({_fmt(n[0])},{_fmt(n[1])},{_fmt(n[2])}))") + dn = E(f"DIRECTION('',{r3(n)})") ref = np.array([1.0, 0.0, 0.0]) if abs(n[0]) < 0.9 else np.array([0.0, 1.0, 0.0]) dr = np.cross(n, ref) dr /= np.linalg.norm(dr) - dd = E(f"DIRECTION('',({_fmt(dr[0])},{_fmt(dr[1])},{_fmt(dr[2])}))") + dd = E(f"DIRECTION('',{r3(dr)})") ax = E(f'AXIS2_PLACEMENT_3D(\'\',{pa},{dn},{dd})') pl = E(f'PLANE(\'\',{ax})') - faces.append(E(f'ADVANCED_FACE(\'\',({fb}),{pl},.T.)')) + faces.append(E(f"FACE_SURFACE('',({fb}),{pl},.T.)")) if faces: shell = E(f"OPEN_SHELL('',({','.join(map(str, faces))}))") sbsm = E(f'SHELL_BASED_SURFACE_MODEL(\'\',({shell}))') rep = E(f'MANIFOLD_SURFACE_SHAPE_REPRESENTATION(\'\',({sbsm}),{ctx})') # attach the component color to every face - col = E(f"COLOUR_RGB('',{p.color[0]:.4g},{p.color[1]:.4g},{p.color[2]:.4g})") + col = E(f"COLOUR_RGB('',{','.join(r(x) for x in p.color)})") fasc = E(f'FILL_AREA_STYLE_COLOUR(\'\',{col})') fas = E(f'FILL_AREA_STYLE(\'\',({fasc}))') ssf = E(f'SURFACE_STYLE_FILL_AREA({fas})') @@ -840,7 +883,7 @@ def pt(p): f.write('HEADER;\n') f.write("FILE_DESCRIPTION(('MCDISPLAY geometry export'),'2;1');\n") f.write(f"FILE_NAME('{Path(path).name}','{__date__}',('mcdisplay')," - f"('McCode'),'mcdisplay dependency-free','');\n") + f"('McCode'),'mcdisplay dependency-free','McCode','');\n") f.write("FILE_SCHEMA(('AUTOMOTIVE_DESIGN { 1 0 10303 214 1 1 1 1 }'));\n") f.write('ENDSEC;\nDATA;\n') f.write('\n'.join(entities) + '\n') @@ -858,7 +901,6 @@ def export_vtkjs(parts, path): connectivity), plus a per-component RGB point data array for the cycling colors. ''' import base64 - import json import zipfile def b64(arr, dtype): @@ -956,6 +998,7 @@ def main(instr=None, dirname=None, **kwds): fmt = kwds.get('format', 'stl') segments = int(kwds.pop('segments', 64) or 64) # tessellation resolution + tubes = kwds.pop('tubes', False) # output directory if dirname is None: @@ -964,7 +1007,7 @@ def main(instr=None, dirname=None, **kwds): dirname = str(p.parent.joinpath(f"{p.stem}_{dt.strftime(dt.now(), '%Y%m%d_%H%M%S')}")) from mccodelib.mcdisplayutils import McDisplayReader - reader = McDisplayReader(dir=dirname, instr=instr, **kwds) + reader = McDisplayReader(dir=dirname, instr=instr, trace=2, **kwds) instrument = reader.read_instrument() # The reader fails if this is done before reader.read_instrument() above ?! @@ -973,6 +1016,8 @@ def main(instr=None, dirname=None, **kwds): # Convert the instrument into tessellated parts, then save parts = instrument_to_parts(instrument, segments) + if tubes: + lines_to_tubes(parts) output = str(root.joinpath(f"{Path(instr).stem}{EXPORTERS[fmt][1]}")) EXPORTERS[fmt][0](parts, output) print(f"Wrote {output} ({len(parts)} components)") @@ -999,6 +1044,7 @@ def output_format(astring): parser.add_argument('--dirname', help='output directory name override') parser.add_argument('-n', '--ncount', dest='n', type=float, default=0, help='Number of particles to simulate') parser.add_argument('-s', '--segments', dest='segments', type=int, default=16, help='Tessellation resolution (circle segments); lower = smaller/lighter files, e.g. 32 or 16. Default: 16') + parser.add_argument('--tubes', action='store_true', help='Draw lines (multiline etc.) as solid cylinders with spheres at the corners, as in mcdisplay-plopp') parser.add_argument('-f', '--format', dest='format', type=output_format, default='glb', help='Output CAD file format (STEP, STL, VRML, GLTF, GLB, OBJ, PLY, XML/X3D, VTK and VTKJS). Default: GLB') args, unknown = parser.parse_known_args() From 810ae94ac237b6528fbab984e5fcf0d32e5d627a Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Sun, 18 Oct 2026 18:25:15 +0200 Subject: [PATCH 5/8] mcdisplay-cad: export particle rays, 12-digit STEP coordinates - with -n > 0, read the traced rays and export them as one grey 'rays' part (instrument coordinates) in every output format; --tubes applies to them too - STEP: write REALs with 12 significant digits; at 6, short ray segments far from the origin collapsed into degenerate faces and OpenCASCADE dropped them - --tubes help: note that tubed rays can make files very large Co-Authored-By: Claude Opus 5.5 --- tools/Python/mcdisplay/cad/mcdisplay.py | 21 ++++++++++++++++++--- 1 file changed, 18 insertions(+), 3 deletions(-) diff --git a/tools/Python/mcdisplay/cad/mcdisplay.py b/tools/Python/mcdisplay/cad/mcdisplay.py index c44f5253ea..8ad3d929fa 100755 --- a/tools/Python/mcdisplay/cad/mcdisplay.py +++ b/tools/Python/mcdisplay/cad/mcdisplay.py @@ -370,6 +370,18 @@ def component_to_part(comp, color, segments=64): return part.transformed(comp.m4) +def rays_to_part(inst, bundle): + ''' one grey part holding every particle ray as a polyline, in instrument coordinates ''' + m4 = {c.name: c.m4 for c in inst.components} + part = Part('rays', (0.75, 0.75, 0.75)) + for story in bundle.rays: + pts = [transform_points(m4[g.compname], [e.args[:3] for e in g.events]) + for g in story.groups if g.compname in m4 and g.events] + if pts: + part.add_polyline(np.vstack(pts)) + return part + + def instrument_to_parts(inst, segments=64): ''' convert components to parts, assigning a cycling color per component ''' parts = [] @@ -792,7 +804,7 @@ def E(text): def r(x): # STEP REAL needs a decimal point: 0 -> 0., 1E-05 -> 1.E-05 - m, e, x = f'{x:.6G}'.partition('E') + m, e, x = f'{x:.12G}'.partition('E') return (m if '.' in m else m + '.') + e + x def r3(p): @@ -1009,6 +1021,7 @@ def main(instr=None, dirname=None, **kwds): from mccodelib.mcdisplayutils import McDisplayReader reader = McDisplayReader(dir=dirname, instr=instr, trace=2, **kwds) instrument = reader.read_instrument() + bundle = reader.read_particles() # None unless -n > 0 # The reader fails if this is done before reader.read_instrument() above ?! root = Path(dirname) @@ -1016,6 +1029,8 @@ def main(instr=None, dirname=None, **kwds): # Convert the instrument into tessellated parts, then save parts = instrument_to_parts(instrument, segments) + if bundle and bundle.rays: + parts.append(rays_to_part(instrument, bundle)) if tubes: lines_to_tubes(parts) output = str(root.joinpath(f"{Path(instr).stem}{EXPORTERS[fmt][1]}")) @@ -1042,9 +1057,9 @@ def output_format(astring): parser, prefix = make_common_parser(__file__, __doc__) parser.add_argument('--dirname', help='output directory name override') - parser.add_argument('-n', '--ncount', dest='n', type=float, default=0, help='Number of particles to simulate') + parser.add_argument('-n', '--ncount', dest='n', type=float, default=0, help='Number of particle rays to trace and include in the output (default 0: geometry only)') parser.add_argument('-s', '--segments', dest='segments', type=int, default=16, help='Tessellation resolution (circle segments); lower = smaller/lighter files, e.g. 32 or 16. Default: 16') - parser.add_argument('--tubes', action='store_true', help='Draw lines (multiline etc.) as solid cylinders with spheres at the corners, as in mcdisplay-plopp') + parser.add_argument('--tubes', action='store_true', help='Draw lines (multiline etc.) as solid cylinders with spheres at the corners, as in mcdisplay-plopp. NB: with -n > 0 the particle rays are tubed too, which can make files very large (e.g. ~150 MB STEP for 100 rays in PSI_DMC)') parser.add_argument('-f', '--format', dest='format', type=output_format, default='glb', help='Output CAD file format (STEP, STL, VRML, GLTF, GLB, OBJ, PLY, XML/X3D, VTK and VTKJS). Default: GLB') args, unknown = parser.parse_known_args() From c48d87820bd403b04c6b294cd6b51546bb57b98d Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Sun, 18 Oct 2026 18:38:15 +0200 Subject: [PATCH 6/8] mcdisplay-cad: add -t/--trace switch (default 2) --trace=1 gives the legacy wireframe drawing; the help notes that STL may then be empty and that many viewers ignore lines in OBJ/PLY (use --tubes). Co-Authored-By: Claude Opus 5.5 --- tools/Python/mcdisplay/cad/mcdisplay.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/tools/Python/mcdisplay/cad/mcdisplay.py b/tools/Python/mcdisplay/cad/mcdisplay.py index 8ad3d929fa..e7b163dc21 100755 --- a/tools/Python/mcdisplay/cad/mcdisplay.py +++ b/tools/Python/mcdisplay/cad/mcdisplay.py @@ -1019,7 +1019,8 @@ def main(instr=None, dirname=None, **kwds): dirname = str(p.parent.joinpath(f"{p.stem}_{dt.strftime(dt.now(), '%Y%m%d_%H%M%S')}")) from mccodelib.mcdisplayutils import McDisplayReader - reader = McDisplayReader(dir=dirname, instr=instr, trace=2, **kwds) + kwds.setdefault('trace', 2) + reader = McDisplayReader(dir=dirname, instr=instr, **kwds) instrument = reader.read_instrument() bundle = reader.read_particles() # None unless -n > 0 @@ -1059,6 +1060,7 @@ def output_format(astring): parser.add_argument('--dirname', help='output directory name override') parser.add_argument('-n', '--ncount', dest='n', type=float, default=0, help='Number of particle rays to trace and include in the output (default 0: geometry only)') parser.add_argument('-s', '--segments', dest='segments', type=int, default=16, help='Tessellation resolution (circle segments); lower = smaller/lighter files, e.g. 32 or 16. Default: 16') + parser.add_argument('-t', '--trace', dest='trace', type=int, default=2, help='mcrun --trace level: 2 (default) gives solid primitives, 1 the legacy wireframe drawing. NB: with 1, surface-only formats (STL) may come out empty and many viewers ignore lines in OBJ/PLY; add --tubes to get solids') parser.add_argument('--tubes', action='store_true', help='Draw lines (multiline etc.) as solid cylinders with spheres at the corners, as in mcdisplay-plopp. NB: with -n > 0 the particle rays are tubed too, which can make files very large (e.g. ~150 MB STEP for 100 rays in PSI_DMC)') parser.add_argument('-f', '--format', dest='format', type=output_format, default='glb', help='Output CAD file format (STEP, STL, VRML, GLTF, GLB, OBJ, PLY, XML/X3D, VTK and VTKJS). Default: GLB') From ed0188963500c417e989fd2039cc39bac809b7e6 Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Sun, 18 Oct 2026 18:46:20 +0200 Subject: [PATCH 7/8] mcdisplay-cad: fix scrambled sphere tessellation, add geometry self-check column_stack on the 2-D meshgrid arrays interleaved x/y/z across points, so every sphere (e.g. the Text.comp stroke joints) came out as spikes in all output formats. Flatten the grids first. test_mcdisplay.py: assert-based check of the tessellated primitives (radii, extents, closed surfaces, polyhedron JSON, tubes), as test_mcdisplopp.py. Co-Authored-By: Claude Opus 5.5 --- tools/Python/mcdisplay/cad/mcdisplay.py | 4 +- tools/Python/mcdisplay/cad/test_mcdisplay.py | 41 ++++++++++++++++++++ 2 files changed, 43 insertions(+), 2 deletions(-) create mode 100644 tools/Python/mcdisplay/cad/test_mcdisplay.py diff --git a/tools/Python/mcdisplay/cad/mcdisplay.py b/tools/Python/mcdisplay/cad/mcdisplay.py index e7b163dc21..b0e6f9fb81 100755 --- a/tools/Python/mcdisplay/cad/mcdisplay.py +++ b/tools/Python/mcdisplay/cad/mcdisplay.py @@ -272,10 +272,10 @@ def tessellate(d, segments=64): nu, nv = segments, max(8, segments // 2) u = np.linspace(0, 2 * np.pi, nu, endpoint=False) v = np.linspace(0, np.pi, nv) - uu, vv = np.meshgrid(u, v) + uu, vv = (g.ravel() for g in np.meshgrid(u, v)) # row-major: index i*nu + j pts = np.column_stack((c[0] + r * np.cos(uu) * np.sin(vv), c[1] + r * np.sin(uu) * np.sin(vv), - c[2] + r * np.cos(vv))).reshape(-1, 3) + c[2] + r * np.cos(vv))) tris = [] for i in range(nv - 1): for j in range(nu): diff --git a/tools/Python/mcdisplay/cad/test_mcdisplay.py b/tools/Python/mcdisplay/cad/test_mcdisplay.py new file mode 100644 index 0000000000..7913439bd7 --- /dev/null +++ b/tools/Python/mcdisplay/cad/test_mcdisplay.py @@ -0,0 +1,41 @@ +''' geometry self-check of the mcdisplay-cad tessellation; run: python3 test_mcdisplay.py (needs mccodelib on PYTHONPATH) ''' +import json +from collections import Counter +import numpy as np +from mccodelib.instrgeom import DrawSphere, DrawCylinder, DrawCone, DrawBox, DrawDisc, DrawPolyhedron +from mcdisplay import tessellate, _tubes + + +def closed(m): # welded by position, every edge shared by exactly two triangles (degenerate pole/cap triangles dropped) + f = np.unique(np.round(m.verts, 9), axis=0, return_inverse=True)[1].reshape(-1)[m.tris] + f = f[(f[:, 0] != f[:, 1]) & (f[:, 1] != f[:, 2]) & (f[:, 0] != f[:, 2])] + return set(Counter(tuple(sorted(e)) for t in f for e in ((t[0], t[1]), (t[1], t[2]), (t[2], t[0]))).values()) == {2} + + +m, _ = tessellate(DrawSphere([1, 2, 3, 0.5]), 16) +assert np.allclose(np.linalg.norm(m.verts - (1, 2, 3), axis=1), 0.5) and closed(m) + +a = np.array([1, 2, 3]) / np.sqrt(14) +m, _ = tessellate(DrawCylinder([1, 2, 3, 0.2, 1.0, 0, *a]), 16) # tilted axis, as in Text.comp strokes +off = m.verts - (1, 2, 3) +h = off @ a +rad = np.linalg.norm(off - np.outer(h, a), axis=1) +assert np.allclose(np.abs(h), 0.5) and (np.isclose(rad, 0.2) | np.isclose(rad, 0)).all() and closed(m) + +m, _ = tessellate(DrawCone([0, 0, 0, 0.3, 1.0, 0, 0, 1]), 16) +assert np.isclose(m.verts[:, 2].max(), 1.0) and np.allclose(np.linalg.norm(m.verts[:, :2], axis=1).max(), 0.3) + +m, _ = tessellate(DrawBox([1, 2, 3, 0.2, 0.4, 0.6, 0, 0, 1, 0]), 16) +assert np.allclose(m.verts.min(0), (0.9, 1.8, 2.7)) and np.allclose(m.verts.max(0), (1.1, 2.2, 3.3)) and closed(m) + +m, _ = tessellate(DrawDisc([0, 0, 1, 0.5, 0, 0, 1]), 16) +assert np.allclose(m.verts[:, 2], 1) and np.isclose(np.linalg.norm(m.verts[:, :2], axis=1).max(), 0.5) + +js = json.dumps({'vertices': [[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1]], + 'faces': [{'face': [0, 2, 1]}, {'face': [0, 1, 3]}, {'face': [0, 3, 2]}, {'face': [1, 2, 3]}]}) +m, _ = tessellate(DrawPolyhedron([js]), 16) # tetrahedron, as emitted at --trace=2 +assert len(m.tris) == 4 and closed(m) + +m = _tubes(np.array([(0, 0, 0), (1, 0, 0), (1, 1, 0.)]), r=0.1) # 2 segments + 1 corner sphere +assert len(m.tris) == 2 * 12 + 20 and m.tris.max() == len(m.verts) - 1 +print('ok') From 962b25d641ce11a0c43477d214b89cf269a92564 Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Sun, 18 Oct 2026 18:51:25 +0200 Subject: [PATCH 8/8] mcdisplay-cad: default to --trace=1 for speed (-t 2 for solid primitives) Co-Authored-By: Claude Opus 5.5 --- tools/Python/mcdisplay/cad/mcdisplay.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/tools/Python/mcdisplay/cad/mcdisplay.py b/tools/Python/mcdisplay/cad/mcdisplay.py index b0e6f9fb81..b3603bdfc7 100755 --- a/tools/Python/mcdisplay/cad/mcdisplay.py +++ b/tools/Python/mcdisplay/cad/mcdisplay.py @@ -1019,7 +1019,6 @@ def main(instr=None, dirname=None, **kwds): dirname = str(p.parent.joinpath(f"{p.stem}_{dt.strftime(dt.now(), '%Y%m%d_%H%M%S')}")) from mccodelib.mcdisplayutils import McDisplayReader - kwds.setdefault('trace', 2) reader = McDisplayReader(dir=dirname, instr=instr, **kwds) instrument = reader.read_instrument() bundle = reader.read_particles() # None unless -n > 0 @@ -1060,7 +1059,7 @@ def output_format(astring): parser.add_argument('--dirname', help='output directory name override') parser.add_argument('-n', '--ncount', dest='n', type=float, default=0, help='Number of particle rays to trace and include in the output (default 0: geometry only)') parser.add_argument('-s', '--segments', dest='segments', type=int, default=16, help='Tessellation resolution (circle segments); lower = smaller/lighter files, e.g. 32 or 16. Default: 16') - parser.add_argument('-t', '--trace', dest='trace', type=int, default=2, help='mcrun --trace level: 2 (default) gives solid primitives, 1 the legacy wireframe drawing. NB: with 1, surface-only formats (STL) may come out empty and many viewers ignore lines in OBJ/PLY; add --tubes to get solids') + parser.add_argument('-t', '--trace', dest='trace', type=int, default=1, help='mcrun --trace level: 1 (default, fastest) gives the legacy wireframe drawing, 2 solid primitives. NB: with 1, surface-only formats (STL) may come out empty and many viewers ignore lines in OBJ/PLY; add --tubes to get solids') parser.add_argument('--tubes', action='store_true', help='Draw lines (multiline etc.) as solid cylinders with spheres at the corners, as in mcdisplay-plopp. NB: with -n > 0 the particle rays are tubed too, which can make files very large (e.g. ~150 MB STEP for 100 rays in PSI_DMC)') parser.add_argument('-f', '--format', dest='format', type=output_format, default='glb', help='Output CAD file format (STEP, STL, VRML, GLTF, GLB, OBJ, PLY, XML/X3D, VTK and VTKJS). Default: GLB')