Computing optimal triangulations using simulated annealing

Computing optimal triangulations using simulated annealing

Computer Aided North-Holland COMAID Geometric Design 329 10 (1993) 329-345 337 Computing optimal triangulations using simulated annealing Larry...

2MB Sizes 0 Downloads 98 Views

Computer Aided North-Holland

COMAID

Geometric

Design

329

10 (1993) 329-345

337

Computing optimal triangulations using simulated annealing Larry L. Schumaker

*

Department of Mathematics, Vanderbilt UniL>ersity,Nash&k, Presented Received

in Oberwolfach, February 1993

June

TN 37240, USA

1992

Abstract Schumaker, L.L., Computing Design 10 (1993) 329-345.

optimal

triangulations

using simulated

annealing,

Computer

Aided

Geometric

Triangulations play an important role in approximation, CAGD, numerical analysis, and elsewhere. In this paper we are concerned with the problem of constructing triangulations which are optimal in some sense. Our aim is to show how simulated annealing can be used to search for globally optimal triangulations for a wide class of optimality criteria. We also give several examples to illustrate its performance on a variety of problems of interest in CAGD and surface fitting.

1. Introduction Definition 1.1. A collection A = {Tli!, of triangles in the plane is called a triangulation of a set R with polygonal boundary provided that 1. any pair of triangles from A intersect at most at a common vertex or along a common edge, 2. the union R of the triangles of A by I/= {q = (x,, yJ);.

{7;}? is a connected

set in iw*. We denote

the set of vertices

Given 0 and V, in general there may be many different associated triangulations. It is well-known, however, that all triangulations have the same number of triangles and edges (see Remark 8.2). Triangulations play an important role in approximation theory, computer-aided geometric design, and numerical analysis, see e.g. [Schumaker ‘87, ‘931. They are particularly useful as tools for defining spline spaces consisting of piecewise polynomials defined over a triangulation. It is well-known (cf. [Baszenski & Schumaker ‘91; Dyn et al. ‘92; Dyn et al. ‘90a, ‘90b, Quak & Schumaker ‘89, ‘90a, ‘90b, ‘91; Rippa ‘89a, ‘89b]) that the choice of triangulation plays an important role in working with splines, and that it usually pays to try to find a triangulation which is best in some sense. The standard approach to finding best triangulations is to apply an appropriate edge swapping algorithm (see Algorithm 2.3). Unfortunately, for all but one of the commonly used optimality criteria, edge swapping leads to triangulations which are only locally optimal and Correspondence to: L.L. Schumaker, Department of Mathematics, Vanderbilt University, * Supported in part by the National Science Foundation under Grant DMS-9208413. 0167-8396/93/$06.00

0 1993 - Elsevier

Science

Publishers

B.V. All rights reserved

Nashville,

TN 37240, USA.

330

L.L. Schumaker

/ Computing optimal triangulations using simulated annealing

not globally optimal. In this paper we show how simulated annealing can be used to search for globally optimal triangulations for a wide class of optimality criteria. The paper is organized as follows. In Section 2 we discuss a general class of criteria for comparing triangulations, and present a local swap algorithm which can be used with these criteria to adjust a given triangulation. In Section 3 we discuss a version of the classical simulated annealing algorithm which can be applied to the problem of finding optimal triangulations. Sections 4-7 are devoted to applications. In particular, we discuss geometric criteria (the max-min angle and min-max angle criteria) in Sections 4 and 5, data dependent interpolation using linear splines in Section 6, and data dependent least-squares fitting in Section 7. We conclude the paper with remarks and references.

2. Optimal

triangulations

Given a set of vertices I/ and a set 0 with polygonal boundary as in Definition 1.1, let Y be the set of all possible triangulations of J2 with vertices from I/. In this section we discuss a general approach to comparing the various triangulations in 7. Definition 2.1. Let %Ybe a mapping of 7 x 7 into R. Then we say that 1. the triangulation A is better than the triangulation A with respect to 59 provided that @A, A) < 0; 2. the two triangulations are equally good if %?‘(A, A) = 0; 3. A is an optimal triangulation in 7 with respect to 59 provided that ‘Z?(A) A) < 0 for all triangulations A in 7. Since the set 7 is finite, it is clear that for a prescribed E”, there always exists at least one optimal triangulation. For most criteria used in practice, optimal triangulations are not unique. We now introduce a local form of optimulity. Given any interior edge e of a triangulation A, there is a natural associated quadrilateral Q, formed from the two triangles of A which share e. We say that such an edge is swappuble provided that Q, is convex and no three of its vertices are collinear. If an edge e of a triangulation A is swappable, then we can create a new triangulation by actually swapping the edge. In particular, if the quadrilateral Q, has the four vertices uil, ui2, ui3, ui4 in counterclockwise order and e has endpoints oil, ui3, then the swapped edge would have endpoints vi2, vi4. Definition 2.2. Let E’ be a criterion for comparing triangulations. We say that a triangulation A is locally optimal with respect to g provided that it cannot be improved by swapping any swappable edge e of A. It is clear that every optimal triangulation must be locally optimal. can be used to construct locally optimal triangulations. Algorithm 2.3 (Local Swap Algorithm). do until no swap is possible choose a swappable edge e if swapping edge e improves the triangulation,

The following

algorithm

make the swap.

In implementing this algorithm, it is often useful to sort the edges according the associated pair of triangles are with respect to the optimality criteria. Then

to how good at each step,

L.L. Schumaker

/ Computing

optimal triangulations

using simulated

annealing

331

the edges are examined in order until one is found which can be swapped to improve the triangulation. If none can be found, the algorithm stops. Clearly, since there are only a finite number of possible triangulations in the set 7 and the local swap algorithm makes a swap only when it improves the triangulation, the algorithm must terminate in a finite number of steps. It always produces a locally optimal triangulation which is guaranteed to be at least as good if not better than the starting triangulation. In general, locally optimal triangulations are not unique. Indeed, for any given criterion 2?, it is easy to construct an example of a triangulation which is locally optimal, but which contains an edge which can be swapped without changing the value of the criterion. In this connection we have the following definition. Definition 2.4. Suppose A is a triangulation which is obtained We say that the swap is neutral provided that %‘(A, A> = 0.

from

A by swapping

an edge.

Another problem with using the local swap algorithm is that generally it is impossible to determine whether a given triangulation is globally optimal or not (the one exceptional case is when we use the max-min angle criterion discussed in Section 4 below-for that criterion, a locally optimal triangulation is actually globally optimal). This suggests the need for an algorithm to search through the set of all possible triangulations looking for a better triangulation. Generally 7 is very large, and so a well-structured algorithm is needed.

3. Simulated

annealing

Simulated annealing (see e.g. [Press et al. ‘861) is a technique which was developed to help solve large combinatorial optimization problems. It can be applied to problems where one can reach all feasible points by performing a sequence of steps, and where for each step, the amount of change in the quantity to be minimized can be quantified. It is based on probabilistic methods, and is designed to avoid getting stuck at local (non-global) minima. At each stage of the process, both good steps (those which improve some measure of quality) and bad steps (those that do not) are made. The bad steps are chosen randomly, and as one proceeds from stage to stage, the probability that they are made is slowly reduced according to a prescribed annealing schedule. To describe the simulated annealing process as it applies to finding (globally) optimal triangulations as defined in Definition 1.1, suppose that ‘$7 is our criterion for comparing triangulations. Let ntemps, nlimit, and glimit be positive integers, and let t, > t, > . . . > t nt~mps> 0 be given real numbers. Algorithm 3.1 (Simulated Annealing). do k = 1 to ntemps do I = 1 to nlimit while the number of good swaps < glimit choose a random edge e in the triangulation if the edge e can be swapped let A be the current triangulation and A the swapped if d = %?(A, A) < 0, make the swap else choose a random number 0 < f3 < 1 if 8 G epdltk, make the swap endif endif

one

332

L.L. Schumaker / Compuhg optimaltriangulationsusing simulated annealing

Fig. 1. Probability

of making

bad swaps.

The parameter ntemps controls the number of stages in the algorithm. The annealing controls the probability that bad swaps will be made. In schedule t, > t, > . * . > t,,,,p, practice we choose the random numbers 0 from a uniform distribution on [O,l]. Then the probability of making a bad swap at temperature t, for a given value d is Pr(B G e-d/tk) = epdltk, see Fig. 1. The t, are analogous to temperatures in the annealing of metal, and as they get smaller, it becomes less and less likely that bad swaps will occur. The parameter nlimit controls the number of swaps to be attempted at each stage of the algorithm. Finally, the parameter glimit controls the number of good swaps allowed at each temperature. The idea is that if one makes too many good swaps too soon, the algorithm will not find its way out of false valleys and into the one where the global minimum is to be found. While this simulated annealing algorithm provides a good way of searching for a globally optimal triangulation, in general, we cannot guarantee that it will find one. Indeed, as remarked in Section 2, in general we have no way of recognizing an optimal triangulation. Thus we have to proceed as follows: 1. we choose some set of parameters, and run the algorithm, keeping track of the best triangulation found during the run, 2. if we are not satisfied with the result, we adjust the parameters, and run the algorithm again. It is well-known that the behavior of the simulated annealing algorithm is fairly sensitive to the choice of parameters (and in particular, to the choice of the annealing schedule). Based on a large number of simulations involving optimal triangulation, we can make the following general suggestions concerning the choice of parameters for this application: 1. To determine a reasonable value for the mitial temperature, run the local swap algorithm, and keep track of the changes d = ‘%?A, A) that correspond to good swaps. Choose the initial temperature t,, to be about twice the largest of these values. 2. Choose t, = rkt, for k = 1, 2,. . . , where 0 < r < 1 is fixed. This gives a linear annealing schedule where the temperature is reduced by the factor r at each step of the algorithm. A reasonable choice for r is 0.9 or 0.95.

L.L. Schumaker

/ Computing optimal triangulations using simulated annealing

333

3. Choose nlimit and glimit as fixed multiples of the number of edges in the triangulation. A factor of 5 or 10 seems to work well in practice. One way to monitor whether the selection of parameters is reasonable is to look for some kind of “convergence”. Assuming that the rk in the annealing schedule are becoming smaller and smaller, as the annealing progresses it becomes more and more difficult to make bad swaps (and so consequently fewer and fewer good swaps are made). Since d can take on only a finite number of values, when tk is sufficiently small, bad swaps become virtually impossible. Of course, no matter how small t, is, it is always theoretically possible that another bad swap will be made if enough edges are tried. However, in practice, at some point the algorithm stops making swaps-in this case we can say that it has “converged”. Clearly, if the initial temperature is high, the algorithm will need to run longer before converging. On the other hand, a high initial temperature and large values of nlimit and ngood allow the exploration of a wider class of possible triangulations, increasing the chance of finding the global minimum. If the initial temperature is taken to be low, the algorithm will converge quickly, and often will converge to a local optimum rather than a global one. With some user interaction, a suitable initial temperature can usually be found in just a few tries.

4. The max-min

angle criterion

The classic criterion

for comparing

triangulations

is the following

Definition 4.1 (max-min angle criterion). For each triangle T in a triangulation A, let m(T) be m(T). We define @A, 6) = m(A) - m(A). the smallest angle in T, and let m(A) = min,,, Under this criterion, an optimal triangulation is one whose smallest angle is maximal. It is well-known (cf. [Lawson ‘77, Schumaker ‘871) that a triangulation is optimal with respect to the max-min angle criterion if an only if it is a Delaunay triangulation. For given V and 0, there may be more than one Delaunay triangulation. Figure 2 shows a simple example involving four points lying on a circle. In this case there are two different Delaunay triangulations. For both triangulations, the minimum angle is 7 degrees. On the other hand, it is also known that one can always convert any particular Delaunay triangulation corresponding to vertices I/ and polygonal domain fi into any other by a finite series of neutral swaps. To compute an optimal triangulation, we may start with an arbitrary triangulation and apply the Local Swap Algorithm 2.3. In this case, the local swap criterion is equivalent to the following well-known local circle criterion: Suppose (2, is a convex quadrilateral with the four vertices L’;,, ui2, c’,~, ui4 in counterclockwise order, formed by two triangles sharing the edge with endpoints L’~, and ui3. Let C be the circle passing through the first three vertices. Then the max-min angle criterion is satisfied by the quadrilateral if and only if ui4 is outside (or on the boundary) of C.

(a)

Fig. 2. Two Delaunay

triangulations

of four points.

L.L. Schumaker

334

/ Computing optimal triangulations using simulated annealing

It is this observation which provides the connection between max-min angle optimal triangulations and Delaunay triangulations. It also implies (cf. [Lawson ‘77, Schumaker ‘871) that (with respect to this criterion) a locally optimal triangulation (as produced by the local swap algorithm) is in fact optimal. There are of course many other other algorithms which can be used to compute a Delaunay triangulation, some of which may be much faster than Algorithm 2.3, see e.g. [Schumaker ‘871 for a discussion and references. We now consider a more refined way of comparing triangulations in terms of minimum angles. Definition 4.2 (Lexicographical max-min be as in Definition 4.1, and let

angle criterion). For each triangle TEA,

a(A) = (~(A),...,Q(A))

let m(T)

(I)

be the vector obtained by arranging the set {m(T)},,, in increasing order. If (Y(A)= a(A), we define %‘(A, A) = 0. Otherwise, there exists a p (depending on A and A) such that q(A)

=a&,

i= l,...,

p-

1,

but

a,(A)

#(Y@),

(2)

and we define ‘@(A, A) =cx,(A)

-q,(A).

(3)

It is obvious that the lexicographical max-min angle criterion is more discriminating than the classical max-min angle criterion. For example, both of the triangulations in Fig. 2 are optimal with respect to the max-min angle criterion, but not with respect to the lexicographical version of the criterion. Indeed, the minimum angles (in degrees) associated with the two triangles in Fig. 2(a) are (7, 71, while for the triangles in Fig. 2(b) they are (7, 22). Thus the second triangulation is better with respect to the lexicographical max-min angle criterion. The local-swap Algorithm 2.3 is well-suited to working with the lexicographical criterion, since swapping an edge always changes exactly two components in the vector, and the new vector is at least as good as the old lexicographically. It is easy to see that Algorithm 2.3 always produces a triangulation which is globally optimal with respect to the lexicographical max-min angle criterion (indeed, it must be a Delaunay triangulation with all neutral swaps locally improved). Thus, with this criterion there is no need for the annealing algorithm.

5. The min-max

angle criterion

The following alternative to the max-min angle criterion has been suggested & Little ‘841:

in [Barnhill

Definition 5.1 (min-max angle criterion). For each triangle T in a triangulation Ai let M(T) be the largest angle in T, and define M(A) = maxr E B M(T). Define HA, AI = M(A) - M(A). Under this criterion, an optimal triangulation is one whose largest angle is minimal. The lexicographical version of this criterion is Definition

5.2 (Lexicographical

min-max angle criterion). For each triangle T E A let M(T) be

the largest angle in the triangle T, and let a(A) = (q(T),...,aN(T))

(4)

L.L. Schumaker

/ Computing

optimal triangulations

using simulated

annealing

335

in decreasing order. If (Y(A) = a(A), be the vector obtained by arranging the set {M(T)},,, we define $%?A,A> = 0. Otherwise, there exists a p (depending on A and A) such that ai

=aJA),

i= l,...,p-1,

but

a,(A)

#c&),

(5)

and we define %‘(A, A) =&)

-a,(A).

With respect to these min-max angle criteria, a locally optimal triangulation necessarily globally optimal (see [Nielson ‘87, Schumaker ‘871 and Example 5.3.). This that using the local swap algorithm, we may find a local optimum which is not a optimum. We now illustrate how simulated annealing can be used to search for a optimum.

(6) is not means global global

Example 5.3. Let I/ be the following set of six points: (0, 0.3), (0.14, 0.04), (0.25, 0.65), (0.26, 0.55), (0.45, O), (0.58, 0.2). Let 0 be the convex hull of these points.

This example is taken from [Nielson ‘871. The Delaunay triangulation of this set of six points is shown in Fig. 3(a). The maximum angles for each of the six triangles (sorted according to size) are (143.27, 128.17, 115.43, 103.90, 67.97). The largest angle in this triangulation is 143.27 degrees. Fig. 3(b) shows the optimal triangulation with respect to the min-max angle criterion. The corresponding vector of maximum angles is (143.27, 128.17, 104.54, 103.90, 85.42), and the largest angle in this triangulation is also 143.27 degrees. However, lexicographically, the second triangulation is better since the third entry is 104.54 as compared to 115.43. The min-max angle optimal triangulation shown in Fig. 3(b) can be found by applying the local swap algorithm starting with the Delaunay triangulation in Fig. 3(a). Only one swap is needed. On the other hand, the local swap algorithm will not produce the optimal triangulation if we start with the triangulation shown in Fig. 3(c). Indeed, as is easily verified, this triangulation is already locally optimal (swapping any edge leads to a worse triangulation). To illustrate how simulated annealing performs, we applied it to this example starting with the triangulation in Fig. 3(c). Instead of working with the actual angles, it is equivalent to work with their cosines. To minimize the maximum angle, we need to maximize cosine values. Table 1 shows a typical annealing run starting with an initial temperature of 0.1 and reducing it by a factor of 0.95 at each step (remember in this case temperature corresponds to changes in cosine values). We choose nlimit = glimit = 50 and ntemp = 6. The first column shows the temperature at each step, while the second and third columns show the number of good and bad swaps made at that temperature. The fact that at temperatures 0.081 and 0.077 neither good nor bad swaps were made suggests that the algorithm has converged in the sense discussed at the end of Section 3. As can be seen from Table 1, the algorithm made a total of 6 swaps. At this point, there is no guarantee that the resulting triangulation is actually min-max angle optimal (although in fact it is, as can be checked by examining all of the 10 possible triangulations). To increase our confidence that we have found an optimal triangulation, we can rerun the algorithm with different parameters. Table 2 shows what happens if we start with an initial temperature of 0.2 and use nlimit = glimit = 50 and ntemp = 40. This time it took 33 temperature steps before the algorithm converges, involving a total of 220 swaps. The result was again the triangulation in Fig. 2(b). The reason for the large number of steps required for convergence in this case is that the initial temperature has been chosen relatively large.

336

L.L. Schurnaker / Computing optimal triangulations using simulated annealing

Fig. 3. Alternate

triangulations

of the 6 points in Example

5.3.

If the initial temperature is set too low, annealing may not be able to make any swaps at all. For example, with an initial temperature of 0.01 and nlimit = glimit = 50, annealing does not change the initial triangulation shown in Fig. 3(c) (and so does not deliver the optimal triangulation). Example

5.4. Let I/ be the set of 30 points shown in Fig. 4(a).

The Delaunay triangulation associated with this set of points is shown in Fig. 4(a). Applying the local swap algorithm results in a total of 4 swaps. The resulting change in

Table 1 An annealing

run on Example

5.3

temp

ngood

nbad

0.100 0.095 0.090 0.086 0.081 0.077

0 1 0 2 0 0

1 0 1 1 0 0

L.L. Schumaker Table 2 An annealing

run on Example

/ Computing opiimal triangulations using simulated annealing

337

5.3

temp

ngood

nbad

temp

ngood

nbad

0.200 0.190 0.181 0.171 0.163 0.155 0.147 0.140 0.133 0.126 0.120 0.114 0.108 0.103 0.098 0.093 0.088 0.084 0.079 0.075

7 3 7 5 6 5 8 3 5 5 6 2 6 4 4 1 2 0 3 4

7 6 5 5 6 5 8 4 5 4 6 3 5 4 4 1 2 1 2 4

0.072 0.068 0.065 0.061 0.058 0.055 0.053 0.050 0.048 0.045 0.043 0.041 0.039 0.037 0.035 0.033 0.032 0.030 0.028 0.027

5 1 4 3 1 4 2 1 0 0 1 2 1 0 0 0 0 0 0 0

4 1 4 3 1 4 2 1 0 0 1 2 1 0 0 0 0 0 0 0

maximum angles can be seen in Table 3 The largest 8 angles were unchanged, and the ninth largest was reduced from 151.76 to 148.98 degrees. To explore how good the triangulation produced by the local swap algorithm for Example 5.4 is, I ran the simulated annealing algorithm, starting with the Delaunay triangulation. The third set of angles shown in Table 3 result from an annealing run with initial temperature 0.5 and nlitnit = glimit = 800, ntemp = 70. This required a total of 10,150 swaps. The resulting triangulation is an improvement over that delivered by the local swap algorithm in as much as the largest 10 angles are the same, but the 11th largest angle is reduced from 132.61 to 130.38 degrees. Further annealing runs did not produce any better triangulation. This example illustrates the fact that, in general, there is not too much difference between an optimal max-min angle (Delaunay) triangulation and an optimal min-max angle triangulation. (b

Fig. 4. Alternate

triangulations

of the 30 points in Example

5.4.

338

L.L. Schumaker

/ Computing optimal triangulations using simulated annealing

Table 3 Maximum angles for Example 5.4 Angles for the Delaunay triangulation (Fig. 4(a)) 173.70 155.21 125.67 117.08 102.24 94.49 85.05 76.29

168.44 151.76 123.82 111.30 101.84 93.84 83.04 75.65

167.45 150.76 122.02 107.38 101.62 90.61 82.58 74.23

163.96 135.63 119.59 106.81 100.63 90.27 81.51 71.82

161.44 132.61 118.88 106.40 97.32 89.85 81.11 70.59

159.75 130.38 118.32 105.70 96.81 86.38 80.47

158.50 126.23 118.04 103.01 94.91 85.43 76.72

Angles after local swap algorithm applied (Fig. 4(b)) 173.70 155.27 123.32 111.30 101.84 93.84 85.43 80.47

168.44 150.76 122.95 108.28 101.62 93.32 85.05 76.29

167.45 135.63 122.02 107.38 100.63 93.16 83.14 76.65

163.96 132.61 119.59 106.40 97.32 90.61 83.04 74.23

161.44 130.38 118.88 103.01 96.81 90.27 82.58 70.59

159.75 126.23 118.04 102.81 94.91 89.85 81.51

158.50 125.67 114.08 102.24 94.49 86.38 81.11

167.45 135.63 118.88 106.81 96.81 90.61 83.04 75.65

163.96 130.38 118.64 106.40 94.91 90.27 82.58 14.23

161.44 126.23 118.04 103.01 94.49 89.85 81.85 70.59

159.75 123.32 144.08 102.24 93.84 86.38 81.51

158.50 122.95 111.30 101.84 93.16 85.43 81.11

Angles after annealing 173.70 155.27 122.02 108.28 100.63 92.94 85.05 80.47

168.44 150.76 119.59 107.38 97.32 91.44 83.14 76.29

6. Data-dependent

interpolation

Suppose we are given a set of points V= ((xi, y,>}r=r in R2 and measurements zi =f(x;, yi>, i=l >*.., n, on a (possibly unknown) function f defined on the convex hull R of I/. In this section we consider the problem of constructing a linear spline surface which approximates f. Given a triangulation A = {T,, . . . , TN} of 0 with vertices at the points of V, we define the associated space of linear splines by S,O(A)=(SEC~(.):,~,~,,~,~=,

,..., N},

(7)

where 9, is the space of linear polynomials in x and y. Given values zi,. . . , z,, it is well-known that there exists a unique linear spline s which interpolates the data in the sense that s(xi, yi) =zi,

i= l,...,

n.

(8)

Indeed, we can write this spline explicitly as S( X, y) = WlZil

+ W*zi2

+ w3zi3

for(X,

Y) E Tl5

(9)

where the vertices of the triangle TI are (Xii, yil), (xi2, yi2), (xi3, yi3), and where w,, w2, w3 are the barycentric coordinates of (x, y) with respect the vertices of 7’,.

L.L. Schumaker

/ Computing optimal triangulations using simulated annealing

339

By (91, the linear spline interpolating the data zi =f(xj, yj> is uniquely determined once we have chosen the triangulation. But since there are many triangulations associated with V and 0, it was proposed in [Dyn et al. ‘90a, b] to look for a best triangulation, depending on the data. Several reasonable data dependent criteria for comparing triangulations were suggested there. To illustrate the method, we choose as our local criterion the size of the angles between the normal vectors to adjacent triangular patches. Given any pair of adjoining triangles T, and T2 in the triangulation sharing an edge e, let n, and n2 be the unit normal vectors to the associated triangular patches. Then we define a(e) to be the angle between n, and n2. We compute this quantity for each of the E, interior edges in A, and let a(A)

= (ay,...,~,)

(IO)

be the corresponding vector whose components are arranged in increasing order. We now define an optimal triangulation to be one which is corresponds to the smallest vector in the lexicographical order (cf. Definition 2.1). To construct a data dependent interpolating surface using this criterion, we can apply the local swap algorithm starting with some initial triangulation. Example

6.1. Suppose f(X,

we want to interpolate

y) = (y-x2)+=

the function

(A’-X2) zt;e;;;

at 36 points on a uniform grid on the unit square 0 = [0, 11 X [O, 11 using a continuous spline based on a triangulation of R using the grid points as vertices.

(11) linear

Figure 5(b) shows the interpolating surface associated with the Delaunay triangulation of the 36 grid points depicted in Fig. 5(a). Table 4 contains the sorted vector of angles associated with the edges of this triangulation. We now apply the local swap Algorithm to adjust the triangulation in Example 6.1. Figures 5(c) and 5(d) show the resulting triangulation and the associated linear spline interpolant. In Table 4 we give the corresponding angle vector. Note that the largest angles have been significantly reduced. The many 0.0 entries correspond to pairs of triangles which lie in the flat part of the surface where y
7. Data-dependent

least squares

In this section we discuss the problem of constructing a surface to fit an unknown function, assuming that we are given measurements on the function which are subject to error (or that there are a very large number of them). In this case, interpolation as in the previous section is not appropriate, and it makes sense to perform some kind of least squares fit.

340

L.L. Schumaker

/ Computing optimal triangulations using simulated annealing

(a)

Fig. 5. Data dependent interpolation (Example 6.1).

Suppose the measurements on f are given by f(xi, yi) + ei, where P = {(xi, yi)}E, is a set errors. We now consider of points in R2, and .zi > 0, i = 1,. . . , m, are the measurement approximating f using the spaces of linear splines 9’:(A) defined in (7).

L.L. Schumaker Table 4 Angles between

patches

for Example

6.1

Angles between the patches for the Delaynay triangulation (Fig. 5(a)) 38.378 48.155 44.684 39.801 - 39.588 55.077 28.898 28.287 33.786 33.561 33.162 30.470 16.066 15.942 21.958 20.814 17.886 23.633 6.843 0.649 7.294 7.294 9.661 7.294 0.000 0.000 0.000 0.000 0.000 0.458 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 Angles 35.993 12.810 8.076 0.458 0.000 0.000 0.000 0.000

between

30.590 11.929 8.047 0.458 0.000 0.000 0.000 0.000

Angles between 26.070 11.929 5.256 0.448 0.000 0.000 0.000 0.000

the patches

for the locally optimal 26.070 11.310 5.794 0.458 0.000 0.000 0.000

the patches

22.929 10.389 1.652 0.447 0.000 0.000 0.000 0.000

23.610 10.646 5.563 0.448 0.000 0.000 0.000

for the annealed

22.113 10.270 1.480 0.020 0.000 0.000 0.000

341

/ Computing optimal triangulations using simulated annealing

20.049 10.261 1.025 0.000 0.000 0.000 0.000

triangulation 21.831 10.261 3.777 0.020 0.000 0.000 0.000

triangulation 17.257 8.954 0.649 0.000 0.000 0.000 0.000

37.734 27.284 15.642 0.649 0.000 0.000 0.000

35.445 27.284 11.310 0.459 0.000 0.000 0.000

33.992 26.003 10.302 0.458 0.000 0.000 0.000

16.563 8.844 0.649 0.000 0.000 0.000 0.000

16.521 8.707 0.459 0.000 0.000 0.000 0.000

15.793 8.321 0.459 0.000 0.000 0.000 0.000

16.521 7.844 0.459 0.000 0.000 0.000 0.000

13.031 5.962 0.458 0.000 0.000 0.000 0.000

12.505 5.794 0.458 0.000 0.000 0.000 0.000

(Fig. 5(c)) 21.558 9.622 0.649 0.020 0.000 0.000 0.000

(Fig. 5(e)) 16.563 8.321 0.648 0.000 0.000 0.000 0.000

We proceed as follows: 1. choose a vertex set V= {uJ1?i for the triangulation; 2. choose a triangulation A of the set I/; 3. choose IzJi”,,; 4. let s be the corresponding interpolating spline defined in (9). Data dependent least squares was introduced independently in [Quack & Schumaker ‘91, Rippa ‘891. Step 1 can be done in various ways. For example, if 0 is a rectangle, we can take the vertices to be the grid points of a rectangular grid. On the other hand, if we have some a priori knowledge about the behavior of the function f being fit, it may be useful to choose less regularly spaced vertices. Step 3 can also be accomplished in several different ways. For example, we may use either moving least squares (where zi is taken to be the value at (xi, yi) of some local least squares polynomial fit using nearby data), or we can perform global least squares using all the data points at once. If accurate values of f at the points of V are known, we can simply take zi =f(xj, yi> for i = 1,. . , n. Our concern here is with Step 2. We want to choose the triangulation to depend on the data since this can lead to vast improvements in the quality of the fit. In this case, the natural criterion to use in comparing triangulations (see [Quak & Schumaker ‘91, Rippa ‘891) is the following.

342

Definition

L.L. Schumaker

/ Computing optimal triangulations

using simulated annealing

7.1 (Least squares fit criterion). For each admissible triangulation

E(A) =

c[

2,

-S(

Xjt

A, let

Yj)]‘.

Define ‘%?(A,A> = E(A) - E(A). The lexicographical version of this criterion is defined as for the interpolation case in Section 6. To illustrate how annealing can be used in connection with this data dependent least squares method, we can consider the following example where the zi are taken to be the exact values of the function to be fit at the vertices I/. Example 7.2. Let I/ be the set of 36 vertices shown in Fig. 5(a) chosen on a regular grid in the unit square 0. Suppose we are given the values of the function f defined in (11) at the points of V, and also at 400 random data points P in 0. Our aim is to fit f using a linear spline which interpolates f on I/, but now in contrast to Example 6.1, we adjust the triangulation to get the best least-squares fit on the points of P.

The Delaunay triangulation of this data set is shown in Fig. 5(a), and the associated spline fit in Fig. 5(b). Applying the local swap algorithm using the above least-squares criterion, we

Fig. 6. Data-dependent

least-squares

for Example

7.1.

L.L. Schumaker

/ Computing optimal triangulations using simulated annealing

343

get the triangulation shown in Fig. 6(a) and the spline fit shown in Fig. 6(b). Clearly, there is a significant improvement in the fit. For this example, the local swap algorithm required a total of 17 swaps, and reduced the total least-squares error from 0.16 to 0.015294. (a)

Id1

Fig. 7. Delaunay

and data-dependent

triangulation

for Example 7.2.

344

L.L. Schumaker / Computing optimal triangulations using simulated annealing

To see how close the triangulation produced by local swapping is to being optimal, we ran our simulated annealing algorithm with several different choices of parameters. For example, starting with an initial temperature of 0.1 and setting nlimit = glimit = 500 and ntemp = 80, we obtained the triangulation shown in Fig. 6(c). The algorithm converged after a total of 12,618 swaps (although the best triangulation was actually encountered in about half that many swaps). The corresponding spline fit (shown in Fig. 6(d)) has a least-squares error of 0.015279. For this example, the local swap algorithm already produced a very good result. Extensive testing with simulated annealing suggests that this is generally the case for smooth functions. The approach also works even if the function being fit is discontinuous. If it is applied to a fairly fine triangulation, it can provide a very good estimate of where the discontinuity occurs as illustrated in the following example. 7.3. Suppose we want to approximate the function f(x, y> = (y -x”)“, on the unit square R using a set P of 400 random measurements in a. Let I/ be the vertices of a 9 X 9 regular grid on the unit square as shown in Fig. 7(a). We want to fit this function using an interpolating linear spline based on a triangulation with vertices at the points of I/, and using the least-squares criterion to adjust the triangulation. The Delaunay triangulation of the points I/ is shown in Fig. 7(a), an the associated linear spline interpolant in Fig. 7(b). Applying the local swap algorithm using the above least-squares criterion, we get the triangulation shown if Fig. 7(c) and the surface in Fig. 7(d). This required a total of 7 swaps, and resulted in a reduction of the least squares error from 18.00 to 12.12. It is also clear from the figure that the line of discontinuity y =x2 is better identified by the adjusted triangulation. Applying simulated annealing with initial temperature 8, nlimit = glimit = 700, and ntemp = 40, we were able to further reduce the least-squares error to 5.80 with a total of 9070 swaps. Figures 7(e) and (f) show the result. Example

8. Remarks Remark 8.1. We are not assuming in this paper that the set R is the convex hull of the set of vertices I/. In fact, we even allow J2 to have holes. Remark 8.2. If A is a triangulation of a set 0 with h holes and y1vertices. Then the number of triangles is N = 2n - nb + 2h - 2, and the number of interior edges is E, = 3n - 2n, + 3h - 3, where nb denotes the number of boundary vertices of the polygon R. Remark 8.3. It is possible to create data-dependent

interpolating splines with higher smoothness. For example, in [Quak & Schumaker ‘89, ‘90a, ‘90b] C’ cubic splines were used, although in order to keep the process local, this required splitting each triangle in a given triangulation into three subtriangles (the so-called Clough-Tocher split). For this space of fitting functions, it is natural to use thin plate energy in comparing triangulations. In particular, if s is a spline defined on a triangulation A, then its thin plate energy is given by

if1 /,rs:x

+ 2s,2, + s;,,] dx dy.

Remark 8.4. The idea of adjusting a triangulation is also useful for non-planar triangulations. For an example involving 3D reconstruction (as arises from CAT scans), see [Baszenski & Schumaker ‘921. It can also be used on spherical triangulations, see [Lawson ‘841.

L.L. Schumaker

Remark vertices ‘921 that T, by a

/ Computing optimal triangulations using simulated annealing

345

8.5. Let 7 be the set of all possible triangulations associated with a given set of I/ and a domain 0 with polygonal boundary. It has recently been shown in [Dyn et al. given any two triangulations T, and T2 in 7, it is always possible to get from T, to finite sequence of edge swaps.

References Barnhill, R.E. and Little, F.F. (1984), Three- and four-dimensional surfaces, Rocky Mt. J. Math. 14, 777102. Baszenski, G. and Schumaker, L.L. (1991), Use of simulated annealing to construct triangular facet surfaces, in: Laurent, P.-J., Le MChautt, A. and Schumaker, L.L., eds., Curces and Surfaces, Academic Press, New York, 27-32. Dyn, N., Goren, I. and Rippa, S. (1992), Transforming triangulations in polygonal domains, preprint. Dyn, N., Levin, D. and Rippa, S. (1990a) Data dependent triangulations for piecewise linear interpolation, IMA J. Numer. Anal. 10, 137-154. Dyn, N., Levitt, D. and Rippa, S. (1990b), Algorithms for the construction of data dependent triangulations, in: Mason, J.C. and Cox, M.G., eds., Algorithms for Approximation II, Chapman and Hall, London, 185-192. Lawson, CL. (1977), Software for C’ surface interpolation, in: Rice, J.R., ed., Mathematical Software III, Academic Press, New York, 161-194. Lawson, C.L. (1985), Some properties of n-dimensional triangulations, JPL Rpt. 85-42, Cal. Tech. Lawson, C.L. (1984), C’ surface interpolation for scattered data on a sphere, Rocky Mountain J. Math. 14, 177-202. Nielson, G. (19871, An example with a local minimum for the minmax ordering of triangulations, TR-87-014, Arizona State C.S. Dept.. Press, W., Flanery, B., Tenkolksy, S. and Lettlich, W. (1986), Numerical Recepies, Cambridge Univ. Press. Quak, E. and Schumaker, L.L. (19891, C’ Surface fitting using data dependent triangulations, in: Chui, C., Schumaker, L. and Ward, J., eds., Approximation Theory VI, Academic Press, New York, 545-548. Quak, E. and Schumaker, L.L. (1990a), Cubic spline interpolation using data dependent triangulations, Computer Aided Geometric Design 7, 293-301. Quak, E. and Schumaker, L.L. (1990bJ, Calculation of the energy of a piecewise polynomial surface, in: Mason, J.C. and Cox, M.G., eds., Algorithms for Approximation II, Chapman and Hall, London, 134-143. Quak, E. and Schumaker, L.L. (19911, Least squares fitting by linear splines on data dependent triangulations, in: Laurent, P.-J., Le MChautt, A. and Schumaker, L.L., eds. Curves and Surfaces, Academic Press, New York, 387-390. Renka, R.J. (1984), Algorithm 624: Triangulation and interpolation of arbitrarily distributed points in the plane, ACM TOMS 10, 440-442. Rippa, S. (1989a), Piecewise linear interpolation and approximation schemes over data dependent triangulations, Ph.D. thesis, Tel Aviv. Rippa, S. (1989b), Data compression using adaptive piecewise linear approximation, preprint. Schumaker, L.L. (1987), Triangulation methods, in: Chui, C.K., Schumaker, L.L. and Utreras, F., eds., Topics in Multirtariate Approximation, Academic Press, New York, 219-232. Schumaker, L.L. (1993), Triangulations in CAGD, Computer Graphics and Applications 13, 47-52.