From 0f82e6e9584f3f435632d6c05a18cb0fbdba9c87 Mon Sep 17 00:00:00 2001 From: Richard Milewski Date: Wed, 5 Feb 2025 11:59:13 -0800 Subject: [PATCH] Update isosurface.scad --- isosurface.scad | 1484 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 1484 insertions(+) diff --git a/isosurface.scad b/isosurface.scad index 78192971..84f6ae28 100644 --- a/isosurface.scad +++ b/isosurface.scad @@ -1,3 +1,4 @@ +<<<<<<< Updated upstream ///////////////////////////////////////////////////////////////////// // LibFile: isosurface.scad // An isosurface is a three-dimensional surface representing points of a constant @@ -1411,3 +1412,1486 @@ for(i=[0:nballs-1]) let( dumm5 = assert(interact[i] != MB_CUSTOM || (interact[i]==MB_CUSTOM && is_def(fieldfuncs[i])), "\nThe MB_CUSTOM ball_type requires a field_function to be defined.") 0 ) 0 ]; +======= +///////////////////////////////////////////////////////////////////// +// 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. +// . +// An isosurface is a three-dimensional surface representing 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. +// . +// 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. +// . +// 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 +// $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. +// . +// 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. +// +// Includes: +// include +// include +// FileGroup: Advanced Modeling +// FileSummary: Isosurfaces and metaballs. +////////////////////////////////////////////////////////////////////// + + +/* +Lookup Tables for Transvoxel's Modified Marching Cubes + +Adapted for OpenSCAD from https://gist.github.com/dwilliamson/72c60fcd287a94867b4334b42a7888ad + +Unlike the original paper (Marching Cubes: A High Resolution 3D Surface Construction Algorithm), these tables guarantee a closed mesh in which connected components are continuous and free of holes. + +Rotations are prioritized over inversions so that 3 of the 6 cases containing ambiguous faces are never added. 3 extra cases are added as a post-process, overriding inversions through custom-built rotations to eliminate the remaining ambiguities. + +The cube index determines the sequence of edges to split. The index ranges from 0 to 255, representing all possible combinations of the 8 corners of the cube being greater or less than the isosurface threshold. + +For example, a cube with corners 2, 3, and 7 greater than the threshold isovalue would have the index 10000110, an 8-bit binary number with bits 2, 3, and 7 set to 1, corresponding to decimal index 134. After determining the cube's index value this way, the triangulation order is looked up in a table. + +Axes are + z + (top) + | y (back) + | / + |/ + +----- x (right) + +Vertex and edge layout (heavier = and # indicate closer to viewer): + + 3 +----------+ 7 +----10----+ + /: /| /: /| + / : / | 1 2 5 6 + 1 +==========+5 | +=====9====+ | + # 2+ - - - # -+ 6 # +- - 11-# -+ + # / # / 0 3 4 7 + #/ #/ #/ #/ + 0 +==========+ 4 +=====8=====+ + +z changes fastest, then y, then x. + +----------------------------------------------------------- +Addition by Alex: +Vertex and face layout for triangulating one voxel face that corrsesponds to a side of the box bounding all voxels. + + 4(back) + 3 +----------+ 7 + /: 5(top) /| + / : / | + 1 +==========+5 | <-- 3 (side) +0(side) --> # 2+ - - - # -+ 6 + # / # / + #/ 2(bot) #/ + 0 +----------+ 4 + 1(front) +*/ + +/// four indices for each face of the cube, counterclockwise looking from inside out +_MCFaceVertexIndices = [ + [], + [0,2,3,1], // left, x=0 plane + [0,1,5,4], // front, y=0 plane + [0,4,6,2], // bottom, z=0 plane + [4,5,7,6], // right, x=voxsize plane + [2,6,7,3], // back, y=voxsize plane + [1,3,7,5], // top, z=voxsize plane +]; + +/// return an array of face indices in _MCFaceVertexIndices if the voxel at coordinate v0 corresponds to the bounding box. +function _bbox_faces(v0, voxsize, bbox) = let( + a = v0-bbox[0], + bb1 = bbox[1] - [voxsize,voxsize,voxsize], + b = v0-bb1 +) [ + if(a[0]==0) 1, + if(a[1]==0) 2, + if(a[2]==0) 3, + if(b[0]>=0) 4, + if(b[1]>=0) 5, + if(b[2]>=0) 6 +]; +/// End of bounding-box face-clipping stuff. Back to the marching cubes triangulation.... + + +/// Pair of vertex indices for each edge on the voxel +_MCEdgeVertexIndices = [ + [0, 1], + [1, 3], + [3, 2], + [2, 0], + [4, 5], + [5, 7], + [7, 6], + [6, 4], + [0, 4], + [1, 5], + [3, 7], + [2, 6] +]; + +/// For each of the 256 configurations of a marching cube, define a list of triangles, specified as triples of edge indices. +_MCTriangleTable = [ + [], + [3,8,0], + [1,0,9], + [9,1,8,8,1,3], + [3,2,11], + [2,11,0,0,11,8], + [1,0,9,3,2,11], + [11,1,2,11,9,1,11,8,9], + [10,2,1], + [2,1,10,0,3,8], + [0,9,2,2,9,10], + [8,2,3,8,10,2,8,9,10], + [1,10,3,3,10,11], + [10,0,1,10,8,0,10,11,8], + [9,3,0,9,11,3,9,10,11], + [9,10,8,8,10,11], + [7,4,8], + [0,3,4,4,3,7], + [0,9,1,4,8,7], + [1,4,9,1,7,4,1,3,7], + [11,3,2,8,7,4], + [4,11,7,4,2,11,4,0,2], + [3,2,11,0,9,1,4,8,7], + [9,1,4,4,1,7,7,1,2,7,2,11], + [7,4,8,1,10,2], + [7,4,3,3,4,0,10,2,1], + [10,2,9,9,2,0,7,4,8], + [7,4,9,7,9,2,9,10,2,3,7,2], + [1,10,3,3,10,11,4,8,7], + [4,0,7,0,1,10,7,0,10,7,10,11], + [7,4,8,9,3,0,9,11,3,9,10,11], + [7,4,11,4,9,11,9,10,11], + [5,9,4], + [8,0,3,9,4,5], + [1,0,5,5,0,4], + [5,8,4,5,3,8,5,1,3], + [3,2,11,5,9,4], + [2,11,0,0,11,8,5,9,4], + [4,5,0,0,5,1,11,3,2], + [11,8,2,8,4,5,2,8,5,2,5,1], + [5,9,4,1,10,2], + [0,3,8,1,10,2,5,9,4], + [2,5,10,2,4,5,2,0,4], + [4,5,8,8,5,3,3,5,10,3,10,2], + [11,3,10,10,3,1,4,5,9], + [4,5,9,10,0,1,10,8,0,10,11,8], + [4,5,10,4,10,3,10,11,3,0,4,3], + [4,5,8,5,10,8,10,11,8], + [5,9,7,7,9,8], + [3,9,0,3,5,9,3,7,5], + [7,0,8,7,1,0,7,5,1], + [3,7,1,1,7,5], + [5,9,7,7,9,8,2,11,3], + [5,9,0,5,0,11,0,2,11,7,5,11], + [2,11,3,7,0,8,7,1,0,7,5,1], + [2,11,1,11,7,1,7,5,1], + [8,7,9,9,7,5,2,1,10], + [10,2,1,3,9,0,3,5,9,3,7,5], + [2,0,10,0,8,7,10,0,7,10,7,5], + [10,2,5,2,3,5,3,7,5], + [5,9,8,5,8,7,1,10,3,10,11,3], + [1,10,0,0,10,11,0,11,7,0,7,5,0,5,9], + [8,7,0,0,7,5,0,5,10,0,10,11,0,11,3], + [5,11,7,10,11,5], + [11,6,7], + [3,8,0,7,11,6], + [1,0,9,7,11,6], + [9,1,8,8,1,3,6,7,11], + [6,7,2,2,7,3], + [0,7,8,0,6,7,0,2,6], + [6,7,2,2,7,3,9,1,0], + [9,1,2,9,2,7,2,6,7,8,9,7], + [10,2,1,11,6,7], + [2,1,10,3,8,0,7,11,6], + [0,9,2,2,9,10,7,11,6], + [6,7,11,8,2,3,8,10,2,8,9,10], + [7,10,6,7,1,10,7,3,1], + [1,10,0,0,10,8,8,10,6,8,6,7], + [9,10,0,10,6,7,0,10,7,0,7,3], + [6,7,10,7,8,10,8,9,10], + [4,8,6,6,8,11], + [6,3,11,6,0,3,6,4,0], + [11,6,8,8,6,4,1,0,9], + [6,4,11,4,9,1,11,4,1,11,1,3], + [2,8,3,2,4,8,2,6,4], + [0,2,4,4,2,6], + [9,1,0,2,8,3,2,4,8,2,6,4], + [9,1,4,1,2,4,2,6,4], + [4,8,6,6,8,11,1,10,2], + [1,10,2,6,3,11,6,0,3,6,4,0], + [0,9,10,0,10,2,4,8,6,8,11,6], + [11,6,3,3,6,4,3,4,9,3,9,10,3,10,2], + [1,10,6,1,6,8,6,4,8,3,1,8], + [1,10,0,10,6,0,6,4,0], + [0,9,3,3,9,10,3,10,6,3,6,4,3,4,8], + [4,10,6,9,10,4], + [4,5,9,6,7,11], + [7,11,6,8,0,3,9,4,5], + [1,0,5,5,0,4,11,6,7], + [11,6,7,5,8,4,5,3,8,5,1,3], + [3,2,7,7,2,6,9,4,5], + [5,9,4,0,7,8,0,6,7,0,2,6], + [1,0,4,1,4,5,3,2,7,2,6,7], + [4,5,8,8,5,1,8,1,2,8,2,6,8,6,7], + [6,7,11,5,9,4,1,10,2], + [5,9,4,7,11,6,0,3,8,2,1,10], + [7,11,6,2,5,10,2,4,5,2,0,4], + [6,7,11,3,8,4,3,4,5,3,5,2,2,5,10], + [9,4,5,7,10,6,7,1,10,7,3,1], + [5,9,4,8,0,1,8,1,10,8,10,7,7,10,6], + [6,7,10,10,7,3,10,3,0,10,0,4,10,4,5], + [4,5,8,8,5,10,8,10,6,8,6,7], + [9,6,5,9,11,6,9,8,11], + [0,3,9,9,3,5,5,3,11,5,11,6], + [1,0,8,1,8,6,8,11,6,5,1,6], + [11,6,3,6,5,3,5,1,3], + [2,6,3,6,5,9,3,6,9,3,9,8], + [5,9,6,9,0,6,0,2,6], + [3,2,8,8,2,6,8,6,5,8,5,1,8,1,0], + [1,6,5,2,6,1], + [2,1,10,9,6,5,9,11,6,9,8,11], + [2,1,10,5,9,0,5,0,3,5,3,6,6,3,11], + [10,2,5,5,2,0,5,0,8,5,8,11,5,11,6], + [10,2,5,5,2,3,5,3,11,5,11,6], + [5,9,6,6,9,8,6,8,3,6,3,1,6,1,10], + [5,9,6,6,9,0,6,0,1,6,1,10], + [8,3,0,5,10,6], + [6,5,10], + [6,10,5], + [3,8,0,5,6,10], + [9,1,0,10,5,6], + [3,8,1,1,8,9,6,10,5], + [6,10,5,2,11,3], + [8,0,11,11,0,2,5,6,10], + [10,5,6,1,0,9,3,2,11], + [5,6,10,11,1,2,11,9,1,11,8,9], + [2,1,6,6,1,5], + [5,6,1,1,6,2,8,0,3], + [6,9,5,6,0,9,6,2,0], + [8,9,3,9,5,6,3,9,6,3,6,2], + [3,6,11,3,5,6,3,1,5], + [5,6,11,5,11,0,11,8,0,1,5,0], + [0,9,3,3,9,11,11,9,5,11,5,6], + [5,6,9,6,11,9,11,8,9], + [7,4,8,5,6,10], + [0,3,4,4,3,7,10,5,6], + [4,8,7,9,1,0,10,5,6], + [6,10,5,1,4,9,1,7,4,1,3,7], + [11,3,2,7,4,8,5,6,10], + [10,5,6,4,11,7,4,2,11,4,0,2], + [7,4,8,3,2,11,9,1,0,10,5,6], + [10,5,6,7,4,9,7,9,1,7,1,11,11,1,2], + [2,1,6,6,1,5,8,7,4], + [7,4,0,7,0,3,5,6,1,6,2,1], + [8,7,4,6,9,5,6,0,9,6,2,0], + [5,6,9,9,6,2,9,2,3,9,3,7,9,7,4], + [4,8,7,3,6,11,3,5,6,3,1,5], + [7,4,11,11,4,0,11,0,1,11,1,5,11,5,6], + [4,8,7,11,3,0,11,0,9,11,9,6,6,9,5], + [5,6,9,9,6,11,9,11,7,9,7,4], + [9,4,10,10,4,6], + [6,10,4,4,10,9,3,8,0], + [0,10,1,0,6,10,0,4,6], + [3,8,4,3,4,10,4,6,10,1,3,10], + [9,4,10,10,4,6,3,2,11], + [8,0,2,8,2,11,9,4,10,4,6,10], + [11,3,2,0,10,1,0,6,10,0,4,6], + [2,11,1,1,11,8,1,8,4,1,4,6,1,6,10], + [4,1,9,4,2,1,4,6,2], + [3,8,0,4,1,9,4,2,1,4,6,2], + [4,6,0,0,6,2], + [3,8,2,8,4,2,4,6,2], + [3,1,11,1,9,4,11,1,4,11,4,6], + [9,4,1,1,4,6,1,6,11,1,11,8,1,8,0], + [11,3,6,3,0,6,0,4,6], + [8,6,11,4,6,8], + [10,7,6,10,8,7,10,9,8], + [10,9,6,9,0,3,6,9,3,6,3,7], + [8,7,0,0,7,1,1,7,6,1,6,10], + [6,10,7,10,1,7,1,3,7], + [3,2,11,10,7,6,10,8,7,10,9,8], + [6,10,7,7,10,9,7,9,0,7,0,2,7,2,11], + [11,3,2,1,0,8,1,8,7,1,7,10,10,7,6], + [6,10,7,7,10,1,7,1,2,7,2,11], + [8,7,6,8,6,1,6,2,1,9,8,1], + [0,3,9,9,3,7,9,7,6,9,6,2,9,2,1], + [8,7,0,7,6,0,6,2,0], + [7,2,3,6,2,7], + [11,3,6,6,3,1,6,1,9,6,9,8,6,8,7], + [11,7,6,1,9,0], + [11,3,6,6,3,0,6,0,8,6,8,7], + [11,7,6], + [10,5,11,11,5,7], + [10,5,11,11,5,7,0,3,8], + [7,11,5,5,11,10,0,9,1], + [3,8,9,3,9,1,7,11,5,11,10,5], + [5,2,10,5,3,2,5,7,3], + [0,2,8,2,10,5,8,2,5,8,5,7], + [0,9,1,5,2,10,5,3,2,5,7,3], + [10,5,2,2,5,7,2,7,8,2,8,9,2,9,1], + [1,11,2,1,7,11,1,5,7], + [8,0,3,1,11,2,1,7,11,1,5,7], + [0,9,5,0,5,11,5,7,11,2,0,11], + [3,8,2,2,8,9,2,9,5,2,5,7,2,7,11], + [5,7,1,1,7,3], + [8,0,7,0,1,7,1,5,7], + [0,9,3,9,5,3,5,7,3], + [9,7,8,5,7,9], + [8,5,4,8,10,5,8,11,10], + [10,5,4,10,4,3,4,0,3,11,10,3], + [1,0,9,8,5,4,8,10,5,8,11,10], + [9,1,4,4,1,3,4,3,11,4,11,10,4,10,5], + [10,5,2,2,5,3,3,5,4,3,4,8], + [10,5,2,5,4,2,4,0,2], + [9,1,0,3,2,10,3,10,5,3,5,8,8,5,4], + [10,5,2,2,5,4,2,4,9,2,9,1], + [1,5,2,5,4,8,2,5,8,2,8,11], + [2,1,11,11,1,5,11,5,4,11,4,0,11,0,3], + [4,8,5,5,8,11,5,11,2,5,2,0,5,0,9], + [5,4,9,2,3,11], + [4,8,5,8,3,5,3,1,5], + [0,5,4,1,5,0], + [0,9,3,3,9,5,3,5,4,3,4,8], + [5,4,9], + [11,4,7,11,9,4,11,10,9], + [0,3,8,11,4,7,11,9,4,11,10,9], + [0,4,1,4,7,11,1,4,11,1,11,10], + [7,11,4,4,11,10,4,10,1,4,1,3,4,3,8], + [9,4,7,9,7,2,7,3,2,10,9,2], + [8,0,7,7,0,2,7,2,10,7,10,9,7,9,4], + [1,0,10,10,0,4,10,4,7,10,7,3,10,3,2], + [7,8,4,10,1,2], + [9,4,1,1,4,2,2,4,7,2,7,11], + [8,0,3,2,1,9,2,9,4,2,4,11,11,4,7], + [7,11,4,11,2,4,2,0,4], + [3,8,2,2,8,4,2,4,7,2,7,11], + [9,4,1,4,7,1,7,3,1], + [9,4,1,1,4,7,1,7,8,1,8,0], + [3,4,7,0,4,3], + [7,8,4], + [8,11,9,9,11,10], + [0,3,9,3,11,9,11,10,9], + [1,0,10,0,8,10,8,11,10], + [10,3,11,1,3,10], + [3,2,8,2,10,8,10,9,8], + [9,2,10,0,2,9], + [1,0,10,10,0,8,10,8,3,10,3,2], + [2,10,1], + [2,1,11,1,9,11,9,8,11], + [2,1,11,11,1,9,11,9,0,11,0,3], + [11,0,8,2,0,11], + [3,11,2], + [1,8,3,9,8,1], + [1,9,0], + [8,3,0], + [] +]; + +/// Same list as above, but with each row in reverse order. Needed for generating shells (two isosurfaces at slightly different iso values). +/// More efficient just to have a static table than to generate it each time by calling reverse() hundreds of times (although this static table was generated that way). +_MCTriangleTable_reverse = [ + [], + [0,8,3], + [9,0,1], + [3,1,8,8,1,9], + [11,2,3], + [8,11,0,0,11,2], + [11,2,3,9,0,1], + [9,8,11,1,9,11,2,1,11], + [1,2,10], + [8,3,0,10,1,2], + [10,9,2,2,9,0], + [10,9,8,2,10,8,3,2,8], + [11,10,3,3,10,1], + [8,11,10,0,8,10,1,0,10], + [11,10,9,3,11,9,0,3,9], + [11,10,8,8,10,9], + [8,4,7], + [7,3,4,4,3,0], + [7,8,4,1,9,0], + [7,3,1,4,7,1,9,4,1], + [4,7,8,2,3,11], + [2,0,4,11,2,4,7,11,4], + [7,8,4,1,9,0,11,2,3], + [11,2,7,2,1,7,7,1,4,4,1,9], + [2,10,1,8,4,7], + [1,2,10,0,4,3,3,4,7], + [8,4,7,0,2,9,9,2,10], + [2,7,3,2,10,9,2,9,7,9,4,7], + [7,8,4,11,10,3,3,10,1], + [11,10,7,10,0,7,10,1,0,7,0,4], + [11,10,9,3,11,9,0,3,9,8,4,7], + [11,10,9,11,9,4,11,4,7], + [4,9,5], + [5,4,9,3,0,8], + [4,0,5,5,0,1], + [3,1,5,8,3,5,4,8,5], + [4,9,5,11,2,3], + [4,9,5,8,11,0,0,11,2], + [2,3,11,1,5,0,0,5,4], + [1,5,2,5,8,2,5,4,8,2,8,11], + [2,10,1,4,9,5], + [4,9,5,2,10,1,8,3,0], + [4,0,2,5,4,2,10,5,2], + [2,10,3,10,5,3,3,5,8,8,5,4], + [9,5,4,1,3,10,10,3,11], + [8,11,10,0,8,10,1,0,10,9,5,4], + [3,4,0,3,11,10,3,10,4,10,5,4], + [8,11,10,8,10,5,8,5,4], + [8,9,7,7,9,5], + [5,7,3,9,5,3,0,9,3], + [1,5,7,0,1,7,8,0,7], + [5,7,1,1,7,3], + [3,11,2,8,9,7,7,9,5], + [11,5,7,11,2,0,11,0,5,0,9,5], + [1,5,7,0,1,7,8,0,7,3,11,2], + [1,5,7,1,7,11,1,11,2], + [10,1,2,5,7,9,9,7,8], + [5,7,3,9,5,3,0,9,3,1,2,10], + [5,7,10,7,0,10,7,8,0,10,0,2], + [5,7,3,5,3,2,5,2,10], + [3,11,10,3,10,1,7,8,5,8,9,5], + [9,5,0,5,7,0,7,11,0,11,10,0,0,10,1], + [3,11,0,11,10,0,10,5,0,5,7,0,0,7,8], + [5,11,10,7,11,5], + [7,6,11], + [6,11,7,0,8,3], + [6,11,7,9,0,1], + [11,7,6,3,1,8,8,1,9], + [3,7,2,2,7,6], + [6,2,0,7,6,0,8,7,0], + [0,1,9,3,7,2,2,7,6], + [7,9,8,7,6,2,7,2,9,2,1,9], + [7,6,11,1,2,10], + [6,11,7,0,8,3,10,1,2], + [6,11,7,10,9,2,2,9,0], + [10,9,8,2,10,8,3,2,8,11,7,6], + [1,3,7,10,1,7,6,10,7], + [7,6,8,6,10,8,8,10,0,0,10,1], + [3,7,0,7,10,0,7,6,10,0,10,9], + [10,9,8,10,8,7,10,7,6], + [11,8,6,6,8,4], + [0,4,6,3,0,6,11,3,6], + [9,0,1,4,6,8,8,6,11], + [3,1,11,1,4,11,1,9,4,11,4,6], + [4,6,2,8,4,2,3,8,2], + [6,2,4,4,2,0], + [4,6,2,8,4,2,3,8,2,0,1,9], + [4,6,2,4,2,1,4,1,9], + [2,10,1,11,8,6,6,8,4], + [0,4,6,3,0,6,11,3,6,2,10,1], + [6,11,8,6,8,4,2,10,0,10,9,0], + [2,10,3,10,9,3,9,4,3,4,6,3,3,6,11], + [8,1,3,8,4,6,8,6,1,6,10,1], + [0,4,6,0,6,10,0,10,1], + [8,4,3,4,6,3,6,10,3,10,9,3,3,9,0], + [4,10,9,6,10,4], + [11,7,6,9,5,4], + [5,4,9,3,0,8,6,11,7], + [7,6,11,4,0,5,5,0,1], + [3,1,5,8,3,5,4,8,5,7,6,11], + [5,4,9,6,2,7,7,2,3], + [6,2,0,7,6,0,8,7,0,4,9,5], + [7,6,2,7,2,3,5,4,1,4,0,1], + [7,6,8,6,2,8,2,1,8,1,5,8,8,5,4], + [2,10,1,4,9,5,11,7,6], + [10,1,2,8,3,0,6,11,7,4,9,5], + [4,0,2,5,4,2,10,5,2,6,11,7], + [10,5,2,2,5,3,5,4,3,4,8,3,11,7,6], + [1,3,7,10,1,7,6,10,7,5,4,9], + [6,10,7,7,10,8,10,1,8,1,0,8,4,9,5], + [5,4,10,4,0,10,0,3,10,3,7,10,10,7,6], + [7,6,8,6,10,8,10,5,8,8,5,4], + [11,8,9,6,11,9,5,6,9], + [6,11,5,11,3,5,5,3,9,9,3,0], + [6,1,5,6,11,8,6,8,1,8,0,1], + [3,1,5,3,5,6,3,6,11], + [8,9,3,9,6,3,9,5,6,3,6,2], + [6,2,0,6,0,9,6,9,5], + [0,1,8,1,5,8,5,6,8,6,2,8,8,2,3], + [1,6,2,5,6,1], + [11,8,9,6,11,9,5,6,9,10,1,2], + [11,3,6,6,3,5,3,0,5,0,9,5,10,1,2], + [6,11,5,11,8,5,8,0,5,0,2,5,5,2,10], + [6,11,5,11,3,5,3,2,5,5,2,10], + [10,1,6,1,3,6,3,8,6,8,9,6,6,9,5], + [10,1,6,1,0,6,0,9,6,6,9,5], + [6,10,5,0,3,8], + [10,5,6], + [5,10,6], + [10,6,5,0,8,3], + [6,5,10,0,1,9], + [5,10,6,9,8,1,1,8,3], + [3,11,2,5,10,6], + [10,6,5,2,0,11,11,0,8], + [11,2,3,9,0,1,6,5,10], + [9,8,11,1,9,11,2,1,11,10,6,5], + [5,1,6,6,1,2], + [3,0,8,2,6,1,1,6,5], + [0,2,6,9,0,6,5,9,6], + [2,6,3,6,9,3,6,5,9,3,9,8], + [5,1,3,6,5,3,11,6,3], + [0,5,1,0,8,11,0,11,5,11,6,5], + [6,5,11,5,9,11,11,9,3,3,9,0], + [9,8,11,9,11,6,9,6,5], + [10,6,5,8,4,7], + [6,5,10,7,3,4,4,3,0], + [6,5,10,0,1,9,7,8,4], + [7,3,1,4,7,1,9,4,1,5,10,6], + [10,6,5,8,4,7,2,3,11], + [2,0,4,11,2,4,7,11,4,6,5,10], + [6,5,10,0,1,9,11,2,3,8,4,7], + [2,1,11,11,1,7,1,9,7,9,4,7,6,5,10], + [4,7,8,5,1,6,6,1,2], + [1,2,6,1,6,5,3,0,7,0,4,7], + [0,2,6,9,0,6,5,9,6,4,7,8], + [4,7,9,7,3,9,3,2,9,2,6,9,9,6,5], + [5,1,3,6,5,3,11,6,3,7,8,4], + [6,5,11,5,1,11,1,0,11,0,4,11,11,4,7], + [5,9,6,6,9,11,9,0,11,0,3,11,7,8,4], + [4,7,9,7,11,9,11,6,9,9,6,5], + [6,4,10,10,4,9], + [0,8,3,9,10,4,4,10,6], + [6,4,0,10,6,0,1,10,0], + [10,3,1,10,6,4,10,4,3,4,8,3], + [11,2,3,6,4,10,10,4,9], + [10,6,4,10,4,9,11,2,8,2,0,8], + [6,4,0,10,6,0,1,10,0,2,3,11], + [10,6,1,6,4,1,4,8,1,8,11,1,1,11,2], + [2,6,4,1,2,4,9,1,4], + [2,6,4,1,2,4,9,1,4,0,8,3], + [2,6,0,0,6,4], + [2,6,4,2,4,8,2,8,3], + [6,4,11,4,1,11,4,9,1,11,1,3], + [0,8,1,8,11,1,11,6,1,6,4,1,1,4,9], + [6,4,0,6,0,3,6,3,11], + [8,6,4,11,6,8], + [8,9,10,7,8,10,6,7,10], + [7,3,6,3,9,6,3,0,9,6,9,10], + [10,6,1,6,7,1,1,7,0,0,7,8], + [7,3,1,7,1,10,7,10,6], + [8,9,10,7,8,10,6,7,10,11,2,3], + [11,2,7,2,0,7,0,9,7,9,10,7,7,10,6], + [6,7,10,10,7,1,7,8,1,8,0,1,2,3,11], + [11,2,7,2,1,7,1,10,7,7,10,6], + [1,8,9,1,2,6,1,6,8,6,7,8], + [1,2,9,2,6,9,6,7,9,7,3,9,9,3,0], + [0,2,6,0,6,7,0,7,8], + [7,2,6,3,2,7], + [7,8,6,8,9,6,9,1,6,1,3,6,6,3,11], + [0,9,1,6,7,11], + [7,8,6,8,0,6,0,3,6,6,3,11], + [6,7,11], + [7,5,11,11,5,10], + [8,3,0,7,5,11,11,5,10], + [1,9,0,10,11,5,5,11,7], + [5,10,11,5,11,7,1,9,3,9,8,3], + [3,7,5,2,3,5,10,2,5], + [7,5,8,5,2,8,5,10,2,8,2,0], + [3,7,5,2,3,5,10,2,5,1,9,0], + [1,9,2,9,8,2,8,7,2,7,5,2,2,5,10], + [7,5,1,11,7,1,2,11,1], + [7,5,1,11,7,1,2,11,1,3,0,8], + [11,0,2,11,7,5,11,5,0,5,9,0], + [11,7,2,7,5,2,5,9,2,9,8,2,2,8,3], + [3,7,1,1,7,5], + [7,5,1,7,1,0,7,0,8], + [3,7,5,3,5,9,3,9,0], + [9,7,5,8,7,9], + [10,11,8,5,10,8,4,5,8], + [3,10,11,3,0,4,3,4,10,4,5,10], + [10,11,8,5,10,8,4,5,8,9,0,1], + [5,10,4,10,11,4,11,3,4,3,1,4,4,1,9], + [8,4,3,4,5,3,3,5,2,2,5,10], + [2,0,4,2,4,5,2,5,10], + [4,5,8,8,5,3,5,10,3,10,2,3,0,1,9], + [1,9,2,9,4,2,4,5,2,2,5,10], + [11,8,2,8,5,2,8,4,5,2,5,1], + [3,0,11,0,4,11,4,5,11,5,1,11,11,1,2], + [9,0,5,0,2,5,2,11,5,11,8,5,5,8,4], + [11,3,2,9,4,5], + [5,1,3,5,3,8,5,8,4], + [0,5,1,4,5,0], + [8,4,3,4,5,3,5,9,3,3,9,0], + [9,4,5], + [9,10,11,4,9,11,7,4,11], + [9,10,11,4,9,11,7,4,11,8,3,0], + [10,11,1,11,4,1,11,7,4,1,4,0], + [8,3,4,3,1,4,1,10,4,10,11,4,4,11,7], + [2,9,10,2,3,7,2,7,9,7,4,9], + [4,9,7,9,10,7,10,2,7,2,0,7,7,0,8], + [2,3,10,3,7,10,7,4,10,4,0,10,10,0,1], + [2,1,10,4,8,7], + [11,7,2,7,4,2,2,4,1,1,4,9], + [7,4,11,11,4,2,4,9,2,9,1,2,3,0,8], + [4,0,2,4,2,11,4,11,7], + [11,7,2,7,4,2,4,8,2,2,8,3], + [1,3,7,1,7,4,1,4,9], + [0,8,1,8,7,1,7,4,1,1,4,9], + [3,4,0,7,4,3], + [4,8,7], + [10,11,9,9,11,8], + [9,10,11,9,11,3,9,3,0], + [10,11,8,10,8,0,10,0,1], + [10,3,1,11,3,10], + [8,9,10,8,10,2,8,2,3], + [9,2,0,10,2,9], + [2,3,10,3,8,10,8,0,10,10,0,1], + [1,10,2], + [11,8,9,11,9,1,11,1,2], + [3,0,11,0,9,11,9,1,11,11,1,2], + [11,0,2,8,0,11], + [2,11,3], + [1,8,9,3,8,1], + [0,9,1], + [0,3,8], + [] +]; + + +/// _cubindex() - private function, called by _isosurface_cubes() +/// Return the index ID of a voxel depending on the field strength at each corner exceeding isoval. +function _cubeindex(f, isoval) = + (f[0] > isoval ? 1 : 0) + + (f[1] > isoval ? 2 : 0) + + (f[2] > isoval ? 4 : 0) + + (f[3] > isoval ? 8 : 0) + + (f[4] > isoval ? 16 : 0) + + (f[5] > isoval ? 32 : 0) + + (f[6] > isoval ? 64 : 0) + + (f[7] > isoval ? 128 : 0); + + +/// isosurface_cubes() - private function, called by isosurface() +/// This implements a marching cube algorithm, sacrificing some memory in favor of speed. +/// Return a list of voxel cube structures that have one or both surfaces isovalmin or isovalmax intersecting them, and cubes inside the isosurface volume that are at the bounds of the bounding box. +/// The cube structure is: +/// [cubecoord, cubeindex_isomin, cubeindex_isomax, field, bfaces] +/// where +/// cubecoord is the [x,y,z] coordinate of the front left bottom corner of the voxel, +/// cubeindex_isomin and cubeindex_isomax are the index IDs of the voxel corresponding to the min and max iso surface intersections +/// cf is vector containing the 6 field strength values at each corner of the voxel cube +/// bfaces is an array of faces corresponding to the sides of the bounding box - this is empty most of the time; it has data only where the isosurface is clipped by the bounding box. +/// The bounding box 'bbox' is expected to be quantized for the voxel size already. + +function _isosurface_cubes(voxsize, bbox, fieldarray, fieldfunc, isovalmin, isovalmax, closed=true) = let( + // get field intensities + 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]]) + fieldfunc(x,y,z) + ] + ] + ], + nx = len(field)-2, + ny = len(field[0])-2, + nz = len(field[0][0])-2, + v0 = bbox[0] +) [ + for(i=[0:nx]) let(x=v0[0]+voxsize*i) + for(j=[0:ny]) let(y=v0[1]+voxsize*j) + for(k=[0:nz]) let(z=v0[2]+voxsize*k) + let(i1=i+1, j1=j+1, k1=k+1, + cf = [ // cube corner field values + field[i][j][k], + field[i][j][k1], + field[i][j1][k], + field[i][j1][k1], + field[i1][j][k], + field[i1][j][k1], + field[i1][j1][k], + field[i1][j1][k1] + ], + mincf = min(cf), + maxcf = max(cf), + cubecoord = [x,y,z], + bfaces = closed ? _bbox_faces(cubecoord, voxsize, bbox) : [], + cubefound_isomin = (mincf<=isovalmin && isovalmin0 && lenmax>0) let( + + // both min and max surfaces intersect a voxel clipped by bounding box + list = concat( + // min surface + [ for(ei=epathmin) let( + edge = _MCEdgeVertexIndices[ei], + vi0 = edge[0], + vi1 = edge[1], + denom = f[vi1] - f[vi0], + u = abs(denom)<0.0001 ? 0.5 : (isovalmin-f[vi0]) / denom + ) vcube[vi0] + u*(vcube[vi1]-vcube[vi0]) ], + // max surface + [ for(ei=epathmax) let( + edge = _MCEdgeVertexIndices[ei], + vi0 = edge[0], + vi1 = edge[1], + denom = f[vi1] - f[vi0], + u = abs(denom)<0.0001 ? 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( + + // only min surface intersects a voxel clipped by bounding box + list = concat( + [ for(ei=epathmin) let( + edge = _MCEdgeVertexIndices[ei], + vi0 = edge[0], + vi1 = edge[1], + denom = f[vi1] - f[vi0], + u = abs(denom)<0.0001 ? 0.5 : (isovalmin-f[vi0]) / denom + ) vcube[vi0] + u*(vcube[vi1]-vcube[vi0]) ], outfacevertices) + ) for(ls = list) ls + else if(lenmin>0) + + // only min surface intersects a voxel + for(ei=epathmin) let( + edge = _MCEdgeVertexIndices[ei], + vi0 = edge[0], + vi1 = edge[1], + denom = f[vi1] - f[vi0], + u = abs(denom)<0.0001 ? 0.5 : (isovalmin-f[vi0]) / denom + ) vcube[vi0] + u*(vcube[vi1]-vcube[vi0]) + else if(n_outer>0 && lenmax>0) let( + + // only max surface intersects the voxel on the bounding box + list = concat( + [ for(ei=epathmax) let( + edge = _MCEdgeVertexIndices[ei], + vi0 = edge[0], + vi1 = edge[1], + denom = f[vi1] - f[vi0], + u = abs(denom)<0.0001 ? 0.5 : (isovalmax-f[vi0]) / denom + ) vcube[vi0] + u*(vcube[vi1]-vcube[vi0]) ], outfacevertices) + ) for(ls = list) ls + else if(lenmax>0) + + // only max surface intersects the voxel + for(ei=epathmax) let( + edge = _MCEdgeVertexIndices[ei], + vi0 = edge[0], + vi1 = edge[1], + denom = f[vi1] - f[vi0], + u = abs(denom)<0.0001 ? 0.5 : (isovalmax-f[vi0]) / denom + ) vcube[vi0] + u*(vcube[vi1]-vcube[vi0]) + else if(n_outer>0) + + // no surface intersects a voxel clipped by bounding box but the bounding box at this voxel is inside the volume between isomin and isomax + for(ls = outfacevertices) ls +]; + + +/// Generate triangles for the special case of voxel faces clipped by the bounding box +function _bbfacevertices(vcube, f, bbface, isovalmax, isovalmin) = let( + vi = _MCFaceVertexIndices[bbface], + vfc = [ for(i=vi) vcube[i] ], + fld = [ for(i=vi) f[i] ], + pgon = flatten([ + for(i=[0:3]) let( + vi0=vi[i], + vi1=vi[(i+1)%4], + f0 = f[vi0], + f1 = f[vi1], + lowhiorder = (f0=f1) let( + u = abs(denom)<0.0001 ? 0.5 : (isovalmax-f0)/denom + ) vcube[vi0] + u*(vcube[vi1]-vcube[vi0]), + if(fbetweenlow && f0>=f1) let( + u = abs(denom)<0.0001 ? 0.5 : (isovalmin-f0)/denom + ) vcube[vi0] + u*(vcube[vi1]-vcube[vi0]) + + ] + ]), + npgon = len(pgon), + triangles = npgon==0 ? [] : [ + for(i=[1:len(pgon)-2]) [pgon[0], pgon[i], pgon[i+1]] + ]) flatten(triangles); + + +/// _showstats() (Private function) - called by isosurface() and isosurface_array() +/// Display statistics about isosurface +function _showstats(voxelsize, bbox, isoval, cubes, faces) = let( + v = column(cubes, 0), // extract cube vertices + x = column(v,0), // extract x values + y = column(v,1), // extract y values + z = column(v,2), // extract z values + xmin = min(x), + xmax = max(x)+voxelsize, + ymin = min(y), + ymax = max(y)+voxelsize, + zmin = min(z), + zmax = max(z)+voxelsize, + ntri = len(faces), + nvox = len(cubes) +) echo(str("\nIsosurface statistics:\n Outer isovalue = ", isoval, "\n Voxel size = ", voxelsize, + "\n Voxels found containing surface = ", nvox, "\n Triangles = ", ntri, + "\n Voxel bounding box for all data = ", bbox, + "\n Voxel bounding box for isosurface = ", [[xmin,ymin,zmin], [xmax,ymax,zmax]], + "\n")) 0; + + +/// ---------- metaball stuff starts here ---------- + +/// Animated metaball demo made with BOSL2 here: https://imgur.com/a/m29q8Qd + +/// Built-in metaball functions corresponding to each MB_ index. +/// Each function takes three parameters: +/// dv = cartesian distance, a vector [dx,dy,dz] being the distances from the ball center to the volume sample point +/// coeff = the intensity (weight, charge, density, etc.) of the metaball, can be a vector if warranted. +/// additional value or array of values needed by the function. +/// cutoff = radial cutoff distance; effect suppression increases with distance until zero at the cutoff distance, and is zero from that point farther out. Default: INF +/// influence = inverse exponent to 1/r, the higher the influence, the further the "reach" of the metaball. Default: 1 + +/// metaball field function, calling any of the other metaball functions above to accumulate +/// the contribution of each metaball at point xyz + +/* +/// metaball suppression and cutoff function to control shape of the falloff +function mb_shaping(dist, coeff, cutoff, influence) = + sign(influence) * 0.5*(cos(180*(dist/cutoff)^4)+1) + * (coeff/dist)^(1/abs(influence)); + + +function _mb_sphere(dv, coeff, cutoff, influence) = + let(r=norm(dv)) r>=cutoff ? 0 + : mb_shaping(r, coeff, cutoff, influence); + + +//function mb_sphere(coeff=10, cutoff=INF, influence=1) = function (dv) _mb_sphere(dv, coeff, cutoff, influence); +*/ + + +/// metaball cutoff function + +function mb_cutoff(dist, cutoff) = dist>=cutoff ? 0 : 0.5*(cos(180*(dist/cutoff)^4)+1); + +/// metaball sphere + +function _mb_sphere_basic(dv, r, neg) = neg*r/norm(dv); +function _mb_sphere_influence(dv, r, ex, neg) = neg * (r/norm(dv))^ex; +function _mb_sphere_cutoff(dv, r, cutoff, neg) = let(dist=norm(dv)) + neg * mb_cutoff(dist, cutoff) * r/dist; +function _mb_sphere_full(dv, r, ex, cutoff, neg) = let(dist=norm(dv)) + neg * mb_cutoff(dist, cutoff) * (r/dist)^ex; + +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.") + let( + r = get_radius(r=r,d=d), + dummy=assert(is_finite(r) && r>0, "\ninvalid radius or diameter."), + neg = negative ? -1 : 1 + ) + !is_finite(cutoff) && influence==1 ? function(dv) _mb_sphere_basic(dv,r,neg) + : !is_finite(cutoff) ? function(dv) _mb_sphere_influence(dv,r,1/influence, neg) + : influence==1 ? function(dv) _mb_sphere_cutoff(dv,r,cutoff,neg) + : function(dv) _mb_sphere_full(dv,r,1/influence,cutoff,neg); + +// metaball cylinder + +function _mb_cylinder_basic(dv, r, hl, neg) = let( + dist = dv.z<-hl ? norm(dv-[0,0,-hl]) + : dv.z0, "\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 + ) + !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); + +// metaball rounded cube + +function _mb_roundcube_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( + 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( + 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( + 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) = + 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.") + 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); + +// metaball octahedron + +function _mb_octahedron_basic(dv, r, neg) = + let(dist = abs(dv.x) + abs(dv.y) + abs(dv.z)) neg*r/dist; +function _mb_octahedron_influence(dv, r, ex, neg) = + let(dist = abs(dv.x) + abs(dv.y) + abs(dv.z)) neg * (r/dist)^ex; +function _mb_octahedron_cutoff(dv, r, cutoff, neg) = + let(dist = abs(dv.x) + abs(dv.y) + abs(dv.z)) neg * mb_cutoff(dist, cutoff) * r/dist; +function _mb_octahedron_full(dv, r, ex, cutoff, neg) = + let(dist = abs(dv.x) + abs(dv.y) + abs(dv.z)) neg * mb_cutoff(dist, cutoff) * (r/dist)^ex; + +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.") + let( + r = get_radius(r=r,d=d), + dummy=assert(is_finite(r) && r>0, "\ninvalid radius or diameter."), + neg = negative ? -1 : 1 + ) + !is_finite(cutoff) && influence==1 ? function(dv) _mb_octahedron_basic(dv,r,neg) + : !is_finite(cutoff) ? function(dv) _mb_octahedron_influence(dv,r,1/influence, neg) + : influence==1 ? function(dv) _mb_octahedron_cutoff(dv,r,cutoff,neg) + : function(dv) _mb_octahedron_full(dv,r,1/influence,cutoff,neg); + +// torus + +function _mb_torus_basic(dv, rmaj, rmin, neg) = + let(dist = norm([norm([dv.x,dv.y])-rmaj, dv.z])) neg*rmin/dist; +function _mb_torus_influence(dv, rmaj, rmin, ex, neg) = + let(dist = norm([norm([dv.x,dv.y])-rmaj, dv.z])) neg * (rmin/dist)^ex; +function _mb_torus_cutoff(dv, rmaj, rmin, cutoff, neg) = + let(dist = norm([norm([dv.x,dv.y])-rmaj, dv.z])) + neg * mb_cutoff(dist, cutoff) * rmin/dist; +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) = + assert(is_num(cutoff) && cutoff>0, "\ncutoff must be a positive number.") + assert(is_num(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."), + neg = negative ? -1 : 1 + ) + !is_finite(cutoff) && influence==1 ? function(dv) _mb_torus_basic(dv, r_maj, r_min, neg) + : !is_finite(cutoff) ? function(dv) _mb_torus_influence(dv, r_maj, r_min, 1/influence, neg) + : influence==1 ? function(dv) _mb_torus_cutoff(dv, r_maj, r_min, cutoff, neg) + : function(dv) _mb_torus_full(dv, r_maj, r_min, 1/influence, cutoff, neg); + + + +/* hard-edge cylinder +_mb_cylinder = 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&Module: metaballs() +// Synopsis: Creates a model of metaballs within a bounding box. +// SynTags: Geom,VNF +// Topics: Metaballs, Isosurfaces, VNF Generators +// See Also: isosurface_array() +// Usage: As a module +// metaballs(funcs, bounding_box, voxel_size, [isovalue=], [closed=], [show_stats=], ...) [ATTACHMENTS]; +// Usage: As a function +// vnf = metaballs(funcs, bounding_box, voxel_size, [isovalue=], [closed=], [show_stats=]); +// Description: +// [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. +// . +// 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. +// . +// 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. +// . +// 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. +// . +// You can also define your own metaball functions as shown in example 5. +// . +// .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. +// . +// 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 +// . +// 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. +// . +// 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. +// . +// Now for the arguments to this metaball() module or function.... +// 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()`. +// 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 +// --- +// 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 +// 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` +// 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: +// "hull" = Anchors to the virtual convex hull of the shape. +// "intersect" = Anchors to the surface of the shape. +// Named Anchors: +// "origin" = Anchor at the origin, oriented UP. +// 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), +// move([0,0,20]), mb_sphere(5), +// move([20,20,10]), mb_sphere(7) +// ]; +// isovalue = 1; +// voxelsize = 1.5; +// boundingbox = [[-30,-31,-31], [32,31,31]]; +// metaballs(funcs, voxelsize, boundingbox, isovalue); +// Example(3D,NoAxes): A metaball can be negative. In this case we have two metaballs in close proximity, with the small negative metaball creating a dent in the large positive one. The positive metaball is shown transparent, and small spheres show the center of each metaball. The negative metaball isn't visible because its field is negative; the isosurface encloses only field values greater than the isovalue of 1. +// centers = [[-1,0,0], [1.25,0,0]]; +// funcs = [ +// move(centers[0]), mb_sphere(8), +// move(centers[1]), mb_sphere(3, negative=true) +// ]; +// voxelsize = 0.25; +// isovalue = 1; +// boundingbox = [[-7,-6,-6], [3,6,6]]; +// #metaballs(funcs, voxelsize, boundingbox, isovalue); +// 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([10,0,10]), mb_octahedron(5) +// ]; +// voxelsize = 0.5; // a bit slow at this resolution +// boundingbox = [[-12,-9,3], [18,10,32]]; +// metaballs(funcs, voxelsize, boundingbox, isovalue=1); +// Example(3D,NoAxes): Interaction between two tori in different orientations. +// funcs = [ +// move([-10,0,17]), mb_torus(r_maj=6, r_min=2), +// move([7,6,21])*xrot(90), mb_torus(r_maj=7, r_min=3) +// ]; +// voxelsize = 0.5; +// boundingbox = [[-19,-9,9], [18,10,32]]; +// metaballs(funcs, voxelsize, boundingbox, isovalue=1); +// Example(3D,NoAxes,VPD=205,Med): A toy airplane, constructed only from metaball spheres with scaling. The bounding box is used to clip the wingtips, tail, and belly of the fuselage. +// bounding_box = [[-55,-50,-5],[35,50,17]]; +// funcs = [ +// move([-20,0,0])*scale([25,4,4]), mb_sphere(1), // fuselage +// move([30,0,5])*scale([4,0.5,8]), mb_sphere(1), // vertical stabilizer +// move([30,0,0])*scale([4,15,0.5]), mb_sphere(1), // horizontal stabilizer +// move([-15,0,0])*scale([6,45,0.5]), mb_sphere(1) // wing +// ]; +// isovalue = 1; +// voxel_size = 1; +// metaballs(funcs, voxel_size, bounding_box, isovalue); +// Example(3D): Demonstration of a custom metaball function, in this case a sphere with some random noise added to its field potential. The `dv` argument must be first; it is calculated internally as a distance vector from the metaball center to a probe point inside the bounding box, and you convert it to a scalar distance `dist` that is calculated inside your function (`dist` could be a more complicated expression, depending on the shape of the metaball). The call to `mb_cutoff()` at the end handles the cutoff function for the noisy ball consistent with the other internal metaball functions; it requires `dist` and `cutoff` as arguments. You are not required to include the `cutoff` and `influence` arguments in a custom function, but this example shows how. +// function noisy_sphere(dv, r, noise_level, cutoff=INF, influence=1) = +// let( +// noise = rands(0, noise_level, 1)[0], +// dist = norm(dv) + noise +// ) 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), +// ]; +// 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. +// include +// 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]) ]; +// +// funcs = [ +// // vertical cylinder arm +// move([0,0,15]), mb_cylinder(2, 17, 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), +// // ring on top +// move([0,0,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) +// ]; +// 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); + +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); + vnf_polyhedron(vnf, convexity=convexity, cp=cp, anchor=anchor, spin=spin, orient=orient, atype=atype) + children(); +} + +function metaballs(funcs, voxel_size, bounding_box, isovalue=1, closed=true, show_stats=false) = + 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, + // 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)) + ], + + // set up field array + bot = bounding_box[0], + top = bounding_box[1], + halfvox = 0.5*voxel_size, + // 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), + 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 stuff starts here ---------- + +// Function&Module: isosurface() +// Synopsis: Creates a 3D isosurface. +// SynTags: Geom,VNF +// Topics: Isosurfaces, VNF Generators +// Usage: As a module +// isosurface(f, isovalue, voxel_size, bounding_box, [reverse=], [closed=], [show_stats=], ...) [ATTACHMENTS]; +// 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. +// . +// 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. +// . +// 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 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 +// 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 +// --- +// 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 +// 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` +// 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: +// "hull" = Anchors to the virtual convex hull of the shape. +// "intersect" = Anchors to the surface of the shape. +// Named Anchors: +// "origin" = Anchor at the origin, oriented UP. +// Example(3D,ThrownTogether,NoAxes): A gyroid is an isosurface defined by all the zero values of a 3D periodic function. To illustrate what the surface looks like, `closed=false` has been set to expose both sides of the surface. The surface is periodic and tileable along all three axis directions. This a non-manifold surface as displayed, not useful for 3D modeling. This example also demonstrates using an additional parameters in the field function beyond just x,y,z; in this case controls the wavelength of the gyroid. +// function gyroid(x,y,z, wavelength) = let( +// p = 360/wavelength, +// px = p*x, py = p*y, pz = p*z +// ) sin(px)*cos(py) + sin(py)*cos(pz) + sin(pz)*cos(px); +// isovalue = 0; +// bbox = [[-100,-100,-100], [100,100,100]]; +// isosurface(function (x,y,z) gyroid(x,y,z, wavelength=200), +// isovalue, voxel_size=5, bounding_box=bbox, +// closed=false); +// Example(3D,NoAxes): If we remove the `closed` parameter or set it to true, the isosurface algorithm encloses the entire half-space bounded by the "inner" gyroid surface, leaving only the "outer" surface exposed. This is a manifold shape but not what we want if trying to model a gyroid. +// function gyroid(x,y,z, wavelength) = let( +// p = 360/wavelength, +// px = p*x, py = p*y, pz = p*z +// ) sin(px)*cos(py) + sin(py)*cos(pz) + sin(pz)*cos(px); +// isovalue = 0; +// bbox = [[-100,-100,-100], [100,100,100]]; +// isosurface(function (x,y,z) gyroid(x,y,z, wavelength=200), +// isovalue, voxel_size=5, bounding_box=bbox); +// Example(3D,ThrownTogether,NoAxes): To make the gyroid a double-sided surface, we need to specify a small range around zero for `isovalue`. Now we have a double-sided surface although with `closed=false` the edges are not closed where the surface is clipped by the bounding box. +// function gyroid(x,y,z, wavelength) = let( +// p = 360/wavelength, +// px = p*x, py = p*y, pz = p*z +// ) sin(px)*cos(py) + sin(py)*cos(pz) + sin(pz)*cos(px); +// isovalue = [-0.3, 0.3]; +// bbox = [[-100,-100,-100], [100,100,100]]; +// isosurface(function (x,y,z) gyroid(x,y,z, wavelength=200), +// isovalue, voxel_size=5, bounding_box=bbox], +// closed = false); +// Example(3D,ThrownTogether,NoAxes): To make the gyroid a valid manifold 3D object, we remove the `closed` parameter (same as setting `closed=true`), which closes the edges where the surface is clipped by the bounding box. The resulting object can be tiled, the VNF returned by the functional version can be wrapped around an axis using {{vnf_bend()}}, and other operations. +// function gyroid(x,y,z, wavelength) = let( +// p = 360/wavelength, +// px = p*x, py = p*y, pz = p*z +// ) sin(px)*cos(py) + sin(py)*cos(pz) + sin(pz)*cos(px); +// isovalue = [-0.3, 0.3]; +// bbox = [[-100,-100,-100], [100,100,100]]; +// isosurface(function (x,y,z) gyroid(x,y,z, wavelength=200), +// isovalue, voxel_size=5, bounding_box=bbox); +// Example(3D,NoAxes): An approximation of the triply-periodic minimal surface known as [Schwartz P](https://en.wikipedia.org/wiki/Schwarz_minimal_surface). +// function schwartz_p(x,y,z, wavelength) = let( +// p = 360/wavelength, +// px = p*x, py = p*y, pz = p*z +// ) cos(px) + cos(py) + cos(pz); +// isovalue = [-0.2, 0.2]; +// bbox = [[-100,-100,-100], [100,100,100]]; +// isosurface(function (x,y,z) schwartz_p(x,y,z, 100), +// isovalue, voxel_size=4, bounding_box=bbox); +// Example(3D,NoAxes): Another approximation of the triply-periodic minimal surface known as [Neovius](https://en.wikipedia.org/wiki/Neovius_surface). +// function neovius(x,y,z, wavelength) = let( +// p = 360/wavelength, +// px = p*x, py = p*y, pz = p*z +// ) 3*(cos(px) + cos(py) + cos(pz)) + 4*cos(px)*cos(py)*cos(pz); +// isovalue = [-0.3, 0.3]; +// bbox = [[-100,-100,-100], [100,100,100]]; +// isosurface(function (x,y,z) neovius(x,y,z,200), +// isovalue, voxel_size=4, bounding_box=bbox); +// Example(3D): Using an array for the `f` argument instead of a function literal. {{metaballs()}} also makes use of this feature, calculating the the 3D grid first. +// field = [ +// repeat(0,[6,6]), +// [ [0,1,2,2,1,0], +// [1,2,3,3,2,1], +// [2,3,4,4,3,2], +// [2,3,4,4,3,2], +// [1,2,3,3,2,1], +// [0,1,2,2,1,0] +// ], +// [ [0,0,0,0,0,0], +// [0,0,1,1,0,0], +// [0,2,3,3,2,0], +// [0,2,3,3,2,0], +// [0,0,1,1,0,0], +// [0,0,0,0,0,0] +// ], +// [ [0,0,0,0,0,0], +// [0,0,0,0,0,0], +// [0,1,2,2,1,0], +// [0,1,2,2,1,0], +// [0,0,0,0,0,0], +// [0,0,0,0,0,0] +// ], +// repeat(0,[6,6]) +// ]; +// rotate([0,-90,180]) +// isosurface_array(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, bounding_box, voxel_size, reverse, closed, show_stats, _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) = + 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, + 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=f, 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 +) [trianglepoints, faces]; +>>>>>>> Stashed changes