Download or view Plot3D.frink in plain text format
/*
This is a simple but rather interesting program that graphs equations in
3 dimensions.
You enter equations in terms of x, y, and z something like one of the
following:
y = sin[x] + cos[z]
x^2 + y^2 + z^2 = 81 (This is a sphere)
y cos[x] = x sin[z]
This version of the program can also graph INEQUALITIES, which have
less-than or greater-than symbols instead of just equals.
Inequalities are important for graphing infinitely-thin objects and making
them printable and sliceable. For example, you might need to convert:
x^2 + y^2 - z^2 = 81
into
abs[x^2 + y^2 - z^2 - 81] <= 2
to give the walls some thickness and give your slicer a chance to print it
successfully. Also increasing the value of "res" below will make the
voxels in the .obj file larger.
Here is an egg. Modify the 1.5 and 1.1 to change its aspect ratio:
(x^2 + y^2 + z^2)^2 <= -6 (1.5 z^3 + (1.5 - 1.1) z (y^2 + x^2))
Here is a Klein bottle:
(x^2 + y^2 + z^2 + 2y - 1)((x^2 + y^2 + z^2 - 2y - 1)^2 - 8 z^2) + 16 x z (x^2 + y^2 + z^2 - 2y - 1) = 0
Here is a solid heart:
320 (-x^2 z^3 - 9 y^2 z^3 / 80 + (x^2 + 9 y^2 / 4 + z^2 - 1)^3) <= 0
This uses a recursive method to subdivide and test cuboids.
You can also use logical relations like AND and OR to combine multiple
shapes.
*/
class Plot3D
{
var xmin = -10
var xmax = 10
var ymin = -10
var ymax = 10
var zmin = -10
var zmax = 10
// Change the doublings to vary the number of voxels. This is the number
// of doublings, so if the number is 10 we have 2^10=1024 doublings for
// a resolution of 1024x1024x1024. (That is over a billion pixels! Don't
// be surprised that graphing at that resolution takes a long time!)
// Be warned that increasing the doublings by 1 makes 8 times as many voxels!
var doublings = 8
// Resolution at which to render the .obj file
var resolution = 254/in
/** Plots the function (specified as a string) and returns a VoxelArray. */
plot[func] :=
{
hasInequality = false
certEq = undef
certFunc = func
// If there's an inequality, let's make a test equation to see if we can
// fill in an entire cuboid using the "CERTAINLY" comparators.
if func =~ %r/([<>]|!=)/
{
hasInequality = true
certFunc =~ %s/<=/ CLE /g // Replace <= with certainly less than or equals
certFunc =~ %s/>=/ CGE /g // Replace >= with certainly greater than or equals
certFunc =~ %s/</ CLT /g // Replace < with certainly less than
certFunc =~ %s/>/ CGT /g // Replace > with certainly greater than
certFunc =~ %s/!=/ CNE /g // Replace = with certainly not equals
certFunc =~ %s/==/ CEQ /g // Replace == with certainly equals
certFunc =~ %s/=/ CEQ /g // Replace = with certainly equals
certEq = parseToExpression[certFunc]
}
// These replacements turn normal comparator and equality tests into
// "POSSIBLY EQUALS" tests.
func =~ %s/<=/ PLE /g // Replace <= with possibly less than or equals
func =~ %s/>=/ PGE /g // Replace >= with possibly greater than or equals
func =~ %s/</ PLT /g // Replace < with possibly less than
func =~ %s/>/ PGT /g // Replace > with possibly greater than
func =~ %s/!=/ PNE /g // Replace = with possibly not equals
func =~ %s/==/ PEQ /g // Replace == with possibly equals
func =~ %s/=/ PEQ /g // Replace = with possibly equals
eq = parseToExpression[func]
numVoxels = 2^doublings
// Scale factors on each axis
sx = numVoxels / (xmax-xmin)
sy = numVoxels / (ymax-ymin)
sz = numVoxels / (zmax-zmin)
v = callJava["frink.graphics.VoxelArray", "construct", [xmin sx, xmax sx, ymin sy, ymax sy, zmin sz, zmax sz, false]]
// Perform the recursive testing of the volume
testCube[xmin, xmax, ymin, ymax, zmin, zmax, v, eq, certEq, doublings, sx, sy, sz]
return v
}
/** Shows an already-rendered VoxelArray on the screen. */
show[v] :=
{
v.projectX[undef].show["X"]
v.projectY[undef].show["Y"]
v.projectZ[undef].show["Z"]
}
getSetBounds[v] :=
{
// To convert from voxel coordinate v to original coordinate x,
// x = xmin + v (xmax - xmin) / (vmax - vmin)
// inverse:
// v = x (vmax - vmin) / (xmax - xmin)
// println["Center of mass (in voxel coords): " + v.centerOfMass[].toString[]]
min = v.getMinimumSetPoints[]
max = v.getMaximumSetPoints[]
println["min is " + min.toString[]]
println["max is " + max.toString[]]
minx1 = min.x (xmax-xmin) / (v.getMaxX[] - v.getMinX[])
maxx1 = max.x (xmax-xmin) / (v.getMaxX[] - v.getMinX[])
miny1 = min.y (ymax-ymin) / (v.getMaxY[] - v.getMinY[])
maxy1 = max.y (ymax-ymin) / (v.getMaxY[] - v.getMinY[])
minz1 = min.z (zmax-zmin) / (v.getMaxZ[] - v.getMinZ[])
maxz1 = max.z (zmax-zmin) / (v.getMaxZ[] - v.getMinZ[])
return [minx1, maxx1, miny1, maxy1, minz1, maxz1]
}
/** Write a rendered VoxelArray to a Wavefront .obj file. */
writeObj[v, name] :=
{
var shortname = name
if [shortname] = name =~ %r/(.*)\.obj$/i
filename = name
else
filename = "$name.obj"
print["Writing $filename..."]
w = new Writer[filename]
w.println[v.toObjFormat[shortname, 1 / (resolution mm)]]
w.close[]
println["done."]
}
// Recursive function to test an interval containing the specified bounds.
// If no possible solution exists, the recursion halts. If the entire cube
// is filled with a "certainly" equation, it is filled and recursion halts.
// If only a possible solution exists, this breaks it down into 8 sub-cuboids
// and tests each of them recursively.
// level is the maximum number of levels to split, so the total
// resolution of the final graph will be 2^level.
testCube[x1, x2, y1, y2, z1, z2, v, eq, certEq, level, sx, sy, sz] :=
{
nextLevel = level - 1
x = new interval[x1, x2]
y = new interval[y1, y2]
z = new interval[z1, z2]
// Test the cuboid. If it possibly contains solutions, recursively
// subdivide.
result = eval[eq]
// If there is a possible solution or an error
if result or result==undef
{
// If we should still recurse
if (nextLevel >= 0)
{
// Do we have inequalities and a CERTAINLY test?
if (certEq != undef) AND (eval[certEq] == true)
{
// If the entire cuboid is a solution, then fill the rectangle
// and stop further recursion on this cuboid.
v.fillCube[x1 sx, x2 sx, y1 sy, y2 sy, z1 sz, z2 sz, true]
return
}
// Further subdivide the cuboid into 8 octants and recursively
// test them all
cx = (x1 + x2)/2
cy = (y1 + y2)/2
cz = (z1 + z2)/2
testCube[x1, cx, y1, cy, z1, cz, v, eq, certEq, nextLevel, sx,sy,sz]
testCube[cx, x2, y1, cy, z1, cz, v, eq, certEq, nextLevel, sx,sy,sz]
testCube[x1, cx, cy, y2, z1, cz, v, eq, certEq, nextLevel, sx,sy,sz]
testCube[cx, x2, cy, y2, z1, cz, v, eq, certEq, nextLevel, sx,sy,sz]
testCube[x1, cx, y1, cy, cz, z2, v, eq, certEq, nextLevel, sx,sy,sz]
testCube[cx, x2, y1, cy, cz, z2, v, eq, certEq, nextLevel, sx,sy,sz]
testCube[x1, cx, cy, y2, cz, z2, v, eq, certEq, nextLevel, sx,sy,sz]
testCube[cx, x2, cy, y2, cz, z2, v, eq, certEq, nextLevel, sx,sy,sz]
} else
if (result) // Valid point at lowest level; fill it
v.fillCube[x1 sx, x2 sx, y1 sy, y2 sy, z1 sz, z2 sz,true]
}
}
/** Sets the bounds to be plotted. */
setBounds[xmin, xmax, ymin, ymax, zmin, zmax] :=
{
this.xmin = xmin
this.xmax = xmax
this.ymin = ymin
this.ymax = ymax
this.zmin = zmin
this.zmax = zmax
}
/** Sets the voxel resolution for printing in either
pixels (dimensionless) / length
or
length.
for example:
setResolution[254/in]
*/
setResolution[res] :=
{
if res conforms (1/mm)
{
resolution = res
return
}
if res conforms (mm)
{
resolution = 1/res
return
}
println["Plot3D.setResolution: Invalid resolution $res"]
}
/** Sets the number of doublings. The number of voxels rendered will be
2^doublings on each side. */
setDoublings[d] :=
{
doublings = d
}
}
Download or view Plot3D.frink in plain text format
This is a program written in the programming language Frink.
For more information, view the Frink
Documentation or see More Sample Frink Programs.
Alan Eliasen was born 19975 days, 5 hours, 28 minutes ago.