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));
}
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.
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.)
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://members.xoom.com/msleonov/pbcomp.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
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.
First find the intersection between the ray and the plane in which the triangle is situated. Then see if the point of intersection is inside the triangle. Details may be found in [O'Rourke (C)] pp.226-238, whose code is at http://cs.smith.edu/ orourke/ . Efficient code complete with statistical tests is described in the Mo:ller- Trumbore paper in J. Graphics Tools (C code downloadable from there): http://www.acm.org/jgt/papers/MollerTrumbore97/ See also the full paper: http://www.Graphics.Cornell.EDU/pubs/1997/MT97.html
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]
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]
See the work of Henrik Wann Jensen at http://www.gk.dtu.dk/ hwj/papers/gi98.html | [Link stale -- can anyone update? -jor/1Dec99]
@inproceedings{j-rcnls-96
, author = "Henrik Wann Jensen"
, title = "Rendering Caustics on Non-Lambertian Surfaces"
, booktitle = "Proc. Graphics Interface '96"
, pages = "116--121"
, location = "Toronto"
, year = 1996
}
Some older references:
An expensive answer: @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", }
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] James Sharman's description: http://www.exaflop.org/docs/marchcubes
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 http://www.patents.ibm.com/details?patent_number=4710876
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 http://www.patents.ibm.com/details?patent_number=4885688
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).
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, Graphics 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/
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.
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 Objects Rendered 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
[This list is in need of updating. -jor]
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.
For unshaded images (i.e., without gradient and light computations),
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
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.
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. What's listed below shows the main ideas; the full code is at http://graphics.lcs.mit.edu/ seth/geomlib/linelinecp.c
// 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 =
static PLANE Ac, *ac = , 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.magic-software.com
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 (C), 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 .
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:
o Not every polyhedron may be partitioned by diagonals into tetrahedra. A counterexample is due to Scho:nhardt [O'Rourke (A), p.254].
o 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).]
o 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]
o 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).]
o Barry Joe's geompack package, at ftp://ftp.cs.ualberta.ca/pub/geompack, includes a 3D convex decomposition Fortran program.
o There seems to be no other publicly available code.
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 |
- 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
Let the two triangles be T1 and T2. If T1 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/
This is a difficult problem. There are two main variants: (1) when the points are organized into parallel slices through the object; (2) when the points are unorganized. For (1), see this survey: D. Meyers, S. Skinner, K. Sloan. "Surfaces from Contours" ACM Trans. Graph. 11(3) Jul 1992, 228--258. Code is available at http://www-sop.inria.fr/prisme/personnel/geiger/nuages.html For (2), see this paper: H. Hoppe, T. DeRose, T. Duchamp, J. McDonald, W. Stuetzle "Surface reconstruction from unorganized points" Proc. SIGGRAPH '92, 71--78. and P. Kumar's collection of links at http://members.tripod.com/ GeomWiz/www.sites.html