## Section 5. 3D computations

(C) 1998 Joseph O'Rourke.

### Subject 5.01: How do I rotate a 3D point?

Assuming you want to rotate vectors around the origin of your coordinate system. (If you want to rotate around some other point, subtract its coordinates from the point you are rotating, do the rotation, and then add back what you subtracted.) In 3-D, you need not only an angle, but also an axis. (In higher dimensions it gets much worse, very quickly.) Actually, you need 3 independent numbers, and these come in a variety of flavors. The flavor I recommend is unit quaternions: 4 numbers that square and add up to +1. You can write these as [(x,y,z),w], with 4 real numbers, or [v,w], with v, a 3-D vector pointing along the axis. The concept of an axis is unique to 3-D. It is a line through the origin containing all the points which do not move during the rotation. So we know if we are turning forwards or back, we use a vector pointing out along the line. Suppose you want to use unit vector u as your axis, and rotate by 2t degrees. (Yes, that's twice t.) Make a quaternion [u sin t, cos t]. You can use the quaternion -- call it q -- directly on a vector v with quaternion multiplication, q v q^-1, or just convert the quaternion to a 3x3 matrix M. If the components of q are {(x,y,z),w], then you want the matrix

```M = {{1-2(yy+zz), 2(xy-wz), 2(xz+wy)},
{ 2(xy+wz),1-2(xx+zz), 2(yz-wx)},
{ 2(xz-wy), 2(yz+wx),1-2(xx+yy)}}.```

Rotations, translations, and much more are explained in all basic computer graphics texts. Quaternions are covered briefly in [Foley], and more extensively in several Graphics Gems, and the SIGGRAPH 85 proceedings.

```/* The following routine converts an angle and a unit axis vector
* to a matrix, returning the corresponding unit quaternion at no
* extra cost. It is written in such a way as to allow both fixed
* point and floating point versions to be created by appropriate
* definitions of FPOINT, ANGLE, VECTOR, QUAT, MATRIX, MUL, HALF,
* TWICE, COS, SIN, ONE, and ZERO.
* The following is an example of floating point definitions.
#define FPOINT double
#define ANGLE FPOINT
#define VECTOR QUAT
typedef struct {FPOINT x,y,z,w;} QUAT;
enum Indices {X,Y,Z,W};
typedef FPOINT MATRIX[4][4];
#define MUL(a,b) ((a)*(b))
#define HALF(a) ((a)*0.5)
#define TWICE(a) ((a)*2.0)
#define COS cos
#define SIN sin
#define ONE 1.0
#define ZERO 0.0
*/
QUAT MatrixFromAxisAngle(VECTOR axis, ANGLE theta, MATRIX m)
{
QUAT q;
ANGLE halfTheta = HALF(theta);
FPOINT cosHalfTheta = COS(halfTheta);
FPOINT sinHalfTheta = SIN(halfTheta);
FPOINT xs, ys, zs, wx, wy, wz, xx, xy, xz, yy, yz, zz;
q.x = MUL(axis.x,sinHalfTheta);
q.y = MUL(axis.y,sinHalfTheta);
q.z = MUL(axis.z,sinHalfTheta);
q.w = cosHalfTheta;
xs = TWICE(q.x);  ys = TWICE(q.y);  zs = TWICE(q.z);
wx = MUL(q.w,xs); wy = MUL(q.w,ys); wz = MUL(q.w,zs);
xx = MUL(q.x,xs); xy = MUL(q.x,ys); xz = MUL(q.x,zs);
yy = MUL(q.y,ys); yz = MUL(q.y,zs); zz = MUL(q.z,zs);
m[X][X] = ONE - (yy + zz); m[X][Y] = xy - wz; m[X][Z] = xz + wy;
m[Y][X] = xy + wz; m[Y][Y] = ONE - (xx + zz); m[Y][Z] = yz - wx;
m[Z][X] = xz - wy; m[Z][Y] = yz + wx; m[Z][Z] = ONE - (xx + yy);
/* Fill in remainder of 4x4 homogeneous transform matrix. */
m[W][X] = m[W][Y] = m[W][Z] = m[X][W] = m[Y][W] = m[Z][W] = ZERO;
m[W][W] = ONE;
return (q);
}
/* The routine just given, MatrixFromAxisAngle, performs rotation about
* an axis passing through the origin, so only a unit vector was needed
* in addition to the angle. To rotate about an axis not containing the
* origin, a point on the axis is also needed, as in the following. For
* mathematical purity, the type POINT is used, but may be defined as:
#define POINT VECTOR
*/
QUAT MatrixFromAnyAxisAngle(POINT o, VECTOR axis, ANGLE theta, MATRIX m)
{
QUAT q;
q = MatrixFromAxisAngle(axis,theta,m);
m[X][W] = o.x-(MUL(m[X][X],o.x)+MUL(m[X][Y],o.y)+MUL(m[X][Z],o.z));
m[Y][W] = o.y-(MUL(m[Y][X],o.x)+MUL(m[Y][Y],o.y)+MUL(m[Y][Z],o.z));
m[Z][W] = o.x-(MUL(m[Z][X],o.x)+MUL(m[Z][Y],o.y)+MUL(m[Z][Z],o.z));
return (q);
}
/* An axis can also be specified by a pair of points, with the direction
* along the line assumed from the ordering of the points. Although both
* the previous routines assumed the axis vector was unit length without
* checking, this routine must cope with a more delicate possibility. If
* the two points are identical, or even nearly so, the axis is unknown.
* For now the auxiliary routine which makes a unit vector ignores that.
* It needs definitions like the following for floating point.
#define SQRT sqrt
#define RECIPROCAL(a) (1.0/(a))
*/
VECTOR Normalize(VECTOR v)
{
VECTOR u;
FPOINT norm = MUL(v.x,v.x)+MUL(v.y,v.y)+MUL(v.z,v.z);
/* Better to test for (near-)zero norm before taking reciprocal. */
FPOINT scl = RECIPROCAL(SQRT(norm));
u.x = MUL(v.x,scl); u.y = MUL(v.y,scl); u.z = MUL(v.z,scl);
return (u);
}
QUAT MatrixFromPointsAngle(POINT o, POINT p, ANGLE theta, MATRIX m)
{
QUAT q;
VECTOR diff, axis;
diff.x = p.x-o.x; diff.y = p.y-o.y; diff.z = p.z-o.z;
axis = Normalize(diff);
return (MatrixFromAnyAxisAngle(o,axis,theta,m));
}```

### Subject 5.02: What is ARCBALL and where is the source?

Arcball is a general purpose 3-D rotation controller described by Ken Shoemake in the Graphics Interface '92 Proceedings. It features good behavior, easy implementation, cheap execution, and optional axis constraints. A Macintosh demo and electronic version of the original paper (Microsoft Word format) may be obtained from ftp://ftp.cis.upenn.edu/pub/graphics.

Complete source code appears in Graphics Gems IV pp. 175-192. A fairly complete sketch of the code appeared in the original article, in Graphics Interface 92 Proceedings, available from Morgan Kaufmann.

### Subject 5.03: How do I clip a polygon against a rectangle?

This is the Sutherland-Hodgman algorithm, an old favorite that should be covered in any textbook. Any of the references listed in the FAQ should have it. According to Vatti (q.v.) "This [algorithm] produces degenerate edges in certain concave / self intersecting polygons that need to be removed as a special extension to the main algorithm" (though this is not a problem if all you are doing with the results is scan converting them.)

### Subject 5.04: How do I clip a polygon against another polygon?

Klamer Schutte, klamer@ph.tn.tudelft.nl has developed and implemented some code in C++ to perform clipping of two possibly concave 2-D polygons. A description can be found at: http://www.ph.tn.tudelft.nl:/People/klamer/clippoly_entry.html To compile the source code you will need a C++ compiler with templates, such as g++. The source code is available at: ftp://ftp.ph.tn.tudelft.nl/pub/klamer/clippoly.tar.gz. See also http://www.propro.ru/leonov/clipdoc.html, which extends the above to permit holes.

Alan Murta released a polygon clipper library (in C) which uses a modified version of the Vatti algorithm: http://www.cs.man.ac.uk/aig/staff/alan/software/index.html

References:

Weiler, K. "Polygon Comparison Using a Graph Representation", SIGGRAPH 80 pg. 10-18

Vatti, Bala R. "A Generic Solution to Polygon Clipping", Communications of the ACM, July 1992, Vol 35, No. 7, pg. 57-63

### Subject 5.05: How do I find the intersection of a line and a plane?

If the plane is defined as:

`a*x + b*y + c*z + d = 0`

and the line is defined as:

```x = x1 + (x2 - x1)*t = x1 + i*t
y = y1 + (y2 - y1)*t = y1 + j*t
z = z1 + (z2 - z1)*t = z1 + k*t```

Then just substitute these into the plane equation. You end up with:

`t = - (a*x1 + b*y1 + c*z1 + d)/(a*i + b*j + c*k)`

When the denominator is zero, the line is contained in the plane if the numerator is also zero (the point at t=0 satisfies the plane equation), otherwise the line is parallel to the plane.

### Subject 5.06: How do I determine the intersection between a ray and a polygon

First find the intersection between the ray and the plane in which the polygon is situated. Then see if the point of intersection is in the polygon.

Reference:

[Glassner:RayTracing]

### Subject 5.07: How do I determine the intersection between a ray and a sphere

Given a ray defined as:

```x = x1 + (x2 - x1)*t = x1 + i*t
y = y1 + (y2 - y1)*t = y1 + j*t
z = z1 + (z2 - z1)*t = z1 + k*t```

and a sphere defined as:

`(x - l)**2 + (y - m)**2 + (z - n)**2 = r**2`

Substituting in gives the quadratic equation:

`a*t**2 + b*t + c = 0`

where:

```a = i**2 + j**2 + k**2
b = 2*i*(x1 - l) + 2*j*(y1 - m) + 2*k*(z1 - n)
c = (x1-l)**2 + (y1-m)**2 + (z1-n)**2 - r**2;```

If the discriminant of this equation is less than 0, the line does not intersect the sphere. If it is zero, the line is tangential to the sphere and if it is greater than zero it intersects at two points. Solving the equation and substituting the values of t into the ray equation will give you the points.

Reference:

[Glassner:RayTracing]

### Subject 5.08: How do I find the intersection of a ray and a bezier surface?

Joy I.K. and Bhetanabhotla M.N., "Ray tracing parametric surfaces utilizing numeric techniques and ray coherence", Computer Graphics, 16, 1986, 279-286

Joy and Bhetanabhotla is only one approach of three major method classes: tessellation, direct computation, and a hybrid of the two. Tessellation is relatively easy (you intersect the polygons making up the tessellation) and has no numerical problems, but can be inaccurate; direct is cheap on memory, but very expensive computationally, and can (and usually does) suffer from precision problems within the root solver; hybrids try to blend the two.

Reference:

[Glassner:RayTracing]

### Subject 5.09: How do I ray trace caustics?

There is a discussion in [Watt:Animation], but I don't know how good it is.

It's hard to get right.

One right (but incredibly expensive) answer is:

@InProceedings{mitchell-1992-illumination,
author = "Don P. Mitchell and Pat Hanrahan",
title = "Illumination From Curved Reflectors",
year = "1992",
month = "July",
volume = "26",
booktitle = "Computer Graphics (SIGGRAPH '92 Proceedings)",
pages = "283--291",
keywords = "caustics, interval arithmetic, ray tracing",
editor = "Edwin E. Catmull",
}

A cheat:

@Article{inakage-1986-caustics,
author = "Masa Inakage",
title = "Caustics and Specular Reflection Models for
Spherical Objects and Lenses ",
pages = "379--383",
journal = "The Visual Computer",
volume = "2",
number = "6",
year = "1986",
keywords = "ray tracing effects",
}

very specialized:

@Article{yuan-1988-gemstone,
author = "Ying Yuan and Tosiyasu L. Kunii and Naota
Inamato and Lining Sun ",
title = "Gemstone Fire: Adaptive Dispersive Ray Tracing
of Polyhedrons ",
year = "1988",
month = "November",
journal = "The Visual Computer",
volume = "4",
number = "5",
pages = "259--70",
keywords = "caustics",
}

### Subject 5.10: What is the marching cubes algorithm?

The marching cubes algorithm is used in volume rendering to construct an isosurface from a 3D field of values.

The 2D analog would be to take an image, and for each pixel, set it to black if the value is below some threshold, and set it to white if it's above the threshold. Then smooth the jagged black outlines by skinning them with lines.

The marching cubes algorithm tests the corner of each cube (or voxel) in the scalar field as being either above or below a given threshold. This yields a collection of boxes with classified corners. Since there are eight corners with one of two states, there are 256 different possible combinations for each cube. Then, for each cube, you replace the cube with a surface that meets the classification of the cube. For example, the following are some 2D examples showing the cubes and their associated surface.

The result of the marching cubes algorithm is a smooth surface that approximates the isosurface that is constant along a given threshold. This is useful for displaying a volume of oil in a geological volume, for example.

References:
"Marching Cubes: A High Resolution 3D Surface Construction Algorithm",
William E. Lorensen and Harvey E. Cline,
Computer Graphics (Proceedings of SIGGRAPH '87), Vol. 21, No. 4, pp. 163-169.

[Watt:Animation] pp. 302-305 and 313-321
[Schroeder]

### Subject 5.11: What is the status of the patent on the "marching cubes" algorithm?

United States Patent Number: 4,710,876
Date of Patent: Dec. 1, 1987
Inventors: Harvey E. Cline, William E. Lorensen
Assignee: General Electric Company
Title: "System and Method for the Display of Surface Structures Contained Within the Interior Region of a Solid Body"
Filed: Jun. 5, 1985

United States Patent Number: 4,885,688
Date of Patent: Dec. 5, 1989
Inventor: Carl R. Crawford
Assignee: General Electric Company
Title: "Minimization of Directed Points Generated in Three-Dimensional Dividing Cubes Method"
Filed: Nov. 25, 1987

### Subject 5.12: How do I do a hidden surface test (backface culling) with 3d points?

Just define all points of all polygons in clockwise order.

```c = (x3*((z1*y2)-(y1*z2))+
(y3*((x1*z2)-(z1*x2))+
(z3*((y1*x2)-(x1*y2))+```
`x1,y1,z1, x2,y2,z2, x3,y3,z3 = the first three points of the polygon.`
`If c is positive the polygon is visible. If c is negative the polygon is invisible (or the other way).`

### Subject 5.13: Where can I find algorithms for 3D collision detection?

Check out "proxima", from Purdue, which is a C++ library for 3D collision detection for arbitrary polyhedra. It's a nice system; the algorithms are sophisticated, but the code is of modest size.

ftp://ftp.cs.purdue.edu/pub/pse/PROX/prox9.1.tar.Z

For information about the I_COLLIDE 3D collision detection system http://www.cs.unc.edu/~geom/I_COLLIDE.html

"Fast Collision Detection of Moving Convex Polyhedra", Rich Rabbitz,
[Gems IV], pages 83-109, includes source in C.

SOLID: "a library for collision detection of three-dimensional objects undergoing rigid motion and deformation. SOLID is designed to be used in interactive 3D graphics applications, and is especially suited for collision detection of objects and worlds described in VRML. Written in standard C++, compiles under GNU g++ version 2.8.1 and Visual C++ 5.0."  See:    http://www.win.tue.nl/cs/tt/gino/solid/

### Subject 5.14: How do I perform basic viewing in 3d?

We describe the shape and position of objects using numbers, usually (x,y,z) coordinates. For example, the corners of a cube might be {(0,0,0), (1,0,0), (1,1,0), (0,1,0), (0,0,1), (1,0,1), (1,1,1), (0,1,1)}. A deep understanding of what we are saying with these numbers requires mathematical study, but I will try to keep this simple. At the least, we must understand that we have designated some point in space as the origin--coordinates
(0,0,0)--and marked off lines in 3 mutually perpendicular directions using equally spaced units to give us (x,y,z) values. It might be helpful to know if we are using 1 to mean 1 foot, 1 meter, or 1 parsec; the numbers alone do not tell us.

A picture on a screen is two steps removed from the 3D world it depicts. First, it is a 2D projection; and second, it is a finite resolution approximation to the infinitely precise projection. I will ignore the approximation (sampling) for now. To know what the projection looks like, we need to know where our viewpoint is, and where the plane of the projection is, both in the 3D world. Think of it as looking out a window into a scene. As artists discovered some 500 years ago, each point in the world appears to be at a point on the window. If you move your head or look out a different window, everything changes. When the mathematicians understood what the artists were doing, they invented perspective geometry.

If your viewpoint is at the origin--(0,0,0)--and the window sits parallel to the x-y plane but at z=1, projection is no more than (x,y,z) in the world appears at (x/z,y/z,1) on the plane. Distant points will have large z values, causing them to shrink in the picture. That's perspective.

The trick is to take an arbitrary viewpoint and plane, and transform the world so we have the simple viewing situation. There are two steps: move the viewpoint to the origin, then move the viewplane to the z=1 plane. If the viewpoint is at (vx,vy,vz), transform every point by the translation (x,y,z) --> (x-vx,y-vy,z-vz). This includes the viewpoint and the viewplane. Now we need to rotate so that the z axis points straight at the viewplane, then scale so it is 1 unit away.

After all this, we may find ourselves looking out upside- down. It is traditional to specify some direction in the world or viewplane as "up", and rotate so the positive y axis points that way (as nearly as possible if it's a world vector). Finally, we have acted so far as if the window was the entire plane instead of a limited portal. A final shift and scale transforms coordinates in the plane to coordinates on the screen, so that a rectangular region of interest (our "window") in the plane fills a rectangular region of the screen (our "canvas" if you like).

I have left out details of how you define and perform the rotation of the viewplane, but I'm sure someone else will be happy to supply those if there is demand. It requires knowing how to describe a plane, and how to rotate vectors to line up. Neither is difficult, but this is already using a lot of net space. One further practical difficulty is the need to clip away parts of the world behind us, so -(x,y,z) doesn't pop up at (x/z,y/z,1). (Notice the mathematics of projection alone would allow that!) But all the viewing transformations can be done using translation, rotation, scale, and a final perspective divide. If a 4x4 homogeneous matrix is used, it can represent everything needed, which saves a lot of work.

### Subject 5.15: How Do I optimize a 3D polygon mesh

References:
"Mesh Optimization" Hoppe, DeRose Duchamp, McDonald, Stuetzle, ACM COMPUTER GRAPHICS Proceedings, Annual Conference Series, 1993.

"Re-Tiling Polygonal Surfaces", Greg Turk, ACM Computer Graphics, 26, 2, July 1992

"Decimation of Triangle Meshes", Schroeder, Zarge, Lorensen, ACM Computer Graphics, 26, 2 July 1992

"Simplification of ObjectsRendered by Polygonal Approximations", Michael J. DeHaemer, Jr. and Michael J. Zyda, Computer & Graphics, Vol. 15, No. 2, 1991, Great Britain: Pergamon Press, pp. 175-184.

"Topological Refinement Procedures for Triangular Finite Element Procedures", S. A. Cannan, S. N. Muthukrishnan and R. P. Phillips, Engineering With Computers, No. 12, p. 243-255, 1996.

"Progressive Meshes", Hoppe, SIGGRAPH 96, http://research.microsoft.com/~hoppe/siggraph96pm/paper.htm

### Subject 5.16: How can I perform volume rendering?

Two principal methods can be used:

• Ray casting or front-to-back, where the volume is behind the projection plane. A ray is projected from each point in the projection plane through the volume. The ray accumulates the properties of each voxel it passes through.

• Object order or back-to-front, where the projection plane is behind the volume. Each slice of the volume is projected on the projection plane, from the farest plane to the nearest plane.

You can also use the marching-cubes algorithm, if the extraction of surfaces from the data set is easy to do (see Subject 5.10).

Here is one algorithm to do front-to-back volume rendering:

Set up a projection plane as a plane tangent to a sphere that encloses the volume. From each pixel of the projection plane, cast a ray through the volume by using a Bresenham 3D algorithm. The ray accumulates properties from each voxel intersected, stopping when the ray exits the volume. The pixel value on the projection plane determines the final color of the ray.

if C is the ray color t the ray transparency
C' the new ray color t' the new ray transparency
Cv the voxel color tv the voxel transparency
then:
C' = C + t*Cv and t' = t * tv
with initial values: C = 0.0 and t = 1.0

An alternate version: instead of C' = C + t*Cv , use :

C' = C + t*Cv*(1-tv)^p with p a float variable. Sometimes this gives the best results. When the ray has accumulated transparency, if it becomes negligible (e.g., t<0.1), the process can stop and proceed to the next ray.

References:

Bresenham 3D:
- http://www.sct.gu.edu.au/~anthony/info/graphics/bresenham.procs
- [Gems IV] p. 366
Volume rendering:
- [Watt:Animation] pp. 297-321
- IEEE Computer Graphics and application
Vol. 10, Nb. 2, March 1990 - pp. 24-53
- "Volume Visualization"
Arie Kaufman - Ed. IEEE Computer Society Press Tutorial
- "Efficient Ray Tracing of Volume Data"
Marc Levoy - ACM Transactions on Graphics, Vol. 9, Nb 3, July 1990

### Subject 5.17: Where can I get the spline description of the famous teapot etc?

See the Standard Procedural Databases software, whose official distribution site is http://www.acm.org/tog/resources/SPD/ This database contains much useful 3D code, including spline surface tessellation, for the teapot.

### Subject 5.18: How can the distance between two lines in space be computed?

The following is robust C code from Seth Teller that computes the "points of closest approach" on two 3D lines. It also classifies the lines as parallel, intersecting, or (the generic case) skew.

```    // computes pB ON line B closest to line A
// computes pA ON line A closest to line B
// return: 0 if parallel; 1 if coincident; 2 if generic (i.e., skew)
int
line_line_closest_points3d (
POINT *pA, POINT *pB,			// computed points
const POINT *a, const VECTOR *adir,		// line A, point-normal form
const POINT *b, const VECTOR *bdir )	// line B, point-normal form
{
static VECTOR Cdir, *cdir = &Cdir;
static PLANE Ac, *ac = &Ac, Bc, *bc = &Bc;

// connecting line is perpendicular to both
vcross ( cdir, adir, bdir );

// check for near-parallel lines
if ( !vnorm( cdir ) )   {
*pA = *a;	// all points are closest
*pB = *b;
return 0;	// degenerate: lines parallel
}

// form plane containing line A, parallel to cdir
plane_from_two_vectors_and_point ( ac, cdir, adir, a );

// form plane containing line B, parallel to cdir
plane_from_two_vectors_and_point ( bc, cdir, bdir, b );

// closest point on A is line A ^ bc
intersect_line_plane ( pA, a, adir, bc );

// closest point on B is line B ^ ac
intersect_line_plane ( pB, b, bdir, ac );

// distinguish intersecting, skew lines
if ( edist( pA, pB ) < 1.0E-5F )
return 1;     // coincident: lines intersect
else
return 2;     // distinct: lines skew
}```

Also Dave Eberly has code for computing distance between various geometric primitives, including MinLineLine(), at http://www.cs.unc.edu/~eberly/gr_dist.htm.

### Subject 5.19: How can I compute the volume of a polyhedron?

Assume that the surface is closed, every face is a triangle, and the vertices of each triangle are oriented ccw from the outside. Let Volume(t,p) be the signed volume of the tetrahedron formed by a point p and a triangle t. This may be computed by a 4x4 determinant, as in [O'Rourke, p.26].

Choose an arbitrary point p (e.g., the origin), and compute the sum Volume(t_i,p) for every triangle t_i of the surface. That is the volume of the object. The justification for this claim is nontrivial, but is essentially the same as the justification for the computation of the area of a polygon (Subject 2.01).

C Code available at http://cs.smith.edu/~orourke/ and http://www.acm.org/jgt/papers/Mirtich96/ .

For computing the volumes of n-d convex polytopes, there is a C implementation by Bueeler and Enge of various algorithms available at|   http://www.Mathpool.Uni-Augsburg.DE/~enge/Volumen.html.

### Subject 5.20: How can I decompose a polyhedron into convex pieces?

Usually this problem is interpreted as seeking a collection of pairwise disjoint convex polyhedra whose set union is the original polyhedron P. The following facts are known about this difficult problem:

• Not every polyhedron may be partitioned by diagonals into tetrahedra. A counterexample is due to Scho:nhardt [O'Rourke (A), p.254].
• Determining whether a polyhedron may be so partitioned is NP-hard, a result due to Seidel & Ruppert [Disc. Comput. Geom. + 7(3) 227-254 (1992).]
• Removing the restriction to diagonals, i.e., permitting so-called Steiner points, there are polyhedra of n vertices that require n^2 convex pieces in any decomposition. This was established by Chazelle [SIAM J. Comput. 13: 488-507 (1984)]. See also [O'Rourke (A), p.256]
• An algorithm of Palios & Chazelle guarantees at most O(n+r^2) pieces, where r is the number of reflex edges (i.e., grooves). [Disc. Comput. Geom. 5:505-526 (1990).]
• Barry Joe's geompack package, at ftp://ftp.cs.ualberta.ca/pub/geompack, includes a 3D convex decomposition Fortran program.
• There seems to be no other publicly available code.

### Subject 5.21: How can the circumsphere of a tetrahedron be computed?

Let a, b, c, and d be the corners of the tetrahedron, with ax, ay, and az the components of a, and likewise for b, c, and d.  Let |a| denote the Euclidean norm of a, and let a x b denote the cross product of a and b.  Then the center m and radius r of the circumsphere are given by

```      |                                                                       |
| |d-a|^2 [(b-a)x(c-a)] + |c-a|^2 [(d-a)x(b-a)] + |b-a|^2 [(c-a)x(d-a)] |
|                                                                       |
r= -------------------------------------------------------------------------
| bx-ax  by-ay  bz-az |
2 | cx-ax  cy-ay  cz-az |
| dx-ax  dy-ay  dz-az |
```

and

```          |d-a|^2 [(b-a)x(c-a)] + |c-a|^2 [(d-a)x(b-a)] + |b-a|^2 [(c-a)x(d-a)]
m= a + ---------------------------------------------------------------------
| bx-ax  by-ay  bz-az |
2 | cx-ax  cy-ay  cz-az |
| dx-ax  dy-ay  dz-az |```

#### Some notes on stability:

• Note that the expression for r is purely a function of differences between coordinates.  The advantage is that the relative error incurred in the  computation of r is also a function of the differences between the vertices, and is not influenced by the absolute coordinates of the  vertices.  In most applications, vertices are usually nearer to each other than to the origin, so this property is advantageous.

Similarly, the formula for m incurs roundoff error proportional to the differences between vertices, but not proportional to the absolute coordinates of the vertices until the final addition.
• These expressions are unstable in only one case:  if the denominator is close to zero.  This instability, which arises if the tetrahedron is nearly degenerate, is unavoidable.  Depending on your application, you   may want to use exact arithmetic to compute the value of the determinant.

See
http://www.geom.umn.edu/software/cglist/alg.html
or
http://www.cs.cmu.edu/~quake/robust.html
•

### Subject 5.22: How do I determine if two triangles in 3D intersect?

Let the two triangles be T1 and T2. If T2 lies strictly to one side of the plane containing T2, or T2 lies strictly to one side of the plane containing T1, the triangles do not intersect. Otherwise, compute the line of intersection L between the planes. Let Ik be the interval (Tk inter L), k=1,2. Either interval may be empty. T1 and T2 intersect iff I1 and I2 overlap.

This method is decribed in Tomas Mo:ller, "A fast triangle-triangle intersection test," J. Graphics Tools 2(2) 25-30 1997. C code at http://www.acm.org/jgt/papers/Moller97/.