From 4f3fb1b708f4eac0ebd87f4fdd40bd4e1bc4d93c Mon Sep 17 00:00:00 2001 From: Yves Chemisky Date: Tue, 24 Feb 2026 14:22:58 +0100 Subject: [PATCH 01/12] Export ImplicitShape and batch_smooth_union Add the implicit_shape module to the shape package and re-export ImplicitShape, the implicit_shape symbol, and batch_smooth_union at the microgen top-level. Update imports and __all__ in microgen/__init__.py and microgen/shape/__init__.py so implicit shapes and the batch smoothing helper are publicly available. --- microgen/__init__.py | 6 ++++++ microgen/shape/__init__.py | 6 +++++- 2 files changed, 11 insertions(+), 1 deletion(-) diff --git a/microgen/__init__.py b/microgen/__init__.py index e5f66939..76ff6ec6 100644 --- a/microgen/__init__.py +++ b/microgen/__init__.py @@ -47,7 +47,9 @@ Ellipsoid, ExtrudedPolygon, FaceCenteredCubic, + ImplicitShape, Infill, + batch_smooth_union, NormedDistance, Octahedron, OctetTruss, @@ -60,6 +62,7 @@ TruncatedCube, TruncatedCuboctahedron, TruncatedOctahedron, + implicit_shape, new_geometry, newGeometry, surface_functions, @@ -82,6 +85,7 @@ "Ellipsoid", "ExtrudedPolygon", "FaceCenteredCubic", + "ImplicitShape", "Infill", "NormedDistance", "Mmg", @@ -101,6 +105,8 @@ "CylindricalTpms", "SphericalTpms", "check_if_only_linear_tetrahedral", + "implicit_shape", + "batch_smooth_union", "cut_phase_by_shape_list", "cutPhaseByShapeList", "cut_phases", diff --git a/microgen/shape/__init__.py b/microgen/shape/__init__.py index d7fc11c3..67ce83ac 100644 --- a/microgen/shape/__init__.py +++ b/microgen/shape/__init__.py @@ -21,12 +21,13 @@ from typing import TYPE_CHECKING, Callable, Literal, Sequence, Tuple -from . import surface_functions +from . import implicit_shape, surface_functions from .box import Box from .capsule import Capsule from .cylinder import Cylinder from .ellipsoid import Ellipsoid from .extruded_polygon import ExtrudedPolygon +from .implicit_shape import ImplicitShape, batch_smooth_union from .polyhedron import Polyhedron from .shape import Shape from .sphere import Sphere @@ -160,7 +161,9 @@ def __init__(self: ShapeError, shape: str) -> None: "Ellipsoid", "ExtrudedPolygon", "FaceCenteredCubic", + "ImplicitShape", "Infill", + "batch_smooth_union", "NormedDistance", "Octahedron", "OctetTruss", @@ -174,6 +177,7 @@ def __init__(self: ShapeError, shape: str) -> None: "TruncatedCube", "TruncatedCuboctahedron", "TruncatedOctahedron", + "implicit_shape", "new_geometry", "newGeometry", "surface_functions", From 0d4f1f11ebe165feaf53b854892740c98d19baf2 Mon Sep 17 00:00:00 2001 From: Yves Chemisky Date: Tue, 24 Feb 2026 14:23:12 +0100 Subject: [PATCH 02/12] Add implicit shapes (F-rep) documentation Add documentation for the implicit shapes (F-rep) feature: a new examples page (docs/examples/implicit_shapes.rst) with included example scripts and images, and a comprehensive tutorials section (docs/tutorials.rst) demonstrating primitives, booleans, smooth blending, batch union, transforms, utilities, and TPMS composition with runnable Jupyter cells. Also expose microgen.shape.implicit_shape in the API docs (docs/microgen.shape.rst) and add the new examples page to the docs index (docs/index.rst). This surfaces the new implicit_shape API and usage patterns for users. --- docs/examples/implicit_shapes.rst | 28 ++++++ docs/index.rst | 1 + docs/microgen.shape.rst | 5 ++ docs/tutorials.rst | 140 ++++++++++++++++++++++++++++++ 4 files changed, 174 insertions(+) create mode 100644 docs/examples/implicit_shapes.rst diff --git a/docs/examples/implicit_shapes.rst b/docs/examples/implicit_shapes.rst new file mode 100644 index 00000000..652bb82c --- /dev/null +++ b/docs/examples/implicit_shapes.rst @@ -0,0 +1,28 @@ +.. _RST ImplicitShapes: + +Implicit Shapes (F-rep) +======================= + +Primitives and Boolean Operations +---------------------------------- + +.. include:: ../../examples/ImplicitShapes/primitives_and_booleans.py + :literal: + +.. image:: ../_static/examples/implicit_booleans.png + +Implicit Lattices +----------------- + +.. include:: ../../examples/Lattices/implicit_lattice.py + :literal: + +.. image:: ../_static/examples/implicit_lattice.png + +Trabecular Bone +--------------- + +.. include:: ../../examples/ImplicitShapes/trabecular_bone.py + :literal: + +.. image:: ../_static/examples/trabecular_bone.png diff --git a/docs/index.rst b/docs/index.rst index 42652012..c2d7f58a 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -39,6 +39,7 @@ Table of contents examples/basic_shapes examples/lattices + examples/implicit_shapes examples/tpms examples/3d_operations examples/mesh diff --git a/docs/microgen.shape.rst b/docs/microgen.shape.rst index 0a0cc052..e4dbf9cf 100644 --- a/docs/microgen.shape.rst +++ b/docs/microgen.shape.rst @@ -52,6 +52,11 @@ microgen.shape package :undoc-members: :show-inheritance: +.. automodule:: microgen.shape.implicit_shape + :members: + :undoc-members: + :show-inheritance: + .. automodule:: microgen.shape.surface_functions :members: :undoc-members: diff --git a/docs/tutorials.rst b/docs/tutorials.rst index 92022584..4c956b5f 100644 --- a/docs/tutorials.rst +++ b/docs/tutorials.rst @@ -260,6 +260,146 @@ Create TPMS on a cylindrical coordinate system: cylindrical_gyroid.sheet.plot(color='white') +Implicit Shapes (F-rep) +----------------------- + +Microgen includes an implicit (F-rep) modeling framework where shapes are +defined by scalar fields ``f(x, y, z)`` — the surface is at ``f = 0`` and the +interior is ``f < 0``. This enables smooth blending, exact booleans, and +composability that is hard to achieve with boundary meshes alone. + +Primitives +^^^^^^^^^^ + +Create basic implicit primitives and visualize them with marching cubes: + +.. jupyter-execute:: + + from microgen.shape.implicit_shape import implicit_sphere, implicit_box + + sphere = implicit_sphere(center=(0, 0, 0), radius=0.4) + sphere.generate_vtk(resolution=60).plot(color='white') + +.. jupyter-execute:: + + box = implicit_box(center=(0, 0, 0), half_extents=(0.3, 0.3, 0.3)) + box.generate_vtk(resolution=60).plot(color='white') + + +Boolean Operations +^^^^^^^^^^^^^^^^^^ + +Combine implicit shapes with standard boolean operators: + +.. jupyter-execute:: + + from microgen.shape.implicit_shape import implicit_sphere, implicit_box + + sphere = implicit_sphere(center=(0, 0, 0), radius=0.4) + box = implicit_box(center=(0, 0, 0), half_extents=(0.3, 0.3, 0.3)) + + union = sphere | box + union.generate_vtk(resolution=60).plot(color='steelblue') + +.. jupyter-execute:: + + intersection = sphere & box + intersection.generate_vtk(resolution=60).plot(color='coral') + +.. jupyter-execute:: + + difference = sphere - box + difference.generate_vtk(resolution=60).plot(color='mediumseagreen') + + +Smooth Blending +^^^^^^^^^^^^^^^ + +Smooth union produces organic-looking junctions controlled by the parameter *k*: + +.. jupyter-execute:: + + smooth = sphere.smooth_union(box, k=0.3) + smooth.generate_vtk(resolution=60).plot(color='orchid') + + +Batch Smooth Union +^^^^^^^^^^^^^^^^^^ + +When combining many primitives, use ``batch_smooth_union`` — it evaluates all +fields in a flat loop, avoiding recursion-depth limits: + +.. jupyter-execute:: + + import numpy as np + from microgen.shape.implicit_shape import implicit_sphere, batch_smooth_union + + rng = np.random.default_rng(42) + spheres = [ + implicit_sphere(center=tuple(rng.uniform(-0.5, 0.5, 3)), radius=0.15) + for _ in range(20) + ] + combined = batch_smooth_union(spheres, k=0.1) + combined.generate_vtk(bounds=(-0.8, 0.8) * 3, resolution=80).plot(color='coral') + + +Transforms +^^^^^^^^^^ + +Translate, rotate, and scale implicit shapes: + +.. jupyter-execute:: + + from microgen.shape.implicit_shape import implicit_box + + box = implicit_box(center=(0, 0, 0), half_extents=(0.3, 0.2, 0.1)) + rotated = box.rotate(angles=(0, 0, 45)) + rotated.generate_vtk(resolution=60).plot(color='white') + + +Utilities: Shell and Repeat +^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +``shell`` hollows a shape; ``repeat`` tiles it periodically: + +.. jupyter-execute:: + + from microgen.shape.implicit_shape import implicit_sphere, shell + + sphere = implicit_sphere(center=(0, 0, 0), radius=0.4) + hollow = shell(sphere, thickness=0.05) + hollow.generate_vtk(resolution=80).plot(color='white') + +.. jupyter-execute:: + + from microgen.shape.implicit_shape import implicit_sphere, repeat + + sphere = implicit_sphere(center=(0, 0, 0), radius=0.15) + tiled = repeat(sphere, spacing=(0.5, 0.5, 0.5)) + tiled.generate_vtk(bounds=(-0.8, 0.8) * 3, resolution=80).plot(color='white') + + +TPMS + F-rep Composition +^^^^^^^^^^^^^^^^^^^^^^^^^ + +Combine TPMS implicit fields with F-rep shapes — for example, intersect a +gyroid with an implicit sphere: + +.. jupyter-execute:: + + from microgen.shape.implicit_shape import implicit_sphere, from_field + from microgen.shape.surface_functions import gyroid + import numpy as np + + gyroid_field = from_field( + lambda x, y, z: gyroid(x * 2 * np.pi, y * 2 * np.pi, z * 2 * np.pi), + bounds=(-0.6, 0.6) * 3, + ) + sphere = implicit_sphere(center=(0, 0, 0), radius=0.5) + result = gyroid_field & sphere + result.generate_vtk(bounds=(-0.6, 0.6) * 3, resolution=80).plot(color='white') + + 3D Operations ------------- From ae3b5fb3653f1fc1209ef9f00680b8693eb72b28 Mon Sep 17 00:00:00 2001 From: Yves Chemisky Date: Tue, 24 Feb 2026 14:23:30 +0100 Subject: [PATCH 03/12] Add implicit shapes and lattice examples Add three new example scripts demonstrating the F-rep implicit modeling utilities and comparisons to strut lattices. - examples/ImplicitShapes/primitives_and_booleans.py: Shows primitive implicit shapes (sphere, box, capsule), hard boolean ops (union, intersection, difference) and smooth blending, plus a batch smooth-union demo combining many spheres. - examples/ImplicitShapes/trabecular_bone.py: Generates a Voronoi-based trabecular-bone-like scaffold (with optional Lloyd relaxation), builds capsules and joint spheres from Voronoi edges/vertices, combines them with batch smooth union and intersects with a clip box, then meshes for visualization. - examples/Lattices/implicit_lattice.py: Builds implicit cubic, BCC and octet-truss lattices from capsule/sphere primitives (batch smooth union), compares implicit vs original strut meshes, and sweeps the smoothness parameter k for visual comparison. Examples rely on numpy, pyvista, scipy.spatial.Voronoi and the microgen implicit-shape utilities for construction and visualization. --- .../ImplicitShapes/primitives_and_booleans.py | 75 +++++++ examples/ImplicitShapes/trabecular_bone.py | 163 +++++++++++++++ examples/Lattices/implicit_lattice.py | 187 ++++++++++++++++++ 3 files changed, 425 insertions(+) create mode 100644 examples/ImplicitShapes/primitives_and_booleans.py create mode 100644 examples/ImplicitShapes/trabecular_bone.py create mode 100644 examples/Lattices/implicit_lattice.py diff --git a/examples/ImplicitShapes/primitives_and_booleans.py b/examples/ImplicitShapes/primitives_and_booleans.py new file mode 100644 index 00000000..7be59380 --- /dev/null +++ b/examples/ImplicitShapes/primitives_and_booleans.py @@ -0,0 +1,75 @@ +"""Implicit shapes: primitives, boolean operations, and smooth blending. + +Demonstrates the F-rep (Function Representation) framework in microgen +for creating implicit shapes and combining them with hard and smooth +boolean operations. +""" + +import pyvista as pv + +from microgen.shape.implicit_shape import ( + batch_smooth_union, + implicit_box, + implicit_sphere, +) + +# --------------------------------------------------------------------------- +# Create primitives +# --------------------------------------------------------------------------- + +sphere = implicit_sphere(center=(0, 0, 0), radius=0.45) +box = implicit_box(center=(0, 0, 0), half_extents=(0.3, 0.3, 0.3)) + +BOUNDS = (-0.8, 0.8) * 3 +RES = 80 + +# --------------------------------------------------------------------------- +# Boolean operations +# --------------------------------------------------------------------------- + +union = sphere | box +intersection = sphere & box +difference = sphere - box +smooth = sphere.smooth_union(box, k=0.3) + +# --------------------------------------------------------------------------- +# Visualize in a 2x2 grid +# --------------------------------------------------------------------------- + +plotter = pv.Plotter(shape=(2, 2), window_size=(1200, 1200)) + +ops = [ + (0, 0, "Union (A | B)", union, "steelblue"), + (0, 1, "Intersection (A & B)", intersection, "coral"), + (1, 0, "Difference (A - B)", difference, "mediumseagreen"), + (1, 1, "Smooth Union (k=0.3)", smooth, "orchid"), +] + +for row, col, title, shape, color in ops: + mesh = shape.generate_vtk(bounds=BOUNDS, resolution=RES) + plotter.subplot(row, col) + plotter.add_text(title, font_size=12) + plotter.add_mesh(mesh, color=color, show_edges=False) + +plotter.link_views() +plotter.show() + +# --------------------------------------------------------------------------- +# Batch smooth union — combining many primitives efficiently +# --------------------------------------------------------------------------- + +import numpy as np + +rng = np.random.default_rng(42) +spheres = [ + implicit_sphere(center=tuple(rng.uniform(-0.5, 0.5, 3)), radius=0.15) + for _ in range(20) +] + +combined = batch_smooth_union(spheres, k=0.1) +mesh = combined.generate_vtk(bounds=(-0.8, 0.8) * 3, resolution=100) + +plotter2 = pv.Plotter() +plotter2.add_text("Batch smooth union (20 spheres, k=0.1)", font_size=12) +plotter2.add_mesh(mesh, color="coral", show_edges=False) +plotter2.show() diff --git a/examples/ImplicitShapes/trabecular_bone.py b/examples/ImplicitShapes/trabecular_bone.py new file mode 100644 index 00000000..0b785613 --- /dev/null +++ b/examples/ImplicitShapes/trabecular_bone.py @@ -0,0 +1,163 @@ +"""Generate a trabecular bone-like structure from Voronoi edges. + +Seed points are scattered in a cube, the 3D Voronoi tessellation provides +the connectivity (edges with natural degree-4 nodes), and each edge becomes +an implicit capsule. All primitives are combined with smooth union and +intersected with a box for a clean cross-section. + +Reference: Voronoi tessellation is the standard approach in literature for +synthetic trabecular bone scaffolds — interior Voronoi vertices have exactly +4 edges, matching real trabecular node valence. +""" + +import numpy as np +import pyvista as pv +from scipy.spatial import Voronoi + +from microgen.shape.implicit_shape import ( + batch_smooth_union, + implicit_box, + implicit_capsule, + implicit_sphere, +) + +# --------------------------------------------------------------------------- +# Parameters +# --------------------------------------------------------------------------- + +SEED = 42 +N_POINTS = 150 # Voronoi seed count + +# Domain +DOMAIN_SIZE = 2.0 +HALF = DOMAIN_SIZE / 2.0 + +# Clip box — the clean cut +BOX_HALF = 0.75 + +# Strut geometry (trabecular dimensions) +RADIUS_MEAN = 0.035 # ~Tb.Th / 2 +RADIUS_STD = 0.006 + +# Smooth blending +K = 0.03 + +# Meshing +RESOLUTION = 200 +BOUNDS = (-BOX_HALF * 1.02, BOX_HALF * 1.02) * 3 + +# Lloyd relaxation iterations (regularizes cell sizes while keeping randomness) +LLOYD_ITERATIONS = 2 + +# --------------------------------------------------------------------------- +# Generate seed points + Lloyd relaxation +# --------------------------------------------------------------------------- + +rng = np.random.default_rng(SEED) + +# Scatter seeds in a region larger than the clip box +seeds = rng.uniform(-HALF, HALF, size=(N_POINTS, 3)) + + +def lloyd_relaxation(pts, bounds_min, bounds_max, iterations=1): + """Move each seed toward the centroid of its Voronoi cell.""" + for _ in range(iterations): + vor = Voronoi(pts) + new_pts = np.copy(pts) + for i, region_idx in enumerate(vor.point_region): + region = vor.regions[region_idx] + if -1 in region or len(region) == 0: + continue # skip open/infinite regions + verts = vor.vertices[region] + centroid = verts.mean(axis=0) + centroid = np.clip(centroid, bounds_min, bounds_max) + new_pts[i] = centroid + pts = new_pts + return pts + + +print(f"Lloyd relaxation ({LLOYD_ITERATIONS} iterations) on {N_POINTS} seeds ...") +seeds = lloyd_relaxation(seeds, -HALF, HALF, iterations=LLOYD_ITERATIONS) + +# --------------------------------------------------------------------------- +# Compute Voronoi and extract edges inside the clip box +# --------------------------------------------------------------------------- + +print("Computing Voronoi tessellation ...") +vor = Voronoi(seeds) + +margin = BOX_HALF + 0.15 +edge_set = set() +for ridge in vor.ridge_vertices: + if -1 in ridge: + continue + n = len(ridge) + for i in range(n): + v1, v2 = ridge[i], ridge[(i + 1) % n] + edge_set.add((min(v1, v2), max(v1, v2))) + +# Filter: keep edges near the clip box +edges = [] +for v1, v2 in edge_set: + p1 = vor.vertices[v1] + p2 = vor.vertices[v2] + if np.all(np.abs(p1) < margin) or np.all(np.abs(p2) < margin): + edges.append((v1, v2)) + +print(f" {len(edges)} finite edges near clip box") + +# Collect unique vertex indices used +used_verts = set() +for v1, v2 in edges: + used_verts.add(v1) + used_verts.add(v2) + +print(f" {len(used_verts)} Voronoi vertices (nodes)") + +# --------------------------------------------------------------------------- +# Build implicit primitives +# --------------------------------------------------------------------------- + +print(f"Building {len(edges)} capsules + {len(used_verts)} joint spheres ...") + +radii = np.clip(rng.normal(RADIUS_MEAN, RADIUS_STD, size=len(edges)), 0.015, None) + +primitives = [] + +for idx, (v1, v2) in enumerate(edges): + primitives.append( + implicit_capsule( + start=tuple(vor.vertices[v1]), + end=tuple(vor.vertices[v2]), + radius=float(radii[idx]), + ) + ) + +for vi in used_verts: + primitives.append( + implicit_sphere(center=tuple(vor.vertices[vi]), radius=RADIUS_MEAN * 1.1) + ) + +# --------------------------------------------------------------------------- +# Combine with smooth union + box intersection +# --------------------------------------------------------------------------- + +print(f"Combining {len(primitives)} primitives with batch smooth union (k={K}) ...") + +result = batch_smooth_union(primitives, k=K) + +print(f"Intersecting with box (half-extent={BOX_HALF}) ...") +box = implicit_box(center=(0, 0, 0), half_extents=(BOX_HALF, BOX_HALF, BOX_HALF)) +result = result & box + +# --------------------------------------------------------------------------- +# Mesh and visualize +# --------------------------------------------------------------------------- + +print(f"Meshing (resolution={RESOLUTION}) ...") +mesh = result.generate_vtk(bounds=BOUNDS, resolution=RESOLUTION) + +plotter = pv.Plotter() +plotter.add_text("Trabecular Bone (Voronoi)", font_size=14) +plotter.add_mesh(mesh, color="ivory", show_edges=False) +plotter.show() diff --git a/examples/Lattices/implicit_lattice.py b/examples/Lattices/implicit_lattice.py new file mode 100644 index 00000000..8d4b9e60 --- /dev/null +++ b/examples/Lattices/implicit_lattice.py @@ -0,0 +1,187 @@ +"""Implicit (F-rep) lattice construction and comparison with strut lattices. + +This example builds BCC, Octet-Truss, and Cubic lattices using the F-rep +implicit modeling framework (capsule primitives + smooth union) and compares +them side-by-side with the original strut-based lattice meshes from microgen. + +A second comparison shows the effect of varying the smoothness parameter *k* +on a BCC lattice — from sharp joints (k ≈ 0) to heavily rounded blends. +""" + +from itertools import product + +import numpy as np +import pyvista as pv + +from microgen import BodyCenteredCubic, Cubic, OctetTruss +from microgen.shape.implicit_shape import ( + ImplicitShape, + batch_smooth_union, + implicit_capsule, + implicit_sphere, +) + +# --------------------------------------------------------------------------- +# Helpers +# --------------------------------------------------------------------------- + +CELL = 1.0 # unit cell size +HALF = CELL / 2.0 +STRUT_RADIUS = 0.08 +RESOLUTION = 80 +BOUNDS = (-HALF * 1.05, HALF * 1.05) * 3 # slightly larger than unit cube + + +def _cube_corners(): + return np.array(list(product([-HALF, HALF], repeat=3))) + + +def _face_centers(): + centers = [] + for axis in range(3): + for sign in [-1, 1]: + pt = [0.0, 0.0, 0.0] + pt[axis] = sign * HALF + centers.append(pt) + return np.array(centers) + + +def _build_implicit_lattice( + vertices: np.ndarray, + pairs: list[tuple[int, int]], + radius: float, + k: float = 0.0, + add_joints: bool = True, +) -> ImplicitShape: + """Build an implicit lattice from vertices and connectivity. + + Each strut is an implicit capsule; joints are implicit spheres. + All primitives are combined with ``batch_smooth_union`` to avoid + recursion-depth issues with large primitive counts. + """ + primitives: list[ImplicitShape] = [] + + for i, j in pairs: + primitives.append( + implicit_capsule( + start=tuple(vertices[i]), + end=tuple(vertices[j]), + radius=radius, + ) + ) + + if add_joints: + for v in vertices: + primitives.append(implicit_sphere(center=tuple(v), radius=radius)) + + return batch_smooth_union(primitives, k=k) + + +# --------------------------------------------------------------------------- +# Lattice definitions (same topology as microgen strut lattices) +# --------------------------------------------------------------------------- + + +def cubic_lattice( + radius: float = STRUT_RADIUS, k: float = 0.0, +) -> ImplicitShape: + """Cubic lattice: 8 corner nodes, 12 edge struts.""" + verts = _cube_corners() + pairs = [] + for i in range(len(verts)): + for j in range(i + 1, len(verts)): + if np.linalg.norm(verts[i] - verts[j]) < CELL + 1e-6: + pairs.append((i, j)) + return _build_implicit_lattice(verts, pairs, radius, k) + + +def bcc_lattice( + radius: float = STRUT_RADIUS, k: float = 0.0, +) -> ImplicitShape: + """Body-centered cubic: center node connected to 8 corners.""" + corners = _cube_corners() + verts = np.vstack(([[0.0, 0.0, 0.0]], corners)) + pairs = [(0, i) for i in range(1, 9)] + return _build_implicit_lattice(verts, pairs, radius, k) + + +def octet_truss_lattice( + radius: float = STRUT_RADIUS, k: float = 0.0, +) -> ImplicitShape: + """Octet-truss: 8 corners + 6 face centers, connected at distance √2/2.""" + corners = _cube_corners() + faces = _face_centers() + verts = np.vstack((corners, faces)) + connection_dist = CELL / np.sqrt(2) + 1e-5 + pairs = [] + for i in range(len(verts)): + for j in range(i + 1, len(verts)): + if np.linalg.norm(verts[i] - verts[j]) < connection_dist: + pairs.append((i, j)) + return _build_implicit_lattice(verts, pairs, radius, k) + + +# --------------------------------------------------------------------------- +# 1. Comparison: original strut lattice vs. implicit lattice (k=0 and k>0) +# --------------------------------------------------------------------------- + +print("Building comparison: strut lattices vs implicit lattices ...") + +lattice_defs = [ + ("Cubic", Cubic, cubic_lattice), + ("BCC", BodyCenteredCubic, bcc_lattice), + ("Octet-Truss", OctetTruss, octet_truss_lattice), +] + +plotter = pv.Plotter(shape=(len(lattice_defs), 3), window_size=(1800, 1800)) + +for row, (name, StrutClass, implicit_fn) in enumerate(lattice_defs): + # Column 0: original strut lattice + strut = StrutClass(strut_radius=STRUT_RADIUS) + strut_mesh = strut.generate_vtk() + + plotter.subplot(row, 0) + plotter.add_text(f"{name}\n(strut)", font_size=10) + plotter.add_mesh(strut_mesh, color="steelblue", show_edges=False) + + # Column 1: implicit lattice, hard union (k=0) + imp_hard = implicit_fn(radius=STRUT_RADIUS, k=0.0) + mesh_hard = imp_hard.generate_vtk(bounds=BOUNDS, resolution=RESOLUTION) + + plotter.subplot(row, 1) + plotter.add_text(f"{name}\nimplicit k=0", font_size=10) + plotter.add_mesh(mesh_hard, color="coral", show_edges=False) + + # Column 2: implicit lattice, smooth union (k=0.15) + imp_smooth = implicit_fn(radius=STRUT_RADIUS, k=0.15) + mesh_smooth = imp_smooth.generate_vtk(bounds=BOUNDS, resolution=RESOLUTION) + + plotter.subplot(row, 2) + plotter.add_text(f"{name}\nimplicit k=0.15", font_size=10) + plotter.add_mesh(mesh_smooth, color="mediumseagreen", show_edges=False) + +plotter.link_views() +plotter.show() + + +# --------------------------------------------------------------------------- +# 2. Smoothness sweep on BCC lattice +# --------------------------------------------------------------------------- + +print("Building smoothness sweep on BCC lattice ...") + +k_values = [0.0, 0.05, 0.10, 0.20, 0.35] +n_k = len(k_values) + +plotter2 = pv.Plotter(shape=(1, n_k), window_size=(400 * n_k, 500)) + +for col, k in enumerate(k_values): + bcc = bcc_lattice(radius=STRUT_RADIUS, k=k) + mesh = bcc.generate_vtk(bounds=BOUNDS, resolution=RESOLUTION) + + plotter2.subplot(0, col) + plotter2.add_text(f"k = {k:.2f}", font_size=12) + plotter2.add_mesh(mesh, color="coral", show_edges=False) + +plotter2.link_views() +plotter2.show() From bfcfbccd0f04db67959effbbdf3b46e3cdfbf00c Mon Sep 17 00:00:00 2001 From: Yves Chemisky Date: Tue, 24 Feb 2026 14:57:49 +0100 Subject: [PATCH 04/12] Update microgen.shape.rst --- docs/microgen.shape.rst | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/docs/microgen.shape.rst b/docs/microgen.shape.rst index e4dbf9cf..b2ea6118 100644 --- a/docs/microgen.shape.rst +++ b/docs/microgen.shape.rst @@ -57,6 +57,11 @@ microgen.shape package :undoc-members: :show-inheritance: +.. automodule:: microgen.shape.implicit_basic_factory + :members: + :undoc-members: + :show-inheritance: + .. automodule:: microgen.shape.surface_functions :members: :undoc-members: From 8f2d7967972251c49b230801195ff9e36243b12c Mon Sep 17 00:00:00 2001 From: Yves Chemisky Date: Tue, 24 Feb 2026 14:57:54 +0100 Subject: [PATCH 05/12] Update tutorials.rst --- docs/tutorials.rst | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/docs/tutorials.rst b/docs/tutorials.rst index 4c956b5f..426cc640 100644 --- a/docs/tutorials.rst +++ b/docs/tutorials.rst @@ -268,6 +268,25 @@ defined by scalar fields ``f(x, y, z)`` — the surface is at ``f = 0`` and the interior is ``f < 0``. This enables smooth blending, exact booleans, and composability that is hard to achieve with boundary meshes alone. +Converting Basic Shapes to Implicit +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Basic shapes (``Sphere``, ``Box``, ``Capsule``, ``Cylinder``, ``Ellipsoid``) +can be converted to implicit shapes with ``to_implicit()``, making it easy to +combine parametric geometry with F-rep boolean operations: + +.. jupyter-execute:: + + import microgen + from microgen.shape.implicit_basic_factory import implicit_box + + sphere = microgen.Sphere(center=(0.3, 0, 0), radius=0.4) + sphere_imp = sphere.to_implicit() + box_imp = implicit_box(center=(-0.2, 0, 0), half_extents=(0.3, 0.3, 0.3)) + result = sphere_imp.smooth_union(box_imp, k=0.2) + result.generate_vtk(resolution=60).plot(color='white') + + Primitives ^^^^^^^^^^ From 1ae5f7ab94677bc2367c418add15591a0c4758d8 Mon Sep 17 00:00:00 2001 From: Yves Chemisky Date: Tue, 24 Feb 2026 14:58:08 +0100 Subject: [PATCH 06/12] Export implicit_basic_factory in package and shape Add implicit_basic_factory to the shape package and the top-level microgen package exports. The change imports implicit_basic_factory in microgen/shape/__init__.py and microgen/__init__.py and appends it to their __all__ lists so the factory is publicly available from both microgen.shape and the top-level microgen namespace. --- microgen/__init__.py | 2 ++ microgen/shape/__init__.py | 3 ++- 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/microgen/__init__.py b/microgen/__init__.py index 76ff6ec6..9bc91eae 100644 --- a/microgen/__init__.py +++ b/microgen/__init__.py @@ -62,6 +62,7 @@ TruncatedCube, TruncatedCuboctahedron, TruncatedOctahedron, + implicit_basic_factory, implicit_shape, new_geometry, newGeometry, @@ -105,6 +106,7 @@ "CylindricalTpms", "SphericalTpms", "check_if_only_linear_tetrahedral", + "implicit_basic_factory", "implicit_shape", "batch_smooth_union", "cut_phase_by_shape_list", diff --git a/microgen/shape/__init__.py b/microgen/shape/__init__.py index 67ce83ac..de39e1d5 100644 --- a/microgen/shape/__init__.py +++ b/microgen/shape/__init__.py @@ -21,7 +21,7 @@ from typing import TYPE_CHECKING, Callable, Literal, Sequence, Tuple -from . import implicit_shape, surface_functions +from . import implicit_basic_factory, implicit_shape, surface_functions from .box import Box from .capsule import Capsule from .cylinder import Cylinder @@ -177,6 +177,7 @@ def __init__(self: ShapeError, shape: str) -> None: "TruncatedCube", "TruncatedCuboctahedron", "TruncatedOctahedron", + "implicit_basic_factory", "implicit_shape", "new_geometry", "newGeometry", From 333a7094c76c2803e71651e61d90ae4d7b9fe03c Mon Sep 17 00:00:00 2001 From: Yves Chemisky Date: Tue, 24 Feb 2026 14:58:25 +0100 Subject: [PATCH 07/12] Add implicit shape factory and conversions Introduce a new implicit primitive factory and wire shape-to-implicit conversions. Adds microgen/shape/implicit_basic_factory.py implementing ImplicitShape factories (implicit_sphere, implicit_box, implicit_cylinder, implicit_plane, implicit_torus, implicit_capsule, implicit_ellipsoid) as SDF-based functions with reasonable bounds. Adds to_implicit methods to Box, Capsule, Cylinder and Ellipsoid that build the corresponding ImplicitShape, apply orientation (as ZXZ Euler angles) and translation. Also updates TYPE_CHECKING imports to reference ImplicitShape. Useful for converting existing geometric shapes to implicit fields for sampling/boolean ops. --- microgen/shape/box.py | 10 + microgen/shape/capsule.py | 10 + microgen/shape/cylinder.py | 14 ++ microgen/shape/ellipsoid.py | 9 + microgen/shape/implicit_basic_factory.py | 249 +++++++++++++++++++++++ 5 files changed, 292 insertions(+) create mode 100644 microgen/shape/implicit_basic_factory.py diff --git a/microgen/shape/box.py b/microgen/shape/box.py index 53741cbb..8f63ee1b 100644 --- a/microgen/shape/box.py +++ b/microgen/shape/box.py @@ -19,6 +19,7 @@ if TYPE_CHECKING: from microgen.shape import KwargsGenerateType, Vector3DType + from microgen.shape.implicit_shape import ImplicitShape class Box(Shape): @@ -85,6 +86,15 @@ def generate_vtk( ) return rotate(box, self.center, self.orientation) + def to_implicit(self: Box) -> ImplicitShape: + """Convert this box to an :class:`ImplicitShape`.""" + from .implicit_basic_factory import implicit_box + + half = (self.dim[0] / 2, self.dim[1] / 2, self.dim[2] / 2) + shape = implicit_box(center=(0, 0, 0), half_extents=half) + angles = tuple(self.orientation.as_euler("ZXZ", degrees=True)) + return shape.rotate(angles).translate(tuple(self.center)) + def generateVtk(self: Box, **kwargs: KwargsGenerateType) -> pv.PolyData: # noqa: N802 """Deprecated. Use :meth:`generate_vtk` instead.""" # noqa: D401 return self.generate_vtk(**kwargs) diff --git a/microgen/shape/capsule.py b/microgen/shape/capsule.py index e4db4cb3..589a5760 100644 --- a/microgen/shape/capsule.py +++ b/microgen/shape/capsule.py @@ -18,6 +18,7 @@ if TYPE_CHECKING: from microgen.shape import KwargsGenerateType, Vector3DType + from microgen.shape.implicit_shape import ImplicitShape class Capsule(Shape): @@ -108,6 +109,15 @@ def generate_vtk( capsule = cylinder.boolean_union(sphere_left).boolean_union(sphere_right) return rotate(capsule, self.center, self.orientation) + def to_implicit(self: Capsule) -> ImplicitShape: + """Convert this capsule to an :class:`ImplicitShape`.""" + from .implicit_basic_factory import implicit_capsule + + h = self.height / 2.0 + shape = implicit_capsule(start=(-h, 0, 0), end=(h, 0, 0), radius=self.radius) + angles = tuple(self.orientation.as_euler("ZXZ", degrees=True)) + return shape.rotate(angles).translate(tuple(self.center)) + def generateVtk( # noqa: N802 self: Capsule, resolution: int = 100, diff --git a/microgen/shape/cylinder.py b/microgen/shape/cylinder.py index 07c26bcc..95de53d8 100644 --- a/microgen/shape/cylinder.py +++ b/microgen/shape/cylinder.py @@ -18,6 +18,7 @@ if TYPE_CHECKING: from microgen.shape import KwargsGenerateType, Vector3DType + from microgen.shape.implicit_shape import ImplicitShape class Cylinder(Shape): @@ -72,6 +73,19 @@ def generate_vtk( ) return rotate(cylinder, self.center, self.orientation) + def to_implicit(self: Cylinder) -> ImplicitShape: + """Convert this cylinder to an :class:`ImplicitShape`.""" + from .implicit_basic_factory import implicit_cylinder + + shape = implicit_cylinder( + center=(0, 0, 0), + axis=(1, 0, 0), + radius=self.radius, + height=self.height, + ) + angles = tuple(self.orientation.as_euler("ZXZ", degrees=True)) + return shape.rotate(angles).translate(tuple(self.center)) + def generateVtk( # noqa: N802 self: Cylinder, resolution: int = 100, diff --git a/microgen/shape/ellipsoid.py b/microgen/shape/ellipsoid.py index a902fd14..705e4102 100644 --- a/microgen/shape/ellipsoid.py +++ b/microgen/shape/ellipsoid.py @@ -20,6 +20,7 @@ if TYPE_CHECKING: from microgen.shape import KwargsGenerateType, Vector3DType + from microgen.shape.implicit_shape import ImplicitShape class Ellipsoid(Shape): @@ -89,6 +90,14 @@ def generate_vtk(self: Ellipsoid, **_: KwargsGenerateType) -> pv.PolyData: ellipsoid = sphere.transform(transform_matrix, inplace=False) return rotate(ellipsoid, self.center, self.orientation) + def to_implicit(self: Ellipsoid) -> ImplicitShape: + """Convert this ellipsoid to an :class:`ImplicitShape`.""" + from .implicit_basic_factory import implicit_ellipsoid + + shape = implicit_ellipsoid(center=(0, 0, 0), radii=self.radii) + angles = tuple(self.orientation.as_euler("ZXZ", degrees=True)) + return shape.rotate(angles).translate(tuple(self.center)) + def generateVtk(self: Ellipsoid, **_: KwargsGenerateType) -> pv.PolyData: # noqa: N802 """Deprecated. Use :meth:`generate_vtk` instead.""" # noqa: D401 return self.generate_vtk() diff --git a/microgen/shape/implicit_basic_factory.py b/microgen/shape/implicit_basic_factory.py new file mode 100644 index 00000000..db88ba42 --- /dev/null +++ b/microgen/shape/implicit_basic_factory.py @@ -0,0 +1,249 @@ +"""Implicit primitive factory functions. + +===================================================================== +Implicit Basic Factory (:mod:`microgen.shape.implicit_basic_factory`) +===================================================================== +""" + +from __future__ import annotations + +import numpy as np +import numpy.typing as npt + +from .implicit_shape import ImplicitShape + + +def implicit_sphere( + center: tuple[float, float, float] = (0, 0, 0), + radius: float = 1.0, +) -> ImplicitShape: + """Sphere SDF: ``|p - c| - r``.""" + cx, cy, cz = center + r = radius + margin = r * 1.1 + return ImplicitShape( + func=lambda x, y, z, _cx=cx, _cy=cy, _cz=cz, _r=r: np.sqrt( + (x - _cx) ** 2 + (y - _cy) ** 2 + (z - _cz) ** 2 + ) + - _r, + bounds=( + cx - margin, + cx + margin, + cy - margin, + cy + margin, + cz - margin, + cz + margin, + ), + ) + + +def implicit_box( + center: tuple[float, float, float] = (0, 0, 0), + half_extents: tuple[float, float, float] = (0.5, 0.5, 0.5), +) -> ImplicitShape: + """Axis-aligned box SDF (exact).""" + cx, cy, cz = center + hx, hy, hz = half_extents + margin = max(hx, hy, hz) * 0.1 + + def _box_sdf( + x: npt.NDArray[np.float64], + y: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], + ) -> npt.NDArray[np.float64]: + qx = np.abs(x - cx) - hx + qy = np.abs(y - cy) - hy + qz = np.abs(z - cz) - hz + outside = np.sqrt( + np.maximum(qx, 0.0) ** 2 + + np.maximum(qy, 0.0) ** 2 + + np.maximum(qz, 0.0) ** 2 + ) + inside = np.minimum(np.maximum(qx, np.maximum(qy, qz)), 0.0) + return outside + inside + + return ImplicitShape( + func=_box_sdf, + bounds=( + cx - hx - margin, + cx + hx + margin, + cy - hy - margin, + cy + hy + margin, + cz - hz - margin, + cz + hz + margin, + ), + ) + + +def implicit_cylinder( + center: tuple[float, float, float] = (0, 0, 0), + axis: tuple[float, float, float] = (0, 0, 1), + radius: float = 0.5, + height: float = 1.0, +) -> ImplicitShape: + """Capped cylinder SDF along arbitrary axis.""" + c = np.asarray(center, dtype=np.float64) + a = np.asarray(axis, dtype=np.float64) + a = a / np.linalg.norm(a) + half_h = height / 2.0 + + def _cyl_sdf( + x: npt.NDArray[np.float64], + y: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], + ) -> npt.NDArray[np.float64]: + px, py, pz = x - c[0], y - c[1], z - c[2] + proj = px * a[0] + py * a[1] + pz * a[2] + perp_x = px - proj * a[0] + perp_y = py - proj * a[1] + perp_z = pz - proj * a[2] + d_radial = np.sqrt(perp_x**2 + perp_y**2 + perp_z**2) - radius + d_axial = np.abs(proj) - half_h + outside = np.sqrt( + np.maximum(d_radial, 0.0) ** 2 + np.maximum(d_axial, 0.0) ** 2 + ) + inside = np.minimum(np.maximum(d_radial, d_axial), 0.0) + return outside + inside + + extent = max(radius, half_h) * 1.1 + np.linalg.norm(c) + return ImplicitShape( + func=_cyl_sdf, + bounds=( + c[0] - extent, + c[0] + extent, + c[1] - extent, + c[1] + extent, + c[2] - extent, + c[2] + extent, + ), + ) + + +def implicit_plane( + point: tuple[float, float, float] = (0, 0, 0), + normal: tuple[float, float, float] = (0, 0, 1), +) -> ImplicitShape: + """Half-space SDF: ``dot(p - point, normal)``.""" + pt = np.asarray(point, dtype=np.float64) + n = np.asarray(normal, dtype=np.float64) + n = n / np.linalg.norm(n) + + return ImplicitShape( + func=lambda x, y, z: (x - pt[0]) * n[0] + + (y - pt[1]) * n[1] + + (z - pt[2]) * n[2], + bounds=None, + ) + + +def implicit_torus( + center: tuple[float, float, float] = (0, 0, 0), + major_r: float = 1.0, + minor_r: float = 0.25, +) -> ImplicitShape: + """Torus SDF in XZ plane.""" + cx, cy, cz = center + + def _torus_sdf( + x: npt.NDArray[np.float64], + y: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], + ) -> npt.NDArray[np.float64]: + px, py, pz = x - cx, y - cy, z - cz + q_xz = np.sqrt(px**2 + pz**2) - major_r + return np.sqrt(q_xz**2 + py**2) - minor_r + + margin = (major_r + minor_r) * 1.1 + return ImplicitShape( + func=_torus_sdf, + bounds=( + cx - margin, + cx + margin, + cy - margin, + cy + margin, + cz - margin, + cz + margin, + ), + ) + + +def implicit_capsule( + start: tuple[float, float, float] = (0, 0, -0.5), + end: tuple[float, float, float] = (0, 0, 0.5), + radius: float = 0.25, +) -> ImplicitShape: + """Capsule SDF (line segment distance minus radius).""" + a = np.asarray(start, dtype=np.float64) + b = np.asarray(end, dtype=np.float64) + ab = b - a + ab_dot = float(np.dot(ab, ab)) + + if ab_dot < 1e-30: + # Degenerate capsule (start == end): fall back to sphere + return implicit_sphere(center=tuple(a), radius=radius) + + def _capsule_sdf( + x: npt.NDArray[np.float64], + y: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], + ) -> npt.NDArray[np.float64]: + px, py, pz = x - a[0], y - a[1], z - a[2] + t = np.clip((px * ab[0] + py * ab[1] + pz * ab[2]) / ab_dot, 0.0, 1.0) + dx = px - t * ab[0] + dy = py - t * ab[1] + dz = pz - t * ab[2] + return np.sqrt(dx**2 + dy**2 + dz**2) - radius + + mid = (a + b) / 2.0 + extent = np.linalg.norm(ab) / 2.0 + radius + margin = extent * 1.1 + return ImplicitShape( + func=_capsule_sdf, + bounds=( + mid[0] - margin, + mid[0] + margin, + mid[1] - margin, + mid[1] + margin, + mid[2] - margin, + mid[2] + margin, + ), + ) + + +def implicit_ellipsoid( + center: tuple[float, float, float] = (0, 0, 0), + radii: tuple[float, float, float] = (1.0, 0.5, 0.25), +) -> ImplicitShape: + """Ellipsoid implicit field: ``sqrt((x/rx)^2 + (y/ry)^2 + (z/rz)^2) - 1``. + + This is an approximate SDF but has the correct zero-level set. + """ + cx, cy, cz = center + rx, ry, rz = radii + margin = max(rx, ry, rz) * 1.1 + + def _ellipsoid_field( + x: npt.NDArray[np.float64], + y: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], + ) -> npt.NDArray[np.float64]: + return ( + np.sqrt( + ((x - cx) / rx) ** 2 + + ((y - cy) / ry) ** 2 + + ((z - cz) / rz) ** 2 + ) + - 1.0 + ) + + return ImplicitShape( + func=_ellipsoid_field, + bounds=( + cx - margin, + cx + margin, + cy - margin, + cy + margin, + cz - margin, + cz + margin, + ), + ) From 133c9fb2ac02f7e2827503dc23b316c034874034 Mon Sep 17 00:00:00 2001 From: Yves Chemisky Date: Tue, 24 Feb 2026 14:58:29 +0100 Subject: [PATCH 08/12] Create implicit_shape.py --- microgen/shape/implicit_shape.py | 556 +++++++++++++++++++++++++++++++ 1 file changed, 556 insertions(+) create mode 100644 microgen/shape/implicit_shape.py diff --git a/microgen/shape/implicit_shape.py b/microgen/shape/implicit_shape.py new file mode 100644 index 00000000..9ff341a8 --- /dev/null +++ b/microgen/shape/implicit_shape.py @@ -0,0 +1,556 @@ +"""F-rep Implicit Modeling. + +==================================================== +Implicit Shape (:mod:`microgen.shape.implicit_shape`) +==================================================== +""" + +from __future__ import annotations + +from typing import TYPE_CHECKING, Callable + +import numpy as np +import numpy.typing as npt +import pyvista as pv +from scipy.spatial.transform import Rotation + +from microgen.operations import rotate as rotate_mesh + +from .shape import Shape + +if TYPE_CHECKING: + from microgen.shape import KwargsGenerateType, Vector3DType + +Field = Callable[ + [npt.NDArray[np.float64], npt.NDArray[np.float64], npt.NDArray[np.float64]], + npt.NDArray[np.float64], +] + +BoundsType = tuple[float, float, float, float, float, float] + + +def _smooth_min( + a: npt.NDArray[np.float64], + b: npt.NDArray[np.float64], + k: float, +) -> npt.NDArray[np.float64]: + """Smooth minimum (Inigo Quilez cubic polynomial).""" + if k <= 0: + return np.minimum(a, b) + h = np.maximum(k - np.abs(a - b), 0.0) / k + return np.minimum(a, b) - h * h * h * k / 6.0 + + +def _smooth_max( + a: npt.NDArray[np.float64], + b: npt.NDArray[np.float64], + k: float, +) -> npt.NDArray[np.float64]: + """Smooth maximum.""" + return -_smooth_min(-a, -b, k) + + +def _merge_bounds( + a: BoundsType | None, + b: BoundsType | None, + mode: str = "union", +) -> BoundsType | None: + """Merge two bounding boxes.""" + if a is None and b is None: + return None + if a is None: + return b + if b is None: + return a + if mode == "union": + return ( + min(a[0], b[0]), + max(a[1], b[1]), + min(a[2], b[2]), + max(a[3], b[3]), + min(a[4], b[4]), + max(a[5], b[5]), + ) + # intersection + return ( + max(a[0], b[0]), + min(a[1], b[1]), + max(a[2], b[2]), + min(a[3], b[3]), + max(a[4], b[4]), + min(a[5], b[5]), + ) + + +class ImplicitShape(Shape): + """Implicit shape defined by a signed distance / scalar field. + + The convention is ``f(x, y, z) < 0`` means inside the shape. + + :param func: scalar field ``(x, y, z) -> array``, or ``None`` for subclasses + :param bounds: ``(xmin, xmax, ymin, ymax, zmin, zmax)`` or ``None`` + :param center: center of the shape + :param orientation: orientation of the shape + """ + + def __init__( + self: ImplicitShape, + func: Field | None = None, + bounds: BoundsType | None = None, + **kwargs: Vector3DType, + ) -> None: + """Initialize the implicit shape.""" + super().__init__(**kwargs) + self._func = func + self._bounds = bounds + + def evaluate( + self: ImplicitShape, + x: npt.NDArray[np.float64], + y: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], + ) -> npt.NDArray[np.float64]: + """Evaluate the scalar field at the given coordinates. + + :param x: x coordinates + :param y: y coordinates + :param z: z coordinates + :return: scalar field values (negative = inside) + """ + if self._func is None: + err_msg = "No scalar field function defined" + raise ValueError(err_msg) + return self._func(x, y, z) + + def generate_vtk( + self: ImplicitShape, + bounds: BoundsType | None = None, + resolution: int = 50, + **_: KwargsGenerateType, + ) -> pv.PolyData: + """Generate a VTK mesh of the implicit shape using the f < 0 convention. + + :param bounds: ``(xmin, xmax, ymin, ymax, zmin, zmax)`` + :param resolution: number of grid points per axis + :return: triangulated surface mesh + """ + bounds = bounds or self._bounds + if bounds is None: + err_msg = ( + "Bounds must be provided either at construction or in generate_vtk()" + ) + raise ValueError(err_msg) + + xmin, xmax, ymin, ymax, zmin, zmax = bounds + xi = np.linspace(xmin, xmax, resolution) + yi = np.linspace(ymin, ymax, resolution) + zi = np.linspace(zmin, zmax, resolution) + x, y, z = np.meshgrid(xi, yi, zi, indexing="ij") + + grid = pv.StructuredGrid(x, y, z) + field = self.evaluate( + x.ravel(order="F"), + y.ravel(order="F"), + z.ravel(order="F"), + ) + grid["implicit"] = field + + clipped = grid.clip_scalar(scalars="implicit", invert=True) + if clipped.n_cells == 0: + return pv.PolyData() + + polydata = clipped.extract_surface().clean().triangulate() + + polydata = rotate_mesh(polydata, center=(0, 0, 0), rotation=self.orientation) + return polydata.translate(xyz=self.center) + + def generate( + self: ImplicitShape, + bounds: BoundsType | None = None, + resolution: int = 50, + **_: KwargsGenerateType, + ) -> "cq.Shape": # noqa: F821 + """Generate a CadQuery Shape from the implicit surface. + + :param bounds: ``(xmin, xmax, ymin, ymax, zmin, zmax)`` + :param resolution: number of grid points per axis + :return: CadQuery Shape + """ + import cadquery as cq + + from .tpms import ShellCreationError + + mesh = self.generate_vtk(bounds=bounds, resolution=resolution) + if mesh.n_cells == 0: + err_msg = "Generated mesh is empty — check bounds and field function" + raise ValueError(err_msg) + + if not mesh.is_all_triangles: + mesh.triangulate(inplace=True) + triangles = mesh.faces.reshape(-1, 4)[:, 1:] + triangles = np.c_[triangles, triangles[:, 0]] + + faces = [] + for tri in triangles: + lines = [ + cq.Edge.makeLine( + cq.Vector(*mesh.points[start]), + cq.Vector(*mesh.points[end]), + ) + for start, end in zip(tri[:], tri[1:]) + ] + wire = cq.Wire.assembleEdges(lines) + faces.append(cq.Face.makeFromWires(wire)) + + try: + shell = cq.Shell.makeShell(faces) + except ValueError as err: + err_msg = ( + "Failed to create the shell, " + "try to increase the resolution or adjust bounds." + ) + raise ShellCreationError(err_msg) from err + + return cq.Shape(shell.wrapped) + + def generateVtk( # noqa: N802 + self: ImplicitShape, + **kwargs: KwargsGenerateType, + ) -> pv.PolyData: + """Deprecated. Use :meth:`generate_vtk` instead.""" # noqa: D401 + return self.generate_vtk(**kwargs) + + # --- Boolean operators --- + + def __or__(self: ImplicitShape, other: ImplicitShape) -> ImplicitShape: + """Union (``a | b``): inside where either field is negative.""" + return union(self, other) + + def __and__(self: ImplicitShape, other: ImplicitShape) -> ImplicitShape: + """Intersection (``a & b``): inside where both fields are negative.""" + return intersection(self, other) + + def __sub__(self: ImplicitShape, other: ImplicitShape) -> ImplicitShape: + """Difference (``a - b``): inside a but not b.""" + return difference(self, other) + + def __invert__(self: ImplicitShape) -> ImplicitShape: + """Complement (``~a``): negate the field.""" + f = self._func + return ImplicitShape( + func=lambda x, y, z, _f=f: -_f(x, y, z), + bounds=self._bounds, + ) + + # --- Smooth booleans --- + + def smooth_union( + self: ImplicitShape, + other: ImplicitShape, + k: float, + ) -> ImplicitShape: + """Smooth union with blending radius *k*.""" + return smooth_union(self, other, k) + + def smooth_intersection( + self: ImplicitShape, + other: ImplicitShape, + k: float, + ) -> ImplicitShape: + """Smooth intersection with blending radius *k*.""" + return smooth_intersection(self, other, k) + + def smooth_difference( + self: ImplicitShape, + other: ImplicitShape, + k: float, + ) -> ImplicitShape: + """Smooth difference with blending radius *k*.""" + return smooth_difference(self, other, k) + + # --- Transforms --- + + def translate( # type: ignore[override] + self: ImplicitShape, + offset: tuple[float, float, float], + ) -> ImplicitShape: + """Return a new shape translated by *offset*.""" + dx, dy, dz = offset + f = self._func + new_bounds = None + if self._bounds is not None: + b = self._bounds + new_bounds = ( + b[0] + dx, + b[1] + dx, + b[2] + dy, + b[3] + dy, + b[4] + dz, + b[5] + dz, + ) + return ImplicitShape( + func=lambda x, y, z, _f=f, _dx=dx, _dy=dy, _dz=dz: _f( + x - _dx, y - _dy, z - _dz + ), + bounds=new_bounds, + ) + + def rotate( # type: ignore[override] + self: ImplicitShape, + angles: tuple[float, float, float], + convention: str = "ZXZ", + ) -> ImplicitShape: + """Return a new shape rotated by Euler *angles* (degrees).""" + rot = Rotation.from_euler(convention, angles, degrees=True) + inv_matrix = rot.inv().as_matrix() + f = self._func + return ImplicitShape( + func=lambda x, y, z, _f=f, _m=inv_matrix: _f( + *(_m @ np.array([x, y, z])) + ), + bounds=self._bounds, # conservative: keep original bounds + ) + + def scale( + self: ImplicitShape, + factor: float, + ) -> ImplicitShape: + """Return a new shape uniformly scaled by *factor*.""" + f = self._func + new_bounds = None + if self._bounds is not None: + b = self._bounds + new_bounds = ( + b[0] * factor, + b[1] * factor, + b[2] * factor, + b[3] * factor, + b[4] * factor, + b[5] * factor, + ) + return ImplicitShape( + func=lambda x, y, z, _f=f, _s=factor: _f(x / _s, y / _s, z / _s) * _s, + bounds=new_bounds, + ) + + +# --------------------------------------------------------------------------- +# Module-level boolean operations +# --------------------------------------------------------------------------- + + +def union(a: ImplicitShape, b: ImplicitShape) -> ImplicitShape: + """Union of two implicit shapes (hard boolean).""" + fa, fb = a._func, b._func # noqa: SLF001 + return ImplicitShape( + func=lambda x, y, z, _fa=fa, _fb=fb: np.minimum(_fa(x, y, z), _fb(x, y, z)), + bounds=_merge_bounds(a._bounds, b._bounds, "union"), # noqa: SLF001 + ) + + +def intersection(a: ImplicitShape, b: ImplicitShape) -> ImplicitShape: + """Intersection of two implicit shapes (hard boolean).""" + fa, fb = a._func, b._func # noqa: SLF001 + return ImplicitShape( + func=lambda x, y, z, _fa=fa, _fb=fb: np.maximum(_fa(x, y, z), _fb(x, y, z)), + bounds=_merge_bounds(a._bounds, b._bounds, "intersection"), # noqa: SLF001 + ) + + +def difference(a: ImplicitShape, b: ImplicitShape) -> ImplicitShape: + """Difference of two implicit shapes (a minus b).""" + fa, fb = a._func, b._func # noqa: SLF001 + return ImplicitShape( + func=lambda x, y, z, _fa=fa, _fb=fb: np.maximum( + _fa(x, y, z), -_fb(x, y, z) + ), + bounds=a._bounds, # noqa: SLF001 + ) + + +def smooth_union( + a: ImplicitShape, + b: ImplicitShape, + k: float, +) -> ImplicitShape: + """Smooth union with blending radius *k*.""" + fa, fb = a._func, b._func # noqa: SLF001 + return ImplicitShape( + func=lambda x, y, z, _fa=fa, _fb=fb, _k=k: _smooth_min( + _fa(x, y, z), _fb(x, y, z), _k + ), + bounds=_merge_bounds(a._bounds, b._bounds, "union"), # noqa: SLF001 + ) + + +def smooth_intersection( + a: ImplicitShape, + b: ImplicitShape, + k: float, +) -> ImplicitShape: + """Smooth intersection with blending radius *k*.""" + fa, fb = a._func, b._func # noqa: SLF001 + return ImplicitShape( + func=lambda x, y, z, _fa=fa, _fb=fb, _k=k: _smooth_max( + _fa(x, y, z), _fb(x, y, z), _k + ), + bounds=_merge_bounds(a._bounds, b._bounds, "intersection"), # noqa: SLF001 + ) + + +def smooth_difference( + a: ImplicitShape, + b: ImplicitShape, + k: float, +) -> ImplicitShape: + """Smooth difference (a minus b) with blending radius *k*.""" + fa, fb = a._func, b._func # noqa: SLF001 + return ImplicitShape( + func=lambda x, y, z, _fa=fa, _fb=fb, _k=k: _smooth_max( + _fa(x, y, z), -_fb(x, y, z), _k + ), + bounds=a._bounds, # noqa: SLF001 + ) + + +def batch_smooth_union( + shapes: list[ImplicitShape], + k: float = 0.0, +) -> ImplicitShape: + """Combine many shapes with smooth union in a flat loop (no recursion). + + This avoids the recursion-depth limit that arises when chaining hundreds + of binary ``smooth_union`` calls, each wrapping the previous in a lambda. + """ + if not shapes: + msg = "batch_smooth_union requires at least one shape" + raise ValueError(msg) + + funcs = [s._func for s in shapes] # noqa: SLF001 + + # Merge all bounds + merged = shapes[0]._bounds # noqa: SLF001 + for s in shapes[1:]: + merged = _merge_bounds(merged, s._bounds, "union") # noqa: SLF001 + + def _batched( + x: npt.NDArray[np.float64], + y: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], + ) -> npt.NDArray[np.float64]: + result = funcs[0](x, y, z) + if k > 0: + for fn in funcs[1:]: + result = _smooth_min(result, fn(x, y, z), k) + else: + for fn in funcs[1:]: + result = np.minimum(result, fn(x, y, z)) + return result + + return ImplicitShape(func=_batched, bounds=merged) + + +# --------------------------------------------------------------------------- +# Primitive factory re-exports (moved to implicit_basic_factory.py) +# --------------------------------------------------------------------------- +from .implicit_basic_factory import ( # noqa: E402 + implicit_box, + implicit_capsule, + implicit_cylinder, + implicit_ellipsoid, + implicit_plane, + implicit_sphere, + implicit_torus, +) + + +# --------------------------------------------------------------------------- +# Utility operations +# --------------------------------------------------------------------------- + + +def shell(shape: ImplicitShape, thickness: float) -> ImplicitShape: + """Hollow shell: ``|f(p)| - thickness / 2``.""" + f = shape._func # noqa: SLF001 + half_t = thickness / 2.0 + return ImplicitShape( + func=lambda x, y, z, _f=f, _ht=half_t: np.abs(_f(x, y, z)) - _ht, + bounds=shape._bounds, # noqa: SLF001 + ) + + +def repeat( + shape: ImplicitShape, + spacing: tuple[float, float, float], + k: float = 0.0, +) -> ImplicitShape: + """Infinite repetition via coordinate modulo. + + :param shape: unit cell shape to tile + :param spacing: ``(sx, sy, sz)`` repetition period per axis + :param k: smooth blending radius across cell boundaries. + When ``k > 0``, the base field is evaluated at the 26 neighboring + periodic images in addition to the current cell and all values + are combined with smooth minimum, so that primitives from adjacent + cells blend seamlessly. When ``k <= 0`` (default), a simple + coordinate-modulo repetition is used (hard tiling). + """ + sx, sy, sz = spacing + f = shape._func # noqa: SLF001 + + if k <= 0: + + def _repeated( + x: npt.NDArray[np.float64], + y: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], + ) -> npt.NDArray[np.float64]: + rx = np.mod(x + sx / 2, sx) - sx / 2 + ry = np.mod(y + sy / 2, sy) - sy / 2 + rz = np.mod(z + sz / 2, sz) - sz / 2 + return f(rx, ry, rz) + + else: + offsets = [ + (dx * sx, dy * sy, dz * sz) + for dx in (-1, 0, 1) + for dy in (-1, 0, 1) + for dz in (-1, 0, 1) + ] + + def _repeated( + x: npt.NDArray[np.float64], + y: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], + ) -> npt.NDArray[np.float64]: + # Wrap to the central cell + cx = np.mod(x + sx / 2, sx) - sx / 2 + cy = np.mod(y + sy / 2, sy) - sy / 2 + cz = np.mod(z + sz / 2, sz) - sz / 2 + # Evaluate at all 27 images, accumulate with smooth min + result = f(cx + offsets[0][0], cy + offsets[0][1], cz + offsets[0][2]) + for ox, oy, oz in offsets[1:]: + result = _smooth_min(result, f(cx + ox, cy + oy, cz + oz), k) + return result + + return ImplicitShape(func=_repeated, bounds=None) + + +def blend( + a: ImplicitShape, + b: ImplicitShape, + factor: float = 0.5, +) -> ImplicitShape: + """Linear interpolation between two fields: ``(1-t)*a + t*b``.""" + fa, fb = a._func, b._func # noqa: SLF001 + t = factor + return ImplicitShape( + func=lambda x, y, z, _fa=fa, _fb=fb, _t=t: (1.0 - _t) * _fa(x, y, z) + + _t * _fb(x, y, z), + bounds=_merge_bounds(a._bounds, b._bounds, "union"), # noqa: SLF001 + ) + + +def from_field(func: Field, bounds: BoundsType | None = None) -> ImplicitShape: + """Wrap any callable ``f(x, y, z) -> scalar`` as an ImplicitShape.""" + return ImplicitShape(func=func, bounds=bounds) From 3a68827443d7eda1a7283d294fd769db3b9bd453 Mon Sep 17 00:00:00 2001 From: Yves Chemisky Date: Tue, 24 Feb 2026 14:58:36 +0100 Subject: [PATCH 09/12] Update sphere.py --- microgen/shape/sphere.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/microgen/shape/sphere.py b/microgen/shape/sphere.py index e41c94bd..947ac555 100644 --- a/microgen/shape/sphere.py +++ b/microgen/shape/sphere.py @@ -17,6 +17,7 @@ if TYPE_CHECKING: from microgen.shape import KwargsGenerateType, Vector3DType + from microgen.shape.implicit_shape import ImplicitShape class Sphere(Shape): @@ -68,6 +69,14 @@ def generate_vtk( phi_resolution=phi_resolution, ) + def to_implicit(self: Sphere) -> ImplicitShape: + """Convert this sphere to an :class:`ImplicitShape`.""" + from .implicit_basic_factory import implicit_sphere + + shape = implicit_sphere(center=(0, 0, 0), radius=self.radius) + angles = tuple(self.orientation.as_euler("ZXZ", degrees=True)) + return shape.rotate(angles).translate(tuple(self.center)) + def generateVtk( # noqa: N802 self: Sphere, theta_resolution: int = 50, From 1203e0fa9edda33fdac0d7d4f8e63e64dc3af09e Mon Sep 17 00:00:00 2001 From: Yves Chemisky Date: Tue, 24 Feb 2026 14:59:25 +0100 Subject: [PATCH 10/12] Tpms: inherit ImplicitShape and set implicit field Change Tpms to inherit from ImplicitShape (import implicit_shape.ImplicitShape) and initialize implicit-interface attributes for compatibility. The constructor now sets self._func as a lambda wrapping surface_function with k_x/k_y/k_z scaling and phase shifts, and computes self._bounds from cell_size and repeat_cell. This enables Tpms instances to be used where an implicit field/bounds are expected. --- microgen/shape/tpms.py | 21 +++++++++++++++++++-- 1 file changed, 19 insertions(+), 2 deletions(-) diff --git a/microgen/shape/tpms.py b/microgen/shape/tpms.py index 0f2e2251..69722e2d 100644 --- a/microgen/shape/tpms.py +++ b/microgen/shape/tpms.py @@ -25,7 +25,7 @@ from microgen.operations import fuseShapes, rotate -from .shape import Shape +from .implicit_shape import ImplicitShape if TYPE_CHECKING: from microgen.shape import KwargsGenerateType, TpmsPartType, Vector3DType @@ -40,7 +40,7 @@ _DIM = 3 -class Tpms(Shape): +class Tpms(ImplicitShape): """Triply Periodical Minimal Surfaces. Class to generate Triply Periodical Minimal Surfaces (TPMS) @@ -149,6 +149,23 @@ def __init__( # noqa: PLR0913 self.offset = offset # call setter self.offset_updated: bool + # Set implicit field for ImplicitShape compatibility + k_x, k_y, k_z = 2.0 * np.pi / self.cell_size + self._func = lambda x, y, z: self.surface_function( + k_x * (x + self.phase_shift[0]), + k_y * (y + self.phase_shift[1]), + k_z * (z + self.phase_shift[2]), + ) + half_dim = 0.5 * self.cell_size * self.repeat_cell + self._bounds = ( + -half_dim[0], + half_dim[0], + -half_dim[1], + half_dim[1], + -half_dim[2], + half_dim[2], + ) + @classmethod def offset_from_density( cls: type[Tpms], From 6700813004722ac22864edf4e5467bf6ec8834f6 Mon Sep 17 00:00:00 2001 From: Yves Chemisky Date: Tue, 24 Feb 2026 15:07:21 +0100 Subject: [PATCH 11/12] Create test_implicit_shape.py --- tests/shapes/test_implicit_shape.py | 486 ++++++++++++++++++++++++++++ 1 file changed, 486 insertions(+) create mode 100644 tests/shapes/test_implicit_shape.py diff --git a/tests/shapes/test_implicit_shape.py b/tests/shapes/test_implicit_shape.py new file mode 100644 index 00000000..3bdfa413 --- /dev/null +++ b/tests/shapes/test_implicit_shape.py @@ -0,0 +1,486 @@ +"""Tests for the F-rep implicit modeling module.""" + +from __future__ import annotations + +import numpy as np +import pytest +import pyvista as pv + +from microgen.shape import implicit_shape, surface_functions, Tpms +from microgen.shape.implicit_shape import ( + ImplicitShape, + blend, + difference, + from_field, + implicit_box, + implicit_capsule, + implicit_cylinder, + implicit_plane, + implicit_sphere, + implicit_torus, + intersection, + repeat, + shell, + smooth_difference, + smooth_intersection, + smooth_union, + union, +) + + +# --------------------------------------------------------------------------- +# Primitive evaluation tests +# --------------------------------------------------------------------------- + + +class TestPrimitives: + """Test that primitives evaluate correctly at known points.""" + + def test_sphere_inside(self): + s = implicit_sphere(center=(0, 0, 0), radius=1.0) + x = np.array([0.0]) + y = np.array([0.0]) + z = np.array([0.0]) + assert s.evaluate(x, y, z)[0] < 0 + + def test_sphere_surface(self): + s = implicit_sphere(center=(0, 0, 0), radius=1.0) + x = np.array([1.0]) + y = np.array([0.0]) + z = np.array([0.0]) + assert abs(s.evaluate(x, y, z)[0]) < 1e-10 + + def test_sphere_outside(self): + s = implicit_sphere(center=(0, 0, 0), radius=1.0) + x = np.array([2.0]) + y = np.array([0.0]) + z = np.array([0.0]) + assert s.evaluate(x, y, z)[0] > 0 + + def test_box_inside(self): + b = implicit_box(center=(0, 0, 0), half_extents=(1, 1, 1)) + assert b.evaluate(np.array([0.0]), np.array([0.0]), np.array([0.0]))[0] < 0 + + def test_box_outside(self): + b = implicit_box(center=(0, 0, 0), half_extents=(1, 1, 1)) + assert b.evaluate(np.array([2.0]), np.array([0.0]), np.array([0.0]))[0] > 0 + + def test_cylinder_inside(self): + c = implicit_cylinder(center=(0, 0, 0), axis=(0, 0, 1), radius=1, height=2) + assert c.evaluate(np.array([0.0]), np.array([0.0]), np.array([0.0]))[0] < 0 + + def test_cylinder_outside(self): + c = implicit_cylinder(center=(0, 0, 0), axis=(0, 0, 1), radius=1, height=2) + assert c.evaluate(np.array([2.0]), np.array([0.0]), np.array([0.0]))[0] > 0 + + def test_plane(self): + p = implicit_plane(point=(0, 0, 0), normal=(0, 0, 1)) + assert p.evaluate(np.array([0.0]), np.array([0.0]), np.array([-1.0]))[0] < 0 + assert p.evaluate(np.array([0.0]), np.array([0.0]), np.array([1.0]))[0] > 0 + + def test_plane_bounds_none(self): + p = implicit_plane() + assert p._bounds is None + + def test_torus_inside(self): + t = implicit_torus(center=(0, 0, 0), major_r=1.0, minor_r=0.25) + # Point on the ring center at (1, 0, 0) should be inside + assert t.evaluate(np.array([1.0]), np.array([0.0]), np.array([0.0]))[0] < 0 + + def test_torus_outside(self): + t = implicit_torus(center=(0, 0, 0), major_r=1.0, minor_r=0.25) + assert t.evaluate(np.array([0.0]), np.array([0.0]), np.array([0.0]))[0] > 0 + + def test_capsule_inside(self): + c = implicit_capsule(start=(0, 0, -1), end=(0, 0, 1), radius=0.5) + assert c.evaluate(np.array([0.0]), np.array([0.0]), np.array([0.0]))[0] < 0 + + def test_capsule_outside(self): + c = implicit_capsule(start=(0, 0, -1), end=(0, 0, 1), radius=0.5) + assert c.evaluate(np.array([2.0]), np.array([0.0]), np.array([0.0]))[0] > 0 + + +# --------------------------------------------------------------------------- +# Boolean operations +# --------------------------------------------------------------------------- + + +class TestBooleans: + """Test hard boolean operations.""" + + def test_union_min(self): + s1 = implicit_sphere(center=(-0.5, 0, 0), radius=1.0) + s2 = implicit_sphere(center=(0.5, 0, 0), radius=1.0) + u = union(s1, s2) + x = np.array([0.0]) + y = np.array([0.0]) + z = np.array([0.0]) + expected = min(s1.evaluate(x, y, z)[0], s2.evaluate(x, y, z)[0]) + assert u.evaluate(x, y, z)[0] == pytest.approx(expected) + + def test_intersection_max(self): + s1 = implicit_sphere(center=(-0.5, 0, 0), radius=1.0) + s2 = implicit_sphere(center=(0.5, 0, 0), radius=1.0) + i = intersection(s1, s2) + x = np.array([0.0]) + y = np.array([0.0]) + z = np.array([0.0]) + expected = max(s1.evaluate(x, y, z)[0], s2.evaluate(x, y, z)[0]) + assert i.evaluate(x, y, z)[0] == pytest.approx(expected) + + def test_difference(self): + s1 = implicit_sphere(center=(0, 0, 0), radius=1.0) + s2 = implicit_sphere(center=(0.5, 0, 0), radius=0.5) + d = difference(s1, s2) + # Inside s1 but outside s2 + x = np.array([-0.5]) + y = np.array([0.0]) + z = np.array([0.0]) + assert d.evaluate(x, y, z)[0] < 0 + + def test_operators(self): + s1 = implicit_sphere(radius=1.0) + s2 = implicit_box(half_extents=(0.5, 0.5, 0.5)) + x = np.array([0.0]) + y = np.array([0.0]) + z = np.array([0.0]) + + u = s1 | s2 + assert u.evaluate(x, y, z)[0] < 0 + + i = s1 & s2 + assert i.evaluate(x, y, z)[0] < 0 + + d = s1 - s2 + # At origin, both are inside, so difference should be max(f_s1, -f_s2) + val = d.evaluate(x, y, z)[0] + assert isinstance(val, (float, np.floating)) + + def test_complement(self): + s = implicit_sphere(radius=1.0) + c = ~s + x = np.array([0.0]) + y = np.array([0.0]) + z = np.array([0.0]) + assert c.evaluate(x, y, z)[0] > 0 # origin was inside, complement is outside + + +# --------------------------------------------------------------------------- +# Smooth booleans +# --------------------------------------------------------------------------- + + +class TestSmoothBooleans: + """Test smooth boolean operations.""" + + def test_smooth_union_blending_zone(self): + s1 = implicit_sphere(center=(-0.5, 0, 0), radius=1.0) + s2 = implicit_sphere(center=(0.5, 0, 0), radius=1.0) + su = smooth_union(s1, s2, k=0.5) + hard = union(s1, s2) + x = np.array([0.0]) + y = np.array([0.0]) + z = np.array([0.0]) + # Smooth union should be <= hard union (more material in blend zone) + assert su.evaluate(x, y, z)[0] <= hard.evaluate(x, y, z)[0] + 1e-10 + + def test_smooth_union_k_zero_equals_hard(self): + s1 = implicit_sphere(center=(-0.5, 0, 0), radius=1.0) + s2 = implicit_sphere(center=(0.5, 0, 0), radius=1.0) + su = smooth_union(s1, s2, k=0.0) + hard = union(s1, s2) + x = np.linspace(-2, 2, 10) + y = np.zeros(10) + z = np.zeros(10) + np.testing.assert_allclose( + su.evaluate(x, y, z), hard.evaluate(x, y, z), atol=1e-10 + ) + + def test_smooth_intersection(self): + s1 = implicit_sphere(radius=1.0) + s2 = implicit_box(half_extents=(0.5, 0.5, 0.5)) + si = smooth_intersection(s1, s2, k=0.3) + x = np.array([0.0]) + y = np.array([0.0]) + z = np.array([0.0]) + assert si.evaluate(x, y, z)[0] < 0 # still inside + + def test_smooth_difference(self): + s1 = implicit_sphere(radius=1.0) + s2 = implicit_sphere(center=(1.0, 0, 0), radius=0.5) + sd = smooth_difference(s1, s2, k=0.2) + x = np.array([-0.5]) + y = np.array([0.0]) + z = np.array([0.0]) + assert sd.evaluate(x, y, z)[0] < 0 # far from subtraction + + def test_smooth_methods(self): + s1 = implicit_sphere(radius=1.0) + s2 = implicit_box(half_extents=(0.5, 0.5, 0.5)) + su = s1.smooth_union(s2, k=0.3) + si = s1.smooth_intersection(s2, k=0.3) + sd = s1.smooth_difference(s2, k=0.3) + x = np.array([0.0]) + y = np.array([0.0]) + z = np.array([0.0]) + assert isinstance(su.evaluate(x, y, z)[0], (float, np.floating)) + assert isinstance(si.evaluate(x, y, z)[0], (float, np.floating)) + assert isinstance(sd.evaluate(x, y, z)[0], (float, np.floating)) + + +# --------------------------------------------------------------------------- +# Transforms +# --------------------------------------------------------------------------- + + +class TestTransforms: + """Test transform operations.""" + + def test_translate(self): + s = implicit_sphere(center=(0, 0, 0), radius=1.0) + st = s.translate((2, 0, 0)) + x = np.array([2.0]) + y = np.array([0.0]) + z = np.array([0.0]) + assert st.evaluate(x, y, z)[0] < 0 # center moved to (2,0,0) + + x0 = np.array([0.0]) + assert st.evaluate(x0, y, z)[0] > 0 # origin is now outside + + def test_rotate_90(self): + # Cylinder along z, rotate 90 degrees around y -> cylinder along x + c = implicit_cylinder(radius=0.5, height=2.0) + cr = c.rotate((0, 90, 0), convention="xyz") + # Point along x axis (was z) should be inside + assert cr.evaluate(np.array([0.5]), np.array([0.0]), np.array([0.0]))[0] < 0 + + def test_scale(self): + s = implicit_sphere(radius=1.0) + ss = s.scale(2.0) + # At distance 1.5, original sphere is outside but scaled should be inside + x = np.array([1.5]) + y = np.array([0.0]) + z = np.array([0.0]) + assert s.evaluate(x, y, z)[0] > 0 + assert ss.evaluate(x, y, z)[0] < 0 + + def test_translate_bounds(self): + s = implicit_sphere(radius=1.0) + st = s.translate((5, 0, 0)) + assert st._bounds is not None + assert st._bounds[0] > 3.0 # xmin shifted + + def test_scale_bounds(self): + s = implicit_sphere(radius=1.0) + ss = s.scale(3.0) + assert ss._bounds is not None + assert ss._bounds[1] > 3.0 # xmax scaled + + +# --------------------------------------------------------------------------- +# generate_vtk +# --------------------------------------------------------------------------- + + +class TestGenerateVtk: + """Test mesh generation.""" + + def test_sphere_mesh(self): + s = implicit_sphere(radius=1.0) + mesh = s.generate_vtk(resolution=30) + assert isinstance(mesh, pv.PolyData) + assert mesh.n_cells > 0 + + def test_box_mesh(self): + b = implicit_box(half_extents=(0.5, 0.5, 0.5)) + mesh = b.generate_vtk(resolution=30) + assert isinstance(mesh, pv.PolyData) + assert mesh.n_cells > 0 + + def test_boolean_mesh(self): + s1 = implicit_sphere(radius=1.0) + s2 = implicit_box(half_extents=(0.6, 0.6, 0.6)) + result = s1 & s2 + mesh = result.generate_vtk(resolution=30) + assert isinstance(mesh, pv.PolyData) + assert mesh.n_cells > 0 + + def test_no_bounds_raises(self): + shape = ImplicitShape(func=lambda x, y, z: x**2 + y**2 + z**2 - 1) + with pytest.raises(ValueError, match="Bounds must be provided"): + shape.generate_vtk() + + def test_explicit_bounds_override(self): + s = implicit_sphere(radius=1.0) + mesh = s.generate_vtk(bounds=(-2, 2, -2, 2, -2, 2), resolution=30) + assert mesh.n_cells > 0 + + def test_generateVtk_deprecated(self): + s = implicit_sphere(radius=1.0) + mesh = s.generateVtk(resolution=20) + assert isinstance(mesh, pv.PolyData) + + +# --------------------------------------------------------------------------- +# Bounds propagation +# --------------------------------------------------------------------------- + + +class TestBoundsPropagation: + """Test bounds merge behavior.""" + + def test_union_expands(self): + s1 = implicit_sphere(center=(-2, 0, 0), radius=1.0) + s2 = implicit_sphere(center=(2, 0, 0), radius=1.0) + u = s1 | s2 + assert u._bounds[0] <= s1._bounds[0] # xmin from s1 + assert u._bounds[0] <= s2._bounds[0] + assert u._bounds[1] >= s1._bounds[1] # xmax includes both + assert u._bounds[1] >= s2._bounds[1] + + def test_intersection_shrinks(self): + s1 = implicit_sphere(center=(-0.5, 0, 0), radius=1.5) + s2 = implicit_sphere(center=(0.5, 0, 0), radius=1.5) + i = s1 & s2 + assert i._bounds[0] >= s1._bounds[0] + assert i._bounds[1] <= s2._bounds[1] + + def test_plane_none_bounds(self): + p = implicit_plane() + s = implicit_sphere(radius=1.0) + u = p | s + # When one is None, result takes the other + assert u._bounds == s._bounds + + +# --------------------------------------------------------------------------- +# Utilities +# --------------------------------------------------------------------------- + + +class TestUtilities: + """Test utility functions.""" + + def test_shell(self): + s = implicit_sphere(radius=1.0) + sh = shell(s, thickness=0.2) + # At radius=1.0 (surface), |f|=0, shell should be inside (-thickness/2) + x = np.array([1.0]) + y = np.array([0.0]) + z = np.array([0.0]) + assert sh.evaluate(x, y, z)[0] < 0 + # At origin, |f|=1.0 >> 0.1, shell should be outside + x0 = np.array([0.0]) + assert sh.evaluate(x0, y, z)[0] > 0 + + def test_repeat(self): + s = implicit_sphere(radius=0.3) + r = repeat(s, spacing=(1.0, 1.0, 1.0)) + # At (1,0,0) should be inside a repeated copy + x = np.array([1.0]) + y = np.array([0.0]) + z = np.array([0.0]) + assert r.evaluate(x, y, z)[0] < 0 + assert r._bounds is None # infinite repetition + + def test_blend(self): + s1 = implicit_sphere(radius=1.0) + s2 = implicit_box(half_extents=(0.5, 0.5, 0.5)) + b = blend(s1, s2, factor=0.5) + x = np.array([0.0]) + y = np.array([0.0]) + z = np.array([0.0]) + expected = 0.5 * s1.evaluate(x, y, z)[0] + 0.5 * s2.evaluate(x, y, z)[0] + assert b.evaluate(x, y, z)[0] == pytest.approx(expected) + + def test_from_field(self): + shape = from_field( + func=lambda x, y, z: x**2 + y**2 + z**2 - 1, + bounds=(-2, 2, -2, 2, -2, 2), + ) + assert shape.evaluate(np.array([0.0]), np.array([0.0]), np.array([0.0]))[0] < 0 + assert shape._bounds == (-2, 2, -2, 2, -2, 2) + + +# --------------------------------------------------------------------------- +# TPMS integration +# --------------------------------------------------------------------------- + + +class TestTpmsIntegration: + """Test that TPMS still works and integrates with ImplicitShape.""" + + def test_isinstance(self): + tpms = Tpms( + surface_function=surface_functions.gyroid, + offset=0.3, + ) + assert isinstance(tpms, ImplicitShape) + + def test_tpms_has_func(self): + tpms = Tpms( + surface_function=surface_functions.gyroid, + offset=0.3, + ) + assert tpms._func is not None + + def test_tpms_has_bounds(self): + tpms = Tpms( + surface_function=surface_functions.gyroid, + offset=0.3, + ) + assert tpms._bounds is not None + + def test_tpms_evaluate(self): + tpms = Tpms( + surface_function=surface_functions.gyroid, + offset=0.3, + ) + x = np.array([0.0]) + y = np.array([0.0]) + z = np.array([0.0]) + val = tpms.evaluate(x, y, z) + assert isinstance(val[0], (float, np.floating)) + + def test_tpms_boolean_with_sphere(self): + tpms = Tpms( + surface_function=surface_functions.gyroid, + offset=0.3, + ) + sphere = implicit_sphere(radius=0.4) + result = tpms & sphere + assert isinstance(result, ImplicitShape) + mesh = result.generate_vtk( + bounds=(-0.5, 0.5, -0.5, 0.5, -0.5, 0.5), + resolution=30, + ) + assert isinstance(mesh, pv.PolyData) + + def test_tpms_generate_vtk_unchanged(self): + """TPMS generate_vtk still works with its own override.""" + tpms = Tpms( + surface_function=surface_functions.gyroid, + offset=0.3, + ) + mesh = tpms.generate_vtk(type_part="sheet") + assert isinstance(mesh, pv.PolyData) + assert mesh.n_cells > 0 + + +# --------------------------------------------------------------------------- +# Error handling +# --------------------------------------------------------------------------- + + +class TestErrorHandling: + """Test error conditions.""" + + def test_no_func_evaluate_raises(self): + shape = ImplicitShape() + with pytest.raises(ValueError, match="No scalar field"): + shape.evaluate(np.array([0.0]), np.array([0.0]), np.array([0.0])) + + def test_no_bounds_generate_vtk_raises(self): + shape = ImplicitShape(func=lambda x, y, z: x) + with pytest.raises(ValueError, match="Bounds must be provided"): + shape.generate_vtk() From 08bc8459c2ab3f85a695f521237c2f4a15308675 Mon Sep 17 00:00:00 2001 From: Yves Chemisky Date: Tue, 24 Feb 2026 23:38:03 +0100 Subject: [PATCH 12/12] revert separate factory --- docs/microgen.shape.rst | 5 - docs/tutorials.rst | 2 +- microgen/__init__.py | 2 - microgen/shape/__init__.py | 3 +- microgen/shape/box.py | 2 +- microgen/shape/capsule.py | 2 +- microgen/shape/cylinder.py | 2 +- microgen/shape/ellipsoid.py | 2 +- microgen/shape/implicit_basic_factory.py | 249 ----------------------- microgen/shape/implicit_shape.py | 247 +++++++++++++++++++++- microgen/shape/sphere.py | 2 +- 11 files changed, 244 insertions(+), 274 deletions(-) delete mode 100644 microgen/shape/implicit_basic_factory.py diff --git a/docs/microgen.shape.rst b/docs/microgen.shape.rst index b2ea6118..e4dbf9cf 100644 --- a/docs/microgen.shape.rst +++ b/docs/microgen.shape.rst @@ -57,11 +57,6 @@ microgen.shape package :undoc-members: :show-inheritance: -.. automodule:: microgen.shape.implicit_basic_factory - :members: - :undoc-members: - :show-inheritance: - .. automodule:: microgen.shape.surface_functions :members: :undoc-members: diff --git a/docs/tutorials.rst b/docs/tutorials.rst index 426cc640..500f05c5 100644 --- a/docs/tutorials.rst +++ b/docs/tutorials.rst @@ -278,7 +278,7 @@ combine parametric geometry with F-rep boolean operations: .. jupyter-execute:: import microgen - from microgen.shape.implicit_basic_factory import implicit_box + from microgen.shape.implicit_shape import implicit_box sphere = microgen.Sphere(center=(0.3, 0, 0), radius=0.4) sphere_imp = sphere.to_implicit() diff --git a/microgen/__init__.py b/microgen/__init__.py index 9bc91eae..76ff6ec6 100644 --- a/microgen/__init__.py +++ b/microgen/__init__.py @@ -62,7 +62,6 @@ TruncatedCube, TruncatedCuboctahedron, TruncatedOctahedron, - implicit_basic_factory, implicit_shape, new_geometry, newGeometry, @@ -106,7 +105,6 @@ "CylindricalTpms", "SphericalTpms", "check_if_only_linear_tetrahedral", - "implicit_basic_factory", "implicit_shape", "batch_smooth_union", "cut_phase_by_shape_list", diff --git a/microgen/shape/__init__.py b/microgen/shape/__init__.py index de39e1d5..67ce83ac 100644 --- a/microgen/shape/__init__.py +++ b/microgen/shape/__init__.py @@ -21,7 +21,7 @@ from typing import TYPE_CHECKING, Callable, Literal, Sequence, Tuple -from . import implicit_basic_factory, implicit_shape, surface_functions +from . import implicit_shape, surface_functions from .box import Box from .capsule import Capsule from .cylinder import Cylinder @@ -177,7 +177,6 @@ def __init__(self: ShapeError, shape: str) -> None: "TruncatedCube", "TruncatedCuboctahedron", "TruncatedOctahedron", - "implicit_basic_factory", "implicit_shape", "new_geometry", "newGeometry", diff --git a/microgen/shape/box.py b/microgen/shape/box.py index 8f63ee1b..bd73c364 100644 --- a/microgen/shape/box.py +++ b/microgen/shape/box.py @@ -88,7 +88,7 @@ def generate_vtk( def to_implicit(self: Box) -> ImplicitShape: """Convert this box to an :class:`ImplicitShape`.""" - from .implicit_basic_factory import implicit_box + from .implicit_shape import implicit_box half = (self.dim[0] / 2, self.dim[1] / 2, self.dim[2] / 2) shape = implicit_box(center=(0, 0, 0), half_extents=half) diff --git a/microgen/shape/capsule.py b/microgen/shape/capsule.py index 589a5760..7d6d156e 100644 --- a/microgen/shape/capsule.py +++ b/microgen/shape/capsule.py @@ -111,7 +111,7 @@ def generate_vtk( def to_implicit(self: Capsule) -> ImplicitShape: """Convert this capsule to an :class:`ImplicitShape`.""" - from .implicit_basic_factory import implicit_capsule + from .implicit_shape import implicit_capsule h = self.height / 2.0 shape = implicit_capsule(start=(-h, 0, 0), end=(h, 0, 0), radius=self.radius) diff --git a/microgen/shape/cylinder.py b/microgen/shape/cylinder.py index 95de53d8..d162edb5 100644 --- a/microgen/shape/cylinder.py +++ b/microgen/shape/cylinder.py @@ -75,7 +75,7 @@ def generate_vtk( def to_implicit(self: Cylinder) -> ImplicitShape: """Convert this cylinder to an :class:`ImplicitShape`.""" - from .implicit_basic_factory import implicit_cylinder + from .implicit_shape import implicit_cylinder shape = implicit_cylinder( center=(0, 0, 0), diff --git a/microgen/shape/ellipsoid.py b/microgen/shape/ellipsoid.py index 705e4102..23789d1e 100644 --- a/microgen/shape/ellipsoid.py +++ b/microgen/shape/ellipsoid.py @@ -92,7 +92,7 @@ def generate_vtk(self: Ellipsoid, **_: KwargsGenerateType) -> pv.PolyData: def to_implicit(self: Ellipsoid) -> ImplicitShape: """Convert this ellipsoid to an :class:`ImplicitShape`.""" - from .implicit_basic_factory import implicit_ellipsoid + from .implicit_shape import implicit_ellipsoid shape = implicit_ellipsoid(center=(0, 0, 0), radii=self.radii) angles = tuple(self.orientation.as_euler("ZXZ", degrees=True)) diff --git a/microgen/shape/implicit_basic_factory.py b/microgen/shape/implicit_basic_factory.py deleted file mode 100644 index db88ba42..00000000 --- a/microgen/shape/implicit_basic_factory.py +++ /dev/null @@ -1,249 +0,0 @@ -"""Implicit primitive factory functions. - -===================================================================== -Implicit Basic Factory (:mod:`microgen.shape.implicit_basic_factory`) -===================================================================== -""" - -from __future__ import annotations - -import numpy as np -import numpy.typing as npt - -from .implicit_shape import ImplicitShape - - -def implicit_sphere( - center: tuple[float, float, float] = (0, 0, 0), - radius: float = 1.0, -) -> ImplicitShape: - """Sphere SDF: ``|p - c| - r``.""" - cx, cy, cz = center - r = radius - margin = r * 1.1 - return ImplicitShape( - func=lambda x, y, z, _cx=cx, _cy=cy, _cz=cz, _r=r: np.sqrt( - (x - _cx) ** 2 + (y - _cy) ** 2 + (z - _cz) ** 2 - ) - - _r, - bounds=( - cx - margin, - cx + margin, - cy - margin, - cy + margin, - cz - margin, - cz + margin, - ), - ) - - -def implicit_box( - center: tuple[float, float, float] = (0, 0, 0), - half_extents: tuple[float, float, float] = (0.5, 0.5, 0.5), -) -> ImplicitShape: - """Axis-aligned box SDF (exact).""" - cx, cy, cz = center - hx, hy, hz = half_extents - margin = max(hx, hy, hz) * 0.1 - - def _box_sdf( - x: npt.NDArray[np.float64], - y: npt.NDArray[np.float64], - z: npt.NDArray[np.float64], - ) -> npt.NDArray[np.float64]: - qx = np.abs(x - cx) - hx - qy = np.abs(y - cy) - hy - qz = np.abs(z - cz) - hz - outside = np.sqrt( - np.maximum(qx, 0.0) ** 2 - + np.maximum(qy, 0.0) ** 2 - + np.maximum(qz, 0.0) ** 2 - ) - inside = np.minimum(np.maximum(qx, np.maximum(qy, qz)), 0.0) - return outside + inside - - return ImplicitShape( - func=_box_sdf, - bounds=( - cx - hx - margin, - cx + hx + margin, - cy - hy - margin, - cy + hy + margin, - cz - hz - margin, - cz + hz + margin, - ), - ) - - -def implicit_cylinder( - center: tuple[float, float, float] = (0, 0, 0), - axis: tuple[float, float, float] = (0, 0, 1), - radius: float = 0.5, - height: float = 1.0, -) -> ImplicitShape: - """Capped cylinder SDF along arbitrary axis.""" - c = np.asarray(center, dtype=np.float64) - a = np.asarray(axis, dtype=np.float64) - a = a / np.linalg.norm(a) - half_h = height / 2.0 - - def _cyl_sdf( - x: npt.NDArray[np.float64], - y: npt.NDArray[np.float64], - z: npt.NDArray[np.float64], - ) -> npt.NDArray[np.float64]: - px, py, pz = x - c[0], y - c[1], z - c[2] - proj = px * a[0] + py * a[1] + pz * a[2] - perp_x = px - proj * a[0] - perp_y = py - proj * a[1] - perp_z = pz - proj * a[2] - d_radial = np.sqrt(perp_x**2 + perp_y**2 + perp_z**2) - radius - d_axial = np.abs(proj) - half_h - outside = np.sqrt( - np.maximum(d_radial, 0.0) ** 2 + np.maximum(d_axial, 0.0) ** 2 - ) - inside = np.minimum(np.maximum(d_radial, d_axial), 0.0) - return outside + inside - - extent = max(radius, half_h) * 1.1 + np.linalg.norm(c) - return ImplicitShape( - func=_cyl_sdf, - bounds=( - c[0] - extent, - c[0] + extent, - c[1] - extent, - c[1] + extent, - c[2] - extent, - c[2] + extent, - ), - ) - - -def implicit_plane( - point: tuple[float, float, float] = (0, 0, 0), - normal: tuple[float, float, float] = (0, 0, 1), -) -> ImplicitShape: - """Half-space SDF: ``dot(p - point, normal)``.""" - pt = np.asarray(point, dtype=np.float64) - n = np.asarray(normal, dtype=np.float64) - n = n / np.linalg.norm(n) - - return ImplicitShape( - func=lambda x, y, z: (x - pt[0]) * n[0] - + (y - pt[1]) * n[1] - + (z - pt[2]) * n[2], - bounds=None, - ) - - -def implicit_torus( - center: tuple[float, float, float] = (0, 0, 0), - major_r: float = 1.0, - minor_r: float = 0.25, -) -> ImplicitShape: - """Torus SDF in XZ plane.""" - cx, cy, cz = center - - def _torus_sdf( - x: npt.NDArray[np.float64], - y: npt.NDArray[np.float64], - z: npt.NDArray[np.float64], - ) -> npt.NDArray[np.float64]: - px, py, pz = x - cx, y - cy, z - cz - q_xz = np.sqrt(px**2 + pz**2) - major_r - return np.sqrt(q_xz**2 + py**2) - minor_r - - margin = (major_r + minor_r) * 1.1 - return ImplicitShape( - func=_torus_sdf, - bounds=( - cx - margin, - cx + margin, - cy - margin, - cy + margin, - cz - margin, - cz + margin, - ), - ) - - -def implicit_capsule( - start: tuple[float, float, float] = (0, 0, -0.5), - end: tuple[float, float, float] = (0, 0, 0.5), - radius: float = 0.25, -) -> ImplicitShape: - """Capsule SDF (line segment distance minus radius).""" - a = np.asarray(start, dtype=np.float64) - b = np.asarray(end, dtype=np.float64) - ab = b - a - ab_dot = float(np.dot(ab, ab)) - - if ab_dot < 1e-30: - # Degenerate capsule (start == end): fall back to sphere - return implicit_sphere(center=tuple(a), radius=radius) - - def _capsule_sdf( - x: npt.NDArray[np.float64], - y: npt.NDArray[np.float64], - z: npt.NDArray[np.float64], - ) -> npt.NDArray[np.float64]: - px, py, pz = x - a[0], y - a[1], z - a[2] - t = np.clip((px * ab[0] + py * ab[1] + pz * ab[2]) / ab_dot, 0.0, 1.0) - dx = px - t * ab[0] - dy = py - t * ab[1] - dz = pz - t * ab[2] - return np.sqrt(dx**2 + dy**2 + dz**2) - radius - - mid = (a + b) / 2.0 - extent = np.linalg.norm(ab) / 2.0 + radius - margin = extent * 1.1 - return ImplicitShape( - func=_capsule_sdf, - bounds=( - mid[0] - margin, - mid[0] + margin, - mid[1] - margin, - mid[1] + margin, - mid[2] - margin, - mid[2] + margin, - ), - ) - - -def implicit_ellipsoid( - center: tuple[float, float, float] = (0, 0, 0), - radii: tuple[float, float, float] = (1.0, 0.5, 0.25), -) -> ImplicitShape: - """Ellipsoid implicit field: ``sqrt((x/rx)^2 + (y/ry)^2 + (z/rz)^2) - 1``. - - This is an approximate SDF but has the correct zero-level set. - """ - cx, cy, cz = center - rx, ry, rz = radii - margin = max(rx, ry, rz) * 1.1 - - def _ellipsoid_field( - x: npt.NDArray[np.float64], - y: npt.NDArray[np.float64], - z: npt.NDArray[np.float64], - ) -> npt.NDArray[np.float64]: - return ( - np.sqrt( - ((x - cx) / rx) ** 2 - + ((y - cy) / ry) ** 2 - + ((z - cz) / rz) ** 2 - ) - - 1.0 - ) - - return ImplicitShape( - func=_ellipsoid_field, - bounds=( - cx - margin, - cx + margin, - cy - margin, - cy + margin, - cz - margin, - cz + margin, - ), - ) diff --git a/microgen/shape/implicit_shape.py b/microgen/shape/implicit_shape.py index 9ff341a8..bd772f9d 100644 --- a/microgen/shape/implicit_shape.py +++ b/microgen/shape/implicit_shape.py @@ -451,17 +451,244 @@ def _batched( # --------------------------------------------------------------------------- -# Primitive factory re-exports (moved to implicit_basic_factory.py) +# Primitive factory functions # --------------------------------------------------------------------------- -from .implicit_basic_factory import ( # noqa: E402 - implicit_box, - implicit_capsule, - implicit_cylinder, - implicit_ellipsoid, - implicit_plane, - implicit_sphere, - implicit_torus, -) + + +def implicit_sphere( + center: tuple[float, float, float] = (0, 0, 0), + radius: float = 1.0, +) -> ImplicitShape: + """Sphere SDF: ``|p - c| - r``.""" + cx, cy, cz = center + r = radius + margin = r * 1.1 + return ImplicitShape( + func=lambda x, y, z, _cx=cx, _cy=cy, _cz=cz, _r=r: np.sqrt( + (x - _cx) ** 2 + (y - _cy) ** 2 + (z - _cz) ** 2 + ) + - _r, + bounds=( + cx - margin, + cx + margin, + cy - margin, + cy + margin, + cz - margin, + cz + margin, + ), + ) + + +def implicit_box( + center: tuple[float, float, float] = (0, 0, 0), + half_extents: tuple[float, float, float] = (0.5, 0.5, 0.5), +) -> ImplicitShape: + """Axis-aligned box SDF (exact).""" + cx, cy, cz = center + hx, hy, hz = half_extents + margin = max(hx, hy, hz) * 0.1 + + def _box_sdf( + x: npt.NDArray[np.float64], + y: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], + ) -> npt.NDArray[np.float64]: + qx = np.abs(x - cx) - hx + qy = np.abs(y - cy) - hy + qz = np.abs(z - cz) - hz + outside = np.sqrt( + np.maximum(qx, 0.0) ** 2 + + np.maximum(qy, 0.0) ** 2 + + np.maximum(qz, 0.0) ** 2 + ) + inside = np.minimum(np.maximum(qx, np.maximum(qy, qz)), 0.0) + return outside + inside + + return ImplicitShape( + func=_box_sdf, + bounds=( + cx - hx - margin, + cx + hx + margin, + cy - hy - margin, + cy + hy + margin, + cz - hz - margin, + cz + hz + margin, + ), + ) + + +def implicit_cylinder( + center: tuple[float, float, float] = (0, 0, 0), + axis: tuple[float, float, float] = (0, 0, 1), + radius: float = 0.5, + height: float = 1.0, +) -> ImplicitShape: + """Capped cylinder SDF along arbitrary axis.""" + c = np.asarray(center, dtype=np.float64) + a = np.asarray(axis, dtype=np.float64) + a = a / np.linalg.norm(a) + half_h = height / 2.0 + + def _cyl_sdf( + x: npt.NDArray[np.float64], + y: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], + ) -> npt.NDArray[np.float64]: + px, py, pz = x - c[0], y - c[1], z - c[2] + proj = px * a[0] + py * a[1] + pz * a[2] + perp_x = px - proj * a[0] + perp_y = py - proj * a[1] + perp_z = pz - proj * a[2] + d_radial = np.sqrt(perp_x**2 + perp_y**2 + perp_z**2) - radius + d_axial = np.abs(proj) - half_h + outside = np.sqrt( + np.maximum(d_radial, 0.0) ** 2 + np.maximum(d_axial, 0.0) ** 2 + ) + inside = np.minimum(np.maximum(d_radial, d_axial), 0.0) + return outside + inside + + extent = max(radius, half_h) * 1.1 + np.linalg.norm(c) + return ImplicitShape( + func=_cyl_sdf, + bounds=( + c[0] - extent, + c[0] + extent, + c[1] - extent, + c[1] + extent, + c[2] - extent, + c[2] + extent, + ), + ) + + +def implicit_plane( + point: tuple[float, float, float] = (0, 0, 0), + normal: tuple[float, float, float] = (0, 0, 1), +) -> ImplicitShape: + """Half-space SDF: ``dot(p - point, normal)``.""" + pt = np.asarray(point, dtype=np.float64) + n = np.asarray(normal, dtype=np.float64) + n = n / np.linalg.norm(n) + + return ImplicitShape( + func=lambda x, y, z: (x - pt[0]) * n[0] + + (y - pt[1]) * n[1] + + (z - pt[2]) * n[2], + bounds=None, + ) + + +def implicit_torus( + center: tuple[float, float, float] = (0, 0, 0), + major_r: float = 1.0, + minor_r: float = 0.25, +) -> ImplicitShape: + """Torus SDF in XZ plane.""" + cx, cy, cz = center + + def _torus_sdf( + x: npt.NDArray[np.float64], + y: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], + ) -> npt.NDArray[np.float64]: + px, py, pz = x - cx, y - cy, z - cz + q_xz = np.sqrt(px**2 + pz**2) - major_r + return np.sqrt(q_xz**2 + py**2) - minor_r + + margin = (major_r + minor_r) * 1.1 + return ImplicitShape( + func=_torus_sdf, + bounds=( + cx - margin, + cx + margin, + cy - margin, + cy + margin, + cz - margin, + cz + margin, + ), + ) + + +def implicit_capsule( + start: tuple[float, float, float] = (0, 0, -0.5), + end: tuple[float, float, float] = (0, 0, 0.5), + radius: float = 0.25, +) -> ImplicitShape: + """Capsule SDF (line segment distance minus radius).""" + a = np.asarray(start, dtype=np.float64) + b = np.asarray(end, dtype=np.float64) + ab = b - a + ab_dot = float(np.dot(ab, ab)) + + if ab_dot < 1e-30: + # Degenerate capsule (start == end): fall back to sphere + return implicit_sphere(center=tuple(a), radius=radius) + + def _capsule_sdf( + x: npt.NDArray[np.float64], + y: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], + ) -> npt.NDArray[np.float64]: + px, py, pz = x - a[0], y - a[1], z - a[2] + t = np.clip((px * ab[0] + py * ab[1] + pz * ab[2]) / ab_dot, 0.0, 1.0) + dx = px - t * ab[0] + dy = py - t * ab[1] + dz = pz - t * ab[2] + return np.sqrt(dx**2 + dy**2 + dz**2) - radius + + mid = (a + b) / 2.0 + extent = np.linalg.norm(ab) / 2.0 + radius + margin = extent * 1.1 + return ImplicitShape( + func=_capsule_sdf, + bounds=( + mid[0] - margin, + mid[0] + margin, + mid[1] - margin, + mid[1] + margin, + mid[2] - margin, + mid[2] + margin, + ), + ) + + +def implicit_ellipsoid( + center: tuple[float, float, float] = (0, 0, 0), + radii: tuple[float, float, float] = (1.0, 0.5, 0.25), +) -> ImplicitShape: + """Ellipsoid implicit field: ``sqrt((x/rx)^2 + (y/ry)^2 + (z/rz)^2) - 1``. + + This is an approximate SDF but has the correct zero-level set. + """ + cx, cy, cz = center + rx, ry, rz = radii + margin = max(rx, ry, rz) * 1.1 + + def _ellipsoid_field( + x: npt.NDArray[np.float64], + y: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], + ) -> npt.NDArray[np.float64]: + return ( + np.sqrt( + ((x - cx) / rx) ** 2 + + ((y - cy) / ry) ** 2 + + ((z - cz) / rz) ** 2 + ) + - 1.0 + ) + + return ImplicitShape( + func=_ellipsoid_field, + bounds=( + cx - margin, + cx + margin, + cy - margin, + cy + margin, + cz - margin, + cz + margin, + ), + ) # --------------------------------------------------------------------------- diff --git a/microgen/shape/sphere.py b/microgen/shape/sphere.py index 947ac555..762fcf66 100644 --- a/microgen/shape/sphere.py +++ b/microgen/shape/sphere.py @@ -71,7 +71,7 @@ def generate_vtk( def to_implicit(self: Sphere) -> ImplicitShape: """Convert this sphere to an :class:`ImplicitShape`.""" - from .implicit_basic_factory import implicit_sphere + from .implicit_shape import implicit_sphere shape = implicit_sphere(center=(0, 0, 0), radius=self.radius) angles = tuple(self.orientation.as_euler("ZXZ", degrees=True))