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";

[lehni]

Uh oh, here we go again : )

Looks like it does find all the intersections... Next I'm looking at resolveCrossings()

[lehni]

Ok, so this intersections get wrongly filtered out as not a crossing:

The problem is in isTouching(): https://github.com/paperjs/paper.js/blob/develop/src/path/CurveLocation.js#L384

Or better, the way it gets used in isCrossing() currently to detect the case where two touchings curves are not an intersection. This is of course wrong...: https://github.com/paperjs/paper.js/blob/develop/src/path/CurveLocation.js#L421

[lehni]

But even when I I deactivate this check, (which doesn't seem to add any issues btw), it still gets wrongly flagged as not a crossing by the code further down. The reason being all the tangents pointing in the same direction. This is where this comment here becomes important:

https://github.com/paperjs/paper.js/blob/develop/src/path/CurveLocation.js#L476

// NOTE: VectorBoolean has code that slowly shifts these points inwards
// until the resulting tangents are not ambiguous. Do we need this too?

The question then also is: Is there a better way to do this? E.g. use a solver to find the first location where the slope changes enough?

[iconexperience]

Isn't this exactly what Cary Clark talks about here?

https://www.youtube.com/watch?v=OmfliNQsk88&feature=youtu.be&t=1372

[lehni]

Yes it appears to be. And it's quite interesting... It sounds like he's calculating the winding contributions differently to what we're doing, based on just the found intersections and the curves going through it, while we're calling getWinding(). I guess his approach is faster.

[lehni]

Oh but after he discussed about propagation, which is what we are doing. I should listen to this properly again.

[iconexperience]

Unfortunately his solution for ordering the curves is for quadratic curves, and he says at the end: "There may be a better one [solution]". So we cannot use this.

[lehni]

Yeah we shouldn't switch now. = ) I also don't feel that our approach is too slow, so that's fine.

[lehni]

But if the two things would happen together (finding intersections and calculating winding contributions) we could probably use the winding numbers to figure out the answer. I think what @kuribas is doing here is smarter in that respect and would help with that: #761 (comment)

But we're too far down our current road probably.

[iconexperience]

Uuuh, I do not want to think about this now :) Just for your information, I found this issue only because I made a mistake when playing with offsetting, which accidentially caused a loop. So this edge case should have never happened (but it did).

[lehni]

hehe, as they always do!

So here is the code mentioned above: https://github.com/adamwulf/vectorboolean/blob/master/VectorBoolean/FBContourOverlap.m#L287

I think we can use the isTouching() check to figure out initially if we need to do this at all or not, nothing else. Then if we need to find non-ambiguous tangents, I wonder if there's a better approach than simply shifting the offsets on the curves (in curve length) away from the intersection point until we find a unambiguous pair of tangents? Do you have a suggestions here?

[lehni]

Ha, interesting commit message here. It looks like we're not the only ones out there : ) adamwulf/vectorboolean@582e65d

[iconexperience]

I am not sure if I fully understand this. So far I tried to stay away from that part of the code. But is the use of getCurvatureAt() of any help here?

Here is an example how curvature could help

[lehni]

Well, it's simple: In this case, both involved curves have collinear tangents at the intersection. The intersection happens to be in existing segments, and the segments are smooth, meaning their in- & out-handles are collinear, too.

So I was wrong about isTouching() being involved, since that only gets called when the intersection is within curves, not in segments at their beginnings / ends.

isCrossing() then gets all the involved tangents and their angles, and compares these to figure out if it is a crossing or not. But since in this example, all the tangents are identical, they are ambiguous.

The only solution that I can think of right now is the same as in VectorBoolean: Slowly moving away until we find unambiguous tangents, meaning they are not collinear anymore, the curves unveil in which direction they are turning after they touch in the intersection. The offset could be a fraction of the total curve length of the shorter curve perhaps, rather than a fixed amount as in their code.

I was wondering if there could be a better solution, at the same time, this happens so rarely that I guess I shouldn't care?

[iconexperience]

I initially forgot the example in my previous post. Wouldn't this be better than slowly moving away from the point?

[lehni]

Oh wow, yes! Let me look into this a bit! : )

[lehni]

Hmm not really... There is no distinction between this sketch and that sketch in terms of curvature values.

[iconexperience]

I wouldn't expect there to be a difference, because in your sketches you only look at the second curve (the one to the right), which is the same.

Look at this example. To be honest, it surprises me a bit that the curvature is different at the same point, depending on which curve you look at.

[lehni]

Yes, you're of course right! I shall give this a go, good thinking.

And I guess that's just how curvature works?

[lehni]

But what about this situation?

They both change signs, and they are actually crossing. Does the amount of curvature then decide? Can we be certain?

[iconexperience]

Doesn't the curvature indicate how fast the curve is crawling away from the tangent and to which direction? I think it should work in your example, but I cannot give you the correct algorithm at the moment.

[lehni]

Yea I was thinking the same. It's one for another day, with a fresher head : ) But definitely solvable.

[iconexperience]

As a first step I have created This little sketch that shows the relation of the curvature at t=0 and the direction in which two curves with identical tangent separate.

The curve with the larger (signed) curvature will always curl away to the right, from the curve with the smaller curvature value. In the example below, the green curve has the larger curvature value.

[iconexperience]

And just as a reminder, the angle in getAngle() is calculated clockwise:

[lehni]

Yes, because: radius = 1 / curvature for the circle that can be fit tightly into the location that the curvature was calculated for: updated sketch

[lehni]

But you're aware that you're comparing signed curvatures, yes?

[iconexperience]

You are correct, the console output is wrong. It should read "red/green curve is to the right". I updated the sketch in my comment above.

[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 #899 (comment): 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: #878 (comment)

[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. #568 (comment)), 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 #899 (comment) or #568 (comment) 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?