 # On the Efficiency of Strategies for Subdividing Polynomial Triangular Surface Patches

In this paper, we investigate the efficiency of various strategies for subdividing polynomial triangular surface patches. We give a simple algorithm performing a regular subdivision in four calls to the standard de Casteljau algorithm (in its subdivision version). A naive version uses twelve calls. We also show that any method for obtaining a regular subdivision using the standard de Casteljau algorithm requires at least 4 calls. Thus, our method is optimal. We give another subdivision algorithm using only three calls to the de Casteljau algorithm. Instead of being regular, the subdivision pattern is diamond-like. Finally, we present a "spider-like" subdivision scheme producing six subtriangles in four calls to the de Casteljau algorithm.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

In this paper, we investigate the efficiency of various strategies for subdividing polynomial triangular surface patches. Subdivision methods based on a version of the de Casteljau algorithm splitting a control net into control subnets (see Farin ) were investigated by Goldman , Boehm and Farin , Böhm , and Seidel  (see also Boehm, Farin, and Kahman , and Filip ). However, except for Böhm , these papers are not particularly concerned with minimizing the number of calls to the standard de Casteljau algorithm. Furthermore, some of these papers (notably Goldman ) use a version of the de Casteljau algorithm computing a -dimensional simplex of polar values, which is more expensive than the standard -dimensional version. In this paper, we give a simple algorithm performing a regular subdivision in four calls to the standard de Casteljau algorithm (in its subdivision version). A naive version uses twelve calls. We also show that any method for obtaining a regular subdivision using the standard de Casteljau algorithm requires at least calls. Thus, our method is optimal. We give another subdivision algorithm using only three calls to the de Casteljau algorithm. Instead of being regular, the subdivision pattern is diamond-like. Finally, we present a “spider-like” subdivision scheme producing six subtriangles in four calls to the de Casteljau algorithm. Some familiarity with affine spaces and affine maps is assumed. Details can be found in Farin , Berger , or Gallier .

## 2 The Polar Form Approach to Polynomial Triangular Surface Patches

The deep reason why polynomial triangular surface patches can be effectively handled in terms of control points is that multivariate polynomials arise from multiaffine symmetric maps (see Ramshaw , Farin [7, 6], Hoschek and Lasser , or Gallier ). Denoting the affine plane as , traditionally, a polynomial surface in is a function , defined such that

for all , where are polynomials in . Given a natural number , if each polynomial has total degree , we say that is a polynomial surface of total degree . The trace of the surface is the set .

Now, given a polynomial surface of total degree in some affine space (typically ), there is unique symmetric and multiaffine map such that

 F(u,v)=f((u,v),…,(u,v)m),

for all . The symmetric and multiaffine map associated with is called the polar form of .

The above result is not hard to prove. Using linearity, it is enough to deal with a single monomial. Given a monomial , with , it is easily shown that the symmetric multiaffine form corresponding to is given

Recall that a map is affine if

 f((1−λ)a+λb)=(1−λ)f(a)+λf(b),

for all , and all . A map is multiaffine if it is affine in each of its arguments. A map is symmetric if it does not depend on the order of its arguments, i.e., , for all , and all permutations .

As an example, consider the following surface known as Enneper’s surface:

We get the polar forms

Furthermore, it turns out that any symmetric multiaffine map is uniquely determined by a family of points (where is any affine space, say ). Let

 Δm={(i,j,k)∈N3 | i+j+k=m}.

The following lemma is easily shown (see Ramshaw  or Gallier ).

###### Lemma 2.1

Given an affine frame in the plane , given a family of points in , there is a unique surface of total degree , defined by a symmetric -affine polar form , such that

 f(r,…,ri,s,…,sj,t,…,tk)=bi,j,k

for all . Furthermore, is given by the expression

where , with , and .

For example, with respect to the standard frame , we obtain the following control points for the Enneper surface:

 f(r,r,r)(23,0,1)f(r,r,t)(23,0,13)f(r,r,s)(23,23,13)f(r,t,t)(13,0,0)f(r,s,t)(13,13,0)f(r,s,s)(23,23,−13)f(t,t,t)(0,0,0)f(s,t,t)(0,13,0)f(s,s,t)(0,23,−13)f(s,s,s)(0,23,−1)

A family of points in is called a (triangular) control net, or Bézier net. Note that the points in

 Δm={(i,j,k)∈N3 | i+j+k=m},

can be thought of as a triangular grid of points in . For example, when , we have the following grid of points:

 500401410302311320203212221230104113122131140005014023032041050

We intentionally let be the row index, starting from the left lower corner, and be the column index, also starting from the left lower corner. The control net can be viewed as an image of the triangular grid in the affine space . It follows from Lemma 2.1 that there is a bijection between polynomial surfaces of degree and control nets . It should also be noted that there are efficient methods for computing control nets from parametric definitions, but this will be published elsewhere.

In the next section, we review a beautiful algorithm to compute a point

on a surface patch using affine interpolation steps, the de Casteljau algorithm.

## 3 The de Casteljau Algorithm for Triangular Patches

In this section, we explain in detail how the de Casteljau algorithm can be used to subdivide a triangular patch into three subpatches. For more details, see Farin [7, 6], Hoschek and Lasser , Risler , or Gallier . In the next section, we will use versions of this algorithm to obtain a triangulation of a surface patch using recursive subdivision.

Given an affine frame , given a triangular control net , recall that in terms of the polar form of the polynomial surface defined by , for every , we have

 bi,j,k=f(r,…,ri,s,…,sj,t,…,tk).

Given in , where , in order to compute , the computation builds a sort of tetrahedron consisting of layers. The base layer consists of the original control points in , which are also denoted as . The other layers are computed in stages, where at stage , , the points are computed such that

 bli,j,k=λbl−1i+1,j,k+μbl−1i,j+1,k+νbl−1i,j,k+1.

During the last stage, the single point is computed. An easy induction shows that

 bli,j,k=f(a,…,al,r,…,ri,s,…,sj,t,…,tk),

where , and thus, .

Assuming that is not on one of the edges of , the crux of the subdivision method is that the three other faces of the tetrahedron of polar values besides the face corresponding to the original control net, yield three control nets

 Nast=(bl0,j,k)(l,j,k)∈∈Δm,

corresponding to the base triangle ,

 Nrat=(bli,0,k)(i,l,k)∈∈Δm,

corresponding to the base triangle , and

 Nrsa=(bli,j,0)(i,j,l)∈∈Δm,

corresponding to the base triangle . If belongs to one of the edges, say , then the triangle is flat, i.e. is not an afine frame, and the net does not define the surface, but instead a curve. However, in such cases, the degenerate net is not needed anyway.

From an implementation point of view, we found it convenient to assume that a triangular net is represented as the list consisting of the concatenation of the rows

 bi,0,m−i,bi,1,m−i−1,…,bi,m−i,0,

i.e.,

 f(r,…,ri,t,…,tm−i),f(r,…,ri,s,t,…,tm−i−1),…,f(r,…,ri,s,…,sm−i−1,t),f(r,…,ri,s,…,sm−i),

where . As a triangle, the net is listed (from top-down) as

The main advantage of this representation is that we can view the net as a two-dimensional array net, such that (with ). In fact, only a triangular portion of this array is filled. This way of representing control nets fits well with the convention that the affine frame is represented as follows:

Instead of simply computing , the de Casteljau algorithm can be easily adapted to output the three nets , , and . The function sdecas3 does that. We also found it convenient to write three distinct functions subdecas3ra, subdecas3sa, and subdecas3ta, computing the control nets with respect to the affine frames , , and . An implementation in Mathematica can be found in Gallier .

## 4 Regular Subdivision Of Triangular Patches

If we want to render a triangular surface patch defined over the affine frame , it seems natural to subdivide into the three subtriangles , , and , where is the center of gravity of the triangle , getting new control nets , and using the functions described earlier, and repeat this process recursively. However, this process does not yield a good triangulation of the surface patch, because no progress is made on the edges , , and , and thus, such a triangulation does not converge to the surface patch. Thus, in order to compute triangulations that converge to the surface patch, we need to subdivide the triangle in such a way that the edges of the affine frame are subdivided. There are many ways of performing such subdivisions, and we will propose a method which has the advantage of yielding a very regular triangulation, and of being very efficient. In fact, we give an optimal method for subdividing an affine frame using four calls to the standard de Casteljau algorithm in its subdivision version. A naive method would require twelve calls.

Goldman  proposed several subdivision algorithms, including one for splitting a triangular patch into four triangular subpatches, but his methods use a generalized version of the de Casteljau algorithm computing a -simplex of polar values. These methods are illustrated graphically in Boehm and Farin . It should be noted that Boehm and Farin do mention that it is possible to compute the control net w.r.t. a new affine frame from the control net w.r.t. an original affine frame in three calls to the standard de Casteljau algorithm. However, they do not explain how to split a triangular patch into four subpatches using four calls to the standard de Casteljau algorithm.

Goldman’s subdivision methods can be justified in a very simple way as shown by Seidel . Given a surface of total degree defined by a triangular control net , w.r.t. the affine frame , for any points (where ), the following -simplex of points where is defined inductively as follows:

where .

If is the polar form of , it is easily shown that

 bl1,…,lni,j,k=f(r,…,ri,s,…,sj,t,…,tk,p1,…,p1l1,…,pn,…,pnln).

For , , as in the standard de Casteljau algorithm. For , is a control net of w.r.t. .

In particular, if are chosen on the edges of , the subnets for the four subpatches are obtained. Goldman observes that some of the nets involved in the computation are trivial, but still, a -simplex of polar values is computed.

It was brought to our attention by Gerald Farin (and it is mentioned in Remark 2 of Seidel’s paper , page 580) that Helmut Prautzsch showed in his dissertation (in German)  that regular subdivision into four subtriangles can be achieved in four calls to the standard de Casteljau algorithm. Prautzsch’s method is briefly described in Böhm , page 348 (figure) and page 349 (in fact, with a typo, one of the barycentric coordinates listed as should be ). The order in which the four patches are obtained is slightly different from ours. Since Prautzsch’s algorithm has not been discussed more extensively in the literature, we feel justified in presenting our method.

The subdivision strategy that we will follow is to divide the affine frame into four subtriangles , , , and , where , , and are the middle points of the sides , and respectively, as shown in the diagram below:

First, we show that any method using the standard version of the de Casteljau algorithm for subdividing a triangular patch into subpatches forming a regular pattern as above requires calls. The crux of the argument is that a call to the de Casteljau algorithm in its subdivision version produces three subpatches containing only one new corner. We want to produce the four subpatches , , , and . After one subdivision step, we have three patches each involving exactly one of . After two subdivision steps, we have six subspatches only two of which involve exactly two of , since we can only subdivide a single patch, and since this patch only has one of . Thus, at least three steps are needed to produce four subpatches involving at least two of . If we produced the patch during the third subdivision step, we would have three patches involving exactly two of , but the subdivision step that produced also produces two patches sharing the same vertex from . However, , , and do not share a vertex from . If was not produced during the third step, at least four steps are needed. Therefore, in all cases, at least four steps are needed to produce the required four subpatches.

We now present our algorithm. The first step is to compute the control net for the affine frame . This can be done using two steps. In the first step, split the triangle into the two triangles and , where is the middle of . Using the function sdecas3 (with ), the nets , , and are obtained, and we throw away (which is degenerate anyway). Then, we split into the two triangles and . For this, we need the barycentric coordinates of with respect to the triangle , which turns out . Using the function sdecas3, the nets , , and are obtained, and we throw away .

We will now compute the net from the net . For this, we need the barycentric coordinates of with respect to the triangle , which turns out to be . Using the function subdecas3sa, the net is obtained.

We can now compute the nets and from the net . For this, we need the barycentric coordinates of with respect to the affine frame which turns out to be . Using the function sdecas3, the snet , , and are obtained, and we throw away .

Finally, we apply to the net to get the net , to to get the net , followed by to the net to get the net , and twice to to get the net ,

Thus, using four calls to the de Casteljau algorithm, we obtained the nets , , , and .

Remarks:

1. For debugging purposes, we assigned different colors to the patches corresponding to , , , and , and we found that they formed a particularly nice pattern under this ordering of the vertices of the triangles. In fact, is blue, is red, is green, and is yellow.

2. In the last step of our algorithm, the subdivision step is performed with respect to a point of barycentric coordinates . One might worry that such a step involving a nonconvex combination is a source of numerical instability. We tested our algorithm on many different examples, and so far, without running into any problem. We also believe that such a nonconvex step is unavoidable if the standard de Casteljau algorithm (building a simplex of polar values of dimension ) is used, but we are unable to prove this.

The subdivision algorithm just presented has been implemented in Mathematica, see Gallier . The subdivision method is illustrated by the following example of a cubic patch specified by the control net

 net = {{0, 0, 0}, {2, 0, 2}, {4, 0, 2}, {6, 0, 0},
{1, 2, 2}, {3, 2, 5}, {5, 2, 2},
{2, 4, 2}, {4, 4, 2},  {3, 6, 0}};


After only three subdivision steps, the triangulation approximates the surface patch very well.

As another example of the use of the above functions, we can display a portion of a well known surface known as the “monkey saddle”, defined by the equations

 x=u,y=v,z=u3−3uv2.

Note that is the real part of the complex number . It is easily shown that the monkey saddle is specified by the following triangular control net monknet over the standard affine frame , where , , and .

monknet = {{0, 0, 0}, {0, 1/3, 0}, {0, 2/3, 0}, {0, 1, 0},
{1/3, 0, 0}, {1/3, 1/3, 0},  {1/3, 2/3, -1},
{2/3, 0, 0}, {2/3, 1/3, 0},  {1, 0, 1}};


We actually display the patch over the rectangle . This can be done by splitting the square into two triangles, and computing control nets with respect to these triangles. This is easy to do, and it is explained for example in Gallier . Subdividing both nets times, we get the following picture.

## 5 A Diamond-Shape Strategy For Subdivision

The strategy of the previous section was to split the affine frame into four congruent subtriangles. We were able to do this using four calls to the de Casteljau algorithm and we showed that it is not possible to do it in fewer calls.

However, it is possible to split the affine frame into four subtriangles using only three calls to the de Casteljau algorithm. The method consists in splitting the triangle into the four subtriangles , , , and :

This can be done by first computing the nets and , which can be done in one call to sdecas3 (dropping ). Next, we split into the two triangles and . For this, we need the barycentric coordinates of with respect to the triangle , which turns out . Using the function sdecas3, the nets , , and are obtained, and we throw away . Finally, we split into the two triangles and . For this, we need the barycentric coordinates of with respect to the triangle , which turns out . Using the function sdecas3, the nets , , and are obtained, and we throw away .

An implementation of the method is given in Gallier .

The result of subdividing two of three times reveals some diamond-shape subdividion patterns. For example, after three iterations, the dome surface is subdivided as follows:

## 6 A Spider-Web Strategy For Subdivision

Is is also possible to split the affine frame into six subtriangles using only four calls to the de Casteljau algorithm.

The triangle is subdivided into the triangles , , , , , and , as shown above, where is the center of gravity. This can be done by first computing the nets , , and , which can be done in one call to sdecas3. We split into the two triangles and using sdecas3 (throwing away ). We split into the two triangles and using sdecas3 (throwing away ). Finally, we split into the two triangles and using sdecas3 (throwing away ).

The result of subdividing recursively yields spider-web like patterns. For example, after three iterations, the dome surface is subdivided as follows:

## 7 Conclusion

We have presented various strategies for subdividing polynomial triangular surface patches. We gave an algorithm performing a regular subdivision in four calls to the standard de Casteljau algorithm, and we showed that this method for obtaining a regular subdivision is optimal. We gave another subdivision algorithm using only three calls to the de Casteljau algorithm. Instead of being regular, the subdivision pattern is diamond-like. Finally, we presented a “spider-web” pattern subdivision scheme producing six subtriangles in four calls to the de Casteljau algorithm. These methods immediately apply to rational surface patches (Gallier ).

An amusing effect is obtained from the regular subdivision scheme if we omit the central triangle . We obtain a “fractalized” representation of the surface patch, in the sense that a Sierpinski gasket pattern is laid onto the patch! It would be interesting to investigate other subdivision strategies and the patterns that they induce, especially if the triangles are colored in various recursive manners.

## References

•  Marcel Berger. Géométrie 1. Nathan, 1990. English edition: Geometry 1, Universitext, Springer Verlag.
•  W. Boehm and G. Farin. Letter to the editor. Computer Aided Design, 15(5):260–261, 1983. Concerning subdivison of Bézier triangles.
•  W. Boehm, G. Farin, and J. Kahman. A survey of curve and surface methods in CAGD. Computer Aided Geometric Design, 1:1–60, 1984.
•  Wolfgang Böhm. Subdividing multivariate splines. Computer Aided Design, 15:345–352, 1983.
•  G. Farin. Triangular Bernstein-Bézier patches. Computer Aided Geometric Design, 3:83–127, 1986.
•  Gerald Farin. NURB Curves and Surfaces, from Projective Geometry to practical use. AK Peters, first edition, 1995.
•  Gerald Farin. Curves and Surfaces for CAGD. Academic Press, fourth edition, 1998.
•  D.J. Filip. Adapative subdivision algorithms for a set of Bézier triangles. Computer Aided Design, 18:74–78, 1986.
•  Jean H. Gallier. Curves and Surfaces In Geometric Modeling: Theory And Algorithms. Morgan Kaufmann, first edition, 1999.
•  Jean H. Gallier. Simple methods for drawing rational surfaces as four or six Bézier patches. Technical report, University of Pennsylvania, Levine Hall, 3330 Walnut, Philadelphia, PA 19104, 2006.
•  R. N. Goldman. Subdivision algorithms for Bézier triangles. Computer Aided Design, 15:159–166, 1983.
•  J. Hoschek and D. Lasser. Computer Aided Geometric Design. AK Peters, first edition, 1993.
•  Helmut Prautzsch. Unterteilungsalgorithmen für multivariate splines. PhD thesis, TU Braunschweig, Braunschweig, Germany, 1984. Dissertation.
•  Lyle Ramshaw. Blossoming: A connect-the-dots approach to splines. Technical report, Digital SRC, Palo Alto, CA 94301, 1987. Report No. 19.
•  J.-J. Risler. Mathematical Methods for CAD. Masson, first edition, 1992.
•  H.P. Seidel. A general subdivision theorem for Bézier triangles. In Tom Lyche and Larry Schumaker, editors, Mathematical Methods in Computer Aided Geometric Design, pages 573–581. Academic Press, Boston, Mass, 1989.