# Scalable vector triangulations

, When I wrote “The tabular method” I drew this figure: it shows all the ways to triangulate an oriented convex polygon, from the unique way to triangulate a triangle up to the 14 ways to triangulate a hexagon. I drew it by hand using Inkscape, but I felt there was something unsatisfactory about that: if I can write a program to count the triangulations, I ought to be able to write a program to draw them!

I use Inkscape for this kind of diagram because its native format is Scalable Vector Graphics. The scalability is important for proofing the diagrams against future increases in graphics resolution. On “retina” (double-resolution) screens, single-resolution bitmap images stand out as jaggy and pixelated against the smoothness of the text rendering. So the generated diagrams need to be in SVG.1

What’s the best way to produce such diagrams from Python? If you look at the Python Package Index, there are tens of packages that mention “SVG” in their description. But I found a very helpful article by Florian Berger2 comparing the two Python packages `pySVG` and `svgwrite`.

Florian wrote, “to someone who needs a module to write SVG quickly, `svgwrite` feels a little more accessible because of the simple import and the good documentation,” and I followed his advice.

But first, generating the triangulations. You’ll recall that counting the triangulations was straightforward once the right recurrence relation had been found:

``````@memoized({2: 1})
def T(n):
"""Return number of ways to triangulate a convex polygon with n sides."""
return sum(T(k + 1) * T(n - k) for k in range(1, n - 1))
``````

The product of the number of triangulations in each component becomes the Cartesian product of the triangulations themselves, and this can be generated using Python’s `itertools.product`:

``````from itertools import product

def triangulations(p):
"""Generate the triangulations of the convex polygon p.
The sequence p consists of the vertices of the polygon.
Each triangulation in the output is a list of triangles, and each
triangle is a tuple of three vertices.

>>> list(triangulations(tuple('abcd')))
[[('a', 'b', 'd'), ('b', 'c', 'd')], [('a', 'b', 'c'), ('a', 'c', 'd')]]
>>> [sum(1 for _ in triangulations(range(i))) for i in range(3, 8)]
[1, 2, 5, 14, 42]

"""
n = len(p)
if n == 2:
yield []
elif n == 3:
yield [p]
else:
for k in range(1, n - 1):
for u, v in product(triangulations(p[:k + 1]), triangulations(p[k:])):
yield u + [(p, p[k], p[-1])] + v
``````

To compute the locations of the vertices in the triangulations, it’ll be convenient to use a vector or geometry package. But the situation with regard to vector packages in Python is as confusing as it was in the case of SVG, except that I can’t find a convenient review or comparison. Here are some quick thoughts:

• `euclid` by Alex Holkner looks like it has solid 2-dimensional and 3-dimensional geometric operations, but it’s not compatible with Python 3 (though it should be straightforward to make it so), it coerces all vector components to floats, and it enforces a distinction between vectors and points.

• `numpy` has fast vector operations, but it doesn’t have convenient functions for geometric operations like rotation.

• `vectors32` by Algis Kabaila supports 3-dimensional vectors and Python 3.2+ only.

• `wasabi.geom` by Daniel Pope has a nice interface, and solid 2-dimensional vector support. But it uses degrees for its angle representation.

It looks as though `wasabi.geom` would be suitable for this use case, but as it happens I have my own `geometry.Vector` class (available here on GitHub).

``````from __future__ import division
from collections import Sequence
from math import acos, atan2, cos, pi, sin, sqrt

__all__ = ['IncompatibleDimensions', 'Vector']

class IncompatibleDimensions(Exception):
pass

class Vector(tuple):
"""Vector is a subclass of tuple representing a vector in two or more
dimensions. You create a Vector by passing an iterable (that
yields the elements of the vector) to the constructor:

>>> Vector(range(5))
Vector(0, 1, 2, 3, 4)

or, for vectors of dimensions other than 1, by passing the
elements as arguments to the constructor:

>>> Vector(1, 2, 3)
Vector(1, 2, 3)

The first three components of a vector are accessible via the x,
y, and z properties:

>>> v = Vector(1, 2, 3)
>>> v.x, v.y, v.z
(1, 2, 3)

The elements must support numeric operations, but otherwise
a Vector is agnostic as to their type:

>>> Vector(1.5, 2.5, 3.5)
Vector(1.5, 2.5, 3.5)
>>> from fractions import Fraction
>>> Vector(Fraction(1, i) for i in (2, 3, 4))
Vector(Fraction(1, 2), Fraction(1, 3), Fraction(1, 4))

Vectors support vector addition and subtraction, and
multiplication and division by scalars:

>>> v, w = Vector(1, 2), Vector(3.5, 4.5)
>>> v + w
Vector(4.5, 6.5)
>>> w - v
Vector(2.5, 2.5)
>>> v * 10
Vector(10, 20)
>>> 100 * v
Vector(100, 200)
>>> -v
Vector(-1, -2)
>>> +w
Vector(3.5, 4.5)
>>> w / 0.5
Vector(7.0, 9.0)
>>> v * 9 // 10
Vector(0, 1)

You can combine vectors and ordinary sequences:

>>> Vector(0, 0) + (1, 2)
Vector(1, 2)

but the dimensions must match:

>>> Vector(0, 0) + (0, 0, 0)
... # doctest: +IGNORE_EXCEPTION_DETAIL
Traceback (most recent call last):
...
IncompatibleDimensions: (2, 3)

The magnitude property of a vector gives its magnitude, as does the
abs() function:

>>> Vector(3, 4).magnitude
5.0
>>> abs(Vector(5, 12))
13.0

Note that the magnitude is always a float because of the square
root. The magnitude_squared property gives the square of the
magnitude, and this varies according to the type of the elements:

>>> Vector(Fraction(1, 3), Fraction(1, 3)).magnitude_squared
Fraction(2, 9)

The angle_to() method computes the unsigned angle between two
vectors:

>>> from math import pi
>>> Vector(1, 0, 0).angle_to((0, 1, 0)) / pi
0.5

The distance() method computes the distance between two points
(expressed as vectors):

>>> Vector(3, -5, 1).distance((1, -2, 7))
7.0

The dot() method computes the dot product.

>>> Vector(1, 2, 3).dot((3, 2, 1))
10

The map() method applies a function to each element in the vector:

>>> Vector(1.5, 2.6, 3.4).map(round)
Vector(2, 3, 3)

The scaled() method returns a vector in the same direction but
with the specified magnitude; the normalized() method returns a unit
vector in the same direction; and the projected() method projects
another vector onto a vector.

>>> Vector(10, 20, 20).scaled(3)
Vector(1.0, 2.0, 2.0)
>>> Vector(0, 0, 5).normalized()
Vector(0.0, 0.0, 1.0)
>>> Vector(2, 0, 0).projected((5, 6, 7))
Vector(5.0, 0.0, 0.0)

Any of the above three methods may raise ZeroDivisionError if
applied to a vector with magnitude zero:

>>> Vector(0, 0).scaled(3)
... # doctest: +IGNORE_EXCEPTION_DETAIL
Traceback (most recent call last):
...
ZeroDivisionError: float division by zero
>>> Vector(0, 0, 0).normalized()
... # doctest: +IGNORE_EXCEPTION_DETAIL
Traceback (most recent call last):
...
ZeroDivisionError: float division by zero
>>> Vector(0, 0, 0, 0).projected((1, 2, 3, 4))
... # doctest: +IGNORE_EXCEPTION_DETAIL
Traceback (most recent call last):
...
ZeroDivisionError: division by zero

Three-dimensional vectors support one extra operation: the cross()
method computes the cross product:

>>> Vector(2, 1, 0).cross((3, 4, 0))
Vector(0, 0, 5)

Two-dimensional vectors support four extra operations.

The angle property is the signed angle (in radians) that the
vector makes with the x-axis:

>>> Vector(1, -1).angle / pi
-0.25

The cross() method returns the signed magnitude of the cross
product between the two vectors (interpreted as if they were
three- dimensional vectors lying in the XY-plane):

>>> Vector(3, 4).cross((2, 1))
-5

The perpendicular() method returns a vector that's perpendicular
to the vector:

>>> Vector(3, 4).perpendicular()
Vector(-4, 3)

The rotated() method rotates a vector through an angle in radians
and returns the result:

>>> Vector(2, 0).rotated(pi / 4)
... # doctest: +ELLIPSIS
Vector(1.414213562373..., 1.414213562373...)

"""
def __new__(cls, *args):
if len(args) == 1: args = args
return super(Vector, cls).__new__(cls, tuple(args))

def __repr__(self):
if len(self) == 1:
fmt = '{0}({1!r})'
else:
fmt = '{0}{1!r}'
return fmt.format(type(self).__name__, tuple(self))

def _check_compatibility(self, other):
if len(self) != len(other):
raise IncompatibleDimensions(len(self), len(other))

def _dimension_error(self, name):
return ValueError('.{0}() is not implemented for {1}-dimensional '
'vectors.'.format(name, len(self)))

if not isinstance(other, Sequence):
return NotImplemented
self._check_compatibility(other)
return Vector(v + w for v, w in zip(self, other))

if not isinstance(other, Sequence):
return NotImplemented
self._check_compatibility(other)
return Vector(w + v for v, w in zip(self, other))

def __sub__(self, other):
if not isinstance(other, Sequence):
return NotImplemented
self._check_compatibility(other)
return Vector(v - w for v, w in zip(self, other))

def __rsub__(self, other):
if not isinstance(other, Sequence):
return NotImplemented
self._check_compatibility(other)
return Vector(w - v for v, w in zip(self, other))

def __mul__(self, s):
return Vector(v * s for v in self)

def __rmul__(self, s):
return Vector(v * s for v in self)

def __div__(self, s):
return Vector(v / s for v in self)

def __truediv__(self, s):
return Vector(v / s for v in self)

def __floordiv__(self, s):
return Vector(v // s for v in self)

def __neg__(self):
return self * -1

def __pos__(self):
return self

def __abs__(self):
return sqrt(self.magnitude_squared)

@property
def angle(self):
"""The signed angle [-pi, pi] between this vector and the x-axis. For
two-dimensional vectors only.

"""
if len(self) == 2:
return atan2(self, self)
else:
raise self._dimension_error('angle')

def angle_to(self, other):
"""Return the unsigned angle [0, pi] to another vector. If either
vector has magnitude zero, raise ZeroDivisionError.

"""
if len(self) in (2, 3):
return atan2(abs(self.cross(other)), self.dot(other))
else:
# We avoid other.magnitude_squared so that this works if
# other is a plain sequence rather than a Vector object.
return acos(self.dot(other)
/ sqrt(self.magnitude_squared
* sum(v * v for v in other)))

def cross(self, other):
"""Return the cross product with another vector. For two-dimensional
and three-dimensional vectors only.

"""
self._check_compatibility(other)
if len(self) == 2:
return self * other - self * other
elif len(self) == 3:
return Vector(self * other - self * other,
self * other - self * other,
self * other - self * other)
else:
raise self._dimension_error('cross')

def distance(self, other):
"""Return the distance to another vector (understanding both vectors
as points).

"""
return abs(self - other)

def dot(self, other):
"""Return the dot product with the other vector."""
self._check_compatibility(other)
return sum(v * w for v, w in zip(self, other))

@property
def is_zero(self):
"""True if this vector has magnitude zero, False otherwise."""
return self.magnitude_squared == 0

@property
def magnitude(self):
"""The magnitude of the vector."""
return abs(self)

@property
def magnitude_squared(self):
"""The squared magnitude of the vector."""
return sum(v * v for v in self)

def map(self, f):
"""Return the vector whose elements are the result of applying the
function f to the elements of this vector.

"""
return Vector(f(v) for v in self)

@property
def non_zero(self):
"""False if this vector has magnitude zero, True otherwise."""
return self.magnitude_squared != 0

def normalized(self):
"""Return a unit vector in the same direction as this vector. If this
has magnitude zero, raise ZeroDivisionError.

"""
return self / abs(self)

def perpendicular(self):
"""Return a vector perpendicular to this vector. For two-dimensional
vectors only.

"""
if len(self) == 2:
return Vector(-self, self)
else:
raise self._dimension_error('perpendicular')

def projected(self, other):
"""Return the projection of another vector onto this vector. If this
vector has magnitude zero, raise ZeroDivisionError.

"""
return self * (self.dot(other) / self.magnitude_squared)

def rotated(self, theta):
origin. For two-dimensional vectors only.

"""
if len(self) == 2:
s, c = sin(theta), cos(theta)
return Vector(self.dot((c, -s)), self.dot((s, c)))
else:
raise self._dimension_error('rotated')

def scaled(self, s):
"""Return a vector of magnitude s in the same direction as this vector.
If this has magnitude zero, raise ZeroDivisionError.

"""
return self * (s / abs(self))

@property
def x(self):
return self

@property
def y(self):
return self

@property
def z(self):
return self
``````

Here’s the code for drawing the triangulations:

``````from itertools import cycle
from math import pi
from svgwrite import Drawing
from vector import Vector

def draw_triangulations(n, filename, radius, columns, edge_style, triangle_styles=()):
"""Draw the triangulations of an oriented n-sided convex polygon and
save it into filename as SVG. Other arguments specify the radius in
pixels of each polygon, the number of columns in the drawing, and
the styles for drawing the edges and the triangles.

"""
polygon = [Vector(1, 0).rotated(i * 2 * pi / n) for i in range(n)]
fill_styles = cycle(triangle_styles)
tt = list(triangulations(polygon))
rows = (len(tt) + columns - 1) // columns
size = (Vector(columns, rows) * 2.1 - Vector(0.1, 0.1)) * radius
drawing = Drawing(filename=filename, size=size)
for i, triangulation in enumerate(tt):
origin = Vector(i % columns, i // columns) * 2.1 + Vector(1, 1)
edges = set()
for triangle in triangulation:
a, b, c = ((v + origin) * radius for v in triangle)
if triangle_styles:
for edge in ((a, b), (b, c), (c, a)):
for edge in edges:
drawing.save()
``````

I’m rubbish when it comes to designing colour palettes or even recognizing good colour schemes, but luckily for me there’s a very useful site `colourlovers.com` where designers submit palettes and people vote on them. The colours in the figure below (showing the 42 ways to triangulate a heptagon) are from this palette by user manekineko.

``````>>> colours = '#BBBB88 #CCC68D #EEDD99 #EEC290 #EEAA88'.split()
>>> triangle_styles = [dict(fill=c) for c in colours]
>>> edge_style = dict(stroke='black', stroke_width=2)
>>> draw_triangulations(7, 'heptagon.svg', 50, 14, edge_style, triangle_styles)
`````` And finally, here are the 429 ways to triangulate a nonagon:

``````>>> draw_triangulations(9, 'nonagon.svg', 50, 33, edge_style, [dict(fill='#eee')])
`````` But the SVG file is rather big: at 1.2 MiB it seems rather large for an image which has 3003 triangles and 6435 lines. See the follow-up for how I was able to reduce it to just over 200 KiB.

1.  SVG still feels a bit cutting-edge, and there are lot of services that don’t support it yet. For example, Google+ ignores SVG images when looking for a representative image for a link, and you can’t upload SVG images to Google+, Facebook, or Twitter. But the image format will get there. It took PNG more than a decade to take over from GIF, after all.

2.  The Python package world is so well populated that there’s a real need for more comparisons along these lines. What’s the best bitmap-handling package? What are the advantages and disadvantages of the many XML-processing tools?