Skip to content

Instantly share code, notes, and snippets.

@wware
Last active August 29, 2015 14:21
Show Gist options
  • Save wware/69b36adbdee2f4a72f43 to your computer and use it in GitHub Desktop.
Save wware/69b36adbdee2f4a72f43 to your computer and use it in GitHub Desktop.
Do some linear algebra and generate a 3D shape. Kind of a skewed truncated pyramid.
import doctest
import numpy
import sys
def is_vector(x):
return isinstance(x, numpy.ndarray) and x.shape == (3,)
def normalize(vec):
length = numpy.linalg.norm(vec)
assert length > 1.e-10
return (1. / length) * vec
def rotate(vec, axis, theta):
axis = normalize(axis)
vpar = numpy.dot(vec, axis) * axis
u = vec - vpar
um = numpy.linalg.norm(u)
uhat = (1. / um) * u
vhat = numpy.cross(axis, uhat)
return vpar + (um * numpy.cos(theta)) * uhat + (um * numpy.sin(theta)) * vhat
class Line:
"""
We can define a line in terms of its direction vector and any point known
to lie on the line. It could also be two points.
"""
def __init__(self, member, direction=None, second=None):
if second is not None and is_vector(second):
assert direction is None
direction = second - member
elif direction is not None and is_vector(direction):
assert second is None
else:
raise Exception
self.direction = direction
self.member = member
def __repr__(self):
return "<Line {0} {1}>".format(self.member, self.direction)
def intersect(self, other):
if isinstance(other, Plane):
return other.intersect(self)
# intersect two lines - raise exception if they are not coplanar
# or if they are parallel and not identical, if identical just return
# either of them
raise Exception("not implemented yet")
class Plane:
"""
We can define a plane in terms of its normal vector and any point known
to lie on the plane. It could also be a line and a point.
"""
def __init__(self, member, normal=None, line=None):
"""
"""
if line is not None and isinstance(line, Line):
assert normal is None
normal = numpy.cross(member - line.member, line.direction)
elif normal is not None and is_vector(normal):
assert line is None
else:
raise Exception
self.normal = normalize(normal)
self.member = member
def intersect(self, other):
if isinstance(other, Line):
x0, n = self.member, self.normal
x1 = other.member
x2 = x1 + other.direction
m = numpy.dot(x0 - x1, n) / numpy.dot(x2 - x1, n)
return m * (x2 - x1) + x1
elif isinstance(other, Plane):
x1, x2 = self.member, other.member
n1, n2 = self.normal, other.normal
direction = normalize(numpy.cross(n1, n2))
n1n2n1 = numpy.cross(direction, n1)
m = numpy.dot(x2 - x1, n2) / numpy.dot(n1n2n1, n2)
x = x1 + m * n1n2n1
return Line(x, direction=direction)
else:
raise Exception
class StlFile:
def __init__(self, vertices, paths):
self.vertices = vertices
self.paths = paths
def dump(self):
print "solid Foo"
for p in self.paths:
points = []
verts = []
for index in p:
point = self.vertices[index]
points.append(point)
verts.append((25.4 * point[0], 25.4 * point[1], 25.4 * point[2]))
n = normalize(numpy.cross(points[1] - points[0], points[2] - points[1]))
v0 = verts[0]
for i in range(1, len(verts) - 1):
# Assume each surface is simple convex for now.
print " facet normal {0:e} {1:e} {2:e}".format(n[0], n[1], n[2])
print " outer loop"
print " vertex {0:e} {1:e} {2:e}".format(v0[0], v0[1], v0[2])
print " vertex {0:e} {1:e} {2:e}".format(verts[i][0], verts[i][1], verts[i][2])
print " vertex {0:e} {1:e} {2:e}".format(verts[i+1][0], verts[i+1][1], verts[i+1][2])
print " endloop"
print " endfacet"
print "endsolid Foo"
# Dimensions of front surface
A = 22
B = 14
C = 5
dx = B / numpy.tan(numpy.pi / 3)
# The front surface
w = numpy.array([0, 0, 0])
x = numpy.array([dx, B, 0])
y = numpy.array([A + dx, B, 0])
z = numpy.array([A, 0, 0])
front_normal = normalize(numpy.cross(w - x, y - x))
theta = (2./3) * numpy.pi
xy_side_normal = rotate(front_normal, y - x, theta)
yz_side_normal = rotate(front_normal, z - y, theta)
zw_side_normal = rotate(front_normal, w - z, theta)
wx_side_normal = rotate(front_normal, x - w, theta)
xy_side_plane = Plane(x, normal=xy_side_normal)
yz_side_plane = Plane(y, normal=yz_side_normal)
zw_side_plane = Plane(z, normal=zw_side_normal)
wx_side_plane = Plane(w, normal=wx_side_normal)
bottom = Plane(numpy.array([0, 0, -C]),
normal=numpy.array([0, 0, 1]))
s = zw_side_plane.intersect(wx_side_plane).intersect(bottom)
t = wx_side_plane.intersect(xy_side_plane).intersect(bottom)
u = xy_side_plane.intersect(yz_side_plane).intersect(bottom)
v = yz_side_plane.intersect(zw_side_plane).intersect(bottom)
stl = StlFile(
(w, x, y, z, s, t, u, v),
(
(0, 2, 1, 3),
(4, 5, 6, 7),
(0, 4, 5, 1),
(1, 5, 6, 2),
(2, 6, 7, 3),
(3, 7, 4, 0)
)
)
stl.dump()
Display the source blob
Display the rendered blob
Raw
solid Foo
facet normal 0.000000e+00 -0.000000e+00 1.000000e+00
outer loop
vertex 0.000000e+00 0.000000e+00 0.000000e+00
vertex 7.641058e+02 3.556000e+02 0.000000e+00
vertex 2.053058e+02 3.556000e+02 0.000000e+00
endloop
endfacet
facet normal 0.000000e+00 -0.000000e+00 1.000000e+00
outer loop
vertex 0.000000e+00 0.000000e+00 0.000000e+00
vertex 2.053058e+02 3.556000e+02 0.000000e+00
vertex 5.588000e+02 0.000000e+00 0.000000e+00
endloop
endfacet
facet normal 0.000000e+00 0.000000e+00 -1.000000e+00
outer loop
vertex -1.270000e+02 -7.332348e+01 -1.270000e+02
vertex 1.629724e+02 4.289235e+02 -1.270000e+02
vertex 8.911058e+02 4.289235e+02 -1.270000e+02
endloop
endfacet
facet normal 0.000000e+00 0.000000e+00 -1.000000e+00
outer loop
vertex -1.270000e+02 -7.332348e+01 -1.270000e+02
vertex 8.911058e+02 4.289235e+02 -1.270000e+02
vertex 6.011333e+02 -7.332348e+01 -1.270000e+02
endloop
endfacet
facet normal 7.500000e-01 -4.330127e-01 -5.000000e-01
outer loop
vertex 0.000000e+00 0.000000e+00 0.000000e+00
vertex -1.270000e+02 -7.332348e+01 -1.270000e+02
vertex 1.629724e+02 4.289235e+02 -1.270000e+02
endloop
endfacet
facet normal 7.500000e-01 -4.330127e-01 -5.000000e-01
outer loop
vertex 0.000000e+00 0.000000e+00 0.000000e+00
vertex 1.629724e+02 4.289235e+02 -1.270000e+02
vertex 2.053058e+02 3.556000e+02 0.000000e+00
endloop
endfacet
facet normal -2.146563e-16 -8.660254e-01 -5.000000e-01
outer loop
vertex 2.053058e+02 3.556000e+02 0.000000e+00
vertex 1.629724e+02 4.289235e+02 -1.270000e+02
vertex 8.911058e+02 4.289235e+02 -1.270000e+02
endloop
endfacet
facet normal -2.146563e-16 -8.660254e-01 -5.000000e-01
outer loop
vertex 2.053058e+02 3.556000e+02 0.000000e+00
vertex 8.911058e+02 4.289235e+02 -1.270000e+02
vertex 7.641058e+02 3.556000e+02 0.000000e+00
endloop
endfacet
facet normal -7.500000e-01 4.330127e-01 -5.000000e-01
outer loop
vertex 7.641058e+02 3.556000e+02 0.000000e+00
vertex 8.911058e+02 4.289235e+02 -1.270000e+02
vertex 6.011333e+02 -7.332348e+01 -1.270000e+02
endloop
endfacet
facet normal -7.500000e-01 4.330127e-01 -5.000000e-01
outer loop
vertex 7.641058e+02 3.556000e+02 0.000000e+00
vertex 6.011333e+02 -7.332348e+01 -1.270000e+02
vertex 5.588000e+02 0.000000e+00 0.000000e+00
endloop
endfacet
facet normal 2.951524e-16 8.660254e-01 -5.000000e-01
outer loop
vertex 5.588000e+02 0.000000e+00 0.000000e+00
vertex 6.011333e+02 -7.332348e+01 -1.270000e+02
vertex -1.270000e+02 -7.332348e+01 -1.270000e+02
endloop
endfacet
facet normal 2.951524e-16 8.660254e-01 -5.000000e-01
outer loop
vertex 5.588000e+02 0.000000e+00 0.000000e+00
vertex -1.270000e+02 -7.332348e+01 -1.270000e+02
vertex 0.000000e+00 0.000000e+00 0.000000e+00
endloop
endfacet
endsolid Foo
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment