Add teapots WIP.teapot

8 files changed, 711 insertions, 8 deletions

diff --git a/dimension/tests/Makefile.am b/dimension/tests/Makefile.am
index 684b2aa..a6c0f15 100644
--- a/dimension/tests/Makefile.am
+++ b/dimension/tests/Makefile.am
@@ -18,7 +18,8 @@
 ###########################################################################
 
 TESTS = demo.dmnsn                                                             \
-        cube.dmnsn
+        cube.dmnsn                                                             \
+        teapot.dmnsn
 TEST_EXTENSIONS = .dmnsn
 DMNSN_LOG_COMPILER = $(top_srcdir)/dimension/dimension
 AM_DMNSN_LOG_FLAGS = --strict -v
diff --git a/dimension/tests/cube.dmnsn b/dimension/tests/cube.dmnsn
index 5f74219..65387ab 100644
--- a/dimension/tests/cube.dmnsn
+++ b/dimension/tests/cube.dmnsn
@@ -42,20 +42,19 @@ objects.append(
   )
 )
 
-class ShinySphere(Sphere):
+class ShinySphere(Teapot):
   def __init__(self, x, y, z, size):
     size -= 1
 
     dx = sin(2*pi*x/size)
     dy = sin(2*pi*y/size)
     dz = sin(2*pi*z/size)
-    Sphere.__init__(
+    center = (5*Vector(x, y, z)/size
+              + Vector(dy + dz, dx + dz, dx + dy)/4
+              - Vector(2.5, 2.5, 2.5))
+    radius = 2/size
+    Teapot.__init__(
       self,
-      center = 5*Vector(x, y, z)/size
-               + Vector(dy + dz, dx + dz, dx + dy)/4
-               - Vector(2.5, 2.5, 2.5),
-      radius = 2/size,
-
       texture = Texture(
         pigment = sRGB(x/size, y/size, z/size),
         finish  = Ambient(sRGB(0.25))
@@ -63,6 +62,7 @@ class ShinySphere(Sphere):
                   + Reflection(0, 0.25)
       )
     )
+    self.scale(radius/2).translate(center)
 
 size = 10
 for x in range(size):
diff --git a/dimension/tests/teapot.dmnsn b/dimension/tests/teapot.dmnsn
new file mode 100644
index 0000000..bcc012a
--- /dev/null
+++ b/dimension/tests/teapot.dmnsn
@@ -0,0 +1,230 @@
+verticies = [
+  Vector( 0.2000,  0.0000, 2.70000), Vector( 0.2000, -0.1120, 2.70000),
+  Vector( 0.1120, -0.2000, 2.70000), Vector( 0.0000, -0.2000, 2.70000),
+  Vector( 1.3375,  0.0000, 2.53125), Vector( 1.3375, -0.7490, 2.53125),
+  Vector( 0.7490, -1.3375, 2.53125), Vector( 0.0000, -1.3375, 2.53125),
+  Vector( 1.4375,  0.0000, 2.53125), Vector( 1.4375, -0.8050, 2.53125),
+  Vector( 0.8050, -1.4375, 2.53125), Vector( 0.0000, -1.4375, 2.53125),
+  Vector( 1.5000,  0.0000, 2.40000), Vector( 1.5000, -0.8400, 2.40000),
+  Vector( 0.8400, -1.5000, 2.40000), Vector( 0.0000, -1.5000, 2.40000),
+  Vector( 1.7500,  0.0000, 1.87500), Vector( 1.7500, -0.9800, 1.87500),
+  Vector( 0.9800, -1.7500, 1.87500), Vector( 0.0000, -1.7500, 1.87500),
+  Vector( 2.0000,  0.0000, 1.35000), Vector( 2.0000, -1.1200, 1.35000),
+  Vector( 1.1200, -2.0000, 1.35000), Vector( 0.0000, -2.0000, 1.35000),
+  Vector( 2.0000,  0.0000, 0.90000), Vector( 2.0000, -1.1200, 0.90000),
+  Vector( 1.1200, -2.0000, 0.90000), Vector( 0.0000, -2.0000, 0.90000),
+  Vector(-2.0000,  0.0000, 0.90000), Vector( 2.0000,  0.0000, 0.45000),
+  Vector( 2.0000, -1.1200, 0.45000), Vector( 1.1200, -2.0000, 0.45000),
+  Vector( 0.0000, -2.0000, 0.45000), Vector( 1.5000,  0.0000, 0.22500),
+  Vector( 1.5000, -0.8400, 0.22500), Vector( 0.8400, -1.5000, 0.22500),
+  Vector( 0.0000, -1.5000, 0.22500), Vector( 1.5000,  0.0000, 0.15000),
+  Vector( 1.5000, -0.8400, 0.15000), Vector( 0.8400, -1.5000, 0.15000),
+  Vector( 0.0000, -1.5000, 0.15000), Vector(-1.6000,  0.0000, 2.02500),
+  Vector(-1.6000, -0.3000, 2.02500), Vector(-1.5000, -0.3000, 2.25000),
+  Vector(-1.5000,  0.0000, 2.25000), Vector(-2.3000,  0.0000, 2.02500),
+  Vector(-2.3000, -0.3000, 2.02500), Vector(-2.5000, -0.3000, 2.25000),
+  Vector(-2.5000,  0.0000, 2.25000), Vector(-2.7000,  0.0000, 2.02500),
+  Vector(-2.7000, -0.3000, 2.02500), Vector(-3.0000, -0.3000, 2.25000),
+  Vector(-3.0000,  0.0000, 2.25000), Vector(-2.7000,  0.0000, 1.80000),
+  Vector(-2.7000, -0.3000, 1.80000), Vector(-3.0000, -0.3000, 1.80000),
+  Vector(-3.0000,  0.0000, 1.80000), Vector(-2.7000,  0.0000, 1.57500),
+  Vector(-2.7000, -0.3000, 1.57500), Vector(-3.0000, -0.3000, 1.35000),
+  Vector(-3.0000,  0.0000, 1.35000), Vector(-2.5000,  0.0000, 1.12500),
+  Vector(-2.5000, -0.3000, 1.12500), Vector(-2.6500, -0.3000, 0.93750),
+  Vector(-2.6500,  0.0000, 0.93750), Vector(-2.0000, -0.3000, 0.90000),
+  Vector(-1.9000, -0.3000, 0.60000), Vector(-1.9000,  0.0000, 0.60000),
+  Vector( 1.7000,  0.0000, 1.42500), Vector( 1.7000, -0.6600, 1.42500),
+  Vector( 1.7000, -0.6600, 0.60000), Vector( 1.7000,  0.0000, 0.60000),
+  Vector( 2.6000,  0.0000, 1.42500), Vector( 2.6000, -0.6600, 1.42500),
+  Vector( 3.1000, -0.6600, 0.82500), Vector( 3.1000,  0.0000, 0.82500),
+  Vector( 2.3000,  0.0000, 2.10000), Vector( 2.3000, -0.2500, 2.10000),
+  Vector( 2.4000, -0.2500, 2.02500), Vector( 2.4000,  0.0000, 2.02500),
+  Vector( 2.7000,  0.0000, 2.40000), Vector( 2.7000, -0.2500, 2.40000),
+  Vector( 3.3000, -0.2500, 2.40000), Vector( 3.3000,  0.0000, 2.40000),
+  Vector( 2.8000,  0.0000, 2.47500), Vector( 2.8000, -0.2500, 2.47500),
+  Vector( 3.5250, -0.2500, 2.49375), Vector( 3.5250,  0.0000, 2.49375),
+  Vector( 2.9000,  0.0000, 2.47500), Vector( 2.9000, -0.1500, 2.47500),
+  Vector( 3.4500, -0.1500, 2.51250), Vector( 3.4500,  0.0000, 2.51250),
+  Vector( 2.8000,  0.0000, 2.40000), Vector( 2.8000, -0.1500, 2.40000),
+  Vector( 3.2000, -0.1500, 2.40000), Vector( 3.2000,  0.0000, 2.40000),
+  Vector( 0.0000,  0.0000, 3.15000), Vector( 0.8000,  0.0000, 3.15000),
+  Vector( 0.8000, -0.4500, 3.15000), Vector( 0.4500, -0.8000, 3.15000),
+  Vector( 0.0000, -0.8000, 3.15000), Vector( 0.0000,  0.0000, 2.85000),
+  Vector( 1.4000,  0.0000, 2.40000), Vector( 1.4000, -0.7840, 2.40000),
+  Vector( 0.7840, -1.4000, 2.40000), Vector( 0.0000, -1.4000, 2.40000),
+  Vector( 0.4000,  0.0000, 2.55000), Vector( 0.4000, -0.2240, 2.55000),
+  Vector( 0.2240, -0.4000, 2.55000), Vector( 0.0000, -0.4000, 2.55000),
+  Vector( 1.3000,  0.0000, 2.55000), Vector( 1.3000, -0.7280, 2.55000),
+  Vector( 0.7280, -1.3000, 2.55000), Vector( 0.0000, -1.3000, 2.55000),
+  Vector( 1.3000,  0.0000, 2.40000), Vector( 1.3000, -0.7280, 2.40000),
+  Vector( 0.7280, -1.3000, 2.40000), Vector( 0.0000, -1.3000, 2.40000),
+  Vector( 0.0000,  0.0000, 0.00000), Vector( 1.4250, -0.7980, 0.00000),
+  Vector( 1.5000,  0.0000, 0.07500), Vector( 1.4250,  0.0000, 0.00000),
+  Vector( 0.7980, -1.4250, 0.00000), Vector( 0.0000, -1.5000, 0.07500),
+  Vector( 0.0000, -1.4250, 0.00000), Vector( 1.5000, -0.8400, 0.07500),
+  Vector( 0.8400, -1.5000, 0.07500)
+]
+
+#print("static const dmnsn_vector dmnsn_teapot_verticies[%d] = {" % \
+#      len(verticies))
+#for i in range(len(verticies)):
+#  v = verticies[i]
+#  v = Vector(v.x, v.z - 1.575, v.y)
+#  print("  { .x = %+.5f, .y = %+.5f, .z = %+.5f }," % (v.x, v.y, v.z))
+#print("};")
+
+ntriangles = 0
+
+B = [
+  lambda x: (1 - x)*(1 -x)*(1 - x),
+  lambda x: 3*x*(1 - x)*(1 - x),
+  lambda x: 3*x*x*(1 - x),
+  lambda x: x*x*x,
+]
+Bprime = [
+  lambda x: -3*(1-x)*(1-x),
+  lambda x: 3*(1-x)*(1-x) - 6*x*(1-x),
+  lambda x: 6*x*(1-x) - 3*x*x,
+  lambda x: 3*x*x,
+]
+
+class _BicubicPatch(Union):
+  def __init__(self, indicies):
+    if len(indicies) != 16:
+      raise ValueError("expected 16 control points")
+
+    self.control_points = []
+    for i in indicies:
+      self.control_points.append(verticies[i])
+
+    resolution = 8
+    triangles = []
+    global ntriangles
+    for i in range(resolution):
+      for j in range(resolution):
+        a = self._eval(i/resolution, j/resolution)
+        an = self._eval_normal(i/resolution, j/resolution)
+        b = self._eval((i + 1)/resolution, j/resolution)
+        bn = self._eval_normal((i + 1)/resolution, j/resolution)
+        c = self._eval((i + 1)/resolution, (j + 1)/resolution)
+        cn = self._eval_normal((i + 1)/resolution, (j + 1)/resolution)
+        d = self._eval(i/resolution, (j + 1)/resolution)
+        dn = self._eval_normal(i/resolution, (j + 1)/resolution)
+        triangles.append(Triangle(a, b, c, an, bn, cn))
+        triangles.append(Triangle(c, d, a, cn, dn, an))
+        ntriangles += 2
+        if ntriangles%1000000 == 0:
+          print(ntriangles)
+
+    Union.__init__(self, triangles)
+
+  def _eval(self, u, v):
+    r = Vector(0)
+    for i in range(4):
+      for j in range(4):
+        r += self.control_points[4*j + i]*B[i](u)*B[j](v)
+    return r
+
+  def _eval_normal(self, u, v):
+    ru = Vector(0)
+    for i in range(4):
+      for j in range(4):
+        ru += self.control_points[4*j + i]*Bprime[i](u)*B[j](v)
+
+    rv = Vector(0)
+    for i in range(4):
+      for j in range(4):
+        rv += self.control_points[4*j + i]*B[i](u)*Bprime[j](v)
+
+    return cross(ru, rv)
+
+class _Teapot(Union):
+  def __init__(self, *args, **kwargs):
+    patches = [
+      self.reflectY(),
+      self.reflectY().scale(1, -1, 1),
+    ]
+    Union.__init__(self, patches, *args, **kwargs)
+    self.rotate(-90*X).translate(-1.575*Y)
+
+  @classmethod
+  def rim(cls):
+    return _BicubicPatch([
+      102, 103, 104, 105, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15
+    ])
+
+  @classmethod
+  def body(cls):
+    body1 = _BicubicPatch([
+      12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27
+    ])
+    body2 = _BicubicPatch([
+      24, 25, 26, 27, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40
+    ])
+    return Union([body1, body2])
+
+  @classmethod
+  def lid(cls):
+    lid1 = _BicubicPatch([
+      96, 96, 96, 96, 97, 98, 99, 100, 101, 101, 101, 101, 0, 1, 2, 3
+    ])
+    lid2 = _BicubicPatch([
+      0, 1, 2, 3, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117
+    ])
+    return Union([lid1, lid2])
+
+  @classmethod
+  def bottom(cls):
+    return _BicubicPatch([
+      118, 118, 118, 118, 124, 122, 119, 121, 123, 126, 125, 120, 40, 39, 38, 37
+    ])
+
+  @classmethod
+  def reflectXY(cls):
+    return Union([
+      cls.rim(),
+      cls.body(),
+      cls.lid(),
+      cls.bottom(),
+    ])
+
+  @classmethod
+  def handle(cls):
+    handle1 = _BicubicPatch([
+      41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56
+    ])
+    handle2 = _BicubicPatch([
+      53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 28, 65, 66, 67
+    ])
+    return Union([handle1, handle2])
+
+  @classmethod
+  def spout(cls):
+    spout1 = _BicubicPatch([
+      68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83
+    ])
+    spout2 = _BicubicPatch([
+      80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95
+    ])
+    return Union([spout1, spout2])
+
+  @classmethod
+  def reflectY(cls):
+    return Union([
+      cls.reflectXY(),
+      cls.reflectXY().scale(-1, 1, 1),
+      cls.handle(),
+      cls.spout(),
+    ])
+
+camera = PerspectiveCamera(
+  location = -5*Z,
+  look_at  = 0,
+)
+
+lights.append(PointLight(location = (-15, 20, -10), color = White))
+
+objects.append(Teapot(pigment = White).rotate(-80*X))
+
+if ntriangles > 0:
+  print(ntriangles)
diff --git a/libdimension-python/dimension.pxd b/libdimension-python/dimension.pxd
index 3f89e71..cddd7ba 100644
--- a/libdimension-python/dimension.pxd
+++ b/libdimension-python/dimension.pxd
@@ -335,6 +335,7 @@ cdef extern from "../libdimension/dimension.h":
   dmnsn_object *dmnsn_new_cube()
   dmnsn_object *dmnsn_new_cone(double r1, double r2, bint open)
   dmnsn_object *dmnsn_new_torus(double major, double minor)
+  dmnsn_object *dmnsn_new_teapot()
 
   dmnsn_object *dmnsn_new_csg_union(dmnsn_array *objects)
   dmnsn_object *dmnsn_new_csg_intersection(dmnsn_object *A, dmnsn_object *B)
diff --git a/libdimension-python/dimension.pyx b/libdimension-python/dimension.pyx
index a34d2eb..704d6fa 100644
--- a/libdimension-python/dimension.pyx
+++ b/libdimension-python/dimension.pyx
@@ -1271,6 +1271,15 @@ cdef class Torus(Object):
     self._object = dmnsn_new_torus(major_radius, minor_radius)
     Object.__init__(self, *args, **kwargs)
 
+cdef class Teapot(Object):
+  """The Utah teapot."""
+  def __init__(self, *args, **kwargs):
+    """
+    Create a Teapot.  Teapot() accepts any arguments that Object() accepts.
+    """
+    self._object = dmnsn_new_teapot()
+    Object.__init__(self, *args, **kwargs)
+
 cdef class Union(Object):
   """A CSG union."""
   def __init__(self, objects, *args, **kwargs):
diff --git a/libdimension/Makefile.am b/libdimension/Makefile.am
index aa5921e..2ee15a9 100644
--- a/libdimension/Makefile.am
+++ b/libdimension/Makefile.am
@@ -108,6 +108,7 @@ libdimension_la_SOURCES = $(nobase_include_HEADERS)                            \
                           scene.c                                              \
                           solid_pigment.c                                      \
                           sphere.c                                             \
+                          teapot.c                                             \
                           texture.c                                            \
                           threads.c                                            \
                           threads.h                                            \
diff --git a/libdimension/dimension/objects.h b/libdimension/dimension/objects.h
index 9ce9e46..02eaba8 100644
--- a/libdimension/dimension/objects.h
+++ b/libdimension/dimension/objects.h
@@ -84,3 +84,9 @@ dmnsn_object *dmnsn_new_cone(double r1, double r2, bool open);
  * @return A torus, centered at the origin and lying in the x-z plane.
  */
 dmnsn_object *dmnsn_new_torus(double major, double minor);
+
+/**
+ * The Utah teapot.
+ * @return A teapot.
+ */
+dmnsn_object *dmnsn_new_teapot(void);
diff --git a/libdimension/teapot.c b/libdimension/teapot.c
new file mode 100644
index 0000000..32c8a16
--- /dev/null
+++ b/libdimension/teapot.c
@@ -0,0 +1,455 @@
+/*************************************************************************
+ * Copyright (C) 2011 Tavian Barnes <tavianator@tavianator.com>          *
+ *                                                                       *
+ * This file is part of The Dimension Library.                           *
+ *                                                                       *
+ * The Dimension Library is free software; you can redistribute it and/  *
+ * or modify it under the terms of the GNU Lesser General Public License *
+ * as published by the Free Software Foundation; either version 3 of the *
+ * License, or (at your option) any later version.                       *
+ *                                                                       *
+ * The Dimension Library is distributed in the hope that it will be      *
+ * useful, but WITHOUT ANY WARRANTY; without even the implied warranty   *
+ * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU  *
+ * Lesser General Public License for more details.                       *
+ *                                                                       *
+ * You should have received a copy of the GNU Lesser General Public      *
+ * License along with this program.  If not, see                         *
+ * <http://www.gnu.org/licenses/>.                                       *
+ *************************************************************************/
+
+/**
+ * @file
+ * Teapots.
+ */
+
+#include "dimension-internal.h"
+
+static const dmnsn_vector dmnsn_teapot_verticies[127] = {
+  { .x = +0.20000, .y = +1.12500, .z = +0.00000 },
+  { .x = +0.20000, .y = +1.12500, .z = -0.11200 },
+  { .x = +0.11200, .y = +1.12500, .z = -0.20000 },
+  { .x = +0.00000, .y = +1.12500, .z = -0.20000 },
+  { .x = +1.33750, .y = +0.95625, .z = +0.00000 },
+  { .x = +1.33750, .y = +0.95625, .z = -0.74900 },
+  { .x = +0.74900, .y = +0.95625, .z = -1.33750 },
+  { .x = +0.00000, .y = +0.95625, .z = -1.33750 },
+  { .x = +1.43750, .y = +0.95625, .z = +0.00000 },
+  { .x = +1.43750, .y = +0.95625, .z = -0.80500 },
+  { .x = +0.80500, .y = +0.95625, .z = -1.43750 },
+  { .x = +0.00000, .y = +0.95625, .z = -1.43750 },
+  { .x = +1.50000, .y = +0.82500, .z = +0.00000 },
+  { .x = +1.50000, .y = +0.82500, .z = -0.84000 },
+  { .x = +0.84000, .y = +0.82500, .z = -1.50000 },
+  { .x = +0.00000, .y = +0.82500, .z = -1.50000 },
+  { .x = +1.75000, .y = +0.30000, .z = +0.00000 },
+  { .x = +1.75000, .y = +0.30000, .z = -0.98000 },
+  { .x = +0.98000, .y = +0.30000, .z = -1.75000 },
+  { .x = +0.00000, .y = +0.30000, .z = -1.75000 },
+  { .x = +2.00000, .y = -0.22500, .z = +0.00000 },
+  { .x = +2.00000, .y = -0.22500, .z = -1.12000 },
+  { .x = +1.12000, .y = -0.22500, .z = -2.00000 },
+  { .x = +0.00000, .y = -0.22500, .z = -2.00000 },
+  { .x = +2.00000, .y = -0.67500, .z = +0.00000 },
+  { .x = +2.00000, .y = -0.67500, .z = -1.12000 },
+  { .x = +1.12000, .y = -0.67500, .z = -2.00000 },
+  { .x = +0.00000, .y = -0.67500, .z = -2.00000 },
+  { .x = -2.00000, .y = -0.67500, .z = +0.00000 },
+  { .x = +2.00000, .y = -1.12500, .z = +0.00000 },
+  { .x = +2.00000, .y = -1.12500, .z = -1.12000 },
+  { .x = +1.12000, .y = -1.12500, .z = -2.00000 },
+  { .x = +0.00000, .y = -1.12500, .z = -2.00000 },
+  { .x = +1.50000, .y = -1.35000, .z = +0.00000 },
+  { .x = +1.50000, .y = -1.35000, .z = -0.84000 },
+  { .x = +0.84000, .y = -1.35000, .z = -1.50000 },
+  { .x = +0.00000, .y = -1.35000, .z = -1.50000 },
+  { .x = +1.50000, .y = -1.42500, .z = +0.00000 },
+  { .x = +1.50000, .y = -1.42500, .z = -0.84000 },
+  { .x = +0.84000, .y = -1.42500, .z = -1.50000 },
+  { .x = +0.00000, .y = -1.42500, .z = -1.50000 },
+  { .x = -1.60000, .y = +0.45000, .z = +0.00000 },
+  { .x = -1.60000, .y = +0.45000, .z = -0.30000 },
+  { .x = -1.50000, .y = +0.67500, .z = -0.30000 },
+  { .x = -1.50000, .y = +0.67500, .z = +0.00000 },
+  { .x = -2.30000, .y = +0.45000, .z = +0.00000 },
+  { .x = -2.30000, .y = +0.45000, .z = -0.30000 },
+  { .x = -2.50000, .y = +0.67500, .z = -0.30000 },
+  { .x = -2.50000, .y = +0.67500, .z = +0.00000 },
+  { .x = -2.70000, .y = +0.45000, .z = +0.00000 },
+  { .x = -2.70000, .y = +0.45000, .z = -0.30000 },
+  { .x = -3.00000, .y = +0.67500, .z = -0.30000 },
+  { .x = -3.00000, .y = +0.67500, .z = +0.00000 },
+  { .x = -2.70000, .y = +0.22500, .z = +0.00000 },
+  { .x = -2.70000, .y = +0.22500, .z = -0.30000 },
+  { .x = -3.00000, .y = +0.22500, .z = -0.30000 },
+  { .x = -3.00000, .y = +0.22500, .z = +0.00000 },
+  { .x = -2.70000, .y = +0.00000, .z = +0.00000 },
+  { .x = -2.70000, .y = +0.00000, .z = -0.30000 },
+  { .x = -3.00000, .y = -0.22500, .z = -0.30000 },
+  { .x = -3.00000, .y = -0.22500, .z = +0.00000 },
+  { .x = -2.50000, .y = -0.45000, .z = +0.00000 },
+  { .x = -2.50000, .y = -0.45000, .z = -0.30000 },
+  { .x = -2.65000, .y = -0.63750, .z = -0.30000 },
+  { .x = -2.65000, .y = -0.63750, .z = +0.00000 },
+  { .x = -2.00000, .y = -0.67500, .z = -0.30000 },
+  { .x = -1.90000, .y = -0.97500, .z = -0.30000 },
+  { .x = -1.90000, .y = -0.97500, .z = +0.00000 },
+  { .x = +1.70000, .y = -0.15000, .z = +0.00000 },
+  { .x = +1.70000, .y = -0.15000, .z = -0.66000 },
+  { .x = +1.70000, .y = -0.97500, .z = -0.66000 },
+  { .x = +1.70000, .y = -0.97500, .z = +0.00000 },
+  { .x = +2.60000, .y = -0.15000, .z = +0.00000 },
+  { .x = +2.60000, .y = -0.15000, .z = -0.66000 },
+  { .x = +3.10000, .y = -0.75000, .z = -0.66000 },
+  { .x = +3.10000, .y = -0.75000, .z = +0.00000 },
+  { .x = +2.30000, .y = +0.52500, .z = +0.00000 },
+  { .x = +2.30000, .y = +0.52500, .z = -0.25000 },
+  { .x = +2.40000, .y = +0.45000, .z = -0.25000 },
+  { .x = +2.40000, .y = +0.45000, .z = +0.00000 },
+  { .x = +2.70000, .y = +0.82500, .z = +0.00000 },
+  { .x = +2.70000, .y = +0.82500, .z = -0.25000 },
+  { .x = +3.30000, .y = +0.82500, .z = -0.25000 },
+  { .x = +3.30000, .y = +0.82500, .z = +0.00000 },
+  { .x = +2.80000, .y = +0.90000, .z = +0.00000 },
+  { .x = +2.80000, .y = +0.90000, .z = -0.25000 },
+  { .x = +3.52500, .y = +0.91875, .z = -0.25000 },
+  { .x = +3.52500, .y = +0.91875, .z = +0.00000 },
+  { .x = +2.90000, .y = +0.90000, .z = +0.00000 },
+  { .x = +2.90000, .y = +0.90000, .z = -0.15000 },
+  { .x = +3.45000, .y = +0.93750, .z = -0.15000 },
+  { .x = +3.45000, .y = +0.93750, .z = +0.00000 },
+  { .x = +2.80000, .y = +0.82500, .z = +0.00000 },
+  { .x = +2.80000, .y = +0.82500, .z = -0.15000 },
+  { .x = +3.20000, .y = +0.82500, .z = -0.15000 },
+  { .x = +3.20000, .y = +0.82500, .z = +0.00000 },
+  { .x = +0.00000, .y = +1.57500, .z = +0.00000 },
+  { .x = +0.80000, .y = +1.57500, .z = +0.00000 },
+  { .x = +0.80000, .y = +1.57500, .z = -0.45000 },
+  { .x = +0.45000, .y = +1.57500, .z = -0.80000 },
+  { .x = +0.00000, .y = +1.57500, .z = -0.80000 },
+  { .x = +0.00000, .y = +1.27500, .z = +0.00000 },
+  { .x = +1.40000, .y = +0.82500, .z = +0.00000 },
+  { .x = +1.40000, .y = +0.82500, .z = -0.78400 },
+  { .x = +0.78400, .y = +0.82500, .z = -1.40000 },
+  { .x = +0.00000, .y = +0.82500, .z = -1.40000 },
+  { .x = +0.40000, .y = +0.97500, .z = +0.00000 },
+  { .x = +0.40000, .y = +0.97500, .z = -0.22400 },
+  { .x = +0.22400, .y = +0.97500, .z = -0.40000 },
+  { .x = +0.00000, .y = +0.97500, .z = -0.40000 },
+  { .x = +1.30000, .y = +0.97500, .z = +0.00000 },
+  { .x = +1.30000, .y = +0.97500, .z = -0.72800 },
+  { .x = +0.72800, .y = +0.97500, .z = -1.30000 },
+  { .x = +0.00000, .y = +0.97500, .z = -1.30000 },
+  { .x = +1.30000, .y = +0.82500, .z = +0.00000 },
+  { .x = +1.30000, .y = +0.82500, .z = -0.72800 },
+  { .x = +0.72800, .y = +0.82500, .z = -1.30000 },
+  { .x = +0.00000, .y = +0.82500, .z = -1.30000 },
+  { .x = +0.00000, .y = -1.57500, .z = +0.00000 },
+  { .x = +1.42500, .y = -1.57500, .z = -0.79800 },
+  { .x = +1.50000, .y = -1.50000, .z = +0.00000 },
+  { .x = +1.42500, .y = -1.57500, .z = +0.00000 },
+  { .x = +0.79800, .y = -1.57500, .z = -1.42500 },
+  { .x = +0.00000, .y = -1.50000, .z = -1.50000 },
+  { .x = +0.00000, .y = -1.57500, .z = -1.42500 },
+  { .x = +1.50000, .y = -1.50000, .z = -0.84000 },
+  { .x = +0.84000, .y = -1.50000, .z = -1.50000 },
+};
+
+#define DMNSN_BEZIER_RESOLUTION 8
+
+#define DMNSN_TEAPOT_NTRIANGLES                                 \
+  (32*2*DMNSN_BEZIER_RESOLUTION*DMNSN_BEZIER_RESOLUTION)
+
+dmnsn_vector dmnsn_teapot_triangles[DMNSN_TEAPOT_NTRIANGLES][6];
+
+static inline double
+dmnsn_Bernstein_value(size_t i, double x)
+{
+  switch (i) {
+  case 0:
+    return (1 - x)*(1 - x)*(1 - x);
+  case 1:
+    return 3.0*x*(1 - x)*(1 - x);
+  case 2:
+    return 3.0*x*x*(1 - x);
+  case 3:
+    return x*x*x;
+  default:
+    dmnsn_unreachable("Wrong Bernstein index.");
+  }
+}
+
+static inline double
+dmnsn_Bernstein_prime(size_t i, double x)
+{
+  switch (i) {
+  case 0:
+    return -3.0*(1 - x)*(1 - x);
+  case 1:
+    return 3.0*(1 - x)*(1 - x) - 6.0*x*(1 - x);
+  case 2:
+    return 6.0*x*(1 - x) - 3.0*x*x;
+  case 3:
+    return 3.0*x*x;
+  default:
+    dmnsn_unreachable("Wrong Bernstein index.");
+  }
+}
+
+static inline dmnsn_vector
+dmnsn_bezier_patch_point(const dmnsn_vector control_points[16],
+                         double u, double v)
+{
+  dmnsn_vector r = dmnsn_zero;
+  for (size_t i = 0; i < 4; ++i) {
+    for (size_t j = 0; j < 4; ++j) {
+      r = dmnsn_vector_add(
+        r,
+        dmnsn_vector_mul(
+          dmnsn_Bernstein_value(i, u)*dmnsn_Bernstein_value(j, v),
+          control_points[4*j + i]
+        )
+      );
+    }
+  }
+  return r;
+}
+
+static inline dmnsn_vector
+dmnsn_bezier_patch_normal(const dmnsn_vector control_points[16],
+                          double u, double v)
+{
+  dmnsn_vector nu = dmnsn_zero, nv = dmnsn_zero;
+  for (size_t i = 0; i < 4; ++i) {
+    for (size_t j = 0; j < 4; ++j) {
+      nu = dmnsn_vector_add(
+        nu,
+        dmnsn_vector_mul(
+          dmnsn_Bernstein_prime(i, u)*dmnsn_Bernstein_value(j, v),
+          control_points[4*j + i]
+        )
+      );
+
+      nv = dmnsn_vector_add(
+        nv,
+        dmnsn_vector_mul(
+          dmnsn_Bernstein_value(i, u)*dmnsn_Bernstein_prime(j, v),
+          control_points[4*j + i]
+        )
+      );
+    }
+  }
+
+  dmnsn_vector n = dmnsn_vector_cross(nu, nv);
+  if (dmnsn_vector_norm(n) < dmnsn_epsilon) {
+    printf("%g\t%g\t" DMNSN_VECTOR_FORMAT "\t" DMNSN_VECTOR_FORMAT "\n",
+           u, v, DMNSN_VECTOR_PRINTF(nu), DMNSN_VECTOR_PRINTF(nv));
+    for (size_t i = 0; i < 16; ++i) {
+      printf("\t" DMNSN_VECTOR_FORMAT "\n", DMNSN_VECTOR_PRINTF(control_points[i]));
+    }
+    nu = dmnsn_zero;
+    for (size_t i = 0; i < 4; ++i) {
+      for (size_t j = 0; j < 4; ++j) {
+        nu = dmnsn_vector_add(
+          nu,
+          dmnsn_vector_mul(
+            dmnsn_Bernstein_prime(i, u)*dmnsn_Bernstein_value(j, v),
+            control_points[4*j + i]
+          )
+        );
+        if (dmnsn_vector_norm(dmnsn_vector_mul(
+            dmnsn_Bernstein_prime(i, u)*dmnsn_Bernstein_value(j, v),
+            control_points[4*j + i]
+                                               )) > dmnsn_epsilon)
+        printf(DMNSN_VECTOR_FORMAT "\n", DMNSN_VECTOR_PRINTF(dmnsn_vector_mul(
+            dmnsn_Bernstein_prime(i, u)*dmnsn_Bernstein_value(j, v),
+            control_points[4*j + i]
+          )));
+      }
+    }
+    dmnsn_assert(false, "asdf.");
+  }
+  return n;
+}
+
+static size_t
+dmnsn_mirror_last_triangle(size_t i, bool mirror_x)
+{
+  return i;
+  for (size_t j = 0; j < 6; ++j) {
+    dmnsn_vector v = dmnsn_teapot_triangles[i - 1][j];
+    if (mirror_x) {
+      v.x = -v.x;
+    } else {
+      v.z = -v.z;
+    }
+
+    dmnsn_teapot_triangles[i][j] = v;
+  }
+
+  return i + 1;
+}
+
+static size_t
+dmnsn_subdivide_patch(const size_t indicies[16], size_t i,
+                      bool mirror_x)
+{
+  dmnsn_vector control_points[16];
+  for (size_t j = 0; j < 16; ++j) {
+    control_points[j] = dmnsn_teapot_verticies[indicies[j]];
+  }
+
+  dmnsn_vector lastrowr[DMNSN_BEZIER_RESOLUTION + 1];
+  dmnsn_vector lastrown[DMNSN_BEZIER_RESOLUTION + 1];
+  for (size_t k = 0; k <= DMNSN_BEZIER_RESOLUTION; ++k) {
+    lastrowr[k] = dmnsn_bezier_patch_point(
+      control_points, 0.0, (double)k/DMNSN_BEZIER_RESOLUTION
+    );
+    lastrown[k] = dmnsn_bezier_patch_normal(
+      control_points, 0.0, (double)k/DMNSN_BEZIER_RESOLUTION
+    );
+  }
+
+  dmnsn_vector ar, an, dr, dn;
+  for (size_t j = 1; j <= DMNSN_BEZIER_RESOLUTION; ++j) {
+    ar = lastrowr[0];
+    an = lastrown[0];
+    dr = dmnsn_bezier_patch_point(
+      control_points, (double)j/DMNSN_BEZIER_RESOLUTION, 0.0
+    );
+    dn = dmnsn_bezier_patch_normal(
+      control_points, (double)j/DMNSN_BEZIER_RESOLUTION, 0.0
+    );
+    lastrowr[0] = dr;
+    lastrown[0] = dn;
+
+    for (size_t k = 1; k <= DMNSN_BEZIER_RESOLUTION; ++k) {
+      dmnsn_vector br = lastrowr[k];
+      dmnsn_vector bn = lastrown[k];
+      dmnsn_vector cr = dmnsn_bezier_patch_point(
+        control_points,
+        (double)j/DMNSN_BEZIER_RESOLUTION, (double)k/DMNSN_BEZIER_RESOLUTION
+      );
+      dmnsn_vector cn = dmnsn_bezier_patch_normal(
+        control_points,
+        (double)j/DMNSN_BEZIER_RESOLUTION, (double)k/DMNSN_BEZIER_RESOLUTION
+      );
+
+      dmnsn_teapot_triangles[i][0] = ar;
+      dmnsn_teapot_triangles[i][1] = br;
+      dmnsn_teapot_triangles[i][2] = cr;
+      dmnsn_teapot_triangles[i][3] = an;
+      dmnsn_teapot_triangles[i][4] = bn;
+      dmnsn_teapot_triangles[i][5] = cn;
+      ++i;
+      i = dmnsn_mirror_last_triangle(i, false);
+      if (mirror_x) {
+        i = dmnsn_mirror_last_triangle(i, true);
+        i = dmnsn_mirror_last_triangle(i, false);
+      }
+
+      dmnsn_teapot_triangles[i][0] = cr;
+      dmnsn_teapot_triangles[i][1] = dr;
+      dmnsn_teapot_triangles[i][2] = ar;
+      dmnsn_teapot_triangles[i][3] = cn;
+      dmnsn_teapot_triangles[i][4] = dn;
+      dmnsn_teapot_triangles[i][5] = an;
+      ++i;
+      i = dmnsn_mirror_last_triangle(i, false);
+      if (mirror_x) {
+        i = dmnsn_mirror_last_triangle(i, true);
+        i = dmnsn_mirror_last_triangle(i, false);
+      }
+
+      ar = br;
+      an = bn;
+      dr = cr;
+      dn = cn;
+      lastrowr[k] = dr;
+      lastrown[k] = dn;
+    }
+  }
+
+  return i;
+}
+
+static const size_t dmnsn_teapot_rim[16] = {
+  102, 103, 104, 105, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15
+};
+static const size_t dmnsn_teapot_body1[16] = {
+  12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27
+};
+static const size_t dmnsn_teapot_body2[16] = {
+  24, 25, 26, 27, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40
+};
+static const size_t dmnsn_teapot_lid1[16] = {
+  96, 96, 96, 96, 97, 98, 99, 100, 101, 101, 101, 101, 0, 1, 2, 3
+};
+static const size_t dmnsn_teapot_lid2[16] = {
+  0, 1, 2, 3, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117
+};
+static const size_t dmnsn_teapot_bottom[16] = {
+  118, 118, 118, 118, 124, 122, 119, 121, 123, 126, 125, 120, 40, 39, 38, 37
+};
+static const size_t dmnsn_teapot_handle1[16] = {
+  41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56
+};
+static const size_t dmnsn_teapot_handle2[16] = {
+  53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 28, 65, 66, 67
+};
+static const size_t dmnsn_teapot_spout1[16] = {
+  68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83
+};
+static const size_t dmnsn_teapot_spout2[16] = {
+  80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95
+};
+
+static void
+dmnsn_teapot_calculate(void)
+{
+  size_t i = 0;
+  i = dmnsn_subdivide_patch(dmnsn_teapot_rim,     i, true);
+  i = dmnsn_subdivide_patch(dmnsn_teapot_body1,   i, true);
+  i = dmnsn_subdivide_patch(dmnsn_teapot_body2,   i, true);
+  i = dmnsn_subdivide_patch(dmnsn_teapot_lid1,    i, true);
+  i = dmnsn_subdivide_patch(dmnsn_teapot_lid2,    i, true);
+  i = dmnsn_subdivide_patch(dmnsn_teapot_bottom,  i, true);
+  i = dmnsn_subdivide_patch(dmnsn_teapot_handle1, i, false);
+  i = dmnsn_subdivide_patch(dmnsn_teapot_handle2, i, false);
+  i = dmnsn_subdivide_patch(dmnsn_teapot_spout1,  i, false);
+  i = dmnsn_subdivide_patch(dmnsn_teapot_spout2,  i, false);
+  //dmnsn_assert(i == DMNSN_TEAPOT_NTRIANGLES, "Wrong number of triangles.");
+}
+
+pthread_mutex_t dmnsn_teapot_mutex = PTHREAD_MUTEX_INITIALIZER;
+
+dmnsn_object *
+dmnsn_new_teapot(void)
+{
+  dmnsn_lock_mutex(&dmnsn_teapot_mutex);
+    static bool verticies_calculated = false;
+    if (!verticies_calculated) {
+      dmnsn_teapot_calculate();
+      verticies_calculated = true;
+    }
+  dmnsn_unlock_mutex(&dmnsn_teapot_mutex);
+
+  dmnsn_array *triangles = dmnsn_new_array(sizeof(dmnsn_object *));
+  for (size_t i = 0; i < DMNSN_TEAPOT_NTRIANGLES; ++i) {
+    dmnsn_vector *vertnorms = dmnsn_teapot_triangles[i];
+    dmnsn_object *triangle = dmnsn_new_triangle(
+      vertnorms[0], vertnorms[1], vertnorms[2],
+      vertnorms[3], vertnorms[4], vertnorms[5]
+    );
+    triangle->texture = dmnsn_new_texture();
+    dmnsn_color randcolor;
+    randcolor.R = drand48();
+    randcolor.G = drand48();
+    randcolor.B = drand48();
+    triangle->texture->pigment = dmnsn_new_solid_pigment(DMNSN_TCOLOR(randcolor));
+    dmnsn_array_push(triangles, &triangle);
+  }
+  dmnsn_object *teapot = dmnsn_new_csg_union(triangles);
+  dmnsn_delete_array(triangles);
+  return teapot;
+}