"""Voronoi ribs on a surface: seeds sampled on a triangle mesh and the cells they own on it.
``surface_seeds`` spreads Poisson-disk seeds over a ``pv.SurfaceMesh``. Each seed owns the
part of the surface closer to it than to any other seed (the restricted Voronoi diagram,
``surface_voronoi``). The boundaries between those cells become ribs: round struts that follow
the surface (``surface_voronoi_ribs``), or walls standing across it (``surface_voronoi_plates``).
"""
# STD
import math
import numbers
# Third-party
import numpy as np
# OpenVCAD
import pyvcad as pv
# Local
from .fields import _coerce_attribute
from .points import StochasticPointSet
from .stochastic import voronoi_plate_lattice
def _seed_array(seeds):
if isinstance(seeds, StochasticPointSet):
seeds = seeds.points
if isinstance(seeds, np.ndarray):
return np.asarray(seeds, dtype=np.float64).reshape(-1, 3)
return np.asarray([
(point.x, point.y, point.z) if isinstance(point, pv.Vec3) else tuple(point)
for point in seeds
], dtype=np.float64).reshape(-1, 3)
def _seed_vectors(seeds):
return [pv.Vec3(*map(float, point)) for point in _seed_array(seeds)]
[docs]
def surface_seeds(mesh, spacing, seed=0, max_count=None, attempts_per_point=30):
"""Sample variable-spacing Poisson-disk seeds on a triangle surface.
Darts land uniformly by area and are kept when their straight-line distance to every kept
seed is at least the mean of the two local spacings, until a long run of darts finds no
room. Each dart needs only its own triangle, so open, non-manifold and noisy scanned meshes
work. Keep the spacing below the thickness of thin walls: distances are straight-line, so
seeds on the two sides of a thin wall would see each other.
Args:
mesh (pv.SurfaceMesh): The surface. Build one from a node with
``pv.SurfaceMesh(node, voxel_size)`` or load a file with ``pv.SurfaceMesh(path)``.
spacing (float | pv.FloatAttribute | callable): Local spacing in mm, positive. A field
is sampled on the surface (with ``d = 0``), so it can grade the cell size.
seed (int): Random seed, default ``0``.
max_count (int | None): Largest number of seeds; ``None`` (default) fills the surface.
attempts_per_point (int): Rejected darts allowed per kept seed before sampling stops.
Returns:
StochasticPointSet: The seeds, with ``local_spacings`` and a ``triangle_ids`` list
giving the triangle each seed lies on.
Example:
>>> mesh = pv.SurfaceMesh(pv.Sphere(pv.Vec3(0, 0, 0), 20.0), 0.5)
>>> seeds = mm.surface_seeds(mesh, spacing=6.0, seed=1)
"""
if isinstance(spacing, numbers.Real) and float(spacing) <= 0.0:
raise ValueError("spacing must be positive")
result = pv.sample_surface_poisson_seeds(
mesh,
_coerce_attribute(spacing),
seed=int(seed),
max_count=0 if max_count is None else int(max_count),
attempts_per_point=int(attempts_per_point),
)
points = StochasticPointSet(
list(result.points),
list(result.spacings),
int(seed),
int(result.candidate_count),
int(result.rejected_count),
)
points.triangle_ids = list(result.triangle_ids)
return points
[docs]
def surface_voronoi(mesh, seeds):
"""Return the Voronoi diagram of seeds restricted to a triangle surface.
Each seed owns the part of the surface closer to it (in a straight line) than to any other
seed. The result's ``segments`` (an ``(N, 2, 3)`` array) are the boundaries between those
cells, one straight piece per triangle crossed, with ``seed_pairs``, ``triangle_ids`` and
each seed's ``cell_areas``. Computed per triangle, exactly, so scanned meshes work.
Args:
mesh (pv.SurfaceMesh): The surface.
seeds (StochasticPointSet | list | numpy.ndarray): Seed points, usually from
:func:`surface_seeds`.
Returns:
pv.SurfaceVoronoi: The restricted diagram.
Example:
>>> diagram = mm.surface_voronoi(mesh, seeds)
>>> diagram.segments.shape
"""
return pv.restricted_voronoi(mesh, _seed_array(seeds))
def _rib_graph(segments, tolerance):
"""Weld segment ends closer than tolerance; return vertices and unique edges."""
keys = {}
vertices = []
edges = set()
for start, end in segments:
ids = []
for point in (start, end):
key = tuple(np.round(point / tolerance).astype(np.int64))
if key not in keys:
keys[key] = len(vertices)
vertices.append(point)
ids.append(keys[key])
if ids[0] != ids[1]:
edges.add((min(ids), max(ids)))
return np.asarray(vertices, dtype=np.float64).reshape(-1, 3), sorted(edges)
def _simplify(vertices, edges, tolerance):
"""Replace chains of two-strut nodes by fewer straight struts within tolerance."""
neighbours = [[] for _ in range(len(vertices))]
for first, second in edges:
neighbours[first].append(second)
neighbours[second].append(first)
used = set()
result = []
def keep_chain(chain):
# Douglas-Peucker on the polyline of vertex indices.
stack = [(0, len(chain) - 1)]
kept = {0, len(chain) - 1}
while stack:
low, high = stack.pop()
if high - low < 2:
continue
start = vertices[chain[low]]
step = vertices[chain[high]] - start
length_squared = float(step @ step)
worst, worst_distance = None, tolerance
for position in range(low + 1, high):
offset = vertices[chain[position]] - start
if length_squared > 0.0:
offset = offset - (offset @ step) / length_squared * step
distance = float(np.linalg.norm(offset))
if distance > worst_distance:
worst, worst_distance = position, distance
if worst is not None:
kept.add(worst)
stack.extend([(low, worst), (worst, high)])
ordered = sorted(kept)
for low, high in zip(ordered, ordered[1:]):
result.append((chain[low], chain[high]))
def walk(start, following):
chain = [start, following]
used.add((min(start, following), max(start, following)))
while len(neighbours[chain[-1]]) == 2 and chain[-1] != start:
current = chain[-1]
nxt = neighbours[current][0] if neighbours[current][0] != chain[-2] else neighbours[current][1]
edge = (min(current, nxt), max(current, nxt))
if edge in used:
break
used.add(edge)
chain.append(nxt)
return chain
for vertex in range(len(vertices)):
if len(neighbours[vertex]) == 2:
continue
for following in neighbours[vertex]:
if (min(vertex, following), max(vertex, following)) not in used:
keep_chain(walk(vertex, following))
# Closed loops made only of two-strut nodes.
for first, second in edges:
if (first, second) not in used:
chain = walk(first, second)
if len(chain) > 3 and chain[-1] == chain[0]:
middle = len(chain) // 2
keep_chain(chain[:middle + 1])
keep_chain(chain[middle:])
else:
keep_chain(chain)
used_vertices = sorted({index for edge in result for index in edge})
remap = {index: position for position, index in enumerate(used_vertices)}
return vertices[used_vertices], sorted({(remap[a], remap[b]) for a, b in result})
[docs]
def surface_voronoi_ribs(
mesh,
seeds,
beam_thickness=1.0,
node_radius=None,
simplify_tolerance=0.05,
):
"""Build round ribs along the Voronoi cell boundaries of seeds on a surface.
The ribs follow the restricted Voronoi diagram (:func:`surface_voronoi`): the struts are
centred on the surface, so half of each rib stands proud of it and half sinks in. Union
them with the part to stiffen its skin, or use them alone as a surface net. Chains of
short pieces are merged into longer straight struts that stay within
``simplify_tolerance`` of the exact boundary.
Args:
mesh (pv.SurfaceMesh): The surface.
seeds (StochasticPointSet | list | numpy.ndarray): Seed points, usually from
:func:`surface_seeds`.
beam_thickness (float | pv.FloatAttribute | callable): Rib diameter in mm. A field is
sampled at both ends of every strut (with ``d = 0``), so the thickness can vary
along the surface.
node_radius (float | pv.FloatAttribute | callable | None): Radius in mm of a ball at
every junction of three or more ribs; ``None`` (default) adds none.
simplify_tolerance (float): Largest distance in mm a merged strut may leave the exact
boundary, default ``0.05``. ``0`` keeps one strut per triangle crossed.
Returns:
pv.BeamLattice: The ribs (new geometry; its distance is exact outside the struts and a
bound inside, like any union of struts).
Example:
>>> seeds = mm.surface_seeds(mesh, spacing=6.0)
>>> ribs = mm.surface_voronoi_ribs(mesh, seeds, beam_thickness=1.2)
"""
diagram = surface_voronoi(mesh, seeds)
segments = diagram.segments
if len(segments) == 0:
raise ValueError("the seeds' cells have no boundaries on this surface; use more seeds")
minimum, maximum = mesh.bounding_box()
size = math.dist((minimum.x, minimum.y, minimum.z), (maximum.x, maximum.y, maximum.z))
vertices, edges = _rib_graph(segments, max(1e-9, 1e-7 * size))
tolerance = float(simplify_tolerance)
if not math.isfinite(tolerance) or tolerance < 0.0:
raise ValueError("simplify_tolerance must be non-negative")
if tolerance > 0.0:
vertices, edges = _simplify(vertices, edges, tolerance)
thickness = _coerce_attribute(beam_thickness)
radii = []
for x, y, z in vertices:
value = float(thickness.sample(float(x), float(y), float(z), 0.0))
if not math.isfinite(value) or value <= 0.0:
raise ValueError("beam_thickness must sample to a positive finite value")
radii.append(0.5 * value)
beams = [pv.Beam(int(a), int(b), r1=radii[a], r2=radii[b]) for a, b in edges]
points = [pv.Vec3(float(x), float(y), float(z)) for x, y, z in vertices]
if node_radius is None:
return pv.BeamLattice(points, beams, max(radii))
degree = np.zeros(len(vertices), dtype=np.int64)
for a, b in edges:
degree[a] += 1
degree[b] += 1
joint = _coerce_attribute(node_radius)
balls = []
for index in np.flatnonzero(degree >= 3):
x, y, z = vertices[index]
balls.append(pv.Ball(int(index), radius=float(joint.sample(float(x), float(y), float(z), 0.0))))
return pv.BeamLattice(points, beams, max(radii), ball_mode=pv.BallMode.Mixed,
ball_radius=max(ball.radius for ball in balls) if balls else 0.0, balls=balls)
[docs]
def surface_voronoi_plates(
body,
seeds,
rib_height,
wall_thickness=1.0,
member_blend_radius=0.0,
surface_tolerance=0.05,
):
"""Build walls along the Voronoi cell boundaries of seeds on the surface of a body.
The walls are the faces of the seeds' 3D Voronoi cells, kept within ``rib_height / 2`` of
the body's surface (inside ``pv.Thicken(rib_height, body)``). For seeds on the surface those
faces cut it exactly along the restricted Voronoi diagram, so every wall is centred on a cell
boundary and stands across the surface. Use the body's mesh for the seeds, for example
``mm.surface_seeds(pv.SurfaceMesh(body, 0.5), 6.0)``.
Args:
body (pv.Node): A closed body; its surface carries the walls.
seeds (StochasticPointSet | list | numpy.ndarray): Seeds on the body's surface.
rib_height (float | pv.FloatAttribute | pv.ScalarParameter | callable): Total wall
height in mm across the surface (half outside, half inside).
wall_thickness (float | pv.FloatAttribute | callable): Wall thickness in mm; a field may
vary it along the surface.
member_blend_radius (float): Fillet radius in mm where walls meet, default ``0``.
surface_tolerance (float): Largest thickness interpolation error in mm.
Returns:
pv.StochasticFaceLattice: The walls.
Example:
>>> walls = mm.surface_voronoi_plates(body, seeds, rib_height=3.0, wall_thickness=0.8)
"""
band = pv.Thicken(rib_height, body)
return voronoi_plate_lattice(
band,
_seed_vectors(seeds),
wall_thickness=wall_thickness,
surface_tolerance=surface_tolerance,
member_blend_radius=member_blend_radius,
)
