Obviously before this is eventually accepted the printfs will have to go -- If the info is highly helpful diagnostically you could turn them into TRACE()es, and if less helpful they could become commented-out TRACE()es (of which there are many in the code).

How does the performance compare to the "old existing" algorithm? If it's good are you keeping is_accurate as a fig leaf and how do you see the risks? This road began because of its shortcomings.

I'm a bit surprised at all the trig. Is that just how the calculation you're working from was written or did you consider and reject replacing a bunch of that with unit vectors?

I have thought that I have already removed the printf's - but obviously in a different version... After the last commit they should not appear anymore. I think there is no need for TRACE - it seems to be stable for merging/simplifying (but I just have tested it with 20 different fonts).

The performance compared to "my" brute force algorithm is much better, I would guess about 1/100 as much loops as in the brute force for the generic case (but I think this doesn't make it exactly 100 times faster - but at least 10 times faster). The performance seems to be comparable to George Williams algorithm. That's why I used it for "simplify" as well. On my computer, the simplify action for a 1600 glyph font still feels instantaneous.

Indeed, I keep is_accurate as a fig leaf: The algorithm does not work out of the box with "stroke path" (causes a segmentation fault that I could not trace). This issue may be connected with a behaviour that may be unpleasant for some type designers: The handles may exceed the Tunni-point (the point where the handles would intersect) after applying the algorithm. This unpleasant behaviour may even (rarely) occur when simplifying. Of course, I could limit the algorithm to a convex control polygon but in my opinion in some situations it may be desirable.

For me, this patch is better than the brute force algorithm and George William's algorithm. But I think somebody of the FontForge team should test this patch and judge himself/herself if it's better for a majority of users.

By the way: Please excuse the introduction of the snake_case for is_accurate. I think, I should convert it to camelCase. Does FontForge have an official code style?

About trigonometry: I don't think that things will speed up if I replace ca=cos(alpha), sa=sin(alpha) by a unit vector (ca,sa). (?)

Indeed, I keep is_accurate as a fig leaf: The algorithm does not work out of the box with "stroke path" (causes a segmentation fault that I could not trace).

Can you give me a reproducible test case for this?

By the way: Please excuse the introduction of the snake_case for is_accurate. I think, I should convert it to camelCase. Does FontForge have an official code style?

We don't have an "official" code style but snake-case for variables is preferred, I think.

About trigonometry: I don't think that things will speed up if I replace ca=cos(alpha), sa=sin(alpha) by a unit vector (ca,sa). (?)

The question was more about why the code is doing anything with angles in the first place. I might be biased but I think of trig as a sort of cheat or work-around for when vectors aren't handy and it has a reputation for being a bit "expensive" computationally.

Try this out for comparison:

diff --git a/fontforge/splinefit.c b/fontforge/splinefit.c
index 5103c6f87..904035335 100644
--- a/fontforge/splinefit.c
+++ b/fontforge/splinefit.c
@@ -33,6 +33,7 @@
 #include "splineorder2.h"
 #include "splineutil.h"
 #include "splineutil2.h"
+#include "utanvec.h"
 
 #include <math.h>
 
@@ -813,17 +814,21 @@ return( SplineMake3(from,to));
     if (is_accurate) { 
        bigreal a,b,alpha,beta,rotation,ca,sa,sasa,cb,sb,sab,f,m,xa,ya,xb,yb,xc,yc,xd,yd;
        int numberOfSolutions;
+    BasePoint alpha_uv = (BasePoint) { BPDot(ftunit, fromunit), BPCross(ftunit, fromunit) };
     alpha = acos(ftunit.x*fromunit.x+ftunit.y*fromunit.y); /* signed angle */
     if ( ftunit.x*fromunit.y-ftunit.y*fromunit.x < 0 ){
                alpha = -alpha; 
        }
+    BasePoint beta_uv = (BasePoint) { BPDot(BPRev(ftunit), tounit), BPCross(ftunit, tounit) };
     beta = acos(-ftunit.x*tounit.x-ftunit.y*tounit.y); /* signed angle */
     if ( ftunit.x*tounit.y-ftunit.y*tounit.x < 0 ){
                beta = -beta;   
        }
+    BasePoint rotation_uv = BPNeg(ftunit);
     rotation = atan2(ftunit.y,ftunit.x);
     ca = cos(-rotation);
        sa = sin(-rotation);
+    printf("rotation_uv: %lf,%lf    ca: %lf    sa: %lf\n", rotation_uv.x, rotation_uv.y, ca, sa);
     SplinePoint *frompoint,*topoint;
     f = 0; /* area */
     m = 0; /* first area moment about y (along x) */
@@ -853,12 +858,18 @@ return( SplineMake3(from,to));
                m = -m;
                f = -f;
        }
+    if ( alpha_uv.y < 0 ) {
+       alpha_uv = BPNeg(alpha_uv);
+       beta_uv = BPNeg(beta_uv);
+    }
        /* start approximation by solving the quartic equation */
        sa = sin(alpha); /* overwriting previous sa */
        ca = cos(alpha); /* overwriting previous ca */
+       printf("alpha_uv: %lf,%lf    ca: %lf    sa:%lf\n", alpha_uv.x, alpha_uv.y, ca, sa);
        sasa = sa*sa;
        sb = sin(beta);
        cb = cos(beta);
+       printf("beta_uv: %lf,%lf    cb: %lf    sb:%lf\n", beta_uv.x, beta_uv.y, cb, sb);
        sab = sin(alpha+beta);
        Polynomial aQuartic;
        aQuartic.length = 5;

(I'm not going to stamp my foot and require a no-trig solution but

    rotation = atan2(ftunit.y,ftunit.x);
    ca = cos(-rotation);
    sa = sin(-rotation);

is a roundabout way of negating a value ... )

This issue may be connected with a behaviour that may be unpleasant for some type designers: The handles may exceed the Tunni-point (the point where the handles would intersect) after applying the algorithm. This unpleasant behaviour may even (rarely) occur when simplifying. Of course, I could limit the algorithm to a convex control polygon but in my opinion in some situations it may be desirable.

The concern about this may not be just a matter of taste. Apparently this exposes a bug in the stroking algorithm, which we should fix (or at least prevent the crash). That demonstrates a general worry: adding such splines to contour geometries could trigger bugs in other downstream spline-processing algorithms.

I don't think we should necessarily prohibit them but at a minimum I'd be inclined to avoid them in Element -> Simplify -> Simplify with a not-checked-by-default checkbox in "Simplify More" and an analogous flag in Contour.Simplify() and its siblings.

Well, the part

rotation = atan2(ftunit.y,ftunit.x);
ca = cos(-rotation);
sa = sin(-rotation);

was unnecessary, indeed. The rotation can be done directly using ftunit. I was not aware of utanvec.h but it makes things much cleaner. I saw there is even BPRot() and BPScale() that would make the notation of the following block easier but not faster nor more compact. I think I will leave it like that:

xa = ((frompoint->me.x-from->me.x)*ftunit.x+(frompoint->me.y-from->me.y)*ftunit.y)/ftlen;
ya = (-(frompoint->me.x-from->me.x)*ftunit.y+(frompoint->me.y-from->me.y)*ftunit.x)/ftlen;
xb = ((frompoint->nextcp.x-from->me.x)*ftunit.x+(frompoint->nextcp.y-from->me.y)*ftunit.y)/ftlen;
yb = (-(frompoint->nextcp.x-from->me.x)*ftunit.y+(frompoint->nextcp.y-from->me.y)*ftunit.x)/ftlen;
xc = ((topoint->prevcp.x-from->me.x)*ftunit.x+(topoint->prevcp.y-from->me.y)*ftunit.y)/ftlen;
yc = (-(topoint->prevcp.x-from->me.x)*ftunit.y+(topoint->prevcp.y-from->me.y)*ftunit.x)/ftlen;
xd = ((topoint->me.x-from->me.x)*ftunit.x+(topoint->me.y-from->me.y)*ftunit.y)/ftlen;
yd = (-(topoint->me.x-from->me.x)*ftunit.y+(topoint->me.y-from->me.y)*ftunit.x)/ftlen;

Thanks for the hint unitdirection = (BPDot(a,b),BPCross(a,b)); I was not aware that one could handle it that elegantly.
I have changed the notation alpha_uv to aunit to be nearer to the notation that George Williams used in the rest of the code (like fromunit, tounit, ftunit).
After the last commit no trigonometric function is used.

About tunni-point-excess: I will make a further commit for it.

The algorithm now does no longer produce handles that exceed the Tunni point (see last commit). It seems to work fine and it may be wise to keep this behaviour anyway. But I was not successful changing all flags is_accurate to true. Any font/glyph that I try to stroke-expand will cause a segmentation fault. May be I am missing something obvious?

I was not aware of utanvec.h but it makes things much cleaner. I saw there is even BPRot() and BPScale()

I added utanvec when I was reimplementing the stroke code. Because the BP stuff is inline it shouldn't have much of an effect on performance -- optimizers are aggressive and capable these days, especially when they have the syntax tree available (in contrast with macros). So optimizing for readability is probably the best heuristic and contributors can use it however they like.

Any font/glyph that I try to stroke-expand will cause a segmentation fault. May be I am missing something obvious?

I think this is fairly simple and needs fixing. If you look at (current) lines 744 through 758 you'll see that the tounit and fromunit calculations handle the cases where to->prev is NULL and from->next is NULL. Normally in calls from Simplify() and such you're potentially replacing a spline or splines but when the stroke code calls this function it's constructing the spline. So it just constructs the to point and sets its slope and lets the fit points do the rest of the work.

The for loop on (current) line 820 is assuming from and to are connected, and is therefore dereferencing from->next, which is NULL.

This seemed to do the trick with some basic testing:

diff --git a/fontforge/splinefit.c b/fontforge/splinefit.c
index fd1215fb1..42c79985d 100644
--- a/fontforge/splinefit.c
+++ b/fontforge/splinefit.c
@@ -817,7 +817,12 @@ return( SplineMake3(from,to));
     SplinePoint *frompoint,*topoint;
     f = 0; /* area */
     m = 0; /* first area moment about y (along x) */
-    for ( frompoint = from, topoint = from->next->to ; ; frompoint = topoint->next->from, topoint = topoint->next->to ) {
+    frompoint = from;
+    if ( from->next==NULL )
+       topoint=to;
+    else
+       topoint=from->next->to;
+    for ( ; ; frompoint = topoint->next->from, topoint = topoint->next->to ) {
                /* normalizing transformation (chord to x axis and length 1) */
                xa = ((frompoint->me.x-from->me.x)*ftunit.x+(frompoint->me.y-from->me.y)*ftunit.y)/ftlen;
                ya = (-(frompoint->me.x-from->me.x)*ftunit.y+(frompoint->me.y-from->me.y)*ftunit.x)/ftlen;
@@ -1273,7 +1278,7 @@ SplinePoint *_ApproximateSplineSetFromGen(SplinePoint *from, SplinePoint *to,
     to->prevcp.x = to->me.x - fp[cnt-1].ut.x;
     to->prevcp.y = to->me.y - fp[cnt-1].ut.y;
     to->noprevcp = false;
-    ApproximateSplineFromPointsSlopes(from,to,fp+1,cnt-2,order2,false);
+    ApproximateSplineFromPointsSlopes(from,to,fp+1,cnt-2,order2,true);
 
     for ( i=0; i<cnt; ++i ) {
        d = SplineMinDistanceToPoint(from->next, &fp[i].p);

what do you think?

I think it works! Thanks for the patch! Now that there are no segmentation faults anymore I have tested the speed: The new algorithm causes stroking to be about 8 times slower than before. I can probably fiddle the algorithm a bit more to be faster (the comparison of the L2-error is somewhat slow, I should take the parametrization time of the fitting points in account as well). This will need some time...
By the way: Shouldn't the accuracy target of the stroking dialog have the default 1 instead of 0.25? I was a bit confused that although I had checked "simplify" in the stroking dialog I could simplify it further after stroking.

At the time I felt like 1em accuracy could lead to some visible distortions. The thing about stroke expansion is that people are pretty good at judging consistent width vs changes in width, so I felt it was OK to be extra picky by default.

I have tested several mechanisms to choose the approximation with the lowest error (those by George Williams as well). But it does not speed up much. And now I suppose that the old algorithm is indeed the better choice for the current stroking algorithm because it calls ApproximateSplineFromPoints() too. If the stroking algorithm stays at it is, I would suggest to merge the current state.
However, it may be worth to change (anyway) a small thing of the stroking algorithm:

The inner path is the result of a stroking without simplifying. The outer path is a circular arc with some intermediate points to show that it would be probably better to shorten the handles by the factor 2/3 when stroking (instead of meeting at the Tunni point). The mathematical reason to do so is that the stroked bezier segments are nearly elliptic, which means nearly parabolic, which means handles = 2/3 of the length to the Tunni point.
If the stroking algorithm changes in this way it could be (I have not tested it), that the new approximation algorithm works better.

@linusromer There are basically four sources of points in the expand stroke system:

1. Points on the source contour (the algorithm doesn't try to fit "through" source points even when they are smooth -- it could but it doesn't.
2. Points on the nib: It doesn't fit "through" these either, although I've considered that as an enhancement, especially as when stroking with a circular nib (the most common case) you get points from the spline approximation of the circle.
3. Points generated by the fitting algorithm.
4. Points related to join geometry (only at joins).

In short, there's nothing in the algorithm that generates co-located control points like those in your sample.

I'm not sure how you generated that curve, and I'm happy to help you look into it, but my initial suspicion is that somewhere in your test code you've converted to a quadratic contour and back to a cubic contour -- that's the easiest way to wind up with contours of that form.

On the more general question I believe any calls from the stroking algorithm to ApproximateSplineFromPoints() are fallbacks in the fitting code. If those fallbacks are still needed it suggests there are "holes" in this new algorithm, and that's something to explore and understand a bit better before merging.

@skef You are right, I was wrong about the standard stroking algorithm. I was misled by the following effect (original S, stroked S using the standard ApproximateSplineFromPointsSlopes() without simplify, stroked S using the new ApproximateSplineFromPointsSlopes() without simplify):

I thought that "without simplify" means that ApproximateSplineFromPointsSlopes() is not used (which is obviously not the case). Please excuse this unnecessary suggestion. I think I have to understand your stroke algorithm deeper to see what goes wrong here.

The simplify option to the stroking algorithm post-processes the contours separately by calling the regular Simplify facility (with particular options). This is to remove potentially "unnecessary" points of types 1 and 2 in my list above.

(If the nib points or source contour points are smooth there's a good chance you won't need the corresponding points on the output contour. On the other hand if you're trying to make a variable font by stroking at different nib magnitudes you probably need all of those...)

@linusromer I think this is good news, by the way.

The place to look for this problem is in _ApproximateSplineSetFromGen(). If its dividing splines this much it doesn't think the splines you're outputting are fitting as well as ones made by the old algorithm. It basically checks all the points to see whether they're within the error tolerance and if any aren't it splits at the worst point and tries to fit splines on either side. So either something is wrong with your "answers" or something is wrong with my tolerance calculations (although it seems to work fine with the old algorithm).

This is good news because there are so many more splines being fit than before, so once the problem is diagnosed I expect the performance will be comparable to the old fitting algorithm.

On the error metrics:

It seems like there must be a maximum "slop" between source t values and candidate t values that could probably be determined experimentally if not analytically. If there is then the mid calculating loop could use 99/t and only check the subset of the 99 calculated points within the slop.

About the idea with the maximal "slop": The slop becomes maximal when one handle is 0 and the other is infinitely long. The maximal slop with "reasonable" handles seems to be about +/- 0.4 when comparing the time to the normed time (normed = parametrized by arc length and divided by the total length). Assuming the maximal slop +/- 0.5 seems to be safe.
Then there is another property that can speed up things: If we have found the nearest approximation point to a given fit point mid[i] at approx[j], we can assume for "reasonable" handles that mid[i] has its nearest counterpart at approx[>=j].
Combining the maximal "slop" and this property, one gets:

for (int i=0; i<cnt; i++) { /* Going through the mid points */
	error = (mid[i].p.x-approx[0].x)*(mid[i].p.x-approx[0].x)
	+(mid[i].p.y-approx[0].y)*(mid[i].p.y-approx[0].y);
	/* Above we have just initialized the error and */
	/* now we are going through the remaining 98 of */
	/* 99 points on the approximate cubic bezier: */
	minIndex = 0;
	for (int j=0; j<fmin(minIndex+50,99); j++) {
		dist = (mid[i].p.x-approx[j].x)*(mid[i].p.x-approx[j].x)
		+(mid[i].p.y-approx[j].y)*(mid[i].p.y-approx[j].y);
		if (dist < error) {
			error = dist;
			minIndex = j;
		}
	}
	errorsum += error;
	if (errorsum > bestError)
		break;
}

And guess what: I have measured the time and it practically does not matter. (Storing minIndex and calculating fmin(minIndex+50,99) seem to compensate the subset optimization.)

Anyway, I have compared the time with the error computation of _ApproximateSplineSetFromGen(), so I have used SplineMinDistanceToPoint() in the following way:

bigreal bestError = 1e30;
bigreal maxerr,d;
SplinePoint fromtemp,totemp;
fromtemp.me = from->me;
totemp.me = to->me;	
Spline *temp;
for (int k=0; k<numberOfSolutions; k++) {
	temp = SplineMake3(&fromtemp,&totemp);
	fromtemp.nextcp.x = from->me.x+ftlen*fromunit.x*abSolutions[k][0];
	fromtemp.nextcp.y = from->me.y+ftlen*fromunit.y*abSolutions[k][0];
	totemp.prevcp.x = to->me.x+ftlen*tounit.x*abSolutions[k][1];
	totemp.prevcp.y = to->me.y+ftlen*tounit.y*abSolutions[k][1];
	maxerr = 0; 
	for ( int i=0; i<cnt; ++i ) {
		d = SplineMinDistanceToPoint(temp,&mid[i].p);
		if ( d>maxerr ) 
			maxerr = d;
	}
	if ( maxerr < bestError ) {
		bestError = maxerr;
		from->nextcp = fromtemp.nextcp;
		to->prevcp = totemp.prevcp;
	}
}

And this was way slower! (About 40 times for maxerr, about 30 times for the summed up error).

Concluding from these (and other experiments) I still favour my current error computation. So I have to check what goes wrong elsewhere with stroking and my implementation.

Going back to the other problem does the version with SplineMinDistanceToPoint() also have the too-many splines problem when stroking? That would be interesting to know because that's the measure being used one level up.

If we eventually get this figured out we could expose an interface (probably with a passed struct or something) to retrieve the worst fitting point, and may be the distance, so that doesn't have to happen twice. That should make the overall system quite a bit faster.

BTW that's not quite the right spline calculation -- SplineMake() or SplineRefigure() need to be run after setting the control points to correctly calculate the coefficients. I just did a bit of testing with this:

diff --git a/fontforge/splinefit.c b/fontforge/splinefit.c
index 33c7d4d12..77160b450 100644
--- a/fontforge/splinefit.c
+++ b/fontforge/splinefit.c
@@ -31,6 +31,7 @@
 #include "fontforge.h"
 #include "splinefit.h"
 #include "splineorder2.h"
+#include "splinerefigure.h"
 #include "splineutil.h"
 #include "splineutil2.h"
 #include "utanvec.h"
@@ -956,51 +957,32 @@ return( SplineMake3(from,to));
                to->prevcp.y = to->me.y+ftlen*tounit.y*abSolutions[0][1];
        } else { /* compare L2 errors to choose the best solution */
                bigreal bestError = 1e30;
-               bigreal t,error,errorsum,dist;
-               BasePoint prevcp,nextcp,coeff1,coeff2,coeff3;
+               bigreal maxerr, d;
+               SplinePoint *fromtemp = SplinePointCreate(from->me.x, from->me.y);
+               SplinePoint *totemp = SplinePointCreate(to->me.x, to->me.y);
+               Spline *temp = SplineMake3(fromtemp, totemp);
                for (int k=0; k<numberOfSolutions; k++) {
-                       nextcp.x = from->me.x+ftlen*fromunit.x*abSolutions[k][0];
-                       nextcp.y = from->me.y+ftlen*fromunit.y*abSolutions[k][0];
-                       prevcp.x = to->me.x+ftlen*tounit.x*abSolutions[k][1];
-                       prevcp.y = to->me.y+ftlen*tounit.y*abSolutions[k][1];
-                       /* Calculate the error of the cubic bezier path from,nextcp,prevcp,to: */
-                       /* In order to do that we calculate 99 points on the bezier path. */
-                       coeff3.x = -from->me.x+3*nextcp.x-3*prevcp.x+to->me.x;
-                       coeff3.y = -from->me.y+3*nextcp.y-3*prevcp.y+to->me.y;
-                       coeff2.x = 3*from->me.x-6*nextcp.x+3*prevcp.x;
-                       coeff2.y = 3*from->me.y-6*nextcp.y+3*prevcp.y;
-                       coeff1.x = -3*from->me.x+3*nextcp.x;
-                       coeff1.y = -3*from->me.y+3*nextcp.y;
-                       BasePoint approx[99];
-                       for (int i=0; i<99; i++) {
-                               t = (i+1)/100.0;
-                               approx[i].x = from->me.x+t*(coeff1.x+t*(coeff2.x+t*coeff3.x));
-                               approx[i].y = from->me.y+t*(coeff1.y+t*(coeff2.y+t*coeff3.y));
-                       }
-                       /* Now we calculate the error by determing the minimal quadratic distance to the mid points. */
-                       errorsum = 0.0;
+                       fromtemp->nextcp.x = from->me.x+ftlen*fromunit.x*abSolutions[k][0];
+                       fromtemp->nextcp.y = from->me.y+ftlen*fromunit.y*abSolutions[k][0];
+                       totemp->prevcp.x = to->me.x+ftlen*tounit.x*abSolutions[k][1];
+                       totemp->prevcp.y = to->me.y+ftlen*tounit.y*abSolutions[k][1];
+                       SplineRefigure3(temp);
+                       maxerr = 0;
                        for (int i=0; i<cnt; i++) { /* Going through the mid points */
-                               error = (mid[i].p.x-approx[0].x)*(mid[i].p.x-approx[0].x)
-                               +(mid[i].p.y-approx[0].y)*(mid[i].p.y-approx[0].y);
-                               /* Above we have just initialized the error and */
-                               /* now we are going through the remaining 98 of */
-                               /* 99 points on the approximate cubic bezier: */
-                               for (int j=1; j<99; j++) {
-                                       dist = (mid[i].p.x-approx[j].x)*(mid[i].p.x-approx[j].x)
-                                       +(mid[i].p.y-approx[j].y)*(mid[i].p.y-approx[j].y);
-                                       if (dist < error)
-                                               error = dist;
-                               }
-                               errorsum += error;
-                               if (errorsum > bestError)
-                                       break;
+                               d = SplineMinDistanceToPoint(temp, &mid[i].p);
+                               if ( d>maxerr )
+                                       maxerr = d;
                        }
-                       if (errorsum < bestError) {
-                               bestError = errorsum;
-                               from->nextcp = nextcp;
-                               to->prevcp = prevcp;
+                       if (maxerr < bestError) {
+                               bestError = maxerr;
+                               from->nextcp = fromtemp->nextcp;
+                               to->prevcp = totemp->prevcp;
                        }
                }
+               printf("Dist: %lf, bestError: %lf\n", BPDist(from->me, to->me), bestError);
+               SplineFree(temp);
+               SplinePointFree(fromtemp);
+               SplinePointFree(totemp);
     } 
     return( SplineMake3(from,to));
        } else { /* original and fast function */

With this very simple single spline:

and a 50 em-unit circle nib at 0 degrees that diff prints:

Dist: 200.841498, bestError: 32.558136
Dist: 115.892443, bestError: 7.359469
Dist: 65.855591, bestError: 1.795036
Dist: 37.082453, bestError: 0.482729
Dist: 20.757076, bestError: 0.137745
Dist: 16.338124, bestError: 0.103985
Dist: 28.872463, bestError: 0.410166
Dist: 16.120087, bestError: 0.119615
Dist: 12.764187, bestError: 0.086315
Dist: 50.983328, bestError: 1.824696
Dist: 28.423327, bestError: 0.516172
Dist: 15.831987, bestError: 0.151382
Dist: 12.610325, bestError: 0.107890
Dist: 22.692114, bestError: 0.395113
Dist: 12.601926, bestError: 0.117874
Dist: 10.104111, bestError: 0.081031
Dist: 95.347409, bestError: 8.000239
Dist: 51.980911, bestError: 2.388849
Dist: 28.674112, bestError: 0.708284
Dist: 15.876248, bestError: 0.212316
Dist: 12.833232, bestError: 0.144283
Dist: 23.528066, bestError: 0.489720
Dist: 12.997261, bestError: 0.149675
Dist: 10.551406, bestError: 0.097623
Dist: 44.707466, bestError: 1.594526
Dist: 24.389037, bestError: 0.496408
Dist: 13.428352, bestError: 0.153448
Dist: 10.981090, bestError: 0.097639
Dist: 20.432920, bestError: 0.309694
Dist: 11.236985, bestError: 0.096306
Dist: 9.205415, bestError: 0.060529
Dist: 146.333312, bestError: 23.799875
Dist: 90.894818, bestError: 5.605748
Dist: 54.968005, bestError: 1.482386
Dist: 31.874513, bestError: 0.413300
Dist: 18.098605, bestError: 0.119903
Dist: 13.786824, bestError: 0.087614
Dist: 23.175440, bestError: 0.327964
Dist: 13.141845, bestError: 0.097345
Dist: 10.043044, bestError: 0.067813
Dist: 36.644712, bestError: 1.301219
Dist: 21.106278, bestError: 0.381694
Dist: 11.974910, bestError: 0.114266
Dist: 9.145421, bestError: 0.078121
Dist: 15.632831, bestError: 0.271532
Dist: 8.803081, bestError: 0.082223
Dist: 6.839331, bestError: 0.054797
Dist: 63.257646, bestError: 5.380805
Dist: 33.667908, bestError: 1.552847
Dist: 18.677151, bestError: 0.461165
Dist: 10.436447, bestError: 0.139449
Dist: 8.263775, bestError: 0.092897
Dist: 15.135292, bestError: 0.315476
Dist: 8.332523, bestError: 0.095994
Dist: 6.816048, bestError: 0.063099
Dist: 30.500877, bestError: 1.091063
Dist: 16.231148, bestError: 0.331017
Dist: 8.833141, bestError: 0.101040
Dist: 7.411607, bestError: 0.065961
Dist: 14.347985, bestError: 0.217792

You can see the "divide and conquer" at the next level up here, with the splines getting shorter. The ~200 is the attempt at the outer offset of the whole spline and the ~146 is the attempt at the inner offset of the whole spline.

The main thing I'm noticing here is that 32.6 em-units seems like a huge error for a relatively boring curve (represented with fit points) with 200 em-units between the endpoints. Same with 146 and 23.8.

Looking at the individual distance values from that 200 em-unit fit:

K: 0, d: 10.813898
K: 0, d: 19.859587
K: 0, d: 26.783539
K: 0, d: 31.174287
K: 0, d: 32.558136
K: 0, d: 30.443860
K: 0, d: 24.430493
K: 0, d: 14.302030
K: 1, d: 10.813898
K: 1, d: 19.859587
K: 1, d: 26.783539
K: 1, d: 31.174287
K: 1, d: 32.558136
K: 1, d: 30.443860
K: 1, d: 24.430493
K: 1, d: 14.302030

So what's it doing? I tried turning the accuracy down to 50 em-units to get the top-level splines and I see this:

After staring at this a bit I tried printing the distance between the control points and corresponding base points and got 1 and 0, which made me suspicious and I think I know roughly what the problem is.

Your calculation of nextcp and prevcp is:

                from->nextcp.x = from->me.x+ftlen*fromunit.x*abSolutions[0][0];
                from->nextcp.y = from->me.y+ftlen*fromunit.y*abSolutions[0][0];
                to->prevcp.x = to->me.x+ftlen*tounit.x*abSolutions[0][1];
                to->prevcp.y = to->me.y+ftlen*tounit.y*abSolutions[0][1];

ftlen is untouched from the calculation above where yours starts. I'm not sure what that's doing here but it won't work with the expand stroke input. All Expand Stroke does is set the slopes with unit tangents (or something like that). It expects the algorithm to use the fit points constrained to those slopes without assuming they have any further meaning.

It seems like the algorithm is making more assumptions than that about the input control point positions but you're the right person to judge that.

Going back to the other problem does the version with SplineMinDistanceToPoint() also have the too-many splines problem when stroking? That would be interesting to know because that's the measure being used one level up.

Yes, in my tests the version with SplineMinDistanceToPoint() has the too-many splines problem as well.

BTW that's not quite the right spline calculation -- SplineMake() or SplineRefigure() need to be run after setting the control points to correctly calculate the coefficients.

I actually already tested it before in the order that you suggest, but strangely the "to.prevcp" point is identical with the "to" point after either applying SplineMake() or SplineRefigure().

I have now again tested the following simple spline

0 14 m 0
 27 14 55 11 86 6 c 4
 125 0 144 0 151 0 c 1024

examplespline.sfd.txt with

bigreal bestError = 1e30;
bigreal maxerr,d;
SplinePoint *fromtemp = SplinePointCreate(from->me.x, from->me.y);
SplinePoint *totemp = SplinePointCreate(to->me.x, to->me.y);
Spline *temp = SplineMake3(fromtemp, totemp);
for (int k=0; k<numberOfSolutions; k++) {
	fromtemp->nextcp.x = from->me.x+ftlen*fromunit.x*abSolutions[k][0];
	fromtemp->nextcp.y = from->me.y+ftlen*fromunit.y*abSolutions[k][0];
	totemp->prevcp.x = to->me.x+ftlen*tounit.x*abSolutions[k][1];
	totemp->prevcp.y = to->me.y+ftlen*tounit.y*abSolutions[k][1]; 
	printf("nextcp (%f,%f)\n",fromtemp->nextcp.x,fromtemp->nextcp.y);
	printf("prevcp (%f,%f)\n",totemp->prevcp.x,totemp->prevcp.y);
	SplineRefigure3(temp);
	printf("nextcp (%f,%f)\n",fromtemp->nextcp.x,fromtemp->nextcp.y);
	printf("prevcp (%f,%f)\n",totemp->prevcp.x,totemp->prevcp.y);
	maxerr = 0; 
	for ( int i=0; i<cnt; ++i ) {
		d = SplineMinDistanceToPoint(temp,&mid[i].p);
		if ( d>maxerr ) 
			maxerr = d;
	}
	printf("solution %i maxerr %f\n",k,maxerr);
	if ( maxerr < bestError ) {
		bestError = maxerr;
		from->nextcp = fromtemp->nextcp;
		to->prevcp = totemp->prevcp;
	}
}
printf("Dist: %lf, bestError: %lf\n", BPDist(from->me, to->me), bestError);
SplineFree(temp);
SplinePointFree(fromtemp);
SplinePointFree(totemp);

and got:

nextcp (60.946940,14.000000)
prevcp (100.541155,0.000000)
nextcp (0.000000,14.000000)
prevcp (151.000000,0.000000)
solution 0 maxerr 1.715541
nextcp (355.669190,14.000000)
prevcp (-194.181095,0.000000)
nextcp (0.000000,14.000000)
prevcp (151.000000,0.000000)
solution 1 maxerr 1.715541
nextcp (208.308065,14.000000)
prevcp (-46.819970,0.000000)
nextcp (0.000000,14.000000)
prevcp (151.000000,0.000000)
solution 2 maxerr 1.715541
nextcp (12.004571,14.000000)
prevcp (149.483524,0.000000)
nextcp (0.000000,14.000000)
prevcp (151.000000,0.000000)
solution 3 maxerr 1.715541
Dist: 151.647618, bestError: 1.715541

As you can see, it is no surprise that the error is the same for every solution, because nextcp and prevcp changed after SplineRefigure3(temp);.

Oh, yeah, my code was only appropriate post- #4125 . Before that you have to manipulate nonextcp and noprevcp yourself (presumably by setting both to false before calling SplineRefigure3()).

Anyway I think the thing to concentrate on here is ftlen. What role is it playing/supposed to play in the calculation?

I think I am beginning to understand now. My assumption that from.me == mid[0].p and to.me == mid[cnt-1].p was wrong and ftlen should be replaced by something like BPDistance(mid[0].p,mid[cnt-1].p). I will look into it tomorrow. Thanks for your patience with me!

Okay, my assumption that mid[0].p == from.me is only true, when merging a path but not when stroking a path. I am still confused because looking at the original code I am basically returning the same as George Williams did:

The original code does some matrix magic and then gets:

to->prevcp.x = to->me.x + tlen*tounit.x; to->prevcp.y = to->me.y + tlen*tounit.y;
from->nextcp.x = from->me.x + flen*fromunit.x; from->nextcp.y = from->me.y + flen*fromunit.y;

and I do some other magic and then instead

from->nextcp.x = from->me.x+ftlen*fromunit.x*abSolutions[0][0];
from->nextcp.y = from->me.y+ftlen*fromunit.y*abSolutions[0][0];
to->prevcp.x = to->me.x+ftlen*tounit.x*abSolutions[0][1];
to->prevcp.y = to->me.y+ftlen*tounit.y*abSolutions[0][1];

where ftlen*abSolutions[][] plays the role of tlen, flen. The reason that I have to use ftlen is because abSolutions[][] refers to the solutions of the normed situation i.e. from.me = (0,0), to.me = (1,0). Hence, I have to rescale the situation by multiplying my solutions by ftlen.

Ah, OK, I've read back through the article and your code and can now offer the primary suggestion to punt on using this for Expand Stroke and the secondary suggestion to perhaps rethink where you're adding this algorithm to the FontForge code.

Discussing the second suggestion should start to clarify the first. The basic issue is that this algorithm isn't really fitting the passed points, it's just using the points as a measure of or check on the fit to the original splines. The easiest way to see this is just to look at the calculation of f and m: Those are the cumulative areas and moments of the passed splines and are not calculated from the fit points.

To use this method with the stroke code we'd either need to: a) come up with a parametric or other closed-form expression of the offset curve from which to calculate f and m or b) generate more fit points (easy enough) and then calculate the area and moment from those points. "b" is probably tractable but I suspect it will wind up being relatively expensive. Given that the stroke system supports arbitrary nibs "a" isn't necessarily impossible but it's very tricky at a minimum.

Accordingly the current algorithm isn't really fitting the points and therefore doesn't obviously belong integrated into ApproximateSplineFromPointsSlopes(). As implemented it's strictly a simplifier: it relies on being passed a spline sequence that it hopes to represent with a single spline. The fit points are only used as a sanity check.

If you wanted you could play around with a linear approximation of f and m from the existing fit points. If you treat the points as connected by lines it shouldn't be too hard to calculate an area (although you would have to be careful about loops and overhangs and such). I dunno about moment -- haven't thought about it.

Just to clarify the scope of the Levien approach a bit:

How does the system generate those stroke fit points? Well, the source curve being stroked takes the form of a cubic spline, so that part is easy. When you're stroking with a circular nib you can just calculate the normal from t and add the radius to it for the offset, and that can (probably?) given a closed form solution easily enough. When stroking with a polygonal nib you don't even fit points -- those are just spline segment offsets connected with lines.

However, the FontForge system in particular doesn't actually stroke with a circular nib even when it claims to -- it strokes with a four-cubic-spline circle approximation (for reasons of generality). There are a lot of heuristics to determine when fitting should start and end so the problem boils down to stroking the source spline with the relevant spline on the nib. The fit points are calculated this way:

1. Pick a t for the source spline.
2. Calculate the position on the spline at t.
3. Calculate the tangent angle of the source spline at t.
4. Find the t' value on the nib spline that has the same tangent angle.
5. Calculate the position on the nib spline at t'.
6. Treat that position as a vector that is added to the source-spline position calculated in step 2

(This sort of makes it sound as if the "center" of the nib is assumed to be at 0,0 as that is kind of what treating a position as a vector amounts to. But it's more that the center of the nib per se is irrelevant -- the offset curve will just be displaced from the source curve by the distance between the nib and 0,0. If you don't want the curve offset centering the nib is your problem.)

It's the transition from 3 to 4 that makes a closed form tricky. You would basically need to reparameterize the nib spline by tangent angle. That's guaranteed to be 1:1 because nibs are required to be convex, but if there's a straightforward equation mapping tangent angle to offset for a cubic spline I haven't found it yet.

a) come up with a parametric or other closed-form expression of the offset curve from which to calculate f and m.

That would be very very hard if not impossible.

b) generate more fit points (easy enough) and then calculate the area and moment from those points. "b" is probably tractable but I suspect it will wind up being relatively expensive.

I agree. Maybe we should use Levien's algorithm exclusively for merging and simplifying (see my last commit)?

There's a remaining wrinkle with this new change. The current stroking algorithm is much, much faster than the previous one but there's nothing dictating any given speed. I was planning to add a knob that lets it use your more accurate fitting algorithm for users who want to reduce the number of fit-generated points.

Unfortunately that algorithm is gone now, which is a shame given the more narrow scope of the new one. Could we get it back, and maybe use an enum (or whatever) to choose between the options?

Adding the mid-old brute force algorithm as a third alternative is surely no problem. But I think it is not necessary, because simplify can be applied after stroking. Wouldn't it be better if your "knob" refers to the exactness of the simplify function?

There's already an accuracy parameter that feeds into both the fitting algorithm and the simplify post-processing (if that's turned on).

I guess that's one way to go. There are some subtle potential advantages to doing it the other way but they can be duplicated with some work.

Taking enum pointtype { pt_curve, pt_corner, pt_tangent, pt_hvcurve }; as antetype, I have added enum mergetype { mt_matrix, mt_levien, mt_bruteforce }; and described it as:

• mt_matrix; original, fast, all-purpose,
• mt_levien; by Raph Levien, fast, accurate, use iff mid is on spline,
• mt_bruteforce; slow, all-purpose, normally more accurate than mt_matrix.

At the moment, mt_bruteforce is not used anywhere but as far as I have tested it should work, too.

Please wait for the merge of this patch. I have discovered a "merge-to-single-cyclic-segment"-issue with the new algorithm that was not around with the old algorithm:

If we merge the two splines by removing the smooth point, we will get a single point instead of a spline. The problem is obvious: Norming the spline such that from == (0,0) and to == (1,0) is impossible. I will have a look into it.

With the last commit, the "merge-to-single-cyclic-segment"-issue is solved.

I once had tested the original _QuarticSolve() versus my primitive newton solver and had found an occurrence where the newton solver found more roots than the _QuarticSolve(). Now I have thoroughly tested both solvers with thousands of random quartics and there was no situation where _QuarticSolve() was worse than the newton solver (but the opposite happened). So I guess my old test was due to a different problem. Hence, I have removed my newton solver and have made _QuarticSolve() from splineuitl.c public.

Sorry -- I'll try to get back to reviewing this shortly.

I lied (obviously). I got pulled into some work-work.

There's another spate of reviewing on the near horizon and I'll help out with that as I can. I fully expect this to be merged before the next release of FontForge.