Triangulations...

... or how I learned that eventually I'll have to stop worrying and love a compiler.

So last week I wrote about how I was re-implementing all the routines needed in Larry's Numerical Approximation with Splines course in Python, because (A) I'm too cheap to buy Matlab, (B) I want to give my programming skills a workout, and (C) the cognitive style of Matlab doesn't agree with me. My examples in that post were pretty sketchy; since I just spent lots of time writing up a pretty nontrivial module, I figured I'd share it to give a better idea of what I mean.

Standard disclaimers: all the code here is my own, and copyright 2013 to me. I am making it available under the MIT license, with the request that I get attribution if it is used or incorporated into any project.

Here's the problem specification:

Definition: a triangulation of a region \(R \subset \mathbb{R}^2\) in the plane is a finite set \(T\) of triangles such that \[ R = \bigcup_{\Delta \in T} \Delta \]and such that the intersection of any two triangles in \(T\) is either the empty set, a single vertex, or a single edge. A triangulation can be viewed as a simplicial complex (whose 1-dimensional reduct is clearly a planar graph), which I will refer to as \[ \langle T^{(0)}, T^{(1)}, T^{(2)} \rangle = \langle V, E, F \rangle \]

Problem(s):

• Given numerical arrays V (type float, size \(|V| \times 2\)) and F (type int, size \(|F| \times 3\)), to represent the triangulation which has a triangle with vertices V[k1,:], V[k2,:], V[k3,:] for each row [k1,k2,k3] of F, if any such triangulation exists. (It is a separate, but worthwhile, problem to think about checking if this data actually does code a triangulation. Even checking whether any of the triangles in F overlap would require a bit of doing. I haven't actually implemented any of this -- you can think about how to use the tools I've defined here to actually perform such a check if you want to.)
• Once such a representation has been made, given a point t = (tx,ty), to decide if \( t \in R \).
• Once we've decided that \(t \in R\), to return some triangle containing \(t\).

Here's the API for the solution to the first part, as it appears in Larry's Matlab suite:

[nb,ne,nt,v1,v2,v3,e1,e2,e3,ie1,ie2,tril,trir,bdy,vadj, [...]
 [...] eadj, adjstart,tadj,tstart,area,TRI] = trilists(V,F)

([...] indicates an implicit linebreak.) See what I mean about how, under this cognitive style, the user must be the smartest player in the whole game? The data above encode indexed lists of vertices, edges, triangles, and other useful data; each of the output variables is a matrix holding some such data. For example, ie1 and ie2 are the \(|E| \times 1\) vectors, holding indices of the initial and terminal vertices of edges \(e \in E\).

Obviously, I have to do better than this, where "better" means "encode a triangulation using one object, not twenty-something".

OK, so we have the problem statement above. Before we dive into writing code, let's plan out what kind of objects we'll want. What do I mean by "what kind"? I mean that we need to answer

What do we want our objects to do, and do efficiently?

Let's start sketching from the top; we'll step down into lower-level objects pretty quickly, then come back up when we have the little guys up to spec. This process might repeat a few times as we find out that we need more structure from our little objects.

Here's my idea for representing a triangulation: it should "be" a graph with additional structure. Why do I say this? Because (A) we have a whole bunch of robust tools and algorithms for dealing with graphs, and (B) the structure of \(T\) is entirely determined by the combination of its graph structure and the embedding into \( \mathbb{R}^2 \).

OK, so how do you implement a graph? My favorite is as a dictionary, with items of the form

T = {vertex: set of neighbors}

which is a variation on the classical adjacency-list representation.

What does this already imply about vertices? Well, like all dictionary keys or set elements, they must be hashable. In particular, we should not make a vertex out of a list or other mutable type. One option would be to have the vertices be the integers \(0, 1, \ldots, |V|-1\), and store a record somewhere of their embedded locations on the plane; but that's not the way an object-oriented person thinks. No, I want a vertex to have its location as its root-level data if possible. OK, we say, let's have a vertex be a tuple v = (vx,vy) where each of vx, vy is a float.

The other consideration that might come to mind is that, when using this kind of representation, you have a duplication of edges: if \(a\) has an edge to \(b\), then \(b\) has an edge to \(a\) too. To eliminate doing some things twice, it's nice to have a way to rank-order the vertices so that we can do some things only if \(a < b\). If we were using integers as vertices, this would be built-in; since we're not, we should code some kind of comparison method. But now we're out of the built-in functionality of tuples: we have to define a Vertex class:

class Vertex(tuple):
    """A vertex object for a triangulation. A Vertex is a 2-tuple
with lexicographic comparison.

Call with one argument v = vertex(t), with t a list or tuple of numbers
the coordinates of v will be t[0] and t[1].

Or call with two arguments v = vertex(vx,vy) where each of vx,vy is numeric.

Vertices are hashable and immutable.
"""
    
    def __new__(cls,x,y = None):
        try:
            retobj = tuple.__new__(cls,(np.float(x[0]),np.float(x[1])))
        except (TypeError,IndexError):
            try:
                retobj = tuple.__new__(cls,[np.float(x),np.float(y)])
            except TypeError:
                raise TypeError('Coordinates must be numeric!')
            return retobj

    def __lt__(self,other):
        return self[0] < other[0] or (self[0] == other[0] and \
self[1] < other[1])

    def __le__(self,other):
        return self[0] < other[0] or (self[0] == other[0] and \
self[1] <= other[1])

It's worth pointing out something here: we only decided to define a new class when we found that we needed functionality beyond what the builtin class provides. By way of comparison: we won't need our edges to do anything; so we won't actually define a corresponding Edge class at all.

We will, however, need our triangles do do things. Specifically, we want them to be able to test for membership. It turns out that the most efficient way to achieve this is using barycentric coordinates (which I'd only ever encountered in the context of algebraic topology -- who knew they were actually useful!)

import numpy as np

class Triangle(tuple):
    """A triangle object for a triangulation."""

    def __new__(cls,A,B,C):
        L = [A,B,C]
        if not all([isinstance(x,Vertex) for x in L]):
            raise Exception('Arguments to the Triangle constructor must be \
instances of the Vertex class!')
        L.sort()
        det = (L[1][0] - L[0][0])*(L[2][1] - L[1][1]) + \
(L[1][0] - L[2][0])*(L[1][1] - L[0][1])
        if det > 0:
            retobj = tuple.__new__(cls,L)
            retobj.det = det
        elif det < 0:
            retobj = tuple.__new__(cls,[L[0],L[2],L[1]])
            retobj.det = -det
        else:
            raise Exception('Triangle is degenerate!')
        
        retobj.M = np.matrix([np.ones( (3,) ), \
[p[0] for p in retobj], [p[1] for p in retobj]])
        return retobj

    def barycoord(self,v):
        return np.linalg.solve(self.M,[1,v[0],v[1]])

    def contains(self,v):
        return all(np.vectorize(lambda t: t >= 0)(self.barycoord(v)))

Of course, when we're all done, the import command should come at the very beginning of the module.

I've set up the Triangle class so that the lex-least vertex always comes first, and then the others come in counter-clockwise order. This isn't strictly necessary; but it does make equality checking much nicer.

The really nice fact about barycentric coordinates which makes this whole module tick, is the following

Fact 1: Let \(\Delta = ABC\) be a triangle, \(P\) a point, and suppose \(P = b_A A + b_B B + b_C C \) in barycentric coordinates. Then \(P\) lies on the same side of segment \( AB \) as \( \Delta \) if and only if \( b_C > 0\). If \( b_C = 0 \), then \(P\) is collinear with \( A,B\).

This is really nice.

OK, so let's go back to thinking about what we'd like a triangulation to be able to do. Well, triangles automatically know what their vertices are; but each vertex should also know what triangles it's on. We'll have to store that data when we create the triangulation object.

Next, our triangulation should know what's on the boundary versus the interior. Here's a tricky way to achieve that. Remember how we put the lexicographic order on vertices? This gives each edge an orientation, if we want to use one; hence each edge can talk about the triangle "on its left" (or right). An edge is a boundary edge precisely if it only has one or the other!

This is all that we need to write the class constructor:

class Triangulation(dict):
    """Call T = Triangulation(vertices,triangles).

Base data: adjacency lists of graph reduct of the triangulation.

vertices should be V x 2 numpy array of xy-positions of vertices.

triangles should a F x 3 numpy array of indices of vertices forming
the triangles of T.
"""

    def __init__(self,vertices,triangles):
        """Constructor"""
        if vertices.ndim != 2 or vertices.shape[1] != 2 or \
triangles.ndim != 2 or triangles.shape[1] != 3:
            raise Exception('Misshapen arguments; \
vertices must be V x 2 and triangles must be F x 3!')

        vertices = np.asarray(vertices)
        triangles = np.asarray(triangles,dtype = np.int)
        
        Pdict = {k:Vertex(vertices[k]) for k in range(len(vertices))}
        dict.__init__(self,[[Pdict[k],set()] for k in range(len(vertices))])

        self.tri_adj = {v:set() for v in self.keys()}
# to be: set of triangles containing vertex v

        self.triangleright = {} 
# to be eventually: (v,w):triangle to the right of edge (v,w), if any.
        self.triangleleft = {} 
# likewise

        for k1,k2,k3 in triangles:
            current_triangle = Triangle(Pdict[k1],Pdict[k2],Pdict[k3])
            v1,v2,v3 = current_triangle
            self[v1].update([v2,v3])
            self[v2].update([v3,v1])
            self[v3].update([v1,v2])
            self.tri_adj[v1].add(current_triangle)
            self.tri_adj[v2].add(current_triangle)
            self.tri_adj[v3].add(current_triangle)
            
            try:
                self.triangleleft[frozenset([v1,v2])]
                raise Exception('Contradictory input:\n\nEdge from vertex {v} \
to vertex {w} has two triangles to its left!'.format(v = v1,w = v2))
            except KeyError:
                self.triangleleft[frozenset([v1,v2])] = current_triangle

            try:
                self.triangleright[frozenset([v1,v3])]
                raise Exception('Contradictory input:\n\nEdge from vertex {v} \
to vertex {w} has two triangles to its right!'.format(v = v1,w = v3))
            except KeyError:
                self.triangleright[frozenset([v1,v3])] = current_triangle

            if v2 < v3:
                # then current_triangle is on the left of {v2,v3}
                try:
                    self.triangleleft[frozenset([v2,v3])]
                    raise Exception('Contradictory input:\n\nEdge from vertex {v} \
to vertex {w} has two triangles to its left!'.format(v = v2,w = v3))
                except KeyError:
                    self.triangleleft[frozenset([v2,v3])] = current_triangle

            else:
                #then current_triangle is on the right of {v3,v2}
                try:
                    self.triangleright[frozenset([v2,v3])]
                    raise Exception('Contradictory input:\n\nEdge from vertex {v} \
to vertex {w} has two triangles to its right!'.format(v = v3,w = v2))
                except KeyError:
                    self.triangleright[frozenset([v2,v3])] = current_triangle

        self.boundaryedges = self.triangleleft.symmetric_difference\
(self.triangleright)
        self.boundaryvertices = set().union(*self.boundaryedges)

Notice that T.triangleleft and T.triangleright are dictionaries keyed by edges; I could have had edges be tuples, but in order to minimize the tedious order-checking, I've put them as frozensets instead.

Now, this wasn't one of the original problems, but I am a visual guy, so I wanted a way to draw pictures of my triangulations (passing plot style arguments through as necessary):

import matplotlib.pyplot as plt

    def draw(self,**kwargs):
        """T.draw() draws the triangulation. Keword arguments are passed to plot.
    """
        try:
            color = kwargs.pop('c')
        except KeyError:
            try:
                color = kwargs.pop('color')
            except KeyError:
                color = 'k'

        wasinteractive = plt.isinteractive()
        plt.ioff()
    
        for v in self:
            for w in self[v]:
                if v < w:
                    plt.plot([v[0],w[0]],[v[1],w[1]],c = color,**kwargs)
    
        if wasinteractive:
            plt.show()
            plt.ion()

This is an instance method of every triangulation object; as you can see, it defaults to drawing the graph in black. (Without that slightly kludgy declaration, every edge would be a different color.) Notice how we've saved ourselves labor by using the ordering on vertices to only draw each edge once. Here's the figure

T.draw(c = 'b')

produced for one of Larry's small sample triangulations:

Remember those bullet points way up at the beginning? We're done with the first one. Another way to spin this is that, so far, all of this is "plumbing": the basic stuff that has to be in place before we can even hope to do anything interesting. That's about to change, right about... now!

Recall that the second problem was to determine (efficiently) whether a point \(t\) lies within the triangulated region or not. One way to do this, of course, would be to run ABC.contains(t) for each triangle ABC. That would, however, be a poor solution to this problem.

However, it took me a fair bit of thinking before I found a criterion which was both correct and more efficient:

Obvious Fact 2: Let \(T\) be any triangulation of the region \(R\), and \(t\) any point. Furthermore, let \(e \in E\) be any closest boundary-edge to \(t\). Then the straight-line path from \(t\) to \(e\) is either entirely inside or entirely outside \(R\).

Nonobvious Fact 3: Let \(e = W -- X\) be as in Fact 2. Since \(e\) is a boundary edge, it subtends only one triangle \( \Delta = UWX \). If the closest point of \(e\) to \(t\) is not an endpoint, then \(t \in R\) if and only if \(t\) and \(U\) are on the same side of \(e\). On the other hand, if the closest point of \(e\) is an endpoint, say \(X\), and \(XY\) is the other boundary edge containing \(X\), then

• if \(t\) and \(Y\) are on the same side of \(WX\), then \(t \notin R\)
• if \(t\) and \(Y\) are on opposite sides of \(WX\), then \(t \in R\)
• if \(WXY\) is a straight line, then \(t \in R \) iff \(t\) and \(U\) are on the same side of \(WXY\)

(You should try to prove this.) This proves the correctness of the following algorithm:

dpp = lambda v,w:(v[0]-w[0])**2 + (v[1]-w[1])**2 #distance squared from v to w

def dpe(t,e):
    """d = dpe(t,e) returns the square of the distance from point t to edge e."""
    v,w = e
    vw2 = dpp(v,w)
    tv2,tw2 = dpp(t,v),dpp(t,w)
    if tv2 == tw2:
        return tv2 - vw2/4
    elif tw2 < tv2:
        v,w,tv2,tw2 = w,v,tw2,tv2 # ensures that v is closer to t than w is

    if tw2 >= tv2 + vw2:
        return tv2
    else:
        return (np.linalg.det([[1,1,1],[t[0],v[0],w[0]],[t[1],v[1],w[1]]])**2)/vw2

    def contains(self,t):
        """Call (retBool, retVertex) = T.contains(t) to test whether t belongs to 
the region covered by the triangles of T.

retBool will be a Boolean, True iff T does contain t.
retVertex will be a closest boundary vertex to t.
"""
        wx = min(self.boundaryedges, key = lambda e: dpe(t,e))
        w,x = wx
        if dpp(t,w) < dpp(t,x):
            x,w = w,x
        dtwx,dtw,dtx,dwx = dpe(t,wx),dpp(t,w),dpp(t,x),dpp(w,x)

        try:
            UWX = self.triangleleft[wx]
        except KeyError:
            UWX = self.triangleright[wx]
        u = next( (u for u in UWX if u not in {w,x}) )

        if dtwx < dtx:
            #then the perpendicular from t to wx cuts it in the middle somewhere
            return (UWX.barycoord(t)[UWX.index(u)] >= 0,x)
        else:
            xy = next( (e for e in self.boundaryedges if x in e and w not in e) )
            y = next( (y for y in xy if y is not x) )
            y_u = UWX.barycoord(y)[UWX.index(u)]
            if y_u < 0:
                return(True,x)
            elif y_u > 0:
                return(False,x)
            else:
                return (UWX.barycoord(t)[UWX.index(u)] >= 0, x)

A note on speed: this method is one of the serious performance bottlenecks in the whole module. This is because the very first thing we do is take a min over the (potentially very large) set of boundary edges. If in practice we can do without this overhead, that could represent a big performance boost.

We've now solved the second problem. The third is just a straightforward application of depth-first search:

    def findtriangle(self,t,guess = None):
        """ABC = T.findtriangle((tx,ty)) returns a triangle such that 
(tx,ty) lies inside ABC. This is not uniquely defined if (tx,ty) lies
 on a triangle boundary, so the method will return the first one it finds.

Currently programmed to return None if (tx,ty) is not inside any
 of the triangles of T.
"""
        #first try all triangles adjacent to the guess
        if guess is not None:
            for ABC in self.tri_adj[guess]:
                if ABC.contains(t):
                    return ABC

        (t_in_region,root) = self.contains(t)
        
        if t_in_region:
            point = root
            DFSparent = {root:None}
            ABCvisited = set()

            while point is not None:
                ABCs = (ABC for ABC in self.tri_adj[point] \
if ABC not in ABCvisited)
                while True:
                    try:
                        ABC = next(ABCs)
                        ABCvisited.add(ABC)
                        if ABC.contains(t):
                            return ABC
                    except StopIteration:
                        try:
                            newpoint = min((v for v in self[point] \
if not v in DFSparent),key = lambda w: dpp(t,w))
                            DFSparent[newpoint] = point
                            point = newpoint
                            break
                        except ValueError:
                            point = DFSparent[point]
                            break

And that's it! We've got a working triangulation module.

The only question now is, how does it perform? And the answer is, a lot more slowly than Matlab, at least if Larry's benchmarks are correct. Which they probably are. The good news is, it looks like the performance hit is a constant multiple (a factor of \(2^9\) or so); so that I've got an algorithm with the right big-O characteristics.

Which brings me back to where I started this post: if I were interested in using this in an application where a thousand calls to T.findtriangle take milliseconds and not seconds, clearly plain old vanilla interpreted Python is not the way to go. And indeed, I'm far from the first person to run into this phenomenon. So many people have, in fact, and yet so popular is Python among people who operate on large amounts of data, that there's a whole project, Cython (Compiled Python) designed towards the end of bridging this gap. I have a feeling the day will come when a time savings factor of \(2^9\) will be worth diving into that stew. 

But that day is not this day.