Inside Boundary Generator: Computational geometry for the GIS programmer

I'm going to take a post or two to dive into the internals of Boundary Generator, my GPL geoprocessing tool that had its "alpha" release last week. The nice thing about open source is that I can point directly at the code as I go along. Up this dreary Saturday afternoon: The underlying computational geometry (some code but mostly math!) behind finding shared boundaries. In my experience, it's fairly rare to get your hands dirty with this when doing GIS development, so hopefully others will find it as interesting as I do. (This is also fairly elementary as far as computational geometry goes; this blog post has nothing that wouldn't be taught in the first few sessions of an undergraduate computer graphics course, or even an expected prerequisite.)

The problem that we need to solve is that of finding overlapping line segments. In addition, our solution should have well-defined tolerances and not be subject to (too much) numerical instability.

Representing geometry

When sitting down to code, the first problem is that of representation: How do we store a line segment? Most peoples' instinctive representation is as a pair of endpoints -- (x0,y0),(x1,y1). But a bit of thought shows that there are many other options. For example, we could store line segments in a more algebraic form, such as the values m, b, x0 and x1 from the following description of a line segment:

• y = mx + b where x0 <= x <= x1

Although this isn't as natural as simply working with endpoints, it has some big advantages over simply storing endpoints for finding overlaps:

• If two line segments have the same slope (m) and y-intercept (b), we can see if they overlap just by checking their x0..x1 ranges for an overlap.
• We can also test for line crossings by solving a system of two equations (y = (m0)(x) + b0 and y = (m1)(x) + b1) for x and y. Boundary Generator doesn't need this, but other geoprocessing functions in the future will.
  

However, this approach is plagued with numerical instability and special cases -- what do you do with lines that are vertical or near-vertical? -- that will quickly turn code into spaghetti as you add special cases to address these problems. It is also difficult to specify overlap tolerances in a meaningful way.

Instead, boundary generator uses a more abstract (and much better) representation that dives a little further (but not too far) past basic algebra: It is a parametric representation that consists of an origin point (o), a vector (v), and a scalar parameter (t). Just as an infinite line may be defined by a slope and y-intercept value, it may instead be defined as an origin and a vector; and just as a line segment can be defined upon an infinite line as a range of x values, a line segment can be defined as a range of t values. The line is all points p which satisfy p = vt + o; the line segment is the portion of the line where 0 <= t <= 1. If the preceeding made sense you may skip the next section or two. :)

Vectors (Go to the code; go to a more in-depth Wikipedia article.)

What is a vector? It is this: a magnitude (length), and a direction. It is typically and most conveniently represented as an (x,y) pair; this defines both direction and length. (The direction is the direction from the origin (0,0) through the point (x,y) and the magnitude is the length of a line segment from the origin to (x,y) ).

Vectors can be added to each other; if we have two vectors a and b, a + b is the vector (ax + bx, ay + by). (They can also be subtracted.) You can also scale a vector by an amount s, changing its length without changing its direction, simply by multiplying both x and y by the scale amount: (sx, sy).

We can also (obviously, since a vector is just an x,y pair) choose to represent a point via a vector simply by interpreting it as the magnitude and direction from the origin.

Line segments (Go to the code.)

Our parametric representation of an infinite line (repeating from before) is: a starting point o, a direction/magnitude v, and a parameter t: vt + o. By varying t, we get all points on the line. When t is 0 it gives us the starting point. We can also test to see if a point p is on the line by trying to solve p = vt +o for t.

And to turn our infinite line into a line segment: simply restrict t to lie between zero and one. Since v already includes a length as well as a direction, the line segment is nicely defined -- the start point is o; the end point is o + v (by setting t to 1.)

This also is a very, very good representation for finding overlapping segments as well as crossing segments. But first, I need to introduce...

The vector dot product (Go to the code; go to the Wikipedia article)

This is the swiss-army-knife of vector operations. It takes in two vectors, and returns (the length of the first vector) * (the length of the second vector) * (the cosine of the angle between them). It is also very easy to calculate. This has a ton of uses, to name a few:

• The dot product of a vector with itself is equal to the square of its length.
• The angle between two vectors a and b, both of magnitude 1, is acos(dot(a,b)).
• If two vectors are perpendicular, their dot product will be zero.

The code uses the dot product in all these ways. For testing whether line segments overlap, we need to ensure they are parallel within a specific tolerance. This can be done by normalizing the vectors (multiplying each by 1/length to preserve its direction while making its length 1), then taking the dot product and inverse cosine; the two are parallel if the angle is within the tolerance to zero or 180 degrees.

The actual boundary generator implementation of AreParallel() uses a bit of algebra during parallel testing to eliminate all the slower and numerically unstable bits. (E.g., replacing divisions w/ multiplications; taking the cos() of the tolerance up front rather than acos() of the dot product).

The distance from a point to a line (Go to the code)

The dot product is also used to help determine the distance from a point to a line. It requires a bit of algebraic reasoning (I put a nicer version of this derivation on mathbin.net, but I don't know how long the link will last):

• The line l is defined by vectors o and v as above; the point is the p.
• If we draw a line from p to the closest point of l, it will be perpendicular to l. Call this closest point q.
• Since q is on l, there is a t value such that q = tv + o.
• Since the line from p to q is perpendicular to l, dot(v, p-q) = 0.
• Need to introduce a pair of identities: dot(a, b+c) = dot(a,b) + dot(a,c) and dot(a, sb) = (s)(dot(a,b))
• Thus, (dot(v,p) - dot(v,q)) = 0
• dot(v,p) - dot(v,(vt+o)) = 0
• dot(v,p) = dot(v,vt+o))
• dot(v,p) = dot(v,vt) + dot(v,o)
• dot(v,p) = t * dot(v,v) + dot(v,o)
• t = (dot(v,p) - dot(v,o)) / dot(v,v)

Once we have t, we can find q easily, then the distance between p and q.

Finally: Finding the overlaps between line segments (Go to the code.)

With all these preliminaries out of the way, finding out if and how two line segments a and b overlap is quite straightforward, and the implementation of LineSeg::Overlay should be easy to read. First, check that angle between a and b's direction vectors is within the angle tolerance. Then, check to see that each endpoint of line b is within the distance tolerance of line b, saving the t values along the way. Finally, if the (min_t_val..max_t_val) range overlaps (0..1), there is an overlap! (If not, b lies upon the same infinite line as a but they don't have anything shared.) Finally, the implementation constructs the overlapping segment as well as the clipped portions of a and/or b that were not overlapping.

There are a few tweaks to keep the implementation is fast and reasonably numerically stable -- it uses squared distances everywhere instead of distances to avoid sqrt(), for example, and ends up only having two division operations involved.

There is also fairly extensive unit testing for all the above code -- essential to keep this type of code from having nasty corner-case errors!

This overlay operation alone would be enough to implement a very slow version of boundary generator: Test every line segment of every polygon against every other line segment of every polygon and record the overlaps. Unfortunately its performance is unacceptable.

A final interesting side note

The representation we chose for lines, and line segments means that everything -- all of the math and all of the code -- remains correct and true, even if we are working in three (or more!) dimensions. All that needs to change is to start using three-dimensional vectors rather than two; everything else is built atop them using vector operations that exist in any number of dimensions.

Next time

Next, when I find the time, will be a description of the quadtree spatial index that Boundary Generator uses to keep things at least somewhat snappy.

Thoughts?

I see this sort of tutorial content all over the place in programming communities, but less often in the geospatial developer world. If you're finding this interesting or helpful, please comment or drop me an email so I keep doing it!

So you'd like to parse GML...

I am a big fan of standards.

However, I'd say that the Open Geospatial Consortium has a serious illness when it comes to defining standards for the industry: W3C-itis. The scope of the GML standard started out huge and became gigantic; the result is plagued with ambiguity and multiple ways of doing things; and is only half-correct on technical issues. (The worst example that I can think of is their botching of spatial references; e.g. axis order.) The OGC has W3C-itis far worse than the W3C themselves.

In short (as far as I can tell by looking at the result and not being involved in the process): the people writing the GML standards and those informing its development were the producers of data, and they took it up as their creed to ensure the standard could represent any and all data they could imagine and do so in a way that deviates as little as possible from their current practices. This has the effect of giving the shaft to any consumer of data, who needs to implement the entire standard in order to be able to meet the ideal of understanding all data produced under it. Which is the entire point of having a data format standard in the first place -- to decouple producers and consumers.

None of this is news and has been complained about loudly elsewhere. And there are now lots of bandaids emerging, e.g. the OGC's "profiles."

What I do have to contribute that's new is the following visualization of the GML schema for a single polygon. This is a railroad diagram, pretty straightforwardly interpreted by following the tracks. (Note that the tool I used sometimes puts sections right-to-left, so follow the arrows.) You'll need to click on it and zoom to 100% for it to be even remotely comprehensible.


A few simple observations and notes:

• Everything is optional. You can have a completely empty polygon with no interior or exterior as far as GML is concerned.
• Yes, there are five different ways to represent the points that make up a LineString for the boundary of a polygon's interior. Oh, and five different ways to represent a polygon boundary itself: as a LinearRing, or as a Ring with Curve, a CompositeCurve, an OrientableCurve or just a plain ol' LineString curveMembers. Nice.
  
• The solid blue boxes are nonterminals ('complex types' in the XML Schema lingo): there are many more tags and choices hidden under them. I couldn't fully expand them without making a ridiculously large image.
• I left off all XML attributes, in particular the grab-bag of attributes that sits on most every element. (srsName, arcRole, ref, etc)
• The tools I used to create the image aren't perfect, so don't use this as a reference to GML. Instead think of it as a nice demonstration of the complexity you face.
  

I produced the image using ebnf2ps, which is awesome; although this is a bit of an abuse of its powers. I used an unfinished homebrew tool, which at some point may be released as open source, to process the schema's XSD files into EBNF. (Its actual purpose is to produce RELAX-NG schemas, which are very close cousins to EBNF.)

Programming functionally: A proto-XPath implementation in Javascript

This has been sitting on my hard drive for a while. I thought I'd clean it up and post it! Well ... at least I've posted it.

---

If you don't know what "functional programming" is, wikipedia has a fairly technical overview. Loosely, "object" is to "object-oriented programming" as "function" is to "functional programming".

I wrote this to give a taste of the kind of power functional programming can bring to bear on a problem, as well as to document the basics of my psuedo-xpath implementation for myself. Some love to hate Javascript, but I think this comes more from the hostile browser environment than the language, which is quite nice. Note that, while I'm heavily using Javascript's functional features, I'm not making any effort to be 'purely functional'. (If I were using a purist style, I would not be using variable assignments and avoiding for/while loops, instead using recursion and if statements to the same effect.)

JS makes it quite easy to pass functions to as arguments, one of the core distinctions of programming in a functional style:

function doSomethingWithPi(arg) {
  arg(3.1415);
}
function eat(x) {
  alert(x + " sure was yummy!");
}
doSomethingWithPi(eat); // ------> "3.1415 sure was yummy!"

Simple enough. Now consider:

function siblingAxis(f,node) {
  var v = node.nextSibling;
  while (v != null) {
    f(v);
    v = v.nextSibling;
  }
}

This function, when given a function & a DOM node, calls that function for every sibling of that node (until out of siblings). If you're familiar with the XPath spec, you already know why the function is called "siblingAxis". Here's another, very similar function -- it calls the provided function once for each child of the given node:

function childAxis(f,node) {
  var v = node.firstChild;
  if (v != null) {
    f(v); // call the function for the first child node
    siblingAxis(f,v);  // do the same for all its siblings...
  }
}

And here's one that takes two functions; the first acts as a test for the given DOM node and the second is only called if the test succeeds with the given node.

function test(test_f, f, node) {
  if (test_f(node)) {
    f(node);
  }
}

Now, if we could wire our test function up to our childAxis function, we would have the beginnings of something fairly powerful. The trick is to curry the function test. It's easier to show what this means rather than try to explain it.

function test(test_f, f) {
  return function(node) {
    if (test_f(node)) {
      f(node);
    }
  }
}

Now, test(test_f,action_f) does not perform its test or action right away. Instead it returns a new function, which expects a node and in turn runs the test against the node, and the other function if it's successful. Javascript implicitly keeps test_f and action_f arguments intact and available in the returned function. This is called a closure in FP jargon. By currying test, we can wire it up with e.g., childAxis:

function isLink(node) {
  return (node.nodeType == 1 && node.localName.toLowerCase() == "a");
}

function doSomethingWithAllLinkChildren(f, node) {
  childAxis(test(isLink,f), node); // easy peasy
}

It only makes sense to do the exact same thing with the axis functions themselves. This means that each axis function takes a "do something with nodes on this axis" function, and returns a function that, from a starting node, actually traverses the given axis and does the something.

function siblingAxis(f) {
  return function(node) {
    var v = node.nextSibling;
    while (v != null) {
      f(v);
      v = v.nextSibling;
    }
  }
}

function childAxis(f) {
  return function(node) {
    var v = node.firstChild;
    if (v != null) {
      f(v); // call the function for the first child node
      siblingAxis(f)(v);  // siblingAxis is curried: siblingAxis(f,v) becomes siblingAxis(f)(v)
    }
  }
}

What have we gained from this? Well, here's a function findem that alerts with the children of each span found within a div child of the given node. In xpath terms, we're after "/div/span/*" (with the starting node treated as the root):

function nameTest(nm) {
  return function(node) {
    return (node.nodeType == 1 && node.localName.toLowerCase() == nm.toLowerCase());
  }
}

function nodealert(node) {
  alert("found " + node.localName);
}

function findem(node) {
  var path = childAxis(test(nameTest("div"),(childAxis(test(nameTest("span"),childAxis(nodealert)))));
  path(node);
}

Let's make another small change to test, currying it again:

function test(test_f) {
  return function(f) {
    return function(node) {
      if (test_f(node)) {
      f(node);
    }
  }
}
}

Now, we can use a handy bit of non-javascript notation and say that f ⋅ g is shorthand for function(x) { f(g(x)) }, and rewrite findem in this style:

function findem(node) {
  var path = childAxis⋅test(nameTest("div"))⋅childAxis⋅test(nameTest("span"))⋅childAxis⋅nodealert;
  path(node);
}

Observe the wonderful correspondence between our functions and the various components of the path they represent ("/div/span/*"):

/     ---> childAxis
div   ---> test(nameTest("div"))
/     ---> childAxis
span  ---> test(nameTest("span"))
/     ---> childAxis
*     ---> (n/a)

At this point, having XPath (or something like it) up and running in javascript is a very straightforward job.

• Each axis (child, sibling, descendant, etc) can be represented by : function someAxis(f) { return function(node) { /* traverse the axis from the node, calling 'f' along the way */ } }


• Each test (both "node tests" and other predicates) can represented as a function that takes a node and returns a boolean, then gets wrapped by test. We'd need to write helper functions such as nameTest above to construct the actual test functions.


• The input path is parsed and translated into a list of axis and test functions. Note that all of them have the form function(f) { return function(node) { /*...*/ } }


• Combine the list of functions into the form a(b(c(d(ff))) where ff is the 'action' function that you want to be called for any nodes returned by the path. This has the effect of calling the outer function(f) { /*...*/} for each node step with f as the next node step, and results in a function awaiting a node...


• Call it with the root node; your ff will be called for every succeeding node!

Click here to see the complete example of the proto-XPath implementation that I've been using along the way. It adds parent,self, and descendant axes; it also insert an implicit test at each step of evaluation to ensure that each node is visited only once. (This prevents repeat visits during paths like "//div//div"). The parsing code is hokey but works.

var DomNav = new Object();

DomNav.siblingAxis = function(f) {
  return function(node) {
    var v = node.nextSibling;
    while (v != null) {
      f(v);
      v = v.nextSibling;
    }
  }
}

DomNav.childAxis = function(f) {
  var n = DomNav.siblingAxis(f);
  return function(node) {
    var v = node.firstChild;
    if (v != null) {
      f(v);
      n(v);
    }
  }
}

DomNav.descendantAxis = function(f) {
  return function(node) {
    var v = node.firstChild;
    while (v != null) {
      f(v);
      var nv = v.firstChild || v.nextSibling;
      while (nv == null) {
        v = v.parentNode;
        if (v == node)
          return;
        nv = v.nextSibling;
      }
      v = nv;
    }
  }
}

DomNav.parentElementAxis = function(f) {
  return function(node) {
    var v = node.parentNode;
    while (v != null && v.nodeType != 1) {
      v = v.parentNode;
    }
    if (v != null)
      f(v);
  }
}

DomNav.selfAxis = function(f) {
  return f;
}

DomNav.test = function(test_f) {
  return function(f) {
    return function(node) {
      if (test_f(node)) {
        f(node);
      }
    }
  }
}

DomNav.nodeNameTest = function(name) {
  return function(node) {
     return (node.nodeType == 1 && node.localName.toLowerCase() == name.toLowerCase());
  }
}

DomNav.nodeIdTest = function(id) {
  return function(node) {
    return (node.nodeType == 1 && node.id == id);
  }
}

DomNav.seenTest = function() {
  // Note: a hash would be faster/correct but dom nodes don't work as dictionary keys. 
  // (it would also work to set a special "seen" property on each dom node itself)
  var seen = [];
  return function(node) {
    if (seen.indexOf(node) >= 0)
      return false;
    seen.push(node);
    return true;
  }
}

DomNav.parse = function(t,f) {
  if (t == null || t == "")
    return function(node) { return; }

  var r = t.split(/(\/\/|\/|\.\.)/); // split by "//", "/", ".."
  var ns = [];
  for (var i = 0; i < r.length; i++) {
  }
    c =" DomNav.test(DomNav.seenTest())" i =" ns.length">= 0; i--) {
    c = DomNav.test(DomNav.seenTest()) (ns[i](c));
  }

  return c;
}

Notes:

• I'm using the term "currying" a bit loosely to refer both to the transformation of a function f(a,b,c) { /* ... */ } into f(a) { return f(b) { return f(c) { /*...*/ } } }, and to calling a function written in this style with partial arguments to return a new function, as in f(a_argument)(b_argument).



• The code I've written is not purely functional, nor even close to it: However, it makes good use of some functional programming features, and is a nice functional/imperative compromise for javascript.



• There is still quite a bit of work to get a full implementation of XPath up and running, as XPath is a large language. However, I believe the remaing work -- if you were interested in a full implementation of XPath -- consists mostly of straightforward extension of the sketch laid out here.