Using Triangle to create voronoi and the result is wrong

Nov 9, 2016 at 6:17 AM
I have encountered problem with using triangle to create voronoi. The results returned by the algorithm is wrong. Please contact me if you have solved this problem.

The input:
22
1 586.3 306.35
2 547.3 306.35
3 567.034 277.849
4 508.248 306.35
5 469.261 306.35
6 547.235 363.506
7 566.8 334.928
8 527.813 334.928
9 488.813 334.928
10 605.787 334.928
11 626.21 306.02
12 645.06 335.06
13 722.8 392.15
14 683.8 392.15
15 666.64 363.11
16 625.3 363.55
17 527.8 277.75
18 488.8 277.75
19 449.8 278.19
20 410.15 277.75
21 391.3 306.35
22 430.3 306.35
Coordinator
Nov 9, 2016 at 10:16 AM
Could you explain, what you think is wrong with the Voronoi diagram?
Nov 10, 2016 at 1:32 AM
Thanks for your reply. The Triangle is really a great work.
You can view this page link: https://github.com/xenovacivus/PathCAM/issues/15 (Sorry, I don't know how to insert a image here.)
However, I have found the wrong place of the code. It is wrong to judge the condition of the line going upwards. Please confirm whether it is correct.
private bool BoxRayIntersection(Point pt, double dx, double dy, out Vertex intersect)
        {
            double x = pt.X;
            double y = pt.Y;

            double t1, x1, y1, t2, x2, y2;

            // Bounding box
            double minX = bounds.Xmin;
            double maxX = bounds.Xmax;
            double minY = bounds.Ymin;
            double maxY = bounds.Ymax;

            // Check if point is inside the bounds
            if (x < minX || x > maxX || y < minY || y > maxY)
            {
                intersect = null;
                return false;
            }

            // Calculate the cut through the vertical boundaries
            if (dx < 0)
            {
                // Line going to the left: intersect with x = minX
                t1 = (minX - x) / dx;
                x1 = minX;
                y1 = y + t1 * dy;
            }
            else if (dx > 0)
            {
                // Line going to the right: intersect with x = maxX
                t1 = (maxX - x) / dx;
                x1 = maxX;
                y1 = y + t1 * dy;
            }
            else
            {
                // Line going straight up or down: no intersection possible
                t1 = double.MaxValue;
                x1 = y1 = 0;
            }

            // Calculate the cut through upper and lower boundaries
            if (dy < 0)
            {
                // Line going downwards: intersect with y = minY
                t2 = (minY - y) / dy;
                x2 = x + t2 * dx;
                y2 = minY;
            }
            //else if (dx > 0) //wrong place
            else if (dy > 0)  //correct code?
            {
                // Line going upwards: intersect with y = maxY
                t2 = (maxY - y) / dy;
                x2 = x + t2 * dx;
                y2 = maxY;
            }
            else
            {
                // Horizontal line: no intersection possible
                t2 = double.MaxValue;
                x2 = y2 = 0;
            }

            if (t1 < t2)
            {
                intersect = new Vertex(x1, y1, -1);
            }
            else
            {
                intersect = new Vertex(x2, y2, -1);
            }

            return true;
        }
Nov 10, 2016 at 2:15 AM
By the way, I'd like to know how to get polygons of the voronoi regions.
My code is as follows. I have encountered with another problem. I wonder if it's wrong to get polygons in my way. o(╯□╰)o
var polygons = new List<DoublePolygon>();
            foreach (var region in voronoi.Regions)
            {
                var ps = new DoublePolygon();
                foreach (var item in region.Vertices)
                {
                    ps.Add(new DoublePoint(item.X, item.Y));
                }
                polygons.Add(new DoublePolygon(ps));
            }
Nov 10, 2016 at 5:50 AM
Another problem is coming.
When the points have too many decimal places, the result is wrong sometimes.
I have a set of points input:
49
1 281.17 324.959999999998
2 133.67 315.809999999998
3 91.95 367.130000000005
4 109.05 320.15
5 240.7 274.940000000002
6 158.29 311.459999999998
7 182.91 307.119999999995
8 181.88 373.259999999998
9 157.26 377.6
10 132.64 381.940000000002
11 108.02 386.280000000005
12 116.57 362.790000000002
13 141.19 358.45
14 165.81 354.109999999998
15 190.42 349.769999999995
16 174.36 330.609999999998
17 149.74 334.95
18 125.12 339.290000000002
19 215.05 345.419999999995
20 198.98 326.269999999995
21 223.6 321.930000000005
22 248.22 317.590000000002
23 256.77 294.1
24 281.39 289.759999999998
25 232.15 298.440000000002
26 207.53 302.780000000005
27 323.62 370.380000000005
28 255.74 360.230000000005
29 231.11 364.569999999995
30 239.66 341.080000000005
31 206.5 368.909999999998
32 222.56 388.069999999995
33 247.18 383.730000000005
34 311.84 326.940000000002
35 312.5 394.190000000002
36 296.43 375.040000000002
37 271.81 379.390000000002
38 304.98 351.55
39 280.36 355.890000000002
40 264.29 336.740000000002
41 272.84 313.25
42 362.77 319.380000000005
43 338.15 323.719999999995
44 371.32 295.880000000005
45 322.08 304.559999999998
46 337.12 389.85
47 345.67 366.359999999998
48 329.6 347.209999999998
49 354.22 342.869999999995

And the result is wrong. Please view the picture:https://cloud.githubusercontent.com/assets/23348727/20165868/ce67164c-a74b-11e6-9ead-6c357fc690bd.jpg
or this page: https://github.com/xenovacivus/PathCAM/issues/15
(The Triangle TestApp you provide will run error. )
Coordinator
Nov 10, 2016 at 1:05 PM
Edited Nov 10, 2016 at 1:28 PM
Ok, you're using beta 3. You should have a look at the latest source code, there have been many improvements since back then. The Voronoi code has been rewritten completely. It now uses a DCEL to represent the Voronoi graph (which gives access to the complete topology). It's in the TriangleNet.Voronoi namespace. The old code can still be found in the TriangleNet.Voronoi.Legacy namespace.

The wrong condition in the intersection code was fixed in changeset 74966

Now, your second example is a completely different case. If you compute the statistics of the mesh, you will find that there's at least one triangle with tiny area (and minimum angle), which means you got a (nearly) flat triangle. Computing the circumcenter (= Voronoi vertex) of that triangle will result in a point with very large (probably useless) coordinates.

Looking at the pictures you have provided on github, you can see that the problem occurs on the right boundary of your mesh. The 5 points seem to be collinear, but due to finite floating point precision, the points might not be exactly collinear (at least one of the points is not on the convex hull but falls slightly inside the mesh). I haven't tested, but I think that's why a flat triangle is created.

This is a general problem when computing Delaunay triangulations. Depending on the problem you are trying to solve, a quick fix might be to perturb the points on the mesh boundary slighty and recompute the triangulation.
Nov 11, 2016 at 1:55 AM
You're awesome. And thanks for your kind help. I will try to achieve it in your solution. You guys are amazing. I'm still in the world of learning to use the Voronoi. Thank you again.
Coordinator
Nov 11, 2016 at 10:11 AM
Depending on what you're doing, there might be better solutions:

Instead of using the built-in Voronoi code, copy the stuff you need and then modify the code to ignore "invalid" triangles.
or
Remove "invalid" triangles after triangulation, reload the modified mesh and then use the built-in Voronoi code.
or
If you know that the point set is actually a polygon, triangulate the polygon. For example, if you export the mesh to Triangle format (a .poly and an .ele file will be created) and look at the poly file, you'll see that segments with vertices [44-42] and [42-46] make up the boundary to the right, so vertices 49 and 47 are the ones that mess up the triangulation. If those were added as segments in the input, the tiny triangles would be removed automatically by Triangle. The input poly file would look like this:
49
1 281.17 324.959999999998
2 133.67 315.809999999998
3 91.95 367.130000000005
4 109.05 320.15
5 240.7 274.940000000002
6 158.29 311.459999999998
7 182.91 307.119999999995
8 181.88 373.259999999998
9 157.26 377.6
10 132.64 381.940000000002
11 108.02 386.280000000005
12 116.57 362.790000000002
13 141.19 358.45
14 165.81 354.109999999998
15 190.42 349.769999999995
16 174.36 330.609999999998
17 149.74 334.95
18 125.12 339.290000000002
19 215.05 345.419999999995
20 198.98 326.269999999995
21 223.6 321.930000000005
22 248.22 317.590000000002
23 256.77 294.1
24 281.39 289.759999999998
25 232.15 298.440000000002
26 207.53 302.780000000005
27 323.62 370.380000000005
28 255.74 360.230000000005
29 231.11 364.569999999995
30 239.66 341.080000000005
31 206.5 368.909999999998
32 222.56 388.069999999995
33 247.18 383.730000000005
34 311.84 326.940000000002
35 312.5 394.190000000002
36 296.43 375.040000000002
37 271.81 379.390000000002
38 304.98 351.55
39 280.36 355.890000000002
40 264.29 336.740000000002
41 272.84 313.25
42 362.77 319.380000000005
43 338.15 323.719999999995
44 371.32 295.880000000005
45 322.08 304.559999999998
46 337.12 389.85
47 345.67 366.359999999998
48 329.6 347.209999999998
49 354.22 342.869999999995
#
# Now the boundary segments.
#
10
1 4 3
2 5 4
3 44 5
4 42 44
5 49 42 # These are the segments that are
6 47 49 # missing when triangulating just
7 46 47 # the point set.
8 35 46
9 11 35
10 3 11
0
Nov 14, 2016 at 1:29 AM
Thank you very much for your kind help and your better solutions. I'm using Triangle to implement a common tool, so I don't know what the input would be like in advance. Whereas, the first and second solution seem to be better. I will try it. Your help was greatly appreciated.
Nov 14, 2016 at 3:01 AM
It worked by using your second solution. But I'm not sure the way I filter the "invalid" triangles is right or wrong. After triangulation, I calculate the circumcenter of the triangle, and record "valid" triangles.
var triangles = new List<ITriangle>();
            foreach (var triangle in mesh.Triangles)
            {
               double bcX = 0, bcY = 0;
                GetTriangleBarycnt(triangle.GetVertex(0).X, triangle.GetVertex(0).Y, triangle.GetVertex(1).X,
                    triangle.GetVertex(1).Y, triangle.GetVertex(2).X, triangle.GetVertex(2).Y, ref bcX, ref bcY);
                if (Math.Abs(bcX) < 10e10 &&Math.Abs(bcY) < 10e10)  //filter condition???
                    triangles.Add(triangle);
            }
            mesh.Load(input, triangles);

            if (mesh.IsPolygon)
            {
                voronoi = new BoundedVoronoi(mesh);
            }
            else
            {
                voronoi = new Voronoi(mesh);
            }
For all that, I still need more tests to verify my program.
I am sorry to have troubled you.
Coordinator
Nov 14, 2016 at 10:08 AM
No trouble at all, you're welcome.

The current code will surely fail if the input vertices are of large magnitude. You could use the mesh.Bounds property to get the approx. center of the mesh and then compare to that position. Multiply with the width/height to get a relative scaling.

If you have a polygon input, you shouldn't remove any triangles. It might give unexpected results.

And I highly recommend to use the latest source code.
Nov 15, 2016 at 8:40 AM
Yes, you got the point. I'm skeptical about the way I removed "invalid" triangles.
And I wonder if I use the latest source code, how can I get the voronoi regions? I need to get polygons from regions.

In the previous version, my code is Like below:
var polygons = new List<DoublePolygon>();
if (voronoi.Regions.Count > 0)
            {
                foreach (var region in voronoi.Regions)
                {
                    var ps = new DoublePolygon();
                    foreach (var item in region.Vertices)
                    {
                        ps.Add(new DoublePoint(item.X, item.Y));
                    }
                    polygons.Add(new DoublePolygon(ps));
                }
            } 
Coordinator
Nov 15, 2016 at 10:26 AM
Moving to the beta 4 code is more about improvements to the triangulation code. The old Voronoi code is still available (see my comment above). If you want to use the new Voronoi diagram using DCELs, the following code should work (not tested):
var voronoi = new StandardVoronoi(mesh);

var polygons = new List<DoublePolygon>();

// DCEL vertex from TriangleNet.Topology.DCEL namespace.
Vertex last;

foreach (var face in voronoi.Faces)
{
    var ps = new DoublePolygon();
    
    foreach (var edge in face.EnumerateEdges())
    {
        ps.Add(new DoublePoint(edge.Origin.X, edge.Origin.Y));

        last = edge.Twin.Origin;
    }

    if (!face.Bounded)
    {
        // If the face is unbounded, we need to manually add the
        // last vertex.
        ps.Add(new DoublePoint(last.X, last.Y));
    }

    polygons.Add(new DoublePolygon(ps));
}
Nov 16, 2016 at 6:34 AM
Thank you for your great patience. You have given me great help. Thanks again. And have a nice day at work.
Nov 16, 2016 at 8:28 AM
I have another idea, just as a discussion. Calculate the area of the triangle, and then ignore triangles with small area values(such as the critical value is 0.001). I think it's much better than calculate the circumcenter of the triangle.
I learned a lot in the process of communication with you. Really thanks!
Coordinator
Nov 16, 2016 at 9:38 AM
Edited Nov 16, 2016 at 9:40 AM
Be careful with that approach. You might remove triangles that are perfectly "valid" (for example two points that are just very close together will make an edge of a thin triangle).

Triangle circumcenters

Imagine to make the triangles even thinner: triangle (1) will stay valid, triangle (2) won't. Both areas will get very small. The second case will usually occur on the boundary, but the first might also appear inside the mesh.

So, if you go with checking the area, also check that it's a boundary triangle (you can do this by checking that at least one of the neighbours is null).
Nov 17, 2016 at 2:16 AM
I overlooked this, and you are always able to hold on to the point.
However, I also need to check the triangles that adjacent to the "invalid" boundary triangles.
I have a situation, that is, an invalid triangle's neighbour is another invalid triangle. It is easy to ignore some of the invalid triangles which are close to each other.