Edge-case bugs in boolean operations

[pearcetm]

@lehni First of all thank you for this great library! It is an integral part of a project I've been working on for quite a while now, osd-paperjs-annotation, which integrates paper.js with OpenSeadragon - there's a demo page here. The project is still an active work-in-progress but it is already proving very useful in a number of other application I'm building.

I've come across a couple of issue with boolean operations (e.g. subtract, intersect). One of them I've tracked down and come up with a fix for, and the other I'm still working on sorting out. Both tend to affect complicated paths like those created by the magic wand tool (which you can see in the demo page for reference).

First issue: float truncation - workaround identified

The first one crops up when large coordinate values are used (in the 10,000-20,000+ range, typically). When the values are this big, the epsilon that gets added to the coordinate can get truncated due to floating point precision, which disrupts the winding calculations. My workaround for this is to test whether the bounds of the paths are large enough to cause this issue, and if so, translate and scale the paths as needed to reduce the size of the coordinates, then do the boolean operation, and then undo the transformation. I've implemented this as a monkey patch, rather than a direct modification of the library code itself at this point. I'm posting the code here in case it is useful to others. If you approve of the approach, I can try putting a PR together to address this more formally within this library itself.

    // Monkey patch the paper.js boolean operations to account for issues with floating point math
    // when large coordinate values are used (1000 is an empiric value that seems to work reliably)
    const funcs = ['unite', 'intersect', 'subtract', 'exclude', 'divide'];
    for(const func of funcs){
        const original = paper.PathItem.prototype[func];
        paper.PathItem.prototype[func] = function(){
            const path = arguments[0],
                    numericThreshold = 1000, // empiric
                    b1 = this.getBounds(),
                    b2 = path.getBounds(),
                    l = Math.min(b1.x, b2.x),
                    r = Math.max(b1.x + b1.width, b2.x + b2.width),
                    t = Math.min(b1.y, b2.y),
                    b = Math.max(b1.y + b1.height, b2.y + b2.height);

            if(l > -numericThreshold &&
                r < numericThreshold &&
                t > -numericThreshold &&
                b < numericThreshold ){
                // Our bounds are within the limit: no need to translate or scale, just call the original function
                return original.apply(this, arguments);
            }
            // One or more of our bounds is out of range
            // Calculate whether we need to scale or just translate
            const w = r - l,
                    h = b - t,
                    scaleX = Math.pow(2, Math.ceil(Math.log2(w/(2*numericThreshold)))),
                    scaleY = Math.pow(2, Math.ceil(Math.log2(h/(2*numericThreshold)))),
                    scale = Math.max(scaleX, scaleY),
                    center = new paper.Point((l + r)/2, (t + b)/2),
                    offset = new paper.Point(-Math.round(center.x), -Math.round(center.y));            
            
            if(scale > 1){
                // we need to scale the path(s) to make them fit within our numeric bounds
                this.scale(1/scale, center);
                if(path !== this){
                    path.scale(1/scale, center);
                }
            }

            // translate the path(s) by the offset
            this.translate(offset);
            if(path !== this){
                path.translate(offset);
            }

            const result = original.apply(this, arguments);

            // restore the path(s)
            this.translate(offset.multiply(-1));
            result.translate(offset.multiply(-1));
            if(path !== this){
                path.translate(offset.multiply(-1));
            }

            if(scale > 1){
                // reset the scale back to the original values
                this.scale(scale, center);
                result.scale(scale, center);
                if(path !== this){
                    path.scale(scale, center);
                }
            }
            
            return result;
        }
    }

Second issue: Many overlapping segments - unresolved

The second issue, which I'm still working on tracking down, comes up when two paths have a lot of overlapping segments - like when PathItem.subtract is used to subtract one path from another, the borders align perfectly. The problem is when the result of the subtract operation and the original path are then tested against each other.

For example, with paths A and B, let R = A.subtract(B). R and B are now non-overlapping, with a shared stretch of border. R.subtract(B) should just return R again (no overlap = nothing to subtract), but sometimes the result is entirely empty. It depends a lot on the exact coordinates of the polygon coordinates. See this sketch for an example.

Note that in the screenshot, the log message for the bottom three tests are truncating the offset: it is actually new Point(0.0000001, 0) (not 0,0).

I'm looking for any pointers on where the issue might be arising, and how to best test for these cases. If I can identify the cause I'll also try to put in a PR for this.

Thanks again for the great library!

[pearcetm]

I've made some progress narrowing this down further. The problem appears to arise when the two curves contain are two nearly identical points, but which are not exactly identical (most commonly this would be due to floating point precision in an earlier calculation). This can lead to super-short segments/curves being added to a path during getCurveIntersections. Because of accumulating precision errors while calculating the intersection, the location of new point that is added can actually be less precise than the original points (this point may or may not be relevant). But the two segments are now so close together that they're less than the epsilon value used in the winding calculations which can mess things up.

Part of this calculation depends on Curve.getTimeOf to return a "time" value along the curve used to determine the "best" point for the intersection. If the "intersection" falls within the epsilon tolerance of an endpoint of the curve, it should be taken to be the same as the endpoint (i.e. time equal to 0 or 1). If not, cubic roots are solved and tested against the curve. This check is currently done in two phases: if the point falls outside of the endpoints with epsilon = 1e-12 we proceed to solve the cubic. If a solution is found, it is returned. If NO solution is found, the endpoints are then checked with geomEpsilon = 1e-7.

paper.js/src/path/Curve.js

Lines 701 to 730 in f1f02cc

	getTimeOf: function(v, point) {
	// Before solving cubics, compare the beginning and end of the curve
	// with zero epsilon:
	var p0 = new Point(v[0], v[1]),
	p3 = new Point(v[6], v[7]),
	epsilon = /*#=*/Numerical.EPSILON,
	geomEpsilon = /*#=*/Numerical.GEOMETRIC_EPSILON,
	t = point.isClose(p0, epsilon) ? 0
	: point.isClose(p3, epsilon) ? 1
	: null;
	if (t === null) {
	// Solve the cubic for both x- and y-coordinates and consider all
	// solutions, testing with the larger / looser geometric epsilon.
	var coords = [point.x, point.y],
	roots = [];
	for (var c = 0; c < 2; c++) {
	var count = Curve.solveCubic(v, c, coords[c], roots, 0, 1);
	for (var i = 0; i < count; i++) {
	var u = roots[i];
	if (point.isClose(Curve.getPoint(v, u), geomEpsilon))
	return u;
	}
	}
	}
	// Since we're comparing with geometric epsilon for any other t along
	// the curve, do so as well now for the beginning and end of the curve.
	return point.isClose(p0, geomEpsilon) ? 0
	: point.isClose(p3, geomEpsilon) ? 1
	: null;
	},

My solution is to do the initial check with the geomEpsilon instead of plain epsilon. This way we shouldn't end up with time values within geomEpsilon of the endpoints, which was possible before (and which was causing the issues described above):

paper.Curve.getTimeOf = function(v, point){

    // Before solving cubics, compare the beginning and end of the curve
    // with zero epsilon:
    var p0 = new paper.Point(v[0], v[1]),
        p3 = new paper.Point(v[6], v[7]),
        geomEpsilon = 1e-7,
        t = point.isClose(p0, geomEpsilon) ? 0
            : point.isClose(p3, geomEpsilon) ? 1
            : null;
    if (t === null) {
        // Solve the cubic for both x- and y-coordinates and consider all
        // solutions, testing with the larger / looser geometric epsilon.
        var coords = [point.x, point.y],
            roots = [];
        for (var c = 0; c < 2; c++) {
            var count = paper.Curve.solveCubic(v, c, coords[c], roots, 0, 1);
            for (var i = 0; i < count; i++) {
                var u = roots[i];
                if (point.isClose(paper.Curve.getPoint(v, u), geomEpsilon))
                    return u;
            }
        }
    }
    
    
    return t;
    
}

@lehni Thoughts on this change? It appears to be working well so far for my application, though I'll continue to test. @iconexperience I see you've done a lot of great work on related issues - any chance you could take a look?

[pearcetm]

I think issue 2 may be related to #1261, and the solution is similar to that proposed by @iconexperience in #1261 (comment). @lehni responded that looking initially at points within GEOMETRIC_EPSILON would cause problems in some cases (#1261 (comment)). However, I'm not sure if that discussion related to the same place within the code that I'm proposing to check the endpoints here (Curve.getTimeOf), or if it was elsewhere within the actual boolean operations, and I'm not sure if that concern would still exist given the other changes since the prior issue.