Determining polygon adjacency

[Judas1012]

Good day, gurus... I'd like to pick some brains for a moment and see if anyone has a better algorithm than the one I've come up with.

In a nutshell, I'm trying to determine whether two polygons are adjacent to each other -- meaning that they share an edge, including situations where the shared edge is just slightly "off" by some threshhold.

Background: I'm working with pretty much a "standard" definition of a polygon.... just a simple list of points, and the edges of the polygon are just "
adjacent" points in the list of points... so, for example, a square would be something like:

0,0,0,1,1,1,1,0

The first point is (0,0) with an edge to (0,1). The next edge is from (0,1) to (1,1), next edge is from (1,1) to (1,0), and finally, the last edge is from (1,0) back to the original (0,0) point. So, it's a square "drawn clockwise". Now, the polygons' points I'm working with are defined by latitude and longitude, so I'm using floating points to represent them.

At any rate, here's the pseudo-algorithm I've come up with:

1. Compute slopes of edges of polygon A
2. Compare slopes of edges of polygon A to slopes of edges of polygon B
3. If two slopes are similar, compare Y-intercepts of the edge in A to the edge in B
4. If the Y-intercepts are similar, then we know that not only is edge A parallel to edge B, but both edges also lie on the same line (because of the similar Y-intercepts)
5. Next, we need to see if the lines actually overlap, which would make the two polygons "adjacent." To do this, I simply use a little point-length test that takes 3 arguments -- the endpoints of one edge in one polygon and a single point from the edge in polygon B as determined by 1-4. If the length from A1 to A2 is equal to the length from A1 to B1 + B1 to A2, then the point is on the line.
6. I then repeat this for all four points, for a total of four "point-on-line" tests... from A1 to A2 testing B1, from A1 to A2 testing B2, from B1 to B2 testing A1 and from B1 to B2 testing A2. (A* are endpoints of an edge from A, likewise for B* points)

Is there a standard way to test for 2D polygon adjacency? Mind you, an adjacent polygon doesn't necessarily share an exact edge -- the edge in B may be much shorter than the edge in A I'm testing against, kind of like the following diagrams (all considered "adjacent" if you imagine the lines overlapping each other):

A: o------------------------------------o
B:                     o-------------o

A: o--------------------------------o
B:                             o------------------------------o

A: o---------------------------------o
B: o---------------------------------o

A:                   o--------o
B: o---------------------------------------o

This scenario would not be considered adjacent:

A: o----------------------------o
B:                              o---------------------------o

...because they only share a single point... they actually need to share a portion of a line segment, however small, to be considered adjacent.

Any ideas?

[JohnnyCoder]

Hello-

Dont know if there are any standard 2D adjacency algorithms... Your algorithm seems pretty sound, and the only thing I noticed (which you may recognize already) is you need to be careful here with lines with infinite slopes (ie, vertical lines). You would need a special test for those, since the slope would be a divide by zero and y-intercept is infinity. Also, I dont know how much precision counts in your application, but certainly a loss of precision will occur as your lines approach vertical.

Just a thought...

[Judas1012]

Hello-

Dont know if there are any standard 2D adjacency algorithms... Your algorithm seems pretty sound, and the only thing I noticed (which you may recognize already) is you need to be careful here with lines with infinite slopes (ie, vertical lines). You would need a special test for those, since the slope would be a divide by zero and y-intercept is infinity. Also, I dont know how much precision counts in your application, but certainly a loss of precision will occur as your lines approach vertical.

Just a thought...

Yep, I took care of the vertical line problem by assigning a "hival" value to any slope that is 0... so, the slope of a vertical line is 2,000,000,000 -- an improbable slope. I know this isn't the greatest programming style, but it works for this scenario -- the only time we use slope is to compare two lines and to compute the y-intercept... of course, hival will throw off the y-intercept computation as well (since there is no y-intercept for a vertical line), so there's another test for that.

Precision problems are what we're running into. We're using blockgroup data from the US government, and our problems are when we're computing adjacency of blockgroups (the polygons) across state lines -- for example, a polygon that is in Mississippi is not reported as adjacent to a polygon in Alabama, even though by visual inspection they're as adjacent as adjacent can get. Within a state is not a problem -- it works perfectly.

I'm thinking this is because blockgroup boundaries are managed individually by each state, and that one state's definition of a blockgroup bordering a state line is just slightly different from the adjacent state's definition, although by both visual inspection and inspection of the actual points/edges the algorithm should work.

I'm making some modifications to the algorithm, we'll see what turns up... thanks for the help, and if anyone else has any suggestions/algorithms, I'd love to hear them!

[RoboTact]

Sort all edges of polygons. In sorting function taking couples of edges, first rule is by direction, and second rule (apllied if directions are the same) is by shift. This way you'll get sets of edges lying on the same line in O(n*log(n)). In such set you can make task one-dimensional and find overlappings in linear (of the number of overlappings) time - just sort ends of line segments and support set of segments owning current position while moving along the line). Resulted complexity is O(n*log(n)+m), n - number of edges, m -number of couples of adjecent edjes. To prevent problems with extremal angles of edges don't divide enything - use vector and scalar products and compare things according to them.

[Judas1012]

Wow -- that worked absolutely great with a few tweaks. With a little pre-processing (computing all slopes of all edges of all polygons of blockgroups in the US took a little over 15 minutes), now we have an indexed table that can be seeked in no time flat. Simply a matter of comparing faces of polygons with sorted (indexed) polygon-slope list.

Time went from about 20 minutes to process one "pair" of states down to about 1 minute flat.

Thanks, RoboTact -- very appreciated.

[RoboTact]

Then please write your limits (number of edges, number of adjacent edges). Because number of adjacent edges is very probably negligible, so complexity is O(n*log(n)), which would work for about 1 second with number of edges up to ~50000. For 1 minute it would take about 5000000 (five millions) edges. And more - depending on your 'real-world' polygons, process can be optimizet to near linear complexity, if you would cover set of edges by certainly non-adjacent sets.

[Judas1012]

In the United States, there are approximately 212,000 different blockgroups, each with between 4 and 500 edges (depending on the complexity/population of the blockgroup).

The slowdown, I believe, comes from the fact that we are using the standard "line test" to determine whether two similarly-sloped edges are in fact adjacent to each other, and not just parallel. We use the distance formula ( sqrt(x^2 + y^2) ) four times: we test one point from one line with the endpoints of the other line, checking to make sure that the distance from endpoint 1 to endpoint 2 is the same as the distance from endpoint 1 to the "test" point + the distance from the "test" point to endpoint 2. Then we repeat 3 more times with different combinations of the points. If we end up with two or more cases where the lengths match, then the lines are adjacent. I think it's here that the code bogs down, since it's using the square root function of floating point numbers (latitudes and longitudes).

We tried normalizing the lines to the x- or y-axis depending on their slope, and then just compared distances of x values (getting around using the square root function) but this proved to be unreliable.

We're working on it, and have it running extremely fast now -- we just got some accuracy issues to work out. Any more suggestions on speeding this up would be greatly appreciated!

[RoboTact]

So, you need to consider the lattice and group edges by belonging to certain cells. Run the algorithm on edges within a numbers of cells, overlapping so that you'll not lose adgecent groups. In can save some log(n) multiplier, adding a portion of edges in double considering.

[Judas1012]

Just a follow-up to our final algorithm.

We found major problems using the slope of the line to see if two lines were parallel -- first was considering infinite slopes, and the next took a bit of math to figure out. A 45-degree line has a slope of 1. That means that everything less than 45-degrees is less than one, and everything above is more than 1. We found that comparison of slopes wasn't effective because we couldn't normalize the "difference" in slopes using the slope alone... for example, a slope of 450 and a slope of 440 may be "similar enough" to consider for our algorithm (a difference of 10), but a difference of 10 for a slope less than 1 just doesn't work.

We ended up "normalizing" our coordinate plane to the southernmost point of the line, then measuring the angle of the line from the x-axis. If a slope is positive, it has an angle from 0 to 90. If negative, 90 to 180. Luckily, the polygon "definitions" we are using are standardized so that every polygon is "drawn" clockwise, eliminating the need to take care of the "180 is the same as 0" special-case.

We ditched the y-intercept method altogether, instead option for a "bounding box" solution. If one line falls within the "bounding box" (a "box" drawn with the two vertices of the lines occupying opposite corners of the box) of the other line, or vice-versa, we then check using the point-on-a-line-four-times method described in a previous post.

The y-intercept test was failing because our coordinates we're checking are so large (latitude/longitude of the U.S.) -- the lines whose slopes are very large have a much higher change in y-intercept than lines with very small slopes, making "normalizing" these values difficult (much like the slopes -- a variation of 10 may be small for certain slopes, ridiculously large for others).

After we determine that the angles are similar enough to consider, and after making sure the "bounding box" conditions are good enough to consider, we then pass the two lines through the "point-on-a-line-test-four-times", which is the exact science to the whole thing. It is a slower function, so narrowing down the field of possible adjacent lines is crucial, and something we feel we accomplished well enough considering time constraints.

Processing several thousand polygons with an average of 100 edges each takes about two minutes. Not too bad. Processing the entire United States should take about 2 days by my calculations (280,000 polygons, narrowed down by distance from each other). Not bad at all.

Step one in the ultimate plan to produce a table of adjacent blockgroups (polygons) spanning the entire United States is underway... thanks to all who helped!

[blevy787]

I realize that this thread is a bit dated and I am new to this forum. First, thanks for all the very useful foregoing discussion. Let me spec my approach and then ask a question.

Approach:
I have very irregular shaped polygons in my data set. I used the border boxes as a first-pass filter. The box of two polygons must intersect using Matlab function inpolygon. Intersection is defined by sharing an overlapping border between two boxes or have a common border. Then, for the subset of polygons with border box overlap/intersection, I test for polygon point intersection by inspecting the subset of polygons with inpolygon. This gives me the 'neighboring' polygon list for each polygon entry. I fix any points that are inadvertently inside another polygon. Last step is to drop any non-reachable polygons (no neighbors) for a directed graph.

Question: any suggestions for fixing the case there there are two polygons that have intersecting/adjacent bounding boxes BUT which have some finite separation between the polygon points s.t. the x,y points of the ith polygon have no intersection with the points of the jth polygon? This is can happen with lousy map quality.

Thanks.

BSL