search

gwlucastrig / Tinfour

Delaunay and Constrained Delaunay Triangulations in Java, providing high-performance utilities for modeling surfaces with support for Lidar LAS files, Digital Elevation Models (DEM), finite element analysis, path planning, natural neighbor interpolation, and other applications of Triangulated Irregular Networks (TIN)
Apache License 2.0
153 stars 34 forks source link

Natural Neighbor Interpolator fails at borders of problematic TIN #53

Closed gwlucastrig closed 3 years ago

gwlucastrig commented 3 years ago

The Natural Neighbor interpolation can fail near the borders of a mesh with problematic positioning of vertices near the borders of the TIN. For example, the following test set was contributed by a user who was having problems with interpolation. It gives rise to a mesh which has very thin triangles located along three edges of the mesh. Those triangles lead to the computation of circumcircles with very large radius values. In these cases, numerical issues lead to incorrect interpolations.

index, x,y,z
0,  120.83333333333333, 398.6666666666667,  28.0
1,  204.16666666666666, 351.3333333333333,  23.0
2,  292.49999999999994, 398.66666666666663, 18.0
3,  375.8333333333333,  351.3333333333333,   4.0
4,  464.16666666666663, 398.66666666666663,  0.0
5,  204.16666666666666, 256.6666666666667,  12.0
6,  287.5,              209.33333333333331,  5.0
7,  375.8333333333333,  256.66666666666663,  9.0
8,  287.5,              114.66666666666666,  0.0

The following images show the layout of the vertices and the resulting interpolation. Vertices 2 and 5 are just inside the line segments that would be formed by Vertices 0,4 and Vertices 0,8. The resulting triangles 2,4,0 and 5,8,0 have very small areas and are almost degenerate. In the picture below, the triangles are so flattened that their three points appear to be lying on a single line. Numerical issues result.

NNI_Problem_Wireframe

Here are the color-coded results from interpolation. Numerical issues lead to erratic behavior. NNI_Problem

gwlucastrig commented 3 years ago

A code change for the Natural Neighbor Interpolator class has been posted to the Tinfour code base. This change corrects the behavior of the interpolator. This issue may now be closed.

NNI_Fix

gwlucastrig commented 3 years ago

In a related development, I noticed that the getArea() method for the SimpleTriangle class (see Triangle Collector Techniques ) was computing a small-magnitude negative area for triangle 5,8,0. This glitch violates the rule that all triangles produced by Tinfour are oriented counterclockwise and have a positive area.

This unfortunate result was due to numeric round-off issues arising from the nearly-degenerate geometry of the triangle. However, when I used the extended-precision logic from the GeometricOperations class, I confirmed that the triangle was, in fact, constructed correctly and did have a positive area.

I am looking at the best way to incorporate extended-precision logic into the SimpleTriangle class. The extended-precision logic is slower than simple floating-point operations but does produce more accurate results.

This change will be integrated into the upcoming 2.1.5 release of Tinfour. Once it is addressed, I will be closing this issue.

gwlucastrig commented 3 years ago

I have updated the area calculation for the SimpleTriangle class (now not quite so simple as it was).

Previously, triangle 5,8,0 was computed to have an area of -1.818989e-12. While the formula in the code was correct, numerical issues resulted in the incorrect value. The current version uses extended-precision arithmetic to find an area of 1.752672e-13

As you can see, the new computation produces a more reasonable value and confirms that the triangle has a positive area. Of course, it is possible that if vertex 5 were moved arbitrarily close to line segment 0,8 it could result in a calculation that exceeded the limits of the extended-precision. Computational geometry is notorious for numerical issues. Currently, Tinfour does include logic to tweak the coordinates of a vertex if it gets too close to a line segment and to insert it into the line rather than building a triangle. In this case, vertex 5 wasn't quite close enough to trigger the adjustment. While this situation did not affect the integrity of the TIN, it did challenge the ordinary-precision area computation.

Jonathan Shewchuk has published algorithms and C-language code for various triangle-geometry computations that work under all conditions. Mr. Shewchuk is pretty much the leading authority on Delaunay Triangulations, an excellent writer, and a pretty fair coder. If you are interested in the Delaunay, his work is well worth considering.

gwlucastrig commented 3 years ago

The code changes for this issue are now complete.

I would like to thank Michael Kemper for identifying this problem and providing such an excellent set of test data. Those 9 vertices in his sample set exercise the Tinfour logic with admirable thoroughness.

gwlucastrig commented 3 years ago

This issue in now fully addressed. Fixes are included in the 2.1.5 release which was issued 10 Jan 2021.