Computing Curve{Surface Incidence Constraints Eciently Paul Michalik Paul.Michalik@prakinf.tu-ilmenau.de Beat Bruderliny bdb@prakinf.tu-ilmenau.de Technical University of Ilmenau Institute for Applied Computer Science Computer Graphics Program Postfach 100565 D-98684 Ilmenau Germany January 19, 1999 Abstract We investigate methods of using constraint{based modeling in a free{form curve and surface enviroment. In this work we concentrate on a problem of maintaining the curve-surface incidence relation while the curve is edited. We formulate the relation between the degrees of freedom (DOF) and parameters (control points of the surface and curve resp.) as an explicit functional prescription. We show how the polynomial composition algorithm based on blossoming and a data structures for ecient storage of intermediate results, signi cantly speeds up the computation. We also discuss how to insert new DOFs, if the incidence condition cannot be satis ed, and sketch how a common solution space for several constraints of this type can be found. Keywords: Incidence, Constraints, Free{form surfaces and curves, Algorithm, Blossoming, B{spline algebra. 1 Introduction The relational modeling paradigm has achieved much progress in recent years especially in 2{D. Editing ge y Tel. ++49 3677 691211 Tel. ++49 3677 692785 ometric or other data while some set of constraints is satis ed has shown to be a very useful and intuitive method for design and modeling. However, keeping the constraints satis ed even in conceptually simple cases such as constrained distances between points or angles between lines is a difcult task especially in 3{D. The situation becomes even more complicated if geometries non{linear in their nature, e.g. polynomial curves or surfaces are involved. Sets of highly non{linear equations arise immediately and many unpredictable solutions might exist. Interactive response of the system, the most popular feature of constraint solving systems, becomes seemingly impossible. Two main approaches for free{form geometries to circumvent the diculties have been proposed. Either only polynomials of low degree and/or of a special type are used to make solving of non{linear systems easier [7], or as in [14, 1, 5] the non{linear part is kept constant which makes the problem linear. Both approaches, of course, have their limitations, but in practical applications they seem to work quite well. The second approach is one of the most wide spread, because of the numerical ecency (only system(s) of linear equations have to be solved). It only depends on application what kind of constraints (formulated as linear equation in the DOFs) are needed. The basic mathematics regarding this approach will be reviewed brie y at the beginning of section 3. A more powerful but also more complicated method is formulating the constraints as non{linear system of equations, but at the cost of additional diculties mentioned above. Nevertheless in e.g. [12] this path has been pursued to constrain the global curvature of a free{form surface. 2 The problem Our work solves a practical reverse engeneering problem formulated by developers of CAD components. Consider the scheme in Fig. 1: A 3{D scanner measures unstructured points which are interpolated by bi{polynomial parametric surface patches. The surface model generated this way is usually of good global quality but in some details does not meet designer's intent. Recomputing the model starting with point data is not a well controllable process. Often some heuristics have to be used (e.g. when de ning the parametrization of the patch), which makes the outcome of the interpolation hard to predict. Experience shows, that from designer's point of view, it is often more plausible to change the model locally, inside restricted areas, near certain points or characteristic curves. We suppose the user marks such points or curves on the surface which should be edited to meet certain design criteria. If a parametrization of these curves is de ned, these curves can be represented exactly as parametric curves: C (t) = S (u(t); v(t)). These curves may be modi ed in 3{D while the initial incidence relation is kept. Thus, the design parameters of the curves de ne parameters of the model. The surface coecients are the degrees of freedom. 3 The method Our method is inspired by [1]. In this article the surface curves are kept xed in space and the surface is deformed by applying forces at marked points, while keeping the curve{surface incidence. In our application the 3{D curve on the surface is seen as the independent (known) function f (t) = C (t) and the uv{ curve inserted into the surface equation as the model function F(t) = S (u(t); v(t)). The problem is formulated as a continuous function approximation. The independent function f is to be expressed in terms of the model function F, such that norm (f; F) = kF ? f k stays minimal. In [1] an error functional is introduced which is to be minimized over the entire domain of the Surface model B−spline 3−D Scanner Real Object Define surface curve Manipulate the surface curves, the incidence will be satisfied by appropriate change of the surfaces Figure 1: Design process curve: Z a b (t) dt = Z b1 a 2 (1) [S (u(t); v(t)) ? C (t)]2 dt = min In context of polynomial curves and bi{polynomial surfaces de ned on bases ,  respectively, the formulae become: S (u(t); v(t)) = and m X n X xij i (u(t))j (v(t)) i=0 j =0 C (t) = q X yk p (t) p=0 The control points of the surface xij are chosen as the parameters to be optimized. After applying the (t) global minima conditions @ @x = 0; 8i; j we get a system of m
n equations in m
n unknowns each of them having a form: ij m X n X i=0 j =0 xij z Z bh a a(8 i;j ) ;kl }| ?
{ ?
i ij t kl t dt ? q X p=0 yp Z bh |a (2) ?
i p (t) kl t dt = 0 {z q kl } where ?
?
ij = i u(t) j v(t)
?
? kl = k u(t) l v(t) Since the parametrization is assumed to be xed, xij and yp may be placed in front of the integral. The equations are linear in the unknowns xij . The condition 1 initially satis es C (t) = S (u(t); v(t)), but it also enforces the best possible approximation of this state if C (t) is changed. In contrast to [1] we are looking for an expression: F : yp ?! xij After separating the knowns and unknowns in equation 2 we obtain an linear equation system in matrix form, using the braced terms a( i;j);kl and qkl from eq. 2: A
x=Q
y yielding after some Gaussian elimination steps (carried out simultaneously on A and Q) and solving for x: A
x = Q
y =) x = Q
y ? N
p (3) Full elimination is in general not possible (and not wanted!), the surface retains some unconstrained degrees of freedom, which occur as singularities in the matrix A. In eq. 3 N is a column matrix representing the nullspace of A and p is a vector of undetermined degrees of freedom { the free variables. 8 0 0 The method, so far, su ers from several complications. What follows is a short discussion of those, sections 5{ 7 deal with proposals how to remedy at least some of them. Before solving equation 2, entries of the matrices A and Q must be computed. The coecients for one of the m
n rows (di erentiation with respect to kl{th unknown) are the terms a( i;j);kl and qkl from equation 2. This means one has to substitute the parametric representa?
tion of the curve u(t); v(t) into the surface basis functions i (u)j (v), multiply the resulting polynomials and then evaluate the integrals. This has to be done (m
n)2 times. Without knowing speci c properties of the basis functions used, the time needed for computation grows beyond acceptable limits. 8 In many cases the degrees of freedom xij of the model function (surface) do not suce for each constellation of yp . This means the surface can not react appropriately if the curve is pulled \too far away" from its initial position. As a consequence, either the error (t) becomes too large, or the problem is ill{conditionded, the existing DOFs are very sensitive to changes and destroy the \good" overall shape of the model, laboriously derived in the reverse engeneering process. The authors of [1] already hinted this problem, but we would like to discuss it in more depth.  Multiple contraints If multiple constraints of some type are de ned, the maintenance of the system becomes quite difcult. In general a graph structure with nodes beeing the constrained objects and arcs the relations among them arises from such problems (see e.g. [3]). Nevertheless, we will try to sketch how this special types of constraints can be satis ed simultaneously in section 7. 0 4 Drawbacks of the method  Ecency  Lack of DOFs 5 Eciency If Bezier or B{spline basis functions are used, the computation can be made much more ecient. In the following, when we talk about basis function, we mean d , B{spline basis, and denote them by the usual Bi; lower index: index, knot sequence, upper index: degree. We will show how ecient B{splines product computations [11, 4], blossom based polynomial composition [2, 9], local support properties of B{splines basis functions (see e.g. [6, 15]) and the combinatorial issues attached can be used to improve the run{time behavior for evaluation of the elements of A and Q. The terms a, q in eq. 2 contain composition of the uv{curve with the 2{dimensional basis functions of the surface. A surface basis function itself can be seen as special case of a bi{polynomial surface with all coecient beeing zero except the one with indices ij equal to one. This point of view allows considerable acceleration of the computation and is also more robust { the zero coecients do not need to be taken into account. In [2] an algorithm for polynomial composition using the blossoming technique [13] is introduced, further analysis of the algorithm can be found in [9]. In those articles multiane properties of a blossom are used to formulate an ecient algorithm to compose Bezier polynomials of appropriate dimension. Com- position of B{splines is done via basis conversion (B{ spline to Bezier basis) and a reconstruction of the B{ spline form of the result. This is not very advantageous in our application, since the extraction of the knowns and unknowns in formula 2 would become quite complicated { we perform the composition directly on B{splines. ?
A polynomial composition C (t) = S u(t); v(t) of a k with p B{spline curve u(t); v(t) de ned on basis Bp; control points (ui ; vi ), and a B{spline surface S (u; v) of degree du , dv in either directions can be written in terms of blossoms as: ?
S u(t); v(t) = p X i0 p X


p X j0 idu


? (4) Bik0 ;


Bikdu ;
p X jdv Bjk0 ;


Bjkdv ;
f ui0 ;


; uidu ; vj0 ;


; vjdv
We skip the use of special notation for multiindices and hyperindices [2, 9] for the moment; i or j is just an integer symbol for which 8i; j 2 f0;


; p ? 1g holds. Suppose we are able to construct the products of the B{spline basis functions again as a B{spline curve, i.e. express it as a linear combination of the \product" B{spline basis functions (di erent than the original one!). Then doing all the iterations and collecting the coecients of the appropriate \product" basis function results in the control points of the B{spline representation of the composed curve. Naive evaluation of the above formulae would of course be wrong. Polynomials on B{spline basis can not be handled like \normal" ones, some further work has to be done. The interesting thing is, as authors of [2] already pointed out, that the terms in eq. 4 only depend on the outer function via blossoms. We will show, how the blossoms can be computed in a coordinate{independent way and that the intermediate results for the blossoms and products have to be computed and stored only once(!) and then can be reused for composing each \basis surface". 5.1 Indexing notation It would be extremely inecient to iterate over all the i and j in eq. 4. Just consider for example if k = 2; p = 11; du = dv = 3 the pdu +dv = 1; 771; 561 necessary iterations. In [2] a di erent indexing scheme was proposed for Bezier geometries. It is applicable (though with slight changes) to B{splines directly. Because of the combinatorial invariance of the permutations of i and j values, it suces to iterate over the \statistics" of them: Suppose there is an ordered set of integers I : fi0 ;


; ip?1g from which d; (d  p) have to be chosen. Assume that the permutations of entries do not matter (as they don't in case of the products) and that each k + 1 consecutive entries build a \segment", such that only entries inside this segment can be chosen at once (this also applies to the products for degree k, see cond. 5 below). Then a multi{index for set I , bandwidth k +1, depth d and at segment s can be de ned as I ~{s;k;d = [# (is ) ; # (is+1 ) ;


; # (is+k )]. Furthermore P I for each ~{s;k;d the relation sr+=ks # (ir ) = d holds. An example: Suppose we have the following constellation (for d = du + dv = 5) of the indices iI and jI in formula 4: [0; 0; 1; 1; 3]. This combination can then be coded as [#(0); #(1); #(2); #(3)] = [2; 2; 0; 1]. One can imagine the data structure as an \scan array" of bandwidth k + 1 moving along all ordered entries in the set I . Authors of [2] de ne the relation ~{curr ! ~{succ, which allows to iterate (segment by segment) directly over these \statistical" vectors. 5.2 The products How to express a product of piecewise B{spline polynomials as a B{spline again has been discussed in [4] where linear interpolation of product values and/or conversion to Bezier basis was used. In [11] a direct algebraic algorithm was introduced, using a method similar to knot insertion algorithm. In both cases, a B{spline space (dennoted by Bld; ) is found, on which the product can be based. In [4] the basis functions of the product space are evaluated at the Greville abscisass of  which yields a banded square matrix, denoted by A . Values of the operands of the product (in our case the values of the basis functions in each product term) are also evaluated at the Greville abscisass and multiplied (yielding vector b ). The coecients of the product B{spline (x ) can then be found by solving A
x = b . Additionally, it is known a priori which of the entries of b are zero. They don't need to be computed at all and are marked in each b and x . 5.2.1 Interpolating the products Taking a closer look at eq. 4 tells us, that the matrix A only needs to be evaluated once! All products can be based on the same basis, the B{spline form of the products can then be computed as xi = A?1
bi . The values of b , of course, have to be computed separately for each product. i 5.2.2 The combinatorics LetQus aggregate the product terms in eq. 4 and call it B(di;j ) 2 (mn); . The question is, how many unique non{zero products of this kind exist? This is a simple exercise in combinatorics. The product does not depend on permutations of the indices. If p is the number of control points of a uv{B{spline ?
of degree k and d = du + dv , there will be k+dd?1 unique products necessary for the ? rst segment. For all the next segments there
are k+dd?2 products, which were already computed during processing the previous segment. Assuming p ? k + 1 overall segments for the uv{B{spline, we can then say about number of products: # Y   B(di;j ) 2 (mn); =   k+d?1 + d     k+d?1 k+d?2 ? (p ? k) d d The set of all indices from which we choose the single d{vector reduces from p to each k + 1 (recall the \bandwidth" above?) immediate neighbors because of the local support property of the basis functions: Y 8 > > < (8i; j 2 ~{ ) ? = 0 max min (8i; j 2 ~| ) > k > : 6= 0  B(di;j ) 2 (mn); > otherwise (5) It can be seen, that the degree of the resulting curve is k
du + k
dv = k (du + dv ) (compare e.g. [8]). 5.3 The blossoms The blossom of a B{spline of degree d is only well dened if its arguments are all in d +1 neighbor intervals (see [13] where several kinds of \overloading" of blossoms arguments are profoundly discussed). We pick the most \strict" one { tame overloading. The arguments of the blossom are the control points of the uv-B{spline. It turns out that the above condition is NOT satis ed i the product terms are zero! One should not forget, that we have inserted knots into the uv{curve at the points of intersection with all knot lines, so that the structure of the uv{space and the parameter space of the uv{curve \agree" (s. [2]). We need to compute only blossoms for those permutations of arguments for which the products are non zero as well. In this case, they will always be well de ned even under the strict tame overloading. 5.3.1 Generalized basis functions As already mentioned, it would be of great advantage if we were able to compute the blossoms in a coordinate independent way. For this purpose, a slight variation of the Cox{deBoor algorithm is used. One only has to insert a di erent parameter value (each argument of the blossom) in each deBoor step. If the coecients arising in each step are collected underway, what we end up with, is the well known Knot insertion algorithm. The cocients are sometimes refered to as linear splines (e.g. in [15] or [11]) and denoted by . We, however, call them generalized basis functions. At this point, we can express a blossom value as a linear combination of control points: ?
fs u0 ;


; udu ; v0 ;


; vdv = +dv suX +du svX i=su j =sv pij i j For further discussion and algorithm see [10]. 5.3.2 The combinatorics As already mentioned, the total number of blossoms needed, grows with the number of non{zero products. However, a blossom of a two{dimensional simploid (bi{polynomial surface) is only invariant under permutation of it's arguments inside either parametric direction. The number of unique blossoms in, for instance, the u{direction can then be expressed as:   # f (u1 ;


; ud ;


) =   k + du ? 1 + u du  k + du ? 1 du    ? k + ddu ? 2 (p ? k) u As in the case of the products of unique
? the number index vectors decreases by k+dd ?2 in all segments but the rst one. It is also possible (as discussed in [2]) to reduce the number of evaluations of the {factors, by reusing those computed on further stages for the same combinations of arguments. u u 5.4 Multi{index tree Upon analyzing the mathematics so far, we can think about the data structure in which the products and blossoms (the coecients) will be stored. We must be able to insert and retrieve the precomputed entries fast, if possible without any kind of preprocessing, sorting or ordering. [ 2, 1, 0, 1 ] ≡ [ 2, 2, 3, 5 ] 2 Because the parametrization of the resulting curve S (u(t); v(t) has to re ect the structure of the domain space of the surface and the uv{curve, it also has to be re ned at parameters of the intersection points with those new knot lines. S( u( t ), v( t ) ) 0 1 2 ... p 0 1 2 ... p 0 1 2 ... p d ... ... 0 1 2 ... p ... ... 2 3 4 ... S( u, v ) ... 4 5 6 ... u( t ), v( t ) tu e Figure 2: Multi{index tree access scheme We propose a tree structure organized as in Fig. 2. Each node (except of a leaf node) consists of an array of dimension p with links to it's children nodes. As a query a multi{index ~{, as described above, is passed to the tree. Depending on the entries of the query, the tree is traversed, at each node the appropriate \turn" is taken. After d steps always a leaf node is reached, where the data (product, blossom coecients, or both) are stored. The depth of the tree never exceeds d. In Fig. 2 an example multi{index for d = 4; k = 3; s = 2 is shown. The solid arrows mark the \way" along the tree for this special entry. Using the algorithm described above and the tree structure, we are able to compute all products and all blossom coecients very quickly. For example computation/storage of all 196 unique products for k = 2; p = 11; du = dv = 3 takes ca. 0:15 seconds on a SGI{O2, which is an improvement factor of ca. 9000 of estimated time of naive approach as sketched in section 5. 6 Explosion of DOFs To get a better response of the surface (and a smaller error) for ill conditioned problems new DOFs have to be inserted. In case of a Bezier or B{spline basis this does not cause any problems. Insertion of new knot lines in the parameter space of the surface in either directions generates new control points (new DOFs). vi a tu b tv ui tv Figure 3: Insertion of new knot lines also re nes the curve The problem is that this of course generates new yp s (i.e. design parameters). The consequence is that though we have more DOFs, we also have more parameters for which again more DOFs may be required! This might cause a vicious circle in the design scheme. Fig. 3 demonstrates the situation. The way out is a hierarchical structure for the yp . Whenever the response of the surface is not satisfactory, insertion of new DOFs is initiated. The number of parameters stays constant from users point of view. Internally the \old" and the \new" y are stored in a hierarchical data structure (see Fig 4). The user is only allowed to operate on the original parameters (y), the new ones y are dependent on the y. 0 7 Multiple constraints In practice, the de nition of more than one constraint may be required. For example it might be very useful to x surface normals along the curve or the position of the boundary curves. Writing down the conditions for constrained normals or any di erential values is Acknowledgement y2 y’2 y0 This work was supported in part by a grant of the Ministry of Science and Culture of Thuringia (TMWFK) Germany. y’3 y3 y’0 y’1 y1 Figure 4: Hierarchical structure for design parameters a nice exercise in algebra, basically the same as for \simple" curve{surface incidence. In [5] some of those cases were realized for Bezier surfaces. Our current work goes in the following direction: Suppose each constraint generates a solution space as eq. 3. One can try to nd the intersection of the solution spaces for each constraint. Then conditions on the DOFs and parameters can be formulated, which enforce the satisfaction of de ned constraints whenever possible. An analogy from 3-D space: Imagine two solution spaces - planes in parametric form. Their line of intersection (if any) would be the space where common solution exists. 8 Conclusion We have combined and extended several methods to solve a practical problem of de ning free form surfaces from incident curves. We have shown a method how to speed up the involved computations and proposed how new DOFs should be inserted without getting caught in a vicious circle using special data structures and the blossoming technique. Future work will look for methods where the surface has to be re ned (at what position the new DOFs should be inserted). We are also working on the storage/management/solving of several constraints of this (and similar) type for one and more surfaces. References [1] G. Celniker and W. Welch. Linear constraints for deformable b-spline surfaces. Computer Graphics, 25(2):171{174, March 1992. [2] T. DeRose, R. Goldman, H. Hagen, and S. Mann. Functional composition algorithm via blossoming. ACM Transactions on Graphics, 0(2), April 1993. [3] U. Doering, P. Michalik, and B. Bruderlin. A constraint{based shape modeling system. In Geometric Constraint Solving & Applications. Springer Verlag, June 1998. [4] G. Elber. Free form surface analysis using a hybrid of symbolic and numerical computations. PhD thesis, University of Utah, 1992. [5] G. Elber and E. Cohen. Filleting and rounding using trimmed tensor product surfaces. In The fourth ACM/IEEE Symposium on Solid Modeling and Applications, pages 201{216, May 1997. [6] G. Farin. Curves and Surfaces for CAGD. Computer Science and Scienti c Computing. Academic Press, Inc., 3. edition, 1992. [7] C. M. Ho mann and J. Peters. Tschirnhaus cubics analyzed for a constraint solver. In M. Dahlen, T. Lyche, and L. L. Schumaker, editors, Mathematical Methods in Computer Aided Geometric Design, volume III. Vanderbilt Press, 1995. [8] J. Hoschek and D. Lasser. Fundamentals of Computer Aided Geometric Design. A K Peters, Ltd., 1989. [9] S. Mann and W. Liu. An analysis of polynomial composition algorithms. Technical Report CS{95{ 24, University of Waterloo, Computer Science Department, 1995. [10] P. Michalik. Objektorientierte Implementierung von parametrischen Kurven und Flachen aufbauend auf dem YART Graphikkern. Master's thesis, Technische Unversitat Ilmenau, 1995. [11] K. Mrken. Some identities for products and degree raising of splines. Constructive Approximation, 7:195{ 208, 1991. [12] M. Plavnik and G. Elber. Surface design using global constraints on total curvature. To retrieve at ftp://ftp.cs.technion.ac.il/pub/misc/gershon/papers, February 1998. [13] L. Ramshaw. Blossoming: A connect{the{dots approach to splines. Technical Report 19, Digital Systems Reaearch Center, Palo Alto CA, June 1987. [14] A. Rappoport, Y. Hel-Or, and M. Werman. Interactive design of smooth objects with probabilistic point constraints. ACM Transactions on Graphics, 13(2):156{176, April 1994. [15] J. J. Risler. Mathematical Methods for CAD. Press Syndicate of the University of Cambridge, 1991.