mis-statement on Gauss Quadrature

[hydrodog]

It would be better to define Gauss Quadrature in a separate article frankly, and link to it. But given that you are trying to given an abridged version here, the statement that Guass Quadrature works better "because you have more points in the steep part of the function" or something like that is not true.

Trapezoidal uses two samples, one on each end of the interval. That's two degrees of freedom, so trapezoidal can perfectly represent a polynomial of zeroth order (constant) or 1st order (linear):

A = w0f(0) + w1f(1)

A Gauss 2nd order method has 4 degrees of freedom -- where to take the samples, and the weights.
That is: w0, w1, x0, x1

The convention is to think of Gauss Quadrature on the interval [-1,1].
The 2nd order method can therefore be perfectly accurate for any polynomial up to cubic:

P3(x) = a + bx + cx^2 + dx^3

Since cubics can represent a wide range of shapes (as you show so ably in your article) they can accurately integrate a wide variety of functions, anything that looks like a polynomial.

The method of proving this is really elegant. First approximate f(x) = 1
You know the area A=2 because it's on the interval [-1,1]
You could argue by symmetry, not knowing any better w0 = w1 which is true, but you don't have to:
w0f(x0) + w1f(x1) = 2. You know for this function that f(x0)=f(x1) = 1, so you have:

w0 + w1 = 2

Next, consider f(x) = x. f(-1) = -1, and f(1) = 1
A = 0 because the function is symmetric.

w0* f(x0) + w1*f(x1) = 0

So:
w0*f(x0) = -w1f(x1)

And on it goes, you can fit all 4 parameters. You will get an error approximating any polynomial higher than 3, or in your case, you are taking square root.

For asymptotes, polynomials are not a great approximation, so gauss quadrature will fail like anything else. You will have to use adaptive quadrature, try taking more samples near the pole, or do a transformation to avoid the problem in the first place. But that's not your problem here.

The error term for Gauss 2nd order is O(h^4).

A 3rd order Gaussian has error O(h^6).

For this problem, I was suggesting Romberg which is even better than Gauss Quadrature. Gauss is good when you have a pole on the end of the interval, because you don't have to sample at the end. But trapezoidal is much simpler, so watch what Romberg can do:

Romberg is black magic. We know that trapezoidal has error O(h^2). The Euler-Maclaurin theorem states that there are only even error terms. I am no mathematician, so I take that on faith (but it works).

So if you numerically integrate with trapezoidal:

S = h(1/2(f0 + fn) + (f1+f2 + ... + fn-1)) + O(h^2) + O(h^4) + O(h^6) ...

Exactly split each interval in half. That means you have twice as many slices, and the new h = h/2
your new s does not have to re-sum everything:

S1 = (1/2(f0 + fn) + (f1+f2 + ... + fn-1)) + O(h^2) + O(h^4) + O(h^6) ...
notice, I leave off the h (you can multiply by h later)

S2 = S1 + (f1.5 + f2.5 + f3.5 + ...fn-.5) + O(h2^2) + O(h2^4) + ...
So you just take your old sum, add in all the new points half way between.
You do not know what the errors are, they are O(h^2) but the new one has h/2 so:

R1 = (4S2 - S1) / 3
cancels your h^2 error term, leaving only the h^4 error term.
If you need to do it again, double the number of samples again (I will call it S4)
R2 = (4S4-S2)/3

Now with two Romberg estimates, you can cancel the next term:

Q1 = (16R2 - R1) / 15

Q1 now has only O(h^6) error. That's pretty amazing considering you could do 2, 4, and 8 samples and have final error O(that final h^6)
You can't go much beyond this, as little roundoff errors will mean that the coefficients between the two are not exactly the same, but the above achieves the equivalent of gauss quadrature 3rd order with far less computation. It's a pretty cool trick!

[Pomax]

I'm not sure I understand what prompted this issue... while I fully agree that the Legendre-Gauss quadrature procedure deserves its own article (written by not-me, ideally! =D), your initial comment suggests I tried to explain the actual reasons of why the LG quadrature works, which I in fact intentionally don't explain at all. I do link to a video lecture series on it, which is quite accessible, and I explain why we're picking it for our arc length computation in terms programmers can accept (it's efficient code) but I don't explain the maths behind it or why that math would even work (better) compared to other possible approaches.

So if I can pull this back a little to the actual article itself, this is the text associated with raising the LG quadrature as a thing we'll be using:

So we turn to numerical approaches again. The method we'll look at here is the Gauss quadrature. This approximation is a really neat trick, because for any nth degree polynomial it finds approximated values for an integral really efficiently. Explaining this procedure in length is way beyond the scope of this page, so if you're interested in finding out why it works, I can recommend the University of South Florida video lecture on the procedure, linked in this very paragraph.

This paragraph does not make any claims on how it works, it only explains that this is the numerical solution we're going to be using, because it'll be efficient.

The article also contains the following text in two later paragraphs:

So, the trick is to come up with useful rectangular strips. A naive way is to simply create n strips, all with the same width, but there is a far better way using special values for C and f(t) depending on the value of n, which indicates how many strips we'll use, and it's called the Legendre-Gauss quadrature.

and

This approach uses strips that are not spaced evenly, but instead spaces them in a special way that works remarkably well. If you look at the earlier sinoid graphic, you could imagine that we could probably get a result similar to the one with 99 strips if we used fewer strips, but spaced them so that the steeper the curve is, the thinner we make the strip, and conversely, the flatter the curve is (especially near the tops of the function), the wider we make the strip. That's akin to how the Legendre values work.

As far as I can tell, neither of these paragraphs try to explain the procedure, merely stating the properties of numeric integral computation in general, and what properties turn that into specifically the Legendre-Gauss way of doing so.

Can you explain where in the text you feel it is trying to make claims it works "better" than some other approach? You give the anecdote that the article claims Guass Quadrature works better "because you have more points in the steep part of the function" or something like that but I can't seem to find any text that would imply this, so it would be most helpful if you can highlight the part(s) of the article that you feel make this claim, be it explicit or implied.

[hydrodog]

This approximation is a really neat trick, because for any nth degree polynomial it finds approximated values for an integral really efficiently.

No, Gauss Quadrature of degree n finds an exact value of a polynomial of any lower degree. Your function is not a polynomial, which is why it is not exact.

so that the steeper the curve is, the thinner we make the strip, and conversely, the flatter the curve is (especially near the tops of the function), the wider we make the strip. That's akin to how the Legendre values work.

I don't think this is at all how the Legendre values work. If you try to integrate a function with a pole, you need more samples around the pole, but for the 45 degree angle sinusoid (is it really sinoid in British English, or is that a mistake?) this statement makes no sense. I tried to show you that the reason the fit for the Gauss method is at the roots of the Legendre polynomial has to do with the fact that you are fitting a perfect solution to polynomials up to degree n, allowing the position of the sample to vary, and that is simply the answer that pops out -- it has nothing at all to do with the slope of the particular function you are using. consider integrating: ``` f(x) = x^2 + 2x + 1 ``` from 0 to 2. The analytical answer is ``` x^3/3 + x^2 + x | = 8/3 + 4 + 2 = 8 1/3 ``` With Gauss 2nd order, you pick the midpoint (x=1), h=1, and the two sample points are: ``` x1 = mid - sqrt(1.0/3), x2 = mid + sqrt(1.0/3) ``` The answer `A = (f(x1) + f(x2))*h` should be accurate to 15 places, no iterations. If you do the same equation where the slope is quite high, say [10,13] the answer is still perfect. But the function: `sqrt((x^2+2x+1)^2 + (3x^2 -4x + 5)^2)` is first of all a polynomial of degree 4, which you are then taking theroot of. I would guess it's going to be very close to polynomial behavior of degree 2, but I'm not a good enough mathematician to know, so I have to try this out. Anyway, I think I pointed out a better way to go, much simpler. I don't have great test code yet but this quick hack indicates the approach. t varies from 0 to 1 (Here I called i X, sorry). For h=t/8 = .125, Romberg should have error .125^6 = 2^-18. That's nearly 6 digit accuracy! If that isn't accurate enough on your interval, then doing two halves ought to do it. Remember that while the function may not be analytically solvable, it's relatively close in behavior to a 2nd order polynomial. I am an efficiency nut so I precomputed the coefficients and then just pass them into a function to calculate the "f'(t)^2" The numerical integration part uses exactly 8 function invokations, and you don't need those odd coefficients and weights from Gauss Quadrature ``` fprime(e, f, g, h, i, x) { // p'(x) = 3ax^2 + 2bx + c // p'^2(x) = 3acx^2 + 2bcx + c^2 + 6abx^3 + 4b^2x^2 + 2bcx + 9a^2x^4 + 6abx^3 + 3acx^2 // = 9a^2x^4 + 12abx^3+(6ac+4b^2)x^2+4bcx+c^2 return (((e*x+f)*x+g)*x+h)*x+i; } // for cubic spline segments where parameter t ranges from 0 to 1 public double arclen() { e = 9*a*a, f = 12*a*b, g = 6*a*c+4*b, h = 4*b*c, i=c*c; S0 = (fprime(e,f,g,h,i,0) + fprime(e,f,g,h,i,1)) *.5; S1 = S0 + fprime(e,f,g,h,i,0.5); S2 = S1 + fprime(e,f,g,h,i,0.25) + fprime(e,f,g,h,i,0.75); S4 = S2 + fprime(e,f,g,h,i,0.125) + fprime(e,f,g,h,i,0.375) + fprime(e,f,g,h,i,.625) + fprime(e,f,g,h,i,.875); R1 = (4*S2 - S1) / 3; R2 = (4*S4 - S2) / 3; Q1 = (16*R2 - R1) / 15; return Q1 * 0.125; } The above code is not yet correct. I just did arclen for a cubic polynomial squared when what you need is arclen for: ``` f(x) = sqrt( x'(t)^2 + y'(t)^2 + z'(t)^2) ``` If I use gauss 3rd order on each half of the interval, I need 6 function invocations with an error of .25^6 which is not as good. Comparable code for Gauss 5th order uses only 5 function invokations with an error estimate of 9.7 x 10 ^-4 ``` // for cubic spline segments where parameter t ranges from 0 to 1 // uses Gauss Quadrature, order 5 mid = 0.5; c1 = Math.sqrt(5 - 2 * Math.sqrt(10.0/7)) / 3; c2 = Math.sqrt(5 + 2 * Math.sqrt(10.0/7)) / 3; w0 =128.0/225; w1 = (322 + 13*Math.sqrt(70)) / 900; w2 = (322 - 13*Math.sqrt(70)) / 900; x0 = 0.5; x1 = x0 - c1; x2 = x0 + c1; x3 = x0 - c2; x4 = x0 + c2; arclen2() { e = 9*a*a; f = 12*a*b; g = 6*a*c+4*b; h = 4*b*c; i=c*c; h2 = 0.5; return (w0*fprime(e,f,g,h,i,x0) + w1*(fprime(e,f,g,h,i,x1) + fprime(e,f,g,h,i, x2)) + w2*(fprime(e,f,g,h,i,x3) + fprime(e,f,g,h,i,x4))) *h2; } ``` I guess the two are roughly comparable for this problem. You can always go up to a higher order method in Gauss Quadrature, but you cannot easily double the number of points. I therefore usually prefer trapezoidal because if we need better accuracy, just interleave the points and calculate the next Romberg value which will be 64 times more accurate. With Gauss it's not trivial to improve the answer without throwing out the previous function invokations. With trapezoid it's much easier.

[behdad]

FWIW, I also was very confused when I read the Guass Quadrature section and the explanation for why it's superior. Took me quite a while of reading Wikipedia to understand what it's doing and appreciate its beauty. The way I appreciate it is this: a Guass Quadrature of sampling n points, gives the exact integral solution for any polynomial of up to degree 2n-1. Since the sampling points are chosen such that all polynomials of degree up to 2n-1 passing through those points have the same integral, this relieves one of first having to approximate the function-to-be-integrated by a polynomial and then integrating it. So, by just sampling 5 points, one gets the exact area of the best-fitting 9th degree polynomial approximation going through those 5 points. The fact that there is a rigorous explanation to the solution is the beauty IMO.

Anyway, for cubic and quadratic curve length I've been using these https://github.com/fonttools/fonttools/blob/master/Lib/fontTools/pens/perimeterPen.py :

	def _addQuadraticQuadrature(self, c0, c1, c2):
		# Approximate length of quadratic Bezier curve using Gauss-Legendre quadrature
		# with n=3 points.
		#
		# This, essentially, approximates the length-of-derivative function
		# to be integrated with the best-matching fifth-degree polynomial
		# approximation of it.
		#
		#https://en.wikipedia.org/wiki/Gaussian_quadrature#Gauss.E2.80.93Legendre_quadrature

		# abs(BezierCurveC[2].diff(t).subs({t:T})) for T in sorted(.5, .5±sqrt(3/5)/2),
		# weighted 5/18, 8/18, 5/18 respectively.
		v0 = abs(-0.492943519233745*c0 + 0.430331482911935*c1 + 0.0626120363218102*c2)
		v1 = abs(c2-c0)*0.4444444444444444
		v2 = abs(-0.0626120363218102*c0 - 0.430331482911935*c1 + 0.492943519233745*c2)

		self.value += v0 + v1 + v2

	def _addCubicQuadrature(self, c0, c1, c2, c3):
		# Approximate length of cubic Bezier curve using Gauss-Lobatto quadrature
		# with n=5 points.
		#
		# This, essentially, approximates the length-of-derivative function
		# to be integrated with the best-matching seventh-degree polynomial
		# approximation of it.
		#
		# https://en.wikipedia.org/wiki/Gaussian_quadrature#Gauss.E2.80.93Lobatto_rules

		# abs(BezierCurveC[3].diff(t).subs({t:T})) for T in sorted(0, .5±(3/7)**.5/2, .5, 1),
		# weighted 1/20, 49/180, 32/90, 49/180, 1/20 respectively.
		v0 = abs(c1-c0)*.15
		v1 = abs(-0.558983582205757*c0 + 0.325650248872424*c1 + 0.208983582205757*c2 + 0.024349751127576*c3)
		v2 = abs(c3-c0+c2-c1)*0.26666666666666666
		v3 = abs(-0.024349751127576*c0 - 0.208983582205757*c1 - 0.325650248872424*c2 + 0.558983582205757*c3)
		v4 = abs(c3-c2)*.15

		self.value += v0 + v1 + v2 + v3 + v4

(I'm using complex type to represent points.)

[Pomax]

Do I not link to the most excellent video lecture for it in the bezier article, only on the dedicated lookup-values page? If so, I should change that. Those lectures are short and sweet, and can explain things infinitely better than I can do in text on the Bezier page itself, which is not the place to explain this in any level of detail.

[behdad]

Do I not link to the most excellent video lecture for it in the bezier article, only on the dedicated lookup-values page? If so, I should change that. Those lectures are short and sweet, and can explain things infinitely better than I can do in text on the Bezier page itself, which is not the place to explain this in any level of detail.

I also got confused by the original poster's point, which is, this is wrong:

"This approach uses strips that are not spaced evenly, but instead spaces them in a special way that works remarkably well. If you look at the earlier sinoid graphic, you could imagine that we could probably get a result similar to the one with 99 strips if we used fewer strips, but spaced them so that the steeper the curve is, the thinner we make the strip, and conversely, the flatter the curve is (especially near the tops of the function), the wider we make the strip. That's akin to how the Legendre values work."

[Pomax]

ah. yeah that's fair, and worth changing the phrasing of to remove any suggestion that "what would be more useful" is literally what LGQ does.

[Pomax]

updated the article text to no longer talk about LGQ being in any way related to the steepness of the curve.

[raphlinus]

Sorry to add to a closed issue (let me know if there's a better place to continue this discussion). LGQ is indeed amazing, and I have evidence. The plots below are the log10 of the error of measuring a quadratic Bezier with endpoints (-1, 0), (x, y), (1, 0), where x, y are the parameters of the plot. The first plot is L0 from "Adaptive subdivision and the length and energy of Bézier curves" by Jens Gravesen. To recap briefly, this is 1/3 the sum of the lengths of the endpoints to the control point + 2/3 the length of the chord. The second plot is the 3-point LGQ as posted by @behdad above.

You can see some of the magic here. The cost to compute the two is about the same, they're both three "hypot" operations plus a handful of adds and multiplies.

The only problem I'm having is: the Gravesen paper gives an error bound, which is very useful for guiding adaptive subdivision. I'm sure there's a good way to compute an error bound for LGQ, but it's not immediately obvious to me. I figure maybe some math expert can help?

[behdad]

Right. I have not been able to derive the error either, and am interested in as well.

[behdad]

I'm also curious to see measurements for cubics.

[behdad]

You can see some of the magic here.

Can you explain the graphs please?

[raphlinus]

The x and y coordinates of the plot represent the position of the control point of the quadratic Béziers, where the endpoints are (-1, 0) and (1, 0). This captures the space of all possible quadratic Béziers, modulo scaling, rotation, and translation, all of which both metrics are invariant to. The color is the number of decimal digits of precision.

[waruyama]

Sorry to add something to a long closed issue, but I actually implemented Romberg's method to approximate the arc length of a cubic Bézier curve.

TL;DR
Gauss Quadrature is better than Romberg's method

For reference, here is my implementation, which is basically the code from Wikipedia converted to Javascript:

function getLength(v, t0, t1) {   
  // Determine arc length using Romberg's method
  // see https://en.wikipedia.org/wiki/Romberg%27s_method

  let x0 = v[0];
  let y0 = v[1];
  let x1 = v[2];
  let y1 = v[3];
  let x2 = v[4];
  let y2 = v[5];
  let x3 = v[6];
  let y3 = v[7];

  // Get coefficients for first derivative
  let ax = 3 * (x3 - x0) - 9 * (x2 - x1);
  let ay = 3 * (y3 - y0) - 9 * (y2 - y1);
  let bx = 6 * (x2 - 2 * x1 + x0);
  let by = 6 * (y2 - 2 * y1 + y0);
  let cx = 3 * (x1 - x0);
  let cy = 3 * (y1 - y0);
  // Make derivative function
  let f = t => Math.hypot((ax * t + bx) * t + cx, (ay * t + by) * t + cy);
  
  let maxSteps = 5;
  let epsilon = 1e-12;
  let h = t1 - t0;

  let Rp = new Array(maxSteps);
  let Rc = new Array(maxSteps);
  Rp[0] = (f(t0) + f(t1)) * h * 0.5;

  for (let i = 1; i < maxSteps; i++) {
    h /= 2;
    let c = 0;
    let ep = 1 << (i - 1);
    for (let j = 1; j <= ep; j++) {
      c += f(t0 + (2 * j - 1) * h);
    }
    Rc[0] = h*c + 0.5*Rp[0]; //R(i,0)
    for (let j = 1; j <= i; j++) {
      let n_k = Math.pow(4, j);
      Rc[j] = (n_k*Rc[j-1] - Rp[j-1])/(n_k-1); // compute R(i,j)
    }
    if (i > 1 && Math.abs(Rp[i - 1] - Rc[i]) < epsilon) {
      return Rc[i - 1];
    }
    let tmp = Rp;
    Rp = Rc;
    Rc = tmp;
  }
  return Rp[maxSteps-1];
}

My function uses maxSteps=5, which results in 16 calls to the derivative function. This is the same number of calls as Paper.js performs in Curve.getLength() using Gauss Quadrature, so the results can be directly compared (Bezierjs uses 24 calls, so it is difficult to compare to Romberg's method). The precision of Gauss Quadrature in Paper.js is in most cases better, sometimes by order of magnitude, than my implementation using Romberg's method.

What I like about Romberg's method is that we can specify a tolerance and stop calculation if the tolerance is reached. Also, there is no need to store the weird precomputed weights, but this is purely a code aesthetic issue.

[devshgraphicsprogramming]

What I like about Romberg's method is that we can specify a tolerance and stop calculation if the tolerance is reached. Also, there is no need to store the weird precomputed weights, but this is purely a code aesthetic issue.

Memory > ALU, you might find that if you do a lot of coefficients/weights the time to fetch them vs compute them might turn against the fetching approach.