info/BezierInfo-2
1
0
Fork 0
mirror of https://github.com/Pomax/BezierInfo-2.git synced 2025-02-22 16:50:13 +01:00
Code Issues Releases Wiki Activity
BezierInfo-2/index.html
Pomax b4feeba376 let's do this
2020-08-08 16:35:45 -07:00

7859 lines
352 KiB
HTML
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

<!DOCTYPE html>
<html lang="en-GB">
<head>
<meta charset="utf-8" />
<meta name="viewport" content="width=device-width, initial-scale=1" />
<title>A Primer on Bézier Curves</title>
<link rel="shortcut icon" href="favicon.png" type="image/png" />
<!-- And a slew of SEO related meta elements, because being discoverable is important -->
<meta
name="description"
content="A detailed explanation of Bézier curves, and how to do the many things that we commonly want to do with them."
/>
<!-- opengraph information -->
<meta property="og:title" content="A Primer on Bézier Curves" />
<meta
property="og:image"
content="https://pomax.github.io/bezierinfo/images/og-image.png"
/>
<meta property="og:type" content="text" />
<meta property="og:url" content="https://pomax.github.io/bezierinfo" />
<meta
property="og:description"
content="A detailed explanation of Bézier curves, and how to do the many things that we commonly want to do with them."
/>
<meta property="og:locale" content="en-GB}}" />
<meta property="og:type" content="article" />
<meta property="og:published_time" content="2013-06-13 12:00:00" />
<meta property="og:author" content="Mike 'Pomax' Kamermans" />
<meta property="og:section" content="Bézier Curves" />
<meta property="og:tag" content="Bézier Curves" />
<!-- twitter card information -->
<meta name="twitter:card" content="summary" />
<meta name="twitter:site" content="@TheRealPomax" />
<meta name="twitter:creator" content="@TheRealPomax" />
<meta
name="twitter:image"
content="https://pomax.github.io/bezierinfo/images/og-image.png"
/>
<meta name="twitter:url" content="https://pomax.github.io/bezierinfo" />
<meta
name="twitter:description"
content="A detailed explanation of Bézier curves, and how to do the many things that we commonly want to do with them."
/>
<!-- my own referral/page hit tracker, because Google knows enough -->
<script src="./lib/site/referrer.js" type="module" async></script>
<!-- the part that makes interactive graphics work: an HTML5 <graphics-element> custom element -->
<script
src="./lib/custom-element/graphics-element.js"
type="module"
async
defer
></script>
<link rel="stylesheet" href="./lib/custom-element/graphics-element.css" />
<!-- page styling -->
<link rel="stylesheet" href="placeholder-style.css" />
</head>
<body>
<header>
<h1>A Primer on Bézier Curves</h1>
<h2>
A free, online book for when you really need to know how to do Bézier
things.
</h2>
<span>Read this in your own language:</span>
<ul>
<li><a href="./index.html">English</a></li>
<li><a href="./ja-JP/index.html">日本語</a></li>
<li><a href="./zh-CN/index.html">中文</a></li>
</ul>
<span
>Don't see your language listed?
<a href="https://github.com/Pomax/BezierInfo-2/wiki/localize"
>Help translate this content!</a
></span
>
<nav aria-labelledby="toc">
<h1 id="toc">Table of Contents</h1>
<ol>
<li><a href="#introduction">A lightning introduction</a></li>
<li><a href="#whatis">So what makes a Bézier Curve?</a></li>
<li><a href="#explanation">The mathematics of Bézier curves</a></li>
<li><a href="#control">Controlling Bézier curvatures</a></li>
<li>
<a href="#weightcontrol"
>Controlling Bézier curvatures, part 2: Rational Béziers</a
>
</li>
<li><a href="#extended">The Bézier interval [0,1]</a></li>
<li><a href="#matrix">Bézier curvatures as matrix operations</a></li>
<li><a href="#decasteljau">de Casteljau's algorithm</a></li>
<li><a href="#flattening">Simplified drawing</a></li>
<li><a href="#splitting">Splitting curves</a></li>
<li><a href="#matrixsplit">Splitting curves using matrices</a></li>
<li><a href="#reordering">Lowering and elevating curve order</a></li>
<li><a href="#derivatives">Derivatives</a></li>
<li><a href="#pointvectors">Tangents and normals</a></li>
<li><a href="#pointvectors3d">Working with 3D normals</a></li>
<li><a href="#components">Component functions</a></li>
<li><a href="#extremities">Finding extremities: root finding</a></li>
<li><a href="#boundingbox">Bounding boxes</a></li>
<li><a href="#aligning">Aligning curves</a></li>
<li><a href="#tightbounds">Tight boxes</a></li>
<li><a href="#inflections">Curve inflections</a></li>
<li><a href="#canonical">Canonical form (for cubic curves)</a></li>
<li><a href="#yforx">Finding Y, given X</a></li>
<li><a href="#arclength">Arc length</a></li>
<li><a href="#arclengthapprox">Approximated arc length</a></li>
<li><a href="#curvature">Curvature of a curve</a></li>
<li>
<a href="#tracing">Tracing a curve at fixed distance intervals</a>
</li>
<li><a href="#intersections">Intersections</a></li>
<li><a href="#curveintersection">Curve/curve intersection</a></li>
<li><a href="#abc">The projection identity</a></li>
<li><a href="#moulding">Manipulating a curve</a></li>
<li><a href="#pointcurves">Creating a curve from three points</a></li>
<li><a href="#curvefitting">Curve fitting</a></li>
<li>
<a href="#catmullconv">Bézier curves and Catmull-Rom curves</a>
</li>
<li>
<a href="#catmullmoulding"
>Creating a Catmull-Rom curve from three points</a
>
</li>
<li><a href="#polybezier">Forming poly-Bézier curves</a></li>
<li><a href="#shapes">Boolean shape operations</a></li>
<li>
<a href="#projections">Projecting a point onto a Bézier curve</a>
</li>
<li><a href="#offsetting">Curve offsetting</a></li>
<li><a href="#graduatedoffset">Graduated curve offsetting</a></li>
<li><a href="#circles">Circles and quadratic Bézier curves</a></li>
<li><a href="#circles_cubic">Circles and cubic Bézier curves</a></li>
<li>
<a href="#arcapproximation"
>Approximating Bézier curves with circular arcs</a
>
</li>
<li><a href="#bsplines">B-Splines</a></li>
<li><a href="#comments">Comments and questions</a></li>
</ol>
</nav>
</header>
<main>
<section id="preface">
<section id="preface">
<h1>Preface</h1>
<p>
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.
</p>
<p>
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.
</p>
<p>
They're named after
<a href="https://en.wikipedia.org/wiki/Pierre_B%C3%A9zier"
>Pierre Bézier</a
>, 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
<a href="https://en.wikipedia.org/wiki/Paul_de_Casteljau"
>Paul de Casteljau</a
>
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
<a href="https://en.wikipedia.org/wiki/Sergei_Natanovich_Bernstein"
>Sergei Natanovich Bernstein</a
>, whose publications on them date back at least as far as 1912.
</p>
<p>
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.
</p>
<p>
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 when you do with these
curves? That's what this page is for. Prepare to be mathed!
</p>
<div class="print">
<h2>PS: buy me a coffee?</h2>
<p>
If you enjoyed this book enough to print it out, you might be
wondering if there is some way to give something back. To answer
that question: you can always buy me a coffee, however-much a
coffee is where you live. Or, if you want to pay a fair price for
this book, you can buy me a really expensive coffee =)
</p>
<p>
This book has grown over the years from a short primer to a 100+
print-page-equivalent ebook on the subject of Bézier curves, 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! Please visit
<a href="https://pomax.github.io/bezierinfo"
>https://pomax.github.io/bezierinfo</a
>
and click on the link in the "Help support the book" preface
section to donate some coffee money.
</p>
</div>
<p>
—Pomax (or in the tweetworld,
<a href="https://twitter.com/TheRealPomax">@TheRealPomax</a>)
</p>
<div class="note">
<h2>Virtually all Bézier graphics are interactive.</h2>
<p>
This page uses interactive examples, relying heavily on
<a href="http://pomax.github.io/bezierjs">Bezier.js</a>, as well
as maths formulae which are typeset into SVG using the
<a href="https://ctan.org/pkg/xetex">XeLaTeX</a> typesetting
system and
<a href="https://github.com/dawbarton/pdf2svg">pdf2svg</a> by
<a href="http://www.cityinthesky.co.uk/">David Barton</a>.
</p>
<h2>This book is open source.</h2>
<p>
This book is an open source software project, and lives on two
github repositories. The first is
<a href="https://github.com/pomax/bezierinfo"
>https://github.com/pomax/bezierinfo</a
>
and is the purely-for-presentation version you are viewing right
now. The other repository is
<a href="https://github.com/pomax/BezierInfo-2"
>https://github.com/pomax/BezierInfo-2</a
>, which is the development version, housing all the HTML,
JavaScript, and CSS. You can fork either of these, and pretty much
do with them as you please, except for passing it off as your own
work wholesale, of course =)
</p>
<h2>How complicated is the maths going to be?</h2>
<p>
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 <em>far</em> 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.
</p>
<h2>Questions, comments:</h2>
<p>
If you have suggestions for new sections, hit up the
<a href="https://github.com/pomax/BezierInfo-2/issues"
>Github issue tracker</a
>
(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.
</p>
<h2>Help support the book!</h2>
<p>
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
<a href="https://patreon.com/bezierinfo">Patreon page</a> for this
book, or if you prefer to make a one-time donation, head on over
to the
<a
href="https://www.paypal.com/cgi-bin/webscr?cmd=_s-xclick&hosted_button_id=QPRDLNGDANJSW"
>buy Pomax a coffee</a
>
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!
</p>
</div>
</section>
</section>
<section id="changelog">
<h1>What's new?</h1>
<p>
This primer is a living document, and so depending on when you last
look at it, there may be new content. Click the following link to
expand this section to have a look at what got added, when.
</p>
<section>
<h2>August 2020</h2>
<ul>
<li>
Completely redid all the code for the site, except for chapter
text: the Primer is now finally a normal web page again.
</li>
</ul>
<h2>June 2020</h2>
<ul>
<li>Added automatic CI/CD using Github Actions</li>
</ul>
<h2>January 2020</h2>
<ul>
<li>Added reset buttons to all graphics</li>
<li>Updated to preface to correctly describe the on-page maths</li>
<li>
Fixed the Catmull-Rom section because it had glaring maths errors
</li>
</ul>
<h2>August 2019</h2>
<ul>
<li>Added a section on (plain) rational Bezier curves</li>
<li>Improved the Graphic component to allow for sliders</li>
</ul>
<h2>December 2018</h2>
<ul>
<li>Added a section on curvature and calculating kappa.</li>
<li>
Added a Patreon page! Head on over to
https://patreon.com/bezierinfo to help support this site!
</li>
</ul>
<h2>August 2018</h2>
<ul>
<li>
Added a section on finding a curve's y, if all you have is the x
coordinate.
</li>
</ul>
<h2>July 2018</h2>
<ul>
<li>
Rewrote the 3D normals section, implementing and explaining
Rotation Minimising Frames.
</li>
<li>
Updated the section on curve order raising/lowering, showing how
to get a least-squares optimized lower order curve.
</li>
<li>
(Finally) updated 'npm test' so that it automatically rebuilds
when files are changed while the dev server is running.
</li>
</ul>
<h2>June 2018</h2>
<ul>
<li>Added a section on direct curve fitting.</li>
<li>Added source links for all graphics.</li>
<li>Added this "What's new?" section.</li>
</ul>
<h2>April 2017</h2>
<ul>
<li>Added a section on 3d normals.</li>
<li>
Added live-updating for the social link buttons, so they always
link to the specific section you're reading.
</li>
</ul>
<h2>February 2017</h2>
<ul>
<li>Finished rewriting the entire codebase for localization.</li>
</ul>
<h2>January 2016</h2>
<ul>
<li>Added a section to explain the Bezier interval.</li>
<li>Rewrote the Primer as a React application.</li>
</ul>
<h2>December 2015</h2>
<ul>
<li>
Set up the split repository between BezierInfo-2 as development
repository, and bezierinfo as live page.
</li>
<li>
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.
</li>
</ul>
<h2>May 2015</h2>
<ul>
<li>
Switched over to pure JS rather than
Processing-through-Processing.js
</li>
<li>
Added Cardano's algorithm for finding the roots of a cubic
polynomial.
</li>
</ul>
<h2>April 2015</h2>
<ul>
<li>Added a section on arc length approximations.</li>
</ul>
<h2>February 2015</h2>
<ul>
<li>Added a section on the canonical cubic Bezier form.</li>
</ul>
<h2>November 2014</h2>
<ul>
<li>Switched to HTTPS.</li>
</ul>
<h2>July 2014</h2>
<ul>
<li>Added the section on arc approximation.</li>
</ul>
<h2>April 2014</h2>
<ul>
<li>Added the section on Catmull-Rom fitting.</li>
</ul>
<h2>November 2013</h2>
<ul>
<li>Added the section on Catmull-Rom / Bezier conversion.</li>
<li>Added the section on Bezier cuves as matrices</li>
</ul>
<h2>April 2013</h2>
<ul>
<li>Added a section on poly-Beziers.</li>
<li>Added a section on boolean shape operations.</li>
</ul>
<h2>March 2013</h2>
<ul>
<li>First drastic rewrite</li>
<li>Added sections on circle approximations.</li>
<li>Added a section on projecting a point onto a curve.</li>
<li>Added a section on tangents and normals.</li>
<li>Added Legendre-Gauss numerical data tables.</li>
</ul>
<h2>October 2011</h2>
<ul>
<li>
First commit for the https://pomax.github.io/bezierinfo site,
based on the pre-Primer webpage that covered the basics of Bezier
curves in HTML with Processing.js examples.
</li>
</ul>
</section>
</section>
<section id="chapters">
<section id="introduction">
<h1>A lightning introduction</h1>
<p>
Let's start with the good stuff: when we're talking about Bézier
curves, we're talking about the things that you can see in the
following graphics. They run from some start point to some end
point, with their curvature influenced by one or more "intermediate"
control points. Now, because all the graphics on this page are
interactive, go manipulate those curves a bit: click-drag the
points, and see how their shape changes based on what you do.
</p>
<div class="figure">
<graphics-element
title="A quadratic Bézier curve"
width="275"
height="275"
src="./chapters/introduction/quadratic.js"
>
<fallback-image>
<img
width="275px"
height="275px"
src="./images/chapters/introduction/quadratic.png"
loading="lazy"
/>
Scripts are disabled. Showing fallback image.
</fallback-image></graphics-element
>
<graphics-element
title="A cubic Bézier curve"
width="275"
height="275"
src="./chapters/introduction/cubic.js"
>
<fallback-image>
<img
width="275px"
height="275px"
src="./images/chapters/introduction/cubic.png"
loading="lazy"
/>
Scripts are disabled. Showing fallback image.
</fallback-image></graphics-element
>
</div>
<p>
These curves are used a lot in computer aided design and computer
aided manufacturing (CAD/CAM) applications, as well as in graphic
design programs like Adobe Illustrator and Photoshop, Inkscape,
GIMP, etc. and in graphic technologies like scalable vector graphics
(SVG) and OpenType fonts (TTF/OTF). A lot of things use Bézier
curves, so if you want to learn more about them... prepare to get
your learn on!
</p>
</section>
<section id="whatis">
<h1>So what makes a Bézier Curve?</h1>
<p>
Playing with the points for curves may have given you a feel for how
Bézier curves behave, but what <em>are</em> 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:
</p>
<p>
Bézier curves are the result of
<a href="https://en.wikipedia.org/wiki/Linear_interpolation"
>linear interpolations</a
>. 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".
</p>
<p>
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:
</p>
<img
class="LaTeX SVG"
src="images/latex/b5aa26284ba3df74970a95cb047a841d.svg"
width="501px"
height="103px"
loading="lazy"
/>
<p>
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:
</p>
<graphics-element
title="Linear Interpolation leading to Bézier curves"
width="825"
height="275"
src="./chapters/whatis/interpolation.js"
>
<fallback-image>
<img
width="825px"
height="275px"
src="./images/chapters/whatis/interpolation.png"
loading="lazy"
/>
Scripts are disabled. Showing fallback image.
</fallback-image></graphics-element
>
<p>And that brings us to the complicated maths: calculus.</p>
<p>
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).
</p>
<p>
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.
</p>
</section>
<section id="explanation">
<h1>The mathematics of Bézier curves</h1>
<p>
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 <strong>single</strong> 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.
</p>
<p>
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 <i>x</i>, to some other value, using some
kind of number manipulation:
</p>
<img
class="LaTeX SVG"
src="images/latex/1caef9931f954e32eae5067b732c1018.svg"
width="96px"
height="17px"
loading="lazy"
/>
<p>
The notation <i>f(x)</i> is the standard way to show that it's a
function (by convention called <i>f</i> if we're only listing one)
and its output changes based on one variable (in this case,
<i>x</i>). Change <i>x</i>, and the output for <i>f(x)</i> changes.
</p>
<p>
So far, so good. Now, let's look at parametric functions, and how
they cheat. Let's take the following two functions:
</p>
<img
class="LaTeX SVG"
src="images/latex/0f5cffd58e864fec6739a57664eb8cbd.svg"
width="93px"
height="36px"
loading="lazy"
/>
<p>
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 <i>a</i>, we're not going to
change the output value for <i>f(b)</i>, since <i>a</i> 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:
</p>
<img
class="LaTeX SVG"
src="images/latex/066a910ae6aba69c40a338320759cdd1.svg"
width="100px"
height="40px"
loading="lazy"
/>
<p>
Multiple functions, but only one variable. If we change the value
for <i>t</i>, we change the outcome of both
<i>f<sub>a</sub>(t)</i> and <i>f<sub>b</sub>(t)</i>. You might
wonder how that's useful, and the answer is actually pretty simple:
if we change the labels <i>f<sub>a</sub>(t)</i> and
<i>f<sub>b</sub>(t)</i> with what we usually mean with them for
parametric curves, things might be a lot more obvious:
</p>
<img
class="LaTeX SVG"
src="images/latex/4cf6fb369841e2c5d36e5567a8db4306.svg"
width="77px"
height="40px"
loading="lazy"
/>
<p>
There we go. <i>x</i>/<i>y</i> coordinates, linked through some
mystery value <i>t</i>.
</p>
<p>
So, parametric curves don't define a <i>y</i> coordinate in terms of
an <i>x</i> coordinate, like normal functions do, but they instead
link the values to a "control" variable. If we vary the value of
<i>t</i>, then with every change we get <strong>two</strong> values,
which we can use as (<i>x</i>,<i>y</i>) coordinates in a graph. The
above set of functions, for instance, generates points on a circle:
We can range <i>t</i> from negative to positive infinity, and the
resulting (<i>x</i>,<i>y</i>) coordinates will always lie on a
circle with radius 1 around the origin (0,0). If we plot it for
<i>t</i> from 0 to 5, we get this (use your up and down arrow keys
to change the plot end value):
</p>
<graphics-element
title="A (partial) circle: x=sin(t), y=cos(t)"
width="275"
height="275"
src="./chapters/explanation/circle.js"
>
<fallback-image>
<img
width="275px"
height="275px"
src="./images/chapters/explanation/circle.png"
loading="lazy"
/>
Scripts are disabled. Showing fallback image.
</fallback-image></graphics-element
>
<p>
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
<i>x</i> and <i>y</i> 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 <i>x</i> and <i>y</i> outputs. So
what are binomial polynomials?
</p>
<p>
You may remember polynomials from high school. They're those sums
that look like this:
</p>
<img
class="LaTeX SVG"
src="images/latex/bb06cb82d372f822a7b35e661502bd72.svg"
width="213px"
height="20px"
loading="lazy"
/>
<p>
If the highest order term they have is <i>x³</i>, they're called
"cubic" polynomials; if it's <i>x²</i>, it's a "square" polynomial;
if it's just <i>x</i>, it's a line (and if there aren't even any
terms with <i>x</i> it's not a polynomial!)
</p>
<p>
Bézier curves are polynomials of <i>t</i>, rather than <i>x</i>,
with the value for <i>t</i> being fixed between 0 and 1, with
coefficients <i>a</i>, <i>b</i> etc. taking the "binomial" form,
which sounds fancy but is actually a pretty simple description for
mixing values:
</p>
<img
class="LaTeX SVG"
src="images/latex/2adc12d0cff01d40d9e1702014a7dc19.svg"
width="367px"
height="64px"
loading="lazy"
/>
<p>
I know what you're thinking: that doesn't look too simple! But if we
remove <i>t</i> and add in "times one", things suddenly look pretty
easy. Check out these binomial terms:
</p>
<img
class="LaTeX SVG"
src="images/latex/9c18f76e76cf684ecd217ad8facc2e93.svg"
width="184px"
height="87px"
loading="lazy"
/>
<p>
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
<a href="https://en.wikipedia.org/wiki/Pascal%27s_triangle"
>Pascal's triangle</a
>. Now <i>that's</i> easy to remember.
</p>
<p>
There's an equally simple way to figure out how the polynomial terms
work: if we rename <i>(1-t)</i> to <i>a</i> and <i>t</i> to
<i>b</i>, and remove the weights for a moment, we get this:
</p>
<img
class="LaTeX SVG"
src="images/latex/449850dead8abbdd11cd4aec1bac082e.svg"
width="903px"
height="63px"
loading="lazy"
/>
<p>
It's basically just a sum of "every combination of <i>a</i> and
<i>b</i>", progressively replacing <i>a</i>'s with <i>b</i>'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:
</p>
<img
class="LaTeX SVG"
src="images/latex/9a6d17c362980775f1425d0d2ad9a36a.svg"
width="300px"
height="55px"
loading="lazy"
/>
<p>
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 ...=&lt;value&gt; and
ending at the value listed on top of the Σ).
</p>
<div class="howtocode">
<h3>How to implement the basis function</h3>
<p>
We could naively implement the basis function as a mathematical
construct, using the function as our guide, like this:
</p>
<pre><code>function Bezier(n,t):
sum = 0
for(k=0; k&lt;n; k++):
sum += n!/(k!*(n-k)!) * (1-t)^(n-k) * t^(k)
return sum</code></pre>
<p>
I say we could, because we're not going to: the factorial function
is <em>incredibly</em> 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.
</p>
<p>
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:
</p>
<pre><code>lut = [ [1], // n=0
[1,1], // n=1
[1,2,1], // n=2
[1,3,3,1], // n=3
[1,4,6,4,1], // n=4
[1,5,10,10,5,1], // n=5
[1,6,15,20,15,6,1]] // n=6
function binomial(n,k):
while(n &gt;= lut.length):
s = lut.length
nextRow = new array(size=s+1)
nextRow[0] = 1
for(i=1, prev=s-1; i&lt;s; i++):
nextRow[i] = lut[prev][i-1] + lut[prev][i]
nextRow[s] = 1
lut.add(nextRow)
return lut[n][k]</code></pre>
<p>
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 <i>n/k</i> pair is requested that isn't in
the LUT yet, we expand it first. Our basis function now looks like
this:
</p>
<pre><code>function Bezier(n,t):
sum = 0
for(k=0; k&lt;=n; k++):
sum += binomial(n,k) * (1-t)^(n-k) * t^(k)
return sum</code></pre>
<p>
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:
</p>
<pre><code>function Bezier(2,t):
t2 = t * t
mt = 1-t
mt2 = mt * mt
return mt2 + 2*mt*t + t2
function Bezier(3,t):
t2 = t * t
t3 = t2 * t
mt = 1-t
mt2 = mt * mt
mt3 = mt2 * mt
return mt3 + 3*mt2*t + 3*mt*t2 + t3</code></pre>
<p>And now we know how to program the basis function. Excellent.</p>
</div>
<p>
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.
</p>
</section>
<section id="control">
<h1>Controlling Bézier curvatures</h1>
<p>
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.
</p>
<p>
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 or
click-drag to see the interpolation percentages for each
curve-defining point at a specific <i>t</i> value.
</p>
<div className="figure">
<Graphic
inline="{true}"
title="Quadratic interpolations"
draw="{this.drawQuadraticLerp}"
/>
<Graphic
inline="{true}"
title="Cubic interpolations"
draw="{this.drawCubicLerp}"
/>
<Graphic
inline="{true}"
title="15th degree interpolations"
draw="{this.draw15thLerp}"
/>
</div>
<p>
Also shown is the interpolation function for a 15<sup>th</sup> 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.
</p>
<p>
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:
</p>
<img
class="LaTeX SVG"
src="images/latex/14cb9fbbaae9e7d87ae6bef3ea7a782e.svg"
width="352px"
height="55px"
loading="lazy"
/>
<p>
That looks complicated, but as it so happens, the "weights" are
actually just the coordinate values we want our curve to have: for
an <i>n<sup>th</sup></i> order curve, w<sub>0</sub> is our start
coordinate, w<sub>n</sub> is our last coordinate, and everything in
between is a controlling coordinate. Say we want a cubic curve that
starts at (120,160), is controlled by (35,200) and (220,260) and
ends at (220,40), we use this Bézier curve:
</p>
<img
class="LaTeX SVG"
src="images/latex/e4759ccbc6fd2b21dd5ac7e36c9399fd.svg"
width="692px"
height="40px"
loading="lazy"
/>
<p>Which gives us the curve we saw at the top of the article:</p>
<Graphic
title="Our cubic Bézier curve"
setup="{this.drawCubic}"
draw="{this.drawCurve}"
/>
<p>
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.
</p>
<div className="howtocode">
<h3>How to implement the weighted basis function</h3>
<p>
Given that we already know how to implement basis function, adding
in the control points is remarkably easy:
</p>
<pre><code>function Bezier(n,t,w[]):
sum = 0
for(k=0; k&lt;=n; k++):
sum += w[k] * binomial(n,k) * (1-t)^(n-k) * t^(k)
return sum</code></pre>
<p>And now for the extremely optimized versions:</p>
<pre><code>function Bezier(2,t,w[]):
t2 = t * t
mt = 1-t
mt2 = mt * mt
return w[0]*mt2 + w[1]*2*mt*t + w[2]*t2
function Bezier(3,t,w[]):
t2 = t * t
t3 = t2 * t
mt = 1-t
mt2 = mt * mt
mt3 = mt2 * mt
return w[0]*mt3 + 3*w[1]*mt2*t + 3*w[2]*mt*t2 + w[3]*t3</code></pre>
<p>And now we know how to program the weighted basis function.</p>
</div>
</section>
<section id="weightcontrol">
<h1>Controlling Bézier curvatures, part 2: Rational Béziers</h1>
<p>
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.
</p>
<p>
Adding these ratio values to the regular Bézier curve function is
fairly easy. Where the regular function is the following:
</p>
<img
class="LaTeX SVG"
src="images/latex/02457b19087540dfb144978419524a85.svg"
width="276px"
height="41px"
loading="lazy"
/>
<p>The function for rational Bézier curves has two more terms:</p>
<img
class="LaTeX SVG"
src="images/latex/d0c52423259ba91a5426324338c83306.svg"
width="587px"
height="44px"
loading="lazy"
/>
<p>
In this, the first new term represents the ratio for each
coordinate, as a multiplication factor. If our ratio values are [1,
0.5, 0.5, 1] then <code>ratio<sub>0</sub> = 1</code>,
<code>ratio<sub>1</sub> = 0.5</code>, and so on. The second term
then corrects for all those multiplications by dividing the total
value by the "basis" value of the Bézier curve, i.e. the value we
get by computing the curve without any weighting at (but
<em>with</em> ratios):
</p>
<img
class="LaTeX SVG"
src="images/latex/03df4af5636ff82d3e964de0b9727f23.svg"
width="271px"
height="41px"
loading="lazy"
/>
<p>So what does this actually do?</p>
<p>
Let's look at the effect of "rationalising" our Bézier curves by
interacting with the curve and ratios. The following graphic shows
the curve from the previous section, enriched with ratio factors:
the closer to zero we set one or more terms, the less relative
influence the associated coordinates exert on the curve (and of
course the higher we set them, the more influence they have).
</p>
<p>
&lt;Graphic title="Our rational cubic Bézier curve"
setup={this.drawCubic} draw={this.drawCurve} sliders={[ {min:0,
max:2, value:1, step:0.01, label:<code>ratio 1</code>}, {min:0,
max:2, value:1, step:0.01, label:<code>ratio 2</code>}, {min:0,
max:2, value:1, step:0.01, label:<code>ratio 3</code>}, {min:0,
max:2, value:1, step:0.01, label:<code>ratio 4</code>} ]}
setSliders={this.setRatio} onSlide={this.changeRatio}
context={this}/&gt;
</p>
<p>
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
<em>relative to all other points</em>.
</p>
<div className="howtocode">
<h3>How to implement rational curves</h3>
<p>
Extending the code of the previous section to include ratios is
almost trivial:
</p>
<pre><code>function RationalBezier(2,t,w[],r[]):
t2 = t * t
mt = 1-t
mt2 = mt * mt
f = [
r[0] * mt2,
2 * r[1] * mt * t,
r[2] * t2
]
basis = f[0] + f[1] + f[2]
return (f[0] * w[0] + f[1] * w[1] + f[2] * w[2])/basis
function RationalBezier(3,t,w[],r[]):
t2 = t * t
t3 = t2 * t
mt = 1-t
mt2 = mt * mt
mt3 = mt2 * mt
f = [
r[0] * mt3,
3 * r[1] * mt2 * t,
3 * r[2] * mt * t2,
r[3] * t3
]
basis = f[0] + f[1] + f[2] + f[3]
return (f[0] * w[0] + f[1] * w[1] + f[2] * w[2] + f[3] * w[3])/basis</code></pre>
<p>And that's all we have to do.</p>
</div>
</section>
<section id="extended">
<h1>The Bézier interval [0,1]</h1>
<p>
Now that we know the mathematics behind Bézier curves, there's one
curious thing that you may have noticed: they always run from
<code>t=0</code> to <code>t=1</code>. Why that particular interval?
</p>
<p>
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:
</p>
<img
class="LaTeX SVG"
src="images/latex/b80a1cac1f9ec476d6f6646ce0e154e7.svg"
width="215px"
height="16px"
loading="lazy"
/>
<p>
The obvious start and end values here need to be
<code>a=1, b=0</code>, so that the mixed value is 100% value 1, and
0% value 2, and <code>a=0, b=1</code>, 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 <strong>and</strong> 100% value 2. In principle that's fine,
but for Bézier curves we always want mixed values
<em>between</em> 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:
</p>
<img
class="LaTeX SVG"
src="images/latex/d930dea961b40f4810708bd6746221a2.svg"
width="215px"
height="16px"
loading="lazy"
/>
<p>
With this we can guarantee that we never sum above 100%. By
restricting <code>a</code> 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.
</p>
<p>
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".
</p>
<p>
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?
</p>
<p>
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 <code>t</code> 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.
</p>
<Graphic
title="Quadratic infinite interval Bézier curve"
setup="{this.setupQuadratic}"
draw="{this.draw}"
/>
<Graphic
title="Cubic infinite interval Bézier curve"
setup="{this.setupCubic}"
draw="{this.draw}"
/>
<p>
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
<a href="http://levien.com/phd/phd.html">"Spiro" curve</a>, which is
a curve based on part of a
<a href="https://en.wikipedia.org/wiki/Euler_spiral"
>Cornu Spiral, also known as Euler's Spiral</a
>. It's a very aesthetically pleasing curve and you'll find it in
quite a few graphics packages like
<a href="https://fontforge.github.io">FontForge</a> and
<a href="https://inkscape.org">Inkscape</a>. It has even been used
in font design, for example for the Inconsolata typeface.
</p>
</section>
<section id="matrix">
<h1>Bézier curvatures as matrix operations</h1>
<p>
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<sub>...</sub>
to refer to coordinate values "in one or more dimensions":
</p>
<img
class="LaTeX SVG"
src="images/latex/5aea6d4d5855135051715fb1cc0ec531.svg"
width="468px"
height="20px"
loading="lazy"
/>
<p>Disregarding our actual coordinates for a moment, we have:</p>
<img
class="LaTeX SVG"
src="images/latex/e524525c62234ce616a1e51c9848c169.svg"
width="353px"
height="19px"
loading="lazy"
/>
<p>We can write this as a sum of four expressions:</p>
<img
class="LaTeX SVG"
src="images/latex/c1f815481ad5132bebc1b1f0a3edf20f.svg"
width="140px"
height="75px"
loading="lazy"
/>
<p>And we can expand these expressions:</p>
<img
class="LaTeX SVG"
src="images/latex/9bc905d79bb22580b8c1cd75a791db73.svg"
width="397px"
height="75px"
loading="lazy"
/>
<p>Furthermore, we can make all the 1 and 0 factors explicit:</p>
<img
class="LaTeX SVG"
src="images/latex/e0d89b48cd11a726c00a2f689d48d57c.svg"
width="217px"
height="75px"
loading="lazy"
/>
<p>
And <em>that</em>, we can view as a series of four matrix
operations:
</p>
<img
class="LaTeX SVG"
src="images/latex/6da69918482a0b6b84d90a72dbeae9dd.svg"
width="607px"
height="72px"
loading="lazy"
/>
<p>If we compact this into a single matrix operation, we get:</p>
<img
class="LaTeX SVG"
src="images/latex/e94ae04eb5732c05d38fa1c97a2a25b0.svg"
width="227px"
height="72px"
loading="lazy"
/>
<p>
This kind of polynomial basis representation is generally written
with the bases in increasing order, which means we need to flip our
<code>t</code> matrix horizontally, and our big "mixing" matrix
upside down:
</p>
<img
class="LaTeX SVG"
src="images/latex/24bdad213879407a35b23c18394293aa.svg"
width="227px"
height="72px"
loading="lazy"
/>
<p>
And then finally, we can add in our original coordinates as a single
third matrix:
</p>
<img
class="LaTeX SVG"
src="images/latex/009c671bc526b5d75c30411c3c3a7e91.svg"
width="323px"
height="73px"
loading="lazy"
/>
<p>
We can perform the same trick for the quadratic curve, in which case
we end up with:
</p>
<img
class="LaTeX SVG"
src="images/latex/77a11d65d7cffc4b84a85c4bec837792.svg"
width="263px"
height="55px"
loading="lazy"
/>
<p>
If we plug in a <code>t</code> 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.
</p>
<p>
<strong>So: why would we bother with matrices?</strong> Matrix
representations allow us to discover things about functions that
would otherwise be hard to tell. It turns out that the curves form
<a href="https://en.wikipedia.org/wiki/Triangular_matrix"
>triangular matrices</a
>, 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
href="https://en.wikipedia.org/wiki/Invertible_matrix#The_invertible_matrix_theorem"
>a ton of properties</a
>
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
<em>immediately</em> 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.
</p>
<p>
So for now, just remember that we can represent curves this way, and
let's move on.
</p>
</section>
<section id="decasteljau">
<h1>de Casteljau's algorithm</h1>
<p>
If we want to draw Bézier curves, we can run through all values of
<code>t</code> from 0 to 1 and then compute the weighted basis
function at each value, getting the <code>x/y</code> values we need
to plot. Unfortunately, the more complex the curve gets, the more
expensive this computation becomes. Instead, we can use
<em>de Casteljau's algorithm</em> 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.
</p>
<p>
Rather than using our calculus function to find
<code>x/y</code> values for <code>t</code>, let's do this instead:
</p>
<ul>
<li>
treat <code>t</code> as a ratio (which it is). t=0 is 0% along a
line, t=1 is 100% along a line.
</li>
<li>
Take all lines between the curve's defining points. For an order
<code>n</code> curve, that's <code>n</code> lines.
</li>
<li>
Place markers along each of these line, at distance
<code>t</code>. So if <code>t</code> is 0.2, place the mark at 20%
from the start, 80% from the end.
</li>
<li>
Now form lines between <code>those</code> points. This gives
<code>n-1</code> lines.
</li>
<li>
Place markers along each of these line at distance <code>t</code>.
</li>
<li>
Form lines between <code>those</code> points. This'll be
<code>n-2</code> lines.
</li>
<li>Place markers, form lines, place markers, etc.</li>
<li>
Repeat this until you have only one line left. The point
<code>t</code> on that line coincides with the original curve
point at <code>t</code>.
</li>
</ul>
<div className="howtocode">
<h3>How to implement de Casteljau's algorithm</h3>
<p>
Let's just use the algorithm we just specified, and implement
that:
</p>
<pre><code>function drawCurve(points[], t):
if(points.length==1):
draw(points[0])
else:
newpoints=array(points.size-1)
for(i=0; i&lt;newpoints.length; i++):
newpoints[i] = (1-t) * points[i] + t * points[i+1]
drawCurve(newpoints, t)</code></pre>
<p>
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
<code>x</code> and <code>y</code> values:
</p>
<pre><code>function drawCurve(points[], t):
if(points.length==1):
draw(points[0])
else:
newpoints=array(points.size-1)
for(i=0; i&lt;newpoints.length; i++):
x = (1-t) * points[i].x + t * points[i+1].x
y = (1-t) * points[i].y + t * points[i+1].y
newpoints[i] = new point(x,y)
drawCurve(newpoints, t)</code></pre>
<p>
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 <i>t</i> ratios (i.e. the "markers"
outlined in the above algorithm), and then call the draw function
for this new list.
</p>
</div>
<p>
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.
</p>
<p>
&lt;Graphictitle="Traversing a curve using de Casteljau's algorithm"
setup={this.setup} draw={this.draw}/&gt;
</p>
</section>
<section id="flattening">
<h1>Simplified drawing</h1>
<p>
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.
</p>
<p>
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 for doing true intersection
detection, or curvature alignment.
</p>
<Graphic
title="Flattening a quadratic curve"
setup="{this.setupQuadratic}"
draw="{this.drawFlattened}"
onKeyDown="{this.onKeyDown}"
/>
<Graphic
title="Flattening a cubic curve"
setup="{this.setupCubic}"
draw="{this.drawFlattened}"
onKeyDown="{this.onKeyDown}"
/>
<p>
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.
</p>
<div className="howtocode">
<h3>How to implement curve flattening</h3>
<p>
Let's just use the algorithm we just specified, and implement
that:
</p>
<pre><code>function flattenCurve(curve, segmentCount):
step = 1/segmentCount;
coordinates = [curve.getXValue(0), curve.getYValue(0)]
for(i=1; i &lt;= segmentCount; i++):
t = i*step;
coordinates.push[curve.getXValue(t), curve.getYValue(t)]
return coordinates;</code></pre>
<p>
And done, that's the algorithm implemented. That just leaves
drawing the resulting "curve" as a sequence of lines:
</p>
<pre><code>function drawFlattenedCurve(curve, segmentCount):
coordinates = flattenCurve(curve, segmentCount)
coord = coordinates[0], _coord;
for(i=1; i &lt; coordinates.length; i++):
_coord = coordinates[i]
line(coord, _coord)
coord = _coord</code></pre>
<p>
We start with the first coordinate as reference point, and then
just draw lines between each point and its next point.
</p>
</div>
</section>
<section id="splitting">
<h1>Splitting curves</h1>
<p>
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 <code>t</code>, the procedure
gives us all the points we need to split a curve at that
<code>t</code> 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.
</p>
<Graphic
title="Splitting a curve"
setup="{this.setupCubic}"
draw="{this.drawSplit}"
/>
<div className="howtocode">
<h3>implementing curve splitting</h3>
<p>
We can implement curve splitting by bolting some extra logging
onto the de Casteljau function:
</p>
<pre><code>left=[]
right=[]
function drawCurve(points[], t):
if(points.length==1):
left.add(points[0])
right.add(points[0])
draw(points[0])
else:
newpoints=array(points.size-1)
for(i=0; i&lt;newpoints.length; i++):
if(i==0):
left.add(points[i])
if(i==newpoints.length-1):
right.add(points[i+1])
newpoints[i] = (1-t) * points[i] + t * points[i+1]
drawCurve(newpoints, t)</code></pre>
<p>
After running this function for some value <code>t</code>, the
<code>left</code> and <code>right</code> arrays will contain all
the coordinates for two new curves - one to the "left" of our
<code>t</code> 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.
</p>
</div>
<p>
This is best illustrated with an animated graphic (click to
play/pause):
</p>
<Graphic
title="Bézier curve splitting"
setup="{this.setupCubic}"
draw="{this.drawAnimated}"
onClick="{this.togglePlay}"
/>
</section>
<section id="matrixsplit">
<h1>Splitting curves using matrices</h1>
<p>
Another way to split curves is to exploit the matrix representation
of a Bézier curve. In <a href="#matrix">the section on matrices</a>,
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)
</p>
<img
class="LaTeX SVG"
src="images/latex/77a11d65d7cffc4b84a85c4bec837792.svg"
width="263px"
height="55px"
loading="lazy"
/>
<p>and</p>
<img
class="LaTeX SVG"
src="images/latex/c58330e12d25c678b593ddbd4afa7c52.svg"
width="323px"
height="73px"
loading="lazy"
/>
<p>
Let's say we want to split the curve at some point
<code>t = z</code>, 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:
</p>
<img
class="LaTeX SVG"
src="images/latex/278b67e9b908f4abcf2e9d069a6b29a4.svg"
width="648px"
height="55px"
loading="lazy"
/>
<p>and</p>
<img
class="LaTeX SVG"
src="images/latex/dbdbbe9aed4dacb1c1c5ae29b4371870.svg"
width="805px"
height="75px"
loading="lazy"
/>
<p>
If we could compact these matrices back to the form
<strong>[t values] · [Bézier matrix] · [column matrix]</strong>,
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 <code>t = 0</code> to <code>t = z</code>. 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!).
</p>
<div className="note">
<h2>Deriving new hull coordinates</h2>
<p>
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:
</p>
<img
class="LaTeX SVG"
src="images/latex/c79b607a92c42789fde57c6a8c4259fd.svg"
width="348px"
height="55px"
loading="lazy"
/>
<img
class="LaTeX SVG"
src="images/latex/8fb4faa046191480e89052102ecd3678.svg"
width="247px"
height="55px"
loading="lazy"
/>
<img
class="LaTeX SVG"
src="images/latex/4063d3462c179e91bb5f97c5e763560a.svg"
width="247px"
height="55px"
loading="lazy"
/>
<img
class="LaTeX SVG"
src="images/latex/a56f198daab08d20ef666599af14f513.svg"
width="252px"
height="55px"
loading="lazy"
/>
<p>
We can do this because [<em>M · M<sup>-1</sup></em
>] 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 [<em>M · M<sup>-1</sup></em
>] has no effect on the total formula, but it does allow us to
change the matrix sequence [<em>something · M</em>] to a sequence
[<em>M · something</em>], and that makes a world of difference: if
we know what [<em>M<sup>-1</sup> · Z · M</em>] is, we can apply
that to our coordinates, and be left with a proper matrix
representation of a quadratic Bézier curve (which is [<em
>T · M · P</em
>]), with a new set of coordinates that represent the curve from
<em>t = 0</em> to <em>t = z</em>. So let's get computing:
</p>
<img
class="LaTeX SVG"
src="images/latex/7d629178a5fb985a35770002d1912535.svg"
width="627px"
height="56px"
loading="lazy"
/>
<p>Excellent! Now we can form our new quadratic curve:</p>
<img
class="LaTeX SVG"
src="images/latex/b5cf45e4b34fdd18f599b79549844d45.svg"
width="417px"
height="55px"
loading="lazy"
/>
<img
class="LaTeX SVG"
src="images/latex/859b7bc7b78e8e297ae5fddd9be40ab7.svg"
width="479px"
height="55px"
loading="lazy"
/>
<img
class="LaTeX SVG"
src="images/latex/abb3edce2229312f351d81092ba2145b.svg"
width="492px"
height="55px"
loading="lazy"
/>
<p>
<strong><em>Brilliant</em></strong
>: if we want a subcurve from <code>t = 0</code> to
<code>t = z</code>, 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
<a href="https://en.wikipedia.org/wiki/Bernstein_polynomial"
>Bernstein polynomial</a
>
of degree two. These new coordinates are actually really easy to
compute directly!
</p>
<p>
Of course, that's only one of the two curves. Getting the section
from <code>t = z</code> to <code>t = 1</code> requires doing this
again. We first observe that in the previous calculation, we
actually evaluated the general interval [0,<code>z</code>]. We
were able to write it down in a more simple form because of the
zero, but what we <em>actually</em> evaluated, making the zero
explicit, was:
</p>
<img
class="LaTeX SVG"
src="images/latex/2f2bec1e77039a40c31220f5bf83641a.svg"
width="381px"
height="55px"
loading="lazy"
/>
<img
class="LaTeX SVG"
src="images/latex/daaae36f13bb97f2a7ac21eec6903755.svg"
width="313px"
height="55px"
loading="lazy"
/>
<p>
If we want the interval [<em>z</em>,1], we will be evaluating this
instead:
</p>
<img
class="LaTeX SVG"
src="images/latex/0f84dbf6e3ea7db732ceb9d71caf9b22.svg"
width="461px"
height="55px"
loading="lazy"
/>
<img
class="LaTeX SVG"
src="images/latex/a34473afe7a4160b45ce0f2a770fad99.svg"
width="412px"
height="57px"
loading="lazy"
/>
<p>
We're going to do the same trick of multiplying by the identity
matrix, to turn <code>[something · M]</code> into
<code>[M · something]</code>:
</p>
<img
class="LaTeX SVG"
src="images/latex/567c29ee78b49c700f54b17780682543.svg"
width="729px"
height="57px"
loading="lazy"
/>
<p>So, our final second curve looks like:</p>
<img
class="LaTeX SVG"
src="images/latex/19049f556723a4f2d985a631a91ae290.svg"
width="421px"
height="55px"
loading="lazy"
/>
<img
class="LaTeX SVG"
src="images/latex/e9f64464287d3d5c6a4cbe64e21746c8.svg"
width="473px"
height="57px"
loading="lazy"
/>
<img
class="LaTeX SVG"
src="images/latex/f2695b6d6417c60343b4934dae8118f8.svg"
width="492px"
height="57px"
loading="lazy"
/>
<p>
<strong><em>Nice</em></strong
>. 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 <em>also</em> really
easy to compute directly!
</p>
</div>
<p>
So, using linear algebra rather than de Casteljau's algorithm, we
have determined that, for any quadratic curve split at some value
<code>t = z</code>, we get two subcurves that are described as
Bézier curves with simple-to-derive coordinates:
</p>
<img
class="LaTeX SVG"
src="images/latex/e16eba6dfb9f0b8d1abc3e1cd3ba63a2.svg"
width="576px"
height="55px"
loading="lazy"
/>
<p>and</p>
<img
class="LaTeX SVG"
src="images/latex/6aeb749eb26f5a9199c1b16d7d421dc0.svg"
width="571px"
height="57px"
loading="lazy"
/>
<p>
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:
</p>
<img
class="LaTeX SVG"
src="images/latex/5e3fae45d325d0f0681731fb606b6fbc.svg"
width="841px"
height="75px"
loading="lazy"
/>
<p>and</p>
<img
class="LaTeX SVG"
src="images/latex/0d2e895e767c4cecb0fccafee1273152.svg"
width="837px"
height="77px"
loading="lazy"
/>
<p>
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 <strong><em>Q</em></strong> 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" <strong><em>Q'</em></strong
>.
</p>
<p>
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.
</p>
</section>
<section id="reordering">
<h1>Lowering and elevating curve order</h1>
<p>
One interesting property of Bézier curves is that an
<em>n<sup>th</sup></em> order curve can always be perfectly
represented by an <em>(n+1)<sup>th</sup></em> order curve, by giving
the higher-order curve specific control points.
</p>
<p>
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/3<sup>rd</sup> start + 2/3<sup>rd</sup> control"
and "2/3<sup>rd</sup> control + 1/3<sup>rd</sup> end". Now we have
exactly the same curve as before, except represented as a cubic
curve rather than a quadratic curve.
</p>
<p>
The general rule for raising an <em>n<sup>th</sup></em> order curve
to an <em>(n+1)<sup>th</sup></em> 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):
</p>
<img
class="LaTeX SVG"
src="images/latex/faf29599c9307f930ec28065c96fde2a.svg"
width="768px"
height="61px"
loading="lazy"
/>
<p>
However, this rule also has as direct consequence that you
<strong>cannot</strong> generally safely lower a curve from
<em>n<sup>th</sup></em> order to <em>(n-1)<sup>th</sup></em> 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.
</p>
<p>
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
<a
href="http://www.sirver.net/blog/2011/08/23/degree-reduction-of-bezier-curves"
>Sirver's Castle</a
>). 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!
</p>
<p>
We start by taking the standard Bézier function, and condensing it a
little:
</p>
<img
class="LaTeX SVG"
src="images/latex/1244a85c1f9044b6f77cb709c682159c.svg"
width="408px"
height="41px"
loading="lazy"
/>
<p>
Then, we apply one of those silly (actually, super useful) calculus
tricks: since our <code>t</code> value is always between zero and
one (inclusive), we know that <code>(1-t)</code> plus
<code>t</code> always sums to 1. As such, we can express any value
as a sum of <code>t</code> and <code>1-t</code>:
</p>
<img
class="LaTeX SVG"
src="images/latex/b2fda1dcce5bb13317aa42ebf5e7ea6c.svg"
width="379px"
height="16px"
loading="lazy"
/>
<p>
So, with that seemingly trivial observation, we rewrite that Bézier
function by splitting it up into a sum of a <code>(1-t)</code> and
<code>t</code> component:
</p>
<img
class="LaTeX SVG"
src="images/latex/41e184228d85023abdadd6ce2acb54c7.svg"
width="316px"
height="67px"
loading="lazy"
/>
<p>
So far so good. Now, to see why we did this, let's write out the
<code>(1-t)</code> and <code>t</code> parts, and see what that gives
us. I promise, it's about to make sense. We start with
<code>(1-t)</code>:
</p>
<img
class="LaTeX SVG"
src="images/latex/4debbed5922d2bd84fd322c616872d20.svg"
width="387px"
height="160px"
loading="lazy"
/>
<p>
So by using this seemingly silly trick, we can suddenly express part
of our n<sup>th</sup> order Bézier function in terms of an (n+1)<sup
>th</sup
>
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
<code>t</code> part, but that's not a problem:
</p>
<img
class="LaTeX SVG"
src="images/latex/483c89c8726f7fd0dca0b7de339b04bd.svg"
width="471px"
height="159px"
loading="lazy"
/>
<p>
So, with both of those changed from an order
<code>n</code> expression to an order <code>(n+1)</code> expression,
we can put them back together again. Now, where the order
<code>n</code> function had a summation from 0 to <code>n</code>,
the order <code>n+1</code> function uses a summation from 0 to
<code>n+1</code>, 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.
</p>
<p>Let's do this:</p>
<img
class="LaTeX SVG"
src="images/latex/dd8d8d98f66ce9f51b95cbf48225e97b.svg"
width="465px"
height="257px"
loading="lazy"
/>
<p>
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:
</p>
<img
class="LaTeX SVG"
src="images/latex/773fdc86b686647c823b4f499aca3a35.svg"
width="71px"
height="16px"
loading="lazy"
/>
<p>
where the matrix <strong>M</strong> is an <code>n+1</code> by
<code>n</code> matrix, and looks like:
</p>
<img
class="LaTeX SVG"
src="images/latex/7a9120997e4a4855ecda435553a7bbdf.svg"
width="336px"
height="187px"
loading="lazy"
/>
<p>
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.
</p>
<p>Not too bad!</p>
<p>
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
<a href="http://mathworld.wolfram.com/NormalEquation.html"
>the normal equation</a
>. 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 <strong>A x = b</strong>, we
can use it. And as it so happens, that's exactly what we're dealing
with, so:
</p>
<img
class="LaTeX SVG"
src="images/latex/d52f60b331c1b8d6733eb5217adfbc4d.svg"
width="272px"
height="116px"
loading="lazy"
/>
<p>The steps taken here are:</p>
<ol>
<li>
We have a function in a form that the normal equation can be used
with, so
</li>
<li>apply the normal equation!</li>
<li>
Then, we want to end up with just B<sub>n</sub> 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
<a href="https://en.wikipedia.org/wiki/Identity_matrix"
>identity matrix</a
>.
</li>
<li>
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
<strong>I</strong>. So we substitute the mess with
<strong>I</strong>, and then
</li>
<li>
because multiplication with the identity matrix does nothing (like
multiplying by 1 does nothing in regular algebra), we just drop
it.
</li>
</ol>
<p>
And we're done: we now have an expression that lets us approximate
an <code>n+1</code><sup>th</sup> order curve with a lower
<code>n</code><sup>th</sup> 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.
</p>
<p>
&lt;Graphic title={"A " + this.getOrder() + " order Bézier curve"}
setup={this.setup} draw={this.draw} onKeyDown={this.props.onKeyDown}
onMouseMove={this.onMouseMove} /&gt;
</p>
</section>
<section id="derivatives">
<h1>Derivatives</h1>
<p>
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.
</p>
<p>
First, let's look at the derivative rule for Bézier curves, which
is:
</p>
<img
class="LaTeX SVG"
src="images/latex/eb4442acc5bc17f4649eb04b2953ed9b.svg"
width="333px"
height="44px"
loading="lazy"
/>
<p>
which we can also write (observing that <i>b</i> in this formula is
the same as our <i>w</i> weights, and that <i>n</i> times a
summation is the same as a summation where each term is multiplied
by <i>n</i>) as:
</p>
<img
class="LaTeX SVG"
src="images/latex/2622790efa97f1915e7998787d8ce977.svg"
width="343px"
height="44px"
loading="lazy"
/>
<p>
Or, in plain text: the derivative of an n<sup>th</sup> degree Bézier
curve is an (n-1)<sup>th</sup> degree Bézier curve, with one fewer
term, and new weights w'<sub>0</sub>...w'<sub>n-1</sub> derived from
the original weights as n(w<sub>i+1</sub> - w<sub>i</sub>). So for a
3<sup>rd</sup> degree curve, with four weights, the derivative has
three new weights: w'<sub>0</sub> = 3(w<sub>1</sub>-w<sub>0</sub>),
w'<sub>1</sub> = 3(w<sub>2</sub>-w<sub>1</sub>) and w'<sub>2</sub> =
3(w<sub>3</sub>-w<sub>2</sub>).
</p>
<div className="note">
<h3>"Slow down, why is that true?"</h3>
<p>
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:
</p>
<img
class="LaTeX SVG"
src="images/latex/e755c2adfec5d266c50e064407ca369b.svg"
width="209px"
height="36px"
loading="lazy"
/>
<p>
Applying the
<a href="http://en.wikipedia.org/wiki/Product_rule">product</a>
and
<a href="http://en.wikipedia.org/wiki/Chain_rule">chain</a> rules
gives us:
</p>
<img
class="LaTeX SVG"
src="images/latex/95a0cd4cc919a3fd5b192ffeb00c231e.svg"
width="412px"
height="28px"
loading="lazy"
/>
<p>Which is hard to work with, so let's expand that properly:</p>
<img
class="LaTeX SVG"
src="images/latex/bd3c740be364071c86ccf42b99d5eba4.svg"
width="344px"
height="27px"
loading="lazy"
/>
<p>
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 <i>n-1</i> and <i>k-1</i>, we'll be on the right
track.
</p>
<img
class="LaTeX SVG"
src="images/latex/6770214cceeb0e13e371bd908867751f.svg"
width="545px"
height="76px"
loading="lazy"
/>
<p>
And that's the first part done: the two components inside the
parentheses are actually regular, lower-order Bézier expressions:
</p>
<img
class="LaTeX SVG"
src="images/latex/fb823558e99662b24d46ae55ac93ce38.svg"
width="533px"
height="48px"
loading="lazy"
/>
<p>
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:
</p>
<img
class="LaTeX SVG"
src="images/latex/28991bba7c13698619f36b6261d91d68.svg"
width="527px"
height="112px"
loading="lazy"
/>
<p>
If we expand this (with some color to show how terms line up), and
reorder the terms by increasing values for <i>k</i> we see the
following:
</p>
<img
class="LaTeX SVG"
src="images/latex/183ac1365b1075166b7066b78a17ad74.svg"
width="447px"
height="109px"
loading="lazy"
/>
<p>
Two of these terms fall way: the first term falls away because
there is no -1<sup>st</sup> 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
<i>B<sub>n-1,n</sub></i
>. This term would have a binomial coefficient of [<i>i</i> choose
<i>i+1</i>], 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:
</p>
<img
class="LaTeX SVG"
src="images/latex/950811d35cebee2a39d881e33cbfb1de.svg"
width="447px"
height="71px"
loading="lazy"
/>
<p>And that's just a summation of lower order curves:</p>
<img
class="LaTeX SVG"
src="images/latex/894ccddbfcec09dd3d035777e7d23245.svg"
width="807px"
height="36px"
loading="lazy"
/>
<p>We can rewrite this as a normal summation, and we're done:</p>
<img
class="LaTeX SVG"
src="images/latex/514090a0fd6c64b7d85a9dc5721a0fa6.svg"
width="545px"
height="51px"
loading="lazy"
/>
</div>
<p>
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:
</p>
<img
class="LaTeX SVG"
src="images/latex/14cb9fbbaae9e7d87ae6bef3ea7a782e.svg"
width="352px"
height="55px"
loading="lazy"
/>
<img
class="LaTeX SVG"
src="images/latex/c010c0df4bb911b84da6e9d379617e4b.svg"
width="533px"
height="59px"
loading="lazy"
/>
<p>
What are the differences? In terms of the actual Bézier curve,
virtually nothing! We lowered the order (rather than <i>n</i>, it's
now <i>n-1</i>), 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:
</p>
<img
class="LaTeX SVG"
src="images/latex/03967e3ecdbff78684995ca9c22a6106.svg"
width="523px"
height="73px"
loading="lazy"
/>
<p>
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
<i>k = 0</i>, 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 <i>n<sup>th</sup></i> order curve
has <i>n-1</i> (meaningful) derivatives, with any further derivative
being zero.
</p>
</section>
<section id="pointvectors">
<h1>Tangents and normals</h1>
<p>
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:
</p>
<img
class="LaTeX SVG"
src="images/latex/d236b7b2ad46c8ced1b43bb2a496379a.svg"
width="141px"
height="40px"
loading="lazy"
/>
<p>
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:
</p>
<img
class="LaTeX SVG"
src="images/latex/6101b2f8b69ebabba4a2c88456a32aa0.svg"
width="251px"
height="25px"
loading="lazy"
/>
<img
class="LaTeX SVG"
src="images/latex/009715fce01e46e7c07f87a8192a8c62.svg"
width="289px"
height="57px"
loading="lazy"
/>
<p>
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 <i>t</i>? In that case we want the
<em>normal</em> 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:
</p>
<img
class="LaTeX SVG"
src="images/latex/2dd2f89d1c762991a86526490a3deef6.svg"
width="339px"
height="60px"
loading="lazy"
/>
<div className="note">
<p>
Rotating coordinates is actually very easy, if you know the rule
for it. You might find it explained as "applying a
<a href="https://en.wikipedia.org/wiki/Rotation_matrix"
>rotation matrix</a
>, 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 cirle by a quarter turn, and done.
</p>
<p>
To turn any point <i>(x,y)</i> into a rotated point
<i>(x',y')</i> (over 0,0) by some angle φ, we apply this nice and
easy computation:
</p>
<img
class="LaTeX SVG"
src="images/latex/deec095950fcd1f9c980be76a7093fe6.svg"
width="183px"
height="37px"
loading="lazy"
/>
<p>
Which is the "long" version of the following matrix
transformation:
</p>
<img
class="LaTeX SVG"
src="images/latex/2a55cb2d23c25408aa10cfd8db13278b.svg"
width="211px"
height="40px"
loading="lazy"
/>
<p>
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 <em>sin</em> and <em>cos</em> for these angles are,
respectively: 0 and 1, -1 and 0, and 0 and -1.
</p>
<p>
But <strong><em>why</em></strong> does this work? Why this matrix
multiplication?
<a
href="http://en.wikipedia.org/wiki/Rotation_matrix#Decomposition_into_shears"
>Wikipedia</a
>
(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.
<a href="http://datagenetics.com/blog/august32013/index.html"
>DataGenetics</a
>
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.
</p>
</div>
<p>
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 <em>t</em>-intervals, not spaced equidistant).
</p>
<div className="figure">
<Graphic
title="Quadratic Bézier tangents and normals"
inline="{true}"
setup="{this.setupQuadratic}"
draw="{this.draw}"
/>
<Graphic
title="Cubic Bézier tangents and normals"
inline="{true}"
setup="{this.setupCubic}"
draw="{this.draw}"
/>
</div>
</section>
<section id="pointvectors3d">
<h1>Working with 3D normals</h1>
<p>
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.
</p>
<p>
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.
</p>
<p>
The "naïve" approach is to construct what is known as the
<a
href="https://en.wikipedia.org/wiki/Frenet%E2%80%93Serret_formulas"
>Frenet normal</a
>, 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".
</p>
<p>
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".
</p>
<p>
So let's do that! And in a twist surprise, we can do this in four
lines:
</p>
<ul>
<li><strong>a</strong> = normalize(B'(t))</li>
<li><strong>b</strong> = normalize(<strong>a</strong> + B''(t))</li>
<li>
<strong>r</strong> = normalize(<strong>b</strong> ×
<strong>a</strong>)
</li>
<li>
<strong>normal</strong> = normalize(<strong>r</strong> ×
<strong>a</strong>)
</li>
</ul>
<p>Let's unpack that a little:</p>
<ul>
<li>
We start by taking the
<a href="https://en.wikipedia.org/wiki/Unit_vector"
>normalized vector</a
>
for the derivative at some point on the curve. We normalize it so
the maths is less work. Less work is good.
</li>
<li>
Then, we compute <strong>b</strong> 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.
</li>
<li>
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
<a href="https://en.wikipedia.org/wiki/Cross_product"
>cross product</a
>. (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.
</li>
<li>
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.
</li>
</ul>
<p>
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:
</p>
<Graphic
title="Some known and unknown vectors"
setup="{this.setup}"
draw="{this.drawFrenetVectors}"
/>
<p>
However, if you've played with that graphic a bit, you might have
noticed something odd. The normal seems to "suddenly twist around"
around between t=0.5 and t=0.75 - why is doing that?
</p>
<p>
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.
</p>
<p>Thankfully, Frenet normals are not our only option.</p>
<p>
Another option is to take a slightly more algorithmic approach and
compute a form of
<a
href="https://www.microsoft.com/en-us/research/wp-content/uploads/2016/12/Computation-of-rotation-minimizing-frames.pdf"
>Rotation Minimising Frame</a
>
(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.
</p>
<p>
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.
</p>
<p>
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:
</p>
<ul>
<li>
Take a point on the curve for which we know the RM frame already,
</li>
<li>
take a next point on the curve for which we don't know the RM
frame yet, and
</li>
<li>
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".
</li>
<li>
This gives the next point a tangent vector that essentially in the
opposite direction of what it should be in, and a normal that's
slightly off-kilter, so:
</li>
<li>
reflect the vectors of our "mirrored frame" a second time, but
this time using the plane through the next point itself as
"mirror".
</li>
<li>
Done: the tangent and normal have been fixed, and we have a good
looking frame to work with.
</li>
</ul>
<p>So, let's write some code for that!</p>
<div className="howtocode">
<h3>Implementing Rotation Minimising Frames</h3>
<p>
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:
</p>
<pre><code>{
o: origin of all vectors, i.e. the on-curve point,
t: tangent vector,
r: rotational axis vector,
n: normal vector
}</code></pre>
<p>
Then, we can write a function that generates a sequence of RM
frames in the following manner:
</p>
<pre><code>generateRMFrames(steps) -&gt; frames:
step = 1.0/steps
// Start off with the standard tangent/axis/normal frame
// associated with the curve at t=0:
frames.add(getFrenetFrame(0))
// start constructing RM frames:
for t0 = 0, t0 &lt; 1.0, t0 += step:
// start with the previous, known frame
x0 = frames.last
// get the next frame: we're going to keep its position and tangent,
// but we're going to recompute the axis and normal.
t1 = t0 + step
x1 = { o: getPoint(t1), t: getDerivative(t) }
// First we reflect x0's tangent and axis of rotation onto x1,
// through/ the plane of reflection at the point between x0 x1
v1 = x1.o - x0.o
c1 = v1 · v1
riL = x0.r - v1 * 2/c1 * v1 · x0.r
tiL = x0.t - v1 * 2/c1 * v1 · x0.t
// note that v1 is a vector, but 2/c1 and (v1 · ...) are just
// plain numbers, so we're just scaling v1 by some constant.
// Then we reflect a second time, over a plane at x1, so that
// the frame tangent is aligned with the curve tangent again:
v2 = x1.t - tiL
c2 = v2 · v2
// and we're done here:
x1.r = riL - v2 * 2/c2 * v2 · riL
x1.n = x1.r × x1.t
frames.add(x1)</code></pre>
<p>
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.
</p>
</div>
<p>
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:
</p>
<Graphic
title="Æsthetically much better 3D curve normals"
setup="{this.setup}"
draw="{this.drawRMFNormals}"
/>
<p>Now that looks much better!</p>
<p>
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 <strong>(0,0,1) × tangent</strong>, 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 quite useful.
</p>
</section>
<section id="components">
<h1>Component functions</h1>
<p>
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.
</p>
<p>
The solution: compute the derivative for each axis separately, and
then fit them back together in the same way we do for the original.
</p>
<p>
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 <i>t</i> (between 0 and 1
inclusive), and the y-axis value, respectively.
</p>
<p>
If you move points in a curve sideways, you should only see the
middle graph change; likely, moving points vertically should only
show a change in the right graph.
</p>
<Graphic
title="Quadratic Bézier curve components"
setup="{this.setupQuadratic}"
draw="{this.draw}"
/>
<Graphic
title="Cubic Bézier curve components"
setup="{this.setupCubic}"
draw="{this.draw}"
/>
</section>
<section id="extremities">
<h1>Finding extremities: root finding</h1>
<p>
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
equations B'(t) = 0 and B''(t) = 0. That said, in the case of
quadratic curves there is no B''(t), so we only need to compute
B'(t) = 0. So, how do we compute the first and second derivatives?
Fairly easily, actually, until our derivatives are 4th order or
higher... then things get really hard. But let's start simple:
</p>
<h3>Quadratic curves: linear derivatives.</h3>
<p>
Finding the solution for "where is this line 0" should be trivial:
</p>
<img
class="LaTeX SVG"
src="images/latex/ec28870926b6ed08625d942b578e6fbe.svg"
width="132px"
height="104px"
loading="lazy"
/>
<p>
Done. And quadratic curves have no meaningful second derivative, so
we're <em>really</em> done.
</p>
<h3>Cubic curves: the quadratic formula.</h3>
<p>
The derivative of a cubic curve is a quadratic curve, and finding
the roots for a quadratic Bézier curve means we can apply the
<a href="https://en.wikipedia.org/wiki/Quadratic_formula"
>Quadratic formula</a
>. If you've seen it before, you'll remember it, and if you haven't,
it looks like this:
</p>
<img
class="LaTeX SVG"
src="images/latex/0ec5cc72a428d75defb480530b50d720.svg"
width="411px"
height="40px"
loading="lazy"
/>
<p>
So, if we can express a 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?
</p>
<p>
First we turn our cubic Bézier function into a quadratic one, by
following the rule mentioned at the end of the
<a href="#derivatives">derivatives section</a>:
</p>
<img
class="LaTeX SVG"
src="images/latex/e06ec558d99b53e559d24524f4201951.svg"
width="537px"
height="36px"
loading="lazy"
/>
<p>
And then, using these <em>v</em> values, we can find out what our
<em>a</em>, <em>b</em>, and <em>c</em> should be:
</p>
<img
class="LaTeX SVG"
src="images/latex/ddc6f99a543afad25c55cf16b9deeed9.svg"
width="315px"
height="119px"
loading="lazy"
/>
<p>
This gives us thee coefficients <em>a</em>, <em>b</em>, and
<em>c</em> that are expressed in terms of <em>v</em> values, where
the <em>v</em> values are just convenient expressions of our
original <em>p</em> values, so we can do some trivial substitution
to get:
</p>
<img
class="LaTeX SVG"
src="images/latex/d9e66caeb45b6643112ce3d971b17e5b.svg"
width="308px"
height="63px"
loading="lazy"
/>
<p>
Easy-peasy. We can now almost trivially find the roots by plugging
those values into the quadratic formula. We also note that the
second derivative of a cubic curve means computing the first
derivative of a quadratic curve, and we just saw how to do that in
the section above.
</p>
<h3>Quartic curves: Cardano's algorithm.</h3>
<p>
Quartic—fourth degree—curves have a cubic function as derivative.
Now, cubic functions are a bit of a problem because they're really
hard to solve. But, way back in the 16<sup>th</sup> century,
<a href="https://en.wikipedia.org/wiki/Gerolamo_Cardano"
>Gerolamo Cardano</a
>
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 "easy", and then the only hard part is figuring out how to go
from that form to the generic form. So:
</p>
<img
class="LaTeX SVG"
src="images/latex/997a8cc704c0ab0e364cb8b532df90b0.svg"
width="253px"
height="44px"
loading="lazy"
/>
<p>
This is easier because for the "easier formula" we can use
<a href="http://www.wolframalpha.com/input/?i=t%5E3+%2B+pt+%2B+q"
>regular calculus</a
>
to find the roots. (As a cubic function, however, it can have up to
three roots, but two of those can be complex. For the purpose of
Bézier curve extremities, we can completely ignore those complex
roots, since our <em>t</em> is a plain real number from 0 to 1.)
</p>
<p>
So, the trick is to figure out how to turn the first formula into
the second formula, and to then work out the maths that gives us the
roots. This is explained in detail over at
<a
href="https://trans4mind.com/personal_development/mathematics/polynomials/cubicAlgebra.htm"
>Ken J. Ward's page</a
>
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.
</p>
<div className="howtocode">
<h3>Implementing Cardano's algorithm for finding all real roots</h3>
<p>
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
<a href="https://en.wikipedia.org/wiki/Complex_number"
>"complex" number</a
>
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!
</p>
<pre><code>// A helper function to filter for values in the [0,1] interval:
function accept(t) {
return 0&lt;=t && t &lt;=1;
}
// A real-cuberoots-only function:
function cuberoot(v) {
if(v&lt;0) return -pow(-v,1/3);
return pow(v,1/3);
}
// Now then: given cubic coordinates {pa, pb, pc, pd} find all roots.
function getCubicRoots(pa, pb, pc, pd) {
var a = (3*pa - 6*pb + 3*pc),
b = (-3*pa + 3*pb),
c = pa,
d = (-pa + 3*pb - 3*pc + pd);
// do a check to see whether we even need cubic solving:
if (approximately(d,0)) {
// this is not a cubic curve.
if (approximately(a,0)) {
// in fact, this is not a quadratic curve either.
if (approximately(b,0)) {
// in fact in fact, there are no solutions.
return [];
}
// linear solution
return [-c / b].filter(accept);
}
// quadratic solution
var q = sqrt(b*b - 4*a*c), 2a = 2*a;
return [(q-b)/2a, (-b-q)/2a].filter(accept)
}
// at this point, we know we need a cubic solution.
a /= d;
b /= d;
c /= d;
var p = (3*b - a*a)/3,
p3 = p/3,
q = (2*a*a*a - 9*a*b + 27*c)/27,
q2 = q/2,
discriminant = q2*q2 + p3*p3*p3;
// and some variables we're going to use later on:
var u1, v1, root1, root2, root3;
// three possible real roots:
if (discriminant &lt; 0) {
var mp3 = -p/3,
mp33 = mp3*mp3*mp3,
r = sqrt( mp33 ),
t = -q / (2*r),
cosphi = t&lt;-1 ? -1 : t&gt;1 ? 1 : t,
phi = acos(cosphi),
crtr = cuberoot(r),
t1 = 2*crtr;
root1 = t1 * cos(phi/3) - a/3;
root2 = t1 * cos((phi+2*pi)/3) - a/3;
root3 = t1 * cos((phi+4*pi)/3) - a/3;
return [root1, root2, root3].filter(accept);
}
// three real roots, but two of them are equal:
if(discriminant === 0) {
u1 = q2 &lt; 0 ? cuberoot(-q2) : -cuberoot(q2);
root1 = 2*u1 - a/3;
root2 = -u1 - a/3;
return [root1, root2].filter(accept);
}
// one real root, two complex roots
var sd = sqrt(discriminant);
u1 = cuberoot(sd - q2);
v1 = cuberoot(sd + q2);
root1 = u1 - v1 - a/3;
return [root1].filter(accept);
}</code></pre>
</div>
<p>
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
reduce 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.
</p>
<h3>Quintic and higher order curves: finding numerical solutions</h3>
<p>
The problem with this is that as the order of the curve goes up, we
can't actually solve those equations the normal way. We can't take
the function, and then work out what the solutions are. Not to
mention that even solving a third order derivative (for a fourth
order curve) is already a royal pain in the backside. We need a
better solution. We need numerical approaches.
</p>
<p>
That's a fancy word for saying "rather than solve the function,
treat the problem as a sequence of identical operations, the
performing of which gets us closer and closer to the real answer".
As it turns out, there is a really nice numerical root-finding
algorithm, called the
<a href="http://en.wikipedia.org/wiki/Newton-Raphson"
>Newton-Raphson</a
>
root finding method (yes, after
<em
><a href="https://en.wikipedia.org/wiki/Isaac_Newton">that</a></em
>
Newton), which we can make use of.
</p>
<p>
The Newton-Raphson approach consists of picking a value
<em>t</em> (any value will do), and getting the corresponding value
of the function at that <em>t</em> value. For normal functions, we
can treat that value as a height. If the height is zero, we're done,
we have found a root. If it's not, we take the tangent of the curve
at that point, and extend it until it passes the x-axis, which will
be at some new point <em>t</em>. We then repeat the procedure with
this new value, and we keep doing this until we find our root.
</p>
<p>
Mathematically, this means that for some <em>t</em>, at step
<em>n=1</em>, we perform the following calculation until
<em>f<sub>y</sub></em
>(<em>t</em>) is zero, so that the next <em>t</em> is the same as
the one we already have:
</p>
<img
class="LaTeX SVG"
src="images/latex/04336edba7697ca91861821e32fc14be.svg"
width="124px"
height="45px"
loading="lazy"
/>
<p>
(The Wikipedia article has a decent animation for this process, so
I'm not adding a sketch for that here unless there are requests for
it)
</p>
<p>
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?
</p>
<p>
As it turns out, Newton-Raphson is so blindingly fast, so we could
get away with just not picking: we simply run the algorithm from
<em>t=0</em> to <em>t=1</em> at small steps (say,
1/200<sup>th</sup>) 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 200<sup>th</sup> order Bézier curve is going to go
wrong, but that's okay: there is no reason, ever, to 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, that's pretty much as high as
you'll ever need to go.
</p>
<h3>In conclusion:</h3>
<p>
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:
</p>
<Graphic
title="Quadratic Bézier curve extremities"
setup="{this.setupQuadratic}"
draw="{this.draw}"
/>
<Graphic
title="Cubic Bézier curve extremities"
setup="{this.setupCubic}"
draw="{this.draw}"
/>
</section>
<section id="boundingbox">
<h1>Bounding boxes</h1>
<p>
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:
</p>
<p><em>Computing the bounding box for a Bézier curve</em>:</p>
<ol>
<li>
Find all <em>t</em> value(s) for the curve derivative's x- and
y-roots.
</li>
<li>
Discard any <em>t</em> value that's lower than 0 or higher than 1,
because Bézier curves only use the interval [0,1].
</li>
<li>
Determine the lowest and highest value when plugging the values
<em>t=0</em>, <em>t=1</em> 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.
</li>
</ol>
<p>
Applying this approach to our previous root finding, we get the
following bounding boxes (with all curve extremity points shown on
the curve):
</p>
<Graphic
title="Quadratic Bézier bounding box"
setup="{this.setupQuadratic}"
draw="{this.draw}"
/>
<Graphic
title="Cubic Bézier bounding box"
setup="{this.setupCubic}"
draw="{this.draw}"
/>
<p>
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.
</p>
</section>
<section id="aligning">
<h1>Aligning curves</h1>
<p>
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.
</p>
<p>
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 <em>n</em> term polynomial functions into <em>n-1</em> 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 <em>n-2</em> term function.
For instance, if we have a cubic curve such as this:
</p>
<img
class="LaTeX SVG"
src="images/latex/eb879c2e965ff4d7c61b70c0c6e5bf87.svg"
width="671px"
height="40px"
loading="lazy"
/>
<p>
Then translating it so that the first coordinate lies on (0,0),
moving all <em>x</em> coordinates by -120, and all
<em>y</em> coordinates by -160, gives us:
</p>
<img
class="LaTeX SVG"
src="images/latex/d6438a6372e9f910b0a3059ebfacccc6.svg"
width="655px"
height="40px"
loading="lazy"
/>
<p>
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:
</p>
<img
class="LaTeX SVG"
src="images/latex/6afb072a451e06ddf4b1ff1eb8948e8a.svg"
width="647px"
height="40px"
loading="lazy"
/>
<p>If we drop all the zero-terms, this gives us:</p>
<img
class="LaTeX SVG"
src="images/latex/1d44794e3aac3ecb881b559eb8a71fab.svg"
width="528px"
height="40px"
loading="lazy"
/>
<p>
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:
</p>
<Graphic
title="Aligning a quadratic curve"
setup="{this.setupQuadratic}"
draw="{this.draw}"
/>
<Graphic
title="Aligning a cubic curve"
setup="{this.setupCubic}"
draw="{this.draw}"
/>
</section>
<section id="tightbounds">
<h1>Tight boxes</h1>
<p>
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:
</p>
<Graphic
title="Aligning a quadratic curve"
setup="{this.setupQuadratic}"
draw="{this.draw}"
/>
<Graphic
title="Aligning a cubic curve"
setup="{this.setupCubic}"
draw="{this.draw}"
/>
<p>
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. If there is high demand for it,
I'll add a section on how to precisely compute the best fit bounding
box, but the maths is fairly grueling and just not really worth
spending time on.
</p>
</section>
<section id="inflections">
<h1>Curve inflections</h1>
<p>
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.
</p>
<p>
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.
</p>
<p>What we need to do is solve a simple equation:</p>
<img
class="LaTeX SVG"
src="images/latex/bafdb6583323bda71d9a15c02d1fdec2.svg"
width="59px"
height="16px"
loading="lazy"
/>
<p>
What we're saying here is that given the curvature function
<em>C(t)</em>, we want to know for which values of <em>t</em> 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 <em>C(t)</em> look like? Actually something that seems not too
hard:
</p>
<img
class="LaTeX SVG"
src="images/latex/2029bca9f4fa15739553636af99b70a8.svg"
width="385px"
height="21px"
loading="lazy"
/>
<p>
The function <em>C(t)</em> 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.
</p>
<p>
However as we've seen in the section on aligning, aligning lets us
simplify things <em>a lot</em>, by completely removing the
contributions of the first coordinate from most mathematical
evaluations, and removing the last <em>y</em> coordinate as well by
virtue of the last point lying on the x-axis. So, while we can
evaluate <em>C(t) = 0</em> for our curve, it'll be much easier to
first axis-align the curve and <em>then</em> evaluating the
curvature function.
</p>
<div className="note">
<h3>Let's derive the full formula anyway</h3>
<p>
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:
</p>
<img
class="LaTeX SVG"
src="images/latex/4d78ebcf8626f777725d67d3672fa480.svg"
width="601px"
height="71px"
loading="lazy"
/>
<p>And of course the same functions for <em>y</em>:</p>
<img
class="LaTeX SVG"
src="images/latex/97b34ad5920612574d1b2a1a9d22d571.svg"
width="399px"
height="69px"
loading="lazy"
/>
<p>
Asking a computer to now compose the <em>C(t)</em> function for us
(and to expand it to a readable form of simple terms) gives us
this rather overly complicated set of arithmetic expressions:
</p>
<img
class="LaTeX SVG"
src="images/latex/b2433959e1f451fa3bf238fc37e04527.svg"
width="552px"
height="97px"
loading="lazy"
/>
<p>
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.
</p>
</div>
<p>
Aligning our curve so that three of the eight coefficients become
zero, we end up with the following simple term function for
<em>C(t)</em>:
</p>
<img
class="LaTeX SVG"
src="images/latex/4dbe6398d0075b5b9ef39458ef620616.svg"
width="560px"
height="21px"
loading="lazy"
/>
<p>
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:
</p>
<img
class="LaTeX SVG"
src="images/latex/d7a657089da19f032dd3b3e1d9ed1d89.svg"
width="509px"
height="73px"
loading="lazy"
/>
<p>
This is a plain quadratic curve, and we know how to solve
<em>C(t) = 0</em>; we use the quadratic formula:
</p>
<img
class="LaTeX SVG"
src="images/latex/b94df866222bed63d123df6b839a4d14.svg"
width="435px"
height="55px"
loading="lazy"
/>
<p>
We can easily compute this value <em>if</em> the discriminator isn't
a negative number (because we only want real roots, not complex
roots), and <em>if</em> <em>x</em> is not zero, because divisions by
zero are rather useless.
</p>
<p>
Taking that into account, we compute <em>t</em>, we disregard any
<em>t</em> value that isn't in the Bézier interval [0,1], and we now
know at which <em>t</em> value(s) our curve will inflect.
</p>
<Graphic
title="Finding cubic Bézier curve inflections"
setup="{this.setupCubic}"
draw="{this.draw}"
/>
</section>
<section id="canonical">
<h1>Canonical form (for cubic curves)</h1>
<p>
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?
</p>
<p>
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
href="http://graphics.pixar.com/people/derose/publications/CubicClassification/paper.pdf"
>"A Geometric Characterization of Parametric Cubic curves"</a
>. 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?
</p>
<p>
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:
</p>
<Graphic
static="{true}"
title="The canonical curve map"
setup="{this.setup}"
draw="{this.drawBase}"
/>
<p>
This is a fairly funky image, so let's see how it breaks down. We
see the three fixed points at (0,0), (0,1) and (1,1), and then the
fourth point is somewhere. Depending on where it is, our curve will
have certain features. Namely, if the fourth point is...
</p>
<ol>
<li>
<p>
anywhere on and in the red zone, the curve will either be
self-intersecting (yielding a loop), or it will have a sharp
discontinuity (yielding a cusp). Anywhere inside the red zone,
this will be a loop. We won't know <em>where</em> that loop is
(in terms of <em>t</em> values), but we are guaranteed that
there is one.
</p>
</li>
<li>
<p>
on the left (red) edge, the curve will have a cusp. We again
don't know <em>where</em>, just that it has one. This edge is
described by the function:
</p>
<img
class="LaTeX SVG"
src="images/latex/63ccae0ebe0ca70dc2afb507ab32e4bd.svg"
width="180px"
height="37px"
loading="lazy"
/>
</li>
<li>
<p>
on the lower right (pink) edge, the curve will have a loop at
t=1, so we know the end coordinate of the curve also lies
<em>on</em> the curve. This edge is described by the function:
</p>
<img
class="LaTeX SVG"
src="images/latex/add5f7fb210a306fe9ff933113f6fb91.svg"
width="231px"
height="39px"
loading="lazy"
/>
</li>
<li>
<p>
on the top (blue) edge, the curve will have a loop at t=0, so we
know the start coordinate of the curve also lies <em>on</em> the
curve. This edge is described by the function:
</p>
<img
class="LaTeX SVG"
src="images/latex/4230e959138d8400e04abf316360009a.svg"
width="153px"
height="37px"
loading="lazy"
/>
</li>
<li>
<p>
inside the green zone, the curve will have a single inflection,
switching concave/convex once.
</p>
</li>
<li>
<p>
between the red and green zones, the curve has two inflections,
meaning its curvature switches between concave/convex form
twice.
</p>
</li>
<li>
<p>
anywhere on the right of the red zone, the curve will have no
inflections. It'll just be a well-behaved arch.
</p>
</li>
</ol>
<p>
Of course, this map is fairly small, but the regions extend to
infinity, with well defined boundaries.
</p>
<div className="note">
<h3>Wait, where do those lines come from?</h3>
<p>
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.
</p>
<p>
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.
</p>
<p>
For the full details, head over to the paper and read through
sections 3 and 4. If you still remember your high school
precalculus, you can probably follow along with this paper,
although you might have to read it a few times before all the bits
"click".
</p>
</div>
<p>
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.
</p>
<p>
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 four 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).
</p>
<p>
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
<em>z</em> 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 <em>z</em> coordinate that is always 1, then we can suddenly
add arbitrary values. For example:
</p>
<img
class="LaTeX SVG"
src="images/latex/d089cc0687982a3302249bb82af3fc16.svg"
width="472px"
height="55px"
loading="lazy"
/>
<p>
Sweet! <em>z</em> 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:
</p>
<img
class="LaTeX SVG"
src="images/latex/058fa85ac31eb666857a860fdedd79df.svg"
width="455px"
height="57px"
loading="lazy"
/>
<p>
Running all our coordinates through this transformation gives a new
set of coordinates, let's call those <strong>U</strong>, 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
<em>x=0</em> line, so what we want is a transformation matrix that,
when we run it, subtracts <em>x</em> from whatever <em>x</em> we
currently have. This is called
<a href="https://en.wikipedia.org/wiki/Shear_matrix">shearing</a>,
and the typical x-shear matrix and its transformation looks like
this:
</p>
<img
class="LaTeX SVG"
src="images/latex/10025fdab2b3fd20f5d389cbe7e3e3ce.svg"
width="195px"
height="53px"
loading="lazy"
/>
<p>
So we want some shearing value that, when multiplied by <em>y</em>,
yields <em>-x</em>, so our x coordinate becomes zero. That value is
simply <em>-x/y</em>, because *-x/y * y = -x*. Done:
</p>
<img
class="LaTeX SVG"
src="images/latex/20684d22b3ddc52fd6abde8ce56608a9.svg"
width="133px"
height="67px"
loading="lazy"
/>
<p>
Now, running this on all our points generates a new set of
coordinates, let's call those <strong>V</strong>, 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
<a href="https://en.wikipedia.org/wiki/Scaling_%28geometry%29"
>do some scaling</a
>
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/point3<sub>x</sub>, and y-scaling by
point2<sub>y</sub>. This is really easy:
</p>
<img
class="LaTeX SVG"
src="images/latex/8cbef24b8c3b26f9daf2f89d27d36e95.svg"
width="137px"
height="71px"
loading="lazy"
/>
<p>
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 <code>y/x</code> rather than <code>x/y</code> 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:
</p>
<img
class="LaTeX SVG"
src="images/latex/0430e8c7f7d4ec80e6527f96f3d56e5c.svg"
width="140px"
height="65px"
loading="lazy"
/>
<p>
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:
</p>
<img
class="LaTeX SVG"
src="images/latex/2f85d84f0e3dd14cc25e48583aed3822.svg"
width="455px"
height="91px"
loading="lazy"
/>
<p>
That looks very complex, but notice that every coordinate value is
being offset by the initial translation, and a lot of terms in there
repeat: it's pretty easy to calculate this fast, since there's so
much we can cache and reuse while we compute this mapped coordinate!
</p>
<p>
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?
</p>
<img
class="LaTeX SVG"
src="images/latex/83262761bb7fa9b832fe483ded436973.svg"
width="353px"
height="59px"
loading="lazy"
/>
<p>
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:
</p>
<img
class="LaTeX SVG"
src="images/latex/f3261ad2802d980ebe6e35b272375700.svg"
width="408px"
height="40px"
loading="lazy"
/>
<p>
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!
</p>
<div className="note">
<h3>How do you track all that?</h3>
<p>
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
<a href="http://www.wolfram.com/mathematica">Mathematica</a>.
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,
<a
href="http://pomax.github.io/gh-weblog-2/downloads/canonical-curve.nb"
>here's</a
>
the Mathematica notebook if you want to see how this works for
yourself.
</p>
<p>
Now, I know, you're thinking "but Mathematica is super expensive!"
and that's true, it's
<a href="http://www.wolfram.com/mathematica-home-edition"
>$295 for home use</a
>, but it's <strong>also</strong>
<a href="http://www.wolfram.com/raspberry-pi"
>free when you buy a $35 raspberry pi</a
>. 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
<em>do</em>, 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.
</p>
</div>
<p>
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:
</p>
<Graphic
title="A cubic curve mapped to canonical form"
setup="{this.setup}"
draw="{this.draw}"
/>
</section>
<section id="yforx">
<h1>Finding Y, given X</h1>
<p>
One common task that pops up in things like CSS work, or parametric
equalisers, or image leveling, or any other number of applications
where Bezier 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 <em>is</em> possible and as long as you have
some code in place to help, it's not a lot of a work either.
</p>
<p>
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 Bezier 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 mouse over the graphic, you will see a red line drawn
that corresponds to the <code>x</code> coordinate under your cursor:
vertical in the left graphic, horizontal in the right.
</p>
<Graphic
title="Finding t, given x=x(t). Left: our curve, right: the x=x(t) function"
setup="{this.tforx.setup}"
draw="{this.tforx.draw}"
onMouseMove="{this.onMouseMove}"
/>
<p>
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). So... can we just use root finding
to get the "t" value associated with an "x" value? To which the
answer should unsurprisingly be " yes, we can". Sure, we'll need to
compute cubic roots, but we've already seen how to do that in the
section on <a href="#extremities">finding extremities</a>, and we're
not dealing with a complicated case: we <em>know</em> there is only
root, so let's go and find it!
</p>
<p>First, let's look at the function for x(t):</p>
<img
class="LaTeX SVG"
src="images/latex/9df91c28af38c1ba2e2d38d2714c9446.svg"
width="335px"
height="19px"
loading="lazy"
/>
<p>
We can rewrite this to a plain polynomial form, by just fully
writing out the expansion and then collecting the polynomial
factors, as:
</p>
<img
class="LaTeX SVG"
src="images/latex/9ab2b830fe7fb73350c19bde04e9441b.svg"
width="445px"
height="19px"
loading="lazy"
/>
<p>
Nothing special here: that's a standard cubic polynomial in "power"
form (i.e. all the terms are ordered by their power of
<code>t</code>). So, given that <code>a</code>, <code>b</code>,
<code>c</code>, <code>d</code>, <em>and</em> <code>x(t)</code> are
all known constants, we can trivially rewrite this (by moving the
<code>x(t)</code> across the equal sign) as:
</p>
<img
class="LaTeX SVG"
src="images/latex/de3bd3e271d72194c730d0ae44f031a8.svg"
width="465px"
height="19px"
loading="lazy"
/>
<p>
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: we
know everything except <code>t</code>, we just need some
mathematical insight to tell us how to do this.
</p>
<p>
...Of course "just" is not the right qualifier here, there is
nothing "just" about finding the roots of a cubic function, but
thankfully we've already covered the tool to do this in the
<a href="#extremities">root finding</a> section: Gerolano Cardano's
solution for cubic roots. Of course, we still need to be a bit
careful, because cubic roots are complicated things: you can get up
to three roots back, even though we only "want" one root. In our
case, only one will be both a real number (as opposed to a complex
number) <em>and</em> lie in the <code>t</code> interval [0,1], so we
need to filter for that:
</p>
<pre><code>double x = some value we know!
double[] roots = getTforX(x);
double t;
if (roots.length &gt; 0) {
for (double _t: roots) {
if (_t&lt;0 || _t&gt;1) continue;
t = _t;
break;
}
}</code></pre>
<p>
And that's it, we're done: we now have the <code>t</code> value
corresponding to our <code>x</code>, and we can just evaluate our
curve for that <code>t</code> value to find a coordinate that has
our original, known <code>x</code>, and our unknown
<code>y</code> value. Mouse over the following graphic, which
performs this root finding based on a knowledge of which "x" value
we're interested in.
</p>
<Graphic
title="Finding y(t), by finding t, given x=x(t)"
setup="{this.yforx.setup}"
draw="{this.yforx.draw}"
onMouseMove="{this.onMouseMove}"
/>
</section>
<section id="arclength">
<h1>Arc length</h1>
<p>
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 <em>f<sub>x</sub>(t)</em> and
<em>f<sub>y</sub>(t)</em>, then the length of the curve, measured
from start point to some point <em>t = z</em>, is computed using the
following seemingly straight forward (if a bit overwhelming)
formula:
</p>
<img
class="LaTeX SVG"
src="images/latex/cb24cda7f7f4bbf3be7104c460e0ec9f.svg"
width="140px"
height="33px"
loading="lazy"
/>
<p>or, more commonly written using Leibnitz notation as:</p>
<img
class="LaTeX SVG"
src="images/latex/d3003177813309f88f58a1f515f5df9f.svg"
width="245px"
height="35px"
loading="lazy"
/>
<p>
This formula says that the length of a parametric curve is in fact
equal to the <strong>area</strong> 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
<a
href="http://www.wolframalpha.com/input/?i=antiderivative+for+sqrt((2*(1-t)*t*B+%2B+t%5E2*C)%27%5E2+%2B+(2*(1-t)*t*E)%27%5E2)&incParTime=true"
>unwieldy computation</a
>, 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:
<strong
><em
>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"</em
></strong
>.
</p>
<p>
Seriously:
<a href="https://en.wikipedia.org/wiki/Abel%E2%80%93Ruffini_theorem"
>It cannot be done</a
>.
</p>
<p>
So we turn to numerical approaches again. The method we'll look at
here is the
<a
href="http://www.youtube.com/watch?v=unWguclP-Ds&feature=BFa&list=PLC8FC40C714F5E60F&index=1"
>Gauss quadrature</a
>. This approximation is a really neat trick, because for any
<em>n<sup>th</sup></em> 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:
</p>
<img
class="LaTeX SVG"
src="images/latex/e168758d35b8f6781617eda5a32b20bf.svg"
width="636px"
height="71px"
loading="lazy"
/>
<p>
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:
</p>
<div className="figure">
<Graphic
inline="{true}"
static="{true}"
title="A function's approximated integral"
setup="{this.setup}"
draw="{this.drawCoarseIntegral}"
/>
<Graphic
inline="{true}"
static="{true}"
title="A better approximation"
setup="{this.setup}"
draw="{this.drawFineIntegral}"
/>
<Graphic
inline="{true}"
static="{true}"
title="An even better approximation"
setup="{this.setup}"
draw="{this.drawSuperFineIntegral}"
/>
</div>
<p>
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.
</p>
<p>
So, the trick is to come up with useful rectangular strips. A naive
way is to simply create <em>n</em> strips, all with the same width,
but there is a far better way using special values for
<em>C</em> and <em>f(t)</em> depending on the value of <em>n</em>,
which indicates how many strips we'll use, and it's called the
Legendre-Gauss quadrature.
</p>
<p>
This approach uses strips that are <em>not</em> 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.
</p>
<div className="note">
<p>
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 <em>z</em>). 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:
</p>
<img
class="LaTeX SVG"
src="images/latex/5509919419288129322cfbd4c60d0a4f.svg"
width="341px"
height="72px"
loading="lazy"
/>
<p>
That may look a bit more complicated, but the fraction involving
<em>z</em> is a fixed number, so the summation, and the evaluation
of the <em>f(t)</em> values are still pretty simple.
</p>
<p>
So, what do we need to perform this calculation? For one, we'll
need an explicit formula for <em>f(t)</em>, because that
derivative notation is handy on paper, but not when we have to
implement it. We'll also need to know what these
<em>C<sub>i</sub></em> and <em>t<sub>i</sub></em> values should
be. Luckily, that's less work because there are actually many
tables available that give these values, for any <em>n</em>, so if
we want to approximate our integral with only two terms (which is
a bit low, really) then
<a href="legendre-gauss.html">these tables</a> would tell us that
for <em>n=2</em> we must use the following values:
</p>
<img
class="LaTeX SVG"
src="images/latex/d0d93f1cc26b560309dade1f1aa012f2.svg"
width="63px"
height="93px"
loading="lazy"
/>
<p>
Which means that in order for us to approximate the integral, we
must plug these values into the approximate function, which gives
us:
</p>
<img
class="LaTeX SVG"
src="images/latex/e96dd431f6ef9433ccf25909dddd5bca.svg"
width="476px"
height="44px"
loading="lazy"
/>
<p>
We can program that pretty easily, provided we have that
<em>f(t)</em> available, which we do, as we know the full
description for the Bézier curve functions B<sub>x</sub>(t) and
B<sub>y</sub>(t).
</p>
</div>
<p>
If we use the Legendre-Gauss values for our <em>C</em> values
(thickness for each strip) and <em>t</em> 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?
</p>
<Graphic
title="Arc length for a Bézier curve"
setup="{this.setupCurve}"
draw="{this.drawCurve}"
/>
</section>
<section id="arclengthapprox">
<h1>Approximated arc length</h1>
<p>
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.
</p>
<p>
If we combine the work done in the previous sections on curve
flattening and arc length computation, we can implement these with
minimal effort:
</p>
<Graphic
title="Approximate quadratic curve arc length"
setup="{this.setupQuadratic}"
draw="{this.draw}"
onKeyDown="{this.props.onKeyDown}"
/>
<Graphic
title="Approximate cubic curve arc length"
setup="{this.setupCubic}"
draw="{this.draw}"
onKeyDown="{this.props.onKeyDown}"
/>
<p>
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 may notice that the error in length is actually pretty
significant, even if the percentage is fairly low: if the number of
segments used yields an error of 0.1% or higher, the flattened curve
already looks fairly obviously flattened. And of course, the longer
the curve, the more significant the error will be.
</p>
</section>
<section id="curvature">
<h1>Curvature of a curve</h1>
<p>
Imagine 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 leaving one curve and entering the next, 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
<a href="https://en.wikipedia.org/wiki/Curvature"
><em>curvature</em></a
>
at the transition from one curve to the next "looks good". So, we
could have them share an end coordinate, and then ensure that the
derivatives for both curves match at that coordinate, and at a
casual glance, that seems the perfect solution: if we make the
derivatives match, then both the "direction" in which we travel from
one curve to the next is the same, and the "speed" at which we
travel the curve will be the same.
</p>
<p>Problem solved!</p>
<p>
But, if we think about this a little more, this cannot possible
work, because of something that you may have noticed in the section
on <a href="#reordering">reordering curves</a>: what a curve looks
like, and the function that draws that curve, are not in some kind
of universal, fixed, one-to-one relation. If we have some quadratic
curve, then simply by raising the curve order we can get
corresponding cubic, quartic, and higher and higher mathematical
expressions that all draw the <em>exact same curve</em> but with
wildly different derivatives. So: if we want to make a transition
from one curve to the next look good, and we want to use the
derivative, then we suddenly need to answer the question: "Which
derivative?".
</p>
<p>
How would you even decide? What makes the cubic derivatives better
or less suited than, say, quintic derivatives? Wouldn't it be nicer
if we could use something that was inherent to the curve, without
being tied to the functions that yield that curve? And (of course)
as it turns out, there is a way to define curvature in such a way
that it only relies on what the curve actually looks like, and given
where this section is in the larger body of this Primer, it should
hopefully not be surprising that we thing we can use to define
curvature is the thing we talked about in the previous section: arc
length.
</p>
<p>
Intuitively, this should make sense, even if we have no idea what
the maths would look like: if we travel some fixed distance along
some curve, then the point at that distance is simply the point at
that distance. It doesn't matter what function we used to draw the
curve: once we know what the curve looks like, the function(s) used
to draw it become irrelevant: a point a third along the full
distance of the curve is simply the point a third along the distance
of the curve.
</p>
<p>
You might think that in order to find the curvature of a curve, we
now need to find and then solve the arc length function, and that
would be a problem because we just saw that there is no way to
actually do that: don't worry, we don't. We do need to know the
<em>form</em> of the arc length function, which we saw above, but
it's not the thing we're actually interested in, and we're going to
be rewriting it in a way that makes most of the crazy complex things
about it just... disappear.
</p>
<p>
In fact, after
<a href="http://mathworld.wolfram.com/Curvature.html"
>running through the steps necessary</a
>
to determine what we're left with if we use the arclength function's
derivative (with another run-through of the maths
<a href="https://math.stackexchange.com/a/275324/71940">here</a>),
rather than the curve's original function's derivative, then the
integral disappears entirely (because of the
<a
href="https://en.wikipedia.org/wiki/Fundamental_theorem_of_calculus"
>fundamental therem of calculus</a
>), and we're left with 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.
</p>
<p>
Let me just highlight that before we move on: we calculate the
curvature of a curve using the arc length function derivative,
because the original function's derivative is entirely unreliable,
and in doing so we end up with a formula that expresses curvature in
terms of the original function's derivatives.
</p>
<p><em>That's crazy!</em></p>
<p>
But, that's what makes maths such an interesting thing: it can show
you that all your assumptions are completely wrong, only to then go
"but actually, you were on the right track all along, here: ..."
with a solution that is so easy to work with as to almost seem
mundane. So: enough of all this text, how do we calculate curvature?
What is the function for κ? Concisely, the function is this:
</p>
<img
class="LaTeX SVG"
src="images/latex/d9c893051586eb8d9de51c0ae1ef8fae.svg"
width="113px"
height="47px"
loading="lazy"
/>
<p>
Which is really just a "short form" that glosses over the fact that
we're dealing with functions:
</p>
<img
class="LaTeX SVG"
src="images/latex/828333034b4fed8e248683760d6bc6f4.svg"
width="239px"
height="55px"
loading="lazy"
/>
<p>
And while that's a litte more verbose, it's still just as simple to
work with as the first function: the curvature at some point on any
(and this cannot be overstated: <em>any</em>) 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
<a href="#derivatives"
>we know what the first and second derivatives are</a
>, and so evaluating this function for any <strong>t</strong> value
is just a matter of basic arithematics.
</p>
<div className="howtocode">
<h3>Implement the kappa function</h3>
<p>In fact, let's just implement it right now:</p>
<pre><code>function kappa(t, B):
d = B.getDerivative()
dd = d.getDerivative()
dx = d.getX(t)
dy = d.getY(t)
ddx = dd.getX(t)
ddy = dd.getY(t)
numerator = dx * ddy - ddx * dy
denominator = pow(dx*dx + dy*dy, 1.5)
return numerator / denominator</code></pre>
<p>That was easy!</p>
<p>
In fact, it stays easy because 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 relation in this entire primer:
</p>
<img
class="LaTeX SVG"
src="images/latex/6ed4fd2ead35c57984caddf9fe375a5f.svg"
width="81px"
height="37px"
loading="lazy"
/>
<p>So that's a rather convenient fact to know, too.</p>
</div>
<p>
So 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 formulat 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 looking
transition we could ask for.
</p>
<Graphic
title="Matching curvatures for a quadratic and cubic Bézier curve"
setup="{this.setup}"
draw="{this.draw}"
/>
<p>
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. In your
code you typically solve this by matching absolute values, but
that's not super easy to program visually... however, we
<em>can</em> just show the curvature on both sides of the curve,
making lining things up a bit easier:
</p>
<Graphic
title="(Easier) curvature matching for a quadratic and cubic Bézier curve"
setup="{this.setup}"
draw="{this.drawOmni}"
/>
</section>
<section id="tracing">
<h1>Tracing a curve at fixed distance intervals</h1>
<p>
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.
</p>
<p>Now you have a problem.</p>
<p>
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 <em>t</em> 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 "<em>t</em> value", by plotting them against one another.
</p>
<p>
The following graphic shows a particularly illustrative curve, and
it's length-to-t plot. For linear traversal, this line needs to be
straight, running from (0,0) to (length,1). This is, it's safe to
say, not what we'll see, we'll see something wobbly instead. To make
matters even worse, the length-to-<em>t</em> 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
once at best, the plot shows that the distance function along the
curve is able to switch forms three times (to see this, try creating
an S curve with the start/end close together, but the control points
far apart).
</p>
<Graphic
title="The t-for-distance function"
setup="{this.setup}"
draw="{this.plotOnly}"
/>
<p>
We see a function that might be invertible, but we won't be able to
do so, symbolically. You may remember from the section on arc length
that we cannot actually compute the true arc length function as an
expression of <em>t</em>, which means we also can't compute the true
inverted function that gives <em>t</em> as an expression of length.
So how do we fix this?
</p>
<p>
One way is to do what the graphic does: simply run through the
curve, determine its <em>t</em>-for-length values as a set of
discrete values at some high resolution (the graphic uses 100
discrete points), and then use those as a basis for finding an
appropriate <em>t</em> value, given a distance along the curve. This
works quite well, actually, and is fairly fast.
</p>
<p>
We can use some colour to show the difference between distance-based
and time based intervals: the following graph is similar to the
previous one, except it segments the curve in terms of
equal-distance intervals. This shows as regular colour intervals
going down the graph, but the mapping to <em>t</em> values is not
linear, so there will be (highly) irregular intervals along the
horizontal axis. It also shows the curve in an alternating colouring
based on the t-for-distance values we find our LUT:
</p>
<Graphic
title="Fixed-interval coloring a curve"
setup="{this.setup}"
draw="{this.drawColoured}"
onKeyDown="{this.props.onKeyDown}"
/>
<p>
Use your up and down arrow keys to increase or decrease the number
of equidistant segments used to colour the curve.
</p>
<p>
However, are there better ways? One such way is discussed in "<a
href="http://www.geometrictools.com/Documentation/MovingAlongCurveSpecifiedSpeed.pdf"
>Moving Along a Curve with Specified Speed</a
>" 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).
</p>
</section>
<section id="intersections">
<h1>Intersections</h1>
<p>
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 <em>might</em> 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.
</p>
<p>
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.
</p>
<h3>Line-line intersections</h3>
<p>
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
<a
href="http://www.topcoder.com/tc?module=Static&d1=tutorials&d2=geometry2#line_line_intersection"
>top coder</a
>
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.
</p>
<p>
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).
</p>
<Graphic
title="Line/line intersections"
setup="{this.setupLines}"
draw="{this.drawLineIntersection}"
/>
<div className="howtocode">
<h3>Implementing line-line intersections</h3>
<p>
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:
</p>
<pre><code>lli8 = function(x1,y1,x2,y2,x3,y3,x4,y4):
var nx=(x1*y2-y1*x2)*(x3-x4)-(x1-x2)*(x3*y4-y3*x4),
ny=(x1*y2-y1*x2)*(y3-y4)-(y1-y2)*(x3*y4-y3*x4),
d=(x1-x2)*(y3-y4)-(y1-y2)*(x3-x4);
if d=0:
return false
return point(nx/d, ny/d)
lli4 = function(p1, p2, p3, p4):
var x1 = p1.x, y1 = p1.y,
x2 = p2.x, y2 = p2.y,
x3 = p3.x, y3 = p3.y,
x4 = p4.x, y4 = p4.y;
return lli8(x1,y1,x2,y2,x3,y3,x4,y4)
lli = function(line1, line2):
return lli4(line1.p1, line1.p2, line2.p1, line2.p2)</code></pre>
</div>
<h3>What about curve-line intersections?</h3>
<p>
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.
</p>
<Graphic
title="Quadratic curve/line intersections"
setup="{this.setupQuadratic}"
draw="{this.draw}"
/>
<Graphic
title="Cubic curve/line intersections"
setup="{this.setupCubic}"
draw="{this.draw}"
/>
<p>
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.
</p>
</section>
<section id="curveintersection">
<h1>Curve/curve intersection</h1>
<p>
Using de Casteljau's algorithm to split the curve we can now
implement curve/curve intersection finding using a "divide and
conquer" technique:
</p>
<ul>
<li>
Take two curves <em>C<sub>1</sub></em> and <em>C<sub>2</sub></em
>, and treat them as a pair.
</li>
<li>
If their bounding boxes overlap, split up each curve into two
sub-curves
</li>
<li>
With <em>C<sub>1.1</sub></em
>, <em>C<sub>1.2</sub></em
>, <em>C<sub>2.1</sub></em> and <em>C<sub>2.2</sub></em
>, form four new pairs (<em>C<sub>1.1</sub></em
>,<em>C<sub>2.1</sub></em
>), (<em>C<sub>1.1</sub></em
>, <em>C<sub>2.2</sub></em
>), (<em>C<sub>1.2</sub></em
>,<em>C<sub>2.1</sub></em
>), and (<em>C<sub>1.2</sub></em
>,<em>C<sub>2.2</sub></em
>).
</li>
<li>
For each pair, check whether their bounding boxes overlap.
<ul>
<li>
If their bounding boxes do not overlap, discard the pair, as
there is no intersection between this pair of curves.
</li>
<li>
If there <em>is</em> overlap, rerun all steps for this pair.
</li>
</ul>
</li>
<li>
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"
<code>t</code> value (we can either throw all but one away, we can
average the cluster's <code>t</code> values, or you can do
something even more creative).
</li>
</ul>
<p>
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.
</p>
<p>
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. The algorithm
resets once it's found a solution, 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!)
</p>
<Graphic
title="Curve/curve intersections"
setup="{this.setup}"
draw="{this.draw}"
>
<button onClick="{this.stepUp}">advance one step</button>
</Graphic>
<p>
Self-intersection is dealt with in the same way, except we turn a
curve into two or more curves first based on the inflection points.
We then form all possible curve pairs with the resultant segments,
and run exactly the same algorithm. All non-overlapping curve pairs
will be removed after the first iteration, and the remaining steps
home in on the curve's self-intersection points.
</p>
</section>
<section id="abc">
<h1>The projection identity</h1>
<p>
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
"mould" a curve, by picking it up at some point, and dragging that
point around to change the curve's shape.
</p>
<p>
How does that work? Succinctly: we run de Casteljau's algorithm in
reverse!
</p>
<p>
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 <em>t</em> 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.
</p>
<p>
So let's use graphics instead of text to see where this "A" is,
because text only gets us so far: in the following graphic, click
anywhere on the curves to see the identity information that we'll be
using to run de Casteljau in reverse (you can manipulate the curve
even after picking a point. Note the "ratio" value when you do so:
does it change?):
</p>
<div className="figure">
<Graphic
inline="{true}"
title="Projections in a quadratic Bézier curve"
setup="{this.setupQuadratic}"
draw="{this.draw}"
onClick="{this.onClick}"
/>
<Graphic
inline="{true}"
title="Projections in a cubic Bézier curve"
setup="{this.setupCubic}"
draw="{this.draw}"
onClick="{this.onClick}"
/>
</div>
<p>Clicking anywhere on the curves shows us three things:</p>
<ol>
<li>our on-curve point; let's call that <b>B</b>,</li>
<li>
a point at the tip of B's "hat", on de Casteljau step up; let's
call that <b>A</b>, and
</li>
<li>
a point that we get by projecting B onto the start--end baseline;
let's call that <b>C</b>.
</li>
</ol>
<p>
These three values A, B, and C hide an important identity formula
for quadratic and cubic Bézier curves: for any point on the curve
with some <em>t</em> value, the ratio distance of C along the
baseline is fixed: if some <em>t</em> 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
<em>t</em> value, C will <em>always</em> 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. The following
function stays true:
</p>
<img
class="LaTeX SVG"
src="images/latex/34fe255294faf45ab02128f7997b92ce.svg"
width="197px"
height="16px"
loading="lazy"
/>
<p>So that just leaves finding A.</p>
<div className="note">
<p>
While that relation is fixed, the function <em>u(t)</em> differs
depending on whether we're working with quadratic or cubic curves:
</p>
<img
class="LaTeX SVG"
src="images/latex/62f2f984e43a22a6b4bda4d399dedfc6.svg"
width="197px"
height="87px"
loading="lazy"
/>
<p>
So, if we know the start and end coordinates, and we know the
<em>t</em> value, we know C:
</p>
<div className="figure">
<Graphic
inline="{true}"
title="Quadratic value of C for t"
draw="{this.drawQCT}"
onMouseMove="{this.setCT}"
/>
<Graphic
inline="{true}"
title="Cubic value of C for t"
draw="{this.drawCCT}"
onMouseMove="{this.setCT}"
/>
</div>
<p>
Mouse-over the graphs to see the expression for C, given the
<em>t</em> value at the mouse pointer.
</p>
</div>
<p>
There's also another important bit of information that is inherent
to the ABC values: while the distances between A and B, and B and C,
are dynamic (based on where we put B), the <em>ratio</em> between
the two distances is stable. Given some <em>t</em> value, the
following always holds:
</p>
<img
class="LaTeX SVG"
src="images/latex/385d1fd4aecbd2066e6e284a84408be6.svg"
width="251px"
height="39px"
loading="lazy"
/>
<p>
This leads to a pretty powerful bit of knowledge: merely by knowing
the <em>t</em> value of some on curve point, we know where C has to
be (as per the above note), and because we know B and C, and thus
have the distance between them, we know where A has to be:
</p>
<img
class="LaTeX SVG"
src="images/latex/12aaf0d7fd20b3c551a0ec76b18bd7d2.svg"
width="217px"
height="37px"
loading="lazy"
/>
<p>And that's it, all values found.</p>
<div className="note">
<p>
Much like the <em>u(t)</em> function in the above note, the
<em>ratio(t)</em> function depends on whether we're looking at
quadratic or cubic curves. Their form is intrinsically related to
the <em>u(t)</em> function in that they both come rolling out of
the same function evaluation, explained over on
<a
href="http://mathoverflow.net/questions/122257/finding-the-formula-for-B%C3%A9zier-curve-ratios-hull-point-point-baseline"
>MathOverflow</a
>
by Boris Zbarsky and myself. The ratio functions are the "s(t)"
functions from the answers there, while the "u(t)" functions have
the same name both here and on MathOverflow.
</p>
<img
class="LaTeX SVG"
src="images/latex/059000c5c8a37dcc8d7fa04154a05df3.svg"
width="245px"
height="41px"
loading="lazy"
/>
<img
class="LaTeX SVG"
src="images/latex/b4987e9b77b0df604238b88596c5f7c3.svg"
width="223px"
height="41px"
loading="lazy"
/>
<p>
Unfortunately, this trick only works for quadratic and cubic
curves. Once we hit higher order curves, things become a lot less
predictable; the "fixed point <em>C</em>" is no longer fixed,
moving around as we move the control points, and projections of
<em>B</em> onto the line between start and end may actually lie on
that line before the start, or after the end, and there are no
simple ratios that we can exploit.
</p>
</div>
<p>
So: if we know B and its corresponding <em>t</em> value, then we
know all the ABC values, which —together with a start and end
coordinate— gives us the necessary information to reconstruct a
curve's "de Casteljau skeleton", which means that two points and a
value between 0 and 1, we can come up with a curve. And that opens
up possibilities: curve manipulation by dragging an on-curve point,
as well as curve fitting of "a bunch of coordinates". These are
useful things, and we'll look at both in the next sections.
</p>
</section>
<section id="moulding">
<h1>Manipulating a curve</h1>
<p>
Armed with knowledge of the "ABC" relation, we can now update a
curve interactively, by letting people click anywhere on the curve,
find the <em>t</em>-value matching that coordinate, and then letting
them drag that point around. With every drag update we'll have a new
point "B", which we can combine with the fixed point "C" to find our
new point A. Once we have those, we can reconstruct the de Casteljau
skeleton and thus construct a new curve with the same start/end
points as the original curve, passing through the user-selected
point B, with correct new control points.
</p>
<Graphic
title="Moulding a quadratic Bézier curve"
setup="{this.setupQuadratic}"
draw="{this.drawMould}"
onClick="{this.placeMouldPoint}"
onMouseDown="{this.markQB}"
onMouseDrag="{this.dragQB}"
onMouseUp="{this.saveCurve}"
/>
<p>
<strong>Click-dragging the curve itself</strong> shows what we're
using to compute the new coordinates: while dragging you will see
the original point B and its corresponding <i>t</i>-value, the
original point C for that <i>t</i>-value, as well as the new point
B' based on the mouse cursor. Since we know the <i>t</i>-value for
this configuration, we can compute the ABC ratio for this
configuration, and we know that our new point A' should like at a
distance:
</p>
<img
class="LaTeX SVG"
src="images/latex/7bba0a4fd605e023cda922de125b3e32.svg"
width="221px"
height="36px"
loading="lazy"
/>
<p>
For quadratic curves, this means we're done, since the new point A'
is equivalent to the new quadratic control point. For cubic curves,
we need to do a little more work:
</p>
<Graphic
title="Moulding a cubic Bézier curve"
setup="{this.setupCubic}"
draw="{this.drawMould}"
onClick="{this.placeMouldPoint}"
onMouseDown="{this.markCB}"
onMouseDrag="{this.dragCB}"
onMouseUp="{this.saveCurve}"
/>
<p>
To help understand what's going on, the cubic graphic shows the full
de Casteljau construction "hull" when repositioning point B. We
compute A' in exactly the same way as before, but we also record the
final strut line that forms B in the original curve. Given A', B',
and the endpoints e1 and e2 of the strut line relative to B', we can
now compute where the new control points should be. Remember that B'
lies on line e1--e2 at a distance <i>t</i>, because that's how
Bézier curves work. In the same manner, we know the distance A--e1
is only line-interval [0,t] of the full segment, and A--e2 is only
line-interval [t,1], so constructing the new control points is
fairly easy.
</p>
<p>First, we construct the one-level-of-de-Casteljau-up points:</p>
<img
class="LaTeX SVG"
src="images/latex/524206c49f317d27d8e07a310b24a7a3.svg"
width="132px"
height="75px"
loading="lazy"
/>
<p>And then we can compute the new control points:</p>
<img
class="LaTeX SVG"
src="images/latex/94f61d17f896aebddcf5a7c676aee7d1.svg"
width="156px"
height="75px"
loading="lazy"
/>
<p>And that's cubic curve manipulation.</p>
</section>
<section id="pointcurves">
<h1>Creating a curve from three points</h1>
<p>
Given the preceding section on curve manipulation, we can also
generate quadratic and cubic curves from any three points. However,
unlike circle-fitting, which requires just three points, Bézier
curve fitting requires three points, as well as a <em>t</em> value,
so we can figure out where point 'C' needs to be.
</p>
<p>
The following graphic lets you place three points, and will use the
preceding sections on the ABC ratio and curve construction to form a
quadratic curve through them. You can move the points you've placed
around by click-dragging, or try a new curve by drawing new points
with pure clicks. (There's some freedom here, so for illustrative
purposes we clamped <em>t</em> to simply be 0.5, lets us bypass some
maths, since a <em>t</em> value of 0.5 always puts C in the middle
of the start--end line segment)
</p>
<Graphic
title="Fitting a quadratic Bézier curve"
setup="{this.setup}"
draw="{this.drawQuadratic}"
onClick="{this.onClick}"
/>
<p>
For cubic curves we also need some values to construct the "de
Casteljau line through B" with, and that gives us quite a bit of
choice. Since we've clamped <em>t</em> to 0.5, we'll set up a line
through B parallel to the line start--end, with a length that is
proportional to the length of the line B--C: the further away from
the baseline B is, the wider its construction line will be, and so
the more "bulby" the curve will look. This still gives us some
freedom in terms of exactly how to scale the length of the
construction line as we move B closer or further away from the
baseline, so I simply picked some values that sort-of-kind-of look
right in that if a circle through (start,B,end) forms a perfect
hemisphere, the cubic curve constructed forms something close to a
hemisphere, too, and if the points lie on a line, then the curve
constructed has the control points very close to B, while still
lying between B and the correct curve end point:
</p>
<Graphic
title="Fitting a cubic Bézier curve"
setup="{this.setup}"
draw="{this.drawCubic}"
onClick="{this.onClick}"
/>
<p>
In each graphic, the blue parts are the values that we "just have"
simply by setting up our three points, combined with our decision on
which <em>t</em> value to use (and construction line orientation and
length for cubic curves). There are of course many ways to determine
a combination of <em>t</em> and tangent values that lead to a more
"æsthetic" curve, but this will be left as an exercise to the
reader, since there are many, and æsthetics are often quite
personal.
</p>
</section>
<section id="curvefitting">
<h1>Curve fitting</h1>
<p>
Given the previous section, one question you might have is "what if
I don't want to guess <code>t</code> 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
<code>t</code> values all on its own?
</p>
<p>
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
<a href="https://en.wikipedia.org/wiki/Least_squares"
>least squares</a
>
<a href="https://en.wikipedia.org/wiki/Polynomial_regression"
>polynomial regression</a
>. 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 <code>nth</code> order curve, what are the
coordinates we can find such that our curve is "off" by the least
amount?".
</p>
<p>
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
<em>squared distance</em> 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 <em>well</em> beyond the scope of this
material.
</p>
<p>
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
<a
href="https://web.archive.org/web/20180403213813/http://jimherold.com/2012/04/20/least-squares-bezier-fit/"
>"Least Squares Bézier Fit"</a
>
article, and end with some nice interactive graphics for doing some
curve fitting.
</p>
<p>
Before we begin, we're going to use the curve in matrix form. In the
<a href="#matrix">section on matrices</a>, 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.
</p>
<p>
As such, the first step in the process is expressing our Bézier
curve as powers/coefficients/coordinate matrix
<strong>T x M x C</strong>, by expanding the Bézier functions.
</p>
<div className="note">
<h2>Revisiting the matrix representation</h2>
<p>
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:
</p>
<img
class="LaTeX SVG"
src="images/latex/bd8e8e294eec10d2bf6ef857c7c0c2c2.svg"
width="295px"
height="43px"
loading="lazy"
/>
<p>
And then we (trivially) rearrange the terms across multiple lines:
</p>
<img
class="LaTeX SVG"
src="images/latex/2b8334727d3b004c6e87263fec6b32b7.svg"
width="216px"
height="64px"
loading="lazy"
/>
<p>
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").
</p>
<p>
With that arrangement, we can easily decompose this as a matrix
multiplication:
</p>
<img
class="LaTeX SVG"
src="images/latex/94acb5850778dcb16c2ba3cfa676f537.svg"
width="572px"
height="53px"
loading="lazy"
/>
<p>
We can do the same for the cubic curve, of course. We know the
base function for cubics:
</p>
<img
class="LaTeX SVG"
src="images/latex/7eada6f12045423de24d9a2ab8e293b1.svg"
width="355px"
height="19px"
loading="lazy"
/>
<p>So we write out the expansion and rearrange:</p>
<img
class="LaTeX SVG"
src="images/latex/875ca8eea72e727ccb881b4c0b6a3224.svg"
width="248px"
height="87px"
loading="lazy"
/>
<p>Which we can then decompose:</p>
<img
class="LaTeX SVG"
src="images/latex/03ec73258d5c95eed39a2ea8665e0b07.svg"
width="404px"
height="72px"
loading="lazy"
/>
<p>
And, of course, we can do this for quartic curves too (skipping
the expansion step):
</p>
<img
class="LaTeX SVG"
src="images/latex/9151c0fdf9689ee598a2d029ab2ffe34.svg"
width="491px"
height="92px"
loading="lazy"
/>
<p>
And so and on so on. Now, let's see how to use these
<strong>T</strong>, <strong>M</strong>, and <strong>C</strong>, to
do some curve fitting.
</p>
</div>
<p>
Let's get started: we're going to assume we picked the right order
curve: for <code>n</code> points we're fitting an <code>n-1</code
><sup>th</sup> order curve, so we "start" with a vector
<strong>P</strong> that represents the coordinates we already know,
and for which we want to do curve fitting:
</p>
<img
class="LaTeX SVG"
src="images/latex/2bef3da3828d63d690460ce9947dbde2.svg"
width="63px"
height="73px"
loading="lazy"
/>
<p>
Next, we need to figure out appropriate <code>t</code> 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
<code>t</code> values), but in this case let's look at two "obvious"
choices:
</p>
<ol>
<li>equally spaced <code>t</code> values, and</li>
<li>
<code>t</code> values that align with distance along the polygon.
</li>
</ol>
<p>
The first one is really simple: if we have <code>n</code> points,
then we'll just assign each point <code>i</code> a
<code>t</code> value of <code>(i-1)/(n-1)</code>. So if we have four
points, the first point will have <code>t=(1-1)/(4-1)=0/3</code>,
the second point will have <code>t=(2-1)/(4-1)=1/3</code>, the third
point will have <code>t=2/3</code>, and the last point will be
<code>t=1</code>. We're just straight up spacing the
<code>t</code> values to match the number of points we have.
</p>
<p>
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
<code>t</code> is equal to the distance along the polygon, scaled to
the [0,1] domain.
</p>
<p>
To get these values, we first compute the general "distance along
the polygon" matrix:
</p>
<img
class="LaTeX SVG"
src="images/latex/78b8ba1aba2e4c9ad3f7890299c90152.svg"
width="395px"
height="40px"
loading="lazy"
/>
<p>
Where <code>length()</code> 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 <code>i=1</code> and
<code>i=n</code> fall in the [0,1] interval, so we need to scale all
values down by whatever the total length of the polygon is:
</p>
<img
class="LaTeX SVG"
src="images/latex/08f4beaebf83dca594ad125bdca7e436.svg"
width="272px"
height="55px"
loading="lazy"
/>
<p>
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.
</p>
<p>
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 <code>t</code> value", which we'll
square to get rid of any pesky negative signs:
</p>
<img
class="LaTeX SVG"
src="images/latex/7e5d59272621baf942bc722208ce70c2.svg"
width="177px"
height="23px"
loading="lazy"
/>
<p>
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:
</p>
<img
class="LaTeX SVG"
src="images/latex/ab334858d3fa309cc1a5ba535a2ca168.svg"
width="195px"
height="41px"
loading="lazy"
/>
<p>
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 p<sub>i</sub> with the
full <strong>P</strong> coordinate matrix, and we can replace
Bézier(s<sub>i</sub>) with the matrix representation
<strong>T x M x C</strong> we talked about before, which gives us:
</p>
<img
class="LaTeX SVG"
src="images/latex/2d42758fba3370f52191306752c2705c.svg"
width="141px"
height="21px"
loading="lazy"
/>
<p>
In which we can replace the rather cumbersome "squaring" operation
with a more conventional matrix equivalent:
</p>
<img
class="LaTeX SVG"
src="images/latex/5f7fcb86ae1c19612b9fe02e23229e31.svg"
width="225px"
height="21px"
loading="lazy"
/>
<p>
Here, the letter <code>T</code> is used instead of the number 2, to
represent the
<a href="https://en.wikipedia.org/wiki/Transpose"
>matrix transpose</a
>; 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).
</p>
<p>
This leaves one problem: <strong>T</strong> isn't actually the
matrix we want: we don't want symbolic <code>t</code> values, we
want the actual numerical values that we computed for
<strong>S</strong>, so we need to form a new matrix, which we'll
call 𝕋, that makes use of those, and then use that 𝕋 instead of
<strong>T</strong> in our error function:
</p>
<img
class="LaTeX SVG"
src="images/latex/6202d7bd150c852b432d807c40fb1647.svg"
width="201px"
height="96px"
loading="lazy"
/>
<p>
Which, because of the first and last values in <strong>S</strong>,
means:
</p>
<img
class="LaTeX SVG"
src="images/latex/4ffad56e281ee79d0688e93033429f0a.svg"
width="212px"
height="92px"
loading="lazy"
/>
<p>
Now we can properly write out the error function as matrix
operations:
</p>
<img
class="LaTeX SVG"
src="images/latex/8d09f2be2c6db79ee966f170ffc25815.svg"
width="231px"
height="21px"
loading="lazy"
/>
<p>
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:
</p>
<img
class="LaTeX SVG"
src="images/latex/283bc9e8fe59a78d3c74860f62a66ecb.svg"
width="197px"
height="36px"
loading="lazy"
/>
<div className="note">
## Where did this derivative come from?
<p>
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.
</p>
<p>
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
<a href="https://math.stackexchange.com/questions/2825438"
>posted a question</a
>
to
<a href="https://math.stackexchange.com">Math.stackexchange</a>,
and received a answer that goes into way more detail than I had
hoped to receive.
</p>
<p>
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.
</p>
</div>
<p>
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 <strong>C</strong>:
</p>
<img
class="LaTeX SVG"
src="images/latex/d84d1c71a3ce1918f53eaf8f9fe98ac4.svg"
width="168px"
height="27px"
loading="lazy"
/>
<p>
Here, the "to the power negative one" is the notation for the
<a href="https://en.wikipedia.org/wiki/Invertible_matrix"
>matrix inverse</a
>. But that's all we have to do: we're done. Starting with
<strong>P</strong> and inventing some <code>t</code> values based on
the polygon the coordinates in <strong>P</strong> define, we can
compute the corresponding Bézier coordinates <strong>C</strong> that
specify a curve that goes through our points. Or, if it can't go
through them exactly, as near as possible.
</p>
<p>
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 prespecified
coordinates.
</p>
<p>
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 Bezier curve of the appropriate
order. Four points? Cubic Bezier. 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 you hit 10<sup>th</sup> or higher order curves. But it
might not!
</p>
<div className="figure">
<Graphic
title="Fitting a Bézier curve"
setup="{this.setup}"
draw="{this.draw}"
onClick="{this.onClick}"
>
<button
onClick="{this.toggle}"
style="position: absolute; right: 0;"
>
toggle
</button>
<SliderSet ref={ set => (this.sliders=set) }
onChange={this.processTimeUpdate} />
</Graphic>
</div>
<p>
You'll note there is a convenient "toggle" buttons that lets you
toggle between equidistance <code>t</code> values, and distance
ratio along the polygon. Arguably more interesting is that once you
have points to abstract a curve, you also get
<em>direct control</em> 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.
</p>
</section>
<section id="catmullconv">
<h1>Bézier curves and Catmull-Rom curves</h1>
<p>
Taking an excursion to different splines, the other common design
curve is the
<a
href="https://en.wikipedia.org/wiki/Cubic_Hermite_spline#Catmull.E2.80.93Rom_spline"
>Catmull-Rom spline</a
>. Now, a Catmull-Rom spline is a form of
<a href="https://en.wikipedia.org/wiki/Cubic_Hermite_spline"
>cubic Hermite spline</a
>, and as it so happens the cubic Bézier curve is <em>also</em> a
cubic Hermite spline, so maybe... maybe we can convert one into the
other, and back, with some simple substitutions?
</p>
<p>
Unlike Bézier curves, Catmull-Rom splines pass through each point
used to define the curve, except the first and last, which makes
sense if you read the "natural language" description for how a
Catmull-Rom spline works: a Catmull-Rom spline is a curve that, at
each point P<sub>x</sub>, has a tangent along the line P<sub
>x-1</sub
>
to P<sub>x+1</sub>. The curve runs from points P<sub>2</sub> to
P<sub>n-1</sub>, and has a "tension" that determines how fast the
curve passes through each point. The lower the tension, the faster
the curve goes through each point, and the bigger its local tangent
is.
</p>
<p>
I'll be showing the conversion to and from Catmull-Rom curves for
the tension that the Processing language uses for its Catmull-Rom
algorithm.
</p>
<p>
We start with showing the Catmull-Rom matrix form, which looks
similar to the Bézier matrix form, with slightly different values in
the matrix:
</p>
<img
class="LaTeX SVG"
src="images/latex/9215d05705c8e8a7ebd718ae6f690371.svg"
width="409px"
height="75px"
loading="lazy"
/>
<p>
However, there's something funny going on here: the coordinate
column matrix looks weird. The reason is that Catmull-Rom curves are
actually curve segments that are described by two coordinate points,
and two tangents; the curve starts at coordinate V1, and ends at
coordinate V2, with the curve "departing" V1 with a tangent vector
V'1 and "arriving" at V2 with tangent vector V'2.
</p>
<p>
This is not particularly useful if we want to draw Catmull-Rom
curves in the same way we draw Bézier curves, i.e. by providing four
points. However, we can fairly easily go from the former to the
latter, but it's going to require some linear algebra, so if you
just want to know how to convert between the two coordinate systems:
skip the following bit.
</p>
<p>
But... if you want to know <em>why</em> that conversion works, let's
do some maths!
</p>
<div className="note">
<h2>Deriving the conversion formulae</h2>
<p>
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.
</p>
<p>
So, let's start with the first, where we want to satisfy the
following equality:
</p>
<img
class="LaTeX SVG"
src="images/latex/1811b59c5ab9233f08590396e5d03303.svg"
width="187px"
height="83px"
loading="lazy"
/>
<p>
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
<em>T</em> 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.
</p>
<p>Computing <em>T</em> is really more "arranging the numbers":</p>
<img
class="LaTeX SVG"
src="images/latex/8d1a91ddd04d2cc02ec2d9372760d1b1.svg"
width="591px"
height="83px"
loading="lazy"
/>
<p>Thus:</p>
<img
class="LaTeX SVG"
src="images/latex/ffa21dac1f629c75d68dc3f8f8b77b92.svg"
width="143px"
height="81px"
loading="lazy"
/>
<p>
However, we're not <em>quite</em> done, because Catmull-Rom curves
have a parameter called "tension", written as τ ("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 <em>T</em>:
</p>
<img
class="LaTeX SVG"
src="images/latex/f6be211b2d3b4e96553fc2340222264a.svg"
width="285px"
height="84px"
loading="lazy"
/>
<p>
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:
</p>
<img
class="LaTeX SVG"
src="images/latex/cc1e2ff43350c32f0ae9ba9a7652b8fb.svg"
width="409px"
height="75px"
loading="lazy"
/>
<p>
Replace point/tangent vector with the expression for
all-coordinates:
</p>
<img
class="LaTeX SVG"
src="images/latex/f08e34395ce2812276fd70548f805041.svg"
width="549px"
height="81px"
loading="lazy"
/>
<p>and merge the matrices:</p>
<img
class="LaTeX SVG"
src="images/latex/f2b2a16a41d134ce0dfd544ab77ff25e.svg"
width="455px"
height="84px"
loading="lazy"
/>
<p>
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:
</p>
<img
class="LaTeX SVG"
src="images/latex/8f56909fcb62b8eef18b9b9559575c13.svg"
width="353px"
height="73px"
loading="lazy"
/>
<p>
So, if we want to express a Catmull-Rom curve using a Bézier
curve, we'll need to turn this Catmull-Rom bit:
</p>
<img
class="LaTeX SVG"
src="images/latex/b21386f86bef8894f108c5441dad10de.svg"
width="227px"
height="84px"
loading="lazy"
/>
<p>Into something that looks like this:</p>
<img
class="LaTeX SVG"
src="images/latex/78ac9df086ec19147414359369b563fc.svg"
width="167px"
height="73px"
loading="lazy"
/>
<p>
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:
</p>
<img
class="LaTeX SVG"
src="images/latex/a47b072a325812ac4f0ff52c22792588.svg"
width="440px"
height="84px"
loading="lazy"
/>
<p>
Then we remove the coordinate vector from both sides without
affecting the equality:
</p>
<img
class="LaTeX SVG"
src="images/latex/841fb6a2a035c9bcf5a2d46f2a67709b.svg"
width="353px"
height="84px"
loading="lazy"
/>
<p>
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:
</p>
<img
class="LaTeX SVG"
src="images/latex/cbdd46d5e2e1a6202ef46fb03711ebe4.svg"
width="657px"
height="88px"
loading="lazy"
/>
<p>
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:
</p>
<img
class="LaTeX SVG"
src="images/latex/3ea54fe939d076f8db605c5b480e7db0.svg"
width="369px"
height="88px"
loading="lazy"
/>
<p>
And now we're <em>basically</em> done. We just multiply those two
matrices and we know what <em>V</em> is:
</p>
<img
class="LaTeX SVG"
src="images/latex/169fd85a95e4d16fe289a75583017a11.svg"
width="161px"
height="77px"
loading="lazy"
/>
<p>
We now have the final piece of our function puzzle. Let's run
through each step.
</p>
<ol>
<li>Start with the Catmull-Rom function:</li>
</ol>
<img
class="LaTeX SVG"
src="images/latex/cc1e2ff43350c32f0ae9ba9a7652b8fb.svg"
width="409px"
height="75px"
loading="lazy"
/>
<ol start="2">
<li>rewrite to pure coordinate form:</li>
</ol>
<img
class="LaTeX SVG"
src="images/latex/f814bb8d627f9c8f33b347c1cf13d4c7.svg"
width="324px"
height="84px"
loading="lazy"
/>
<ol start="3">
<li>rewrite for "normal" coordinate vector:</li>
</ol>
<img
class="LaTeX SVG"
src="images/latex/5f2750de827497375d9a915f96686885.svg"
width="441px"
height="81px"
loading="lazy"
/>
<ol start="4">
<li>merge the inner matrices:</li>
</ol>
<img
class="LaTeX SVG"
src="images/latex/79e333cd0c569657eea033b04fb5e61b.svg"
width="348px"
height="84px"
loading="lazy"
/>
<ol start="5">
<li>rewrite for Bézier matrix form:</li>
</ol>
<img
class="LaTeX SVG"
src="images/latex/1b8a782f7540503d38067317e4cd00b0.svg"
width="431px"
height="77px"
loading="lazy"
/>
<ol start="6">
<li>
and transform the coordinates so we have a "pure" Bézier
expression:
</li>
</ol>
<img
class="LaTeX SVG"
src="images/latex/e3d30ab368dcead1411532ce3814d3f3.svg"
width="348px"
height="81px"
loading="lazy"
/>
<p>
And we're done: we finally know how to convert these two curves!
</p>
</div>
<p>
If we have a Catmull-Rom curve defined by four coordinates P<sub
>1</sub
>
through P<sub>4</sub>, then we can draw that curve using a Bézier
curve that has the vector:
</p>
<img
class="LaTeX SVG"
src="images/latex/ba31c32eba62f1e3b15066cd5ddda597.svg"
width="249px"
height="85px"
loading="lazy"
/>
<p>
Similarly, if we have a Bézier curve defined by four coordinates
P<sub>1</sub> through P<sub>4</sub>, we can draw that using a
standard tension Catmull-Rom curve with the following coordinate
values:
</p>
<img
class="LaTeX SVG"
src="images/latex/eae7f01976e511ee38b08b6edc8765d2.svg"
width="284px"
height="77px"
loading="lazy"
/>
<p>
or, if your API requires specifying Catmull-Rom curves using "point
+ tangent" form:
</p>
<img
class="LaTeX SVG"
src="images/latex/26363fc09f8cf2d41ea5b4256656bb6d.svg"
width="284px"
height="77px"
loading="lazy"
/>
</section>
<section id="catmullmoulding">
<h1>Creating a Catmull-Rom curve from three points</h1>
<p>
Now, we saw how to fit a Bézier curve to three points, but if
Catmull-Rom curves go through points, why can't we just use those to
do curve fitting, instead?
</p>
<p>
As a matter of fact, we can, but there's a difference between the
kind of curve fitting we did in the previous section, and the kind
of curve fitting that we can do with Catmull-Rom curves. In the
previous section we came up with a single curve that goes through
three points. There was a decent amount of maths and computation
involved, and the end result was three or four coordinates that
described a single curve, depending on whether we were fitting a
quadratic or cubic curve.
</p>
<p>
Using Catmull-Rom curves, we need virtually no computation, but even
though we end up with one Catmull-Rom curve of <i>n</i> points, in
order to draw the equivalent curve using cubic Bézier curves we need
a massive <i>3n-2</i> points (and that's without double-counting
points that are shared by consecutive cubic curves).
</p>
<p>
In the following graphic, on the left we see three points that we
want to draw a Catmull-Rom curve through (which we can move around
freely, by the way), with in the second panel some of the
"interesting" Catmull-Rom information: in black there's the baseline
start--end, which will act as tangent orientation for the curve at
point p2. We also see a virtual point p0 and p4, which are initially
just point p2 reflected over the baseline. However, by using the up
and down cursor key we can offset these points parallel to the
baseline. Why would we want to do this? Because the line p0--p2 acts
as departure tangent at p1, and the line p2--p4 acts as arrival
tangent at p3. Play around with the graphic a bit to get an idea of
what all of that meant:
</p>
<Graphic
title="Catmull-Rom curve fitting"
setup="{this.setup}"
draw="{this.draw}"
onKeyDown="{this.props.onKeyDown}"
/>
<p>
As should be obvious by now, Catmull-Rom curves are great for
"fitting a curvature to some points", but if we want to convert that
curve to Bézier form we're going to end up with a lot of separate
(but visually joined) Bézier curves. Depending on what we want to
do, that'll be either unnecessary work, or exactly what we want:
which it is depends entirely on you.
</p>
</section>
<section id="polybezier">
<h1>Forming poly-Bézier curves</h1>
<p>
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:
</p>
<ol>
<li>
the end point of each section is the starting point of the
following section, and
</li>
<li>the derivatives across that dual point line up.</li>
</ol>
<p>
Unless you want sharp corners, of course. Then you don't even need
2.
</p>
<p>
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:
</p>
<Graphic
title="Unlinked quadratic poly-Bézier"
setup="{this.setupQuadratic}"
draw="{this.draw}"
/>
<Graphic
title="Unlinked cubic poly-Bézier"
setup="{this.setupCubic}"
draw="{this.draw}"
/>
<p>
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:
</p>
<img
class="LaTeX SVG"
src="images/latex/408dd95905a5f001179c4da6051e49c5.svg"
width="124px"
height="17px"
loading="lazy"
/>
<p>
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:
</p>
<img
class="LaTeX SVG"
src="images/latex/8c1b570b3efdfbbc39ddedb4adcaaff6.svg"
width="304px"
height="40px"
loading="lazy"
/>
<p>
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...
</p>
<Graphic
title="Connected quadratic poly-Bézier"
setup="{this.setupQuadratic}"
draw="{this.draw}"
onMouseMove="{this.linkDerivatives}"
/>
<Graphic
title="Connected cubic poly-Bézier"
setup="{this.setupCubic}"
draw="{this.draw}"
onMouseMove="{this.linkDerivatives}"
/>
<p>
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's positions is different depending on whether
it's the reflection of its left or right neighbouring control point:
we can't even form a proper rule-conforming curve! 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.
</p>
<p>So: let's relax the requirement a little.</p>
<p>
We can change the constraint so that we still preserve the
<em>angle</em> 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
<em>vector length</em>. Doing so will give us a much more useful
kind of poly-Bézier curve:
</p>
<Graphic
title="Angularly connected quadratic poly-Bézier"
setup="{this.setupQuadratic}"
draw="{this.draw}"
onMouseMove="{this.linkDirection}"
/>
<Graphic
title="Angularly connected cubic poly-Bézier"
setup="{this.setupCubic}"
draw="{this.draw}"
onMouseMove="{this.linkDirection}"
/>
<p>
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...
</p>
<p>
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.
</p>
<Graphic
title="Standard connected quadratic poly-Bézier"
setup="{this.setupQuadratic}"
draw="{this.draw}"
onMouseDown="{this.bufferPoints}"
onMouseMove="{this.modelCurve}"
/>
<Graphic
title="Standard connected cubic poly-Bézier"
setup="{this.setupCubic}"
draw="{this.draw}"
onMouseDown="{this.bufferPoints}"
onMouseMove="{this.modelCurve}"
/>
<p>
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...
</p>
<p>
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.
</p>
</section>
<section id="shapes">
<h1>Boolean shape operations</h1>
<p>
We can apply the topics covered so far in this primer to effect
boolean shape operations: getting the union, intersection, or
exclusion, between two or more shapes that involve Bézier curves.
For simplicity (well... sort of, more homogeneity), we'll be looking
at poly-Bézier shapes only, but a shape that consists of a mix of
lines and Bézier curves is technically a simplification. (Although
it does mean we need to write a definition for the class of shapes
that mix lines and Bézier curves. Since poly-Bézier curves are a
superset, we'll be using those in the following examples.)
</p>
<p>
The procedure for performing boolean operations consists, broadly,
of four steps:
</p>
<ol>
<li>Find the intersection points between both shapes,</li>
<li>
cut up the shapes into multiple sections between these
intersections,
</li>
<li>
discard any section that isn't part of the desired operation's
resultant shape, and
</li>
<li>link up the remaining sections to form the new shape.</li>
</ol>
<p>
Finding all intersections between two poly-Bézier curves, or any
poly-line-section shape, is similar to the iterative algorithm
discussed in the section on curve/curve intersection. For each
segment in the poly-Bézier curve, we check whether its bounding box
overlaps with any of the segment bounding boxes in the other
poly-Bézier curve. If so, we run normal intersection detection.
</p>
<p>
After finding all intersection points, we split up our poly-Bézier
curves, and make sure to record which of the newly formed
poly-Bézier curves might potentially link up at the points we split
the originals up at. This will let us quickly glue poly-Bézier
curves back together after the next step.
</p>
<p>
Once we have all the new poly-Bézier curves, we run the first step
of the desired boolean operation.
</p>
<ul>
<li>
Union: discard all poly-Bézier curves that lie "inside" our union
of our shapes. E.g. if we want the union of two overlapping
circles, the resulting shape is the outline.
</li>
<li>
Intersection: discard all poly-Bézier curves that lie "outside"
the intersection of the two shapes. E.g. if we want the
intersection of two overlapping circles, the resulting shape is
the tapered ellipse where they overlap.
</li>
<li>
Exclusion: none of the sections are discarded, but we will need to
link the shapes back up in a special way. Flip any section that
would qualify for removal under UNION rules.
</li>
</ul>
<table className="sketch">
<tbody>
<tr>
<td className="labeled-image">
<img src="images/op_base.gif" height="169" />
Two overlapping shapes.
</td>
<td className="labeled-image">
<img src="images/op_union.gif" height="169" />
The unified region.
</td>
<td className="labeled-image">
<img src="images/op_intersection.gif" height="169" />
Their intersection.
</td>
<td className="labeled-image">
<img src="images/op_exclusion.gif" height="169" />
Their exclusion regions.
</td>
</tr>
</tbody>
</table>
<p>
The main complication in the outlined procedure here is determining
how sections qualify in terms of being "inside" and "outside" of our
shapes. For this, we need to be able to perform point-in-shape
detection, for which we'll use a classic algorithm: getting the
"crossing number" by using ray casting, and then testing for
"insidedness" by applying the
<a
href="http://folk.uio.no/bjornw/doc/bifrost-ref/bifrost-ref-12.html"
>even-odd rule</a
>: For any point and any shape, we can cast a ray from our point, to
some point that we know lies outside of the shape (such as a corner
of our drawing surface). We then count how many times that line
crosses our shape (remember that we can perform line/curve
intersection detection quite easily). If the number of times it
crosses the shape's outline is even, the point did not actually lie
inside our shape. If the number of intersections is odd, our point
did lie inside out shape. With that knowledge, we can decide whether
to treat a section that such a point lies on "needs removal" (under
union rules), "needs preserving" (under intersection rules), or
"needs flipping" (under exclusion rules).
</p>
<p>
These operations are expensive, and implementing your own code for
this is generally a bad idea if there is already a geometry package
available for your language of choice. In this case, for JavaScript
the most excellent <a href="http://paperjs.org">Paper.js</a> already
comes with all the code in place to perform efficient boolean shape
operations, so rather that implement an inferior version here, I can
strongly recommend the Paper.js library if you intend to do any
boolean shape work.
</p>
<p>
The following graphic shows Paper.js doing its thing for two shapes:
one static, and one that is linked to your mouse pointer. If you
move the mouse around, you'll see how the shape intersections are
resolved. The base shapes are outlined in blue, and the boolean
result is coloured red.
</p>
<Graphic
title="Boolean shape operations with Paper.js"
paperjs="{true}"
setup="{this.setup}"
draw="{this.draw}"
onMouseMove="{this.onMouseMove}"
>
<button onclick="() => this.setMode(mode)">mode</button>
</Graphic>
</section>
<section id="projections">
<h1>Projecting a point onto a Bézier curve</h1>
<p>
Say we have a Bézier curve and some point, not on the curve, of
which we want to know which <code>t</code> value on the curve gives
us an on-curve point closest to our off-curve point. Or: say we want
to find the projection of a random point onto a curve. How do we do
that?
</p>
<p>
If the Bézier curve is of low enough order, we might be able to
<a
href="https://web.archive.org/web/20140713004709/http://jazzros.blogspot.com/2011/03/projecting-point-on-bezier-curve.html"
>work out the maths for how to do this</a
>, and get a perfect <code>t</code> 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
<code>t</code> value using a
<a href="https://en.wikipedia.org/wiki/Binary_search_algorithm"
>binary search</a
>. First, we do a coarse distance-check based on
<code>t</code> values associated with the curve's "to draw"
coordinates (using a lookup table, or LUT). This is pretty fast.
Then we run this algorithm:
</p>
<ol>
<li>
with the <code>t</code> value we found, start with some small
interval around <code>t</code> (1/length_of_LUT on either side is
a reasonable start),
</li>
<li>
if the distance to <code>t ± interval/2</code> is larger than the
distance to <code>t</code>, try again with the interval reduced to
half its original length.
</li>
<li>
if the distance to <code>t ± interval/2</code> is smaller than the
distance to <code>t</code>, replace <code>t</code> with the
smaller-distance value.
</li>
<li>
after reducing the interval, or changing <code>t</code>, go back
to step 1.
</li>
</ol>
<p>
We keep repeating this process until the interval is small enough to
claim the difference in precision found is irrelevant for the
purpose we're trying to find <code>t</code> for. In this case, I'm
arbitrarily fixing it at 0.0001.
</p>
<p>
The following graphic demonstrates the result of this procedure.
Simply move the cursor around, and if it does not lie on top of the
curve, you will see a line that projects the cursor onto the curve
based on an iteratively found "ideal" <code>t</code> value.
</p>
<Graphic
title="Projecting a point onto a Bézier curve"
setup="{this.setup}"
draw="{this.draw}"
onMouseMove="{this.onMouseMove}"
/>
</section>
<section id="offsetting">
<h1>Curve offsetting</h1>
<p>
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.
</p>
<p>
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?
</p>
<p>
Bottom line: <strong>you can't</strong>. 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.
</p>
<div className="note">
<h3>"What do you mean, 'you can't'? Prove it."</h3>
<p>
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
<strong>has</strong> to be a polynomial. Let's look at why this
is:
</p>
<p>
From a mathematical point of view, an offset curve
<code>O(t)</code> is a curve such that, given our original curve
<code>B(t)</code>, any point on <code>O(t)</code> is a fixed
distance <code>d</code> away from coordinate <code>B(t)</code>. So
let's math that:
</p>
<img
class="LaTeX SVG"
src="images/latex/1d4be24e5896dce3c16c8e71f9cc8881.svg"
width="108px"
height="16px"
loading="lazy"
/>
<p>
However, we're working in 2D, and <code>d</code> is a single
value, so we want to turn it into a vector. If we want a point
distance <code>d</code> "away" from the curve
<code>B(t)</code> then what we really mean is that we want a point
at <code>d</code> times the "normal vector" from point
<code>B(t)</code>, where the "normal" is a vector that runs
perpendicular ("at a right angle") to the tangent at
<code>B(t)</code>. Easy enough:
</p>
<img
class="LaTeX SVG"
src="images/latex/5bfee4f2ae27304475673d0596e42f9a.svg"
width="151px"
height="16px"
loading="lazy"
/>
<p>
Now this still isn't very useful unless we know what the formula
for <code>N(t)</code> is, so let's find out.
<code>N(t)</code> runs perpendicular to the original curve
tangent, and we know that the tangent is simply
<code>B'(t)</code>, so we could just rotate that 90 degrees and be
done with it. However, we need to ensure that
<code>N(t)</code> has the same magnitude for every <code>t</code>,
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 <code>N(t)</code> always has length 1, which we can
achieve by dividing <code>B'(t)</code> by its magnitude:
</p>
<img
class="LaTeX SVG"
src="images/latex/fa6c243de2aa78b7451e0086848dfdfc.svg"
width="120px"
height="40px"
loading="lazy"
/>
<p>
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 <code>a</code> to end <code>b</code>,
we must use the formula we saw earlier. Noting that "length" is
usually denoted with double vertical bars:
</p>
<img
class="LaTeX SVG"
src="images/latex/b262e50c085815421d94e120fc17f1c8.svg"
width="169px"
height="36px"
loading="lazy"
/>
<p>
So if we want the length of the tangent, we plug in
<code>B'(t)</code>, with <code>t = 0</code> as start and
<code>t = 1</code> as end:
</p>
<img
class="LaTeX SVG"
src="images/latex/1d586b939b44ff9bdb42562a12ac2779.svg"
width="209px"
height="36px"
loading="lazy"
/>
<p>
And that's where things go wrong. It doesn't even really matter
what the second derivative for <code>B(t)</code> is, that square
root is screwing everything up, because it turns our nice
polynomials into things that are no longer polynomials.
</p>
<p>
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 <strong>base</strong> 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.
</p>
<p>
So, since the tangent length isn't a polynomial, the normalised
tangent won't be a polynomial either, which means
<code>N(t)</code> won't be a polynomial, which means that
<code>d</code> times <code>N(t)</code> won't be a polynomial,
which means that, ultimately, <code>O(t)</code> won't be a
polynomial, which means that even if we can determine the function
for <code>O(t)</code> just fine (and that's far from trivial!), it
simply cannot be represented as a Bézier curve.
</p>
<p>
And that's one reason why Bézier curves are tricky: there are
actually a <em>lot</em> 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.
</p>
</div>
<p>
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 <code>t=0.5</code> 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).
</p>
<p>
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
<code>t=0.5</code> 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.
</p>
<p>
The following graphics show off curve offsetting, and you can use
your up and down arrow keys 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.
</p>
<Graphic
title="Offsetting a quadratic Bézier curve"
setup="{this.setupQuadratic}"
draw="{this.draw}"
onKeyDown="{this.props.onKeyDown}"
/>
<Graphic
title="Offsetting a cubic Bézier curve"
setup="{this.setupCubic}"
draw="{this.draw}"
onKeyDown="{this.props.onKeyDown}"
/>
<p>
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.
</p>
</section>
<section id="graduatedoffset">
<h1>Graduated curve offsetting</h1>
<p>
What if we want to do graduated offsetting, starting at some
distance <code>s</code> but ending at some other distance
<code>e</code>? 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.
</p>
<p>
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
<code>S</code> (and if the current sub-curve is the first one,
<code>S</code> is zero), and we call the full length of our original
curve <code>L</code>, then we get the following graduation values:
</p>
<ul>
<li>
start: map <code>S</code> from interval (<code>0,L</code>) to
interval <code>(s,e)</code>
</li>
<li>
c1: <code>map(&lt;strong&gt;S+d1&lt;/strong&gt;, 0,L, s,e)</code>,
d1 = distance along curve to projection of c1
</li>
<li>
c2: <code>map(&lt;strong&gt;S+d2&lt;/strong&gt;, 0,L, s,e)</code>,
d2 = distance along curve to projection of c2
</li>
<li>...</li>
<li>
end:
<code
>map(&lt;strong&gt;S+length(subcurve)&lt;/strong&gt;, 0,L,
s,e)</code
>
</li>
</ul>
<p>
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):
</p>
<Graphic
title="Offsetting a quadratic Bézier curve"
setup="{this.setupQuadratic}"
draw="{this.draw}"
onKeyDown="{this.props.onKeyDown}"
/>
<Graphic
title="Offsetting a cubic Bézier curve"
setup="{this.setupCubic}"
draw="{this.draw}"
onKeyDown="{this.props.onKeyDown}"
/>
</section>
<section id="circles">
<h1>Circles and quadratic Bézier curves</h1>
<p>
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?
</p>
<p>You approximate.</p>
<p>
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.
</p>
<p>
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.
</p>
<p>
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:
</p>
<Graphic
title="Quadratic Bézier arc approximation"
setup="{this.setup}"
draw="{this.draw}"
onMouseMove="{this.onMouseMove}"
/>
<p>
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.
</p>
<p>
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 <em>φ</em>:
</p>
<img
class="LaTeX SVG"
src="images/latex/8374c4190d6213b0ac0621481afaa754.svg"
width="175px"
height="40px"
loading="lazy"
/>
<p>
What we want to find is the intersection of the tangents, so we want
a point C such that:
</p>
<img
class="LaTeX SVG"
src="images/latex/a127f926eced2751a09c54bf7c361b4a.svg"
width="284px"
height="40px"
loading="lazy"
/>
<p>
i.e. we want a point that lies on the vertical line through S (at
some distance <em>a</em> from S) and also lies on the tangent line
through E (at some distance <em>b</em> from E). Solving this gives
us:
</p>
<img
class="LaTeX SVG"
src="images/latex/b5d864e9ed0c44c56d454fbaa4218d5e.svg"
width="219px"
height="40px"
loading="lazy"
/>
<p>First we solve for <em>b</em>:</p>
<img
class="LaTeX SVG"
src="images/latex/9e4d886c372f916f6511c41245ceee39.svg"
width="560px"
height="17px"
loading="lazy"
/>
<p>which yields:</p>
<img
class="LaTeX SVG"
src="images/latex/c22f6d343ee0cce7bff6a617c946ca17.svg"
width="101px"
height="39px"
loading="lazy"
/>
<p>which we can then substitute in the expression for <em>a</em>:</p>
<img
class="LaTeX SVG"
src="images/latex/ef3ab62bb896019c6157c85aae5d1ed3.svg"
width="231px"
height="195px"
loading="lazy"
/>
<p>
A quick check shows that plugging these values for <em>a</em> and
<em>b</em> into the expressions for C<sub>x</sub> and C<sub>y</sub>
give the same x/y coordinates for both "<em>a</em> away from A" and
"<em>b</em> 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
<em>t=0.5</em> (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:
</p>
<img
class="LaTeX SVG"
src="images/latex/fe32474b4616ee9478e1308308f1b6bf.svg"
width="188px"
height="32px"
loading="lazy"
/>
<p>
We compute T, observing that if <em>t=0.5</em>, the polynomial
values (1-t)², 2(1-t)t, and t² are 0.25, 0.5, and 0.25 respectively:
</p>
<img
class="LaTeX SVG"
src="images/latex/e1059e611aa1e51db41f9ce0b4ebb95a.svg"
width="252px"
height="35px"
loading="lazy"
/>
<p>Which, worked out for the x and y components, gives:</p>
<img
class="LaTeX SVG"
src="images/latex/7754bc3c96ae3c90162fec3bd46bedff.svg"
width="408px"
height="77px"
loading="lazy"
/>
<p>
And the distance between these two is the standard Euclidean
distance:
</p>
<img
class="LaTeX SVG"
src="images/latex/adbd056f4b8fcd05b1d4f2fce27d7657.svg"
width="399px"
height="153px"
loading="lazy"
/>
<p>
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?
</p>
<table>
<tbody>
<tr>
<td>
<img src="images/arc-q-pi.gif" height="190" />
plotted for 0 ≤ φ ≤ π:
</td>
<td>
<img src="images/arc-q-pi2.gif" height="187" />
plotted for 0 ≤ φ ≤ ½π:
</td>
<td>
{ this.props.showhref ?
"http://www.wolframalpha.com/input/?i=plot+sqrt%28%281%2F4+*+%28sin%28x%29+%2B+2tan%28x%2F2%29%29+-+sin%28x%2F2%29%29%5E2+%2B+%282sin%5E4%28x%2F4%29%29%5E2%29+for+0+%3C%3D+x+%3C%3D+pi%2F4"
: null }
<img src="images/arc-q-pi4.gif" height="174" />
plotted for 0 ≤ φ ≤ ¼π:
</td>
</tr>
</tbody>
</table>
<p>
We now see why the eighth circle arc looks decent, but the quarter
circle arc doesn't: an error of roughly 0.06 at <em>t=0.5</em> 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!
</p>
<p>
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:
</p>
<img
class="LaTeX SVG"
src="images/latex/df87674db0f31fc3944aaeb6b890e196.svg"
width="247px"
height="53px"
loading="lazy"
/>
<p>
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).
</p>
<p>
The bottom line?
<strong>Quadratic curves are kind of lousy</strong> 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.
</p>
</section>
<section id="circles_cubic">
<h1>Circles and cubic Bézier curves</h1>
<p>
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?
</p>
<p>
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.
</p>
<p>
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.
</p>
<p>
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:
</p>
<Graphic
title="Cubic Bézier arc approximation"
setup="{this.setup}"
draw="{this.draw}"
onMouseMove="{this.onMouseMove}"
/>
<p>
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.
</p>
<p>
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!
</p>
<p>
So, maths time again: how okay is "okay"? Let's apply some more
maths to find out.
</p>
<p>
Unlike for the quadratic curve, we can't use <i>t=0.5</i> as our
reference point because by its very nature it's one of the three
points that are actually guaranteed to lie on the circular curve.
Instead, we need a different <i>t</i> value. If we run some analysis
on the curve we find that the actual <i>t</i> value at which the
curve is furthest from what it should be is 0.211325 (rounded), but
we don't know "why", since finding this value involves root-finding,
and is nearly impossible to do symbolically without pages and pages
of math just to express one of the possible solutions.
</p>
<p>
So instead of walking you through the derivation for that value,
let's simply take that <i>t</i> value and see what the error is for
circular arcs with an angle ranging from 0 to 2π:
</p>
<table>
<tbody>
<tr>
<td>
<img src="images/arc-c-2pi.gif" height="187" />
plotted for 0 ≤ φ ≤ 2π:
</td>
<td>
<img src="images/arc-c-pi.gif" height="187" />
plotted for 0 ≤ φ ≤ π:
</td>
<td>
<img src="images/arc-c-pi2.gif" height="187" />
plotted for 0 ≤ φ ≤ ½π:
</td>
</tr>
</tbody>
</table>
<p>
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!
</p>
<p>
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.
</p>
<p>
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:
</p>
<img
class="LaTeX SVG"
src="images/latex/a4f0dafbfe80c88723c3cc22277a9682.svg"
width="175px"
height="40px"
loading="lazy"
/>
<p>
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:
</p>
<img
class="LaTeX SVG"
src="images/latex/dfb83eec053c30e0a41b0a52aba24cd4.svg"
width="113px"
height="40px"
loading="lazy"
/>
<p>where "a" is some scaling factor, and:</p>
<img
class="LaTeX SVG"
src="images/latex/e75a848f5f8aead495e35175e2955e06.svg"
width="163px"
height="40px"
loading="lazy"
/>
<p>where "b" is also some scaling factor.</p>
<p>
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!
</p>
<div className="note">
<h2>Let's do this thing.</h2>
<p>
Unlike for the quadratic case, we need some more information in
order to compute <i>a</i> and <i>b</i>, since they're no longer
dependent variables. First, we observe that the curve is
symmetrical, so whatever values we end up finding for C<sub
>1</sub
>
will apply to C<sub>2</sub> as well (rotated along its tangent),
so we'll focus on finding the location of C<sub>1</sub> only. So
here's where we do something that you might not expect: we're
going to ignore for a moment, because we're going to have a much
easier time if we just solve this problem with geometry first,
then move to calculus to solve a much simpler problem.
</p>
<p>
If we look at the triangle that is formed between our starting
point, or initial guess C<sub>1</sub> and our real C<sub>1</sub>,
there's something funny going on: if we treat the line
{start,guess} as our opposite side, the line {guess,real} as our
adjacent side, with {start,real} our hypotenuse, then the angle
for the corner hypotenuse/adjacent is half that of the arc we're
covering. Try it: if you place the end point at a quarter circle
(pi/2, or 90 degrees), the angle in our triangle is half a quarter
(pi/4, or 45 degrees). With that knowledge, and a knowledge of
what the length of any of our lines segments are (as a function),
we can determine where our control points are, and thus have
everything we need to find the error distance function. Of the
three lines, the one we can easiest determine is {start,guess}, so
let's find out what the guessed control point is. Again
geometrically, because we have the benefit of an on-curve
<i>t=0.5</i> value.
</p>
<p>
The distance from our guessed point to the start point is exactly
the same as the projection distance we looked at earlier. Using
<i>t=0.5</i> as our point "B" in the "A,B,C" projection, then we
know the length of the line segment {C,A}, since it's d<sub
>1</sub
>
= {A,B} + d<sub>2</sub> = {B,C}:
</p>
<img
class="LaTeX SVG"
src="images/latex/3189cac1ddac07c1487e1e51740ecc88.svg"
width="397px"
height="40px"
loading="lazy"
/>
<p>
So that just leaves us to find the distance from <i>t=0.5</i> to
the baseline for an arbitrary angle φ, which is the distance from
the centre of the circle to our <i>t=0.5</i> point, minus the
distance from the centre to the line that runs from start point to
end point. The first is the same as the point P we found for the
quadratic curve:
</p>
<img
class="LaTeX SVG"
src="images/latex/fe32474b4616ee9478e1308308f1b6bf.svg"
width="188px"
height="32px"
loading="lazy"
/>
<p>
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:
</p>
<img
class="LaTeX SVG"
src="images/latex/0364731626a530c8a9b30f424ada53c5.svg"
width="261px"
height="67px"
loading="lazy"
/>
<p>
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}:
</p>
<img
class="LaTeX SVG"
src="images/latex/ee08d86b7497c7ab042ee899bf15d453.svg"
width="397px"
height="48px"
loading="lazy"
/>
<img
class="LaTeX SVG"
src="images/latex/195790bae7de813aec342ea82b5d8781.svg"
width="211px"
height="47px"
loading="lazy"
/>
<p>
Which means we can now determine the distance {start,guessed},
which is the same as the distance {C,A}, and use that to determine
the vertical distance from our start point to our C<sub>1</sub>:
</p>
<img
class="LaTeX SVG"
src="images/latex/178a838274748439778e2a29f5a27d0b.svg"
width="252px"
height="56px"
loading="lazy"
/>
<p>
And after this tedious detour to find the coordinate for
C<sub>1</sub>, we can find C<sub>2</sub> fairly simply, since it's
lies at distance -C<sub>1y</sub> along the end point's tangent:
</p>
<img
class="LaTeX SVG"
src="images/latex/49dbf244d50c787a4ab18694488d9b69.svg"
width="524px"
height="79px"
loading="lazy"
/>
<p>
And that's it, we have all four points now for an approximation of
an arbitrary circular arc with angle φ.
</p>
</div>
<p>
So, to recap, given an angle φ, the new control coordinates are:
</p>
<img
class="LaTeX SVG"
src="images/latex/e2258660a796dcd6189a6f5e14326dad.svg"
width="205px"
height="40px"
loading="lazy"
/>
<p>and</p>
<img
class="LaTeX SVG"
src="images/latex/acbc5efb06bc34571ccc0322376e0b9b.svg"
width="321px"
height="40px"
loading="lazy"
/>
<p>
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:
</p>
<img
class="LaTeX SVG"
src="images/latex/877f9c217c51c0087be751a7580ed459.svg"
width="412px"
height="33px"
loading="lazy"
/>
<p>
Which, in decimal values, rounded to six significant digits, is:
</p>
<img
class="LaTeX SVG"
src="images/latex/05d36e051a38905dcb81e65db8261f24.svg"
width="412px"
height="16px"
loading="lazy"
/>
<p>
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:
</p>
<Graphic
title="Cubic Bézier circle approximation"
draw="{this.drawCircle}"
static="{true}"
/>
</section>
<section id="arcapproximation">
<h1>Approximating Bézier curves with circular arcs</h1>
<p>
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.
</p>
<p>
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.
</p>
<p>
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.
</p>
<p>
So: step 1, how do we find a circle through three points? That part
is actually really simple. You may remember (if you ever learned
it!) that a line between two points on a circle is called a
<a href="https://en.wikipedia.org/wiki/Chord_%28geometry%29"
>chord</a
>, and 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.
</p>
<p>
So: if we have have three points, we have three (different) chords,
and consequently, three (different) lines that go from those chords
through the center of the circle. So we find the centers of the
chords, find the perpendicular lines, find the intersection of those
lines, and thus find the center of the circle.
</p>
<p>
The following graphic shows this procedure with a different colour
for each chord and its associated perpendicular through the center.
You can move the points around as much as you like, those lines will
always meet!
</p>
<Graphic
title="Finding a circle through three points"
setup="{this.setupCircle}"
draw="{this.drawCircle}"
/>
<p>
So, with the procedure on how to find a circle through three points,
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, over the
remaining point, to the end point.
</p>
<p>
So how can we convert a Bézier curve into a (sequence of) circular
arc(s)?
</p>
<ul>
<li>Start at <em>t=0</em></li>
<li>
Pick two points further down the curve at some value
<em>m = t + n</em> and <em>e = t + 2n</em>
</li>
<li>Find the arc that these points define</li>
<li>
Determine how close the found arc is to the curve:
<ul>
<li>
Pick two additional points <em>e1 = t + n/2</em> and
<em>e2 = t + n + n/2</em>.
</li>
<li>
These points, if the arc is a good approximation of the curve
interval chosen, should lie <em>on</em> the 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.
</li>
<li>
For point points, determine the (absolute) error between the
radius of the circle, and the <em>actual</em> distance from
the center of the circle to the point on the curve.
</li>
<li>
If this error is too high, we consider the arc bad, and try a
smaller interval.
</li>
</ul>
</li>
</ul>
<p>
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 <em>t=0.5</em>, and then we start performing a
<a href="https://en.wikipedia.org/wiki/Binary_search_algorithm"
>Binary Search</a
>.
</p>
<ol>
<li>We start with {0, 0.5, 1}</li>
<li>
That'll fail, so we retry with the interval halved: {0, 0.25, 0.5}
<ul>
<li>
If that arc's good, we move back up by half distance: {0,
0.375, 0.75}.
</li>
<li>
However, if the arc was still bad, we move <em>down</em> by
half the distance: {0, 0.125, 0.25}.
</li>
</ul>
</li>
<li>
We keep doing this over and over until we have two arcs found 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
former.
</li>
</ol>
<p>
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
<em>combined</em> 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.
</p>
<Graphic
title="Arc approximation of a Bézier curve"
setup="{this.setupCubic}"
draw="{this.drawSingleArc}"
onKeyDown="{this.props.onKeyDown}"
/>
<p>
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
<em>t=1</em>, 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:
</p>
<Graphic
title="Arc approximation of a Bézier curve"
setup="{this.setupCubic}"
draw="{this.drawArcs}"
onKeyDown="{this.props.onKeyDown}"
/>
<p>
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!
</p>
</section>
<section id="bsplines">
<h1>B-Splines</h1>
<p>
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.
</p>
<p>
First off: B-Splines are
<a href="https://en.wikipedia.org/wiki/Piecewise"
>piecewise polynomial interpolation curves</a
>, 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.
</p>
<p>
What do they look like? They look like this! .. okay that's an empty
graph, but simply click to place some point, with the stipulation
that you need at least four point to see any curve. More than four
points simply draws a longer B-Spline curve:
</p>
<BSplineGraphic sketch="{this.basicSketch}" />
<p>
The important part to notice here is that we are
<strong>not</strong> 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
<em>each possible curve involving four consecutive points</em>, 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.
</p>
<p>Consider the difference to be this:</p>
<ul>
<li>
for Bézier curves, the curve is defined as an interpolation of
points, but:
</li>
<li>
for B-Splines, the curve is defined as an interpolation of
<em>curves</em>.
</li>
</ul>
<p>
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.
</p>
<h2>How to compute a B-Spline curve: some maths</h2>
<p>
Given a B-Spline of degree <code>d</code> and thus order
<code>k=d+1</code> (so a quadratic B-Spline is degree 2 and order 3,
a cubic B-Spline is degree 3 and order 4, etc) and
<code>n</code> control points
<code>P&lt;sub&gt;0&lt;/sub&gt;</code> through
<code>P&lt;sub&gt;n-1&lt;/sub&gt;</code>, we can compute a point on
the curve for some value <code>t</code> 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:
</p>
<img
class="LaTeX SVG"
src="images/latex/0f3451c711c0fe5d0b018aa4aa77d855.svg"
width="169px"
height="41px"
loading="lazy"
/>
<p>
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 <em>N(...)</em> function, subscripted with an obvious
parameter <code>i</code>, which comes from our summation, and some
magical parameter <code>k</code>. So we need to know two things: 1.
what does N(t) do, and 2. what is that <code>k</code>? Let's cover
both, in reverse order.
</p>
<p>
The parameter <code>k</code> represents the "knot interval" over
which a section of curve is defined. As we learned earlier, a
B-Spline curve is itself an interpoliation 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 <code>d</code> B-Spline with <code>n</code> control
point gives us <code>d + n + 1</code> knots, defining
<code>d + n</code> intervals along the curve, and it is these
intervals that the above <code>k</code> subscript to the N()
function applies to.
</p>
<p>Then the N() function itself. What does it look like?</p>
<img
class="LaTeX SVG"
src="images/latex/cf45d1ea00d4866abc8a058b130299b4.svg"
width="559px"
height="43px"
loading="lazy"
/>
<p>
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
<code>i</code> goes up, and <code>k</code> 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
<code>i</code>/<code>k</code> values:
</p>
<img
class="LaTeX SVG"
src="images/latex/adac18ea69cc58e01c8d5e15498e4aa6.svg"
width="240px"
height="40px"
loading="lazy"
/>
<p>
And this function finally has a straight up evaluation: if a
<code>t</code> value lies within a knot-specific interval once we
reach a <code>k=1</code> value, it "counts", otherwise it doesn't.
We did cheat a little, though, because for all these values we need
to scale our <code>t</code> value first, so that it lies in the
interval bounded by <code>knots[d]</code> and <code>knots[n]</code>,
which are the start point and end point where curvature is
controlled by exactly <code>order</code> 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 <code>t</code> from [the
interval 0,1] to the interval [4,8], and then use that value in the
functions above, instead.
</p>
<h2>Can we simplify that?</h2>
<p>We can, yes.</p>
<p>
People far smarter than us have looked at this work, and two in
particular —
<a href="http://www.npl.co.uk/people/maurice-cox">Maurice Cox</a>
and
<a href="https://en.wikipedia.org/wiki/Carl_R._de_Boor"
>Carl de Boor</a
>
— came to a mathematically pleasing solution: to compute a point
P(t), we can compute this point by evaluating <em>d(t)</em> on a
curve section between knots <em>i</em> and <em>i+1</em>:
</p>
<img
class="LaTeX SVG"
src="images/latex/763838ea6f9e6c6aa63ea5f9c6d9542f.svg"
width="281px"
height="21px"
loading="lazy"
/>
<p>
This is another recursive function, with <em>k</em> values
decreasing from the curve order to 1, and the value
<em>α</em> (alpha) defined by:
</p>
<img
class="LaTeX SVG"
src="images/latex/892209dad8fd1f839470dd061e870913.svg"
width="255px"
height="39px"
loading="lazy"
/>
<p>
That looks complicated, but it's not. Computing alpha is just a
fraction involving known, plain numbers and once we have our alpha
value, computing (1-alpha) is literally just "computing one minus
alpha". Computing this d() function is thus simply a matter of
"computing simple arithmetics but with recursion", which might be
computationally expensive because we're doing "a lot of" steps, but
is also computationally cheap because each step only involves very
simple maths. Of course as before the recursion has to stop:
</p>
<img
class="LaTeX SVG"
src="images/latex/4c8f9814c50c708757eeb5a68afabb7f.svg"
width="368px"
height="40px"
loading="lazy"
/>
<p>
So, we see two stopping conditions: either <code>i</code> becomes 0,
in which case d() is zero, or <code>k</code> becomes zero, in which
case we get the same "either 1 or 0" that we saw in the N() function
above.
</p>
<p>
Thanks to Cox and de Boor, we can compute points on a B-Spline
pretty easily: we just need to compute a triangle of interconnected
values. For instance, d() for i=3, k=3 yields the following
triangle:
</p>
<img
class="LaTeX SVG"
src="images/latex/7962d6fea86da6f53a7269fba30f0138.svg"
width="417px"
height="231px"
loading="lazy"
/>
<p>
That is, we compute d(3,3) as a mixture of d(2,3) and d(2,2): d(3,3)
= a(3,3) x d(2,3) + (1-a(3,3)) x d(2,2)... and we simply keep
expanding our triangle until we reach the terminating function
parameters. Done deal!
</p>
<p>
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
<code>d(..., k)</code> 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: that's
pretty essential!
</p>
<p>
If we run this computation "down", starting at d(3,3), then without
special code in place we would be computing quite a few terms
multiple times at each step. On the other hand, we can also start
with that last "column", we can generate the terminating d() values
first, then compute the a() constants, perform our multiplications,
generate the previous step's d() values, compute their a()
constants, do the multiplications, etc. until we end up all the way
back at the top. If we run our computation this way, we don't need
any explicit caching, we can just "recycle" the list of numbers we
start with and simply update them as we move up the triangle. So,
let's implement that!
</p>
<h2>
Cool, cool... but I don't know what to do with that information
</h2>
<p>
I know, this is pretty mathy, so let's have a look at what happens
when we change parameters here. We can't change the maths for the
interpolation functions, so that gives us only one way to control
what happens here: the knot vector itself. As such, let's look at
the graph that shows the interpolation functions for a cubic
B-Spline with seven points with a uniform knot vector (so we see
seven identical functions), representing how much each point
(represented by one function each) influences the total curvature,
given our knot values. And, because exploration is the key to
discovery, let's make the knot vector a thing we can actually
manipulate. Normally a proper knot vector has a constraint that any
value is strictly equal to, or larger than the previous ones, but
screw it this is programming, let's ignore that hard restriction and
just mess with the knots however we like.
</p>
<div className="two-column">
<KnotController ref="interpolation-graph" />
<BSplineGraphic sketch={this.interpolationGraph} controller={(owner,
knots) => this.bindKnots(owner, knots, "interpolation-graph")}/>
</div>
<p>
Changing the values in the knot vector changes how much each point
influences the total curvature (with some clever knot value
manipulation, we can even make the influence of certain points
disappear entirely!), so we can see that while the control points
define the hull inside of which we're going to be drawing a curve,
it is actually the knot vector that determines the actual
<em>shape</em> of the curve inside that hull.
</p>
<p>
After reading the rest of this section you may want to come back
here to try some specific knot vectors, and see if the resulting
interpolation landscape makes sense given what you will now think
should happen!
</p>
<h2>Running the computation</h2>
<p>
Unlike the de Casteljau algorithm, where the <code>t</code> 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
<a
href="http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/de-Boor.html"
>this Michigan Tech</a
>
page, but an easier to read version is implemented by
<a
href="https://github.com/thibauts/b-spline/blob/master/index.js#L59-L71"
>b-spline.js</a
>, so we'll look at its code.
</p>
<p>
Given an input value <code>t</code>, 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
<code>s</code> that this mapped <code>t</code> value lies on:
</p>
<pre><code>for(s=domain[0]; s &lt; domain[1]; s++) {
if(knots[s] &lt;= t && t &lt;= knots[s+1]) break;
}</code></pre>
<p>
after running this code, <code>s</code> 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):
</p>
<pre><code>let v = copy of control points
for(let L = 1; L &lt;= order; L++) {
for(let i=s; i &gt; s + L - order; i--) {
let numerator = t - knots[i]
let denominator = knots[i - L + order] - knots[i]
let alpha = numerator / denominator
let v[i] = alpha * v[i] + (1-alpha) * v[i-1]
}
}</code></pre>
<p>
(A nice bit of behaviour in this code is that we work the
interpolation "backwards", starting at <code>i=s</code> at each
level of the interpolation, and we stop when
<code>i = s - order + level</code>, so we always end up with a value
for <code>i</code> such that those <code>v[i-1]</code> don't try to
use an array index that doesn't exist)
</p>
<h2>Open vs. closed paths</h2>
<p>
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 <em>the same point</em>. However, because B-Splines
are an interpolation of curves, not just point, we can't simply make
the first and last point the same, we need to link a few point
point: for an order <code>d</code> B-Spline, we need to make the
last <code>d</code> point the same as the first
<code>d</code> points. And the easiest way to do this is to simply
append <code>points.splice(0,d)</code> to <code>points</code>. Done!
</p>
<p>
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 <code>points[0]</code> and <code>points[n-k]</code> etc. are the
same coordinate, and manipulating one will equally manipulate the
other, but programming generally makes this really easy by storing
references to coordinates (or other linked values such as coordinate
weights, discussed in the NURBS section) rather than separate
coordinate objects.
</p>
<h2>Manipulating the curve through the knot vector</h2>
<p>
The most important thing to understand when it comes to B-Splines is
that they work <em>because</em> 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
<em>values</em> 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:
</p>
<ol>
<li>
we can use a uniform knot vector, with equally spaced intervals,
</li>
<li>
we can use a non-uniform knot vector, without enforcing equally
spaced intervals,
</li>
<li>
we can collapse sequential knots to the same value, locally
lowering curve complexity using "null" intervals, and
</li>
<li>
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.
</li>
</ol>
<h3>Uniform B-Splines</h3>
<p>
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
<em>the same intervals</em>, or even [0,2,3,...,18,20,22], which
also defines <em>the same intervals</em>, just scaled by a constant
factor, which becomes normalised during interpolation and so does
not contribute to the curvature.
</p>
<div className="two-column">
<KnotController ref="uniform-spline" />
<BSplineGraphic sketch={this.uniformBSpline} controller={(owner,
knots) => this.bindKnots(owner, knots, "uniform-spline")}/>
</div>
<p>
This is an important point: the intervals that the knot vector
defines are <em>relative</em> 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.
</p>
<p>
The problem with uniform knot vectors is that, as we need
<code>order</code> 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...
</p>
<h3>Reducing local curve complexity by collapsing intervals</h3>
<p>
By collapsing knot intervals by making two or more consecutive knots
have the same value, we can 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
<code>order</code> knots creates a situation where all continuity is
lost and the curve "kinks".
</p>
<div className="two-column">
<KnotController ref="center-cut-bspline" />
<BSplineGraphic sketch={this.centerCutBSpline} controller={(owner,
knots) => this.bindKnots(owner, knots, "center-cut-bspline")}/>
</div>
<h3>Open-Uniform B-Splines</h3>
<p>
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:
</p>
<p>
For any curve of degree <code>D</code> with control points
<code>N</code>, we can define a knot vector of length
<code>N+D+1</code> in which the values <code>0 ... D+1</code> are
the same, the values <code>D+1 ... N+1</code> follow the "uniform"
pattern, and the values <code>N+1 ... N+D+1</code> 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.
</p>
<div className="two-column">
<KnotController ref="open-uniform-bspline" />
<BSplineGraphic sketch={this.openUniformBSpline} controller={(owner,
knots) => this.bindKnots(owner, knots, "open-uniform-bspline")}/>
</div>
<h3>Non-uniform B-Splines</h3>
<p>
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,
and that is that any value <code>knots[k+1]</code> should be equal
to, or greater than <code>knots[k]</code>.
</p>
<h2>One last thing: Rational B-Splines</h2>
<p>
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 close to that point the
spline curve will lie, a bit like turning up the gravity of a
control point.
</p>
<div className="two-column">
{ // <KnotController ref="rational-uniform-bspline" />
}
<WeightController ref="rational-uniform-bspline-weights" />
<BSplineGraphic scrolling={true}
sketch={this.rationalUniformBSpline} controller={(owner, knots,
weights, closed) => { // this.bindKnots(owner, knots,
"rational-uniform-bspline"); this.bindWeights(owner, weights,
closed, "rational-uniform-bspline-weights"); }} />
</div>
<p>
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 are an
important type of curve in computer-facilitated design, used a lot
in 3D modelling (as NURBS surfaces) as well as in
arbitrary-precision 2D design due to the level of control a NURBS
curve offers designers.
</p>
<p>
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.
</p>
<h2>Extending our implementation to cover rational splines</h2>
<p>
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.
</p>
<p>
For example, a 2D point <code>(x,y)</code> with weight
<code>w</code> becomes a 3D point <code>(w * x, w * y, w)</code>.
</p>
<p>
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.
</p>
<p>
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 <code>w'</code> and divide all the
regular coordinate dimensions by it, then throw away the weight
information.
</p>
<p>
Based on our previous example, we take the final 3D point
<code>(x', y', w')</code>, which we then turn back into a 2D point
by computing <code>(x'/w', y'/w')</code>. And that's it, we're done!
</p>
</section>
<section id="comments">
<script>
/* ----------------------------------------------------------------------------- *
*
* PLEASE DO NOT LOCALISE THIS FILE
*
* I can't respond to questions that aren't asked in English, so this is one of
* the few cases where there is a content.en-GB.md but you shouldn't change it.
*
* ----------------------------------------------------------------------------- */
</script>
<h1>Comments and questions</h1>
<p>
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
<a href="https://patreon.com/bezierinfo">Patreon page</a> for this
book, or if you prefer to make a one-time donation, head on over to
the
<a
href="https://www.paypal.com/cgi-bin/webscr?cmd=_s-xclick&hosted_button_id=QPRDLNGDANJSW"
>buy Pomax a coffee</a
>
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.
</p>
<p>With that said, on to the comments!</p>
<div id="disqus_thread" />
</section>
</section>
</main>
</body>
</html>
Reference in New Issue View Git Blame Copy Permalink