delaunay sweep circle initialization

Question

I recently came accross an interesting "sweep circle" algorithm for delaunay triangulation [1].

It's basic idea is seeding with a trivial triangulation (3 points), then gradually extending it outwards (points are sorted by their distance from a chosen origin - the sweep-circle center) while preserving the legality of the triangulation by flipping triangles if needed, always starting from the new one, and continueing to neighbouring triangles of flipped ones.

I've implemented it, and it works nicely, but I still struggle with making the seeding part flawless.

The thing is :

• In order to ensure that each new point is outside of the current envelope ("frontier" in the article), you have to make sure that the swept points are closer to the origin than the unswept points.
• In order to find where to connect the new point, you have to have the origin inside the envelope (so angular sort of the envelope and binary search will do).

Starting at some iteration this will natuarally happen, but the problem is making it happen for the seed points.

So what I aim at is finding an origin which fulfills :

• Be inside the seed triangulation.
• Have the seed points the closest points to it (have the sweep circle strictly containing the seed points empty of other points).

The article suggests taking some average point as the origin, and start with the 3 points closest to it, but this clearly fails if these 3 points are on the same side of the chosen origin, so they form a triangle that is not containing it.

A different but somewhat similar algorithm [2] is suggesting choosing an arbitrary point, then taking the point closest to it, and the point that creates the smallest circumcircle with the other two, and fix the origin as the center of their circumcircle, but that's clearly not working either, unless these 3 points happen to form an acute triangle (is this circle even ensured to be empty of other points?).

Currently I use a similar idea - arbitrary point + 2 closest points to it, as the first 3 points, but I fix the origin as some point inside the formed triangle (e.g their average), which basically works, but can fail in principal.

I'm looking for another simple way to do it, or a way to detect when the proposed method is not working and fix it in a clever way involving the point that "ruined it", before the iterations begin and without having to deal with more than 2 seed triangles (4 points), or having to make many iterations or guesses at init stage.

I know I can take the first point and shift it by some small $\epsilon$ towards inside the triangle, and it should work, as long as it's smaller than say the minimum distance bewteen 2 points, but I'm looking for a clever solution, and I am especially interested in a solution that can chose an integer origin for a set of integer points, for then the whole algorithm can be implemented using only integral arithmetics (faster, less error prone).

Bottom line :

I think it all comes down to this :

Given a set of n points in the plane, find a point O which lies strictly inside the triangle formed by the 3 points of the set closest to it (to O).

But I guess that put in the context, there might be solutions to implementing a working version of the algorithm that bypass this question.

References :

1. A faster circle-sweep Delaunay triangulation algorithm by Ahmad Biniaz and Gholamhossein Dastghaibyfard (2011)
2. S-hull: a fast radial sweep-hull routine for Delaunay triangulation by D. A. Sinclair (2016)

Answer 1

I just realized that my question is irrelevant to making the algorithm work. I mixed the radial and angular atributes of the points and mistakably thought them to need to be based at the same origin, but that's just confusing semantics. Actually, the radial origin decides on the order of insertion of points and if set correctly, ensures no new point will apear inside an existing triangle, and the angular origin only helps to find the projection of a new point on the envelope (frontier / convex-hull in the case of [2]), so the angular origin indeed needs to be inside the seed envelope, but the radial origin need not. Moving the angular-origin inside the seed envelope will affect the actual progression of the algorithm, but will not affect correctness.

This makes things a lot easier - taking arbitrary point of the set, or somewhere in the plane and find the 2/3 points closest to it will work in general position, but it may result in a colinear triplet. Following [2] and taking a point from the set, its nearest neighbour, and a 3rd point that creates the smallest circumcircle with them will always work as this is a delaunay triangle (nearest-neighbours are always delaunay edges, and the 2 triangles based on this edge have to create an empty circumcircle, so one of them has to be the smallest) which makes the circumcircle empty of other points, so using its origin as the radial origin is correct.

To do this for integers, we may take the second approach and round the result (checking for super-edge-cases), or take the first approach which is more straight forward for integer radial origin, and restart if we find a colinear triplet, or when the triangle is too obtuse to contain an integer point in its interior. The last condition is only for making the algorithm fully integral (e.g for machines wihout Floating Point Arithmetics), but using precalculated float pseudo-angles [3] (avoiding use of slow atan2) and comparing them later in the binary search (tree), might be faster than using fully integral orientation test by demand.

3. pseudo-angles