mapbox / delaunator

An incredibly fast JavaScript library for Delaunay triangulation of 2D points
https://mapbox.github.io/delaunator/
ISC License
2.24k stars 139 forks source link

Triangulation is not robust (fails for certain floating point input) #61

Closed MK-3PP closed 3 years ago

MK-3PP commented 4 years ago

Note: A more thorough introduction is here in the original issue at delaunator-cpp .

I am processing point clouds around the size of 100 k to 1.000 k Points. The points' outer hull is a near perfect rectangle and their distribution is grid like. In pre-processing the point cloud gets slightly tilted towards a plane for error compensation.

While I have not yet observed the issue below in our raw data, occasionally there was a number of big triangles (hundreds in a cloud of a million points) spanning across the whole pre-processed (tilted) point cloud. To reproduce this, I wrote a test. It generates simple square grid, rotates it by a few degrees and then calls delaunator to triangulate.

I have ported my test from C++ to JS by modifying your demo (Forgive me, it is my first JS ever). Erroneous grids get dumped to console. You can rotate the grid by moving the mouse left or right on the canvas. Uncomment the "Test loop" bit to automatically iterate over 10.000 test cases. One grid that fails validation is at 302.328° rotation (Firefox 74.0, 64-bit, Win10):

i accidentally the hull

points = [[-3.0381276552207055, 10.481881920449052], [-2.1931270567413446, 11.016647278872279], [-1.3481264582619854, 11.551412637295508], [-0.5031258597826245, 12.086177995718735], [0.3418747386967347, 12.620943354141964], [1.1868753371760938, 13.155708712565193], [-2.5033622967974765, 9.63688132196969], [-1.6583616983181173, 10.171646680392918], [-0.8133610998387582, 10.706412038816147], [0.03163949864060278, 11.241177397239376], [0.8766400971199619, 11.775942755662605], [1.721640695599322, 12.310708114085832], [-1.9685969383742474, 8.791880723490332], [-1.1235963398948883, 9.326646081913559], [-0.2785957414155291, 9.861411440336788], [0.5664048570638318, 10.396176798760017], [1.411405455543191, 10.930942157183246], [2.25640605402255, 11.465707515606473], [-1.4338315799510184, 7.9468801250109715], [-0.5888309814716592, 8.4816454834342], [0.2561696170076999, 9.016410841857429], [1.10117021548706, 9.551176200280658], [1.94617081396642, 10.085941558703885], [2.791171412445779, 10.620706917127112], [-0.8990662215277911, 7.1018795265316115], [-0.05406562304843021, 7.6366448849548405], [0.7909349754309281, 8.17141024337807], [1.635935573910288, 8.706175601801297], [2.4809361723896473, 9.240940960224526], [3.3259367708690073, 9.775706318647753], [-0.3643008631045621, 6.256878928052252], [0.48069973537479704, 6.7916442864754805], [1.3257003338541562, 7.326409644898709], [2.1707009323335162, 7.861175003321938], [3.0157015308128763, 8.395940361745165], [3.8607021292922354, 8.930705720168394]];

Although the error is similar, it does occur under different rotation angles under C++ and JS. Also, @abellgithub had at first not been able to recreate the error under OSX. I suppose the difference lies in different function implementations in the standard libraries, e.g. for sine and cosine that lead to different rounding errors.

mourner commented 4 years ago

This seems to be a robustness issue related to orientation checks. If I use the robust orient2d check from robust-predicates, the test case works properly — so it seems that #51 wasn't enough, and it's time to start using the fully robust check. Maybe we could get away with just orient2d for now (since incircle is pretty big so will increase the bundle quite a bit).

For the C++ version, I guess you can use the original Shewchuk's predicates (written in C).

Samfox2 commented 4 years ago

@MK-3PP I am facing the same problem. Could you share your c++ code to see if it is also solved with my point cloud?

MK-3PP commented 4 years ago

@Samfox2 I posted my test code in the opening post here. But validation depends on prior knowledge about the point cloud. In this case: regular 1-by-1 grid means that triangles are halved 1-by-1 squares. So it will most likely not apply to your point cloud.

If you meant a method for visual inspection see here. I just save the mesh as text file and display it in Octave.

Samfox2 commented 4 years ago

Great thanks. My point cloud is not regular but I am facing the same problem - about 20 from roughly 1 million triangles show the same behavior: occuring of big overlapping triangles

mcwhittemore commented 3 years ago

@mourner I think I'm running into this bug writing using delaunator to recursively break a polygon up into small triangles. Here is a gist showing what I'm doing. Would it be helpful if I PR'd this as a failing test case?

mourner commented 3 years ago

@mcwhittemore can you try the v5 branch I just pushed? I switched the orientation check to the one from robust-predicates there. The only gotcha is that it's now native ESM, so you'd have to also switch to Node v12+, "type": "module" (or .mjs extension) and import instead of require.

mcwhittemore commented 3 years ago

@mourner The result changed but the effect stayed the same using 9b457a9db823c1dd23dcd63bee4d80b2978d71d7. Here is the update code with the output I'm seeing.. I'll try to PR a test vs the v5 branch to better show what is happening.

mourner commented 3 years ago

Fixed by #68.