"""Deterministic world-space point sampling inside prepared implicit hosts."""
# STD
import math
import numbers
# Third-party
import numpy as np
# OpenVCAD
import pyvcad as pv
# Local
from .fields import _coerce_attribute
[docs]
class StochasticPointSet:
"""Inspectable output from :func:`sample_points_in_host`."""
def __init__(self, points, local_spacings, seed, candidate_count, rejected_count):
#: Accepted points as a list of ``pv.Vec3``.
self.points = points
#: Local target spacing in mm at each accepted point.
self.local_spacings = local_spacings
#: Random seed the points were drawn with.
self.seed = seed
#: Number of random candidates drawn before relaxation.
self.candidate_count = candidate_count
#: Candidates that were not kept (``candidate_count - len(points)``).
self.rejected_count = rejected_count
def __len__(self):
return len(self.points)
def __iter__(self):
return iter(self.points)
def __getitem__(self, index):
return self.points[index]
def _field(value, name):
if isinstance(value, numbers.Real) and float(value) <= 0.0:
raise ValueError(f"{name} must be positive")
return _coerce_attribute(value)
def _sample_field(field, point, signed_distance, name):
value = float(field.sample(
float(point[0]), float(point[1]), float(point[2]), float(signed_distance)
))
if not math.isfinite(value) or value <= 0.0:
raise ValueError(f"{name} must sample to a positive finite value")
return value
def _inside_distance(host, point, boundary_clearance):
distance = host.evaluate(float(point[0]), float(point[1]), float(point[2]))
if distance is None or not math.isfinite(float(distance)):
return None
if float(distance) > -boundary_clearance:
return None
return float(distance)
def _is_separated(candidate, spacing, points, spacings, ignored_index=None):
if not points:
return True
coordinates = np.asarray(points, dtype=np.float64)
deltas = coordinates - candidate
squared = np.einsum("ij,ij->i", deltas, deltas)
required = 0.5 * (np.asarray(spacings, dtype=np.float64) + spacing)
if ignored_index is not None:
squared[ignored_index] = math.inf
return bool(np.all(squared >= required * required))
[docs]
def sample_points_in_host(
host,
spacing=None,
density=None,
seed=0,
max_count=None,
relaxation=0,
boundary_clearance=0.0,
attempts_per_point=80,
):
"""Sample deterministic, variably spaced world-space points in a prepared host.
Exactly one of ``spacing`` or ``density`` is required. Each may be a positive
number, ``FloatAttribute``, or scalar callable accepted by OpenVCAD. Density is
interpreted as target points per cubic world unit and converted locally to
``density ** (-1/3)`` spacing. Two samples are accepted when their separation is
at least the mean of their local spacings.
The host must already have been prepared. A point is inside when its host SDF is
at most ``-boundary_clearance``. Rejection sampling therefore remains correct for
disconnected and non-convex hosts, but a loose bounding box can be expensive.
``relaxation`` applies deterministic separation-preserving jitter passes.
"""
if (spacing is None) == (density is None):
raise ValueError("provide exactly one of spacing or density")
if max_count is not None and int(max_count) < 1:
raise ValueError("max_count must be at least one")
if int(relaxation) < 0:
raise ValueError("relaxation must be non-negative")
if float(boundary_clearance) < 0.0:
raise ValueError("boundary_clearance must be non-negative")
if int(attempts_per_point) < 1:
raise ValueError("attempts_per_point must be at least one")
minimum, maximum = host.bounding_box()
lower = np.array([minimum.x, minimum.y, minimum.z], dtype=np.float64)
upper = np.array([maximum.x, maximum.y, maximum.z], dtype=np.float64)
extent = upper - lower
if not np.all(np.isfinite(extent)) or np.any(extent <= 0.0):
raise ValueError("prepared host must have a finite, non-empty bounding box")
source_name = "spacing" if spacing is not None else "density"
source_field = _field(spacing if spacing is not None else density, source_name)
rng = np.random.default_rng(int(seed))
if max_count is None:
probe = 0.5 * (lower + upper)
probe_distance = _inside_distance(host, probe, float(boundary_clearance))
for _ in range(4096):
if probe_distance is not None:
break
probe = rng.uniform(lower, upper)
probe_distance = _inside_distance(host, probe, float(boundary_clearance))
if probe_distance is None:
raise RuntimeError(
"could not find an interior host point; check preparation, bbox, and clearance"
)
source_value = _sample_field(source_field, probe, probe_distance, source_name)
estimate_spacing = source_value if spacing is not None else source_value ** (-1.0 / 3.0)
max_count = max(4, min(10000, int(np.prod(extent) / estimate_spacing ** 3 * 2.0)))
else:
max_count = int(max_count)
accepted = []
local_spacings = []
candidate_count = 0
consecutive_rejections = 0
rejection_limit = max(200, int(attempts_per_point) * max(4, max_count))
while len(accepted) < max_count and consecutive_rejections < rejection_limit:
candidate_count += 1
candidate = rng.uniform(lower, upper)
signed_distance = _inside_distance(host, candidate, float(boundary_clearance))
if signed_distance is None:
consecutive_rejections += 1
continue
source_value = _sample_field(source_field, candidate, signed_distance, source_name)
local_spacing = source_value if spacing is not None else source_value ** (-1.0 / 3.0)
if _is_separated(candidate, local_spacing, accepted, local_spacings):
accepted.append(candidate)
local_spacings.append(local_spacing)
consecutive_rejections = 0
else:
consecutive_rejections += 1
for _ in range(int(relaxation)):
for index in rng.permutation(len(accepted)):
candidate = accepted[index] + rng.normal(0.0, 0.12 * local_spacings[index], 3)
signed_distance = _inside_distance(host, candidate, float(boundary_clearance))
if signed_distance is None:
continue
source_value = _sample_field(source_field, candidate, signed_distance, source_name)
local_spacing = source_value if spacing is not None else source_value ** (-1.0 / 3.0)
if _is_separated(
candidate, local_spacing, accepted, local_spacings, ignored_index=index
):
accepted[index] = candidate
local_spacings[index] = local_spacing
points = [pv.Vec3(float(point[0]), float(point[1]), float(point[2])) for point in accepted]
return StochasticPointSet(
points,
[float(value) for value in local_spacings],
int(seed),
candidate_count,
candidate_count - len(points),
)
