Function Based 3-D Solid Modeller
Aaron Schilling, Fareed Faghih, Marcus Lew
Math480B Final Project -- 6/4/2010
3D Model of a Jack in a Cube
3D Modelling Software in action
Self-Replicating RepRap 3D Printer
Commercial Z406 3D printer
Example 3D scan and 3D print
{{{id=1|
%auto
#Define Variables
x,y,z = var('x, y, z')
#Transformations: Mathematical functions to change the geometrical size or location of solids. Includes translations, scaling, reflection and rotation.
#Translations: Movement in the specific direction
def TransX(a,tx) : return a.substitute(x=(x-tx))
def TransY(a,ty) : return a.substitute(y=(y-ty))
def TransZ(a,tz) : return a.substitute(z=(z-tz))
#Scaling: Size factor change in specific direction
def ScaleX(a,sx) : return a.substitute(x=(x/sx))
def ScaleY(a,sy) : return a.substitute(y=(y/sy))
def ScaleZ(a,sz) : return a.substitute(z=(z/sz))
def ScaleXYZ(a,s) : return a.substitute(x=(x/s),y=(y/s),z=(z/s))
#Reflection: Reflects across specific plane
def ReflectXY(a) : return a.substitute(z=(-z))
def ReflectZX(a) : return a.substitute(y=(-y))
def ReflectYZ(a) : return a.substitute(x=(-x))
#Rotation: Turning in specific plane
def DegToRad(a) : return (pi*a/180)
def RotateRX(a,r) : return a.substitute(y=(y*cos(r) - z*sin(r)), z=(y*sin(r) + z*cos(r)))
def RotateRY(a,r) : return a.substitute(x=(x*cos(r) + z*sin(r)), z=(-x*sin(r) + z*cos(r)))
def RotateRZ(a,r) : return a.substitute(x=(x*cos(r) + y*sin(r)), y=(-x*sin(r) + y*cos(r)))
def RotateX(a,d) : return RotateRX(a, DegToRad(d))
def RotateY(a,d) : return RotateRY(a, DegToRad(d))
def RotateZ(a,d) : return RotateRZ(a, DegToRad(d))
# Simpler Plot Function: Reduces implicit_plot3d to a simpler more practical plotter
def Plot3d (sld,size,pts) : return implicit_plot3d(sld, (x, -size, size), (y, -size, size), (z, -size, size), plot_points=pts)
# Boolean Operators: For shape interactions
def Union(a,b) : return Min(a,b)
def Intersect(a,b) : return Max(a,b)
def Subtract(a,b) : return Max(a,-b)
def Morph (a,b,t) : return ((b-a)*t+a)
def Blend (a,b,l) : return (a*b)-l
def Combine (a,b,k) : return k*(a+b)
Max(a,b) = (a + b + abs(a - b)) / 2
Min(a,b) = (a + b - abs(a - b)) / 2
# Primitives and Shapes: Defines the specific building blocks which can be pieced together or transformed to form the desired object. Adapted from w3dsurf and wolfram mathworld
def Sphere (a) : return (x^2+y^2+z^2-a^2)
def Cube(h) : return (Max(Max(abs(x/(h/2)), abs(y/(h/2))), abs(z/(h/2))) - 1)
def Prism(a,b,c) : return (Max(Max(abs(x/(a/2)), abs(y/(b/2))), abs(z/(c/2))) - 1)
def Torus(a,b) : return ((a-sqrt(x^2+y^2)^2+z^2-b^2))
def Cyyl(a,b,c) : return Intersect(((x^2/a)+(y^2/b)-1),Close(z,c))
def Cone(a,b,c) : return Intersect(((c*sqrt((x^2/a^2)+(y^2/b^2)))-z),Close(z,c))
def Ellip(a,b,c) : return ((x^2/a^2)+(y^2/b^2)+(z^2/c^2)-1)
def Para (a,b,c) : return (Max((x^2/a^2+y^2/b^2-c*z),z)-1)
def Close(a,b) : return (abs(a/(b/2))-1)
def HSpace (h) : return z-h
def Tetrahedral (a,b,c) : return ((x^2 + y^2 + z^2)^2 + a*x*y*z - b*(x^2 + y^2 + z^2) + c)
def Sinus (r,h) : return (sin(pi*(x^2+y^2-r^2))/2+h*z)
def Bug (a,b,c) : return ((x*cos(b*y)-z*sin(b*y))+a)^2+(y/c)^2+(x*sin(b*y)+z*cos(b*y))^2 - 1
def Drope (h) : return (z - h*x*exp(-x^2-y^2))
def Heart (a,b,c,d,e) : return (a*(x^2)+b*(y^2)+c*(z^2)-1)^3-d*(x^2)*(z^3)-e*(y^2)*(z^3)
def surface_to_stl_filesave(surface,filename):
# Return an STL representation of the surface to a file.
# Original work by Christopher Olah and later modified to write to a
# file.
#
# INPUT:
# - `surface` -- any surface, eg. output of a 3d plot function.
#
# OUTPUT:
# A string that represents the surface in the STL format.
#
# COMMENTS:
# (1) You must view the surface before plotting it.
# Otherwise, this will not work.
# (2) In order to do 3d printing with this, you will need to
# convert it into gcode. Skeinforge is an open-source
# program that can do this.
# (3) The size of the surface is not normalized in export.
# Sage's units will become the units of the STL
# description. These seem to be ~0.05 cm (at least when
# printed using skeinforge -> replicatorg -> hacklab.to's
# cupcake).
# (4) Be aware of the overhang limits of your 3d printer;
# most printers can only handle an overhang of Pi/2 (45*)
# before your model will start drooping.
#
# EXAMPLES:
# sage: x,y,z = var('x,y,z')
# sage: a = implicit_plot3d(x^2+y^2+z^2-9, [x,-5,5], [y,-5,5],[z,
# -5,5])
# sage: surface_to_stl_filesave(a,'file')
#
stl=open(filename + ".stl",'w');
stl.write('solid mathsurface\n');
for i in surface.face_list():
n = ( i[1][1]*i[2][2] - i[2][1]*i[1][2],
i[1][2]*i[2][0] - i[1][0]*i[2][2],
i[1][0]*i[2][1] - i[2][0]*i[1][1] )
abs = (n[0]^2+n[1]^2+n[2]^2)^0.5
n= (n[0]/abs,n[1]/abs,n[2]/abs)
stl.write(' facet normal ' + repr(n[0]) + ' ' + repr(n[1]) + ' ' + repr(n[2]) + '\n');
stl.write(' outer loop\n');
stl.write(' vertex ' + repr(i[0][0]) + ' ' + repr(i[0][1]) + ' ' + repr(i[0][2]) + '\n');
stl.write(' vertex ' + repr(i[1][0]) + ' ' + repr(i[1][1]) + ' ' + repr(i[1][2]) + '\n');
stl.write(' vertex ' + repr(i[2][0]) + ' ' + repr(i[2][1]) + ' ' + repr(i[2][2]) + '\n');
stl.write(' endloop\n');
stl.write(' endfacet\n');
stl.write('endsolid surface\n');
stl.close();
return
# Functions for Curve Followers
# Translation Distances - evaluates the parametric function at t = s to determine the translation distance.
def trans_point(r,s): return f[r](t=s)
# Rotation Angles - takes the tangent line at time t = v and projects it onto the xy, yz, and xz planes when u = 0, 1, or 2 respectively. Then it finds the angle between the projection and the respective x, y and z unit vectors using arctangent of the slope.
def angle(u,v):
if diff(f[u])(t = v) == 0:
return pi/2
else:
return arctan(diff(f[(u+1)%3])(t = v)/diff(f[u])(t = v))
# Sphere Curve Follower - f is a parametric function for a curve such that f = [x(t),y(t),z(t)], n is the number of shapes desired on the curve, (a,b) is the interval for t, and radius is the radius of each sphere.
def curve_fol_sphere(f,n,(a,b),radius):
step = (b - a)/n
x0 = trans_point(0,a)
y0 = trans_point(1,a)
z0 = trans_point(2,a)
combo = TransZ(TransY(TransX(Sphere(radius),x0),y0),z0)
for i in range(1,n):
tnew = a + i*step
x1 = trans_point(0,tnew)
y1 = trans_point(1,tnew)
z1 = trans_point(2,tnew)
next_sphere = TransZ(TransY(TransX(Sphere(radius),x1),y1),z1)
combo = Union(combo, next_sphere)
return combo
# Cylinder Curve Follower - f is a parametric function for a curve such that f = [x(t),y(t),z(t)], n is the number of shapes desired on the curve, (a,b) is the interval for t, and """""size is the size of each cube."""""""
def curve_fol_cyyl(f,n,(a,b),(c,d,e)):
step = (b - a)/n
x0 = trans_point(0,a)
y0 = trans_point(1,a)
z0 = trans_point(2,a)
angxy0 = angle(0,a)
angyz0 = angle(1,a)+pi/2
angxz0 = angle(2,a)
cyll = RotateRZ(RotateRY(RotateRX(Cyyl(c,d,e),angxy0),angyz0),angxz0)
combo = TransZ(TransY(TransX(cyll,x0),y0),z0)
for i in range(1,n):
tnew = a + i*step
x1 = trans_point(0,tnew)
y1 = trans_point(1,tnew)
z1 = trans_point(2,tnew)
angxy1 = angle(0,tnew)
angyz1 = angle(1,tnew)+pi/2
angxz1 = angle(2,tnew)
next_cyyl = RotateRZ(RotateRY(RotateRX(Cyyl(c,d,e),angxy1),angyz1),angxz1)
move_next_cyyl = TransZ(TransY(TransX(next_cyyl,x1),y1),z1)
combo = Union(combo, move_next_cyyl)
return combo
# Box Curve Follower - f is a parametric function for a curve such that f = [x(t),y(t),z(t)], n is the number of shapes desired on the curve, (a,b) is the interval for t, and size is the size of each cube.
def curve_fol_box(f,n,(a,b),size):
step = (b - a)/n
x0 = trans_point(0,a)
y0 = trans_point(1,a)
z0 = trans_point(2,a)
angxy0 = angle(0,a)
angyz0 = angle(1,a)
angxz0 = angle(2,a)
box = RotateRZ(RotateRY(RotateRX(Cube(size),angxy0),angyz0),angxz0)
combo = TransZ(TransY(TransX(box,x0),y0),z0)
for i in range(1,n):
tnew = a + i*step
x1 = trans_point(0,tnew)
y1 = trans_point(1,tnew)
z1 = trans_point(2,tnew)
angxy1 = angle(0,tnew)
angyz1 = angle(1,tnew)
angxz1 = angle(2,tnew)
next_box = RotateRZ(RotateRY(RotateRX(Cube(size),angxy1),angyz1),angxz1)
move_next_box = TransZ(TransY(TransX(next_box,x1),y1),z1)
combo = Union(combo,move_next_box)
return combo
///
}}}
{{{id=92|
# "Drope Block" a sphere with the drope function on 4 of its 6 "sides" - Adapted from Nicholas Lewis
sp = Sphere(4)
dr = TransZ(Drope (4),2.5)
top = dr
side1 = RotateY(top,90)
ts1 = Intersect(top ,side1)
bs2 = RotateY(ts1,180)
tube =Intersect(ts1 ,bs2)
sld = Intersect(tube,sp)
dbplt = Plot3d(sld,5,100)
dbplt
///
}}}
{{{id=94|
# "Sinus Bugs" two interlocking bugs patterned with a sinus - Adapted from Nicholas Lewis
t1 = 5*Bug (-.5,1, 4)
t2 = 5*Bug (.5,1, 4)
s = Sinus (0,5)
sld = Combine(Union(t1,t2),s,1)
plt = Plot3d(sld,6,200)
plt
///
}}}
{{{id=93|
# "Super Wiffle" a ball locked inside a tetrahedron
c = Tetrahedral (18,18,60)
sp = Sphere(1.5)
sld = Union(c,sp)
swplt = Plot3d(sld,10,100)
swplt
///
}}}
{{{id=95|
# "Sinus Sphere" a sphere patterned with a sinus - Adapted from Nicholas Lewis
a = Sinus(0,4)
b = (cos(x) * cos(y)) * 10
c = Sphere(4)
sld = Intersect(b-a+2,c)
Splt = Plot3d(sld,4,100)
Splt
///
}}}
{{{id=96|
#Heart
hearts=Plot3d(Heart(1,9/4,1,1,9/80),1.5,100)
hearts
///
}}}
{{{id=86|
# "Spaceship" a spaceship modeled using prisms, a cone, a cylinder and the sinus sphere
q=TransZ(TransX(RotateRY(Prism(1,.25,3),-pi/9),.75),8)
r=Union(q,RotateRZ(q,pi))
r2=Union(r,RotateRZ(r,pi/2))
spaceship=implicit_plot3d(Union(r2,Union(TransZ(RotateRX(ScaleXYZ(sld,.25),pi),10),Union(TransZ(Cyyl(.35,.35,6),7),Cone(1,1,9)))),(x, -5, 5), (y, -5, 5), (z, -.25, 12), plot_points=100)
spaceship
///
}}}
{{{id=103|
# "Sphere Circle"
t = var('t')
f = [4*sin(t), 4*cos(t), 0*t]
circ=curve_fol_sphere(f,7,(0,2*pi),2)
Plot3d(circ,10,100)
///
}}}
{{{id=105|
# "Stacked Sphere Circle"
t = var('t')
f = [4*sin(t), 4*cos(t), 2*t/t]
circ2=curve_fol_sphere(f,7,(0,2*pi),2)
criclplt=Plot3d(Union(circ,circ2),10,100)
criclplt
///
}}}
{{{id=110|
#stl_outputs
surface_to_stl_filesave(spaceship,'spaceship')
surface_to_stl_filesave(criclplt,'criclplt')
surface_to_stl_filesave(hearts,'hearts')
surface_to_stl_filesave(Splt,'Splt')
surface_to_stl_filesave(swplt,'swplt')
surface_to_stl_filesave(plt,'plt')
surface_to_stl_filesave(dbplt,'dbplt')
///
}}}
{{{id=111|
%auto
# Interactive Curve Follower: Various parameters controlling the creation of a parametric curve on which 3D object will be placed are presented. Once the proper parameters have been selected a 3d model of the figure is displayed.
var('t x y z')
inter_plot = 0
@interact
def curve_fol(shape_type = ("Shape",['Sphere','Cube','Cylinder']), rotshapes = checkbox(True, "Rotate Shapes to Match Curve?"), n = ("# of Shapes on Curve",1), xfun = ("x(t)", t), yfun = ("y(t)", t), zfun = ("z(t)", t), a = ("t start", 0), b = ("t end",1), radius = ("radius", 1), height = ("height", 1), showplot = checkbox(True, "Display Shape?"), xaxes = ("X-axis Size", 10), yaxes = ("Y-axis Size", 10), zaxes = ("Z-axis Size", 10), saveplot = checkbox(False, "Save Plot?"), auto_update = False):
var('t x y z')
g = (xfun,yfun,zfun)
if shape_type == 'Cylinder':
shape_fun = Cyyl
if shape_type == 'Sphere':
shape_fun = Sphere
if shape_type == 'Cube':
shape_fun = Cube
def trans_point2(r,s): return g[r](t=s)
def angle(u,v):
if diff(g[u])(t = v) == 0:
return pi/2
else:
return arctan(diff(g[(u+1)%3])(t = v)/diff(g[u])(t = v))
if shape_type == 'Cylinder':
step = (b - a)/n
x0 = trans_point2(0,a)
y0 = trans_point2(1,a)
z0 = trans_point2(2,a)
angxy0 = angle(0,a)
angyz0 = angle(1,a)+pi/2
angxz0 = angle(2,a)
if rotshapes:
shape1 = RotateRZ(RotateRY(RotateRX(Cyyl(radius,radius,height),angxy0),angyz0),angxz0)
else:
shape1 = Cyyl(radius,radius,height)
combo = TransZ(TransY(TransX(shape1,x0),y0),z0)
for i in range(1,n):
tnew = a + i*step
x1 = trans_point2(0,tnew)
y1 = trans_point2(1,tnew)
z1 = trans_point2(2,tnew)
angxy1 = angle(0,tnew)
angyz1 = angle(1,tnew)+pi/2
angxz1 = angle(2,tnew)
if rotshapes:
shape_next = RotateRZ(RotateRY(RotateRX(Cyyl(radius,radius,height),angxy0),angyz0),angxz0)
else:
shape_next = Cyyl(radius,radius,height)
next_sphere = TransZ(TransY(TransX(shape_next,x1),y1),z1)
combo = Union(combo, next_sphere)
elif shape_type == 'Cube':
step = (b - a)/n
x0 = trans_point2(0,a)
y0 = trans_point2(1,a)
z0 = trans_point2(2,a)
angxy0 = angle(0,a)
angyz0 = angle(1,a)
angxz0 = angle(2,a)
if rotshapes:
shape1 = RotateRZ(RotateRY(RotateRX(Cube(radius),angxy0),angyz0),angxz0)
else:
shape1 = Cube(radius)
combo = TransZ(TransY(TransX(shape1,x0),y0),z0)
for i in range(1,n):
tnew = a + i*step
x1 = trans_point2(0,tnew)
y1 = trans_point2(1,tnew)
z1 = trans_point2(2,tnew)
angxy1 = angle(0,tnew)
angyz1 = angle(1,tnew)
angxz1 = angle(2,tnew)
if rotshapes:
shape_next = RotateRZ(RotateRY(RotateRX(Cube(radius),angxy0),angyz0),angxz0)
else:
shape_next = Cube(radius)
next_sphere = TransZ(TransY(TransX(shape_next,x1),y1),z1)
combo = Union(combo, next_sphere)
else:
step = (b - a)/n
x0 = trans_point2(0,a)
y0 = trans_point2(1,a)
z0 = trans_point2(2,a)
combo = TransZ(TransY(TransX(Sphere(radius),x0),y0),z0)
for i in range(1,n):
tnew = a + i*step
x1 = trans_point2(0,tnew)
y1 = trans_point2(1,tnew)
z1 = trans_point2(2,tnew)
next_sphere = TransZ(TransY(TransX(Sphere(radius),x1),y1),z1)
combo = Union(combo, next_sphere)
inter_plot = implicit_plot3d(combo,(x,-xaxes,xaxes),(y,-yaxes,yaxes),(z,-zaxes,zaxes))
if showplot:
show(inter_plot)
if saveplot:
surface_to_stl_filesave(inter_plot,'inter_plot')
return inter_plot
///
}}}
{{{id=114|
var('t')
f = [t,t*t,t]
v = curve_fol_cyyl(f,3,(0,1),(1,1,1))
Plot3d(v,6,100)
///
}}}