A case that failed in triangulation
pan3rock opened this issue · 20 comments
I encountered some failed cases that said "Duplicate vertex detected: # is a duplicate of #". After checking it, I found it's because of the triangulation at the last iteration in my codes. There is a simple example below. It's kind of like #166. I don't know the reason, is it a bug or due to invalid points?
The vertices are
0 0.000010000000 0.000010000000
1 1.000000000000 0.000000000000
2 0.999990000000 0.099990000000
3 0.000000000000 0.100000000000
4 0.500005000000 0.000005000000
5 0.500000000000 0.050000000000
6 0.499995000000 0.099995000000The generated triangles are
5 0 4
0 5 3
5 2 6
2 3 6
3 5 6
2 5 1
4 1 5The triangle (2, 3, 6) is the red line, a flat triangle.
2 0.999990000000 0.099990000000
3 0.000000000000 0.100000000000
6 0.499995000000 0.099995000000
It's clear these points aren't collinear from their y values! you can round your input points to match with your need.
Yes, it is the duplicate. The points are not exactly on the same line. Exactly means that if you solve orient2d predicate: the determinant is exactly zero. It is easy to see that 0.09999 != 0.1 or 0.099995.
@artem-ogre @Islam0mar The line is slanted.. After checking the area or slope, the three points are collinear.
x1, y1 = 0.999990000000, 0.099990000000 # 2
x2, y2 = 0.000000000000, 0.100000000000 # 3
x3, y3 = 0.499995000000, 0.099995000000 # 6
def collinear(x1, y1, x2, y2, x3, y3):
a = x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2)
if (a == 0):
print("Yes")
else:
print("No")
collinear(x1, y1, x2, y2, x3, y3)
print((y1 - y2)/(x1 - x2) == (y2 - y3)/(x2 - x3))The output is Yes.
This code is likely returning “Yes” because of the floating point rounding errors. This is exactly the problem that is solved by the robust geometric predicates (RGPs). CDT is using RGPs. There are a lot of very good RGP explanations that can be found online.
Thanks. I tried the RGPs provided by predicates. The value of orient2d of the three points is also 0. I'll check it again tomorrow and find if there is another mistake.
What about this one used in CDT?
Line 115 in 58f34da
I used the orient2d in CDT, the results are below:
5 0 4 2.49973e-02
0 5 3 4.99945e-02
5 2 6 2.49973e-02
2 3 6 6.93847e-18
3 5 6 2.49973e-02
2 5 1 4.99945e-02
4 1 5 2.49973e-02
The value for (2, 3, 6) is 6.93847e-18. It may be a precision problem.
I print the orient2d in my codes when it's absolute value < 1.0e-5, it shows:
29 114 132 0.00000e+00
112 57 98 0.00000e+00
51 92 0 0.00000e+00
211 183 223 0.00000e+00
145 127 52 -1.73472e-18
195 236 222 0.00000e+00
212 209 193 0.00000e+00
119 155 100 -2.60209e-18
95 79 77 0.00000e+00
iter_b: 13
29 114 132 0.00000e+00
183 223 211 0.00000e+00
51 92 0 0.00000e+00
212 209 193 0.00000e+00
112 57 98 0.00000e+00
245 244 263 3.46945e-18
127 52 145 -1.73472e-18
195 236 222 0.00000e+00
119 155 100 -2.60209e-18
95 79 77 0.00000e+00
iter_b: 14
127 52 145 -1.73472e-18
112 57 98 0.00000e+00
51 92 0 0.00000e+00
183 223 211 0.00000e+00
212 209 193 0.00000e+00
245 244 263 3.46945e-18
195 236 222 0.00000e+00
132 29 114 -1.73472e-18
119 155 100 -2.60209e-18
77 95 79 0.00000e+00
you will find value of some triangles are at the level of 1.0e-18 and some are zeros. Under my understanding, these triangles should be ignored when generating, but now they exist.
I also computed orient2d(p1, p2, p3), orient2d(p1, p3, p2), orient2d(p3, p2, p1) for those weird triangles
e.g.
15 8 3 0.00000e+00
15 8 3 -1.73472e-18
15 8 3 -8.67362e-19
e.g.
93 175 176 2.16840e-19
93 175 176 -2.16840e-19
93 175 176 -2.16840e-19
It says that those with values of zero will be at levels of <1.0e-18 if sequences are changed.
Are all of these results for adaptive::orient2d from predicates used in CDT?
Yes.
auto v1 = predicates::adaptive::orient2d(x1, y1, x2, y2, x3, y3);In the last examples change of sign is due to changing the order of points (which changes winding). It however should not change from zero to negative in my understanding. What are the coordinates of points (15, 8, 3)?
Was your other claim that orient2d should return zero for (2,3,6) compared to another RGP implementation.
Another claim is that there are triangles in final triangulation for which adaptive::orient2d is zero?
I would appreciate the minimal reproducers. In that case I could take a closer look when I’m back at the keyboard.
A demo is provided. You can observe the change from zero to other values. Another RGP implementation is also included. In the first case, they are the same, but the second is different.
the double precision is 1e-16, is that 1.0e-18 should be regarded as zero? Like my job, the update of new vertices is from some splitted edges (x_new = (x_i + x_j) / 2.0), the location has an error of 1.0e-16, is it possible that's the reason to produce such problem?
git clone https://github.com/pan3rock/CDT-issue-172.git
cd CDT-issue-172
git submodule init
git submodule update
mkdir build
cd build
cmake ..
make -j4
../bin/demo1
my output:
(6, 4, 7): 2.16840e-19
(6, 7, 4): -2.16840e-19
(7, 4, 6): -2.16840e-19
another predicates: 2.16840e-19
another predicates: -2.16840e-19
another predicates: -2.16840e-19
(1, 2, 3): -8.67362e-19
(1, 3, 2): 4.33681e-19
(3, 2, 1): 0.00000e+00
another predicates: -4.33681e-19
another predicates: 4.33681e-19
another predicates: 4.33681e-19
Only exactly zero should be considered as zero. E.g., smallest representable positive double 2.225E-307 should be considered positive.
What happens if you change fmt::print to display the actual value of double?
you mean to print the double value with full precision?
fmt::print("({:d}, {:d}, {:d}): {}\n", v[0], v[1], v[2], v1);output:
(6, 4, 7): 2.168404344971009e-19
(6, 7, 4): -2.168404344971009e-19
(7, 4, 6): -2.168404344971009e-19
another predicates: 2.168404344971009e-19
another predicates: -2.168404344971009e-19
another predicates: -2.168404344971009e-19
(1, 2, 3): -8.673617379884035e-19
(1, 3, 2): 4.336808689942018e-19
(3, 2, 1): 0
another predicates: -4.336808689942018e-19
another predicates: 4.336808689942018e-19
another predicates: 4.336808689942018e-19
(1, 2, 3): -8.673617379884035e-19 (1, 3, 2): 4.336808689942018e-19 (3, 2, 1): 0 another predicates: -4.336808689942018e-19 another predicates: 4.336808689942018e-19 another predicates: 4.336808689942018e-19
Doesn’t look right. I will check (1,2,3)-case next week.
Can't reproduce on gcc and clang:
https://godbolt.org/z/s5MqjrTvf
It turns out that predicate breaks when -Ofast flag is used, works fine with -O3
https://godbolt.org/z/Kx9xrcETs
from Ofast docs:
Disregard strict standards compliance. -Ofast enables all -O3 optimizations. It also enables optimizations that are not valid for all standard-compliant programs. It turns on -ffast-math, -fallow-store-data-races and the Fortran-specific -fstack-arrays, unless -fmax-stack-var-size is specified, and -fno-protect-parens. It turns off -fsemantic-interposition.
ffast-math is not supported by robust predicates.
For the (2, 3, 6) case both predicates agree that points are collinear with O3.
My guess is that if Ofast flag is removed, the triangulation will not contain flat triangles.
Great, you're right. It's Ofast. Thanks a lot.