NASA Technical Reports Server (NTRS) 19960020954: Numerical Conformal Mapping Using Cross-Ratios and Delaunay Triangulation

35 Pages·2.8 MB·English

by  NASA Technical Reports Server (NTRS)

I
Research Institute for Advanced Computer Science
NASA Ames Research Center
I Numerical conformal mapping using cross-ratios and
Delaunay triangulation
Tobin A. Driscoll
Stephen A. Vavasis
RIACS Technical Report 96.05
January 1996
Numerical conformal mapping using cross-ratios and
Delaunay triangulation
Tobin A. Driscoll
Stephen A. Vavasis
The Research Institute for Advanced Computer Science is operated by Universities Space Research
Association, The American City Building, Suite 212, Columbia, MD 21044 [phone]
Work reported herein was supported by NASA via Contract NAS 2-13721 between NASA and the
Universities Space Research Association (USRA). Work performed at the Research Institute for Advanced
Computer Science (RIACS), NASA Ames Research Center, Moffett Field, CA 94035-1000
Numerical conformal mapping using cross-ratios
and Delaunay triangulation
Tobin A. Driscoll* Stephen A. Vavasisi
January 23, 1996
Abstract
We propose a new algorithm for computing the Riemann mapping of the
unit disk to a polygon, also known as the Schwarz-Christoffel transformation.
The new algorithm, CRDT, is based on cross-ratios of the prevertices, and also
on cross-ratios of quadrilaterals in a Delaunay triangulation of the polygon.
The CRDT algorithm produces an accurate representation of the Riemann
mapping even in the presence of arbitrary long, thin regions in the polygon,
d i ea ny previous conformal mapping algorithm. We believe that CRDT can
never fail to converge to the correct Riemann mapping, but the correctness and
convergence proof depend on conjectures that we have so far not been able to
prove. We demonstrate convergence with computational experiments.
The Riemann mapping has applications to problems in two-dimensional
potential theory and to finite-difference mesh generation. We use CRDT to
produce a mapping and solve a boundary value problem on long, thin regions
for which no other algorithm can solve these problems:-- .
1 Conformal mapping
Let P be an open, simply-connected, nonempty subset of the complex plane C that is
not the entire plane. The celebrated Riemann Mapping Theorem [8] states that there
*Center for Applied Mathematics, Cornell University, Ithaca, NY 14853,
driscollQcam.cornell,edu. This work was funded in part by DOE Grant DEFG02-94ER25199
to L. N. Trefethen.
tDepartment of Computer Science, Cornell University, Ithaca, NY 14853,
vavasisQcs.cornell.edu. This author's work was supported in part by an NSF Presiden-
tial Young Investigator grant with matching funds from Xerox, AT&T, Sun and Tektronix. This
work was partially supported by Xerox Corporation, during a visit to the Xerox Palo Alto &search
Center. Copyright (c) 1995 by Xerox Corp. All rights reserved. This work was partially supported
by the Research Institute for Advanced Computer Science (RIACS) during a visit to RIACS,
sponsored by NASA under contract NAS 2-13721.
1
is an analytic function f with a nowhere-vanishing derivative such that f(D)= P,
where D denotes the open unit disk, and such that f is bijective on D. Furthermore,
let zo be an arbitrarily specified point in P and let a be an arbitrary angle in [0,27r).
Then f can be chosen so that f(0) = ~0 and argf'(0) = a. With such a specification,
f is uniquely determined. The point zo is called the conformal center of f.
This mapping f can be used to solve problems in potential theory posed on the
original domain P. The classical example is Laplace's equation Au = 0 with Dirichlet
data. Because f preserves Laplace solutions, the original problem is reduced to solving
Laplace's equation on a disk, which can be done via Fourier transforms.
A second application is finite-difference mesh generation. Because the conformal
mapping f preserves angles, any grid on the disk with orthogonal grid-line inter-
sections is mapped%y f to a 'grid 'on PI'whose grid lines also meet orthogonally.
Orthogonal grid lines simplify the task of discretizing a differential operator via finite
difference approximations.
In the case that P is a simple polygon, the Riemann mapping can be written down
almost in closed form. Let P be a finite polygon whose boundary is piecewise linear,
with no interior angles equal to 0 (no cusps). Let the vertices of P in counterclockwise
. . , . ., . .
order be denoted 21,. zn. Let the interior angles at 21,. Zn be a1,. ,an, and
let pj = aj/n - 1 for each j , SO that pj E (-1,1] for each j . (If pj=l, Zj is the tip of
a slit, and the sides adjacent to zj coincide at least partially.) Then any conformal
+
mapping f : D P has the form 181
where A, B are complex parameters (B # 0), w1,. . ., w n are points in counterclock-
wise order on the boundary of the unit disk, and the integral denotes a complex
contour integral. The points w1, . . . , are called prevertices; they map to the points
Wn .
21, . . . ,2 , under f. Formula '(1)i s dhown .as'the"Sch^.w I ar_z.- Chrilsto~e(lS -C) formula.
The Schwarz-Christoffel formula is not quite in closed form, because there is no
explicit expression for the n+4 real parameters A, B, 01,. . . , O n, where 0, = arg wi. In
practice, these parameters are determined by solving a system of nonlinear equations
derived from geometric constraints. Any particular specification of the unknown
parameters will yield some polygon (possibly covering parts of C more than once)
whose side lengths and orientation can be measured and compared to the desired
image polygon. By using ratios of sides, we can eliminate the affine scaling constants
from the system, and by arbitrarily specifying the three degrees of freedom in the
mapping, we can reduce the size of the square nonlinear system to n - 3 [18].
Actual software packages for S-C mapping like SCPACK 1171 and its cousin, the
[?I,
SC Toolbox for MATLAB solve such nonlinear systems numerically, after apply-
ing a transformation to the primitive variables that eliminates the need for explicit
enforcement of the ordering constraints on the Oj's. But two difliculties limit the
generality of polygons these packages can map:
2
The system of nonlinear equations does not have any special structure that
0
lends itself to easy solution. In fact, the system can have local minima that can
trap nonlinear solvers and prevent convergence entirely [9].
More important, SCPACK and the SC Toolbox cannot generally handle crowd-
0
ing, a phenomenon of conformal mapping that occurs whenever the domain has
any long, narrow channel [14]. The effect of such a channel is to make the pre-
vertex positions badly skewed, to a degree exponential in the aspect ratio of the
narrow region (see Fig. 5). The canonical example of crowding is a rectangle
R of side lengths a and 1, where a >> 1. If the conformal center is chosen at
the center of R and are endpoints of one of the short edges of R, then
- 21, 22
. e2- *& .=-lO(e-*Tia):. Ths; if the’ aspect-ratio-is2 0 ‘to 1, then *atl east 14
significant digits are lost when computing a difference between two 6)j’s.
One partial solution to the problem of crowding in the SC Toolbox is its pro-
vision for mapping elongated domains to rectangles. The elongation in the domain
can be matched by a similar elongation in image rectangle, alleviating the crowding
problem [9, 121. While this technique can be generalized to certain classes of mul-
tiply elongated polygons Ell], the technique becomes much more delicate, and the
fundamental domain is no longer a disk.
We propose a new algorithm that remedies both of these deficiencies. There are
two principal innovations in the new algorithm:
Our algorithm uses as primitive variables certain cross-ratios of the wj’s. Cross-
0
ratios are defined in Section 4. Because cross-ratios are invariant under frac-
tional linear transformations, we can compute many different embeddings of the
wj’s that are all conformally equivalent and hence yield the same polygon. In
particular, when evaluating f for one part of the domain, we can recompute the
wj’s so that no crowding occurs near the-points-where-f needs to be evaluated.
Thus, crowding is no longer a problem.
Our system of equations enforces the constraints that certain absolute cross-
0
ratios come out correctly in the image polygon (rather than enforcing con-
ditions about side lengths and orientations as above). These cross-ratios in
the image polygon appear to be strongly correlated with the corresponding
primitive-variable cross-ratios. The resulting nonlinear system appears to have
a monotonicity property that makes it much easier to solve than the other
formulations.
The CRDT algorithm consists of the following steps:
1. Some of the edges of the polygon P are split, i.e., new vertices are introduced
whose angles are 7r. Let the number of vertices in the split polygon be n.
2. A Delaunay triangulation of P with the new vertices is computed.
3
3. A solver is called for an implicitly specified (n- 3) x (n- 3) system of nonlinear
-
equations. The variables of these equations are the n 3 cross-ratios of the wi’s
associated with the Delaunay triangulation. Each of the n-3 equations enforces
a constraint that a crossratio in the image polygon comes out correctly.
CRDT stands for “crossratios of the Delaunay triangulation.”
The remainder of the paper is organized as follows. In Section 2 we provide the
definition and basic properties of the constrained Delaunay triangulation. In Section 3
we describe the splitting step of our algorithm. In Section 4 we define crossratios
-
and establish some of their properties. In particular, we show that n 3 cross-
ratios uniquely determine the image polygon under (1) up to similarity transform. In
Section- 5 wep resent-$hewhole ~orithm-and.,detai~-on.htoow co mpute the forward
mapping. In Section 6 we explain how CRDT circumvents crowding. In Section 7
we discuss solvers for the nonlinear system and report on experiments with various
polygonal domains. In Section 8 we describe how to use CRDT in two applications.
We know of no other algorithm that could duplicate the results of Section 8.
2 The Delaunay triangulation
Let P be a simple polygon. A triangulation of the polygon is a division of P into
nondegenerate triangles such that (a) the intersection of any two triangles is either a
common edge of the two, a common vertex, or empty; (b) the union of the triangles
is P; (6) all three vertices of every triangle are also vertices of P; and conversely (d)
every vertex of P is also a vertex of a triangle.
It is well known that for any n-vertex simple polygon P, there exists a triangu-
lation of P, and furthermore, any triangulation of P has exactly n - 2 triangles. A
triangulation edge that is not a polygon edge is called a diagonal. It is also known
that any triangulation of P hzis’edctly n -’Q-diStihct diagonals: . .
-
Let the dual of a triangulation be the following (n 2)-node graph. The graph
has one node for each triangle, and an edge between two nodes if their corresponding
triangles have a common diagonal. Thus, the dual graph has exactly n-2 vertices and
exactly n-3 edges; each edge is in correspondence with a diagonal of the triangulation.
It is well-known that the dual graph of a triangulation of a simple polygon is always
a tree. See Fig. 1 for an example of a triangulation and its dual.
Among all triangulations of P, there is a distinguished triangulation known as
the constrained Delaunay triangulation or just the Delaunay triangulation [l].T his
is a triangulation with the following (defining) property: If d is a diagonal of the
triangulation, let Q(d) be the quadrilateral associated with d, that is, the union of
the two triangles on either side of d. Then the sum of the two opposite interior angles
of Q(d)t hat are split by d is at least
7r.
It can be proved that such a triangulation always exists and, if no four points
of P are cocircular, that it is unique. Moreover, there is a simple algorithm that
4
Figure 1: The Delaunay triangulation of a 7-sided polygon. The heavy solid segments
are boundary edges, and the light solid segments are diagonals of the triangulation.
The dual graph has five nodes (circles) and four edges (dashed lines). The dual is
abstract; the geometry shown here is for convenience.
converges to the Delaunay triangulation in O(n2)s teps: First, compute an arbitrary
triangulation of P. Find a diagonal d such that the condition in the last paragraph
is violated. Then “flip” d-i.e., delete d and replace it with the other diagonal of
Q(d), and replace the two triangles formerly adjacent to d with the two triangles
thus formed. Repeat until no flips are possible. This is the algorithm used in our
computational tests, although- a’m ore eEcient O(nl og n) algorithm for constrained
Delaunay triangulation was developed by Chew [3].
3 Splitting edges
The first step of our algorithm is to split some edges of the polygon. Splitting an
edge means replacing it by several smaller edges whose union is equal to the original
edge. These new edges are joined by vertices whose angles are Notice that this
7r.
operation does not affect the Schwarz-Christoffel formula (1);a vertex whose interior
angle is has its ,8 exponent equal to 0 in that formula.
7r
The purpose of splitting is to make sure each individual quadrilateral in the De-
launay triangulation is well-conditioned. By “well-conditioned” we mean that the
prevertices of the quadrilateral are not too crowded in some valid arrangement of
the S-C prevertices. In particular, we want to avoid quadrilaterals that are long and
5
narrow with the long edges equal to polygon edges, because the prevertices of such a
quadrilateral will be crowded on the unit circle. (A long and narrow quadrilateral is
acceptable provided that the polygon is “fat” around it. See Fig. 2.)
\ /
. ,
. A L
/ \
Figure 2: On the left is an octagon with a long, narrow quadrilateral in its triangula-
tion (heavy outline). This quadrilateral would have to have its sides split because of
crowding: in a mapping with conformal center at the center of the quadrilateral, the
short edges of the quadrilateral are exponentially contracted in the preimage on the
disk. In contrast, for the polygon on the right no splitting is necessary, because the
polygon is “fat” around the quadrilateral. That is, the quadrilateral can be enclosed
in a disk that is mostly interior to the polygon.
The splitting procedure has two phases. First, for every vertex with an interior
II
angle of 7r/4 or less, we chop off the corner at v as follows. Find the largest isoceles
triangle T that can be formed by u and its two adjacent edges such that T is contained
in P, and introduce new vertices along the. two edges-that are adjacent to v at the
midpoints of the two sides of T.A fter this split, the two edges adjacent to v are said
to be protected; that is, we do not allow them to be split during the second phase.
Let P‘ denote the polygon obtained after this first part of the splitting procedure is
complete.
The second phase of the edge-splitting procedure is iterative and generates a
sequence of poIygons, each of which is a subdivision of its predecessor, starting with
P’. Let e be an unprotected edge of some polygon occurring in the iteration. Let Z(e)
be its length. Let d(e) be the smallest distance from e to any foreign vertex, where
“foreign” means a vertex other than the endpoints of e, and distance is measured
geodesically, i.e., along the shortest piecewise linear path that remains inside the
polygon. (It turns out that d(e) can be determined efficiently given a triangulation
of the polygon.) Then we say e is ill-separated if
6
-.
- .
_ .
“
Figure 3: On the left is the Delaunay triangulation of a polygon. The middle shows
the Delaunay triangulation after the sharp corners have been chopped in the first
splitting phase. On the right is the subsequent result of the second phase, in which
narrow regions are subdivided and the Delaunay triangulation is recomputed.
At each iteration we identify all ill-separated edges and split them into three equal
pieces. We repeat this until all edges are well-separated. See Fig. 3 for an example of
both phases of the splitting process. The splitting of edges and protecting of sharp
angles is reminiscent of techniques previously introduced in the finite-element mesh
generation literature; see for example [2, 4, 151. In the mesh generation literature,
the purpose of these techniques is similar to our own purpose, namely, to prevent
the occurrence of poorly shaped triangles that could arise in a triangulation of the
original (unsplit) polygon. The main difference is that finite-element mesh generation
subdivides the interior of the domain as well as its boundary and thus would avoid
both kinds of long, skinny quadrilaterals illustrated in Fig. 2.
We do not try to prove that the splits computed by this procedure are “effective”
for our algorithm, because we do not yet have an a priori characterization of well-
conditioned quadrilaterals. We do prove, however, that the second phase of the
splitting procedure described in this section always terminates after a finite number
of steps. Let r(e) be the geodesic distance from edge e to the closest foreign edge.
(A foreign edge is one that is not adjacent to e.) Let ro be the minimum of r(e) over
all unprotected edges of P‘. Notice that there can be no edge shorter than ro in P‘,
for if there were an edge eo = (v1,v2) of length shorter than TO, let el be the other
edge whose endpoint is and let be the other edge whose endpoint is Then
v1 e2 v2.
one checks that dist(e1, 5 < ro, contradicting the choice of ro.
e2) Z(e0)
We claim that the splitting procedure above never produces an edge shorter than
TO. To see this, let e = (v1,v2) be an unprotected edge of a polygon at some inter-
mediate stage of the above algorithm whose length is less than 3ro. We must argue
that we could never split e. By induction, let us assume that no edges up to now are
shorter than ro. Let eo denote the original edge of P‘ that contains e. Let v be the
7
foreign vertex closest to e, i.e., dist(v,e) = d(e). There are three cases: v lies on an
edge that was foreign to eo in P'; it lies on an edge that was adjacent to eo; or it lies
on eo itself.
In the first case, we know that dist(e0,v) 2 ro and hence dist(e,v) 2 ro. But
Z(e) < 3r0, so (2) is not satisfied and e is not split.
In the second case, let e; be the original edge that contains v, so that eo and e; are
adjacent; say their common point is v'. Since eo cannot be protected by assumption,
the interior angle at v' greater than n/4. By assumption, the distance from to v is
21'
at least ro. Therefore, some simple trigonometry shows that the distance from v to eo
is greater than To/@. Thus, &(e)> ro/@ whereas Z(e) < 3r0, so (2) is not satisfied.
In the third case is collinear with e, and its distance from e again must be at
21
least ro; so the same reGoningshows that (2) is not satisfied.
4 Cross-ratios and embeddings
Let a, b, c,d be four distinct points in the complex plane such that the order abcd forms
a quadrilateral with counterclockwise vertex order and such that ac is an interior
diagonal of the quadrilateral. We define the cross-ratio of these points to be
Note the identity p(a,b,c,d) = p(c,d,a,b ). Thus, the cross-ratio depends on the
quadrilateral and the diagonal, but not on which endpoint of the diagonal we start
at.
In general, the cross-ratio is a complex number, but there is an important special
case when it is real.
, . . - _ . -~ .
, ~
Lemma 1 Let a,b , c, d be four distinct points on a circle in counterclockwise order.
Then p(a, b, c, d) is a negative real number.
Proof. The angle of the quadrilateral ubcd at a and the angle at c are inscribed
in complementary arcs of the circle, so the sum of these angles must be n. A quick
diagram shows that (d-a)/@-a) has its arg equal to the angle at a, and (b-c)/(d-c)
has its arg equal to the angle at c. Therefore, the arg of the cross-ratio, which is the
sum of these args, is n. I
As mentioned in the introduction, the n-3 primitive real variables of the nonlinear
system arise from n - 3 cross-ratios of prevertices. The preceding lemma confirms
that these variables are indeed real. However, cross-ratios are not quite suitable as
variables, because we would have to impose side constraints that they be negative.
Instead, our unconstrained primitive variables are logarithms of the negatives of the
cross-ratios (see (3) below).
8