http://www.perlmonks.org?node_id=695003

Monks! Want to see some pretty colors? Check out a demo of my new module, Math::Geometry::Voronoi:

http://sam.tregar.com/voronoi.cgi

I'm using this at my day job to make overlays for Google maps. So far it's been pretty great - producing very natural looking maps. (If that link is slow or broken try again in a bit - it's a CGI on a shared server with fairly low limits.)

Want to help me make it better? I tried (and failed) to write code to join polygons when mapping multiple points in a set. In a Voronoi diagram each point generates a polygon. If I color all the polygons from a set the same color I get a map that looks right, but the code would run a lot faster if I could draw just a few bigger polygons instead of hundreds of little ones.

The code I wrote (which I never released because it didn't work) walked through the set of polygons looking for common edges, combining polygons along each edge. That worked fine in early rounds, joining big polygons with small ones. The problem is that after a while you'll reach a case where a polygon needs to be joined with itself - that is, it has an internal edge. I couldn't find a way to reliably remove the internal edges (which may be long strings of multiple edges) without sometimes botching the job and killing the polygon.

Another useful note - the input polygons produced by voronoi mapping are convex, but the output usually won't be. So just computing the convex hull won't work.

Mathmonks, I call on you! The first one to turn in a working implementation gets to be a co-maintainer of the module on CPAN. Think of the glory! And the recriminations when the code doesn't compile on an ancient AIX machine! Forget that last one!

-sam

Replies are listed 'Best First'.
Re: Better maps with Math::Geometry::Voronoi, and a Challenge for Math Monks (holes)
by tye (Sage) on Jul 01, 2008 at 19:52 UTC

    To solve the general case requires that you allow for more than just polygons because you can end up with a topologically interesting shape, that is, a shape with one or more holes in it.

    I have usually seen this handled by having a clock-wise polygon around the outside and zero or more counter-clock-wise polygons describing the holes. And often those are simplified by having them joined by lines contained within the shape.

    So your best bet is probably to only remove a "shared" edge that appears twice in one polygon if both instances of the edge are adjacent in the list. This will remove useless internal lines (and each such removal can make a new removal possible) while preventing you from splitting your shape into two shapes (the outer edge and the hole).

    The next trick is to flag shared edges that you end up being unable to remove so that you don't draw lines for those. You fill the polygon like normal but you need to just skip over self-shared edges when you go to draw a line around the border of the polygon. Even better is to draw the "internal" (self-shared) edges using the "fill" color so you don't end up with unshaded pixels when they fall too precisely onto such a line. If you draw the internal / self-shared edges before the other edges, then that is best of all (to make sure you don't have single-pixel gaps in your border line).

    - tye        

      Hmmmm, holes. Yes, that would be a problem. I don't think Google maps can draw that kind of shape, so I'd probably have to slice it into two C-shaped polygons. Tough to do in a generic way I bet.

      Good advice otherwise though. Thanks.

      -sam

Re: Better maps with Math::Geometry::Voronoi, and a Challenge for Math Monks
by BrowserUk (Patriarch) on Jul 02, 2008 at 00:14 UTC

    I realise that Voronoi diagrams have their own (many) uses, and that you cannot talk about your work application, but...is there any way that you can describe or analogise an application for deriving a single encompassing polygon from a set of points?

    I'm having trouble envisaging the application that starts with a set of points, needs to construct a single polygon that "encompasses" them all, and for which the simplest or best mechanism is to go through the construction of a voronoi diagram and then coalesce the polygons.

    I worked on an application years ago, something to do with x-ray microscopy, that given a very large set of points needed to calculate the smallest area that encompassed all the points. That too did not want a convex hull solution. I'm trying to remember the details (and name) of the algorithm we ended up using, but I remember that it was fairly simple and very much faster than the first couple of attempts. And it definitely didn't calculate a voronoi diagram at any stage. I think it used vector math?


    Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
    "Science is about questioning the status quo. Questioning authority".
    In the absence of evidence, opinion is indistinguishable from prejudice.
      The classic Voronoi Diagram use-case is phone-booths on a campus. Let's say you need to use the phone and you want to know which phone is closest to your current location. A map that shows you this in the simplest way is a Voronoi diagram where the points are phone-booths and the polygons are the area around each phone where that phone is the closest one. Find out which region you are in and you know which phone to head for, no measuring required.

      Generally speaking I think Voronoi diagrams are useful any time you have a collection of points and you want to visualize them as geometric areas instead of clumps of points.

      My use-case isn't so different from the campus phone maps, but I don't want to get into specifics. Our competitors are everywhere!

      -sam

        I get the idea of the Voronoi diagrams. The bit that confuses me is the conversion of the cells whihc are the useful bits into a single polygon...but that said, I just thought of an application: showing (for example) cell phone coverage.


        Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
        "Science is about questioning the status quo. Questioning authority".
        In the absence of evidence, opinion is indistinguishable from prejudice.

      Given a large number of points and wanting the smallest area that encompasses all of the points but not wanting the convex hull, the answer would be a shape with an area of zero (assuming the at-least-not-clearly-stated additional requirement of a single connected area, a series of line segments can suffice).

      So I think there must be another requirement hiding in there that I'm not inferring. :)

      - tye        

        Hm. Good point. I wasn't involved in the sourcing of the data, and don't really know what it represented. The description I was given, as best I remember it, was that the points represented scintillations on the memrane that forms the cell walls of an amoeba-like thing. The creature can stretch bits of itself out into limb-like appendages, much like you see here.

        The problem is that being 3 dimensional, you get scinitallations from all over the outer membrane, not just around the 2D edges. So reconstructing a top-down, 2D representation of the thing, required connecting the dots that formed the outer edge of it, whilst ignoring those that came from the top of it. If you plot the points, your eyes/brain allow you to "see" the outline, and people would spend hours manually tracing the outines in GDDM so that they could extract the creatures image from the background before processing it further.

        Effectively it was an "edge tracing algorithm", but as the entire image was constructed of dots, and the background had plenty of them also, it was rather more complex. First you had to low-pass filter the points to exclude noise and as much of the background as possible without lossing too much definition in the creature. then attempt to draw a single line around the thing. As the points were digital on a defined grid, you end up with a complex polygon.

        Maybe that "single line" is the extra qualification you are after?


        Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
        "Science is about questioning the status quo. Questioning authority".
        In the absence of evidence, opinion is indistinguishable from prejudice.
Re: Better maps with Math::Geometry::Voronoi, and a Challenge for Math Monks
by hossman (Prior) on Jul 01, 2008 at 21:55 UTC

    Do you have any demos of using this on an actual map? (instead of just random points in a blank space)

    FWIW: more fun with Voronoi diagrams. They look really cool when you color the polygons using gradients between the points and the polygons.

    Which reminds me: Is there a way to tell Math::Geometry::Voronoi that you want the coordinate space to wrap?

      Map demos - no, I don't. So far all the work stuff is still private. I think I could post some screenshots, but I'd have to clear it up the line.

      Wrapping space - no, there isn't. I think that's the kind of thing you'd have to hack into the underlying way-math-heavy C code. Patches welcome!

      -sam

Re: Better maps with Math::Geometry::Voronoi, and a Challenge for Math Monks
by salva (Canon) on Jul 02, 2008 at 08:05 UTC
    Supposing the polygon descriptions obtained from M::G::Voronoi doesn't return different overlapping edges (that I think it doesn't):
    1. select all the points/polygons you want to merge: @points
    2. calculate the set of all the edges involved: @edges
    3. discard from @edges all the edges that have both side points belonging to @points
    4. reconstruct polygons from the remaining edges:
      1. calculate the set of all the vertices in @edges: @vertices
      2. set @polygon = (); $direction = vector (0, -1)
      3. get the leftmost $vertex from @vertices
      4. select from @edges the ones that have $vertex at some extreme: @departing
      5. for my $dep (@departing) { calculate the angle in the range (-π,π) between $direction and the vector ($dep->the_other_vertex - $vertex) } and select the edge $this_edge with minimum angle (minimum = most negative, not minimum in abs value!)
      6. push $polygon, $this_edge
      7. set $direction = previously calculated vector ($this_edge->the_other_vertex - $vertex)
      8. set $vertex = $this_edge->the_other_vertex
      9. repeat from 4.4 until $vertex is equal again to the starting vertex found in 4.3
      10. look for holes inside @polygon (left as an exercise for the reader :-) )
      11. remove the edges inside @polygon from @edges
      12. repeat from 4.1 until @edges is empty.
      Huh, sounds promising. A few questions:

      Re #5: Can you explain how I might implement this? My math skills are failing me here. Code would help.

      Re #10: Hah! If that's the case then I guess it's game over for me. Just solving the easy part isn't going to help me much. It seems to me that without this step I'm highly likely to produce invalid polygons that won't render.

      Re #11: How do you do that? I think if I could solve this one then my original approach might have worked.

      -sam

        Re #5: Can you explain how I might implement this? My math skills are failing me here. Code would help.
        You can compute the angle as ...
        use Math::Trig qw(pi); $angle = atan2($v2->{y}, $v2->{x}) - atan2($v1->{y}, $v1->{x}); # and normalize the angle to be in the range [-pi, pi] my $rounds = floor(($angle + pi) / (2 * pi)); $angle -= 2 * pi * $rounds;
        Re #10: Hah! If that's the case then I guess it's game over for me. Just solving the easy part isn't going to help me much. It seems to me that without this step I'm highly likely to produce invalid polygons that won't render.
        It's easy (at least conceptually!):

        At this point, @polygon contains a set of edges forming a closed polygon, all you have to do is to select all the edges inside from @edges and form (negative) polygons with them again. You will have to recursively call the procedure from 4.1 and take into account that this time you get polygons representing holes (that can contain holes, that are holes in holes and so "positive" polygons and so on).

        To find the set of edges inside @polygon from @edges, the easiest way is probably to recursively select all the edges connected to the edges/points in @polygon but that are not in @polygon.

        Re #11: How do you do that? I think if I could solve this one then my original approach might have worked.
        That's trivial (maybe some point on my previous post was confusing): @polygon contains a subset of @edges forming a closed polygon, you only have to ...
        my %polygon = map { $_ => 1 } @polygon, @polygon_holes; @edges = grep { !$polygon{$_} } @edges
Re: Better maps with Math::Geometry::Voronoi, and a Challenge for Math Monks (Delaunay graph)
by blokhead (Monsignor) on Jul 02, 2008 at 12:28 UTC
    The Voronoi diagram of a set of points is the dual graph of the Delaunay triangulation of the points. The Voronoi regions for p,q share an edge if and only if p,q have a Delaunay edge between them.

    This means that questions about Voronoi region adjacency can be more conveniently framed as questions about adjacency/connectivity in a plain graph, if you have both structures available as well as the correspondence between the two. In terms of your sets of points, if some Delaunay edge has its two endpoints in the same set/grouping, then you would have colored those two Voronoi regions the same, and you can throw out the dual edge (this is the edge in the Voronoi diagram that the two regions share).

    Following this Delaunay-adjacency idea through, you can compute in the Delaunay graph the "contiguous" components (there's probably a better name for this -- I just made it up). In other words, partition the Delaunay graph into maximal connected subgraphs whose vertices are all in the same set (starting at each unseen vertex, traverse with BFS/DFS only to neighbors in the same set as the root). These partitions will correspond exactly to the polygons you want. You can take each contiguous component, and traverse around its perimeter, accumulating the border Voronoi edges. If removing a particular component disconnects the Delaunay graph into k components, then the corresponding compound Voronoi polygon will have k-1 holes, so you'd need to traverse those perimeters as well (although the perimeters are shared by other components and traced out there too).

    blokhead

      This sounds correct to me. You're definitely right about Delauney graphs and you can get that data from my module by looking at which points correspond to an edge. The problem for me is that I don't know if I can actually translate this solution into code! Maybe I'll give it a try and see what I can come up with.

      -sam

Re: Better maps with Math::Geometry::Voronoi, and a Challenge for Math Monks
by moritz (Cardinal) on Jul 01, 2008 at 19:08 UTC
    I tried (and failed) to write code to join polygons when mapping multiple points in a set.

    If you want help, your best bet is to share the code you have so far, and your test cases as well. This mightily reduces the entry barrier.

      I doubt it. I think my approach was fundamentally wrong, not buggy. I don't think it can work the way I tried it.

      As far as test cases - that's a good idea. I'll see what I can do. At the very least I should be able to whip up a visual test.

      -sam

Re: Better maps with Math::Geometry::Voronoi, (Working* code)
by BrowserUk (Patriarch) on Jul 03, 2008 at 07:50 UTC

    *almost

    There are are still occasional "extra" points being included, but I believe these are being output (wrongly) from M::G::V, as with the duplicated points I noted here. And once that problem has been solved, the spurious edges that show up occasionally will go away.

    How it works

    If you look at the example (requires FireFox) that samtregar posted, and consider just the edges that form the overall outline of the shaded polygons. In the set of all the edges of all the small polygons, all the "internal edges", that is edges that are not part of the bounding polygon, will appear in two of the smaller polygons. But those edges making up the bounding polygon will only appear in one.

    By processing all the edges from all the small polygons through a filter selecting just those edges that only appear once in the total set of edges, you will be left with just those edges that make up the bounding polygon. So, the problem then becomes one of simply (!) re-constructing the boundary polygon from that disordered set of edges.

    The exclamation mark is because although the final cut of the reconstruction (reordering) process is actually quite simple, arriving at the algorithm was anything but.

    The problem on a simple scale, is that you end up with an array of arrays of arrays like this:

    [ [[177, 518], [205, 353]], [[259, 614], [177, 518]], [[205, 353], [383, 489]], [[259, 614], [380, 498]], [[367, 274], [415, 456]], [[464, 246], [367, 274]], [[380, 498], [486, 622]], [[383, 489], [415, 456]], [[452, 160], [464, 246]], [[490, 135], [452, 160]], [[486, 622], [579, 606]], [[490, 135], [586, 244]], [[579, 606], [632, 643]], [[586, 244], [717, 227]], [[632, 643], [719, 549]], [[717, 227], [746, 269]], [[719, 549], [775, 555]], [[746, 269], [856, 312]], [[775, 555], [856, 312]], ]

    And you need to re-order them so that the (right-most) end of one edge mates up with the start (left-most) end of the next. And do this all the way down the list until the end of the last edge matches the start of the first. You can then take the set of starts (or ends) and the first (or last) point from the other end, and you end up with an ordered set of points to construct the closed polygon.

    The problem is that the end of the first point could match up with any of the subsequent edges. Start or end. Having tried half a dozen combinations of: sorting and then swapping; or swapping then sorting; or sorting, swapping some and then resorting (rinse, repeat); I finally arrived at the algorithm by working through the above dataset manually and noting what I did. I then reproduced that in code. The result is the following subroutine:

    sub edges2poly { my $edges = shift; ## Ref to the AoAoA above. ## For each edge in the list (which is reordered as this loop prog +resses) for my $i ( 0 .. $#$edges - 1 ) { ## Stringyfy the end-point my $target = "@{ $edges->[ $i ][ 1 ] }"; ## And scan the remain edges for my $j ( $i + 1 .. $#$edges ) { ## Comparing the target, ## first against the start if( $target eq "@{ $edges->[ $j ][ 0 ] }" ) { ## found the next edge correctly oriented. last if $j == $i + 1; ## nowt to do if its already the + next edge } ## and then the end of each elsif( $target eq "@{ $edges->[ $j ][ 1 ] }" ) { ## Found it, but its ends need swapping @{ $edges->[ $j ] } = reverse @{ $edges->[ $j ] }; last if $j == $i + 1; ## nowt more to do if its the ne +xt edge } else { next; ## try the next edge } ## If we got here, we found it (and possible swapped then +ends, ## But it is in the wrong place in the list, so swap the e +dges. @{ $edges }[ $i+1, $j ] = @{ $edges }[ $j, $i+1 ]; last; ## And skip to the next target } } ## return an AoA of just the required points for the bounding poly +gon. return( $edges->[ 0 ][ 0 ], map{ $_->[ 1 ] } @$edges ); }

    The POC code uses GD to plot the points (in red), the Voronoi diagram (in green), and the boundary polygon (in blue). (Update: Added the pre-filtered polygons in grey to show that incomplete polygons are being removed prior to the 1-edging process).

    There are three operational modes: -PAT=[square|diamond|random]. The square and diamond options show up the duplicate vertices problem noted elsewhere, (but I filter these out) and show that the algorthm works.

    The random (default) pattern takes another option (-N=nnn) which is the number of randomly generated points passed to M::G::V to start the thing off. This demonstrates that the algorithm is very fast, but also shows up the fact that my filtering process is not succeeding in filtering all the spurious elements.

    If you run with the default (-PAT=random -N=20), then most of the time everything is fine. If you increase that to -N=100, it will sometimes be fine, but will often show the effects of the spurious points. Run with -N=1000 and it will (mostly) draw the full bounding polygon, but usually with a few "extra" edges.

    Update Added a couple of extra test modes: -PAT=hex|hex2 which more clearly demonstrate the "spurious points" problem.

    (If you run this on linux, you'll need to tweak the line:system 'voronoi.png'; to cause the program to display the resulting .png)

    The code: ('scuse the mass of debug comments)


    Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
    "Science is about questioning the status quo. Questioning authority".
    In the absence of evidence, opinion is indistinguishable from prejudice.
      Awesome, thanks! I'll give this a try and let you know how it goes.

      One question - does it deal with shapes that have holes in them? Unfortunately I think tye is right, that will happen. Perhaps it can explain some of the spurious edges? I'm skeptical of blaming them on the underlying Voronoi code - that's some very old, well tested code and it doesn't show any anomalies when you graph the results without joining. Happy to be proven wrong, of course, but I'd need proof based on the original data to believe it.

      -sam

        does it deal with shapes that have holes in them?

        I think it should, though it may need some modification to isolate the two (or more) polygons. I'll need to construct a known testcase and see what happens.

        As for the source of the errors. The presence of 'holes' could explain the anomolies I'm seeing. But either way, the problem is exacerbated by the presence of duplicate points that are being returned in individual polys. I've tried printing the values from the simple 'squares' testcase I posted earlier, to the fullest precision Perl can give me, and this isn't a case of points that differ by infinitesimal amounts being mapped to the same pixel. The values as returned bu M::G::V are all exact integer values.

        I'm currently filtering each of the individual polys for adjacent duplicate vertices, but that won't handle the case of there being 'extra points' in the solution. If you run the code above with -PAT=hex2 and look at the grey (unfiltered) polys on the left-hand edge, you'll see that there are various edges showing up that do not obviously result from the normal dividors from the points.

        In particular, if you look at the poly in the lower left corner, the extreme left-hand edge appears to be parallel to the left edge of the coordinate space (white box). There is no way you should be able to generate a 'parallel normal' from points wholy contained within the white coordinate space--but there it is.

        I'm not accusing the underlying libraries of having errors, I just don't know any way to explain the presence of that edge (and others) based on what I know of the algorithms involved. Now I'm wondering what hoops would be involved in compiling the underlying library on Win32? (Did you post a pointer to teh source code anywhere?)

        I had to make a couple of minor edits to memory.c (in myalloc() to get it to compile with MSC. Apparently cl doesn't know how to calculate sizeof( void* )? I just switched the two occurances to char* until I got a chance to investigate that, work out the correct solution and produce a patch for you. But I do not see that it could have any affect on the math.

        I'll play and try and come up with a simple dataset that produces a 'hole' to see what happens. If you have any ideas on how to produce such a dataset, please let me know. I'm having trouble conceiving of how that can arise right now?


        Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
        "Science is about questioning the status quo. Questioning authority".
        In the absence of evidence, opinion is indistinguishable from prejudice.
        One question - does it deal with shapes that have holes in them?

        You know, the more I think about this, the more I am convinced (on the basis of my own brand of logic rather than any real mathematic understanding), that it isn't possible to have a 'hole' in a Voronoi diagram?

        If someone knows better, please demonstrate?


        Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
        "Science is about questioning the status quo. Questioning authority".
        In the absence of evidence, opinion is indistinguishable from prejudice.
Re: Better maps with Math::Geometry::Voronoi, and a Challenge for Math Monks
by roboticus (Chancellor) on Jul 02, 2008 at 15:02 UTC
    samtregar:

    I've some ideas on how to add polygon creation to the code that actually creates the vertex / line / edge list, but I don't really have time to code it up at the moment. But if I had to do it in a separate routine, then here's what I'd do:

    • First, I'd scan the point list, and for each "match", I'd mark the point with the higher index as "replaced by" the other.

    • Next, edit your list of lines to replace the replaced points with the replaced-by points.

    • Finally, build your polygon edge lists by scanning through the lines list. If a line has the same point on both sides, skip it. Otherwise, add the matching edge to both polygons.

    Note: This method uses integer comparisons only, so you needn't worry about float mismatches...

    Suppose, for example, you have:

    my @points = ([1.0, 4.0], [2.3, 0.7], [3.6, 6.0], [4.0, 3.6]); my @vertices=([0.0, 0.0], [0.5, 7.0], [2.0, 2.2], [2.2, 5.3], [4.6, -0.2], [5.2, 7.0]); my @edges=( [0, 0, 1], # Left bound for polygon on P0 [1, 0, 2], # P0 | P1 [2, 1, 3], # P0 | P2 [3, 2, 3], # P0 | P3 [4, 0, 4], # Lower bound for P1 [5, 2, 4], # P1 | P3 [6, 3, 5], # P2 | P3 [7, 1, 5], # Upper bound for P2 [8, 4, 5] ); # Right bound for P3 my @lines=( [a, b, c, 0, -1], [a, b, c, 0, 1], [a, b, c, 0, 2], [a, b, c, 0, 3], [a, b, c, 1, -1], [a, b, c, 1, 3], [a, b, c, 2, 3], [a, b, c, 2, -1], [a, b, c, 3, -1]);

    Note: This is a hand-constructed dataset, which is definitely *not* a voronoi diagram. (I don't have the voronoi package installed on this machine. (In retrospect, it would've been faster for me to install & run it and build a *real* dataset... oh, well.)

    If points 0 and 3 are the same color, while the other two differ, then when you scan your edge list, you'll mark edge/line 3 as "removed", and edit the list to replace all references of point 3 with point 0, giving:

    my @points = ([1.0, 4.0], [2.3, 0.7], [3.6, 6.0], [4.0, 3.6, 0]); # Matches 0 my @edges=( [0, 0, 1], [1, 0, 2], [2, 1, 3], [3, 2, 3, 'removed'], [4, 0, 4], [5, 2, 4], [6, 3, 5], [7, 1, 5], [8, 4, 5] ); my @lines=( [a, b, c, 0, -1], [a, b, c, 0, 1], [a, b, c, 0, 2], [a, b, c, 0, 0, 3], # removed! (same on both sides) [a, b, c, 1, -1], [a, b, c, 1, 0, 3], # P3 replaced by P0 [a, b, c, 2, 0, 3], # P3 replaced by P0 [a, b, c, 2, -1], [a, b, c, 0, -1]); my @polygons=( [0, 1, 2, 5, 6, 8], # Poly 0 edges [1, 4, 5], # Poly 1 edges [2, 6, 7] ); # Poly 2 edges

    As I mentioned above, I was thinking I'd change the C code to create the polygons in the same process as computing the lines/edges/vertices. The method I was going to use was to add a column to @points for "match" determination. (I.e., your code would do the first step before constructing the diagram.) Then it would construct the vertices & edges normally, but during line construction perform the replacement on the fly. (You could have two pairs of point references, one as is computed currently and the other holding the "edited" point references, in case you wanted the original regions *and* the combined ones...)

    ...roboticus
      This sounds very similar to the approach I tried, which worked a lot of the time and then failed on perverse shapes. How does your approach deal with shapes containing holes?

      -sam

        samtregar:

        If you collect the edges together into continuous boundaries, then you'd get the expected external border. Then you'd wind up with more edges. Collecting a continuous list of those would yield an internal polygon. Rinse, lather, repeat until you exhaust the list of perimeter segments. Would that be acceptable to you, or do you need something different? If you have a particular way you want to represent them, let me know, and I'll see what I can do...

        ...roboticus
Re: Better maps with Math::Geometry::Voronoi, and a Challenge for Math Monks
by BrowserUk (Patriarch) on Jul 02, 2008 at 18:34 UTC

    One thing that may be confusing your attempts to join the polygons, is that the polygons returned from M::G::V->polygons can contain duplicate points. Are you aware of this?

    By way of example, the following creates a regular grid of 9x7 equi-spaced points and the uses M::G::V to compute and return the polys. The result is 35 (7x5) "square" polygons. But each of the polys return has 6 rather than 4 points? With 2 points being duplicated in each poly. This doesn't affect the drawing of the polys, but it gets confusing when trying to manipulate the edges.

    #! perl -slw use strict; use Data::Dump qw[ pp ]; $Data::Dump::MAX_WIDTH = 200; use Math::Geometry::Voronoi; my @points = map{ my $y = $_; map [ $_, $y ], map $_ * 100, 1 .. 9; } map $_ * 100, 1 .. 7; my $geo = Math::Geometry::Voronoi->new( points => \@points ); $geo->compute; my @geoPolys = $geo->polygons; pp \@geoPolys;

    Produces:

    C:\test>voronoi [ ## {**** duplicates ****} {**** duplicates ****} [10, [250, 250], [250, 250], [250, 150], [150, 150], [150, 150], [15 +0, 250]], [11, [350, 250], [350, 250], [350, 150], [250, 150], [250, 150], [25 +0, 250]], [12, [450, 250], [450, 250], [450, 150], [350, 150], [350, 150], [35 +0, 250]], [13, [550, 250], [550, 250], [550, 150], [450, 150], [450, 150], [45 +0, 250]], [14, [650, 250], [650, 250], [650, 150], [550, 150], [550, 150], [55 +0, 250]], [15, [750, 250], [750, 250], [750, 150], [650, 150], [650, 150], [65 +0, 250]], [16, [850, 250], [850, 250], [850, 150], [750, 150], [750, 150], [75 +0, 250]], [19, [250, 350], [250, 350], [250, 250], [150, 250], [150, 250], [15 +0, 350]], [20, [350, 350], [350, 350], [350, 250], [250, 250], [250, 250], [25 +0, 350]], [21, [450, 350], [450, 350], [450, 250], [350, 250], [350, 250], [35 +0, 350]], [22, [550, 350], [550, 350], [550, 250], [450, 250], [450, 250], [45 +0, 350]], [23, [650, 350], [650, 350], [650, 250], [550, 250], [550, 250], [55 +0, 350]], [24, [750, 350], [750, 350], [750, 250], [650, 250], [650, 250], [65 +0, 350]], [25, [850, 350], [850, 350], [850, 250], [750, 250], [750, 250], [75 +0, 350]], [28, [250, 450], [250, 450], [250, 350], [150, 350], [150, 350], [15 +0, 450]], [29, [350, 450], [350, 450], [350, 350], [250, 350], [250, 350], [25 +0, 450]], [30, [450, 450], [450, 450], [450, 350], [350, 350], [350, 350], [35 +0, 450]], [31, [550, 450], [550, 450], [550, 350], [450, 350], [450, 350], [45 +0, 450]], [32, [650, 450], [650, 450], [650, 350], [550, 350], [550, 350], [55 +0, 450]], [33, [750, 450], [750, 450], [750, 350], [650, 350], [650, 350], [65 +0, 450]], [34, [850, 450], [850, 450], [850, 350], [750, 350], [750, 350], [75 +0, 450]], [37, [250, 550], [250, 550], [250, 450], [150, 450], [150, 450], [15 +0, 550]], [38, [350, 550], [350, 550], [350, 450], [250, 450], [250, 450], [25 +0, 550]], [39, [450, 550], [450, 550], [450, 450], [350, 450], [350, 450], [35 +0, 550]], [40, [550, 550], [550, 550], [550, 450], [450, 450], [450, 450], [45 +0, 550]], [41, [650, 550], [650, 550], [650, 450], [550, 450], [550, 450], [55 +0, 550]], [42, [750, 550], [750, 550], [750, 450], [650, 450], [650, 450], [65 +0, 550]], [43, [850, 550], [850, 550], [850, 450], [750, 450], [750, 450], [75 +0, 550]], [46, [250, 650], [250, 650], [250, 550], [150, 550], [150, 550], [15 +0, 650]], [47, [350, 650], [350, 650], [350, 550], [250, 550], [250, 550], [25 +0, 650]], [48, [450, 650], [450, 650], [450, 550], [350, 550], [350, 550], [35 +0, 650]], [49, [550, 650], [550, 650], [550, 550], [450, 550], [450, 550], [45 +0, 650]], [50, [650, 650], [650, 650], [650, 550], [550, 550], [550, 550], [55 +0, 650]], [51, [750, 650], [750, 650], [750, 550], [650, 550], [650, 550], [65 +0, 650]], [52, [850, 650], [850, 650], [850, 550], [750, 550], [750, 550], [75 +0, 650]], ]

    Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
    "Science is about questioning the status quo. Questioning authority".
    In the absence of evidence, opinion is indistinguishable from prejudice.
      No I didn't know that. Should be pretty easy to filter out, I'd imagine. The underlying algorithm isn't mine, so I can't say for sure why that might happen. My guess would be very small edges collapsing to a zero-width edge. By producing squares you might very well have constructed a worst-case scenario for the algorithm.

      -sam

      BrowserUk:

      I used your code and examined the results of all the output arrays. From what I can see, it appears that the extra vertices are due to extra lines. Not extra in the sense that they're not valid, just in that they're unnecessary/unwanted. In the case of a square of four points (Pll, Pul, Plr, and Prr where u=upper, l=lower/left, r=right), you expect a line between Pll and Pul, you also expect a line between Pll and Plr. But a line between Pll and Prr is also valid. Since they intersect at a single point you don't need the last one. I suspect that the code is either experiencing just enough rounding error to insert the extra points, or that the code simply has all those three lines as boundaries, and it simply computes duplicates as it walks through the list of lines.

      ...roboticus
        . In the case of a square of four points (Pll, Pul, Plr, and Prr where u=upper, l=lower/left, r=right), you expect a line between Pll and Pul, you also expect a line between Pll and Plr. But a line between Pll and Prr is also valid.

        No, that cannot happen. A square can only arise from a cross pattern of 5 dots. And a diagonal line connecting the opposite corners would pass through the central dot and could not be 'normal' to any divisor between that dot and any other.


        Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
        "Science is about questioning the status quo. Questioning authority".
        In the absence of evidence, opinion is indistinguishable from prejudice.
Re: Better maps with Math::Geometry::Voronoi, and a Challenge for Math Monks
by toma (Vicar) on Jul 03, 2008 at 04:49 UTC
    I suggest using Triangle (offsite link) to mesh your data, then load the resultant undirected graph into the Graph module, then do whatever work needs to be done. Triangle will handle whatever challenges your data has, unless it doesn't, in which case you can clean up the data on the way in to Triangle. If you can't clean it up, you will have a research project on your hands.

    This is an example where floating point numbers can be a problem. For example, consider three lines that intersect at a single point. Due to floating point errors, this tends to generate a tiny triangle instead of a single point intersection. Issues such as this are covered in the Triangle documentation.

    It should work perfectly the first time! - toma
      Interesting - I hadn't run into Triangle before. I'll dig into the docs. But I do have one question if you have the time - what does it mean to "mesh" my data? If I pass Triangle a set of polygons that should all be considered a single shape, what would I get back?

      Luckily for me my data is on an integer grid so I think I'll be spared the problems floats can cause.

      -sam

        "Meshing" is the process of splitting up your geometry into fundamental shapes such as triangles or squares or rectangles.

        You can pass triangle a set of polygons and get back a bunch of triangles that correspond to a single shape. For example, you can make a very good polygon merge program by starting with Triangle. It handles holes, non-convex shapes, and other problems without difficulty.

        You'll have to look at Triangle to see what all it does. It is a large topic. I use it to find messed up geometries in large collections of polygons. For example, my code detects self-intersecting polygons (that is, 'bow ties'). There are about a half-dozen such checks that I do on each polygon. The polygons can have hundreds of points or more. My code typically runs one polygon at a time over a collection of about 60,000 polygons. My perl code creates Triangle input files, calls Triangle from backticks, and then processes the Triangle output file. The Graph module is handy for loading the Triangle output data.

        It should work perfectly the first time! - toma
Re: Better maps with Math::Geometry::Voronoi, and a Challenge for Math Monks
by BrowserUk (Patriarch) on Jul 05, 2008 at 06:20 UTC

    Sam. I've finally come up with an example that I believe proves that M::G::V is producing wrong output. Whether the underlying library, or something happens on the way through I cannot identify.

    If you run the following code, it produces 10 polygons, including the triangle that appears as an asymmetry at roughly 11 o'clock, in the png produced.

    However, when I feed the same points into this voronoi applet, only 9 polygons and no asymmetry result. ( You have to feed the points into that applet by clicking, so I've arranged the input data in roughly the correct (physical) pattern in the code below.)

    If you compare the results from the two implementations I think you'll agree that M::G::V is producing an extraneous polygon.

    #! perl -slw use strict; use GD; use Data::Dump qw[ pp ]; $Data::Dump::MAX_WIDTH = 200; use Math::Geometry::Voronoi; my @points = ( [150,150], [500,150], [850,150], [500,300], [429,329], [429,471], [400,400], [500,400], [600,400], [571,329], [571,471], [500,500], [150,650], [500,650], [850,650], ); my $geo = Math::Geometry::Voronoi->new( points => \@points ); $geo->compute; my @geoPolys = $geo->polygons; pp \@geoPolys; my @gdPolys = map { gdPolyFromPoints( @{ $_ }[ 1 .. $#$_ ] ); } @geoPo +lys; my $img = GD::Image->new( 1000, 800, 1 ); $img->filledRectangle( 0, 0, 1000, 800, 0x00ffffff ); my $color = 0xA0; for ( @gdPolys ) { $img->filledPolygon( $_, rgb2n( $color, $color, $color ) ); $color -= 0x10; } open IMG, '>:raw', 'voronoi-t.png' or die $!; print IMG $img->png; close IMG; system 'voronoi-t.png'; sub rgb2n{ unpack 'N', pack 'CCCC', 0x40, @_ } sub gdPolyFromPoints { return unless @_; my $p = new GD::Polygon; $p->addPt( @$_ ) for @_; return $p; }

    Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
    "Science is about questioning the status quo. Questioning authority".
    In the absence of evidence, opinion is indistinguishable from prejudice.
      Thanks, this is very useful. Have you looked at whether the extra poly is present in the lines() and edges()? That's the data that comes straight from the C implementation. If not, no problem, that's the first thing I'll look at.

      -sam

        No. Lines() and Edges() return 36 elements each. These are the 30 edges that make up the 9 polygons, plus six lines off the edge of the coordinate space, which are not part of any complete polygon.

        So the polygons() method must be manufacturing the extra poly.


        Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
        "Science is about questioning the status quo. Questioning authority".
        In the absence of evidence, opinion is indistinguishable from prejudice.