Friday, May 01, 2009

engineering a weighted straight skeleton algorithm



This post details an implementation of the Straight Skeleton in the method of Felkel ( pdf | Felkel, Obdrzalek | 1998 ), modified to construct weighted skeletons. Source code is available here. [edit: At the time I wrote the post, I didn't understand a few things about some of the strange cases that occur in the weighted straight skeleton - you can read more about them, here.]

It contributes a method for dealing with multiple co-sited "split events" in a floating point environment, as well as a resolution strategy for horizontal roof edges, this leads to a degenerate case in which the weighted straight skeleton is ambiguous (detailed here). In this post I concentrate on showing an algorithm exists, rather than finding it's speed.

The straight skeleton takes a simple polygon as input and outputs a graph connecting the input vertices as output.

As defined by A Novel Type of Skeleton for Polygons ( doi | pdf | Aichholzer, Alberts, Aurenhammer, Gärtner | 1995 ) the straight skeleton shrinks in the edges of a polygon, tracing out the lines made by the corners of the polygon. These traced lines form the skeleton.

Strictly the skeleton is a 2D construction, however here I regard as a 3D "roof" - given the floorplan of a house, the skeleton is the gutters and crest lines of a roof when viewed from above (4, figure below).

Some properties of the skeleton include:
  • Unique for a given input polygon
  • Related to the medial axis and voronoi tesselations for convex shapes
  • The skeleton of convex shapes is a directed, acyclic graph (DAG) with the input corners as leaves
  • Each face is monotonic in a direction perpendicular to it's defining edge
Overview

The algorithm uses a sweep plane - staring from the input polygon as the lowest set of points, intersections between planes are considered in height order (perpendicularly to the input polygon's plane). The planes are defined by the gradient associated with each edge of the input polygon, in the unweighed case there are all assumed to be the same value.

examining_collisions

There are collisions that occur between faces that aren't valid. These must be ignored.

dont happen

After all collisions have been accounted for, the algorithm is complete.

Loops of corners

The primary structure used is a set of corners. Each corner has pointers to:
  • the adjacent (next and previous) input edges (I assume a counter clockwise ordering of points).
  • the adjacent corners
A partially evaluated skeleton is shown in the following diagram.

skeleton_pointer_config


Via the normal graphic conventions, holes are represented by backwards (clockwise) loops of corners. There are additional structures that list all input edges, current edges and current corners. Each edge keeps track of the skeleton's output edges.

As the algorithm runs there is one loop of corner pointers for every volume that intersects the current sweep plane. If a skeleton produces a roof with two peaks, there would be two loops of corner points at the corresponding heights. Here it is important to note that an edge can then belong to two different corners. This is the cause of much grief later on.

Collisions between three edges

When the planes defined by edges collide to a point we have to modify the data structures accordingly. Often we will remove old corners, insert new corners or add output edges to the polygon. There are only so many ways that three edges can collide:

three way


In 1) two adjacent edges meet another edge in it's centre, 2) shows three consecutive edges colliding in 3) we see a variation where three consecutive edges that form a loop collide. There are two other important cases:

green_red


When there is only one plane at a point, it is just a plane(!?). When there are two, it is just an edge (above, green). These aren't registered as events and no action is taken. 1,2 & 3) covered all possible topologies except three unconnected edges colliding (red, above). The arrangement is impossible (in the simple unweighted straight skeleton case) as there must be something between the edges that also collides at that point, at which point it is a degenerate case. I'll get onto these in a bit.

Each plane defined by each input edge cannot collide with just any old other edge-plane. Eppstein (Raising Roofs, Crashing Cycles, and Playing Pool (doi | web | Eppstein, Erickson | 1999 )) gives more ideas on a fast approach, but to just understand correctness we can start by noting that when two planes collide it forms an edge. We need three planes to collide to form a point; Before three planes can collide, two of these planes must already be adjacent. (Unlike Felkel's approach I'm less interested in the "bisectors" defined by each corner's pair of edges - this quickly becomes misleading in the weighted case). Using the above list of possible intersections, we can see we should only be testing for collisions between pairs of adjacent edges and one other edge.

This digram gives an overview of the data structure additions for each three plane case (the manipulations implemented are detailed in the following section):

pointer twiddling


These three types of collisions are sufficient the calculate the straight skeleton of almost all random input polygons (see following video or image). The other cases are degenerate and are discussed in the next section.


rect26802

Collisions must occur on the faces

One way to view the algorithm, is that after the first event, we calculate a "new input polygon" at the given height and start again - so we might see a the algorithms executing on a pile of contours, one for every event, fig. 1 below.

redo


So far the collisions have been presented as infinite planes, however we are really colliding a set of three sided polygons or faces. For example the back edge (bold, 2) in the above figure changes it's topology (red, 2) after some of the events. After event b (view the above larger) the face belongs to two volumes, and hence gives two polygons to collide against. This means that for all of the three faces involved in a collision we must ensure that a collision occurs within it's bounds.

[diagram showing what happens when we don't bother with these checks]

There are two obvious places to perform this check:
  • Store a list of events between polygons - After an event, find new events and remove old events involving each edge or bisector that's changed.
  • Store a list of events between planes (planes don't have bounds). Before evaluating an event, ensure it occurs within all three faces. This is the simpler method I implemented.
One trick here (as Felkel describes) is to realize that if two adjacent edges are involved, then the collision is within the bounds for those faces. The faces are adjacent on a surface of the skeleton volume, so they intersect only along this line (and any intersection must also occur on this line). These adjacent edges form Felkel's "bisectors", and as Eppstein explains reduces this to a ray-casting problem that has very fast known solutions. Each of the three types of collisions identified involves a pair of adjacent edges colliding.

These "three sided" faces that we're colliding are sometime triangles, but sometimes they are unbounded, when the adjacent edge move away from each other. (They are triangles with infinite area!)

Floating point calculations sometimes miss collisions that we expect, topologically, to happen. However we can increase the chances of getting a topologically sound result if we increase the margin (epsilon-error range) until we catch all collisions, even if we pick up some additional errors. There are two places where we need to be careful - in checking that a collision really happens on an edge's face and in finding events that happen "at the same height" or "in the same place". I use a cylinder shaped volume to accumulate events that occur at a given height range and within a certain distance to a centre point.

Collisions between many edges


In a degenerate case, more than three planes can meet at one point. This happens in the center of regular polygons, as well as examples such as the following image (where five or eight faces meet single points). Here I present a resolution algorithm that processes pairs of adjacent corners, for any number of edges meeting at a point.

rect27677


It is important that these occasions are detected and treated with respect, as we assume that the linked-list of corners and their locations form a simple polygon (when projected down onto the input plane). This takes a little patience when working with floating point calculations. Luckily we can generalise some rules for handling point collisions, whether three, or any number of planes meet at a point.

The collision resolution technique takes a set of planes that meet at a point as input. It processes these into lists of chains of adjacent corners that reference these edges. While constructing these chains we perform the face collisions to determine which (if any) corner references the correct portion of the edge.

Chains consist of ordered lists of edges meeting at a point. If the edges are adjacent, they are in the same chain.

Remember that one edge can be referenced by more that one corner, when there are are more than one "peaks" in the output roof shape.

pairwise


In figure 1 we see the input structure - a chain of edges abc and their associated corners (dark blue circles). We now process each pair of edges that are adjacent (orange lines) as in figures 2 and 3:
  • We remove each of these corners from the master-list. (note that inter-corner pointers, green arrows, don't need to be adjusted).
  • We add in output edges to the skeleton.
Then we process the non-adjacent edge pairs as in figure 4, adding in a new input corner and finally adjusting the inter-corner pointers, (cutting out the internal corners). However non-adjacent edge pairs take a little more thought as multiple chains can collide at one point.

pairwise_intra

As the above shows, once the inter-chain edges are processed (figure 2, above) we can then move onto intra-chain collisions (figure 3). We still process corners by pairs, but this time we process the three pairs that are non-adjacent, the red, yellow and orange edges in figure 3. We order these around the collision point (in the unweighted case this can be done by the angle any corner on the chain makes with the collision point, but see later comments on the weighted case), and process the pair made by first corner in the first chain, with the last corner of the next.

For each of these pairs we add one additional corner that is adjacent to each edge. This corner is at the location of the collision, and cuts out the inter-chain corners from any edge loops.


Horizontal roof edges

When evaluating rectilinear shapes, such as the above cross-shape, we often see horizontal skeleton lines. This is another degenerate case that we often see along side the >3 planes at a point case.

We can note that before processing events at one height, there should be no adjacent edges that form a horizontal edge in the skeleton. Also that after processing there should be no such edges.

Because every horizontal skeleton edge must have an event at it's start and end, we can be sure that if we process all the events at one height together, we will have closed all horizontal edges. [edit: this assumption is wrong - see later post on ambiguous cases]

As all events co-heighted occur simultaneously there is no interaction between events and the algorithm is to be deterministic. There is no way to sort co-heighted events, so we let them occur in an arbitrary order.

After one event occurs we may be in the situation where two consecutive, adjacent edges form a horizontal edge. Because this edge has no height, there is no way to determine which events it causes can occur first:

horiz bisector barf

In the above figure, there is a collision at x, then at the same height, faces a,b,c &d collide at y. Face a should be rejected because we're outside the bisector. However as the edge forms a horizontal skeleton line, the bisector is poorly defined (green boundaries), so we would have to accept that all four faces collide, resulting in deformed output.

We do not track events created by consecutive horizontal edges, safe in the knowledge that after we've dealt with all events at the same height there won't be be any left. In a similar way, we can't ascertain whether collisions between planes occur within the current face boundaries.

For this reason all face-collision checks for events at one height are performed before any are evaluated. I do this by caching a list of chains. This brings with it the problem that the corner used in a chain may have been superseded, or "cut off" from it's intended collision (see figure below) by an earlier event at the same height.

The resolution strategy here depends on the chain type:
  • If the chain has more than one corner in it, it represents the intersection of two faces, so there will always be a corner lying on this intersection whose adjacencies will lead to the correct corner.
  • If the chain is of length 1, we must resort to more geometric calculations to project all referenced corners onto the input edge, and choose the preceding corner.
validate chains

In the above figure 1, faces a,b,c,d and e collide at point x. We find their positions in the given volumes by locating the corresponding corners. As co-heighted events occur in an arbitrary order, by the time figure 2 occurs the corners (u & v) are no longer valid.

In the case of u (or any chain with 2+ edges), we can can use the following corner's previous corner to find u'. This is because the corner references two edges used in a collision, so must take part it (and horizontal events are independent - they can't effect each other's edges or validity).

This trick won't work to find v' as it is part of a chain of length 1. We project all corners that reference the edge onto the collision line. We then use the point (v') that is just before the projected collision as the previous corner. We ignore any points that reference the edge before the found corner. This is valid as all these points lie on the edge's plane, and all the points occur at the same height, so their order is defined by their projection. There is no question that the event occurs as there is no interaction between co-heighted events.

Negative edge gradients

By changing the weight or gradient of an edge, we can drastically change the output shape. There is even the possibility of giving the edges a negative (outwards) gradient, to create "overhangs".

No modifications are required to the given algorithm to allow these gradients.

It is fun to note that negating the gradient (reversed polarity result in the figure below) produces the same shape as reversing the polarity (loop direction) of the input, but with a different final meaning - no longer is the result the top of the roof-shape (a surface whose normals point upwards), it becomes the ceiling of the building from the inside (a surface whose normals point downwards).

reversed

When edges have negative gradients, it becomes much harder to guarantee termination of the skeleton - it may define an infinite volume! Sometimes it will terminate, others not. Determining which polygons will terminate and which won't would be a good future research direction. Of course it is always possible to specify an arbitrarily high horizontal plane, which terminates the algorithm.

The following animation places several weighted straight skeletons with varying gradients on top of each other to create some primitive architecture.



Determining the order of chains

Determining the order of corners within a chain can be accomplished by following adjacency pointers to other corners. But between chains it becomes harder. In the weighted skeleton case we must project the corners to the collision height, before determining the order.


chain_order

Angles around the intersection point are less than a full circle. If they weren't they would have overlapped collided at a lower height.

Adjacent parallel edges

If two skeleton edges are parallel (here I mean edges going in the same direction, not opposites), the intersection of their planes with a third is no longer a point (it is either a line or not an intersection at all). This leaves the straight skeleton poorly defined. This condition can occur in two places
  1. the input polygon
  2. when intermediate edges are removed to being two parallel edges together
This also causes much "instability" or chaotic behavior in the algorithm. In the unweighted case it is possible to define the collision of two consecutive parallel edge-planes as the gradient vector starting from the intersection point. This is not possible if the edges have different weights, see the following video as the weight changes:

There has to be an artificial resolution method for parallel consecutive edges. The simplest is to remove one of the offending edges - this choice can be arbitrary and is justified by the limiting case (as we approach Pi degrees from either side, we resolve to one edge or the other). This works in the first case, parallel consecutives in the input polygon.

two become one

In the above figure edge a and b are parallel and have the same weight, they are not consecutive at the start, but become consecutive after an event. Theoretically any line that lies on a and b's face could become the bisector after the collision, however it must be consistently chosen throughout the algorithm. As I'm interested in clean 3D solids, the outcome in the figure is most desirable. By merging the two faces together we can see that the input to the algorithm shouldn't have been the two separate edges a and b, but rather one edge referenced twice. This is the approach I take.

On the subject of chaotic elements of the straight skeleton, another situation deserves to be pointed out. If concave faces' bisectors pass each other ("split events" in Felkel's terminology), there is often a sudden shift in the topology of the skeleton.



Adjacent and horizontal faces form an ambiguous case

When the skeleton is weighted (different gradients on each edge), there is an ill defined condition:

ambiguous

In figure 1, we see the input, a weighted skeleton where edge a has a lower gradient than edge b. These two edges meet to form a horizontal edge. At this point there are two resolution possibilities that make topological sense, 2a, letting the lower gradient take priority and 2b, letting the higher gradient take priority.

[edit: this is wrong, either 2a or 2b, above, are equally valid] A sensible implementation should choose to give the steeper edge priority (figure 2a) for two reasons:
  1. If we look at the limiting case as edges a and b become parallel, as we reduce the angle between them, we get output topology favoring the steeper edge.
  2. Topology 2b leads to the situation where one horizontal event causes another, breaking earlier assumptions in the implementation. This would mean that we could no longer order events by height.
However this is not consistent with the concept of the skeleton being a "wavefront".

Aichholzer et al. do not describe the weighted skeleton, and so do not come across this situation in their text. I'm not aware of a description of the weighted skeleton that avoids this ambiguity.

Building faces

Once the edges of the skeleton are known, it is desirable to be able to define faces (we might want to tile our roof). There is one edge for each face.

The implemented algorithm keeps a single graph (a hash-map) of all edges in the skeleton, and for each input edge there is a list of corners that once referenced that edge. This information is sufficient to re-construct each face:

edges_one

As there may be parallel edges with the same gradient, the algorithm needs to be a little smarter:

edges_two

15 comments:

  1. I admire you for having detailed information and posting informative post. You are truly a master in this craft.

    ReplyDelete
  2. Anonymous09:19

    Hi,

    I'm trying to implement Felkel's algorithm and I don't understand how it is supposed to handle this case :

    http://www.freeimagehosting.net/image.php?8d9fd46df5.png

    During the initialization, vertices A,B,C,D create intersection X and vertices B,C,D,E create intersection Y.

    If X is the nearest (not very accurate on the image...), a new vertex is created and a new intersection is computed : Z.

    How am I supposed to know that Y shouldn't be picked up at the next iteration ? It's at the top of the priority queue.

    The only possibility I see is that ONE of its parents (C) is unactive, but the algorithm says that both C and D should be unactive for Y to be discarded.

    Thanks in advance !

    ReplyDelete
  3. So basically the edge BC can't be involved in collision at BCDE because it it removed in the collision ABCD.

    If X happens before Y (I think you mean that), the next collision will be the edges (ab), (cd) and (de) colliding.

    ReplyDelete
  4. Anonymous11:46

    Makes sense. Correct me if I'm wrong, but there is nowhere in the paper's algorithm where such a check is performed ?!

    Thanks for the quick answer anyway :-)

    ReplyDelete
  5. It'll be to do with automatic removal from the LAV/SLAV and the check that the point is still in one of these.

    ReplyDelete
  6. Anonymous07:34

    Hey Twak,

    thanks for the great work and time you obviously invested in describing your algorithm. I´m also working on an implementation of Felk`s basic algorithm, but I´m a little unsatisfied with it, especially as I´d like to change my system to a weighted solution. This is why I´m really interested in your solution. Though, as I´ve dealed quite a lot with Felk and Eppstein, I´m always searching for "Edge", "Split"- or "Vertex"-events and I don´t quite get, where these fit in your algorithm. Could you give me a short step-by-step-description what to do in every iteration? And am I getting this right, that you´re only dealing with planes and their intersections in your method? This seems to be much less error-prone than calculating bisectors, rotated-planes and so on...

    Thanks in advance for your help,
    patrick

    ReplyDelete
  7. This comment has been removed by the author.

    ReplyDelete
  8. Do you have an updated full functioning version of your code at:
    http://campskeleton.googlecode.com/svn/trunk/ ??

    I'm not succeeding in compiling it with Eclipse Juno, which outputs several warnings regarding generic types and generic methods. I would really like see it working. Indeed, I would like to get a binary .jar from it, so that I could call it from within Matlab.

    I appreciate any help.

    Regards,

    JLeandro - www.ime.usp.br - Brazil

    ReplyDelete
  9. Hi, I was wondering if you have an example of an implementation of the offset in a weighted straight skeleton.

    ReplyDelete
    Replies
    1. Hi Yushan, the code is over here https://github.com/twak/campskeleton

      Delete
    2. For the Offsetskeleton class there is a function called shrink which creates a consistent offset for all edges of the shape. I was wondering if there's a way that I can offset the edges at different variables. What if I only wanted to shrink one edge by 5, and another edge by 1, and another edge by 2? I tried changing the angles to change the skeleton shape but I've been looking through the code and I'm not sure how exactly the offset of the straight skeleton would work in this case.

      Delete
    3. You can use the example here: https://github.com/twak/campskeleton/blob/wiki/PageName.md

      But create a different Machine with a different angle (in radians) for each edge.

      Can you ask future questions on github? Thanks!

      Delete
    4. Hi, sorry I'm not sure where on github I can ask questions.
      But basically I want to get the points for the new offset shape in a weighted skeleton. Is there an example of doing getting it, or how would I go about doing it?

      Delete
  10. I have recently implemented the straight skeleton algorithm.
    Thank you very much for your description of how to deal with singular cases!

    As the first step, I implemented the typical algorithm for nonsingular cases, with edge events and split events. In some scientific paper I saw a comment saying that almost all singular cases can be solved by handling events one-by-one. However, no matter how hard I tried, my algorithm could not solve all the test cases I created by hand.

    Then I implemented handling of all simultaneous events at once by using these "edge chains", and now all the cases work perfectly!

    ReplyDelete