The problem of curve fitting is fundamental to font technology, as we want to make Béziers which most closely resemble the "true" shape of the glyph. Font tools need to apply curve fitting to simplify outlines, apply transformations such as offset curve, delete a smooth on-curve point, and other applications. However, the problem is surprisingly tricky and there is no definitive solution in the literature. I recently worked out a much better solution and wrote about it in my blog:
https://raphlinus.github.io/curves/2021/03/11/bezier-fitting.htmlThe main reason it's so tricky is that C-shaped Bézier curves come in triplets of similar shape, causing local minima so that an optimizer may fall into one of these minima, even if one of the other two has lower error. My blog post gives a relatively simple formula for finding all three minima.
One subtle detail is that the goal of optimizing Bézier curves for shape accuracy is in tension with making them interpolation compatible. To me, this is a strike against cubic Béziers as a fundamental shape representation for fonts, and I continue to explore techniques that have smoothness guarantees by construction. But hopefully that's a topic for a future blog post.
Comments
As far as I understand, the first moment is determined only for the approximation curve without any connection to the source curve (delta1 may depend on the area of the source curve, but delta0 is independent). How can it be a measure for the error then?
At some point I'll release code (still working on it). I could also polish up my Python notebook a bit and release that, but at present it's in a pretty messy state. In any case, I'd be thrilled to see it start to be applied in type tools.
Choosing the solution close to the other moment probably works well, but I haven't tested it carefully. In my notebook, and likely in the version I implement in kurbo, I'm computing an L2 error norm based on the normalized arc length parameterization of the two curves.
Btw, while experimenting with the curve algorithms I found that a great way to visualize before/after difference is to make a closed path with the original and modified curve:
That is possible to do with 2 overlapping colored contours one of which is semi-transparent. Can be done in modern tools, but the b/w option is much easier to do (duplicate segment, drag back to auto-join at the ends).
12*ca*((((cb*(30*f*cb-sb)-15*f)*ca+2*sa-cb*sa*(cb+30*f*sb))*ca+cb*(sb-15*f*cb))*ca-sa*cb*cb),
12*ca*((((70*m+15*f)*sb*sb+cb*(9*sb-70*cb*m-5*cb*f))*ca-5*sa*sb*(3*sb-4*cb*(7*m+f)))*ca-cb*(9*sb-70*cb*m-5*cb*f)),
16*ca*(((12*sa-5*ca*(42*m-17*f))*sb-70*cb*(3*m-f)*sa-75*ca*cb*f*f)*sb-75*cb*cb*f*f*sa),
80*ca*sb*(42*sb*m-25*f*(sb-cb*f))]
12*ca^3*(2*sa-15*ca*f),
12*ca^2*(ca*(70*m+15*f)-15*sa),
16*ca*(12*sa-5*ca*(42*m-17*f)),
80*ca*(42*m-25*f)]
So after I either clean up my Python notebook so I can release it, or get far enough into my Rust implementation to try your two-bezier input, I'd be able to compare more directly. In the meantime, I would be slightly suspicious that the area and moment of your input are exactly correct, as the asymmetry in output resembles what would happen if that were the case. In particular, if the input curve is symmetrical, the x-moment should be half the area; when I perturb the area slightly but change x-moment to match, both my numerical solver and your quartic agree on the result, which is symmetrical.
Mateusz: that's a very complex question. The previous blog post had an O(n^5) fit from Euler spirals (not based on moments or areas) that is not a problem for interpolation compatibility, but I think an O(n^6) solution will inevitably hit different branches.
More seriously, I'm pleased and impressed you were able to get this to work. Your derivation is better than the direction I was headed - I was dividing by sin(θ₀+θ₁) to resolve the δ₀δ₁ terms, and was planning to special-case the case where that approaches zero. That seems to cancel out, and your coefficients seem pretty solid.
As I mentioned, the remaining case that I'm finding tricky is the "near misses." What's really bothering me there is that if we compute the roots of the derivative of that polynomial (ie solve the resulting cubic) we should get pretty good results, but it's not symmetric to a left-right swap. This is why I'm still leaning to numerically finding the zero-crossing of the derivative of the x-moment. I can think of a couple ways to do that, but haven't yet settled on one.
This feels like a real open source collaboration!
- One should norm θ₀ to being non-negative (what you may have done anyway but I did not see any mathematical necessity). The numerical computations unfortunately seem to depend on the signs...
- The case θ₁ = 180° - θ₀ leads to bad approximations and should be handled separately. This case happens quite often in type design: Think of the outer outline of an /o that consists of two horizontal and two vertical on-curve points. (In the end one probably want to add extrema again.)
- There are situations (e.g. happening on the /o of the bold Euler Fraktur eufb10.pfb) where both the zeroes of the x-moment and the zeroes of its derivative are suboptimal (because δ1 is negative for them)! E.g. the x-moment has the function `f(x):=0.113110x^4-0.334960x^3+0.302089x^2-0.140921x+0.048084`(where x stands for δ0) and `δ1(x)=(20*0.024743-6*x*sin(0.116270))/(3*(2*sin(0.220310)-x*sin(0.116270+0.220310)))` and it looks like this (green = x-moment, orange = δ1):

I have solved the 3. problem by additionally considering cases where δ1 is approximately 0 (exactely 0 does not preserve the direction, I took 0.01 instead).Looks impressive. ~45/60 reduction in points. If you count the pixels in the 3rd, rightmost image, you get the difference of the areas. Can also be done by xor of the two rendered images. Would be interesting how large the relative area-difference and the relative thickness difference is (pixel_difference / em_size).
In theory 1 unit / 2000 M-size = 0.0005 is smaller than visual resolution 0.15 mm / 250 mm reading distance = 0.0006.
The interesting part for me is a threshold, which we can use in software.
Here is a showcase between https://fonts.google.com/specimen/Castoro?preview.text_type=custom&preview.text=iftſtfiſißſsſijA and Castoro downloaded from github/TiroTypeworks TTF-Version.
Left is Google rendered by Chrome right is Tiro rendered in Mac TextEdit. Look without cheating first at size 14 for the difference, then next sizes.
For the readers not detecting it: The sequence longs_i is different between the two versions of the font.
Newton solvers are tempting due to their high rate of convergence, and, in this case, ease of computing the derivative analytically. However, they go unstable when there are multiple roots (or nearly so), and that's also a case that can happen here. I've been exploring the ITP method and have it implemented in kurbo now, as it's extremely robust numerically as well as converging reasonably fast. I'm still chipping away at the kurbo implementation of curve fitting (I really want to get it perfect) and haven't made final decisions, but am leaning strongly toward ITP as the main solving technique.
Let me see if I can explain the left-right swap issue more clearly. If you're only looking at zeros of the x-moment error, those don't move if you multiply it by a factor of (aθ₀ + b)², which is what we do to get the quartic polynomial. Another way of phrasing this is that the x-moment error can be expressed as a rational polynomial with a quartic on top and that quadratic on bottom. And thus the zeros also don't move if you swap θ₀ and θ₁. But if you look at extrema of the x-moment error vs extrema of just the quartic (the latter of which is two different versions depending on whether it's expressed in terms of δ₀ or δ₁), all three are in different places. Yes, you could make it symmetrical by computing both of the latter extrema and choosing the minimal error, but don't you think it's more satisfying to pick the extremum of the x-moment itself?
ya < 0 < yb
quite difficult for a general quartic?Thank you for pointing this out. Indeed, I have (falsely) assumed that the zeroes of the derivative of the quartic are the zeroes of the derivative of the moment. But of course you are right. However, the derivative of the moment ends up in a quartic as well.
The quartic for the roots of the derivative of the moment is:
where ca = cos(θ₀), sa = sin(θ₀), cb = cos(θ₁), sb = sin(θ₁), sab = sin(θ₀ + θ₁), f = area, a = δ₀.This quartic equation reduces to a linear for the special cases θ₀ + θ₁ = 180° or θ₀ + θ₁ = 0°:
θ₁ = 12.5°
area = 0.025
moment = 0.013
zeroes (θ₀,θ₁):
normal zero (-0.238080,0.432158)
normal zero (1.808710,1.556856)
derivative zero (0.381235,0.260242)
derivative zero (0.497571,0.199395)
θ₀ = 12.5°
θ₁ = 6.5°
area = 0.025
moment = 0.013
zeroes (θ₀,θ₁):
normal zero (0.501812,-0.802088)
normal zero (1.420217,1.899029)
derivative zero (0.260242,0.381235)
derivative zero (0.199395,0.497571)
I cannot understand how anybody could disagree with that... Let me assure that both statements are correct
Regarding poor approximations when the sum of the angles is 180°, I think in this case the equation becomes quadratic and it's best to just set the x^3 and x^4 terms to zero rather than try to run them through a quartic solver with epsilon values - my current solver code also doesn't handle this well, though I'm looking at ways of making it more robust. Also I can heartily recommend https://cristiano-de-michele.netlify.app/publication/orellana-2020/orellana-2020.pdf as a good approach to quartic solving.
I will post another update and some code before too long. The main reason I'm looking at this now is that I think it will be the basis of a most excellent solution to stroke offset.