An arc of up to 90 degrees can be approximated by a single cubic Bezier spline reasonably well. To find an appropriate approximating spline we set the start and end points of the spline to match the arc and then can find the middle two points using the following formula, where A is the start angle and B is the end angle:h = 4/3 * tan ((B - A)/4)
pt1 = {center.x + r*cos(A), center.y + r*sin(A)}
pt2 = {center.x + r*cos(A) - h*r*sin(A), center.y + r*sin(A) + h*r*cos(A)}
pt3 = {center.x + r*cos(A) + h*r*sin(B), center.y + r*sin(A) - h*r*cos(B)}
pt4 = {center.x + r*cos(B), center.y + r*sin(B)}
This approach works quite well and is currently used by cairo. However, the problem with this approach is that it requires converting to polar coordinates and then back to Cartesian ones. This conversion is slow because it requires the use of trigonometric functions.
Fortunately, this paper gives a solution that avoids the conversion to polar coordinates:
ax = pt1.x - center.x
ay = pt1.y - center.y
bx = pt4.x - center.x
by = pt4.y - center.y
q1 = ax*ax + ay*ay
q2 = q1 + ax*bx + ay*by
k2 = 4/3 (sqrt(2*q1*q2) - q2) / (ax*by - ay*bx)
pt2.x = pt1.x - k2*ay
pt2.y = pt1.y + k2*ax
pt3.x = pt4.x + k2*bx
pt3.y = pt4.y - k2*by
Unfortunately, there is no explanation, derivation or proof for this formula. Even more problematic, the formula for k2 becomes less stable as pt1 approaches pt4, becoming NaN when pt1 equals pt4.
Therefore, the homework problem has two parts:
1. Give an explanation or derivation for the formula for k2 provided above.
2. Provide a similar formulation for pt2 and pt3 that doesn't degenerate as pt1 and pt4 become close or prove that one doesn't exist.
The best answer will be credited in the new cairo stroking code.
17 comments:
I spent a bunch of time staring at this problem for rounded rectangle drawing in Mozilla's Thebes wrapper. You might find the big long comment starting at http://mxr.mozilla.org/mozilla-central/source/gfx/thebes/src/gfxContext.cpp#806 helpful, and the papers cited there (especially the one by Michael Goldapp; doi:10.1016/0167-8396(91)90007-X). Unfortunately, it's not available online for free. And I have no idea about the numerical stability when the two arc endpoints are close.
Incidentally, it would be great if Cairo had direct support for SVG elliptical arcs as path segments...
First off, your formulas for pt2 and pt3 don't match the paper's. The paper has, for example:
pt2.x = pt1.x + center.x - k2*pt1.y
This is clearly bogus, if I understand how this stuff is supposed to work, since that gives the wrong direction at pt1. Your formula for pt2 is more likely to be correct (though your formula for pt3 seems to have bx and by switched).
OK. So the formula for k2 there gives exactly the same value as the formula for h, up to sign. I didn't check the signs, but they should be checkable without too much trouble.
Proof: The 4/3 is the same in both formulas and can be dropped for purposes of comparing the two values. If R is the radius of the circle, then we have:
q1 = R^2
q2 = R^2 + R^2 cos (B - A)
(since q1 - q1 is just the dot product of (ax,ay) and (bx,by), and the two vectors both have length R and have angle B-A between them.
|ax*by - ay*bx| is just the length of the cross product of the two vectors considered as three-dimensional vectors, and has magnitude R^2*|sin(B-A)|. I'm going to drop the absolute value; this is where you'd need to check your signs.
Now we have (ignoring the 4/3):
k2 = (sqrt(2*R^2*R^2(1 + cos(B-A))) - R^2 - R^2*cos(B-A))/(R^2*sin(B-A).
All the R^2 go away, and we are left with:
k2 = (sqrt(2(1+cos(B-A)) - 1 - cos(B-A))/sin(B-A)
Now we use the identity: 1 + cos(x) = 2cos^2(x/2) to get:
k2 = (2cos((B-A)/2) - 2cos^2((B-A)/2)) / sin(B-A)
and then use sin(2x) = 2sin(x)cos(x) to get:
k2 = (1 - cos((B-A)/2))/(sin((B-A)/2)
Now use 1 - cos(x) = 2sin^2(x/2) and the sin(2x) formula again to get:
k2 = sin((B-A)/4)/cos((B-A)/4)
and we're done.
At least an explanation, right? ;)
As a note, you could turn my explanation into a derivation by simply running the steps backwards, as long as |(B-A/2)| < pi/2 and B-A != 0.
I'll try the challenge...but is there a reason you're using cubic splines? Jim Blinn's excellent "How many ways can you draw a circle?" gives a great selection of circle (and arc) drawing algorithms, most of them for vector displays...perfect if you need to decompose the arc into line segments.
Jeff, what's the actual input into the algorithm here? Is it the center and the two arc endpoints? Or the two angles A and B? Or something else?
rethinking this, maybe you use cubic splines because you can transform just the 4 points? if you do want to check out Blinn's article: page 5 and 6 have two of the better algorithms, but...
Google Books: How Many Ways Can You Draw A Circle?
Boris, the cairo api has
void cairo_arc(cairo_t *cr, double xc, double yc, double radius, double angle1, double angle2);
Boris,
The input is the center and the two end points for the arc.
Brendan,
This code will not necessarily be used for implementing cairo_arc because that already use polar coordinates. The Cartesian version is for round joins and caps.
@Brendan: We convert arcs to cubic bezier splines because that's what's supported by cairo's path format. There's been talk of keeping arcs as arcs instead of approximating them by splines. However, the additional curve type adds complexity.
Interestingly enough, Postscript had support for arcs but this was removed in PDF.
@Zack,
I've quickly skimmed through the paper by Goldapp. He gives the tan(a/4) formulation in the "Hermite-type approximation" section. The other formulations, although potentially more accurate, seemed more complicated than they were worth.
Also, as far as I could tell, the formulations he gives all use polar coordinates?
@Jeff,
You can derive from the "Hermite type approximation" a constant α and and some formulas that let you calculate any given quarter-circle without any polar coordinates. A diagram would probably help: see here Basically, given the center point C and the desired radius, you can derive the other three points of the bounding square, A, B, D. B and D are the Bezier endpoints, and the constant α lets you calculate the control points P and Q trivially - the notation α(D→A) means "form the vector from D to A, then multiply it by α".
You then have to cut up the Bezier if you don't want the whole quarter-circle, but I believe that can be done with no trig either.
Last night Boris suggested bisecting the angle twice and then computing the tan of it. This turns into the following, which seems to work pretty well:
vector_t mid;
mid.x = a.x + b.x;
mid.y = a.y + b.y;
double l = sqrt(mid.x*mid.x + mid.y*mid.y);
mid.x /= l;
mid.y /= l;
vector_t mid2;
mid2.x = a.x + mid.x;
mid2.y = a.y + mid.y;
k2 = (4./3.)*(-a.y*mid.x + a.x*mid.y)/(a.x*mid2.x + a.y*mid2.y);
Where 'a' and 'b' are the vectors from the center to the start and end points respectively
mid2, not mid, in the numerator in the k2 calculation there at the end, right?
Yep, it should be mid2. I put that error in when manually inlining dot() and perp() to post the solution here. Interestingly enough it still passed the tests that I was running.
It would be very interesting to see the performance numbers comparing this algorithm to the existing code.
A quick and dirty comparison puts the polar coordinate version at about 1050 cycles and the new Cartesian one at about 360 cycles. So just under 3x faster.
Post a Comment