Frequently Asked Questions
Section 6. Geometric Structures and Mathematics
(C) 1998 Joseph O'Rourke.
For 2-d Delaunay triangulation, try Shewchuk's triangle program. It includes options for constrained triangulation and quality mesh generation. It uses exact arithmetic.
The Delaunay triangulation is equivalent to computing the convex hull of the points lifted to a paraboloid. For n-d Delaunay triangulation try Clarkson's hull program (exact arithmetic) or Barber and Huhdanpaa's Qhull program (floating point arithmetic). The hull program also computes Voronoi volumes and alpha shapes. The Qhull program also computes 2-d Voronoi diagrams and n-d Voronoi vertices. The output of both programs may be visualized with Geomview.
There are many other codes for Delaunay triangulation and Voronoi diagrams. See Amenta's list of computational geometry software.
The Delaunay triangulation satisfies the following property: the circumcircle of each triangle is empty. The Voronoi diagram is the closest-point map, i.e., each Voronoi cell identifies the points that are closest to an input site. The Voronoi diagram is the dual of the Delaunay triangulation. Both structures are defined for general dimension. Delaunay triangulation is an important part of mesh generation.
There is a FAQ of polyhedral computation explaining how to compute n-d Delaunay triangulation and n-d Voronoi diagram using a convex hull code, and how to use the linear programming technique to determine the Voronoi cells adjacent to a given Voronoi cell efficiently for large scale or higher dimensional cases.
Avis' lrs code uses the same file formats as cdd. It uses exact arithmetic and is useful for problems with very large output size, since it does not require storing the output.
Polyhedral Computation FAQ: http://www.ifor.math.ethz.ch/ifor/staff/fukuda/fukuda.html
[O' Rourke] pp. 168-204
[Preparata & Shamos] pp. 204-225
Chew, L. P. 1987. "Constrained Delaunay
Triangulations," Proc. Third
Annual ACM Symposium on Computational Geometry.
Chew, L. P. 1989. "Constrained Delaunay
4:97-108. (UPDATED VERSION OF CHEW 1987.)
Fang, T-P. and L. A. Piegl. 1994.
"Algorithm for Constrained Delaunay
Triangulation," The Visual Computer 10:255-265. (RECOMMENDED!)
Frey, W. H. 1987. "Selective Refinement:
A New Strategy for Automatic
Node Placement in Graded Triangular Meshes," International Journal for
Numerical Methods in Engineering 24:2183-2200.
Guibas, L. and J. Stolfi. 1985.
"Primitives for the Manipulation of
General Subdivisions and the Computation of Voronoi Diagrams," ACM
Transactions on Graphics 4(2):74-123.
Karasick, M., D. Lieber, and L. R. Nackman.
1991. "Efficient Delaunay
Triangulation Using Rational Arithmetic," ACM Transactions on Graphics
Lischinski, D. 1994. "Incremental
Delaunay Triangulation," Graphics
Gems IV, P. S. Heckbert, Ed. Cambridge, MA: Academic Press Professional.
(INCLUDES C++ SOURCE CODE -- THANK YOU, DANI!)
Schuierer, S. 1989. "Delaunay
Triangulation and the Radiosity
Approach," Proc. Eurographics '89, W. Hansmann, F. R. A. Hopgood, and
W. Strasser, Eds. Elsevier Science Publishers, 345-353.
Subramanian, G., V. V. S. Raveendra, and M.
G. Kamath. 1994. "Robust
Boundary Triangulation and Delaunay Triangulation of Arbitrary
Triangular Domains," International Journal for Numerical Methods in
Watson, D. F. and G. M. Philip. 1984.
Triangulations," Computer Vision, Graphics, and Image Processing
For n-d convex hulls, try Clarkson's hull program (exact arithmetic) or Barber and Huhdanpaa's Qhull program (floating point arithmetic). Qhull handles numeric precision problems by merging facets. The output of both programs may be visualized with Geomview.
In higher dimensions, the number of facets may be much smaller than the number of lower-dimensional features. When this is the case, Fukuda's cdd program is much faster than Qhull or hull.
There are many other codes for convex hulls. See Amenta's list of computational geometry software.
Barber & http://www.geom.umn.edu/locate/qhull Huhdanpaa ftp://geom.umn.edu/pub/software/qhull.tar.Z
Clarkson: http://cm.bell-labs.com/netlib/voronoi/hull.html / ftp://cm.bell-labs.com/netlib/voronoi/hull.shar.Z
Geomview: http://www.geom.umn.edu/locate/geomview / ftp://geom.umn.edu/pub/software/geomview/
Fukuda: http://www.ifor.math.ethz.ch/staff/fukuda/cdd_home/cdd.html / ftp://ftp.ifor.math.ethz.ch/pub/fukuda/cdd/
C code for Graham's algorithm on p.80-96.
Three-dimensional convex hull discussed in Chapter 4 (p.113-67).
C code for the incremental algorithm on p.130-60.
[Preparata & Shamos] pp. 95-184
For n-d halfspace intersection, try Fukuda's cdd program or Barber and Huhdanpaa's Qhull program. Both use floating point arithmetic. Fukuda includes code for exact arithmetic. Qhull handles numeric precision problems by merging facets.
Qhull computes halfspace intersection by computing a convex hull. The intersection of halfspaces about the origin is equivalent to the convex hull of the halfspace coefficients divided by their offsets.
[Preparata & Shamos] pp. 315-320
Let p1, p2, p3 be the three vertices
(corners) of a closed triangle T (in 2D or 3D or any
dimension). Let t1, t2, t3 be real numbers that sum to 1,
with each between 0 and 1: t1 + t2 + t3 = 1, 0 <= ti <=
1. Then the point p = t1*p1 + t2*p2 + t3*p3 lies in the plane
of T and is inside T. The numbers (t1,t2,t3) are called the
barycentric coordinates of p with respect to T. They uniquely
identify p. They can be viewed as masses placed at the
vertices whose center of gravity is p. For example, let
p1=(0,0), p2=(1,0), p3=(3,2). The barycentric coordinates
(1/2,0,1/2) specify the point
p = (0,0)/2 + 0*(1,0) + (3,2)/2 = (3/2,1), the midpoint of the p1-p3 edge. If p is joined to the three vertices, T is partitioned into three triangles. Call them T1, T2, T3, with Ti not incident to pi. The areas of these triangles Ti are proportional to the barycentric coordinates ti of p.
[Coxeter, Intro. to Geometry, p.217].
Use barycentric coordinates (see item 53) . Let A, B, C be the three vertices of your triangle. Any point P inside can be expressed uniquely as P = aA + bB + cC, where a+b+c=1 and a,b,c are each >= 0. Knowing a and b permits you to calculate c=1-a-b. So if you can generate two random numbers a and b, each in [0,1], such that their sum <=1, you've got a random point in your triangle. Generate random a and b independently and uniformly in [0,1] (just divide the standard C rand() by its max value to get such a random number.) If a+b>1, replace a by 1-a, b by 1-b. Let c=1-a-b. Then aA + bB + cC is uniformly distributed in triangle ABC: the reflection step a=1-a; b=1-b gives a point (a,b) uniformly distributed in the triangle (0,0)(1,0)(0,1), which is then mapped affinely to ABC. Now you have barycentric coordinates a,b,c. Compute your point P = aA + bB + cC.
[Gems I], Turk, pp. 24-28.
"Evenly distributed" doesn't have a
single definition. There is a sense in which only the five
Platonic solids achieve regular tesselations, as they are the
only ones whose faces are regular and equal, with each vertex
incident to the same number of faces. But generally
"even distribution" focusses not so much on the
induced tesselation, as it does on the distances and
arrangement of the points/vertices. For example, eight points
can be placed
on the sphere further away from one another than is achieved by the vertices of an inscribed cube: start with an inscribed cube, twist the top face 45 degrees about the north pole, and then move the top and bottom points along the surface towards the equator a bit.
The various definitions of "evenly distributed" lead into moderately complex mathematics. This topic happens to be a FAQ on sci.math as well as on comp.graphics.algorithms; a much more extensive and rigorous discussion written by Dave Rusin can be found at: http://www.math.niu.edu/~rusin/papers/known-math/spheres/sphere.faq
A simple method of tesselating the sphere is repeated subdivision. An example starts with a unit octahedron. At each stage, every triangle in the tesselation is divided into 4 smaller triangles by creating 3 new vertices halfway along the edges. The new vertices are normalized, "pushing" them out to lie on the sphere. After N steps of subdivision, there will be 8 * 4^N triangles in the tesselation.
A simple example of tesselation by
subdivision is available at
One frequently used definition of "evenness" is a distribution that minimizes a 1/r potential energy function of all the points, so that in a global sense points are as "far away" from each other as possible. Starting from an arbitrary distribution, we can either use numerical minimization algorithms or point repulsion, in which all the points are considered to repel each other according to a 1/r^2 force law and dynamics are simulated. The algorithm is run until some convergence criterion is satisfied, and the resulting distribution is our approximation.
Jon Leech developed code to do just this. It is available at ftp://ftp.cs.unc.edu/pub/users/leech/points.tar.gz. See his README file for more information. His distribution includes sample output files for various n < 300, which may be used off-the-shelf if that is all you need.
Another method that is simpler than the above, but attains less uniformity, is based on spacing the points along a spiral that encircles the sphere. Code available at ftp://grendel.csc.smith.edu/pub/code/sphere.tar.gz (4K)
Data on various polyhedra is available at
For the icosahedron, the twelve vertices are:
(+- 1, 0, +-t), (0, +-t, +-1), and (+-t, +-1, 0)
where t = (sqrt(5)-1)/2, the golden ratio.
Reference: "Beyond the Third Dimension" by Banchoff, p.168.
There are several methods. Perhaps the easiest to understand is the "rejection method": generate random points in an origin- centered cube with opposite corners (r,r,r) and (-r,-r,-r). Reject any point p that falls outside of the sphere of radius r. Scale the vector to lie on the surface. Because the cube to sphere volume ratio is pi/6, the average number of iterations before an acceptable vector is found is 6/pi = 1.90986. This essentially doubles the effort, and makes this method slower than the "trig method." A timing comparison conducted by Ken Sloan showed that the trig method runs in about 2/3's the time of the rejection method. He found that methods based on the use of normal distributions are twice as slow as the rejection method.
The trig method. This method works only in 3-space, but it is very fast. It depends on the slightly counterintuitive fact (see proof below) that each of the three coordinates is uniformly distributed on [-1,1] (but the three are not independent, obviously). Therefore, it suffices to choose one axis (Z, say) and generate a uniformly distributed value on that axis. This constrains the chosen point to lie on a circle parallel to the X-Y plane, and the obvious trig method may be used to obtain the remaining coordinates.
(a) Choose z uniformly distributed in [-1,1]. (b) Choose t uniformly distributed on [0, 2*pi). (c) Let r = sqrt(1-z^2). (d) Let x = r * cos(t). (e) Let y = r * sin(t).
This method uses uniform deviates (faster to generate than normal deviates), and no set of coordinates is ever rejected.
Here is a proof of correctness for the fact that the z-coordinate is uniformly distributed. The proof also works for the x- and y- coordinates, but note that this works only in 3-space.
The area of a surface of revolution in 3-space is given by
A = 2 * pi * int_a^b f(x) * sqrt(1 + [f'(x}]^2) dx
Consider a zone of a sphere of radius R. Since we are integrating in the z direction, we have
f(z) = sqrt(R^2 - z^2) f'(z) = -z / sqrt(R^2-z^2)1 + [f'(z)]^2 = r^2 / (R^2-z^2)A = 2 * pi * int_a^b sqrt(R^2-z^2) * R/sqrt(R^2-z^2) dz = 2 * pi * R int_a^b dz = 2 * pi * R * (b-a) = 2 * pi * R * h,
where h = b-a is the vertical height of the zone. Notice how the integrand reduces to a constant. The density is therefore uniform.
Code available at ftp://grendel.csc.smith.edu/pub/code/sphere.tar.gz (4K)
A common convention is to write ü as "ue", so you'll also see "Pluecker". Lines in 3D can easily be given by listing the coordinates of two distinct points, for a total of six numbers. Or, they can be given as the coordinates of two distinct planes, eight numbers. What's wrong with these? Nothing; but we can do better. Plücker coordinates are, in a sense, halfway between these extremes, and can trivially generate either. Neither extreme is as efficient as Plücker coordinates for computations.
When Plücker coordinates generalize to Grassmann coordinates, as laid out beautifully in [Hodge], Chapter VII, the determinant definition is clearly the one to use. But 3D lines can use a simpler definition. Take two distinct points on a line, say
P = (Px,Py,Pz)
Q = (Qx,Qy,Qz)
Think of these as vectors from the origin, if you like. The Plücker coordinates for the line are essentially
U = P - Q
V = P x Q
Except for a scale factor, which we ignore, U and V do not depend on the specific P and Q! Cross products are perpendicular to their factors, so we always have U.V = 0. In [Stolfi] lines have orientation, so are the same only if their Plücker coordinates are related by a positive scale factor.
As determinants of homogeneous coordinates, begin with the 4x2 matrix
[ Px Qx ] row x [ Py Qy ] row y [ Pz Qz ] row z [ 1 1 ] row w
Define Plücker coordinate Gij as the determinant of rows i and j, in that order. Notice that Giw = Pi - Qi, which is Ui. Now let (i,j,k) be a cyclic permutation of (x,y,z), namely (x,y,z) or (y,z,x) or (z,x,y), and notice that Gij = Vk. Determinants are anti-symmetric in the rows, so Gij = -Gji. Thus all possible Plücker line coordinates are either zero (if i=j) or components of U or V, perhaps with a sign change. Taking the w component of a vector as 0, the determinant form will operate just as well on a point P and vector U as on two points. We can also begin with a 2x4 matrix whose rows are the coefficients of homogeneous plane equations, E.P=0, from which come dual coordinates G'ij. Now if (h,i,j,k) is an even permutation of (x,y,z,w), then Ghi = G'jk. (Just swap U and V!)
Got Plücker, want points? No problem. At least one component of U is non-zero, say Ui, which is Giw. Create homogeneous points Pj = Gjw + Gij, and Qj = Gij. (Don't expect the P and Q that we started with, and don't expect w=1.) Want plane equations? Let (i,j,k,w) be an even permutation of (x,y,z,w), so G'jk = Giw. Then create Eh = G'hk, and Fh = G'jh.
Example: Begin with P = (2,4,8) and Q = (2,3,5). Then U = (0,1,3) and V = (-4,6,-2). The direct determinant forms are Gxw=0, Gyw=1, Gzw=3, Gyz=-4, Gzx=6, Gxy=-2, and the dual forms are G'yz=0, G'zx=1, G'xy=3, G'xw=-4, G'yw=6, G'zw=-2. Take Uz = Gzw = G'xy = 3 as a suitable non-zero element. Then two planes meeting in the line are
(G'xy G'yy G'zy G'wy).P = 0
(G'xx G'xy G'xz G'xw).P = 0
That is, a point P is on the line if it satisfies both these equations:
3 Px + 0 Py + 0 Pz - 6 Pw = 0
0 Px + 3 Py - 1 Pz - 4 Pw = 0
We can also easily determine if two lines meet, or if not, how they pass. If U1 and V1 are the coordinates of line 1, U2 and V2, of line 2, we look at the sign of U1.V2 + V1.U2. If it's zero, they meet. The determinant form reveals even permutations of (x,y,z,w):
G1xw G2yz + G1yw G2zx + G1zw G2xy + G1yz G2xw + G1zx p2yw + G1xy p2zw
Two oriented lines L1 and L2 can interact in three
L1 might intersect L2, L1 might go clockwise around L2, or L1 might go counterclockwise around L2. Here are some examples:
The first and second rows are just different views of the same lines, once from the "front" and once from the "back." Here's what they might look like if you look straight down line L2 (shown here as a dot).
The Plücker coordinates of L1 and L2 give you a quick way to test which of the three it is.
cw: U1.V2 + V1.U2 < 0
ccw: U1.V2 + V1.U2 > 0
thru: U1.V2 + V1.U2 = 0
So why is this useful? Suppose you want to test if a ray intersects a triangle in 3-space. One way to do this is to represent the ray and the edges of the triangle with Plücker coordinates. The ray hits the triangle if and only if it hits one of the triangle's edges, or it's "clockwise" from all three edges, or it's "counterclockwise" from all three edges. For example...
Using Plücker coordinates, ray shooting tests like this take only a few lines of code.
Grassmann coordinates allow similar methods to be used for points, lines, planes, and so on, and in a space of any dimension. Consult [Hodge] for a thorough discussion of the theory, [Stolfi] for a little theory with a concise implementation for low dimensions, and this article for different but equally lucid description:
J. Erickson, Plücker coordinates, Ray Tracing News 10(3)
One of the best sources on Plücker coordinates is [Stolfi].
Plücker coordinates are an incredibly useful way of specifying lines and rays in 3-dimensional space by six-dimensional vectors. Suppose L is the oriented line passing through two points with homogeneous coordinates [a,b,c,d] and [w,x,y,z] -- or with standard Cartesian coordinates (b/a,c/a,d/a) and (x/w,y/w,z/w) -- in that order. Then the Plücker coordinates of L are
[a*x-b*w, a*y-c*w, b*y-c*x, a*z-d*w, b*z-d*x, c*z-d*y].
These coordinates are determinants of 2x2 submatrices (i.e., pairs of columns) of the matrix
[ a b c d ]
[ w x y z ]
listed in the following order: 12, 13, 23, 14, 24, 34. The exact order doesn't matter, as long as you match up the right pairs in the formulas: 12 with 34, 13 with 24, and 14 with 23. For the formulas to work, you have to subtract instead of adding whenever you use the 13 coordinate, "because 13 is unlucky".
Plücker coordinates are homogeneous -- multiplying all six coordinates by any real number gives you new Plücker coordinates for the same line. Also, not every set of six numbers is the Plücker coordinates of a line. Plücker coordinates [L1, L2, L3, L4, L5, L6] always satisfy the following equation:
L1*L6 - L2*L5 + L3*L4 = 0.
Two oriented lines L and M can interact in three different ways: L might intersect M, L might go clockwise around M, or L might go counterclockwise around M. Here are some examples:
The first and second rows are just different views of the same lines, once from the "front" and once from the "back." Here's what they might look like if you look straight down line M (shown here as a dot).
The Plücker coordinates of L and M give you a quick way to test which of the three it is.
cw: L1*M6 - L2*M5 + L3*M4 + L4*M3 - L5*M2 + L6*M1 < 0 ccw: L1*M6 - L2*M5 + L3*M4 + L4*M3 - L5*M2 + L6*M1 > 0 thru: L1*M6 - L2*M5 + L3*M4 + L4*M3 - L5*M2 + L6*M1 = 0
So why is this useful? Well, suppose you want to test if a ray intersects a triangle in 3-space. One way to do this is to represent the ray and the edges of the triangle with Plücker coordiantes. The ray hits the triangle if and only if it hits one of the triangle's edges, or it's "clockwise" from all three edges, or it's "counterclockwise" from all three edges. For example in this picture, the ray is oriented counterclockwise from all three edges, so it must intersect the triangle.
So if you use Plücker coordinates, ray shooting tests like this take only about ten lines of code.