Improve classification and handling of cubic curves

[lehni]

As @deanm suggested here: https://github.com/paperjs/paper.js/issues/1074#issuecomment-239798954:

We could benefit from a better system of classification of cubic curves. This issue is created to track this progress separately.

It sounds like following Loop and Blinn should get us there quickly:

http://research.microsoft.com/en-us/um/people/cloop/LoopBlinn05.pdf

[lehni]

Some more information by @deanm here: https://github.com/paperjs/paper.js/issues/773#issuecomment-212840326

[deanm]

I've only implemented the Loop Blinn analysis, but I just wanted to mention (not recommending) another paper that presents an alternative basis that gives you the same sort of classification:

https://arxiv.org/pdf/1605.08669.pdf

If nothing else, page 16 presents some nice examples of the different cases.

[lehni]

@hkrish just told me about this interesting paper here that we should also look into:

Stephen Vincent. Fast Detection of the Geometric Form of Two-Dimensional Cubic Bézier Curves. Journal of Graphics Tools, 7(3):43-51, 2002 (see https://github.com/rougier/gl-bezier/blob/master/bezier_type.py)

[deanm]

Yes, I know the paper, it's from Adobe... It predates Loop and Blinn by a few years. In my opinion Loop and Blinn's analysis is definitely simpler, probably faster, etc. If I remember correctly they are more or less achieving the same thing, what Blinn did was look at the problem by extending it to homogenous coordinates which I think makes everything work out in a way that is simpler, and more general (less cases). The Loop Blinn classification is extremely simple to implement. Also I'm not sure that Vincent's approach gives an easy way to calculate the parameters for a loop double point, for example?

[deanm]

Here is some simple JavaScript for the Loop Blinn classification. Depending on how you plan to use it, there should probably be a few numerical adjustments, basically rounding towards zero in a few places depending what you want out of the analysis. If you need the double point parameter for a loop / cusp, I can also give you that code, it can be directly calculated from d1/d2/d3.

// int bez3_classify_cubic(x0, y0, x1, y1, x2, y2, x3, y3)
//
// Return the type of cubic Bézier via discriminant classification.  See:
//   Resolution Independent Curve Rendering using Programmable Graphics Hardware
//   GPU Gems 3 chapter 25.
//
// NOTE: If, for example, a curve is classified as loop, it doesn't mean that
//       is has a loop over the Bézier interval of 0 .. 1, it just means that
//       the underlying polynomial is of the "loop type".  So this characterizes
//       the type of the cubic polynomial, not necessarily the type of Bézier.
//
// Returns:
//   0 - Point       (all points equal)
//   1 - Line        (d1 == d2 == d3 == 0)
//   2 - Quadratic   (d1 == d2 == 0)
//   3 - Serpentine  (discr > 0)
//   4 - Cusp        (discr == 0)
//   5 - Loop        (discr < 0)
function bez3_classify_cubic(x0, y0, x1, y1, x2, y2, x3, y3) {
  if (x0 === x1 && x0 === x2 && x0 === x3 && y0 === y1 && y0 === y2 && y0 === y3)
    return 0;  // Point

  // Calculate coeffecients of I(s, t), roots of which are inflection points.
  var a1 = y0*(x2 - x3) + x3*y2 - x2*y3 + x0*(y3 - y2),
      a2 = x0*y3 - x3*y0 + y1*(x3 - x0) + x1*(y0 - y3),
      a3 = x1*y0 - x0*y1 + x2*(y1 - y0) + y2*(x0 - x1);  // Error in GPU Gems.

  var d3 = 3*a3, d2 = d3 - a2, d1 = d2 - a2 + a1;

  if (d1 === 0 && d2 === 0) return d3 === 0 ? 1 : 2;  // Line or quadratic.

  // Just for numerics to keep the squares from getting huge.  Theoretically
  // optional.  Could also scale the input points instead.
  var max = Math.max(Math.abs(d3), Math.abs(d2), Math.abs(d1));
  d1 /= max; d2 /= max; d3 /= max;

  var discr = d1*d1 * (3*d2*d2 - 4*d1*d3);
  return discr > 0 ? 3 : discr === 0 ? 4 : 5;
}

[lehni]

@deanm thanks! I was just looking at the WebKit implementation and saw how they were handling rounding. Very interesting. It looks liek Vincent's approach has the advantage of finding the amount of inflection points?

[deanm]

@lehni I posted a table on one of these bugs before, but here it is again:

                        A           B           C           D
                       Func        Cusp        Loop

  Inflection Points     1           -           -           2
              Cusps     -           1           -           -
 Self-Intersections     -           -           1           -
 Curvature Extremum     2          0/2         1/3         3/5

Those categories (A, B, C, D) map to a Loop and Blinn result. So you know the number of inflection points from Loop and Blinn. In the end that is exactly what Loop and Blinn is doing, it's looking at the roots of the equation of the inflection points, the discriminant tells you how many roots are real and whether there are multiplicities, etc.

[deanm]

I hope I don't make any mistakes, but off the top of my head, the mapping between the above table and the cases in Loop and Blinn's 2005 paper. "A" is 3b, "B" is 3a, "C" is 2, and "D" is 1.

I think it's possible the code I posted above doesn't discern between 3a and 3b, but that is a very simple thing to add.

[lehni]

@deanm thank you so much! I shall start putting this all together and see where I get : ) Definitely a whole lot of help.

[deanm]

Btw, again I hope I don't make a mistake, but Loop and Blinn call their case 3b a cusp, but with the cusp at infinity, so I think by a practical definition it is not really a cusp type. Other literature refers to the same case as a serpentine, the important distinction for you probably is there are two cases still, A and D, 1 or 2 inflection points.

[lehni]

I have programmed this as an interactive sketch now, with further distinction between single and double inflections. You can use the mouse to change the curve and see its classification in the log:

Sketch

This shall serve as a playground for further discussion.

[lehni]

@deanm I am getting a whole lot of loops for very normal looking curves, even the initial one. Isn't this wrong?

[lehni]

For the record, here the comparison with your method:

Updated Sketch

[deanm]

Hey Lehni,

Make sure to read the comment at the top of the function : )

// NOTE: If, for example, a curve is classified as loop, it doesn't mean that
//       is has a loop over the Bézier interval of 0 .. 1, it just means that
//       the underlying polynomial is of the "loop type".  So this characterizes
//       the type of the cubic polynomial, not necessarily the type of Bézier.

The classification is for the entire curve, not just the interval from t=0..1. If you extended your curve it (should) form a loop.

[deanm]

And the caterpillar becomes a butterfly... ; )

[lehni]

Yes I was thinking that... But couldn't imagine that particular cure to behave this way. Thanks for illustrating it that beautifully : )

[lehni]

I've quickly ported over the code from https://github.com/rougier/gl-bezier/blob/master/bezier_type.py as well, and it does yield more fine-grained results, which I think are more useful:

Updated Sketch

It also seems that the 'double' / 'single' distinction is more correct there.

[lehni]

And I just added display of inflection points, so it is easier to verify if the results are correct:

Updated Sketch

[deanm]

What do you mean by more fine-grained results? I probably should have been a lot more clear that the code I posted was just an example of Loop and Blinn's analysis, and how straight-forward the math is, but it's not at all meant to be a final implementation, and there are a lot of further things you can get out from Loop and Blinn's approach. I think it's effectively the same concept as Vincent, it's just Loop and Blinn approach the math a bit differently but basically it's the same idea. In the case of Loop and Blinn, once you have d1/d2/d3, it's direct to compute the inflection, double point, cusp, etc parameter values. Once you have the parameter values you can just check them if they are between 0 .. 1 and you know if they are part of the segment or not.

I think the confusion between "double" and "single" is again just about the parameter range. Loop and Blinn analysis is always on the entire underlying curve. So it tells you there are, for example, 2 inflection points, then you can actually check the t values for those inflection points and decide if they are on the segment or not. So I don't think Loop and Blinn is not "correct" in this sense it's just using a different definition. This is because Vincent's analysis is sort of adhoc and more based on the control polygon, where Loop and Blinn's analysis is more on the polynomial of the underlying curve (irrespective of the parameterization).

[lehni]

Apologies for being vague. What I meant is that I can get concrete information about the curve between 0..1 more quickly with Vincent's approach. It is true that I can get the t values for the inflections or the double point more easily from Loop and Blinn, but the calculations involved aren't all that simple, and probably not much faster than simply solving for quadratic roots with our current solver already... I shall do some more tests today.

[lehni]

@deanm I have implemented your suggestion now, and you are right, it does work really beautifully and the computation should be fairly cheap. This works really well in all my tests so far:

Updated Sketch

[lehni]

@deanm thank you for all the suggestions and help! This worked out beautifully.