Rendering Bézier Forms via Forward Differencing
-- Hin Jang
Forthcoming game engines will render, as its base primitive, curves and
non-planar surfaces. One representation is the Bézier form, named after French
mathematician Pierre Bézier. Herein is a very brief review of these entities primarily
used in computer-aided geometric design for the past thirty years. Following this is an
introduction to the method of evaluation known as forward differencing and its application
to rendering curves and surfaces.
A Bézier curve of degree n is defined as
n
----
p(t) = \ pi Bi,n(t)
/
----
i=0
0 <= t <= 1
The points pi form the Bézier control polygon. The
Berstein polynomials are
n!
Bi,n(t) = ---------- (1 - t)n - i ti
i!(n - i)!
With cubic curves (n = 3), for example,
B0,3 = -t3 + 3t2 - 3t + 1
B1,3 = 3t3 - 6t2 + 3t
B2,3 = -3t3 + 3t2
B3,3 = t3
The expression for p(t) is then
p(t) = (1 - t)3p0 + 3t(1 - t)2p1 + 3t2(1 - t)p2 + t3p3
A tensor product Bézier surface p(s, t) of degree m
by n is defined as
m n
---- ----
p(s, t) = \ \ pij Bi,m(s) Bj,n(t)
/ /
---- ----
i=0 j=0
0 <= s <= 1
0 <= t <= 1
The points pi j = (x, y, z)ij form a Bézier
control net, the vertices of the characteristic polyhedron that surrounds the surface.
The Berstein polynomials Bi, m(s) and Bj, n(t) are defined as for
Bézier curves. For a bicubic surface, the expression is
p(s, t) = SPT
where
S = [(1 - s)3 3s(1 - s)2 3s2(1 - s) s3]
| p00 p01 p02 p03 |
P = | p10 p11 p12 p13 |
| p20 p21 p22 p23 |
| p30 p31 p32 p33 |
| (1 - t)3 |
T = | 3t(1 - t)2 |
| 3t2(1 - t) |
| t3 |
Varying the entries of the control net P deforms the underlying
surface, and thus allowing a high degree of shape modification.
—
Foward Differencing
This method of evaulating a polynomial was known since the time of
Isacc Newton and is well-suited for a variety of rendering contexts. Suppose
y = f(x)
represents a function that is differentiable in an open interval
containing x. The differential of y is
dy = f'(x)dx
where dx is the differential of x. if f(x) is known, then
f(x + k) can be approximated via differentials where
f(x + k) ~= f(x) + dy = f(x) + f'(x)dx
k > 0
Within the context of computer graphics programming, the forward
difference Df(t) of a function f(t) is
Df(t) = f(t + k) - f(t)
so
f(t + k) = f(t) + Df(t)
in iterative terms,
fn+1 = fn + Dfn
Evaluating successive function values using the forward difference(s)
is known as forward differencing.
—
Quadratic Bézier Curves
Bézier curves of degree two (i.e., quadratic) are the simplest in the
family of parametric curves. Level-of-detail control involves, partially, the
determination of the ideal resolution at which the following equation is evaluated for a
given scene
p(t) = (1 - t)2p0 + 2t(1 - t)p1 + t2p2
as derived from the general form
n
----
p(t) = \ pi Bi,n(t)
/
----
i=0
0 <= t <= 1
with blending functions
n!
Bi,n(t) = ---------- (1 - t)n - i ti
i!(n - i)!
Forward differencing or some hybrid of such a method will most likely
be used evaluate the network of Bézier curves that populate a virtual environment. For a
quadratic Bézier curve p(t), then,
p(t) = (1 - t)2p0 + 2t(1 - t)p1 + t2p2
= (t2 - 2t + 1)p0 + (-2t2 + 2t)p1 + t2p2
Dp(t) = p(t + k) - p(t)
= ((t + k)2 - 2(t + k) + 1)p0 + (-2(t + k)2 + 2(t + k))p1 +
(t + k)2p2 - [(t2 - 2t + 1)p0 + (-2t2 + 2t)p1 + t2p2]
= (2tk + k2 - 2k)p0 + (-4tk - 2k2 + 2k)p1 + (2tk + k2)p2
D2p(t) = Dp(t + k) - Dp(t)
= (2(t + k)k + k2 - 2k)p0 + (-4(t + k)k - 2k2 + 2k)p1 +
(2(t + k)k + k2)p2 - [(2tk + k2 - 2k)p0 + (-4tk - 2k2 + 2k)p1 +
(2tk + k2)p2]
= 2k2p0 - 4k2p1 + 2k2p2
at initial condition t = 0 with k = STEP = 0.1
p(0) = p0
Dp(0) = (k2 - 2k)p0 + (-2k2 + 2k)p1 + k2p2
= -0.19p0 + 0.18p1 + 0.01p2
D2p(0) = 2k2p0 - 4k2p1 + 2k2p2
= 0.02p0 - 0.04p1 + 0.02p2
Thus, an implementation of forward differencing to render p(t)
would be
#include <stdio.h>
#define STEP 0.1 // range = (0, 1]
typedef struct point Point;
struct point {
float x, y, z;
};
//
// entries are valid for STEP equal to 0.1
static float M[3][3] = {
{ 1.00, 0.00, 0.00 },
{ -0.19, 0.18, 0.01 },
{ 0.02, -0.04, 0.02 },
};
int main(int args, char *argv[]) {
FILE *fp;
Point p[3]; // control points
int i;
if (args < 2) {
fprintf(stderr, "usage: %s <filename>", argv[0]);
exit(1);
}
if ((fp = fopen(argv[1], "r")) == NULL) {
fprintf(stderr, "%s not found\n", argv[1]);
exit(1);
}
i = 0;
//
// the call to fscanf() assumes the datafile consists of
// three rows of three floating point numbers
while (fscanf(fp, "%f %f %f", &p[i].x, &p[i].y, &p[i].z) != EOF)
i++;
Render_Quadratic_BezierCurve(p);
return 0;
}
//
// replace fprintf() with additional code or appropriate
// function calls
void Render_Quadratic_BezierCurve(Point p[3]) {
float x, y, z, dx, dy, dz, ddx, ddy, ddz, t;
x = p[0].x;
y = p[0].y;
z = p[0].z;
dx = M[1][0]*p[0].x + M[1][1]*p[1].x + M[1][2]*p[2].x;
dy = M[1][0]*p[0].y + M[1][1]*p[1].y + M[1][2]*p[2].y;
dz = M[1][0]*p[0].z + M[1][1]*p[1].z + M[1][2]*p[2].z;
ddx = M[2][0]*p[0].x + M[2][1]*p[1].x + M[2][2]*p[2].x;
ddy = M[2][0]*p[0].y + M[2][1]*p[1].y + M[2][2]*p[2].y;
ddz = M[2][0]*p[0].z + M[2][1]*p[1].z + M[2][2]*p[2].z;
fprintf(stdout, "(%f, %f, %f)\n", x, y, z);
for (t = 0.0; t <= 1.0; t += STEP) {
x += dx; dx += ddx;
y += dy; dy += ddy;
z += dz; dz += ddz;
fprintf(stdout, "(%f, %f, %f)\n", x, y, z);
}
}
Although efficient, Render_Quadratic_BezierCurve( ) can be
optimised further. Regardless of the geometry of the curve, p(t) is always
evaluated at fixed increments of size STEP. ideally, however, we want
|| p(t) - p(t + k) || >> 0
Satisfying this condition for all sample points for all quadratic
Bézier curves p(t) requires that k be a function of curvature. in other
words, the curve is sampled at non-fixed intervals of the paramteric variable. an adaptive
foward differencing algorithm evaluates p(t) at non-fixed intervals. To
appreciate the efficiency of forward differencing consider the following table for one
hundred samples (k = 0.01) of the quadratic Bézier curve p(t).
| summary of
arithmetic operations |
| method |
setup |
rendering |
|
|
multiply |
add |
multiply |
add |
| direct evaluation |
0 |
0 |
1300 |
700 |
| forward differencing |
18 |
12 |
0 |
600 |
The computational effort for transformation, rasterisation, and the
loop increment and testing is not included in these figures. In general, forward
differencing performs dn additions per sample for a d-dimensional Bézier
curve of degree n.
—
Quadratic Bézier Surfaces
There are several ways to render a Bézier surface q(s, t), all
of which are analogous to the one-parameter case: a simple curve. The setup, however, is
more complex. Recall that the forward difference equations for the quadratic curve p(t)
at initial condition t = 0 with k = STEP = 0.1 are
p(0) = p0
Dp(0) = (k2 - 2k)p0 + (-2k2 + 2k)p1 + k2p2
= -0.19p0 + 0.18p1 + 0.01p2
D2p(0) = 2k2p0 - 4k2p1 + 2k2p2
= 0.02p0 - 0.04p1 + 0.02p2
or, as a matrix product
| 1.00 0.00 0.00 | | p0 |
D = | -0.19 0.18 0.01 | | p1 |
| 0.02 -0.04 0.02 | | p2 |
The second and third rows of D are the forward differences. In
the case of the quadratic surface q(s, t) with the same initial condition and step
size in s and t
| 1.00 0.00 0.00 | | p00 p01 p02 | | 1.00 -0.19 0.02 |
D = | -0.19 0.18 0.01 | | p10 p11 p12 | | 0.00 0.18 -0.04 |
| 0.02 -0.04 0.02 | | p20 p12 p22 | | 0.00 0.01 0.02 |
= S . P . T
where P is the matrix of control points.
The first row of matrix D now has the values p(0), Dp(0),
D2p(0) in t. the other rows of D contain the first and
second forward differences in s of the forward differences in the first row. The
matrix T is the transpose of S. with matrix D defined as such, an
algorithm to render q(s, t) is
k = 1.0 / (maximum number of steps along t = 0 and s = 0)
X = x-coordinates of the control points
Y = y-coordinates of the control points
Z = z-coordinates of the control points
px = S * X * T
py = S * Y * T
pz = S * Z * T
//
// points along constant t, step size = k
for (0.0 <= s <= 1.0) {
Xd = forward differences for x-coordinate = first row of px
Yd = forward differences for y-coordinate = first row of py
Zd = forward differences for z-coordinate = first row of pz
//
// points along constant s, step size = k
for (0.0 <= t <= 1.0) {
generate points using differences Xd, Yd and Zd
}
update first and second rows of px
update first and second rows of py
update first and second rows of pz
}
An actual implementation of this algorithm is shown below. The context
in which the code was written values modularity higher than efficieny. As such, for
real-time applications, a degree of optimisation is recommended. Note also that the order
in which the points are generated does not lead to an efficient means to send triangle
strips to the graphics subsystem. Modification of the inner loop is therefore required.
//
// written by hin: 19990114
#include <stdio.h>
#define STEP 0.1
//
// entries are valid for STEP = 0.1
static float S[3][3] = {
{ 1.00, 0.00, 0.00 },
{ -0.19, 0.18, 0.01 },
{ 0.02, -0.04, 0.02 },
};
//
// transpose of matrix S
static float T[3][3] = {
{ 1.00, -0.19, 0.02 },
{ 0.00, 0.18, -0.04 },
{ 0.00, 0.01, 0.02 },
};
int main(int args, char *argv[]) {
FILE *fp;
float X[3][3], Y[3][3], Z[3][3]; // control points
int i, j;
if (args < 2) {
fprintf(stderr, "usage: %s <filename>", argv[0]);
exit(1);
}
if ((fp = fopen(argv[1], "r")) == NULL) {
fprintf(stderr, "%s not found\n", argv[1]);
exit(1);
}
//
// the call to fscanf() assumes the datafile consists of nine
// (x, y, z) triplets in row major order
for (j = 0; j < 3; j++)
for (i = 0; i < 3; i++)
(void)fscanf(fp, "%f %f %f", &X[i][j], &Y[i][j], &Z[i][j]);
Render_Quadratic_BezierSurface(X, Y, Z);
return 0;
}
//
// replace fprintf() with additional code or appropriate
// function calls.
void Render_Quadratic_BezierSurface(float X[][3],
float Y[][3],
float Z[][3]) {
register float x, y, z, dx, dy, dz, ddx, ddy, ddz, s, t;
float xp[3][3], yp[3][3], zp[3][3], tp[3][3];
MatrixMultiply(X, T, tp);
MatrixMultiply(S, tp, xp);
MatrixMultiply(Y, T, tp);
MatrixMultiply(S, tp, yp);
MatrixMultiply(Z, T, tp);
MatrixMultiply(S, tp, zp);
for (s = 0.0; s <= 1.0; s += STEP) {
x = xp[0][0];
y = yp[0][0];
z = zp[0][0];
dx = xp[0][1];
dy = yp[0][1];
dz = zp[0][1];
ddx = xp[0][2];
ddy = yp[0][2];
ddz = zp[0][2];
for (t = 0.0; t <= 1.0; t += STEP) {
// note: any system that renders a surface is in
// fact sending triangle vertex data to the
// graphics subsystem. the reason for this
// is purely a matter of resources. to
// faithfully sample a continuous function on a
// closed interval would require extremely high
// sample densities and hence large data volumes
// and processing capacity. it is therefore
// advantageous to reduce the sample density and
// compensate for the artifacts in favour of
// interactive rates of display.
fprintf(stdout, "(%f, %f, %f)\n", x, y, z);
x += dx; dx += ddx;
y += dy; dy += ddy;
z += dz; dz += ddz;
}
AddRows(xp);
AddRows(yp);
AddRows(zp);
}
}
//
// c = a * b
void MatrixMultiply(float a[][3], float b[][3], float c[][3]) {
c[0][0] = a[0][0]*b[0][0] + a[0][1]*b[1][0] + a[0][2]*b[2][0];
c[0][1] = a[0][0]*b[0][1] + a[0][1]*b[1][1] + a[0][2]*b[2][1];
c[0][2] = a[0][0]*b[0][2] + a[0][1]*b[1][2] + a[0][2]*b[2][2];
c[1][0] = a[1][0]*b[0][0] + a[1][1]*b[1][0] + a[1][2]*b[2][0];
c[1][1] = a[1][0]*b[0][1] + a[1][1]*b[1][1] + a[1][2]*b[2][1];
c[1][2] = a[1][0]*b[0][2] + a[1][1]*b[1][2] + a[1][2]*b[2][2];
c[2][0] = a[2][0]*b[0][0] + a[2][1]*b[1][0] + a[2][2]*b[2][0];
c[2][1] = a[2][0]*b[0][1] + a[2][1]*b[1][1] + a[2][2]*b[2][1];
c[2][2] = a[2][0]*b[0][2] + a[2][1]*b[1][2] + a[2][2]*b[2][2];
}
//
// row0 += row1
// row1 += row2
void AddRows(float m[][3]) {
m[0][0] += m[1][0];
m[0][1] += m[1][1];
m[0][2] += m[1][2];
m[1][0] += m[2][0];
m[1][1] += m[2][1];
m[1][2] += m[2][2];
}
further reading
Bernstein, S., Démonstration du théoreme de Weierstrass fondeé sur
le calcul des probabilités. Harkov Soobs. Matem ob-va, 13:1-2, 1912
Bézier, P., Procédé de définition numérique des courbes et
surfaces non mathématiques. Automatisme, XIII(5):189-196, 1968
Bézier, P., Mathematical and Practical Possibilities of UNISURF. In R.
Barnhill and R. Riesenfeld, editors, Computer Aided Geometric Design, Academic
Press, 127-152, 1974
Bézier, P., Essay de définition numérique des courbes et des
surfaces expérimentales. Ph.D. dissertation, University of Paris VI, France, 1977
Boehm, W., "Generating the Bézier Points of B-spline Curves and
Surfaces," Computer Aided Design, 13(6):365-366, 1981
Farin, G., Curves and Surfaces for Computer-Aided Geometric Design:
A Practical Guide, Fourth Edition, Academic Press, San Diego, 1997
Foley, J.D., A. van Dam, S.K. Feiner, and J.F. Hughes, Computer
Graphics Principles and Practice, Second Edition, Addison-Wesley, Reading, 488-491,
1990
Klassen, R.V., "Integer Forward Differencing of Cubic Polynomials:
Analysis and Algorithms," ACM Transactions on Graphics, 10(2):152-181, April
1991
Klassen, R.V., "Exact Integer Hybrid Subdivision and Forward
Differencing of Cubics," ACM Transactions on Graphics, 13(3):240-255, July
1994
Lien, S.L., M. Shantz, and V. Pratt, "Adaptive Forward
Differencing for Rendering Curves and Surfaces," Computer Graphics, SIGGRAPH 1987
Proceedings, 21(4):111-117
Rappoport, A., "Renderings Curves and Surfaces with Hybrid
Subdivision and Forward Differencing," ACM Transcations on Graphics,
10(4):323-341, October 1991
Shantz, M., and S.L. Chang, "Rendering Trimmed NURBS with Adaptive
Forward Differencing," Computer Graphics, SIGGRAPH 1988 Proceedings,
22(4):189-196
Shao, L., and H. Zhou, "Curve Fitting with Bézier Cubics," Graphical
Models and Image Processing, 58(3):223-232, May 1996
Vaishnav, H., and A. Rockwood, "Calculating Offsets of a Bézier
Curve," ACM Proceedings on the Second Symposium on Solid Modeling, 491-492,
1993
Hin Jang (c) March 9, 1999 - permission to copy all or part of
the material contained herein is granted, upon approval from the author via e-mail, provided the copies are not made or
distributed for direct commerical advantage, the author's name and the title of the
article and its date appear, and notice is given that copying is by permission of the
author.
|
