+float Squared3xCurveArea(const float *a, const float *control, const float *b, int components)
+{
+#if 0
+ // mimicing the old behaviour with the new code...
+
+ float deviation;
+ float quartercurvearea = 0;
+ int c;
+ for (c = 0;c < components;c++)
+ {
+ deviation = control[c] * 0.5f - a[c] * 0.25f - b[c] * 0.25f;
+ quartercurvearea += deviation*deviation;
+ }
+
+ // But as the new code now works on the squared 2x curve area, let's scale the value
+ return quartercurvearea * quartercurvearea * 64.0;
+
+#else
+ // ideally, we'd like the area between the spline a->control->b and the line a->b.
+ // but as this is hard to calculate, let's calculate an upper bound of it:
+ // the area of the triangle a->control->b->a.
+ //
+ // one can prove that the area of a quadratic spline = 2/3 * the area of
+ // the triangle of its control points!
+ // to do it, first prove it for the spline through (0,0), (1,1), (2,0)
+ // (which is a parabola) and then note that moving the control point
+ // left/right is just shearing and keeps the area of both the spline and
+ // the triangle invariant.
+ //
+ // why are we going for the spline area anyway?
+ // we know that:
+ //
+ // the area between the spline and the line a->b is a measure of the
+ // error of approximation of the spline by the line.
+ //
+ // also, on circle-like or parabola-like curves, you easily get that the
+ // double amount of line approximation segments reduces the error to its quarter
+ // (also, easy to prove for splines by doing it for one specific one, and using
+ // affine transforms to get all other splines)
+ //
+ // so...
+ //
+ // let's calculate the area! but we have to avoid the cross product, as
+ // components is not necessarily 3
+ //
+ // the area of a triangle spanned by vectors a and b is
+ //
+ // 0.5 * |a| |b| sin gamma
+ //
+ // now, cos gamma is
+ //
+ // a.b / (|a| |b|)
+ //
+ // so the area is
+ //
+ // 0.5 * sqrt(|a|^2 |b|^2 - (a.b)^2)
+ int c;
+ float aa = 0, bb = 0, ab = 0;
+ for (c = 0;c < components;c++)
+ {
+ float xa = a[c] - control[c];
+ float xb = b[c] - control[c];
+ aa += xa * xa;
+ ab += xa * xb;
+ bb += xb * xb;
+ }
+ // area is 0.5 * sqrt(aa*bb - ab*ab)
+ // 2x TRIANGLE area is sqrt(aa*bb - ab*ab)
+ // 3x CURVE area is sqrt(aa*bb - ab*ab)
+ return aa * bb - ab * ab;
+#endif
+}
+