A sweeping surface S(u; v), where (u; v) 2 u0; u1] v0; v1] ..... We require the calling function to specify u0 and some vector v not in the direction of the tangent ...
Computing and approximating sweeping surfaces based on rotation minimizing frames Tim Poston
Shiaofen Fang
Wayne Lawton
CIeMed, Institute of Systems Science
National University of Singapore, Singapore 0511
ABSTRACT The sweeping surface, or generalized cylinder, is important in graphics, computer vision and elsewhere. A rotation minimizing frame, avoid the unnecessary twist given it by the Frenet frame. We give a quadratically convergent algorithm for computing it along a 3D trajectory curve, with an error bound. This is then used to create NURBS-approximated sweeping surfaces, to aid its incorporation into existing geometric modeling systems.
Keywords { surface modeling, sweeping, approximation, NURBS.
1 Introduction Sweeping a 2D cross-section curve C along a 3D trajectory curve A (perhaps changing C's size and shape en route), is a powerful, intuitive way to represent and construct 3D geometric shapes in computer graphics, animation, CAD/CAM and computer vision2,7,8,10. The Frenet frame, extensively used in the graphics and CAD/CAM communities in de ning it1,2,6,10, creates needless twist, avoided (Section 2) is by the rotation minimizing frame; originally introduced by Shani and Ballard8, and characterized as an ODE solution by Klok6. This slides the cross-section along the trajectory without rotation within the cross-section plane. It has not been widely accepted in graphics or CAD/CAM, partly because it has no closed formula for the frame at arbitrary points of A. Section 3 derives an ecient, quadratically convergent numerical scheme to nd it, with an error bound. Section 4 treats the NURBS surface approximation of the sweeping surface analogous to that for a Frenet-based one1 .
2 SWEEPING SURFACE AND FRENET FRAMES A sweeping surface S(u; v), where (u; v) 2 [u0; u1] [v0 ; v1], based on a 3D trajectory curve A(u) with unit tangent vector T(u), is de ned as S(u; v) = A(u) + c1(v)X(u) + c2(v)Y(u) , where (1)
F (u) = (X(u); Y(u); T(u)) is a dynamic orthonormal frame along the trajectory.
c(v) = (c1 (v); c2 (v)) is a 2D cross-section curve de ned in the abstract XY plane, mapped successively to each plane (T(u))? normal to the trajectory at A(u). h(u) is a scalar function that changes the size of the cross-section along A.
Figure 1. (a) The ribbon of normals N along c(s) = (s + s3 ; 2s; 2s4 + s), with a quartic term for spatial generality. (A cubic curve with an in ection is necessarily planar.) Any sweeping surface based on N is discontinuous. (b) Continuous but violently changing N along c(s) = (s + s3 ; 2s + s2=10; 2s4 + s).
If A is at least G1 continuous, T is uniquely de ned at each A(u), but the X and Y axes can rotate about T in the plane (A(u))? , producing di erent S(u; v), which thus depends strongly on the choice of F . Parametrizing A(u) by arc length s : a(s) = A(u), where u0(s) = 1=kA0(u(s))k, the Frenet frame FFr = (N(s); B(s); T(s)) is 0 (s) T(s) = a0(s) ; N(s) = kTT0(s) k ; B(s) = T(s) N(s) : For a C 2 (hence G2) trajectory curve with non-vanishing curvature, this dynamic frame is de ned and computable in closed form. It has been the most popular choice in the geometric modeling and graphics communities, in de ning sweeping surfaces. It is an unsatisfactory choice, for two reasons.
?
Figure 2. Let At(u) = t(cos u; sin u; 0) + (1 ? t) cos u2 ; sin u2 ; sin u2 ; gures (a,b,c) correspond to t = 0:7; 0:25; 0:15: The curve At (u) + Nt (u) in (a) has linking number l = 1 with At(u); the corresponding curve in (c) is unlinked (l = 0) with At. The integer value of l is continuous|hence locally constant|on the space of in ection-free G2 loops (shrinking N if necessary to avoid a loop piercing its ribbon of normals), so any smooth change At (u) from A0 (u) to A1 (u) must include an in ected At. In an animation sequence, in ection can be topologically unavoidable.
The trajectory A(u) may not be C 2 everywhere: C 1 and even G1 curves are common in graphics and modeling:
FFr is not even de ned where a curve is less than G2. Second, FFr has an intrinsic twist, forcing S(u; v) to turn
(unnecessarily and unnaturally) along the trajectory. This is measured by the torsion (s) = (dN=ds) B of the curve. For small a00 (for instance, on a cubic spline with a000 constant and a00 passing near 0), both a00(s) 0 00 000 ? ka00(s)k0 a00(s) 0 N (s) = ka00(s)k = ? ka (s)k a (s) ka00(s)k2
and the torsion (s) can become arbitrarily large, In the worst case ka00(s)k vanishes to give a point of in ection, where both N and B change signs, causing weird surface behavior. This can happen for seemingly well de ned C 2 curves (Fig.1). It is not enough to say6 that \for a `generic' non-planar segment the curvature is in fact non-zero." Almost any small perturbation removes a point of in ection from a space curve, but (Fig.1b) the discontinuity merely becomes a large continuous lurch; and the genericity result fails for animated curves (Fig.2).
3 LEAST ROTATION FRAME
Figure 3. (a) If we parameterize the outer curve|two arcs of unequal radius|by arc length, it is C 1 ; but equal spacing along it corresponds to di erent spacings along the two arcs that we nd by going unit distance along the inward normals, giving a G1 but ? 1 2 2 2 not C sweeping curve. (b) The C piecewise cubic A(u) = u; u ; u (u + juj) has discontinuous torsion, generating a curve of Frenet normals which changes direction at u = 0; swept tubes are not even G1 .
3.1 Parallel transport and least rotation frame The condition that the (X; Y) not twist around T is clearest via the curve that T(u) traces on the unit sphere S 2 . Along this, we require vectors in the tangent plane PT(u) S 2 to rotate by the minimum necessary to remain in it. This, equivalent to the non-twisting we require for (X; Y) at A(u), is the criterion for parallel transport3 of tangent vectors and tangent spaces along curves in a surface in Rn. It can be expressed, for a vector rolling with PT(u)S 2 , by the conditions (i) that T = 0 and (ii) that 0 (u) ? PT(u) S 2 so that 0 (u) is a scalar multiple (say, by (u) ) of T(u). Substituting this in ( T)0 = 0 from (i) xes = ?T0 =T T ?T0 =1, giving 0 (u) = ?
?T0(u) (u) T(u) ;
(2)
the ODE for parallel transport along the curve T(u) in S 2 : linear in , with non-constant coecients. Its solutions 1(u), 2(u) for an arbitrary initial orthonormal basis 1 (u0), 2 (u0 ) for PT(u0) S 2 provide an orthonormal basis for each PT(u) S 2 , and hence for each normal plane at A(u). By (2), if A(u) is not C 2 , then (u) is not C 1. This is inevitable, even in 2D, where all reasonable frames coincide (Fig.3a). However, since for a G1 curve the direction of T is continuous, so also by (2) is the direction of 0(u); that is, is G1. Given a G1 sweeping curve, this results in a G1 surface. In contrast, the Frenet frame's use of the second derivative makes it discontinuous where A is not G2; and even for B-splines the Frenet frame does not give G1 curves (Fig.3b).
3.2 Computing the least rotation frame This frame requires an arbitrary initial choice of ( 1 (u0); 2(u0 ). Finding it at another A(u) entails solving (2), possible in closed form only for a few special curves. However, solving numerically for the frames at u1; u2; u3; : : :
costs no more than the adaptive steps required by the Frenet frame in Fig 2. One could solve (2) by Euler, Runge-Kutta, etc., methods: we use a scheme tailored to the relevant geometry. The unique smallest-angle rotation that carries T(u) into T(u + u) 6= T(u) is given by Ru;u
2 a2 ab ac 3 2 cos ?c b 3 1 ? cos = 4 c cos ?a 5 + a2 + b2 + c2 4 ab b2 bc 5 ac bc c2 ?b a cos
(3)
where Q (a; b; c) = T(u) T(u + u); cos = T(u) T(u + u). If u = (u1 ? u0)=(n + 1), the product Rn = ni=0 Ru0 +iu;u is a rotation carrying T(u0 ) exactly into T(u1 ) and N(u0) into N(u1 ), so that any curve-discretization error is purely a mis-rotation within N. To compute Rn uses 27n oating multiplies8; nding 1(u1 ) = Ru0 +n; ( (Ru0; ( 1 (u0 ))) ) needs only 9n, after which 2 can be found with 6n more as T 1 . However, the rotational viewpoint makes it simpler to analyse the discretization error. For any surface , parallel transport around a closed curve t in gives3 a net rotation of the initial tangent plane equal to the integral of the Gaussian curvature K over the region R enclosed by t. Since for S 2 K 1, this reduces to the area of R. In transport around a wholly right-angled triangle on S 2 , with corners at the N.Pole and 90 apart on the Equator, a net turn of =2 is geometrically obvious. So is the 2 0 case of a plane circular A with the corresponding t a great circle. Thus we can analyse errors in solving (3) by considering the area between the true curve T(u) and the curve assumed in a numerical approximation. Assume T(u1 ) = T(u0 ) , or make it so by adding to A a circular arc, on which (3) is exact transport. Then Rn rotates PT(u0 ) by the area on S 2 contained within the geodesic polygon T(u0 )T(u0 +u) T(u0 +nu)T(u0), not by that enclosed by T .
Figure 4. The error
P
n
(shaded areas) in taking the area inside T as that inside a polygon shrinks with 1=n2 for A smooth.
For a(s) a helix with curvature and torsion , the error in approximately solving (2) from a(se) to a(sb) of helix by (3) is easily seen (using the relation between enclosed area and rotation by parallel transport) as 2 3 E( ; ; ) = ? 2 arcsin p sin2 2 + ? 12 ; 1 ? sin p2 2 where = + ; = ; = ; = bs ? es ; = 2
Only for helices are and constant; however, (se)2 (se)3 =12 is the exact cubic term in the expansion in of the error for an arbitrary curve, giving a direct way to determine the number n of equal length arcs in which to divide a curve of length L to ensure that total accumulated angular error is less than " : n
s
j j2L3 12"
(4)
For a cubic polynomial segment, with a(s) = c0 + s c1 + s2 =2 c2 + s3 =3 c3 ; a0 (s) = c1 + s c2 + s2 c3 , we need not explicitly compute the non-polynomial quantities and . For any curve, consider a small piece of the
curve like that from T(se) to T(sb) in Fig.4. Setting q = t(se) t(sb), the shaded area A is bounded by the product of arcsinkqk kqk and the maximum (over ) of kqqk T(s), and thus by
(T (e) T (b)) a e+2b
s s max j(T (es) T (bs)) a (s)j = (numerator (5) max j ( T ( e s) T (b s)) T(s)j quadratic) min k a ( s ) k min k a ( s ) k es+bs (a (es) a (bs)) a 2 c1 c2 ) c3 (es ? bs)3 = C (es ? bs)3 ; say, = 4(min (6) 3 min ka (s)k3 ka (s)k giving an explicit bound C(s1 ? s0 )3 =n2 for the total error. The local canonical form for the curve a(s) on page 27 of do Carmo3 yields an exact expression (c1 c2) c3 = ? (s)2 (s)ka0 (s)k6 =12 for the triple product, which 0
0
0
0
0
s s
0
0
0
0
implies that the general upper bound (6) is four times the actual error for the helix.
3.3 Initial frames and closed curves The circle's worth of possible minimal rotation frame systems f 1 (u); 2 (u)g for A(u) corresponds to the 2 of angles through which a given f 1 (u0 ); 2 (u0)g at some A(u0) can be rotated. One angle is enough to change any least-rotation frame to any other, so a single angle should apparently suce as a label. Ideally, this should be translation invariant; the frame corresponding to for curve a(u) should have the same 1 (u), 2 (u) as the frame for k+a(u), with k any constant vector. However, suppose we had such a scheme: then, for every curve we could choose the frame system corresponding to = 0 . In particular, for every straight A we would have a vector 1 orthogonal to the unit vector t along L. By translation invariance, 1 would depend only on t. That is to say, for every point t on the sphere S 2 we would get a vector 1 in the tangent plane at t to S 2 . But the well-known Hairy Ball Theorem (Spanier9 , Lemma 11, x 4.7) states that such a system of tangent vectors cannot be continuous. Thus, if we add continuity to the requirement of translation invariance, no single-angle frame speci cation can exist. We require the calling function to specify u0 and some vector v not in the direction of the tangent vector t(u0); if the surface is known never to be vertical, v = (0; 0; 1) suces. De ning a vector within the initialization as the component ? ( t)t of in the plane normal to t, we can set 1(u0 ) = =kk, 2(u0 ) = 1 (u0) t.
Figure 5. A closed-curve sweeping surface with mismatching ends.
For a closed curve A, a minimal-rotation 1 (u), 2 (u) does not give 1 (u1), 2(u1 ) equal to 1 (u0), 2 (u0), unless the area enclosed by T(u) is some 2n (Fig.5). The Frenet frame must match, at G2 non-in ection points; but it is precisely for closed curves that it gives most trouble (Fig.2). there is often an unavoidable topological obstruction to the existence of any continuous movement between given curves, that would keep the Frenet frame always well de ned. Either Frenet or pure minimal-rotation frames thus involve us in orientation discontinuity around an evolving loop. Starting from a minimal-rotation (X(u); Y(u)), we can at least keep the rotation under control. Any moving
frame for the (A(u))? can be written, for some function w(u), as
Fe = Xe (u); Ye (u) = (X(u) cos w + Y(u) sin w; ?X(u) sin w + bY (u) cos w ) :
?
(7)
Fixing = arctan2 X(u0) Y(u1 ); X(u0) X(u1 ) and setting w(u) = (u ? u0 ) =(u1 ? u0 ) gives the unique continuous frame that matches F at A(u0 ) and minimizes the twist energy integral E(Fe) = =
Z s1 s0
2 Xe 0 Ye ds
(8)
Among rotating frames continuous around A, it is thus minimally rotating as measured by (8): it is most geometrically natural, and the twist looks most evenly distributed, when the parameter u is arc length. This solves the continuity problem for a single curve, but because of the necessary discontinuity of arctan2( ) it cannot be applied independently to an evolving sequence of sweep trajectories. Each time the total area in S 2 enclosed (with multiplicities) by T(u) passed through an odd multiple of , the value ? of would jump by 2, between . This would not spoil the continuity at the join, but the value of w (u0 + u1)=2 would jump by , turning that section of the surface exactly around between one image and the next. It is thus necessary to enforce continuity by adding a varying multiple 2n of a full rotation: for each successive curve Ai+1 , supposed close to Ai , nd for the minimal-rotation frame (Xi+1 ; Yi+1 ) i+1 := arctan2 (Xi+1 (u0) Yi+1(u1 ); Xi+1(u0 ) Xi+1(u1 )) if ( i+1 ? i) > ; n = n ? 1 if ( i+1 ? i) < ; n = n + 1 : The nal result, with a shared-out twist of last + 2n, can be quite strongly twisted (Fig.6), but this is topologically inevitable. Certain curves are indeed twistier than others. If they are modelled in an easily-bent material (such as a rubber band or a catheter), their twist-energy-minimising `preference' for positions that match the twist in their relaxed state strongly a ects their equilibrium shapes.
Figure 6. Two `least twist for the family' closed sweeping surfaces with ends. The second surface inherits an extra 2 twist.
4 SURFACE APPROXIMATION & SYSTEM INTEGRATION A sweeping surface (1) cannot, in general, be represented as a tensor-product spline surface such as NURBS, the CAD industry standard for surfaces. An approximation is necessary, so as to incorporate the sweeping surface into existing NURBS-based modeling systems for applications. The sweeping surface derived above forms part of a geometric modeler4 developed at CIeMed1, which uniformly uses a special C 1 B spline representation (a subclass of NURBS) for curve and surface models. A sweeping surface also needs to be approximated by such
1 The Centre for Information-Enhanced Medicine, a research partnership between the Institute of Systems Science of the National University of Singapore and the Medical School of Johns Hopkins University
a model. We brie y outline the curve and surface representation of the CIeMed geometric modeler, and its approach to approximating sweeping surfaces with C 1 B spline surfaces.
4.1 Curve and surface representation A C 1 B spline curve (t), i.e., a non-uniform cubic B spline curve5 with knots t1 ; : : :; tn, of multiplicity 2 for all inner 4 and four for the ends (Fig.7). has a computationally ecient form: for a given t 2 [ti; ti+1],
(t) = bi (t ? ti ) = bi (s) = [s3; s2 ; s; 1]Mi( )[C2i?1; C2i; C2i+1; C2i+2]T (9) where s = t ? ti 2 [0; li], li = ti+1 ? ti , and we de ne t0 = t1 ; tn+1 = tn.
0? 1 2l ?1 +3l ? 3l +2l +1 BB l2(l ?1 +l ) l33((ll ?1 ++2l l) ) l3 (l +l +1 ) 3 B l (l ?31 +l ) ? l2 (l??11 +l ) l2 Mi ( ) = B BB 3 0 B@ ? l ?13+l (l ?1 +l ) l i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
l2i (li +li+1 )
i
i
0 0
0
0
i
?1 (li?1 +li ) i
li?1 +li
1
i
i
i
i
li
i
i
1 CC CC CC CA
A C 1 B spline surface is the tensor product of C 1 curves. Let surface S (u; v) have U knots = fu1; : : :; ung, V knots = fv1 ; : : :; vm g and control points fPij g; (i = 1 2n; j = 1 2m). Its B Spline form (the tensor product of (9) applied to u and v) is, for each [ui; ui+1] [vj ; vj +1]: S (u; v) = bij (u ? ui; v ? vj ) = bij (s; t) = [s3 ; s2; s; 1]Mi()BMjT ()[t3; t2; t; 1]T (10) where s = u ? ui ; t = v ? vj , and 0C 1 2i?1;2j ?1 C2i?1;2j C2i?1;2j +1 C2i?1;2j +2 B C2i;2j?1 C2i;2j C2i;2j+1 C2i;2j+2 CC B=B @ C2i+1;2j?1 C2i+1;2j C2i+1;2j+1 C2i+1;2j+2 A C2i+2;2j ?1 C2i+2;2j C2i+2;2j +1 C2i+2;2j +2 P2
P3
P6
P7
P1 P10 P4
P5 P8
P9
Figure 7. A C 1 curve (t) with ten control points and ve knots.
4.2 Sweeping surface construction To get a C 1 B Spline surface approximating (to within total angular error ") a numerical S(u; v) given as
Cross section: a set of contour points on the cross section contour curve, centred on (0; 0). Trajectory: a set of centreline points on the trajectory curve. Scaling factors at each centreline point, that multiply the contour points there,
we 1. Fit a C 1 curve A(u) to t the centreline points.
Fit the contour points by a cross section curve c with 2D control points: fP1;1; P1;2; : : :; P1;2mg. Interpolate the scaling factors between the centreline points. Use " and (4) to x the number n of sampling points, on the trajectory for the rotation minimizing frames. Compute the frames, using (3). Transform c to the cross sections at all n sampling points, scaled appropriately, to generate a cross section curves, with control points fPi;1; Pi;2; : : :; Pi;2mg for each i = 1; : : :; n, 6. Fit fP1;j ; P2;j ; : : :; Pn;j g, for j = 1; : : :; 2m, to a C 1 B spline curve with control points fC1j ; C2j ; : : :; C2n;j g. Then (9) with fCij g (i = 1; : : :; 2n; j = 1; : : :; 2m) interpolates the cross section curves. 2. 3. 4. 5.
In conclusion, this achieves the results claimed above; an ecient, untwisted sweeping surface in industrystandard NURBS form.
5 REFERENCES [1] M. Bloomenthal and R. Riesenfeld. Approximation of sweep surfaces by tensor product NURBS. In Proc. Curves and Surfaces in Computer Vision and Graphics, SPIE 1610, Boston, Massachusetts, November 1991. [2] W. F. Bronsvoort and F. Klok. Ray tracing generalized cylinders. ACM Trans. on Graphics, 4(4):291{302, October 1985. [3] Manfredo P. do Carmo. Di erential Geometry of Curves and Surfaces. Prentice-Hall, Inc., Englewood Cli s, New Jersey, 1976. [4] Shiaofen Fang, Wieslaw Nowinski, Lakshmipathy Jagannathan, and Rajagopalan Srinivasan. Geometric model reconstruction from cross sections for branched structures in medical application. June 1995. [5] Gerald Farin. Curves and Surfaces for CAGD, A Practical Guide. Academic Press, 3rd edition, 1993. [6] F. Klok. Two moving coordinate frames for sweeping along a 3D trajectory. Computer Aided Geometric Design, 3:217{229, 1986. [7] S. K. Semwal and J. Hallauer. Biomedical modeling: Implementing line-of-action algorithm for human muscles and bones using generalized cylinder. Computers and Graphics, 18(1):105{112, 1994. [8] U. Shani and D. Ballard. Splines as embeddings for generalized cylinders. Computer Vision, Graphics, and Image Processing, 27:129{156, 1984. [9] Edwin H. Spanier. Algebraic Topology. McGraw-Hill Book Company, 1966. [10] W. P. Wang and K. K. Wang. Geometric modeling for swept volume of moving solids. IEEE Computer Graphics and Application, 6:8{17, 1986.
