mirror of
https://github.com/Pomax/BezierInfo-2.git
synced 2025-08-27 10:15:05 +02:00
full circle chapters rewrite
This commit is contained in:
@@ -1,16 +1,10 @@
|
||||
let guess, w, h, pad = 75, r;
|
||||
let curve, w, h, pad = 75, r;
|
||||
|
||||
setup() {
|
||||
w = this.width;
|
||||
h = this.height;
|
||||
r = w/2 - pad;
|
||||
guess = new Bezier(this,
|
||||
{ x: w - pad, y: h/2},
|
||||
{ x: w - pad, y: h/2},
|
||||
{ x: w - pad, y: h/2},
|
||||
{ x: w - pad, y: h/2}
|
||||
);
|
||||
setSlider(`.slide-control`, `angle`, -PI/4, v => this.updateCurve(v));
|
||||
w = this.width;
|
||||
h = this.height;
|
||||
r = w/2 - pad;
|
||||
setSlider(`.slide-control`, `angle`, 1.4);
|
||||
}
|
||||
|
||||
draw() {
|
||||
@@ -28,58 +22,24 @@ draw() {
|
||||
let a = this.angle;
|
||||
wedge(w/2, h/2, r, a < 0 ? a : 0, a < 0 ? 0 : a);
|
||||
|
||||
guess.drawSkeleton(`lightblue`);
|
||||
guess.drawCurve(`lightblue`);
|
||||
|
||||
let real = this.getRealCurve(guess.points[0], guess.points[3], this.angle);
|
||||
real.drawSkeleton();
|
||||
real.drawCurve();
|
||||
|
||||
setColor(`black`);
|
||||
real.points.forEach(p => {
|
||||
this.updateCurve(this.angle);
|
||||
curve.drawSkeleton();
|
||||
curve.drawCurve();
|
||||
curve.points.forEach(p => {
|
||||
circle(p.x, p.y, 2);
|
||||
text(`(${p.x|0},${p.y|0})`, p.x+5, p.y);
|
||||
});
|
||||
}
|
||||
|
||||
updateCurve(a) {
|
||||
let angle = -a;
|
||||
|
||||
const S = guess.points[0],
|
||||
B = {
|
||||
x: w/2 + r * cos(angle/2),
|
||||
y: h/2 + r * sin(angle/2)
|
||||
},
|
||||
E = guess.points[3] = {
|
||||
x: w/2 + r * cos(angle),
|
||||
y: h/2 + r * sin(angle)
|
||||
};
|
||||
|
||||
guess = this.guessCurve(S,B,E);
|
||||
|
||||
return angle;
|
||||
}
|
||||
|
||||
guessCurve(S, B, E) {
|
||||
const C = { x: (S.x + E.x)/2, y: (S.y + E.y)/2 }, // we know we're working with t=0.5
|
||||
A = { x: B.x + (B.x-C.x)/3, y: B.y + (B.y-C.y)/3 }, // cubic ratio at t=0.5 is 1/3
|
||||
bx = (E.x-S.x)/4,
|
||||
by = (E.y-S.y)/4,
|
||||
e1 = { x: B.x - bx, y: B.y - by },
|
||||
e2 = { x: B.x + bx, y: B.y + by },
|
||||
v1 = { x: A.x + (e1.x-A.x)*2, y: A.y + (e1.y-A.y)*2 },
|
||||
v2 = { x: A.x + (e2.x-A.x)*2, y: A.y + (e2.y-A.y)*2 },
|
||||
C1 = { x: S.x + (v1.x-S.x)*2, y: S.y + (v1.y-S.y)*2 },
|
||||
C2 = { x: E.x + (v2.x-E.x)*2, y: E.y + (v2.y-E.y)*2 };
|
||||
return new Bezier(this, [S, C1, C2, E]);
|
||||
}
|
||||
|
||||
getRealCurve(S, E, angle) {
|
||||
updateCurve(angle) {
|
||||
const f = 4/3 * tan(angle/4);
|
||||
const S = { x: w/2 + r, y: h/2 };
|
||||
const C1 = { x: w/2 + r, y: h/2 + r * f };
|
||||
const C2 = {
|
||||
x: w/2 + r * (cos(angle) + f * sin(angle)),
|
||||
y: h/2 + r * (sin(angle) - f * cos(angle))
|
||||
};
|
||||
return new Bezier(this, [S, C1, C2, E]);
|
||||
const E = { x: w/2 + r * cos(angle), y: h/2 + r * sin(angle) };
|
||||
curve = new Bezier(this, [S, C1, C2, E]);
|
||||
}
|
||||
|
119
docs/chapters/circles_cubic/arc-basics.js
Normal file
119
docs/chapters/circles_cubic/arc-basics.js
Normal file
@@ -0,0 +1,119 @@
|
||||
let curve, w, h, pad = 75, r;
|
||||
|
||||
setup() {
|
||||
w = this.width;
|
||||
h = this.height;
|
||||
r = w - 2 * pad;
|
||||
curve = new Bezier(this,
|
||||
// { x: r, y: h/2},
|
||||
{ x: r, y: h/2},
|
||||
{ x: r, y: h/2},
|
||||
{ x: r, y: h/2}
|
||||
);
|
||||
setSlider(`.slide-control`, `angle`, 1.57, v => this.updateCurve2(v));
|
||||
}
|
||||
|
||||
draw() {
|
||||
clear();
|
||||
|
||||
translate(1.5 * pad,pad);
|
||||
setColor(`grey`);
|
||||
line(0,-h,0,h);
|
||||
line(-w,0,w,0);
|
||||
curve.drawSkeleton(`grey`);
|
||||
let p = curve.drawStruts(0.5, `grey`, false);
|
||||
curve.drawCurve(`black`);
|
||||
// this.drawCubicInformation(p);
|
||||
this.drawQuadraticInformation(p);
|
||||
}
|
||||
|
||||
drawQuadraticInformation(p) {
|
||||
noFill();
|
||||
setStroke(`black`)
|
||||
let a = this.angle;
|
||||
let ar = 20
|
||||
arc(0,0,ar,0,a);
|
||||
setStroke(`grey`);
|
||||
let br = 4 * ar;
|
||||
arc(0,0,br,0,a/2);
|
||||
line(0,0, p[2].x, p[2].y);
|
||||
line(0,0, p[1].x, p[1].y);
|
||||
|
||||
setFill(`black`);
|
||||
setTextStroke(`white`, 3);
|
||||
text(`P1`, p[0]);
|
||||
text(`P2`, p[1]);
|
||||
text(`P3`, p[2]);
|
||||
text(`B`, p[5]);
|
||||
|
||||
text(`θ`, 1.2 * ar * cos(a / 3), 1.2 * ar * sin(a / 3));
|
||||
text(`θ/2`, 1.2 * br * cos(a / 4), 1.2 * br * sin(a / 4));
|
||||
}
|
||||
|
||||
|
||||
drawCubicInformation(p) {
|
||||
noFill();
|
||||
setStroke(`black`)
|
||||
let a = this.angle;
|
||||
let ar = 20
|
||||
arc(0,0,ar,0,a);
|
||||
setStroke(`grey`);
|
||||
let br = 4 * ar;
|
||||
arc(0,0,br,0,a/2);
|
||||
line(0,0, p[3].x, p[3].y);
|
||||
line(0,0, p[9].x, p[9].y);
|
||||
|
||||
setFill(`black`);
|
||||
setTextStroke(`white`, 3);
|
||||
text(`P1`, p[0]);
|
||||
text(`P2`, p[1]);
|
||||
text(`P3`, p[2]);
|
||||
text(`P4`, p[3]);
|
||||
text(`A`, p[5]);
|
||||
text(`B`, p[9]);
|
||||
text(`e1`, p[7]);
|
||||
text(`e2`, p[8]);
|
||||
|
||||
text(`θ`, 1.2 * ar * cos(a / 3), 1.2 * ar * sin(a / 3));
|
||||
text(`θ/2`, 1.2 * br * cos(a / 4), 1.2 * br * sin(a / 4));
|
||||
}
|
||||
|
||||
updateCurve3(a) {
|
||||
this.angle = a;
|
||||
let p = curve.points;
|
||||
let d = r * 4/3 * tan(a/4);
|
||||
|
||||
// start point
|
||||
p[0].x = r;
|
||||
p[0].y = 0;
|
||||
|
||||
// end point
|
||||
p[3].x = r * cos(a);
|
||||
p[3].y = r * sin(a);
|
||||
|
||||
// first control point
|
||||
p[1].x = r;
|
||||
p[1].y = d;
|
||||
|
||||
// second control point
|
||||
p[2].x = p[3].x + d * sin(a);
|
||||
p[2].y = p[3].y - d * cos(a);
|
||||
}
|
||||
|
||||
updateCurve2(a) {
|
||||
this.angle = a;
|
||||
let p = curve.points;
|
||||
let d = r * tan(a/2);
|
||||
|
||||
// start point
|
||||
p[0].x = r;
|
||||
p[0].y = 0;
|
||||
|
||||
// control point
|
||||
p[1].x = r;
|
||||
p[1].y = d;
|
||||
|
||||
// end point
|
||||
p[2].x = r * cos(a);
|
||||
p[2].y = r * sin(a);
|
||||
}
|
@@ -1,185 +1,309 @@
|
||||
# Circles and cubic Bézier curves
|
||||
# Circular arcs and cubic Béziers
|
||||
|
||||
In the previous section we tried to approximate a circular arc with a quadratic curve, and it mostly made us unhappy. Cubic curves are much better suited to this task, so what do we need to do?
|
||||
|
||||
For cubic curves, we basically want the curve to pass through three points on the circle: the start point, the mid point at "angle/2", and the end point at "angle". We then also need to make sure the control points are such that the start and end tangent lines line up with the circle's tangent lines at the start and end point.
|
||||
|
||||
The first thing we can do is "guess" what the curve should look like, based on the previously outlined curve-through-three-points procedure. This will give use a curve with correct start, mid and end points, but possibly incorrect derivatives at the start and end, because the control points might not be in the right spot. We can then slide the control points along the lines that connect them to their respective end point, until they effect the corrected derivative at the start and end points. However, if you look back at the section on fitting curves through three points, the rules used were such that they optimized for a near perfect hemisphere, so using the same guess won't be all that useful: guessing the solution based on knowing the solution is not really guessing.
|
||||
|
||||
So have a graphical look at a "bad" guess versus the true fit, where we'll be using the bad guess and the description in the second paragraph to derive the maths for the true fit:
|
||||
Let's look at approximating circles and circular arcs using cubic Béziers. How much better is that?
|
||||
|
||||
<graphics-element title="Cubic Bézier arc approximation" width="400" height="400" src="./arc-approximation.js">
|
||||
<input type="range" min="-3.1415" max="3.1415" step="0.01" value="-0.7854" class="slide-control">
|
||||
<input type="range" min="-3.1415" max="3.1415" step="0.01" value="1.4" class="slide-control">
|
||||
</graphics-element>
|
||||
|
||||
We see two curves here; in blue, our "guessed" curve and its control points, and in grey/black, the true curve fit, with proper control points that were shifted in, along line between our guessed control points, such that the derivatives at the start and end points are correct.
|
||||
At cursory glance, a fair bit better, but let's find out _how much_ better by looking at how to construct the Bézier curve.
|
||||
|
||||
We can already see that cubic curves are a lot better than quadratic curves, and don't look all that wrong until we go well past a quarter circle; ⅜th starts to hint at problems, and half a circle has an obvious "gap" between the real circle and the cubic approximation. Anything past that just looks plain ridiculous... but quarter curves actually look pretty okay!
|
||||

|
||||
|
||||
So, maths time again: how okay is "okay"? Let's apply some more maths to find out.
|
||||
The start and end points are trivial, but the mid point requires a bit of work, but it's mostly basic trigonometry once we know the angle θ for our circular arc: if we scale our circular arc to a unit circle, we can always start our arc, with radius 1, at (1,0) and then given our arc angle θ, we also know that the circular arc has length θ (because unit circles are nice that way). We also know our end point, because that's just (cos(θ), sin(θ)), and so the challenge is to figure out what control points we need in order for the curve at _t_=0.5 to exactly touch the circular arc at the angle θ/2:
|
||||
|
||||
Unlike for the quadratic curve, we can't use <i>t=0.5</i> as our reference point because by its very nature it's one of the three points that are actually guaranteed to lie on the circular curve. Instead, we need a different <i>t</i> value. If we run some analysis on the curve we find that the actual <i>t</i> value at which the curve is furthest from what it should be is 0.211325 (rounded), but we don't know "why", since finding this value involves root-finding, and is nearly impossible to do symbolically without pages and pages of math just to express one of the possible solutions.
|
||||
|
||||
So instead of walking you through the derivation for that value, let's simply take that <i>t</i> value and see what the error is for circular arcs with an angle ranging from 0 to 2π:
|
||||
|
||||
<table><tbody><tr><td>
|
||||
<img src="images/arc-c-2pi.gif" height="187"/>
|
||||
plotted for 0 ≤ φ ≤ 2π:
|
||||
</td><td>
|
||||
<img src="images/arc-c-pi.gif" height="187"/>
|
||||
plotted for 0 ≤ φ ≤ π:
|
||||
</td><td>
|
||||
<img src="images/arc-c-pi2.gif" height="187"/>
|
||||
plotted for 0 ≤ φ ≤ ½π:
|
||||
</td></tr></tbody></table>
|
||||
|
||||
We see that cubic Bézier curves are much better when it comes to approximating circular arcs, with an error of less than 0.027 at the two "bulge" points for a quarter circle (which had an error of 0.06 for quadratic curves at the mid point), and an error near 0.001 for an eighth of a circle, so we're getting less than half the error for a quarter circle, or: at a slightly lower error, we're getting twice the arc. This makes cubic curves quite useful!
|
||||
|
||||
In fact, the precision of a cubic curve at a quarter circle is considered "good enough" by so many people that it's generally considered "just fine" to use four cubic Bézier curves to fake a full circle when no circle primitives are available; generally, people won't notice that it's not a real circle unless you also happen to overlay an actual circle, so that the difference becomes obvious.
|
||||
|
||||
So with the error analysis out of the way, how do we actually compute the coordinates needed to get that "true fit" cubic curve? The first observation is that we already know the start and end points, because they're the same as for the quadratic attempt:
|
||||
|
||||
\[ S = \begin{pmatrix} 1 \\ 0 \end{pmatrix} ~, ~\ E = \begin{pmatrix} cos(φ) \\ sin(φ) \end{pmatrix} \]
|
||||
|
||||
But we now need to find two control points, rather than one. If we want the derivatives at the start and end point to match the circle, then the first control point can only lie somewhere on the vertical line through S, and the second control point can only lie somewhere on the line tangent to point E, which means:
|
||||
So let's again formally describe this:
|
||||
|
||||
\[
|
||||
C_1 = S + a \cdot \begin{pmatrix} 0 \\ 1 \end{pmatrix}
|
||||
\begin{aligned}
|
||||
P_1 &= (1, 0) \\
|
||||
P_2 &= (1, c) \\
|
||||
P_3 &= P_4 + k \cdot (sin(θ), -cos(θ)) \\
|
||||
P_4 &= (cos(θ), sin(θ))
|
||||
\end{aligned}
|
||||
\]
|
||||
|
||||
where "a" is some scaling factor we'll need to find the expression for, and:
|
||||
Only P<sub>3</sub> isn't quite straight-forward here, and its description is based on the fact that the triangle (origin, P<sub>4</sub>, P<sub>3</sub>) is a right angled triangle, with the distance between the origin and P<sub>4</sub> being 1 (because we're working with a unit circle), and the distance between P<sub>4</sub> and P<sub>3</sub> being _c , so that we can represent P<sub>3</sub> as "The point P<sub>4</sub> plus the vector from the origin to P<sub>4</sub> but then rotated a quarter circle, counter-clockwise, and scaled by _c_".
|
||||
|
||||
With that, we can determine the _y_-coordinates for A, B, e<sub>1</sub>, and e<sub>2</sub>, after which we have all the information we need to determine what the value of _k_ is. We can find these values by using (no surprise here) linear interpolation between known points, as A is midway between P<sub>2</sub> and P<sub>3</sub>, e<sub>1</sub> is between A and "midway between P<sub>1</sub> and P<sub>2</sub>" (which is "half height" P<sub>2</sub>), and so forth:
|
||||
|
||||
\[
|
||||
C_2 = E + a \cdot \begin{pmatrix} sin(φ) \\ cos(φ) \end{pmatrix}
|
||||
\begin{aligned}
|
||||
A_y &= \frac{P_{2_y} + P_{3_y}}{2} = \frac{k + sin(θ) - k \cdot cos(θ)}{2} \\
|
||||
e_{1_y} &= \frac{A_y + \frac{1}{2}P_{2_y}}{2} = \frac{\frac{k + sin(θ) - k \cdot cos(θ)}{2} + \frac{}{2}}{2} = \frac{2k + sin(θ) + k \cdot cos(θ)}{4} \\
|
||||
e_{2_y} &= \frac{A_y + \textit{mid}(P_4, P_3)}{2} = \frac{A_y + sin(θ) - \frac{k}{2} cos(θ)}{2} = \frac{k + 3sin(θ) 2k \cdot cos(θ)}{4} \\
|
||||
B_y &= \frac{e_{1_y} + e_{2_y}}{2} = \frac{3k + 4sin(θ) - 3k \cdot cos(θ)}{8}
|
||||
\end{aligned}
|
||||
\]
|
||||
|
||||
using the same scaling factor, because circular arcs are symmetrical, so our approximation will need to be symmetrical, too.
|
||||
Which now gives us two identities for B, because in addition to determining B through linear interpolation, we also know that B's _y_ coordinate is just _sin(θ/2)_: we started this exercise by saying we were going to approximate the circular arc using a Bézier curve that had its midpoint, which is point B, touching the unit circle at the arc's half-angle, by definition making B the point at (cos(θ/2), sin(θ/2)).
|
||||
|
||||
Starting with this information, we slowly maths our way to success, but I won't lie: the maths for this is pretty trig-heavy, and it's easy to get lost if you remember (or know!) some of the core trigonometric identities, so if you just want to see the final result just skip past the next section!
|
||||
This means we can equate the two identities we now have for B<sub>y</sub> and solve for _k_.
|
||||
|
||||
<div class="note">
|
||||
|
||||
## Let's do this thing.
|
||||
## Deriving _k_
|
||||
|
||||
Unlike for the quadratic case, we need some more information in order to compute <i>a</i> and <i>b</i>, since they're no longer dependent variables. First, we observe that the curve is symmetrical, so whatever values we end up finding for C<sub>1</sub> will apply to C<sub>2</sub> as well (rotated along its tangent), so we'll focus on finding the location of C<sub>1</sub> only. So here's where we do something that you might not expect: we're going to ignore for a moment, because we're going to have a much easier time if we just solve this problem with geometry first, then move to calculus to solve a much simpler problem.
|
||||
|
||||
If we look at the triangle that is formed between our starting point, or initial guess C<sub>1</sub> and our real C<sub>1</sub>, there's something funny going on: if we treat the line {start,guess} as our opposite side, the line {guess,real} as our adjacent side, with {start,real} our hypotenuse, then the angle for the corner hypotenuse/adjacent is half that of the arc we're covering. Try it: if you place the end point at a quarter circle (pi/2, or 90 degrees), the angle in our triangle is half a quarter (pi/4, or 45 degrees). With that knowledge, and a knowledge of what the length of any of our lines segments are (as a function), we can determine where our control points are, and thus have everything we need to find the error distance function. Of the three lines, the one we can easiest determine is {start,guess}, so let's find out what the guessed control point is. Again geometrically, because we have the benefit of an on-curve <i>t=0.5</i> value.
|
||||
|
||||
The distance from our guessed point to the start point is exactly the same as the projection distance we looked at earlier. Using <i>t=0.5</i> as our point "B" in the "A,B,C" projection, then we know the length of the line segment {C,A}, since it's d<sub>1</sub> = {A,B} + d<sub>2</sub> = {B,C}:
|
||||
Solving for _k_ is fairly straight forward, but it's a fair few steps, and if you just the immediate result: using a tool like [Wolfram Alpha](https://www.wolframalpha.com/) is definitely the way to go. That said, let's get going:
|
||||
|
||||
\[
|
||||
||{A,C}|| = d_2 + d_1 = d_2 + d_2 \cdot \textit{ratio}_3 \left(\frac{1}{2}\right) = d_2 + \frac{1}{3}d_2 = \frac{4}{3}d_2
|
||||
\begin{aligned}
|
||||
\frac{3c + 4sin(θ)) - 3k \cdot cos(θ)}{8} &= sin(\frac{θ}{2}) \\
|
||||
3k + 4sin(θ)) - 3k \cdot cos(θ) &= 8sin\left(\frac{θ}{2}\right) \\
|
||||
3k - 3k \cdot cos(θ) &= 8sin\left(\frac{θ}{2}\right) - 4sin(θ) \\
|
||||
3k (1 - cos(θ)) &= 4 \left ( 2sin\left(\frac{θ}{2} \right) - sin(θ) \right ) \\
|
||||
3k &= 4 \cdot \frac{2sin(\frac{θ}{2}) - sin(θ)}{1 - cos(θ)} \\
|
||||
k &= \frac{4}{3} \cdot \frac {2sin\left(\frac{θ}{2}\right) - sin(θ)}{1 - cos(θ)}
|
||||
\end{aligned}
|
||||
\]
|
||||
|
||||
So that just leaves us to find the distance from <i>t=0.5</i> to the baseline for an arbitrary angle φ, which is the distance from the centre of the circle to our <i>t=0.5</i> point, minus the distance from the centre to the line that runs from start point to end point. The first is the same as the point P we found for the quadratic curve:
|
||||
And finally, we can take further advantage of several trigonometric identities to _drastically_ simplify our formula for _k_:
|
||||
|
||||
\[
|
||||
P_x = cos(\frac{φ}{2}) ~, ~\ P_y = sin(\frac{φ}{2})
|
||||
\begin{aligned}
|
||||
k &= \frac{4}{3} \cdot \frac {2sin\left(\frac{θ}{2}\right) - sin(θ)}{1 - cos(θ)}\\
|
||||
k &= \frac{4}{3} \cdot \left ( \frac {2sin\left(\frac{θ}{2}\right)}{1 - cos(θ)} - \frac {sin(θ)}{1 - cos(θ)} \right )\\
|
||||
k &= \frac{4}{3} \cdot \left (csc\left(\frac{θ}{2}\right) - cot\left(\frac{θ}{2}\right) \right )\\
|
||||
k &= \frac{4}{3} \cdot tan\left ( \frac{θ}{4} \right )\\
|
||||
\end{aligned}
|
||||
\]
|
||||
|
||||
And the distance from the origin to the line start/end is another application of angles, since the triangle {origin,start,C} has known angles, and two known sides. We can find the length of the line {origin,C}, which lets us trivially compute the coordinate for C:
|
||||
|
||||
\[
|
||||
\begin{array}{l}
|
||||
l = cos(\frac{φ}{2}) ~, \\
|
||||
\left\{\begin{array}{l}
|
||||
C_x = l \cdot cos\left(\frac{φ}{2}\right) = cos^2\left(\frac{φ}{2}\right)~, \\
|
||||
C_y = l \cdot sin\left(\frac{φ}{2}\right) = cos(\frac{φ}{2}) \cdot sin\left(\frac{φ}{2}\right)~, \\
|
||||
\end{array}\right.
|
||||
\end{array}
|
||||
\]
|
||||
|
||||
With the coordinate C, and knowledge of coordinate B, we can determine coordinate A, and get a vector that is identical to the vector {start,guess}:
|
||||
|
||||
\[
|
||||
\left\{\begin{array}{l}
|
||||
B_x - C_x = cos\left(\frac{φ}{2}\right) - cos^2\left(\frac{φ}{2}\right) \\
|
||||
B_y - C_y = sin\left(\frac{φ}{2}\right) - cos(\frac{φ}{2}) \cdot sin\left(\frac{φ}{2}\right)
|
||||
= sin\left(\frac{φ}{2}\right) - \frac{sin(φ)}{2}
|
||||
\end{array}\right.
|
||||
\]
|
||||
|
||||
\[
|
||||
\left\{\begin{array}{l}
|
||||
\vec{v}_x = \{C,A\}_x = \frac{4}{3} \cdot (B_x - C_x) \\
|
||||
\vec{v}_y = \{C,A\}_y = \frac{4}{3} \cdot (B_y - C_y)
|
||||
\end{array}\right.
|
||||
\]
|
||||
|
||||
Which means we can now determine the distance {start,guessed}, which is the same as the distance {C,A}, and use that to determine the vertical distance from our start point to our C<sub>1</sub>:
|
||||
|
||||
\[
|
||||
\left\{\begin{array}{l}
|
||||
C_{1x} = 1 \\
|
||||
C_{1y} = \frac{d}{sin\left(\frac{φ}{2}\right)}
|
||||
= \frac{\sqrt{\vec{v}^2_x + \vec{v}^2_y}}{sin\left(\frac{φ}{2}\right)}
|
||||
= \frac{4}{3} tan \left( \frac{φ}{4} \right)
|
||||
\end{array}\right.
|
||||
\]
|
||||
|
||||
And after this tedious detour to find the coordinate for C<sub>1</sub>, we can find C<sub>2</sub> fairly simply, since it's lies at distance -C<sub>1y</sub> along the end point's tangent:
|
||||
|
||||
\[
|
||||
\begin{array}{l}
|
||||
E'_x = -sin(φ) ~,\\
|
||||
E'_y = cos(φ) ~, \\
|
||||
||E'|| = \sqrt{ (-sin(φ))^2 + cos^2(φ)} = 1 ~, \\
|
||||
\\
|
||||
\left\{\begin{array}{l}
|
||||
C_2x = E_x - C_{1y} \cdot \frac{E_x'}{||E'||}
|
||||
= cos(φ) + C_{1y} \cdot sin(φ)
|
||||
= cos(φ) + \frac{4}{3} tan \left( \frac{φ}{4} \right) \cdot sin(φ) \\
|
||||
C_2y = E_y - C_{1y} \cdot \frac{E_y'}{||E'||}
|
||||
= sin(φ) - C_{1y} \cdot cos(φ)
|
||||
= sin(φ) - \frac{4}{3} tan \left( \frac{φ}{4} \right) \cdot cos(φ)
|
||||
\end{array}\right.
|
||||
\end{array}
|
||||
\]
|
||||
|
||||
And that's it, we have all four points now for an approximation of an arbitrary circular arc with angle φ.
|
||||
And we're done.
|
||||
|
||||
</div>
|
||||
|
||||
So, to recap, given an angle φ, the new control coordinates are:
|
||||
So, the distance of our control points to the start/end points can be expressed as a number that we get from an almost trivial expression involving the circular arc's angle:
|
||||
|
||||
\[
|
||||
C_1 = \left [ \begin{matrix}
|
||||
1 \\
|
||||
f
|
||||
\end{matrix} \right ]~,~\textit{with}~f = \frac{4}{3} tan \left( \frac{φ}{4} \right)
|
||||
k = f(θ) = \frac{4}{3} tan\left (\frac{θ}{4} \right )
|
||||
\]
|
||||
|
||||
and
|
||||
Which means that for any circular arc with angle θ and radius _r_, our Bézier approximation based on three points of incidence is:
|
||||
|
||||
\[
|
||||
C_2 = \left [ \begin{matrix}
|
||||
cos(φ) + f \cdot sin(φ) \\
|
||||
sin(φ) - f \cdot cos(φ)
|
||||
\end{matrix} \right ]~,~\textit{with}~f = \frac{4}{3} tan \left( \frac{φ}{4} \right)
|
||||
\begin{aligned}
|
||||
\textit{start} &= (r,~0) \\
|
||||
\textit{control}_{~1} &= (r,~k) \\
|
||||
\textit{control}_{~2} &= r\cdot(cos(θ) + k \cdot sin(θ), sin(θ) - k \cdot cos(θ)) \\
|
||||
\textit{end} &= r \cdot (cos(θ),~sin(θ))
|
||||
\end{aligned}
|
||||
\]
|
||||
|
||||
And, because the "quarter curve" special case comes up so incredibly often, let's look at what these new control points mean for the curve coordinates of a quarter curve, by simply filling in φ = π/2:
|
||||
Which also gives us the commonly found value of 0.55228 for quarter circles, based on them having an angle of half π:
|
||||
|
||||
\[
|
||||
\begin{array}{l}
|
||||
S = (1, 0) ~, \
|
||||
C_1 = \left ( 1, 4 \frac{\sqrt{2}-1}{3} \right ) ~, \
|
||||
C_2 = \left ( 4 \frac{\sqrt{2}-1}{3} , 1 \right ) ~, \
|
||||
E = (0, 1)
|
||||
\end{array}
|
||||
f\left ( \frac{\pi}{2} \right ) = \frac{4}{3} \cdot tan\left(\frac{\pi}{8}\right) = \frac{4}{3}(\sqrt{2}-1)\approx 0.55228474983[...]
|
||||
\]
|
||||
|
||||
Which, in decimal values, rounded to six significant digits, is:
|
||||
And thus giving us the following Bézier coordinates for a quarter circle of radius _r_:
|
||||
|
||||
\[
|
||||
\begin{array}{l}
|
||||
S = (1, 0) ~, \
|
||||
C_1 = (1, 0.55228) ~, \
|
||||
C_2 = (0.55228 , 1) ~, \
|
||||
E = (0, 1)
|
||||
\end{array}
|
||||
\begin{aligned}
|
||||
\textit{start} &= (r,~0) \\
|
||||
\textit{control}_{~1} &= (r,~0.55228 \cdot r) \\
|
||||
\textit{control}_{~2} &= (0.55228 \cdot r,~r) \\
|
||||
\textit{end} &= (0,~r)
|
||||
\end{aligned}
|
||||
\]
|
||||
|
||||
Of course, this is for a circle with radius 1, so if you have a different radius circle, simply multiply the coordinate by the radius you need. And then finally, forming a full curve is now a simple a matter of mirroring these coordinates about the origin:
|
||||
<div class="note">
|
||||
|
||||
<graphics-element title="Cubic Bézier circle approximation" width="340" height="300" src="./circle.js"></graphics-element>
|
||||
## So, how accurate is this?
|
||||
|
||||
Unlike for the quadratic curve, we can't use <i>t=0.5</i> as our reference point because by its very nature it's one of the three points that are actually guaranteed to be on the circular arc itself. Instead, we need a different <i>t</i> value that will give us the maximum deflection - there are two possible choices (as our curve is still strictly "overshoots" the circular arc, and it's symmetrical) but rather than trying to use calculus to find the perfect _t_ value—which we could! the maths is perfectly reasonable as long as we get to use computers—we can also just perform a binary search for the biggest deflection and not bother with all this maths stuff.
|
||||
|
||||
So let's do that instead: we can run a maximum deflection check that just runs through _t_ from 0 to 1 at some coarse interval, finds a _t_ value that has "the highest deflection of the bunch", then reruns the same check with a much smaller interval around that _t_ value, repeating as many times as necessary to get us an arbitrarily precise value of _t_:
|
||||
|
||||
```
|
||||
getMostWrongT(radius, bezier, start, end, epsilon=1e-15):
|
||||
if end-start < epsilon:
|
||||
return (start+end)/2
|
||||
worst_t = 0
|
||||
max = 0
|
||||
stepsize = (end-start)/10
|
||||
for t=start to end, using stepsize:
|
||||
p = bezier.get(t)
|
||||
diff = p.magnitude() - radius
|
||||
if diff > max:
|
||||
worst_t = t
|
||||
max = diff
|
||||
return getMostWrongT(radius, bezier, worst_t - stepsize, worst_t + stepsize)
|
||||
```
|
||||
|
||||
Plus, how often do you get to write a function with that name?
|
||||
|
||||
Using this code, we find that our _t_ values are approximately 0.211325 and 0.788675, so let's pick the lower of the two and see what the maximum deflection is across our domain of angles, with the original quadratic error show in green (rocketing off to infinity first, and then coming back down as we approach 2π)
|
||||
|
||||
<table><tbody><tr><td>
|
||||
<img src="images/chapter-assets/circles/image-20210417173811587.png" width="95%"/>
|
||||
</td><td>
|
||||
<img src="images/chapter-assets/circles/image-20210417174019035.png" width="95%"/>
|
||||
</td><td>
|
||||
<img src="images/chapter-assets/circles/image-20210417174100036.png" width="95%"/>
|
||||
</td></tr>
|
||||
<tr><td>
|
||||
error plotted for 0 ≤ φ ≤ 2π
|
||||
</td><td>
|
||||
error plotted for 0 ≤ φ ≤ π
|
||||
</td><td>
|
||||
error plotted for 0 ≤ φ ≤ ½π
|
||||
</td></tr>
|
||||
</tbody></table>
|
||||
|
||||
That last image is probably not quite clear enough: the cubic approximation of a quarter circle is so incredibly much better that we can't even really see it at the same scale of our quadratic curve. Let's scale the y-axis a little, and try that again:
|
||||
|
||||
<p style="text-align: center"><img src="images/chapter-assets/circles/image-20210417174215876.png" height="350px"></p>
|
||||
|
||||
Yeah... the error of a cubic approximation for a quarter circle turns out to be _two orders of magnitude_ better. At approximately 0.00027 (or: just shy of being 2.7 pixels off for a circle with a radius of 10,000 pixels) the increase in precision over quadratic curves is quite spectacular - certainly good enough that no one in their right mind should ever use quadratic curves.
|
||||
|
||||
</div>
|
||||
|
||||
So that's it, kappa is _4/3 · tan(θ/4)_ , we're done! ...or are we?
|
||||
|
||||
## Can we do better?
|
||||
|
||||
Technically: yes, we can. But I'm going to prefix this section with "we can, and we should investigate that possibility, but let me warn you up front that the result is _only_ better if we're going to hard-code the values". We're about to get into the weeds and the standard three-points-of-incidence value is so good already that for most applications, trying to do better won't make any sense at all.
|
||||
|
||||
So with that said: what we calculated above is an _upper bound_ for a best fit Bézier curve for a circular arc: anywhere we don't touch the circular arc in our approximation, we've "overshot" the arc. What if we dropped our value for _k_ just a little, so that the curve starts out as an over-estimation, but then crosses the circular arc, yielding an region of underestimation, and then crosses the circular arc again, with another region of overestimation. This might give us a lower overall error, so let's see what we can do.
|
||||
|
||||
First, let's express the total error (given circular arc angle θ, and some _k_) using standard calculus notation:
|
||||
|
||||
\[
|
||||
\textit{erf}~(θ, k) = \int_0^1{\left \| \sqrt{B_x(t,θ,k)^2 + B_y(t,θ,k)^2} - r \right \|dt}
|
||||
\]
|
||||
|
||||
This says that the error function for a given angle and value of _k_ is equal to the "infinite" sum of differences between our curve and the circular arc, as we run _t_ from 0 to 1, using an infinitely small step size. between subsequent _t_ values.
|
||||
|
||||
Now, since we want to find the minimal error, that means we want to know where along this function things go from "error is getting progressively less" to "error is increasing again", which means we want to know where its derivative is zero, which as mathematical expression looks like:
|
||||
|
||||
\[
|
||||
\left ( \int_0^1{\left \| \sqrt{B_x^2 + B_y^2} - r \right \|dt} \right ) \frac{d}{dt} = 0
|
||||
\]
|
||||
|
||||
And here we have the most direct application of the [Fundamental Theorem of Calculus](https://en.wikipedia.org/wiki/Fundamental_theorem_of_calculus): the derivative and integral are each other's inverse operations, so they cancel out, leaving us with our original function:
|
||||
|
||||
\[
|
||||
\left \| \sqrt{B_x^2 + B_y^2} - r \right \| = 0, ~~ t \in [0,1]
|
||||
\]
|
||||
|
||||
And now we just solve for that... oh wait. We've seen this before. In order to solve this, we'd end up needing to solve this:
|
||||
|
||||
\[
|
||||
B_x^2 + B_y^2 = r
|
||||
\]
|
||||
|
||||
And both of those terms on the left of the equal sign are 6<sup>th</sup> degree polynomials, which means—as we've covered in the section on arc lengths—[there is no symbolic solution for this equasion](https://en.wikipedia.org/wiki/Abel%E2%80%93Ruffini_theorem). Instead, we'll have to use a numerical approach to find the solutions here, so... to the computer!
|
||||
|
||||
<div class="note">
|
||||
|
||||
## Iterating on a solution
|
||||
|
||||
By which I really mean "to the binary search algorithm", because we're dealing with a reasonably well behaved function: depending on the value for _k_ , we're either going to end up with a Bézier curve that's on average "not at distance _r_ from the arc's center", "exactly distance _r_ from the arc's center", or "more than distance _r_ from the arc's center", so we can just binary search our way to the most accurate value for _c_ that gets us that middle case.
|
||||
|
||||
First our setup, where we determine our upper and lower bounds, before entering our binary search:
|
||||
|
||||
```
|
||||
findBest(radius, angle, points[]):
|
||||
lowerBound = 0
|
||||
upperBound = 4.0/3.0 * Math.tan(abs(angle) / 4)
|
||||
return binarySearch(radius, angle, points, lowerBound, upperBound)
|
||||
```
|
||||
|
||||
And then the binary search algorithm, which can be found in pretty much any CS textbook, as well as more online articles, tutorials, and blog posts than you can ever read in a life time:
|
||||
|
||||
```
|
||||
binarySearch(radius, angle, points[], lowerBound, upperBound, epsilon=1e-15):
|
||||
value = (upperBound + lowerBound)/2
|
||||
|
||||
if (upperBound - lowerBound < epsilon) return value
|
||||
|
||||
// recompute the control points, based on our current "value"
|
||||
d = (points[3].y < 0 ? -1 : 1) * value * radius
|
||||
points[1] = new Point(radius, d)
|
||||
points[2] = new Point(
|
||||
points[3].x + d * sin(angle)
|
||||
points[3].y - d * cos(angle)
|
||||
)
|
||||
|
||||
if radialError(radius, points) > 0:
|
||||
// our bezier curve is longer than we want it to be: reduce the upper bound
|
||||
return binarySearch(radius, angle, points, lowerBound, value)
|
||||
else:
|
||||
// our bezier curve is shorter than we want it to be: increase the lower bound
|
||||
return binarySearch(radius, angle, points, value, upperBound)
|
||||
```
|
||||
|
||||
Using the following `radialError` function, which samples the curve's approximation of the circular arc over several points (although the first and last point will never contribute anything, so we skip them):
|
||||
|
||||
```
|
||||
radialError(radius, points[]):
|
||||
err = 0
|
||||
steps = 5.0
|
||||
for (int i=1; i<steps; i++):
|
||||
Point p = getOnCurvePoint(points, i/steps)
|
||||
err += p.magnitude()/radius - 1
|
||||
return err
|
||||
```
|
||||
|
||||
In this, `getOnCurvePoint` is just the standard Bézier evaluation function, yielding a point. Treating that point as a vector, we can get its length to the origin using a `magnitude` call.
|
||||
|
||||
## Examining the result
|
||||
|
||||
Running the above code we can get a list of _k_ values associated with a list of angles θ from 0 to π, and we can use that to, for each angle, plot what the difference between the circular arc and the Bézier approximation looks like:
|
||||
|
||||

|
||||
|
||||
Here we see the difference between an arc and its Bézier approximation plotted as we run _t_ from 0 to 1. Just by looking at the plot we can tell that there is maximum deflection at _t_ = 0.5, so let's plot the maximum deflection "function", for angles from 0 to θ:
|
||||
|
||||
In fact, let's plot the maximum deflections for both approaches as a functions over θ:
|
||||
|
||||
<table><tbody><tr><td>
|
||||
<img src="images/chapter-assets/circles/image-20210418111929371.png" width="95%"/>
|
||||
</td><td>
|
||||
<img src="images/chapter-assets/circles/image-20210418112008676.png" width="95%"/>
|
||||
</td><td>
|
||||
<img src="images/chapter-assets/circles/image-20210418112038613.png" width="95%"/>
|
||||
</td></tr>
|
||||
<tr><td>
|
||||
max deflection using unit scale
|
||||
</td><td>
|
||||
max deflection at 10x scale
|
||||
</td><td>
|
||||
max deflection at 100x scale
|
||||
</td></tr>
|
||||
</tbody></table>
|
||||
|
||||
That doesn't actually appear to be all that much better, so let's look at some numbers, to see what the improvement actually is:
|
||||
|
||||
| angle | "improved" deflection | "upper bound" deflection | difference |
|
||||
| ----- | --------------------- | ------------------------ | --------------------- |
|
||||
| 1/8 π | 6.202833502388927E-8 | 6.657161222278773E-8 | 4.5432771988984655E-9 |
|
||||
| 1/4 π | 3.978021202111215E-6 | 4.246252911066506E-6 | 2.68231708955291E-7 |
|
||||
| 3/8 π | 4.547652269037972E-5 | 4.8397483513262785E-5 | 2.9209608228830675E-6 |
|
||||
| 1/2 π | 2.569196199214696E-4 | 2.7251652752280364E-4 | 1.559690760133403E-5 |
|
||||
| 5/8 π | 9.877526288810667E-4 | 0.0010444175859711802 | 5.666495709011343E-5 |
|
||||
| 3/4 π | 0.00298164978679627 | 0.0031455628414580605 | 1.6391305466179062E-4 |
|
||||
| 7/8 π | 0.0076323182807019885 | 0.008047777909948373 | 4.1545962924638413E-4 |
|
||||
| π | 0.017362185964043708 | 0.018349016519545902 | 9.86830555502194E-4 |
|
||||
|
||||
As we can see, the increase in precision is not particularly big: for a quarter circle (π/2) the traditional _k_ will be off by 2.75 pixels on a circle with radius 10,000 pixels, whereas this "better" fit will be off by 2.56 pixels. And while that's certainly an almost 10% improvement, it's also nowhere near enough of an improvement to make a discernible difference.
|
||||
|
||||
</div>
|
||||
|
||||
At this point it should be clear that while, yes, there are improvement to be had, they're essentially insignificant while also being _much_ more computationally expensive.
|
||||
|
||||
## TL;DR: just tell me which value I should be using
|
||||
|
||||
It depends on what we need to do. If we just want the best value for quarter circles, and we're going to hard code the value for _k_, then there is no reason to hard-code the constant `k=4/3*tan(pi/8)` when you can just as easily hard-code the constant as `k=0.551784777779014` instead.
|
||||
|
||||
**If you need "the" value for quarter circles, use 0.551785 instead of 0.55228**
|
||||
|
||||
However, for dynamic arc approximation, in code that tries to fit circular paths using Bézier paths instead, it should be fairly obvious that the simple function involving a tangent computation, two divisions, and one multiplication, is vastly more performant than running all the code we ended writing just to get a 25% lower error value, and most certainly worth preferring over getting the "more accurate" value.
|
||||
|
||||
**If you need to fit Béziers to circular arcs on the fly, use `4/3 * tan(θ/4)`**
|
||||
|
||||
However, always remember that if you're writing for humans, you can typically use the best of both worlds: as the user interacts with their curves, you should draw _their curves_ instead of drawing approximations of them. If they need to draw circles or circular arcs, draw those, and only approximate them with a Bézier curve when the data needs to be exported to a format that doesn't support those. Ideally with a preview mechanism that highlights where the errors will be, and how large they will be.
|
||||
|
||||
**If you're writing code for graphics design by humans, use circular arcs for circular arcs**
|
||||
|
||||
And that's it. We have pretty well exhausted this subject. There are different metrics we could use to find "different best _k_ values", like trying to match arc length (e.g. when we're optimizing for material cost), or minimizing the area between the circular arc and the Bézier curve (e.g. when we're optimizing for inking), or minimizing the rate of change of the Bézier's curvature (e.g. when we're optimizing for curve traversal) and they all yield values that are so similar that it's almost certainly not worth it. (For instance, for quarter circle approximations those values are 0.551777, 0.5533344, and 0.552184 respectively. Much like the 0.551785 we get from minimizing the maximum deflection, none of these values are significantly better enough to prefer them over the upper bound value).
|
Reference in New Issue
Block a user