diff --git a/.gitignore b/.gitignore
index e43b0f9..263bf47 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1 +1,4 @@
 .DS_Store
+__pycache__
+*.pyc
+.*.swp
diff --git a/csg/core.py b/csg/core.py
index 355d671..d4ec318 100644
--- a/csg/core.py
+++ b/csg/core.py
@@ -1,7 +1,11 @@
 import math
 import operator
-from csg.geom import *
 from functools import reduce
+from csg.geom import *
+from csg.stl import (
+    save_polys_to_stl_file,
+    read_polys_from_stl_file,
+)
 
 class CSG(object):
     """
@@ -215,6 +219,25 @@ def saveVTK(self, filename):
                     f.write('{0} '.format(index))
                 f.write('\n')
 
+    def saveSTL(self, filename, binary=True):
+        """
+        Save polygons in STL file.
+        binary - if true (default), file is written in binary STL format.  Otherwise ASCII STL format.
+        """
+        save_polys_to_stl_file(
+            self.toPolygons(),
+            filename,
+            binary
+        )
+
+    @classmethod
+    def readSTL(cls, filename):
+        """
+        Reads polygons from an STL file and generates a CSG object from them.
+        """
+        polygons = read_polys_from_stl_file(filename)
+        return CSG.fromPolygons(polygons)
+
     def union(self, csg):
         """
         Return a new CSG solid representing space in either this solid or in the
diff --git a/csg/stl.py b/csg/stl.py
new file mode 100644
index 0000000..8d6f26d
--- /dev/null
+++ b/csg/stl.py
@@ -0,0 +1,182 @@
+import struct
+from csg.geom import Polygon, Vertex
+
+
+class StlEndOfFileException(Exception):
+    """Exception class for reaching the end of the STL file while reading."""
+    pass
+
+
+class StlMalformedLineException(Exception):
+    """Exception class for malformed lines in the STL file being read."""
+    pass
+
+
+def _float_fmt(val):
+    """
+    Make a short, clean, string representation of a float value.
+    """
+    s = ("%.6f" % val).rstrip('0').rstrip('.')
+    return '0' if s == '-0' else s
+
+
+def _stl_write_facet(poly, f, binary=True):
+    """
+    Writes a single triangle facet to the given STL file stream.
+    binary - Save in binary format if True, else ASCII format.
+    """
+    norm = poly.plane.normal
+    v0, v1, v2 = poly.vertices
+    if binary:
+        data = struct.pack(
+            '<3f 3f 3f 3f H',
+            norm[0], norm[1], norm[2],
+            v0.pos[0], v0.pos[1], v0.pos[2],
+            v1.pos[0], v1.pos[1], v1.pos[2],
+            v2.pos[0], v2.pos[1], v2.pos[2],
+            0
+        )
+        f.write(data)
+    else:
+        v0 = " ".join(_float_fmt(x) for x in v0.pos)
+        v1 = " ".join(_float_fmt(x) for x in v1.pos)
+        v2 = " ".join(_float_fmt(x) for x in v2.pos)
+        norm = " ".join(_float_fmt(x) for x in norm)
+        vfmt = (
+            "  facet normal {norm}\n"
+            "    outer loop\n"
+            "      vertex {v0}\n"
+            "      vertex {v1}\n"
+            "      vertex {v2}\n"
+            "    endloop\n"
+            "  endfacet\n"
+        )
+        data = vfmt.format(norm=norm, v0=v0, v1=v1, v2=v2)
+        f.write(bytes(data, encoding='ascii'))
+
+
+def save_polys_to_stl_file(polys, filename, binary=True):
+    """
+    Save polygons in STL file.
+    polys - list of Polygons.
+    filename - Name fo the STL file to save to.
+    binary - if true (default), file is written in binary STL format.  Otherwise ASCII STL format.
+    """
+    # Convert all polygons to triangles.
+    tris = []
+    for poly in polys:
+        vlen = len(poly.vertices)
+        for n in range(1,vlen-1):
+            tris.append(
+                Polygon([
+                    poly.vertices[0],
+                    poly.vertices[n%vlen],
+                    poly.vertices[(n+1)%vlen],
+                ])
+            )
+    if binary:
+        with open(filename, 'wb') as f:
+            f.write(b'%-80s' % b'Binary STL Model')
+            f.write(struct.pack('<I', len(tris)))
+            for tri in tris:
+                _stl_write_facet(tri, f, binary=binary)
+    else:
+        with open(filename, 'wb') as f:
+            f.write(b"solid Model\n")
+            for tri in tris:
+                _stl_write_facet(tri, f, binary=binary)
+            f.write(b"endsolid Model\n")
+
+
+def _read_ascii_line(f, watchwords=None):
+    """
+    Reads a single line from an ASCII STL file stream and checks for required keywords.
+    Returns array of float values from the read line.
+    Throws StlEndOfFileException if 'endsolid' line is read.
+    Throws StlMalformedLineException if keywords are not found.
+    """
+    line = f.readline(1024).decode("ascii")
+    if line == "":
+        raise StlEndOfFileException()
+    words = line.strip(' \t\n\r').lower().split()
+    if words[0] == 'endsolid':
+        raise StlEndOfFileException()
+    argstart = 0
+    if watchwords:
+        watchwords = watchwords.lower().split()
+        argstart = len(watchwords)
+        for i in range(argstart):
+            if words[i] != watchwords[i]:
+                raise StlMalformedLineException()
+    return [float(val) for val in words[argstart:]]
+
+
+def _read_ascii_facet(f):
+    """
+    Load a single facet triangle from the ASCII STL file stream.
+    Skips corrupted facets if it can.
+    Returns a Polygon.
+    Throws StlEndOfFileException if EOF is reached.
+    """
+    while True:
+        try:
+            normal = _read_ascii_line(f, watchwords='facet normal')
+            _read_ascii_line(f, watchwords='outer loop')
+            v0 = _read_ascii_line(f, watchwords='vertex')
+            v1 = _read_ascii_line(f, watchwords='vertex')
+            v2 = _read_ascii_line(f, watchwords='vertex')
+            _read_ascii_line(f, watchwords='endloop')
+            _read_ascii_line(f, watchwords='endfacet')
+            if v0 == v1 or v1 == v2 or v2 == v0:
+                continue  # zero area facet.  Skip to next facet.
+        except StlEndOfFileException:
+            return None
+        except StlMalformedLineException:
+            continue  # Skip to next facet.
+        v0 = Vertex(v0)
+        v1 = Vertex(v1)
+        v2 = Vertex(v2)
+        return Polygon([v0, v1, v2])
+
+
+def _read_binary_facet(f):
+    """
+    Load a single facet triangle from the binary STL file stream.
+    Returns a Polygon.
+    """
+    data = struct.unpack('<3f 3f 3f 3f H', f.read(4*4*3+2))
+    normal = data[0:3]
+    v0 = Vertex(data[3:6])
+    v1 = Vertex(data[6:9])
+    v2 = Vertex(data[9:12])
+    return Polygon([v0, v1, v2])
+
+
+def read_polys_from_stl_file(filename):
+    """
+    Read array of triangle polygons from an STL file.
+    filename - Name fo the STL file to read from.
+    """
+    polygons = []
+    with open(filename, 'rb') as f:
+        line = f.readline(80)
+        if line == "":
+            return  # End of file.
+        line = line.lstrip()
+        if line[0:6].lower() == b"solid ":
+            # Reading ASCII STL file.
+            while True:
+                poly = _read_ascii_facet(f)
+                if not poly:
+                    break
+                polygons.append(poly)
+        else:
+            # Reading Binary STL file.
+            chunk = f.read(4)
+            facets = struct.unpack('<I', chunk)[0]
+            for n in range(facets):
+                poly = _read_binary_facet(f)
+                if not poly:
+                    break
+                polygons.append(poly)
+    return polygons
diff --git a/examples/pyopengl.py b/examples/pyopengl.py
index 28b65f4..5c0a72f 100644
--- a/examples/pyopengl.py
+++ b/examples/pyopengl.py
@@ -1,5 +1,6 @@
 import sys
 import os
+from collections import namedtuple
 
 from OpenGL.GL import *
 from OpenGL.GLUT import *
@@ -12,28 +13,49 @@
 
 from optparse import OptionParser
 
-light_ambient = [0.3, 0.3, 0.3, 1.0]
-light_diffuse = [0.7, 0.7, 0.7, 1.0]  # Red diffuse light
-light_position = [100.0, 100.0, 100.0, 0.0]  # Infinite light location.
+Light = namedtuple('Light', ['num', 'ambient', 'diffuse', 'position'])
+
+lights = [
+    Light(
+        num=GL_LIGHT0,
+        ambient=[0.3, 0.3, 0.3, 1.0],
+        diffuse=[0.4, 0.4, 0.4, 1.0],
+        position=[100.0, 100.0, 100.0, 0.0], # up and right, towards viewer
+    ),
+    Light(
+        num=GL_LIGHT1,
+        ambient=[0.3, 0.3, 0.3, 1.0],
+        diffuse=[0.4, 0.4, 0.4, 1.0],
+        position=[-100.0, 100.0, 100.0, 0.0], # up and left, towards viewer
+    ),
+    Light(
+        num=GL_LIGHT2,
+        ambient=[0.3, 0.3, 0.3, 1.0],
+        diffuse=[0.4, 0.4, 0.4, 1.0],
+        position=[0.0, -100.0, 100.0, 0.0], # down and center, towards user
+    ),
+]
+
 
 rot = 0.0
 
 class TestRenderable(object):
-    def __init__(self, operation):
+    def __init__(self, operation, wireframe):
         self.faces = []
         self.normals = []
         self.vertices = []
         self.colors = []
         self.vnormals = []
         self.list = -1
-        
+        self.wireframe = wireframe
+
         a = CSG.cube()
         b = CSG.cylinder(radius=0.5, start=[0., -2., 0.], end=[0., 2., 0.])
         for p in a.polygons:
             p.shared = [1.0, 0.0, 0.0, 1.0]
         for p in b.polygons:
             p.shared = [0.0, 1.0, 0.0, 1.0]
-            
+
         recursionlimit = sys.getrecursionlimit()
         sys.setrecursionlimit(10000)
         try:
@@ -48,7 +70,7 @@ def __init__(self, operation):
         except RuntimeError as e:
             raise RuntimeError(e)
         sys.setrecursionlimit(recursionlimit)
-        
+
         for polygon in polygons:
             n = polygon.plane.normal
             indices = []
@@ -63,7 +85,7 @@ def __init__(self, operation):
             self.faces.append(indices)
             self.normals.append([n.x, n.y, n.z])
             self.colors.append(polygon.shared)
-        
+
         # setup vertex-normals
         ns = []
         for vns in self.vnormals:
@@ -73,38 +95,40 @@ def __init__(self, operation):
             n = n.dividedBy(len(vns))
             ns.append([a for a in n])
         self.vnormals = ns
-        
+
     def render(self):
+        if self.wireframe:
+            glPolygonMode(GL_FRONT_AND_BACK, GL_LINE)
         if self.list < 0:
             self.list = glGenLists(1)
             glNewList(self.list, GL_COMPILE)
-            
+
             for n, f in enumerate(self.faces):
-                glMaterialfv(GL_FRONT, GL_DIFFUSE, self.colors[n])
-                glMaterialfv(GL_FRONT, GL_SPECULAR, self.colors[n])
+                colors = self.colors[n]
+                glMaterialfv(GL_FRONT, GL_DIFFUSE, colors)
+                glMaterialfv(GL_FRONT, GL_SPECULAR, colors)
                 glMaterialf(GL_FRONT, GL_SHININESS, 50.0)
-                glColor4fv(self.colors[n])
-            
-                glBegin(GL_POLYGON)
-                if self.colors[n][0] > 0:
-                    glNormal3fv(self.normals[n])
+                glColor4fv(colors)
 
+                glBegin(GL_POLYGON)
+                glNormal3fv(self.normals[n])
                 for i in f:
-                    if self.colors[n][1] > 0:
-                        glNormal3fv(self.vnormals[i])
+                    # Disabled vertex normals to make faces clearer.
+                    # if sum(x*x for x in self.vnormals[i]) > 1e-4:
+                    #     glNormal3fv(self.vnormals[i])
                     glVertex3fv(self.vertices[i])
                 glEnd()
             glEndList()
         glCallList(self.list)
-        
+
 renderable = None
 
 def init():
-    # Enable a single OpenGL light.
-    glLightfv(GL_LIGHT0, GL_AMBIENT, light_ambient)
-    glLightfv(GL_LIGHT0, GL_DIFFUSE, light_diffuse)
-    glLightfv(GL_LIGHT0, GL_POSITION, light_position)
-    glEnable(GL_LIGHT0);
+    for light in lights:
+        glLightfv(light.num, GL_AMBIENT, light.ambient)
+        glLightfv(light.num, GL_DIFFUSE, light.diffuse)
+        glLightfv(light.num, GL_POSITION, light.position)
+        glEnable(light.num);
     glEnable(GL_LIGHTING);
 
     # Use depth buffering for hidden surface elimination.
@@ -115,38 +139,39 @@ def init():
     gluPerspective(40.0, 640./480., 1.0, 10.0);
     glMatrixMode(GL_MODELVIEW);
     gluLookAt(0.0, 0.0, 5.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.)
-    
+
 def display():
     global rot
     glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
-    
+
     glPushMatrix()
     glTranslatef(0.0, 0.0, -1.0);
-    glRotatef(rot, 1.0, 0.0, 0.0);
-    glRotatef(rot, 0.0, 0.0, 1.0);
+    glRotatef(rot*7, 1.0, 0.0, 0.0);
+    glRotatef(rot*13, 0.0, 0.0, 1.0);
     rot += 0.1
-    
+
     renderable.render()
-    
+
     glPopMatrix()
     glFlush()
     glutSwapBuffers()
     glutPostRedisplay()
-    
+
 if __name__ == '__main__':
     parser = OptionParser()
     parser.add_option('-o', '--operation', dest='operation',
         type='str', default='subtract')
+    parser.add_option('-w', '--wireframe', action='store_true')
     (options, args) = parser.parse_args()
 
-    renderable = TestRenderable(options.operation)
-    
+    renderable = TestRenderable(options.operation, options.wireframe)
+
     glutInit()
     glutInitWindowSize(640,480)
     glutCreateWindow("CSG Test")
     glutInitDisplayMode(GLUT_DEPTH | GLUT_DOUBLE | GLUT_RGBA)
     glutDisplayFunc(display)
-     
+
     init()
 
     glutMainLoop()
diff --git a/examples/stlglshow.py b/examples/stlglshow.py
new file mode 100644
index 0000000..147db25
--- /dev/null
+++ b/examples/stlglshow.py
@@ -0,0 +1,165 @@
+import sys
+import os
+from collections import namedtuple
+
+from OpenGL.GL import *
+from OpenGL.GLUT import *
+from OpenGL.GLU import *
+
+sys.path.insert(0, os.getcwd())
+
+from csg.core import CSG
+from csg.geom import Vertex, Vector
+
+from optparse import OptionParser
+
+Light = namedtuple('Light', ['num', 'ambient', 'diffuse', 'position'])
+
+lights = [
+    Light(
+        num=GL_LIGHT0,
+        ambient=[0.3, 0.3, 0.3, 1.0],
+        diffuse=[0.4, 0.4, 0.4, 1.0],
+        position=[1000.0, 1000.0, 1000.0, 0.0], # up and right, towards viewer
+    ),
+    Light(
+        num=GL_LIGHT1,
+        ambient=[0.3, 0.3, 0.3, 1.0],
+        diffuse=[0.4, 0.4, 0.4, 1.0],
+        position=[-1000.0, 1000.0, 1000.0, 0.0], # up and left, towards viewer
+    ),
+    Light(
+        num=GL_LIGHT2,
+        ambient=[0.3, 0.3, 0.3, 1.0],
+        diffuse=[0.4, 0.4, 0.4, 1.0],
+        position=[0.0, -1000.0, 1000.0, 0.0], # down and center, towards user
+    ),
+]
+
+
+rot = 0.0
+
+class TestRenderable(object):
+    def __init__(self, filename, wireframe):
+        self.faces = []
+        self.normals = []
+        self.vertices = []
+        self.colors = []
+        self.vnormals = []
+        self.list = -1
+        self.wireframe = wireframe
+
+        recursionlimit = sys.getrecursionlimit()
+        sys.setrecursionlimit(10000)
+        try:
+            obj = CSG.readSTL(filename)
+            polygons = obj.toPolygons()
+        except RuntimeError as e:
+            raise RuntimeError(e)
+        sys.setrecursionlimit(recursionlimit)
+
+        for p in polygons:
+            p.shared = [0.0, 1.0, 0.0, 1.0]
+
+        for polygon in polygons:
+            n = polygon.plane.normal
+            indices = []
+            for v in polygon.vertices:
+                pos = [v.pos.x, v.pos.y, v.pos.z]
+                if not pos in self.vertices:
+                    self.vertices.append(pos)
+                    self.vnormals.append([])
+                index = self.vertices.index(pos)
+                indices.append(index)
+                self.vnormals[index].append(v.normal)
+            self.faces.append(indices)
+            self.normals.append([n.x, n.y, n.z])
+            self.colors.append(polygon.shared)
+
+        # setup vertex-normals
+        ns = []
+        for vns in self.vnormals:
+            n = Vector(0.0, 0.0, 0.0)
+            for vn in vns:
+                n = n.plus(vn)
+            n = n.dividedBy(len(vns))
+            ns.append([a for a in n])
+        self.vnormals = ns
+
+    def render(self):
+        if self.wireframe:
+            glPolygonMode(GL_FRONT_AND_BACK, GL_LINE)
+        if self.list < 0:
+            self.list = glGenLists(1)
+            glNewList(self.list, GL_COMPILE)
+
+            for n, f in enumerate(self.faces):
+                colors = self.colors[n]
+                glMaterialfv(GL_FRONT, GL_DIFFUSE, colors)
+                glMaterialfv(GL_FRONT, GL_SPECULAR, colors)
+                glMaterialf(GL_FRONT, GL_SHININESS, 50.0)
+                glColor4fv(colors)
+
+                glBegin(GL_POLYGON)
+                glNormal3fv(self.normals[n])
+                for i in f:
+                    # Disabled vertex normals to make faces clearer.
+                    # if sum(x*x for x in self.vnormals[i]) > 1e-4:
+                    #     glNormal3fv(self.vnormals[i])
+                    glVertex3fv(self.vertices[i])
+                glEnd()
+            glEndList()
+        glCallList(self.list)
+
+renderable = None
+
+def init():
+    for light in lights:
+        glLightfv(light.num, GL_AMBIENT, light.ambient)
+        glLightfv(light.num, GL_DIFFUSE, light.diffuse)
+        glLightfv(light.num, GL_POSITION, light.position)
+        glEnable(light.num);
+    glEnable(GL_LIGHTING);
+
+    # Use depth buffering for hidden surface elimination.
+    glEnable(GL_DEPTH_TEST);
+
+    # Setup the view of the cube.
+    glMatrixMode(GL_PROJECTION);
+    gluPerspective(40.0, 640./480., 1.0, 1000.0);
+    glMatrixMode(GL_MODELVIEW);
+    gluLookAt(0.0, 0.0, 200.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.)
+
+def display():
+    global rot
+    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
+
+    glPushMatrix()
+    glTranslatef(0.0, 0.0, -1.0);
+    glRotatef(rot*7, 1.0, 0.0, 0.0);
+    glRotatef(rot*13, 0.0, 0.0, 1.0);
+    rot += 0.1
+
+    renderable.render()
+
+    glPopMatrix()
+    glFlush()
+    glutSwapBuffers()
+    glutPostRedisplay()
+
+if __name__ == '__main__':
+    parser = OptionParser()
+    parser.add_option('-w', '--wireframe', action='store_true')
+    (options, args) = parser.parse_args()
+
+    renderable = TestRenderable(args[0], options.wireframe)
+
+    glutInit()
+    glutInitWindowSize(640,480)
+    glutCreateWindow("STL Show")
+    glutInitDisplayMode(GLUT_DEPTH | GLUT_DOUBLE | GLUT_RGBA)
+    glutDisplayFunc(display)
+
+    init()
+
+    glutMainLoop()
diff --git a/examples/stlsave.py b/examples/stlsave.py
new file mode 100644
index 0000000..5f474ae
--- /dev/null
+++ b/examples/stlsave.py
@@ -0,0 +1,42 @@
+import sys
+import os
+
+from optparse import OptionParser
+
+sys.path.insert(0, os.getcwd())
+
+from csg.core import CSG
+
+
+def render_to_stl_files(operation):
+    a = CSG.cube()
+    b = CSG.cylinder(radius=0.5, start=[0., -2., 0.], end=[0., 2., 0.])
+
+    recursionlimit = sys.getrecursionlimit()
+    sys.setrecursionlimit(10000)
+    try:
+        if not operation:
+            raise Exception('Unknown operation: \'%s\'' % operation)
+        elif 'subtract'.startswith(operation):
+            result = a.subtract(b)
+        elif 'union'.startswith(operation):
+            result = a.union(b)
+        elif 'intersection'.startswith(operation):
+            result = a.intersect(b)
+        else:
+            raise Exception('Unknown operation: \'%s\'' % operation)
+    except RuntimeError as e:
+        raise RuntimeError(e)
+    sys.setrecursionlimit(recursionlimit)
+
+    result.saveSTL('ascii.stl', binary=False)
+    result.saveSTL('binary.stl', binary=True)
+
+
+if __name__ == '__main__':
+    parser = OptionParser()
+    parser.add_option('-o', '--operation', dest='operation',
+        type='str', default='subtract')
+    (options, args) = parser.parse_args()
+    render_to_stl_files(options.operation)
+