Infinite loop due to reaching outer half-edge

[xionluhnis]

I've been looking for a solution as simple as the one this library is providing to get constrained Delaunay where the constraint is basically a polygon outline.

However, is the current implementation making specific assumptions about the point distribution? e.g. that there are no points inside of the constrained shape?

If not, then there may some missing edge cases. Here is a test case that fails because of an infinite loop:

const Delaunator = require('delaunator');
const Constrainautor = require('@kninnug/constrainautor');

const positions = [[2, 1], [3, 1], [4, 1], [5, 1], [6, 1], [7, 1], [8, 1], [9, 1], [10, 1], [11, 1], [12, 1], [13, 1], [2, 2], [3, 2], [4, 2], [5, 2], [6, 2], [7, 2], [8, 2], [9, 2], [10, 2], [11, 2], [12, 2], [13, 2], [2, 3], [3, 3], [4, 3], [5, 3], [6, 3], [7, 3], [8, 3], [9, 3], [10, 3], [11, 3], [12, 3], [13, 3], [2, 4], [3, 4], [4, 4], [5, 4], [6, 4], [7, 4], [8, 4], [9, 4], [10, 4], [11, 4], [12, 4], [13, 4], [2, 5], [3, 5], [4, 5], [5, 5], [6, 5], [7, 5], [8, 5], [9, 5], [10, 5], [11, 5], [12, 5], [13, 5], [2, 6], [3, 6], [4, 6], [5, 6], [6, 6], [7, 6], [8, 6], [9, 6], [10, 6], [11, 6], [12, 6], [13, 6], [2, 7], [3, 7], [4, 7], [5, 7], [6, 7], [7, 7], [8, 7], [9, 7], [10, 7], [11, 7], [12, 7], [13, 7], [2, 8], [3, 8], [4, 8], [5, 8], [6, 8], [7, 8], [8, 8], [9, 8], [10, 8], [11, 8], [12, 8], [13, 8], [1, 9], [2, 9], [3, 9], [4, 9], [5, 9], [6, 9], [7, 9], [8, 9], [9, 9], [10, 9], [11, 9], [12, 9], [13, 9], [1, 10], [2, 10], [3, 10], [4, 10], [5, 10], [6, 10], [7, 10], [8, 10], [9, 10], [10, 10], [11, 10], [12, 10], [13, 10], [1, 11], [2, 11], [3, 11], [4, 11], [5, 11], [6, 11], [7, 11], [8, 11], [9, 11], [10, 11], [11, 11], [12, 11], [13, 11], [1, 12], [2, 12], [3, 12], [4, 12], [5, 12], [6, 12], [7, 12], [8, 12], [9, 12], [10, 12], [11, 12], [12, 12], [1, 13], [2, 13], [3, 13], [4, 13], [5, 13], [6, 13], [7, 13], [8, 13], [9, 13], [10, 13], [11, 13], [12, 13], [1, 14], [2, 14], [3, 14], [4, 14], [5, 14], [6, 14], [7, 14], [8, 14], [9, 14], [10, 14], [11, 14], [12, 14], [0, 14.404910396560176], [0.1097833356951295, 13.512342151010488], [0.219566671390259, 12.6197739054608], [0.3293500070853885, 11.727205659911116], [0.439133342780518, 10.834637414361424], [0.5489166784756471, 9.942069168811738], [0.6587000141707766, 9.04950092326205], [0.768483349865906, 8.156932677712362], [0.8782666855610356, 7.264364432162673], [0.988050021256165, 6.371796186612986], [1.0978333569512946, 5.479227941063298], [1.2076166926464236, 4.58665969551361], [1.3174000283415532, 3.6940914499639224], [1.4271833640366824, 2.8015232044142344], [1.536966699731812, 1.908954958864547], [1.6467500354269418, 1.0163867133148592], [1.7565333711220712, 0.12381846776517148], [2.625357004835959, 0.18948424425307334], [3.49774548578131, 0.17614930575222743], [4.377480950347606, 0], [5.03667324191489, 0.03551984647853645], [5.695865533482174, 0.0710396929570729], [6.2976567830833465, 0.22587193892879634], [6.899448032684518, 0.38070418490052144], [7.003061369329261, 0.37532483863500393], [7.5514736580060084, 0.3124430512698645], [8.099885946682761, 0.2495612639047251], [8.717896070125907, 0.3415446219466246], [9.33590619356905, 0.4335279799885241], [10.26668585466961, 0.4655189291422352], [11.197465515770173, 0.4975098782959462], [11.844773441110322, 0.4704867438747797], [12.492081366450472, 0.44346360945361146], [13.099038695146895, 0.5505512714034877], [13.705996023843316, 0.6576389333533622], [14, 0.7212017967478266], [13.9812640368062, 1.5897138741882504], [13.9625280736124, 2.4582259516286733], [13.943792110418599, 3.3267380290690935], [13.9250561472248, 4.195250106509518], [13.906320184030996, 5.0637621839499385], [13.887584220837196, 5.932274261390362], [13.868848257643396, 6.800786338830786], [13.850112294449596, 7.669298416271207], [13.831376331255797, 8.53781049371163], [13.812640368061993, 9.406322571152053], [13.607759847190183, 10.205419346701756], [13.402879326318374, 11.004516122251466], [13.19799880544656, 11.803612897801173], [12.99311828457475, 12.60270967335088], [12.788237763702941, 13.401806448900587], [12.58335724283113, 14.200903224450293], [12.378476721959318, 15], [11.426286204885525, 14.954223876658476], [10.474095687811731, 14.908447753316947], [9.521905170737938, 14.862671629975427], [8.569714653664144, 14.8168955066339], [7.617524136590351, 14.771119383292376], [6.665333619516556, 14.725343259950849], [5.7131431024427615, 14.679567136609325], [4.760952585368969, 14.6337910132678], [3.808762068295175, 14.588014889926276], [2.8565715512213803, 14.542238766584749], [1.9043810341475864, 14.496462643243225], [0.9521905170737933, 14.450686519901703]];

const constraints = [[171, 172], [172, 173], [173, 174], [174, 175], [175, 176], [176, 177], [177, 178], [178, 179], [179, 180], [180, 181], [181, 182], [182, 183], [183, 184], [184, 185], [185, 186], [186, 187], [187, 188], [188, 189], [189, 190], [190, 191], [191, 192], [192, 193], [193, 194], [194, 195], [195, 196], [196, 197], [197, 198], [198, 199], [199, 200], [200, 201], [201, 202], [202, 203], [203, 204], [204, 205], [205, 206], [206, 207], [207, 208], [208, 209], [209, 210], [210, 211], [211, 212], [212, 213], [213, 214], [214, 215], [215, 216], [216, 217], [217, 218], [218, 219], [219, 220], [220, 221], [221, 222], [222, 223], [223, 224], [224, 225], [225, 226], [226, 227], [227, 228], [228, 229], [229, 230], [230, 231], [231, 232], [232, 233], [233, 234], [234, 235], [235, 171]];

const del = Delaunator.from(positions);
const con = new Constrainautor(del);
con.constrainAll(constraints);

const { triangles } = del;
console.log('Triangles:');
for(const tri of triangles)
  console.log('[' + tri.join(',') + ']');

In the above, the regular integer coordinates are basically an internal grid, and all the irregular ones are points around those.
The constraints are basically on those irregular points that form an outer polygon.
Several edge constraints are part of the convex hull (and thus with outer half-edge -1). But the samples you showed at https://observablehq.com/@fil/polygon-triangulation seem to have such conditions (except they seem to never have internal points).

Debugging what's happening in the code, it fails at constrainOne(segP1, segP2) with segP1=173 and segP2=174, which is the third constraint edge.
There, it ends up with conEdge=1311, and it enters your posterior loop:

while(edg !== -1){
	// edg is the intersecting half-edge in the triangle we came from
	// adj is now the opposite half-edge in the adjacent triangle, which
	// is away from segP1.
	const adj = del.halfedges[edg], // cross diagonal
		bot = prevEdge(edg),
		top = prevEdge(adj),
		rgt = nextEdge(adj);

with adj=-1, and then the rest goes crazy because we're then indexing the typed arrays at negative indices, which gives undefined, and eventually NaN with the rest of the operations ... and we're stuck.

I'm not clear exactly what that part of the algorithm is, but I wonder if the loop should not check that adj is not -1.

[xionluhnis]

Debugging further, it seems that the constrained segment that makes things fail is near a region with successive boundary edges, for which Delaunator creates a very thin triangle overlapping one constraint edge. This leads to the intersection call:

intersectSegments(p1x, p1y, p2x, p2y, p3x, p3y, p4x, p4y)
with

p1x: 0.219566671390259
p1y: 12.6197739054608
p2x: 0.3293500070853885
p2y: 11.727205659911116
p3x: 0.5489166784756471
p3y: 9.942069168811738
p4x: 0.1097833356951295
p4y: 13.512342151010488

visualized there: https://www.geogebra.org/calculator/zayn3rkp

This returns as intersecting, and then the opposite edge opp from segP1 is used for the further propagation, i.e. your conEdge. This ends up not having an adjacent opposite half-edge, because it's on the convex hull, and then we go into the troubles.

Note that the constraint is p1-p2, whereas p3-p4 is not a constraint. It's constructed by Delaunator as an edge of an outer triangle.

Those very thin triangles don't matter in my code, because I rule out all very thin triangles in the post-processing I do (since I don't need triangles outside my enclosing polygon). But I wonder if you assumed to keep the convex hull by construction (whereas I do not).

I will check if I can do that processing before using your code (those triangles should all be outside the enclosing constraint polygon by construction because my points are constructed to be far enough from each others).
Though it would be great if such thin triangles were not breaking Constrainautor::constrainOne (assuming it's an issue with the thin triangle, which it may not be ... maybe other cases trigger entering the second loop with a half-edge that doesn't have an adjacent opposite half-edge.
Actually, being able to remove those thin triangles implies also fixing the constrained edges that are not currently correctly set by flipping other edges, so it's not a solution for me because it means implementing what you did, but while discarding the impact of those thin triangles...

[kninnug]

Thank you for finding this issue.

However, is the current implementation making specific assumptions about the point distribution? e.g. that there are no points inside of the constrained shape?

This should not be an issue. The algorithm doesn't care about the shape that the constrained edges forms, other than that they may not intersect each other. Constraining edges on the hull is also not (supposed to be) a problem.

with adj=-1, and then the rest goes crazy because we're then indexing the typed arrays at negative indices, which gives undefined, and eventually NaN with the rest of the operations ... and we're stuck.

It is assumed there that adj can never be -1, since that would mean the constraining segment goes to a point outside the hull, or that there is a hole in the triangulation. Constrainautor does assume that the hull is convex, and that the triangulation has no holes. This is what Delaunator generally guarantees of its output, but I'll add it to the documentation. I can add a check & throw for it as well, but that will probably not solve your problem.

Culling triangles on the hull won't work either, since that will break the Constrainautor too.

It seems that indeed Constrainautor.intersectSegments is not robust enough and breaks on smallish line segments. You can try a more robust intersect test implementation, for example robust-segment-intersect by Mikola Lysenko by overriding Constrainautor#intersectSegments:

import robustIntersect from 'robust-segment-intersect';
class RobustConstrainautor extends Constrainautor {
	intersectSegments(p1, p2, p3, p4){
		// ignore when the segments share an end-point
		if(p1 === p3 || p1 === p4 || p2 === p3 || p2 === p4){
			return false;
		}
		
		const pts = this.del.coords;
		
		return robustIntersect(
			[pts[p1 * 2], pts[p1 * 2 + 1]],
			[pts[p2 * 2], pts[p2 * 2 + 1]],
			[pts[p3 * 2], pts[p3 * 2 + 1]],
			[pts[p4 * 2], pts[p4 * 2 + 1]]
		);
	}
}

const del = Delaunator.from(...),
	con = new RobustConstrainautor(del);
con.constrainAll(...);

The above code solves the issue (on my end), but I am somewhat hesitant to switch over to that implementation as it will incur a dependency on several packages, and has a small but measurable performance impact.

[xionluhnis]

Thanks for confirming the problem and finding a solution!
I understand that the choice of segment intersection algorithm is a tricky one and your option of switching to an external one while extending the class sounds good to me.

While the "throw" may not solve my problem directly, having a guard in the code that throws an error upon detecting such invalid case would be good. It's definitely better than an infinite loop when debugging.