-
- 
-
- 
Preface
-- In order to draw things in 2D, we usually rely on lines, which typically get classified into two categories: straight lines, and curves. The - first of these are as easy to draw as they are easy to make a computer draw. Give a computer the first and last point in the line, and BAM! - straight line. No questions asked. -
-- Curves, however, are a much bigger problem. While we can draw curves with ridiculous ease freehand, computers are a bit handicapped in that - they can't draw curves unless there is a mathematical function that describes how it should be drawn. In fact, they even need this for - straight lines, but the function is ridiculously easy, so we tend to ignore that as far as computers are concerned; all lines are - "functions", regardless of whether they're straight or curves. However, that does mean that we need to come up with fast-to-compute - functions that lead to nice looking curves on a computer. There are a number of these, and in this article we'll focus on a particular - function that has received quite a bit of attention and is used in pretty much anything that can draw curves: Bézier curves. -
-- They're named after Pierre Bézier, who is principally responsible for making - them known to the world as a curve well-suited for design work (publishing his investigations in 1962 while working for Renault), although - he was not the first, or only one, to "invent" these type of curves. One might be tempted to say that the mathematician - Paul de Casteljau was first, as he began investigating the nature of these - curves in 1959 while working at Citroën, and came up with a really elegant way of figuring out how to draw them. However, de Casteljau did - not publish his work, making the question "who was first" hard to answer in any absolute sense. Or is it? Bézier curves are, at their core, - "Bernstein polynomials", a family of mathematical functions investigated by - Sergei Natanovich Bernstein, whose publications on them date back at - least as far as 1912. -
-- Anyway, that's mostly trivia, what you are more likely to care about is that these curves are handy: you can link up multiple Bézier curves - so that the combination looks like a single curve. If you've ever drawn Photoshop "paths" or worked with vector drawing programs like Flash, - Illustrator or Inkscape, those curves you've been drawing are Bézier curves. -
-- But what if you need to program them yourself? What are the pitfalls? How do you draw them? What are the bounding boxes, how do you - determine intersections, how can you extrude a curve, in short: how do you do everything that you might want to do with these curves? That's - what this page is for. Prepare to be mathed! -
-Virtually all Bézier graphics are interactive.
-- This page uses interactive examples, relying heavily on Bezier.js, as well as maths - formulae which are typeset into SVG using the XeLaTeX typesetting system and - pdf2svg by David Barton. -
-This book is open source.
-- This book is an open source software project, and lives on two github repositories. The first is - https://github.com/pomax/bezierinfo and is the purely-for-presentation version you are - viewing right now. The other repository is https://github.com/pomax/BezierInfo-2, - which is the development version, housing all the code that gets turned into the web version, and is also where you should file - issues if you find bugs or have ideas on what to change or add to the primer. -
-How complicated is the maths going to be?
-- Most of the mathematics in this Primer are early high school maths. If you understand basic arithmetic, and you know how to read English, - you should be able to get by just fine. There will at times be far more complicated maths, but if you don't feel like digesting - them, you can safely skip over them by either skipping over the "detail boxes" in section or by just jumping to the end of a section with - maths that looks too involving. The end of sections typically simply list the conclusions so you can just work with those values directly. -
-What language is all this example code in?
-- There are way too many programming languages to favour one of all others, soo all the example code in this Primer uses a form of - pseudo-code that uses a syntax that's close enough to, but not actually, modern scripting languages like JS, Python, etc. That means you - won't be able to copy-paste any of it without giving it any thought, but that's intentional: if you're reading this primer, presumably you - want to learn, and you don't learn by copy-pasting. You learn by doing things yourself, making mistakes, and then fixing - those mistakes. Now, of course, I didn't intentionally add errors in the example code just to trick you into making mistakes (that would - be horrible!) but I did intentionally keep the code from favouring one programming language over another. Don't worry though, if - you know even a single procedural programming language, you should be able to read the examples without any difficulties. -
-Questions, comments:
-- If you have suggestions for new sections, hit up the Github issue tracker (also - reachable from the repo linked to in the upper right). If you have questions about the material, there's currently no comment section - while I'm doing the rewrite, but you can use the issue tracker for that as well. Once the rewrite is done, I'll add a general comment - section back in, and maybe a more topical "select this section of text and hit the 'question' button to ask a question about it" system. - We'll see. -
-Help support the book!
-- If you enjoyed this book, or you simply found it useful for something you were trying to get done, and you were wondering how to let me - know you appreciated this book, you have two options: you can either head on over to the - Patreon page for this book, or if you prefer to make a one-time donation, head on over to - the buy Pomax a coffee page. This - work has grown from a small primer to a 100-plus print-page-equivalent reader on the subject of Bézier curves over the years, and a lot of - coffee went into the making of it. I don't regret a minute I spent on writing it, but I can always do with some more coffee to keep on - writing! -
-Що нового?
-- Цей підручник постійно розвививається, тож залежно від того, коли ви востаннє його переглядали, тут можуть бути оновлення. Перейдіть за - цим посиланням, щоб побачити, що було додано. (Також доступний RSS-канал) -
- - - -August-September 2020
--
-
-
-
- Completely overhauled the site: the Primer is now a normal web page that works fine with JS disabled, but obviously better with JS - turned on. -
-
-
June 2020
--
-
Added automatic CI/CD using Github Actions
-
January 2020
--
-
Added reset buttons to all graphics
- Updated to preface to correctly describe the on-page maths
- Fixed the Catmull-Rom section because it had glaring maths errors
-
August 2019
--
-
Added a section on (plain) rational Bezier curves
- Improved the Graphic component to allow for sliders
-
December 2018
--
-
Added a section on curvature and calculating kappa.
- -
-
- Added a Patreon page! Head on over to patreon.com/bezierinfo to help support this - site! -
-
-
August 2018
--
-
Added a section on finding a curve's y, if all you have is the x coordinate.
-
July 2018
--
-
Rewrote the 3D normals section, implementing and explaining Rotation Minimising Frames.
- Updated the section on curve order raising/lowering, showing how to get a least-squares optimized lower order curve.
- -
-
(Finally) updated 'npm test' so that it automatically rebuilds when files are changed while the dev server is running.
-
-
June 2018
--
-
Added a section on direct curve fitting.
- Added source links for all graphics.
- Added this "What's new?" section.
-
April 2017
--
-
Added a section on 3d normals.
- Added live-updating for the social link buttons, so they always link to the specific section you're reading.
-
February 2017
--
-
Finished rewriting the entire codebase for localization.
-
January 2016
--
-
Added a section to explain the Bezier interval.
- Rewrote the Primer as a React application.
-
December 2015
--
-
Set up the split repository between BezierInfo-2 as development repository, and bezierinfo as live page.
- -
-
- Removed the need for client-side LaTeX parsing entirely, so the site doesn't take a full minute or more to load all the graphics. -
-
-
May 2015
--
-
Switched over to pure JS rather than Processing-through-Processing.js
- Added Cardano's algorithm for finding the roots of a cubic polynomial.
-
April 2015
--
-
Added a section on arc length approximations.
-
February 2015
--
-
Added a section on the canonical cubic Bezier form.
-
November 2014
--
-
Switched to HTTPS.
-
July 2014
--
-
Added the section on arc approximation.
-
April 2014
--
-
Added the section on Catmull-Rom fitting.
-
November 2013
--
-
Added the section on Catmull-Rom / Bezier conversion.
- Added the section on Bezier cuves as matrices.
-
April 2013
--
-
Added a section on poly-Beziers.
- Added a section on boolean shape operations.
-
March 2013
--
-
First drastic rewrite.
- Added sections on circle approximations.
- Added a section on projecting a point onto a curve.
- Added a section on tangents and normals.
- Added Legendre-Gauss numerical data tables.
-
October 2011
--
-
-
-
- First commit for the bezierinfo site, based on the pre-Primer webpage that covered - the basics of Bezier curves in HTML with Processing.js examples. -
-
-
- - Короткий вступ -
-- Давайте розпочнемо з добрих новин: криві Безьє, про які ми говоритимемо, ви зможете побачити далі на графіках. Ці криві розпочинаються у - якійсь певній точці, і закінчуються у якійсь певній точці. Їх кривизна залежить від однієї або кількох "проміжних" контрольних точок. - Зараз, оскільки всі графіки на цій сторінці інтерактивні, поекспериментуйте трохи з цими кривими. Клікніть на точку мишкою й потягніть - - так ви зможете відчути, як форма кривої змінюється в залежності від ваших дій. -
-
-
- 
-
- - Ці криві інтенсивно використовуються у системах автоматизованого проектування та виробництва (CAD/CAM), а також у програмах для графічного - дизайну, таких як Adobe Illustrator, Photoshop, Inkscape, GIMP, тощо. Також криві Безьє використовуються у графічних технологіях, таких як - масштабована векторна графіка (SVG) та шрифти OpenType (TTF/OTF). Криві Безьє використовуються багато де, тому якщо хочете дізнатись про - них більше, приготуйтесь трохи повчитися! -
-- - So what makes a Bézier Curve? -
-- Playing with the points for curves may have given you a feel for how Bézier curves behave, but what are Bézier curves, really? - There are two ways to explain what a Bézier curve is, and they turn out to be the entirely equivalent, but one of them uses complicated - maths, and the other uses really simple maths. So... let's start with the simple explanation: -
-- Bézier curves are the result of linear interpolations. That sounds - complicated but you've been doing linear interpolation since you were very young: any time you had to point at something between two other - things, you've been applying linear interpolation. It's simply "picking a point between two points". -
-- If we know the distance between those two points, and we want a new point that is, say, 20% the distance away from the first point (and - thus 80% the distance away from the second point) then we can compute that really easily: -
- -
- - 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 - those lines gives us two points, between which we can again perform linear interpolation, yielding a single point. And that point —and all - points we can form in this way for all ratios taken together— form our Bézier curve: -
-
-
- 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 - of the fascinating parts about Bézier curves is that they can both be described in terms of polynomial functions, as well as in terms of - very simple interpolations of interpolations of [...]. That, in turn, means we can look at what these curves can do based on both "real - maths" (by examining the functions, their derivatives, and all that stuff), as well as by looking at the "mechanical" composition (which - tells us, for instance, that a curve will never extend beyond the points we used to construct it). -
-- So let's start looking at Bézier curves a bit more in depth: their mathematical expressions, the properties we can derive from them, and - the various things we can do to, and with, Bézier curves. -
-- - The mathematics of Bézier curves -
-- Bézier curves are a form of "parametric" function. Mathematically speaking, parametric functions are cheats: a "function" is actually a - well defined term representing a mapping from any number of inputs to a single output. Numbers go in, a single number - comes out. Change the numbers that go in, and the number that comes out is still a single number. -
-- Parametric functions cheat. They basically say "alright, well, we want multiple values coming out, so we'll just use more than one - 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 - function. Parametric functions cheat by changing that. In a parametric function all the different functions share a variable, like this: -
- -
- - Multiple functions, but only one variable. If we change the value for t, we change the outcome of both fa(t) and - fb(t). You might wonder how that's useful, and the answer is actually pretty simple: if we change the labels - 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 - link the values to a "control" variable. If we vary the value of t, then with every change we get two values, - which we can use as (x,y) coordinates in a graph. The above set of functions, for instance, generates points on a circle: We - can range t from negative to positive infinity, and the resulting (x,y) coordinates will always lie on a circle with - radius 1 around the origin (0,0). If we plot it for t from 0 to 5, we get this: -
-
-
- - 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 - sine, the other a cosine); but Bézier curves use the "binomial polynomial" for both the x and y outputs. So what are - binomial polynomials? -
-You may remember polynomials from high school. They're those sums that look like this:
- -
- - If the highest order term they have is x³, they're called "cubic" polynomials; if it's x², it's a "square" polynomial; if - it's just x, it's a line (and if there aren't even any terms with x it's not a polynomial!) -
-- Bézier curves are polynomials of t, rather than x, with the value for t being fixed between 0 and 1, with - coefficients a, b etc. taking the "binomial" form, which sounds fancy but is actually a pretty simple description for mixing - values: -
- -
- - I know what you're thinking: that doesn't look too simple! But if we remove t and add in "times one", things suddenly look pretty - easy. Check out these binomial terms: -
- -
- - Notice that 2 is the same as 1+1, and 3 is 2+1 and 1+2, and 6 is 3+3... As you can see, each time we go up a dimension, we simply start - and end with 1, and everything in between is just "the two numbers above it, added together", giving us a simple number sequence known as - Pascal's triangle. Now that's easy to remember. -
-- There's an equally simple way to figure out how the polynomial terms work: if we rename (1-t) to a and t to b, - and remove the weights for a moment, we get this: -
- -
- - It's basically just a sum of "every combination of a and b", progressively replacing a's with b's after every - + sign. So that's actually pretty simple too. So now you know binomial polynomials, and just for completeness I'm going to show you the - generic function for this: -
- -
- - And that's the full description for Bézier curves. Σ in this function indicates that this is a series of additions (using the variable - listed below the Σ, starting at ...=<value> and ending at the value listed on top of the Σ). -
-How to implement the basis function
-We could naively implement the basis function as a mathematical construct, using the function as our guide, like this:
- -| 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 - [1,3,3,1], and so on, with each next row fitting 1 more number than the previous row, starting and ending with "1", with all the numbers - in between being the sum of the previous row's elements on either side "above" the one we're computing. -
-- 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: -
- -| 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: -
- -| 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: -
- -| 1 | -- - | -
| 2 | -|
| 3 | -|
| 4 | -|
| 5 | -|
| 6 | -|
| 7 | -|
| 8 | -|
| 9 | -|
| 10 | -|
| 11 | -|
| 12 | -|
| 13 | -
And now we know how to program the basis function. Excellent.
-So, now we know what the basis function looks like, time to add in the magic that makes Bézier curves so special: control points.
-- - Controlling Bézier curvatures -
-- Bézier curves are, like all "splines", interpolation functions. This means that they take a set of points, and generate values somewhere - "between" those points. (One of the consequences of this is that you'll never be able to generate a point that lies outside the outline - for the control points, commonly called the "hull" for the curve. Useful information!). In fact, we can visualize how each point - contributes to the value generated by the function, so we can see which points are important, where, in the curve. -
-- The following graphs show the interpolation functions for quadratic and cubic curves, with "S" being the strength of a point's - contribution to the total sum of the Bézier function. Click-and-drag to see the interpolation percentages for each curve-defining point at - a specific t value. -
-
-
- 
-
- 
-
- - Also shown is the interpolation function for a 15th order Bézier function. As you can see, the start and end point contribute - considerably more to the curve's shape than any other point in the control point set. -
-- If we want to change the curve, we need to change the weights of each point, effectively changing the interpolations. The way to do this - 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. -
-How to implement the weighted basis function
-Given that we already know how to implement basis function, adding in the control points is remarkably easy:
- -| 1 | -- - | -
| 2 | -|
| 3 | -|
| 4 | -|
| 5 | -
And now for the optimized versions:
- -| 1 | -- - | -
| 2 | -|
| 3 | -|
| 4 | -|
| 5 | -|
| 6 | -|
| 7 | -|
| 8 | -|
| 9 | -|
| 10 | -|
| 11 | -|
| 12 | -|
| 13 | -
And now we know how to program the weighted basis function.
-- - Controlling Bézier curvatures, part 2: Rational Béziers -
-- We can further control Bézier curves by "rationalising" them: that is, adding a "ratio" value in addition to the weight value discussed in - the previous section, thereby gaining control over "how strongly" each coordinate influences the curve. -
-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:
- -
-
- In this, the first new term represents an additional weight for each coordinate. For example, if our ratio values are [1, 0.5, 0.5, 1]
- then ratio0 = 1, ratio1 = 0.5, and so on, and is effectively identical as if we were just
- using different weight. So far, nothing too special.
-
- However, the second new term is what makes the difference: every point on the curve isn't just a "double weighted" point, it is a - fraction of the "doubly weighted" value we compute by introducing that ratio. When computing points on the curve, we compute the - "normal" Bézier value and then divide that by the Bézier value for the curve that only uses ratios, not weights. -
-- This does something unexpected: it turns our polynomial into something that isn't a polynomial anymore. It is now a kind of curve - that is a super class of the polynomials, and can do some really cool things that Bézier curves can't do "on their own", such as perfectly - describing circles (which we'll see in a later section is literally impossible using standard Bézier curves). -
-- But the best way to show what this does is to do literally that: let's look at the effect of "rationalising" our Bézier curves using an - interactive graphic for a rationalised curves. The following graphic shows the Bézier curve from the previous section, "enriched" with - ratio factors for each coordinate. The closer to zero we set one or more terms, the less relative influence the associated coordinate - exerts on the curve (and of course the higher we set them, the more influence they have). Try to change the values and see how it affects - what gets drawn: -
-
-
- - 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 - with gravity, if the relative strengths stay the same, nothing really changes. The values define each coordinate's influence - relative to all other points. -
-How to implement rational curves
-Extending the code of the previous section to include ratios is almost trivial:
- -| 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 Bézier interval [0,1] -
-
- Now that we know the mathematics behind Bézier curves, there's one curious thing that you may have noticed: they always run from
- t=0 to t=1. Why that particular interval?
-
- It all has to do with how we run from "the start" of our curve to "the end" of our curve. If we have a value that is a mixture of two - other values, then the general formula for this is: -
- -
-
- 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:
- if they are, then we could just pick whatever values we like, and end up with a mixed value that is, for example, 100% value 1
- and 100% value 2. In principle that's fine, but for Bézier curves we always want mixed values between the start
- 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.
-
- But... what if we use this form, which is based on the assumption that we will only ever use values between 0 and 1, and instead use - values outside of that interval? Do things go horribly wrong? Well... not really, but we get to "see more". -
-- In the case of Bézier curves, extending the interval simply makes our curve "keep going". Bézier curves are simply segments of some - polynomial curve, so if we pick a wider interval we simply get to see more of the curve. So what do they look like? -
-
- The following two graphics show you Bézier curves rendered "the usual way", as well as the curves they "lie on" if we were to extend the
- t values much further. As you can see, there's a lot more "shape" hidden in the rest of the curve, and we can model those
- parts by moving the curve points around.
-
-
- 
-
- - In fact, there are curves used in graphics design and computer modelling that do the opposite of Bézier curves; rather than fixing the - interval, and giving you freedom to choose the coordinates, they fix the coordinates, but give you freedom over the interval. A great - example of this is the "Spiro" curve, which is a curve based on part of a - Cornu Spiral, also known as Euler's Spiral. It's a very aesthetically pleasing - curve and you'll find it in quite a few graphics packages like FontForge and - Inkscape. It has even been used in font design, for example for the Inconsolata typeface. -
-- - Bézier curvatures as matrix operations -
-- We can also represent Bézier curves as matrix operations, by expressing the Bézier formula as a polynomial basis function and a - 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.
-
- So: why would we bother with matrices? Matrix representations allow us to discover things about functions that would - otherwise be hard to tell. It turns out that the curves form - triangular matrices, and they have a determinant equal to the product of the - actual coordinates we use for our curve. It's also invertible, which means there's - a ton of properties that are all satisfied. Of - course, the main question is "why is this useful to us now?", and the answer to that is that it's not immediately useful, but - you'll be seeing some instances where certain curve properties can be either computed via function manipulation, or via clever use of - matrices, and sometimes the matrix approach can be (drastically) faster. -
-So for now, just remember that we can represent curves this way, and let's move on.
-- - de Casteljau's algorithm -
-
- If we want to draw Bézier curves, we can run through all values of t from 0 to 1 and then compute the weighted basis function
- at each value, getting the x/y values we need to plot. Unfortunately, the more complex the curve gets, the more expensive
- this computation becomes. Instead, we can use de Casteljau's algorithm to draw curves. This is a geometric approach to curve
- drawing, and it's really easy to implement. So easy, in fact, you can do it by hand with a pencil and ruler.
-
Rather than using our calculus function to find x/y values for t, let's do this instead:
-
-
- treat
tas a ratio (which it is). t=0 is 0% along a line, t=1 is 100% along a line.
- - Take all lines between the curve's defining points. For an order
ncurve, that'snlines.
- -
- Place markers along each of these line, at distance
t. So iftis 0.2, place the mark at 20% from the start, - 80% from the end. -
- - Now form lines between
thosepoints. This givesn-1lines.
- - Place markers along each of these line at distance
t.
- - Form lines between
thosepoints. This'll ben-2lines.
- - Place markers, form lines, place markers, etc. -
-
- Repeat this until you have only one line left. The point
ton that line coincides with the original curve point at -t. -
-
- To see this in action, mouse-over the following sketch. Moving the mouse changes which curve point is explicitly evaluated using de - Casteljau's algorithm, moving the cursor left-to-right (or, of course, right-to-left), shows you how a curve is generated using this - approach. -
-
-
- How to implement de Casteljau's algorithm
-Let's just use the algorithm we just specified, and implement that:
- -| 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:
-
| 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 - new list. -
-- - Simplified drawing -
-- We can also simplify the drawing process by "sampling" the curve at certain points, and then joining those points up with straight lines, - a process known as "flattening", as we are reducing a curve to a simple sequence of straight, "flat" lines. -
-- We can do this is by saying "we want X segments", and then sampling the curve at intervals that are spaced such that we end up with the - number of segments we wanted. The advantage of this method is that it's fast: instead of evaluating 100 or even 1000 curve coordinates, we - can sample a much lower number and still end up with a curve that sort-of-kind-of looks good enough. The disadvantage of course is that we - lose the precision of working with "the real curve", so we usually can't use the flattened for doing true intersection detection, or - curvature alignment. -
-
-
- 
-
- - Try clicking on the sketch and using your up and down arrow keys to lower the number of segments for both the quadratic and cubic curve. - You'll notice that for certain curvatures, a low number of segments works quite well, but for more complex curvatures (try this for the - cubic curve), a higher number is required to capture the curvature changes properly. -
-How to implement curve flattening
-Let's just use the algorithm we just specified, and implement that:
- -| 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:
- -| 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.
-- - Splitting curves -
-
- Using de Casteljau's algorithm, we can also find all the points we need to split up a Bézier curve into two, smaller curves, which taken
- together form the original curve. When we construct de Casteljau's skeleton for some value t, the procedure gives us all the
- points we need to split a curve at that t value: one curve is defined by all the inside skeleton points found prior to our
- on-curve point, with the other curve being defined by all the inside skeleton points after our on-curve point.
-
-
- implementing curve splitting
-We can implement curve splitting by bolting some extra logging onto the de Casteljau function:
- -| 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
- same order as the original curve, and can be overlaid exactly on the original curve.
-
- - Splitting curves using matrices -
-- Another way to split curves is to exploit the matrix representation of a Bézier curve. In the section on matrices, - we saw that we can represent curves as matrix multiplications. Specifically, we saw these two forms for the quadratic and cubic curves - respectively: (we'll reverse the Bézier coefficients vector for legibility) -
- -
- and
- -
-
- Let's say we want to split the curve at some point t = z, forming two new (obviously smaller) Bézier curves. To find the
- coordinates for these two Bézier curves, we can use the matrix representation and some linear algebra. First, we separate out the actual
- "point on the curve" information into a new matrix multiplication:
-
- and
- -
-
- If we could compact these matrices back to the form [t values] · [Bézier matrix] · [column matrix], with the first two
- staying the same, then that column matrix on the right would be the coordinates of a new Bézier curve that describes the first segment,
- from t = 0 to t = z. As it turns out, we can do this quite easily, by exploiting some simple rules of linear
- algebra (and if you don't care about the derivations, just skip to the end of the box for the results!).
-
Deriving new hull coordinates
-- Deriving the two segments upon splitting a curve takes a few steps, and the higher the curve order, the more work it is, so let's look - at the quadratic curve first: -
- -
-
- 
-
- 
-
- 
- - We can do this because [M · M-1] is the identity matrix. It's a bit like multiplying something by x/x in calculus: it doesn't do anything to the function, but it does - allow you to rewrite it to something that may be easier to work with, or can be broken up differently. In the same way, multiplying our - matrix by [M · M-1] has no effect on the total formula, but it does allow us to change the matrix sequence [something · M] to a sequence [M · something], and that makes a world of difference: if we know what [M-1 · Z · M] is, we can apply that to our coordinates, - and be left with a proper matrix representation of a quadratic Bézier curve (which is [T · M · P]), with a new set of - coordinates that represent the curve from t = 0 to t = z. So let's get computing: -
- -
- Excellent! Now we can form our new quadratic curve:
- -
-
- 
-
- 
-
- Brilliant: if we want a subcurve from t = 0 to t = z, we can keep the first coordinate the same (which makes sense),
- our control point becomes a z-ratio mixture of the original control point and the start point, and the new end point is a mixture that
- looks oddly similar to a Bernstein polynomial of degree two. These new
- coordinates are actually really easy to compute directly!
-
- Of course, that's only one of the two curves. Getting the section from t = z to t = 1 requires doing this
- again. We first observe that in the previous calculation, we actually evaluated the general interval [0,z]. We were able to
- write it down in a more simple form because of the zero, but what we actually evaluated, making the zero explicit, was:
-
-
- 
- If we want the interval [z,1], we will be evaluating this instead:
- -
-
- 
-
- We're going to do the same trick of multiplying by the identity matrix, to turn [something · M] into
- [M · something]:
-
- So, our final second curve looks like:
- -
-
- 
-
- 
- - Nice. We see the same thing as before: we can keep the last coordinate the same (which makes sense); our control point becomes a z-ratio - mixture of the original control point and the end point, and the new start point is a mixture that looks oddly similar to a bernstein - polynomial of degree two, except this time it uses (z-1) rather than (1-z). These new coordinates are also really easy to - compute directly! -
-
- So, using linear algebra rather than de Casteljau's algorithm, we have determined that, for any quadratic curve split at some value
- t = z, we get two subcurves that are described as Bézier curves with simple-to-derive coordinates:
-
- and
- -
- - We can do the same for cubic curves. However, I'll spare you the actual derivation (don't let that stop you from writing that out - yourself, though) and simply show you the resulting new coordinate sets: -
- -
- and
- -
- - So, looking at our matrices, did we really need to compute the second segment matrix? No, we didn't. Actually having one segment's matrix - means we implicitly have the other: push the values of each row in the matrix Q to the right, with zeroes - getting pushed off the right edge and appearing back on the left, and then flip the matrix vertically. Presto, you just "calculated" - Q'. -
-- Implementing curve splitting this way requires less recursion, and is just straight arithmetic with cached values, so can be cheaper on - systems where recursion is expensive. If you're doing computation with devices that are good at matrix multiplication, chopping up a - Bézier curve with this method will be a lot faster than applying de Casteljau. -
-- - Lowering and elevating curve order -
-- One interesting property of Bézier curves is that an nth order curve can always be perfectly represented by an - (n+1)th order curve, by giving the higher-order curve specific control points. -
-- If we have a curve with three points, then we can create a curve with four points that exactly reproduces the original curve. First, we - give it the same start and end points, and for its two control points we pick "1/3rd start + 2/3rd control" and - "2/3rd control + 1/3rd end". Now we have exactly the same curve as before, except represented as a cubic curve - rather than a quadratic curve. -
-- The general rule for raising an 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 - try to, but the resulting curve will not be identical to the original, and may in fact look completely different. -
-- However, there is a surprisingly good way to ensure that a lower order curve looks "as close as reasonably possible" to the original - curve: we can optimise the "least-squares distance" between the original curve and the lower order curve, in a single operation (also - explained over on Sirver's Castle). However, to - use it, we'll need to do some calculus work and then switch over to linear algebra. As mentioned in the section on matrix representations, - 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):
-
-
- So by using this seemingly silly trick, we can suddenly express part of our nth order Bézier function in terms of an (n+1)th
- order Bézier function. And that sounds a lot like raising the curve order! Of course we need to be able to repeat that trick for the
- t part, but that's not a problem:
-
-
- So, with both of those changed from an order n expression to an order (n+1) expression, we can put them back
- together again. Now, where the order n function had a summation from 0 to n, the order n+1 function
- uses a summation from 0 to n+1, but this shouldn't be a problem as long as we can add some new terms that "contribute
- nothing". In the next section on derivatives, there is a discussion about why "higher terms than there is a binomial for" and "lower than
- zero terms" both "contribute nothing". So as long as we can add terms that have the same form as the terms we need, we can just include
- them in the summation, they'll sit there and do nothing, and the resulting function stays identical to the lower order curve.
-
Let's do this:
- -
- - 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:
- - That might look unwieldy, but it's really just a mostly-zeroes matrix, with a very simply fraction on the diagonal, and an even simpler - fraction to the left of it. Multiplying a list of coordinates with this matrix means we can plug the resulting transformed coordinates - into the one-order-higher function and get an identical looking curve. -
-Not too bad!
-- Equally interesting, though, is that with this matrix operation established, we can now use an incredibly powerful and ridiculously simple - way to find out a "best fit" way to reverse the operation, called - the normal equation. What it does is minimize the sum of the square - differences between one set of values and another set of values. Specifically, if we can express that as some function - A x = b, we can use it. And as it so happens, that's exactly what we're dealing with, so: -
- -
- The steps taken here are:
--
-
- We have a function in a form that the normal equation can be used with, so -
- apply the normal equation! -
- - Then, we want to end up with just Bn on the left, so we start by left-multiply both sides such that we'll end up with lots of - stuff on the left that simplified to "a factor 1", which in matrix maths is the - identity matrix. - -
- - In fact, by left-multiplying with the inverse of what was already there, we've effectively "nullified" (but really, one-inified) that - big, unwieldy block into the identity matrix I. So we substitute the mess with I, and then - -
- - because multiplication with the identity matrix does nothing (like multiplying by 1 does nothing in regular algebra), we just drop it. - -
- And we're done: we now have an expression that lets us approximate an n+1th order curve with a lower nth order curve. It won't be an exact fit, but it's definitely a best approximation. So, let's implement these rules for
- raising and lowering curve order to a (semi) random curve, using the following graphic. Select the sketch, which has movable control
- points, and press your up and down arrow keys to raise or lower the curve order.
-
-
- - - Derivatives -
-- There's a number of useful things that you can do with Bézier curves based on their derivative, and one of the more amusing observations - about Bézier curves is that their derivatives are, in fact, also Bézier curves. In fact, the differentiation of a Bézier curve is - relatively straightforward, although we do need a bit of math. -
-First, let's look at the derivative rule for Bézier curves, which is:
- -
- - which we can also write (observing that b in this formula is the same as our w weights, and that n times a summation - is the same as a summation where each term is multiplied by n) as: -
- -
- - Or, in plain text: the derivative of an nth degree Bézier curve is an (n-1)th degree Bézier curve, with one fewer - term, and new weights w'0...w'n-1 derived from the original weights as n(wi+1 - wi). So for a - 3rd degree curve, with four weights, the derivative has three new weights: w'0 = 3(w1-w0), - w'1 = 3(w2-w1) and w'2 = 3(w3-w2). -
-"Slow down, why is that true?"
-- Sometimes just being told "this is the derivative" is nice, but you might want to see why this is indeed the case. As such, let's have a - look at the proof for this derivative. First off, the weights are independent of the full Bézier function, so the derivative involves - only the derivative of the polynomial basis function. So, let's find that: -
- -
- - Applying the product and - chain rules gives us: -
- -
- Which is hard to work with, so let's expand that properly:
- -
- - Now, the trick is to turn this expression into something that has binomial coefficients again, so we want to end up with things that - look like "x! over y!(x-y)!". If we can do that in a way that involves terms of n-1 and k-1, we'll be on the right track. -
- -
- And that's the first part done: the two components inside the parentheses are actually regular, lower-order Bézier expressions:
- -
- - Now to apply this to our weighted Bézier curves. We'll write out the plain curve formula that we saw earlier, and then work our way - through to its derivative: -
- -
- - If we expand this (with some color to show how terms line up), and reorder the terms by increasing values for k we see the - following: -
- -
- - Two of these terms fall way: the first term falls away because there is no -1st term in a summation. As such, it always - contributes "nothing", so we can safely completely ignore it for the purpose of finding the derivative function. The other term is the - very last term in this expansion: one involving Bn-1,n. This term would have a binomial coefficient of [i choose i+1], which is a non-existent binomial coefficient. Again, - this term would contribute "nothing", so we can ignore it, too. This means we're left with: -
- -
- And that's just a summation of lower order curves:
- -
- We can rewrite this as a normal summation, and we're done:
- -
- - Let's rewrite that in a form similar to our original formula, so we can see the difference. We will first list our original formula for - Bézier curves, and then the derivative: -
- -
-
- 
- - What are the differences? In terms of the actual Bézier curve, virtually nothing! We lowered the order (rather than n, it's now - n-1), but it's still the same Bézier function. The only real difference is in how the weights change when we derive the curve's - function. If we have four points A, B, C, and D, then the derivative will have three points, the second derivative two, and the third - derivative one: -
- -
- - We can keep performing this trick for as long as we have more than one weight. Once we have one weight left, the next step will see - k = 0, and the result of our "Bézier function" summation is zero, because we're not adding anything at all. As such, a quadratic - curve has no second derivative, a cubic curve has no third derivative, and generalized: an nth order curve has - n-1 (meaningful) derivatives, with any further derivative being zero. -
-- - Tangents and normals -
-- If you want to move objects along a curve, or "away from" a curve, the two vectors you're most interested in are the tangent vector and - normal vector for curve points. These are actually really easy to find. For moving and orienting along a curve, we use the tangent, which - indicates the direction of travel at specific points, and is literally just the first derivative of our curve: -
- -
- - This gives us the directional vector we want. We can normalize it to give us uniform directional vectors (having a length of 1.0) at each - point, and then do whatever it is we want to do based on those directions: -
- -
- - The tangent is very useful for moving along a line, but what if we want to move away from the curve instead, perpendicular to the curve at - some point t? In that case we want the normal vector. This vector runs at a right angle to the direction of the curve, and - is typically of length 1.0, so all we have to do is rotate the normalized directional vector and we're done: -
- -
- - Rotating coordinates is actually very easy, if you know the rule for it. You might find it explained as "applying a - rotation matrix, which is what we'll look at here, too. Essentially, the - idea is to take the circles over which we can rotate, and simply "sliding the coordinates" over these circles by the desired angle. If - we want a quarter circle turn, we take the coordinate, slide it along the circle by a quarter turn, and done. -
-- To turn any point (x,y) into a rotated point (x',y') (over 0,0) by some angle φ, we apply this nice and easy computation: -
- -
- Which is the "long" version of the following matrix transformation:
- -
- - And that's all we need to rotate any coordinate. Note that for quarter, half, and three-quarter turns these functions become even - easier, since sin and cos for these angles are, respectively: 0 and 1, -1 and 0, and 0 and -1. -
-- But why does this work? Why this matrix multiplication? - Wikipedia (technically, Thomas Herter and Klaus - Lott) tells us that a rotation matrix can be treated as a sequence of three (elementary) shear operations. When we combine this into a - single matrix operation (because all matrix multiplications can be collapsed), we get the matrix that you see above. - DataGenetics have an excellent article about this very thing: it's - really quite cool, and I strongly recommend taking a quick break from this primer to read that article. -
-- The following two graphics show the tangent and normal along a quadratic and cubic curve, with the direction vector coloured blue, and the - normal vector coloured red (the markers are spaced out evenly as t-intervals, not spaced equidistant). -
-
-
- 
-
- - - Working with 3D normals -
-- Before we move on to the next section we need to spend a little bit of time on the difference between 2D and 3D. While for many things - this difference is irrelevant and the procedures are identical (for instance, getting the 3D tangent is just doing what we do for 2D, but - for x, y, and z, instead of just for x and y), when it comes to normals things are a little more complex, and thus more work. Mind you, - it's not "super hard", but there are more steps involved and we should have a look at those. -
-- Getting normals in 3D is in principle the same as in 2D: we take the normalised tangent vector, and then rotate it by a quarter turn. - However, this is where things get that little more complex: we can turn in quite a few directions, since "the normal" in 3D is a plane, - not a single vector, so we basically need to define what "the" normal is in the 3D case. -
-- The "naïve" approach is to construct what is known as the - Frenet normal, where we follow a simple recipe that works in - many cases (but does super bizarre things in some others). The idea is that even though there are infinitely many vectors that are - perpendicular to the tangent (i.e. make a 90 degree angle with it), the tangent itself sort of lies on its own plane already: since each - point on the curve (no matter how closely spaced) has its own tangent vector, we can say that each point lies in the same plane as the - local tangent, as well as the tangents "right next to it". -
-- Even if that difference in tangent vectors is minute, "any difference" is all we need to find out what that plane is - or rather, what the - vector perpendicular to that plane is. Which is what we need: if we can calculate that vector, and we have the tangent vector that we know - lies on a plane, then we can rotate the tangent vector over the perpendicular, and presto. We have computed the normal using the same - logic we used for the 2D case: "just rotate it 90 degrees". -
-So let's do that! And in a twist surprise, we can do this in four lines:
--
-
- a = normalize(B'(t)) -
- b = normalize(a + B''(t)) -
- r = normalize(b × a) -
- normal = normalize(r × a) -
Let's unpack that a little:
--
-
- - We start by taking the normalized vector for the derivative at some point on the - curve. We normalize it so the maths is less work. Less work is good. - -
- - Then, we compute b which represents what a next point's tangent would be if the curve stopped changing at our point and - just had the same derivative and second derivative from that point on. - -
- - This lets us find two vectors (the derivative, and the second derivative added to the derivative) that lie on the same plane, which - means we can use them to compute a vector perpendicular to that plane, using an elementary vector operation called the - cross product. (Note that while that operation uses the × operator, it's most - definitely not a multiplication!) The result of that gives us a vector that we can use as the "axis of rotation" for turning the tangent - a quarter circle to get our normal, just like we did in the 2D case. - -
- - Since the cross product lets us find a vector that is perpendicular to some plane defined by two other vectors, and since the normal - vector should be perpendicular to the plane that the tangent and the axis of rotation lie in, we can use the cross product a second - time, and immediately get our normal vector. - -
- And then we're done, we found "the" normal vector for a 3D curve. Let's see what that looks like for a sample curve, shall we? You can - move your cursor across the graphic from left to right, to show the normal at a point with a t value that is based on your cursor - position: all the way on the left is 0, all the way on the right = 1, midway is t=0.5, etc: -
-
-
- - 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? -
-- As it turns out, it's doing that because that's how the maths works, and that's the problem with Frenet normals: while they are - "mathematically correct", they are "practically problematic", and so for any kind of graphics work what we really want is a way to compute - normals that just... look good. -
-Thankfully, Frenet normals are not our only option.
-- Another option is to take a slightly more algorithmic approach and compute a form of - Rotation Minimising Frame - (also known as "parallel transport frame" or "Bishop frame") instead, where a "frame" is a set made up of the tangent, the rotational - axis, and the normal vector, centered on an on-curve point. -
-- These type of frames are computed based on "the previous frame", so we cannot simply compute these "on demand" for single points, as we - could for Frenet frames; we have to compute them for the entire curve. Thankfully, the procedure is pretty simple, and can be performed at - the same time that you're building lookup tables for your curve. -
-- The idea is to take a starting "tangent/rotation axis/normal" frame at t=0, and then compute what the next frame "should" look like by - applying some rules that yield a good looking next frame. In the case of the RMF paper linked above, those rules are: -
--
-
- Take a point on the curve for which we know the RM frame already, -
- take a next point on the curve for which we don't know the RM frame yet, and -
- - reflect the known frame onto the next point, by treating the plane through the curve at the point exactly between the next and previous - points as a "mirror". - -
- - This gives the next point a tangent vector that's essentially pointing in the opposite direction of what it should be, and a normal - that's slightly off-kilter, so: - -
- - reflect the vectors of our "mirrored frame" a second time, but this time using the plane through the "next point" itself as "mirror". - -
- Done: the tangent and normal have been fixed, and we have a good looking frame to work with. -
So, let's write some code for that!
-Implementing Rotation Minimising Frames
-- 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:
- -| 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. -
-- Speaking of better looking, what does this actually look like? Let's revisit that earlier curve, but this time use rotation minimising - frames rather than Frenet frames: -
-
-
- 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 - the z-axis as our initial axis of rotation, so that our initial normal is (0,0,1) × tangent, and then take things from - there, but having that initial "mathematically correct" frame so that the initial normal seems to line up based on the curve's orientation - in 3D space is just nice. -
-- - Component functions -
-- One of the first things people run into when they start using Bézier curves in their own programs is "I know how to draw the curve, but - how do I determine the bounding box?". It's actually reasonably straightforward to do so, but it requires having some knowledge on - exploiting math to get the values we need. For bounding boxes, we aren't actually interested in the curve itself, but only in its - "extremities": the minimum and maximum values the curve has for its x- and y-axis values. If you remember your calculus (provided you ever - took calculus, otherwise it's going to be hard to remember) we can determine function extremities using the first derivative of that - function, but this poses a problem, since our function is parametric: every axis has its own function. -
-- The solution: compute the derivative for each axis separately, and then fit them back together in the same way we do for the original. -
-- Let's look at how a parametric Bézier curve "splits up" into two normal functions, one for the x-axis and one for the y-axis. Note the - leftmost figure is again an interactive curve, without labeled axes (you get coordinates in the graph instead). The center and rightmost - figures are the component functions for computing the x-axis value, given a value for t (between 0 and 1 inclusive), and the y-axis - value, respectively. -
-- If you move points in a curve sideways, you should only see the middle graph change; likewise, moving points vertically should only show a - change in the right graph. -
-
- -
- - - Finding extremities: root finding -
-- Now that we understand (well, superficially anyway) the component functions, we can find the extremities of our Bézier curve by finding - maxima and minima on the component functions, by solving the equation B'(t) = 0. We've already seen that the derivative of a Bézier curve - is a simpler Bézier curve, but how do we solve the equality? Fairly easily, actually, until our derivatives are 4th order or higher... - then things get really hard. But let's start simple: -
-Quadratic curves: linear derivatives.
-
- The derivative of a quadratic Bézier curve is a linear Bézier curve, interpolating between just two terms, which means finding the
- solution for "where is this line 0" is effectively trivial by rewriting it to a function of t and solving. First we turn our
- cubic Bézier function into a quadratic one, by following the rule mentioned at the end of the
- derivatives section:
-
- And then we turn this into our solution for t using basic arithmetics:
- Done.
-
- Although with the caveat that if b-a is zero, there
- is no solution and we probably shouldn't try to perform that division.
-
Cubic curves: the quadratic formula.
-- The derivative of a cubic Bézier curve is a quadratic Bézier curve, and finding the roots for a quadratic polynomial means we can apply - the Quadratic formula. If you've seen it before, you'll remember it, and if - you haven't, it looks like this: -
- -
- - So, if we can rewrite the Bézier component function as a plain polynomial, we're done: we just plug in the values into the quadratic - formula, check if that square root is negative or not (if it is, there are no roots) and then just compute the two values that come out - (because of that plus/minus sign we get two). Any value between 0 and 1 is a root that matters for Bézier curves, anything below or above - that is irrelevant (because Bézier curves are only defined over the interval [0,1]). So, how do we convert? -
-- First we turn our cubic Bézier function into a quadratic one, by following the rule mentioned at the end of the - derivatives section: -
- -
- And then, using these v values, we can find out what our a, b, and c should be:
- -
-
- This gives us three coefficients {a, b, c} that are expressed in terms of v values, where the v values are
- expressions of our original coordinate values, so we can do some substitution to get:
-
- Easy-peasy. We can now almost trivially find the roots by plugging those values into the quadratic formula.
-- And as a cubic curve, there is also a meaningful second derivative, which we can compute by simple taking the derivative of the - derivative. -
-Quartic curves: Cardano's algorithm.
-- We haven't really looked at them before now, but the next step up would be a Quartic curve, a fourth degree Bézier curve. As expected, - these have a derivative that is a cubic function, and now things get much harder. Cubic functions don't have a "simple" rule to find their - roots, like the quadratic formula, and instead require quite a bit of rewriting to a form that we can even start to try to solve. -
-- Back in the 16th century, before Bézier curves were a thing, and even before calculus itself was a thing, - Gerolamo Cardano figured out that even if the general cubic function is - really hard to solve, it can be rewritten to a form for which finding the roots is "easier" (even if not "easy"): -
- -
-
- We can see that the easier formula only has two constants, rather than four, and only two expressions involving t, rather
- than three: this makes things considerably easier to solve because it lets us use
- regular calculus to find the values that satisfy the equation.
-
- Now, there is one small hitch: as a cubic function, the solutions may be - complex numbers rather than plain numbers... And Cardano realised this, - centuries before complex numbers were a well-understood and established part of number theory. His interpretation of them was "these - numbers are impossible but that's okay because they disappear again in later steps", allowing him to not think about them too much, but we - have it even easier: as we're trying to find the roots for display purposes, we don't even care about complex numbers: we're - going to simplify Cardano's approach just that tiny bit further by throwing away any solution that's not a plain number. -
-- So, how do we rewrite the hard formula into the easier formula? This is explained in detail over at - Ken J. Ward's page for solving the - cubic equation, so instead of showing the maths, I'm simply going to show the programming code for solving the cubic equation, with the - complex roots getting totally ignored, but if you're interested you should definitely head over to Ken's page and give the procedure a - read-through. -
-Implementing Cardano's algorithm for finding all real roots
-- The "real roots" part is fairly important, because while you cannot take a square, cube, etc. root of a negative number in the "real" - number space (denoted with ℝ), this is perfectly fine in the - "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! -
- -| 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 | -
- And that's it. The maths is complicated, but the code is pretty much just "follow the maths, while caching as many values as we can to - prevent recomputing things as much as possible" and now we have a way to find all roots for a cubic function and can just move on with - using that to find extremities of our curves. -
-- And of course, as a quartic curve also has meaningful second and third derivatives, we can quite easily compute those by using the - derivative of the derivative (of the derivative), just as for cubic curves. -
-Quintic and higher order curves: finding numerical solutions
-- And this is where thing stop, because we cannot find the roots for polynomials of degree 5 or higher using algebra (a fact known - as the Abel–Ruffini theorem). Instead, for occasions like these, - where algebra simply cannot yield an answer, we turn to numerical analysis. -
-- That's a fancy term for saying "rather than trying to find exact answers by manipulating symbols, find approximate answers by describing - the underlying process as a combination of steps, each of which can be assigned a number via symbolic manipulation". For example, - trying to mathematically compute how much water fits in a completely crazy three dimensional shape is very hard, even if it got you the - perfect, precise answer. A much easier approach, which would be less perfect but still entirely useful, would be to just grab a buck and - start filling the shape until it was full: just count the number of buckets of water you used. And if we want a more precise answer, we - can use smaller buckets. -
-
- So that's what we're going to do here, too: we're going to treat the problem as a sequence of steps, and the smaller we can make each
- step, the closer we'll get to that "perfect, precise" answer. And as it turns out, there is a really nice numerical root-finding
- algorithm, called the Newton-Raphson root finding method (yes, after
- that Newton), which we can make use of. The Newton-Raphson approach
- consists of taking our impossible-to-solve function f(x), picking some initial value x (literally any value will
- do), and calculating f(x). We can think of that value as the "height" of the function at x. If that height is
- zero, we're done, we have found a root. If it isn't, we calculate the tangent line at f(x) and calculate at which
- x value its height is zero (which we've already seen is very easy). That will give us a new x and we
- repeat the process until we find a root.
-
- Mathematically, this means that for some x, at step n=1, we perform the following calculation until
- fy(x) is zero, so that the next t is the same as the one we already have:
-
- (The Wikipedia article has a decent animation for this process, so I will not add a graphic for that here)
-- Now, this works well only if we can pick good starting points, and our curve is - continuously differentiable and doesn't have - oscillations. Glossing over the exact meaning of those terms, the - curves we're dealing with conform to those constraints, so as long as we pick good starting points, this will work. So the question is: - which starting points do we pick? -
-- As it turns out, Newton-Raphson is so blindingly fast that we could get away with just not picking: we simply run the algorithm from - t=0 to t=1 at small steps (say, 1/200th) and the result will be all the roots we want. Of course, this may - pose problems for high order Bézier curves: 200 steps for a 200th order Bézier curve is going to go wrong, but that's okay: - there is no reason (at least, none that I know of) to ever use Bézier curves of crazy high orders. You might use a fifth order - curve to get the "nicest still remotely workable" approximation of a full circle with a single Bézier curve, but that's pretty much as - high as you'll ever need to go. -
-In conclusion:
-- So now that we know how to do root finding, we can determine the first and second derivative roots for our Bézier curves, and show those - roots overlaid on the previous graphics. For the quadratic curve, that means just the first derivative, in red: -
-
- And for cubic curves, that means first and second derivatives, in red and purple respectively:
-
- - - Bounding boxes -
-- If we have the extremities, and the start/end points, a simple for-loop that tests for min/max values for x and y means we have the four - values we need to box in our curve: -
-Computing the bounding box for a Bézier curve:
--
-
- Find all t value(s) for the curve derivative's x- and y-roots. -
- Discard any t value that's lower than 0 or higher than 1, because Bézier curves only use the interval [0,1]. -
- - Determine the lowest and highest value when plugging the values t=0, t=1 and each of the found roots into the original - functions: the lowest value is the lower bound, and the highest value is the upper bound for the bounding box we want to construct. - -
- Applying this approach to our previous root finding, we get the following - axis-aligned bounding boxes (with all curve extremity points - shown on the curve): -
-
-
- 
-
- - We can construct even nicer boxes by aligning them along our curve, rather than along the x- and y-axis, but in order to do so we first - need to look at how aligning works. -
-- - Aligning curves -
-- While there are an incredible number of curves we can define by varying the x- and y-coordinates for the control points, not all curves - are actually distinct. For instance, if we define a curve, and then rotate it 90 degrees, it's still the same curve, and we'll find its - extremities in the same spots, just at different draw coordinates. As such, one way to make sure we're working with a "unique" curve is to - "axis-align" it. -
-- Aligning also simplifies a curve's functions. We can translate (move) the curve so that the first point lies on (0,0), which turns our - n term polynomial functions into n-1 term functions. The order stays the same, but we have less terms. Then, we can - rotate the curves so that the last point always lies on the x-axis, too, making its coordinate (...,0). This further simplifies the - function for the y-component to an 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: -
-
- -
- - - Tight bounding boxes -
-- With our knowledge of bounding boxes, and curve alignment, We can now form the "tight" bounding box for curves. We first align our curve, - recording the translation we performed, "T", and the rotation angle we used, "R". We then determine the aligned curve's normal bounding - box. Once we have that, we can map that bounding box back to our original curve by rotating it by -R, and then translating it by -T. -
-We now have nice tight bounding boxes for our curves:
-
-
- 
-
- - These are, strictly speaking, not necessarily the tightest possible bounding boxes. It is possible to compute the optimal bounding box by - determining which spanning lines we need to effect a minimal box area, but because of the parametric nature of Bézier curves this is - actually a rather costly operation, and the gain in bounding precision is often not worth it. -
-- - Curve inflections -
-- Now that we know how to align a curve, there's one more thing we can calculate: inflection points. Imagine we have a variable size circle - that we can slide up against our curve. We place it against the curve and adjust its radius so that where it touches the curve, the - curvatures of the curve and the circle are the same, and then we start to slide the circle along the curve - for quadratic curves, we can - always do this without the circle behaving oddly: we might have to change the radius of the circle as we slide it along, but it'll always - sit against the same side of the curve. -
-- But what happens with cubic curves? Imagine we have an S curve and we place our circle at the start of the curve, and start sliding it - along. For a while we can simply adjust the radius and things will be fine, but once we get to the midpoint of that S, something odd - happens: the circle "flips" from one side of the curve to the other side, in order for the curvatures to keep matching. This is called an - 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 - circle being on the other side of the curve. So what does C(t) look like? Actually something that seems not too hard: -
- -
- - The function C(t) is the cross product between the first and second derivative functions for the parametric dimensions of our - curve. And, as already shown, derivatives of Bézier curves are just simpler Bézier curves, with very easy to compute new coefficients, so - this should be pretty easy. -
-- However as we've seen in the section on aligning, aligning lets us simplify things a lot, by completely removing the - contributions of the first coordinate from most mathematical evaluations, and removing the last y coordinate as well by virtue of - the last point lying on the x-axis. So, while we can evaluate C(t) = 0 for our curve, it'll be much easier to first axis-align - the curve and then evaluating the curvature function. -
-Let's derive the full formula anyway
-- Of course, before we do our aligned check, let's see what happens if we compute the curvature function without axis-aligning. We start - with the first and second derivatives, given our basis functions: -
- -
- And of course the same functions for y:
- -
- - Asking a computer to now compose the C(t) function for us (and to expand it to a readable form of simple terms) gives us this - rather overly complicated set of arithmetic expressions: -
- -
- - That is... unwieldy. So, we note that there are a lot of terms that involve multiplications involving x1, y1, and y4, which would all - disappear if we axis-align our curve, which is why aligning is a great idea. -
-
- Aligning our curve so that three of the eight coefficients become zero, and observing that scale does not affect finding
- t values, we end up with the following simple term function for C(t):
-
- - That's a lot easier to work with: we see a fair number of terms that we can compute and then cache, giving us the following - simplification: -
- -
- This is a plain quadratic curve, and we know how to solve C(t) = 0; we use the quadratic formula:
- -
- - We can easily compute this value if the discriminator isn't a negative number (because we only want real roots, not complex - roots), and if x is not zero, because divisions by zero are rather useless. -
-- Taking that into account, we compute t, we disregard any t value that isn't in the Bézier interval [0,1], and we now - know at which t value(s) our curve will inflect. -
-
-
- - - The canonical form (for cubic curves) -
-- While quadratic curves are relatively simple curves to analyze, the same cannot be said of the cubic curve. As a curvature is controlled - by more than one control point, it exhibits all kinds of features like loops, cusps, odd colinear features, and as many as two inflection - points because the curvature can change direction up to three times. Now, knowing what kind of curve we're dealing with means that some - algorithms can be run more efficiently than if we have to implement them as generic solvers, so is there a way to determine the curve type - without lots of work? -
-- As it so happens, the answer is yes, and the solution we're going to look at was presented by Maureen C. Stone from Xerox PARC and Tony D. - deRose from the University of Washington in their joint paper - "A Geometric Characterization of Parametric Cubic curves". It was published in 1989, and defines curves as having a "canonical" form (i.e. a form that all curves can be reduced to) from which we - can immediately tell what features a curve will have. So how does it work? -
-- The first observation that makes things work is that if we have a cubic curve with four points, we can apply a linear transformation to - these points such that three of the points end up on (0,0), (0,1) and (1,1), with the last point then being "somewhere". After applying - that transformation, the location of that last point can then tell us what kind of curve we're dealing with. Specifically, we see the - following breakdown: -
-
- 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 - have, if the fourth point is in/on that region or boundary. Specifically, if the fourth point is... -
--
-
-
-
- ...anywhere inside the red zone, but not on its boundaries, the curve will either be self-intersecting (yielding a loop). We won't - know where it self-intersects (in terms of t values), but we are guaranteed that it does. -
-
- -
-
- ...on the left (red) edge of the red zone, the curve will have a cusp. We again don't know where, but we know there is one. - This edge is described by the function: -
- -
-
- -
-
- ...on the almost circular, lower right (pink) edge, the curve's end point touches the curve, forming a loop. This edge is described by - the function: -
- -
-
- -
-
...on the top (blue) edge, the curve's start point touches the curve, forming a loop. This edge is described by the function:
- -
-
- -
-
...inside the lower (green) zone, past
-y=1, the curve will have a single inflection (switching concave/convex once).
- -
-
- ...between the left and lower boundaries (below the cusp line but above the single-inflection line), the curve will have two - inflections (switching from concave to convex and then back again, or from convex to concave and then back again). -
-
- ...anywhere on the right of self-intersection zone, the curve will have no inflections. It'll just be a simple arch.
-
Of course, this map is fairly small, but the regions extend to infinity, with well defined boundaries.
-Wait, where do those lines come from?
-- Without repeating the paper mentioned at the top of this section, the loop-boundaries come from rewriting the curve into canonical form, - and then solving the formulae for which constraints must hold for which possible curve properties. In the paper these functions yield - formulae for where you will find cusp points, or loops where we know t=0 or t=1, but those functions are derived for the full cubic - expression, meaning they apply to t=-∞ to t=∞... For Bézier curves we only care about the "clipped interval" t=0 to t=1, so some of the - properties that apply when you look at the curve over an infinite interval simply don't apply to the Bézier curve interval. -
-- The right bound for the loop region, indicating where the curve switches from "having inflections" to "having a loop", for the general - cubic curve, is actually mirrored over x=1, but for Bézier curves this right half doesn't apply, so we don't need to pay attention to - it. Similarly, the boundaries for t=0 and t=1 loops are also nice clean curves but get "cut off" when we only look at what the general - curve does over the interval t=0 to t=1. -
-- For the full details, head over to the paper and read through sections 3 and 4. If you still remember your high school pre-calculus, you - can probably follow along with this paper, although you might have to read it a few times before all the bits "click". -
-- So now the question becomes: how do we manipulate our curve so that it fits this canonical form, with three fixed points, and one "free" - point? Enter linear algebra. Don't worry, I'll be doing all the math for you, as well as show you what the effect is on our curves, but - basically we're going to be using linear algebra, rather than calculus, because "it's way easier". Sometimes a calculus approach is very - hard to work with, when the equivalent geometrical solution is super obvious. -
-- The approach is going to start with a curve that doesn't have all-colinear points (so we need to make sure the points don't all fall on a - straight line), and then applying three graphics operations that you will probably have heard of: translation (moving all points by some - fixed x- and y-distance), scaling (multiplying all points by some x and y scale factor), and shearing (an operation that turns rectangles - into parallelograms). -
-- Step 1: we translate any curve by -p1.x and -p1.y, so that the curve starts at (0,0). We're going to make use of an interesting trick - here, by pretending our 2D coordinates are 3D, with the z coordinate simply always being 1. This is an old trick in graphics to - overcome the limitations of 2D transformations: without it, we can only turn (x,y) coordinates into new coordinates of the form (ax + by, - 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 - x=0 line, so what we want is a transformation matrix that, when we run it, subtracts x from whatever x we - 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 - do some scaling to make sure it ends up at (0,1). Additionally, we want - 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
- right place, and y-scaling would move it out of (0,1) again, so our only option is to y-shear point three, just like how we x-sheared
- point 2 earlier. In this case, we do the same trick, but with y/x rather than x/y because we're not x-shearing
- 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 - pull this apart a little and end up with an easy-to-calculate bit of code! -
-- 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! -
-How do you track all that?
-- Doing maths can be a pain, so whenever possible, I like to make computers do the work for me. Especially for things like this, I simply - use Mathematica. Tracking all this math by hand is insane, and we invented computers, - literally, to do this for us. I have no reason to use pen and paper when I can write out what I want to do in a program, and have the - program do the math for me. And real math, too, with symbols, not with numbers. In fact, - here's the Mathematica notebook if you want to see how - this works for yourself. -
-- Now, I know, you're thinking "but Mathematica is super expensive!" and that's true, it's - $344 for home use, up from $295 when I original wrote this, but it's - also free when you buy a $35 raspberry pi. Obviously, I bought a - raspberry pi, and I encourage you to do the same. With that, as long as you know what you want to do, Mathematica can just do - it for you. And we don't have to be geniuses to work out what the maths looks like. That's what we have computers for. -
-- So, let's write up a sketch that'll show us the canonical form for any curve drawn in blue, overlaid on our canonical map, so that we can - immediately tell which features our curve must have, based on where the fourth coordinate is located on the map: -
-
- - - Finding Y, given X -
-- One common task that pops up in things like CSS work, or parametric equalizers, or image leveling, or any other number of applications - where Bézier curves are used as control curves in a way that there is really only ever one "y" value associated with one "x" value, you - might want to cut out the middle man, as it were, and compute "y" directly based on "x". After all, the function looks simple enough, - finding the "y" value should be simple too, right? Unfortunately, not really. However, it is possible and as long as you have - some code in place to help, it's not a lot of a work either. -
-
- We'll be tackling this problem in two stages: the first, which is the hard part, is figuring out which "t" value belongs to any given "x"
- value. For instance, have a look at the following graphic. On the left we have a Bézier curve that looks for all intents and purposes like
- it fits our criteria: every "x" has one and only one associated "y" value. On the right we see the function for just the "x" values:
- that's a cubic curve, but not a really crazy cubic curve. If you move the graphic's slider, you will see a red line drawn that corresponds
- to the x coordinate: this is a vertical line in the left graphic, and a horizontal line on the right.
-
-
- - 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 - we've already seen how to calculate roots, so let's just run cubic root finding - and not even the complicated - 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: -
- -| 1 | -- - | -
| 2 | -|
| 3 | -|
| 4 | -|
| 5 | -|
| 6 | -|
| 7 | -|
| 8 | -|
| 9 | -|
| 10 | -
- So the procedure is fairly straight forward: pick an x, find the associated 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
- x. Move the slider for the following graphic to see this in action:
-
-
- - - Arc length -
-- How long is a Bézier curve? As it turns out, that's not actually an easy question, because the answer requires maths that —much like root - finding— cannot generally be solved the traditional way. If we have a parametric curve with fx(t) and - 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, - it's far from simple... cutting straight to after the chase is over: for quadratic curves, this formula generates an - unwieldy computation, and we're simply not going to implement things that way. For cubic Bézier curves, things get even more fun, because there is no "closed - form" solution, meaning that due to the way calculus works, there is no generic formula that allows you to calculate the arc length. Let - me just repeat this, because it's fairly crucial: - for cubic and higher Bézier curves, there is no way to solve this function if you want to use it "for all possible coordinates". -
-Seriously: It cannot be done.
-- So we turn to numerical approaches again. The method we'll look at here is the - Gauss quadrature. This approximation - is a really neat trick, because for any nth degree polynomial it finds approximated values for an integral really - efficiently. Explaining this procedure in length is way beyond the scope of this page, so if you're interested in finding out why it - 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: -
- -
- - 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 - more strips we use (and of course the more we use, the thinner they get) the closer we get to the true area under the curve, and thus the - better the approximation: -
-
-
- 
-
- 
-
- - Now, infinitely many terms to sum and infinitely thin rectangles are not something that computers can work with, so instead we're going to - approximate the infinite summation by using a sum of a finite number of "just thin" rectangular strips. As long as we use a high enough - number of thin enough rectangular strips, this will give us an approximation that is pretty close to what the real value is. -
-- So, the trick is to come up with useful rectangular strips. A naive way is to simply create n strips, all with the same width, - but there is a far better way using special values for C and f(t) depending on the value of n, which indicates - how many strips we'll use, and it's called the Legendre-Gauss quadrature. -
-- This approach uses strips that are not spaced evenly, but instead spaces them in a special way based on describing the function - as a polynomial (the more strips, the more accurate the polynomial), and then computing the exact integral for that polynomial. We're - essentially performing arc length computation on a flattened curve, but flattening it based on the intervals dictated by the - Legendre-Gauss solution. -
-- Note that one requirement for the approach we'll use is that the integral must run from -1 to 1. That's no good, because we're dealing - with Bézier curves, and the length of a section of curve applies to values which run from 0 to "some value smaller than or equal to 1" - (let's call that value z). Thankfully, we can quite easily transform any integral interval to any other integral interval, by - shifting and scaling the inputs. Doing so, we get the following: -
- -
- - That may look a bit more complicated, but the fraction involving z is a fixed number, so the summation, and the evaluation of - the f(t) values are still pretty simple. -
-- So, what do we need to perform this calculation? For one, we'll need an explicit formula for f(t), because that derivative - notation is handy on paper, but not when we have to implement it. We'll also need to know what these Ci and - ti values should be. Luckily, that's less work because there are actually many tables available that give these - 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: -
- -
- - We can program that pretty easily, provided we have that f(t) available, which we do, as we know the full description for the - Bézier curve functions Bx(t) and By(t). -
-- If we use the Legendre-Gauss values for our C values (thickness for each strip) and t values (location of each strip), - we can determine the approximate length of a Bézier curve by computing the Legendre-Gauss sum. The following graphic shows a cubic curve, - with its computed lengths; Go ahead and change the curve, to see how its length changes. One thing worth trying is to see if you can make - a straight line, and see if the length matches what you'd expect. What if you form a line with the control points on the outside, and the - start/end points on the inside? -
-
-
- - - Approximated arc length -
-- Sometimes, we don't actually need the precision of a true arc length, and we can get away with simply computing the approximate arc length - instead. The by far fastest way to do this is to flatten the curve and then simply calculate the linear distance from point to point. This - will come with an error, but this can be made arbitrarily small by increasing the segment count. -
-- If we combine the work done in the previous sections on curve flattening and arc length computation, we can implement these with minimal - effort: -
-
-
- 
-
- - You may notice that even though the error in length is actually pretty significant in absolute terms, even at a low number of segments we - get a length that agrees with the true length when it comes to just the integer part of the arc length. Quite often, approximations can - drastically speed things up! -
-- - Curvature of a curve -
-- If we have two curves, and we want to line them in up in a way that "looks right", what would we use as metric to let a computer decide - what "looks right" means? -
-- For instance, we can start by ensuring that the two curves share an end coordinate, so that there is no "gap" between the end of one and - the start of the next curve, but that won't guarantee that things look right: both curves can be going in wildly different directions, and - the resulting joined geometry will have a corner in it, rather than a smooth transition from one curve to the next. -
-- What we want is to ensure that the curvature at the transition from one curve to the - next "looks good". So, we start with a shared coordinate, and then also require that derivatives for both curves match at that coordinate. - That way, we're assured that their tangents line up, which must mean the curve transition is perfectly smooth. We can even make the - second, third, etc. derivatives match up for better and better transitions. -
-Problem solved!
-- However, there's a problem with this approach: if we think about this a little more, we realise that "what a curve looks like" and its - derivative values are pretty much entirely unrelated. After all, the section on reordering curves showed us that - the same looking curve can have an infinite number of curve expressions of arbitrarily high Bézier degree, and each of those will have - wildly different derivative values. -
-
- So what we really want is some kind of expression that's not based on any particular expression of t, but is based on
- something that is invariant to the kind of function(s) we use to draw our curve. And the prime candidate for this is our curve
- expression, reparameterised for distance: no matter what order of Bézier curve we use, if we were able to rewrite it as a function of
- distance-along-the-curve, all those different degree Bézier functions would end up being the same function for "coordinate at
- some distance D along the curve".
-
We've seen this before... that's the arc length function.
-- So you might think that in order to find the curvature of a curve, we now need to solve the arc length function itself, and that this - would be quite a problem because we just saw that there is no way to actually do that. Thankfully, we don't. We only need to know the - form of the arc length function, which we saw above and is fairly simple, rather than needing to solve the arc length - function. If we start with the arc length expression and the - run through the steps necessary to determine its derivative (with an - alternative, shorter demonstration of how to do this found - over on Stackexchange), then the - integral that was giving us so much problems in solving the arc length function disappears entirely (because of the - fundamental theorem of calculus), and what we're left with us - some surprisingly simple maths that relates curvature (denoted as κ, "kappa") to—and this is the truly surprising bit—a specific - combination of derivatives of our original function. -
-Let me highlight what just happened, because it's pretty special:
--
-
- we wanted to make curves line up, and initially thought to match the curves' derivatives, but -
- that turned out to be a really bad choice, so instead -
- we picked a function that is basically impossible to work with, and then worked with that, which -
- gives us a simple formula that is and expression using the curves' derivatives. -
That's crazy!
-- But that's also one of the things that makes maths so powerful: even if your initial ideas are off the mark, you might be much closer than - you thought you were, and the journey from "thinking we're completely wrong" to "actually being remarkably close to being right" is where - we can find a lot of insight. -
-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 little 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 - looks oddly similar to the standard Euclidean distance function. And nothing in these functions is hard to calculate either: for Bézier - curves, simply knowing our curve coordinates means we know what the first and second derivatives are, and so - evaluating this function for any t value is just a matter of basic arithematics. -
-In fact, let's just implement it right now:
- -| 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 anyway) -
-- With all of that covered, let's line up some curves! The following graphic gives you two curves that look identical, but use quadratic and - cubic functions, respectively. As you can see, despite their derivatives being necessarily different, their curvature (thanks to being - derived based on maths that "ignores" specific function derivative, and instead gives a formula that smooths out any differences) is - exactly the same. And because of that, we can put them together such that the point where they overlap has the same curvature for both - curves, giving us the smoothest transition. -
-
- - 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: -
-
-
- - - Tracing a curve at fixed distance intervals -
-- Say you want to draw a curve with a dashed line, rather than a solid line, or you want to move something along the curve at fixed distance - intervals over time, like a train along a track, and you want to use Bézier curves. -
-Now you have a problem.
-
- The reason you have a problem is that Bézier curves are parametric functions with non-linear behaviour, whereas moving a train along a
- track is about as close to a practical example of linear behaviour as you can get. The problem we're faced with is that we can't just pick
- t values at some fixed interval and expect the Bézier functions to generate points that are spaced a fixed distance apart. In
- fact, let's look at the relation between "distance long a curve" and "t value", by plotting them against one another.
-
- The following graphic shows a particularly illustrative curve, and it's distance-for-t plot. For linear traversal, this line needs to be - straight, running from (0,0) to (length,1). That is, it's safe to say, not what we'll see: we'll see something very wobbly, instead. To - make matters even worse, the distance-for-t function is also of a much higher order than our curve is: while the curve we're using for - this exercise is a cubic curve, which can switch concave/convex form twice at best, the distance function is our old friend the arc length - function, which can have more inflection points. -
-
-
- 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
- then we find "the closest t value that matches some required distance" using those values instead. If we have a low number of
- points sampled, we can then even refine which t value "should" work for our desired distance by interpolating between two
- points, but if we have a high enough number of samples, we don't even need to bother.
-
- So let's do exactly that: the following graph is similar to the previous one, showing how we would have to "chop up" our distance-for-t
- curve in order to get regularly spaced points on the curve. It also shows what using those t values on the real curve looks
- like, by coloring each section of curve between two distance markers differently:
-
-
- 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 "Moving Along a Curve with Specified Speed" by David Eberly of Geometric Tools, LLC, but basically because we have no explicit length function (or rather, one we don't have to - constantly compute for different intervals), you may simply be better off with a traditional lookup table (LUT). -
-- - Intersections -
-- Let's look at some more things we will want to do with Bézier curves. Almost immediately after figuring out how to get bounding boxes to - work, people tend to run into the problem that even though the minimal bounding box (based on rotation) is tight, it's not sufficient to - perform true collision detection. It's a good first step to make sure there might be a collision (if there is no bounding box - overlap, there can't be one), but in order to do real collision detection we need to know whether or not there's an intersection on the - actual curve. -
-- We'll do this in steps, because it's a bit of a journey to get to curve/curve intersection checking. First, let's start simple, by - implementing a line-line intersection checker. While we can solve this the traditional calculus way (determine the functions for both - lines, then compute the intersection by equating them and solving for two unknowns), linear algebra actually offers a nicer solution. -
-Line-line intersections
-- If we have two line segments with two coordinates each, segments A-B and C-D, we can find the intersection of the lines these segments are - an intervals on by linear algebra, using the procedure outlined in this - top coder - article. Of course, we need to make sure that the intersection isn't just on the lines our line segments lie on, but actually on our line - segments themselves. So after we find the intersection, we need to verify that it lies without the bounds of our original line segments. -
-- The following graphic implements this intersection detection, showing a red point for an intersection on the lines our segments lie on - (thus being a virtual intersection point), and a green point for an intersection that lies on both segments (being a real intersection - point). -
-
-
- Implementing line-line intersections
-- Let's have a look at how to implement a line-line intersection checking function. The basics are covered in the article mentioned above, - 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: -
- -| 1 | -- - | -
| 2 | -|
| 3 | -|
| 4 | -|
| 5 | -|
| 6 | -|
| 7 | -|
| 8 | -|
| 9 | -|
| 10 | -|
| 11 | -|
| 12 | -|
| 13 | -|
| 14 | -|
| 15 | -|
| 16 | -|
| 17 | -
What about curve-line intersections?
-- Curve/line intersection is more work, but we've already seen the techniques we need to use in order to perform it: first we - translate/rotate both the line and curve together, in such a way that the line coincides with the x-axis. This will position the curve in - a way that makes it cross the line at points where its y-function is zero. By doing this, the problem of finding intersections between a - curve and a line has now become the problem of performing root finding on our translated/rotated curve, as we already covered in the - section on finding extremities. -
-
-
- 
-
- - Curve/curve intersection, however, is more complicated. Since we have no straight line to align to, we can't simply align one of the - curves and be left with a simple procedure. Instead, we'll need to apply two techniques we've met before: de Casteljau's algorithm, and - curve splitting. -
-- - Curve/curve intersection -
-- Using de Casteljau's algorithm to split the curve we can now implement curve/curve intersection finding using a "divide and conquer" - technique: -
--
-
- - Take two curves C1 and C2, and treat them as a pair. - -
- If their bounding boxes overlap, split up each curve into two sub-curves -
- - With C1.1, C1.2, C2.1 and C2.2, form four new pairs (C1.1,C2.1), (C1.1, C2.2), (C1.2,C2.1), and (C1.2,C2.2). - -
-
- For each pair, check whether their bounding boxes overlap.
-
-
-
- If their bounding boxes do not overlap, discard the pair, as there is no intersection between this pair of curves. -
- If there is overlap, rerun all steps for this pair. -
- -
- Once the sub-curves we form are so small that they effectively occupy sub-pixel areas, we consider an intersection found, noting that we
- might have a cluster of multiple intersections at the sub-pixel level, out of which we pick one to act as "found"
tvalue - (we can either throw all but one away, we can average the cluster'stvalues, or you can do something even more creative). -
-
- This algorithm will start with a single pair, "balloon" until it runs in parallel for a large number of potential sub-pairs, and then - taper back down as it homes in on intersection coordinates, ending up with as many pairs as there are intersections. -
-- The following graphic applies this algorithm to a pair of cubic curves, one step at a time, so you can see the algorithm in action. Click - the button to run a single step in the algorithm, after setting up your curves in some creative arrangement. You can also change the value - that is used in step 5 to determine whether the curves are small enough. Manipulating the curves or changing the threshold will reset the - algorithm, so you can try this with lots of different curves. -
-(can you find the configuration that yields the maximum number of intersections between two cubic curves? Nine intersections!)
-
-
-
- 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
- t=0.5 first, and then running the exact same algorithm as above, with all non-overlapping curve pairs getting removed at each
- iteration, and each successive step homing in on the curve's self-intersection points.
-
- - The projection identity -
-- De Casteljau's algorithm is the pivotal algorithm when it comes to Bézier curves. You can use it not just to split curves, but also to - draw them efficiently (especially for high-order Bézier curves), as well as to come up with curves based on three points and a tangent. - Particularly this last thing is really useful because it lets us "mold" a curve, by picking it up at some point, and dragging that point - around to change the curve's shape. -
-How does that work? Succinctly: we run de Casteljau's algorithm in reverse!
-- In order to run de Casteljau's algorithm in reverse, we need a few basic things: a start and end point, a point on the curve that want to - be moving around, which has an associated t value, and a point we've not explicitly talked about before, and as far as I know has - no explicit name, but lives one iteration higher in the de Casteljau process then our on-curve point does. I like to call it "A" for - reasons that will become obvious. -
-
- So let's use graphics instead of text to see where this "A" is, because text only gets us so far: move the sliders for the following
- graphics to see what, given specific t value, our A coordinate is. As well as some other coordinates, which
- taken together let us derive a value that the graphics call "ratio": if you move the curve's points around, A, B, and C will move, what
- happens to that value?
-
-
- 
-
- So these graphics show us several things:
--
-
- a point at the tip of the curve construction's "hat": let's call that
A, as well as
- - our on-curve point give our chosen
tvalue: let's call thatB, and finally,
- -
- a point that we get by projecting A, through B, onto the line between the curve's start and end points: let's call that
C. -
- -
- for both quadratic and cubic curves, two points
e1ande2, which represent the single-to-last step in de - Casteljau's algorithm: in the last step, we findBat(1-t) * e1 + t * e2. -
- -
- for cubic curves, also the points
v1andv2, which together withArepresent the first step in de - Casteljau's algorithm: in the next step, we finde1ande2. -
-
- These three values A, B, and C allow us to derive an important identity formula for quadratic and cubic Bézier curves: for any point on
- the curve with some t value, the ratio of distances from A to B and B to C is fixed: if some t value sets up a C
- that is 20% away from the start and 80% away from the end, then it doesn't matter where the start, end, or control points are;
- for that t value, C will always lie at 20% from the start and 80% from the end point. Go ahead, pick an
- on-curve point in either graphic and then move all the other points around: if you only move the control points, start and end won't move,
- and so neither will C, and if you move either start or end point, C will move but its relative position will not change.
-
- So, how can we compute 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:
-
-
- If we can figure out what the function 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.
- Running through the maths
- (with thanks to Boris Zbarsky) shows us the following two formulae:
-
- And
- -
-
- So, if we know the start and end coordinates, and we know the t value, we know C, without having to calculate the
- A or even B coordinates. In fact, we can do the same for the ratio function: as another function of
- t, we technically don't need to know what A or B or C are, we can express it was a
- pure function of t, too.
-
We start by observing that, given A, B, and C, the following always holds:
- 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
- t value, we can construct 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:
-
-
- With 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
- e2. For cubic curves, there is no single pair of points that can act as e1 and e2: as long as the
- 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:
-
- And then reverse engineer the curve's control control points:
- -
-
- So: if we have a curve's start and end point, then for any 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
- to reconstruct a curve's "de Casteljau skeleton". Which means that we can now do several things: we can "fit" curves using only three
- points, which means we can also "mold" curves by moving an on-curve point but leaving its start and end point, and then reconstructing the
- curve based on where we moved the on-curve point to. These are very useful things, and we'll look at both in the next few sections.
-
- - Creating a curve from three points -
-- Given the preceding section, you might be wondering if we can use that knowledge to just "create" curves by placing some points and having - the computer do the rest, to which the answer is: that's exactly what we can now do! -
-
- For quadratic curves, things are pretty easy. Technically, we'll need a t value in order to compute the ratio function used
- in computing the ABC coordinates, but we can just as easily approximate one by treating the distance between the start and
- B point, and B and end point as a ratio, using
-
-
- With this code in place, creating a quadratic curve from three points is literally just computing the ABC values, and using
- A as our curve's control point:
-
-
- - 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 - you ever learned it!) that a line between two points on a circle is called a - chord, and that one property of chords is that the line from the center - of any chord, perpendicular to that chord, passes through the center of the circle. -
-- That means that if we have have three points on a circle, we have three (different) chords, and consequently, three (different) lines that - go from those chords through the center of the circle: if we find two of those lines, then their intersection will be our circle's center, - and the circle's radius will—by definition!—be the distance from the center to any of our three points: -
-
-
-
- With that covered, we now also know the tangent line to our point 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
- e1 and e2 on that tangent, so that we can construct a new cubic curve hull. We use the approach as we did for
- quadratic curves to automatically determine a reasonable t value, and then our e1 and
- e2 coordinates must obey the standard de Casteljau rule for linear interpolation:
-
-
- Where d is the total length of the line segment from e1 to e2. So how long do we make that? There
- are again all kinds of approaches we can take, and a simple-but-effective one is to set the length of that segment to "one third the
- length of the baseline". This forces e1 and e2 to always be the "linear curve" distance apart, which means if we
- place our three points on a line, it will actually look like a line. Nice! The last thing we'll need to do is make sure to flip
- the sign of d depending on which side of the baseline our B is located, so we don't up creating a funky curve
- with a loop in it. To do this, we can use the atan2 function:
-
-
- This angle φ will be between 0 and π if B is "above" the baseline (rotating all three points so that the start is on the left
- and the end is the right), so we can use a relatively straight forward check to make sure we're using the correct sign for our value
- d:
-
- 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 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;
- more on that later), we're only trying to come up with some reasonable e1 and
- e2 points so we can construct a new cubic curve... so now that we have those: let's see what kind of cubic curve that gives
- us:
-
-
- That looks perfectly serviceable!
-- 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" - curves, where we can reshape a curve by dragging a point on the curve around while leaving the start and end fixed, effectively molding - the shape as if it were clay or the like. We'll see the last tool we need to do that in the next section, and then we'll look at - implementing curve molding in the section after that, so read on! -
-- - Projecting a point onto a Bézier curve -
-- Before we can move on to actual curve molding, it'll be good if know how to actually be able to find "some point on the curve" that we're - trying to click on. After all, if all we have is our Bézier coordinates, that is not in itself enough to figure out which point on the - curve our cursor will be closest to. So, how do we project points onto a curve? -
-
- If the Bézier curve is of low enough order, we might be able to
- work out the maths for how to do this, and get a perfect t value back, but in general this is an incredibly hard problem and the easiest solution is, really, a
- numerical approach again. We'll be finding our ideal t value using a
- binary search. First, we do a coarse distance-check based on
- t values associated with the curve's "to draw" coordinates (using a lookup table, or LUT). This is pretty fast:
-
| 1 | -- - | -
| 2 | -|
| 3 | -|
| 4 | -|
| 5 | -|
| 6 | -|
| 7 | -
- After this runs, we know that 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
- better guess than both LUT[i-1] and LUT[i+1], but there might be an even better projection somewhere else between those two
- values, so that's what we're going to be testing for, using a variation of the binary search.
-
-
-
-
- we start with our point
p, and thetvaluest1=LUT[i-1].tandt2=LUT[i+1].t, which - span an intervalv = t2-t1. -
- - - we test this interval in five spots: the start, middle, and end (which we already have), and the two points in between the middle and - start/end points - -
-
- we then check which of these five points is the closest to our original point
p, and then repeat step 1 with the points - before and after the closest point we just found. -
-
- This makes the interval we check smaller and smaller at each iteration, and we can keep running the three steps until the interval becomes - so small as to lead to distances that are, for all intents and purposes, the same for all points. -
-- So, let's see that in action: in this case, I'm going to arbitrarily say that if we're going to run the loop until the interval is smaller - than 0.001, and show you what that means for projecting your mouse cursor or finger tip onto a rather complex Bézier curve (which, of - course, you can reshape as you like). Also shown are the original three points that our coarse check finds. -
-
- - - Molding a curve -
-- Armed with knowledge of the "ABC" relation, point-on-curve projection, and guestimating reasonable looking helper values for cubic curve - construction, we can finally cover curve molding: updating a curve's shape interactively, by dragging points on the curve around. -
-
- For quadratic curve, this is a really simple trick: we project our cursor onto the curve, which gives us a t value and
- initial B coordinate. We don't even need the latter: with our t value and "wherever the cursor is" as target
- B, we can compute the associated C:
-
- And then the associated A:
- 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 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
- creation, we were free to pick an appropriate t value ourselves, which allowed us to find appropriate e1 and
- e2 coordinates. That's great, but when we're curve molding we don't have that luxury: whatever point we decide to start
- moving around already has its own t value, and its own e1 and e2 values, and those may not make
- sense for the rest of the curve.
-
- For example, let's see what happens if we just "go with what we get" when we pick a point and start moving it around, preserving its
- t value and e1/e2 coordinates:
-
- - 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. -
-
- One way to combat this might be to combine the above approach with the approach from the
- creating curves section: generate both the "unchanged t/e1/e2" curve, as
- well as the "idealized" curve through the start/cursor/end points, with idealized t value, and then interpolating between
- those two curves:
-
-
-
- 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 t/e1/e2" closer to the original point, and bias
- towards "idealized" form the further away we move our point, with anything that's further than our falloff distance simply
- being the idealized curve. We don't even try to interpolate at that point.
-
- A more advanced way to try to smooth things out is to implement continuous molding, where we constantly update the curve as we
- move around, and constantly change what our B point is, based on constantly projecting the cursor on the curve
- as we're updating it - this is, you won't be surprised to learn, tricky, and beyond the scope of this section: interpolation
- (with a reasonable distance) will do for now!
-
- - Curve fitting -
-
- Given the previous section, one question you might have is "what if I don't want to guess t values?". After all, plenty of
- graphics packages do automated curve fitting, so how can we implement that in a way that just finds us reasonable t values
- all on its own?
-
- And really this is just a variation on the question "how do I get the curve through these X points?", so let's look at that. Specifically,
- let's look at the answer: "curve fitting". This is in fact a rather rich field in geometry, applying to anything from data modelling to
- path abstraction to "drawing", so there's a fair number of ways to do curve fitting, but we'll look at one of the most common approaches:
- something called a least squares
- polynomial regression. In this approach, we look at the number of points
- we have in our data set, roughly determine what would be an appropriate order for a curve that would fit these points, and then tackle the
- question "given that we want an nth order curve, what are the coordinates we can find such that our curve is "off" by the
- least amount?".
-
- Now, there are many ways to determine how "off" points are from the curve, which is where that "least squares" term comes in. The most - common tool in the toolbox is to minimise the squared distance between each point we have, and the corresponding point on the - curve we end up "inventing". A curve with a snug fit will have zero distance between those two, and a bad fit will have non-zero distances - between every such pair. It's a workable metric. You might wonder why we'd need to square, rather than just ensure that distance is a - positive value (so that the total error is easy to compute by just summing distances) and the answer really is "because it tends to be a - little better". There's lots of literature on the web if you want to deep-dive the specific merits of least squared error metrics versus - least absolute error metrics, but those are well beyond the scope of this material. -
-- So let's look at what we end up with in terms of curve fitting if we start with the idea of performing least squares Bézier fitting. We're - going to follow a procedure similar to the one described by Jim Herold over on his - "Least Squares Bézier Fit" - article, and end with some nice interactive graphics for doing some curve fitting. -
-- Before we begin, we're going to use the curve in matrix form. In the section on matrices, I mentioned that some - things are easier if we use the matrix representation of a Bézier curve rather than its calculus form, and this is one of those things. -
-- As such, the first step in the process is expressing our Bézier curve as powers/coefficients/coordinate matrix T x M x C, - by expanding the Bézier functions. -
-Revisiting the matrix representation
-- Rewriting Bézier functions to matrix form is fairly easy, if you first expand the function, and then arrange them into a multiple line - form, where each line corresponds to a power of t, and each column is for a specific coefficient. First, we expand the function: -
- -
- And then we (trivially) rearrange the terms across multiple lines:
- -
- - This rearrangement has "factors of t" at each row (the first row is t⁰, i.e. "1", the second row is t¹, i.e. "t", the third row is t²) - and "coefficient" at each column (the first column is all terms involving "a", the second all terms involving "b", the third all terms - involving "c"). -
-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):
- -
- - And so and on so on. Now, let's see how to use these T, M, and C, to do some curve - fitting. -
-
- Let's get started: we're going to assume we picked the right order curve: for n points we're fitting an n-1th order curve, so we "start" with a vector P that represents the coordinates we already know, and for which
- we want to do curve fitting:
-
-
- Next, we need to figure out appropriate t values for each point in the curve, because we need something that lets us tie "the
- actual coordinate" to "some point on the curve". There's a fair number of different ways to do this (and a large part of optimizing "the
- perfect fit" is about picking appropriate t values), but in this case let's look at two "obvious" choices:
-
-
-
- equally spaced
tvalues, and
- tvalues that align with distance along the polygon.
-
- The first one is really simple: if we have n points, then we'll just assign each point i a t value
- of (i-1)/(n-1). So if we have four points, the first point will have t=(1-1)/(4-1)=0/3, the second point will
- have t=(2-1)/(4-1)=1/3, the third point will have t=2/3, and the last point will be t=1. We're just
- straight up spacing the t values to match the number of points we have.
-
- The second one is a little more interesting: since we're doing polynomial regression, we might as well exploit the fact that our base
- coordinates just constitute a collection of line segments. At the first point, we're fixing t=0, and the last point, we want t=1, and
- anywhere in between we're simply going to say that t is equal to the distance along the polygon, scaled to the [0,1] domain.
-
To get these values, we first compute the general "distance along the polygon" matrix:
- -
-
- Where length() is literally just that: the length of the line segment between the point we're looking at, and the previous
- point. This isn't quite enough, of course: we still need to make sure that all the values between i=1 and
- i=n fall in the [0,1] interval, so we need to scale all values down by whatever the total length of the polygon is:
-
- - And now we can move on to the actual "curve fitting" part: what we want is a function that lets us compute "ideal" control point values - such that if we build a Bézier curve with them, that curve passes through all our original points. Or, failing that, have an overall error - distance that is as close to zero as we can get it. So, let's write out what the error distance looks like. -
-
- As mentioned before, this function is really just "the distance between the actual coordinate, and the coordinate that the curve evaluates
- to for the associated t value", which we'll square to get rid of any pesky negative signs:
-
- - Since this function only deals with individual coordinates, we'll need to sum over all coordinates in order to get the full error - function. So, we literally just do that; the total error function is simply the sum of all these individual errors: -
- -
- - And here's the trick that justifies using matrices: while we can work with individual values using calculus, with matrices we can compute - as many values as we make our matrices big, all at the "same time", We can replace the individual terms pi with the full - P coordinate matrix, and we can replace Bézier(si) with the matrix representation - T x M x C we talked about before, which gives us: -
- -
- In which we can replace the rather cumbersome "squaring" operation with a more conventional matrix equivalent:
- -
-
- Here, the letter T is used instead of the number 2, to represent the
- matrix transpose; each row in the original matrix becomes a column in the transposed
- matrix instead (row one becomes column one, row two becomes column two, and so on).
-
- This leaves one problem: T isn't actually the matrix we want: we don't want symbolic t values, we want the
- actual numerical values that we computed for S, so we need to form a new matrix, which we'll call 𝕋, that makes use of
- those, and then use that 𝕋 instead of T in our error function:
-
- Which, because of the first and last values in S, means:
- -
- Now we can properly write out the error function as matrix operations:
- -
- - So, we have our error function: we now need to figure out the expression for where that function has minimal value, e.g. where the error - between the true coordinates and the coordinates generated by the curve fitting is smallest. Like in standard calculus, this requires - taking the derivative, and determining where that derivative is zero: -
- -
- Where did this derivative come from?
-- That... is a good question. In fact, when trying to run through this approach, I ran into the same question! And you know what? I - straight up had no idea. I'm decent enough at calculus, I'm decent enough at linear algebra, and I just don't know. -
-- So I did what I always do when I don't understand something: I asked someone to help me understand how things work. In this specific - case, I - posted a question - to Math.stackexchange, and received a answer that goes into way more detail than I had - hoped to receive. -
-- Is that answer useful to you? Probably: no. At least, not unless you like understanding maths on a recreational level. And I do mean - maths in general, not just basic algebra. But it does help in giving us a reference in case you ever wonder "Hang on. Why was that - true?". There are answers. They might just require some time to come to understand. -
-- Now, given the above derivative, we can rearrange the terms (following the rules of matrix algebra) so that we end up with an expression - for C: -
- -
-
- Here, the "to the power negative one" is the notation for the
- matrix inverse. But that's all we have to do: we're done. Starting with
- P and inventing some t values based on the polygon the coordinates in P define, we can
- compute the corresponding Bézier coordinates C that specify a curve that goes through our points. Or, if it can't go
- through them exactly, as near as possible.
-
- So before we try that out, how much code is involved in implementing this? Honestly, that answer depends on how much you're going to be - writing yourself. If you already have a matrix maths library available, then really not that much code at all. On the other hand, if you - are writing this from scratch, you're going to have to write some utility functions for doing your matrix work for you, so it's really - anywhere from 50 lines of code to maybe 200 lines of code. Not a bad price to pay for being able to fit curves to pre-specified - coordinates. -
-- So let's try it out! The following graphic lets you place points, and will start computing exact-fit curves once you've placed at least - three. You can click for more points, and the code will simply try to compute an exact fit using a Bézier curve of the appropriate order. - Four points? Cubic Bézier. Five points? Quartic. And so on. Of course, this does break down at some point: depending on where you place - your points, it might become mighty hard for the fitter to find an exact fit, and things might actually start looking horribly off once - there's enough points for compound - floating point rounding errors to start making a - difference (which is around 10~11 points). -
-
-
-
- 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
- direct control over the time values through sliders for each, because if the time values are our degree of freedom, you should be
- able to freely manipulate them and see what the effect on your curve is.
-
- - Bézier curves and Catmull-Rom curves -
-- Taking an excursion to different splines, the other common design curve is the - Catmull-Rom spline, which unlike Bézier curves - pass through each control point, so they offer a kind of "built-in" curve fitting. -
-- In fact, let's start with just playing with one: the following graphic has a predefined curve that you manipulate the points for, and lets - you add points by clicking/tapping the background, as well as let you control "how fast" the curve passes through its point using the - tension slider. The tenser the curve, the more the curve tends towards straight lines from one point to the next. -
-
-
- - 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 - curve built up of lots of identically-computed pieces, similar to if you just took a whole bunch of Bézier curves, placed them end to end, - and lined up their control points so that things look like a single curve. For a Catmull-Rom curve, each "piece" between two points is - defined by the point's coordinates, and the tangent for those points, the latter of which - can trivially be derived from knowing the - previous and next point: -
- -
- - One downside of this is that—as you may have noticed from the graphic—the first and last point of the overall curve don't actually join up - with the rest of the curve: they don't have a previous/next point respectively, and so there is no way to calculate what their tangent - should be. Which also makes it rather tricky to fit a Catmull-Rom curve to three points like we were able to do for Bézier curves. More on - that in the next section. -
-- In fact, before we move on, let's look at how to actually draw the basic form of these curves (I say basic, because there are a number of - variations that make things - considerable more - complex): -
- -| 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 convert - one to the other and back, and the maths for doing so is surprisingly simple! -
-The main difference between Catmull-Rom curves and Bézier curves is "what the points mean":
--
-
- - A cubic Bézier curve is defined by a start point, a control point that implies the tangent at the start, a control point that implies - the tangent at the end, and an end point, plus a characterizing matrix that we can multiply by that point vector to get on-curve - coordinates. - -
- - A Catmull-Rom curve is defined by a start point, a tangent that for that starting point, an end point, and a tangent for that end point, - plus a characteristic matrix that we can multiple by the point vector to get on-curve coordinates. - -
- Those are very similar, so let's see exactly how similar they are. We've already see the matrix form for Bézier curves, - so how different is the matrix form for Catmull-Rom curves?: -
- -
- - That's pretty dang similar. So the question is: how can we convert that expression with Catmull-Rom matrix and vector into an expression - of the Bézier matrix and vector? The short answer is of course "by using linear algebra", but the longer answer is the rest of this - section, and involves some maths that you may not even care for: if you just want to know the (incredibly simple) conversions between the - two curve forms, feel free to skip to the end of the following explanation, but if you want to how we can get one from the - other... let's get mathing! -
-Deriving the conversion formulae
-- In order to convert between Catmull-Rom curves and Bézier curves, we need to know two things. Firstly, how to express the Catmull-Rom - curve using a "set of four coordinates", rather than a mix of coordinates and tangents, and secondly, how to convert those Catmull-Rom - coordinates to and from Bézier form. -
-- We start with the first part, to figure out how we can go from Catmull-Rom V coordinates to Bézier - P coordinates, by applying "some matrix T". We don't know what that T is yet, but - we'll get to that: -
- -
- - So, this mapping says that in order to map a Catmull-Rom "point + tangent" vector to something based on an "all coordinates" vector, we - need to determine the mapping matrix such that applying T yields P2 as start point, P3 as end point, and two tangents based on - the lines between P1 and P3, and P2 nd P4, respectively. -
-Computing T is really more "arranging the numbers":
- -
- Thus:
- -
- - However, we're not quite done, because Catmull-Rom curves have that "tension" parameter, written as τ (a lowercase"tau"), which - is a scaling factor for the tangent vectors: the bigger the tension, the smaller the tangents, and the smaller the tension, the bigger - the tangents. As such, the tension factor goes in the denominator for the tangents, and before we continue, let's add that tension - factor into both our coordinate vector representation, and mapping matrix T: -
- -
- - With the mapping matrix properly done, let's rewrite the "point + tangent" Catmull-Rom matrix form to a matrix form in terms of four - coordinates, and see what we end up with: -
- -
- 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:
- -
- - A matrix times its inverse is the matrix equivalent of 1, and because "something times 1" is the same as "something", so we can just - outright remove any matrix/inverse pair: -
- -
- 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.
--
-
- Start with the Catmull-Rom function: -
- -
-
- rewrite to pure coordinate form: -
- -
-
- rewrite for "normal" coordinate vector: -
- -
-
- merge the inner matrices: -
- -
-
- rewrite for Bézier matrix form: -
- -
-
- and transform the coordinates so we have a "pure" Bézier expression: -
- And we're done: we finally know how to convert these two curves!
-- If we have a Catmull-Rom curve defined by four coordinates P1 through P4, then we can draw that curve using a Bézier - curve that has the vector: -
- -
- - Similarly, if we have a Bézier curve defined by four coordinates P1 through P4, we can draw that using a standard - tension Catmull-Rom curve with the following coordinate values: -
- -
- Or, if your API allows you to specify Catmull-Rom curves using plain coordinates:
- -
- - - Creating a Catmull-Rom curve from three points -
-- Much shorter than the previous section: we saw that Catmull-Rom curves need at least 4 points to draw anything sensible, so how do we - create a Catmull-Rom curve from three points? -
-Short and sweet: we don't.
-- We run through the maths that lets us create a cubic Bézier curve, and then convert its coordinates to - Catmull-Rom form using the conversion formulae we saw above. -
-- - Forming poly-Bézier curves -
-- Much like lines can be chained together to form polygons, Bézier curves can be chained together to form poly-Béziers, and the only trick - required is to make sure that: -
--
-
- the end point of each section is the starting point of the following section, and -
- the derivatives across that dual point line up. -
Unless you want sharp corners, of course. Then you don't even need 2.
-- We'll cover three forms of poly-Bézier curves in this section. First, we'll look at the kind that just follows point 1. where the end - point of a segment is the same point as the start point of the next segment. This leads to poly-Béziers that are pretty hard to work with, - but they're the easiest to implement: -
-
-
- 
-
- - 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... -
-
-
- 
-
- - 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 - control point cannot satisfy the constraint that it's the reflection of its two neighbouring control points... This means that we cannot - use quadratic poly-Béziers for anything other than really, really simple shapes. And even then, they're probably the wrong choice. Cubic - curves are pretty decent, but the fact that the derivatives are linked means we can't manipulate curves as well as we might if we relaxed - the constraints a little. -
-So: let's relax the requirement a little.
-- We can change the constraint so that we still preserve the angle of the derivatives across sections (so transitions from one - section to the next will still look natural), but give up the requirement that they should also have the same vector length. - Doing so will give us a much more useful kind of poly-Bézier curve: -
-
-
- 
-
- - 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 - work. Quadratic curves really aren't very useful to work with... -
-- Finally, we also want to make sure that moving the on-curve coordinates preserves the relative positions of the associated control points. - With that, we get to the kind of curve control that you might be familiar with from applications like Photoshop, Inkscape, Blender, etc. -
-
-
- 
-
- - 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, - for instance. There is no way to preserve correct control points without a kink at the bottom point. Quadratic curves: just not that - good... -
-- A final improvement is to offer fine-level control over which points behave which, so that you can have "kinks" or individually controlled - segments when you need them, with nicely well-behaved curves for the rest of the path. Implementing that, is left as an exercise for the - reader. -
-- - Curve offsetting -
-- Perhaps you're like me, and you've been writing various small programs that use Bézier curves in some way or another, and at some point - you make the step to implementing path extrusion. But you don't want to do it pixel based; you want to stay in the vector world. You find - that extruding lines is relatively easy, and tracing outlines is coming along nicely (although junction caps and fillets are a bit of a - hassle), and then you decide to do things properly and add Bézier curves to the mix. Now you have a problem. -
-- Unlike lines, you can't simply extrude a Bézier curve by taking a copy and moving it around, because of the curvatures; rather than a - uniform thickness, you get an extrusion that looks too thin in places, if you're lucky, but more likely will self-intersect. The trick, - then, is to scale the curve, rather than simply copying it. But how do you scale a Bézier curve? -
-- Bottom line: you can't. So you cheat. We're not going to do true curve scaling, or rather curve offsetting, because - that's impossible. Instead we're going to try to generate 'looks good enough' offset curves. -
-"What do you mean, you can't? Prove it."
-- First off, when I say "you can't," what I really mean is "you can't offset a Bézier curve with another Bézier curve", not even by using - a really high order curve. You can find the function that describes the offset curve, but it won't be a polynomial, and as such it - cannot be represented as a Bézier curve, which has to be a polynomial. Let's look at why this is: -
-
- From a mathematical point of view, an offset curve O(t) is a curve such that, given our original curve B(t),
- any point on O(t) is a fixed distance d away from coordinate B(t). So let's math that:
-
-
- However, we're working in 2D, and d is a single value, so we want to turn it into a vector. If we want a point distance
- d "away" from the curve B(t) then what we really mean is that we want a point at d times the
- "normal vector" from point B(t), where the "normal" is a vector that runs perpendicular ("at a right angle") to the tangent
- at B(t). Easy enough:
-
-
- Now this still isn't very useful unless we know what the formula for N(t) is, so let's find out. N(t) runs
- perpendicular to the original curve tangent, and we know that the tangent is simply B'(t), so we could just rotate that 90
- degrees and be done with it. However, we need to ensure that N(t) has the same magnitude for every t, or the
- offset curve won't be at a uniform distance, thus not being an offset curve at all. The easiest way to guarantee this is to make sure
- N(t) always has length 1, which we can achieve by dividing B'(t) by its magnitude:
-
-
- Determining the length requires computing an arc length, and this is where things get Tricky with a capital T. First off, to compute arc
- length from some start a to end b, we must use the formula we saw earlier. Noting that "length" is usually
- denoted with double vertical bars:
-
-
- So if we want the length of the tangent, we plug in B'(t), with t = 0 as start and t = 1 as end:
-
-
- And that's where things go wrong. It doesn't even really matter what the second derivative for B(t) is, that square root is
- screwing everything up, because it turns our nice polynomials into things that are no longer polynomials.
-
- There is a small class of polynomials where the square root is also a polynomial, but they're utterly useless to us: any polynomial with - unweighted binomial coefficients has a square root that is also a polynomial. Now, you might think that Bézier curves are just fine - because they do, but they don't; remember that only the base function has binomial coefficients. That's before we - factor in our coordinates, which turn it into a non-binomial polygon. The only way to make sure the functions stay binomial is to make - all our coordinates have the same value. And that's not a curve, that's a point. We can already create offset curves for points, we call - them circles, and they have much simpler functions than Bézier curves. -
-
- So, since the tangent length isn't a polynomial, the normalised tangent won't be a polynomial either, which means
- N(t) won't be a polynomial, which means that d times N(t) won't be a polynomial, which means
- that, ultimately, O(t) won't be a polynomial, which means that even if we can determine the function for
- O(t) just fine (and that's far from trivial!), it simply cannot be represented as a Bézier curve.
-
- And that's one reason why Bézier curves are tricky: there are actually a lot of curves that cannot be represented as a Bézier - curve at all. They can't even model their own offset curves. They're weird that way. So how do all those other programs do it? Well, - much like we're about to do, they cheat. We're going to approximate an offset curve in a way that will look relatively close to what the - real offset curve would look like, if we could compute it. -
-
- So, you cannot offset a Bézier curve perfectly with another Bézier curve, no matter how high-order you make that other Bézier curve.
- However, we can chop up a curve into "safe" sub-curves (where "safe" means that all the control points are always on a single side of the
- baseline, and the midpoint of the curve at t=0.5 is roughly in the center of the polygon defined by the curve coordinates)
- and then point-scale each sub-curve with respect to its scaling origin (which is the intersection of the point normals at the start and
- end points).
-
- A good way to do this reduction is to first find the curve's extreme points, as explained in the earlier section on curve extremities, and
- use these as initial splitting points. After this initial split, we can check each individual segment to see if it's "safe enough" based
- on where the center of the curve is. If the on-curve point for t=0.5 is too far off from the center, we simply split the
- segment down the middle. Generally this is more than enough to end up with safe segments.
-
- The following graphics show off curve offsetting, and you can use the slider to control the distance at which the curve gets offset. The - curve first gets reduced to safe segments, each of which is then offset at the desired distance. Especially for simple curves, - particularly easily set up for quadratic curves, no reduction is necessary, but the more twisty the curve gets, the more the curve needs - to be reduced in order to get segments that can safely be scaled. -
-
-
- 
-
- - 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 - is large enough, this may still lead to incorrect offsets. -
-- - Graduated curve offsetting -
-
- What if we want to do graduated offsetting, starting at some distance s but ending at some other distance e?
- Well, if we can compute the length of a curve (which we can if we use the Legendre-Gauss quadrature approach) then we can also determine
- how far "along the line" any point on the curve is. With that knowledge, we can offset a curve so that its offset curve is not uniformly
- wide, but graduated between with two different offset widths at the start and end.
-
- Like normal offsetting we cut up our curve in sub-curves, and then check at which distance along the original curve each sub-curve starts
- and ends, as well as to which point on the curve each of the control points map. This gives us the distance-along-the-curve for each
- interesting point in the sub-curve. If we call the total length of all sub-curves seen prior to seeing "the current" sub-curve
- S (and if the current sub-curve is the first one, S is zero), and we call the full length of our original curve
- L, then we get the following graduation values:
-
-
-
- start: map
Sfrom interval (0,L) to interval(s,e)
- - c1:
map(<strong>S+d1</strong>, 0,L, s,e), d1 = distance along curve to projection of c1
- - c2:
map(<strong>S+d2</strong>, 0,L, s,e), d2 = distance along curve to projection of c2
- - ... -
- end:
map(<strong>S+length(subcurve)</strong>, 0,L, s,e)
-
- At each of the relevant points (start, end, and the projections of the control points onto the curve) we know the curve's normal, so - offsetting is simply a matter of taking our original point, and moving it along the normal vector by the offset distance for each point. - Doing so will give us the following result (these have with a starting width of 0, and an end width of 40 pixels, but can be controlled - with your up and down arrow keys): -
-
-
- 
-
- - - Circles and quadratic Bézier curves -
-- Circles and Bézier curves are very different beasts, and circles are infinitely easier to work with than Bézier curves. Their formula is - much simpler, and they can be drawn more efficiently. But, sometimes you don't have the luxury of using circles, or ellipses, or arcs. - Sometimes, all you have are Bézier curves. For instance, if you're doing font design, fonts have no concept of geometric shapes, they only - know straight lines, and Bézier curves. OpenType fonts with TrueType outlines only know quadratic Bézier curves, and OpenType fonts with - Type 2 outlines only know cubic Bézier curves. So how do you draw a circle, or an ellipse, or an arc? -
-You approximate.
-- We already know that Bézier curves cannot model all curves that we can think of, and this includes perfect circles, as well as ellipses, - and their arc counterparts. However, we can certainly approximate them to a degree that is visually acceptable. Quadratic and cubic curves - offer us different curvature control, so in order to approximate a circle we will first need to figure out what the error is if we try to - approximate arcs of increasing degree with quadratic and cubic curves, and where the coordinates even lie. -
-- Since arcs are mid-point-symmetrical, we need the control points to set up a symmetrical curve. For quadratic curves this means that the - control point will be somewhere on a line that intersects the baseline at a right angle. And we don't get any choice on where that will - be, since the derivatives at the start and end point have to line up, so our control point will lie at the intersection of the tangents at - the start and end point. -
-- First, let's try to fit the quadratic curve onto a circular arc. In the following sketch you can move the mouse around over a unit circle, - to see how well, or poorly, a quadratic curve can approximate the arc from (1,0) to where your mouse cursor is: -
-
-
- - 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 - our on-curve coordinates are with respect to a circular arc, given a specific start and end angle. We'll be looking at how much space - there is between the circular arc, and the quadratic curve's midpoint. -
-- 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? -
-
- 
- plotted for 0 ≤ φ ≤ π:
- |
-
- 
- plotted for 0 ≤ φ ≤ ½π:
- |
-
-
- 
-
- plotted for 0 ≤ φ ≤ ¼π:
- |
-
- We now see why the eighth circle arc looks decent, but the quarter circle arc doesn't: an error of roughly 0.06 at t=0.5 means - we're 6% off the mark... we will already be off by one pixel on a circle with pixel radius 17. Any decent sized quarter circle arc, say - with radius 100px, will be way off if approximated by a quadratic curve! For the eighth circle arc, however, the error is only roughly - 0.003, or 0.3%, which explains why it looks so close to the actual eighth circle arc. In fact, if we want a truly tiny error, like 0.001, - we'll have to contend with an angle of (rounded) 0.593667, which equates to roughly 34 degrees. We'd need 11 quadratic curves to form a - full circle with that precision! (technically, 10 and ten seventeenth, but we can't do partial curves, so we have to round up). That's a - whole lot of curves just to get a shape that can be drawn using a simple function! -
-- 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, - 1.038, 0.594 and 0.3356; in degrees, that means we can cover roughly 100 degrees (requiring four curves), 59.5 degrees (requiring six - curves), 34 degrees (requiring 11 curves), and 19.2 degrees (requiring a whopping nineteen curves). -
-- The bottom line? Quadratic curves are kind of lousy if you want circular (or elliptical, which are circles that have been - squashed in one dimension) curves. We can do better, even if it's just by raising the order of our curve once. So let's try the same thing - for cubic curves. -
-- - Circles and cubic Bézier curves -
-- In the previous section we tried to approximate a circular arc with a quadratic curve, and it mostly made us unhappy. Cubic curves are - much better suited to this task, so what do we need to do? -
-- For cubic curves, we basically want the curve to pass through three points on the circle: the start point, the mid point at "angle/2", and - the end point at "angle". We then also need to make sure the control points are such that the start and end tangent lines line up with the - circle's tangent lines at the start and end point. -
-- The first thing we can do is "guess" what the curve should look like, based on the previously outlined curve-through-three-points - procedure. This will give use a curve with correct start, mid and end points, but possibly incorrect derivatives at the start and end, - because the control points might not be in the right spot. We can then slide the control points along the lines that connect them to their - respective end point, until they effect the corrected derivative at the start and end points. However, if you look back at the section on - fitting curves through three points, the rules used were such that they optimized for a near perfect hemisphere, so using the same guess - won't be all that useful: guessing the solution based on knowing the solution is not really guessing. -
-- So have a graphical look at a "bad" guess versus the true fit, where we'll be using the bad guess and the description in the second - paragraph to derive the maths for the true fit: -
-
-
- - We see two curves here; in blue, our "guessed" curve and its control points, and in grey/black, the true curve fit, with proper control - points that were shifted in, along line between our guessed control points, such that the derivatives at the start and end points are - correct. -
-- We can already see that cubic curves are a lot better than quadratic curves, and don't look all that wrong until we go well past a quarter - circle; ⅜th starts to hint at problems, and half a circle has an obvious "gap" between the real circle and the cubic approximation. - Anything past that just looks plain ridiculous... but quarter curves actually look pretty okay! -
-So, maths time again: how okay is "okay"? Let's apply some more maths to find out.
-- Unlike for the quadratic curve, we can't use t=0.5 as our reference point because by its very nature it's one of the three points - that are actually guaranteed to lie on the circular curve. Instead, we need a different t value. If we run some analysis on the - curve we find that the actual t value at which the curve is furthest from what it should be is 0.211325 (rounded), but we don't - know "why", since finding this value involves root-finding, and is nearly impossible to do symbolically without pages and pages of math - just to express one of the possible solutions. -
-- So instead of walking you through the derivation for that value, let's simply take that t value and see what the error is for - circular arcs with an angle ranging from 0 to 2π: -
-
- 
- plotted for 0 ≤ φ ≤ 2π:
- |
-
- 
- plotted for 0 ≤ φ ≤ π:
- |
-
- 
- plotted for 0 ≤ φ ≤ ½π:
- |
-
- We see that cubic Bézier curves are much better when it comes to approximating circular arcs, with an error of less than 0.027 at the two - "bulge" points for a quarter circle (which had an error of 0.06 for quadratic curves at the mid point), and an error near 0.001 for an - eighth of a circle, so we're getting less than half the error for a quarter circle, or: at a slightly lower error, we're getting twice the - arc. This makes cubic curves quite useful! -
-- In fact, the precision of a cubic curve at a quarter circle is considered "good enough" by so many people that it's generally considered - "just fine" to use four cubic Bézier curves to fake a full circle when no circle primitives are available; generally, people won't notice - that it's not a real circle unless you also happen to overlay an actual circle, so that the difference becomes obvious. -
-- So with the error analysis out of the way, how do we actually compute the coordinates needed to get that "true fit" cubic curve? The first - observation is that we already know the start and end points, because they're the same as for the quadratic attempt: -
- -
- - But we now need to find two control points, rather than one. If we want the derivatives at the start and end point to match the circle, - then the first control point can only lie somewhere on the vertical line through S, and the second control point can only lie somewhere on - the line tangent to point E, which means: -
- -
- where "a" is some scaling factor we'll need to find the expression for, and:
- -
- using the same scaling factor, because circular arcs are symmetrical, so our approximation will need to be symmetrical, too.
-- Starting with this information, we slowly maths our way to success, but I won't lie: the maths for this is pretty trig-heavy, and it's - easy to get lost if you remember (or know!) some of the core trigonometric identities, so if you just want to see the final result just - skip past the next section! -
-Let's do this thing.
-- Unlike for the quadratic case, we need some more information in order to compute a and b, since they're no longer - dependent variables. First, we observe that the curve is symmetrical, so whatever values we end up finding for C1 will apply - to C2 as well (rotated along its tangent), so we'll focus on finding the location of C1 only. So here's where we - do something that you might not expect: we're going to ignore for a moment, because we're going to have a much easier time if we just - solve this problem with geometry first, then move to calculus to solve a much simpler problem. -
-- If we look at the triangle that is formed between our starting point, or initial guess C1 and our real C1, there's - something funny going on: if we treat the line {start,guess} as our opposite side, the line {guess,real} as our adjacent side, with - {start,real} our hypotenuse, then the angle for the corner hypotenuse/adjacent is half that of the arc we're covering. Try it: if you - place the end point at a quarter circle (pi/2, or 90 degrees), the angle in our triangle is half a quarter (pi/4, or 45 degrees). With - that knowledge, and a knowledge of what the length of any of our lines segments are (as a function), we can determine where our control - points are, and thus have everything we need to find the error distance function. Of the three lines, the one we can easiest determine - is {start,guess}, so let's find out what the guessed control point is. Again geometrically, because we have the benefit of an on-curve - t=0.5 value. -
-- The distance from our guessed point to the start point is exactly the same as the projection distance we looked at earlier. Using - t=0.5 as our point "B" in the "A,B,C" projection, then we know the length of the line segment {C,A}, since it's d1 = - {A,B} + d2 = {B,C}: -
- -
- - So that just leaves us to find the distance from t=0.5 to the baseline for an arbitrary angle φ, which is the distance from the - centre of the circle to our t=0.5 point, minus the distance from the centre to the line that runs from start point to end point. - The first is the same as the point P we found for the quadratic curve: -
- -
- - And the distance from the origin to the line start/end is another application of angles, since the triangle {origin,start,C} has known - angles, and two known sides. We can find the length of the line {origin,C}, which lets us trivially compute the coordinate for C: -
- -
- - With the coordinate C, and knowledge of coordinate B, we can determine coordinate A, and get a vector that is identical to the vector - {start,guess}: -
- -
-
- 
- - Which means we can now determine the distance {start,guessed}, which is the same as the distance {C,A}, and use that to determine the - vertical distance from our start point to our C1: -
- -
- - And after this tedious detour to find the coordinate for C1, we can find C2 fairly simply, since it's lies at - distance -C1y along the end point's tangent: -
- -
- 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
- -
- - And, because the "quarter curve" special case comes up so incredibly often, let's look at what these new control points mean for the curve - coordinates of a quarter curve, by simply filling in φ = π/2: -
- -
- Which, in decimal values, rounded to six significant digits, is:
- -
- - Of course, this is for a circle with radius 1, so if you have a different radius circle, simply multiply the coordinate by the radius you - need. And then finally, forming a full curve is now a simple a matter of mirroring these coordinates about the origin: -
-
- - - Approximating Bézier curves with circular arcs -
-- Let's look at doing the exact opposite of the previous section: rather than approximating circular arc using Bézier curves, let's - approximate Bézier curves using circular arcs. -
-- We already saw in the section on circle approximation that this will never yield a perfect equivalent, but sometimes you need circular - arcs, such as when you're working with fabrication machinery, or simple vector languages that understand lines and circles, but not much - else. -
-- The approach is fairly simple: pick a starting point on the curve, and pick two points that are further along the curve. Determine the - circle that goes through those three points, and see if it fits the part of the curve we're trying to approximate. Decent fit? Try spacing - the points further apart. Bad fit? Try spacing the points closer together. Keep doing this until you've found the "good approximation/bad - approximation" boundary, record the "good" arc, and then move the starting point up to overlap the end point we previously found. Rinse - and repeat until we've covered the entire curve. -
-- We already saw how to fit a circle through three points in the section on creating a curve from three points, - and finding the arc through those points is straight-forward: pick one of the three points as start point, pick another as an end point, - and the arc has to necessarily go from the start point, to the end point, over the remaining point. -
-So, how can we convert a Bézier curve into a (sequence of) circular arc(s)?
--
-
- Start at
t=0
- - Pick two points further down the curve at some value
m = t + nande = t + 2n
- - Find the arc that these points define -
-
- Determine how close the found arc is to the curve:
-
-
-
- Pick two additional points
e1 = t + n/2ande2 = t + n + n/2.
- -
- These points, if the arc is a good approximation of the curve interval chosen, should lie
onthe circle, so their - distance to the center of the circle should be the same as the distance from any of the three other points to the center. -
- -
- For point points, determine the (absolute) error between the radius of the circle, and the
actualdistance from the - center of the circle to the point on the curve. -
- - If this error is too high, we consider the arc bad, and try a smaller interval. -
- - Pick two additional points
- The result of this is shown in the next graphic: we start at a guaranteed failure: s=0, e=1. That's the entire curve. The midpoint is
- simply at t=0.5, and then we start performing a
- binary search.
-
-
-
- We start with
low=0,mid=0.5andhigh=1
- -
- That'll fail, so we retry with the interval halved:
{0, 0.25, 0.5}--
-
- If that arc's good, we move back up by half distance:
{0, 0.375, 0.75}.
- - However, if the arc was still bad, we move down by half the distance:
{0, 0.125, 0.25}.
-
- - If that arc's good, we move back up by half distance:
- - We keep doing this over and over until we have two arcs, in sequence, of which the first arc is good, and the second arc is bad. When we - find that pair, we've found the boundary between a good approximation and a bad approximation, and we pick the good arc. - -
- The following graphic shows the result of this approach, with a default error threshold of 0.5, meaning that if an arc is off by a - combined half pixel over both verification points, then we treat the arc as bad. This is an extremely simple error policy, but - already works really well. Note that the graphic is still interactive, and you can use your up and down arrow keys keys to increase or - decrease the error threshold, to see what the effect of a smaller or larger error threshold 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 - point we are done. Again, the following graphic allows for up and down arrow key input to increase or decrease the error threshold, so you - can see how picking a different threshold changes the number of arcs that are necessary to reasonably approximate a curve: -
-
-
- - 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 - can determine whether a coordinate lies "on" your curve really easily (simply compute the distance to each circular arc center, and if any - of those are close to the arc radii, at an angle between the arc start and end, bingo, this point can be treated as lying "on the curve"). - Another benefit is that this approximation is "linear": you can almost trivially travel along the arcs at fixed speed. You can also - trivially compute the arc length of the approximated curve (it's a bit like curve flattening). The only thing to bear in mind is that this - is a lossy equivalence: things that you compute based on the approximation are guaranteed "off" by some small value, and depending on how - much precision you need, arc approximation is either going to be super useful, or completely useless. It's up to you to decide which, - based on your application! -
-- - B-Splines -
-- No discussion on Bézier curves is complete without also giving mention of that other beast in the curve design space: B-Splines. Easily - confused to mean Bézier splines, that's not actually what they are; they are "basis function" splines, which makes a lot of difference, - and we'll be looking at those differences in this section. We're not going to dive as deep into B-Splines as we have for Bézier curves - (that would be an entire primer on its own) but we'll be looking at how B-Splines work, what kind of maths is involved in computing them, - and how to draw them based on a number of parameters that you can pick for individual B-Splines. -
-- First off: B-Splines are piecewise, - polynomial interpolation curves, where the "single curve" is built by - performing polynomial interpolation over a set of points, using a sliding window of a fixed number of points. For instance, a "cubic" - B-Spline defined by twelve points will have its curve built by evaluating the polynomial interpolation of four points, and the curve can - be treated as a lot of different sections, each controlled by four points at a time, such that the full curve consists of smoothly - connected sections defined by points {1,2,3,4}, {2,3,4,5}, ..., {8,9,10,11}, and finally {9,10,11,12}, for eight sections. -
-- What do they look like? They look like this! Tap on the graphic to add more points, and move points around to see how they map to the - spline curve drawn. -
-
- - 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 - poly-Bézier curve is actually impossible, because we start with a four point curve, and then add three more points for each section that - follows, so we can only have 4, 7, 10, 13, 16, etc. point Poly-Béziers. Similarly, while Catmull-Rom curves can grow by adding single - points, this addition of a single point introduces three implicit Bézier points. Cubic B-Splines, on the other hand, are smooth - interpolations of each possible curve involving four consecutive points, such that at any point along the curve except for our - start and end points, our on-curve coordinate is defined by four control points. -
-Consider the difference to be this:
--
-
- for Bézier curves, the curve is defined as an interpolation of points, but: -
- for B-Splines, the curve is defined as an interpolation of curves. -
In fact, let's look at that again, but this time with the base curves shown, too. Each consecutive four points define one curve:
-
-
- - 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. -
-How to compute a B-Spline curve: some maths
-
- Given a B-Spline of degree d and thus order k=d+1 (so a quadratic B-Spline is degree 2 and order 3, a cubic
- B-Spline is degree 3 and order 4, etc) and n control points P0 through Pn-1, 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
- i, which comes from our summation, and some magical parameter k. So we need to know two things: 1. what does
- N(t) do, and 2. what is that k? Let's cover both, in reverse order.
-
- The parameter k represents the "knot interval" over which a section of curve is defined. As we learned earlier, a B-Spline
- curve is itself an interpolation of curves, and we can treat each transition where a control point starts or stops influencing the total
- curvature as a "knot on the curve". Doing so for a degree d B-Spline with n control point gives us
- d + n + 1 knots, defining d + n intervals along the curve, and it is these intervals that the above
- 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
- our t value first, so that it lies in the interval bounded by knots[d] and knots[n], which are the
- start point and end point where curvature is controlled by exactly order control points. For instance, for degree 3 (=order
- 4) and 7 control points, with knot vector [1,2,3,4,5,6,7,8,9,10,11], we map t from [the interval 0,1] to the interval [4,8],
- and then use that value in the functions above, instead.
-
Can we simplify that?
-We can, yes.
-
- People far smarter than us have looked at this work, and two in particular —
- Maurice Cox and
- 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
- matter of computing pretty simple arithmetical statements, with some caching of results so we can refer to them as we recurve. While the
- 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.
-
- Thanks to Cox and de Boor, we can compute points on a B-Spline pretty easily using the same kind of linear interpolation we saw in de
- 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
- themselves a mixture of etc. etc. We simply keep expanding our terms until we reach the stop conditions, and then sum everything back up.
- It's really quite elegant.
-
- One thing we need to keep in mind is that we're working with a spline that is constrained by its control points, so even though the
- d(..., k) values are zero or one at the lowest level, they are really "zero or one, times their respective control point", so
- in the next section you'll see the algorithm for running through the computation in a way that starts with a copy of the control point
- vector and then works its way up to that single point, rather than first starting "on the left", working our way "to the right" and then
- summing back up "to the left". We can just start on the right and work our way left immediately.
-
Running the computation
-
- Unlike the de Casteljau algorithm, where the t value stays the same at every iteration, for B-Splines that is not the case,
- and so we end having to (for each point we evaluate) run a fairly involving bit of recursive computation. The algorithm is discussed on
- this Michigan Tech page, but an easier to read version
- is implemented by b-spline.js, so we'll look at its code.
-
- Given an input value t, we first map the input to a value from the domain [0,1] to the domain
- [knots[degree], knots[knots.length - 1 - degree]. Then, we find the section number s that this mapped
- t value lies on:
-
| 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):
-
| 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
- v[i-1] don't try to use an array index that doesn't exist)
-
Open vs. closed paths
-
- Much like poly-Béziers, B-Splines can be either open, running from the first point to the last point, or closed, where the first and last
- point are the same coordinate. However, because B-Splines are an interpolation of curves, not just points, we can't simply make the first
- and last point the same, we need to link as many points as are necessary to form "a curve" that the spline performs interpolation with. As
- such, for an order d B-Spline, we need to make the first and last d points the same. This is of course hardly
- more work than before (simply append points.splice(0,d) to points) but it's important to remember that you need
- more than just a single point.
-
- Of course if we want to manipulate these kind of curves we need to make sure to mark them as "closed" so that we know the coordinate for
- points[0] and points[n-k] etc. don't just happen to have the same x/y values, but really are the same
- coordinate, so that manipulating one will equally manipulate the other, but programming generally makes this really easy by storing
- references to points, rather than copies (or other linked values such as coordinate weights, discussed in the NURBS section) rather than
- separate coordinate objects.
-
Manipulating the curve through the knot vector
-- The most important thing to understand when it comes to B-Splines is that they work because of the concept of a knot vector. As - mentioned above, knots represent "where individual control points start/stop influencing the curve", but we never looked at the - values that go in the knot vector. If you look back at the N() and a() functions, you see that interpolations are based on - intervals in the knot vector, rather than the actual values in the knot vector, and we can exploit this to do some pretty interesting - things with clever manipulation of the knot vector. Specifically there are four things we can do that are worth looking at: -
--
-
- we can use a uniform knot vector, with equally spaced intervals, -
- we can use a non-uniform knot vector, without enforcing equally spaced intervals, -
- we can collapse sequential knots to the same value, locally lowering curve complexity using "null" intervals, and -
- - we can form a special case non-uniform vector, by combining (1) and (3) to for a vector with collapsed start and end knots, with a - uniform vector in between. - -
Uniform B-Splines
-- The most straightforward type of B-Spline is the uniform spline. In a uniform spline, the knots are distributed uniformly over the entire - curve interval. For instance, if we have a knot vector of length twelve, then a uniform knot vector would be [0,1,2,3,...,9,10,11]. Or - [4,5,6,...,13,14,15], which defines the same intervals, or even [0,2,3,...,18,20,22], which also defines - the same intervals, just scaled by a constant factor, which becomes normalised during interpolation and so does not contribute to - the curvature. -
-
-
- - 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. -
-
- The problem with uniform knot vectors is that, as we need order control points before we have any curve with which we can
- perform interpolation, the curve does not "start" at the first point, nor "ends" at the last point. Instead there are "gaps". We can get
- rid of these, by being clever about how we apply the following uniformity-breaking approach instead...
-
Reducing local curve complexity by collapsing intervals
-
- Collapsing knot intervals, by making two or more consecutive knots have the same value, allows us to reduce the curve complexity in the
- sections that are affected by the knots involved. This can have drastic effects: for every interval collapse, the curve order goes down,
- and curve continuity goes down, to the point where collapsing order knots creates a situation where all continuity is lost
- and the curve "kinks".
-
-
- 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 - curve not starting and ending where we'd kind of like it to: -
-
- For any curve of degree D with control points N, we can define a knot vector of length N+D+1 in
- which the values 0 ... D+1 are the same, the values D+1 ... N+1 follow the "uniform" pattern, and the values
- N+1 ... N+D+1 are the same again. For example, a cubic B-Spline with 7 control points can have a knot vector
- [0,0,0,0,1,2,3,4,4,4,4], or it might have the "identical" knot vector [0,0,0,0,2,4,6,8,8,8,8], etc. Again, it is the relative differences
- that determine the curve shape.
-
-
- 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
- pick specific knot intervals, there is nothing particularly interesting going on. There is one constraint to the knot vector, other than
- that any value knots[k+1] should be greater than or equal to knots[k].
-
One last thing: Rational B-Splines
-- While it is true that this section on B-Splines is running quite long already, there is one more thing we need to talk about, and that's - "Rational" splines, where the rationality applies to the "ratio", or relative weights, of the control points themselves. By introducing a - ratio vector with weights to apply to each control point, we greatly increase our influence over the final curve shape: the more weight a - control point carries, the closer to that point the spline curve will lie, a bit like turning up the gravity of a control point, just like - for rational Bézier curves. -
-
-
- - 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, - the capital S actually just stands for "spline", but a lot of people mistakenly treat it as if it is, so now you know better). NURBS is an - important type of curve in computer-facilitated design, used a lot in 3D modelling (typically as NURBS surfaces) as well as in - arbitrary-precision 2D design due to the level of control a NURBS curve offers designers. -
-- While a true non-uniform rational B-Spline would be hard to work with, when we talk about NURBS we typically mean the Open-Uniform - Rational B-Spline, or OURBS, but that doesn't roll off the tongue nearly as nicely, and so remember that when people talk about NURBS, - they typically mean open-uniform, which has the useful property of starting the curve at the first control point, and ending it at the - last. -
-Extending our implementation to cover rational splines
-- The algorithm for working with Rational B-Splines is virtually identical to the regular algorithm, and the extension to work in the - control point weights is fairly simple: we extend each control point from a point in its original number of dimensions (2D, 3D, etc.) to - one dimension higher, scaling the original dimensions by the control point's weight, and then assigning that weight as its value for the - extended dimension. -
-For example, a 2D point (x,y) with weight w becomes a 3D point (w * x, w * y, w).
- We then run the same algorithm as before, which will automatically perform weight interpolation in addition to regular coordinate - interpolation, because all we've done is pretended we have coordinates in a higher dimension. The algorithm doesn't really care about how - many dimensions it needs to interpolate. -
-
- In order to recover our "real" curve point, we take the final result of the point generation algorithm, and "unweigh" it: we take the
- final point's derived weight w' and divide all the regular coordinate dimensions by it, then throw away the weight
- information.
-
- Based on our previous example, we take the final 3D point (x', y', w'), which we then turn back into a 2D point by computing
- (x'/w', y'/w'). And that's it, we're done!
-
-previous
- Comments and questions
-
-- First off, if you enjoyed this book, or you simply found it useful for something you were trying to get done, and you were wondering how - to let me know you appreciated this book, you have two options: you can either head on over to the - Patreon page for this book, or if you prefer to make a one-time donation, head on over to - the buy Pomax a coffee page. This - work has grown from a small primer to a 70-plus print-page-equivalent reader on the subject of Bézier curves over the years, and a lot of - coffee went into the making of it. I don't regret a minute I spent on writing it, but I can always do with some more coffee to keep on - writing. -
-With that said, on to the comments!
- -