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]);
[download]

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
[download]

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...)