isosurface.scad
An isosurface is a three-dimensional surface representing points of a constant value (e.g. density pressure, temperature, electric field strength, density) in a 3D volume. It is essentially a 3D cross-section of a 4-dimensional function. An isosurface may be represented generally by any function of three variables, that is, the function returns a single value based on [x,y,z] inputs. The isosurface is defined by all return values equal to a constant isovalue.
A gryoid (often used as a volume infill pattern in FDM 3D printing) is an exmaple of an isosurface that is unbounded and periodic in all three dimensions. Other typical examples in 3D graphics are metaballs (also known as "blobby objects"), which are bounded and closed organic-looking surfaces that meld together when in close proximity.
To use, add the following lines to the beginning of your file:
include <BOSL2/std.scad>
include <BOSL2/isosurface.scad>
-
isosurface()– Creates a 3D isosurface. [Geom] [VNF] -
isosurface_array()– Creates a 3D isosurface from a 3D array of densities. [Geom] [VNF] -
metaballs()– Creates a model of metaballs within a bounding box. [Geom] [VNF]
Synopsis: Creates a 3D isosurface. [Geom] [VNF]
Topics: Isosurfaces, VNF Generators
Usage: As a module
- isosurface(voxel_size, bounding_box, isovalue, field_function, [additional=], [reverse=], [close_clip=], [show_stats=]);
Usage: As a function
- vnf = isosurface(voxel_size, bounding_box, isovalue, field_function, [additional=], [close_clip=], [show_stats=]);
Description:
When called as a function, returns a VNF structure (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 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:
| By Position | What it does |
|---|---|
voxel_size |
The size (scalar) of the voxel cube that determines the resolution of the surface. |
bounding_box |
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. |
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. |
field_function |
A function literal taking as input an [x,y,z] coordinate and optional additional parameters, and returns a single value. |
| By Name | What it does |
|---|---|
additional |
A single value, or an array of optional additional parameters that may be required by the field function. It is your responsibility to create a function literal compatible with these inputs. If additional is not set, only the [x,y,z] parameter is passed to the function; no additional parameters are passed. Default: undef |
reverse |
When true, reverses the orientation of the facets in the mesh. Default: false |
close_clip |
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 |
Example 1: A gyroid is an isosurface defined by all the zero values of a 3D periodic function. To illustrate what the surface looks like, close_clip=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 the use of the additional parameter, which in this case controls the wavelength of the gyroid.
include <BOSL2/std.scad>
include <BOSL2/isosurface.scad>
gyroid = function (xyz, wavelength) let(
p = 360/wavelength,
px = p*xyz[0],
py = p*xyz[1],
pz = p*xyz[2]
) sin(px)*cos(py) + sin(py)*cos(pz) + sin(pz)*cos(px);
bbox = [[-100,-100,-100], [100,100,100]];
isosurface(voxel_size=5, bounding_box=bbox, isovalue=0,
field_function=gyroid, additional=200, close_clip=false);
Example 2: If we remove the close_clip 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.
include <BOSL2/std.scad>
include <BOSL2/isosurface.scad>
gyroid = function (xyz, wavelength) let(
p = 360/wavelength,
px = p*xyz[0],
py = p*xyz[1],
pz = p*xyz[2]
) sin(px)*cos(py) + sin(py)*cos(pz) + sin(pz)*cos(px);
bbox = [[-100,-100,-100], [100,100,100]];
isosurface(voxel_size=5, bounding_box=bbox, isovalue=0,
field_function=gyroid, additional=200);
Example 3: 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 close_clip=false the edges are not closed where the surface is clipped by the bounding box.
include <BOSL2/std.scad>
include <BOSL2/isosurface.scad>
gyroid = function (xyz, wavelength) let(
p = 360/wavelength,
px = p*xyz[0],
py = p*xyz[1],
pz = p*xyz[2]
) sin(px)*cos(py) + sin(py)*cos(pz) + sin(pz)*cos(px);
bbox = [[-100,-100,-100], [100,100,100]];
isosurface(voxel_size=5, bounding_box=bbox, isovalue=[-0.3, 0.3],
field_function=gyroid, additional=200, close_clip=false);
Example 4: To make the gyroid a valid manifold 3D object, we remove the close_clip parameter (same as setting close_clip=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.
include <BOSL2/std.scad>
include <BOSL2/isosurface.scad>
gyroid = function (xyz, wavelength) let(
p = 360/wavelength,
px = p*xyz[0],
py = p*xyz[1],
pz = p*xyz[2]
) sin(px)*cos(py) + sin(py)*cos(pz) + sin(pz)*cos(px);
bbox = [[-100,-100,-100], [100,100,100]];
isosurface(voxel_size=5, bounding_box=bbox, isovalue=[-0.3, 0.3],
field_function=gyroid, additional=200);
Example 5: An approximation of the triply-periodic minimal surface known as Schwartz P.
include <BOSL2/std.scad>
include <BOSL2/isosurface.scad>
schwartz_p = function (xyz, wavelength) let(
p = 360/wavelength,
px = p*xyz[0],
py = p*xyz[1],
pz = p*xyz[2]
) cos(px) + cos(py) + cos(pz);
bbox = [[-100,-100,-100], [100,100,100]];
isosurface(voxel_size=4, bounding_box=bbox, isovalue=[-0.2,0.2],
field_function=schwartz_p, additional=100);
Example 6: Another approximation of the triply-periodic minimal surface known as Neovius.
include <BOSL2/std.scad>
include <BOSL2/isosurface.scad>
neovius = function (xyz, wavelength) let(
p = 360/wavelength,
px = p*xyz[0],
py = p*xyz[1],
pz = p*xyz[2]
) 3*(cos(px) + cos(py) + cos(pz)) + 4*cos(px)*cos(py)*cos(pz);
bbox = [[-100,-100,-100], [100,100,100]];
isosurface(voxel_size=4, bounding_box=bbox, isovalue=[-0.3,0.3],
field_function=neovius, additional=200);
Synopsis: Creates a 3D isosurface from a 3D array of densities. [Geom] [VNF]
Topics: Isosurfaces, VNF Generators
Usage: As a module
- isosurface_array(voxel_size, isovalue, fields, [origin=], [reverse=], [close_clip=], [show_stats=]);
Usage: As a function
- vnf = isosurface_array(voxel_size, isovalue, fields, [origin=], [reverse=], [close_clip=], [show_stats=]);
Description:
When called as a function, returns a VNF structure (list of triangles and
faces) representing a 3D isosurface within the passed array at a single isovalue or
range of isovalues.
When called as a module, displays the isosurface within the passed array 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.
Use this when you already have a 3D array of intensity or density data, for example like what you may get from a CT scan.
By default, the returned VNF structure occupies a volume with its origin at [0,0,0]
extending in the positive x, y, and z directions by multiples of voxel_size.
This origin can be overridden by the origin parameter.
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 at the outer limits
of the fields array 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:
| By Position | What it does |
|---|---|
voxel_size |
The size (scalar) of the voxel cube that determines the resolution of the surface. |
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 surfaces (such as gyroids) that have both sides exposed within the bounding box. |
fields |
3D array of field intesities. This array should be organized so that the indices are in order of x, y, and z when the array is referenced; that is, fields[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. |
| By Name | What it does |
|---|---|
origin |
Origin in 3D space corresponding to fields[0][0][0]. The bounding box of the isosurface extends from this origin by multiples of voxel_size according to the size of the fields array. Default: [0,0,0] |
reverse |
When true, reverses the orientation of the facets in the mesh. Default: false |
close_clip |
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, exposes 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 |
Example 1:
include <BOSL2/std.scad>
include <BOSL2/isosurface.scad>
fields = [
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(voxel_size=10,
isovalue=0.5, fields=fields);
Synopsis: Creates a model of metaballs within a bounding box. [Geom] [VNF]
Topics: Metaballs, Isosurfaces, VNF Generators
See Also: isosurface_array()
Usage: As a module
- metaballs(voxel_size, bounding_box, isovalue, ball_centers, [ball_sizes=], [ball_type=], [rotation=], [field_function=], [radial_cutoff=], [close_clip=], [show_stats=]);
Usage: As a function
- vnf = metaballs(voxel_size, bounding_box, isovalue, ball_centers, [ball_sizes=], [ball_type=], [rotation=], [field_function=], [radial_cutoff=], [close_clip=], [show_stats=]);
Description:
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 "charge" of each ball and their distance from one another. If you consider a "ball" to be a point charge in 3D space, the electric field surrounding that charge decreases in intensity with distance from the charge. The metaball is the isosurface corresponding to all value where the electric field intensity is a constant value. A stronger charge results in a stronger the electric field, and correspondingly a larger metaball. Fields from two charges add together, changing the shape of the two corresponding metaballs when they are in close proximity.
In physics, the electric field intensity falls off as an inverse-square relationship
with distance; that is, the field is proportional to field_function parameter.
Six shapes of fields around each metaball center are possible. You can specify different types for each metaball in the list, and you can also specify your own custom field equation. The five types are:
-
MB_SPHERE- the standard spherical metaball with a 1/r field strength falloff. -
MB_ELLIPSOID- an ellipsoid-shaped field that requires specifying a [x,y,z] vector for the charge, representing field strength in each of the x, y, and z directions -
MB_ROUNDCUBE- a cube-shaped metaball with corners that get more rounded with size. The squareness can be controlled with a value between 0 (spherical) or 1 (cubical) in theadditionalparameter, and defaults to 0.5 if omitted. -
MB_CUBE- a cube-shaped metaball with sharp edges and corners, resulting from using Chebyshev distance rather than Euclidean distance calculations. -
MB_OCTAHEDRON- an octahedron-shaped metaball with sharp edges and corners, resulting from using taxicab distance rather than Euclidean distance calculations. -
MB_TORUS- a toroidal field oriented perpendicular to the x, y, or z axis. Thechargeis a two-element vector determining the major and minor diameters, and theadditionalparamater sets the axis directions for each ball center (defaults to [0,0,1] if not set). -
MB_CUSTOM- your own custom field definition, requiring you to set thefield_functionparameter to your own function literal. If eitherMB_ELLIPSOIDorMB_TORUSoccur in the list, the list of charges must be explicitly defined rather than supplying a single value for all.
Arguments:
| By Position | What it does |
|---|---|
voxel_size |
The size (scalar) of the voxel cube that determines the resolution of the metaball surface. |
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 of the metaballs. |
ball_centers |
an array of 3D points specifying each of the metaball centers. |
| By Name | What it does |
|---|---|
charge |
a single value, or an array of values corresponding to ball_centers, specifying the charge intensity of each ball center. Default: 10 |
ball_type |
shape of field that falls off from the metaball center. Can be one of MB_SPHERE, MB_ELLIPSOID, MB_ROUNDCUBE, MB_CUBE, MB_OCTAHEDRON, MB_TORUS, or MB_CUSTOM. This may be an array of values corresponding to each ball. Where this value is MB_CUSTOM, the corresponding array element in field_function must also be set. Default: _MB_SPHERE
|
rotation |
A vector [x_rotation, y_rotation, z_rotation], or list of vectors for each ball, specifying the rotation angle in degrees around the x, y, and z axes. This is meaningless for _MB_SPHERE but allows you to orient the other metaball types. Default: undef |
field_function |
A single function literal or array of function literals that return a single field value from one metaball, and takes as inputs a 3D distance vector, a single charge or list of charges, and a single additional parameter or list of parameters (that third parameter must exist in the function definition even if it isn't used). If the corresponding ball_type parameter is not MB_CUSTOM, then the function specified in ball_type is used instead; only where ball_type is MB_CUSTOM does this custom field function get invoked. Default: undef |
additional |
A single value, or a list of optional additional parameters that may be required by the field function. If you make a custom function, it is your responsibility to create a function literal compatible with these inputs. Nothing is passed to the function literal if additional is not set. This parameter must be specified as an entire list for all metaballs if MB_ELLIPSOID or MB_TORUS is included in ball_type. Default: undef for ball_type=CUSTOM
|
radial_cutoff |
Maximum radial distance of a metaball's influence. This isn't a sharp cutoff; rather, the suppression increases with distance until the influence is zero at the radial_cutoff distance. Can be a single value or an array of values corresponding to each ball center, but typically it's sufficient to supply a single value approximately the average separation of each ball, so each ball mostly acts on its nearest neighbors. Default: INF |
close_clip |
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 |
Example 1: A group of five spherical metaballs with different charges. The parameter show_stats=true (not shown here) was used to find a compact bounding box for this figure.
include <BOSL2/std.scad>
include <BOSL2/isosurface.scad>
centers = [[-20,-20,-20], [-0,-20,-20],
[0,0,0], [0,0,20], [20,20,10] ];
charges = [5, 4, 3, 5, 7];
type = MB_SPHERE;
isovalue = 1;
voxelsize = 1.5;
boundingbox = [[-30,-31,-31], [32,31,31]];
metaballs(voxelsize, boundingbox, isovalue=isovalue,
ball_centers=centers, charge=charges, ball_type=type);
Example 2: A metaball can have negative charge. 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.
include <BOSL2/std.scad>
include <BOSL2/isosurface.scad>
centers = [[-1,0,0], [1.25,0,0]];
charges = [8, -3];
type = MB_SPHERE;
voxelsize = 0.25;
isovalue = 1;
boundingbox = [[-7,-6,-6], [3,6,6]];
#metaballs(voxelsize, boundingbox, isovalue=isovalue,
ball_centers=centers, charge=charges, ball_type=type);
color("green") for(c=centers) translate(c) sphere(d=1, $fn=16);
Example 3: A cube, a rounded cube, and an octahedron interacting.
include <BOSL2/std.scad>
include <BOSL2/isosurface.scad>
centers = [[-7,-3,27], [7,5,21], [10,0,10]];
charge = 5;
type = [MB_CUBE, MB_ROUNDCUBE, MB_OCTAHEDRON];
voxelsize = 0.4; // a bit slow at this resolution
isovalue = 1;
boundingbox = [[-13,-9,3], [16,11,33]];
metaballs(voxelsize, boundingbox, isovalue=isovalue,
ball_centers=centers, charge=charge, ball_type=type);
Example 4: Interaction between two torus-shaped fields in different orientations.
include <BOSL2/std.scad>
include <BOSL2/isosurface.scad>
centers = [[-10,0,17], [7,6,21]];
charges = [[6,2], [7,3]];
type = MB_TORUS;
axis_orient = [[0,0,1], [0,1,0]];
voxelsize = 0.5;
isovalue = 1;
boundingbox = [[-19,-9,9], [18,10,32]];
metaballs(voxelsize, boundingbox, isovalue=isovalue,
ball_centers=centers, charge=charges, ball_type=type,
additional=axis_orient);
Example 5: Demonstration of a custom metaball function, in this case a sphere with some random noise added to its electric field.
include <BOSL2/std.scad>
include <BOSL2/isosurface.scad>
noisy_sphere = function (cdist, charge, additional,
rotation_matrix_unused, rcutoff=INF)
let(
r = norm(cdist) + rands(0, 0.2, 1)[0],
suppress = let(a=min(r,rcutoff)/rcutoff) 1-a*a
) r==0 ? 1000*charge : suppress*charge / r;
centers = [[-9,0,0], [9,0,0]];
charge = 5;
type = [MB_SPHERE, MB_CUSTOM];
fieldfuncs = [undef, noisy_sphere];
voxelsize = 0.4;
boundingbox = [[-16,-8,-8], [16,8,8]];
metaballs(voxelsize, boundingbox, isovalue=1,
ball_centers=centers, charge=charge, ball_type=type,
field_function=fieldfuncs);
Example 6: 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 <BOSL2/std.scad>
include <BOSL2/isosurface.scad>
ztheta = 90-acos(-1/3);
cz = cos(ztheta);
sz = sin(ztheta);
type = [
MB_ELLIPSOID, MB_ELLIPSOID,
MB_ELLIPSOID, MB_ELLIPSOID,
MB_TORUS, MB_SPHERE, MB_SPHERE, MB_SPHERE
];
centers = [
[0,0,20], [20*cz,0,20*sz],
zrot(120, p=[20*cz,0,20*sz]),
zrot(-120, p=[20*cz,0,20*sz]),
[0,0,35], [32*cz,0,32*sz],
zrot(120, p=[32*cz,0,32*sz]),
zrot(-120, p=[32*cz,0,32*sz])];
cutoff = 40; // extent of influence of each ball
rotation = [
[0,90,0], [0,-ztheta,0], [0,-ztheta,120], [0,-ztheta,-120],
[0,0,0], undef, undef, undef];
axis = [
undef, undef, undef, undef,
[0,1,0], undef, undef, undef
];
charge = [
[6,2,2], [7,2,2], [7,2,2], [7,2,2],
[8,2], 5, 5, 5
];
voxelsize = 1;
isovalue = 1;
boundingbox = [[-23,-36,-15], [39,36,46]];
// useful to save as VNF for copies and manipulations
vnf = metaballs(voxelsize, boundingbox, isovalue=isovalue, ball_centers=centers,
charge=charge, ball_type=type, additional=axis, rotation=rotation,
radial_cutoff=cutoff);
vnf_polyhedron(vnf);