"""Random-to-regular Voronoi transitions by seed regularization.
A Voronoi diagram of a *regular* point set is itself a familiar strut lattice:
========== ========================================= ==========================
Sites Voronoi cells Matching graph unit cell
========== ========================================= ==========================
``"sc"`` cubes (sites at cell centres) ``"cubic"``
``"bcc"`` truncated octahedra (Kelvin cells) ``"kelvin"``
``"fcc"`` rhombic dodecahedra :func:`voronoi_unit_cell`
========== ========================================= ==========================
So a stochastic Voronoi lattice turns into a periodic strut lattice when its seeds move
from random positions onto such sites. The helpers below build those seeds from a
regularity field (0 = random, 1 = regular) and feed them to one ``StochasticGraphLattice``:
the random foam and the regular lattice are the same Voronoi graph, so they are connected
node to node by construction. Placing the sites through a ``CellMap`` makes the regular
region follow a curved part and line up with a mapped graph lattice of the matching cell.
"""
# STD
import math
import numbers
# Third-party
import numpy as np
# OpenVCAD
import pyvcad as pv
# Local
from .fields import _coerce_attribute
from .maps import _as_maps
from .points import StochasticPointSet
from .stochastic import voronoi_graph_lattice
SITE_KINDS = ("sc", "bcc", "fcc")
# Site offsets inside one unit cell, in cell units, chosen so the Voronoi cells line up with
# the cell boundaries of a graph lattice on the same grid.
_SITE_OFFSETS = {
"sc": ((0.5, 0.5, 0.5),),
"bcc": ((0.0, 0.0, 0.0), (0.5, 0.5, 0.5)),
"fcc": ((0.0, 0.0, 0.0), (0.5, 0.5, 0.0), (0.5, 0.0, 0.5), (0.0, 0.5, 0.5)),
}
# Nearest-neighbour distance of each site set for a unit cubic cell.
_NEAREST_NEIGHBOR = {"sc": 1.0, "bcc": math.sqrt(3.0) / 2.0, "fcc": 1.0 / math.sqrt(2.0)}
# Saturated Poisson-disk samples from pv.sample_poisson_seeds hold about this many points
# per spacing^3 (measured in a large box); used to match random and regular seed density.
_POISSON_DENSITY = 0.60
def _kind(value):
name = str(value).lower()
if name not in _SITE_OFFSETS:
raise ValueError("kind must be one of: " + ", ".join(SITE_KINDS))
return name
def _vec3_tuple(value):
if isinstance(value, pv.Vec3):
return (float(value.x), float(value.y), float(value.z))
if isinstance(value, numbers.Real):
return (float(value),) * 3
values = tuple(float(item) for item in value)
if len(values) != 3:
raise ValueError("expected a number or three coordinates")
return values
def _point_list(points):
if isinstance(points, StochasticPointSet):
return list(points.points)
if isinstance(points, TransitionSeeds):
return list(points.points)
if isinstance(points, np.ndarray):
return [pv.Vec3(float(x), float(y), float(z)) for x, y, z in points.reshape(-1, 3)]
return [point if isinstance(point, pv.Vec3) else pv.Vec3(*_vec3_tuple(point)) for point in points]
def _as_array(points):
return np.asarray([(p.x, p.y, p.z) for p in points], dtype=np.float64).reshape(-1, 3)
[docs]
def voronoi_unit_cell(kind):
"""Return the graph unit cell whose tiling equals the Voronoi edges of ``kind`` sites.
``"sc"`` gives the cubic cell, ``"bcc"`` the Kelvin cell (truncated octahedra) and
``"fcc"`` the rhombic-dodecahedron cell. The edges are in the same origin-centred,
normalized ``[-0.5, 0.5]`` convention as the built-in cells, so the result can be passed
to :func:`pyvcad_metamaterials.graph_lattice` with the same CellMap used for
:func:`cell_map_sites`. It is :func:`periodic_voronoi` of the sites, so other seed sets
work the same way.
"""
from .tessellations import periodic_voronoi
return periodic_voronoi(_SITE_OFFSETS[_kind(kind)]).edges()
[docs]
def lattice_sites(kind, bounds, cell_size, origin=None, margin=0.0):
"""Return regular ``"sc"``, ``"bcc"`` or ``"fcc"`` sites inside an axis-aligned box.
``bounds`` is a ``(minimum, maximum)`` pair; ``cell_size`` is a number or three numbers in
mm; ``origin`` is a lattice corner (defaults to ``minimum``). Sites within ``margin`` of
the box are included. To line the sites up with a rectangular map's graph lattice, use
the map's minimum corner as ``origin`` and its cell size.
"""
name = _kind(kind)
minimum = np.asarray(_vec3_tuple(bounds[0]))
maximum = np.asarray(_vec3_tuple(bounds[1]))
size = np.asarray(_vec3_tuple(cell_size))
if np.any(size <= 0.0) or not np.all(np.isfinite(size)):
raise ValueError("cell_size must be positive and finite")
if np.any(maximum <= minimum):
raise ValueError("bounds maximum must exceed minimum on every axis")
start = minimum if origin is None else np.asarray(_vec3_tuple(origin))
lower = minimum - float(margin)
upper = maximum + float(margin)
first = np.floor((lower - start) / size) - 1
last = np.ceil((upper - start) / size) + 1
axes = [np.arange(first[axis], last[axis] + 1) for axis in range(3)]
cells = np.stack(np.meshgrid(*axes, indexing="ij"), axis=-1).reshape(-1, 3)
sites = []
tolerance = 1e-9 * float(np.max(size))
for offset in _SITE_OFFSETS[name]:
positions = start + (cells + np.asarray(offset)) * size
inside = np.all((positions >= lower - tolerance) & (positions <= upper + tolerance), axis=1)
sites.append(positions[inside])
positions = np.concatenate(sites, axis=0)
order = np.lexsort((positions[:, 0], positions[:, 1], positions[:, 2]))
return [pv.Vec3(float(x), float(y), float(z)) for x, y, z in positions[order]]
def _extended_forward(cell_map, logical, counts):
"""Map a logical point, continuing the map linearly beyond its logical box."""
clamped = [min(max(logical[axis], 0.0), float(counts[axis])) for axis in range(3)]
base = cell_map.forward(pv.Vec3(*clamped))
if all(clamped[axis] == logical[axis] for axis in range(3)):
return base
du, dv, dw = cell_map.jacobian(pv.Vec3(*clamped))
step = [logical[axis] - clamped[axis] for axis in range(3)]
return pv.Vec3(
base.x + du.x * step[0] + dv.x * step[1] + dw.x * step[2],
base.y + du.y * step[0] + dv.y * step[1] + dw.y * step[2],
base.z + du.z * step[0] + dv.z * step[1] + dw.z * step[2],
)
[docs]
def cell_map_sites(cell_map, kind, margin=0.0):
"""Return ``kind`` sites placed through a CellMap ("align to graph points").
Sites are generated in logical cell coordinates and mapped to world space, so the
regular region of a transition follows the map, and its Voronoi cells line up with a
graph lattice of :func:`voronoi_unit_cell` on the same map. Sites in trimmed-away
elements are skipped, periodic seams are not duplicated, and a
``pv.CellMapCollection`` contributes the sites of every map.
``margin`` (in cells) also returns sites up to that far outside the map's logical box
on non-periodic axes, continuing the map linearly from its boundary. These outside
sites bound the Voronoi cells at the map surface (see :func:`voronoi_transition`).
Voronoi cells are only similarity-invariant: the alignment is exact for cubic cells on a
rectangular map and approximate where the map stretches or bends cells. Keep cells near
cubic (similar size along U, V and W).
"""
name = _kind(kind)
extra = float(margin)
if not math.isfinite(extra) or extra < 0.0:
raise ValueError("margin must be non-negative and finite")
result = []
for item in _as_maps(cell_map):
counts = item.cells
periodic = item.periodic_axes
for offset in _SITE_OFFSETS[name]:
ranges = []
for axis in range(3):
if periodic[axis]:
values = [index + offset[axis] for index in range(counts[axis])]
else:
low = -int(math.ceil(extra)) - 1
high = counts[axis] + int(math.ceil(extra)) + 1
values = [index + offset[axis] for index in range(low, high + 1)]
values = [value for value in values
if -extra - 1e-9 <= value <= counts[axis] + extra + 1e-9]
ranges.append(values)
for w in ranges[2]:
for v in ranges[1]:
for u in ranges[0]:
try:
result.append(_extended_forward(item, (u, v, w), counts))
except (RuntimeError, ValueError):
continue
return result
def _nearest_neighbor_distance(kind, cell_size=None, cell_map=None):
if cell_size is not None:
size = np.asarray(_vec3_tuple(cell_size))
return _NEAREST_NEIGHBOR[kind] * float(np.prod(size) ** (1.0 / 3.0))
lengths = []
for item in _as_maps(cell_map):
counts = item.cells
for u in np.linspace(0.25, counts[0] - 0.25, 5):
for v in np.linspace(0.25, counts[1] - 0.25, 5):
for w in np.linspace(0.25, counts[2] - 0.25, 3):
try:
du, dv, dw = item.jacobian(pv.Vec3(float(u), float(v), float(w)))
except (RuntimeError, ValueError):
continue
volume = abs(np.linalg.det(np.asarray([_vec3_tuple(du), _vec3_tuple(dv),
_vec3_tuple(dw)])))
if volume > 0.0 and math.isfinite(volume):
lengths.append(volume ** (1.0 / 3.0))
if not lengths:
raise ValueError("could not measure the CellMap's cell size")
return _NEAREST_NEIGHBOR[kind] * float(np.median(lengths))
[docs]
def matched_spacing(kind, cell_size):
"""Return the Poisson-disk spacing whose seed density equals ``kind`` sites.
``cell_size`` is the cubic cell edge in mm (or three edges). Random seeds sampled with
this spacing make Voronoi cells of the same average size as the regular lattice, so the
cell size stays continuous across a transition.
"""
name = _kind(kind)
size = np.asarray(_vec3_tuple(cell_size))
site_density = len(_SITE_OFFSETS[name]) / float(np.prod(size))
return (_POISSON_DENSITY / site_density) ** (1.0 / 3.0)
[docs]
def random_seeds(host, spacing, seed=0, max_count=None, boundary_clearance=0.0,
attempts_per_point=30):
"""Sample Poisson-disk seeds in a prepared host quickly (native, ~10^5 points).
``spacing`` is a positive number, ``FloatAttribute`` or scalar callable (mm). Two seeds
are accepted together when their distance is at least the mean of their local spacings,
the same rule as :func:`sample_points_in_host`, but the native sampler grows the point
set front by front instead of testing every pair, so it suits large hosts. Results are
deterministic for a given ``seed`` but differ from :func:`sample_points_in_host`.
Returns a :class:`StochasticPointSet`.
"""
if isinstance(spacing, numbers.Real) and float(spacing) <= 0.0:
raise ValueError("spacing must be positive")
result = pv.sample_poisson_seeds(
host,
_coerce_attribute(spacing),
seed=int(seed),
max_count=0 if max_count is None else int(max_count),
boundary_clearance=float(boundary_clearance),
attempts_per_point=int(attempts_per_point),
)
return StochasticPointSet(
list(result.points),
list(result.spacings),
int(seed),
int(result.candidate_count),
int(result.rejected_count),
)
[docs]
class TransitionSeeds:
"""Seeds from :func:`transition_seeds` with provenance for inspection.
Attributes:
points (list[pv.Vec3]): Seeds to pass to a Voronoi lattice.
regularity (list[float]): Regularity weight that placed each seed.
sites (list[pv.Vec3]): Regular sites inside the host.
random_points (list[pv.Vec3]): Random seeds used (empty for ``method="jitter"``).
regular_index, random_index (list[int]): Source site / random seed of each point or -1.
matched_count (int): Random seeds paired with a site (``method="match"``).
dropped_count (int): Seeds removed by the minimum-separation check.
nearest_neighbor (float): Site nearest-neighbour distance in mm.
"""
def __init__(self, result, sites, random_points, nearest_neighbor, method):
self.points = list(result.points)
self.regularity = list(result.regularity)
self.regular_index = list(result.regular_index)
self.random_index = list(result.random_index)
self.matched_count = int(result.matched_count)
self.dropped_count = int(result.separation_dropped_count)
self.sites = sites
self.random_points = random_points
self.nearest_neighbor = float(nearest_neighbor)
self.method = method
def __len__(self):
return len(self.points)
def __iter__(self):
return iter(self.points)
def __getitem__(self, index):
return self.points[index]
[docs]
def transition_seeds(
host,
regularity,
kind="bcc",
cell_size=None,
cell_map=None,
origin=None,
method="match",
random_points=None,
random_spacing=None,
jitter=0.5,
match_radius=None,
minimum_separation=None,
seed=0,
boundary_clearance=0.0,
site_margin=0.0,
):
"""Build Voronoi seeds that change from random to a regular lattice.
``regularity`` (number, ``FloatAttribute`` or callable) is 0 where the foam should be
random and 1 where it should be the regular lattice; ramps from :func:`ramp` work well.
Where it cannot be evaluated (for example outside a CellMap domain) it counts as 0.
The regular sites are ``kind`` (``"sc"``, ``"bcc"`` or ``"fcc"``) sites either on an
axis-aligned grid with ``cell_size`` (and optional lattice ``origin``) or, with
``cell_map``, placed through the map (see :func:`cell_map_sites`). Only sites inside the
prepared ``host`` are used.
``method="match"`` (default) pairs each site with a nearby random seed and slides the
seed onto the site by the local regularity, so the seed count and cell size stay
continuous. The random seeds are ``random_points`` (for example a graded
:func:`sample_points_in_host` result) or, by default, a Poisson-disk sample whose density
matches the sites (``random_spacing`` overrides it). ``method="jitter"`` instead displaces
every site by up to ``(1 - regularity) * jitter * nearest_neighbour`` (a perturbed
lattice), which needs no random set.
``match_radius`` (default 0.75 x the site nearest-neighbour distance) limits pairing and
``minimum_separation`` (default 0.3 x) removes seeds that would make sliver cells.
``site_margin`` (cells) extends CellMap sites past the map boundary (see
:func:`cell_map_sites`). Returns :class:`TransitionSeeds`.
"""
name = _kind(kind)
if method not in ("match", "jitter"):
raise ValueError("method must be 'match' or 'jitter'")
if (cell_size is None) == (cell_map is None):
raise ValueError("provide exactly one of cell_size or cell_map")
weight = _coerce_attribute(regularity).with_fallback(0.0)
neighbor = _nearest_neighbor_distance(name, cell_size=cell_size, cell_map=cell_map)
if cell_map is not None:
candidates = cell_map_sites(cell_map, name, margin=site_margin)
else:
minimum, maximum = host.bounding_box()
candidates = lattice_sites(name, (minimum, maximum), cell_size, origin=origin)
kept = pv.select_points_in_host(host, candidates, float(boundary_clearance))
sites = [candidates[index] for index in kept]
if not sites:
raise ValueError("no regular sites lie inside the host; check cell_size, origin and the host")
radius = 0.75 * neighbor if match_radius is None else float(match_radius)
separation = 0.3 * neighbor if minimum_separation is None else float(minimum_separation)
if method == "jitter":
result = pv.jitter_seed_sites(
host, sites, weight, float(jitter) * neighbor, seed=int(seed),
minimum_separation=separation, boundary_clearance=float(boundary_clearance))
return TransitionSeeds(result, sites, [], neighbor, method)
if random_points is None:
if random_spacing is None:
if cell_size is not None:
random_spacing = matched_spacing(name, cell_size)
else:
cube = neighbor / _NEAREST_NEIGHBOR[name]
random_spacing = matched_spacing(name, cube)
random_points = random_seeds(
host, random_spacing, seed=seed, boundary_clearance=boundary_clearance).points
else:
random_points = _point_list(random_points)
result = pv.blend_seed_sets(
host, random_points, sites, weight, radius,
minimum_separation=separation, boundary_clearance=float(boundary_clearance))
return TransitionSeeds(result, sites, random_points, neighbor, method)
[docs]
def voronoi_transition(
host,
regularity,
kind="bcc",
cell_size=None,
cell_map=None,
beam_thickness=1.0,
node_radius=None,
seeds=None,
ghost_layer=None,
boundary="open",
member_blend_radius=0.0,
field_tolerance=0.05,
boundary_tolerance=0.1,
prune_loose=True,
**seed_options
):
"""One call: a Voronoi strut lattice that turns from random foam into a regular lattice.
Builds :func:`transition_seeds` (pass seed options such as ``method``, ``seed`` or
``random_points`` as keywords) and a Voronoi graph lattice from them. Where
``regularity`` is 1 the struts are the edges of cubes (``"sc"``), Kelvin cells
(``"bcc"``) or rhombic dodecahedra (``"fcc"``). Both regions belong to one Voronoi graph,
so they share nodes and the transition itself cannot pinch apart.
``ghost_layer`` (mm, default one site spacing) also places seeds in a layer just outside
the host; the host still clips the geometry. Without that layer the outermost Voronoi
cells are unbounded, so regular cells at the surface lose the struts that lie on it.
``prune_loose=True`` (default) drops strut pieces that the host cuts off at both ends
(chords through a corner or a concave part of the host), which would otherwise print as
loose fragments. With precomputed ``seeds`` no ghost layer is added unless requested.
The host must be prepared. Check the result with
:func:`pyvcad_metamaterials.check_connectivity`; a very thin host can still split.
"""
if seeds is not None and seed_options:
raise ValueError("seed options cannot be combined with precomputed seeds")
if seeds is None:
name = _kind(kind)
if (cell_size is None) == (cell_map is None):
raise ValueError("provide exactly one of cell_size or cell_map")
neighbor = _nearest_neighbor_distance(name, cell_size=cell_size, cell_map=cell_map)
ghost = neighbor if ghost_layer is None else float(ghost_layer)
else:
ghost = 0.0 if ghost_layer is None else float(ghost_layer)
if not math.isfinite(ghost) or ghost < 0.0:
raise ValueError("ghost_layer must be non-negative and finite")
work_host = host
if ghost > 0.0:
work_host = pv.Offset(ghost, host)
step = 0.25 * ghost
work_host.prepare(pv.Vec3(step, step, step), 2.0 * ghost)
if seeds is None:
if cell_map is not None and ghost > 0.0:
seed_options.setdefault("site_margin", 1.0)
seeds = transition_seeds(work_host, regularity, kind=kind, cell_size=cell_size,
cell_map=cell_map, **seed_options)
return voronoi_graph_lattice(
host,
_point_list(seeds),
beam_thickness=beam_thickness,
node_radius=node_radius,
boundary=boundary,
field_tolerance=field_tolerance,
boundary_tolerance=boundary_tolerance,
member_blend_radius=member_blend_radius,
prune_loose=prune_loose,
allow_outside_seeds=True,
)
[docs]
def voronoi_tpms_transition(
host,
region,
tpms,
half_width,
regularity=0.0,
kind="bcc",
cell_size=None,
cell_map=None,
beam_thickness=1.0,
check_voxel_size=None,
**options
):
"""Join a Voronoi foam inside ``region`` to a TPMS outside it through an overlap band.
A TPMS sheet and a strut foam have nothing in common to blend, so this keeps both
structures in a ``2 * half_width`` band around the region surface and lets them
intersect (see :func:`pyvcad_metamaterials.overlap`). The foam is built by
:func:`voronoi_transition` on ``host`` grown into the band, with loose strut pieces
pruned against exactly that boundary, so the band edge leaves no loose fragments.
``regularity`` (0 = random foam) and ``kind``/``cell_size``/``cell_map`` set the foam's
seeds; use the TPMS's own ``cell_map`` so the foam cells match the TPMS cell size.
``beam_thickness`` and the TPMS wall thickness may be fields, for example ramps that
thin each structure as it crosses the band. Extra ``options`` go to
:func:`voronoi_transition`.
Choose ``half_width`` between a quarter and a half of the cell size: narrower bands give
few contacts, much wider ones double the material. With ``check_voxel_size`` the result
is checked with :func:`pyvcad_metamaterials.require_connected` (raises ``ValueError``
when it is not one piece). Returns the joined node, clipped to ``host``.
"""
from .connectivity import require_connected
from .mixing import overlap
width = float(half_width)
if not math.isfinite(width) or width <= 0.0:
raise ValueError("half_width must be positive and finite")
foam_host = pv.Intersection(host, pv.Offset(width, region))
step = 0.25 * width
foam_host.prepare(pv.Vec3(step, step, step), 2.0 * width)
foam = voronoi_transition(
foam_host,
regularity,
kind=kind,
cell_size=cell_size,
cell_map=cell_map,
beam_thickness=beam_thickness,
**options
)
result = pv.Intersection(overlap(foam, tpms, region, width), host)
if check_voxel_size is not None:
require_connected(result, check_voxel_size)
return result
[docs]
def voronoi_foam(
host,
spacing,
seed=0,
beam_thickness=1.0,
node_radius=None,
ghost_layer=None,
prune_loose=True,
member_blend_radius=0.0,
boundary="open",
field_tolerance=0.05,
boundary_tolerance=0.1,
):
"""One call: a random Voronoi strut foam that fills a prepared host in one piece.
Seeds come from :func:`random_seeds` with ``spacing`` (mm; a number, ``FloatAttribute``
or scalar callable), so a spacing field grades the cell size, for example smaller cells
near a head with ``mm.ramp(0, 15, coordinates=mm.distance_field(head))``. Seeds are also
placed in a ``ghost_layer`` just outside the host (default: the typical local spacing)
so the cells at the surface are closed, and ``prune_loose`` drops strut pieces the host
cuts off at both ends. A spacing field must therefore be defined a little outside the
host as well. The host can be any prepared node, including ``mm.domain(cell_map)`` for
a conformal liner volume. Returns ``voronoi_graph_lattice(host, ...)``.
"""
if ghost_layer is None:
probe = pv.sample_poisson_seeds(host, _coerce_attribute(spacing), seed=int(seed), max_count=64)
if len(probe.spacings) == 0:
raise ValueError("could not place seeds in the host; check preparation and spacing")
ghost = float(np.median(np.asarray(probe.spacings)))
else:
ghost = float(ghost_layer)
if not math.isfinite(ghost) or ghost < 0.0:
raise ValueError("ghost_layer must be non-negative and finite")
work_host = host
if ghost > 0.0:
work_host = pv.Offset(ghost, host)
step = 0.25 * ghost
work_host.prepare(pv.Vec3(step, step, step), 2.0 * ghost)
seeds = random_seeds(work_host, spacing, seed=seed)
return voronoi_graph_lattice(
host,
seeds,
beam_thickness=beam_thickness,
node_radius=node_radius,
boundary=boundary,
field_tolerance=field_tolerance,
boundary_tolerance=boundary_tolerance,
member_blend_radius=member_blend_radius,
prune_loose=prune_loose,
allow_outside_seeds=True,
)
