Revise triangulation

- Avoid double-splitting a single point
- Check for mid-line intersections
- Avoid splitting zero-area polygons
- Keep full z-list for intersection checks

Fixes https://gitlab.com/kicad/code/kicad/-/issues/16241

libs/kimath/include/geometry/polygon_triangulation.h

@@ -66,6 +66,8 @@ public:
     bool TesselatePolygon( const SHAPE_LINE_CHAIN& aPoly )
     {
         m_bbox = aPoly.BBox();
+        m_prefactor_x = 32767.0 / m_bbox.GetWidth();
+        m_prefactor_y = 32767.0 / m_bbox.GetHeight();
         m_result.Clear();
 
         if( !m_bbox.GetWidth() || !m_bbox.GetHeight() )
@@ -271,8 +273,8 @@ private:
      */
     int32_t zOrder( const double aX, const double aY ) const
     {
-        int32_t x = static_cast<int32_t>( 32767.0 * ( aX - m_bbox.GetX() ) / m_bbox.GetWidth() );
+        int32_t x = static_cast<int32_t>( m_prefactor_x * ( aX - m_bbox.GetX() ) );
-        int32_t y = static_cast<int32_t>( 32767.0 * ( aY - m_bbox.GetY() ) / m_bbox.GetHeight() );
+        int32_t y = static_cast<int32_t>( m_prefactor_y * ( aY - m_bbox.GetY() ) );
 
         x = ( x | ( x << 8 ) ) & 0x00FF00FF;
         x = ( x | ( x << 4 ) ) & 0x0F0F0F0F;
@@ -301,7 +303,7 @@ private:
 
         while( p != aStart )
         {
-            if( area( p->prev, p, p->next ) == 0.0 )
+            if( *p == *( p->next ) || area( p->prev, p, p->next ) == 0.0 )
             {
                 p = p->prev;
                 p->next->remove();
@@ -310,6 +312,7 @@ private:
                 if( p == p->next )
                     break;
             }
+
             p = p->next;
         };
 
@@ -436,26 +439,38 @@ private:
                 continue;
             }
 
-            Vertex* nextNext = next->next;
+            Vertex* p = next;
+            bool removed = false;
 
-            if( *prev != *nextNext && intersects( prev, aPoint, next, nextNext ) &&
+            do
-                    locallyInside( prev, nextNext ) &&
-                    locallyInside( nextNext, prev ) )
             {
-                m_result.AddTriangle( prev->i, aPoint->i, nextNext->i );
+                Vertex* nextNext = p->next->next;
+                prev = p->prev;
 
-                // remove two nodes involved
+                if( *prev != *nextNext && intersects( prev, p, p->next, nextNext ) &&
-                next->remove();
+                        locallyInside( prev, nextNext ) &&
-                aPoint->remove();
+                        locallyInside( nextNext, prev ) )
+                {
+                    m_result.AddTriangle( prev->i, p->i, nextNext->i );
 
-                aPoint = nextNext;
+                    // remove two nodes involved
-                stop = nextNext;
+                    p->next->remove();
+                    p->remove();
 
-                continue;
+                    next = nextNext;
-            }
+                    p = nextNext;
+                    removed = true;
+                }
+
+                p = p->next;
+
+            } while ( p != next );
 
             aPoint = next;
 
+            if( removed )
+                continue;
+
             /*
              * We've searched the entire polygon for available ears and there are still
              * un-sliced nodes remaining.
@@ -472,7 +487,9 @@ private:
                 }
 
                 // If we don't have any NULL triangles left, cut the polygon in two and try again
-                splitPolygon( aPoint );
+                if ( !splitPolygon( aPoint ) )
+                    return false;
+
                 break;
             }
         }
@@ -557,7 +574,7 @@ private:
      * independently.  This is assured to generate at least one new ear if the
      * split is successful
      */
-    void splitPolygon( Vertex* start )
+    bool splitPolygon( Vertex* start )
     {
         Vertex* origPoly = start;
 
@@ -565,19 +582,29 @@ private:
         {
             Vertex* marker = origPoly->next->next;
 
+            if( m_splits.count( origPoly ) )
+            {
+                origPoly = origPoly->next;
+                continue;
+            }
+
             while( marker != origPoly->prev )
             {
+                if( m_splits.count( marker ) )
+                {
+                    marker = marker->next;
+                    continue;
+                }
+
                 // Find a diagonal line that is wholly enclosed by the polygon interior
                 if( origPoly->i != marker->i && goodSplit( origPoly, marker ) )
                 {
                     Vertex* newPoly = origPoly->split( marker );
 
-                    origPoly->updateList();
+                    m_splits.insert( origPoly );
-                    newPoly->updateList();
+                    m_splits.insert( marker );
 
-                    earcutList( origPoly );
+                    return ( earcutList( origPoly ) && earcutList( newPoly ) );
-                    earcutList( newPoly );
-                    return;
                 }
 
                 marker = marker->next;
@@ -585,6 +612,8 @@ private:
 
             origPoly = origPoly->next;
         } while( origPoly != start );
+
+        return false;
     }
 
     /**
@@ -597,10 +626,15 @@ private:
      */
     bool goodSplit( const Vertex* a, const Vertex* b ) const
     {
-        return a->next->i != b->i &&
+        bool a_on_edge = ( a->nextZ && *a == *a->nextZ ) || ( a->prevZ && *a == *a->prevZ );
-               a->prev->i != b->i &&
+        bool b_on_edge = ( b->nextZ && *b == *b->nextZ ) || ( b->prevZ && *b == *b->prevZ );
-               !intersectsPolygon( a, b ) &&
+        bool no_intersect = a->next->i != b->i && a->prev->i != b->i && !intersectsPolygon( a, b );
-               locallyInside( a, b );
+        bool local_split = locallyInside( a, b ) && locallyInside( b, a ) && middleInside( a, b );
+        bool same_dir = area( a->prev, a, b->prev ) != 0.0 || area( a, b->prev, b ) != 0.0;
+        bool has_len = ( *a == *b ) && area( a->prev, a, a->next ) > 0 && area( b->prev, b, b->next ) > 0;
+
+
+        return no_intersect && local_split && ( same_dir || has_len ) && !a_on_edge && !b_on_edge;
     }
 
     /**
@@ -611,6 +645,22 @@ private:
         return ( q->y - p->y ) * ( r->x - q->x ) - ( q->x - p->x ) * ( r->y - q->y );
     }
 
+    constexpr int sign( double aVal ) const
+    {
+        return ( aVal > 0 ) - ( aVal < 0 );
+    }
+
+    /**
+     * If p, q, and r are collinear and r lies between p and q, then return true.
+    */
+    constexpr bool overlapping( const Vertex* p, const Vertex* q, const Vertex* r ) const
+    {
+        return q->x <= std::max( p->x, r->x ) &&
+               q->x >= std::min( p->x, r->x ) &&
+               q->y <= std::max( p->y, r->y ) &&
+               q->y >= std::min( p->y, r->y );
+    }
+
     /**
      * Check for intersection between two segments, end points included.
      *
@@ -621,8 +671,28 @@ private:
         if( ( *p1 == *q1 && *p2 == *q2 ) || ( *p1 == *q2 && *p2 == *q1 ) )
             return true;
 
-        return ( area( p1, q1, p2 ) > 0 ) != ( area( p1, q1, q2 ) > 0 )
+        int sign1 = sign( area( p1, q1, p2 ) );
-                && ( area( p2, q2, p1 ) > 0 ) != ( area( p2, q2, q1 ) > 0 );
+        int sign2 = sign( area( p1, q1, q2 ) );
+        int sign3 = sign( area( p2, q2, p1 ) );
+        int sign4 = sign( area( p2, q2, q1 ) );
+
+        if( sign1 != sign2 && sign3 != sign4 )
+            return true;
+
+        if( sign1 == 0 && overlapping( p1, p2, q1 ) )
+            return true;
+
+        if( sign2 == 0 && overlapping( p1, q2, q1 ) )
+            return true;
+
+        if( sign3 == 0 && overlapping( p2, p1, q2 ) )
+            return true;
+
+        if( sign4 == 0 && overlapping( p2, q1, q2 ) )
+            return true;
+
+
+        return false;
     }
 
     /**
@@ -666,6 +736,28 @@ private:
             return area( a, b, a->prev ) < 0 || area( a, a->next, b ) < 0;
     }
 
+    /**
+     * Check to see if the segment halfway point between a and b is inside the polygon
+    */
+    bool middleInside( const Vertex* a, const Vertex* b ) const
+    {
+        const Vertex* p = a;
+        bool          inside = false;
+        double        px = ( a->x + b->x ) / 2;
+        double        py = ( a->y + b->y ) / 2;
+
+        do
+        {
+            if( ( ( p->y > py ) != ( p->next->y > py ) )
+                && ( px < ( p->next->x - p->x ) * ( py - p->y ) / ( p->next->y - p->y ) + p->x ) )
+                inside = !inside;
+
+            p = p->next;
+        } while( p != a );
+
+        return inside;
+    }
+
     /**
      * Create an entry in the vertices lookup and optionally inserts the newly created vertex
      * into an existing linked list.
@@ -696,7 +788,10 @@ private:
 
 private:
     BOX2I                                 m_bbox;
+    double                                m_prefactor_x;
+    double                                m_prefactor_y;
     std::deque<Vertex>                    m_vertices;
+    std::set<Vertex*>                     m_splits;
     SHAPE_POLY_SET::TRIANGULATED_POLYGON& m_result;
 };