openscad/BOSL2
1
0
Fork 0
mirror of https://github.com/revarbat/BOSL2.git synced 2025-08-18 18:11:16 +02:00
Code Issues Projects Releases Wiki Activity

Update isosurface.scad

Browse Source
This commit is contained in:
Richard Milewski
2025-02-09 13:38:10 -08:00
parent 6631e14cd3
commit 65d5471441
1 changed files with 431 additions and 201 deletions
Download Patch File Download Diff File Expand all files Collapse all files

632
isosurface.scad
View File

@@ -1,33 +1,39 @@
/////////////////////////////////////////////////////////////////////
// LibFile: isosurface.scad
// [metaballs](https://en.wikipedia.org/wiki/Metaballs) (also known as "blobby objects"),
// are bounded and closed organic-looking surfaces that smoothly meld together when in close proximity.
// Metaballs are a specific example of isosurfaces.
// are bounded and closed organic surfaces that smoothly blend together.
// Metaballs are one specific kind of isosurface.
// .
// An isosurface is a three-dimensional surface representing points of a constant
// value (e.g. pressure, temperature, electric potential, density) in a
// An isosurface, or implicit surface, is a three-dimensional surface representing all points of a
// constant value (e.g. pressure, temperature, electric potential, density) in a
// 3D volume. It's the 3D version of a 2D contour; in fact, any 2D cross-section of an
// isosurface *is* a 2D contour.
// isosurface **is** a 2D contour.
// .
// For computer-aided design, isosurfaces of abstract functions can generate complex
// curved surfaces and organic-looking shapes. Metaballs, for example, are typically generated by
// point sources (the metaball centers) that contribute values to all the points in a 3D volume depending
// on distance to the metaball centers. The combined contributions from all metaballs in the volume
// result variying values throughout the volume. Surfaces that connect all points of equal value in the
// volume that appear as blobby shapes that blend together.
// For computer-aided design, isosurfaces of abstract functions can generate complex curved surfaces
// and organic shapes. For example, spherical metaballs can be formulated using a set of point
// centers that define the metaballs locations. For a single metaball, a function is defined each
// each point in a 3D volume based on the distance from that point to the metaball center. The
// combined contributions from all the metaballs results in a function that varies in a complicated
// way across the volume. When two metaballs are far apart, they appear simply as spheres, but when
// they are close together they enlarge amd reach toward each other and meld together in a smooth
// fashion. The resulting metaball model appears as smoothly blended blobby shapes. The
// implementation below provides metaballs of a variety of types including spheres, cuboids and
// cylinders (cones), with optional parameters to adjust the influence of one metaball on others,
// and the cutoff distance where the metaball's influence stops.
// .
// In general temrs, an isosurface may be represented by any function of three variables,
// that is, the isosurface of a function $f(x,y,z)$ is the set of points where
// An isosurface can be defined using any function of three variables:
// the isosurface of a function $f(x,y,z)$ is the set of points where
// $f(x,y,z)=c$ for some constant value $c$. The constant $c$ is referred to as the "isovalue".
// We use the term "isovalue" in the context of metaballs also, although you may find other applications
// referring to "threshold", which means the same thing.
// Changing the isovalue will tend to grow or shrink the isosurface, depending on how the function is
// defined. Since metaballs are isosurfaces, they also have an isovalue. The isovalue is also known
// as the "threshold".
// .
// Some isosurface functions are unbounded, extending infinitely in all directions. A familiar example may
// be a [gryoid](https://en.wikipedia.org/wiki/Gyroid), which is often used as a volume infill pattern in
// [fused filament fabrication](https://en.wikipedia.org/wiki/Fused_filament_fabrication)). The gyroid
// isosurface is unbounded and periodic in all three dimensions.
// .
// Below are modules and functions to create manifold 3D models of metaballs and other isosurfaces.
// This file provides modules and functions to create a VNF using metaballs, or from general isosurfaces.
//
// Includes:
// include <BOSL2/std.scad>
@@ -682,9 +688,9 @@ function _isosurface_cubes(voxsize, bbox, fieldarray, fieldfunc, isovalmin, isov
field = is_def(fieldarray)
? fieldarray
: let(v = bbox[0], hv = 0.5*voxsize, b1 = bbox[1]+[hv,hv,hv]) [
for(x=[v[0]:voxsize:b1[0]]) [
for(y=[v[1]:voxsize:b1[1]]) [
for(z=[v[2]:voxsize:b1[2]])
for(x=[v.x:voxsize:b1.x]) [
for(y=[v.y:voxsize:b1.y]) [
for(z=[v.z:voxsize:b1.z])
fieldfunc(x,y,z)
]
]
@@ -765,7 +771,7 @@ function _isosurface_triangles(cubelist, cubesize, isovalmin, isovalmax, tritabl
vi0 = edge[0],
vi1 = edge[1],
denom = f[vi1] - f[vi0],
u = abs(denom)<0.0001 ? 0.5 : (isovalmin-f[vi0]) / denom
u = abs(denom)<0.0001 || denom==-INF ? 0.5 : (isovalmin-f[vi0]) / denom
) vcube[vi0] + u*(vcube[vi1]-vcube[vi0]) ],
// max surface
[ for(ei=epathmax) let(
@@ -773,7 +779,7 @@ function _isosurface_triangles(cubelist, cubesize, isovalmin, isovalmax, tritabl
vi0 = edge[0],
vi1 = edge[1],
denom = f[vi1] - f[vi0],
u = abs(denom)<0.0001 ? 0.5 : (isovalmax-f[vi0]) / denom
u = abs(denom)<0.0001 || denom==-INF ? 0.5 : (isovalmax-f[vi0]) / denom
) vcube[vi0] + u*(vcube[vi1]-vcube[vi0]) ], outfacevertices)
) for(ls = list) ls
else if(n_outer>0 && lenmin>0) let(
@@ -785,7 +791,7 @@ function _isosurface_triangles(cubelist, cubesize, isovalmin, isovalmax, tritabl
vi0 = edge[0],
vi1 = edge[1],
denom = f[vi1] - f[vi0],
u = abs(denom)<0.0001 ? 0.5 : (isovalmin-f[vi0]) / denom
u = abs(denom)<0.0001 || denom==-INF ? 0.5 : (isovalmin-f[vi0]) / denom
) vcube[vi0] + u*(vcube[vi1]-vcube[vi0]) ], outfacevertices)
) for(ls = list) ls
else if(lenmin>0)
@@ -796,7 +802,7 @@ function _isosurface_triangles(cubelist, cubesize, isovalmin, isovalmax, tritabl
vi0 = edge[0],
vi1 = edge[1],
denom = f[vi1] - f[vi0],
u = abs(denom)<0.0001 ? 0.5 : (isovalmin-f[vi0]) / denom
u = abs(denom)<0.0001 || denom==-INF ? 0.5 : (isovalmin-f[vi0]) / denom
) vcube[vi0] + u*(vcube[vi1]-vcube[vi0])
else if(n_outer>0 && lenmax>0) let(
@@ -807,7 +813,7 @@ function _isosurface_triangles(cubelist, cubesize, isovalmin, isovalmax, tritabl
vi0 = edge[0],
vi1 = edge[1],
denom = f[vi1] - f[vi0],
u = abs(denom)<0.0001 ? 0.5 : (isovalmax-f[vi0]) / denom
u = abs(denom)<0.0001 || denom==-INF ? 0.5 : (isovalmax-f[vi0]) / denom
) vcube[vi0] + u*(vcube[vi1]-vcube[vi0]) ], outfacevertices)
) for(ls = list) ls
else if(lenmax>0)
@@ -818,7 +824,7 @@ function _isosurface_triangles(cubelist, cubesize, isovalmin, isovalmax, tritabl
vi0 = edge[0],
vi1 = edge[1],
denom = f[vi1] - f[vi0],
u = abs(denom)<0.0001 ? 0.5 : (isovalmax-f[vi0]) / denom
u = abs(denom)<0.0001 || denom==-INF ? 0.5 : (isovalmax-f[vi0]) / denom
) vcube[vi0] + u*(vcube[vi1]-vcube[vi0])
else if(n_outer>0)
@@ -843,7 +849,7 @@ function _bbfacevertices(vcube, f, bbface, isovalmax, isovalmin) = let(
fmax = max(f0, f1),
fbetweenlow = (fmin <= isovalmin && isovalmin <= fmax),
fbetweenhigh = (fmin <= isovalmax && isovalmax <= fmax),
denom = f1 - f0
denom = f1==INF || f0==INF ? 0 : f1-f0
) [
if(isovalmin <= f0 && f0 <= isovalmax) vcube[vi0],
if(fbetweenlow && f0<=f1) let(
@@ -935,7 +941,7 @@ function _mb_sphere_full(dv, r, ex, cutoff, neg) = let(dist=norm(dv))
function mb_sphere(r, cutoff=INF, influence=1, negative=false, d) =
assert(is_num(cutoff) && cutoff>0, "\ncutoff must be a positive number.")
assert(is_num(influence) && influence>0, "\ninfluence must be a positive number.")
assert(is_finite(influence) && influence>0, "\ninfluence must be a positive number.")
let(
r = get_radius(r=r,d=d),
dummy=assert(is_finite(r) && r>0, "\ninvalid radius or diameter."),
@@ -946,71 +952,150 @@ function mb_sphere(r, cutoff=INF, influence=1, negative=false, d) =
: influence==1 ? function(dv) _mb_sphere_cutoff(dv,r,cutoff,neg)
: function(dv) _mb_sphere_full(dv,r,1/influence,cutoff,neg);
// metaball cylinder
// metaball capsule (round-ended cylinder)
function _mb_cylinder_basic(dv, r, hl, neg) = let(
function _mb_capsule_basic(dv, hl, r, neg) = let(
dist = dv.z<-hl ? norm(dv-[0,0,-hl])
: dv.z<hl ? norm([dv.x,dv.y]) : norm(dv-[0,0,hl])
: dv.z<=hl ? norm([dv.x,dv.y]) : norm(dv-[0,0,hl])
) neg*r/dist;
function _mb_cylinder_influence(dv, r, hl, ex, neg) = let(
function _mb_capsule_influence(dv, hl, r, ex, neg) = let(
dist = dv.z<-hl ? norm(dv-[0,0,-hl])
: dv.z<hl ? norm([dv.x,dv.y]) : norm(dv-[0,0,hl])
: dv.z<=hl ? norm([dv.x,dv.y]) : norm(dv-[0,0,hl])
) neg * (r/dist)^ex;
function _mb_cylinder_cutoff(dv, r, hl, cutoff, neg) = let(
function _mb_capsule_cutoff(dv, hl, r, cutoff, neg) = let(
dist = dv.z<-hl ? norm(dv-[0,0,-hl])
: dv.z<hl ? norm([dv.x,dv.y]) : norm(dv-[0,0,hl])
) neg * mb_cutoff(dist, cutoff) * r/dist;
function _mb_cylinder_full(dv, r, hl, ex, cutoff, neg) = let(
function _mb_capsule_full(dv, hl, r, ex, cutoff, neg) = let(
dist = dv.z<-hl ? norm(dv-[0,0,-hl])
: dv.z<hl ? norm([dv.x,dv.y]) : norm(dv-[0,0,hl])
) neg * mb_cutoff(dist, cutoff) * (r/dist)^ex;
function mb_cylinder(r, length, cutoff=INF, influence=1, negative=false, d) =
assert(is_num(cutoff) && cutoff>0, "\ncutoff must be a positive number.")
assert(is_num(influence) && influence>0, "\ninfluence must be a positive number.")
assert(is_num(length) && length>0, "\nlength must be a positive number.")
let(
r = get_radius(r=r,d=d),
dummy=assert(is_finite(r) && r>0, "\ninvalid radius or diameter."),
neg = negative ? -1 : 1
function mb_capsule(h, r, cutoff=INF, influence=1, negative=false, d,l,height,length) =
assert(is_num(cutoff) && cutoff>0, "\ncutoff must be a positive number.")
assert(is_finite(influence) && influence>0, "\ninfluence must be a positive number.")
let(
h = one_defined([h,l,height,length],"h,l,height,length"),
dum1 = assert(is_finite(h) && h>0, "\ncylinder height must be a positive number."),
r = get_radius(r=r,d=d),
dum2 = assert(is_finite(r) && r>0, "\ninvalid radius or diameter."),
sh = h-2*r, // straight side length
dum3 = assert(sh>0, "\nTotal length must accommodate rounded ends of cylinder."),
neg = negative ? -1 : 1
)
!is_finite(cutoff) && influence==1 ? function(dv) _mb_cylinder_basic(dv,r,length/2,neg)
: !is_finite(cutoff) ? function(dv) _mb_cylinder_influence(dv,r,length/2,1/influence, neg)
: influence==1 ? function(dv)_mb_cylinder_cutoff(dv,r,length/2,cutoff,neg)
: function (dv) _mb_cylinder_full(dv, r, length/2, 1/influence, cutoff, neg);
!is_finite(cutoff) && influence==1 ? function(dv) _mb_capsule_basic(dv,sh/2,r,neg)
: !is_finite(cutoff) ? function(dv) _mb_capsule_influence(dv,sh/2,r,1/influence, neg)
: influence==1 ? function(dv) _mb_capsule_cutoff(dv,sh/2,r,cutoff,neg)
: function (dv) _mb_capsule_full(dv, sh/2, r, 1/influence, cutoff, neg);
// metaball connector cylinder - calls mb_capsule* functions after transform
function mb_connector(p1, p2, r, cutoff=INF, influence=1, negative=false, d) =
assert(is_num(cutoff) && cutoff>0, "\ncutoff must be a positive number.")
assert(is_finite(influence) && influence>0, "\ninfluence must be a positive number.")
let(
dum1 = assert(is_matrix([p1],1,3), "\nConnector start point p1 must be a 3D coordinate."),
dum2 = assert(is_matrix([p2],1,3), "\nConnector end point p2 must be a 3D coordinate."),
r = get_radius(r=r,d=d),
dum3 = assert(is_finite(r) && r>0, "\ninvalid radius or diameter."),
neg = negative ? -1 : 1,
dc = p2-p1, // center-to-center distance
midpt = reverse(-0.5*(p1+p2)),
h = norm(dc)/2, // center-to-center length (cylinder height)
transform = submatrix(down(h)*rot(from=dc,to=UP)*move(-p1) ,[0:2], [0:3])
)
!is_finite(cutoff) && influence==1 ? function(dv)
let(newdv = transform * [each dv,1])
_mb_capsule_basic(newdv,h,r,neg)
: !is_finite(cutoff) ? function(dv)
let(newdv = transform * [each dv,1])
_mb_capsule_influence(newdv,h,r,1/influence, neg)
: influence==1 ? function(dv)
let(newdv = transform * [each dv,1])
_mb_capsule_cutoff(newdv,h,r,cutoff,neg)
: function (dv)
let(newdv = transform * [each dv,1])
_mb_capsule_full(newdv, h, r, 1/influence, cutoff, neg);
// metaball rounded cylinder / cone
function revsurf(dv, path, coef, exp, cutoff) =
let(
pt = [norm([dv.x,dv.y]), dv.z],
segs = pair(path),
dist = min([for(seg=segs)
let(
c=seg[1]-seg[0],
s0 = seg[0]-pt,
t = -s0*c/(c*c)
)
t<0 ? norm(s0)
: t>1 ? norm(seg[1]-pt)
: norm(s0+t*c)]),
inside_check = [for(seg=segs)
if (cross(seg[1]-seg[0], pt-seg[0]) > EPSILON) 1]
)
mb_cutoff(dist, cutoff) * (inside_check==[]
? (coef*(1+dist))^exp : (coef/(1+dist))^exp );
function revsurf_rounded(path, rounding, cutoff, influence) =
let(shifted = offset(path, delta=-rounding, closed=false))
function (dv) revsurf(dv, shifted, 1+rounding, 1/influence, cutoff);
function mb_cyl(h,r,rounding=0,r1,r2,l,height,length,d1,d2,d, cutoff=INF, influence=1, negative=false) =
let(
r1 = get_radius(r1=r1,r=r, d1=d1, d=d),
r2 = get_radius(r1=r2,r=r, d1=d2, d=d),
h = first_defined([h,l,height,length],"h,l,height,length")
)
assert(is_finite(rounding) && rounding>=0, "rounding must be a nonnegative number")
assert(is_finite(r1) && r1>0, "r/r1/d/d1 must be a positive number")
assert(is_finite(r2) && r2>0, "r/r2/d/d2 must be a positive number")
let(
vang = atan2(r1-r2,h),
facelen = adj_ang_to_hyp(h, abs(vang)),
roundlen1 = rounding/tan(45-vang/2),
roundlen2 = rounding/tan(45+vang/2),
sides = [[0,h/2],[r2,h/2],[r1,-h/2], [0, -h/2]],
)
assert(roundlen1 <= r1, "size of rounding is larger than the r1 radius of the cylinder/cone")
assert(roundlen2 <= r2, "size of rounding is larger than the r2 radius of the cylinder/cone")
assert(roundlen1+roundlen2 <=facelen, "Roundings don't fit on the edge length of the cylinder/cone")
revsurf_rounded(sides, rounding, influence);
// metaball rounded cube
function _mb_roundcube_basic(dv, siz, xp, neg) = let(
function _mb_cuboid_basic(dv, siz, xp, neg) = let(
dist = xp >= 1100 ? max(v_abs(dv))
: (abs(dv.x)^xp + abs(dv.y)^xp + abs(dv.z)^xp) ^ (1/xp)
) neg*siz/dist;
function _mb_roundcube_influence(dv, siz, xp, ex, neg) = let(
function _mb_cuboid_influence(dv, siz, xp, ex, neg) = let(
dist = xp >= 1100 ? max(v_abs(dv))
:(abs(dv.x)^xp + abs(dv.y)^xp + abs(dv.z)^xp) ^ (1/xp)
) neg * (siz/dist)^ex;
function _mb_roundcube_cutoff(dv, siz, xp, cutoff, neg) = let(
function _mb_cuboid_cutoff(dv, siz, xp, cutoff, neg) = let(
dist = xp >= 1100 ? max(v_abs(dv))
: (abs(dv.x)^xp + abs(dv.y)^xp + abs(dv.z)^xp) ^ (1/xp)
) neg * mb_cutoff(dist, cutoff) * siz/dist;
function _mb_roundcube_full(dv, siz, xp, ex, cutoff, neg) = let(
function _mb_cuboid_full(dv, siz, xp, ex, cutoff, neg) = let(
dist = xp >= 1100 ? max(v_abs(dv))
:(abs(dv.x)^xp + abs(dv.y)^xp + abs(dv.z)^xp) ^ (1/xp)
) neg * mb_cutoff(dist, cutoff) * (siz/dist)^ex;
function mb_roundcube(size, squareness=0.5, cutoff=INF, influence=1, negative=false) =
function mb_cuboid(size, squareness=0.5, cutoff=INF, influence=1, negative=false) =
assert(is_num(cutoff) && cutoff>0, "\ncutoff must be a positive number.")
assert(is_num(influence) && influence>0, "\ninfluence must be a positive number.")
assert(is_num(size) && size>0, "\nsize must be a positive number.")
assert(is_finite(influence) && influence>0, "\ninfluence must be a positive number.")
assert(is_finite(size) && size>0, "\nsize must be a positive number.")
let(
dummy=assert(is_num(size) && size>0, "\ninvalid size."),
xp = _squircle_se_exponent(squareness),
neg = negative ? -1 : 1
)
!is_finite(cutoff) && influence==1 ? function(dv) _mb_roundcube_basic(dv, size/2, xp, neg)
: !is_finite(cutoff) ? function(dv) _mb_roundcube_influence(dv, size/2, xp, 1/influence, neg)
: influence==1 ? function(dv) _mb_roundcube_cutoff(dv, size/2, xp, cutoff, neg)
: function (dv) _mb_roundcube_full(dv, size/2, xp, 1/influence, cutoff, neg);
!is_finite(cutoff) && influence==1 ? function(dv) _mb_cuboid_basic(dv, size/2, xp, neg)
: !is_finite(cutoff) ? function(dv) _mb_cuboid_influence(dv, size/2, xp, 1/influence, neg)
: influence==1 ? function(dv) _mb_cuboid_cutoff(dv, size/2, xp, cutoff, neg)
: function (dv) _mb_cuboid_full(dv, size/2, xp, 1/influence, cutoff, neg);
// metaball octahedron
@@ -1025,7 +1110,7 @@ function _mb_octahedron_full(dv, r, ex, cutoff, neg) =
function mb_octahedron(r, cutoff=INF, influence=1, negative=false, d) =
assert(is_num(cutoff) && cutoff>0, "\ncutoff must be a positive number.")
assert(is_num(influence) && influence>0, "\ninfluence must be a positive number.")
assert(is_finite(influence) && is_num(influence) && influence>0, "\ninfluence must be a positive number.")
let(
r = get_radius(r=r,d=d),
dummy=assert(is_finite(r) && r>0, "\ninvalid radius or diameter."),
@@ -1049,14 +1134,23 @@ function _mb_torus_full(dv, rmaj, rmin, ex, cutoff, neg) =
let(dist = norm([norm([dv.x,dv.y])-rmaj, dv.z]))
neg * mb_cutoff(dist, cutoff) * (rmin/dist)^ex;
function mb_torus(r_maj, r_min, cutoff=INF, influence=1, negative=false, d_maj, d_min) =
function mb_torus(r_maj, r_min, cutoff=INF, influence=1, negative=false, d_maj, d_min, or,od,ir,id) =
assert(is_num(cutoff) && cutoff>0, "\ncutoff must be a positive number.")
assert(is_num(influence) && influence>0, "\ninfluence must be a positive number.")
assert(is_finite(influence) && influence>0, "\ninfluence must be a positive number.")
let(
r_maj = get_radius(r=r_maj,d=d_maj),
r_min = get_radius(r=r_min,d=d_min),
dum1=assert(is_finite(r_maj) && r_maj>0, "\ninvalid major radius or diameter."),
dum2=assert(is_finite(r_min) && r_min>0, "\ninvalid minor radius or diameter."),
_or = get_radius(r=or, d=od, dflt=undef),
_ir = get_radius(r=ir, d=id, dflt=undef),
_r_maj = get_radius(r=r_maj, d=d_maj, dflt=undef),
_r_min = get_radius(r=r_min, d=d_min, dflt=undef),
r_maj = is_finite(_r_maj)? _r_maj :
is_finite(_ir) && is_finite(_or)? (_or + _ir)/2 :
is_finite(_ir) && is_finite(_r_min)? (_ir + _r_min) :
is_finite(_or) && is_finite(_r_min)? (_or - _r_min) :
assert(false, "Bad major size parameter."),
r_min = is_finite(_r_min)? _r_min :
is_finite(_ir)? (maj_rad - _ir) :
is_finite(_or)? (_or - maj_rad) :
assert(false, "\nBad minor size parameter."),
neg = negative ? -1 : 1
)
!is_finite(cutoff) && influence==1 ? function(dv) _mb_torus_basic(dv, r_maj, r_min, neg)
@@ -1067,18 +1161,18 @@ function mb_torus(r_maj, r_min, cutoff=INF, influence=1, negative=false, d_maj,
/* hard-edge cylinder
_mb_cylinder = function (dv, coeff, length, cutoff, influence)
_mb_capsule = function (dv, coeff, length, cutoff, influence)
let(
dist = max(abs(dv.z*2*coeff/length), norm([dv.x,dv.y])),
suppress = let(a = min(r,cutoff)/cutoff) 1-a*a,
) suppress*coeff / dist;
function mb_cylinder(coeff=10, length=10, cutoff=INF, influence=1) = function (dv) _mb_cylinder(dv, coeff, length, cutoff, influence);
function mb_capsule(coeff=10, length=10, cutoff=INF, influence=1) = function (dv) _mb_capsule(dv, coeff, length, cutoff, influence);
*/
// Function&Module: metaballs()
// Synopsis: Creates a model of metaballs within a bounding box.
// Synopsis: Creates a group of 3D metaballs (smoothly connected blobs).
// SynTags: Geom,VNF
// Topics: Metaballs, Isosurfaces, VNF Generators
// See Also: isosurface_array()
@@ -1087,80 +1181,104 @@ function mb_cylinder(coeff=10, length=10, cutoff=INF, influence=1) = function (d
// Usage: As a function
// vnf = metaballs(funcs, bounding_box, voxel_size, [isovalue=], [closed=], [show_stats=]);
// Description:
// ![Metaball andimation](./images/metaball_demo.png)
// .
// [Metaballs](https://en.wikipedia.org/wiki/Metaballs), also known as "blobby objects",
// are organic-looking ball-shaped blobs that meld together when in close proximity.
// The melding property is determined by an interaction formula based on the coefficient
// weight (which can be thought of as a charge, strength, density, or intensity) of
// each ball and their distances from one another.
// can produce soothly varying blobs and organic forms. You create metaballs by placing metaball
// objects at different locations. These objects have a basic size and shape when placed in
// isolation, but if another metaball object is nearby, the two objects interact, growing larger
// and melding together. The closer the objects are, the more they blend and meld.
// .
// One analagous way to think of metaballs is, consider each "ball" to be a point-light source in
// a dark room. Pick an illumination value, and every point in the volume of the room with
// that intensity of illumination defines the isosurface, which would be a sphere around a
// single source, or a blob surrounding multiple points because the illumination is additive between them.
// The simplest metaball specification is a list of alternating transformation matrices and
// metaball functions: `[trans0, func0, trans1, func1, ... ]`. The transformation matrix
// you supply can be constructed using the usual transformation commands such as {{up()}},
// {{right()}}, {{back()}}, {{move()}}, {{scale()}}, {{rot()}} and so on. You can multiply
// the transformations together, similar to how the transformations can be applied
// to regular objects in OpenSCAD. For exmaple, to transform an object in regular OpenSCAD you might
// write `up(5) xrot(25) zrot(45) scale(4)`. You would provide that transformation
// as the transformation matrix `up(5) * xrot(25) * zrot(45) * scale(4)`. You can use
// scaling to produce an ellipse from a sphere, and you can even use {{skew()}} if desired.
// When no transformation is needed, give `IDENT` as the transformation.
// .
// Regardless of how you think of it (charge, light, heat, pressure), a stronger metaball
// intensity results in stronger "field" values around the metaball, and correspondingly a
// larger metaball due to the isosurface of a particular value being farther away.
// A metaball is basically a contour surface; that is, a 3D version of a 2D contour.
// You can create metaballs in a variety of standard shapes using the predefined functions
// listed below. If you wish, you can also create custom metaball shapes using your own functions
// (see Example 11). Three parameters are available on all of the built-in metaballs to control the
// interaction of the metaballs with each other: `cutoff`, `influence`, and `negative'.
// .
// Most implementations of metaballs instead use a simple inverse relationship proportional to $1/d$ to
// control how the contributions from each metaball fall off with distance. That
// is the default falloff used for the field types available here. The optional `influence`
// argument is a reciprocal exponent on $d$, defaulting to 1. It controls how much the metaball influences
// others at a given distance. If you set `influence=0.5`, the reciprocal is 2, so you get a $1/d^2$ falloff.
// The `cutoff` parameter specifies the distance beyond which the metaball has no interaction with
// other balls. When you apply `cutoff`, a smooth suppression factor begins begins decreasing the
// interaction strength at half the cutoff distance and reduces the interaction to zero at the cutoff.
// .
// You can also define your own metaball functions as shown in example 5.
// The `influence` parameter adjusts the strength of the interaction metaball objects have with each
// other. If you increase `influence` from its default of 1, the metaball interacts with other
// metaballs at a longer range, and surrounding balls grow bigger. The metaball with larger
// influence can also grow bigger because it couples more strongly with other nearby balls, but it
// can also remain nearly unchanged while influencing others when `isovalue` is greater than 1.
// Decreasing influence has the reverse effect. Small changes in influence can have a large
// effect. For example, setting `influence=2` dramatically increases the interactions at longer
// distances, and you may want to set the `cutoff` argument to limit the range influence.
// .
// The `negative` parameter, if set to `true`, creates a negative metaball, which can create
// hollows or dents in other metaballs, or swallow other metaballs entirely, making them disappear.
// Negative metaballs are always below the isovalue, so they are never directly visible;
// only their effects are visible. See Example 7.
// .
// For complicated metaball assemblies you may wish to repeat a structure in different locations or
// otherwise transformed. Nesting metaball specifications are supported:
// Instad of specifying a transform and function, you specify a transform and then a metaball
// specification. For example, you could set `finger=[t0,f0,t1,f1,t2,f2]` and then set
// `hand=[u0,finger,u1,finger,...]` and then invoke metaball with `[s0, hand]`.
// Basically, any list of metaballs is, itself, a metaball. It can be used in place of a function in another list.
// This is a powerful technique that lets you make groups of metaballs that you can use as individual
// metaballs in other groups, and can make your code compact and simpler to understand. See Example 13.
// .
// .h3 Built-in metaball functions
// Various shapes of metaball field density functions are built into this library. You can specify different
// ones for each metaball in the list, and you can also specify your own custom function.
// Several metaballs are defined for you to use in your models.
// All of the built-in metaballs take positional and named parameters that specify the size of the
// metaball (e.g. radius, height). The size arguments are the same as those for the regular objects
// of the same type (e.g. a sphere accepts both `r` for radius and the named parameter `d=` for
// diameter). The size parameters always specify the size of the metaball **in isolation** with
// `isovalue=1`. The metaballs can grow much bigger than their specified sizes when they interact
// with each other. The metaballs also grow bigger than their specified sizes, even in isolation,
// if `isovalue<1` and smaller than their specified sizes if `isovalue>1`.
// .
// All built-in functions have these arguments in common:
// * The first argument is always a size (such as a radius or diameter) that determines the size of
// the metaball. Other size arguments may follow as appropriate. All the
// metaball functions are designed so that an isolated metaball with `isovalue=1` appears with a
// radius or size approximately equal to this coefficient, but the metaball can get
// significantly larger when other metaballs are in the bounding box, depending on proximity.
// * `cutoff` - specifies the distance beyond which the metaball has no influence.
// A smooth suppression factor is applied to the metaball's influence on others, starting at half
// the cutoff distance, suppressing the influence to zero at the cutoff distance. Default: INF
// * `influence` - determines the extent of influence of the metaball. This is an inverse
// distance relationship proportional to $1/d$ where $d$ is distance. The `influence` argument is the
// reciprocal of the exponent; for example, If `influence=0.5`, you get an inverse-square falloff
// $1/d^2$, resulting in less influence at a given distance than the default `influence=1`. Setting
// `influence=2` results in a gentle $1/\sqrt d$ falloff, dramatically increasing
// the influence at distances, and you may want to set the `cutoff` argument to limit that influence.
// * `negative` - when true, causes the metaball to have a negative influence on its surroundings. A
// negative metaball can create hollows or dents in other metaballs, or swallow other metaballs
// entirely, making them disappear if the metaball's negative influence is large. A negative
// metaball is never visible directly, only its effect is visible, because the isosurface surrounds
// only field values greater than the isovalue (see Example 2 below). Default: false
// All of the built-in functions all accept these named arguments, which are not repeated in the list below:
// * `cutoff` - specifies the distance beyond which the metaball has no influence. Default: INF
// * `influence` - a positive number that determines how much influence the metaball has on others. Default: 1
// * `negative` - when true, causes the metaball to have a negative influence on its surroundings. Default: false
// .
// These are the built-in metaball functions. Arguments with default values are optional:
// * `mb_sphere(r/d, cutoff=INF, influence=1, negative=false)` - the standard spherical metaball with a $1/d$ falloff when `influence=1`. The `r` or `d` argument controls the radius or diameter, respectively. For a spherical metaball by itself, you get a sphere of radius `r` at `isovalue=1`. You can create an ellipsoid using `scale()` as the last transformation element in the metaball `funcs` array.
// * `mb_cylinder(r/d, length, cutoff=INF, influence=1, negative=false)` - a cylindrical-shaped field with rounded ends of radius `r` or diameter `d`, useful as a connector. For a single cylindrical metaball by itself at `isovalue=1`, you get a cylinder of radius `r` and straight-side length of `length`, but it grows when other metaballs are nearby.
// * `mb_roundcube(size, squareness=0.5, cutoff=INF, influence=1, negative=false)` - a cuboid metaball with rounded corners that get more rounded at farther distances, depending on isovalue and influence from other metaballs. The corner sharpness is controlled by the `squareness` parameter ranging from 0 (spherical) to 1 (cubical), and defaults to 0.5 if omitted. By itself with `isovalue=1`, you get a rounded cube having `size` distance from side to side.
// * `mb_octahedron(r/d, cutoff=INF, influence=1, negative=false)` - an octahedron-shaped metaball with sharp edges and corners, resulting from using [taxicab distance](https://en.wikipedia.org/wiki/Taxicab_geometry) rather than Euclidean distance calculations. By itself with `isovalue=1` you get an octahedron with tip radius of `r` or tip-to-top diameter `d`.
// * `mb_torus(r_maj/d_maj, r_min/d_min, cutoff=INF, influence=1, negative=false)` - a toroidal field oriented perpendicular to the z axis. The arguments `r_maj` and `r_min` control the major and minor radii; otherwise `d_maj` and `d_min` control the major and minor diameters.
// The built-in metaball functions are:
// * `mb_sphere(r|d=)` - spherical metaball, with radius r or diameter d. You can create an ellipsoid using `scale()` as the last transformation entry of the metaball `funcs` array.
// * `mb_cuboid(size, [squareness=0.5])` - cuboid metaball with rounded edges and corners. The corner sharpness is controlled by the `squareness` parameter ranging from 0 (spherical) to 1 (cubical), and defaults to 0.5. The `size` specifies the dimensions of the cuboid shape between the face centers, but except when `squareness=1`, the faces are always a little bit curved.
// * `mb_cyl(h|l|height|length=, [r|d=], [r1=|d1=], [r2=|d2=], [rounding=0])` - vertical cylinder or cone metaball with the same dimenional arguments as {{cyl()}}. The rounding is the same at both ends. If you just want a cylindrical shape, consider `mb_capsule()`, which has a faster execution time.
// * `mb_capsule(h|l|height|length=, r|d=)` - vertical cylinder of radius `r` or diameter `d` with hemispherical caps. The height or length specifies the **total** height including the rounded ends.
// * `mb_connector(p1, p2, r|d=)` - cylinder of radius `r` or diameter `d` with hemispherical caps (like mb_capsule), but specified to connect point `p1` to point `p2` (where `p1` and `p2` are 3D vectors). The specified points are at the ends of the straight portion of the shape (the centers of the two capping hemispheres).
// * `mb_octahedron(r|d=])` - octahedral metaball with sharp edges and corners. The `r` parameter specifies the distance from center to tip. The vertex parameter specifies the distance between the two opposite tips.
// * `mb_torus([r_maj|d_maj=], [r_min|d_min=],[or=|od=],[ir=|id=])` - torus metaball oriented perpendicular to the z axis. You can specify the torus dimensions using the same arguments as {{torus()}}; that is, major radius (or diameter) with `r_maj` or `d_maj`, and minor radius and diameter using `r_min` or `d_min`. Alternatively you can give the inner radius or diameter with `ir` or `id` and the outer radius or diameter with `or` or `od`.
// .
// Your own custom function must be written as a [function literal](https://en.wikibooks.org/wiki/OpenSCAD_User_Manual/User-Defined_Functions_and_Modules#Function_literals)
// and take `dv` as the first argument with a size as the second argument. `dv` is passed to your
// function as a 3D distance vector from the ball center to the point in the bounding box volume for
// which to calculate the field intensity. The function must return a single number such that higher
// values are enclosed by the metaball and lower values are outside the metaball.
// In this case, if you have written `my_func()`, the array element you initialize must appear
// as `function (dv) my_func(dv, ...)`. See Example 5 below.
// .h3 Metaball functions and user defined functions
// Each metaball function is defined as a function of a 3-vector that gives the value of the metaball function
// for that point in space. As is common in metaball implementations, we define the built-in metaballs using an
// inverse relationship where the metaball functions fall off as $1/d$. The spherical metaball therefore has
// a simple basic definition as `f(v) = 1/norm(v)`. Note that with this framework, `f(v)>=c` defines a bounded
// object. Increasing the isovalue shrinks the object, and decreasing the isovalue grows the object.
// .
// Now for the arguments to this metaball() module or function....
// In order to adjust interaction strength, the influence parameter applies an exponent, so if `influence=a` then the
// decay becomes $1/d^(1/a)$. This means, for example, that if you set influence to 2 you get a $1/d^2$ falloff.
// Changing this exponent changes how the balls interact.
// .
// You can pass a custom function as a [function literal](https://en.wikibooks.org/wiki/OpenSCAD_User_Manual/User-Defined_Functions_and_Modules#Function_literals)
// that takes a single argument (a 3-vector) and returns a single numerical value.
// The returned value should define a function where in isovalue range [c,INF] defines a bounded object. See Example 11.
// Arguments:
// funcs = a 1-D list of transform and function pairs in the form `[trans0, func0, trans1, func1, ...]`, with one pair for each metaball. The transform should be at least a position such as `move([x,y,z])` to specify the location of the metaball center, but you can also include rotations, such as `move([x,y,z])*rot([ax,ay,az])`. You can multiply together any of BOSL2's affine operations like {{xrot()}}, {{scale()}}, and {{skew()}}. This is useful for orienting non-spherical metaballs. The priority order of the transforms is right to left, that is, `move([4,5,6])*rot([45,0,90])` does the rotation first, and then the move, similar to normal OpenSCAD syntax `translate([4,5,6]) rotate([45,0,90]) children()`.
// funcs = a 1-D list of transform and function pairs in the form `[trans0, func0, trans1, func1, ...]`, with one pair for each metaball. The transform should be at least a position such as `move([x,y,z])` to specify the location of the metaball center, but you can also include rotations, such as `move([x,y,z])*rot([ax,ay,az])`. You can multiply together any of BOSL2's affine operations like {{xrot()}}, {{scale()}}, and {{skew()}}. This is useful for orienting non-spherical metaballs. The priority order of the transforms is right to left, that is, `move([4,5,6])*rot([45,0,90])` does the rotation first, and then the move, similar to normal OpenSCAD syntax `translate([4,5,6]) rotate([45,0,90]) children()`. The transform `IDENT` may be used if you don't want to specify a transform, resulting in a metaball positioned at the origin. **Any function in the list may, itself, be another list of metaballs instead of a function name.** See Example 13 for a demonstration.
// voxel_size = The size (scalar) of the voxel cube that determines the resolution of the metaball surface. **Start with a larger size for experimenting, and refine it gradually.** A small voxel size can significantly slow down processing time, especially with a large `bounding_box`.
// bounding_box = A pair of 3D points `[[xmin,ymin,zmin], [xmax,ymax,zmax]]`, specifying the minimum and maximum box corner coordinates. The voxels needn't fit perfectly inside the bounding box.
// isovalue = A scalar value specifying the isosurface value (threshold value) of the metaballs. At the default value of 1.0, the internal metaball functions are designd so the coefficient corresponds to the radial size of the metaball, when rendered in isolation with no other metaballs. Default: 1.0
// isovalue = A scalar value specifying the isosurface value (threshold value) of the metaballs. At the default value of 1.0, the internal metaball functions are designd so the size arguments correspond to the size parameter (such as radius) of the metaball, when rendered in isolation with no other metaballs. Default: 1.0
// ---
// closed = When true, maintains a manifold surface where the bounding box clips it (there is a negligible speed penalty in doing this). When false, the bounding box clips the surface, exposing the back sides of facets. Setting this to false can be useful with OpenSCAD's "View > Thrown together" menu option to distinguish inside from outside. Default: true
// show_stats = If true, display statistics about the metaball isosurface in the console window. Besides the number of voxels found to contain the surface, and the number of triangles making up the surface, this is useful for getting information about a smaller bounding box possible, to improve speed for subsequent renders. Enabling this parameter has a speed penalty. Default: false
// convexity = Max number of times a line could intersect a wall of the shape. Affects preview only. Default: 6
// show_stats = If true, display statistics about the metaball isosurface in the console window. Besides the number of voxels found to contain the surface, and the number of triangles making up the surface, this is useful for getting information about a smaller bounding box possible, to improve speed for subsequent renders. Enabling this parameter has a small speed penalty. Default: false
// convexity = Maximum number of times a line could intersect a wall of the shape. Affects preview only. Default: 6
// cp = (Module only) Center point for determining intersection anchors or centering the shape. Determines the base of the anchor vector. Can be "centroid", "mean", "box" or a 3D point. Default: "centroid"
// anchor = (Module only) Translate so anchor point is at origin (0,0,0). See [anchor](attachments.scad#subsection-anchor). Default: `"origin"`
// spin = (Module only) Rotate this many degrees around the Z axis after anchor. See [spin](attachments.scad#subsection-spin). Default: `0`
@@ -1171,11 +1289,58 @@ function mb_cylinder(coeff=10, length=10, cutoff=INF, influence=1) = function (d
// "intersect" = Anchors to the surface of the shape.
// Named Anchors:
// "origin" = Anchor at the origin, oriented UP.
// Example(3D,VPD=110): This is the first of a series of five examples demonstrating the differnent types of metaball interactions. We start with two spheres 30 units apart. Each would have a radius of 10 in isolation, but because they are influencing their surroundings, each sphere mutually contributes to the size of the other. The sum of contributions between the spheres add up so that a surface plotted around the region exceeding the threshold defined by `isovalue=1` looks like a peanut shape surrounding the two spheres.
// funcs = [
// left(15), mb_sphere(10),
// right(15), mb_sphere(10)
// ];
// voxelsize = 1;
// boundingbox = [[-30,-19,-19], [30,19,19]];
// metaballs(funcs, voxelsize, boundingbox,
// isovalue=1);
// Example(3D,VPD=110): Adding a cutoff of 25 to the left sphere causes its influence to disappear completely 25 units away (which is the center of the right sphere). The left sphere is bigger because it still receives the full influence of the right sphere, but the right sphere is smaller because the left sphere has no contribution past 25 units. Setting cutoff too small can remove the interactions of one metaball from all other metaballs, leaving that metaball alone by itself.
// funcs = [
// left(15), mb_sphere(10, cutoff=25),
// right(15), mb_sphere(10)
// ];
// voxelsize = 1;
// boundingbox = [[-30,-19,-19], [30,19,19]];
// metaballs(funcs, voxelsize, boundingbox,
// isovalue=1);
// Example(3D,VPD=110): Here, the left sphere has less influence in addition to a cutoff. Setting `influence=0.5` results in a steeper falloff of contribution from the left sphere. Each sphere has a different size and shape due to unequal contributions based on distance.
// funcs = [
// left(15), mb_sphere(10, influence=0.5, cutoff=25),
// right(15), mb_sphere(10)
// ];
// voxelsize = 1;
// boundingbox = [[-30,-19,-19], [30,19,19]];
// metaballs(funcs, voxelsize, boundingbox,
// isovalue=1);
// Example(3D,VPD=110): In this example, we have two size-10 spheres as before and one tiny sphere of 1.5 units radius offset a bit on the y axis. With an isovalue of 1, this figure would appear similar to Example 1 above, but here the isovalue has been set to 2, causing the surface to shrink around a smaller volume values greater than 2. Remember, higher isovalue thresholds cause metaballs to shrink.
// funcs = [
// left(15), mb_sphere(10),
// right(15), mb_sphere(10),
// fwd(15), mb_sphere(1.5)
// ];
// voxelsize = 1;
// boundingbox = [[-30,-19,-19], [30,19,19]];
// metaballs(funcs, voxelsize, boundingbox,
// isovalue=2);
// Example(3D,VPD=110): The influence of the tiny sphere has been set quite high, to 10. Notice that the tiny sphere does not change in size, but its contribution to its surroundings has dramatically increased, causing the two other spheres to grow and meld into each other. The `influence` argument on a small metaball affects its surroundings more than itself.
// funcs = [
// move([-15,0,0]), mb_sphere(10),
// move([15,0,0]), mb_sphere(10),
// move([0,-15,0]), mb_sphere(1.5, influence=10)
// ];
// voxelsize = 1;
// boundingbox = [[-30,-19,-19], [30,19,19]];
// metaballs(funcs, voxelsize, boundingbox,
// isovalue=2);
// Example(3D,NoAxes): A group of five spherical metaballs with different sizes. The parameter `show_stats=true` (not shown here) was used to find a compact bounding box for this figure.
// funcs = [ // spheres of different sizes
// move([-20,-20,-20]), mb_sphere(5),
// move([0,-20,-20]), mb_sphere(4),
// move([0,0,0]), mb_sphere(3),
// IDENT, mb_sphere(3),
// move([0,0,20]), mb_sphere(5),
// move([20,20,10]), mb_sphere(7)
// ];
@@ -1196,8 +1361,8 @@ function mb_cylinder(coeff=10, length=10, cutoff=INF, influence=1) = function (d
// color("green") move_copies(centers) sphere(d=1, $fn=16);
// Example(3D,NoAxes): A cube, a rounded cube, and an octahedron interacting. Because the surface is generated through cubical voxels, voxel corners are always cut off, resulting in difficulty resolving some sharp edges.
// funcs = [
// move([-7,-3,27])*zrot(55), mb_roundcube(6, squareness=1),
// move([5,5,21]), mb_roundcube(5),
// move([-7,-3,27])*zrot(55), mb_cuboid(6, squareness=1),
// move([5,5,21]), mb_cuboid(5),
// move([10,0,10]), mb_octahedron(5)
// ];
// voxelsize = 0.5; // a bit slow at this resolution
@@ -1230,36 +1395,55 @@ function mb_cylinder(coeff=10, length=10, cutoff=INF, influence=1) = function (d
// ) mb_cutoff(dist,cutoff) * (r/dist)^(1/influence);
//
// funcs = [
// move([-9,0,0]), mb_sphere(5),
// move([9,0,0]), function (dv) noisy_sphere(dv, 5, 0.2),
// left(9), mb_sphere(5),
// right(9), function (dv) noisy_sphere(dv, 5, 0.2),
// ];
// voxelsize = 0.4;
// boundingbox = [[-16,-8,-8], [16,8,8]];
// metaballs(funcs, voxelsize, boundingbox, isovalue=1);
// Example(3D,Med,NoAxes,VPR=[55,0,0],VPD=200,VPT=[7,2,2]): A complex example using ellipsoids, spheres, and a torus to make a tetrahedral object with rounded feet and a ring on top. The bottoms of the feet are flattened by limiting the minimum z value of the bounding box. The center of the object is thick due to the contributions of four ellipsoids converging. Designing an object like this using metaballs requires trial and error with low-resolution renders.
// Example(3D,Med,NoAxes,VPR=[55,0,0],VPD=200,VPT=[7,2,2]): A complex example using ellipsoids, a cylinder, spheres, and a torus to make a tetrahedral object with rounded feet and a ring on top. The bottoms of the feet are flattened by limiting the minimum z value of the bounding box. The center of the object is thick due to the contributions of four ellipsoids converging. Designing an object like this using metaballs requires trial and error with low-resolution renders.
// include <BOSL2/polyhedra.scad>
// tetpts = zrot(15, p = 22 * regular_polyhedron_info("vertices", "tetrahedron"));
// tetxform = [ for(pt = tetpts) move(pt)*rot(from=RIGHT, to=pt)*scale([7,1.5,1.5]) ];
// tettransform = [ for(pt = tetpts) move(pt)*rot(from=RIGHT, to=pt)*scale([7,1.5,1.5]) ];
//
// funcs = [
// // vertical cylinder arm
// move([0,0,15]), mb_cylinder(2, 17, influence=0.8),
// up(15), mb_capsule(17, 2, influence=0.8),
// // ellipsoid arms
// tetxform[0], mb_sphere(1, cutoff=30),
// tetxform[1], mb_sphere(1, cutoff=30),
// tetxform[2], mb_sphere(1, cutoff=30),
// for(i=[0:2]) each [tettransform[i], mb_sphere(1, cutoff=30)],
// // ring on top
// move([0,0,35])*xrot(90), mb_torus(r_maj=8, r_min=2.5, cutoff=35),
// up(35)*xrot(90), mb_torus(r_maj=8, r_min=2.5, cutoff=35),
// // feet
// move(2.2*tetpts[0]), mb_sphere(5, cutoff=30),
// move(2.2*tetpts[1]), mb_sphere(5, cutoff=30),
// move(2.2*tetpts[2]), mb_sphere(5, cutoff=30)
// for(i=[0:2]) each [move(2.2*tetpts[i]), mb_sphere(5, cutoff=30)],
// ];
// voxelsize = 1;
// boundingbox = [[-22,-32,-13], [36,32,46]];
// // useful to save as VNF for copies and manipulations
// vnf = metaballs(funcs, voxelsize, boundingbox, isovalue=1);
// vnf_polyhedron(vnf);
// Example(3D,Med,NoAxes): Example of grouping metaballs together and nesting them in lists of other metaballs. Here, just one finger is defined, and a thumb is defined from one less joint in the finger. Individual fingers are grouped together with different positions and scaling, along with the thumb. Finally, this group of all fingers is used to combine with a rounded cuboid, with a slight dent subtracted to hollow out the palm, to make the hand.
// joints = [[0,0,1], [0,0,85], [0,-5,125], [0,-16,157], [0,-30,178]];
// finger = [
// for(i=[0:3]) each [IDENT, mb_connector(joints[i], joints[i+1], 6+i/4, influence=.5)]
// ];
// thumb = [
// for(i=[0:2]) each [IDENT, mb_connector(joints[i], joints[i+1], 6+i, influence=.4)]
// ];
// allfingers = [
// left(25)*zrot(5)*yrot(-50)*scale([1,1,0.6])*zrot(30), thumb,
// left(15)*yrot(-9)*scale([1,1,0.9]), finger,
// IDENT, finger,
// right(15)*yrot(8)*scale([1,1,0.92]), finger,
// right(30)*yrot(17)*scale([0.9,0.9,0.75]), finger
// ];
// hand = [
// IDENT, allfingers,
// scale([1,0.25,1.4]), mb_cuboid(75, squareness=0.3, cutoff=80),
// move([-10,-100,55])*yrot(10)*scale([2,2,1]),
// mb_sphere(r=15, cutoff=50, influence=2, negative=true)
// ];
// bbox = [[-100,-40,-10], [76,18,186]];
// metaballs(hand, 2, bbox, isovalue=1, show_stats=true);
module metaballs(funcs, voxel_size, bounding_box, isovalue=1, closed=true, convexity=6, cp="centroid", anchor="origin", spin=0, orient=UP, atype="hull", show_stats=false) {
vnf = metaballs(funcs, voxel_size, bounding_box, isovalue, closed, show_stats);
@@ -1271,36 +1455,67 @@ function metaballs(funcs, voxel_size, bounding_box, isovalue=1, closed=true, sho
assert(all_defined([funcs, isovalue, bounding_box, voxel_size]), "\nThe parameters funcs, isovalue, bounding_box, and voxel_size must all be defined.")
assert(len(funcs)%2==0, "\nThe funcs parameter must be an even-length list of alternating transforms and functions")
let(
nballs = len(funcs)/2,
funclist = _mb_unwind_list(funcs),
nballs = len(funclist)/2,
dummycheck = [
for(i=[0:len(funcs)/2-1]) let(j=2*i)
assert(is_matrix(funcs[j],4,4), str("\nfuncs entry at position ", j, " must be a 4×4 matrix."))
assert(is_function(funcs[j+1]) || is_list(funcs[j+1]), str("\nfuncs entry at position ", j+1, " must be a function literal or a metaball list.")) 0
],
// set up transformation matrices in advance
transmatrix = [
for(i=[0:nballs-1])
let(j=2*i)
assert(is_matrix(funcs[j],4,4), str("\nfuncs entry at position ", j, " must be a 4×4 matrix."))
assert(is_function(funcs[j+1]), str("\nfuncs entry at position ", j+1, " must be a function literal."))
transpose(select(matrix_inverse(funcs[j]), 0,2))
transpose(select(matrix_inverse(funclist[j]), 0,2))
],
// set up field array
bot = bounding_box[0],
top = bounding_box[1],
// new bounding box centered around original, forced to integer multiples of voxel size
halfvox = 0.5*voxel_size,
bbcenter = mean(bounding_box),
bbnums = v_ceil((bounding_box[1]-bounding_box[0]) / voxel_size),
newbbox = [bbcenter - halfvox*bbnums, bbcenter + halfvox*bbnums],
// set up field array
bot = newbbox[0],
top = newbbox[1],
// accumulate metaball contributions using matrices rather than sums
xset = [bot.x:voxel_size:top.x+halfvox],
yset = list([bot.y:voxel_size:top.y+halfvox]),
zset = list([bot.z:voxel_size:top.z+halfvox]),
allpts = [for(x=xset, y=yset, z=zset) [x,y,z,1]],
trans_pts = [for(i=[0:nballs-1]) allpts*transmatrix[i]],
allvals = [for(i=[0:nballs-1]) [for(pt=trans_pts[i]) funcs[2*i+1](pt)]],
total = _sum(allvals,allvals[0]*0),
allvals = [for(i=[0:nballs-1]) [for(pt=trans_pts[i]) funclist[2*i+1](pt)]],
//total = _sum(allvals,allvals[0]*EPSILON),
total = _sum(slice(allvals,1,-1), allvals[0]),
fieldarray = list_to_matrix(list_to_matrix(total,len(zset)),len(yset))
) isosurface(fieldarray, isovalue, voxel_size, closed=closed, show_stats=show_stats, _origin=bounding_box[0]);
) isosurface(fieldarray, isovalue, voxel_size, closed=closed, show_stats=show_stats, _mb_origin=newbbox[0]);
function _mb_unwind_list(list, parent_trans=[IDENT]) =
let(
dum1 = assert(is_list(list), "\nDid not find valid list of metaballs."),
n=len(list),
dum2 = assert(n%2==0, "\nList of metaballs must have an even number of elements with alternating transforms and functions/lists.")
) [
for(i=[0:2:n-1])
let(
dum = assert(is_matrix(list[i],4,4), str("\nInvalid 4×4 transformation matrix found at position ",i,".")),
trans = parent_trans[0] * list[i],
j=i+1
) if(is_function(list[j]))
each [trans, list[j]]
else if (is_list(list[j]))
each _mb_unwind_list(list[j], [trans])
else
assert(false, str("\nExpected function literal or list at position ",j,"."))
];
/// ---------- isosurface stuff starts here ----------
// Function&Module: isosurface()
// Synopsis: Creates a 3D isosurface.
// Synopsis: Creates a 3D isosurface (a 3D contour) from a function or array of values.
// SynTags: Geom,VNF
// Topics: Isosurfaces, VNF Generators
// Usage: As a module
@@ -1308,52 +1523,65 @@ function metaballs(funcs, voxel_size, bounding_box, isovalue=1, closed=true, sho
// Usage: As a function
// vnf = isosurface(f, isovalue, voxel_size, bounding_box, [reverse=], [closed=], [show_stats=]);
// Description:
// When called as a function, returns a [VNF structure](vnf.scad) (list of triangles and faces) representing a 3D isosurface within the specified bounding box at a single isovalue or range of isovalues.
// When called as a module, displays the isosurface within the specified bounding box at a single isovalue or range of isovalues. This module just passes the parameters to the function, and then calls {{vnf_polyhedron()}} to display the isosurface.
// Computes a [VNF structure](vnf.scad) of a 3D isosurface within a bounded box at a single
// isovalue or range of isovalues.
// The isosurface of a function $f(x,y,z)$ is the set of points where $f(x,y,z)=c$ for some
// constant isovalue, $c$.
// To provide a function you supply a [function literal](https://en.wikibooks.org/wiki/OpenSCAD_User_Manual/User-Defined_Functions_and_Modules#Function_literals)
// taking three parameters as input to define the grid coordinate location (e.g. `x,y,z`) and
// returning a single numerical value.
// You can also define an isosurface using a 3D array of values instead of a function, in which
// case the isosurface is the set of points where the array is equal to the isovalue. The array
// indices are in the order `[x][y][z]`.
// .
// A [marching cubes](https://en.wikipedia.org/wiki/Marching_cubes) algorithm is used
// to identify an envelope containing the isosurface within the bounding box. The surface
// intersecttion with a voxel cube is then triangulated to form a surface fragment, which is
// combined with all other surface fragments. Ambiguities in triangulating the surfaces
// in certain voxel cube configurations are resolved so that all triangular facets are
// properly oriented with no holes in the surface. If a side of the bounding box clips
// the isosurface, this clipped area is filled in so that the surface remains manifold.
// The VNF that is computed has the isosurface as its bounding surface, with all the points where
// $f(x,y,z)>c$ on the interior side of the surface.
// When the isovalue is a range, `[c1, c2]`, then the resulting VNF has two bounding surfaces
// corresponding to `c1` and `c2`, and the interior of the object are the points with intermediate
// isovalues; this generally produces a shell object that has an inside and outside surface. The
// range can start at `-INF` or end at `INF`. A single isovalue `c` is equivalent to `[c,INF]`.
// .
// Be mindful of how you set `voxel_size` and `bounding_box`. For example a voxel size
// of 1 unit with a bounding box volume of 200×200×200 may be noticeably slow,
// requiring calculation and storage of 8,000,000 field values, and more processing
// and memory to generate the triangulated mesh. On the other hand, a voxel size of 5
// in a 100×100×100 bounding box requires only 8,000 field values and the mesh
// generates fairly quickly, just a handful of seconds. A good rule is to keep the
// number of field values below 10,000 for preview, and adjust the voxel size
// smaller for final rendering. If the isosurface fits completely within the bounding
// box, you can call {{pointlist_bounds()}} on `vnf[0]` returned from the
// `isosurface()` function to get an idea of a more optimal smaller bounding box to use,
// possibly allowing increasing resolution by decresing the voxel size. You can also set
// the parameter `show_stats=true` to get the bounds of the voxels containing the surface.
// The isosurface is evaluated over a bounding box which is divided into voxels of the specified
// `voxel_size`. Smaller voxels produce a finer, smoother result at the expense of execution time.
// If the voxel size doesn't exactly divide your specified bounding box, then the bounding box is
// enlarged enough to contain all whole voxels, and centered on your requested box.
// If the bounding box clips the isosurface, then if `closed=true` (default) the clipped area is
// filled in to produce a closed surface. Setting `closed=false` causes the VNF to end at the
// bounding box, resulting in a non-manifold shape with holes, exposing the inside of the object.
// .
// The `voxel_size` and `bounding_box` parameters affect the run time, which can be long.
// A voxel size of 1 with a bounding box volume of 200×200×200 may be slow because it requires the
// calculation and storage of 8,000,000 function values, and more processing and memory to generate
// the triangulated mesh. On the other hand, a voxel size of 5 over a 100×100×100 bounding box
// requires only 8,000 function values and a modest computation time. A good rule is to keep the
// number of voxels below 10,000 for preview, and adjust the voxel size smaller for final
// rendering. A bounding box that is larger than your isosurface wastes time computing function
// values that are not needed. If the isosurface fits completely within the bounding box, you can
// call {{pointlist_bounds()}} on `vnf[0]` returned from the `isosurface()` function to get an
// idea of a the optimal bounding box to use. You may be able to decrease run time, or keep the
// same run time but increase the resolution. You can also set the parameter `show_stats=true` to
// get the bounds of the voxels containing the surface.
// .
// The point list in the VNF structure contains many duplicated points. This is not a
// problem for rendering the shape, but if you want to eliminate these, you can pass
// the structure to {{vnf_merge_points()}}. Additionally, flat surfaces (often
// resulting from clipping by the bounding box) are triangulated at the voxel size
// resolution, and these can be unified into a single face by passing the vnf
// structure to {{vnf_unify_faces()}}. These steps can be expensive for execution time
// structure to {{vnf_unify_faces()}}. These steps can be computationall expensive
// and are not normally necessary.
// Arguments:
// f = The isosurface function. Can be a [function literal](https://en.wikibooks.org/wiki/OpenSCAD_User_Manual/User-Defined_Functions_and_Modules#Function_literals) taking as input the `x,y,z` grid coordinates and returning a single value, or a 3D array (all points in the grid precomputed).
// **As a function literal:** Say you have you created your own function, `my_func(x,y,z,a,b,c)` (call it whatever you want), which depends on x, y, z, and additional parameters a, b, c, and returns a single value. In the parameter list to `isosurface()`, you would set the `f` parameter to `function (x,y,z) my_func(x,y,z,a,b,c)`.
// **As an array:** The array you pass in should be organized so that the indices are in order of `[x][y][z]` when the array is referenced; that is, `f[x_index][y_index][z_index]` has `z_index` changing most rapidly as the array is traversed. If you organize the array differently, you may have to perform a `rotate()` or `mirror()` operation on the final result to orient it properly.
// isovalue = As a scalar, specifies the output value of `field_function` corresponding to the isosurface. As a vector `[min_isovalue, max_isovalue]`, specifies the range of isovalues around which to generate a surface. For closed surfaces, a single value results in a closed volume, and a range results in a shell (with an inside and outside surface) enclosing a volume. A range must be specified for infinite-extent surfaces (such as gyroids) to create a manifold shape within the bounding box.
// voxel_size = The size (scalar) of the voxel cube that determines the resolution of the surface.
// bounding_box = Applicable only when `f` is a function literal. This is a pair of 3D points `[[xmin,ymin,zmin], [xmax,ymax,zmax]]`, specifying the minimum and maximum corner coordinates of the bounding box. You don't have ensure that the voxels fit perfectly inside the bounding box. While the voxel at the minimum bounding box corner is aligned on that corner, the last voxel at the maximum box corner may extend a bit beyond it. Default: undef
// f = The isosurface function or array.
// isovalue = a scalar giving the isovalue parameter or a 2-vector giving an isovalue range.
// voxel_size = scalar size of the voxel cube that is used to sample the surface.
// bounding_box = When `f` is a function, a pair of 3D points `[[xmin,ymin,zmin], [xmax,ymax,zmax]]`, specifying the minimum and maximum corner coordinates of the bounding box. The actual bounding box enlarged if necessary to make the voxels fit perfectly, and centered around your requested box.
// ---
// reverse = When true, reverses the orientation of the facets in the mesh. Default: false
// closed = When true, maintains a manifold surface where the bounding box clips it (there is a negligible speed penalty in doing this). When false, the bounding box clips the surface, exposing the back sides of facets. Setting this to false can be useful with OpenSCAD's "View > Thrown Together" menu option to distinguish inside from outside. Default: true
// show_stats = If true, display statistics about the isosurface in the console window. Besides the number of voxels found to contain the surface, and the number of triangles making up the surface, this is useful for getting information about a smaller bounding box possible for the isosurface, to improve speed for subsequent renders. Enabling this parameter has a speed penalty. Default: false
// convexity = Max number of times a line could intersect a wall of the shape. Affects preview only. Default: 6
// closed = When true, close the surface if it intersects the bounding box by adding a closing face. When false, do not add a closing face and instead produce a non-manfold VNF that has holes. Default: true
// reverse = When true, reverses the orientation of the VNF faces. Default: false
// show_stats = If true, display statistics in the console window about the isosurface: number of voxels that contain the surface, number of triangles, bounding box of the voxels, and voxel-rounded bounding box of the surface, which may help you reduce your bounding box to improve speed. Enabling this parameter has a slight speed penalty. Default: false
// convexity = Maximum number of times a line could intersect a wall of the shape. Affects preview only. Default: 6
// cp = (Module only) Center point for determining intersection anchors or centering the shape. Determines the base of the anchor vector. Can be "centroid", "mean", "box" or a 3D point. Default: "centroid"
// anchor = (Module only) Translate so anchor point is at origin (0,0,0). See [anchor](attachments.scad#subsection-anchor). Default: `"origin"`
// spin = (Module only) Rotate this many degrees around the Z axis after anchor. See [spin](attachments.scad#subsection-spin). Default: `0`
// anchor = (Module only) Translate so anchor point is at origin (0,0,0). See [anchor](attachments.scad#subsection-anchor). Default: `"origin"`
// spin = (Module only) Rotate this many degrees around the Z axis after anchor. See [spin](attachments.scad#subsection-spin). Default: `0`
// orient = (Module only) Vector to rotate top toward, after spin. See [orient](attachments.scad#subsection-orient). Default: `UP`
// atype = (Module only) Select "hull" or "intersect" anchor type. Default: "hull"
// Anchor Types:
@@ -1447,35 +1675,37 @@ function metaballs(funcs, voxel_size, bounding_box, isovalue=1, closed=true, sho
// isosurface(field, isovalue=0.5,
// voxel_size=10);
module isosurface(f, isovalue, voxel_size, bounding_box, reverse=false, closed=true, convexity=6, cp="centroid", anchor="origin", spin=0, orient=UP, atype="hull", show_stats=false, _origin=undef) {
vnf = isosurface(f, isovalue, voxel_size, bounding_box, reverse, closed, show_stats, _origin);
module isosurface(f, isovalue, voxel_size, bounding_box, reverse=false, closed=true, convexity=6, cp="centroid", anchor="origin", spin=0, orient=UP, atype="hull", show_stats=false, _mb_origin=undef) {
vnf = isosurface(f, isovalue, voxel_size, bounding_box, reverse, closed, show_stats, _mb_origin);
vnf_polyhedron(vnf, convexity=convexity, cp=cp, anchor=anchor, spin=spin, orient=orient, atype=atype)
children();
}
function isosurface(f, isovalue, voxel_size, bounding_box, reverse=false, closed=true, show_stats=false, _origin=undef) =
function isosurface(f, isovalue, voxel_size, bounding_box, reverse=false, closed=true, show_stats=false, _mb_origin=undef) =
assert(all_defined([f, isovalue, voxel_size]), "\nThe parameters f, isovalue, and bounding_box must all be defined.")
assert((is_function(f) && is_def(bounding_box)) || (is_list(f) && is_undef(bounding_box)),
"\nbounding_box must be passed if f is a function, and cannot be passed if f is an array.")
let(
isovalmin = is_list(isovalue) ? isovalue[0] : isovalue,
isovalmax = is_list(isovalue) ? isovalue[1] : INF,
dum1 = assert(isovalmin < isovalmax, str("\nBad isovalue range (", isovalmin, ", >= ", isovalmax, "), should be expressed as [min_value, max_value].")),
hv = 0.5*voxel_size,
bbox = is_function(f)
? let( // new bounding box quantized for voxel_size
hv = 0.5*voxel_size,
bbn = (bounding_box[1]-bounding_box[0]+[hv,hv,hv]) / voxel_size,
bbsize = [round(bbn[0]), round(bbn[1]), round(bbn[2])] * voxel_size
) [bounding_box[0], bounding_box[0]+bbsize]
: let(
nx = len(f)-1,
ny = len(f[0])-1,
nz = len(f[0][0])-1
) is_def(_origin) ? [_origin, _origin+voxel_size*[nx,ny,nz]]
: [-0.5*voxel_size*[nx,ny,nz], 0.5*voxel_size*[nx, ny, nz]],
cubes = _isosurface_cubes(voxel_size, bbox, fieldarray=is_function(f)?undef:f, fieldfunc=is_function(f)?f:undef, isovalmin=isovalmin, isovalmax=isovalmax, closed=closed),
? let( // new bounding box quantized for voxel_size, centered around original box
bbcenter = mean(bounding_box),
bbn = v_ceil((bounding_box[1]-bounding_box[0]) / voxel_size)
) [bbcenter - hv*bbn, bbcenter + hv*bbn]
: let( // new bounding box, either centered on origin or using metaball origin
dims = list_shape(f) - [1,1,1]
) is_def(_mb_origin)
? [_mb_origin, _mb_origin+voxel_size*dims] // metaball bounding box
: [-hv*dims, hv*dims], // centered bounding box
cubes = _isosurface_cubes(voxel_size, bbox,
fieldarray=is_function(f)?undef:f, fieldfunc=is_function(f)?f:undef,
isovalmin=isovalmin, isovalmax=isovalmax, closed=closed),
tritablemin = reverse ? _MCTriangleTable_reverse : _MCTriangleTable,
tritablemax = reverse ? _MCTriangleTable : _MCTriangleTable_reverse,
trianglepoints = _isosurface_triangles(cubes, voxel_size, isovalmin, isovalmax, tritablemin, tritablemax),
faces = [ for(i=[0:3:len(trianglepoints)-1]) [i,i+1,i+2] ],
dummy = show_stats ? _showstats(voxel_size, bbox, isovalmin, cubes, faces) : 0
dum2 = show_stats ? _showstats(voxel_size, bbox, isovalmin, cubes, faces) : 0
) [trianglepoints, faces];