voronoi issue

[jianchunhuang]

wo,

Following your suggestion, I found out the problem. I did use options.ConformingDelaunay = True. however I also smoothed the mesh. After I removed the smooth process, I got the correct voronoi diagram.

Looks like that all these infinite loop of voronoi faces are in the faces with circumcenter outside of the domain. Does it mean that HandleCase2 subroutine in BoundedVoronoi has a bug? If I use HandleCase1(edge, v1, v2) for both dir<=0 and dir>0 I got the following diagram.

Mesh is attached in the end of this discussion.

Thank you so much for your help.

Victor

P.S. Ele file 60 3 0 0 0 16 1 1 16 0 15 2 1 18 2 3 18 1 16 4 16 17 18 5 15 14 17 6 20 17 23 7 20 18 17 8 24 23 17 9 25 19 20 10 20 23 25 11 17 16 15 12 18 20 22 13 2 22 3 14 2 18 22 15 3 28 4 16 28 3 22 17 4 28 5 18 22 19 28 19 19 36 29 20 22 20 19 21 5 29 6 22 5 28 29 23 27 6 29 24 27 29 36 25 29 28 19 26 19 25 36 27 14 24 17 28 24 13 31 29 13 24 14 30 21 37 25 31 25 23 21 32 31 21 24 33 33 21 31 34 13 30 31 35 21 23 24 36 30 12 11 37 12 30 13 38 33 11 10 39 33 31 11 40 37 21 33 41 11 31 30 42 26 33 10 43 34 35 38 44 27 36 35 45 27 34 7 46 34 27 35 47 7 34 32 48 27 7 6 49 35 37 26 50 35 25 37 51 38 26 9 52 26 38 35 53 32 9 8 54 9 32 38 55 32 8 7 56 26 10 9 57 38 32 34 58 37 33 26 59 35 36 25

Poly file 39 2 0 1 0 0 0 1 1 0 2.5 1 2 0 5 1 3 0 7.5 1 4 0 10 1 5 2.5 10 1 6 5 10 1 7 7.5 10 1 8 10 10 1 9 10 7.5 1 10 10 5 1 11 10 2.5 1 12 10 0 1 13 7.5 0 1 14 5 0 1 15 2.5 0 1 16 1.51861710066147 1.48499150296617 0 17 3.56297273334277 1.83491922832983 0 18 1.67057105736466 3.51907249206639 0 19 3.2793394746019 6.48331740223419 0 20 3.19059981615234 4.38058385471337 0 21 6.48879234800104 3.55014378244063 0 22 1.56286531666954 5.9099858328593 0 23 4.73424030064428 3.19990159051327 0 24 6.04333884466454 1.57501028782374 0 25 5.00137615979257 5.14466239115557 0 26 8.55019125365695 5.85133438478974 0 27 5.79648073566084 8.64052238997417 0 28 1.73029750210372 8.17265903863546 0 29 3.74163838264543 8.48343553638656 0 30 8.6847111437244 1.14254889518071 0 31 7.86887243569513 2.50728128405131 0 32 8.95829235800975 8.9663320711651 0 33 8.57914874589765 4.05379433283979 0 34 7.35531613365615 8.36417098203868 0 35 6.69904322370934 6.79279909176506 0 36 4.99761592056865 6.98587498139638 0 37 6.87432030993352 5.16074406499738 0 38 8.49764585032612 7.30159943893413 0 16 1 0 0 1 1 1 1 2 1 2 2 3 1 3 3 4 1 4 4 5 1 5 5 6 1 6 6 7 1 7 7 8 1 8 8 9 1 9 9 10 1 10 10 11 1 11 11 12 1 12 12 13 1 13 13 14 1 14 14 15 1 15 15 0 1 0

[wo80]

Well, the thing about DCELs is: it's so easy to miss topology changes.

I found two glitches in HandleCase2:

1. The newly generated vertex has no leaving edge ... D'oh!
2. If edge.twin.face was closed before, the next pointer wasn't updated.

Please replace HandleCase2 with the following code, and let me know, if it fixes the problem! Leave this issue open until the code gets merged.

private void HandleCase2(HalfEdge edge, TVertex v1, TVertex v2)
{
    // The vertices of the infinite edge.
    var p1 = (Point)edge.origin;
    var p2 = (Point)edge.twin.origin;

    // The two edges leaving p1, pointing into the mesh.
    var e1 = edge.twin.next;
    var e2 = e1.twin.next;

    // Check if the neighboring cell was closed before.
    if (edge.twin.id != edge.twin.face.edge.id)
    {
        edge.twin.face.edge.next = e1;
    }
    else
    {
        // If the cell isn't closed yet, make sure to update the faces edge pointer.
        e1.face.edge = e1;
    }

    // Find the two intersections with boundary edge.
    IntersectionHelper.IntersectSegments(v1, v2, e1.origin, e1.twin.origin, ref p2);
    IntersectionHelper.IntersectSegments(v1, v2, e2.origin, e2.twin.origin, ref p1);

    // The infinite edge will now lie on the boundary. Update pointers:
    e1.twin.next = edge.twin;
    edge.twin.next = e2;
    edge.twin.face = e2.face;

    e1.origin = edge.twin.origin;

    // Dissolve edge from other edges (origin and face stay the same).
    edge.twin.twin = null;
    edge.twin = null;

    // Close the cell.
    var gen = factory.CreateVertex(v1.x, v1.y);
    var he = factory.CreateHalfEdge(gen, edge.face);

    gen.leaving = he;

    edge.next = he;
    he.next = edge.face.edge;
    e2.twin.next = edge;

    // Let the face edge point to the edge leaving at generator.
    edge.face.edge = he;

    base.edges.Add(he);

    he.id = base.edges.Count;

    gen.id = offset++;
    base.vertices.Add(gen);
}

[jianchunhuang]

That is great. It does fix my problem.
Many many thanks. Victor