From ba58f3c9881bcb5c61d8325f7e7d225ff9d27ab0 Mon Sep 17 00:00:00 2001
From: Carlos Scheidegger <cscheid@research.att.com>
Date: Thu, 20 Feb 2014 17:06:30 -0500
Subject: [PATCH 01/16] moving development to 500lines repo

---
 rasterizer/README.md              |  49 ++++++++++++++
 rasterizer/rasterizer/__init__.py |   8 +++
 rasterizer/rasterizer/color.py    |  26 ++++++++
 rasterizer/rasterizer/csg.py      |  45 +++++++++++++
 rasterizer/rasterizer/ellipse.py  |  67 +++++++++++++++++++
 rasterizer/rasterizer/geometry.py | 106 ++++++++++++++++++++++++++++++
 rasterizer/rasterizer/image.py    |  24 +++++++
 rasterizer/rasterizer/poly.py     |  46 +++++++++++++
 rasterizer/rasterizer/scene.py    |  24 +++++++
 rasterizer/rasterizer/shape.py    |  47 +++++++++++++
 rasterizer/test_rasterizer.py     |  69 +++++++++++++++++++
 11 files changed, 511 insertions(+)
 create mode 100644 rasterizer/README.md
 create mode 100644 rasterizer/rasterizer/__init__.py
 create mode 100644 rasterizer/rasterizer/color.py
 create mode 100644 rasterizer/rasterizer/csg.py
 create mode 100644 rasterizer/rasterizer/ellipse.py
 create mode 100644 rasterizer/rasterizer/geometry.py
 create mode 100644 rasterizer/rasterizer/image.py
 create mode 100644 rasterizer/rasterizer/poly.py
 create mode 100644 rasterizer/rasterizer/scene.py
 create mode 100644 rasterizer/rasterizer/shape.py
 create mode 100644 rasterizer/test_rasterizer.py

diff --git a/rasterizer/README.md b/rasterizer/README.md
new file mode 100644
index 000000000..ab163378f
--- /dev/null
+++ b/rasterizer/README.md
@@ -0,0 +1,49 @@
+# `tiny_gfx`: A tiny rasterizer
+
+A *rasterizer* is a piece of software that converts arbitrary shapes
+into *raster images*, which is just a funny name for rectangular grids
+of pixels. A rasterizer, in some way or another, is at the heart of
+pretty much every modern display technology today, from computer
+displays to e-ink displays to printers both 2D and 3D (the graphics in
+laser shows are the most notable exception).
+
+In this chapter, I will teach you a little about how rasterizers work,
+by describing `tiny_gfx`, a simple rasterizer in pure Python. Along
+the way we will pick up some techniques that show up repeatedly in
+computer graphics code. Having a bit of mathematical background on
+linear algebra will help, but I hope the presentation will be
+self-contained.
+
+`tiny_gfx` is not practical in many ways. Besides being slow (see
+below), the main shortcoming in `tiny_gfx` is that shapes are all of a
+single solid color. Still, for 500 lines of code, `tiny_gfx` has a
+relatively large set of features:
+
+- alpha blending for semi-transparent colors
+- polygonal shapes (possibly concave, but not self-intersecting)
+  (TODO: right now only convex polygons are supported)
+- circles and ellipses
+- transformations of shapes
+- boolean operations on shapes (union, intersection, subtraction)
+- antialiasing
+- empty space skipping and fast rasterization of runs for general
+  shapes
+
+Maybe most interestingly, shapes in `tiny_gfx` are extensible. New
+shapes are easily added, and they compose well with the other parts of
+the code.
+
+## A performance caveat
+
+Rasterizers are so central to display technology that their
+performance can make or break a piece of software, and these days the
+fastest rasterizers are all based in hardware. Your videogame graphics
+card (and even the cheapest smartphones these days) is rasterizing
+polygons in highly parallel processors; 192 cores is a typical
+number. It should be no surprise, then, that the rasterizer we will
+see here is slow: if CPU-intensive tasks in Python run around 50 times
+slower than heavily-optimized, low-level code, and if a graphics
+driver has around 200 cores at its disposal, a slowdown of 10,000
+times should not be surprising. In reality, `tiny_gfx` is closer to
+1,000,000 times slower than the special-purpose graphics rasterizer
+from the laptop in which I'm developing it.
diff --git a/rasterizer/rasterizer/__init__.py b/rasterizer/rasterizer/__init__.py
new file mode 100644
index 000000000..2c142778f
--- /dev/null
+++ b/rasterizer/rasterizer/__init__.py
@@ -0,0 +1,8 @@
+from shape import *
+from geometry import *
+from poly import *
+from color import *
+from ellipse import *
+from image import *
+from scene import *
+from csg import *
diff --git a/rasterizer/rasterizer/color.py b/rasterizer/rasterizer/color.py
new file mode 100644
index 000000000..7d328cc83
--- /dev/null
+++ b/rasterizer/rasterizer/color.py
@@ -0,0 +1,26 @@
+class Color:
+    def __init__(self, r=0, g=0, b=0, a=1, rgb=None):
+        self.rgb = rgb or (r, g, b)
+        self.a = a
+    def draw(self, o):
+        if self.a == o.a == 0.0:
+            return
+        if o.a == 1.0:
+            self.rgb = o.rgb
+            self.a = 1
+        else:
+            u = 1.0 - o.a
+            self.rgb = (u * self.rgb[0] + o.a * o.rgb[0],
+                        u * self.rgb[1] + o.a * o.rgb[1],
+                        u * self.rgb[2] + o.a * o.rgb[2])
+            self.a = 1.0 - (1.0 - self.a) * (1.0 - o.a)
+    def fainter(self, k):
+        return Color(rgb=self.rgb, a=self.a*k)
+    def as_ppm(self):
+        def byte(v):
+            return int(v ** (1.0 / 2.2) * 255)
+        return "%c%c%c" % (byte(self.rgb[0] * self.a),
+                           byte(self.rgb[1] * self.a),
+                           byte(self.rgb[2] * self.a))
+    def __repr__(self):
+        return "[" + str(self.rgb) + "," + str(self.a) + "]"
diff --git a/rasterizer/rasterizer/csg.py b/rasterizer/rasterizer/csg.py
new file mode 100644
index 000000000..47dc3795f
--- /dev/null
+++ b/rasterizer/rasterizer/csg.py
@@ -0,0 +1,45 @@
+from shape import Shape
+from geometry import *
+
+class CSG(Shape):
+    def __init__(self, v1, v2, color=None):
+        Shape.__init__(self, color or v1.color or v2.color)
+        self.v1 = v1
+        self.v2 = v2
+    def transform(self, t):
+        return self.__class__(self.v1.transform(t), self.v2.transform(t),
+                              color=self.color)
+
+class Union(CSG):
+    def __init__(self, v1, v2, color=None):
+        CSG.__init__(self, v1, v2, color=color)
+        self.bound = Vector.union(v1.bound.low, v1.bound.high,
+                                  v2.bound.low, v2.bound.high)
+    def contains(self, p):
+        return self.v1.contains(p) or self.v2.contains(p)
+    def signed_distance_bound(self, p):
+        b1 = self.v1.signed_distance_bound(p)
+        b2 = self.v2.signed_distance_bound(p)
+        return b1 if b1 > b2 else b2
+
+class Intersection(CSG):
+    def __init__(self, v1, v2, color=None):
+        CSG.__init__(self, v1, v2, color=color)
+        self.bound = v1.bound.intersection(v2.bound)
+    def contains(self, p):
+        return self.v1.contains(p) and self.v2.contains(p)
+    def signed_distance_bound(self, p):
+        b1 = self.v1.signed_distance_bound(p)
+        b2 = self.v2.signed_distance_bound(p)
+        return b1 if b1 < b2 else b2
+
+class Subtraction(CSG):
+    def __init__(self, v1, v2, color=None):
+        CSG.__init__(self, v1, v2, color=color)
+        self.bound = self.v1.bound
+    def contains(self, p):
+        return self.v1.contains(p) and not self.v2.contains(p)
+    def signed_distance_bound(self, p):
+        b1 = self.v1.signed_distance_bound(p)
+        b2 = -self.v2.signed_distance_bound(p)
+        return b1 if b1 < b2 else b2
diff --git a/rasterizer/rasterizer/ellipse.py b/rasterizer/rasterizer/ellipse.py
new file mode 100644
index 000000000..4037c1327
--- /dev/null
+++ b/rasterizer/rasterizer/ellipse.py
@@ -0,0 +1,67 @@
+from shape import Shape
+from geometry import Vector, Transform, quadratic, scale, translate
+from color import Color
+
+class Ellipse(Shape):
+    def __init__(self, a=1.0, b=1.0, c=0.0, d=0.0, e=0.0, f=-1.0, color=None):
+        Shape.__init__(self, color)
+        self.a = a
+        self.b = b
+        self.c = c
+        self.d = d
+        self.e = e
+        self.f = f
+        t = Transform(2 * a, c, 0, c, 2 * b, 0)
+        self.center = t.inverse() * Vector(-d, -e)
+        l1, l0 = quadratic(1, 2 * (-a - b), 4 * a * b - c * c)
+        v = t.eigv()
+        axes = [v[0] * ((l0 / 2) ** -0.5), v[1] * ((l1 / 2) ** -0.5)]
+        self.bound = Vector.union(self.center - axes[0] - axes[1],
+                                  self.center - axes[0] + axes[1],
+                                  self.center + axes[0] - axes[1],
+                                  self.center + axes[0] + axes[1])
+    def value(self, p):
+        return self.a*p.x*p.x + self.b*p.y*p.y + self.c*p.x*p.y \
+               + self.d*p.x + self.e*p.y + self.f
+    def contains(self, p):
+        return self.value(p) < 0
+    def transform(self, transform):
+        i = transform.inverse()
+        ((m00, m01, m02), (m10, m11, m12),_) = i.m
+        aa = self.a*m00*m00 + self.b*m10*m10 + self.c*m00*m10
+        bb = self.a*m01*m01 + self.b*m11*m11 + self.c*m01*m11
+        cc = 2*self.a*m00*m01 + 2*self.b*m10*m11 \
+             + self.c*(m00*m11 + m01*m10)
+        dd = 2*self.a*m00 * m02 + 2*self.b*m10*m12 \
+             + self.c*(m00*m12 + m02*m10) + self.d*m00 + self.e*m10
+        ee = 2*self.a*m10*m02 + 2*self.b*m11*m12 \
+             + self.c*(m01*m12 + m02*m11) + self.d*m01 + self.e*m11
+        ff = self.a*m02*m02 + self.b*m12*m12 + self.c*m02*m12 \
+             + self.d*m02 + self.e*m12 + self.f
+        return Ellipse(aa, bb, cc, dd, ee, ff, color=self.color)
+    def signed_distance_bound(self, p):
+        def sgn(x):
+            return 0 if x == 0 else x / abs(x)
+        v = -sgn(self.value(p))
+        c = self.center
+        pc = p - c
+        u2 = self.a*pc.x**2 + self.b*pc.y**2 + self.c*pc.x*pc.y
+        u1 = 2*self.a*c.x*pc.x + 2*self.b*c.y*pc.y \
+             + self.c*c.y*pc.x + self.c*c.x*pc.y + self.d*pc.x \
+             + self.e*pc.y
+        u0 = self.a*c.x**2 + self.b*c.y**2 + self.c*c.x*c.y \
+             + self.d*c.x + self.e*c.y + self.f
+        sols = quadratic(u2, u1, u0)
+        crossings = c+pc*sols[0], c+pc*sols[1]
+        if (p - crossings[0]).length() < (p - crossings[1]).length():
+            surface_pt = crossings[0]
+        else:
+            surface_pt = crossings[1]
+        d = Vector(2*self.a*surface_pt.x + self.c*surface_pt.y + self.d,
+                   2*self.b*surface_pt.y + self.c*surface_pt.x + self.e)
+        return v * abs(d.dot(p - surface_pt) / d.length())
+
+def Circle(center, radius, color=None):
+    return Ellipse(color=color).transform(
+        scale(radius, radius)).transform(
+        translate(center.x, center.y))
diff --git a/rasterizer/rasterizer/geometry.py b/rasterizer/rasterizer/geometry.py
new file mode 100644
index 000000000..35f7d5c49
--- /dev/null
+++ b/rasterizer/rasterizer/geometry.py
@@ -0,0 +1,106 @@
+import random
+import math
+from itertools import product
+
+def quadratic(a, b, c):
+    d = (b * b - 4 * a * c) ** 0.5
+    if b >= 0:
+        return (-b - d) / (2 * a), (2 * c) / (-b - d)
+    else:
+        return (2 * c) / (-b + d), (-b + d) / (2 * a)
+
+class Vector:
+    def __init__(self, *args):
+        self.x, self.y = args
+    def __add__(self, o):
+        return Vector(self.x + o.x, self.y + o.y)
+    def __sub__(self, o):
+        return Vector(self.x - o.x, self.y - o.y)
+    def __mul__(self, k):
+        return Vector(self.x * k, self.y * k)
+    def dot(self, o):
+        return self.x * o.x + self.y * o.y
+    def min(self, o):
+        return Vector(min(self.x, o.x), min(self.y, o.y))
+    def max(self, o):
+        return Vector(max(self.x, o.x), max(self.y, o.y))
+    def union(*args):
+        return AABox(reduce(Vector.min, args), reduce(Vector.max, args))
+    def length(self):
+        return (self.x * self.x + self.y * self.y) ** 0.5
+
+class AABox:
+    def __init__(self, p1, p2):
+        self.low = p1.min(p2)
+        self.high = p1.max(p2)
+    def midpoint(self):
+        return (self.low + self.high) * 0.5
+    def size(self):
+        return self.high - self.low
+    def contains(self, p):
+        return self.low.x <= p.x <= self.high.x and \
+               self.low.y <= p.y <= self.high.y
+    def overlaps(self, r):
+        return not (r.low.x >= self.high.x or r.high.x <= self.low.x or
+                    r.low.y >= self.high.y or r.high.y <= self.low.y)
+    def intersection(self, other):
+        return AABox(self.low.max(other.low), self.high.min(other.high))
+
+class HalfPlane:
+    def __init__(self, p1, p2):
+        self.v = Vector(-p2.y + p1.y, p2.x - p1.x)
+        l = self.v.length()
+        self.c = -self.v.dot(p1) / l
+        self.v = self.v * (1.0 / l)
+    def signed_distance(self, p):
+        return self.v.dot(p) + self.c
+
+class Transform:
+    def __init__(self, m11, m12, tx, m21, m22, ty):
+        self.m = [[m11, m12, tx],
+                  [m21, m22, ty],
+                  [0, 0, 1]]
+    def __mul__(self, other): # ugly
+        if isinstance(other, Transform):
+            t = [[0] * 3 for i in xrange(3)]
+            for i, j, k in product(xrange(3), repeat=3):
+                t[i][j] += self.m[i][k] * other.m[k][j]
+            return Transform(t[0][0], t[0][1], t[0][2],
+                             t[1][0], t[1][1], t[1][2])
+        else:
+            nx = self.m[0][0] * other.x + self.m[0][1] * other.y + self.m[0][2]
+            ny = self.m[1][0] * other.x + self.m[1][1] * other.y + self.m[1][2]
+            return Vector(nx, ny)
+    def det(s):
+        return s.m[0][0] * s.m[1][1] - s.m[0][1] * s.m[1][0]
+    def inverse(self):
+        d = 1.0 / self.det()
+        return Transform(d * self.m[1][1], -d * self.m[0][1], -self.m[0][2],
+                         -d * self.m[1][0], d * self.m[0][0], -self.m[1][2])
+    def eigv(self): # power iteration, ignores translation, assumes SPD
+        a = Vector(random.random(), random.random())
+        last = Vector(0,0)
+        t = Transform(self.m[0][0], self.m[0][1], 0,
+                      self.m[1][0], self.m[1][1], 0)
+        while (a - last).length() > 1e-6:
+            last = a
+            a = t * a
+            a = a * (1.0 / a.length())
+        return a, Vector(-a.y, a.x)
+
+def identity():
+    return Transform(1, 0, 0, 0, 1, 0)
+
+def rotate(theta):
+    s = math.sin(theta)
+    c = math.cos(theta)
+    return Transform(c, -s, 0, s, c, 0)
+
+def translate(tx, ty):
+    return Transform(1, 0, tx, 0, 1, ty)
+
+def scale(x, y):
+    return Transform(x, 0, 0, 0, y, 0)
+
+def around(v, t):
+    return translate(v.x, v.y) * t * translate(-v.x, -v.y)
diff --git a/rasterizer/rasterizer/image.py b/rasterizer/rasterizer/image.py
new file mode 100644
index 000000000..9c6de0b11
--- /dev/null
+++ b/rasterizer/rasterizer/image.py
@@ -0,0 +1,24 @@
+from color import Color
+from geometry import AABox, Vector
+
+class PPMImage:
+    def __init__(self, resolution, bg=Color()):
+        self.resolution = resolution
+        self.pixels = []
+        for i in xrange(self.resolution):
+            lst = []
+            for j in xrange(self.resolution):
+                lst.append(Color(rgb=bg.rgb, a=bg.a))
+            self.pixels.append(lst)
+    def bounds(self):
+        return AABox(Vector(0,0), Vector(1,1))
+    def __getitem__(self, a):
+        return self.pixels[a.y][a.x]
+    def __setitem__(self, a, color):
+        self.pixels[a.y][a.x] = color
+    def write_ppm(self, out):
+        n = self.resolution
+        out.write("P6\n%s\n%s\n255\n" % (n,n))
+        for y in xrange(n-1, -1, -1):
+            for x in xrange(n):
+                out.write(self.pixels[y][x].as_ppm())
diff --git a/rasterizer/rasterizer/poly.py b/rasterizer/rasterizer/poly.py
new file mode 100644
index 000000000..5a409cbd8
--- /dev/null
+++ b/rasterizer/rasterizer/poly.py
@@ -0,0 +1,46 @@
+from shape import Shape
+from geometry import HalfPlane, Vector
+
+class Poly(Shape): # a *convex* poly, in ccw order, with no repeating vertices
+    def __init__(self, ps, color=None):
+        Shape.__init__(self, color)
+        mn = min(enumerate(ps), key=lambda (i,v): (v.y, v.x))[0]
+        self.vs = list(ps[mn:]) + list(ps[:mn])
+        self.bound = Vector.union(*self.vs)
+        self.half_planes = []
+        for i in xrange(len(self.vs)):
+            h = HalfPlane(self.vs[i], self.vs[(i+1) % len(self.vs)])
+            self.half_planes.append(h)
+    def signed_distance_bound(self, p):
+        plane = self.half_planes[0]
+        min_inside = 1e30
+        max_outside = -1e30
+        for plane in self.half_planes:
+            d = plane.signed_distance(p)
+            if d <= 0 and d > max_outside:
+                max_outside = d
+            if d >= 0 and d < min_inside:
+                min_inside = d
+        return max_outside if max_outside <> -1e30 else min_inside
+    def contains(self, p):
+        for plane in self.half_planes:
+            if plane.signed_distance(p) < 0:
+                return False
+        return True
+    def transform(self, xform):
+        return Poly(list(xform * v for v in self.vs), color=self.color)
+
+Triangle, Quad = Poly, Poly
+
+def Rectangle(v1, v2, color=None):
+    return Quad([Vector(min(v1.x, v2.x), min(v1.y, v2.y)),
+                 Vector(max(v1.x, v2.x), min(v1.y, v2.y)),
+                 Vector(max(v1.x, v2.x), max(v1.y, v2.y)),
+                 Vector(min(v1.x, v2.x), max(v1.y, v2.y))],
+                color=color)
+
+def LineSegment(v1, v2, thickness, color=None):
+    d = v2 - v1
+    d.x, d.y = -d.y, d.x
+    d *= thickness / d.length() / 2
+    return Quad([v1 + d, v1 - d, v2 - d, v2 + d], color=color)
diff --git a/rasterizer/rasterizer/scene.py b/rasterizer/rasterizer/scene.py
new file mode 100644
index 000000000..8c31235e3
--- /dev/null
+++ b/rasterizer/rasterizer/scene.py
@@ -0,0 +1,24 @@
+from geometry import identity
+from shape import Shape
+
+class Scene:
+    def __init__(self, nodes=None, transform=None):
+        if transform is None:
+            transform = identity()
+        if nodes is None:
+            nodes = []
+        self.transform = transform
+        self.nodes = nodes
+    def add(self, node):
+        self.nodes.append(node)
+    def draw(self, image):
+        for shape in self.traverse(identity()):
+            shape.draw(image)
+    def traverse(self, xform):
+        this_xform = xform * self.transform
+        for node in self.nodes:
+            if isinstance(node, Scene):
+                for n in node.traverse(this_xform):
+                    yield n
+            elif isinstance(node, Shape):
+                yield node.transform(this_xform)
diff --git a/rasterizer/rasterizer/shape.py b/rasterizer/rasterizer/shape.py
new file mode 100644
index 000000000..505d45b19
--- /dev/null
+++ b/rasterizer/rasterizer/shape.py
@@ -0,0 +1,47 @@
+from color import Color
+from itertools import product
+from geometry import Vector
+import random
+
+class Shape:
+    def __init__(self, color=None):
+        self.color = color if color is not None else Color()
+    def draw(self, image, super_sampling = 6):
+        color = self.color
+        r = float(image.resolution)
+        jitter = [Vector((x + random.random()) / super_sampling / r,
+                         (y + random.random()) / super_sampling / r)
+                  for (x, y) in product(xrange(super_sampling), repeat=2)]
+        lj = len(jitter)
+        total_pixels = 0
+        tb = self.bound
+        if not tb.overlaps(image.bounds()):
+            return
+        l_x = max(0, int(tb.low.x * r))
+        l_y = max(0, int(tb.low.y * r))
+        h_x = min(r-1, int(tb.high.x * r))
+        h_y = min(r-1, int(tb.high.y * r))
+        corners = [(0, 0), (1.0/r, 0), (0, 1.0/r), (1.0/r, 1.0/r)]
+        for y in xrange(l_y, int(h_y+1)):
+            x = l_x
+            while x <= h_x:
+                corner = Vector(x / r, y / r)
+                b = self.signed_distance_bound(corner) * r
+                if b > 1.414:
+                    steps = int(b - 0.414)
+                    for x_ in xrange(x, min(x + steps, int(h_x+1))):
+                        image.pixels[y][x_].draw(color)
+                    x += steps
+                    total_pixels += min(x + steps, int(h_x+1)) - x
+                    continue
+                elif b < -1.414:
+                    steps = int(-b - 0.414)
+                    x += steps
+                    continue
+                total_pixels += 1
+                coverage = 0
+                for j in jitter:
+                    if self.contains(corner + j):
+                        coverage += 1.0
+                image.pixels[y][x].draw(color.fainter(coverage / lj))
+                x += 1
diff --git a/rasterizer/test_rasterizer.py b/rasterizer/test_rasterizer.py
new file mode 100644
index 000000000..a425113d1
--- /dev/null
+++ b/rasterizer/test_rasterizer.py
@@ -0,0 +1,69 @@
+from rasterizer import *
+import sys
+import cProfile
+import random
+import math
+
+def do_it():
+    f = open(sys.argv[1], 'w')
+    i = PPMImage(512, Color(1,1,1,1))
+    s = Scene()
+    s3 = Scene()
+    s3.add(Union(Circle(Vector(0.3, 0.1), 0.1),
+                 Circle(Vector(0.35, 0.1), 0.1),
+                 Color(0,0,0,0.5)))
+    s3.add(Intersection(Circle(Vector(0.3, 0.3), 0.1),
+                        Circle(Vector(0.35, 0.3), 0.1),
+                        Color(0,0.5,0,1)))
+    s3.add(Subtraction(Circle(Vector(0.3, 0.5), 0.1),
+                       Circle(Vector(0.35, 0.5), 0.1),
+                       Color(0,0,0.5,1)))
+    s3.add(Subtraction(Circle(Vector(0.35, 0.7), 0.1),
+                       Circle(Vector(0.3, 0.7), 0.1),
+                       Color(0,0.5,0.5,1)))
+    s2 = Scene([
+        LineSegment(Vector(0.0,0), Vector(0.0,1), 0.01, Color(1,0,0,1)),
+        LineSegment(Vector(0.1,0), Vector(0.1,1), 0.01, Color(1,0,0,1)),
+        LineSegment(Vector(0.2,0), Vector(0.2,1), 0.01, Color(1,0,0,1)),
+        LineSegment(Vector(0.3,0), Vector(0.3,1), 0.01, Color(1,0,0,1)),
+        LineSegment(Vector(0.4,0), Vector(0.4,1), 0.01, Color(1,0,0,1)),
+        LineSegment(Vector(0.5,0), Vector(0.5,1), 0.01, Color(1,0,0,1)),
+        LineSegment(Vector(0.6,0), Vector(0.6,1), 0.01, Color(1,0,0,1)),
+        LineSegment(Vector(0.2,0.1), Vector(0.7, 0.4), 0.01, Color(0,0.5,0,1))])
+    for shape in [
+        Poly([Vector(0.2,0.1), Vector(0.9,0.3), Vector(0.1,0.4)],
+             Color(1,0,0,0.5)),
+        Poly([Vector(0.5,0.2), Vector(1,0.4), Vector(0.2,0.7)],
+             Color(0,1,0,0.5)),
+        Rectangle(Vector(0.1,0.7), Vector(0.6,0.8),
+                  Color(0,0.5,0.8,0.5)),
+        Poly([Vector(-1, 0.6), Vector(0.2, 0.8), Vector(-2, 0.7)],
+             Color(1,0,1,0.9)),
+        Circle(Vector(0.5,0.9), 0.2, Color(0,0,0,1)),
+        Circle(Vector(0.9, 0.5), 0.2,
+               Color(0,1,1,0.5)).transform(around(Vector(0.9, 0.5), scale(1.0, 0.5))),
+        s2,
+        Scene([s2], around(Vector(0.5, 0.5), rotate(math.radians(90)))),
+        Scene([s3], translate(0.5, 0))
+        ]:
+        s.add(shape)
+    s.draw(i)
+    i.write_ppm(f)
+    f.close()
+
+def test_ellipse():
+    Ellipse().transform(
+        scale(1.5, 1)).transform(
+        translate(0, 2)).transform(
+        rotate(math.radians(45))).transform(
+        scale(-1, 2))
+
+def test_eigv():
+    v = Transform(2, 0, 0, 0, 1, 0).eigv()
+    assert(abs(v[0].dot(Vector(1,0))) > 0.9999)
+    assert(abs(v[1].dot(Vector(0,1))) > 0.9999)
+
+if __name__ == '__main__':
+    test_ellipse()
+    test_eigv()
+    do_it()

From 4e3902898a7f6592988776e5c3e45dc2e727fe6c Mon Sep 17 00:00:00 2001
From: Carlos Scheidegger <cscheid@research.att.com>
Date: Thu, 20 Feb 2014 17:11:18 -0500
Subject: [PATCH 02/16] added author line to README.md

---
 README.md | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/README.md b/README.md
index 93aa1b760..05e8d27ca 100644
--- a/README.md
+++ b/README.md
@@ -113,4 +113,12 @@ Contributors
     <td>@jaaaarel</td>
     <td>taavi.burns@points.com</td>
   </tr>
+  <tr>
+    <td>Carlos Scheidegger</td>
+    <td>AT&amp;T Research</td>
+    <td>rasterizer</td>
+    <td>@cscheid</td>
+    <td>@cscheid</td>
+    <td>carlos.scheidegger@gmail.com</td>
+  </tr>
 </table>

From 86bf52ee6aaa4637975ae8cdb785ebf5f6eb9342 Mon Sep 17 00:00:00 2001
From: Carlos Scheidegger <cscheid@research.att.com>
Date: Fri, 21 Feb 2014 22:22:48 -0500
Subject: [PATCH 03/16] some cleanup and documentation

---
 rasterizer/README.md              |   2 +
 rasterizer/doc/README.md          | 141 ++++++++++++++++++++++++++++++
 rasterizer/rasterizer/ellipse.md  |  67 ++++++++++++++
 rasterizer/rasterizer/ellipse.py  |  22 ++---
 rasterizer/rasterizer/geometry.py |  22 ++---
 rasterizer/rasterizer/poly.py     |  13 ++-
 rasterizer/rasterizer/shape.py    |  14 ++-
 rasterizer/test_rasterizer.py     |  31 ++++---
 rasterizer/todo.org               |   6 ++
 9 files changed, 272 insertions(+), 46 deletions(-)
 create mode 100644 rasterizer/doc/README.md
 create mode 100644 rasterizer/rasterizer/ellipse.md
 create mode 100644 rasterizer/todo.org

diff --git a/rasterizer/README.md b/rasterizer/README.md
index ab163378f..05f257e52 100644
--- a/rasterizer/README.md
+++ b/rasterizer/README.md
@@ -33,6 +33,8 @@ Maybe most interestingly, shapes in `tiny_gfx` are extensible. New
 shapes are easily added, and they compose well with the other parts of
 the code.
 
+A description of the current version of the code is in `doc/README.md`.
+
 ## A performance caveat
 
 Rasterizers are so central to display technology that their
diff --git a/rasterizer/doc/README.md b/rasterizer/doc/README.md
new file mode 100644
index 000000000..9e7b5ffde
--- /dev/null
+++ b/rasterizer/doc/README.md
@@ -0,0 +1,141 @@
+# Overview
+
+By far, the most important class in this project is `Shape`, in
+`shape.py`.
+
+
+## `shape.py`
+
+`Shape` is the main class of `tiny_gfx`. The class itself contains a
+very basic constructor, and the main code for the scanline rasterizer
+itself, in `Shape.draw`.
+
+Concrete subclasses of `Shape` need to implement two methods:
+`contains` and `signed_distance_bound`. Both methods take a point in
+the plane. `contains` should return true if the shape contains the
+point (no surprise there). `signed_distance_bound` is the main method
+used by `tiny_gfx` to *accelerate* the rendering. The idea is that
+`signed_distance_bound` should return a real value that describes a
+*distance certificate*. If this value (call it `r`) is negative, then the shape is
+promising that not only is the point is not in the shape, but no
+points inside a ball of radius `-r` around the point is in the shape as
+well. If the value is positive, then the shape is
+promising that not only is the point in the shape, but a ball of
+radius `r` is entirely contained in the shape as well. If the shape
+cannot guarantee anything either way, it should return zero. In other
+words, `signed_distance_bound` is used by `Shape.draw` to skip large
+numbers of checks against `contains`, both inside and outside the
+shape. It's ok if `signed_distance_bound` is conservative. In fact,
+always returning zero would be correct behavior. It would just be
+much slower.
+
+Finally, concrete subclasses of `Shape` must have a `self.bound` field
+that stores an `AABox` object that represents an axis-aligned bounding
+box of the object. The bounding box does not need to be exact, but it
+must be *conservative*: `self.contains(p)` must imply
+`self.bounds.contains(p)`.  The smaller the area of this box, the more
+efficient rasterization will be.
+
+`Shape.draw` works by traversing candidate image pixels inside the
+shape bounding box in row order. Anti-aliasing is performed by
+jittered sampling (see Figure 5.16
+[here](http://www.cs.utah.edu/~shirley/papers/thesis/ch5b.pdf)). Note
+that if `Shape.draw` knows that the pixel is entirely inside the image
+(because `r` from above is greater than the length of a pixel
+diagonal), then no anti-aliasing tests are needed. This allows the
+default setting for super-sampling to be to take 36 samples inside a
+pixel with no large loss in performance. By the same token, if `r` is
+a negative number, then `Shape.draw` skips those pixels with no
+further checks. Anti-aliasing is performed by counting the number of
+points that pass the `contains` call, and updating the image pixel.
+
+
+
+## `geometry.py`
+
+This file contains geometry classes generally needed for the
+rasterizer.
+
+* `Vector`, a 2D vector with your basic operator overloading and
+  methods. In this code we use this class to store both points and
+  vectors. There are reasons why this is a bad idea, but for sake of
+  simplicity and brevity, we do it.
+
+* `AABox`, a 2D axis-aligned bounding box
+
+* `HalfPlane`, a class that models one-half of the 2D plane by a
+  linear equation
+  
+* `Transform`, a class for 2D affine transformations of `Vector`s
+
+* utility functions to build transforms (which should perhaps be
+  `classmethod`s of Transform, except that using them leads to long,
+  unreadable lines for things that should have short names)
+  
+
+## `color.py`
+
+* `Color`: Self-contained RGBA class in plain old dumb RGBA.
+
+## `csg.py`
+
+This file contains classes used for Boolean operations with
+shapes. (CSG stands for Constructive Solid Geometry, the three-letter
+acronym used in graphics for the idea). The base class is `CSG`,
+and there's a class for each supported Boolean operation: `Union`,
+`Intersection` and `Subtraction`.
+
+If the `signed_distance_bound` and `contains` ideas above makes sense,
+then the code for the three classes should be self-explanatory.
+
+
+## `ellipse.py`
+
+`Ellipse` is the most complicated `Shape` available, and is presented
+as an example of the generality of the approach used here. The
+internal representation for an ellipse is the implicit equation form,
+which defines the ellipse as the set of points for which a certain
+function is less than 0. In this case, the function is a quadratic
+polynomial in x and y.
+
+`Ellipse.contains` simply evaluates the implicit equation.
+
+The code for `Ellipse.signed_distance_bound` is actually relatively
+clever and non-trivial. This flavor of geometric insights is prevalent
+in graphics, so I wanted to give an actual example of it in the
+code. However, it takes a diagram to show why it works, so I don't
+really expect you to understand it without comments.
+
+I have to draw a diagram to explain how it
+works, so if you run into trouble on that one, send me a message and
+I'll move it to the top of my priority queue.
+
+There's more description of the code for `Ellipse` under `rasterizer/ellipse.md`.
+
+## `image.py`
+
+`PPMImage`: Simple, self-contained class that stores a two-dimensional array of
+`Color` pixels, and writes them as a
+[PPM file](http://netpbm.sourceforge.net/doc/ppm.html).
+
+
+## `poly.py`
+
+`ConvexPoly` represents a convex polygon by a set of
+half-planes. `Shape.contains` simply tests all half-planes
+corresponding to each edge of the polygon, and
+`Shape.signed_distance_bound` takes the most conservative value across
+all of the shape half-planes. This actually gives values of 0 for
+points on the "infinite line" spanned by polygon edges, but that's
+fine because the result needs only be conservative.
+
+
+## `scene.py`
+
+`Scene` stores a hierarchical scene graph, where each node is either a
+`Shape` or a `Scene` itself. In addition, each `Scene` object carries
+an affine transform that's applied to every element under it. By
+having different scenes holding the same object lists with different
+transformations, it becomes easy to express scenes with repeated
+elements.
+
diff --git a/rasterizer/rasterizer/ellipse.md b/rasterizer/rasterizer/ellipse.md
new file mode 100644
index 000000000..d535a110f
--- /dev/null
+++ b/rasterizer/rasterizer/ellipse.md
@@ -0,0 +1,67 @@
+# The ellipse class
+
+The ellipse class is the most complicated shape of the system. I
+include it to show how complex shapes can be incorporated in the
+tiny_gfx. This section will require more mathematical background than
+the other sections, but not more than calculus.
+
+We use an *implicit equation* for the ellipse. This is a function of
+the input coordinates x and y, which will return a positive value if
+the point lies outside the shape. Any ellipse can be represented by
+the following polynomial in x and y:
+
+v(x, y) = a x^2 + b y^2 + c x y + d x + e y + f
+
+Most of the class will 
+
+## Finding the center of the ellipse
+
+We will need the center of the ellipse to compute the signed distance
+bound below. By a happy accident, the gradient of the implicit ellipse
+function is an affine transformation of the input vector, and the
+center of the ellipse is the only point in the function where the
+gradient is zero. To find the input vector whose gradient is zero,
+then, we invert the transformation and multiply by Vector(0, 0).
+
+In a real system, you wouldn't solve a system Ax=b by explicitly
+finding the inverse: it is slower than necessary and numerically
+unstable. We use the inverse here for convenience (we already need the
+inverse of a transformation in order to transform the ellipse, see
+below) and to fit in the 500 line budget.
+
+## Finding the bounding box of the ellipse
+
+Our shape class needs a bounding box. Take the partial with respect to
+x to find where the ellipse boundary is vertical, df/dx(x, y) = 0, and
+v(x, y) = 0, so that the point is on the boundary. The solutions will
+give the leftmost and rightmost points of the ellipse. solving the
+system of equations df/dy(x, y) and v(x, y) = 0 gives the top and
+bottom.
+
+## Transforming an ellipse
+
+Say we're given an affine transformation T(x, y) = M . (x, y) = (x', y') that we
+need to apply to the ellipse to get a new ellipse. This means that the
+implicit equation for the new ellipse needs to be written as a
+function of x' and y', but our current equation is a polynomial in x
+and y. To get the new equation, we need to write x and y themselves as
+functions of x' and y'
+
+T(x, y) = (x', y')
+
+## Bounding the distance to an ellipse
+
+The idea is that the distance from a point to an ellipse is not
+exactly trivial to compute, and requires numerical optimization. We
+don't want to waste lines of code on a numerical optimizer and an
+ugly(ier?) function. Instead, we can get a good bound quickly by:
+
+* computing the intersection between the point `p` and the ellipse
+center (call this intersection `i`)
+* finding a line that goes through `i` in the direction of the normal
+ of the ellipse (call this line `n`)
+* project the `p` onto the `l` (call it `s`)
+* return the distance between `l` and `s`
+
+The proof that this is a lower bound of the distance to an ellipse is
+easy, but it takes a diagram to make it obvious. I'll make one soon.
diff --git a/rasterizer/rasterizer/ellipse.py b/rasterizer/rasterizer/ellipse.py
index 4037c1327..3f05a748c 100644
--- a/rasterizer/rasterizer/ellipse.py
+++ b/rasterizer/rasterizer/ellipse.py
@@ -1,25 +1,27 @@
 from shape import Shape
-from geometry import Vector, Transform, quadratic, scale, translate
+from geometry import Vector, Transform, quadratic, scale, translate, AABox
 from color import Color
+import sys
 
 class Ellipse(Shape):
     def __init__(self, a=1.0, b=1.0, c=0.0, d=0.0, e=0.0, f=-1.0, color=None):
         Shape.__init__(self, color)
+        if c*c - 4*a*b >= 0:
+            raise Exception("Not an ellipse")
         self.a = a
         self.b = b
         self.c = c
         self.d = d
         self.e = e
         self.f = f
-        t = Transform(2 * a, c, 0, c, 2 * b, 0)
-        self.center = t.inverse() * Vector(-d, -e)
-        l1, l0 = quadratic(1, 2 * (-a - b), 4 * a * b - c * c)
-        v = t.eigv()
-        axes = [v[0] * ((l0 / 2) ** -0.5), v[1] * ((l1 / 2) ** -0.5)]
-        self.bound = Vector.union(self.center - axes[0] - axes[1],
-                                  self.center - axes[0] + axes[1],
-                                  self.center + axes[0] - axes[1],
-                                  self.center + axes[0] + axes[1])
+        g = Transform(2*a, c, d, c, 2*b, e)
+        self.center = g.inverse() * Vector(0, 0)
+        y1, y2 = quadratic(b-c*c/4*a, e-c*d/2*a, f-d*d/4*a)
+        x1, x2 = quadratic(a-c*c/4*b, d-c*e/2*b, f-e*e/4*b)
+        self.bound = Vector.union(Vector(-(d + c*y1)/2*a, y1),
+                                  Vector(-(d + c*y2)/2*a, y2),
+                                  Vector(x1, -(e + c*x1)/2*b),
+                                  Vector(x2, -(e + c*x2)/2*b))
     def value(self, p):
         return self.a*p.x*p.x + self.b*p.y*p.y + self.c*p.x*p.y \
                + self.d*p.x + self.e*p.y + self.f
diff --git a/rasterizer/rasterizer/geometry.py b/rasterizer/rasterizer/geometry.py
index 35f7d5c49..7952b0043 100644
--- a/rasterizer/rasterizer/geometry.py
+++ b/rasterizer/rasterizer/geometry.py
@@ -28,6 +28,8 @@ def union(*args):
         return AABox(reduce(Vector.min, args), reduce(Vector.max, args))
     def length(self):
         return (self.x * self.x + self.y * self.y) ** 0.5
+    def __repr__(self):
+        return "[%.3f %.3f]" % (self.x, self.y)
 
 class AABox:
     def __init__(self, p1, p2):
@@ -62,7 +64,7 @@ def __init__(self, m11, m12, tx, m21, m22, ty):
                   [0, 0, 1]]
     def __mul__(self, other): # ugly
         if isinstance(other, Transform):
-            t = [[0] * 3 for i in xrange(3)]
+            t = [[0.0] * 3 for i in xrange(3)]
             for i, j, k in product(xrange(3), repeat=3):
                 t[i][j] += self.m[i][k] * other.m[k][j]
             return Transform(t[0][0], t[0][1], t[0][2],
@@ -75,18 +77,12 @@ def det(s):
         return s.m[0][0] * s.m[1][1] - s.m[0][1] * s.m[1][0]
     def inverse(self):
         d = 1.0 / self.det()
-        return Transform(d * self.m[1][1], -d * self.m[0][1], -self.m[0][2],
-                         -d * self.m[1][0], d * self.m[0][0], -self.m[1][2])
-    def eigv(self): # power iteration, ignores translation, assumes SPD
-        a = Vector(random.random(), random.random())
-        last = Vector(0,0)
-        t = Transform(self.m[0][0], self.m[0][1], 0,
-                      self.m[1][0], self.m[1][1], 0)
-        while (a - last).length() > 1e-6:
-            last = a
-            a = t * a
-            a = a * (1.0 / a.length())
-        return a, Vector(-a.y, a.x)
+        t = Transform(d * self.m[1][1], -d * self.m[0][1], 0,
+                      -d * self.m[1][0], d * self.m[0][0], 0)
+        v = t * Vector(self.m[0][2], self.m[1][2])
+        t.m[0][2] = -v.x
+        t.m[1][2] = -v.y
+        return t
 
 def identity():
     return Transform(1, 0, 0, 0, 1, 0)
diff --git a/rasterizer/rasterizer/poly.py b/rasterizer/rasterizer/poly.py
index 5a409cbd8..5c62b0088 100644
--- a/rasterizer/rasterizer/poly.py
+++ b/rasterizer/rasterizer/poly.py
@@ -1,7 +1,7 @@
 from shape import Shape
 from geometry import HalfPlane, Vector
 
-class Poly(Shape): # a *convex* poly, in ccw order, with no repeating vertices
+class ConvexPoly(Shape): # a *convex* poly, in ccw order, with no repeating vertices
     def __init__(self, ps, color=None):
         Shape.__init__(self, color)
         mn = min(enumerate(ps), key=lambda (i,v): (v.y, v.x))[0]
@@ -28,9 +28,9 @@ def contains(self, p):
                 return False
         return True
     def transform(self, xform):
-        return Poly(list(xform * v for v in self.vs), color=self.color)
+        return ConvexPoly(list(xform * v for v in self.vs), color=self.color)
 
-Triangle, Quad = Poly, Poly
+Triangle, Quad = ConvexPoly, ConvexPoly
 
 def Rectangle(v1, v2, color=None):
     return Quad([Vector(min(v1.x, v2.x), min(v1.y, v2.y)),
@@ -44,3 +44,10 @@ def LineSegment(v1, v2, thickness, color=None):
     d.x, d.y = -d.y, d.x
     d *= thickness / d.length() / 2
     return Quad([v1 + d, v1 - d, v2 - d, v2 + d], color=color)
+
+# A simple polygon
+# It will work by triangulating a polygon (via simple, slow ear-clipping)
+# and then rendering it by a CSG union of the triangles
+class Poly(Shape):
+    def __init__(self, *args, **kwargs):
+        raise Exception("Unimplemented")
diff --git a/rasterizer/rasterizer/shape.py b/rasterizer/rasterizer/shape.py
index 505d45b19..54160360c 100644
--- a/rasterizer/rasterizer/shape.py
+++ b/rasterizer/rasterizer/shape.py
@@ -7,6 +7,8 @@ class Shape:
     def __init__(self, color=None):
         self.color = color if color is not None else Color()
     def draw(self, image, super_sampling = 6):
+        if not self.bound.overlaps(image.bounds()):
+            return
         color = self.color
         r = float(image.resolution)
         jitter = [Vector((x + random.random()) / super_sampling / r,
@@ -14,14 +16,10 @@ def draw(self, image, super_sampling = 6):
                   for (x, y) in product(xrange(super_sampling), repeat=2)]
         lj = len(jitter)
         total_pixels = 0
-        tb = self.bound
-        if not tb.overlaps(image.bounds()):
-            return
-        l_x = max(0, int(tb.low.x * r))
-        l_y = max(0, int(tb.low.y * r))
-        h_x = min(r-1, int(tb.high.x * r))
-        h_y = min(r-1, int(tb.high.y * r))
-        corners = [(0, 0), (1.0/r, 0), (0, 1.0/r), (1.0/r, 1.0/r)]
+        l_x = max(0, int(self.bound.low.x * r))
+        l_y = max(0, int(self.bound.low.y * r))
+        h_x = min(r-1, int(self.bound.high.x * r))
+        h_y = min(r-1, int(self.bound.high.y * r))
         for y in xrange(l_y, int(h_y+1)):
             x = l_x
             while x <= h_x:
diff --git a/rasterizer/test_rasterizer.py b/rasterizer/test_rasterizer.py
index a425113d1..2aa16b8fe 100644
--- a/rasterizer/test_rasterizer.py
+++ b/rasterizer/test_rasterizer.py
@@ -31,14 +31,14 @@ def do_it():
         LineSegment(Vector(0.6,0), Vector(0.6,1), 0.01, Color(1,0,0,1)),
         LineSegment(Vector(0.2,0.1), Vector(0.7, 0.4), 0.01, Color(0,0.5,0,1))])
     for shape in [
-        Poly([Vector(0.2,0.1), Vector(0.9,0.3), Vector(0.1,0.4)],
-             Color(1,0,0,0.5)),
-        Poly([Vector(0.5,0.2), Vector(1,0.4), Vector(0.2,0.7)],
-             Color(0,1,0,0.5)),
+        Triangle([Vector(0.2,0.1), Vector(0.9,0.3), Vector(0.1,0.4)],
+                 Color(1,0,0,0.5)),
+        Triangle([Vector(0.5,0.2), Vector(1,0.4), Vector(0.2,0.7)],
+                 Color(0,1,0,0.5)),
         Rectangle(Vector(0.1,0.7), Vector(0.6,0.8),
                   Color(0,0.5,0.8,0.5)),
-        Poly([Vector(-1, 0.6), Vector(0.2, 0.8), Vector(-2, 0.7)],
-             Color(1,0,1,0.9)),
+        Triangle([Vector(-1, 0.6), Vector(0.2, 0.8), Vector(-2, 0.7)],
+                 Color(1,0,1,0.9)),
         Circle(Vector(0.5,0.9), 0.2, Color(0,0,0,1)),
         Circle(Vector(0.9, 0.5), 0.2,
                Color(0,1,1,0.5)).transform(around(Vector(0.9, 0.5), scale(1.0, 0.5))),
@@ -58,12 +58,19 @@ def test_ellipse():
         rotate(math.radians(45))).transform(
         scale(-1, 2))
 
-def test_eigv():
-    v = Transform(2, 0, 0, 0, 1, 0).eigv()
-    assert(abs(v[0].dot(Vector(1,0))) > 0.9999)
-    assert(abs(v[1].dot(Vector(0,1))) > 0.9999)
+def test_inverse():
+    for i in xrange(10000):
+        f = Transform(random.random(), random.random(), random.random(),
+                      random.random(), random.random(), random.random())
+        v = Vector(random.random(), random.random())
+        m = f * f.inverse()
+        if (m * v - v).length() > 1e-4:
+            print >>sys.stderr, "inverse failed!"
+            print f.m
+            print m.m
+            raise Exception("foo")
 
 if __name__ == '__main__':
-    test_ellipse()
-    test_eigv()
+    # test_inverse()
+    # test_ellipse()
     do_it()
diff --git a/rasterizer/todo.org b/rasterizer/todo.org
new file mode 100644
index 000000000..e7c0d48e5
--- /dev/null
+++ b/rasterizer/todo.org
@@ -0,0 +1,6 @@
+* Refactoring
+** TODO Vector.union becomes AABox.union
+
+* Comments
+** Describe Ellipse.signed_distance_bound
+** Explain why Points and Vectors should be different (but aren't)

From 4691fa15a627eb3f815c84d291d6b790fd37f1e4 Mon Sep 17 00:00:00 2001
From: Carlos Scheidegger <cscheid@research.att.com>
Date: Sun, 2 Mar 2014 22:52:07 -0500
Subject: [PATCH 04/16] make changes suggested by @EkkiD, add some other
 comments, small cosmetic changes

---
 rasterizer/rasterizer/geometry.py | 17 +++++++++++++-
 rasterizer/rasterizer/poly.py     |  3 +--
 rasterizer/rasterizer/shape.py    | 37 ++++++++++++++-----------------
 3 files changed, 34 insertions(+), 23 deletions(-)

diff --git a/rasterizer/rasterizer/geometry.py b/rasterizer/rasterizer/geometry.py
index 7952b0043..e53a2f6f7 100644
--- a/rasterizer/rasterizer/geometry.py
+++ b/rasterizer/rasterizer/geometry.py
@@ -2,6 +2,9 @@
 import math
 from itertools import product
 
+# Solves the quadratic equation ax^2 + bx + c = 0
+# using a variant of the standard quadratic formula
+# that behaves better numerically
 def quadratic(a, b, c):
     d = (b * b - 4 * a * c) ** 0.5
     if b >= 0:
@@ -57,19 +60,28 @@ def __init__(self, p1, p2):
     def signed_distance(self, p):
         return self.v.dot(p) + self.c
 
+# Transform represents an affine transformation
+# of a 2D vector ("affine" because it's a linear transformation
+# combined with a translation)
 class Transform:
     def __init__(self, m11, m12, tx, m21, m22, ty):
         self.m = [[m11, m12, tx],
                   [m21, m22, ty],
                   [0, 0, 1]]
-    def __mul__(self, other): # ugly
+    def __mul__(self, other):
         if isinstance(other, Transform):
+            # if the other element is also a transform,
+            # then return a transform corresponding to the
+            # composition of the two transforms
             t = [[0.0] * 3 for i in xrange(3)]
             for i, j, k in product(xrange(3), repeat=3):
                 t[i][j] += self.m[i][k] * other.m[k][j]
             return Transform(t[0][0], t[0][1], t[0][2],
                              t[1][0], t[1][1], t[1][2])
         else:
+            # if the other element is a vector, then
+            # apply the transformation to the vector via
+            # a matrix-vector multiplication
             nx = self.m[0][0] * other.x + self.m[0][1] * other.y + self.m[0][2]
             ny = self.m[1][0] * other.x + self.m[1][1] * other.y + self.m[1][2]
             return Vector(nx, ny)
@@ -98,5 +110,8 @@ def translate(tx, ty):
 def scale(x, y):
     return Transform(x, 0, 0, 0, y, 0)
 
+# To perform a transformation 'around' some point,
+# we first translate that point to the origin, perform a transformation,
+# then translate back that same amount
 def around(v, t):
     return translate(v.x, v.y) * t * translate(-v.x, -v.y)
diff --git a/rasterizer/rasterizer/poly.py b/rasterizer/rasterizer/poly.py
index 5c62b0088..b0000b095 100644
--- a/rasterizer/rasterizer/poly.py
+++ b/rasterizer/rasterizer/poly.py
@@ -4,8 +4,7 @@
 class ConvexPoly(Shape): # a *convex* poly, in ccw order, with no repeating vertices
     def __init__(self, ps, color=None):
         Shape.__init__(self, color)
-        mn = min(enumerate(ps), key=lambda (i,v): (v.y, v.x))[0]
-        self.vs = list(ps[mn:]) + list(ps[:mn])
+        self.vs = ps
         self.bound = Vector.union(*self.vs)
         self.half_planes = []
         for i in xrange(len(self.vs)):
diff --git a/rasterizer/rasterizer/shape.py b/rasterizer/rasterizer/shape.py
index 54160360c..7868c5808 100644
--- a/rasterizer/rasterizer/shape.py
+++ b/rasterizer/rasterizer/shape.py
@@ -15,31 +15,28 @@ def draw(self, image, super_sampling = 6):
                          (y + random.random()) / super_sampling / r)
                   for (x, y) in product(xrange(super_sampling), repeat=2)]
         lj = len(jitter)
-        total_pixels = 0
-        l_x = max(0, int(self.bound.low.x * r))
-        l_y = max(0, int(self.bound.low.y * r))
-        h_x = min(r-1, int(self.bound.high.x * r))
-        h_y = min(r-1, int(self.bound.high.y * r))
+        l_x = max(int(self.bound.low.x * r), 0)
+        l_y = max(int(self.bound.low.y * r), 0)
+        h_x = min(int(self.bound.high.x * r), r-1)
+        h_y = min(int(self.bound.high.y * r), r-1)
         for y in xrange(l_y, int(h_y+1)):
             x = l_x
             while x <= h_x:
                 corner = Vector(x / r, y / r)
-                b = self.signed_distance_bound(corner) * r
-                if b > 1.414:
-                    steps = int(b - 0.414)
+                b = self.signed_distance_bound(corner)
+                pixel_diameter = (2 ** 0.5) / r
+                if b > pixel_diameter:
+                    steps = int(r * (b - (pixel_diameter - 1.0/r)))
                     for x_ in xrange(x, min(x + steps, int(h_x+1))):
                         image.pixels[y][x_].draw(color)
                     x += steps
-                    total_pixels += min(x + steps, int(h_x+1)) - x
-                    continue
-                elif b < -1.414:
-                    steps = int(-b - 0.414)
+                elif b < -pixel_diameter:
+                    steps = int(r * (-b - (pixel_diameter - 1.0/r)))
                     x += steps
-                    continue
-                total_pixels += 1
-                coverage = 0
-                for j in jitter:
-                    if self.contains(corner + j):
-                        coverage += 1.0
-                image.pixels[y][x].draw(color.fainter(coverage / lj))
-                x += 1
+                else:
+                    coverage = 0
+                    for j in jitter:
+                        if self.contains(corner + j):
+                            coverage += 1.0
+                    image.pixels[y][x].draw(color.fainter(coverage / lj))
+                    x += 1

From 4df3da8f1f71b2ef474a1438df9d5209f719e659 Mon Sep 17 00:00:00 2001
From: Carlos Scheidegger <cscheid@research.att.com>
Date: Mon, 3 Mar 2014 00:00:24 -0500
Subject: [PATCH 05/16] added diagram and explanation of one of the ellipse
 cases

---
 rasterizer/doc/ellipse_1.svg     | 55 ++++++++++++++++++++++++++++++++
 rasterizer/rasterizer/ellipse.md | 46 +++++++++++++++++++-------
 rasterizer/rasterizer/ellipse.py | 11 ++++---
 3 files changed, 96 insertions(+), 16 deletions(-)
 create mode 100644 rasterizer/doc/ellipse_1.svg

diff --git a/rasterizer/doc/ellipse_1.svg b/rasterizer/doc/ellipse_1.svg
new file mode 100644
index 000000000..4a9d368f3
--- /dev/null
+++ b/rasterizer/doc/ellipse_1.svg
@@ -0,0 +1,55 @@
+<!-- Made in Mondrian.io -->
+<svg id="main" width="1000" height="800" version="1.1" xmlns="http://www.w3.org/2000/svg" xmlns:mondrian="http://mondrian.io/xml" viewbox="0 0 1000 800" enable-background="new 0 0 1000 800" style="width: 1000px; height: 800px; -webkit-transform: scale(1.6006652139172843);">
+<ellipse cx="516.3178322529608" cy="376.24088826133703" rx="164.99999999999986" ry="85.4999999999998"  fill="rgba(255, 255, 255, 1)" stroke="rgba(0, 0, 0, 1)" stroke-width="1">
+</ellipse>
+<path stroke="rgba(0, 0, 0, 1)" fill="rgba(255, 255, 255, 1)" stroke-width="1" d="M515.4505155384234,375.04931161960553 L817.4505155384212,113.04931161960516" >
+</path>
+<ellipse cx="515.2815484240302" cy="375.29280959954326" rx="3.559151217920075" ry="3.5591512179200713"  fill="rgba(0, 0, 0, 1)" stroke="rgba(0, 0, 0, 1)" stroke-width="1">
+</ellipse>
+<path fill="rgba(0, 0, 0, 1)" stroke="rgba(0, 0, 0, 1)" stroke-width="1" d="M816.0798056290876,109.22082508844622 C818.0454536637212,109.22082508844652 819.6389568470105,110.82244313651044 819.6389568470078,112.77997630636705 S818.0373387989447,116.33912752428645 816.0798056290887,116.33912752428645 S812.5206544111683,114.73750947622293 812.5206544111683,112.77997630636676 S815.0137215764473,110.86259270409474 816.0798056290887,109.2208250884465z" >
+</path>
+<path fill="rgba(0, 0, 0, 1)" stroke="rgba(0, 0, 0, 1)" stroke-width="1" d="M950.5460781313685,313.91121703294397 C952.5117261660018,313.9112170329449 954.1052293492888,315.5128350810089 954.1052293492855,317.4703682508654 S952.5036113012249,321.02951946878466 950.5460781313689,321.02951946878466 S946.9869269134484,319.4279014207211 946.9869269134484,317.4703682508649 S948.5885449615126,313.9112170329447 950.5460781313689,313.9112170329447z" >
+</path>
+<path fill="rgba(0, 0, 0, 1)" stroke="rgba(0, 0, 0, 1)" stroke-width="1" d="M599.5460781313666,298.91121703294283 C601.5117261660029,298.91121703294505 603.1052293492904,300.51283508100886 603.1052293492864,302.47036825086553 S601.5036113012256,306.0295194687849 599.5460781313698,306.0295194687849 S595.9869269134495,304.4279014207214 595.9869269134495,302.47036825086514 S597.5885449615134,298.911217032945 599.5460781313698,298.911217032945z" >
+</path>
+<path stroke="rgba(0, 0, 0, 1)" fill="rgba(255, 255, 255, 1)" stroke-width="1" d="M598.8491828110405,302.82990763596547 L693.626415809849,34.596601283944736" >
+</path>
+<path stroke="rgba(0, 0, 0, 1)" fill="rgba(255, 255, 255, 1)" stroke-width="1" d="M816.2372618649003,113.84909411145918 L682.6730232777935,68.2268977007501" >
+</path>
+<path stroke="rgba(0, 0, 0, 1)" fill="rgba(0, 0, 0, 1)" stroke-width="1" d="M815.8884371874112,114.4363002808318 L656.2427343264982,329.8873594515336" >
+</path>
+<path stroke="rgba(0, 0, 0, 1)" fill="rgba(0, 0, 0, 1)" stroke-width="1" d="M602.556568762651,303.0442766696096 L757.2574932163678,358.1432360640845" >
+</path>
+<path stroke="rgba(0, 0, 0, 1)" fill="rgba(0, 0, 0, 1)" stroke-width="1" d="M814.8884371874109,115.26188794208673 L728.7080135191308,346.96007616500555" >
+</path>
+<path stroke="rgba(0, 0, 0, 1)" fill="rgba(0, 0, 0, 1)" stroke-width="1" d="M585.9571176757494,298.45357314304226 L590.9018960829482,282.20644409082666" >
+</path>
+<path stroke="rgba(0, 0, 0, 1)" fill="rgba(0, 0, 0, 1)" stroke-width="1" d="M591.0771806773158,282.71126372261324 L604.7511318841725,286.5710304906553" >
+</path>
+<path stroke="rgba(0, 0, 0, 1)" fill="rgba(0, 0, 0, 1)" stroke-width="1" d="M678.0307646199999,79.15214321393834 L690.8516194558027,83.63944240646939" >
+</path>
+<path stroke="rgba(0, 0, 0, 1)" fill="rgba(0, 0, 0, 1)" stroke-width="1" d="M690.3306826209499,83.45566233965721 L694.3853138056319,72.36547953525258" >
+</path>
+<path stroke="rgba(0, 0, 0, 1)" fill="rgba(0, 0, 0, 1)" stroke-width="1" d="M733.481282670545,336.04661940533634 L720.3723459991867,331.43057332903004" >
+</path>
+<path stroke="rgba(0, 0, 0, 1)" fill="rgba(0, 0, 0, 1)" stroke-width="1" d="M720.1765748278801,332.05815529964013 L716.0963995149973,343.17683425854875" >
+</path>
+<text x="828.2093848336237" y="129.6557822902478" font-size="24" font-family="Arial" fill="rgba(0, 0, 0, 1)" stroke="rgba(0, 0, 0, 1)"  style="display: block;" stroke-width="1" transform="rotate(0 834.709 117.656) scale(1 1) translate(0 0)">
+p</text>
+<text x="597.475013084448" y="324.93326663760064" font-size="24" font-family="Arial" fill="rgba(0, 0, 0, 1)" stroke="rgba(0, 0, 0, 1)"  style="display: block;" stroke-width="1" transform="rotate(0 599.975 312.933) scale(1 1) translate(0 0)">
+i</text>
+<text x="620.7601331692273" y="179.93047338238472" font-size="24" font-family="Arial" fill="rgba(0, 0, 0, 1)" stroke="rgba(0, 0, 0, 1)"  style="display: block;" stroke-width="1" transform="rotate(0 627.26 167.93) scale(1 1) translate(0 0)">
+n</text>
+<text x="640.6578186977242" y="345.9801241250146" font-size="24" font-family="Arial" fill="rgba(0, 0, 0, 1)" stroke="rgba(0, 0, 0, 1)"  style="display: block;" stroke-width="1" transform="rotate(0 647.158 333.98) scale(1 1) translate(0 0)">
+e</text>
+<path fill="rgba(0, 0, 0, 1)" stroke="rgba(0, 0, 0, 1)" stroke-width="1" d="M682.3922642189752,65.24754040408655 C684.3579122536091,65.2475404040874 685.9514154368969,66.84915845215116 685.9514154368945,68.80669162200795 S684.3497973888324,72.36584283992732 682.3922642189764,72.36584283992732 S678.833113001056,70.7642247918638 678.833113001056,68.80669162200758 S680.4347310491202,65.24754040408737 682.3922642189764,65.24754040408737z" >
+</path>
+<text x="663.6989971065256" y="65.73629684945303" font-size="24" font-family="Arial" fill="rgba(0, 0, 0, 1)" stroke="rgba(0, 0, 0, 1)"  style="display: block;" stroke-width="1" transform="rotate(0 669.699 53.736) scale(1 1) translate(0 0)">
+s</text>
+<text x="686.6671008250099" y="353.23359511874764" font-size="24" font-family="Arial" fill="rgba(0, 0, 0, 1)" stroke="rgba(0, 0, 0, 1)"  style="display: block;" stroke-width="1" transform="rotate(0 690.167 341.234) scale(1 1) translate(0 0)">
+t</text>
+<text x="718.3472438849867" y="368.2816630722368" font-size="24" font-family="Arial" fill="rgba(0, 0, 0, 1)" stroke="rgba(0, 0, 0, 1)"  style="display: block;" stroke-width="1" transform="rotate(0 726.847 356.282) scale(1 1) translate(0 0)">
+s'</text>
+<text x="651.6040899317613" y="317.6175319896437" font-size="24" font-family="Arial" fill="rgba(0, 0, 0, 1)" stroke="rgba(0, 0, 0, 1)"  style="display: block;" stroke-width="1" transform="rotate(0 660.604 305.618) scale(1 1) translate(0 0)">
+e'</text>
+</svg>
diff --git a/rasterizer/rasterizer/ellipse.md b/rasterizer/rasterizer/ellipse.md
index d535a110f..4508464cd 100644
--- a/rasterizer/rasterizer/ellipse.md
+++ b/rasterizer/rasterizer/ellipse.md
@@ -12,7 +12,11 @@ the following polynomial in x and y:
 
 v(x, y) = a x^2 + b y^2 + c x y + d x + e y + f
 
-Most of the class will 
+The boundary of the ellipse is exactly the set of points where v(x,y)
+= 0. When we evaluate a point x,y and get a value greater than zero,
+we consider that point to be outside the ellipse. Most of the ellipse
+class will involve manipulating this implicit equation.
+
 
 ## Finding the center of the ellipse
 
@@ -29,14 +33,18 @@ unstable. We use the inverse here for convenience (we already need the
 inverse of a transformation in order to transform the ellipse, see
 below) and to fit in the 500 line budget.
 
+
 ## Finding the bounding box of the ellipse
 
-Our shape class needs a bounding box. Take the partial with respect to
-x to find where the ellipse boundary is vertical, df/dx(x, y) = 0, and
-v(x, y) = 0, so that the point is on the boundary. The solutions will
-give the leftmost and rightmost points of the ellipse. solving the
-system of equations df/dy(x, y) and v(x, y) = 0 gives the top and
-bottom.
+Our shape class needs a bounding box. The partial derivative measures
+the rate of change of the function with respect a variable.  When the
+implicit equation of the ellipse is zero, *and* the rate of change of
+the implicit function with respect to x is zero, then ellipse boundary
+is vertical. So we solve the system of equations df/dx(x, y) = 0, v(x,
+y) = 0, and the two solutions are going to be touching the top and
+bottom of the bounding box exactly. Solving the system of equations
+df/dy(x, y) = 0, v(x, y) = 0 gives the left and right points.
+
 
 ## Transforming an ellipse
 
@@ -49,6 +57,10 @@ functions of x' and y'
 
 T(x, y) = (x', y')
 
+This gives 
+
+
+
 ## Bounding the distance to an ellipse
 
 The idea is that the distance from a point to an ellipse is not
@@ -60,8 +72,18 @@ ugly(ier?) function. Instead, we can get a good bound quickly by:
 center (call this intersection `i`)
 * finding a line that goes through `i` in the direction of the normal
  of the ellipse (call this line `n`)
-* project the `p` onto the `l` (call it `s`)
-* return the distance between `l` and `s`
-
-The proof that this is a lower bound of the distance to an ellipse is
-easy, but it takes a diagram to make it obvious. I'll make one soon.
+* project the `p` onto the `i` (call it `s`)
+* return the distance between `i` and `s`
+
+Referring to the [diagram](../doc/ellipse_1.svg), `p` is outside the ellipse,
+`t` is the line tangent to `i`, and `e` is the closest point in the
+ellipse to `p`. Notice that the distance between `s` and `i` is equal
+to the distance between `p` and `s'`, since they form parallel
+lines. Now take the line from `p` to `e`. Notice that it intersects
+the line from `i` to `s'`, at a point we will call `e'`. Now notice
+that `e'`, `s'` and `p` form a right triangle, of which the line from
+`e'` to `s'` is the hypotenuse (and it's *shorter* than the original
+line too). This shortened line, in turn, is yet shorter than the leg
+`s'` to `p`. But length of this leg is precisely the distance from `s`
+to `i`, and so we know that it's a lower bound to the distance from
+`e` to `p`, and we are done.
diff --git a/rasterizer/rasterizer/ellipse.py b/rasterizer/rasterizer/ellipse.py
index 3f05a748c..4c4bb5b07 100644
--- a/rasterizer/rasterizer/ellipse.py
+++ b/rasterizer/rasterizer/ellipse.py
@@ -14,8 +14,8 @@ def __init__(self, a=1.0, b=1.0, c=0.0, d=0.0, e=0.0, f=-1.0, color=None):
         self.d = d
         self.e = e
         self.f = f
-        g = Transform(2*a, c, d, c, 2*b, e)
-        self.center = g.inverse() * Vector(0, 0)
+        gradient = Transform(2*a, c, d, c, 2*b, e)
+        self.center = gradient.inverse() * Vector(0, 0)
         y1, y2 = quadratic(b-c*c/4*a, e-c*d/2*a, f-d*d/4*a)
         x1, x2 = quadratic(a-c*c/4*b, d-c*e/2*b, f-e*e/4*b)
         self.bound = Vector.union(Vector(-(d + c*y1)/2*a, y1),
@@ -44,9 +44,11 @@ def transform(self, transform):
     def signed_distance_bound(self, p):
         def sgn(x):
             return 0 if x == 0 else x / abs(x)
-        v = -sgn(self.value(p))
+        s = -sgn(self.value(p))
         c = self.center
         pc = p - c
+        # to find where the line from p to c intersects
+        # we solve for the value u such that f(c + (p-c) * u) = 0
         u2 = self.a*pc.x**2 + self.b*pc.y**2 + self.c*pc.x*pc.y
         u1 = 2*self.a*c.x*pc.x + 2*self.b*c.y*pc.y \
              + self.c*c.y*pc.x + self.c*c.x*pc.y + self.d*pc.x \
@@ -55,13 +57,14 @@ def sgn(x):
              + self.d*c.x + self.e*c.y + self.f
         sols = quadratic(u2, u1, u0)
         crossings = c+pc*sols[0], c+pc*sols[1]
+        # the surface point we want is the one closest to p
         if (p - crossings[0]).length() < (p - crossings[1]).length():
             surface_pt = crossings[0]
         else:
             surface_pt = crossings[1]
         d = Vector(2*self.a*surface_pt.x + self.c*surface_pt.y + self.d,
                    2*self.b*surface_pt.y + self.c*surface_pt.x + self.e)
-        return v * abs(d.dot(p - surface_pt) / d.length())
+        return s * abs(d.dot(p - surface_pt) / d.length())
 
 def Circle(center, radius, color=None):
     return Ellipse(color=color).transform(

From a3d24bf341693deb5972d91088adb59241f0baf3 Mon Sep 17 00:00:00 2001
From: Carlos Scheidegger <cscheid@research.att.com>
Date: Mon, 3 Mar 2014 02:01:22 -0500
Subject: [PATCH 06/16] fix incorrect implementation of
 Ellipse.signed_distance_bound for internal points

---
 rasterizer/rasterizer/ellipse.py  | 54 ++++++++++++++++++++-----------
 rasterizer/rasterizer/geometry.py |  2 ++
 rasterizer/test_rasterizer.py     | 32 +++++++++++++-----
 3 files changed, 62 insertions(+), 26 deletions(-)

diff --git a/rasterizer/rasterizer/ellipse.py b/rasterizer/rasterizer/ellipse.py
index 4c4bb5b07..ac4d63df7 100644
--- a/rasterizer/rasterizer/ellipse.py
+++ b/rasterizer/rasterizer/ellipse.py
@@ -1,5 +1,5 @@
 from shape import Shape
-from geometry import Vector, Transform, quadratic, scale, translate, AABox
+from geometry import Vector, Transform, quadratic, scale, translate, AABox, HalfPlane
 from color import Color
 import sys
 
@@ -14,8 +14,8 @@ def __init__(self, a=1.0, b=1.0, c=0.0, d=0.0, e=0.0, f=-1.0, color=None):
         self.d = d
         self.e = e
         self.f = f
-        gradient = Transform(2*a, c, d, c, 2*b, e)
-        self.center = gradient.inverse() * Vector(0, 0)
+        self.gradient = Transform(2*a, c, d, c, 2*b, e)
+        self.center = self.gradient.inverse() * Vector(0, 0)
         y1, y2 = quadratic(b-c*c/4*a, e-c*d/2*a, f-d*d/4*a)
         x1, x2 = quadratic(a-c*c/4*b, d-c*e/2*b, f-e*e/4*b)
         self.bound = Vector.union(Vector(-(d + c*y1)/2*a, y1),
@@ -41,14 +41,10 @@ def transform(self, transform):
         ff = self.a*m02*m02 + self.b*m12*m12 + self.c*m02*m12 \
              + self.d*m02 + self.e*m12 + self.f
         return Ellipse(aa, bb, cc, dd, ee, ff, color=self.color)
-    def signed_distance_bound(self, p):
-        def sgn(x):
-            return 0 if x == 0 else x / abs(x)
-        s = -sgn(self.value(p))
-        c = self.center
+    def intersections(self, c, p):
+        # returns the two intersections of the line through c and p
+        # and the ellipse
         pc = p - c
-        # to find where the line from p to c intersects
-        # we solve for the value u such that f(c + (p-c) * u) = 0
         u2 = self.a*pc.x**2 + self.b*pc.y**2 + self.c*pc.x*pc.y
         u1 = 2*self.a*c.x*pc.x + 2*self.b*c.y*pc.y \
              + self.c*c.y*pc.x + self.c*c.x*pc.y + self.d*pc.x \
@@ -56,15 +52,37 @@ def sgn(x):
         u0 = self.a*c.x**2 + self.b*c.y**2 + self.c*c.x*c.y \
              + self.d*c.x + self.e*c.y + self.f
         sols = quadratic(u2, u1, u0)
-        crossings = c+pc*sols[0], c+pc*sols[1]
-        # the surface point we want is the one closest to p
-        if (p - crossings[0]).length() < (p - crossings[1]).length():
-            surface_pt = crossings[0]
+        return c+pc*sols[0], c+pc*sols[1]
+    def signed_distance_bound(self, p):
+        v = self.value(p)
+        if v == 0:
+            return 0
+        elif v < 0:
+            # if inside the ellipse, create an inscribed quadrilateral
+            # that contains the given point and use the minimum distance
+            # from the point to the quadrilateral as a bound. Since
+            # the quadrilateral lies entirely inside the ellipse, the
+            # distance from the point to the ellipse must be smaller.
+            v0, v2 = self.intersections(p, p + Vector(1, 0))
+            v1, v3 = self.intersections(p, p + Vector(0, 1))
+            hps = [HalfPlane(v0,v1), HalfPlane(v1,v2),
+                   HalfPlane(v2,v3), HalfPlane(v3,v0)]
+            return min(*(abs(hp.signed_distance(p))
+                         for hp in hps))
         else:
-            surface_pt = crossings[1]
-        d = Vector(2*self.a*surface_pt.x + self.c*surface_pt.y + self.d,
-                   2*self.b*surface_pt.y + self.c*surface_pt.x + self.e)
-        return s * abs(d.dot(p - surface_pt) / d.length())
+            c = self.center
+            crossings = self.intersections(c, p)
+            # the surface point we want is the one closest to p
+            if (p - crossings[0]).length() < (p - crossings[1]).length():
+                surface_pt = crossings[0]
+            else:
+                surface_pt = crossings[1]
+            # n is the normal at surface_pt
+            n = self.gradient * surface_pt
+            n = n * (1.0 / n.length())
+            # returns the length of the projection of p - surface_pt
+            # along the normal
+            return -abs(n.dot(p - surface_pt))
 
 def Circle(center, radius, color=None):
     return Ellipse(color=color).transform(
diff --git a/rasterizer/rasterizer/geometry.py b/rasterizer/rasterizer/geometry.py
index e53a2f6f7..10aa4989f 100644
--- a/rasterizer/rasterizer/geometry.py
+++ b/rasterizer/rasterizer/geometry.py
@@ -6,6 +6,8 @@
 # using a variant of the standard quadratic formula
 # that behaves better numerically
 def quadratic(a, b, c):
+    if a == 0:
+        return -c/b, -c/b
     d = (b * b - 4 * a * c) ** 0.5
     if b >= 0:
         return (-b - d) / (2 * a), (2 * c) / (-b - d)
diff --git a/rasterizer/test_rasterizer.py b/rasterizer/test_rasterizer.py
index 2aa16b8fe..88e56be19 100644
--- a/rasterizer/test_rasterizer.py
+++ b/rasterizer/test_rasterizer.py
@@ -31,6 +31,7 @@ def do_it():
         LineSegment(Vector(0.6,0), Vector(0.6,1), 0.01, Color(1,0,0,1)),
         LineSegment(Vector(0.2,0.1), Vector(0.7, 0.4), 0.01, Color(0,0.5,0,1))])
     for shape in [
+        Circle(Vector(0.0, 0.0), 0.8, Color(1,0.5,0,1)).transform(scale(0.05,1)),
         Triangle([Vector(0.2,0.1), Vector(0.9,0.3), Vector(0.1,0.4)],
                  Color(1,0,0,0.5)),
         Triangle([Vector(0.5,0.2), Vector(1,0.4), Vector(0.2,0.7)],
@@ -51,12 +52,27 @@ def do_it():
     i.write_ppm(f)
     f.close()
 
-def test_ellipse():
-    Ellipse().transform(
-        scale(1.5, 1)).transform(
-        translate(0, 2)).transform(
-        rotate(math.radians(45))).transform(
-        scale(-1, 2))
+def test_quadratic():
+    for i in xrange(1000):
+        a = random.random()
+        b = random.random()
+        c = random.random()
+        if b * b - 4 * a * c >= 0:
+            v1, v2 = quadratic(a, b, c)
+            if v1 * v1 * a + v1 * b + c > 1e-5:
+                raise Exception("fail")
+            if v2 * v2 * a + v2 * b + c > 1e-5:
+                raise Exception("fail")
+    for i in xrange(1000):
+        a = 0
+        b = random.random()
+        c = random.random()
+        if b * b - 4 * a * c >= 0:
+            v1, v2 = quadratic(a, b, c)
+            if v1 * v1 * a + v1 * b + c > 1e-5:
+                raise Exception("fail")
+            if v2 * v2 * a + v2 * b + c > 1e-5:
+                raise Exception("fail")
 
 def test_inverse():
     for i in xrange(10000):
@@ -71,6 +87,6 @@ def test_inverse():
             raise Exception("foo")
 
 if __name__ == '__main__':
-    # test_inverse()
-    # test_ellipse()
+    test_inverse()
+    test_quadratic()
     do_it()

From b17264cd91ba07e8d89c76dbc1420ddfb0a98e30 Mon Sep 17 00:00:00 2001
From: Carlos Scheidegger <cscheid@research.att.com>
Date: Mon, 3 Mar 2014 02:31:18 -0500
Subject: [PATCH 07/16] added explanation for internal-point bound, simplified
 the code by reducing call to ConvexPoly

---
 rasterizer/rasterizer/ellipse.md | 54 +++++++++++++++++++++++++-------
 rasterizer/rasterizer/ellipse.py | 10 +++---
 2 files changed, 47 insertions(+), 17 deletions(-)

diff --git a/rasterizer/rasterizer/ellipse.md b/rasterizer/rasterizer/ellipse.md
index 4508464cd..fdd2dcc39 100644
--- a/rasterizer/rasterizer/ellipse.md
+++ b/rasterizer/rasterizer/ellipse.md
@@ -48,18 +48,39 @@ df/dy(x, y) = 0, v(x, y) = 0 gives the left and right points.
 
 ## Transforming an ellipse
 
-Say we're given an affine transformation T(x, y) = M . (x, y) = (x', y') that we
-need to apply to the ellipse to get a new ellipse. This means that the
-implicit equation for the new ellipse needs to be written as a
-function of x' and y', but our current equation is a polynomial in x
-and y. To get the new equation, we need to write x and y themselves as
-functions of x' and y'
-
-T(x, y) = (x', y')
-
-This gives 
-
-
+Say we're given an affine transformation T(x, y) = M . (x, y) = (x',
+y') that we need to apply to the ellipse to get a new ellipse. The
+straightforward way is to transform query points every time we want to
+test them against the new ellipse. That is, if we translated an
+ellipse by (10,0), then we need to transform the query points by the
+*inverse* of T. To see this, think of the translated ellipse
+above. Its new center is at (10,0), which means that query points at
+(10,0) need to be transformed to (0,0) on the old coordinate
+system. This approach works, but it only gives us the `contains`
+primitive. Worse yet, affine transformations do not preserve distances
+to shapes, so we cannot hope to make a query to
+`signed_distance_bound` and then invert that result somehow. Instead,
+we will use the inverse transformation to compute entirely new
+coefficients. We write (x, y) = M^-1(x', y'), which gives us an
+expression of the new coordinates in terms of the old ones, and then
+recreate the quadratic polynomial. This is tedious but
+straightforward algebra that goes sort of like this:
+
+x = m00 x' + m01 y' + m02
+y = m10 x' + m11 y' + m12
+(where m00, m01, etc. are the coefficients of the inverse matrix)
+
+f(x, y) = ax^2 + by^2 + cxy + dx + ey + f
+f(x, y) = a(m00 x' + m01 y' + m02)^2 + b(m10 x' + m11 y' + m12)^2
+        + ...
+
+Then you collect all terms that are quadratic in x' and call then a',
+all the terms quadratic in y' and call them b', etc. You end up with
+
+f'(x', y') = a'x'^2 + b'y^2 + c'x'y' + d'x' + e'y' + f'
+
+In the code, the variables a', b', c', etc. are denoted respectively
+as aa, bb, cc, etc.
 
 ## Bounding the distance to an ellipse
 
@@ -87,3 +108,12 @@ line too). This shortened line, in turn, is yet shorter than the leg
 `s'` to `p`. But length of this leg is precisely the distance from `s`
 to `i`, and so we know that it's a lower bound to the distance from
 `e` to `p`, and we are done.
+
+If `p` is inside the ellipse, then the problem is simpler. We first
+look for a convex polygon that inscribes the ellipse and contains
+`p`. This is simple: we simply shoot rays in the horizontal and
+vertical directions from `p` and collect the intersection points of
+those rays and the ellipse. Because ellipses are convex, a polygon
+that connects those intersections will be entirely inside the
+ellipse. Then, we simply return the signed distance bound for
+that polygon.
diff --git a/rasterizer/rasterizer/ellipse.py b/rasterizer/rasterizer/ellipse.py
index ac4d63df7..a72e57915 100644
--- a/rasterizer/rasterizer/ellipse.py
+++ b/rasterizer/rasterizer/ellipse.py
@@ -1,4 +1,5 @@
 from shape import Shape
+from poly import ConvexPoly
 from geometry import Vector, Transform, quadratic, scale, translate, AABox, HalfPlane
 from color import Color
 import sys
@@ -43,7 +44,9 @@ def transform(self, transform):
         return Ellipse(aa, bb, cc, dd, ee, ff, color=self.color)
     def intersections(self, c, p):
         # returns the two intersections of the line through c and p
-        # and the ellipse
+        # and the ellipse. Defining a line as a function of a single
+        # parameter u, x(u) = c.x + u * (p.x - c.x), (and same for y)
+        # this simply solves the quadratic equation f(x(u), y(u)) = 0
         pc = p - c
         u2 = self.a*pc.x**2 + self.b*pc.y**2 + self.c*pc.x*pc.y
         u1 = 2*self.a*c.x*pc.x + 2*self.b*c.y*pc.y \
@@ -65,10 +68,7 @@ def signed_distance_bound(self, p):
             # distance from the point to the ellipse must be smaller.
             v0, v2 = self.intersections(p, p + Vector(1, 0))
             v1, v3 = self.intersections(p, p + Vector(0, 1))
-            hps = [HalfPlane(v0,v1), HalfPlane(v1,v2),
-                   HalfPlane(v2,v3), HalfPlane(v3,v0)]
-            return min(*(abs(hp.signed_distance(p))
-                         for hp in hps))
+            return abs(ConvexPoly([v0,v1,v2,v3]).signed_distance_bound(p))
         else:
             c = self.center
             crossings = self.intersections(c, p)

From 40cef03c1406bd925516763f75c55d1b915256d5 Mon Sep 17 00:00:00 2001
From: Carlos Scheidegger <cscheid@research.att.com>
Date: Thu, 13 Mar 2014 16:08:52 -0400
Subject: [PATCH 08/16] add support for hex specification

---
 rasterizer/rasterizer/color.py | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/rasterizer/rasterizer/color.py b/rasterizer/rasterizer/color.py
index 7d328cc83..9656b7175 100644
--- a/rasterizer/rasterizer/color.py
+++ b/rasterizer/rasterizer/color.py
@@ -24,3 +24,13 @@ def byte(v):
                            byte(self.rgb[2] * self.a))
     def __repr__(self):
         return "[" + str(self.rgb) + "," + str(self.a) + "]"
+    @staticmethod
+    def hex(code, a=1):
+        if len(code) == 4:
+            return Color(int(code[1], 16) / 15.0,
+                         int(code[2], 16) / 15.0,
+                         int(code[3], 16) / 15.0, a)
+        elif len(code) == 7:
+            return Color(int(code[1:3], 16) / 255.0,
+                         int(code[3:5], 16) / 255.0,
+                         int(code[5:7], 16) / 255.0, a)

From 6dae39b06dc6530e4cf6ff1fb1e0bafc1177e4ef Mon Sep 17 00:00:00 2001
From: Carlos Scheidegger <cscheid@research.att.com>
Date: Thu, 13 Mar 2014 16:09:12 -0400
Subject: [PATCH 09/16] fix bug in transform from typo

---
 rasterizer/rasterizer/ellipse.py | 17 +++++++++++++----
 1 file changed, 13 insertions(+), 4 deletions(-)

diff --git a/rasterizer/rasterizer/ellipse.py b/rasterizer/rasterizer/ellipse.py
index a72e57915..925eb9078 100644
--- a/rasterizer/rasterizer/ellipse.py
+++ b/rasterizer/rasterizer/ellipse.py
@@ -23,6 +23,8 @@ def __init__(self, a=1.0, b=1.0, c=0.0, d=0.0, e=0.0, f=-1.0, color=None):
                                   Vector(-(d + c*y2)/2*a, y2),
                                   Vector(x1, -(e + c*x1)/2*b),
                                   Vector(x2, -(e + c*x2)/2*b))
+        if not self.contains(self.center):
+            raise Exception("Internal error, center not inside ellipse")
     def value(self, p):
         return self.a*p.x*p.x + self.b*p.y*p.y + self.c*p.x*p.y \
                + self.d*p.x + self.e*p.y + self.f
@@ -35,9 +37,9 @@ def transform(self, transform):
         bb = self.a*m01*m01 + self.b*m11*m11 + self.c*m01*m11
         cc = 2*self.a*m00*m01 + 2*self.b*m10*m11 \
              + self.c*(m00*m11 + m01*m10)
-        dd = 2*self.a*m00 * m02 + 2*self.b*m10*m12 \
+        dd = 2*self.a*m00*m02 + 2*self.b*m10*m12 \
              + self.c*(m00*m12 + m02*m10) + self.d*m00 + self.e*m10
-        ee = 2*self.a*m10*m02 + 2*self.b*m11*m12 \
+        ee = 2*self.a*m01*m02 + 2*self.b*m11*m12 \
              + self.c*(m01*m12 + m02*m11) + self.d*m01 + self.e*m11
         ff = self.a*m02*m02 + self.b*m12*m12 + self.c*m02*m12 \
              + self.d*m02 + self.e*m12 + self.f
@@ -50,11 +52,18 @@ def intersections(self, c, p):
         pc = p - c
         u2 = self.a*pc.x**2 + self.b*pc.y**2 + self.c*pc.x*pc.y
         u1 = 2*self.a*c.x*pc.x + 2*self.b*c.y*pc.y \
-             + self.c*c.y*pc.x + self.c*c.x*pc.y + self.d*pc.x \
+             + self.c*c.y*pc.x +   self.c*c.x*pc.y + self.d*pc.x \
              + self.e*pc.y
         u0 = self.a*c.x**2 + self.b*c.y**2 + self.c*c.x*c.y \
              + self.d*c.x + self.e*c.y + self.f
-        sols = quadratic(u2, u1, u0)
+        try:
+            sols = quadratic(u2, u1, u0)
+        except:
+            print self.value(c)
+            print self.value(p)
+            print self.a, self.b, self.c, self.d, self.e, self.f
+            print u2, u1, u0, c, p
+            raise
         return c+pc*sols[0], c+pc*sols[1]
     def signed_distance_bound(self, p):
         v = self.value(p)

From 8fe60bd999c11b82bd255e82d95f097b17f3cb56 Mon Sep 17 00:00:00 2001
From: Carlos Scheidegger <cscheid@research.att.com>
Date: Thu, 13 Mar 2014 16:09:28 -0400
Subject: [PATCH 10/16] add examples

---
 rasterizer/rasterizer/examples/__init__.py |  4 ++++
 rasterizer/rasterizer/examples/destijl.py  | 16 ++++++++++++++++
 rasterizer/rasterizer/examples/e1.py       |  6 ++++++
 rasterizer/rasterizer/examples/e2.py       | 15 +++++++++++++++
 rasterizer/rasterizer/examples/e3.py       | 19 +++++++++++++++++++
 rasterizer/run_examples.py                 | 17 +++++++++++++++++
 6 files changed, 77 insertions(+)
 create mode 100644 rasterizer/rasterizer/examples/__init__.py
 create mode 100644 rasterizer/rasterizer/examples/destijl.py
 create mode 100644 rasterizer/rasterizer/examples/e1.py
 create mode 100644 rasterizer/rasterizer/examples/e2.py
 create mode 100644 rasterizer/rasterizer/examples/e3.py
 create mode 100755 rasterizer/run_examples.py

diff --git a/rasterizer/rasterizer/examples/__init__.py b/rasterizer/rasterizer/examples/__init__.py
new file mode 100644
index 000000000..19f529fa8
--- /dev/null
+++ b/rasterizer/rasterizer/examples/__init__.py
@@ -0,0 +1,4 @@
+import e1
+import e2
+import e3
+import destijl
diff --git a/rasterizer/rasterizer/examples/destijl.py b/rasterizer/rasterizer/examples/destijl.py
new file mode 100644
index 000000000..26eb126e4
--- /dev/null
+++ b/rasterizer/rasterizer/examples/destijl.py
@@ -0,0 +1,16 @@
+from .. import *
+
+def painting():
+    scene = Scene()
+    scene.add(Rectangle(Vector(0.95, 0.0), Vector(1.0, 0.175), Color(0.75, 0.55, 0, 1)))
+    scene.add(Rectangle(Vector(0.0, 0.73), Vector(1.0, 0.76), Color(0,0,0,1)))
+    scene.add(Rectangle(Vector(0.94, 0), Vector(0.96, 1), Color(0,0,0,1)))
+    scene.add(Rectangle(Vector(0.95, 0.16), Vector(1, 0.19), Color(0,0,0,1)))
+    scene.add(Rectangle(Vector(0.25, 0.35), Vector(1.0, 1.0), Color(0.75,0,0,1)))
+    scene.add(Rectangle(Vector(0.0, 0.0), Vector(0.25, 0.35), Color(0.05,0,0.85,1)))
+    scene.add(Rectangle(Vector(0.24, 0.0), Vector(0.26, 1.0), Color(0,0,0,1)))
+    scene.add(Rectangle(Vector(0.0, 0.34), Vector(1.0, 0.36), Color(0,0,0,1)))
+    return scene
+
+def run(image):
+    painting().draw(image)
diff --git a/rasterizer/rasterizer/examples/e1.py b/rasterizer/rasterizer/examples/e1.py
new file mode 100644
index 000000000..a19b66782
--- /dev/null
+++ b/rasterizer/rasterizer/examples/e1.py
@@ -0,0 +1,6 @@
+from .. import *
+def run(image):
+    scene = Scene()
+    scene.add(Triangle([Vector(0.5, 0.5), Vector(0.8, 0.5), Vector(0.5, 0.8)],
+                       Color(1,0,0,1)))
+    scene.draw(image)
diff --git a/rasterizer/rasterizer/examples/e2.py b/rasterizer/rasterizer/examples/e2.py
new file mode 100644
index 000000000..7dde53f48
--- /dev/null
+++ b/rasterizer/rasterizer/examples/e2.py
@@ -0,0 +1,15 @@
+import destijl
+from .. import *
+
+def run(image):
+    painting = destijl.painting()
+    frame = Scene([Rectangle(Vector(-0.025,-0.025), Vector(1.025, 1.025), Color(0,0,0,1)),
+                   Rectangle(Vector(0,0), Vector(1, 1), Color(1,1,1,1))])
+    small_painting = Scene([frame, painting], scale(0.45, 0.45))
+    p1 = Scene([small_painting], translate(0.025, 0.025))
+    p2 = Scene([small_painting], translate(0.025, 0.525))
+    p3 = Scene([small_painting], translate(0.525, 0.025))
+    p4 = Scene([small_painting], translate(0.525, 0.525))
+    s = Scene([p1, p2, p3, p4])
+    s.draw(image)
+   
diff --git a/rasterizer/rasterizer/examples/e3.py b/rasterizer/rasterizer/examples/e3.py
new file mode 100644
index 000000000..c2b22e191
--- /dev/null
+++ b/rasterizer/rasterizer/examples/e3.py
@@ -0,0 +1,19 @@
+from .. import *
+import math
+
+def run(image):
+    s = Scene()
+    radius = 0.02
+    angle = 1
+    distance = 0.02
+    colors = [Color.hex("#1f77b4"), Color.hex("#ff7f0e"), Color.hex("#2ca02c"),
+              Color.hex("#d62728"), Color.hex("#9467bd"), Color.hex("#8c564b")]
+    for i in xrange(500):
+        obj = Circle(Vector(0.5, 0.5), radius, colors[i % len(colors)])
+        obj = obj.transform(translate(distance, 0))
+        obj = obj.transform(around(Vector(0.5, 0.5), rotate(angle)))
+        s.add(obj)
+        step = 3 * radius / distance
+        distance = distance + 0.009 * step
+        angle = angle + step
+    s.draw(image)
diff --git a/rasterizer/run_examples.py b/rasterizer/run_examples.py
new file mode 100755
index 000000000..bdc49890c
--- /dev/null
+++ b/rasterizer/run_examples.py
@@ -0,0 +1,17 @@
+#!/usr/bin/env python
+
+import rasterizer.examples as examples
+from rasterizer import *
+import subprocess
+
+def run_example(example, filename):
+    image = PPMImage(512, Color(1, 1, 1, 1))
+    example.run(image)
+    image.write_ppm(open(filename + '.ppm', 'w'))
+    subprocess.call(["convert", filename + '.ppm', filename + '.png'])
+
+run_example(examples.e1, 'e1')
+run_example(examples.destijl, 'destijl')
+run_example(examples.e2, 'e2')
+run_example(examples.e3, 'e3')
+

From 502bc58a8366e964012d6ecebea2571deab37f40 Mon Sep 17 00:00:00 2001
From: Carlos Scheidegger <cscheid@research.att.com>
Date: Wed, 30 Apr 2014 14:32:28 -0400
Subject: [PATCH 11/16] add stub definitions to base abstract class Shape

---
 rasterizer/rasterizer/shape.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/rasterizer/rasterizer/shape.py b/rasterizer/rasterizer/shape.py
index 7868c5808..9c4718c84 100644
--- a/rasterizer/rasterizer/shape.py
+++ b/rasterizer/rasterizer/shape.py
@@ -6,6 +6,11 @@
 class Shape:
     def __init__(self, color=None):
         self.color = color if color is not None else Color()
+        self.bound = None
+    def contains(self, p):
+        raise NotImplementedError("Undefined method")
+    def signed_distance_bound(self, p):
+        raise NotImplementedError("Undefined method")
     def draw(self, image, super_sampling = 6):
         if not self.bound.overlaps(image.bounds()):
             return

From eeaf978c7fa760dfae30b35d8814e86907fffe20 Mon Sep 17 00:00:00 2001
From: Carlos Scheidegger <cscheid@research.att.com>
Date: Fri, 2 May 2014 11:17:37 -0400
Subject: [PATCH 12/16] fix typo spotted by @EkkiD

---
 rasterizer/rasterizer/ellipse.md | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/rasterizer/rasterizer/ellipse.md b/rasterizer/rasterizer/ellipse.md
index fdd2dcc39..d6ced8ad6 100644
--- a/rasterizer/rasterizer/ellipse.md
+++ b/rasterizer/rasterizer/ellipse.md
@@ -103,9 +103,9 @@ to the distance between `p` and `s'`, since they form parallel
 lines. Now take the line from `p` to `e`. Notice that it intersects
 the line from `i` to `s'`, at a point we will call `e'`. Now notice
 that `e'`, `s'` and `p` form a right triangle, of which the line from
-`e'` to `s'` is the hypotenuse (and it's *shorter* than the original
-line too). This shortened line, in turn, is yet shorter than the leg
-`s'` to `p`. But length of this leg is precisely the distance from `s`
+`p` to `e'` is the hypotenuse (and it's *shorter* than the original
+line from `p` to `i`). This shortened line, in turn, is yet shorter than the leg
+`p` to `s'`. But length of this leg is precisely the distance from `s`
 to `i`, and so we know that it's a lower bound to the distance from
 `e` to `p`, and we are done.
 

From 02c66516187dfad321093550743e5173bb763c91 Mon Sep 17 00:00:00 2001
From: Carlos Scheidegger <cscheid@research.att.com>
Date: Fri, 2 May 2014 11:22:52 -0400
Subject: [PATCH 13/16] raise internal error if quadratic() fails due to no
 real roots

---
 rasterizer/rasterizer/ellipse.py | 8 ++------
 1 file changed, 2 insertions(+), 6 deletions(-)

diff --git a/rasterizer/rasterizer/ellipse.py b/rasterizer/rasterizer/ellipse.py
index 925eb9078..3c9bf343f 100644
--- a/rasterizer/rasterizer/ellipse.py
+++ b/rasterizer/rasterizer/ellipse.py
@@ -58,12 +58,8 @@ def intersections(self, c, p):
              + self.d*c.x + self.e*c.y + self.f
         try:
             sols = quadratic(u2, u1, u0)
-        except:
-            print self.value(c)
-            print self.value(p)
-            print self.a, self.b, self.c, self.d, self.e, self.f
-            print u2, u1, u0, c, p
-            raise
+        except ValueError:
+            raise Exception("Internal error, solutions be real numbers")
         return c+pc*sols[0], c+pc*sols[1]
     def signed_distance_bound(self, p):
         v = self.value(p)

From 04ce3bdbd4d90703701410331b0d2637e14b9bc2 Mon Sep 17 00:00:00 2001
From: Carlos Scheidegger <cscheid@research.att.com>
Date: Fri, 2 May 2014 11:29:00 -0400
Subject: [PATCH 14/16] replace Vector.union with AABox.from_vectors as @EkkiD
 suggested

---
 rasterizer/rasterizer/csg.py      | 4 ++--
 rasterizer/rasterizer/ellipse.py  | 8 ++++----
 rasterizer/rasterizer/geometry.py | 5 +++--
 rasterizer/rasterizer/poly.py     | 4 ++--
 4 files changed, 11 insertions(+), 10 deletions(-)

diff --git a/rasterizer/rasterizer/csg.py b/rasterizer/rasterizer/csg.py
index 47dc3795f..1117b70f5 100644
--- a/rasterizer/rasterizer/csg.py
+++ b/rasterizer/rasterizer/csg.py
@@ -13,8 +13,8 @@ def transform(self, t):
 class Union(CSG):
     def __init__(self, v1, v2, color=None):
         CSG.__init__(self, v1, v2, color=color)
-        self.bound = Vector.union(v1.bound.low, v1.bound.high,
-                                  v2.bound.low, v2.bound.high)
+        self.bound = AABox.from_vectors(v1.bound.low, v1.bound.high,
+                                        v2.bound.low, v2.bound.high)
     def contains(self, p):
         return self.v1.contains(p) or self.v2.contains(p)
     def signed_distance_bound(self, p):
diff --git a/rasterizer/rasterizer/ellipse.py b/rasterizer/rasterizer/ellipse.py
index 3c9bf343f..0518ffd55 100644
--- a/rasterizer/rasterizer/ellipse.py
+++ b/rasterizer/rasterizer/ellipse.py
@@ -19,10 +19,10 @@ def __init__(self, a=1.0, b=1.0, c=0.0, d=0.0, e=0.0, f=-1.0, color=None):
         self.center = self.gradient.inverse() * Vector(0, 0)
         y1, y2 = quadratic(b-c*c/4*a, e-c*d/2*a, f-d*d/4*a)
         x1, x2 = quadratic(a-c*c/4*b, d-c*e/2*b, f-e*e/4*b)
-        self.bound = Vector.union(Vector(-(d + c*y1)/2*a, y1),
-                                  Vector(-(d + c*y2)/2*a, y2),
-                                  Vector(x1, -(e + c*x1)/2*b),
-                                  Vector(x2, -(e + c*x2)/2*b))
+        self.bound = AABox.from_vectors(Vector(-(d + c*y1)/2*a, y1),
+                                        Vector(-(d + c*y2)/2*a, y2),
+                                        Vector(x1, -(e + c*x1)/2*b),
+                                        Vector(x2, -(e + c*x2)/2*b))
         if not self.contains(self.center):
             raise Exception("Internal error, center not inside ellipse")
     def value(self, p):
diff --git a/rasterizer/rasterizer/geometry.py b/rasterizer/rasterizer/geometry.py
index 10aa4989f..b3bce9e9c 100644
--- a/rasterizer/rasterizer/geometry.py
+++ b/rasterizer/rasterizer/geometry.py
@@ -29,8 +29,6 @@ def min(self, o):
         return Vector(min(self.x, o.x), min(self.y, o.y))
     def max(self, o):
         return Vector(max(self.x, o.x), max(self.y, o.y))
-    def union(*args):
-        return AABox(reduce(Vector.min, args), reduce(Vector.max, args))
     def length(self):
         return (self.x * self.x + self.y * self.y) ** 0.5
     def __repr__(self):
@@ -52,6 +50,9 @@ def overlaps(self, r):
                     r.low.y >= self.high.y or r.high.y <= self.low.y)
     def intersection(self, other):
         return AABox(self.low.max(other.low), self.high.min(other.high))
+    @staticmethod
+    def from_vectors(*args):
+        return AABox(reduce(Vector.min, args), reduce(Vector.max, args))
 
 class HalfPlane:
     def __init__(self, p1, p2):
diff --git a/rasterizer/rasterizer/poly.py b/rasterizer/rasterizer/poly.py
index b0000b095..d759195dd 100644
--- a/rasterizer/rasterizer/poly.py
+++ b/rasterizer/rasterizer/poly.py
@@ -1,11 +1,11 @@
 from shape import Shape
-from geometry import HalfPlane, Vector
+from geometry import HalfPlane, Vector, AABox
 
 class ConvexPoly(Shape): # a *convex* poly, in ccw order, with no repeating vertices
     def __init__(self, ps, color=None):
         Shape.__init__(self, color)
         self.vs = ps
-        self.bound = Vector.union(*self.vs)
+        self.bound = AABox.from_vectors(*self.vs)
         self.half_planes = []
         for i in xrange(len(self.vs)):
             h = HalfPlane(self.vs[i], self.vs[(i+1) % len(self.vs)])

From c3a506fc37de6a82b347e2a381d236711dcff1ee Mon Sep 17 00:00:00 2001
From: Carlos Scheidegger <cscheid@research.att.com>
Date: Fri, 2 May 2014 11:39:46 -0400
Subject: [PATCH 15/16] For didactic purposes, define a base class SceneObject
 from which Scene and Shape derive.

---
 rasterizer/rasterizer/scene.py | 8 ++++----
 rasterizer/rasterizer/shape.py | 6 +++++-
 2 files changed, 9 insertions(+), 5 deletions(-)

diff --git a/rasterizer/rasterizer/scene.py b/rasterizer/rasterizer/scene.py
index 8c31235e3..ff93c42ca 100644
--- a/rasterizer/rasterizer/scene.py
+++ b/rasterizer/rasterizer/scene.py
@@ -1,7 +1,7 @@
 from geometry import identity
-from shape import Shape
+from shape import Shape, SceneObject
 
-class Scene:
+class Scene(SceneObject):
     def __init__(self, nodes=None, transform=None):
         if transform is None:
             transform = identity()
@@ -12,8 +12,8 @@ def __init__(self, nodes=None, transform=None):
     def add(self, node):
         self.nodes.append(node)
     def draw(self, image):
-        for shape in self.traverse(identity()):
-            shape.draw(image)
+        for scene_object in self.traverse(identity()):
+            scene_object.draw(image)
     def traverse(self, xform):
         this_xform = xform * self.transform
         for node in self.nodes:
diff --git a/rasterizer/rasterizer/shape.py b/rasterizer/rasterizer/shape.py
index 9c4718c84..faf735301 100644
--- a/rasterizer/rasterizer/shape.py
+++ b/rasterizer/rasterizer/shape.py
@@ -3,7 +3,11 @@
 from geometry import Vector
 import random
 
-class Shape:
+class SceneObject:
+    def draw(self, image):
+        raise NotImplementedError("Undefined method")
+
+class Shape(SceneObject):
     def __init__(self, color=None):
         self.color = color if color is not None else Color()
         self.bound = None

From 9df7ffe6a572e8539c496069d2a1b63c9e534f75 Mon Sep 17 00:00:00 2001
From: Michael DiBernardo <mikedebo@gmail.com>
Date: Mon, 12 May 2014 07:03:37 -0400
Subject: [PATCH 16/16] Restructure contributor table.

- Add 'Online' column to provide space for contributors to represent
  themselves however they would like.
- Move Twitter handles into 'Online' column
---
 README.md | 132 ++++++++++++++++++++++++++++++++++--------------------
 1 file changed, 83 insertions(+), 49 deletions(-)

diff --git a/README.md b/README.md
index 954ea1a1f..91f78a58d 100644
--- a/README.md
+++ b/README.md
@@ -1,40 +1,33 @@
 *500 Lines or Less*
 ===================
 
-This is the source for the book *500 Lines or Less*,
-the fourth in the [Architecture of Open Source Applications](http://aosabook.org) series.
-As with other books in the series,
-all written material will be covered by the Creative Commons - Attribution license,
-and all code by the MIT License:
-please see the [license description](LICENSE.md) for details.
-In addition,
-all royalties from paid-for versions will all go to Amnesty International.
+This is the source for the book *500 Lines or Less*, the fourth in the
+[Architecture of Open Source Applications](http://aosabook.org) series.  As
+with other books in the series, all written material will be covered by the
+Creative Commons - Attribution license, and all code by the MIT License: please
+see the [license description](LICENSE.md) for details.  In addition, all
+royalties from paid-for versions will all go to Amnesty International.
 
 Mission
 -------
 
-Every architect studies family homes,
-apartments,
-schools,
-and other common types of buildings during her training.
-Equally,
-every programmer ought to know how a compiler turns text into instructions,
-how a spreadsheet updates cells,
-and how a browser decides what to put where when it renders a page.
-This book's goal is to give readers that broad-ranging overview,
-and while doing so,
-to help them understand how software designers think.
+Every architect studies family homes, apartments, schools, and other common
+types of buildings during her training.  Equally, every programmer ought to
+know how a compiler turns text into instructions, how a spreadsheet updates
+cells, and how a browser decides what to put where when it renders a page.
+This book's goal is to give readers that broad-ranging overview, and while
+doing so, to help them understand how software designers think.
 
-Contributions should not focus on the details of one algorithm
-or on the features of a particular language.
-Instead,
-they should discuss the decisions and tradeoffs software architects make when crafting an application:
+Contributions should not focus on the details of one algorithm or on the
+features of a particular language.  Instead, they should discuss the decisions
+and tradeoffs software architects make when crafting an application:
 
-*   Why divide the application into these particular modules with these particular interfaces?
+*   Why divide the application into these particular modules with these
+    particular interfaces?
 *   Why use inheritance here and composition there?
 *   Why multi-thread this but not that?
-*   When should the application allow for or rely on plugins,
-    and how should they be configured and loaded?
+*   When should the application allow for or rely on plugins, and how should
+    they be configured and loaded?
 
 Contribution Guidelines
 -----------------------
@@ -61,96 +54,137 @@ Contributors
     <th>Name</th>
     <th>Affiliation</th>
     <th>Project</th>
+    <th>Online</th>
     <th>GitHub</th>
-    <th>Twitter</th>
     <th>Email (if you choose)</th>
   </tr>
   <tr>
     <td>Mike DiBernardo</td>
     <td>freelance</td>
     <td>editorial</td>
-    <td>@MichaelDiBernardo</td>
-    <td></td>
+    <td>
+        <ul>
+            <li><a href="https://twitter.com/mdibernardo">@mdibernardo</a></li>
+            <li><a href="http://mikedebo.ca">mikedebo.ca</a></li>
+        </ul>
+    </td>
+    <td><a href="/MichaelDiBernardo">MichaelDiBernardo</a></td>
     <td>mikedebo@gmail.com</td>
   </tr>
   <tr>
     <td>Dustin Mitchell</td>
     <td>Mozilla</td>
     <td>cluster</td>
-    <td>@djmitche</td>
     <td>&nbsp;</td>
+    <td><a href="/djmitche">djmitche</a></td>
     <td>dustin@mozila.com</td>
   </tr>
   <tr>
     <td>Audrey Tang</td>
     <td>g0v.tw, Socialtext, Apple</td>
     <td>spreadsheet</td>
-    <td>@audreyt</td>
-    <td>@audreyt</td>
+    <td>
+        <ul>
+            <li><a href="https://twitter.com/audreyt">@audreyt</a></li>
+        </ul>
+    </td>
+    <td><a href="/audreyt">audreyt</a></td>
     <td>audreyt@audreyt.org</td>
   </tr>
   <tr>
     <td>Greg Wilson</td>
     <td>Mozilla</td>
     <td>web-server</td>
-    <td>@gvwilson</td>
-    <td>@gvwilson</td>
+    <td>
+        <ul>
+            <li><a href="https://twitter.com/gvwilson">@gvwilson</a></li>
+        </ul>
+    </td>
+    <td><a href="/gvwilson">gvwilson</a></td>
     <td>gvwilson@third-bit.com</td>
   </tr>
   <tr>
     <td>Kresten Krab Thorup</td>
     <td>Trifork</td>
     <td>torrent client</td>
-    <td>@krestenkrab</td>
-    <td>@drkrab</td>
+    <td>
+        <ul>
+            <li><a href="https://twitter.com/drkrab">@drkrab</a></li>
+        </ul>
+    </td>
+    <td><a href="/krestenkrab">krestenkrab</a></td>
     <td>krab@trifork.com</td>
   </tr>
   <tr>
     <td>Taavi Burns</td>
     <td>Points.com</td>
     <td>data-store</td>
-    <td>@taavi</td>
-    <td>@jaaaarel</td>
+    <td>
+        <ul>
+            <li><a href="https://twitter.com/jaaaarel">@jaaaarel</a></li>
+        </ul>
+    </td>
+    <td><a href="/taavi">taavi</a></td>
     <td>taavi.burns@points.com</td>
   </tr>
   <tr>
     <td>Guido van Rossum</td>
     <td>Dropbox</td>
     <td>crawler</td>
-    <td>@gvanrossum</td>
-    <td>@gvanrossum</td>
+    <td>
+        <ul>
+            <li><a href="https://twitter.com/gvanrossum">@gvanrossum</a></li>
+        </ul>
+    </td>
+    <td><a href="/gvanrossum">gvanrossum</a></td>
     <td>guido@python.org</td>
   </tr>
   <tr>
     <td>Erick Dransch</td>
     <td>Upverter</td>
     <td>Modeller</td>
-    <td>@EkkiD</td>
-    <td>@ErickDransch</td>
+    <td>
+        <ul>
+            <li><a href="https://twitter.com/ErickDransch">@ErickDransch</a></li>
+        </ul>
+    </td>
+    <td><a href="/EkkiD">EkkiD</a></td>
     <td>erick.dransch@upverter.com</td>
   </tr>
   <tr>
     <td>Sarah Mei</td>
     <td>Ministry of Velocity</td>
     <td>testing framework</td>
-    <td><a href='http://github.com/sarahmei'>@sarahmei</a></td>
-    <td><a href='http://twitter.com/sarahmei'>@sarahmei</a></td>
-    <td></td>
+    <td>
+        <ul>
+            <li><a href="https://twitter.com/sarahmei">@sarahmei</a></li>
+        </ul>
+    </td>
+    <td><a href="/sarahmei">sarahmei</a></td>
+    <td>&nbsp;</td>
   </tr>
   <tr>
     <td>Leah Hanson</td>
     <td>Google</td>
     <td>static analysis</td>
-    <td>@astrieanna</td>
-    <td>@astrieanna</td>
+    <td>
+        <ul>
+            <li><a href="https://twitter.com/astrieanna">@astrieanna</a></li>
+        </ul>
+    </td>
+    <td><a href="/astrieanna">astrieanna</a></td>
     <td>leah.a.hanson@gmail.com</td>
   </tr>
   <tr>
     <td>Christian Muise</td>
     <td>University of Melbourne</td>
     <td>flow-shop</td>
-    <td>@haz</td>
-    <td>@cjmuise</td>
+    <td>
+        <ul>
+            <li><a href="https://twitter.com/cjmuise">@cjmuise</a></li>
+        </ul>
+    </td>
+    <td><a href="/haz">haz</a></td>
     <td>christian.muise@gmail.com</td>
   </tr>
 </table>