path.unite() fails to resolve intersections

[iconexperience]

In this example calling unite() on a path fails to resolve the self intersections of the path. The result is a path that still has self intersections.

var path = new Path({segments:[[349.31714952528665, 176.97192499841802, 2.0667935742653185, -10.850610233997372, 0, 0], [308.23418424120047, 394.39737504104323, 0, 0, -0.3852007263769224, 2.0386172397749647], [273.60284923271973, 415.7138856106701, 9.194683617471242, 1.7373535056759692, 0, 0], [340.6967045975081, 152.3115389894092, -4.006231967752797, -3.6687868382265094, 5.694026309937783, 5.214418174126138], [349.72123243229134, 170.85618880187394, -0.36146393999291604, -6.612268801346318, 0.1110858449889065, 2.0320961120115726], [349.31714952528654, 176.97192499841861, 0.3852007263769224, -2.038617239775249, -2.0667935742656027, 10.850610233997315], [333.4126425432937, 153.58289639999975, -10.850610233997315, -2.0667935742653754, 10.850610233997372, 2.0667935742653754]], closed:true});
var path2 = path.unite();
console.log(path2.className);
path2.translate(200, 0);
path.strokeColor = "green";
path2.strokeColor = "red";

Haha, so just one more recursion is needed? Pfff! Are you happy otherwise with the fine-tuned values in addCurveIntersections()?

[iconexperience]

What now? What should the new value for recursions be? 28?

Also, while we are at this we can change the split condition from oldTDiff > 0.5 && tDiff > 0.5 to tDiff > 0.8. This is the recommendation in the original algorithm and does not seem to make any difference in our tests. The parameter oldTDiff is then not needed anymore.

[lehni]

Yes that's what I meant with my question (tDiff > 0.8)... I remember you saying something about that elsewhere, and shall change it right now. Regarding the recursion: What if we set it to 26, and keep testing this first if we encounter new situations? Or should we go 32 right away? I love powers of two.

[iconexperience]

Powers of two are nice. But I prefer 26, because you have to keep in mind that we have a binary tree here. If we have a freak curve pair that fills the whole tree, changing the recursion to 32 would mean 256 times more calls to addCurveIntersections(). This would mean 4,294,967,296 calls as a theoretical limit. Although I have not encountered such a curve pair anymore, one case would probably kill most browsers.

So let's set it to 26.

[lehni]

Yes I was wondering about exactly that. 26 it is! We don't need more crashing web browsers.

[iconexperience]

This was actually quite an interesting issue, I learned a lot.

[lehni]

Super!

[lehni]

For reference, here a situation that I thought could still make problems due to exactly identical tangents, but it doesn't because we're using tMin / tMax!

[iconexperience]

Yep, using tMin/tMax was a clever move. And if we ever run into issues that require curvature, it's all prepared.

[lehni]

It would have been clever if I knew why it solved a few edge cases. I guess I had a hunch. But now we know! : )

[iconexperience]

Here are some metrics for fat line clipping on this issue.

• The function addCurveIntersections() is called 126 times for the intersection. This is far less than I would have expected.
• The curve gets split 38 times

At least for this intersection the fear that we will end up with a huge call tree is unfounded. That's a great relief.

[iconexperience]

Another thing that I noticed is that there seem to be many cases that reach the maximum recursion number without any curve splitting at all. These cases do not yield in an intersection and occur if two curves are connected tangentially at an end point, like in this example:

These curve pairs do not seem to be a problem, but it would certainly be nice to recognize them earlier so we do not have to go to the maximum number of recusions if there is no intersection anyway.

Let's keep this in mind as a potential performance enhancement.

[iconexperience]

Here are some more metrics. I ran my large test case and measured the number of calls to addCurveIntersections().

• In only 36 cases (of about 300,00) addCurveIntersections() was called more than 100 times.
• The maximum number of calls was 672, the second largest number of calls was 256.

Overall this is another great relief. The call tree does not grow as much as I had suspected.

[lehni]

Very interesting! Regarding https://github.com/paperjs/paper.js/issues/899#issuecomment-171200668: The handle bounds of these curves touch merely. If we exclude touch, which we can because we handle start- and end-points separately, we should get rid of all of these right away, correct?

[iconexperience]

Correct, but how do you exclude touch in fat line clipping while ensuring that any other intersections will be found? That is the big question.

If we find an easy way it seems like we can also raise the maximum recursions value to 32. I just ran another test with 1,000,000 curve pairs and maximum recursions set to 32. The number of calls to addCurveIntersections() did not increase, I still have only one case with more than 500 calls.

The problem with raising the max recursions value is that there are many cases of tangentially touching curves (just think of an ellipse). Since these cases always reach maximum recursion without finding an intersection, raising the number of recursion directly affects performance.

[lehni]

Hmm but that suggestion would kill this again: https://github.com/paperjs/paper.js/issues/878#issuecomment-168153188

[iconexperience]

Another option would be to wait a certain number of recursions (until the curves are almost stright) and then check if the handle bounds touch only in one point. I think this can be done after just a few recursions.

[lehni]

I just switched the code linked above to compare against epsilon without the equal sign in <= and >= checks. Now the easy solution could be to introduce a value boundsEpsilon = 0, and detect the case where the bounding boxes have either no with or no height, in which case it is set to epsilon.

[iconexperience]

Sorry, I do not understand :)

[lehni]

Let me just test it : )

[lehni]

Something like this, I was thinking:

    var epsilon = /*#=*/Numerical.GEOMETRIC_EPSILON,
        boundsEpsilon = 0,
        c1p1x = v1[0], c1p1y = v1[1],
        c1p2x = v1[6], c1p2y = v1[7],
        c2p1x = v2[0], c2p1y = v2[1],
        c2p2x = v2[6], c2p2y = v2[7],
        // 's' stands for scaled handles...
        c1s1x = (3 * v1[2] + c1p1x) / 4,
        c1s1y = (3 * v1[3] + c1p1y) / 4,
        c1s2x = (3 * v1[4] + c1p2x) / 4,
        c1s2y = (3 * v1[5] + c1p2y) / 4,
        c2s1x = (3 * v2[2] + c2p1x) / 4,
        c2s1y = (3 * v2[3] + c2p1y) / 4,
        c2s2x = (3 * v2[4] + c2p2x) / 4,
        c2s2y = (3 * v2[5] + c2p2y) / 4,
        min = Math.min,
        max = Math.max,
        c1x1 = min(c1p1x, c1s1x, c1s2x, c1p2x),
        c1x2 = max(c1p1x, c1s1x, c1s2x, c1p2x),
        c2x1 = min(c2p1x, c2s1x, c2s2x, c2p2x),
        c2x2 = max(c2p1x, c2s1x, c2s2x, c2p2x),
        c1y1 = min(c1p1y, c1s1y, c1s2y, c1p2y),
        c1y2 = max(c1p1y, c1s1y, c1s2y, c1p2y),
        c2y1 = min(c2p1y, c2s1y, c2s2y, c2p2y),
        c2y2 = max(c2p1y, c2s1y, c2s2y, c2p2y);
    if (c1x1 === c1x2 || c1y1 === c1y2 || c2x1 === c2x2 || c2y1 === c2y2)
        boundsEpsilon = epsilon;
    if (!(  c1x2 + boundsEpsilon > c2x1 &&
            c1x1 - boundsEpsilon < c2x2 &&
            c1y2 + boundsEpsilon > c2y1 &&
            c1y1 - boundsEpsilon < c2y2))
        return locations;

[lehni]

Oh this doesn't work for edge cases yet (e.g. https://github.com/paperjs/paper.js/issues/568#issuecomment-63792068), because this is only the bounds check to be used before calling add*Intersections(), the code that checks start- end points needs to be performed still...

[lehni]

This works for me:

_getIntersections: function(v1, v2, c1, c2, locations, param) {
    if (!v2) {
        // If v2 is not provided, search for a self-intersection on v1.
        return Curve._getSelfIntersection(v1, c1, locations, param);
    }
    // Avoid checking curves if completely out of control bounds. As
    // a little optimization, we can scale the handles with 0.75
    // before calculating the control bounds and still be sure that
    // the curve is fully contained.
    var epsilon = /*#=*/Numerical.GEOMETRIC_EPSILON,
        c1p1x = v1[0], c1p1y = v1[1],
        c1p2x = v1[6], c1p2y = v1[7],
        c2p1x = v2[0], c2p1y = v2[1],
        c2p2x = v2[6], c2p2y = v2[7],
        // 's' stands for scaled handles...
        c1s1x = (3 * v1[2] + c1p1x) / 4,
        c1s1y = (3 * v1[3] + c1p1y) / 4,
        c1s2x = (3 * v1[4] + c1p2x) / 4,
        c1s2y = (3 * v1[5] + c1p2y) / 4,
        c2s1x = (3 * v2[2] + c2p1x) / 4,
        c2s1y = (3 * v2[3] + c2p1y) / 4,
        c2s2x = (3 * v2[4] + c2p2x) / 4,
        c2s2y = (3 * v2[5] + c2p2y) / 4,
        min = Math.min,
        max = Math.max,
        c1x1 = min(c1p1x, c1s1x, c1s2x, c1p2x),
        c1x2 = max(c1p1x, c1s1x, c1s2x, c1p2x),
        c2x1 = min(c2p1x, c2s1x, c2s2x, c2p2x),
        c2x2 = max(c2p1x, c2s1x, c2s2x, c2p2x),
        c1y1 = min(c1p1y, c1s1y, c1s2y, c1p2y),
        c1y2 = max(c1p1y, c1s1y, c1s2y, c1p2y),
        c2y1 = min(c2p1y, c2s1y, c2s2y, c2p2y),
        c2y2 = max(c2p1y, c2s1y, c2s2y, c2p2y);

    function checkBounds(epsilon) {
        return  c1x2 + epsilon > c2x1 && c1x1 - epsilon < c2x2 &&
                c1y2 + epsilon > c2y1 && c1y1 - epsilon < c2y2;
    }

    if (!checkBounds(epsilon))
        return locations;
    // Now detect and handle overlaps:
    var overlaps = Curve.getOverlaps(v1, v2);
    if (overlaps) {
        for (var i = 0; i < 2; i++) {
            var overlap = overlaps[i];
            addLocation(locations, param,
                v1, c1, overlap[0], null,
                v2, c2, overlap[1], null, true);
        }
        return locations;
    }

    var straight1 = Curve.isStraight(v1),
        straight2 = Curve.isStraight(v2),
        straight = straight1 && straight2,
        before = locations.length;
    // If both bounding boxes have with and height, perform a stricter
    // bounds check, excluding merely touching at their borders.
    if ((c1x1 === c1x2 || c1y1 === c1y2 ||
         c2x1 === c2x2 || c2y1 === c2y2) || checkBounds(0)) {
        // Determine the correct intersection method based on whether
        // one or curves are straight lines:
        (straight
            ? addLineIntersection
            : straight1 || straight2
                ? addCurveLineIntersections
                : addCurveIntersections)(
                    v1, v2, c1, c2, locations, param,
                    // Define the defaults for these parameters of
                    // addCurveIntersections():
                    // tMin, tMax, uMin, uMax, reverse, recursion
                    0, 1, 0, 1, 0, 0);
    }
...

[lehni]

The crucial part being:

    // If both bounding boxes have with and height, perform a stricter
    // bounds check, excluding merely touching at their borders.
    if ((c1x1 === c1x2 || c1y1 === c1y2 ||
         c2x1 === c2x2 || c2y1 === c2y2) || checkBounds(0)) {
        // Determine the correct intersection method based on whether
        // one or curves are straight lines:

The question is: Are these the right conditions for when to not perform this additional, more strict check?

[iconexperience]

So if you have two vertical lines with a distance of 1e-15, no intersections will be found?

[lehni]

In that situation, the same intersections as currently are found, because the first call to checkBounds() is the same as right now, and the 2nd, more strict call is only performed if both bounding boxes have area > 0, excluding cases as the one in https://github.com/paperjs/paper.js/issues/899#issuecomment-171200668 or https://github.com/paperjs/paper.js/issues/568#issuecomment-63792068 before calling the fat-line clipping code. Afterwards, the beginning and end points are checked against each other regardless of the outcome of this stricter check.

[iconexperience]

But it would not work for these cases, right?

[lehni]

Yeah you're right : /

[lehni]

And for some strange reason that I don't quite comprehend it has lead to slowdowns in the BooleanOperations.html example.

[iconexperience]

I think I found a very nice way to solve this. If we change the code in Curve.js (lines 1517 - 1526) from

            } else if (tDiff > 0) { // Iterate
                addCurveIntersections(v2, v1Clip, c2, c1, locations, param,
                        uMin, uMax, tMinNew, tMaxNew, !reverse, recursion);
            } else {
                // Curve 1 has converged to a point. Since we cannot construct a
                // fat-line from a point, we dismiss this clipping so we can
                // continue clipping curve 2.
                addCurveIntersections(v2, v1, c2, c1, locations, param,
                        uMin, uMax, tMin, tMax, !reverse, recursion);
            }

to

            } else { // Iterate
                addCurveIntersections(v2, v1Clip, c2, c1, locations, param,
                        uMin, uMax, tMinNew, tMaxNew, !reverse, recursion);
            }

everything seems to work. In my big test case this reduces the calls to addCurveIntersecticions() from 4940290 to 1013953! I was probably thinking too complicated when I created #863.

Could you check how this affects performance?

[lehni]

Wowowow! Will check right now : )

[lehni]

Very nice. In the particular case in BooleanOperations.html I've been testing with, it's a little faster now, maybe 5%. I also ran all the usual tests and couldn't find any problem with this.

[lehni]

Some nice code simplifications resulted from this again.

[iconexperience]

If I now increase maximum recursion to 32, the number of calls rises marginally from 1013953 to 1019560, but one new glitch (7751) appears in my test case. This was probably hidden before by a missed intersection. I can investige this later, maybe tonight.

So fat line clipping seems to be very robust, with a clean code. The two remaining minor issues that I see are

• Should we further raise the maximum recursions? Are ther cases where this is necessary? Can this lead to a massive call tree?
• Are there any cases that need special handling due to the cange in e0d2d0d? I am specifically thinking of curves with vertical handles, because if we call getSignedDistance() with a vector of length zero, it returns the distance as if the line was vertical.

[iconexperience]

I have started to investigate the curve pairs that create many calls to addCurveIntersecticions().

Initially I found this curve pair which creates 672 calls to addCurveIntersecticions(). As you can see, curve2 is almost identical to the first part of curve1, but of course not completely, because otherwise our function getOverlaps() would have filtered this curve pair.

So it seems like the problematic curve pairs are the ones where one curve is almost identical to part of the other curve. Therefore I decided to generate these curves artifically. I take one curve and create another curve that is exactly part of the first curve. Then I move one end point of the second curve just slightly, so that getOverlaps() does not filter these curves out. Here is the code how I do it:

        var crv1 = CubicBezierUtils.createRandomCurve(100);
        var crv2 = crv1.clone().divide(Math.random(), true);
        crv2.segment2.point.x += 1.0001 * Numerical.GEOMETRIC_EPSILON;
        crv2.segment2.point.y += 1.0001 * Numerical.GEOMETRIC_EPSILON;

If you call crv1.getIntersections(crv2) on these two curves, things start to get very interesting. These curve pairs usually create more than 50,000 calls to addCurveIntersecticions(), in single cases up to 300,000 calls. And there are cases when more than the theoretically possible 9 intersections are found. One case had 21 intersections.

Another thing that I found is that the number of calls to crv1.getIntersections(crv2) and the maximum number of found intersections rises with larger curve sizes.

Now this sounds quite grim, but let's keep in mind that these curve pairs should be the worst of the worst cases that we can encounter, and I have never seen any similar pair in a real case scenario.

But I want to look further into these edge cases and see if there is a way to recognize and handle them in a better way.

[lehni]

Interesting indeed! I think this deserves a new issue. Could you open one with the above content copied to it? Does this mean that getOverlaps()is too strict? If more than 9 overlaps are found, it should mean that the curve are actually identical, no? Otherwise, the epsilons in addCurveIntersecticions() are too permissive?