- Email: [email protected]
A modified Clough-Tocher interpolant
Computer Aided Geometric North-Holland
A modified
Presented
19
Design 2 (1985) 19-27
at Oberwolfach
Clough-Tocher
15 November
interpolant
1984
Abstract. Clough-Tocher
schemes are used to create surfaces that interpolate to positional and gradient values that may be arbitrarily located. A standard and a modified version of the CloughkTocher mterpolation scheme are described in term5 of Bernstein~BCzier polynomials over triangles. Both method, are compared by a method that simulates reflection lines. Keywords. Surfaces. triangular
interpolants,
Bernstein-Bkzier
polynomials.
reflectlon
line\
1. Introduction Clough-Tocher interpolants were originally invented as a tool for the finite element method [Strang, Fix ‘731, but they have also turned out to be a useful tool in CAGD, namely in the field. of scattered data interpolation [Lawson ‘771. A Clough-Tocher interpolant produces a C’ piecewise polynomial surface, defined over a triangulation of given data sites. Recently the original interpolant has been generalized [Alfeld ‘85. Sablonniere ‘851, but in this paper we present a modification rather than a generalization. It is based on increasing the smoothness of the overall interpolant by an appropriate specification of the available degrees of freedom. The relevant smoothness conditions are conveniently formulated in terms of Bernstein-Bezier polynomials defined over triangles. The modified method is then compared to the standard method, and to this end a (surface interrogation) method is introduced that simulates ‘reflection lines’.
2. Continuity
conditions
for triangular BCzier patches
For the theory of BernsteinBCzier polynomials defined over triangles. the reader is referred to [Farin ‘80,‘83,‘85,BBhm/Farin/Kahmann ‘841. An nth degree Bernstein-Bezier polynomial defined over a triangle is of the form
,+J+h=" where U, u, M:are barycentric coordinates of the domain triangle and the b,,,,l, are called Bezier ordinates of p, see Fig. 1. The Beziei ordinates form a triangular net, the socalled “control net”. The Bezier ordinates of a boundary curve are the corresponding boundary Bezier points of the net. Similarly, a cross-boundary derivative, evaluated along that boundary, depends only on the first two rows of BCzier ordinates parallel to that boundary. Consider now two adjacent triangles and a Bezier net defined over each of them. We ask for the conditions that the nets must satisfy in order for the two patches defined by them to be C’ or C2. 0167-8396/M/$3.30
8) 1985, Elsewer Science Publishers
B.V. (North-Holland)
Fig. 1. Bernstein-Bkzier one triangle.
Fig. 2a. Adjacent
triangular
Fig. 2b. Adjacent
patches
triangular
wjith C’ continuity
patches
with C2 continuity
polynomial
over
G. Fmn
VLL
/ A modified Clough - Tocher inrerpoiani
Fig. 3. Three triangle.
21
C’ control
nets over split
For the C’ case, Fig. 2a illustrates the answer: The condition for C’ continuity is that the shaded pairs of triangles he coplunur (Farin ‘80). For C2 continuity, the shaded pairs of triangles in Fig. 2b are constructed to be coplanar. Each of the marked points will in general correspond to two different z-values. two from each control net. The condition for C2 continuity is thut these values coincide. These ‘auxiliary points’ are analogous to the auxiliary points used to define C2 conditions for Bezier curves [Bohm/ Farin/Kahmann ‘841. Those points could also be interpreted as B-spline coefficients. It is an interesting (and hence open) question whether the auxiliary points from Fig. 2 play a similar role for triangular splines.
3. A parameter count As an application of the above continuity conditions, we consider three Bernstein-Bezier triangles as in Fig. 3. We require any two adjacent patches to be C’; i.e. we consider C’ piecewise cubits defined over a triangulation that is obtained by splitting a triangle (called ‘macro-triangle’) into three triangles (called ‘mini-triangles’) at its centroid. These C’ cubits form a linear space, and we are interested in its dimension (i.e. the degrees of freedom present in the three C’ BCzier nets from Fig. 3). To count these degrees of freedom, one observes that the centroid and the three Bezier points next to it must be coplanar, thus providing three degrees of freedom. Moving out from the centroid one more step, we find six more Bezier points. Taking into account that they are related by the C’ conditions. one finds that they add three more degrees of freedom. Moving one step outwards again, we now have to consider the nine boundary Bezier points. They provide six additional degrees of freedom; the linear space of C’ cubits ooer the split triangle therefore has dimension twelve. The above counting procedure also leads to an interesting result: C’ cubits ouer the split triangle are C2 at the centroid [Farin ‘80, Alfeld ‘851. For a proof, strip the piecewise cubits from Fig. 3 of the boundary Bezier points, leaving a C’ piecewise quadratic function. The above counting of degrees of freedom asserted that such a function has six degrees of freedom. This is the dimension of the linear space of (global) quadratics defined over the macro-triangle.
G. Farm / A modified Clough - Tocher inierpohnr
22
They form a subspace of C’ piecewise quadratics. Since both spaces are of dimension six, they must be equal. Thus C’ piecewise quadratics over the split triangle are actually global quadratics, i.e. instead of being just C’, they are also C2. The considered piecewise quadratic and the original piecewise cubic agree in all derivatives up to order two at the centroid. Hence C’ cubits over the split triangle are actually C2 at the centroid.
4. The Clough-Tocher
interpolant
Suppose z-values and gradients are given over a set of triangulated data points, and one wants a C’ piecewise cubic surface that interpolates to this data. As an ‘initial guess’, the following piecewise cubic interpolation scheme is useful: for each triangle, define a cubic Bernstein-Bezier cubic. Its nine boundary Bezier ordinates (all the h,,,,I except 6,,,,,) are readily determined from the given data by univariate cubic Hermite interpolation. The remaining term h,,,,, is given by b 1.1.1= i(bz.o,l + b,,,., + b,.,,, + bo.l.2 +
b2.1.0
+
b,.,.,,)
-Theo
+
h,,.,
+
h.u.3).
(2)
the nine-parumeter interThis choice of b,,,,, ensures quadratic precision of this interpolant, [Farin ‘831. A surface composed from nine-parameter interpolants will however fail to be C’ in general. It must be modified to ensure the desired overall C’ smoothness. This is done by splitting each triangle in the triangulation into three mini-triangles and consequently splitting each nineparameter cubic into three ‘mini-cubits’ (for the algorithm to do this, see, e.g., [Bohm/Farin/ Kahmann ‘841). The subdivided surface now has enough degrees of freedom (twelve per subdivided macro-triangle instead of ten before the subdivision) to allow C’ continuity of the overall interpolant. In Fig. 4, let C,, P2, P, and C,, P,, P2 be the vertices of two adjacent microtriangles, coordinates of the opposite respectively. We can express C, and C, in terms of barycentric triangle:
polunt
c,=K2+6P,+ivP2,
(3)
c, = UC, + UP2+ WP, . While the C’ continuity
conditions
(4) from Section
2 are fulfilled
for the triangle
pairs formed
Pl
p2
Fig. 4. Adjacent C’ triangular defined over mini-triangles.
patches
23
G. Furin / A modified Clough - Tocher rnrerpolunt
by al, b,, subtriangles.
u2.
cl and
a3.
b,,
The C’ condition
by the middle a4. c2, they are in general violated that has to be met can be expressed as
pair
J = vu3 + wa2 + ux.
of
(5)
The standard way to guarantee (5) is as follows: let L denote any direction not parallel to P,. P2. Then the directional derivative of the mini-cubic p, defined over C,, P2. P, is a univariate quadratic Bezier polynomial with Bezier ordinates 3(l,b,
+ j2u2 + [,a,),
3(1,x + 12a, + /,a,).
3(l,bz + /,a4 + laj).
Here, I,, I,, I, denote the barycentric representation of the direction L in terms of the triangle C,, P2. P,. One way to fix the unknown x is to demand that the quadratic directional derivative be linear instead of quadratic. This condition can be expressed as [see Bohm/Farin/ Kahmann ‘841: (I,b,
+ 12a2 + /,a,)
+ 2( I,x + I,u,
+ /,a,)
+ (1,b,
+ /,a,
+ /,a,)
= 0.
(6)
It is now possible to determine x from (6). The remaining unknown y may be found from (5) or by an analogous procedure for the mini-cubic p2, defined over triangle C,, P,, P2. It is essential that L denote the same direction for both p, and p2. The standard way to achieve this is to make L perpendicular to P,P,. The obvious drawback is that the resulting interpolant is not affinely invariant. An alternative is to set L = Cl - C,, but this does not seem to have been considered in the literature. The above method is then applied to all pairs of adjacent macro-triangles. Returning to Fig. 3, we see that we thus determine all Bezier ordinates marked by full circles. The remaining ones, marked by open circles, have to be computed in order to ensure internal C’ continuity. This is accomplished by applying the Cl-conditions four times [see e.g. Farin ‘831.
5. A smoother interpolant The ‘linear cross boundary derivative’ condition (6) appears to be somewhat arbitrary, and in this section we will develop a geometrically more meaningful alternative. The C’ condition (5) is not sufficient to determine both x and JJ, and consequently an additional condition is needed. We have already (in Section 3) observed that C’ cubits over the split triangle enjoy an extra degree of continuity at the centroid, and thus define a surface that is smoother than could be expected ab initio. Therefore it seems reasonable to come as close to C2 continuity as possible across the edges of adjacent macro-triangles. In terms of Fig. 4, we then have to minimize the error in the C2 equations ud, + c’x + wb, = iie, + DC, + ky,
(7a)
ud, + ub, + wx = iie, + 6.y + ha,.
We therefore have the constrained least squares problem: solve (5) subject error in eqs. (7). We use standard least squares methods to obtain y=(-us,-
ua,,r,
z - u s2 - rXal l >/D
with s1 = 2( UT, - PVr2), sz = -2(
rYr, + ilrz),
r, = fie, + DC, - US, - wb, , r, = ite, + ivu, - us, - ub 23
(7b)
to minimizing
the
G. Furin
/ A modified
Clough - Tocher mterpoluni
r3 = vu3 t- wuz, u,., = 2( vz + w’), al.2
=
u 2,2 =
-quits
WC),
2( ti2 + C’))
D = -2ua,.,
The remaining
+
unknown
- u2,2u2 - a ,,,.
x is then found
from (5).
6. A comparison between both interpolants Let us now compare some properties of the interpolants from Section 4 and 5: Support. If the prescribed data values at only one data site are changed, the standard interpolant from Section 4 will change over all triangles having that data site as a vertex (Fig. 5a). The new interpolant from Section 5 has twice as many triangles in its support (Fig. 5b). This means that the effect of sharp changes in the data is distributed over a larger region. As a consequence, the new interpolant is not as much affected by ‘bad’ gradient data (e.g. from poor gradient estimation) as much as the standard interpolant. Visual impression. When we tested the interpolants on sample functions. their contour lines (computed by routines from [Petersen ‘841) looked very similar. However, when the more sophisticated tool of reflection line simulation (see next section) was employed, drastic differences became apparent, see Fig. 6. The test function was f( x, y) = (x - o.3)3 + X( 4’ - o.3)2 - 0.1X. the exact gradients were supplied. and the light source was assumed to be parallel to the y-axis. The data points consisted of the corners of the unit square with the additional two points (0.4, 0.7) and (0.6, 0.6). The reflection lines of the new scheme are close to the theoretically correct ellipses. Note that the standard scheme produces several families of concentric closed reflection lines (where one would be correct). The boundaries of the domain triangles are much more apparent in the
Fig. 5a. Support
of standard
Clough-Tocher
scheme
Fig. 5b. Support
of proposed
Clough-Tocher
scheme
G. Form / A modrfied Clough
Fig. 6a. Reflection lines of standard Clough-Tocher scheme.
Tocher mterpolunt
25
Fig. 6b. Reflection lines of proposed Clough-Tocher scheme.
reflection lines of the standard surface than in those of the one produced by the new scheme. general, the new scheme seems to produce visually more pleasing (smoother) surfaces.
In
7. Reflection lines It is often not sufficient to simply look at the 2D display of an object that has been generated by some CAGD method in order to judge its visual quality. A typical example is the curvature plot of a curve ~ it reveals curve properties on the screen that could otherwise only be detected from precision plots. For surfaces. one such ‘interrogation method’ is the use of reflection line simulation [Klass ‘801, a method that is used in the automobile industry to judge the esthetic quality of a car body. A car is brought to a show room to judge the reflection pattern of a series of parallel fuorescent light bulbs. If this reflection pattern is ‘nice’, then the esthetic quality of the car body is considered satisfactory. In the design phase, it is helpful to simulate these reflection patterns in order to detect surface irregularities before the surface is actually manufactured. Here we try to employ this concept to the surfaces described in the previous sections. Since we are only considering surfaces of the form z =f(x, y), the general (parametric) case can be substantially simplified. Let us now define how to simulate the show room environment for our surfaces: we want to plot the reflection line pattern with respect to some assumed fluorscent light bulb as perceived by some assumed observer. Thus in addition to the surface under consideration, assume a light source that is idealized by straight line parallel to the xy-plane. Also, assume an observer in a direction perpendicular to the light source. For simplicity, assume that both the observer and the light source are infinitely distant from the surface, see Fig. 7. Let us now characterize all points p on the surface that the observer perceives as reflecting light from the light source. We shall call all those points a ‘reflection line’ of the surface with respect to the assumed observer and the assumed light source. Let d = (d,, d,, 0) be a direction of unit length in the xy-plane perpendicular to the light source. Let n be the surface normal at p and let ti be the projection of n into the plane containing p and d perpendicular to the xy-plane. Then a reflection line is the collection of all points p for which the projected normal 8: fi = (nTd)d + (0, 0, ni)T
26
G. Farm
/ A modified
Clough - Tocher rnterpolunr
d
Fig. 7. Environment simulation.
for
reflection
line
i-:angle with d. This may be expressed
forms a constant
as
const. = dTii = dT((nTd)d + (0, 0, )13)-r) = dTn. Since the surface normal d,f,
is given by nT = (-f,,
+ dzf,. = const.
-f, , l), the above equations
are equivalent
to (9)
In other words, the desired reflection line is a contour of the directional derivative of f with respect to the direction d. To obtain several reflection lines, we simply contour the directional derivative surface at several values, thereby simulating several parallel light sources of different heights over the X, y-plane. Note that if f is a surface consisting of polynomial segments of degree n defined over triangles, the directional derivative of f is a piecewise polynomial surface of degree n - 1. This property is not shared by tensor product surfaces; their directional derivatives are in general polynomials of higher order, and obtaining reflection lines of them is thus computationally involved.
Acknowledgments This research was supported by DOE contract DE-AC02-82ER12046 to the University of Utah and by the University of Utah Research Committee. The paper has benefitted from discussions with the members of the Math CAGD group, especially with R.E. Barnhill, and from a suggestion by W. Bohm.
References Alfeld. P. (1985). A bivariate C2 Clough-Tocher scheme. Computer Aided Geometric Design 2, to appear. Barnhill, R.E. and Farin, G. (1981), C’ quintic interpolation over triangles: Two explicit representations, Int. J. Num. Methods in Engineering 17, 1763-1778. Biihm, W., Farin, G.. Kahmann, J. (1984), A survey of curve and surface methods in Computer Aided Geometric Design, Computer Aided Geometric Design 1, l-60. Farin. G. (1980) Bkzier polynomials over triangles and the construction of piecewise C’ polynomials, TR/91. Mathematics Department, Brunei University, Uxbridge, Middlesex, England. Farin. G. (1983). Smooth interpolation to scattered 3D data. in: R.E. Barnhill and W. Biihm. eds., Surfaces UI CAGD, North-Holland, Amsterdam, 43-63.
G. Furin
/ A modified
Clough
Tocher rnierpolunt
27
Farin. G. (1985). A survey of BernsteinBerier trtangular interpolants. in preparation. Klass. R. (1980) Correction of local surface irregularities using reflection lines, Computer Aided Design 12. 73377. Lawson. C.L. (1977). Software for C’ surface interpolation. in: J. Rice. ed.. Murhemcr/iccr/ Softwrrre 111. Academic Press, New York, 161-194. Petersen, C.S (1984). Adaptive contouring of three dimensional surfaces, Computer Aided Geometrtc Design 1, 61-70. Sahlonnibe, P. (1984) BernateinBezier methods for the construction of bivariate spline approximants, in this volume. Strang. G. and Fix. G. (1973). An Ana!\xis o/the Ftnite Element Method. Prentice-Hall, New York.
A modified
Presented
19
Design 2 (1985) 19-27
at Oberwolfach
Clough-Tocher
15 November
interpolant
1984
Abstract. Clough-Tocher
schemes are used to create surfaces that interpolate to positional and gradient values that may be arbitrarily located. A standard and a modified version of the CloughkTocher mterpolation scheme are described in term5 of Bernstein~BCzier polynomials over triangles. Both method, are compared by a method that simulates reflection lines. Keywords. Surfaces. triangular
interpolants,
Bernstein-Bkzier
polynomials.
reflectlon
line\
1. Introduction Clough-Tocher interpolants were originally invented as a tool for the finite element method [Strang, Fix ‘731, but they have also turned out to be a useful tool in CAGD, namely in the field. of scattered data interpolation [Lawson ‘771. A Clough-Tocher interpolant produces a C’ piecewise polynomial surface, defined over a triangulation of given data sites. Recently the original interpolant has been generalized [Alfeld ‘85. Sablonniere ‘851, but in this paper we present a modification rather than a generalization. It is based on increasing the smoothness of the overall interpolant by an appropriate specification of the available degrees of freedom. The relevant smoothness conditions are conveniently formulated in terms of Bernstein-Bezier polynomials defined over triangles. The modified method is then compared to the standard method, and to this end a (surface interrogation) method is introduced that simulates ‘reflection lines’.
2. Continuity
conditions
for triangular BCzier patches
For the theory of BernsteinBCzier polynomials defined over triangles. the reader is referred to [Farin ‘80,‘83,‘85,BBhm/Farin/Kahmann ‘841. An nth degree Bernstein-Bezier polynomial defined over a triangle is of the form
,+J+h=" where U, u, M:are barycentric coordinates of the domain triangle and the b,,,,l, are called Bezier ordinates of p, see Fig. 1. The Beziei ordinates form a triangular net, the socalled “control net”. The Bezier ordinates of a boundary curve are the corresponding boundary Bezier points of the net. Similarly, a cross-boundary derivative, evaluated along that boundary, depends only on the first two rows of BCzier ordinates parallel to that boundary. Consider now two adjacent triangles and a Bezier net defined over each of them. We ask for the conditions that the nets must satisfy in order for the two patches defined by them to be C’ or C2. 0167-8396/M/$3.30
8) 1985, Elsewer Science Publishers
B.V. (North-Holland)
Fig. 1. Bernstein-Bkzier one triangle.
Fig. 2a. Adjacent
triangular
Fig. 2b. Adjacent
patches
triangular
wjith C’ continuity
patches
with C2 continuity
polynomial
over
G. Fmn
VLL
/ A modified Clough - Tocher inrerpoiani
Fig. 3. Three triangle.
21
C’ control
nets over split
For the C’ case, Fig. 2a illustrates the answer: The condition for C’ continuity is that the shaded pairs of triangles he coplunur (Farin ‘80). For C2 continuity, the shaded pairs of triangles in Fig. 2b are constructed to be coplanar. Each of the marked points will in general correspond to two different z-values. two from each control net. The condition for C2 continuity is thut these values coincide. These ‘auxiliary points’ are analogous to the auxiliary points used to define C2 conditions for Bezier curves [Bohm/ Farin/Kahmann ‘841. Those points could also be interpreted as B-spline coefficients. It is an interesting (and hence open) question whether the auxiliary points from Fig. 2 play a similar role for triangular splines.
3. A parameter count As an application of the above continuity conditions, we consider three Bernstein-Bezier triangles as in Fig. 3. We require any two adjacent patches to be C’; i.e. we consider C’ piecewise cubits defined over a triangulation that is obtained by splitting a triangle (called ‘macro-triangle’) into three triangles (called ‘mini-triangles’) at its centroid. These C’ cubits form a linear space, and we are interested in its dimension (i.e. the degrees of freedom present in the three C’ BCzier nets from Fig. 3). To count these degrees of freedom, one observes that the centroid and the three Bezier points next to it must be coplanar, thus providing three degrees of freedom. Moving out from the centroid one more step, we find six more Bezier points. Taking into account that they are related by the C’ conditions. one finds that they add three more degrees of freedom. Moving one step outwards again, we now have to consider the nine boundary Bezier points. They provide six additional degrees of freedom; the linear space of C’ cubits ooer the split triangle therefore has dimension twelve. The above counting procedure also leads to an interesting result: C’ cubits ouer the split triangle are C2 at the centroid [Farin ‘80, Alfeld ‘851. For a proof, strip the piecewise cubits from Fig. 3 of the boundary Bezier points, leaving a C’ piecewise quadratic function. The above counting of degrees of freedom asserted that such a function has six degrees of freedom. This is the dimension of the linear space of (global) quadratics defined over the macro-triangle.
G. Farm / A modified Clough - Tocher inierpohnr
22
They form a subspace of C’ piecewise quadratics. Since both spaces are of dimension six, they must be equal. Thus C’ piecewise quadratics over the split triangle are actually global quadratics, i.e. instead of being just C’, they are also C2. The considered piecewise quadratic and the original piecewise cubic agree in all derivatives up to order two at the centroid. Hence C’ cubits over the split triangle are actually C2 at the centroid.
4. The Clough-Tocher
interpolant
Suppose z-values and gradients are given over a set of triangulated data points, and one wants a C’ piecewise cubic surface that interpolates to this data. As an ‘initial guess’, the following piecewise cubic interpolation scheme is useful: for each triangle, define a cubic Bernstein-Bezier cubic. Its nine boundary Bezier ordinates (all the h,,,,I except 6,,,,,) are readily determined from the given data by univariate cubic Hermite interpolation. The remaining term h,,,,, is given by b 1.1.1= i(bz.o,l + b,,,., + b,.,,, + bo.l.2 +
b2.1.0
+
b,.,.,,)
-Theo
+
h,,.,
+
h.u.3).
(2)
the nine-parumeter interThis choice of b,,,,, ensures quadratic precision of this interpolant, [Farin ‘831. A surface composed from nine-parameter interpolants will however fail to be C’ in general. It must be modified to ensure the desired overall C’ smoothness. This is done by splitting each triangle in the triangulation into three mini-triangles and consequently splitting each nineparameter cubic into three ‘mini-cubits’ (for the algorithm to do this, see, e.g., [Bohm/Farin/ Kahmann ‘841). The subdivided surface now has enough degrees of freedom (twelve per subdivided macro-triangle instead of ten before the subdivision) to allow C’ continuity of the overall interpolant. In Fig. 4, let C,, P2, P, and C,, P,, P2 be the vertices of two adjacent microtriangles, coordinates of the opposite respectively. We can express C, and C, in terms of barycentric triangle:
polunt
c,=K2+6P,+ivP2,
(3)
c, = UC, + UP2+ WP, . While the C’ continuity
conditions
(4) from Section
2 are fulfilled
for the triangle
pairs formed
Pl
p2
Fig. 4. Adjacent C’ triangular defined over mini-triangles.
patches
23
G. Furin / A modified Clough - Tocher rnrerpolunt
by al, b,, subtriangles.
u2.
cl and
a3.
b,,
The C’ condition
by the middle a4. c2, they are in general violated that has to be met can be expressed as
pair
J = vu3 + wa2 + ux.
of
(5)
The standard way to guarantee (5) is as follows: let L denote any direction not parallel to P,. P2. Then the directional derivative of the mini-cubic p, defined over C,, P2. P, is a univariate quadratic Bezier polynomial with Bezier ordinates 3(l,b,
+ j2u2 + [,a,),
3(1,x + 12a, + /,a,).
3(l,bz + /,a4 + laj).
Here, I,, I,, I, denote the barycentric representation of the direction L in terms of the triangle C,, P2. P,. One way to fix the unknown x is to demand that the quadratic directional derivative be linear instead of quadratic. This condition can be expressed as [see Bohm/Farin/ Kahmann ‘841: (I,b,
+ 12a2 + /,a,)
+ 2( I,x + I,u,
+ /,a,)
+ (1,b,
+ /,a,
+ /,a,)
= 0.
(6)
It is now possible to determine x from (6). The remaining unknown y may be found from (5) or by an analogous procedure for the mini-cubic p2, defined over triangle C,, P,, P2. It is essential that L denote the same direction for both p, and p2. The standard way to achieve this is to make L perpendicular to P,P,. The obvious drawback is that the resulting interpolant is not affinely invariant. An alternative is to set L = Cl - C,, but this does not seem to have been considered in the literature. The above method is then applied to all pairs of adjacent macro-triangles. Returning to Fig. 3, we see that we thus determine all Bezier ordinates marked by full circles. The remaining ones, marked by open circles, have to be computed in order to ensure internal C’ continuity. This is accomplished by applying the Cl-conditions four times [see e.g. Farin ‘831.
5. A smoother interpolant The ‘linear cross boundary derivative’ condition (6) appears to be somewhat arbitrary, and in this section we will develop a geometrically more meaningful alternative. The C’ condition (5) is not sufficient to determine both x and JJ, and consequently an additional condition is needed. We have already (in Section 3) observed that C’ cubits over the split triangle enjoy an extra degree of continuity at the centroid, and thus define a surface that is smoother than could be expected ab initio. Therefore it seems reasonable to come as close to C2 continuity as possible across the edges of adjacent macro-triangles. In terms of Fig. 4, we then have to minimize the error in the C2 equations ud, + c’x + wb, = iie, + DC, + ky,
(7a)
ud, + ub, + wx = iie, + 6.y + ha,.
We therefore have the constrained least squares problem: solve (5) subject error in eqs. (7). We use standard least squares methods to obtain y=(-us,-
ua,,r,
z - u s2 - rXal l >/D
with s1 = 2( UT, - PVr2), sz = -2(
rYr, + ilrz),
r, = fie, + DC, - US, - wb, , r, = ite, + ivu, - us, - ub 23
(7b)
to minimizing
the
G. Furin
/ A modified
Clough - Tocher mterpoluni
r3 = vu3 t- wuz, u,., = 2( vz + w’), al.2
=
u 2,2 =
-quits
WC),
2( ti2 + C’))
D = -2ua,.,
The remaining
+
unknown
- u2,2u2 - a ,,,.
x is then found
from (5).
6. A comparison between both interpolants Let us now compare some properties of the interpolants from Section 4 and 5: Support. If the prescribed data values at only one data site are changed, the standard interpolant from Section 4 will change over all triangles having that data site as a vertex (Fig. 5a). The new interpolant from Section 5 has twice as many triangles in its support (Fig. 5b). This means that the effect of sharp changes in the data is distributed over a larger region. As a consequence, the new interpolant is not as much affected by ‘bad’ gradient data (e.g. from poor gradient estimation) as much as the standard interpolant. Visual impression. When we tested the interpolants on sample functions. their contour lines (computed by routines from [Petersen ‘841) looked very similar. However, when the more sophisticated tool of reflection line simulation (see next section) was employed, drastic differences became apparent, see Fig. 6. The test function was f( x, y) = (x - o.3)3 + X( 4’ - o.3)2 - 0.1X. the exact gradients were supplied. and the light source was assumed to be parallel to the y-axis. The data points consisted of the corners of the unit square with the additional two points (0.4, 0.7) and (0.6, 0.6). The reflection lines of the new scheme are close to the theoretically correct ellipses. Note that the standard scheme produces several families of concentric closed reflection lines (where one would be correct). The boundaries of the domain triangles are much more apparent in the
Fig. 5a. Support
of standard
Clough-Tocher
scheme
Fig. 5b. Support
of proposed
Clough-Tocher
scheme
G. Form / A modrfied Clough
Fig. 6a. Reflection lines of standard Clough-Tocher scheme.
Tocher mterpolunt
25
Fig. 6b. Reflection lines of proposed Clough-Tocher scheme.
reflection lines of the standard surface than in those of the one produced by the new scheme. general, the new scheme seems to produce visually more pleasing (smoother) surfaces.
In
7. Reflection lines It is often not sufficient to simply look at the 2D display of an object that has been generated by some CAGD method in order to judge its visual quality. A typical example is the curvature plot of a curve ~ it reveals curve properties on the screen that could otherwise only be detected from precision plots. For surfaces. one such ‘interrogation method’ is the use of reflection line simulation [Klass ‘801, a method that is used in the automobile industry to judge the esthetic quality of a car body. A car is brought to a show room to judge the reflection pattern of a series of parallel fuorescent light bulbs. If this reflection pattern is ‘nice’, then the esthetic quality of the car body is considered satisfactory. In the design phase, it is helpful to simulate these reflection patterns in order to detect surface irregularities before the surface is actually manufactured. Here we try to employ this concept to the surfaces described in the previous sections. Since we are only considering surfaces of the form z =f(x, y), the general (parametric) case can be substantially simplified. Let us now define how to simulate the show room environment for our surfaces: we want to plot the reflection line pattern with respect to some assumed fluorscent light bulb as perceived by some assumed observer. Thus in addition to the surface under consideration, assume a light source that is idealized by straight line parallel to the xy-plane. Also, assume an observer in a direction perpendicular to the light source. For simplicity, assume that both the observer and the light source are infinitely distant from the surface, see Fig. 7. Let us now characterize all points p on the surface that the observer perceives as reflecting light from the light source. We shall call all those points a ‘reflection line’ of the surface with respect to the assumed observer and the assumed light source. Let d = (d,, d,, 0) be a direction of unit length in the xy-plane perpendicular to the light source. Let n be the surface normal at p and let ti be the projection of n into the plane containing p and d perpendicular to the xy-plane. Then a reflection line is the collection of all points p for which the projected normal 8: fi = (nTd)d + (0, 0, ni)T
26
G. Farm
/ A modified
Clough - Tocher rnterpolunr
d
Fig. 7. Environment simulation.
for
reflection
line
i-:angle with d. This may be expressed
forms a constant
as
const. = dTii = dT((nTd)d + (0, 0, )13)-r) = dTn. Since the surface normal d,f,
is given by nT = (-f,,
+ dzf,. = const.
-f, , l), the above equations
are equivalent
to (9)
In other words, the desired reflection line is a contour of the directional derivative of f with respect to the direction d. To obtain several reflection lines, we simply contour the directional derivative surface at several values, thereby simulating several parallel light sources of different heights over the X, y-plane. Note that if f is a surface consisting of polynomial segments of degree n defined over triangles, the directional derivative of f is a piecewise polynomial surface of degree n - 1. This property is not shared by tensor product surfaces; their directional derivatives are in general polynomials of higher order, and obtaining reflection lines of them is thus computationally involved.
Acknowledgments This research was supported by DOE contract DE-AC02-82ER12046 to the University of Utah and by the University of Utah Research Committee. The paper has benefitted from discussions with the members of the Math CAGD group, especially with R.E. Barnhill, and from a suggestion by W. Bohm.
References Alfeld. P. (1985). A bivariate C2 Clough-Tocher scheme. Computer Aided Geometric Design 2, to appear. Barnhill, R.E. and Farin, G. (1981), C’ quintic interpolation over triangles: Two explicit representations, Int. J. Num. Methods in Engineering 17, 1763-1778. Biihm, W., Farin, G.. Kahmann, J. (1984), A survey of curve and surface methods in Computer Aided Geometric Design, Computer Aided Geometric Design 1, l-60. Farin. G. (1980) Bkzier polynomials over triangles and the construction of piecewise C’ polynomials, TR/91. Mathematics Department, Brunei University, Uxbridge, Middlesex, England. Farin. G. (1983). Smooth interpolation to scattered 3D data. in: R.E. Barnhill and W. Biihm. eds., Surfaces UI CAGD, North-Holland, Amsterdam, 43-63.
G. Furin
/ A modified
Clough
Tocher rnierpolunt
27
Farin. G. (1985). A survey of BernsteinBerier trtangular interpolants. in preparation. Klass. R. (1980) Correction of local surface irregularities using reflection lines, Computer Aided Design 12. 73377. Lawson. C.L. (1977). Software for C’ surface interpolation. in: J. Rice. ed.. Murhemcr/iccr/ Softwrrre 111. Academic Press, New York, 161-194. Petersen, C.S (1984). Adaptive contouring of three dimensional surfaces, Computer Aided Geometrtc Design 1, 61-70. Sahlonnibe, P. (1984) BernateinBezier methods for the construction of bivariate spline approximants, in this volume. Strang. G. and Fix. G. (1973). An Ana!\xis o/the Ftnite Element Method. Prentice-Hall, New York.