triangulation failed for surfaces built with non-convex polygons

[videlec]

See eg #238.

from flatsurf import Polygon, MutableOrientedSimilaritySurface
K = QuadraticField(3)
a = K.gen()
# build the vectors
l = [(1, 0), (1, 2), (a, 11), (a, 1), (1, 4), (1, 6), (a, 3),
     (a, 5), (1, 8), (1, 6), (a, 9), (a, 7), (1, 10), (1, 0)]
vecs = []
for m, e in l:
    v = vector(K, [m * cos(2*pi*e/12), m * sin(2*pi*e/12)])
    vecs.append(v)
p = Polygon(edges=vecs)

from collections import defaultdict
d = defaultdict(list)
for i, e in enumerate(p.edges()):
    e.set_immutable()
    d[e].append(i)

Sbase = MutableOrientedSimilaritySurface(K)
_ = Sbase.add_polygon(p)
for v in list(d):
    if v in d:
        indices = d[v]
        v_op = -v
        v_op.set_immutable()
        opposite_indices = d[v_op]
        assert len(indices) == len(opposite_indices), (len(indices), len(opposite_indices))
        if len(indices) == 1:
            del d[v]
            del d[v_op]
            Sbase.glue((0, indices[0]), (0, opposite_indices[0]))

(i0, j0), (i1, j1) = d.values()

S1 = MutableOrientedSimilaritySurface.from_surface(Sbase)
S1.glue((0, i0), (0, i1))
S1.glue((0, j0), (0, j1))
S1.set_immutable()

S1.triangulate()
ValueError: polygon has negative area; probably the vertices are not in counter-clockwise order

Even though triangulation of non-convex polygons seem to work:

from flatsurf import Polygon

K = QuadraticField(3)
a = K.gen()

# build the vectors
l = [(1, 0), (1, 2), (a, 11), (a, 1), (1, 4), (1, 6), (a, 3),
     (a, 5), (1, 8), (1, 6), (a, 9), (a, 7), (1, 10), (1, 0)]
vecs = []
for m, e in l:
    v = vector(K, [m * cos(2*pi*e/12), m * sin(2*pi*e/12)])
    vecs.append(v)
p = Polygon(edges=vecs)
p.triangulation()

[saraedum]

The polygon in question is Polygon(vertices=[(0, 0), (-3, 0), (-3/2, 1/2*a)]). It's created from "edges" [(-3, 0), (3/2, 1/2*a), (3/2, -1/2*a)].

This happens while triangulating the polygon Polygon(vertices=[(0, 0), (3/2, 1/2*a), (3, 0), (9/2, 1/2*a), (4, a), (3, a), (3, 2*a), (3/2, 5/2*a), (1, 2*a), (0, 2*a), (0, a), (-3/2, 1/2*a), (-1, 0)]), aka

We try to create a polygon from the vertices 0, 1, 2. And that's indeed not convex.

[saraedum]

It's a bit bizarre that this is not using the triangulation of the polygon but that surfaces are implementing their own triangulation.

So, triangulating an infinite surface uses polygon triangulation but finite surfaces implement it on their own 🤷

[saraedum]

Anyway, this is the fix:

diff --git a/flatsurf/geometry/surface.py b/flatsurf/geometry/surface.py
index 096c6ee..2e26f70 100644
--- a/flatsurf/geometry/surface.py
+++ b/flatsurf/geometry/surface.py
@@ -2056,7 +2056,7 @@ class MutableOrientedSimilaritySurface(
             for i in range(n):
                 e1 = poly.edge(i)
                 e2 = poly.edge((i + 1) % n)
-                if ccw(e1, e2) != 0:
+                if ccw(e1, e2) > 0:
                     # This is in case the polygon is a triangle with subdivided edge.
                     e3 = poly.edge((i + 2) % n)
                     if ccw(e1 + e2, e3) != 0: