Algorithm to determine the intersection curves between bezier surfaces by the solution of multivariable polynomial system and the differential marching method

Algorithm to Determine the Intersection Curves between Bezier Surfaces by the Solution of Multivariable Polynomial System and the Differential Marching Method

Mário Carneiro Faustini

Marcos Sales Guerra Tsuzuki

Universidade de São Paulo, Escola Politécnica, Departamento de Engenharia Mecânica, Av. Prof. Mello Moraes, 2231 - 05508-900, São Paulo - SP, Brasil

The determination of the intersection curve between Bézier Surfaces may be seen as the composition of two separated problems: determining initial points and tracing the intersection curve from these points. The Bézier Surface is represented by a parametric function (polynomial with two variables) that maps a point in the tridimensional space from the bidimensional parametric space. In this article, it is proposed an algorithm to determine the initial points of the intersection curve of Bézier Surfaces, based on the solution of polynomial systems with the Projected Polyhedral Method, followed by a method for tracing the intersection curves (Marching Method with differential equations). In order to allow the use of the Projected Polyhedral Method, the equations of the system must be represented in terms of the Bernstein basis, and towards this goal it is proposed a robust and reliable algorithm to exactly transform a multivariable polynomial in terms of power basis to a polynomial written in terms of Bernstein basis .

Keywords: geometric modeling, parametric surfaces, intersection curves, multivariable polynomial systems

Introduction

In the development of solid modelers supporting the geometric representation of surfaces, it is crucial that the problem of determining the intersection curve between surfaces is well algorithmically resolved. The applications of such an algorithm would be the generation of the solid resulting from boolean operations, the definition of the toolpath for sculptured surface machines, the program definitions to stereolithography machines and smoothing edges of solid models.

This problem was initially addressed by Hosaka (1992), where it is shown that the intersection problem may be split in two different and sequential problems: determining initial points of the intersection curve and, then, tracing it from these points. This way, it is possible to define an algorithm to each application domain. However, some quite common situations couldn’t be included in the application domain of existing algorithms. For instance, with the use of these previous methods sometimes it isn’t possible to determine the intersection curve’s initial points if there are singularities^* on the curve.

Several authors studied the determination of initial points for tracing intersection curves (Kriezis et al., 1992; Hu et al., 1997; Grandine & Klein, 1997). However, as presented in (Faustini, 1999), comparing the approaches of (Hu et al., 1997) and (Grandine & Klein, 1997), we have that the first one demands a larger processing time, as long as it determines a higher number of initial points. The second approach, on the other hand, finds a more concise number of initial points, thus demanding a smaller processing time. In both approaches, the determination of initial points faces the task of solving multivariable polynomial systems. To deal with this task, it is used the Projected Polyhedral Method (Sherbrooke & Patrikalakis, 1993). Nevertheless, in order to use this method it is necessary to convert multivariable polynomials from power basis to Bernstein basis, what represents the new contribution from this work (Faustini, 1999). This way, once the initial points are determined, it remains the problem of tracing the intersection curves. For that, it is used the proposal presented by (Hu et al., 1997) that determines the curve by differentiating the vectorial equations of the surfaces at a given intersection point. There is a singular case that occurs when both surfaces have a common tangent plane, and this problem is solved using the oriented distance function presented in Wang et al. (1991) and Kriezis et al. (1992) which allows to process these singular situations.

All the algorithms presented in this article were implemented and tested, and some results are shown at the last section.

Bézier Surfaces

The Bézier surfaces are widely used in CAD applications. Such parametric surfaces are defined by the following expression (Hosaka, 1992):

(1)

where P_ij are the control points of the surface and and are Bernstein polynomials, defined by:

(2)

where is the binomial coefficient:

(3)

The Bernstein polynomials have, among others, the following main properties:

(4)

(5)

The control points are disposed in a rectangular net of n+1 by m+1 points. This way, the position of a surface point in the xyz space is given by:

(6)

where each coordinate is defined by a polynomial in u and v. The surface degree (and therefore, the degree of the polynomials that describe each coordinate) is given by m and n. For a matter of simplicity, the surfaces presented in this work are all 3^rd degree Bézier surfaces (m=n=3). Nevertheless, as it will become clear, the presented method isn’t limited by the surface’s degree.

Figure 1 shows a 3^rd degree Bézier surface and its respective control points. It is interesting to notice that the surface points are obtained from the u and v parameters, so that it can be mapped by the parametric domain. In other words, for each point of the parametric uv domain there is a correspondent point on the surface in the xyz space.

Problem Definition

In this section, the problem of determining the intersection curve between two surfaces is defined in the following way: given F(u, v) and G(s, t) surfaces parametrized within the domain [0, 1]¹, it is desired to obtain all the pairs (u, v) and (s, t) so that

(7)

The vectorial expression above encloses a system of 3 equations and 4 unknowns, which makes a direct solution quite complex. In order to trace the intersection curves, some marching method (Hu et al., 1997) may be used to generate the curves from some discrete points on them (initial points). Obtaining these points is the main task of this work.

The solution of Eq. (7) may consist of none, one or more than one isolated curves, as shown by Figure 2. To use the marching method, it is necessary to have at least one initial point for each isolated intersection curve.

Determining the Initial Points

The following approach, presented by Grandine & Klein (1997), intends to find the initial points for the intersection curves. Given the surfaces F(u, v) and G(s, t), it is needed to obtain the set u, v, s and t so that:

(8)

or, if we reparametrize by the variable t representing the arc length of the intersection curve, it becomes:

(9)

In order to determine the points of the intersection curves in which they cross the borders of the parametric domains of F and G, we may substitute each parameter for zero or one in Eq. (8) and solve the remaining 3x3 system. For u = 0 or 1 and v = 0 or 1, we obtain initial points for the intersection curves that cross the borders of the parametric domain of F and that may correspond to any point inside the parametric domain of G. In the other hand, for S=0 or 1 and t = 0 or 1, we obtain the initial points for curves that may begin and end anywhere inside the domain of F, but that for sure are limited by the borders of the parametric domain of G. This way, the cases of the border points of F and G may be processed separately.

Supposing that all the zeros on the borders of the parametric domain of F have been determined (u = 0 and 1, v = 0 and 1 in Eq. (8)). The resulting points could be those seen in Figure 3.

So, the possible topologies of the intersection curves limited by the points in Figure 4 are the ones shown in Figure 4. Other than those, we may have intersection curves that cross no border of any parametric domain. These curves are called intersection loops, and they may be in any number and at any position inside the parametric domain.

In order to determine these loops, we "scan" the parametric domain of the surface F with parallel lines inclined of q in relation to the direction of parameter v, as shown in Figure 5.

This way, the vector (- sin q , cos q ) is parallel to the scanning lines, with 0 £ q £ p /2. This allows us to determine the points where the intersection curves change their direction from these lines and return towards them. These are called turning points, and they are characterized by the fact that the tangent vector of the curve on this point is parallel to the direction of the scanning lines. In the example of the intersection curves of Figure 6, it is shown how the scanning lines "find" the turning points both for opened curves and loops, detecting the presence of the last ones and providing their initial points.

A turning point is, then, given by:

(10)

what means that (u’(t ) , v’(t )) is perpendicular to (cos q , sin q ).

Differentiating Eq. (9):

(11)

So, we have that the vector (u’(t ) , v’(t ) , s’(t ) , t’(t )) is inside the null space of the 3x4 matrix [F_uF_v -G_s -G_t]. This null space is given by multiples of the vector:

(12)

From Eqs. (10), (11) and (12), it can be seen that the turning points are the solutions of the following 4x4 system:

(13)

The solution set for the system above will also contain the points where one intersection curve crosses another (called critical points) in addition to the turning points.

It is interesting to notice that there are some special cases where the algorithm doesn’t work efficiently. For instance, the situation in which an intersection curve consists of a straight line within the parametric domain of F inclined by the same angle q (from the direction of the parameter v) chosen to the method, then all the points of this curve will represent solutions of the system of Eq. (13). The same limitation happens when there is an intersection curve (open or loop) of singular points, in this case all of them will be solutions of Eq. (13).

Finally, the main idea of this method is to find initial points for each of the intersection curves. The solution of the original 3 equations and 4 unknowns system in Eq. (8) is simplified with the addition of a 4^th equation, originating the 4x4 system of Eq. (13). It is interesting to notice that this system fundamentally deals with intersection loops, since initial points for opened intersection curves are given by the determination of the points where the intersection curves cross the borders of the parametric domains of the surfaces, as explained in the beginning of this section.

Therefore, it must be chosen a robust and reliable numerical method to solve the multivariable polynomial systems that arise, since their solutions will represent the desired initial points which will be classified in lists that will define the relations among them, allowing a complete topological analysis of the intersection curves as seen through the parametric domains. In order to fulfill this task, it is used the Projected Polyhedral Method, explained in the next section.

Projected Polyhedral Method

Solving multivariable polynomial systems means, in algebraic terms, determining all the n-uples x = (x ₁ , x ₂ , ... , x _n) so that

(14)

with x_i Î [0,1]^** and f _k a polynomial with n variables.

The projected polyhedral method was fully presented in (Sherbrooke & Patrikalakis, 1993). However, as it is a complex algorithm, we will present a brief description of the method. Its first step is to transform each equation f_k to the Bernstein basis, what means:

(15)

where d_i^(k) is the degree of the variable x_i in f_k .

Now, let’s redefine the problem of Eq. (14) as the task of determining the intersection of the graphs^*** of each f_k (each graph can be seen as an hypersurface in Âⁿ⁺¹) and the hyperplane x_n+1 = 0. Each graph of f_k is given by

(16)

So that the system of Eq. (14) becomes:

(17)

Since each graph F_k is also a parametric hypersurface, then its control points would be:

(18)

Therefore a new way to look at the problem of solving the system of Eq. (14) would be to work out the equivalent system of Eq. (17). And this means to determine the intersection of the hypersurfaces defined by each F_k. An initial approach to this problem would be determining the intersection among the convex-hulls that are formed by the control points of theses hypersurfaces in the space  ⁿ⁺¹. Moreover, it is possible to convert this problem of dimension (n + 1) into n bidimensional problems by projecting the control points of each F_k on every plane of the space  ⁿ⁺¹ formed by one coordinate related to the variable x_j and the last coordinate of Eq. (16), related to the original function f_k (in other words, one plane for each variable of the initial system of Eq. (14)). This way, the control points of each and every F_k are projected on n planes, so that for the plane j, of the variable x_j, their coordinates would be

(19)

Then it is formed, in each plane, the bidimensional convex-hull of the projected control points of each and every F_k. And within each plane, it is determined the intersection of all convex-hulls with the segment [0,1] (that is related to the domain of the variable x_j) of the abscissa of the plane. The result may be one segment, one point or null (because all the projected convex-hulls are convex polygons). If null, then there is no value for the variable of this plane to be considered as a solution of the system. If the result is a point, then its value (a real within [0,1]) is a solution itself for the related variable in the system of Eq. (14). And if the result is a segment, then it is still possible that inside of it there are one or more points of solution to the related variable in the system.

This means that, for each variable x_i, we will possibly have decreased the range of search from the initial variable’s domain [0,1] down to a segment [a_i, b_i]. So, we can reduce the box of solutions from the initial full domain [0,1]ⁿ down to the newly found box of solution:

(20)

At this point, for each k, we define a new function f_k’ so that:

(21)

This function will be mapped within the domain [0,1]ⁿ in the same way that f_k was mapped by the domain of Eq. (14).

The whole process is then iteratively repeated until we find boxes of solution (like that of Eq. (20)) with all sides small enough to be seen as isolated roots for the variables of the system. In other words: in S (Eq. (20)), " i so that 1£ i £ n, it must be true that (b_i - a_i) £ tolerance. This tolerance defines the precision of the solutions.

However, after processing the intersection of the convex-hulls within the planes, it is possible that the value (b_i - a_i) of some variables in S (Eq. (20)) may have not significantly decreased from 1. This happens when there is more than one solution to the system of Eq. (14). In such case, we must split in two the segment of the variable that wasn’t considerably reduced. So on, we deal with the two resulting boxes as independent and separated problems. For instance, if we obtained the box S = [a₁, b₁] x [a₂, b₂] x...x [a_n, b_n] and (b₂ - a₂) > criteria (» 1), then two solution sub-boxes will be separately considered at next iteration:

and

It is interesting to notice that such subdivision is done recursively until a sub-box with only one solution is isolated.

Basis Conversion of Multivariable Polynomials From Power To Bernstein Basis

Such basis conversion is a fundamental step of the Projected Polyhedral Method, as seen in the previous section, and its evaluation represents the main problem found when implementing such algorithm. A generic approach may be proposed in pursue of this task (Faustini, 1999). This way, the following polynomial, function of the parameters u₁, u₂ until u_n with the maximum degree of each of them as being, respectively, d₁, d₂ until d_n:

(22)

may be rewritten in terms of the Bernstein basis:

(23)

Then it is possible to design an algorithm for each case, when needed, from the following generic element of the matrix of weights W (as in Eq. (15)):

(24)

Tracing the Intersection Curves: The Marching Method with Differential Equations

Once the initial points of the intersection curves are obtained, each curve C(w), parametrized by the arc length w, is traced using the Marching Method proposed by (Hu et al., 1997).

The marching direction is given by the tangent vector of C(w), which is perpendicular to the normal vector of both surfaces at the given intersection point:

(25)

Differentiating the surfaces in w, we obtain:

(26)

Manipulating:

(27)

Above are the equations from which the points of the intersection curves may be successively determined using a numerical integration method for the problem of the initial value with the given system of ordinary differential equations of Eq. (27) Hu et al. did not explain which method was used to integrate the ordinary differential equations, as a first implementation we decided to use the 4^th order Runge-Kutta Method to accomplish this task.

So, the Marching algorithm has the following steps:

Obtaining initial points for the algorithm, which means at least one point for each intersection curve. This is done using the procedures presented in the previous sections.

For each of these points, are obtained the derived equations of the surface parameters with respect to the arc length of the intersection curve using Eq. (27).

This system of ordinary differential equations is integrated using a numerical method (4^th order Runge-Kutta, for instance), providing this way the next point of the intersection curve.

Steps 2 and 3 are repeated until the traced curve reaches a turning point or an initial border points (previously determined).

Special Cases: Situations of indetermination and the oriented distance function

However, some special cases may be found when tracing the curves, and they result in the indetermination of Eq. (25) because both surfaces have parallel normal vectors at the singular curve point. To deal with this problem, we use the oriented distance function as an alternative method to obtain the derived equations usually determined by Eq. (27).

The real oriented distance function f between the surface G and a point moving on the surface F, presented in (Wang et al., 1991; Kriezis et al., 1992), is defined by:

(28)

where g is the  ^3® [0,1]² function that provides the values of the parameters s and t of G, that define the point on this surface that is closer to a given point of  ³, G(g(F(u,v))) is a point on the surface G which is closer to the point defined by F(u,v), and N^G[g(F(u,v))] is the normal unitary vector of the surface G at that same point. Thus, |f (u,v)| denotes the actual euclidean distance between the point F(u,v) and G(g(F(u,v))), since the vectors N^G[g(F(u,v))] and [F(u,v) - G(g(F(u,v)))] are colinear (colinearity condition). In this case, G(g(F(u,v))) is the orthogonal projection of the point F(u,v) over the surface G, so that the first and second order derived equations of f may be explicitly determined. Thus, the intersection between the surfaces is equivalent to the set of points that satisfy the implicit equation f (u,v) = 0.

Using the first and second order derived equations of f , we may develop the expressions to obtain the values of the derived equations of Eq. (27) at the singular points.

So, when, during the Marching Method presented by Eq. (27), the algorithm finds a point that defines an indetermination in Eq. (25), then it uses an alternative approach by solving the following equation:

(29)

where H is the Hessian of f , which means:

(30)

The second order derived equations of f are given by:

(31)

(32)

(33)

The partial derived equations of the normal vectors with crossed parameters presented above are given by:

(34)

(35)

where L is the orthogonal projection tensor of surface G(s,t), given by:

(36)

The direct partial derived equations of the normal vectors are:

(37)

(38)

If we analyze the sign of the determinant of H, we face three possibilities:

if det[H] > 0 then we have an isolated point of tangency between the surfaces;

if det[H] = 0 then we have a cusp point (Figure 7 (a));

se det[H] < 0 the we have a saddle point (Figure 7 (b)).

This way, a general solution of Eq. (29) (when det[H] £ 0) is:

(39)

and to obtain l ₁we impose that .

At this point, the derived equations of s and t are determined by:

(40)

where K is the first fundamental tensor of G(s,t), given by:

(41)

and E is defined by:

(42)

There is also the case of surfaces that are tangent along a whole intersection curve. In this situation the expressions in Eq. (27) may never be used, they are undefined along the whole curve.

In other words, det[H] = 0 along a whole curve. With this, we obtain:

(43)

or

(44)

After determining l ₁ (imposing that ) and the derived equations of s and t (with the same procedure presented in Eq. (40)), we have all the derived equations needed for the Marching Method.

Results

Some examples are now presented. The algorithms to generate them were implemented in C++ programming language and executed with a Pentium 200MHz computer. All the surfaces presented in both Figure 8 and Figure 11 are 3^rd degree Bézier surfaces. In Figure 9 and Figure 12 we have initial points determined for the intersection curves, and in Figure 10 and Figure 13 the intersection curves themselves found using the algorithms presented in this article.

It is important to comment the example of surfaces from Figure 12. In Figure 13, are shown the initial points determined for the intersection curves. These should be (within a parametric domain): two border points, one critical point by the half of the interval domain of the ordinate and two points where the scanning lines (inclined of 150° referring to the direction of the ordinate) are tangent to the intersection curves. However, only one of these last two points was found with a given precision in the Projected Polyhedral Method. In order to determine also the other one, it is necessary to increase the precision in the referred method (what leads also to a significant increase in processing time). This problem, of numerical nature, is still being addressed. Nevertheless, since the missing point was redundant for that configuration of curves (as was also redundant the turning point that was found), it was possible to correctly trace the curves (as seen in Figure 13).

Conclusions

In this work we presented a complete proposal of determining the intersection curves between two Bézier Surfaces, which integrated several approaches from the literature. The main contribution of this work is a method to convert multivariable polynomials from power basis into equivalent ones in Bernstein basis, which is a fundamental step in the implementation of the Projected Polyhedral Method for solving multivariable nonlinear polynomial systems.

However, it still remains the limitation that happens when all the points along an intersection curve are singular points. Then, the algorithm determines a very large (and unnecessary) number of initial points (this number varies accordingly to the precision desired), in both the method proposed by Grandine & Klein (1997) and the one shown in (Hu et al., 1997). Thus, this is a situation that still compromises the efficiency of the algorithm.

Acknowledgements

We would like to thanks FAPESP for supporting the first author (proc. 97/02939-8) and CNPq for supporting the second author (proc. 300.224/96-6).

Manuscript received: May 1999, Technical Editor: Paulo Eigi Miyagi.

** or any other finite domain [a_i, b_i] which must be, then, mapped in the domain [0,1] through the variable transformation x_i = a_i + x_i . (b_i - a_i) in every equation f _k .

*** Be f:Âⁿ®Â a function of x. Then, its graph is the function F: Âⁿ®Âⁿ⁺¹ defined by: F(x) = (x, f(x)).

• Faustini, M. C., 1999, Proposta de Algoritmo para a Determinação da Intersecção entre Superfícies NURBS, Tese de Mestrado, Dep. de Eng. Mecânica, Escola Politécnica, USP.
• Grandine, T. A. & Klein IV, F. W., 1997, A new approach to the surface intersection problem, Computer Aided Geometric Design, No. 14, pp. 111-134.
• Hosaka, M., 1992, Modeling of Curves and Surfaces in CAD/CAM, Springer-Verlag, Berlim.
• Hu, C. Y., Maekawa, T., Patrikalakis, N. M. and Ye, X., 1997, Robust interval algorithm for surface intersections, Computer-Aided Design, Vol. 29, No. 9, pp.617-627.
• Kriezis, G. A., Patrikalakis, N. M. and Wolter, F. E., 1992, Topological and differential-equation methods for surface intersections, Computer-Aided Design, Vol. 24, No. 1, pp.41-55.
• Sherbrooke, E. C. & Patrikalakis N. M., 1993, Computation of the solutions of nonlinear polynomial systems, Computer Aided Geometric Design No. 10, pp. 379-405.
• Wang, Y., Gursoz, E. L., Chen, J. M., Prinz, F. B. and Patrikalakis, N. M., 1991, Intersection of Parametric Surfaces For Next Generation Geometric Modelers in Product Modeling for Computer-Aided Design and Manufacturing, Elsevier Science Publishers B. V., New York, pp. 75-96.

* Singularities are defined as intersection points for which the two respective normal vectors, referring to each surface, are parallel.

Singularities are defined as intersection points for which the two respective normal vectors, referring to each surface, are parallel.

Publication Dates

• Publication in this collection
  15 Dec 2000
• Date of issue
  2000

History

• Received
  May 1999