Efficient triangulation of weakly simple polygon with known, ordered boundary points

Question

I have a point set P which represents the (ordered) boundary points of a planar polygon, without any points on the interior. The polygon may be concave, and is at minimum weakly simple. It may be possible to constrain all polygons to be simple. My goal is to triangulate the planar polygon faster than O(n^2).

The size of set P ranges from several hundreds of points to 1-2 million. For reasons specific to how these polygons are generated, it is difficult to simply down-sample and get a smaller set P.

I'm working in Python -- I don't really have any restrictions about which packages I use, or what licenses they have.

Approaches I've tried

My "stable" software is built around an earclipping algorithm. This worked very well for small examples, with hundreds or thousands of points, but has excessive runtime for large polygons. It does not leverage knowledge of the boundary and boundary ordering (and to be nitpicky, produces very ugly meshes).

More recently I have tried using Triangle, through some python bindings, with the constrained-Delaunay algorithm. Even when supplying the boundary edges, constrained Delaunay run on concave shapes returns triangles inside the polygon boundary cavities. Note: I do not require a (quasi-) Delaunay triangulation, I am looking at Triangle because it's compact, well-regarded, and fast. A simple solution to the concavity-filling is to check whether the generated tris are inside the boundary (for example, computing the centroid of each tri and running a point-in-polygon algorithm). However, simple point-in-polygon algorithms such as ray-casting run in O(n^2) time because for n points in P there are n boundary edges which must be checked for ~n generated tris.

Specific forms of the question

Is there a (different? than earclipping or constrained delaunay) algorithm which takes advantage of my highly-specific case, or
is there a more efficient way of checking whether generated tris are inside my boundary?

Based on some other SO posts I've linked to below, it sounds like one or both of the above are satisfied by Triangle and I simply am not taking full advantage of it (possibly also a limitation of the python bindings?). I have functionally-nonexistent knowledge of C; digging through the Triangle code has been fruitless thus far.

1. For the fastest triangulations

There is 2 kind of algorithms (as I know):

those who directly triangulates
those who split the given non-convex outline in smaller convex outlines (or at least monotone), and only then triangulate them (this is easy because this is convex)

Both are based on what's called sweepline algorithms which is in short:

choose a direction in your plane (whatever it is)
sort your points against that direction
consider matching points in the order they come from that sorted list

[ ] - A publication for the 1st kind

[ ] - this wikipedia page can also be usefull.

2. for the nicest triangulations

There is multiple ways but the best I know is triangulating by reducing the outline:

choose one point you will place a triangle with its direct neighboors
create an triangle and remove the latest point from the outline
do it again until there is no more outline
have some heuristic to know each iteration which triangle to make in priority

This is what's done in madcad.triangulation, Its complexity is O(n*k) with k = number of concave points

from madcad.triangulation import *

outline = Wire([ ... your points ... ])

triangulation(outline)  # this is triangulating the outline, it even works with holes in the surface
triangulation_outline(outline)  # this function is a subpart of the latest, working only on outlines with no holes

conclusion

I depends on what you need :) I personally abandonned the fastest triangulations for the most of my use cases, because of the extreemly poor quality of the triangles in case of smooth outlines.

If you need to perform some computation on the resulting mesh, the very thin triangles will lead to significan loss of precision using floating points.