add section: area under a bézier curve

[jamesgk]

@behdad and I recently came up with a relatively simple way to calculate the area under a cubic (or by simple extension, quadratic) bézier curve. He thinks this would be a great place to describe it.

Basically the equation is this (latex source in alt text):

We derive this equation first by taking the definite integral

Note that this yields a signed area, which in Processing's coordinate system is positive if the area lies to the left of the curve and negative if the area is to the right (reverse of a typical Cartesian system).

We know:

and this equation's derivative:

and so the integral comes out to:

We then subtract excess area between x0 and x3:

or:

Then, factoring out 3 from the resulting coefficients gives us our final equation for A. Here is a method of your BezierCurve class demonstrating this calculation:

  float getArea() {
    if (points.length < 3) return 0;
    float x0 = points[0].x, y0 = points[0].y,
      x1, y1, x2, y2, x3, y3;
    if (points.length == 4) {
      x1 = points[1].x; y1 = points[1].y;
      x2 = points[2].x; y2 = points[2].y;
      x3 = points[3].x; y3 = points[3].y;
    } else if (points.length == 3) {
      x1 = points[0].x * 1 / 3 + points[1].x * 2 / 3;
      y1 = points[0].y * 1 / 3 + points[1].y * 2 / 3;
      x2 = points[2].x * 1 / 3 + points[1].x * 2 / 3;
      y2 = points[2].y * 1 / 3 + points[1].y * 2 / 3;
      x3 = points[2].x; y3 = points[2].y;
    } else {
      return 0; // just return zero if order > 3
    }
    return (
        x0 * (      - 2*y1 -   y2 + 3*y3)
      + x1 * ( 2*y0        -   y2 -   y3)
      + x2 * (   y0 +   y1        - 2*y3)
      + x3 * (-3*y0 +   y1 + 2*y2       )
    ) * 3 / 20;
  }

Hopefully you think this is useful, and I didn't make any mistakes! Let us know what you think.

[behdad]

Thanks James. Note that the final getArea() is not "area under the curve anymore". That would have been the initial integral, before subtracting the area under the line segment p0-p3. After the subtraction this because the signed area in the closed shaped formed by the curve and the line segment p3-p0. We can call this the area of the curve. Ie, if the curve is a degenerate line segment, the area is 0.

[Pomax]

yeah, "area under a curve" for parametric curves is kind of hard due to the absence of a strict "above" or "below", but getting the area between the curve and its start/end line segment is still pretty useful! Nice stuff, it's definitely something worth adding to the article!

You might be able to simplify the formulae by axis aligning the curve first, so that some of the terms fall away (x₀, y₀ already fall away after a translation, y₃ would fall away after a rotation but then requires introducing potential imprecisions, so might not be worth it)

[jamesgk]

In fact, this only calculates the "area of a curve" as I would intuitively understand it so long as the curve has no inflection points. Area to the left of the curve has different sign as area to the right, so they cancel and e.g. an S-shaped curve will have an "area" of zero. But maybe "area of a curve" is still an accurate-enough name.

We discussed the idea of aligning the curve before calculating the integral, but I don't think it actually simplifies the formulae per se considering the additional calculations needed to transform the points in the curve. Though I suppose that's pretty subjective.

In any case, would you like us to create a pull request with this code? I'd be happy to but I thought you might want to write the actual section in the primer yourself? But I'd be happy to do that as well.

[behdad]

Write. It's still the "signed area".

[Pomax]

as a lossless operation the translation by -(x,y)₀ seems a step that generally makes sense (if only to homogenise with most other papers I've read that try to use canonical curve forms). But yes, if you want to do a PR with the code I'd be more than happy to merge it in!

I'm also working on a newer pure JS rather than Processing bezier library that I want to swithc over to probably over the holidays too, so if you'd rather file it against that over on https://github.com/Pomax/bezierjs then that'd be great, too. I'm happy to port if you don't, though!

If you do, though, the current active branch is the browserify one, because I started writing bezier.js as a single file and ran into the "what am I doing" wall, it's been split up with a browserify build step instead, but I've not merged that to master yet.

[jamesgk]

Even more interesting -- doing this calculation with a quadratic curve and translating by -(x0,y0) yields:

(x2 * y1 - x1 * y2) / 3

or, one third the cross product of (x1,y1) and (x2,y2), i.e. one third the area of the parallelogram with these vectors as sides, i.e. two thirds the area of the triangle made by the control points! Which is a nice relationship. Behdad and I are thinking that there may be a similarly nice way to relate the area of a cubic curve to the polygons made by its control points, but we haven't looked into it too deeply.

I'm going to upload the code for the quadratic area, if you care to pull it into bezierjs. I might also make a pull request for bezierinfo, though I was having a little trouble visualizing the signed area with the Processing library.

[Pomax]

bezierjs is pure javascript though, not Processing =)

[jamesgk]

Right, the trouble I was having was with bezierinfo. I can't remember what the problem was off the top of my head, but perhaps it was related to the issue you recently posted. I'll come back to it eventually.

[Pomax]

cool. I'm still trying to figure out a good way to migrate bezierinfo to bezierjs instead - with the recent developments in React, it feels like I can significantly improve things by exploiting React as "section elements" including per-section native code.

[deanm]

I just stumbled across this discussion, and just in case anyone is interested I have some old notes on the topic:

http://deanm.github.io/graphics/2016/03/30/CurveArea.html

It perhaps provides some more references, the cubic equations were worked out long ago (in a metafont discussion from 2000), and there is a good writeup of it by Jackowski (https://tug.org/TUGboat/tb33-1/tb103jackowski.pdf).

The main point of the above article was trying to tackle the equation for rational quadratic Béziers, I haven't been able to find anywhere that has been done before. It's a bit tricky but should be usable.

[behdad]

@deanm Thanks for the writeup. Saved me writing that in the future!

I ended up using sympy to implement Green's theorem and apply it to find mean, variance, and covariance as well: https://github.com/behdad/fonttools/blob/5cd0a55635062f382ff86ccea6e42bbf9b8b030a/Snippets/symfont.py

[jamiemlaw]

Behdad and I are thinking that there may be a similarly nice way to relate the area of a cubic curve to the polygons made by its control points, but we haven't looked into it too deeply.

Well, it might be ten years later, but I was doing some research today and I think I found the formula you were looking for. For a cubic Bézier defined by control points p0, p1, p2 and p3, you take the area of polygon p0 p1 p2 p3 + area of polygon p0 p1 p3 + area of polygon p0 p2 p3, and multiply everything by 0.3.

Alternatively, you take twice the area of polygon p0 p1 p2 p3, add the area of polygon p0 p2 p1 p3, and multiply that by 0.3.

[behdad]

Thanks. How did you arrive at that?

[jamiemlaw]

In terms of how I got there, I admit it was mainly just trial and error. I was convinced there was a natural way to express a cubic Bézier's area in terms of its control points, and I was pushing numbers around until the problem relented and the symmetry popped out.

However, here's how I can demonstrate the result:

Via the shoelace formula, the area of polygon p0 p1 p2 p3 is

$$ x_0(y_3 - y_1) + x_1(y_0 - y_2) + x_2(y_1 - y_3) + x_3(y_2 - y_0) $$

or

$$ \begin{pmatrix} &x_0(&&-y_1&&+y_3&) \ +&x_1(&+y_0&&-y_2&&) \ +&x_2(&&+y_1&&-y_3&) \ +&x_3(&-y_0&&+y_2&&) \end{pmatrix} $$

polygon p0 p1 p3 is

$$ x_0(y_3 - y_1) + x_1(y_0 - y_3) + x_3(y_1 - y_0) $$

or

$$ \begin{pmatrix} &x_0(&&-y_1&&+y_3&) \ +&x_1(&+y_0&&&-y_3&) \ +&x_2(&&&&&) \ +&x_3(&-y_0&+y_1&&&) \end{pmatrix} $$

and polygon p0 p2 p3 is

$$ x_0(y_3 - y_2) + x_2(y_0 - y_3) + x_3(y_2 - y_0) $$

or

$$ \begin{pmatrix} &x_0(&&&-y_2&+y_3&) \ +&x_1(&&&&&) \ +&x_2(&+y_0&&&-y_3&) \ +&x_3(&-y_0&&+y_2&&) \end{pmatrix} $$

Summing these three areas together gives the familiar

$$ \begin{pmatrix} &x_0(&&-2y_1&-y_2&+3y_3&) \ +&x_1(&2y_0&&-y_2&-y_3&) \ +&x_2(&y_0&y_1&&-2y_3&) \ +&x_3(&-3y_0&y_1&+2y_2&&) \ \end{pmatrix} $$

[behdad]

Beautiful.

[jamesgk]

Nice! I like the alternative formulation as well; it's pretty funny how that works out.

[rice7th]

@jamiemlaw And for quadratic beziers how would that work?

[jamiemlaw]

@rice7th For quadratics, there’s a different formula. @jamesgk covered it earlier in the thread. It’s just two thirds the area of the triangle formed by the three control points.

[rice7th]

@jamiemlaw Oh I may have missed it then! Thank you so much!