Estimating Tessellation Parameter Intervals for Rational
Curves and Surfaces
Jianmin Zheng
Department of Applied Mathematics
Zhejiang University
Hangzhou, 310027, P.R.CHINA
and
Thomas W. Sederberg
Department of Computer Science
Brigham Young University
Provo, Utah 84602
This paper presents a method for determining a priori a constant parameter interval with which a rational curve or surface can be
tessellated such that the deviation of the curve or surface from its piecewise linear approximation is within a speci¯ed tolerance. The
parameter interval is estimated based on information about the second order derivatives in the homogeneous coordinates, instead
of using a±ne coordinates directly. This new step size can be found with roughly the same amount of computation as the step
size presented in [Cheng 1992], though it can be proven to always be larger than Cheng's step size. In fact, numerical experiments
show the new step is typically orders of magnitude larger than the step size in [Cheng 1992]. Furthermore, for rational cubic and
quartic curves, the new step size is generally twice as large as the step size found by computing bounds on the Bernstein polynomial
coe±cients of the second derivatives function.
Categories and Subject Descriptors: I.3.5 [Computer Graphics]: Computational Geometry and Object Modelling|geometric
algorithms; J.6 [Computer Applications]: Computer-Aided Engineering|computer-aided design
General Terms: Algorithms
Additional Key Words and Phrases: Rational curves and surfaces, tessellation, °atness, derivative bounds, step size, projection
distance
1. INTRODUCTION
In computer graphics and modeling, parametric curves and surfaces are often tessellated into piecewise linear
segments for rendering [Lane and Carpenter 1979; Rappoport 1991; Abi-Ezzi and Wozny 1991; Klassen 1994]
mesh generation [Sheng and Hirsch 1992; Piegl and Richard 1995] and surface/surface intersection [Wang 1984;
Filip et al. 1986]. In such applications, the approximation error is typically taken as the maximal distance
between the original and the approximating segments.
Many criteria and methods have been developed for the task of assuring that a piecewise linear tessellation
satis¯es an error bound [Lane and Riesenfeld 1980; Schaback 1993; Elber 1996; Tookey and Cripps 1997]. One
popular approach is to determine a global parameter interval or step size ± that is valid over the entire domain
Name: Jianmin Zheng
A±liation: Zhejiang University
Address: Hangzhou, 310027, P.R.CHINA
Name: Thomas W. Sederberg
A±liation: Brigham Young University
Address: Provo, Utah 84602; Email: tom@byu.edu
Permission to make digital or hard copies of part or all of this work for personal or classroom use is granted without fee provided
that copies are not made or distributed for pro¯t or direct commercial advantage and that copies show this notice on the ¯rst page
or initial screen of a display along with the full citation. Copyrights for components of this work owned by others than ACM must
be honored. Abstracting with credit is permitted. To copy otherwise, to republish, to post on servers, to redistribute to lists, or to
use any component of this work in other works, requires prior speci¯c permission and/or a fee. Permissions may be requested from
Publications Dept, ACM Inc., 1515 Broadway, New York, NY 10036 USA, fax +1 (212) 869-0481, or permissions@acm.org.
2 ¢
[Wang 1984; Filip et al. 1986; Cheng 1992; Sheng and Hirsch 1992; Klassen 1994]. Then the tessellation can
be generated either by recursively splitting the curve or surface in half to a depth computed from ±, or by just
sampling the curve or surface at points separated by a parameter interval ±. Such an approach lends itself
well to forward di®erencing. A step size should ideally be as large as possible without violating the prescribed
tolerance. The calculation of such a step size is usually not easy for rational curves and surfaces.
The general approach to obtaining such a step size is to establish a relatively simple relationship between the
step size and a bound on the maximal deviation. Once the bound is found, the step size follows immediately.
Based on theorems from approximation theory [de Boor 1978], Wang proposed using bounds on second order
derivatives of the curve or surface to estimate deviation from the linear approximation [Wang 1984]. This
approach was later enhanced by Filip et al [Filip et al. 1986]. For a C2 continuous parametric curve r(t); t 2 [®; ¯],
a bound on the displacement of r(t) from its linear approximation L(t) which interpolates r(®) and r(¯) is given
by
sup
t2[®;¯] kr(t) ¡ L(t)k ·
1
8±2 sup
t2[®;¯] kr00(t)k; ±= ¯ ¡ ® (1)
For a C2 continuous parametric surface r(u; v) de¯ned over a right triangle T with two vertical edge lengths of
±u and ±v, a bound on the deviation of r(u; v) from the base triangle L(u; v) which linearly interpolates r(u; v)
at three vertices of T is
sup
(u;v)2T kr(u; v) ¡ L(u; v)k ·
1
8
¡
Duu±2u
+ 2Duv±u±v + Dvv±2
v
¢
(2)
where
Duu = sup
(u;v)2T
°°°°
@2r(u; v)
@u2
°°°°
; Duv = sup
(u;v)2T
°°°°
@2r(u; v)
@u@v
°°°°
; Dvv = sup
(u;v)2T
°°°°
@2r(u; v)
@v2
°°°°
(3)
In these formulas, we need to compute the sup's of the second order (partial) derivatives over the domain of
the curve or surface. For polynomial curves and surfaces, the convex hull property of the control points of the
Bezier representation o®ers a straightforward way to compute upper bounds. For rational curves and surfaces,
estimating the sup's of their second order (partial) derivatives is much more expensive. The rational case is
important not only because it is more powerful in modeling than the polynomial scheme and because NURBS
have become an industry standard in CAD/CAM systems, but also because perspective transformation changes
a polynomial curve or surface to a rational curve or surface. (While these formulae apply to any C2 curve or
surface, including C2 NURBS, they are more commonly applied to individual B¶ezier curve segments or surface
patches).
To date, not many papers have been published that address the rational case. Several papers have been
written that study e±cient ways to obtain bounds on ¯rst derivatives for rational curves and surfaces [Floater
1992; Hermann 1992; Wang et al. 1997] but these techniques are still somewhat costly. Experience with the
¯rst derivative suggests that global bounds on second derivatives are either very expensive to compute, or very
conservative. Cheng proposed a new method that does not require bounds on the second derivatives [Cheng
1992]. He instead treats a rational surface as a polynomial surface in projective 4D space and shows that, given
a tolerance ² and a rational surface r(u; v) = R(u; v)=w(u; v), we only need to compute a bound for its associated
homogeneous surface S(u; v) = (R(u; v); w(u; v)) with a new tolerance ²H, which guarantees the ²-°atness of
r(u; v). The new tolerance is
²H =
8>>>><
>>>>:
inf jw(u;v)j²
(1+(sup jjr(u;v)jj¡²)2)1=2 ; 0 · ² · sup jjr(u; v)jj
inf jw(u; v)j²; sup jjr(u; v)jj < ² · 2 sup jjr(u; v)jj
+1; 2 sup jjr(u; v)jj < ²
(4)
Thus the problem of evaluating the bounds on r(u; v)'s second order derivatives turns to estimating S(u; v)'s,
which is much easier. However, numerical experiments show that this estimation is often too loose.
In this paper, we adhere to Cheng's idea that using homogeneous coordinates for rational curves and surfaces
makes derivatives simple, and hence we attempt to bound the maximum deviation based on the associated
homogeneous coordinates. Section 2 derives a simple bound expressed by the homogeneous coordinates for
¢ 3
the distance of two projection points. Then the new step size formulas are developed for rational curves and
surfaces in Sections 3 and 4 respectively. It is found that the new bound is tighter than Cheng's. When all the
weights of the rational curves or surfaces are the same, the bound reduces to that of polynomial B¶ezier curves or
surfaces. It is expected that the improved bound will lead to more e±cient algorithms in graphics and modeling
applications.
2. PROJECTION DISTANCE ESTIMATION
Our basic strategy is to perform evaluation in homogeneous space. This requires us to approximate the e®ect of
the perspective projection e±ciently. The following theorems provide relationships between homogeneous and
Euclidean spaces for this purpose. Throughout the paper, En stands for n-dimensional Euclidean space.
Theorem 1. Let Q1 = (R1; w1) and Q2 = (R2; w2) with w1; w2 > 0 be two points in En+1, whose
corresponding projections are r1 = R1=w1; r2 = R2=w2 2 En, and let r be any number not smaller than
max(kr1k; kr2k). Then
kr1 ¡ r2k · kR1 ¡ R2k + (r ¡ kr1 ¡ r2k)jw1 ¡ w2j
min(w1; w2)
(5)
Proof: Without loss of generality, we assume that w1 · w2. Then
kr1 ¡ r2k = kR1
w1 ¡ R2
w2k · kR1
w1 ¡ R1
w2 k + kR1
w2 ¡ R2
w2 k
· kR1
w1 kjw1¡w2j w2
+ kR1¡R2k w2 · r jw1¡w2j w2
+ kR1¡R2k w2
= rjw1¡w2j+kR1¡R2k w1+w2¡w1
= rjw1¡w2j+kR1¡R2k min(w1;w2)+jw1¡w2j
Thus we obtain
(min(w1; w2) + jw1 ¡ w2j)kr1 ¡ r2k · rjw1 ¡ w2j + kR1 ¡ R2k
i.e.
min(w1; w2)kr1 ¡ r2k · kR1 ¡ R2k + (r ¡ kr1 ¡ r2k)jw1 ¡ w2j
Therefore
kr1 ¡ r2k · kR1 ¡ R2k + (r ¡ kr1 ¡ r2k)jw1 ¡ w2j
min(w1; w2)
Theorem 1 has an intuitive geometric interpretation. See Figure 1 for an example in the case of n = 1.
The multiplication of w1 and the Euclidean distance of r1 and r2 equals the length of Q1I, which is the sum
of kQ1Hk and kHIk. Note that kHIk = kQ2Hk= tan ®, kQ1Hk is kR1 ¡ R2k and kQ2Hk is just jw1 ¡ w2j.
Therefore the inequality (5) is equivalent to tan ® = kOCk=kr2Ck ¸ 1=(r ¡ kr1 ¡ r2k). A similar but looser
estimation was given by Tiller [Tiller 1992]and Schaback [Schaback 1993].
Theorem 2. Let Q1 = (R1; w1) and Q2 = (R2; w2) be two points in En+1, and let ² and r be two positive
numbers. If Q1 and Q2 satisfy
w1; w2 ¸ 0; max(k
R1
w1 k; k
R2
w2 k) · r (6)
and
kR1 ¡ R2k + (r ¡ ²)jw1 ¡ w2j
min(w1; w2) · ² (7)
then the distance of their projections r1 = R1=w1 and r2 = R2=w2 is less than or equal to ²:
kr1 ¡ r2k · ² (8)
4 ¢
a r
1
C
I
H
O
w1
w2
r1 r2
Q2
Q1
R1 R2
Fig. 1. Geometric interpretation for the case of n=1.
Proof: Let ´ = kr1 ¡ r2k. Then
´ · kR1¡R2k+(r¡´)jw1¡w2j min(w1;w2) = kR1¡R2k+(r¡²)jw1¡w2j min(w1;w2) + (²¡´)jw1¡w2j min(w1;w2)
· ² + (² ¡ ´) jw1¡w2j min(w1;w2)
Thus (² ¡ ´)
³
1 + jw1¡w2j min(w1;w2)
´
¸ 0, which arrives at the conclusion: ² ¸ ´ = kr1 ¡ r2k.
This theorem shows that the closeness requirement of two points can be imposed on a simple expression
represented by the homogeneous coordinates.
3. STEP SIZE FOR RATIONAL CURVES
3.1 General theory
Since curves have simpler equations, we discuss the curve case ¯rst. For a rational curve r(t) = R(t)
w(t) ; t 2 [0; 1],
the tessellation problem can be formulated: given a tolerance ², ¯nd a maximum parameter step ± such that
sup
t2[®;¯] jjr(t) ¡ L(t)jj · ² (9)
where ®; ¯ 2 [0; 1] and ± = ¯ ¡ ®, L(t) is the line segment connecting r(®) and r(¯) parameterized as follows:
L(t) =
Ln(t)
Ld(t)
=
[(¯ ¡ t)R(®) + (t ¡ ®)R(¯)]=(¯ ¡ ®)
[(¯ ¡ t)w(®) + (t ¡ ®)w(¯)]=(¯ ¡ ®)
(10)
In practice, we typically content ourselves with ¯nding a step size for which (9) is an upper bound. Surprisingly,
situations exist for which a step size ±1 < ± can actually violate the ²-error, even though a step size of ± does
not. Figure 2 illustrates such a case.
In light of this problem, we modify our de¯nition of ± as follows:
± = supf½j sup
t2[®;®+´] jjr(t) ¡ L(t)jj · ² for any positive ´ · ½; ®; ® + ½ 2 [0; 1]g
With Theorem 2, we need to compute a bound for kR(t) ¡ Ln(t)k + (r ¡ ²)jw(t) ¡ Ld(t)j. The tighter
the bound, the larger (and the more economical) the step size ±. The formula (1) can be simply used here
by ¯nding the bounds on kR00(t)k and jw00(t)j, and combining them. However, this is di®erent from taking
¢ 5
>e
<=e
r(d)
r(d1)
r(0)
Fig. 2. Line r(0)r(±1) is not within ²-°atness while r(0)r(±) is.
kR00(t)k + (r ¡ ²)jw00(t)j as a whole. In general, the latter is smaller than the former since the maxima for
kR00(t)k and jw00(t)j may occur at di®erent parameter values. In the following we show that the proof of
Theorem 2 in [Filip et al. 1986] can be generalized for this.
Lemma 1. Let Q(t) = (R(t); w(t)) 2 En+1 be a C2 vector-valued function de¯ned over [®; ¯] with Q(®) =
Q(¯) = 0. Let k be any nonnegative number. Then
sup
t2[®;¯]
(kR(t)k + kjw(t)j) ·
(¯ ¡ ®)2
8
sup
t2[®;¯]
(kR00(t)k + kjw00(t)j) (11)
Proof: Let d(t) = kR(t)k + kjw(t)j. Then d(t) is a continuous function in [®; ¯] with d(®) = d(¯) = 0 and
d(t) ¸ 0, which means d(t) can attain its maximum value at some » 2 (®; ¯). It is obvious that d(t) is C2
continuous in [®; ¯] possibly except those points where w(t) = 0 or R(t) = 0. Now we discuss the problem in
three cases:
case 1 . w(») = 0.
In this case, d(») = kR(»)k. From (1),
kR(»)k ·
(¯ ¡ ®)2
8
sup
t2[®;¯] kR00(t)k
Let ´ be so that kR00(´)k = sup
t2[®;¯] kR00(t)k. We have
d(») = kR(»)k · (¯¡®)2
8 kR00(´)k
· (¯¡®)2
8 (kR00(´)k + kjw00(´)j) · (¯¡®)2
8 sup
t2[®;¯]
(kR00(´)k + kjw00(´)j)
case 2 . R(») = 0.
Then d(») = kjw(»)j. By a like argument as in case 1, we obtain
d(») =
(¯ ¡ ®)2
8 k sup
t2[®;¯] jw00(t)j ·
(¯ ¡ ®)2
8
sup
t2[®;¯]
(kR00(t)k + kjw00(t)j)
case 3 . w(») 6= 0; R(») 6= 0.
Thus d(t) is C2 at t = » and d0(») = 0, i.e.
d0(») =
R
kRk
R0(») + k
w
jwj
w0(») = 0 (12)
Now we express Q(t) in its Taylor expansion about »:
(R(t); w(t)) = (R(»); w(»)) + (R0(»); w0(»))(t ¡ ») +
Z t
»
(R00(s); w00(s))(t ¡ s)ds (13)
6 ¢
Without loss of generality, we assume that » 2 [®+¯
2 ; ¯). Then substituting t = ¯ into the equation (13) yields
0 = (R(¯); w(¯)) = (R(»); w(»)) + (R0(»); w0(»))(¯ ¡ ») +
Z ¯
»
(R00(s); w00(s))(¯ ¡ s)ds (14)
Taking the dot product of equation (14) with ( R(»)
kR(»)k
; k w(»)
jw(»)j
), we obtain
0 = kR(»)k + kjw(»)j+0+
Z ¯
»
µ
R(»)R00(s)
kR(»)k
+ k
w(»)w00(s)
jw(»)j
¶
(¯ ¡ s)ds
Thus
d(») = kR(»)k + kjw(»)j =
¯¯¯
R ¯
» (R(»)R00(s)
kR(»)k
+ k w(»)w00(s)
jw(»)j
)(¯ ¡ s)ds
¯¯¯
·
R ¯
» (jR(»)R00(s)
kR(»)k j + kjw(»)w00(s)
jw(»)j j)(¯ ¡ s)ds ·
R ¯
» (kR00(s)k + kjw00(s)j)(¯ ¡ s)ds
· sup
t2[®;¯]
(kR00(t)k + kjw00(t)j)
R ¯
» (¯ ¡ s)ds = sup
t2[®;¯]
(kR00(t)k + kjw00(t)j) (¯¡»)2
2
· (¯¡®)2
8 sup
t2[®;¯]
(kR00(t)k + kjw00(t)j)
Now we are ready to derive the formula for the step size.
Theorem 3. Given a C2 rational curve r(t) = R(t)
w(t) ; t 2 [®; ¯] with w(t) > 0, and r is a number satisfying
r ¸ sup
t2[®;¯] kr(t)k = sup
t2[®;¯] kR(t)
w(t) k. Let L(t) = Ln(t)
Ld(t) be de¯ned as in (10), then sup
t2[®;¯] kr(t) ¡ L(t)k · ² if
± = ¯ ¡ ® ·
8>>>>>><
>>>>>>:
q
8 inftfw(t)g²
supt(kR00(t)k+(r¡²)jw00(t)j); ² < r
q
8 inftfw(t)g²
supt kR00(t)k
; r· ² < 2r
1; 2r · ²
(15)
Proof: First let us consider the case of 2r · ². Since
L(t) =
R(®)(¯ ¡ t) + R(¯)(t ¡ ®)
w(®)(¯ ¡ t) + w(¯)(t ¡ ®)
=
R(®)
w(®)
(¯ ¡ t)w(®)
(¯ ¡ t)w(®) + (t ¡ ®)w(¯)
+
R(¯)
w(¯)
(t ¡ ®)w(¯)
(¯ ¡ t)w(®) + (t ¡ ®)w(¯) ;
kL(t)k · max(kR(®)
w(®) k; kR(¯)
w(¯) k) · r. Then for any t 2 [®; ¯], kr(t) ¡ L(t)k · kr(t)k + kL(t)k · r + r = 2r · ².
Second, in the case of r · ² < 2r, applying (1) to R(t), we get
sup
t2[®;¯] kR(t) ¡ Ln(t)k ·
(¯ ¡ ®)2
8
sup
t2[®;¯] kR00(t)k
As Ld(t) is a convex combination of w(®) and w(¯), the value of Ld(t) is between w(®) and w(¯). Hence
sup
t2[®;¯]
kR(t) ¡ Ln(t)k + (r ¡ ²)jw(t) ¡ Ld(t)j
min(w(t); Ld(t)) · sup
t2[®;¯]
kR(t) ¡ Ln(t)k
min(w(t); Ld(t)) ·
±2
8
sup
t2[®;¯] kR00(t)k
inf
t2[®;¯]
w(t) · ²
By Theorem 2, sup
t2[®;¯] kr(t) ¡ L(t)k · ².
Finally, for the case of ² < r, r ¡ ² > 0. Note that R(®) ¡ Ln(®) = R(¯) ¡ Ln(¯) = 0 and w(®) ¡ Ld(®) =
w(¯) ¡ Ld(¯) = 0. From Lemma 1,
sup
t2[®;¯]
(kR(t) ¡ Ln(t)k + (r ¡ ²)jw(t) ¡ Ld(t)j) ·
(¯ ¡ ®)2
8
sup
t2[®;¯]
(kR00(t)k + (r ¡ ²)jw00(t)j)
¢ 7
So
sup
t2[®;¯]
kR(t) ¡ Ln(t)k + (r ¡ ²)jw(t) ¡ Ld(t)j
min(w(t); Ld(t)) ·
sup
t2[®;¯]
(kR00(t)k + (r ¡ ²)jw00(t)j)
inft w(t)
(¯ ¡ ®)2
8 · ²
Thus the proof is completed by Theorem 2.
3.2 Rational Bezier curves
For a rational B¶ezier curve r(t) de¯ned by
r(t) =
R(t)
w(t)
=
Pn
i=0
PiwiBn
i (t)
Pn
i=0
wiBn
i (t)
; t 2 [0; 1] (16)
R(t) and w(t) are in B¶ezier form. The estimation for the sup's of their second order derivatives and for inffw(t)g is straightforward due to the convex-hull property of polynomials in Bernstein form. De¯ne ¢ to be a forward
di®erence operator so that ¢Pi = Pi+1 ¡ Pi. We have
Corollary 1. Let r(t) be a rational B¶ezier curve (16), r = max
i kPik and w = min
i fwig > 0. For any ² > 0,
the deviation of the curve segment over [®; ¯] µ [0; 1] from its fractional linear approximation is within ², i.e.
sup
t2[®;¯] kr(t) ¡ L(t)k · ² if ± = ¯ ¡ ® satis¯es
± ·
8>>>>>>>><
>>>>>>>>:
r
8w²
n(n¡1) max
i
(k¢2(wiPi)k+(r¡²)j¢2wij); ² < r
r
8w²
n(n¡1) max
i k¢2(wiPi)k
; r· ² < 2r
1; 2r · ²
(17)
Proof: The second order derivative of a B¶ezier curve is also a B¶ezier curve with degree lower by 2, e.g.
(R00(t); w00(t)) = n(n ¡ 1)
nX¡2
i=0
¢2(wiPi; wi)Bn¡2
i (t) (18)
The conclusion follows immediately from the inequalities
kR00(t)k · n(n ¡ 1)
nX¡2
i=0
k¢2(wiPi)kBn¡2
i (t) · n(n ¡ 1) max
i k¢2(wiPi)k (19)
and
kR00(t)k + (r ¡ ²)jw00(t)j · n(n ¡ 1)
nP¡2
i=0
(k¢2(wiPi)k + (r ¡ ²)j¢2wij)Bn¡2
i (t)
· n(n ¡ 1) max
i
(k¢2(wiPi)k + (r ¡ ²)j¢2wij)
(20)
It is possible for one of the wi to be zero, but still have a well-de¯ned curve. In this case, a better estimation
for inffw(t)g can be obtained by ¯rst degree elevating or subdividing the polynomial w(t) and then taking
minfwig of the resulting weights.
Examples.
(1) Quadratic rational B¶ezier curve: the step size for the case of ² < r is
± =
s
4 minfw0; w1; w2g²
kw2P2 ¡ 2w1P1 + w0P0k + (r ¡ ²)jw2 ¡ 2w1 + w0j
8 ¢
(2) Cubic rational B¶ezier curve: the step size for the case of ² < r is
± =
s
4 minfw0; w1; w2; w3g²
3 maxfk¢2(w0P0)k + (r ¡ ²)j¢2w0j; k¢2(w1P1)k + (r ¡ ²)j¢2w1jg
Remarks.
(1) Here we show that the new step size ±n given by (17) is larger than Cheng's. For a rational B¶ezier curve,
Cheng's step size is
±c =
s
8²H
n(n ¡ 1) maxi jj¢2(wiPi; wi)jj
and ²H is de¯ned in (4). When ² < r,
(8w²=n(n ¡ 1)±2
c )2 = (1 + (r ¡ ²)2) max
i
(k¢2(wiPi)k2 + j¢2wij2)
= max
i
(k¢2(wiPi)k2 + (r ¡ ²)2j¢2wij2 + (r ¡ ²)2k¢2(wiPi)k2 + j¢2wij2)
¸ max
i
(k¢2(wiPi)k2 + (r ¡ ²)2j¢2wij2 + 2(r ¡ ²)k¢2(wiPi)kj¢2wij)
= max
i
(k¢2(wiPi)k + (r ¡ ²)j¢2wij)2 = (8w²=n(n ¡ 1)±2n
)2
so ±n ¸ ±c. The other cases are trivial.
(2) When all weights are identical, the rational B¶ezier curve reduces to a polynomial curve
r(t) =
Xn
i=0
PiBn
i (t); t 2 [0; 1]:
In this case ¢2wi = 0, so the new step size is the same as that for a B¶ezier curve, which was given in [Filip
et al. 1986], i.e.
± ·
s
8²
n(n ¡ 1) maxi jj¢2Pijj
(21)
However Cheng's estimation, in general, does not reduce to | and is more conservative than | (21).
3.3 Further improvements
Further improvements without much additional computation are possible. The ¯rst observation is that the
bound for the second derivative curve (18) might be tightened. Farin points out that a rational curve can be
bounded in a tighter convex hull generated by the end points and so-called \weight points" [Farin 1993]. The
weight points are obtained after an additional subdivision without doing the full work of subdividing. In our
case, the weight points are just the average of the adjacent homogeneous points. That is
Wi = (WRi;Wwi) = n(n ¡ 1)
¢2(wi¡1Pi¡1; wi¡1)+¢2(wiPi; wi)
2 ; i = 1¢ ¢ ¢ n ¡ 2 (22)
If we also denote the end points by W0 = n(n ¡ 1)¢2(w0P0; w0) and Wn¡1 = n(n ¡ 1)¢2(wn¡2Pn¡2; wn¡2),
then the curve (18) lies in the convex hull that includes only Wi (i = 0; ¢ ¢ ¢ ; n¡1). Therefore for any t 2 [0; 1],
there exist nonnegative coe±cients ¸i(t) such that
(R00(t); w00(t)) =
nX¡1
i=0
¸i(t)Wi; ¸0(t) + ¸1(t) + ¢ ¢ ¢ + ¸n¡1(t) = 1
¢ 9
Therefore
kR00(t)k ·
nP¡1
i=0 kWRik¸i(t) · max
i kWRik
= n(n ¡ 1) max
i fk¢2(w0P0)k; k¢2(wi¡1Pi¡1)+¢2(wiPi)k 2 ; k¢2(wn¡2Pn¡2)kg
(23)
and
kR00(t)k + (r ¡ ²)jw00(t)j ·
nP¡1
i=0
(kWRik + (r ¡ ²)jWwij)¸i(t)
· max
i
(kWRik + (r ¡ ²)jWwij)
= n(n ¡ 1) max
i fk¢2(w0P0)k + (r ¡ ²)j¢2w0j;
k¢2(wi¡1Pi¡1)+¢2(wiPi)k 2 + (r ¡ ²) j¢2wi¡1+¢2wij 2 ;
k¢2(wn¡2Pn¡2)k + (r ¡ ²)j¢2wn¡2jg
(24)
Usually this gives a tighter bound. The improvement depends on the distribution of the control points of curve
(18). If the adjacent control points are close, the averages formed will not be much less in magnitude than the
original ones. For a cubic rational curve, (18) is a linear polynomial, in which case there is no improvement in
using the weight-point method.
The second observation is that translation does not change the shape of the curve and the translated curve
should have the same step size as the original curve. However the formula (15) or (17) is dependent on the
a±ne coordinate system. This suggests that we move the original point of the a±ne coordinate system to make
the value of r = max
i kPik as small as possible. One simple way is to ¯nd the min-max bounding box of the
rational curve and to choose the center of the box as the origin of a local coordinate system, or we can ¯nd a
bounding circle or sphere and then use its center as the origin. The method proposed above can be applied to
these control points in their local coordinate system, which also means that the tessellation will be invariant
to translation. It should be pointed out that this is just a heuristic approach. When the control points are
moved, the values of k¢2(wiPi)k are also changed. So sometimes the result might be worse. But the numerical
examples indicate that in most cases doing a certain translation gives a better step size, especially when the
control points are located very asymmetrically around the original point.
Summarizing, we have the following algorithm:
Input: the control points Pi and weights wi, the error tolerance ².
Output: a global step size
Procedure:
Step 1. Find a min-max bounding box of the curve, and compute its center point, denoted by C.
Step 2. Translate the control points by C: Pi Ã Pi ¡ C.
Step 3. Compute the B¶ezier representation (18) for the second derivative of the homogeneous B¶ezier curve.
Step 4. Compute sup(kR00(t)k + (r ¡ ²)jw00(t)j) or supkR00(t)k by formula (24) or (23).
Step 5. compute the step size ± by (15).
10 ¢
3.4 Numerical experiments
We present a few planar curve examples. The data de¯ning the rational B¶ezier curves are generated randomly.
Ci(t) represents the curve of degree i.
C1(t) : (2; 5; 5:6); (1; 8; 0:7)
C2(t) : (¡3;¡10; 0:96); (6; 8; 2:3); (2; 4; 0:63)
C3(t) : (19; 61; 0:08); (¡61; 52; 0:5); (17; 55; 1); (49;¡20; 0:4)
C4(t) : (1; 5; 6:1); (7; 7; 0:39); (¡8;¡10; 18:4); (¡1;¡10; 1:1); (¡6;¡3; 0:03)
C5(t) : (53;¡6; 0:7); (¡7; 66; 1:8); (¡64;¡46; 147); (¡71; 43; 6:6); (97;¡68; 4); (¡66; 57; 0:7)
C6(t) : (36;¡23; 1:7); (48; 54; 0:8); (14;¡13; 0:2); (64; 13; 1); (¡68; 54; 1:4); (43;¡1; 0:4); (34; 92; 0:2)
C7(t) : (9; 9; 1:5); (¡4; 0; 3:1); (5; 0; 3:3); (¡7; 0; 2:7); (10;¡10; 2:6); (2;¡3; 0:7);
(¡4; 9; 1:1); (1;¡5; 1:3)
C8(t) : (3; 2; 0:2); (¡5; 8; 1:6); (1; 4; 0:8); (7; 10; 0:8); (1;¡8; 1:1); (¡7; 5; 1:3); (0; 10; 0:4);
(4; 1; 0:7); (6;¡1; 2:3)
In each tuple, the ¯rst two numbers stand for the control point's x; y-components, and the third for the weight.
Di®erent methods are used to compute the step size ± with the tolerance ² = 0:1. The results are listed in Table 1.
Here, we test the new methods with just using formula (17), or using weight-point method, or translating the
control points ¯rst and then using (17), or translating the control points ¯rst and then using weight-point
method. They are respectively denoted by \new", \new-w", \new-t" and \new-tw". Cheng's method is labelled
\cheng". A well-known method for ¯nding a step in a parametric curve is
± =
r
8²
D
(25)
where D is an upper bound on the magnitude of the second derivative vector of the curve. Variations of this
step size are discussed in [Wang et al. 1997; Filip et al. 1986]. For a rational curve of degree n, the second
derivative vector is degree 3n. For comparison with the new method, we computed in two di®erent ways a step
size using the method in (25). First, we expressed the second derivative in rational Bernstein form and took
D to be the largest distance from the origin to any control point. This method is labelled \approx-D" in the
table. We also went to the expense of computing the tight upper bound on the second derivative curve. This
amounted to ¯nding the global maximum of a degree 6n polynomial over the [0; 1] interval, a computation that
is far too costly for practical use, but which is of interest because it yields the largest step size of which (25) is
capable. The step size produced in this way is labelled \exact-D" in the table.
± C1(t) C2(t) C3(t) C4(t) C5(t) C6(t) C7(t) C8(t)
cheng 1.0 0.017 0.0003 0.0001 0.00008 0.00018 0.0027 0.0012
new 1.0 0.0549 0.0075 0.0015 0.0007 0.0035 0.0108 0.0072
new-w 1.0 0.0549 0.0075 0.0019 0.00098 0.0048 0.0132 0.0072
new-t 1.0 0.0551 0.0081 0.0016 0.0007 0.0037 0.0111 0.0075
new-tw 1.0 0.0551 0.0081 0.0021 0.001 0.0052 0.0132 0.0075
approx-D 0.048 0.0335 0.0039 0.0015 0.001 0.0052 0.0121 0.0032
exact-D 0.085 0.0457 0.0039 0.0019 0.0013 0.0058 0.0130 0.0036
Table 1. The step sizes computed by di®erent methods with ² = 0:1.
We also ran several sets of numerical comparisons on several hundred randomly generated test cases of degree
three and four. In one batch of tests, the (x; y) coordinates of the control points were generated as random
numbers in the interval ([¡1000; 1000]; [¡1000; 1000]). We also ran tests in the intervals ([¡100; 100]; [¡100; 100])
and ([¡10; 10]; [¡10; 10]). Control point weights were generated as the ratio of two random numbers in the
intervals [1; 10000], and ² was ¯xed at 0:1. In each batch, we ran several hundred test cases, computed the ratio
between the \new-tw" step size, and those produced by Cheng's method and by the \approx-D" method. The
results are displayed in Table 2.
Two things are noteworthy in this table. First, it is somewhat surprising that \new-tw" gives on average twice
as large of a step size as \approx-D". This can be explained. In the case of a polynomial B¶ezier curve, it is easily
shown that the two methods will give exactly the same step size. However, if the weights vary signi¯cantly,
¢ 11
Interval Degree new¡tw
cheng
new¡tw
approx¡D
([¡1000; 1000]; [¡1000; 1000]) 3 66.4 2.03
([¡1000; 1000]; [¡1000; 1000]) 4 72.6 1.95
([¡100; 100]; [¡100; 100]) 3 20.2 1.92
([¡10; 10]; [¡10; 10]) 3 6.8 2.56
Table 2. Average ratio of step sizes with ² = 0:1.
\new-tw" is often able to provide a larger step size than \approx-D" because it is ¯nding the distance between
a rational curve segment and a fractional-linearly parametrized line segment. By contrast, \approx-D" ¯nds the
distance between the rational curve segment and a linearly-parametrized line segment.
The second curiosity in Table 2 is that Cheng's method appears to improve as the size of the control point
coordinates decreases! Inspection of Cheng's formulae con¯rms that this should indeed occur. Furthermore, it
can be seen that Cheng's method will produce a di®erent step size even if the control points and ² are scaled the
same amount! This suggests that Cheng's method might be able to produce a larger step size if we ¯rst scale
the control points and ² by some constant, and that there might be an optimal scale factor that will produce
an optimally large step size. Closer study reveals that, in fact, our step size is exactly the one produced by
Cheng's method after applying the optimal scaling!
4. STEP SIZES FOR RATIONAL SURFACES
The development for rational surfaces parallels that for rational curves. In this section we just list the results
without proofs except the following lemma. Of course, when tessellating two adjacent patches, it is a good idea
to compute step sizes for all four boundary curves independently from the step size for the patch interior to
avoid cracks (see [Rockwood et al. 1989]).
Lemma 2. Let Q(u; v) = (R(u; v); w(u; v)) be a C2 vector-valued function de¯ned over a right triangle T
with vertices A;B and C of the form B = A + (±u; 0) and C = A + (0; ±v), and vanish at these three vertices
Q(A) = Q(B) = Q(C) = 0. Then for any nonnegative number k,
sup
(u;v)2T
(kR(u; v)k + kjw(u; v)j) · (Duu±2u
+ 2Duv±u±v + Dvv±2
v)=8 (26)
where
Duu = sup
(u;v)2T
(kR00 uu(u; v)k + kjw00 uu(u; v)j)
Duv = sup
(u;v)2T
(kR00 uv(u; v)k + kjw00 uv(u; v)j)
Dvv = sup
(u;v)2T
(kR00 vv(u; v)k + kjw00 vv(u; v)j)
(27)
Proof: This lemma can be proven in a way similar to that in [Filip et al. 1986]. Let D(u; v) = kR(u; v)k +
kjw(u; v)j. D(u; v) is a continuous function and D(u; v) ¸ 0; D(A) = D(B) = D(C) = 0. Then D(u; v) must
get its maximum value at some point P0 = (u0; v0) in T. Suppose P0 lies in region T1 (see Figure 3). The other
cases can be proven similarly.
Let V = P0 ¡ A = (l cos µ; l sin µ) with l = kP0 ¡ Ak and µ the angle between AB and V . Now consider the
curve g(t) from A to P0 on Q(u; v). Let g(t) = Q(A+tV ) and d(t) = D(A+tV ) = kR(A+tV )k+kjw(A+tV )j.
Then g(0) = Q(A) and g(1) = Q(P0). So d(0) = D(A) = 0 and d(1) = D(P0) = sup
(u;v)2TfD(u; v)g. Like the
argument in Lemma 1, here we also consider three cases:
case 1 . w(P0) = 0. Applying (2) gives
D(u; v) · d(P0) = kR(P0)k · 1
8
¡
sup kR00 uuk±2u
+ 2 sup kR00 uvk±u±v + sup kR00 vvk±2
v
¢
· 1
8
¡
Duu±2u
+ 2Duv±u±v + Dvv±2
v
¢
12 ¢
q
du
dv
V
A B
C
P0 T0
T2 T1
T3
Fig. 3. The region of parametric triangle T.
case 2 . R(P0) = 0. Just as in case 1,
D(u; v) · d(P0) = kjw(P0)j · 1
8
¡
sup jw00 uuj±2u
+ 2 sup jw00 uvj±u±v + sup jw00 vvj±2
v
¢
· 1
8
¡
Duu±2u
+ 2Duv±u±v + Dvv±2
v
¢
case 3 . w(P0) 6= 0; R(P0) 6= 0.
Thus d(t) is C2 at t = 1 and d0(1) = 0, i.e.
d0(1) =
R(P0)
kR(P0)k
(R0u(P0)l cos µ + R0v(P0)l sin µ) + k
w(P0)
jw(P0)j
(w0u
(P0)l cos µ + w0v
(P0)l sin µ) = 0 (28)
Writing g(t) in its Taylor expansion at t = 1:
g(t) = g(1) + g0(1)(t ¡ 1) +
Z t
1
g00(s)(t ¡ s)ds (29)
Letting t = 0 gives
g(0) = g(1) ¡ g0(1) +
Z 1
0
g00(s)sds
or
0 = (R(A); w(A) = (R(P0); w(P0)) ¡ (R0u(P0); w0u
(P0))l cos µ ¡ (R0v(P0); w0v
(P0))l sin µ + I (30)
where
I =
Z 1
0
µ
(@2R(s)
@u2 ;
@2w(s)
@u2 )l2 cos2 µ + 2(@2R(s)
@u@v
;
@2w(s)
@u@v
)l2 cos µ sin µ + (@2R(s)
@v2 ;
@2w(s)
@v2 )l2 sin2 µ
¶
sds
Taking the dot product of equation (30) with ( R(P0)
kR(P0)k
; k w(P0)
jw(P0)j
) yields
0 = kR(P0)k + kjw(P0)j ¡ 0 ¡ 0 + I ¢ (
R(P0)
kR(P0)k
; k
w(P0)
jw(P0)j
)
Therefore we get
d(P0) · jI ¢ ( R(P0)
kR(P0)k
; k w(P0)
jw(P0)j
)j
·
R 1
0
£
(kR00 uu(s)k + kjw00 uu(s)j)±2u
=4 + 2(kR00 uvk + kjw00 uvj)±u±v=4 + (kR00 vvk + kjw00 vvj)±2
v=4
¤
sds
· 1
4 (Duu±2u
+ 2Duv±u±v + Dvv±2
v)
R 1
0 sds · 1
8 (Duu±2u
+ 2Duv±u±v + Dvv±2
v)
¢ 13
From Lemma 2 and Theorem 2, we can derive the required step sizes for rational surfaces:
Theorem 4. Let T ½ E2 be a right triangle with vertices A;B and C : B = A + (±u; 0); C = A + (0; ±v).
r(u; v) = R(u;v)
w(u;v) ; (w(u; v) > 0) : T¡ > E3 is a C2 rational surface. L(u; v) is a fractional linearly parameterized
triangle with r(A) = L(A); r(B) = L(B) and r(C) = L(C). Then sup
(u;v)2T kr(u; v) ¡ L(u; v)k · ² if
Duu±2u
+ 2Duv±u±v + Dvv±2
v · 8² inf
(u;v)2T
w(u; v) (31)
where
r ¸ sup
T kr(u; v)k
Duu =
8>><
>>:
sup
T
(kR00 uu(u; v)k + (r ¡ ²)jw00 uu(u; v)j); ² < r
sup
T kR00 uu(u; v)k; r· ² < 2r
0; 2r · ²
Duv =
8>><
>>:
sup
T
(kR00 uv(u; v)k + (r ¡ ²)jw00 uv(u; v)j); ² < r
sup
T kR00 uv(u; v)k; r· ² < 2r
0; 2r · ²
Dvv =
8>><
>>:
sup
T
(kR00 vv(u; v)k + (r ¡ ²)jw00 vv(u; v)j); ² < r
sup
T kR00 vv(u; v)k; r· ² < 2r
0; 2r · ²
(32)
Note that formula (31) has two unknowns ±u and ±v. It requires a second condition to determine them. In
case Duu = 0 or Dvv = 0, which implies the surface is fractional linear in the corresponding direction, we let
±u = 1 and get
±v =
p
D2
uv + 8Dvv² inffw(u; v)g ¡ Duv
Dvv
or let ±v = 1 and get
±u =
p
D2
uv + 8Duu² inffw(u; v)g ¡ Duv
Duu
Otherwise, set ±u=±v =
p
Dvv=Duu and we have
±u =
s
4Dvv² inffw(u; v)g
DuuDvv + DuvpDuuDvv
; ±v =
s
4Duu² inffw(u; v)g
DuuDvv + DuvpDuuDvv
This setting will minimize the number of triangles for tessellating the surface, as indicated in [Abi-Ezzi and
Wozny 1991]. A special case is Duu = Dvv = 0 and Duv 6= 0. We simply let ±u = ±v, and thus
±u = ±v =
s
4² inffw(u; v)g
Duv
For a rational B¶ezier surface expressed by
r(u; v) =
R(u; v)
w(u; v)
=
Pn
i=0
mP
j=0
wijPijBn
i (u)Bm
j (v)
Pn
i=0
mP
j=0
wijBn
i (u)Bm
j (v)
; u;v 2 [0; 1] (33)
14 ¢
R(u; v) and w(u; v) are in B¶ezier form, and the evaluation of Duu;Duv and Dvv is easy.
Duu = n(n ¡ 1) max
0·i·n¡2
0·j·m
fkwi+2;jPi+2;j ¡ 2wi+1;jPi+1;j + wi;jPi;jk + (r ¡ ²)jwi+2;j ¡ 2wi+1;j + wi;jjg
Duv = nm max
0·i·n¡1
0·j·m¡1
fkwi+1;j+1Pi+1;j+1 ¡ wi+1;jPi+1;j ¡ wi;j+1Pi;j+1 + wi;jPi;jk
+ (r ¡ ²)jwi+1;j+1 ¡ wi+1;j ¡ wi;j+1 + wi;jjg
Dvv = m(m ¡ 1) max
0·j·m¡2
0·i·n
fkwi;j+2Pi;j+2 ¡ 2wi;j+1Pi;j+1 + wi;jPi;jk + (r ¡ ²)jwi;j+2 ¡ 2wi;j+1 + wi;jjg
where r can be chosen to be r = max
0·i·n
0·j·m
kPijk.
It can be shown that the above estimation improves Cheng's result for surfaces. Also further improvements
may be possible by bounding the second order derivatives with the \weight points", or by transforming the
control points with a certain translation.
5. CONCLUSION
We have presented an approach to computing the tessellation step sizes for rational curves and surfaces. The
method is derived in the homogeneous coordinates, so it is numerically accessible. The new formulas make a
substantial improvement over previous known results. This is achieved by:
¢ e±ciently estimating the e®ect of perspective transformation;
¢ treating kR(t) ¡ Ln(t)k + kjw(t) ¡ Ld(t)j as a whole;
¢ performing more precise estimations.
Our approach is developed based on the parametric displacement of a curve or surface r from its fractional
linearly parameterized interpolant L. Future work could consider the geometric displacement, i.e., the perpendicular
distance of r from L, for the global step size estimation. Reparameterization of the curve does not
change the shape. In particular, the reparameterization by a MÄobius transformation [Farouki 1997] might yield
a larger step size. How to ¯nd an MÄobius transformation that optimizes step size warrants further study.
ACKNOWLEDGMENTS
The authors are supported in part by NSF grant number CCR-9712407. Jianmin Zheng is also supported by
Doctoral Science Foundation of China.
REFERENCES
Abi-Ezzi, S. S. and Wozny, M. J. 1991. Factoring a homogeneous transformation for a more e±cient graphics pipeline.
Computers & Graphics 15, 2, 249{258.
Cheng, F. April 1992. Estimating subdivision depths for rational curves and surfaces. ACM Transactions on Graphics 11, 2,
140{151. ISSN 0730-0301.
de Boor, C. 1978. A Practical Guide to Splines. Springer.
Elber, G. 1996. Error bounded piecewise linear approximation of freeform surfaces. Computer-aided Design 28, 1, 51{57.
Farin, G. April 1993. Tighter convex hulls for rational B¶ezier curves. Computer Aided Geometric Design 10, 2, 123{126.
Farouki, R. 1997. Optimal parameterizations. Computer Aided Geometric Design 14, 153{168.
Filip, D., Magedson, R., and Markot, R. 1986. Surface algorithms using bounds on derivatives. Computer Aided Geometric
Design 3, 4, 295{311.
Floater, M. 1992. Derivatives of rational B¶ezier curves. Computer Aided Geometric Design 9, 3, 161{174.
Hermann, T. 1992. On a tolerance problem of parametric curves and surfaces. Computer Aided Geometric Design 9, 2,
109{118.
Klassen, R. V. July 1994. Exact integer hybrid subdivision and forward di®erencing of cubics. ACM Transactions on
Graphics 13, 3, 240{255. ISSN 0730-0301.
Lane, J. and Carpenter, L. November 1979. A generalized scan line algorithm for the computer display of parametrically
de¯ned surfaces. Computer Graphics and Image Processing 11, 290|297.
Lane, J. and Riesenfeld, R. 1980. A theoretical development for the computer generation and display of piecewise polynomial
surfaces. IEEE Trans. Pattern Analysis Machine Intell. 2, 1, 35{46.
Piegl, L. A. and Richard, A. M. 1995. Tessellating trimmed nurbs surfaces. Computer-aided Design 27, 1, 16|26.
¢ 15
Rappoport, A. October 1991. Rendering curves and surfaces with hybrid subdivision and forward di®erencing. ACM Transactions
on Graphics 10, 4, 323{341. ISSN 0730-0301.
Rockwood, A. P., Heaton, K., and Davis, T. July 1989. Real-time rendering of trimmed surfaces. Computer Graphics
(Proceedings of SIGGRAPH 89) 23, 3, 107{116. Held in Boston, Massachusetts.
Schaback, R. February 1993. Error estimates for approximations from control nets. Computer Aided Geometric Design 10, 1,
57{66.
Sheng, X. and Hirsch, B. E. 1992. Triangulation of trimmed surfaces in parametric space. Computer-Aided Design 14,
437{444.
Tiller, W. 1992. Knot removal algorithms for nurbs curves and surfaces. Computer-Aided Design 24, 445{453.
Tookey, R. and Cripps, R. 1997. Improved surface bounds based on derivatives. Computer Aided Geometric Design 14, 8,
787{791. ISSN 0167-8396.
Wang, G.-J., Sederberg, T., and Saito, T. 1997. Partial derivatives of rational B¶ezier surfaces. Computer Aided Geometric
Design 14, 4, 377{381.
Wang, G.-Z. 1984. The subdivision method for ¯nding the intersection between two B¶ezier curves or surfaces. Zhejiang
University Journal: Special Issue on Computational Geometry 14, 108{119.