-
- 
+ 
So let's look at that in action: the following graphic is interactive in that you can use your up and down arrow keys to increase or decrease the interpolation ratio, to see what happens. We start with three points, which gives us two lines. Linear interpolation over @@ -514,7 +514,6 @@ -
And that brings us to the complicated maths: calculus.
While it doesn't look like that's what we've just done, we actually just drew a quadratic curve, in steps, rather than in a single go. One @@ -543,13 +542,13 @@ function". An illustration: Let's say we have a function that maps some value, let's call it x, to some other value, using some kind of number manipulation:
-
+ 
The notation f(x) is the standard way to show that it's a function (by convention called f if we're only listing one) and its output changes based on one variable (in this case, x). Change x, and the output for f(x) changes.
So far, so good. Now, let's look at parametric functions, and how they cheat. Let's take the following two functions:
-
+ 
There's nothing really remarkable about them, they're just a sine and cosine function, but you'll notice the inputs have different names. If we change the value for a, we're not going to change the output value for f(b), since a isn't used in that @@ -557,8 +556,8 @@
@@ -568,7 +567,7 @@
fa(t) and fb(t) with what we usually mean with them for parametric curves, things might be a lot more
obvious:
- 
+ 
There we go. x/y coordinates, linked through some mystery value t.
So, parametric curves don't define a y coordinate in terms of an x coordinate, like normal functions do, but they instead @@ -585,7 +584,6 @@ -
Bézier curves are just one out of the many classes of parametric functions, and are characterised by using the same base function for all of the output values. In the example we saw above, the x and y values were generated by different functions (one uses a @@ -595,9 +593,9 @@
You may remember polynomials from high school. They're those sums that look like this:
@@ -611,8 +609,8 @@
@@ -622,9 +620,9 @@
@@ -638,8 +636,8 @@
@@ -650,9 +648,9 @@
@@ -662,11 +660,34 @@
We could naively implement the basis function as a mathematical construct, using the function as our guide, like this:
-function Bezier(n,t):
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+
I say we could, because we're not going to: the factorial function is incredibly expensive. And, as we can see from the above
explanation, we can actually create Pascal's triangle quite easily without it: just start at [1], then [1,1], then [1,2,1], then
@@ -677,7 +698,13 @@
We can generate this as a list of lists lightning fast, and then never have to compute the binomial terms because we have a lookup
table:
- lut = [ [1], // n=0
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+ 8
+
+
+ 9
+
+
+ 10
+
+
+ 11
+
+
+ 12
+
+
+ 13
+
+
+ 14
+
+
+ 15
+
+
+ 16
+
+
+ 17
+
+
+ 18
+
+
+
So what's going on here? First, we declare a lookup table with a size that's reasonably large enough to accommodate most lookups. Then,
we declare a function to get us the values we need, and we make sure that if an n/k pair is requested that isn't in the LUT yet,
we expand it first. Our basis function now looks like this:
- function Bezier(n,t):
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+
Perfect. Of course, we can optimize further. For most computer graphics purposes, we don't need arbitrary curves (although we will also
provide code for arbitrary curves in this primer); we need quadratic and cubic curves, and that means we can drastically simplify the
code:
- function Bezier(2,t):
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+ 8
+
+
+ 9
+
+
+ 10
+
+
+ 11
+
+
+ 12
+
+
+ 13
+
+
+
And now we know how to program the basis function. Excellent.
+
+
+
-
-
+
+
+
@@ -782,14 +939,14 @@ function Bezier(3,t): is about as straightforward as possible: just multiply each point with a value that changes its strength. These values are conventionally called "weights", and we can add them to our original Bézier function:
-
+ 
That looks complicated, but as it so happens, the "weights" are actually just the coordinate values we want our curve to have: for an nth order curve, w0 is our start coordinate, wn is our last coordinate, and everything in between is a controlling coordinate. Say we want a cubic curve that starts at (110,150), is controlled by (25,190) and (210,250) and ends at (210,30), we use this Bézier curve:
-
+ 
Which gives us the curve we saw at the top of the article:
What else can we do with Bézier curves? Quite a lot, actually. The rest of this article covers a multitude of possible operations and algorithms that we can apply, and the tasks they achieve. @@ -806,13 +962,42 @@ function Bezier(3,t):
Given that we already know how to implement basis function, adding in the control points is remarkably easy:
-function Bezier(n,t,w[]):
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+
And now for the optimized versions:
- function Bezier(2,t,w[]):
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+ 8
+
+
+ 9
+
+
+ 10
+
+
+ 11
+
+
+ 12
+
+
+ 13
+
+
+
And now we know how to program the weighted basis function.
Adding these ratio values to the regular Bézier curve function is fairly easy. Where the regular function is the following:
The function for rational Bézier curves has two more terms:
@@ -886,7 +1112,6 @@ function Bezier(3,t,w[]): -
You can think of the ratio values as each coordinate's "gravity": the higher the gravity, the closer to that coordinate the curve will want to be. You'll also notice that if you simply increase or decrease all the ratios by the same amount, nothing changes... much like @@ -896,7 +1121,13 @@ function Bezier(3,t,w[]):
Extending the code of the previous section to include ratios is almost trivial:
-function RationalBezier(2,t,w[],r[]):
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+ 8
+
+
+ 9
+
+
+ 10
+
+
+ 11
+
+
+ 12
+
+
+ 13
+
+
+ 14
+
+
+ 15
+
+
+ 16
+
+
+ 17
+
+
+ 18
+
+
+ 19
+
+
+ 20
+
+
+ 21
+
+
+ 22
+
+
+ 23
+
+
+ 24
+
+
+ 25
+
+
+ 26
+
+
+
And that's all we have to do.
+ 
The obvious start and end values here need to be a=1, b=0, so that the mixed value is 100% value 1, and 0% value 2, and
a=0, b=1, so that the mixed value is 0% value 1 and 100% value 2. Additionally, we don't want "a" and "b" to be independent:
@@ -947,7 +1258,7 @@ function RationalBezier(3,t,w[],r[]):
and end point, so we need to make sure we can never set "a" and "b" to some values that lead to a mix value that sums to more than 100%.
And that's easy:
+ 
With this we can guarantee that we never sum above 100%. By restricting a to values in the interval [0,1], we will always be
somewhere between our two values (inclusively), and we will always sum to a 100% mix.
@@ -1013,28 +1324,28 @@ function RationalBezier(3,t,w[],r[]):
coefficients matrix, and the actual coordinates as a matrix. Let's look at what this means for the cubic curve, using P... to
refer to coordinate values "in one or more dimensions":
+ 
Disregarding our actual coordinates for a moment, we have:
-
+ 
We can write this as a sum of four expressions:
-
+ 
And we can expand these expressions:
-
+ 
Furthermore, we can make all the 1 and 0 factors explicit:
-
+ 
And that, we can view as a series of four matrix operations:
-
+ 
If we compact this into a single matrix operation, we get:
-
+ 
This kind of polynomial basis representation is generally written with the bases in increasing order, which means we need to flip our
t matrix horizontally, and our big "mixing" matrix upside down:
+ 
And then finally, we can add in our original coordinates as a single third matrix:
-
+ 
We can perform the same trick for the quadratic curve, in which case we end up with:
-
+ 
If we plug in a t value, and then multiply the matrices, we will get exactly the same values as when we evaluate the original
polynomial function, or as when we evaluate the curve using progressive linear interpolation.
@@ -1097,32 +1408,101 @@ function RationalBezier(3,t,w[],r[]):
-
Let's just use the algorithm we just specified, and implement that:
-function drawCurve(points[], t):
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+ 8
+
+
+
And done, that's the algorithm implemented. Except usually you don't get the luxury of overloading the "+" operator, so let's also give
the code for when you need to work with x and y values:
- function drawCurve(points[], t):
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+ 8
+
+
+ 9
+
+
+ 10
+
+
+
So what does this do? This draws a point, if the passed list of points is only 1 point long. Otherwise it will create a new list of
points that sit at the t ratios (i.e. the "markers" outlined in the above algorithm), and then call the draw function for this
@@ -1162,14 +1542,16 @@ function RationalBezier(3,t,w[],r[]):
-
-
- Scripts are disabled. Showing fallback image.
- 
-
-
-
-
+
+
+
+ Scripts are disabled. Showing fallback image.
+
+
+
+
+
+
@@ -1180,21 +1562,79 @@ function RationalBezier(3,t,w[],r[]):
Let's just use the algorithm we just specified, and implement that:
-function flattenCurve(curve, segmentCount):
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+
And done, that's the algorithm implemented. That just leaves drawing the resulting "curve" as a sequence of lines:
- function drawFlattenedCurve(curve, segmentCount):
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+
We start with the first coordinate as reference point, and then just draw lines between each point and its next point.
We can implement curve splitting by bolting some extra logging onto the de Casteljau function:
-left=[]
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+ 8
+
+
+ 9
+
+
+ 10
+
+
+ 11
+
+
+ 12
+
+
+ 13
+
+
+ 14
+
+
+ 15
+
+
+ 16
+
+
+
After running this function for some value t, the left and right arrays will contain all the
coordinates for two new curves - one to the "left" of our t value, the other on the "right". These new curves will have the
@@ -1256,7 +1751,7 @@ function drawCurve(points[], t):
and
and
Excellent! Now we can form our new quadratic curve:
If we want the interval [z,1], we will be evaluating this instead:
So, our final second curve looks like:
and
and
nth order curve to an (n+1)th order curve is as follows (observing
that the start and end weights are the same as the start and end weights for the old curve):
- 
+ 
However, this rule also has as direct consequence that you cannot generally safely lower a curve from
nth order to (n-1)th order, because the control points cannot be "pulled apart" cleanly. We can
@@ -1536,27 +2031,27 @@ function drawCurve(points[], t):
some things can be done much more easily with matrices than with calculus functions, and this is one of those things. So... let's go!
We start by taking the standard Bézier function, and condensing it a little:
- 
+ 
Then, we apply one of those silly (actually, super useful) calculus tricks: since our t value is always between zero and one
(inclusive), we know that (1-t) plus t always sums to 1. As such, we can express any value as a sum of
t and 1-t:
- 
+ 
So, with that seemingly trivial observation, we rewrite that Bézier function by splitting it up into a sum of a (1-t) and
t component:
- 
+ 
So far so good. Now, to see why we did this, let's write out the (1-t) and t parts, and see what that gives us.
I promise, it's about to make sense. We start with (1-t):
@@ -1568,9 +2063,9 @@ function drawCurve(points[], t):
@@ -1584,8 +2079,8 @@ function drawCurve(points[], t):
Let's do this:
@@ -1593,13 +2088,13 @@ function drawCurve(points[], t):
And this is where we switch over from calculus to linear algebra, and matrices: we can now express this relation between Bézier(n,t) and
Bézier(n+1,t) as a very simple matrix multiplication:
- 
+ 
where the matrix M is an n+1 by n matrix, and looks like:
@@ -1617,9 +2112,9 @@ function drawCurve(points[], t):
The steps taken here are:
@@ -1668,7 +2163,7 @@ function drawCurve(points[], t):
First, let's look at the derivative rule for Bézier curves, which is:
Which is hard to work with, so let's expand that properly:
And that's the first part done: the two components inside the parentheses are actually regular, lower-order Bézier expressions:
And that's just a summation of lower order curves:
We can rewrite this as a normal summation, and we're done:
@@ -1857,9 +2352,9 @@ function drawCurve(points[], t):
@@ -1869,9 +2364,9 @@ function drawCurve(points[], t):
@@ -1886,16 +2381,16 @@ function drawCurve(points[], t):
Which is the "long" version of the following matrix transformation:
@@ -2019,7 +2514,6 @@ function drawCurve(points[], t):
-
However, if you've played with that graphic a bit, you might have noticed something odd. The normal seems to "suddenly twist around the
curve" between t=0.65 and t=0.75... Why is it doing that?
@@ -2070,14 +2564,46 @@ function drawCurve(points[], t):
We first assume we have a function for calculating the Frenet frame at a point, which we already discussed above, inn a way that it
yields a frame with properties:
- {
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+
Then, we can write a function that generates a sequence of RM frames in the following manner:
- generateRMFrames(steps) -> frames:
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+ 8
+
+
+ 9
+
+
+ 10
+
+
+ 11
+
+
+ 12
+
+
+ 13
+
+
+ 14
+
+
+ 15
+
+
+ 16
+
+
+ 17
+
+
+ 18
+
+
+ 19
+
+
+ 20
+
+
+ 21
+
+
+ 22
+
+
+ 23
+
+
+ 24
+
+
+ 25
+
+
+ 26
+
+
+ 27
+
+
+ 28
+
+
+ 29
+
+
+ 30
+
+
+ 31
+
+
+ 32
+
+
+ 33
+
+
+ 34
+
+
+ 35
+
+
+ 36
+
+
+
Ignoring comments, this is certainly more code than when we were just computing a single Frenet frame, but it's not a crazy amount more
code to get much better looking normals.
@@ -2131,7 +2767,6 @@ function drawCurve(points[], t):
-
That looks so much better!
For those reading along with the code: we don't even strictly speaking need a Frenet frame to start with: we could, for instance, treat
@@ -2178,7 +2813,6 @@ function drawCurve(points[], t):
-
And then we turn this into our solution for t using basic arithmetics:
@@ -2239,9 +2873,9 @@ function drawCurve(points[], t):
@@ -2256,17 +2890,17 @@ function drawCurve(points[], t):
And then, using these v values, we can find out what our a, b, and c should be:
@@ -2275,9 +2909,9 @@ function drawCurve(points[], t):
Easy-peasy. We can now almost trivially find the roots by plugging those values into the quadratic formula.
@@ -2298,9 +2932,9 @@ function drawCurve(points[], t):
@@ -2331,14 +2965,20 @@ function drawCurve(points[], t):
"complex" number space (denoted with ℂ). And, as it so happens, Cardano is
also attributed as the first mathematician in history to have made use of complex numbers in his calculations. For this very algorithm!
- // A helper function to filter for values in the [0,1] interval:
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+ 8
+
+
+ 9
+
+
+ 10
+
+
+ 11
+
+
+ 12
+
+
+ 13
+
+
+ 14
+
+
+ 15
+
+
+ 16
+
+
+ 17
+
+
+ 18
+
+
+ 19
+
+
+ 20
+
+
+ 21
+
+
+ 22
+
+
+ 23
+
+
+ 24
+
+
+ 25
+
+
+ 26
+
+
+ 27
+
+
+ 28
+
+
+ 29
+
+
+ 30
+
+
+ 31
+
+
+ 32
+
+
+ 33
+
+
+ 34
+
+
+ 35
+
+
+ 36
+
+
+ 37
+
+
+ 38
+
+
+ 39
+
+
+ 40
+
+
+ 41
+
+
+ 42
+
+
+ 43
+
+
+ 44
+
+
+ 45
+
+
+ 46
+
+
+ 47
+
+
+ 48
+
+
+ 49
+
+
+ 50
+
+
+ 51
+
+
+ 52
+
+
+ 53
+
+
+ 54
+
+
+ 55
+
+
+ 56
+
+
+ 57
+
+
+ 58
+
+
+ 59
+
+
+ 60
+
+
+ 61
+
+
+ 62
+
+
+ 63
+
+
+ 64
+
+
+ 65
+
+
+ 66
+
+
+ 67
+
+
+ 68
+
+
+ 69
+
+
+ 70
+
+
+ 71
+
+
+ 72
+
+
+ 73
+
+
+ 74
+
+
+ 75
+
+
+ 76
+
+
+ 77
+
+
+ 78
+
+
+ 79
+
+
+ 80
+
+
+ 81
+
+
@@ -2454,9 +3338,9 @@ function getCubicRoots(pa, pb, pc, pd) {
(The Wikipedia article has a decent animation for this process, so I will not add a graphic for that here)
@@ -2492,7 +3376,6 @@ function getCubicRoots(pa, pb, pc, pd) {
-
And for cubic curves, that means first and second derivatives, in red and purple respectively:
n-2 term function. For instance, if we have a cubic curve such as this:
- 
+ 
Then translating it so that the first coordinate lies on (0,0), moving all x coordinates by -120, and all y coordinates
by -160, gives us:
- 
+ 
If we then rotate the curve so that its end point lies on the x-axis, the coordinates (integer-rounded for illustrative purposes here)
become:
- 
+ 
If we drop all the zero-terms, this gives us:
- 
+ 
We can see that our original curve definition has been simplified considerably. The following graphics illustrate the result of aligning
our example curves to the x-axis, with the cubic case using the coordinates that were just used in the example formulae:
@@ -2598,7 +3481,6 @@ function getCubicRoots(pa, pb, pc, pd) {
-
@@ -2666,7 +3548,7 @@ function getCubicRoots(pa, pb, pc, pd) {
inflection, and we can find out where those happen relatively easily.
What we need to do is solve a simple equation:
- 
+ 
What we're saying here is that given the curvature function C(t), we want to know for which values of t this function is
zero, meaning there is no "curvature", which will be exactly at the point between our circle being on one side of the curve, and our
@@ -2674,9 +3556,9 @@ function getCubicRoots(pa, pb, pc, pd) {
@@ -2698,16 +3580,16 @@ function getCubicRoots(pa, pb, pc, pd) {
And of course the same functions for y:
@@ -2717,9 +3599,9 @@ function getCubicRoots(pa, pb, pc, pd) {
@@ -2734,9 +3616,9 @@ function getCubicRoots(pa, pb, pc, pd) {
@@ -2745,17 +3627,17 @@ function getCubicRoots(pa, pb, pc, pd) {
This is a plain quadratic curve, and we know how to solve C(t) = 0; we use the quadratic formula:
@@ -2806,7 +3688,6 @@ function getCubicRoots(pa, pb, pc, pd) {
-
This is a fairly funky image, so let's see what the various parts of it mean...
We see the three fixed points at (0,0), (0,1) and (1,1). The various regions and boundaries indicate what property the original curve will
@@ -2826,9 +3707,9 @@ function getCubicRoots(pa, pb, pc, pd) {
@@ -2839,9 +3720,9 @@ function getCubicRoots(pa, pb, pc, pd) {
@@ -2849,8 +3730,8 @@ function getCubicRoots(pa, pb, pc, pd) {
...on the top (blue) edge, the curve's start point touches the curve, forming a loop. This edge is described by the function:
@@ -2907,12 +3788,12 @@ function getCubicRoots(pa, pb, pc, pd) {
cx + dy), which means we can't do translation, since that requires we end up with some kind of (x + a, y + b). If we add a bogus
z coordinate that is always 1, then we can suddenly add arbitrary values. For example:
- 
+ 
Sweet! z stays 1, so we can effectively ignore it entirely, but we added some plain values to our x and y coordinates. So, if we
want to subtract p1.x and p1.y, we use:
- 
+ 
Running all our coordinates through this transformation gives a new set of coordinates, let's call those U, where the
first coordinate lies on (0,0), and the rest is still somewhat free. Our next job is to make sure point 2 ends up lying on the
@@ -2920,12 +3801,12 @@ function getCubicRoots(pa, pb, pc, pd) {
currently have. This is called shearing, and the typical x-shear matrix and its
transformation looks like this:
- 
+ 
So we want some shearing value that, when multiplied by y, yields -x, so our x coordinate becomes zero. That value is
simply -x/y, because *-x/y * y = -x*. Done:
- 
+ 
Now, running this on all our points generates a new set of coordinates, let's call those V, which now have point 1 on
(0,0) and point 2 on (0, some-value), and we wanted it at (0,1), so we need to
@@ -2933,7 +3814,7 @@ function getCubicRoots(pa, pb, pc, pd) {
point 3 to end up on (1,1), so we can also scale x to make sure its x-coordinate will be 1 after we run the transform. That means we'll be
x-scaling by 1/point3x, and y-scaling by point2y. This is really easy:
- 
+ 
Then, finally, this generates a new set of coordinates, let's call those W, of which point 1 lies on (0,0), point 2 lies on (0,1), and
point three lies on (1, ...) so all that's left is to make sure point 3 ends up at (1,1) - but we can't scale! Point 2 is already in the
@@ -2942,13 +3823,13 @@ function getCubicRoots(pa, pb, pc, pd) {
but y-shearing. Additionally, we don't actually want to end up at zero (which is what we did before) so we need to shear towards an
offset, in this case 1:
- 
+ 
And this generates our final set of four coordinates. Of these, we already know that points 1 through 3 are (0,0), (0,1) and (1,1), and
only the last coordinate is "free". In fact, given any four starting coordinates, the resulting "transformation mapped" coordinate will
be:
- 
+ 
Okay, well, that looks plain ridiculous, but: notice that every coordinate value is being offset by the initial translation, and also
notice that a lot of terms in that expression are repeated. Even though the maths looks crazy as a single expression, we can just
@@ -2958,13 +3839,13 @@ function getCubicRoots(pa, pb, pc, pd) {
First, let's just do that translation step as a "preprocessing" operation so we don't have to subtract the values all the time. What does
that leave?
- 
+ 
Suddenly things look a lot simpler: the mapped x is fairly straight forward to compute, and we see that the mapped y actually contains the
mapped x in its entirety, so we'll have that part already available when we need to evaluate it. In fact, let's pull out all those common
factors to see just how simple this is:
- 
+ 
That's kind of super-simple to write out in code, I think you'll agree. Coding math tends to be easier than the formulae initially make it
look!
@@ -3031,7 +3912,6 @@ function getCubicRoots(pa, pb, pc, pd) {
-
Now, if you look more closely at that right graphic, you'll notice something interesting: if we treat the red line as "the x axis", then
the point where the function crosses our line is really just a root for the cubic function x(t) through a shifted "x-axis"... and
@@ -3039,23 +3919,29 @@ function getCubicRoots(pa, pb, pc, pd) {
cubic case either: because of the kind of curve we're starting with, we know there is only root, simplifying the code we need!
First, let's look at the function for x(t):
- 
+ 
We can rewrite this to a plain polynomial form, by just fully writing out the expansion and then collecting the polynomial factors, as:
- 
+ 
Nothing special here: that's a standard cubic polynomial in "power" form (i.e. all the terms are ordered by their power of
t). So, given that a, b, c, d, and x(t) are all
known constants, we can trivially rewrite this (by moving the x(t) across the equal sign) as:
- 
+ 
You might be wondering "where did all the other 'minus x' for all the other values a, b, c, and d go?" and the answer there is that they
all cancel out, so the only one we actually need to subtract is the one at the end. Handy! So now we just solve this equation using
Cardano's algorithm, and we're left with some rather short code:
- // prepare our values for root finding:
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+ 8
+
+
+ 9
+
+
+ 10
+
+
+
So the procedure is fairly straight forward: pick an x, find the associted t value, evaluate our curve
for that t value, which gives us the curve's {x,y} coordinate, which means we know y for this
@@ -3090,9 +4008,9 @@ y = curve.get(t).y
fy(t), then the length of the curve, measured from start point to some point t = z, is computed using the
following seemingly straight forward (if a bit overwhelming) formula:
- 
+ 
or, more commonly written using Leibnitz notation as:
- 
+ 
This formula says that the length of a parametric curve is in fact equal to the area underneath a function that looks a
remarkable amount like Pythagoras' rule for computing the diagonal of a straight angled triangle. This sounds pretty simple, right? Sadly,
@@ -3118,7 +4036,7 @@ y = curve.get(t).y
+ 
In plain text: an integral function can always be treated as the sum of an (infinite) number of (infinitely thin) rectangular strips sitting "under" the function's plotted graph. To illustrate this idea, the following graph shows the integral for a sinusoid function. The @@ -3180,9 +4098,9 @@ y = curve.get(t).y
@@ -3196,15 +4114,15 @@ y = curve.get(t).y values, for any n, so if we want to approximate our integral with only two terms (which is a bit low, really) then these tables would tell us that for n=2 we must use the following values:
-
+ 
Which means that in order for us to approximate the integral, we must plug these values into the approximate function, which gives us:
@@ -3257,7 +4175,6 @@ y = curve.get(t).y
-
So, what does the function look like? This:
-
+ 
Which is really just a "short form" that glosses over the fact that we're dealing with functions of t, so let's expand that a
tiny bit:
+ 
And while that's a litte more verbose, it's still just as simple to work with as the first function: the curvature at some point on any (and this cannot be overstated: any) curve is a ratio between the first and second derivative cross product, and something that @@ -3356,13 +4273,42 @@ y = curve.get(t).y evaluating this function for any t value is just a matter of basic arithematics.
In fact, let's just implement it right now:
-function kappa(t, B):
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+
That was easy! (Well okay, that "not a number" value will need to be taken into account by downstream code, but that's a reality of
programming anwyay)
@@ -3385,14 +4331,13 @@ y = curve.get(t).y
-
One thing you may have noticed in this sketch is that sometimes the curvature looks fine, but seems to be pointing in the wrong direction, making it hard to line up the curves properly. A way around that, of course, is to show the curvature on both sides of the curve, so let's just do that. But let's take it one step further: we can also compute the associated "radius of curvature", which gives us the implicit circle that "fits" the curve's curvature at any point, using what is possibly the simplest bit of maths found in this entire primer:
-
+ 
So let's revisit the previous graphic with the curvature visualised on both sides of our curves, as well as showing the circle that "fits"
our curve at some point that we can control by using a slider:
@@ -3441,7 +4386,6 @@ y = curve.get(t).y
-
So, how do we "cut up" the arc length function at regular intervals, when we can't really work with it? We basically cheat: we run through
the curve using t values, determine the distance-for-this-t-value at each point we generate during the run, and
@@ -3462,7 +4406,6 @@ y = curve.get(t).y
-
Use the slider to increase or decrease the number of equidistant segments used to colour the curve.
However, are there better ways? One such way is discussed in "
-
@@ -3519,7 +4461,13 @@ y = curve.get(t).y
but sometimes you need more function signatures, because you might not want to call your function with eight distinct parameters. Maybe
you're using point structs for the line. Let's get coding:
Finding self-intersections is effectively the same procedure, except that we're starting with a single curve, so we need to turn that into
two separate curves first. This is trivially achieved by splitting at an inflection point, or if there are none, just splitting at
@@ -3736,7 +4735,7 @@ lli = function(line1, line2):
So, how can we compute
If we can figure out what the function And
So, if we know the start and end coordinates, and we know the t value, we know C, without having to calculate the
We start by observing that, given Working out the maths for this, we see the following two formulae for quadratic and cubic curves: And
Which now leaves us with some powerful tools: given thee points (start, end, and "some point on the curve"), as well as a
With And then reverse engineer the curve's control control points:
So: if we have a curve's start and end point, then for any
@@ -3817,7 +4816,6 @@ lli = function(line1, line2):
-
For cubic curves we need to do a little more work, but really only just a little. We're first going to assume that a decent curve through
the three points should approximate a circular arc, which first requires knowing how to fit a circle to three points. You may remember (if
@@ -3837,7 +4835,6 @@ lli = function(line1, line2):
-
With that covered, we now also know the tangent line to our point
@@ -3862,9 +4859,9 @@ lli = function(line1, line2):
@@ -3874,9 +4871,9 @@ lli = function(line1, line2):
The result of this approach looks as follows:
It is important to remember that even though we're using a circular arc to come up with decent That looks perfectly servicable!
Of course, we can take this one step further: we can't just "create" curves, we also have (almost!) all the tools available to "mold"
@@ -3936,13 +4931,42 @@ lli = function(line1, line2):
binary search. First, we do a coarse distance-check based on
After this runs, we know that And then the associated And we're done, because that's our new quadratic control point!
As before, cubic curves are a bit more work, because while it's easy to find our initial
That looks reasonable, close to the original point, but the further we drag our point, the less "useful" things become. Especially if we
drag our point across the baseline, rather than turning into a nice curve.
@@ -4054,7 +5076,6 @@ for (coordinate, index) in LUT:
-
The slide controls the "falloff distance" relative to where the original point on the curve is, so that as we drag our point around, it
interpolates with a bias towards "preserving And then we (trivially) rearrange the terms across multiple lines:
@@ -4142,40 +5163,40 @@ for (coordinate, index) in LUT:
With that arrangement, we can easily decompose this as a matrix multiplication: We can do the same for the cubic curve, of course. We know the base function for cubics: So we write out the expansion and rearrange: Which we can then decompose: And, of course, we can do this for quartic curves too (skipping the expansion step):
Implementing line-line intersections
lli8 = function(x1,y1,x2,y2,x3,y3,x4,y4):
+
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+ 8
+
+
+ 9
+
+
+ 10
+
+
+ 11
+
+
+ 12
+
+
+ 13
+
+
+ 14
+
+
+ 15
+
+
+ 16
+
+
+ 17
+ What about curve-line intersections?
@@ -3643,7 +4643,6 @@ lli = function(line1, line2):
-
C? We start with our observation that C always lies somewhere between the start and ends
points, so logically C will have a function that interpolates between those two coordinates:
+ 
u(t) looks like, we'll be done. Although we do need to remember that this
u(t) will have a different for depending on whether we're working with quadratic or cubic curves.
@@ -3745,9 +4744,9 @@ lli = function(line1, line2):
>
(with thanks to Boris Zbarsky) shows us the following two formulae:
+ 
+ 
A or even B coordinates. In fact, we can do the same for the ratio function: as another function of
@@ -3755,18 +4754,18 @@ lli = function(line1, line2):
pure function of t, too.
A, B, and C, the following always holds:
+ 
+ 
+ 
t value, we can contruct curves: we can compute C using the start and end points, and our
u(t) function, and once we have C, we can use our on-curve point (B) and the
ratio(t) function to find A:
+ 
A found, finding e1 and e2 for quadratic curves is a matter of running the linear
interpolation with t between start and A to yield e1, and between A and end to yield
@@ -3774,9 +4773,9 @@ lli = function(line1, line2):
distance ratio between e1 to B and B to e2 is the Bézier ratio (1-t):t,
we can reverse engineer v1 and v2:
+ 
+ 
t value we implicitly know all the ABC values, which (combined
with an educated guess on appropriate e1 and e2 coordinates for cubic curves) gives us the necessary information
@@ -3801,9 +4800,9 @@ lli = function(line1, line2):
B, because the tangent to any point on the circle is a line
through that point, perpendicular to the line from that point to the center. That just leaves marking appropriate points
@@ -3847,9 +4844,9 @@ lli = function(line1, line2):
e1 and e2 terms,
we're not trying to perfectly create a circular arc with a cubic curve (which is good, because we can't;
@@ -3908,7 +4904,6 @@ lli = function(line1, line2):
-
t values associated with the curve's "to draw" coordinates (using a lookup table, or LUT). This is pretty fast:
p = some point to project onto the curve
+
+
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+ LUT[i] is the coordinate on the curve in our LUT that is closest to the point we want
to project, so that's a pretty good initial guess as to what the best projection onto our curve is. To refine it, we note that LUT[i] is a
@@ -3993,9 +5017,9 @@ for (coordinate, index) in LUT:
initial B coordinate. We don't even need the latter: with our t value and "whever the cursor is" as target
B, we can compute the associated C:
+ 
A:
+ 
t value and ABC values, getting
those all-important e1 and e2 coordinates is going to pose a bit of a problem... in the section on curve
@@ -4028,7 +5051,6 @@ for (coordinate, index) in LUT:
-
t/e1/e2" closer to the original point, and bias
@@ -4121,17 +5142,17 @@ for (coordinate, index) in LUT:
@@ -4192,8 +5213,8 @@ for (coordinate, index) in LUT:
@@ -4220,8 +5241,8 @@ for (coordinate, index) in LUT:
To get these values, we first compute the general "distance along the polygon" matrix:
@@ -4232,8 +5253,8 @@ for (coordinate, index) in LUT:
@@ -4248,8 +5269,8 @@ for (coordinate, index) in LUT:
@@ -4259,8 +5280,8 @@ for (coordinate, index) in LUT:
@@ -4272,17 +5293,17 @@ for (coordinate, index) in LUT:
In which we can replace the rather cumbersome "squaring" operation with a more conventional matrix equivalent:
@@ -4297,25 +5318,25 @@ for (coordinate, index) in LUT:
Which, because of the first and last values in S, means:
Now we can properly write out the error function as matrix operations:
@@ -4325,9 +5346,9 @@ for (coordinate, index) in LUT:
@@ -4395,7 +5416,6 @@ for (coordinate, index) in LUT: -
You'll note there is a convenient "toggle" buttons that lets you toggle between equidistant t values, and distance ratio
along the polygon formed by the points. Arguably more interesting is that once you have points to abstract a curve, you also get
@@ -4426,7 +5446,6 @@ for (coordinate, index) in LUT:
-
Now, it may look like Catmull-Rom curves are very different from Bézier curves, because these curves can get very long indeed, but what looks like a single Catmull-Rom curve is actually a spline: a single @@ -4438,9 +5457,9 @@ for (coordinate, index) in LUT:
@@ -4455,7 +5474,13 @@ for (coordinate, index) in LUT: considerable more complex):
-tension = some value greater than 0, defaulting to 1
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+ 8
+
+
+ 9
+
+
+ 10
+
+
+ 11
+
+
+ 12
+
+
+ 13
+
+
+ 14
+
+
+ 15
+
+
+ 16
+
+
+ 17
+
+
+ 18
+
+
+ 19
+
+
+
Now, since a Catmull-Rom curve is a form of cubic Hermite spline, and as
cubic Bézier curves are also a form of cubic Hermite spline, we run into an interesting bit of maths programming: we can
@@ -4497,8 +5581,8 @@ for p = 1 to points.length-3 (inclusive):
@@ -4523,9 +5607,9 @@ for p = 1 to points.length-3 (inclusive):
@@ -4536,17 +5620,17 @@ for p = 1 to points.length-3 (inclusive):
Computing T is really more "arranging the numbers":
Thus:
@@ -4557,9 +5641,9 @@ for p = 1 to points.length-3 (inclusive):
@@ -4568,73 +5652,73 @@ for p = 1 to points.length-3 (inclusive):
Replace point/tangent vector with the expression for all-coordinates:
and merge the matrices:
This looks a lot like the Bézier matrix form, which as we saw in the chapter on Bézier curves, should look like this:
So, if we want to express a Catmull-Rom curve using a Bézier curve, we'll need to turn this Catmull-Rom bit:
Into something that looks like this:
And the way we do that is with a fairly straight forward bit of matrix rewriting. We start with the equality we need to ensure:
Then we remove the coordinate vector from both sides without affecting the equality:
Then we can "get rid of" the Bézier matrix on the right by left-multiply both with the inverse of the Bézier matrix:
@@ -4643,17 +5727,17 @@ for p = 1 to points.length-3 (inclusive):
And now we're basically done. We just multiply those two matrices and we know what V is:
We now have the final piece of our function puzzle. Let's run through each step.
@@ -4662,8 +5746,8 @@ for p = 1 to points.length-3 (inclusive):
@@ -4672,9 +5756,9 @@ for p = 1 to points.length-3 (inclusive):
@@ -4682,9 +5766,9 @@ for p = 1 to points.length-3 (inclusive):
@@ -4692,9 +5776,9 @@ for p = 1 to points.length-3 (inclusive):
@@ -4702,9 +5786,9 @@ for p = 1 to points.length-3 (inclusive):
@@ -4712,9 +5796,9 @@ for p = 1 to points.length-3 (inclusive):
And we're done: we finally know how to convert these two curves!
@@ -4726,9 +5810,9 @@ for p = 1 to points.length-3 (inclusive):
@@ -4737,17 +5821,17 @@ for p = 1 to points.length-3 (inclusive):
Or, if your API allows you to specify Catmull-Rom curves using plain coordinates:
@@ -4813,20 +5897,19 @@ for p = 1 to points.length-3 (inclusive):
-
Dragging the control points around only affects the curve segments that the control point belongs to, and moving an on-curve point leaves
the control points where they are, which is not the most useful for practical modelling purposes. So, let's add in the logic we need to
make things a little better. We'll start by linking up control points by ensuring that the "incoming" derivative at an on-curve point is
the same as it's "outgoing" derivative:
- 
+ 
We can effect this quite easily, because we know that the vector from a curve's last control point to its last on-curve point is equal to
the derivative vector. If we want to ensure that the first control point of the next curve matches that, all we have to do is mirror that
last control point through the last on-curve point. And mirroring any point A through any point B is really simple:
- 
+ 
So let's implement that and see what it gets us. The following two graphics show a quadratic and a cubic poly-Bézier curve again, but this
time moving the control points around moves others, too. However, you might see something unexpected going on for quadratic curves...
@@ -4859,7 +5942,6 @@ for p = 1 to points.length-3 (inclusive):
-
As you can see, quadratic curves are particularly ill-suited for poly-Bézier curves, as all the control points are effectively linked.
Move one of them, and you move all of them. Not only that, but if we move the on-curve points, it's possible to get a situation where a
@@ -4902,7 +5984,6 @@ for p = 1 to points.length-3 (inclusive):
-
Cubic curves are now better behaved when it comes to dragging control points around, but the quadratic poly-Bézier still has the problem
that moving one control points will move the control points and may ending up defining "the next" control point in a way that doesn't
@@ -4940,7 +6021,6 @@ for p = 1 to points.length-3 (inclusive):
-
Again, we see that cubic curves are now rather nice to work with, but quadratic curves have a new, very serious problem: we can move an
on-curve point in such a way that we can't compute what needs to "happen next". Move the top point down, below the left and right points,
@@ -4986,9 +6066,9 @@ for p = 1 to points.length-3 (inclusive):
@@ -4999,9 +6079,9 @@ for p = 1 to points.length-3 (inclusive):
@@ -5013,8 +6093,8 @@ for p = 1 to points.length-3 (inclusive):
@@ -5025,9 +6105,9 @@ for p = 1 to points.length-3 (inclusive):
@@ -5035,9 +6115,9 @@ for p = 1 to points.length-3 (inclusive):
@@ -5099,7 +6179,6 @@ for p = 1 to points.length-3 (inclusive):
-
-
You may notice that this may still lead to small 'jumps' in the sub-curves when moving the curve around. This is caused by the fact that
we're still performing a naive form of offsetting, moving the control points the same distance as the start and end points. If the curve
@@ -5166,7 +6244,6 @@ for p = 1 to points.length-3 (inclusive):
-
-
As you can see, things go horribly wrong quite quickly; even trying to approximate a quarter circle using a quadratic curve is a bad idea.
An eighth of a turns might look okay, but how okay is okay? Let's apply some maths and find out. What we're interested in is how far off
@@ -5230,33 +6306,33 @@ for p = 1 to points.length-3 (inclusive):
We start out with our start and end point, and for convenience we will place them on a unit circle (a circle around 0,0 with radius 1), at
some angle φ:
- 
+ 
What we want to find is the intersection of the tangents, so we want a point C such that:
- 
+ 
i.e. we want a point that lies on the vertical line through S (at some distance a from S) and also lies on the tangent line
through E (at some distance b from E). Solving this gives us:
- 
+ 
First we solve for b:
- 
+ 
which yields:
- 
+ 
which we can then substitute in the expression for a:
- 
+ 
A quick check shows that plugging these values for a and b into the expressions for Cx and Cy give
the same x/y coordinates for both "a away from A" and "b away from B", so let's continue: now that we know the
coordinate values for C, we know where our on-curve point T for t=0.5 (or angle φ/2) is, because we can just evaluate the Bézier
polynomial, and we know where the circle arc's actual point P is for angle φ/2:
- 
+ 
We compute T, observing that if t=0.5, the polynomial values (1-t)², 2(1-t)t, and t² are 0.25, 0.5, and 0.25 respectively:
- 
+ 
Which, worked out for the x and y components, gives:
- 
+ 
And the distance between these two is the standard Euclidean distance:
- 
+ 
So, what does this distance function look like when we plot it for a number of ranges for the angle φ, such as a half circle, quarter
circle and eighth circle?
@@ -5297,7 +6373,7 @@ for p = 1 to points.length-3 (inclusive):
In fact, let's flip the function around, so that if we plug in the precision error, labelled ε, we get back the maximum angle for that
precision:
- 
+ 
And frankly, things are starting to look a bit ridiculous at this point, we're doing way more maths than we've ever done, but thankfully
this is as far as we need the maths to take us: If we plug in the precisions 0.1, 0.01, 0.001 and 0.0001 we get the radians values 1.748,
@@ -5344,7 +6420,6 @@ for p = 1 to points.length-3 (inclusive):
-
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
@@ -5403,8 +6478,8 @@ for p = 1 to points.length-3 (inclusive):
@@ -5415,16 +6490,16 @@ for p = 1 to points.length-3 (inclusive):
where "a" is some scaling factor we'll need to find the expression for, and:
@@ -5460,9 +6535,9 @@ for p = 1 to points.length-3 (inclusive):
@@ -5472,7 +6547,7 @@ for p = 1 to points.length-3 (inclusive):
@@ -5494,16 +6569,16 @@ for p = 1 to points.length-3 (inclusive):
@@ -5512,9 +6587,9 @@ for p = 1 to points.length-3 (inclusive):
@@ -5523,9 +6598,9 @@ for p = 1 to points.length-3 (inclusive):
And that's it, we have all four points now for an approximation of an arbitrary circular arc with angle φ.
@@ -5534,16 +6609,16 @@ for p = 1 to points.length-3 (inclusive):
So, to recap, given an angle φ, the new control coordinates are:
and
@@ -5553,16 +6628,16 @@ for p = 1 to points.length-3 (inclusive):
Which, in decimal values, rounded to six significant digits, is:
@@ -5657,7 +6732,6 @@ for p = 1 to points.length-3 (inclusive):
-
With that in place, all that's left now is to "restart" the procedure by treating the found arc's end point as the new to-be-determined
arc's starting point, and using points further down the curve. We keep trying this until the found end point is for t=1, at which
@@ -5672,7 +6746,6 @@ for p = 1 to points.length-3 (inclusive):
-
So... what is this good for? Obviously, if you're working with technologies that can't do curves, but can do lines and circles, then the
answer is pretty straightforward, but what else? There are some reasons why you might need this technique: using circular arcs means you
@@ -5715,7 +6788,6 @@ for p = 1 to points.length-3 (inclusive):
-
The important part to notice here is that we are not doing the same thing with B-Splines that we do for poly-Béziers or
Catmull-Rom curves: both of the latter simply define new sections as literally "new sections based on new points", so a 12 point cubic
@@ -5739,7 +6811,6 @@ for p = 1 to points.length-3 (inclusive):
-
In order to make this interpolation of curves work, the maths is necessarily more complex than the maths for Bézier curves, so let's have
a look at how things work.
@@ -5751,7 +6822,7 @@ for p = 1 to points.length-3 (inclusive):
>, we can compute a point on the curve for some value t in the interval [0,1] (where 0 is the start of the curve, and 1 the
end, just like for Bézier curves), by evaluating the following function:
- 
+ 
Which, honestly, doesn't tell us all that much. All we can see is that a point on a B-Spline curve is defined as "a mix of all the control
points, weighted somehow", where the weighting is achieved through the N(...) function, subscripted with an obvious parameter
@@ -5766,14 +6837,14 @@ for p = 1 to points.length-3 (inclusive):
k subscript to the N() function applies to.
Then the N() function itself. What does it look like?
- 
+ 
So this is where we see the interpolation: N(t) for an (i,k) pair (that is, for a step in the above summation, on a specific
knot interval) is a mix between N(t) for (i,k-1) and N(t) for (i+1,k-1), so we see that this is a recursive
iteration where i goes up, and k goes down, so it seem reasonable to expect that this recursion has to stop at
some point; obviously, it does, and specifically it does so for the following i/k values:
- 
+ 
And this function finally has a straight up evaluation: if a t value lies within a knot-specific interval once we reach a
k=1 value, it "counts", otherwise it doesn't. We did cheat a little, though, because for all these values we need to scale
@@ -5790,12 +6861,12 @@ for p = 1 to points.length-3 (inclusive):
Carl de Boor — came to a mathematically pleasing solution: to compute a point
P(t), we can compute this point by evaluating d(t) on a curve section between knots i and i+1:
- 
+ 
This is another recursive function, with k values decreasing from the curve order to 1, and the value α (alpha) defined
by:
- 
+ 
That looks complicated, but it's not. Computing alpha is just a fraction involving known, plain numbers. And, once we have our alpha
value, we also have (1-alpha) because it's a trivial subtraction. Computing the d() function is thus mostly a
@@ -5803,7 +6874,7 @@ for p = 1 to points.length-3 (inclusive):
recursion might see computationally expensive, the total algorithm is cheap, as each step only involves very simple maths.
Of course, the recursion does need a stop condition:
- 
+ 
So, we actually see two stopping conditions: either i becomes 0, in which case d() is zero, or
k becomes zero, in which case we get the same "either 1 or 0" that we saw in the N() function above.
@@ -5813,7 +6884,7 @@ for p = 1 to points.length-3 (inclusive):
Casteljau's algorithm. For instance, if we write out d() for i=3 and k=3, we get the following
recursion diagram:
- 
+ 
That is, we compute d(3,3) as a mixture of d(2,3) and d(2,2), where those two are themselves a
mixture of d(1,3) and d(1,2), and d(1,2) and d(1,1), respectively, which are
@@ -5839,23 +6910,78 @@ for p = 1 to points.length-3 (inclusive):
[knots[degree], knots[knots.length - 1 - degree]. Then, we find the section number s that this mapped
t value lies on:
- for(s=domain[0]; s < domain[1]; s++) {
- if(knots[s] <= t && t <= knots[s+1]) break;
-}
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+
after running this code, s is the index for the section the point will lie on. We then run the algorithm mentioned on the MU
page (updated to use this description's variable names):
- let v = copy of control points
-for(let L = 1; L <= order; L++) {
- for(let i=s; i > s + L - order; i--) {
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+ 8
+
+
+ 9
+
+
+ 10
+
+
+
(A nice bit of behaviour in this code is that we work the interpolation "backwards", starting at i=s at each level of the
interpolation, and we stop when i = s - order + level, so we always end up with a value for i such that those
@@ -5910,7 +7036,6 @@ for(let L = 1; L <= order; L++) {
-
This is an important point: the intervals that the knot vector defines are relative intervals, so it doesn't matter if every
interval is size 1, or size 100 - the relative differences between the intervals is what shapes any particular curve.
@@ -5935,7 +7060,6 @@ for(let L = 1; L <= order; L++) {
-
Open-Uniform B-Splines
By combining knot interval collapsing at the start and end of the curve, with uniform knots in between, we can overcome the problem of the
@@ -5956,7 +7080,6 @@ for(let L = 1; L <= order; L++) {
-
Non-uniform B-Splines
This is essentially the "free form" version of a B-Spline, and also the least interesting to look at, as without any specific reason to
@@ -5979,7 +7102,6 @@ for(let L = 1; L <= order; L++) {
-
Of course this brings us to the final topic that any text on B-Splines must touch on before calling it a day: the
NURBS, or Non-Uniform Rational B-Spline (NURBS is not a plural,
diff --git a/docs/ja-JP/index.html b/docs/ja-JP/index.html
index cf759e48..fd72865b 100644
--- a/docs/ja-JP/index.html
+++ b/docs/ja-JP/index.html
@@ -33,7 +33,7 @@
-
+
@@ -469,7 +469,7 @@
例えば、2点間の距離がわかっているとして、一方の点から距離の20%だけ離れた（すなわち、もう一方の点から80%離れた）新しい点を求めたい場合、次のようにとても簡単に計算できます。
- 
+ 
では、実際に見てみましょう。下の図はインタラクティブになっています。上下キーで補間の比率が増減しますので、どうなるか確かめてみましょう。最初に3点があり、それを結んで2本の直線が引かれています。この直線の上でそれぞれ線形補間を行うと、2つの点が得られます。この2点の間でさらに線形補間を行うと、1つの点を得ることができます。そして、あらゆる比率に対して同様に点を求め、それをすべて集めると、このようにベジエ曲線ができるのです。
@@ -481,7 +481,6 @@
-
また、これが複雑な方の数学につながっていきます。微積分です。
いま上で行ったものとは似つかないように思えますが、実はあれは2次曲線を描いていたのです。ただし一発で描きあげるのではなく、手順を追って描いていきました。ベジエ曲線は多項式関数で表現できますが、その一方で、とても単純な補間の補間の補間の……というふうにも説明できます。これがベジエ曲線のおもしろいところです。これはまた、ベジェ曲線は「本当の数学」で見る（関数を調べたり微分を調べたり、あらゆる方法で）ことも可能ですし、「機械的」な組み立て操作として見る（例えば、ベジエ曲線は組み立てに使う点の間からは決してはみ出ないということがわかります）ことも可能だということを意味しています。
@@ -498,27 +497,27 @@
ベジエ曲線は「パラメトリック」関数の一種です。数学的に言えば、パラメトリック関数というのはインチキです。というのも、「関数」はきっちり定義された用語であり、いくつかの入力を1つの出力に対応させる写像を表すものだからです。いくつかの数値を入れると、1つの数値が出てきます。入れる数値が変わっても、出てくる数値はやはり1つだけです。パラメトリック関数はインチキです。基本的には「じゃあわかった、値を複数個出したいから、関数を複数個使うことにするよ」ということです。例として、ある値xに何らかの操作を行い、別の値へと写す関数があるとします。
-
+
f(x)という記法は、これが関数（1つしかない場合は慣習的にfと呼びます）であり、その出力が1つの変数（この場合はxです）に応じて変化する、ということを示す標準的な方法です。xを変化させると、f(x)の出力が変化します。
ここまでは順調です。では、パラメトリック関数について、これがどうインチキなのかを見てみましょう。以下の2つの関数を考えます。
-
+
注目すべき箇所は特に何もありません。ただの正弦関数と余弦関数です。ただし、入力が別々の名前になっていることに気づくでしょう。仮にaの値を変えたとしても、f(b)の出力の値は変わらないはずです。なぜなら、こちらの関数にはaは使われていないからです。パラメトリック関数は、これを変えてしまうのでインチキなのです。パラメトリック関数においては、どの関数も変数を共有しています。例えば、
複数の関数がありますが、変数は1つだけです。tの値を変えた場合、fa(t)とfb(t)の両方の出力が変わります。これがどのように役に立つのか、疑問に思うかもしれません。しかし、実際には答えは至ってシンプルです。fa(t)とfb(t)のラベルを、パラメトリック曲線の表示によく使われているもので置き換えてやれば、ぐっとはっきりするかと思います。
-
+
きました。x/y座標です。謎の値tを通して繫がっています。
というわけで、普通の関数ではy座標をx座標によって定義しますが、パラメトリック曲線ではそうではなく、座標の値を「制御」変数と結びつけます。tの値を変化させるたびに2つの値が変化するので、これをグラフ上の座標
@@ -532,16 +531,15 @@
-
ベジエ曲線はパラメトリック関数の一種であり、どの次元に対しても同じ基底関数を使うという点で特徴づけられます。先ほどの例では、xの値とyの値とで異なる関数（正弦関数と余弦関数）を使っていましたが、ベジエ曲線ではxとyの両方で「二項係数多項式」を使います。では、二項係数多項式とは何でしょう？
高校で習った、こんな形の多項式を思い出すかもしれません。
@@ -552,7 +550,7 @@
基底関数の実装方法
上で説明した関数を使えば、数学的な組み立て方で、基底関数をナイーブに実装することもできます。
- function Bezier(n,t):
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+
「こともできる」と書いたのは、この方法では実装しない方が良いからです。階乗はとてつもなく重い計算なのです。また、先ほどの説明からわかるように、実際は階乗を使わなくても、かなり簡単にパスカルの三角形を作ることができます。[1]から始めて[1,1]、[1,2,1]、[1,3,3,1]、……としていくだけです。下の段は上の段よりも1つ要素が増え、各段の先頭と末尾は1になります。中間の数はどれも、左右斜め上にある両要素の和になります。
このパスカルの三角形は、「リストのリスト」として瞬時に生成できます。そして、これをルックアップテーブルとして利用すれば、二項係数を計算する必要はまったくなくなります。
- lut = [ [1], // n=0
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+ 8
+
+
+ 9
+
+
+ 10
+
+
+ 11
+
+
+ 12
+
+
+ 13
+
+
+ 14
+
+
+ 15
+
+
+ 16
+
+
+ 17
+
+
+ 18
+
+
+
これはどのように動くのでしょう？最初に、十分に大きなサイズのルックアップテーブルを宣言します。次に、求めたい値を得るための関数を定義します。この関数は、求めたい値のn/kのペアがテーブル中にまだ存在しない場合、先にテーブルを拡張するようになっています。さて、これで基底関数は次のようになりました。
- function Bezier(n,t):
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+
完璧です。もちろん、さらなる最適化を施すこともできます。コンピュータグラフィクス用途ではたいてい、任意の次数の曲線が必要になるわけではありません。2次と3次の曲線だけが必要であれば、以下のようにコードを劇的に単純化することができます（実際、この入門では任意の次数までは扱いませんので、これに似たようなコードが出てきます）。
- function Bezier(2,t):
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+ 8
+
+
+ 9
+
+
+ 10
+
+
+ 11
+
+
+ 12
+
+
+ 13
+
+
+
これで基底関数をプログラムする方法がわかりました。すばらしい。
+
+
+
-
-
+
+
+
@@ -700,11 +857,11 @@ function Bezier(3,t):
曲線を変更したい場合は、各点の重みを変える（実質的には補間率を変える）必要があります。これはとても単純で、寄与の大きさを変えるための値を、各点にただ掛ければいいのです。この値は「重み」と呼ばれていますが、これを元のベジエ関数に組み込めば、次のようになります。
-
+ 
複雑そうに見えますが、運がいいことに「重み」というのは実はただの座標値です。というのはn次の曲線の場合、w0が始点の座標、wnが終点の座標となり、その間はどれも制御点の座標になります。例えば、始点が(120,160)、制御点が(35,200)と(220,260)、終点が(220,40)となる3次ベジエ曲線は、次のようになります。
-
+
この式からは、記事の冒頭に出てきた曲線が得られます。
基底関数の実装方法はすでに知っていますし、これに制御点を組み込むのは非常に簡単です。
-function Bezier(n,t,w[]):
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+
そして、最適化を行ったバージョンは以下のようになります。
- function Bezier(2,t,w[]):
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+ 8
+
+
+ 9
+
+
+ 10
+
+
+ 11
+
+
+ 12
+
+
+ 13
+
+
+
これで、重みつき基底関数をプログラムする方法がわかりました。
Adding these ratio values to the regular Bézier curve function is fairly easy. Where the regular function is the following:
The function for rational Bézier curves has two more terms:
@@ -794,7 +1021,6 @@ function Bezier(3,t,w[]): -
You can think of the ratio values as each coordinate's "gravity": the higher the gravity, the closer to that coordinate the curve will want to be. You'll also notice that if you simply increase or decrease all the ratios by the same amount, nothing changes... much like @@ -804,7 +1030,13 @@ function Bezier(3,t,w[]):
Extending the code of the previous section to include ratios is almost trivial:
-function RationalBezier(2,t,w[],r[]):
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+ 8
+
+
+ 9
+
+
+ 10
+
+
+ 11
+
+
+ 12
+
+
+ 13
+
+
+ 14
+
+
+ 15
+
+
+ 16
+
+
+ 17
+
+
+ 18
+
+
+ 19
+
+
+ 20
+
+
+ 21
+
+
+ 22
+
+
+ 23
+
+
+ 24
+
+
+ 25
+
+
+ 26
+
+
+
And that's all we have to do.
このことは、曲線の「始点」から曲線の「終点」までどうやって動かすか、ということにすべて関係しています。2つの値を混ぜ合わせて1つの値をつくる場合、一般の式は次のようになります。
-
+ 
明らかに、始点ではa=1, b=0とする必要があります。こうすれば、値1が100%、値2が0%で混ぜ合わさるからです。また、終点ではa=0, b=1とする必要があります。こうすれば、値1が0%、値2が100%で混ぜ合わさります。これに加えて、aとbを独立にしておきたくはありません。独立になっている場合、何でも好きな値にすることできますが、こうすると例えば「値1が100%かつ値2が100%」のようなことが可能になってしまいます。これはこれで原則としてはかまいませんが、ベジエ曲線の場合は混ぜ合わさった値が常に始点と終点の間になってほしいのです。というわけで、混ぜ合わせの和が100%を決して超えないように、aとbの値を設定する必要があります。これは次のようにすれば簡単です。
+ 
こうすれば、和が100%を超えることはないと保証できます。aの値を区間[0,1]に制限してしまえば、混ぜ合わさった値は常に2つの値の間のどこか（両端を含む）になり、また和は常に100%になります。
ベジエ曲線は、行列演算の形でも表現することができます。ベジエ曲線の式を多項式基底と係数行列で表し、実際の座標も行列として表現するのです。これがどういうことを意味しているのか、3次ベジエ曲線について見てみましょう。
-
+
実際の座標を一旦無視すると、次のようになります。
-
+
これは、4つの項の和になっています。
-
+
それぞれの項を展開します。
-
+
その上で、係数の0や1もすべて明示的に書けば、このようになります。
-
+
さらに、これは4つの行列演算の和として見ることができます。
-
+
これを1つの行列演算にまとめると、以下のようになります。
-
+
多項式基底をこのような形で表現する場合、通常はその基底を昇冪の順に並べます。したがって、tの行列を左右反転させ、大きな「混合」行列は上下に反転させる必要があります。
+
そして最後に、もともとあった座標を3番目の行列として付け加えます。
-
+
2次ベジエ曲線の場合も同様に変形することができ、最終的には以下のようになります。
-
+
tの値を代入して行列の乗算を行えば、もともとの多項式関数から計算したときの値や、線形補間によって順次求めたときの値と、まったく同じものが得られます。
@@ -969,31 +1281,100 @@ function RationalBezier(3,t,w[],r[]):
-
いま説明したアルゴリズムを実装すると、以下のようになります。
-function drawCurve(points[], t):
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+ 8
+
+
+
これで実装完了です。ただし、+演算子のオーバーロードなどという贅沢品はたいてい無いでしょうから、xやyの値を直接扱う場合のコードも示しておきます。
- function drawCurve(points[], t):
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+ 8
+
+
+ 9
+
+
+ 10
+
+
+
さて、これは何をしているのでしょう？関数に渡す点のリストが長さ1であれば、点を1つ描きます。それ以外であれば、比率tの位置の点（すなわち、さきほどの説明に出てきた「印」）のリストを作り、そしてこの新しいリストを引数にして関数を呼び出します。
@@ -1020,14 +1401,16 @@ function RationalBezier(3,t,w[],r[]):
-
-
- Scripts are disabled. Showing fallback image.
-
-
-
-
-
+
+
+
+ Scripts are disabled. Showing fallback image.
+
+
+
+
+
+
@@ -1036,21 +1419,79 @@ function RationalBezier(3,t,w[],r[]):
上でいま解説したアルゴリズムを使って実装するだけです。
-function flattenCurve(curve, segmentCount):
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+
これで完了です。実装できました。あとは、一連の直線で結果の「曲線」を描画するだけです。
- function drawFlattenedCurve(curve, segmentCount):
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+
先頭の座標を参照点にしてスタートし、あとはそれぞれの点からその次の点へと、直線を引いていくだけです。
ド・カステリョの関数の中に記録する処理を追加すれば、曲線の分割が実装できます。
-left=[]
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+ 8
+
+
+ 9
+
+
+ 10
+
+
+ 11
+
+
+ 12
+
+
+ 13
+
+
+ 14
+
+
+ 15
+
+
+ 16
+
+
+
ある値tに対してこの関数を実行すると、leftとrightに新しい2曲線の座標が入ります。一方はtの「左」側、もう一方は「右」側の曲線です。この2曲線は元の曲線と同じ次数になり、また元の曲線とぴったり重なります。
@@ -1106,7 +1602,7 @@ function drawCurve(points[], t):
ならびに
ならびに
いいですね！これで、新しい2次ベジエ曲線が得られます。
区間[z,1]を求めたい場合は、かわりに次のような計算になります。
先ほどと同じ手法を使い、[なにか · M]を[M · なにか]に変えます。
よって、後半部分の曲線は結局のところ以下のようになります。
および
および
nth order curve to an (n+1)th order curve is as follows (observing
that the start and end weights are the same as the start and end weights for the old curve):
-
+
However, this rule also has as direct consequence that you cannot generally safely lower a curve from
nth order to (n-1)th order, because the control points cannot be "pulled apart" cleanly. We can
@@ -1372,27 +1868,27 @@ function drawCurve(points[], t):
some things can be done much more easily with matrices than with calculus functions, and this is one of those things. So... let's go!
We start by taking the standard Bézier function, and condensing it a little:
-
+
Then, we apply one of those silly (actually, super useful) calculus tricks: since our t value is always between zero and one
(inclusive), we know that (1-t) plus t always sums to 1. As such, we can express any value as a sum of
t and 1-t:
-
+
So, with that seemingly trivial observation, we rewrite that Bézier function by splitting it up into a sum of a (1-t) and
t component:
-
+
So far so good. Now, to see why we did this, let's write out the (1-t) and t parts, and see what that gives us.
I promise, it's about to make sense. We start with (1-t):
@@ -1404,9 +1900,9 @@ function drawCurve(points[], t):
@@ -1420,8 +1916,8 @@ function drawCurve(points[], t):
Let's do this:
@@ -1429,13 +1925,13 @@ function drawCurve(points[], t):
And this is where we switch over from calculus to linear algebra, and matrices: we can now express this relation between Bézier(n,t) and
Bézier(n+1,t) as a very simple matrix multiplication:
-
+
where the matrix M is an n+1 by n matrix, and looks like:
@@ -1453,9 +1949,9 @@ function drawCurve(points[], t):
The steps taken here are:
@@ -1504,7 +2000,7 @@ function drawCurve(points[], t):
First, let's look at the derivative rule for Bézier curves, which is:
Which is hard to work with, so let's expand that properly:
And that's the first part done: the two components inside the parentheses are actually regular, lower-order Bézier expressions:
And that's just a summation of lower order curves:
We can rewrite this as a normal summation, and we're done:
@@ -1693,9 +2189,9 @@ function drawCurve(points[], t):
@@ -1705,9 +2201,9 @@ function drawCurve(points[], t):
@@ -1722,16 +2218,16 @@ function drawCurve(points[], t):
Which is the "long" version of the following matrix transformation:
@@ -1855,7 +2351,6 @@ function drawCurve(points[], t):
-
However, if you've played with that graphic a bit, you might have noticed something odd. The normal seems to "suddenly twist around the
curve" between t=0.65 and t=0.75... Why is it doing that?
@@ -1906,14 +2401,46 @@ function drawCurve(points[], t):
We first assume we have a function for calculating the Frenet frame at a point, which we already discussed above, inn a way that it
yields a frame with properties:
- {
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+
Then, we can write a function that generates a sequence of RM frames in the following manner:
- generateRMFrames(steps) -> frames:
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+ 8
+
+
+ 9
+
+
+ 10
+
+
+ 11
+
+
+ 12
+
+
+ 13
+
+
+ 14
+
+
+ 15
+
+
+ 16
+
+
+ 17
+
+
+ 18
+
+
+ 19
+
+
+ 20
+
+
+ 21
+
+
+ 22
+
+
+ 23
+
+
+ 24
+
+
+ 25
+
+
+ 26
+
+
+ 27
+
+
+ 28
+
+
+ 29
+
+
+ 30
+
+
+ 31
+
+
+ 32
+
+
+ 33
+
+
+ 34
+
+
+ 35
+
+
+ 36
+
+
+
Ignoring comments, this is certainly more code than when we were just computing a single Frenet frame, but it's not a crazy amount more
code to get much better looking normals.
@@ -1967,7 +2604,6 @@ function drawCurve(points[], t):
-
That looks so much better!
For those reading along with the code: we don't even strictly speaking need a Frenet frame to start with: we could, for instance, treat
@@ -2014,7 +2650,6 @@ function drawCurve(points[], t):
-
And then we turn this into our solution for t using basic arithmetics:
@@ -2075,9 +2710,9 @@ function drawCurve(points[], t):
@@ -2092,17 +2727,17 @@ function drawCurve(points[], t):
And then, using these v values, we can find out what our a, b, and c should be:
@@ -2111,9 +2746,9 @@ function drawCurve(points[], t):
Easy-peasy. We can now almost trivially find the roots by plugging those values into the quadratic formula.
@@ -2134,9 +2769,9 @@ function drawCurve(points[], t):
@@ -2167,14 +2802,20 @@ function drawCurve(points[], t):
"complex" number space (denoted with ℂ). And, as it so happens, Cardano is
also attributed as the first mathematician in history to have made use of complex numbers in his calculations. For this very algorithm!
- // A helper function to filter for values in the [0,1] interval:
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+ 8
+
+
+ 9
+
+
+ 10
+
+
+ 11
+
+
+ 12
+
+
+ 13
+
+
+ 14
+
+
+ 15
+
+
+ 16
+
+
+ 17
+
+
+ 18
+
+
+ 19
+
+
+ 20
+
+
+ 21
+
+
+ 22
+
+
+ 23
+
+
+ 24
+
+
+ 25
+
+
+ 26
+
+
+ 27
+
+
+ 28
+
+
+ 29
+
+
+ 30
+
+
+ 31
+
+
+ 32
+
+
+ 33
+
+
+ 34
+
+
+ 35
+
+
+ 36
+
+
+ 37
+
+
+ 38
+
+
+ 39
+
+
+ 40
+
+
+ 41
+
+
+ 42
+
+
+ 43
+
+
+ 44
+
+
+ 45
+
+
+ 46
+
+
+ 47
+
+
+ 48
+
+
+ 49
+
+
+ 50
+
+
+ 51
+
+
+ 52
+
+
+ 53
+
+
+ 54
+
+
+ 55
+
+
+ 56
+
+
+ 57
+
+
+ 58
+
+
+ 59
+
+
+ 60
+
+
+ 61
+
+
+ 62
+
+
+ 63
+
+
+ 64
+
+
+ 65
+
+
+ 66
+
+
+ 67
+
+
+ 68
+
+
+ 69
+
+
+ 70
+
+
+ 71
+
+
+ 72
+
+
+ 73
+
+
+ 74
+
+
+ 75
+
+
+ 76
+
+
+ 77
+
+
+ 78
+
+
+ 79
+
+
+ 80
+
+
+ 81
+
+
@@ -2290,9 +3175,9 @@ function getCubicRoots(pa, pb, pc, pd) {
(The Wikipedia article has a decent animation for this process, so I will not add a graphic for that here)
@@ -2328,7 +3213,6 @@ function getCubicRoots(pa, pb, pc, pd) {
-
And for cubic curves, that means first and second derivatives, in red and purple respectively:
n-2 term function. For instance, if we have a cubic curve such as this:
-
+
Then translating it so that the first coordinate lies on (0,0), moving all x coordinates by -120, and all y coordinates
by -160, gives us:
-
+
If we then rotate the curve so that its end point lies on the x-axis, the coordinates (integer-rounded for illustrative purposes here)
become:
-
+
If we drop all the zero-terms, this gives us:
-
+
We can see that our original curve definition has been simplified considerably. The following graphics illustrate the result of aligning
our example curves to the x-axis, with the cubic case using the coordinates that were just used in the example formulae:
@@ -2434,7 +3318,6 @@ function getCubicRoots(pa, pb, pc, pd) {
-
@@ -2502,7 +3385,7 @@ function getCubicRoots(pa, pb, pc, pd) {
inflection, and we can find out where those happen relatively easily.
What we need to do is solve a simple equation:
-
+
What we're saying here is that given the curvature function C(t), we want to know for which values of t this function is
zero, meaning there is no "curvature", which will be exactly at the point between our circle being on one side of the curve, and our
@@ -2510,9 +3393,9 @@ function getCubicRoots(pa, pb, pc, pd) {
@@ -2534,16 +3417,16 @@ function getCubicRoots(pa, pb, pc, pd) {
And of course the same functions for y:
@@ -2553,9 +3436,9 @@ function getCubicRoots(pa, pb, pc, pd) {
@@ -2570,9 +3453,9 @@ function getCubicRoots(pa, pb, pc, pd) {
@@ -2581,17 +3464,17 @@ function getCubicRoots(pa, pb, pc, pd) {
This is a plain quadratic curve, and we know how to solve C(t) = 0; we use the quadratic formula:
@@ -2642,7 +3525,6 @@ function getCubicRoots(pa, pb, pc, pd) {
-
This is a fairly funky image, so let's see what the various parts of it mean...
We see the three fixed points at (0,0), (0,1) and (1,1). The various regions and boundaries indicate what property the original curve will
@@ -2662,9 +3544,9 @@ function getCubicRoots(pa, pb, pc, pd) {
@@ -2675,9 +3557,9 @@ function getCubicRoots(pa, pb, pc, pd) {
@@ -2685,8 +3567,8 @@ function getCubicRoots(pa, pb, pc, pd) {
...on the top (blue) edge, the curve's start point touches the curve, forming a loop. This edge is described by the function:
@@ -2743,12 +3625,12 @@ function getCubicRoots(pa, pb, pc, pd) {
cx + dy), which means we can't do translation, since that requires we end up with some kind of (x + a, y + b). If we add a bogus
z coordinate that is always 1, then we can suddenly add arbitrary values. For example:
-
+
Sweet! z stays 1, so we can effectively ignore it entirely, but we added some plain values to our x and y coordinates. So, if we
want to subtract p1.x and p1.y, we use:
-
+
Running all our coordinates through this transformation gives a new set of coordinates, let's call those U, where the
first coordinate lies on (0,0), and the rest is still somewhat free. Our next job is to make sure point 2 ends up lying on the
@@ -2756,12 +3638,12 @@ function getCubicRoots(pa, pb, pc, pd) {
currently have. This is called shearing, and the typical x-shear matrix and its
transformation looks like this:
-
+
So we want some shearing value that, when multiplied by y, yields -x, so our x coordinate becomes zero. That value is
simply -x/y, because *-x/y * y = -x*. Done:
-
+
Now, running this on all our points generates a new set of coordinates, let's call those V, which now have point 1 on
(0,0) and point 2 on (0, some-value), and we wanted it at (0,1), so we need to
@@ -2769,7 +3651,7 @@ function getCubicRoots(pa, pb, pc, pd) {
point 3 to end up on (1,1), so we can also scale x to make sure its x-coordinate will be 1 after we run the transform. That means we'll be
x-scaling by 1/point3x, and y-scaling by point2y. This is really easy:
-
+
Then, finally, this generates a new set of coordinates, let's call those W, of which point 1 lies on (0,0), point 2 lies on (0,1), and
point three lies on (1, ...) so all that's left is to make sure point 3 ends up at (1,1) - but we can't scale! Point 2 is already in the
@@ -2778,13 +3660,13 @@ function getCubicRoots(pa, pb, pc, pd) {
but y-shearing. Additionally, we don't actually want to end up at zero (which is what we did before) so we need to shear towards an
offset, in this case 1:
-
+
And this generates our final set of four coordinates. Of these, we already know that points 1 through 3 are (0,0), (0,1) and (1,1), and
only the last coordinate is "free". In fact, given any four starting coordinates, the resulting "transformation mapped" coordinate will
be:
-
+
Okay, well, that looks plain ridiculous, but: notice that every coordinate value is being offset by the initial translation, and also
notice that a lot of terms in that expression are repeated. Even though the maths looks crazy as a single expression, we can just
@@ -2794,13 +3676,13 @@ function getCubicRoots(pa, pb, pc, pd) {
First, let's just do that translation step as a "preprocessing" operation so we don't have to subtract the values all the time. What does
that leave?
-
+
Suddenly things look a lot simpler: the mapped x is fairly straight forward to compute, and we see that the mapped y actually contains the
mapped x in its entirety, so we'll have that part already available when we need to evaluate it. In fact, let's pull out all those common
factors to see just how simple this is:
-
+
That's kind of super-simple to write out in code, I think you'll agree. Coding math tends to be easier than the formulae initially make it
look!
@@ -2867,7 +3749,6 @@ function getCubicRoots(pa, pb, pc, pd) {
-
Now, if you look more closely at that right graphic, you'll notice something interesting: if we treat the red line as "the x axis", then
the point where the function crosses our line is really just a root for the cubic function x(t) through a shifted "x-axis"... and
@@ -2875,23 +3756,29 @@ function getCubicRoots(pa, pb, pc, pd) {
cubic case either: because of the kind of curve we're starting with, we know there is only root, simplifying the code we need!
First, let's look at the function for x(t):
-
+
We can rewrite this to a plain polynomial form, by just fully writing out the expansion and then collecting the polynomial factors, as:
-
+
Nothing special here: that's a standard cubic polynomial in "power" form (i.e. all the terms are ordered by their power of
t). So, given that a, b, c, d, and x(t) are all
known constants, we can trivially rewrite this (by moving the x(t) across the equal sign) as:
-
+
You might be wondering "where did all the other 'minus x' for all the other values a, b, c, and d go?" and the answer there is that they
all cancel out, so the only one we actually need to subtract is the one at the end. Handy! So now we just solve this equation using
Cardano's algorithm, and we're left with some rather short code:
- // prepare our values for root finding:
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+ 8
+
+
+ 9
+
+
+ 10
+
+
+
So the procedure is fairly straight forward: pick an x, find the associted t value, evaluate our curve
for that t value, which gives us the curve's {x,y} coordinate, which means we know y for this
@@ -2926,9 +3845,9 @@ y = curve.get(t).y
fy(t), then the length of the curve, measured from start point to some point t = z, is computed using the
following seemingly straight forward (if a bit overwhelming) formula:
-
+
or, more commonly written using Leibnitz notation as:
-
+
This formula says that the length of a parametric curve is in fact equal to the area underneath a function that looks a
remarkable amount like Pythagoras' rule for computing the diagonal of a straight angled triangle. This sounds pretty simple, right? Sadly,
@@ -2954,7 +3873,7 @@ y = curve.get(t).y
In plain text: an integral function can always be treated as the sum of an (infinite) number of (infinitely thin) rectangular strips sitting "under" the function's plotted graph. To illustrate this idea, the following graph shows the integral for a sinusoid function. The @@ -3016,9 +3935,9 @@ y = curve.get(t).y
@@ -3032,15 +3951,15 @@ y = curve.get(t).y values, for any n, so if we want to approximate our integral with only two terms (which is a bit low, really) then these tables would tell us that for n=2 we must use the following values:
-
+
Which means that in order for us to approximate the integral, we must plug these values into the approximate function, which gives us:
@@ -3093,7 +4012,6 @@ y = curve.get(t).y
-
So, what does the function look like? This:
-
+
Which is really just a "short form" that glosses over the fact that we're dealing with functions of t, so let's expand that a
tiny bit:
+
And while that's a litte more verbose, it's still just as simple to work with as the first function: the curvature at some point on any (and this cannot be overstated: any) curve is a ratio between the first and second derivative cross product, and something that @@ -3192,13 +4110,42 @@ y = curve.get(t).y evaluating this function for any t value is just a matter of basic arithematics.
In fact, let's just implement it right now:
-function kappa(t, B):
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+
That was easy! (Well okay, that "not a number" value will need to be taken into account by downstream code, but that's a reality of
programming anwyay)
@@ -3221,14 +4168,13 @@ y = curve.get(t).y
One thing you may have noticed in this sketch is that sometimes the curvature looks fine, but seems to be pointing in the wrong direction, making it hard to line up the curves properly. A way around that, of course, is to show the curvature on both sides of the curve, so let's just do that. But let's take it one step further: we can also compute the associated "radius of curvature", which gives us the implicit circle that "fits" the curve's curvature at any point, using what is possibly the simplest bit of maths found in this entire primer:
-
+
So let's revisit the previous graphic with the curvature visualised on both sides of our curves, as well as showing the circle that "fits"
our curve at some point that we can control by using a slider:
@@ -3277,7 +4223,6 @@ y = curve.get(t).y
-
So, how do we "cut up" the arc length function at regular intervals, when we can't really work with it? We basically cheat: we run through
the curve using t values, determine the distance-for-this-t-value at each point we generate during the run, and
@@ -3298,7 +4243,6 @@ y = curve.get(t).y
-
Use the slider to increase or decrease the number of equidistant segments used to colour the curve.
However, are there better ways? One such way is discussed in "
-
@@ -3355,7 +4298,13 @@ y = curve.get(t).y
but sometimes you need more function signatures, because you might not want to call your function with eight distinct parameters. Maybe
you're using point structs for the line. Let's get coding:
Finding self-intersections is effectively the same procedure, except that we're starting with a single curve, so we need to turn that into
two separate curves first. This is trivially achieved by splitting at an inflection point, or if there are none, just splitting at
@@ -3572,7 +4572,7 @@ lli = function(line1, line2):
So, how can we compute
If we can figure out what the function And
So, if we know the start and end coordinates, and we know the t value, we know C, without having to calculate the
We start by observing that, given Working out the maths for this, we see the following two formulae for quadratic and cubic curves: And
Which now leaves us with some powerful tools: given thee points (start, end, and "some point on the curve"), as well as a
With And then reverse engineer the curve's control control points:
So: if we have a curve's start and end point, then for any
@@ -3653,7 +4653,6 @@ lli = function(line1, line2):
-
For cubic curves we need to do a little more work, but really only just a little. We're first going to assume that a decent curve through
the three points should approximate a circular arc, which first requires knowing how to fit a circle to three points. You may remember (if
@@ -3673,7 +4672,6 @@ lli = function(line1, line2):
-
With that covered, we now also know the tangent line to our point
@@ -3698,9 +4696,9 @@ lli = function(line1, line2):
@@ -3710,9 +4708,9 @@ lli = function(line1, line2):
The result of this approach looks as follows:
It is important to remember that even though we're using a circular arc to come up with decent That looks perfectly servicable!
Of course, we can take this one step further: we can't just "create" curves, we also have (almost!) all the tools available to "mold"
@@ -3772,13 +4768,42 @@ lli = function(line1, line2):
binary search. First, we do a coarse distance-check based on
After this runs, we know that And then the associated And we're done, because that's our new quadratic control point!
As before, cubic curves are a bit more work, because while it's easy to find our initial
That looks reasonable, close to the original point, but the further we drag our point, the less "useful" things become. Especially if we
drag our point across the baseline, rather than turning into a nice curve.
@@ -3890,7 +4913,6 @@ for (coordinate, index) in LUT:
-
The slide controls the "falloff distance" relative to where the original point on the curve is, so that as we drag our point around, it
interpolates with a bias towards "preserving And then we (trivially) rearrange the terms across multiple lines:
@@ -3978,40 +5000,40 @@ for (coordinate, index) in LUT:
With that arrangement, we can easily decompose this as a matrix multiplication: We can do the same for the cubic curve, of course. We know the base function for cubics: So we write out the expansion and rearrange: Which we can then decompose: And, of course, we can do this for quartic curves too (skipping the expansion step):
Implementing line-line intersections
lli8 = function(x1,y1,x2,y2,x3,y3,x4,y4):
+
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+ 8
+
+
+ 9
+
+
+ 10
+
+
+ 11
+
+
+ 12
+
+
+ 13
+
+
+ 14
+
+
+ 15
+
+
+ 16
+
+
+ 17
+ What about curve-line intersections?
@@ -3479,7 +4480,6 @@ lli = function(line1, line2):
-
C? We start with our observation that C always lies somewhere between the start and ends
points, so logically C will have a function that interpolates between those two coordinates:
+
u(t) looks like, we'll be done. Although we do need to remember that this
u(t) will have a different for depending on whether we're working with quadratic or cubic curves.
@@ -3581,9 +4581,9 @@ lli = function(line1, line2):
>
(with thanks to Boris Zbarsky) shows us the following two formulae:
+
+
A or even B coordinates. In fact, we can do the same for the ratio function: as another function of
@@ -3591,18 +4591,18 @@ lli = function(line1, line2):
pure function of t, too.
A, B, and C, the following always holds:
+
+
+
t value, we can contruct curves: we can compute C using the start and end points, and our
u(t) function, and once we have C, we can use our on-curve point (B) and the
ratio(t) function to find A:
+
A found, finding e1 and e2 for quadratic curves is a matter of running the linear
interpolation with t between start and A to yield e1, and between A and end to yield
@@ -3610,9 +4610,9 @@ lli = function(line1, line2):
distance ratio between e1 to B and B to e2 is the Bézier ratio (1-t):t,
we can reverse engineer v1 and v2:
+
+
t value we implicitly know all the ABC values, which (combined
with an educated guess on appropriate e1 and e2 coordinates for cubic curves) gives us the necessary information
@@ -3637,9 +4637,9 @@ lli = function(line1, line2):
B, because the tangent to any point on the circle is a line
through that point, perpendicular to the line from that point to the center. That just leaves marking appropriate points
@@ -3683,9 +4681,9 @@ lli = function(line1, line2):
e1 and e2 terms,
we're not trying to perfectly create a circular arc with a cubic curve (which is good, because we can't;
@@ -3744,7 +4741,6 @@ lli = function(line1, line2):
-
t values associated with the curve's "to draw" coordinates (using a lookup table, or LUT). This is pretty fast:
p = some point to project onto the curve
+
+
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+ LUT[i] is the coordinate on the curve in our LUT that is closest to the point we want
to project, so that's a pretty good initial guess as to what the best projection onto our curve is. To refine it, we note that LUT[i] is a
@@ -3829,9 +4854,9 @@ for (coordinate, index) in LUT:
initial B coordinate. We don't even need the latter: with our t value and "whever the cursor is" as target
B, we can compute the associated C:
+
A:
+
t value and ABC values, getting
those all-important e1 and e2 coordinates is going to pose a bit of a problem... in the section on curve
@@ -3864,7 +4888,6 @@ for (coordinate, index) in LUT:
-
t/e1/e2" closer to the original point, and bias
@@ -3957,17 +4979,17 @@ for (coordinate, index) in LUT:
@@ -4028,8 +5050,8 @@ for (coordinate, index) in LUT:
@@ -4056,8 +5078,8 @@ for (coordinate, index) in LUT:
To get these values, we first compute the general "distance along the polygon" matrix:
@@ -4068,8 +5090,8 @@ for (coordinate, index) in LUT:
@@ -4084,8 +5106,8 @@ for (coordinate, index) in LUT:
@@ -4095,8 +5117,8 @@ for (coordinate, index) in LUT:
@@ -4108,17 +5130,17 @@ for (coordinate, index) in LUT:
In which we can replace the rather cumbersome "squaring" operation with a more conventional matrix equivalent:
@@ -4133,25 +5155,25 @@ for (coordinate, index) in LUT:
Which, because of the first and last values in S, means:
Now we can properly write out the error function as matrix operations:
@@ -4161,9 +5183,9 @@ for (coordinate, index) in LUT:
@@ -4231,7 +5253,6 @@ for (coordinate, index) in LUT: -
You'll note there is a convenient "toggle" buttons that lets you toggle between equidistant t values, and distance ratio
along the polygon formed by the points. Arguably more interesting is that once you have points to abstract a curve, you also get
@@ -4262,7 +5283,6 @@ for (coordinate, index) in LUT:
-
Now, it may look like Catmull-Rom curves are very different from Bézier curves, because these curves can get very long indeed, but what looks like a single Catmull-Rom curve is actually a spline: a single @@ -4274,9 +5294,9 @@ for (coordinate, index) in LUT:
@@ -4291,7 +5311,13 @@ for (coordinate, index) in LUT: considerable more complex):
-tension = some value greater than 0, defaulting to 1
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+ 8
+
+
+ 9
+
+
+ 10
+
+
+ 11
+
+
+ 12
+
+
+ 13
+
+
+ 14
+
+
+ 15
+
+
+ 16
+
+
+ 17
+
+
+ 18
+
+
+ 19
+
+
+
Now, since a Catmull-Rom curve is a form of cubic Hermite spline, and as
cubic Bézier curves are also a form of cubic Hermite spline, we run into an interesting bit of maths programming: we can
@@ -4333,8 +5418,8 @@ for p = 1 to points.length-3 (inclusive):
@@ -4359,9 +5444,9 @@ for p = 1 to points.length-3 (inclusive):
@@ -4372,17 +5457,17 @@ for p = 1 to points.length-3 (inclusive):
Computing T is really more "arranging the numbers":
Thus:
@@ -4393,9 +5478,9 @@ for p = 1 to points.length-3 (inclusive):
@@ -4404,73 +5489,73 @@ for p = 1 to points.length-3 (inclusive):
Replace point/tangent vector with the expression for all-coordinates:
and merge the matrices:
This looks a lot like the Bézier matrix form, which as we saw in the chapter on Bézier curves, should look like this:
So, if we want to express a Catmull-Rom curve using a Bézier curve, we'll need to turn this Catmull-Rom bit:
Into something that looks like this:
And the way we do that is with a fairly straight forward bit of matrix rewriting. We start with the equality we need to ensure:
Then we remove the coordinate vector from both sides without affecting the equality:
Then we can "get rid of" the Bézier matrix on the right by left-multiply both with the inverse of the Bézier matrix:
@@ -4479,17 +5564,17 @@ for p = 1 to points.length-3 (inclusive):
And now we're basically done. We just multiply those two matrices and we know what V is:
We now have the final piece of our function puzzle. Let's run through each step.
@@ -4498,8 +5583,8 @@ for p = 1 to points.length-3 (inclusive):
@@ -4508,9 +5593,9 @@ for p = 1 to points.length-3 (inclusive):
@@ -4518,9 +5603,9 @@ for p = 1 to points.length-3 (inclusive):
@@ -4528,9 +5613,9 @@ for p = 1 to points.length-3 (inclusive):
@@ -4538,9 +5623,9 @@ for p = 1 to points.length-3 (inclusive):
@@ -4548,9 +5633,9 @@ for p = 1 to points.length-3 (inclusive):
And we're done: we finally know how to convert these two curves!
@@ -4562,9 +5647,9 @@ for p = 1 to points.length-3 (inclusive):
@@ -4573,17 +5658,17 @@ for p = 1 to points.length-3 (inclusive):
Or, if your API allows you to specify Catmull-Rom curves using plain coordinates:
@@ -4649,20 +5734,19 @@ for p = 1 to points.length-3 (inclusive):
-
Dragging the control points around only affects the curve segments that the control point belongs to, and moving an on-curve point leaves
the control points where they are, which is not the most useful for practical modelling purposes. So, let's add in the logic we need to
make things a little better. We'll start by linking up control points by ensuring that the "incoming" derivative at an on-curve point is
the same as it's "outgoing" derivative:
-
+
We can effect this quite easily, because we know that the vector from a curve's last control point to its last on-curve point is equal to
the derivative vector. If we want to ensure that the first control point of the next curve matches that, all we have to do is mirror that
last control point through the last on-curve point. And mirroring any point A through any point B is really simple:
-
+
So let's implement that and see what it gets us. The following two graphics show a quadratic and a cubic poly-Bézier curve again, but this
time moving the control points around moves others, too. However, you might see something unexpected going on for quadratic curves...
@@ -4695,7 +5779,6 @@ for p = 1 to points.length-3 (inclusive):
-
As you can see, quadratic curves are particularly ill-suited for poly-Bézier curves, as all the control points are effectively linked.
Move one of them, and you move all of them. Not only that, but if we move the on-curve points, it's possible to get a situation where a
@@ -4738,7 +5821,6 @@ for p = 1 to points.length-3 (inclusive):
-
Cubic curves are now better behaved when it comes to dragging control points around, but the quadratic poly-Bézier still has the problem
that moving one control points will move the control points and may ending up defining "the next" control point in a way that doesn't
@@ -4776,7 +5858,6 @@ for p = 1 to points.length-3 (inclusive):
-
Again, we see that cubic curves are now rather nice to work with, but quadratic curves have a new, very serious problem: we can move an
on-curve point in such a way that we can't compute what needs to "happen next". Move the top point down, below the left and right points,
@@ -4822,9 +5903,9 @@ for p = 1 to points.length-3 (inclusive):
@@ -4835,9 +5916,9 @@ for p = 1 to points.length-3 (inclusive):
@@ -4849,8 +5930,8 @@ for p = 1 to points.length-3 (inclusive):
@@ -4861,9 +5942,9 @@ for p = 1 to points.length-3 (inclusive):
@@ -4871,9 +5952,9 @@ for p = 1 to points.length-3 (inclusive):
@@ -4935,7 +6016,6 @@ for p = 1 to points.length-3 (inclusive):
-
-
You may notice that this may still lead to small 'jumps' in the sub-curves when moving the curve around. This is caused by the fact that
we're still performing a naive form of offsetting, moving the control points the same distance as the start and end points. If the curve
@@ -5002,7 +6081,6 @@ for p = 1 to points.length-3 (inclusive):
-
-
As you can see, things go horribly wrong quite quickly; even trying to approximate a quarter circle using a quadratic curve is a bad idea.
An eighth of a turns might look okay, but how okay is okay? Let's apply some maths and find out. What we're interested in is how far off
@@ -5066,33 +6143,33 @@ for p = 1 to points.length-3 (inclusive):
We start out with our start and end point, and for convenience we will place them on a unit circle (a circle around 0,0 with radius 1), at
some angle φ:
-
+
What we want to find is the intersection of the tangents, so we want a point C such that:
-
+
i.e. we want a point that lies on the vertical line through S (at some distance a from S) and also lies on the tangent line
through E (at some distance b from E). Solving this gives us:
-
+
First we solve for b:
-
+
which yields:
-
+
which we can then substitute in the expression for a:
-
+
A quick check shows that plugging these values for a and b into the expressions for Cx and Cy give
the same x/y coordinates for both "a away from A" and "b away from B", so let's continue: now that we know the
coordinate values for C, we know where our on-curve point T for t=0.5 (or angle φ/2) is, because we can just evaluate the Bézier
polynomial, and we know where the circle arc's actual point P is for angle φ/2:
-
+
We compute T, observing that if t=0.5, the polynomial values (1-t)², 2(1-t)t, and t² are 0.25, 0.5, and 0.25 respectively:
-
+
Which, worked out for the x and y components, gives:
-
+
And the distance between these two is the standard Euclidean distance:
-
+
So, what does this distance function look like when we plot it for a number of ranges for the angle φ, such as a half circle, quarter
circle and eighth circle?
@@ -5133,7 +6210,7 @@ for p = 1 to points.length-3 (inclusive):
In fact, let's flip the function around, so that if we plug in the precision error, labelled ε, we get back the maximum angle for that
precision:
-
+
And frankly, things are starting to look a bit ridiculous at this point, we're doing way more maths than we've ever done, but thankfully
this is as far as we need the maths to take us: If we plug in the precisions 0.1, 0.01, 0.001 and 0.0001 we get the radians values 1.748,
@@ -5180,7 +6257,6 @@ for p = 1 to points.length-3 (inclusive):
-
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
@@ -5239,8 +6315,8 @@ for p = 1 to points.length-3 (inclusive):
@@ -5251,16 +6327,16 @@ for p = 1 to points.length-3 (inclusive):
where "a" is some scaling factor we'll need to find the expression for, and:
@@ -5296,9 +6372,9 @@ for p = 1 to points.length-3 (inclusive):
@@ -5308,7 +6384,7 @@ for p = 1 to points.length-3 (inclusive):
@@ -5330,16 +6406,16 @@ for p = 1 to points.length-3 (inclusive):
@@ -5348,9 +6424,9 @@ for p = 1 to points.length-3 (inclusive):
@@ -5359,9 +6435,9 @@ for p = 1 to points.length-3 (inclusive):
And that's it, we have all four points now for an approximation of an arbitrary circular arc with angle φ.
@@ -5370,16 +6446,16 @@ for p = 1 to points.length-3 (inclusive):
So, to recap, given an angle φ, the new control coordinates are:
and
@@ -5389,16 +6465,16 @@ for p = 1 to points.length-3 (inclusive):
Which, in decimal values, rounded to six significant digits, is:
@@ -5493,7 +6569,6 @@ for p = 1 to points.length-3 (inclusive):
-
With that in place, all that's left now is to "restart" the procedure by treating the found arc's end point as the new to-be-determined
arc's starting point, and using points further down the curve. We keep trying this until the found end point is for t=1, at which
@@ -5508,7 +6583,6 @@ for p = 1 to points.length-3 (inclusive):
-
So... what is this good for? Obviously, if you're working with technologies that can't do curves, but can do lines and circles, then the
answer is pretty straightforward, but what else? There are some reasons why you might need this technique: using circular arcs means you
@@ -5551,7 +6625,6 @@ for p = 1 to points.length-3 (inclusive):
-
The important part to notice here is that we are not doing the same thing with B-Splines that we do for poly-Béziers or
Catmull-Rom curves: both of the latter simply define new sections as literally "new sections based on new points", so a 12 point cubic
@@ -5575,7 +6648,6 @@ for p = 1 to points.length-3 (inclusive):
-
In order to make this interpolation of curves work, the maths is necessarily more complex than the maths for Bézier curves, so let's have
a look at how things work.
@@ -5587,7 +6659,7 @@ for p = 1 to points.length-3 (inclusive):
>, we can compute a point on the curve for some value t in the interval [0,1] (where 0 is the start of the curve, and 1 the
end, just like for Bézier curves), by evaluating the following function:
-
+
Which, honestly, doesn't tell us all that much. All we can see is that a point on a B-Spline curve is defined as "a mix of all the control
points, weighted somehow", where the weighting is achieved through the N(...) function, subscripted with an obvious parameter
@@ -5602,14 +6674,14 @@ for p = 1 to points.length-3 (inclusive):
k subscript to the N() function applies to.
Then the N() function itself. What does it look like?
-
+
So this is where we see the interpolation: N(t) for an (i,k) pair (that is, for a step in the above summation, on a specific
knot interval) is a mix between N(t) for (i,k-1) and N(t) for (i+1,k-1), so we see that this is a recursive
iteration where i goes up, and k goes down, so it seem reasonable to expect that this recursion has to stop at
some point; obviously, it does, and specifically it does so for the following i/k values:
-
+
And this function finally has a straight up evaluation: if a t value lies within a knot-specific interval once we reach a
k=1 value, it "counts", otherwise it doesn't. We did cheat a little, though, because for all these values we need to scale
@@ -5626,12 +6698,12 @@ for p = 1 to points.length-3 (inclusive):
Carl de Boor — came to a mathematically pleasing solution: to compute a point
P(t), we can compute this point by evaluating d(t) on a curve section between knots i and i+1:
-
+
This is another recursive function, with k values decreasing from the curve order to 1, and the value α (alpha) defined
by:
-
+
That looks complicated, but it's not. Computing alpha is just a fraction involving known, plain numbers. And, once we have our alpha
value, we also have (1-alpha) because it's a trivial subtraction. Computing the d() function is thus mostly a
@@ -5639,7 +6711,7 @@ for p = 1 to points.length-3 (inclusive):
recursion might see computationally expensive, the total algorithm is cheap, as each step only involves very simple maths.
Of course, the recursion does need a stop condition:
-
+
So, we actually see two stopping conditions: either i becomes 0, in which case d() is zero, or
k becomes zero, in which case we get the same "either 1 or 0" that we saw in the N() function above.
@@ -5649,7 +6721,7 @@ for p = 1 to points.length-3 (inclusive):
Casteljau's algorithm. For instance, if we write out d() for i=3 and k=3, we get the following
recursion diagram:
-
+
That is, we compute d(3,3) as a mixture of d(2,3) and d(2,2), where those two are themselves a
mixture of d(1,3) and d(1,2), and d(1,2) and d(1,1), respectively, which are
@@ -5675,23 +6747,78 @@ for p = 1 to points.length-3 (inclusive):
[knots[degree], knots[knots.length - 1 - degree]. Then, we find the section number s that this mapped
t value lies on:
- for(s=domain[0]; s < domain[1]; s++) {
- if(knots[s] <= t && t <= knots[s+1]) break;
-}
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+
after running this code, s is the index for the section the point will lie on. We then run the algorithm mentioned on the MU
page (updated to use this description's variable names):
- let v = copy of control points
-for(let L = 1; L <= order; L++) {
- for(let i=s; i > s + L - order; i--) {
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+ 8
+
+
+ 9
+
+
+ 10
+
+
+
(A nice bit of behaviour in this code is that we work the interpolation "backwards", starting at i=s at each level of the
interpolation, and we stop when i = s - order + level, so we always end up with a value for i such that those
@@ -5746,7 +6873,6 @@ for(let L = 1; L <= order; L++) {
-
This is an important point: the intervals that the knot vector defines are relative intervals, so it doesn't matter if every
interval is size 1, or size 100 - the relative differences between the intervals is what shapes any particular curve.
@@ -5771,7 +6897,6 @@ for(let L = 1; L <= order; L++) {
-
Open-Uniform B-Splines
By combining knot interval collapsing at the start and end of the curve, with uniform knots in between, we can overcome the problem of the
@@ -5792,7 +6917,6 @@ for(let L = 1; L <= order; L++) {
-
Non-uniform B-Splines
This is essentially the "free form" version of a B-Spline, and also the least interesting to look at, as without any specific reason to
@@ -5815,7 +6939,6 @@ for(let L = 1; L <= order; L++) {
-
Of course this brings us to the final topic that any text on B-Splines must touch on before calling it a day: the
NURBS, or Non-Uniform Rational B-Spline (NURBS is not a plural,
diff --git a/docs/news/00-draft.md b/docs/news/00-draft.md
new file mode 100644
index 00000000..04f30b66
--- /dev/null
+++ b/docs/news/00-draft.md
@@ -0,0 +1,30 @@
+# Writing a web page
+
+## Let's talk about React
+
+- React vs HTML: web apps are not web pages
+
+## Let's talk about web pages
+
+- web pages are not web sites
+- lazy loading your assets
+- HTML vs "don't even write HTML"
+- the modern web: distrust as a baseline.
+ - why should people trust your JS?
+ - will it drain their battery?
+ - will it mine crypto?
+ - will it break if someone doesn't trust Google?
+ - etc.
+ - pregenerate everything that can be generated. You're building a web page, not a web page builder. It should work immediately.
+
+## Let's talk about modern JS
+
+- modern JS
+ `const`, `let`, `...`, `??`, `Object.fromEntries`/`Object.entries`, `async`+`await`
+ `async` and `defer` attributes
+ esmodules (remember `type="module"`!)
+ custom elements
+ - not quite react
+ - let's invent a new framework!
+ - mutation observer for attribute and child changes
+
diff --git a/docs/news/01-draft.md b/docs/news/01-draft.md
new file mode 100644
index 00000000..4111a389
--- /dev/null
+++ b/docs/news/01-draft.md
@@ -0,0 +1,19 @@
+# Writing a static site generator
+
+- markdown content
+ - "github flavoured" markdown
+ - latex blocks
+ - custom elements
+ - simple code formatting
+
+- localisation per chapter
+ - `en-GB`, `zh-CN`, `ja-JP`, arbitrary locale code
+
+- nunjucks templates:
+ - main index page
+ - news index + posts pages
+ - RSS feed
+
+- localisation of template strings in JS
+- prettification
+- link checker
diff --git a/docs/news/02-draft.md b/docs/news/02-draft.md
new file mode 100644
index 00000000..4231431c
--- /dev/null
+++ b/docs/news/02-draft.md
@@ -0,0 +1,14 @@
+# Processing LaTeX in Markdown. Using LaTeX
+
+extract blocks, replace with templating identifier
+ - hash filename based on latex code
+ - wrap in (locale-appropriate) LaTeX preamble
+ - send over to `xelatex`
+ - send that over to `pdfcrop`
+ - send that over to `pdf2svg`
+ - send that over to `svgo`
+ - put that in `images/chapters/.../` or `images/news/.../`
+
+replace the template identifier with lazy-loading `![]()
` pointing at the generated `.svg` image.
+
+have JS in place that strips lazy-loading and replaces it with longer-distance IntersectionObserver for better lazy-loading.
diff --git a/docs/news/03-draft.md b/docs/news/03-draft.md
new file mode 100644
index 00000000..dd376b61
--- /dev/null
+++ b/docs/news/03-draft.md
@@ -0,0 +1,16 @@
+# HTML custom elements and reinventing Processing.js
+
+- `<graphics-element>`
+- BaseAPI + GraphicsAPI
+- `<fallback-image>
+ - generate snapshot of first frame
+ - canvas in Node!
+ - image filename hashed based on file-content
+ - store in images/chapter/... or images/news/...
+ - lazy loading for if there is no JS
+- if JS, IntersectionObserver on the graphics element
+
+- Styling: `:defined` to the rescue
+
+- Why not `extends HTMLCanvasElement`?
+ - node + browser
diff --git a/docs/news/2020-09-18.html b/docs/news/2020-09-18.html
index 65f4ed9d..cdb54813 100644
--- a/docs/news/2020-09-18.html
+++ b/docs/news/2020-09-18.html
@@ -27,7 +27,7 @@
-
+
@@ -165,13 +165,76 @@
- 
+ 
- "Lazy loaded everything", so that you get what you need, only when or even just before you need it.
- Localized content based on a simple filenaming scheme.
- Nicely formatted HTML, CSS, and JS thanks to
prettier.
+ - with some code formatting so that there are line numbers without needing JS:
+
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+ 8
+
+
+ 9
+
+
+ 10
+
+
+ 11
+
+
+ 12
+
+
+ 13
+
+
+
+
- Responsive CSS, so the content intelligently reflows where possible.
- A "Live build" setup for working on the content and code.
- Automatic link-checking to make sure none of the links in the Primer lead you to a 404.
diff --git a/docs/news/2020-09-18.md b/docs/news/2020-09-18.md
index f67feb7f..73a4fef8 100644
--- a/docs/news/2020-09-18.md
+++ b/docs/news/2020-09-18.md
@@ -36,6 +36,24 @@ To give a bit of a teaser, some of the things I'll be writing about:
- "Lazy loaded everything", so that you get what you need, only when or even just before you need it.
- Localized content based on a simple filenaming scheme.
- Nicely formatted HTML, CSS, and JS thanks to `prettier`.
+- with some code formatting so that there are line numbers without needing JS:
+
+```
+let curve;
+
+setup() {
+ curve = Bezier.defaultCubic();
+ setMovable(curve.points);
+ noGrid();
+}
+
+draw() {
+ clear(`lightblue`);
+ curve.drawCurve();
+ curve.drawPoints();
+}
+```
+
- Responsive CSS, so the content intelligently reflows where possible.
- A "Live build" setup for working on the content and code.
- Automatic link-checking to make sure none of the links in the Primer lead you to a 404.
diff --git a/docs/news/index.html b/docs/news/index.html
index 2d026af1..40fffda5 100644
--- a/docs/news/index.html
+++ b/docs/news/index.html
@@ -26,7 +26,7 @@
-
+
diff --git a/docs/news/rss.xml b/docs/news/rss.xml
index 83450466..8038efa8 100644
--- a/docs/news/rss.xml
+++ b/docs/news/rss.xml
@@ -6,18 +6,18 @@
News updates for the primer on Bézier Curves by Pomax
en-GB
- Wed Sep 23 2020 17:48:27 +00:00
+ Fri Sep 25 2020 15:42:02 +00:00
https://pomax.github.io/bezierinfo/images/og-image.png
A Primer on Bézier Curves
https://pomax.github.io/bezierinfo
-
-
Rewriting the tech stack
https://pomax.github.io/bezierinfo/news/2020-09-18.html
+
<p>Once upon a time, I needed to draw some Bezier curves because I was trying to create a Japanese kanji composition system that turned strokes into outlines, and that required knowing how to offset Bezier curves and... at the time (2011, time flies) there was no good single source of information for Bezier curves on the web. So I made one. Sure it started small, but it turns out that if you just keep adding bits to something, several years later you have quite the monster, and a single HTML file becomes intractible.</p>
<p>So, in 2016, when <a href="">React.js</a> exploded onto the scene, I rewrote the primer as a React app, and it became a lot easier to maintain. Like, <em>a lot</em> a lot. However, there was a downside: no JS meant no content. Sure, server-side rendering sort of existed, but not really, and because the Primer is hosted through github, there was no "server" to run. Plus, trying to rehydrate an app the size of the Primer from a giant HTML file had truly <em>dire</em> performance.</p>
<p>So I left it a regular React app, and every time I thought "wouldn't it be nice if it was just... a web page again?" the browser landscape just hadn't caught up. Finally, in 2020, things are different: with a global pandemic, and some vacation time, and something random causing me to look up the state of HTML custom elements, everything was pointing at it being time to finally, <em>finally</em>, turn the Primer back into a normal web page.</p>
@@ -44,13 +44,45 @@
<li>Real LaTeX code, that gets compiled into optimized SVG using <code>xelatex</code>, <code>pdfcrop</code>, <code>pdf2svg</code>, and <code>svgo</code>:</li>
</ul>
<blockquote>
-<img class="LaTeX SVG" src="./images/news/2020-09-18.html/24b6ead799277b0d73687d2fc5532a84.svg" width="401px" height="97px" loading="lazy">
+<img class="LaTeX SVG" src="./images/news/2020-09-18.html/15225da473048d8c7b5b473b89de0b66.svg" width="401px" height="97px" loading="lazy">
</blockquote>
<ul>
<li>"Lazy loaded everything", so that you get what you need, only when or even just before you need it.</li>
<li>Localized content based on a simple filenaming scheme.</li>
<li>Nicely formatted HTML, CSS, and JS thanks to <code>prettier</code>.</li>
+<li>with some code formatting so that there are line numbers without needing JS:</li>
+</ul>
+
+ <table class="code"><tr><td>1</td><td rowspan="13">
+ <textarea disabled rows="13" role="doc-example">let curve;
+
+setup() {
+ curve = Bezier.defaultCubic();
+ setMovable(curve.points);
+ noGrid();
+}
+
+draw() {
+ clear(`lightblue`);
+ curve.drawCurve();
+ curve.drawPoints();
+}</textarea>
+ </td></tr>
+<tr><td>2</td></tr>
+<tr><td>3</td></tr>
+<tr><td>4</td></tr>
+<tr><td>5</td></tr>
+<tr><td>6</td></tr>
+<tr><td>7</td></tr>
+<tr><td>8</td></tr>
+<tr><td>9</td></tr>
+<tr><td>10</td></tr>
+<tr><td>11</td></tr>
+<tr><td>12</td></tr>
+<tr><td>13</td></tr></table>
+
+<ul>
<li>Responsive CSS, so the content intelligently reflows where possible.</li>
<li>A "Live build" setup for working on the content and code.</li>
<li>Automatic link-checking to make sure none of the links in the Primer lead you to a 404.</li>
@@ -60,11 +92,10 @@
<p>Enjoy <a href="https://pomax.github.io/bezierinfo">The new Primer on Bézier Curves</a>, and if you find any problems, <a href="https://github.com/Pomax/BezierInfo-2/issues">you know where to go</a>.</p>
<p>See you in the next post!</p>
<p>— Pomax</p>
-
+
+
Thu Sep 17 2020 17:00:00 +00:00
2020-09-18.html
-
-
diff --git a/docs/style.css b/docs/style.css
index 55c3e816..c45ac266 100644
--- a/docs/style.css
+++ b/docs/style.css
@@ -25,6 +25,10 @@
--code-block-color: rgba(0, 0, 0, 0.03);
--heading-color: #cfe6ff;
--heading-text: #333;
+
+ --line-numbers: 2.2em;
+ --side-padding: 2em;
+ --line-height: 16px;
}
body {
@@ -194,7 +198,7 @@ div.note:before {
img.LaTeX.SVG {
margin: 0;
background: white;
- padding: 1em 1em 1em 0;
+ padding: 1em 1em 0.25em 0;
border-radius: 0.3em;
}
@@ -380,3 +384,55 @@ div.github img area {
h5.post-date {
text-align: right;
}
+
+table.code {
+ border-collapse: collapse;
+ margin: 0;
+ margin-left: var(--side-padding);
+ width: calc(100% - 3 * var(--side-padding));
+ border: 0.5em solid #f0f0f0;
+ border-width: 0.5em 0;
+}
+
+table.code td {
+ background: #f0f0f0;
+ border: none;
+ padding: 0;
+ cursor: default;
+}
+
+table.code td,
+table.code td textarea {
+ font-family: Courier;
+ font-size: 16px;
+ white-space: pre;
+}
+
+table.code td textarea {
+ position: absolute;
+ top: 0;
+ left: 0;
+ width: 100%;
+ height: 100%;
+ border: none;
+ resize: none;
+ cursor: text;
+ color: inherit;
+ padding-left: 1em;
+ background: transparent;
+ overflow: hidden;
+}
+
+table.code tr td[rowspan] {
+ position: relative;
+}
+
+table.code tr td:first-child {
+ width: 1.5em;
+ padding-right: 0.25em;
+ text-align: right;
+ user-select: none;
+ color: #999;
+ cursor: default;
+ border-right: 1px solid red;
+}
diff --git a/docs/zh-CN/index.html b/docs/zh-CN/index.html
index a39514bb..8aa3d1b7 100644
--- a/docs/zh-CN/index.html
+++ b/docs/zh-CN/index.html
@@ -33,7 +33,7 @@
-
+
@@ -449,7 +449,7 @@
>的结果。这听起来很复杂，但你在很小的时候就做过线性插值：当你指向两个物体中的另外一个物体时，你就用到了线性插值。它就是很简单的“选出两点之间的一个点”。
如果我们知道两点之间的距离，并想找出离第一个点20%间距的一个新的点(也就是离第二个点80%的间距)，我们可以通过简单的计算来得到：
-
+
让我们来通过实际操作看一下：下面的图形都是可交互的，因此你可以通过上下键来增加或减少插值距离，来观察图形的变化。我们从三个点构成的两条线段开始。通过对各条线段进行线性插值得到两个点，对点之间的线段再进行线性插值，产生一个新的点。最终这些点——所有的点都可以通过选取不同的距离插值产生——构成了贝塞尔曲线
：
@@ -462,6 +462,7 @@
+
这为我们引出了复杂的数学：微积分。
虽然我们刚才好像没有用到这个，我们实际上只是逐步地画了一条二次曲线，而不是一次画好。贝塞尔曲线的一个很棒的特性就是它们可以通过多项式方程表示，也可以用很简单的插值形式表示。因此，反过来说，我们可以基于“真正的数学”（检查方程式，导数之类的东西），也可以通过观察曲线的“机械”构成（比如说，可以得知曲线永远不会延伸超过我们用来构造它的点），来看看这些曲线能够做什么。
@@ -476,26 +477,26 @@
贝塞尔曲线是“参数”方程的一种形式。从数学上讲，参数方程作弊了：“方程”实际上是一个从输入到唯一输出的、良好定义的映射关系。几个输入进来，一个输出返回。改变输入变量，还是只有一个输出值。参数方程在这里作弊了。它们基本上干了这么件事，“好吧，我们想要更多的输出值，所以我们用了多个方程”。举个例子：假如我们有一个方程，通过一些计算，将假设为x的一些值映射到另外的值:
-
+
记号f(x)是表示函数的标准方式（为了方便起见，如果只有一个的话，我们称函数为f），函数的输出根据一个变量（本例中是x）变化。改变x，f(x)的输出值也会变。
到目前没什么问题。现在，让我们来看一下参数方程，以及它们是怎么作弊的。我们取以下两个方程：
-
+
这俩方程没什么让人印象深刻的，只不过是正弦函数和余弦函数，但正如你所见，输入变量有两个不同的名字。如果我们改变了a的值，f(b)的输出不会有变化，因为这个方程没有用到a。参数方程通过改变这点来作弊。在参数方程中，所有不同的方程共用一个变量，如下所示：
多个方程，但只有一个变量。如果我们改变了t的值，fa(t)和fb(t)的输出都会发生变化。你可能会好奇这有什么用，答案其实很简单：对于参数曲线，如果我们用常用的标记来替代fa(t)和fb(t)，看起来就有些明朗了：
-
+
好了，通过一些神秘的t值将x/y坐标系联系起来。
所以，参数曲线不像一般函数那样，通过x坐标来定义y坐标，而是用一个“控制”变量将它们连接起来。如果改变t的值，每次变化时我们都能得到两个值，这可以作为图形中的(x,y)坐标。比如上面的方程组，生成位于一个圆上的点：我们可以使t在正负极值间变化，得到的输出(x,y)都会位于一个以原点(0,0)为中心且半径为1的圆上。如果我们画出t从0到5时的值，将得到如下图像：
@@ -508,16 +509,15 @@
-
贝塞尔曲线是（一种）参数方程，并在它的多个维度上使用相同的基本方程。在上述的例子中x值和y值使用了不同的方程，与此不同的是，贝塞尔曲线的x和y都用了“二项多项式”。那什么是二项多项式呢？
你可能记得高中所学的多项式，看起来像这样：
@@ -528,17 +528,17 @@
我明白你在想什么：这看起来并不简单，但如果我们拿掉t并让系数乘以1，事情就会立马简单很多，看看这些二次项：
@@ -549,8 +549,8 @@
@@ -559,25 +559,54 @@
这就是贝塞尔曲线完整的描述。在这个函数中的Σ表示了这是一系列的加法（用Σ下面的变量，从...=<值>开始，直到Σ上面的数字结束）。
如何实现基本方程
我们可以用之前说过的方程，来简单地实现基本方程作为数学构造，如下：
- function Bezier(n,t):
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+
我说我们“可以用”是因为我们不会这么去做：因为阶乘函数开销非常大。并且，正如我们在上面所看到的，我们不用阶乘也能够很容易地构造出帕斯卡三角形：一开始是[1]，接着是[1,2,1],然后是[1,3,3,1]等等。下一行都比上一行多一个数，首尾都为1，中间的数字是上一行两边元素的和。
我们可以很快的生成这个列表，并在之后使用这个查找表而不用再计算二次多项式的系数：
- lut = [ [1], // n=0
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+ 8
+
+
+ 9
+
+
+ 10
+
+
+ 11
+
+
+ 12
+
+
+ 13
+
+
+ 14
+
+
+ 15
+
+
+ 16
+
+
+ 17
+
+
+ 18
+
+
+
这里做了些什么？首先，我们声明了一个足够大的查找表。然后，我们声明了一个函数来获取我们想要的值，并且确保当一个请求的n/k对不在LUT查找表中时，先将表扩大。我们的基本函数如下所示：
- function Bezier(n,t):
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+
完美。当然我们可以进一步优化。为了大部分的计算机图形学目的，我们不需要任意的曲线。我们需要二次曲线和三次曲线（实际上这篇文章没有涉及任意次的曲线，因此你会在其他地方看到与这些类似的代码），这说明我们可以彻底简化代码:
- function Bezier(2,t):
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+ 8
+
+
+ 9
+
+
+ 10
+
+
+ 11
+
+
+ 12
+
+
+ 13
+
+
+
现在我们知道如何代用码实现基本方程了。很好。
@@ -645,34 +800,38 @@ function Bezier(3,t):
-
-
- Scripts are disabled. Showing fallback image.
-
-
-
-
-
+
+
+
+ Scripts are disabled. Showing fallback image.
+
+
+
+
+
+
-
-
- Scripts are disabled. Showing fallback image.
-
-
-
-
-
+
+
+
+ Scripts are disabled. Showing fallback image.
+
+
+
+
+
+
上面有一张是15th阶的插值方程。如你所见，在所有控制点中，起点和终点对曲线形状的贡献比其他点更大些。
如果我们要改变曲线，就需要改变每个点的权重，有效地改变插值。可以很直接地做到这个：只要用一个值乘以每个点，来改变它的强度。这个值照惯例称为“权重”，我们可以将它加入我们原始的贝塞尔函数：
-
+
看起来很复杂，但实际上“权重”只是我们想让曲线所拥有的坐标值：对于一条nth阶曲线，w0是起始坐标，wn是终点坐标，中间的所有点都是控制点坐标。假设说一条曲线的起点为（120，160），终点为（220，40），并受点（35，200）和点（220，260）的控制，贝塞尔曲线方程就为：
-
+
这就是我们在文章开头看到的曲线：
鉴于我们已经知道怎样实现基本函数，在其加入控制点是非常简单的：
-function Bezier(n,t,w[]):
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+
下面是优化过的版本：
- function Bezier(2,t,w[]):
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+ 8
+
+
+ 9
+
+
+ 10
+
+
+ 11
+
+
+ 12
+
+
+ 13
+
+
+
现在我们知道如何编程实现基本权重函数了。
Adding these ratio values to the regular Bézier curve function is fairly easy. Where the regular function is the following:
The function for rational Bézier curves has two more terms:
@@ -760,7 +989,6 @@ function Bezier(3,t,w[]): -
You can think of the ratio values as each coordinate's "gravity": the higher the gravity, the closer to that coordinate the curve will want to be. You'll also notice that if you simply increase or decrease all the ratios by the same amount, nothing changes... much like @@ -770,7 +998,13 @@ function Bezier(3,t,w[]):
Extending the code of the previous section to include ratios is almost trivial:
-function RationalBezier(2,t,w[],r[]):
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+ 8
+
+
+ 9
+
+
+ 10
+
+
+ 11
+
+
+ 12
+
+
+ 13
+
+
+ 14
+
+
+ 15
+
+
+ 16
+
+
+ 17
+
+
+ 18
+
+
+ 19
+
+
+ 20
+
+
+ 21
+
+
+ 22
+
+
+ 23
+
+
+ 24
+
+
+ 25
+
+
+ 26
+
+
+
And that's all we have to do.
t=0到t=1。为什么是这个特殊区间？
这一切都与我们如何从曲线的“起点”变化到曲线“终点”有关。如果有一个值是另外两个值的混合，一般方程如下：
-
+
很显然，起始值需要a=1, b=0，混合值就为100%的value 1和0%的value 2。终点值需要a=0, b=1，则混合值是0%的value
1和100%的value
2。另外，我们不想让“a”和“b”是互相独立的:如果它们是互相独立的话，我们可以任意选出自己喜欢的值，并得到混合值，比如说100%的value1和100%的value2。原则上这是可以的，但是对于贝塞尔曲线来说，我们通常想要的是起始值和终点值之间的混合值，所以要确保我们不会设置一些“a”和"b"而导致混合值超过100%。这很简单：
+
用这个式子我们可以保证相加的值永远不会超过100%。通过将a限制在区间[0，1],我们将会一直处于这两个值之间（包括这两个端点），并且相加为100%。
通过将贝塞尔公式表示成一个多项式基本方程、系数矩阵以及实际的坐标，我们也可以用矩阵运算来表示贝塞尔。让我们看一下这对三次曲线来说有什么含义：
-
+
暂时不用管我们具体的坐标，现在有：
-
+
可以将它写成四个表达式之和：
-
+
我们可以扩展这些表达式：
-
+
更进一步，我们可以加上所有的1和0系数，以便看得更清楚：
-
+
现在，我们可以将它看作四个矩阵运算：
-
+
如果我们将它压缩到一个矩阵操作里，就能得到：
-
+
这种多项式表达式一般是以递增的顺序来写的，所以我们应该将t矩阵水平翻转，并将大的那个“混合”矩阵上下颠倒：
+
最终，我们可以加入原始的坐标，作为第三个单独矩阵：
-
+
我们可以对二次曲线运用相同的技巧，可以得到：
-
+
如果我们代入t值并乘以矩阵来计算，得到的值与解原始多项式方程或用逐步线性插值计算的结果一样。
因此：为什么我们要用矩阵来计算？ 用矩阵形式来表达曲线可以让我们去探索函数的一些很难被发现的性质。可以证明的是曲线构成了
-
让我们使用刚才描述过的算法，并实现它： 好了，这就是算法的实现。一般来说你不能随意重载“+”操作符，因此我们给出计算
以上算法做了什么？如果参数points列表只有一个点，
就画出一个点。如果有多个点，就生成以t为比例的一系列点（例如，以上算法中的"标记点"），然后为新的点列表调用绘制函数。
@@ -974,14 +1357,16 @@ function RationalBezier(3,t,w[],r[]):
-
+
@@ -990,21 +1375,79 @@ function RationalBezier(3,t,w[],r[]):
让我们来实现刚才简述过的算法： 好了，这就是算法的实现。它基本上是画出一系列的线段来模拟“曲线”。 我们将第一个坐标作为参考点，然后在相邻两个点之间画线。 通过在 de Casteljau 函数里插入一些额外的输出代码，我们就可以实现曲线的分割：
对某个给定
However, this rule also has as direct consequence that you cannot generally safely lower a curve from
nth order to (n-1)th order, because the control points cannot be "pulled apart" cleanly. We can
@@ -1344,27 +1842,27 @@ function drawCurve(points[], t):
some things can be done much more easily with matrices than with calculus functions, and this is one of those things. So... let's go!
We start by taking the standard Bézier function, and condensing it a little:
Then, we apply one of those silly (actually, super useful) calculus tricks: since our
So, with that seemingly trivial observation, we rewrite that Bézier function by splitting it up into a sum of a
So far so good. Now, to see why we did this, let's write out the
@@ -1376,9 +1874,9 @@ function drawCurve(points[], t):
@@ -1392,8 +1890,8 @@ function drawCurve(points[], t):
Let's do this: where the matrix M is an
@@ -1425,9 +1923,9 @@ function drawCurve(points[], t):
The steps taken here are: First, let's look at the derivative rule for Bézier curves, which is:
@@ -1677,9 +2175,9 @@ function drawCurve(points[], t):
Which is the "long" version of the following matrix transformation:
However, if you've played with that graphic a bit, you might have noticed something odd. The normal seems to "suddenly twist around the
curve" between t=0.65 and t=0.75... Why is it doing that?
@@ -1878,14 +2375,46 @@ function drawCurve(points[], t):
We first assume we have a function for calculating the Frenet frame at a point, which we already discussed above, inn a way that it
yields a frame with properties:
Then, we can write a function that generates a sequence of RM frames in the following manner:
Ignoring comments, this is certainly more code than when we were just computing a single Frenet frame, but it's not a crazy amount more
code to get much better looking normals.
@@ -1939,7 +2578,6 @@ function drawCurve(points[], t):
-
That looks so much better!
For those reading along with the code: we don't even strictly speaking need a Frenet frame to start with: we could, for instance, treat
@@ -1986,7 +2624,6 @@ function drawCurve(points[], t):
And then we turn this into our solution for
@@ -2064,17 +2701,17 @@ function drawCurve(points[], t):
And then, using these v values, we can find out what our a, b, and c should be:
@@ -2083,9 +2720,9 @@ function drawCurve(points[], t):
Easy-peasy. We can now almost trivially find the roots by plugging those values into the quadratic formula.
@@ -2139,14 +2776,20 @@ function drawCurve(points[], t):
"complex" number space (denoted with ℂ). And, as it so happens, Cardano is
also attributed as the first mathematician in history to have made use of complex numbers in his calculations. For this very algorithm!
@@ -2262,9 +3149,9 @@ function getCubicRoots(pa, pb, pc, pd) {
(The Wikipedia article has a decent animation for this process, so I will not add a graphic for that here) And for cubic curves, that means first and second derivatives, in red and purple respectively:
Then translating it so that the first coordinate lies on (0,0), moving all x coordinates by -120, and all y coordinates
by -160, gives us:
If we then rotate the curve so that its end point lies on the x-axis, the coordinates (integer-rounded for illustrative purposes here)
become:
If we drop all the zero-terms, this gives us:
We can see that our original curve definition has been simplified considerably. The following graphics illustrate the result of aligning
our example curves to the x-axis, with the cubic case using the coordinates that were just used in the example formulae:
@@ -2406,7 +3292,6 @@ function getCubicRoots(pa, pb, pc, pd) {
What we need to do is solve a simple equation:
What we're saying here is that given the curvature function C(t), we want to know for which values of t this function is
zero, meaning there is no "curvature", which will be exactly at the point between our circle being on one side of the curve, and our
@@ -2482,9 +3367,9 @@ function getCubicRoots(pa, pb, pc, pd) {
@@ -2506,16 +3391,16 @@ function getCubicRoots(pa, pb, pc, pd) {
And of course the same functions for y:
@@ -2542,9 +3427,9 @@ function getCubicRoots(pa, pb, pc, pd) {
@@ -2553,17 +3438,17 @@ function getCubicRoots(pa, pb, pc, pd) {
This is a plain quadratic curve, and we know how to solve C(t) = 0; we use the quadratic formula:
@@ -2614,7 +3499,6 @@ function getCubicRoots(pa, pb, pc, pd) {
This is a fairly funky image, so let's see what the various parts of it mean...
We see the three fixed points at (0,0), (0,1) and (1,1). The various regions and boundaries indicate what property the original curve will
@@ -2634,9 +3518,9 @@ function getCubicRoots(pa, pb, pc, pd) {
...on the top (blue) edge, the curve's start point touches the curve, forming a loop. This edge is described by the function:
Sweet! z stays 1, so we can effectively ignore it entirely, but we added some plain values to our x and y coordinates. So, if we
want to subtract p1.x and p1.y, we use:
Running all our coordinates through this transformation gives a new set of coordinates, let's call those U, where the
first coordinate lies on (0,0), and the rest is still somewhat free. Our next job is to make sure point 2 ends up lying on the
@@ -2728,12 +3612,12 @@ function getCubicRoots(pa, pb, pc, pd) {
currently have. This is called shearing, and the typical x-shear matrix and its
transformation looks like this:
So we want some shearing value that, when multiplied by y, yields -x, so our x coordinate becomes zero. That value is
simply -x/y, because *-x/y * y = -x*. Done:
Now, running this on all our points generates a new set of coordinates, let's call those V, which now have point 1 on
(0,0) and point 2 on (0, some-value), and we wanted it at (0,1), so we need to
@@ -2741,7 +3625,7 @@ function getCubicRoots(pa, pb, pc, pd) {
point 3 to end up on (1,1), so we can also scale x to make sure its x-coordinate will be 1 after we run the transform. That means we'll be
x-scaling by 1/point3x, and y-scaling by point2y. This is really easy:
Then, finally, this generates a new set of coordinates, let's call those W, of which point 1 lies on (0,0), point 2 lies on (0,1), and
point three lies on (1, ...) so all that's left is to make sure point 3 ends up at (1,1) - but we can't scale! Point 2 is already in the
@@ -2750,13 +3634,13 @@ function getCubicRoots(pa, pb, pc, pd) {
but y-shearing. Additionally, we don't actually want to end up at zero (which is what we did before) so we need to shear towards an
offset, in this case 1:
And this generates our final set of four coordinates. Of these, we already know that points 1 through 3 are (0,0), (0,1) and (1,1), and
only the last coordinate is "free". In fact, given any four starting coordinates, the resulting "transformation mapped" coordinate will
be:
Okay, well, that looks plain ridiculous, but: notice that every coordinate value is being offset by the initial translation, and also
notice that a lot of terms in that expression are repeated. Even though the maths looks crazy as a single expression, we can just
@@ -2766,13 +3650,13 @@ function getCubicRoots(pa, pb, pc, pd) {
First, let's just do that translation step as a "preprocessing" operation so we don't have to subtract the values all the time. What does
that leave?
Suddenly things look a lot simpler: the mapped x is fairly straight forward to compute, and we see that the mapped y actually contains the
mapped x in its entirety, so we'll have that part already available when we need to evaluate it. In fact, let's pull out all those common
factors to see just how simple this is:
That's kind of super-simple to write out in code, I think you'll agree. Coding math tends to be easier than the formulae initially make it
look!
@@ -2839,7 +3723,6 @@ function getCubicRoots(pa, pb, pc, pd) {
-
Now, if you look more closely at that right graphic, you'll notice something interesting: if we treat the red line as "the x axis", then
the point where the function crosses our line is really just a root for the cubic function x(t) through a shifted "x-axis"... and
@@ -2847,23 +3730,29 @@ function getCubicRoots(pa, pb, pc, pd) {
cubic case either: because of the kind of curve we're starting with, we know there is only root, simplifying the code we need!
First, let's look at the function for x(t):
We can rewrite this to a plain polynomial form, by just fully writing out the expansion and then collecting the polynomial factors, as:
Nothing special here: that's a standard cubic polynomial in "power" form (i.e. all the terms are ordered by their power of
You might be wondering "where did all the other 'minus x' for all the other values a, b, c, and d go?" and the answer there is that they
all cancel out, so the only one we actually need to subtract is the one at the end. Handy! So now we just solve this equation using
Cardano's algorithm, and we're left with some rather short code:
So the procedure is fairly straight forward: pick an or, more commonly written using Leibnitz notation as:
This formula says that the length of a parametric curve is in fact equal to the area underneath a function that looks a
remarkable amount like Pythagoras' rule for computing the diagonal of a straight angled triangle. This sounds pretty simple, right? Sadly,
@@ -2926,7 +3847,7 @@ y = curve.get(t).y
In plain text: an integral function can always be treated as the sum of an (infinite) number of (infinitely thin) rectangular strips
sitting "under" the function's plotted graph. To illustrate this idea, the following graph shows the integral for a sinusoid function. The
@@ -2988,9 +3909,9 @@ y = curve.get(t).y
@@ -3004,15 +3925,15 @@ y = curve.get(t).y
values, for any n, so if we want to approximate our integral with only two terms (which is a bit low, really) then
these tables would tell us that for n=2 we must use the following values:
Which means that in order for us to approximate the integral, we must plug these values into the approximate function, which gives us:
@@ -3065,7 +3986,6 @@ y = curve.get(t).y
-
So, what does the function look like? This:
Which is really just a "short form" that glosses over the fact that we're dealing with functions of
And while that's a litte more verbose, it's still just as simple to work with as the first function: the curvature at some point on any
(and this cannot be overstated: any) curve is a ratio between the first and second derivative cross product, and something that
@@ -3164,13 +4084,42 @@ y = curve.get(t).y
evaluating this function for any t value is just a matter of basic arithematics.
In fact, let's just implement it right now:
That was easy! (Well okay, that "not a number" value will need to be taken into account by downstream code, but that's a reality of
programming anwyay)
@@ -3193,14 +4142,13 @@ y = curve.get(t).y
One thing you may have noticed in this sketch is that sometimes the curvature looks fine, but seems to be pointing in the wrong direction,
making it hard to line up the curves properly. A way around that, of course, is to show the curvature on both sides of the curve, so let's
just do that. But let's take it one step further: we can also compute the associated "radius of curvature", which gives us the implicit
circle that "fits" the curve's curvature at any point, using what is possibly the simplest bit of maths found in this entire primer:
So let's revisit the previous graphic with the curvature visualised on both sides of our curves, as well as showing the circle that "fits"
our curve at some point that we can control by using a slider:
@@ -3249,7 +4197,6 @@ y = curve.get(t).y
So, how do we "cut up" the arc length function at regular intervals, when we can't really work with it? We basically cheat: we run through
the curve using Use the slider to increase or decrease the number of equidistant segments used to colour the curve.
However, are there better ways? One such way is discussed in "
-
@@ -3327,7 +4272,13 @@ y = curve.get(t).y
but sometimes you need more function signatures, because you might not want to call your function with eight distinct parameters. Maybe
you're using point structs for the line. Let's get coding:
Finding self-intersections is effectively the same procedure, except that we're starting with a single curve, so we need to turn that into
two separate curves first. This is trivially achieved by splitting at an inflection point, or if there are none, just splitting at
@@ -3544,7 +4546,7 @@ lli = function(line1, line2):
So, how can we compute
If we can figure out what the function And
So, if we know the start and end coordinates, and we know the t value, we know C, without having to calculate the
We start by observing that, given Working out the maths for this, we see the following two formulae for quadratic and cubic curves: And
Which now leaves us with some powerful tools: given thee points (start, end, and "some point on the curve"), as well as a
With And then reverse engineer the curve's control control points:
So: if we have a curve's start and end point, then for any
@@ -3625,7 +4627,6 @@ lli = function(line1, line2):
-
For cubic curves we need to do a little more work, but really only just a little. We're first going to assume that a decent curve through
the three points should approximate a circular arc, which first requires knowing how to fit a circle to three points. You may remember (if
@@ -3645,7 +4646,6 @@ lli = function(line1, line2):
-
With that covered, we now also know the tangent line to our point
@@ -3670,9 +4670,9 @@ lli = function(line1, line2):
@@ -3682,9 +4682,9 @@ lli = function(line1, line2):
The result of this approach looks as follows:
It is important to remember that even though we're using a circular arc to come up with decent That looks perfectly servicable!
Of course, we can take this one step further: we can't just "create" curves, we also have (almost!) all the tools available to "mold"
@@ -3744,13 +4742,42 @@ lli = function(line1, line2):
binary search. First, we do a coarse distance-check based on
After this runs, we know that And then the associated And we're done, because that's our new quadratic control point!
As before, cubic curves are a bit more work, because while it's easy to find our initial
That looks reasonable, close to the original point, but the further we drag our point, the less "useful" things become. Especially if we
drag our point across the baseline, rather than turning into a nice curve.
@@ -3862,7 +4887,6 @@ for (coordinate, index) in LUT:
-
The slide controls the "falloff distance" relative to where the original point on the curve is, so that as we drag our point around, it
interpolates with a bias towards "preserving And then we (trivially) rearrange the terms across multiple lines:
@@ -3950,40 +4974,40 @@ for (coordinate, index) in LUT:
With that arrangement, we can easily decompose this as a matrix multiplication: We can do the same for the cubic curve, of course. We know the base function for cubics: So we write out the expansion and rearrange: Which we can then decompose: And, of course, we can do this for quartic curves too (skipping the expansion step): To get these values, we first compute the general "distance along the polygon" matrix: In which we can replace the rather cumbersome "squaring" operation with a more conventional matrix equivalent:
@@ -4105,25 +5129,25 @@ for (coordinate, index) in LUT:
Which, because of the first and last values in S, means: Now we can properly write out the error function as matrix operations:
@@ -4133,9 +5157,9 @@ for (coordinate, index) in LUT:
@@ -4203,7 +5227,6 @@ for (coordinate, index) in LUT:
-
You'll note there is a convenient "toggle" buttons that lets you toggle between equidistant
Now, it may look like Catmull-Rom curves are very different from Bézier curves, because these curves can get very long indeed, but what
looks like a single Catmull-Rom curve is actually a spline: a single
@@ -4246,9 +5268,9 @@ for (coordinate, index) in LUT:
@@ -4263,7 +5285,13 @@ for (coordinate, index) in LUT:
considerable more
complex):
Now, since a Catmull-Rom curve is a form of cubic Hermite spline, and as
cubic Bézier curves are also a form of cubic Hermite spline, we run into an interesting bit of maths programming: we can
@@ -4305,8 +5392,8 @@ for p = 1 to points.length-3 (inclusive):
@@ -4344,17 +5431,17 @@ for p = 1 to points.length-3 (inclusive):
Computing T is really more "arranging the numbers": Thus:
@@ -4365,9 +5452,9 @@ for p = 1 to points.length-3 (inclusive):
@@ -4376,73 +5463,73 @@ for p = 1 to points.length-3 (inclusive):
Replace point/tangent vector with the expression for all-coordinates: and merge the matrices: This looks a lot like the Bézier matrix form, which as we saw in the chapter on Bézier curves, should look like this: So, if we want to express a Catmull-Rom curve using a Bézier curve, we'll need to turn this Catmull-Rom bit: Into something that looks like this: And the way we do that is with a fairly straight forward bit of matrix rewriting. We start with the equality we need to ensure: Then we remove the coordinate vector from both sides without affecting the equality: Then we can "get rid of" the Bézier matrix on the right by left-multiply both with the inverse of the Bézier matrix:
@@ -4451,17 +5538,17 @@ for p = 1 to points.length-3 (inclusive):
And now we're basically done. We just multiply those two matrices and we know what V is: We now have the final piece of our function puzzle. Let's run through each step. And we're done: we finally know how to convert these two curves!
@@ -4545,17 +5632,17 @@ for p = 1 to points.length-3 (inclusive):
Or, if your API allows you to specify Catmull-Rom curves using plain coordinates:
Dragging the control points around only affects the curve segments that the control point belongs to, and moving an on-curve point leaves
the control points where they are, which is not the most useful for practical modelling purposes. So, let's add in the logic we need to
make things a little better. We'll start by linking up control points by ensuring that the "incoming" derivative at an on-curve point is
the same as it's "outgoing" derivative:
We can effect this quite easily, because we know that the vector from a curve's last control point to its last on-curve point is equal to
the derivative vector. If we want to ensure that the first control point of the next curve matches that, all we have to do is mirror that
last control point through the last on-curve point. And mirroring any point A through any point B is really simple:
So let's implement that and see what it gets us. The following two graphics show a quadratic and a cubic poly-Bézier curve again, but this
time moving the control points around moves others, too. However, you might see something unexpected going on for quadratic curves...
@@ -4667,7 +5753,6 @@ for p = 1 to points.length-3 (inclusive):
-
As you can see, quadratic curves are particularly ill-suited for poly-Bézier curves, as all the control points are effectively linked.
Move one of them, and you move all of them. Not only that, but if we move the on-curve points, it's possible to get a situation where a
@@ -4710,7 +5795,6 @@ for p = 1 to points.length-3 (inclusive):
-
Cubic curves are now better behaved when it comes to dragging control points around, but the quadratic poly-Bézier still has the problem
that moving one control points will move the control points and may ending up defining "the next" control point in a way that doesn't
@@ -4748,7 +5832,6 @@ for p = 1 to points.length-3 (inclusive):
-
Again, we see that cubic curves are now rather nice to work with, but quadratic curves have a new, very serious problem: we can move an
on-curve point in such a way that we can't compute what needs to "happen next". Move the top point down, below the left and right points,
@@ -4794,9 +5877,9 @@ for p = 1 to points.length-3 (inclusive):
@@ -4807,9 +5890,9 @@ for p = 1 to points.length-3 (inclusive):
@@ -4821,8 +5904,8 @@ for p = 1 to points.length-3 (inclusive):
@@ -4843,9 +5926,9 @@ for p = 1 to points.length-3 (inclusive):
@@ -4907,7 +5990,6 @@ for p = 1 to points.length-3 (inclusive):
-
You may notice that this may still lead to small 'jumps' in the sub-curves when moving the curve around. This is caused by the fact that
we're still performing a naive form of offsetting, moving the control points the same distance as the start and end points. If the curve
@@ -4974,7 +6055,6 @@ for p = 1 to points.length-3 (inclusive):
-
As you can see, things go horribly wrong quite quickly; even trying to approximate a quarter circle using a quadratic curve is a bad idea.
An eighth of a turns might look okay, but how okay is okay? Let's apply some maths and find out. What we're interested in is how far off
@@ -5038,33 +6117,33 @@ for p = 1 to points.length-3 (inclusive):
We start out with our start and end point, and for convenience we will place them on a unit circle (a circle around 0,0 with radius 1), at
some angle φ:
What we want to find is the intersection of the tangents, so we want a point C such that:
i.e. we want a point that lies on the vertical line through S (at some distance a from S) and also lies on the tangent line
through E (at some distance b from E). Solving this gives us:
First we solve for b: which yields: which we can then substitute in the expression for a:
A quick check shows that plugging these values for a and b into the expressions for Cx and Cy give
the same x/y coordinates for both "a away from A" and "b away from B", so let's continue: now that we know the
coordinate values for C, we know where our on-curve point T for t=0.5 (or angle φ/2) is, because we can just evaluate the Bézier
polynomial, and we know where the circle arc's actual point P is for angle φ/2:
We compute T, observing that if t=0.5, the polynomial values (1-t)², 2(1-t)t, and t² are 0.25, 0.5, and 0.25 respectively: Which, worked out for the x and y components, gives: And the distance between these two is the standard Euclidean distance:
So, what does this distance function look like when we plot it for a number of ranges for the angle φ, such as a half circle, quarter
circle and eighth circle?
@@ -5105,7 +6184,7 @@ for p = 1 to points.length-3 (inclusive):
In fact, let's flip the function around, so that if we plug in the precision error, labelled ε, we get back the maximum angle for that
precision:
And frankly, things are starting to look a bit ridiculous at this point, we're doing way more maths than we've ever done, but thankfully
this is as far as we need the maths to take us: If we plug in the precisions 0.1, 0.01, 0.001 and 0.0001 we get the radians values 1.748,
@@ -5152,7 +6231,6 @@ for p = 1 to points.length-3 (inclusive):
-
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
@@ -5211,8 +6289,8 @@ for p = 1 to points.length-3 (inclusive):
where "a" is some scaling factor we'll need to find the expression for, and:
@@ -5280,7 +6358,7 @@ for p = 1 to points.length-3 (inclusive):
@@ -5302,16 +6380,16 @@ for p = 1 to points.length-3 (inclusive):
@@ -5320,9 +6398,9 @@ for p = 1 to points.length-3 (inclusive):
@@ -5331,9 +6409,9 @@ for p = 1 to points.length-3 (inclusive):
And that's it, we have all four points now for an approximation of an arbitrary circular arc with angle φ. So, to recap, given an angle φ, the new control coordinates are: and Which, in decimal values, rounded to six significant digits, is:
With that in place, all that's left now is to "restart" the procedure by treating the found arc's end point as the new to-be-determined
arc's starting point, and using points further down the curve. We keep trying this until the found end point is for t=1, at which
@@ -5480,7 +6557,6 @@ for p = 1 to points.length-3 (inclusive):
-
So... what is this good for? Obviously, if you're working with technologies that can't do curves, but can do lines and circles, then the
answer is pretty straightforward, but what else? There are some reasons why you might need this technique: using circular arcs means you
@@ -5523,7 +6599,6 @@ for p = 1 to points.length-3 (inclusive):
The important part to notice here is that we are not doing the same thing with B-Splines that we do for poly-Béziers or
Catmull-Rom curves: both of the latter simply define new sections as literally "new sections based on new points", so a 12 point cubic
@@ -5547,7 +6622,6 @@ for p = 1 to points.length-3 (inclusive):
-
In order to make this interpolation of curves work, the maths is necessarily more complex than the maths for Bézier curves, so let's have
a look at how things work.
@@ -5559,7 +6633,7 @@ for p = 1 to points.length-3 (inclusive):
>, we can compute a point on the curve for some value
Which, honestly, doesn't tell us all that much. All we can see is that a point on a B-Spline curve is defined as "a mix of all the control
points, weighted somehow", where the weighting is achieved through the N(...) function, subscripted with an obvious parameter
@@ -5574,14 +6648,14 @@ for p = 1 to points.length-3 (inclusive):
Then the N() function itself. What does it look like?
So this is where we see the interpolation: N(t) for an
And this function finally has a straight up evaluation: if a
This is another recursive function, with k values decreasing from the curve order to 1, and the value α (alpha) defined
by:
That looks complicated, but it's not. Computing alpha is just a fraction involving known, plain numbers. And, once we have our alpha
value, we also have Of course, the recursion does need a stop condition:
So, we actually see two stopping conditions: either
That is, we compute
after running this code,
(A nice bit of behaviour in this code is that we work the interpolation "backwards", starting at
This is an important point: the intervals that the knot vector defines are relative intervals, so it doesn't matter if every
interval is size 1, or size 100 - the relative differences between the intervals is what shapes any particular curve.
@@ -5743,7 +6871,6 @@ for(let L = 1; L <= order; L++) {
-
By combining knot interval collapsing at the start and end of the curve, with uniform knots in between, we can overcome the problem of the
@@ -5764,7 +6891,6 @@ for(let L = 1; L <= order; L++) {
-
This is essentially the "free form" version of a B-Spline, and also the least interesting to look at, as without any specific reason to
@@ -5787,7 +6913,6 @@ for(let L = 1; L <= order; L++) {
-
Of course this brings us to the final topic that any text on B-Splines must touch on before calling it a day: the
NURBS, or Non-Uniform Rational B-Spline (NURBS is not a plural,
diff --git a/src/build/create-news-pages.js b/src/build/create-news-pages.js
index 118af164..a9b77b7c 100644
--- a/src/build/create-news-pages.js
+++ b/src/build/create-news-pages.js
@@ -48,14 +48,14 @@ async function generatePost(file, localeStrings) {
// split off the post's title
let post = await convertMarkDown(
+ data,
{
imagepath: path.join(paths.images, `news`, filename),
modulepubdir: `./news/`,
file: file,
id: filename,
},
- localeStrings,
- data
+ localeStrings
);
const title = post.substring(post.indexOf(`
如何实现de Casteljau算法
function drawCurve(points[], t):
+
+
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+ 8
+ x和y坐标的实现：function drawCurve(points[], t):
+
+
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+ 8
+
+
+ 9
+
+
+ 10
+
-
-
+
+
@@ -1026,11 +1469,16 @@ function RationalBezier(3,t,w[],r[]):
-
如何实现曲线的拉平
function flattenCurve(curve, segmentCount):
+
+
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+ function drawFlattenedCurve(curve, segmentCount):
+
+
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
分割曲线的代码实习
works, I can recommend the University of South Florida video lecture on the procedure, linked in this very paragraph. The general solution
we're looking for is the following:
- left=[]
+
+
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+ 8
+
+
+ 9
+
+
+ 10
+
+
+ 11
+
+
+ 12
+
+
+ 13
+
+
+ 14
+
+
+ 15
+
+
+ 16
+ t 值，该函数执行后，数组 left 和 right 将包含两条曲线的所有点的坐标 --
一条是t值左侧的曲线，一条是t值右侧的曲线， 与原始曲线同序且完全重合。
@@ -1064,7 +1562,7 @@ function drawCurve(points[], t):
and
and
Excellent! Now we can form our new quadratic curve:
If we want the interval [z,1], we will be evaluating this instead:
So, our final second curve looks like:
and
and
nth order curve to an (n+1)th order curve is as follows (observing
that the start and end weights are the same as the start and end weights for the old curve):
-
+
+
t value is always between zero and one
(inclusive), we know that (1-t) plus t always sums to 1. As such, we can express any value as a sum of
t and 1-t:
+
(1-t) and
t component:
+
(1-t) and t parts, and see what that gives us.
I promise, it's about to make sense. We start with (1-t):
@@ -1401,13 +1899,13 @@ function drawCurve(points[], t):
And this is where we switch over from calculus to linear algebra, and matrices: we can now express this relation between Bézier(n,t) and
Bézier(n+1,t) as a very simple matrix multiplication:
-
+
n+1 by n matrix, and looks like:
Which is hard to work with, so let's expand that properly:
And that's the first part done: the two components inside the parentheses are actually regular, lower-order Bézier expressions:
And that's just a summation of lower order curves:
We can rewrite this as a normal summation, and we're done:
@@ -1665,9 +2163,9 @@ function drawCurve(points[], t):
@@ -1694,16 +2192,16 @@ function drawCurve(points[], t):
@@ -1827,7 +2325,6 @@ function drawCurve(points[], t):
-
{
+
+
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+ generateRMFrames(steps) -> frames:
+
+
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+ 8
+
+
+ 9
+
+
+ 10
+
+
+ 11
+
+
+ 12
+
+
+ 13
+
+
+ 14
+
+
+ 15
+
+
+ 16
+
+
+ 17
+
+
+ 18
+
+
+ 19
+
+
+ 20
+
+
+ 21
+
+
+ 22
+
+
+ 23
+
+
+ 24
+
+
+ 25
+
+
+ 26
+
+
+ 27
+
+
+ 28
+
+
+ 29
+
+
+ 30
+
+
+ 31
+
+
+ 32
+
+
+ 33
+
+
+ 34
+
+
+ 35
+
+
+ 36
+
-
t using basic arithmetics:
@@ -2047,9 +2684,9 @@ function drawCurve(points[], t):
// A helper function to filter for values in the [0,1] interval:
+
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+ 8
+
+
+ 9
+
+
+ 10
+
+
+ 11
+
+
+ 12
+
+
+ 13
+
+
+ 14
+
+
+ 15
+
+
+ 16
+
+
+ 17
+
+
+ 18
+
+
+ 19
+
+
+ 20
+
+
+ 21
+
+
+ 22
+
+
+ 23
+
+
+ 24
+
+
+ 25
+
+
+ 26
+
+
+ 27
+
+
+ 28
+
+
+ 29
+
+
+ 30
+
+
+ 31
+
+
+ 32
+
+
+ 33
+
+
+ 34
+
+
+ 35
+
+
+ 36
+
+
+ 37
+
+
+ 38
+
+
+ 39
+
+
+ 40
+
+
+ 41
+
+
+ 42
+
+
+ 43
+
+
+ 44
+
+
+ 45
+
+
+ 46
+
+
+ 47
+
+
+ 48
+
+
+ 49
+
+
+ 50
+
+
+ 51
+
+
+ 52
+
+
+ 53
+
+
+ 54
+
+
+ 55
+
+
+ 56
+
+
+ 57
+
+
+ 58
+
+
+ 59
+
+
+ 60
+
+
+ 61
+
+
+ 62
+
+
+ 63
+
+
+ 64
+
+
+ 65
+
+
+ 66
+
+
+ 67
+
+
+ 68
+
+
+ 69
+
+
+ 70
+
+
+ 71
+
+
+ 72
+
+
+ 73
+
+
+ 74
+
+
+ 75
+
+
+ 76
+
+
+ 77
+
+
+ 78
+
+
+ 79
+
+
+ 80
+
+
+ 81
+
-
+
+
+
+
+
@@ -2525,9 +3410,9 @@ function getCubicRoots(pa, pb, pc, pd) {
@@ -2647,9 +3531,9 @@ function getCubicRoots(pa, pb, pc, pd) {
@@ -2657,8 +3541,8 @@ function getCubicRoots(pa, pb, pc, pd) {
@@ -2715,12 +3599,12 @@ function getCubicRoots(pa, pb, pc, pd) {
cx + dy), which means we can't do translation, since that requires we end up with some kind of (x + a, y + b). If we add a bogus
z coordinate that is always 1, then we can suddenly add arbitrary values. For example:
-
+
+
+
+
+
+
+
+
+
+
+
t). So, given that a, b, c, d, and x(t) are all
known constants, we can trivially rewrite this (by moving the x(t) across the equal sign) as:
+
fy(t), then the length of the curve, measured from start point to some point t = z, is computed using the
following seemingly straight forward (if a bit overwhelming) formula:
- // prepare our values for root finding:
+
+
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+ 8
+
+
+ 9
+
+
+ 10
+ x, find the associted t value, evaluate our curve
for that t value, which gives us the curve's {x,y} coordinate, which means we know y for this
@@ -2898,9 +3819,9 @@ y = curve.get(t).y
+
+
+
+
+
t, so let's expand that a
tiny bit:
+
function kappa(t, B):
+
+
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
-
+
-
t values, determine the distance-for-this-t-value at each point we generate during the run, and
@@ -3270,7 +4217,6 @@ y = curve.get(t).y
-
Implementing line-line intersections
lli8 = function(x1,y1,x2,y2,x3,y3,x4,y4):
+
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+ 8
+
+
+ 9
+
+
+ 10
+
+
+ 11
+
+
+ 12
+
+
+ 13
+
+
+ 14
+
+
+ 15
+
+
+ 16
+
+
+ 17
+ What about curve-line intersections?
@@ -3451,7 +4454,6 @@ lli = function(line1, line2):
-
C? We start with our observation that C always lies somewhere between the start and ends
points, so logically C will have a function that interpolates between those two coordinates:
+
u(t) looks like, we'll be done. Although we do need to remember that this
u(t) will have a different for depending on whether we're working with quadratic or cubic curves.
@@ -3553,9 +4555,9 @@ lli = function(line1, line2):
>
(with thanks to Boris Zbarsky) shows us the following two formulae:
+
+
A or even B coordinates. In fact, we can do the same for the ratio function: as another function of
@@ -3563,18 +4565,18 @@ lli = function(line1, line2):
pure function of t, too.
A, B, and C, the following always holds:
+
+
+
t value, we can contruct curves: we can compute C using the start and end points, and our
u(t) function, and once we have C, we can use our on-curve point (B) and the
ratio(t) function to find A:
+
A found, finding e1 and e2 for quadratic curves is a matter of running the linear
interpolation with t between start and A to yield e1, and between A and end to yield
@@ -3582,9 +4584,9 @@ lli = function(line1, line2):
distance ratio between e1 to B and B to e2 is the Bézier ratio (1-t):t,
we can reverse engineer v1 and v2:
+
+
t value we implicitly know all the ABC values, which (combined
with an educated guess on appropriate e1 and e2 coordinates for cubic curves) gives us the necessary information
@@ -3609,9 +4611,9 @@ lli = function(line1, line2):
B, because the tangent to any point on the circle is a line
through that point, perpendicular to the line from that point to the center. That just leaves marking appropriate points
@@ -3655,9 +4655,9 @@ lli = function(line1, line2):
e1 and e2 terms,
we're not trying to perfectly create a circular arc with a cubic curve (which is good, because we can't;
@@ -3716,7 +4715,6 @@ lli = function(line1, line2):
-
t values associated with the curve's "to draw" coordinates (using a lookup table, or LUT). This is pretty fast:
p = some point to project onto the curve
+
+
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+ LUT[i] is the coordinate on the curve in our LUT that is closest to the point we want
to project, so that's a pretty good initial guess as to what the best projection onto our curve is. To refine it, we note that LUT[i] is a
@@ -3801,9 +4828,9 @@ for (coordinate, index) in LUT:
initial B coordinate. We don't even need the latter: with our t value and "whever the cursor is" as target
B, we can compute the associated C:
+
A:
+
t value and ABC values, getting
those all-important e1 and e2 coordinates is going to pose a bit of a problem... in the section on curve
@@ -3836,7 +4862,6 @@ for (coordinate, index) in LUT:
-
t/e1/e2" closer to the original point, and bias
@@ -3929,17 +4953,17 @@ for (coordinate, index) in LUT:
@@ -4000,8 +5024,8 @@ for (coordinate, index) in LUT:
@@ -4028,8 +5052,8 @@ for (coordinate, index) in LUT:
@@ -4040,8 +5064,8 @@ for (coordinate, index) in LUT:
@@ -4056,8 +5080,8 @@ for (coordinate, index) in LUT:
@@ -4067,8 +5091,8 @@ for (coordinate, index) in LUT:
@@ -4080,17 +5104,17 @@ for (coordinate, index) in LUT:
@@ -4166,9 +5190,9 @@ for (coordinate, index) in LUT:
t values, and distance ratio
along the polygon formed by the points. Arguably more interesting is that once you have points to abstract a curve, you also get
@@ -4234,7 +5257,6 @@ for (coordinate, index) in LUT:
-
tension = some value greater than 0, defaulting to 1
+
+
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+ 8
+
+
+ 9
+
+
+ 10
+
+
+ 11
+
+
+ 12
+
+
+ 13
+
+
+ 14
+
+
+ 15
+
+
+ 16
+
+
+ 17
+
+
+ 18
+
+
+ 19
+
@@ -4331,9 +5418,9 @@ for p = 1 to points.length-3 (inclusive):
@@ -4480,9 +5567,9 @@ for p = 1 to points.length-3 (inclusive):
@@ -4490,9 +5577,9 @@ for p = 1 to points.length-3 (inclusive):
@@ -4500,9 +5587,9 @@ for p = 1 to points.length-3 (inclusive):
@@ -4510,9 +5597,9 @@ for p = 1 to points.length-3 (inclusive):
@@ -4520,9 +5607,9 @@ for p = 1 to points.length-3 (inclusive):
@@ -4621,20 +5708,19 @@ for p = 1 to points.length-3 (inclusive):
-
+
+
@@ -4833,9 +5916,9 @@ for p = 1 to points.length-3 (inclusive):
+
+
+
+
+
+
+
+
+
+
+
@@ -5223,16 +6301,16 @@ for p = 1 to points.length-3 (inclusive):
@@ -5268,9 +6346,9 @@ for p = 1 to points.length-3 (inclusive):
@@ -5361,16 +6439,16 @@ for p = 1 to points.length-3 (inclusive):
@@ -5465,7 +6543,6 @@ for p = 1 to points.length-3 (inclusive):
-
-
t in the interval [0,1] (where 0 is the start of the curve, and 1 the
end, just like for Bézier curves), by evaluating the following function:
+
k subscript to the N() function applies to.
+
(i,k) pair (that is, for a step in the above summation, on a specific
knot interval) is a mix between N(t) for (i,k-1) and N(t) for (i+1,k-1), so we see that this is a recursive
iteration where i goes up, and k goes down, so it seem reasonable to expect that this recursion has to stop at
some point; obviously, it does, and specifically it does so for the following i/k values:
+
t value lies within a knot-specific interval once we reach a
k=1 value, it "counts", otherwise it doesn't. We did cheat a little, though, because for all these values we need to scale
@@ -5598,12 +6672,12 @@ for p = 1 to points.length-3 (inclusive):
Carl de Boor — came to a mathematically pleasing solution: to compute a point
P(t), we can compute this point by evaluating d(t) on a curve section between knots i and i+1:
+
+
(1-alpha) because it's a trivial subtraction. Computing the d() function is thus mostly a
@@ -5611,7 +6685,7 @@ for p = 1 to points.length-3 (inclusive):
recursion might see computationally expensive, the total algorithm is cheap, as each step only involves very simple maths.
+
i becomes 0, in which case d() is zero, or
k becomes zero, in which case we get the same "either 1 or 0" that we saw in the N() function above.
@@ -5621,7 +6695,7 @@ for p = 1 to points.length-3 (inclusive):
Casteljau's algorithm. For instance, if we write out d() for i=3 and k=3, we get the following
recursion diagram:
+
d(3,3) as a mixture of d(2,3) and d(2,2), where those two are themselves a
mixture of d(1,3) and d(1,2), and d(1,2) and d(1,1), respectively, which are
@@ -5647,23 +6721,78 @@ for p = 1 to points.length-3 (inclusive):
[knots[degree], knots[knots.length - 1 - degree]. Then, we find the section number s that this mapped
t value lies on:
+
+ for(s=domain[0]; s < domain[1]; s++) {
- if(knots[s] <= t && t <= knots[s+1]) break;
-}
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+ s is the index for the section the point will lie on. We then run the algorithm mentioned on the MU
page (updated to use this description's variable names):
let v = copy of control points
-for(let L = 1; L <= order; L++) {
- for(let i=s; i > s + L - order; i--) {
+
+
+
+
+ 1
+
+
+
+
+
+ 2
+
+
+ 3
+
+
+ 4
+
+
+ 5
+
+
+ 6
+
+
+ 7
+
+
+ 8
+
+
+ 9
+
+
+ 10
+ i=s at each level of the
interpolation, and we stop when i = s - order + level, so we always end up with a value for i such that those
@@ -5718,7 +6847,6 @@ for(let L = 1; L <= order; L++) {
-
Open-Uniform B-Splines
Non-uniform B-Splines
`) + 4, post.indexOf(`
`));
diff --git a/src/build/get-all-chapter-files.js b/src/build/get-all-chapter-files.js
index 961e9fc9..d6b06aba 100644
--- a/src/build/get-all-chapter-files.js
+++ b/src/build/get-all-chapter-files.js
@@ -2,6 +2,8 @@ import glob from "glob";
import path from "path";
import paths from "../project-paths.js";
+import sectionOrder from "../../docs/chapters/toc.js";
+
function getAllChapterFiles() {
// async, by returning a Promise
return new Promise((resolve, reject) => {
@@ -11,6 +13,10 @@ function getAllChapterFiles() {
const locales = {};
files.forEach((file) => {
+ const dir = path.relative(paths.chapters, path.dirname(file));
+
+ if (!sectionOrder.includes(dir)) return;
+
let locale = file.match(/content\.([^.]+)\.md/)[1];
if (!locales[locale]) {
locales[locale] = [];
diff --git a/src/build/markdown/convert-markdown.js b/src/build/markdown/convert-markdown.js
index a7873981..8698edec 100644
--- a/src/build/markdown/convert-markdown.js
+++ b/src/build/markdown/convert-markdown.js
@@ -1,34 +1,110 @@
-import convert from "./convert.js";
-import latexToSVG from "../latex/latex-to-svg.js";
-import preprocessGraphicsElement from "./preprocess-graphics-element.js";
-import extractLaTeX from "./extract-latex.js";
+import convert from "./util/convert.js";
+
+import latexToSVG from "./processors/latex/latex-to-svg.js";
+import processGraphicsElement from "./processors/graphics-element/process-graphics-element.js";
+import markupCodeBlock from "./processors/code/markup-code-block.js";
+
import nunjucks from "nunjucks";
nunjucks.configure(".", { autoescape: false });
/**
- * ...docs go here...
+ * The code currently has processors for three types of chunks.
*/
-async function convertMarkDown(pathdata, localeStrings, markdown) {
- try {
- markdown = await preprocessGraphicsElement(pathdata, localeStrings, markdown);
- } catch (e) {
- console.error(`Unrecoverable error in ${pathdata.file}.`);
- console.error(e);
- process.exit(1);
+const Processors = [
+ // Classic LaTeX, but only using the modern block syntax.
+ // Ancient $$ is not acceptable.
+ class LaTeXProcessor {
+ static marks = { start: `\\[`, end: `\\]` };
+ static async process(chunk, chunkPos, pathdata, localeStrings) {
+ return await latexToSVG(chunk, pathdata, localeStrings, chunkPos);
+ }
+ },
+ // Rewrite custom element stubs to their full form
+ class GraphicsElementProcessor {
+ static marks = { start: `