# Copyright (c) 2026, pinkfish
#
# Licensed under the BSD 2-Clause License. See the LICENSE file in the project
# root for the full license text.
# SPDX-License-Identifier: BSD-2-Clause
# LibFile: pybosl2/vnf.py
# Minimal pure-Python port of BOSL2's VNF ("Vertices and Faces") structure
# from vnf.scad -- just the pieces the bezier surface functions
# (pybosl2/beziers.py's BezierPatch) need to turn a grid of surface sample
# points into a polyhedron: vnf_vertex_array() (grid -> VNF with the quad
# subdivision styles), vnf_join() (merge VNFs), and rendering to PythonSCAD's
# native polyhedron(). No osuse()/BOSL2 runtime dependency.
#
# A VNF is [vertices, faces]: vertices a list of 3-D points, faces a list of
# index lists (each a polygon into `vertices`). That maps straight onto
# OpenSCAD's polyhedron(points=, faces=). The class carries the pair and, like
# Path2D/Bezier, keeps every operation as a method.
#
# FileSummary: VNF (vertices+faces) surface structure and grid meshing (BOSL2 vnf.scad).
# DocCategory: Paths, regions & surfaces
# FileGroup: BOSL2
"""VNF (vertices+faces) surface structure and grid meshing (BOSL2 vnf.scad)."""
from __future__ import annotations
import math
from collections import defaultdict
from typing import TYPE_CHECKING, Any
import numpy as np
from pybosl2._mctable import CORNER_OFFSETS, EDGE_CORNERS, TRI_TABLE
from pybosl2.bounds import Bounds2D, Bounds3D
from pybosl2.enums import SamplingType, SkinMethod, VNFStyle
if TYPE_CHECKING:
from collections.abc import Callable, Sequence
from numpy.typing import NDArray
from pybosl2.caps import CapSpec, CapsSpec
from pybosl2.isosurface import _MetaballSpec
from pybosl2.path3d import Path3D
from pybosl2.shapes3d import Bosl2Solid
_EPS = 1e-9
def _plane_edge_t(
pt0: list[float] | np.ndarray,
pt1: list[float] | np.ndarray,
a: float,
b: float,
c: float,
d: float,
) -> float:
"""Return parametric t (0..1) where the edge pt0→pt1 crosses the plane A*x+B*y+C*z=D."""
d0 = a * pt0[0] + b * pt0[1] + c * pt0[2] - d
d1 = a * pt1[0] + b * pt1[1] + c * pt1[2] - d
denom = d0 - d1
if abs(denom) < _EPS:
return 0.5
return d0 / denom
def _interpolate(
pt0: list[float] | np.ndarray,
pt1: list[float] | np.ndarray,
t: float,
) -> list[float]:
"""Linear interpolation between two 3-D points."""
return [
pt0[0] + t * (pt1[0] - pt0[0]),
pt0[1] + t * (pt1[1] - pt0[1]),
pt0[2] + t * (pt1[2] - pt0[2]),
]
def _triangle_area(
a: list[float] | np.ndarray,
b: list[float] | np.ndarray,
c: list[float] | np.ndarray,
) -> float:
"""Signed triangle area from three 3-D points (half the cross-product magnitude)."""
u = np.array(b, dtype=float) - np.array(a, dtype=float)
v = np.array(c, dtype=float) - np.array(a, dtype=float)
return float(np.linalg.norm(np.cross(u, v))) * 0.5
def _assemble_edge_paths(
edges: list[tuple[int, int]],
) -> list[list[int]]:
"""Assemble disconnected directed edges into closed loops.
Each edge ``(i, j)`` is treated as a directed connection i→j.
Returns a list of vertex-index paths forming closed polygons.
"""
if not edges:
return []
adj: dict[int, list[int]] = defaultdict(list)
for a, b in edges:
adj[a].append(b)
visited: set[int] = set()
paths: list[list[int]] = []
for start in list(adj):
if start in visited:
continue
path: list[int] = [start]
visited.add(start)
current = start
while True:
next_candidates = list(adj.get(current, []))
if not next_candidates:
break
nxt = next_candidates[0]
if nxt in visited:
if nxt == path[0]:
paths.append(path)
break
path.append(nxt)
visited.add(nxt)
current = nxt
# Handle any remaining edges not in a loop by assembling orphan paths
remaining: set[tuple[int, int]] = set(edges)
path_edges: set[tuple[int, int]] = set()
for p in paths:
for k1 in range(len(p)):
k2 = (k1 + 1) % len(p)
path_edges.add((p[k1], p[k2]))
remaining -= path_edges
return paths
# -- marching-squares lookup table -------------------------------------------
# In the two ambiguous cases with two opposite corners above and the other
# two below the isovalue, it is assumed the high values connect (ridge, not valley).
# This makes the contour compatible with marching cubes at pixel boundaries.
_MSQUARE_SEGMENT_TABLE: list[list[list[int]]] = [
[[], []],
[[0, 3], []],
[[1, 0], []],
[[1, 3], []],
[[3, 2], []],
[[0, 2], []],
[[1, 2], [3, 0]],
[[1, 2], []],
[[2, 1], []],
[[0, 1], [2, 3]],
[[2, 0], []],
[[2, 3], []],
[[3, 1], []],
[[0, 1], []],
[[3, 0], []],
[[], []],
]
_MSQUARE_VERTEX_INDEX_MAP: list[list[float]] = [
[0.0, 0.0],
[0.0, 1.0],
[1.0, 0.0],
[1.0, 1.0],
]
def _msquare_index(fvals: Sequence[float], isovalue: float) -> int:
"""Return 0..15 marching-square case index for 4 corner values."""
idx = 0
for i, v in enumerate(fvals):
if float(v) >= isovalue:
idx |= 1 << i
return idx
def _assemble_partial_paths_2d(
segments: list[list[list[float]]],
closed: bool,
) -> list[list[list[float]]]:
"""Assemble 2-D line segments into connected paths (contour polygons).
Each segment is ``[[x0,y0], [x1,y1]]``. Returns a list of paths,
each a list of ``[x, y]`` points. If *closed* is True, only closed
loops are kept; otherwise dangling paths are also returned open.
"""
graph: dict[tuple[float, float], list[tuple[float, float]]] = defaultdict(list)
for seg in segments:
if len(seg) < 2:
continue
p0 = (float(seg[0][0]), float(seg[0][1]))
p1 = (float(seg[1][0]), float(seg[1][1]))
if abs(p0[0] - p1[0]) < _EPS and abs(p0[1] - p1[1]) < _EPS:
continue
graph[p0].append(p1)
graph[p1].append(p0)
visited: set[tuple[float, float]] = set()
paths: list[list[list[float]]] = []
for start in graph:
if start in visited:
continue
path: list[tuple[float, float]] = [start]
visited.add(start)
curr = start
prev: tuple[float, float] | None = None
while True:
neigh = graph.get(curr, [])
nxt: tuple[float, float] | None = None
if len(path) == 1:
if neigh:
nxt = neigh[0]
else:
for n in neigh:
if n != prev:
nxt = n
break
if nxt is None:
break
if nxt == start:
paths.append([[float(x), float(y)] for x, y in path])
break
if nxt in visited:
break
path.append(nxt)
visited.add(nxt)
prev = curr
curr = nxt
if not closed:
for start2 in graph:
if start2 in visited:
continue
opath: list[tuple[float, float]] = [start2]
visited.add(start2)
cur = start2
while True:
neigh2 = [n for n in graph.get(cur, []) if n not in visited]
if not neigh2:
break
nx = neigh2[0]
opath.append(nx)
visited.add(nx)
cur = nx
if len(opath) >= 2:
paths.append([[float(x), float(y)] for x, y in opath])
return paths
def _marching_squares(
field: np.ndarray,
xs: np.ndarray,
ys: np.ndarray,
isovalue: float,
) -> list[list[list[float]]]:
"""Run marching squares on a 2-D scalar field, returning contour paths."""
nx, ny = field.shape
segments: list[list[list[float]]] = []
_edge_verts = [(0, 1), (1, 3), (2, 3), (0, 2)]
_vert_coords = [(0.0, 0.0), (0.0, 1.0), (1.0, 0.0), (1.0, 1.0)]
for i in range(nx - 1):
dx = xs[i + 1] - xs[i]
for j in range(ny - 1):
dy = ys[j + 1] - ys[j]
fvals = [
float(field[i, j]),
float(field[i, j + 1]),
float(field[i + 1, j]),
float(field[i + 1, j + 1]),
]
idx = _msquare_index(fvals, isovalue)
for edge_group in _MSQUARE_SEGMENT_TABLE[idx]:
if not edge_group:
continue
seg: list[list[float]] = []
for e in edge_group:
va, vb = _edge_verts[e]
fa, fb = fvals[va], fvals[vb]
denom = fb - fa
u = 0.5 if abs(denom) < _EPS else (isovalue - fa) / denom
ca = _vert_coords[va]
cb = _vert_coords[vb]
x = xs[i] + ca[0] * dx + u * (cb[0] - ca[0]) * dx
y = ys[j] + ca[1] * dy + u * (cb[1] - ca[1]) * dy
seg.append([float(x), float(y)])
if len(seg) == 2:
segments.append(seg)
paths = _assemble_partial_paths_2d(segments, closed=True)
for p in paths:
if p and p[0] == p[-1]:
p.pop()
return paths
def _to_grid(points: Any) -> np.ndarray:
"""Convert points (Path3D, list of Path3D, arrays) to a 3-D numpy grid."""
from pybosl2.path3d import Path3D
if isinstance(points, Path3D):
return np.array([list(p) for p in points], dtype=float)
if isinstance(points, list) and points and isinstance(points[0], Path3D):
return np.array([[list(p) for p in row] for row in points], dtype=float)
return np.asarray(points, dtype=float)
# -- marching-cubes helpers ----------------------------------------------------
def _resolve_grid(
bb: Bounds3D,
voxel_size: float | None,
voxel_count: int | None,
exact_bounds: bool,
) -> tuple[Bounds3D, float]:
import math
if voxel_size is None:
w, h, d = bb.max_x - bb.min_x, bb.max_y - bb.min_y, bb.max_z - bb.min_z
voxvol = (w * h * d) / (voxel_count if voxel_count else 22**3)
voxel_size = voxvol ** (1 / 3)
if exact_bounds:
return bb, voxel_size
vs = voxel_size
nx = math.ceil((bb.max_x - bb.min_x) / vs)
ny = math.ceil((bb.max_y - bb.min_y) / vs)
nz = math.ceil((bb.max_z - bb.min_z) / vs)
cx, cy, cz = (bb.min_x + bb.max_x) / 2, (bb.min_y + bb.max_y) / 2, (bb.min_z + bb.max_z) / 2
hx, hy, hz = 0.5 * vs * nx, 0.5 * vs * ny, 0.5 * vs * nz
return Bounds3D(
min_x=cx - hx,
min_y=cy - hy,
min_z=cz - hz,
max_x=cx + hx,
max_y=cy + hy,
max_z=cz + hz,
width=2 * hx,
length=2 * hy,
height=2 * hz,
), voxel_size
def _grid_axes_2d(
bb: Bounds2D,
pixel_size: float,
) -> tuple[np.ndarray, np.ndarray]:
"""Build uniform 2-D grid axes from a bounding box and pixel size."""
xs = np.arange(bb.min_x, bb.max_x + pixel_size * 0.5, pixel_size)
ys = np.arange(bb.min_y, bb.max_y + pixel_size * 0.5, pixel_size)
return xs, ys
def _grid_axes(bb: Bounds3D, voxel_size: float) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
import math
def axis(lo: float, hi: float, step: float) -> np.ndarray:
sides: int = int(math.floor((hi - lo) / step + 0.5)) + 1
return lo + step * np.arange(sides)
return (
axis(bb.min_x, bb.max_x, voxel_size),
axis(bb.min_y, bb.max_y, voxel_size),
axis(bb.min_z, bb.max_z, voxel_size),
)
def _resolve_grid_2d(
bb: Bounds2D,
pixel_size: float | None,
pixel_count: int | None,
exact_bounds: bool,
) -> tuple[Bounds2D, float]:
"""Resolve 2-D grid parameters from a bounding box and optional pixel size/count."""
if pixel_size is None:
w, h = bb.max_x - bb.min_x, bb.max_y - bb.min_y
pixvol = (w * h) / (pixel_count if pixel_count else 32**2)
pixel_size = math.sqrt(pixvol)
if exact_bounds:
return bb, pixel_size
vs = pixel_size
nx = math.ceil((bb.max_x - bb.min_x) / vs)
ny = math.ceil((bb.max_y - bb.min_y) / vs)
cx = (bb.min_x + bb.max_x) / 2
cy = (bb.min_y + bb.max_y) / 2
hx = 0.5 * vs * nx
hy = 0.5 * vs * ny
return Bounds2D(
min_x=cx - hx,
min_y=cy - hy,
max_x=cx + hx,
max_y=cy + hy,
width=2 * hx,
length=2 * hy,
), vs
def _sample_field_2d(
f: np.ndarray | Callable[[np.ndarray], np.ndarray],
xs: np.ndarray,
ys: np.ndarray,
) -> np.ndarray:
"""Sample a 2-D scalar field on a grid, returning a 2-D numpy array."""
if isinstance(f, np.ndarray) or (isinstance(f, (list, tuple)) and not callable(f)):
return np.asarray(f, dtype=float)
gx, gy = np.meshgrid(xs, ys, indexing="ij")
pts = np.stack([gx.ravel(), gy.ravel()], axis=1)
with np.errstate(all="ignore"):
vals = np.asarray(f(pts), dtype=float)
if vals.shape == (len(pts),):
return vals.reshape(gx.shape)
raise TypeError("_sample_field_2d: callable f must accept (N,2) array and return (N,) array.")
[docs]
def contour(
f: np.ndarray | Callable[[np.ndarray], np.ndarray],
isovalue: float,
bounding_box: Bounds2D,
pixel_size: float | None = None,
pixel_count: int | None = None,
closed: bool = True,
exact_bounds: bool = False,
) -> list[list[list[float]]]:
"""Generate 2-D contour paths at a given isovalue from a scalar field.
Uses marching squares on a uniform 2-D grid to trace the contour where
``f(x, y) == isovalue``. Returns a list of closed (or open) polyline
paths, each being a list of ``[x, y]`` points.
Args:
f: A 2-D numpy array or a callable ``(N,2)→(N,)`` or ``(x,y)→float``.
isovalue: Scalar threshold.
bounding_box: A :class:`~pybosl2.bounds.Bounds2D`.
pixel_size: Isotropic pixel size.
pixel_count: Approximate total pixel count (ignored if *pixel_size* given).
closed: If True, return only closed contour loops.
exact_bounds: If True, use *bounding_box* exactly.
Returns:
A list of contour paths, each a list of ``[x, y]`` points.
Examples:
.. pythonscad-example::
import numpy as np
from pybosl2 import contour, Bounds2D
from pybosl2.path2d import Path2D
def field(p):
r = np.hypot(p[:, 0], p[:, 1])
return r
paths = contour(field, 10, Bounds2D(-15, -15, 15, 15, 30, 30), pixel_size=0.5)
Path2D(paths[0]).stroke(width=0.5).linear_extrude(height=2).show()
"""
bb, ps = _resolve_grid_2d(bounding_box, pixel_size, pixel_count, exact_bounds)
xs, ys = _grid_axes_2d(bb, ps)
field_arr = _sample_field_2d(f, xs, ys)
paths = _marching_squares(field_arr, xs, ys, float(isovalue))
if not closed:
return paths
return [p for p in paths if len(p) >= 3 and p[0] != p[-1]]
def _sample_field(
f: np.ndarray | Callable[[np.ndarray], np.ndarray],
xs: np.ndarray,
ys: np.ndarray,
zs: np.ndarray,
) -> np.ndarray:
if isinstance(f, np.ndarray) or (isinstance(f, (list, tuple)) and not callable(f)):
return np.asarray(f, dtype=float)
gx, gy, gz = np.meshgrid(xs, ys, zs, indexing="ij")
pts = np.stack([gx.ravel(), gy.ravel(), gz.ravel()], axis=1)
with np.errstate(all="ignore"):
try:
vals = np.asarray(f(pts), dtype=float)
if vals.shape == (len(pts),):
return vals.reshape(gx.shape)
except (ValueError, TypeError, IndexError, np.exceptions.AxisError):
pass
vals = np.array([float(f(np.array([p[0], p[1], p[2]]))) for p in pts])
return vals.reshape(gx.shape)
def _marching_cubes(
field: np.ndarray,
xs: np.ndarray,
ys: np.ndarray,
zs: np.ndarray,
iso: float,
closed: bool,
) -> tuple[list[list[float]], list[list[int]]]:
if closed:
field = np.pad(field, 1, mode="constant", constant_values=-1e30)
xs = np.concatenate([[xs[0] - (xs[1] - xs[0])], xs, [xs[-1] + (xs[-1] - xs[-2])]])
ys = np.concatenate([[ys[0] - (ys[1] - ys[0])], ys, [ys[-1] + (ys[-1] - ys[-2])]])
zs = np.concatenate([[zs[0] - (zs[1] - zs[0])], zs, [zs[-1] + (zs[-1] - zs[-2])]])
nx, ny, nz = field.shape
coords = (xs, ys, zs)
verts: list[list[float]] = []
faces: list[list[int]] = []
cache: dict[tuple[int, int, int, int, int, int], int] = {}
def corner_pos(ci: int, i: int, j: int, k: int) -> tuple[int, int, int]:
di, dj, dk = CORNER_OFFSETS[ci]
return (i + di, j + dj, k + dk)
def edge_vertex(ca: tuple[int, int, int], cb: tuple[int, int, int]) -> int:
ordered = (ca, cb) if ca < cb else (cb, ca)
key = (ordered[0][0], ordered[0][1], ordered[0][2], ordered[1][0], ordered[1][1], ordered[1][2])
idx = cache.get(key)
if idx is not None:
return idx
ia, ja, ka = ca
ib, jb, kb = cb
va, vb = field[ia, ja, ka], field[ib, jb, kb]
t: float = 0.5 if va == vb else (iso - va) / (vb - va)
pa = np.array([coords[0][ia], coords[1][ja], coords[2][ka]])
pb = np.array([coords[0][ib], coords[1][jb], coords[2][kb]])
idx = len(verts)
verts.append(list(pa + t * (pb - pa)))
cache[key] = idx
return idx
for i in range(nx - 1):
for j in range(ny - 1):
for k in range(nz - 1):
cvals = [
field[i + CORNER_OFFSETS[c][0], j + CORNER_OFFSETS[c][1], k + CORNER_OFFSETS[c][2]]
for c in range(8)
]
cubeindex: int = 0
for c in range(8):
if cvals[c] < iso:
cubeindex |= 1 << c
tris = TRI_TABLE[cubeindex]
if not tris:
continue
for t in range(0, len(tris), 3):
face: list[int] = []
for e in tris[t : t + 3]:
c0, c1 = EDGE_CORNERS[e]
face.append(edge_vertex(corner_pos(c0, i, j, k), corner_pos(c1, i, j, k)))
if face[0] != face[1] and face[1] != face[2] and face[0] != face[2]:
faces.append(face)
return verts, faces
def _count(sides: int, s: int = 0, reverse: bool = False) -> list[int]:
radius = list(range(s, s + sides))
return radius[::-1] if reverse else radius
def _lofttri(
p1: np.ndarray,
p2: np.ndarray,
i1off: int,
i2off: int,
n1: int,
n2: int,
reverse: bool,
trimax: float,
) -> list[list[int]]:
"""Triangulate between two rows (possibly unequal length) by shortest new edge.
BOSL2 _lofttri.
"""
a1 = np.asarray(p1, dtype=float)
a2 = np.asarray(p2, dtype=float)
tris: list[Any] = []
if n1 != n2:
i1 = i2 = 0
tc1 = tc2 = 0
while True:
t1 = i1 + 1 if i1 < n1 else n1
t2 = i2 + 1 if i2 < n2 else n2
d12 = 9e9 if t2 >= n2 else float(np.linalg.norm(a2[t2] - a1[i1]))
d21 = 9e9 if t1 >= n1 else float(np.linalg.norm(a1[t1] - a2[i2]))
userow = (2 if tc1 < trimax else 1) if d12 < d21 else (1 if tc2 < trimax else 2)
newt = (t1 if t1 < n1 else i1) if userow == 1 else (t2 if t2 < n2 else i2)
newofft = i2off + newt if userow == 2 else i1off + newt
tc1n = tc1 + 1 if (d12 < d21 and tc1 < trimax) else 0
tc2n = tc2 + 1 if (d21 < d12 and tc2 < trimax) else 0
triangle = [i1off + i1, i2off + i2, newofft] if reverse else [i2off + i2, i1off + i1, newofft]
if t1 >= n1 and t2 >= n2:
break
tris.append(triangle)
if userow == 1:
i1 = i1 if t1 >= n1 else t1
else:
i2 = i2 if t2 >= n2 else t2
tc1, tc2 = tc1n, tc2n
else:
sides = n1
i = 0
while True:
t = i + 1 if i < sides else sides
if t >= sides:
break
d12 = float(np.linalg.norm(a2[t] - a1[i]))
d21 = float(np.linalg.norm(a1[t] - a2[i]))
third1 = i2off + t if d12 < d21 else i1off + t
third2 = i1off + i if d12 < d21 else i2off + i
if reverse:
tris.append([i1off + i, i2off + i, third1])
tris.append([i2off + t, i1off + t, third2])
else:
tris.append([i2off + i, i1off + i, third1])
tris.append([i1off + t, i2off + t, third2])
i = t
return tris
def _rows_as_float(rows: Any) -> list[list[float]]:
"""Rows of plain floats.
An array converts in one C-level pass; rows that are already plain floats only need copying.
Converting a list of lists *through* numpy is slower than either, so it is not used here.
"""
if rows is None or len(rows) == 0:
return []
if isinstance(rows, np.ndarray):
return rows.astype(float, copy=False).tolist() # type: ignore[no-any-return]
first = rows[0]
if type(first) is list and (not first or type(first[0]) is float):
return [list(r) for r in rows]
return [[float(x) for x in r] for r in rows]
def _rows_as_int(rows: Any) -> list[list[int]]:
"""Rows of plain ints, converted the same way as :func:`_rows_as_float`."""
if rows is None or len(rows) == 0:
return []
if isinstance(rows, np.ndarray):
return rows.astype(np.intp, copy=False).tolist() # type: ignore[no-any-return]
first = rows[0]
if type(first) is list and (not first or type(first[0]) is int):
return [list(r) for r in rows]
return [[int(i) for i in r] for r in rows]
def _norms(vectors: NDArray[np.float64]) -> NDArray[np.float64]:
"""Length of each row vector."""
# via a typed local: numpy's stubs type einsum() as Any, which --strict rejects returning
lengths: NDArray[np.float64] = np.sqrt(np.einsum("ij,ij->i", vectors, vectors))
return lengths
def _face_array(faces: Sequence[Sequence[int]]) -> NDArray[np.intp] | None:
"""Return the faces as an (N, sides) index array, or None when they are not all the same length."""
if not faces:
return None
try:
arr = np.asarray(faces, dtype=np.intp)
except (TypeError, ValueError):
return None
return arr if arr.ndim == 2 else None
def _fan_triangles(faces: Sequence[Sequence[int]]) -> NDArray[np.intp]:
"""Every face fan-triangulated, as an (N, 3) array of vertex indices."""
arr = _face_array(faces)
if arr is not None and arr.shape[1] >= 3:
return np.concatenate(
[np.stack([arr[:, 0], arr[:, k], arr[:, k + 1]], axis=1) for k in range(1, arr.shape[1] - 1)]
)
tris = [(f[0], f[k], f[k + 1]) for f in faces for k in range(1, len(f) - 1)]
return np.asarray(tris, dtype=np.intp) if tris else np.empty((0, 3), dtype=np.intp)
VnfStyle = VNFStyle
[docs]
class VNF:
"""A VNF surface: ``vertices`` (3-D points) plus ``faces`` (index polygons into vertices).
Renders to PythonSCAD's native ``polyhedron`` via :meth:`polyhedron`. Build one from a
rectangular grid of sample points with :meth:`vertex_array`, merge several with
:meth:`union`, or mesh a scalar field with :meth:`from_field` and combine metaball
primitives with :meth:`from_metaballs`.
Args:
vertices: list of [x, y, z] points
faces: list of index lists (each polygon into *vertices*)
Examples:
Meshing a bumpy grid of sample points into a surface and rendering it as a polyhedron:
.. pythonscad-example::
import math
from pybosl2 import VNF
grid = [[[x, y, 4 * math.sin(x / 6) * math.cos(y / 6)] for y in range(0, 60, 4)]
for x in range(0, 60, 4)]
VNF.vertex_array(grid).polyhedron().show()
"""
def __init__(self, vertices: list[list[float]] | None = None, faces: list[list[int]] | None = None) -> None:
"""Initialize the VNF with vertices and faces."""
self.vertices = _rows_as_float(vertices)
self.faces = _rows_as_int(faces)
def __repr__(self) -> str:
"""Return a string representation of the VNF."""
return f"VNF({len(self.vertices)} verts, {len(self.faces)} faces)"
def __bool__(self) -> bool:
"""Check if the VNF has any faces."""
return len(self.faces) > 0
[docs]
def bounds(self) -> Bounds3D:
"""Axis-aligned :class:`~pybosl2.bounds.Bounds3D` of the VNF."""
arr = np.asarray(self.vertices, dtype=float)
mn, mx = arr.min(axis=0), arr.max(axis=0)
return Bounds3D(
min_x=float(mn[0]),
min_y=float(mn[1]),
min_z=float(mn[2]),
max_x=float(mx[0]),
max_y=float(mx[1]),
max_z=float(mx[2]),
width=float(mx[0] - mn[0]),
length=float(mx[1] - mn[1]),
height=float(mx[2] - mn[2]),
)
[docs]
def reverse(self) -> "VNF":
"""Return a copy with every face wound the other way (flips the surface normals)."""
return VNF(self.vertices, [f[::-1] for f in self.faces])
[docs]
def volume(self) -> float:
"""Signed enclosed volume (BOSL2 vnf_volume()); negative when the faces wind inward.
Used to detect and fix inverted meshes (a swept/skinned surface whose winding came out
inside-out): ``vnf if vnf.volume() >= 0 else vnf.reverse()``.
"""
tris = _fan_triangles(self.faces) # every face fan-triangulated, all at once
if not len(tris):
return 0.0
v = np.asarray(self.vertices, dtype=float)
a, b, c = v[tris[:, 0]], v[tris[:, 1]], v[tris[:, 2]]
return float(np.einsum("ij,ij->i", a, np.cross(b, c)).sum()) / 6.0
[docs]
@classmethod
def union(cls, vnfs: list["VNF"]) -> "VNF":
"""Merge a list of VNFs into one, offsetting each VNF's face indices (BOSL2 vnf_join())."""
vnfs = list(vnfs)
if len(vnfs) == 1:
return vnfs[0]
verts: list[Any] = []
faces: list[Any] = []
off = 0
for v in vnfs:
faces.extend([off + j for j in f] for f in v.faces if len(f) >= 3)
verts.extend(v.vertices)
off += len(v.vertices)
return cls(verts, faces)
[docs]
@classmethod
def join(cls, vnfs: list["VNF"]) -> "VNF":
"""Merge multiple VNFs into a single consolidated VNF with shared vertices.
Each input VNF's vertices and faces are copied into a combined vertex array,
with face indices offset appropriately. No deduplication is performed.
Args:
vnfs: A list of :class:`VNF` objects to merge.
Returns:
A new :class:`VNF` containing all vertices and faces from the inputs.
Examples:
.. pythonscad-example::
from pybosl2 import VNF
a = VNF.vertex_array([[ [0,0,0],[1,0,0] ], [ [0,1,0],[1,1,0] ]])
b = VNF.vertex_array([[ [0,0,1],[1,0,1] ], [ [0,1,1],[1,1,1] ]])
VNF.join([a, b]).polyhedron().show()
"""
return cls.union(vnfs)
[docs]
def halfspace(
self,
plane: Sequence[float],
keep: bool = True,
closed: bool = True,
) -> "VNF":
"""Clip a VNF to one side of a plane, optionally closing the cut face.
A plane is defined as ``[A, B, C, D]`` for ``A*x + B*y + C*z = D``.
If *keep* is True, the positive halfspace (``A*x + B*y + C*z > D``)
is retained. If *keep* is False, the negative halfspace is retained.
Args:
plane: Plane equation ``[A, B, C, D]``.
keep: If True, keep the positive halfspace. Defaults to True.
closed: If True, triangulate and close the cut face. Defaults to True.
Returns:
A new :class:`VNF` containing only the requested halfspace.
Raises:
AssertionError: If *plane* does not have exactly 4 elements.
Examples:
.. pythonscad-example::
import numpy as np
from pybosl2 import VNF, Bounds3D
cube_vnf = VNF.from_field(
lambda p: 5 - np.max(np.abs(p), axis=1),
0, Bounds3D(-10,-10,-10,10,10,10,20,20,20), voxel_size=1
)
cut = cube_vnf.halfspace([0, 0, 1, 0], keep=True, closed=True)
cut.polyhedron().show()
"""
assert len(plane) == 4, "halfspace(): plane must be [A, B, C, D]."
a, b, c, d = plane[0], plane[1], plane[2], plane[3]
verts_in = np.asarray(self.vertices, dtype=float)
if len(verts_in) == 0:
return VNF([], [])
n: np.ndarray = np.array([a, b, c], dtype=float)
dists: np.ndarray = verts_in @ n - d
if keep:
inside_mask: np.ndarray = dists >= -_EPS
else:
inside_mask = dists <= _EPS
inside_indices: list[int] = [i for i, m in enumerate(inside_mask) if m]
vertex_map: dict[int, int] = {}
for new_idx, old_idx in enumerate(inside_indices):
vertex_map[old_idx] = new_idx
new_verts: list[list[float]] = [list(verts_in[i]) for i in inside_indices]
new_faces: list[list[int]] = []
cut_edges: list[tuple[int, int]] = []
for face in self.faces:
face_inside: list[bool] = [inside_mask[i] for i in face]
all_in = all(face_inside)
none_in = not any(face_inside)
if all_in:
new_faces.append([vertex_map[i] for i in face])
elif not none_in:
fv = len(new_verts)
clipped: list[int] = []
nv = len(face)
for idx in range(nv):
i0 = face[idx]
i1 = face[(idx + 1) % nv]
v0_in = inside_mask[i0]
v1_in = inside_mask[i1]
if v0_in and v1_in:
if not clipped or clipped[-1] != vertex_map[i0]:
clipped.append(vertex_map[i0])
clipped.append(vertex_map[i1])
elif v0_in and not v1_in:
if not clipped or clipped[-1] != vertex_map[i0]:
clipped.append(vertex_map[i0])
t = _plane_edge_t(self.vertices[i0], self.vertices[i1], a, b, c, d)
pt = _interpolate(self.vertices[i0], self.vertices[i1], t)
new_verts.append(pt)
clipped.append(fv)
cut_edges.append((fv, fv + 1))
fv += 1
elif not v0_in and v1_in:
t = _plane_edge_t(self.vertices[i0], self.vertices[i1], a, b, c, d)
pt = _interpolate(self.vertices[i0], self.vertices[i1], t)
new_verts.append(pt)
clipped.append(fv)
cut_edges.append((fv, fv + 1))
fv += 1
clipped.append(vertex_map[i1])
if len(clipped) >= 3:
# fan-triangulate the clipped polygon
base = clipped[0]
for k in range(1, len(clipped) - 1):
tri = [base, clipped[k], clipped[k + 1]]
if _triangle_area(new_verts[tri[0]], new_verts[tri[1]], new_verts[tri[2]]) > _EPS:
new_faces.append(tri)
if closed and cut_edges:
edge_list = list(cut_edges)
paths: list[list[int]] = _assemble_edge_paths(edge_list)
for path in paths:
if len(path) >= 3:
pbase = path[0]
for k in range(1, len(path) - 1):
tri = [pbase, path[k], path[k + 1]]
if _triangle_area(new_verts[tri[0]], new_verts[tri[1]], new_verts[tri[2]]) > _EPS:
new_faces.append(tri)
return VNF(new_verts, new_faces)
[docs]
def slice(
self,
plane: Sequence[float],
closed: bool = True,
) -> tuple["VNF", "VNF"]:
"""Slice a VNF into two VNFs along a plane, closing both cut faces.
Returns ``(vnf_above, vnf_below)`` where *vnf_above* is the positive
halfspace and *vnf_below* is the negative halfspace.
Args:
plane: Plane equation ``[A, B, C, D]`` for ``A*x + B*y + C*z = D``.
closed: If True, close both cut faces. Defaults to True.
Returns:
A ``(above, below)`` tuple of :class:`VNF` objects.
Examples:
.. pythonscad-example::
import numpy as np
from pybosl2 import VNF, Bounds3D
cube_vnf = VNF.from_field(
lambda p: 5 - np.max(np.abs(p), axis=1),
0, Bounds3D(-10,-10,-10,10,10,10,20,20,20), voxel_size=1
)
above, below = cube_vnf.slice([0, 0, 1, 0], closed=True)
above.polyhedron().show()
"""
above = self.halfspace(plane, keep=True, closed=closed)
below = self.halfspace(plane, keep=False, closed=closed)
return above, below
[docs]
@classmethod
def vertex_array(
cls,
points: Path3D | list[Path3D] | list[list[list[float]]] | list[np.ndarray] | np.ndarray,
caps: "CapsSpec | None" = None,
col_wrap: bool = False,
row_wrap: bool = False,
reverse: bool = False,
style: VNFStyle | VnfStyle = VNFStyle.DEFAULT,
) -> "VNF":
"""Build a VNF from a rectangular grid of 3-D points (BOSL2 vnf_vertex_array()).
Each grid cell becomes triangles (or a quad) chosen by *style*: "default", "alt",
"min_edge", "min_area", "convex", "concave", "quincunx", "quad", "flip1", "flip2".
*col_wrap*/*row_wrap* close the grid into a tube/torus; *caps* closes the
column-wrapped ends with :class:`~pybosl2.caps.CapType` or
:class:`~pybosl2.caps.CapSpec` styles (see :data:`~pybosl2.caps.CapsSpec`).
*reverse* flips face winding. Degenerate (zero-area) faces are dropped.
Args:
points: Input grid points.
caps: Cap specification for both ends: a single :class:`CapType`,
:class:`CapSpec`, or a two-element pair ``[cap_start, cap_end]``.
Pass ``None`` (the default) for no caps.
col_wrap: Close the column direction into a tube.
row_wrap: Close the row direction into a torus.
reverse: Flip face winding.
style: Triangulation method.
"""
from pybosl2.caps import CapType, norm_caps
cap_specs: list[CapSpec] = norm_caps(caps) if caps is not None else norm_caps(CapType.NONE)
cap1, cap2 = cap_specs[0], cap_specs[1]
def _resolve_cap(cap: CapSpec) -> tuple[bool, bool]:
if cap.cap_type == CapType.NONE:
return False, False
if cap.cap_type in (CapType.ROUND, CapType.SPHERE):
return True, True
return True, False
grid = _to_grid(points)
rows = len(grid)
if rows == 0:
return cls([], [])
cols = len(grid[0])
if rows <= 1 or cols <= 1:
return cls([], [])
make_cap1, cap1_round = _resolve_cap(cap1)
make_cap2, cap2_round = _resolve_cap(cap2)
if (make_cap1 or make_cap2) and not col_wrap:
raise AssertionError("col_wrap must be true if caps are requested")
if (make_cap1 or make_cap2) and row_wrap:
raise AssertionError("cannot combine caps with row_wrap")
parr = np.asarray(grid, dtype=float).reshape(rows * cols, -1) # flattened, row-major
pcnt = rows * cols
colcnt = cols - (0 if col_wrap else 1)
rowcnt = rows - (0 if row_wrap else 1)
def idx(r: int, c: int) -> int:
return (r % rows) * cols + (c % cols)
# the four corner indices of every cell, in row-major cell order
cell_r = np.repeat(np.arange(rowcnt), colcnt)
cell_c = np.tile(np.arange(colcnt), rowcnt)
i1 = (cell_r % rows) * cols + (cell_c % cols)
i2 = ((cell_r + 1) % rows) * cols + (cell_c % cols)
i3 = ((cell_r + 1) % rows) * cols + ((cell_c + 1) % cols)
i4 = (cell_r % rows) * cols + ((cell_c + 1) % cols)
verts = parr.tolist()
if style == VNFStyle.QUINCUNX:
centres = (parr[i1] + parr[i2] + parr[i3] + parr[i4]) / 4.0
verts.extend(centres.tolist())
vertsarr = np.asarray(verts, dtype=float)
faces: list[Any] = []
if make_cap1:
if cap1_round:
row0 = parr[:cols]
center: list[float] = list(row0.mean(axis=0))
dome_radius: float = float(max(np.linalg.norm(p[:-1] - center[:-1]) for p in row0))
apex: list[float] = center.copy()
apex[2] -= dome_radius
apex_idx = len(verts)
verts.append(apex)
for i in range(cols):
j = (i + 1) % cols
faces.append([apex_idx, idx(0, i), idx(0, j)] if reverse else [idx(0, i), apex_idx, idx(0, j)])
else:
faces.append(_count(cols, 0, reverse=not reverse))
if make_cap2:
if cap2_round:
row_last = parr[(rows - 1) * cols : rows * cols]
center = list(row_last.mean(axis=0))
dome_radius = float(max(np.linalg.norm(p[:-1] - center[:-1]) for p in row_last))
apex = center.copy()
apex[2] += dome_radius
apex_idx = len(verts)
verts.append(apex)
for i in range(cols):
j = (i + 1) % cols
faces.append(
[apex_idx, idx(rows - 1, i), idx(rows - 1, j)]
if not reverse
else [idx(rows - 1, i), apex_idx, idx(rows - 1, j)]
)
else:
faces.append(_count(cols, (rows - 1) * cols, reverse=reverse))
cells = VNF._cell_faces(style, parr, i1, i2, i3, i4, cell_r, cell_c, pcnt, colcnt, reverse)
faces.extend(VNF._keep_real_faces(cells, vertsarr, reverse))
return cls(verts, faces)
@staticmethod
def _cell_faces(
style: VNFStyle,
parr: NDArray[np.float64],
i1: NDArray[np.intp],
i2: NDArray[np.intp],
i3: NDArray[np.intp],
i4: NDArray[np.intp],
cell_r: NDArray[np.intp],
cell_c: NDArray[np.intp],
pcnt: int,
colcnt: int,
reverse: bool,
) -> NDArray[np.intp]:
"""Which way each grid cell is cut, for every cell at once.
Returns a ``(cells, faces_per_cell, indices_per_face)`` array in the same order the cells
are walked -- row-major, and within a cell in the order the style lays the faces down.
Args:
style: The triangulation style.
parr: The grid points, flattened row-major.
i1: Index of each cell's first corner (and *i2*, *i3*, *i4* the rest, going round).
i2: Index of each cell's second corner.
i3: Index of each cell's third corner.
i4: Index of each cell's fourth corner.
cell_r: Grid row of each cell.
cell_c: Grid column of each cell.
pcnt: How many grid points there are (where the quincunx centres start).
colcnt: How many cells per row.
reverse: Whether the faces are being wound the other way.
Returns:
The face indices per cell.
"""
p1, p2, p3, p4 = parr[i1], parr[i2], parr[i3], parr[i4]
if style == VNFStyle.QUAD:
return np.stack([np.stack([i1, i2, i3, i4], axis=1)], axis=1)
if style == VNFStyle.QUINCUNX:
i5 = pcnt + cell_r * colcnt + cell_c
return np.stack(
[
np.stack([i1, i5, i2], axis=1),
np.stack([i2, i5, i3], axis=1),
np.stack([i3, i5, i4], axis=1),
np.stack([i4, i5, i1], axis=1),
],
axis=1,
)
alt = np.stack([np.stack([i1, i4, i2], axis=1), np.stack([i2, i4, i3], axis=1)], axis=1)
default = np.stack([np.stack([i1, i3, i2], axis=1), np.stack([i1, i4, i3], axis=1)], axis=1)
if style == VNFStyle.MIN_AREA:
area42 = _norms(np.cross(p2 - p1, p4 - p1)) + _norms(np.cross(p4 - p3, p2 - p3))
area13 = _norms(np.cross(p1 - p4, p3 - p4)) + _norms(np.cross(p3 - p2, p1 - p2))
use_alt = area42 < area13 + _EPS
elif style == VNFStyle.MIN_EDGE:
use_alt = _norms(p4 - p2) < _norms(p1 - p3) + _EPS
elif style in (VNFStyle.CONVEX, VNFStyle.CONCAVE):
sides = (-1 if reverse else 1) * np.cross(p2 - p1, p3 - p1)
dot4 = np.einsum("ij,ij->i", sides, p4)
dot1 = np.einsum("ij,ij->i", sides, p1)
use_alt = dot4 > dot1 if style == VNFStyle.CONVEX else dot4 <= dot1
flat = ~np.any(sides, axis=1) # a cell with no normal collapses to one triangle
if flat.any():
cells = np.where(use_alt[:, None, None], alt, default)
cells[flat, 0] = np.stack([i1, i4, i3], axis=1)[flat]
cells[flat, 1] = i1[flat][:, None] # a single repeated corner: dropped as degenerate
return cells
elif style == VNFStyle.ALT:
use_alt = np.ones(len(i1), dtype=bool)
elif style == VNFStyle.FLIP1:
use_alt = (cell_r + cell_c) % 2 == 0
elif style == VNFStyle.FLIP2:
use_alt = (cell_r + cell_c) % 2 == 1
else: # default
use_alt = np.zeros(len(i1), dtype=bool)
return np.where(use_alt[:, None, None], alt, default)
@staticmethod
def _keep_real_faces(cells: NDArray[np.intp], verts: NDArray[np.float64], reverse: bool) -> list[list[int]]:
"""Return the faces of *cells* that enclose an area, wound to match *reverse*.
Args:
cells: Face indices per cell, as :meth:`_cell_faces` returns them.
verts: The vertex positions the indices point into.
reverse: If True, each face comes back wound the other way.
Returns:
The faces, in cell order, with the degenerate ones dropped.
"""
flat = cells.reshape(-1, cells.shape[2])
a, b, c = verts[flat[:, 0]], verts[flat[:, 1]], verts[flat[:, 2]]
keep = _norms(np.cross(b - a, c - a)) > _EPS
kept = flat[keep]
return (kept[:, ::-1] if reverse else kept).tolist() # type: ignore[no-any-return]
[docs]
@classmethod
def tri_array(
cls,
points: list[list[list[float]]],
caps: bool = False,
cap1: bool | None = None,
cap2: bool | None = None,
col_wrap: bool = False,
row_wrap: bool = False,
reverse: bool = False,
limit_bunching: bool = True,
) -> "VNF":
"""Build a VNF from an array of rows whose lengths may differ (BOSL2 vnf_tri_array()).
Triangulates between adjacent rows by repeatedly adding the shortest new edge, so it
meshes triangular / irregular point arrays (what the degenerate bezier patches produce).
"""
if (caps or cap1 or cap2) and row_wrap:
raise AssertionError("cannot combine caps with row_wrap")
plen = len(points)
st = []
for row in points:
row = [list(p) for p in row]
if col_wrap and not np.array_equal(row[0], row[-1]):
row = row + [list(row[0])]
st.append(row)
addcol = (len(st[0]) - len(points[0])) if col_wrap else 0
rowstarts = [len(r) for r in st]
pcumlen = [0]
for n in rowstarts:
pcumlen.append(pcumlen[-1] + n)
capfirst = cap1 if cap1 is not None else (caps if caps is not None else False)
caplast = cap2 if cap2 is not None else (caps if caps is not None else False)
faces: list[Any] = []
if capfirst:
rng = list(range(rowstarts[0] - addcol)) if reverse else list(range(rowstarts[0] - 1 - addcol, -1, -1))
faces.append(rng)
for i in range(plen - 1 + (1 if row_wrap else 0)):
j = (i + 1) % plen
trimax = max(1, abs(len(st[i]) - len(st[j]))) if limit_bunching else float("inf")
faces.extend(
_lofttri(
st[i], # type: ignore[arg-type]
st[j], # type: ignore[arg-type]
pcumlen[i],
pcumlen[j],
rowstarts[i],
rowstarts[j],
reverse,
trimax,
)
)
if caplast:
if reverse:
rng = list(range(pcumlen[plen] - 1 - addcol, pcumlen[plen - 1] - 1, -1))
else:
rng = list(range(pcumlen[plen - 1], pcumlen[plen] - addcol))
faces.append(rng)
verts = [p for row in st for p in row]
return cls(verts, faces)
[docs]
def polyhedron(self) -> Any:
"""Native geometry for this VNF via PythonSCAD's ``polyhedron(points=, faces=)``.
A VNF winds its faces counter-clockwise seen from outside (so :meth:`volume` is positive
for a solid); ``polyhedron()`` wants them the other way round, so each face is reversed
on the way out. Handing them over as-is builds the solid inside out -- it still looks
right on its own, but every union or difference with it then does the opposite of what
it should.
"""
from pythonscad import polyhedron as _polyhedron
# The native polyhedron() rejects an empty point list with a bare "There must at least be
# one point in the polyhedron"; say which VNF that came from instead.
if not self.vertices or not self.faces:
raise ValueError(
f"polyhedron(): this VNF has no geometry to build "
f"({len(self.vertices)} vertices, {len(self.faces)} faces)."
)
return _polyhedron(points=self.vertices, faces=[f[::-1] for f in self.faces], convexity=10)
[docs]
def geometry(self) -> Any:
"""Return the VNF as native polyhedron geometry, matching Path2D/Region's geometry() surface."""
return self.polyhedron()
[docs]
@classmethod
def from_field(
cls,
f: np.ndarray | Path3D | Callable[[np.ndarray], np.ndarray] | Callable[[Path3D], np.ndarray],
isovalue: float,
bounding_box: Bounds3D | float | Sequence[float] | Sequence[Sequence[float]] | None = None,
voxel_size: float | None = None,
voxel_count: int | None = None,
closed: bool = True,
reverse: bool = False,
exact_bounds: bool = False,
) -> "VNF":
"""Mesh a scalar field into a :class:`VNF` via marching cubes.
The solid is the region where ``f >= isovalue``.
Args:
f: A :class:`~pybosl2.path3d.Path3D`, a 3-D numpy array,
a ``(N,3) → (N,)`` callable, or a
``(:class:`~pybosl2.path3d.Path3D`) → (N,)`` callable.
isovalue: Scalar threshold.
bounding_box: A :class:`~pybosl2.bounds.Bounds3D` or ``None``
(auto-computed from array shape when *f* is an array).
voxel_size: Isotropic voxel size.
voxel_count: Approximate total voxel count (ignored if *voxel_size* given).
closed: If True, pad field so mesh closes at bounding-box faces.
reverse: If True, reverse inside/outside sense.
exact_bounds: If True, use *bounding_box* exactly.
Returns:
A :class:`VNF`.
Raises:
NotImplementedError: If *isovalue* is a tuple range; only scalar thresholds are supported.
Examples:
.. pythonscad-example::
import numpy as np
from pybosl2 import VNF, Bounds3D
def field(p):
x, y, z = p[:, 0], p[:, 1], p[:, 2]
return 20 / np.sqrt(x*x + y*y + z*z) + 3 * np.sin(x / 3)
VNF.from_field(
field, 1,
Bounds3D(-30, -30, -30, 30, 30, 30, 60, 60, 60),
voxel_size=2,
).polyhedron().show()
"""
from pybosl2.path3d import Path3D
bb: Bounds3D | None = None
if bounding_box is not None:
if isinstance(bounding_box, Bounds3D):
bb = bounding_box
elif isinstance(bounding_box, (int, float)):
size = float(bounding_box)
bb = Bounds3D(-size / 2, -size / 2, -size / 2, size / 2, size / 2, size / 2, size, size, size)
elif isinstance(bounding_box, (list, tuple, np.ndarray)):
val_list = list(bounding_box)
if len(val_list) == 2 and isinstance(val_list[0], (list, tuple, np.ndarray)):
p1 = [float(x) for x in val_list[0]]
p2 = [float(x) for x in val_list[1]]
bb = Bounds3D(p1[0], p1[1], p1[2], p2[0], p2[1], p2[2], p2[0] - p1[0], p2[1] - p1[1], p2[2] - p1[2])
elif len(val_list) == 3:
val = [float(x) for x in val_list]
bb = Bounds3D(
-val[0] / 2,
-val[1] / 2,
-val[2] / 2,
val[0] / 2,
val[1] / 2,
val[2] / 2,
val[0],
val[1],
val[2],
)
elif len(val_list) == 6:
val = [float(x) for x in val_list]
bb = Bounds3D(
val[0],
val[1],
val[2],
val[3],
val[4],
val[5],
val[3] - val[0],
val[4] - val[1],
val[5] - val[2],
)
if isinstance(f, Path3D):
f = np.asarray(f, dtype=float)
elif callable(f) and not isinstance(f, np.ndarray):
_original = f
def _wrapped(pts: np.ndarray) -> np.ndarray:
try:
return np.asarray(_original(pts), dtype=float) # type: ignore[arg-type]
except (TypeError, ValueError):
return np.asarray(_original(Path3D(pts)), dtype=float) # type: ignore[arg-type]
f = _wrapped
if isinstance(isovalue, tuple):
raise NotImplementedError(
"from_field(): tuple (lo, hi) isovalue ranges are not yet implemented. "
"Use a single float isovalue instead."
)
iso = float(isovalue)
if isinstance(f, np.ndarray) or (isinstance(f, (list, tuple)) and not callable(f)):
field = np.asarray(f, dtype=float)
if bb is None:
vs = voxel_size if voxel_size is not None else 1.0
half = 0.5 * vs * (np.array(field.shape) - 1)
bb = Bounds3D(
min_x=-half[0],
min_y=-half[1],
min_z=-half[2],
max_x=half[0],
max_y=half[1],
max_z=half[2],
width=2 * half[0],
length=2 * half[1],
height=2 * half[2],
)
vs_final = vs
else:
vs_final = (bb.max_x - bb.min_x) / (field.shape[0] - 1)
xs, ys, zs = _grid_axes(bb, vs_final)
else:
assert bb is not None, "from_field(): a callable field needs a bounding_box."
bb, vs_final = _resolve_grid(bb, voxel_size, voxel_count, exact_bounds)
xs, ys, zs = _grid_axes(bb, vs_final)
field = _sample_field(f, xs, ys, zs)
verts, faces = _marching_cubes(field, xs, ys, zs, iso, closed)
vnf = cls(verts, faces)
if len(faces):
vol = vnf.volume()
if (vol < 0) != reverse:
vnf = vnf.reverse()
return vnf
[docs]
@classmethod
def from_skin(
cls,
profiles: Sequence[Sequence[Sequence[float]]],
slices: int,
refine: float = 1.0,
method: SkinMethod = SkinMethod.DIRECT,
sampling: SamplingType | None = None,
caps: "CapsSpec" = "butt",
closed: bool = False,
style: VNFStyle = VNFStyle.MIN_EDGE,
z: Sequence[float] | None = None,
) -> "VNF | Bosl2Solid":
"""Blend a stack of 2-D/3-D profiles into a skinned surface, returning a VNF or Bosl2Solid.
Consecutive profiles are connected vertex-to-vertex; *slices* extra interpolated profiles are
inserted between each pair to smooth the transition.
Args:
profiles: list of >= 2 closed profiles (each a list of points). If 2-D, give matching *z*.
slices: number of interpolated profiles inserted between each pair (int or per-gap list)
refine: subdivide every profile by this factor before skinning (default 1)
method: "direct" (connect vertex i to vertex i) or "reindex" (rotate each profile to
best-align with the previous).
sampling: "length" or "segment" resampling (default "length")
caps: cap the ends; supports decorative cap types
closed: the stack loops back to the first profile (default False)
style: vnf_vertex_array quad-subdivision style
z: per-profile Z heights, required when the profiles are 2-D
Examples:
Skinning a round profile up to a square one (a lofted transition):
.. pythonscad-example::
import math
import numpy as np
from pybosl2 import VNF
from pybosl2.enums import SkinMethod
circle = [[6 * math.cos(t), 6 * math.sin(t)] for t in np.linspace(0, 2 * math.pi, 24, endpoint=False)]
square = [[-8, -8], [8, -8], [8, 8], [-8, 8]]
VNF.from_skin([circle, square], slices=20, method=SkinMethod.REINDEX, z=[0, 25]).polyhedron().show()
"""
from pybosl2.skin import _skin
return _skin(
profiles,
slices,
refine=refine,
method=method,
sampling=sampling,
caps=caps,
closed=closed,
style=style,
z=z,
)
__all__ = [
"VNF",
"VnfStyle",
"contour",
]