Không rõ

3D computations

How do I rotate a 3D point?

Let's assume 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 3D, 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 3D vector pointing along the axis. The concept

of an axis is unique to 3D. 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));

}

What is ARCBALL and where is the source?

Arcball is a general purpose 3D 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.

The original arcball code was written for IRIS GL. A translation

into OpenGL/GLUT, and for IRIS Performer, may be found at:

http://cs.anu.edu.au/people/Hugh.Fisher/3dstuff/

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. See the selected list below.

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.)

It should be noted that the Sutherland-Hodgman algorithm

may be used to clip a polygon against any convex polygon.

Cf. also Subject 5.04.

[Foley, van Dam]: Section 3.14.1 (pp 124 - 126)

[Hearn]: Section 6-8, pp 237 - 242 (with actual C code!)

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 2D

polygons. A description can be found at:

/www.ph.tn.tudelft.nl:/People/klamer/clippoly_entry.html">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://home.attbi.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/sof...ware/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

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.

How do I determine the intersection between a ray and a triangle?

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 ©] 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/

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]

See also the "3D Object Intersection" page, described in Subject 0.05.

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]

See also the "3D Object Intersection" page, described in Subject 0.05.

How do I ray trace caustics?

See the work of Henrik Wann Jensen at

http://graphics.ucsd.edu/~henrik/

@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

}

Metropolis Light Transport handles this phenomenon well:

http://www-graphics.stanford.edu/papers/metro/

Bidirectional path tracing also handles caustics.

http://graphics.stanford.EDU/papers/veach_thesis/ (Chapter 10)

http://www.graphics.cornell.edu/~eric/thesis/

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",

}

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]

For alternatives to the (patented; cf. Subj. 5.11) marching cubes algorithm,

see

http://www.unchainedgeometry.com/jbloom/pa...pers/index.html

under "Implicit Surface Polygonization."

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

http://www.delphion.com/

Type in "4710876" (w/o commas, w/o quotes) into their search engine.

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

Access as above.

For alternative, unpatented algorithms, cf. Subj. 5.10.

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).

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/

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/

SWIFT++: a C++ library for collision detection, exact and approximate

distance computation, and contact determination of three-dimensional

polyhedral objects undergoing rigid motion.

Some preliminary results indicate that it is faster than I-COLLIDE and

V-CLIP, and more robust than I-COLLIDE.

http://www.cs.unc.edu/~geom/SWIFT++

ColDet: C++ library for 3D collison detection. Works on generic

polyhedra, and even polygon soups. Uses bounding box hierarchies

and triangle intersection tests. Released as open source under LGPL.

Tested on Windows, MacOS, and Linux. http://photoneffect.com/coldet/ .

Terdiman's lib, which might need less RAM than the above:

http://www.codercorner.com/Opcode.htm

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

objects 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).

Details of how to define and perform the rotation of the viewplane

have been left out, but see "How can I aim a camera in a specific

direction?" elsewhere in this FAQ. One simple way to designate a

plane is with the point closest to the origin, call it D. Then

a point P is on the plane if D.P = D.D; or using d = ||D|| and

N = D/d, if N.P = d. Aim the camera with N, and scale with d.

A 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!) In

fact ordinarily a clipping box, the "viewing frustum", is used

to eliminate parts of the scene outside the window left or right,

top or bottom, and too close or too far.

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.

How do I optimize/simplify 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 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/siggr...h96pm/paper.htm

Several papers by Michael Garland (quadric-based error metric):

http://graphics.cs.uiuc.edu/~garland/

Demos:

By Stan Melax: http://www.cs.ualberta.ca/~melax/polychop/

By Stefan Krause: http://www.stefan-krause.com [Gnu Open Source]

By "klaudius": http://www.klaudius.free.fr

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.

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/gra...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

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.

How can the distance between two lines in space be computed?

Let x_i be points on the respective lines and n_i unit direction

vectors along the lines. Then the distance is

| (x_1 - x_0)^T (n_1 X n_0) | / || n_1 X n_0 ||.

Often one wants the points of closest approach as well as the distance.

The following is robust C code from Seth Teller that computes the

these points 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 = &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.magic-software.com

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 .

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:

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.

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

http://www.cs.cmu.edu/~quake/robust.html

How do I determine if two triangles in 3D intersect?

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/ . See also

http://www.ce.chalmers.se/staff/tomasm/code/

http://www.magic-software.com/MgcIntersection.html

See also the "3D Object Intersection" page, described in Subject 0.05.

NASA's "Intersect" code will intersect any number of triangulated

surfaces provided that each of the surfaces is both closed and manifold.

http://www.nas.nasa.gov/~aftosmis/cart3d/s...g.html#AuxProgs

Based on "Robust and Efficient Cartesian Mesh Generation for

Component-Based Geometry" AIAA Paper 97-0196. Michael Aftosmis.

How can a 3D surface be reconstructed from a collection of points?

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.

http://www.acm.org/pubs/citations/journals...-3/p228-meyers/

Code (NUAGES) is available at

http://www-sop.inria.fr/prisme/logiciel/nuages.html.en

ftp://ftp-sop.inria.fr/prisme/NUAGES/Nuag...AGES_SRC.tar.gz

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

New code, Cocone, written with CGAL, based on recent work by

N. Amenta, S. Choi, T. K. Dey and N. Leekha:

http://www.cis.ohio-state.edu/~tamaldey/cocone.html

How can I find the smallest sphere enclosing a set of points in 3D?

Although not obvious, the smallest sphere enclosing a set of points

in any dimension can be found by Linear Programming. This was proved

by Emo Welzl in the paper, "Smallest enclosing disks (balls and

ellipsoids)" [Lecture Notes Comput. Sci., 555, Springer-Verlag, 1991,

359--370]. + Code developed by Bernd Gaertner available (GNU

General Public License) at:

http://www.inf.ethz.ch/~gaertner/miniball.html

This code is remarkably efficient: e.g., 2 seconds for 10^6 points in

3D on an Ultra-Sparc II. See also Dave Eberly's direct implementation

of Welzl's algorithm:

http://www.magic-software.com/MgcContainment3D.html

What's the big deal with quaternions?

This could mean "Why do they evoke such heated debate?" or "What are

their virtues?"

The heat of debate is hard to explain, but it's been happening for many

decades. When Gibbs first deprecated the quaternion product and split it

into a cross product and a dot product, one outraged Victorian called

the result a "hermaphrodite monster" -- and that before the Internet's

flame wars. Generally, the quaternion advocates seem to feel the

opponents are lazy or thick-headed, and that deeper understanding of

quaternions would lead to deeper appreciation. The opponents don't

appreciate that attitude, and seem to feel the advocates are snooty or

sheep, and that matrices and such are less abstract and do just fine.

(Advocates of Clifford algebra would claim that both sides are mired in

the past.) Passion aside, quaternions have appropriate uses, as do their

alternatives.

Someone new to the debate first needs to know what quaternions are, and

what they're supposed to be good for. Quaternions are a quadruple of

numbers, used to represent 3D rotations.

q = [x,y,z,w] = [(x,y,z),w] = [V,w]

The "norm" of a quaternion N(q) is conventionally the sum of the squares

of the four components. Some writers confuse this with the magnitude,

which is the square root of the norm. Another common misconception is

that only quaternions of unit norm can be used, those with the sum of

the four squares equal to 1, but that is wrong (though they are

preferred).

[U sin a,cos a] rotates by angle 2a around unit vector U

Popular non-quaternion options are 3x3 special orthogonal matrices (9

numbers with constraints), Euler angles (3 numbers), axis-angle (4

numbers), and angular velocity vectors (3 numbers). None of these

options actually _are_ rotations, which are physical; they _represent_

rotations. The distinction is important because, for example, it is

common to use an axis-angle with an angle greater than 360 degrees to

tell an animation system to spin an object more than a full turn,

something a matrix cannot say. In mathematics, the usual meaning of a

rotation would not allow the multiple spin version, which can lead to

confusing debates.

q and -q represent the same rotation

Two rotations, the physical things, can be applied one after the other.

Assuming the two rotation axes have a least one point in common, the

result will be another rotation. Some rotation representations handle

this gracefully, some don't. For quaternions and matrices, forms of

multiplication are defined such that the product gives the desired

result. For Euler angles especially there is no simple computation.

q = q2 q1 = [V2,w2][V1,w1] = [V2xV1+w2V1+w1V2,w2w1-V2.V1]

Every rotation has a reverse rotation, such that the combination of the

two leaves an object as it was. (Rotations are an algebraic "group".)

Euler angles make reversals difficult to compute. Other representations,

including quaternions, make them simple.

reverse([V,w]) = [V,-w]

q^(-1) = [-V,w]/(V.V+ww)

Two physical rotations are also more or less similar. Unit quaternions

do a particularly good job of representing similar rotations with

similar numbers.

similar(q1,q2) = |q1.q2| = | x1x2+y1y2+z1z2+w1w2 |

Points in space, the physical things, are normally represented as 3 or 4

numbers. The effect of a rotation on a collection of points can be

computed from the representation of the rotation, and here matrices seem

fastest, using three dot products. Using their own product twice,

quaternions are a bit less efficient. (They are usually converted to

matrices at the last minute.)

p2 = q p1 q^(-1)

Sequences of rotations can be interpolated, so that the object being

turned is rotated to specific poses at specific times. This motivated

Ken Shoemake's early use of quaternions in computer graphics, as

published in 1985. He used an analog of linear interpolation (sometimes

called "lerp") that he called "Slerp", and also introduced an analog of

a piecewise Bezier curve. A few years later in some course notes he

described another curve variation he called "Squad", which still seems

to be popular. Later authors have proposed many alternatives.

sin (1-t)A sin tA

Slerp(q1,q2;t) = q1 ---------- + q2 ------, cos A = q1.q2

sin A sin A

Squad(q1,a1,b2,q2;t) = Slerp(Slerp(q1,q2;t),

Slerp(a1,b2;t);

2t(1-t))

Physics simulation, aerospace control, and robotics are examples of

computations which also depend on rotation representation. Constrained

rotations like a wheel on an axle or the elbow bend of a robot typically

use specialized representations, such as an angle alone. In many general

situations, however, quaternions have proved valuable.

2 dq = W q dt, W is the angular velocity vector

User interfaces for 3D rotation also require a representation. Direct

manipulation interfaces typically use angles for jointed figures, but

for freer manipulation may use quaternions, as in Arcball or

through-the-lens camera control. As Shoemake's _full_ Graphics Gems code

for Arcball demonstrates (with the [CAPS LOCK] key), any rotation can be

graphed as an arc on a sphere. (Not to be confused with the quaternion

unit sphere in 4D.) Whether quaternions, or any other representation,

are helpful for numeric presentation and input seems a matter of taste

and circumstance.

q = U2 U1^(-1) = [U1xU2,U1.U2]

Because of their metric properties for representing rotations, unit

quaternions are most common. Advocates frequently point out that it is

far cheaper to normalize the length of a non-zero quaternion than to

bring a matrix back to rotation form. Also Shoemake's later conversion

code cheaply creates a correct rotation matrix from _any_ quaternion

(found with his Euler angle code from Graphics Gems, which does the same

for all 24 variations of that representation).

Normalize(q) = q/Sqrt(q.q)

Comparisons to Euler angles may mention "gimbal lock" (frequently

misspelled) as a disadvantage quaternions avoid. In the physical world

where gyroscopes are mounted on nested pivots, which are called gimbals,

locking is a real problem quaternions cannot help. What's usually meant

is that because the similarity of rotations organizes them somewhat like

a sphere, while similarity of vectors is quite different, an inevitable

misfit plagues Euler angles. This can show up in behavior much like

physical gimbal lock, but also in other ways. The difficulties are

topological, and aiming runs into them as well, even if quaternions are

used. Quaternion authors who propose using curves in the vector space of

quaternion logarithms often risk the misfit unawares. Frankly, you must

pick through the literature carefully, whether informal and online or

refereed and printed, because mistakes are tragically common.

To explore Graphics Gems code, see "Where is all the source?" in this

FAQ. To read more about quaternions, you have many options. Since they

were discovered in 1843 by Hamilton (the same Irish mathematician and

physicist whose name shows up in quantum mechanics), quaternions have

found many friends, as a web search will reveal. Quaternions can be

approached and applied in numerous different ways, so if you keep

looking it's likely you will find something that suits your taste and

your needs.

(Subject 0.04) [Eberly], Ch. 2.

http://www.maths.tcd.ie/pub/HistMath/Peopl...amilton/OnQuat/

Hamilton's original paper. Not easy for today's readers.

K. Shoemake. Animating Rotation with Quaternion Curves.

Proceedings of Siggraph 85.

Original animation method. Covers all the basics.

ftp://ftp.cis.upenn.edu/pub/graphics/shoemake/

Later Shoemake tutorial. More complete than most authors.

Graphics Gems I-V, various authors, various articles.

As usual, a helpful source of code and discussion.

ftp://ftp.netcom.com/pub/hb/hbaker/quaternion/

Henry Baker collects good quaternion stuff. Access spotty.

http://linux.rice.edu/~rahul/hbaker/quaternion/

Henry Baker collection with more reliable access.

http://www.eecs.wsu.edu/~hart/papers/vqr.pdf

Visualizing quaternion rotation. May help build intuition.

http://graphics.stanford.edu/courses/cs348c-95-fall

/software/quatdemo/

The GL code implementing above Hart et al. paper.

http://www.diku.dk/students/myth/quat.html

Mathematical, but not too fast and not too fancy.

http://www.cs.berkeley.edu/~laura/cs184/qu...quaternion.html

Laura Downs covers the fundamentals.

http://graphics.cs.ucdavis.edu/GraphicsNot...tes/Quaternions

/Quaternions.htm

Ken Joy covers the fundamentals.

http://www.gg.caltech.edu/STC/rr_sig97.html

High-tech interpolation method. Demanding but rewarding.

Duff, Tom: Quaternion Splines for Animating Orientation,

Proceedings of the USENIX Association Second Computer Graphics

Workshop (held in Monterey, CA 12-13 Dec. 1985), pp. 54-62.

Subdivision paper in odd place. Author last seen at Pixar.

How can I aim a camera in a specific direction?

What's needed is a method for creating a rotation that turns one unit

vector to line up with another. To aim at an object, you can subtract

the position of the camera from the position of the object to get a

vector which you then normalize. The vector you want to turn is the

camera forward vector, commonly a unit vector along the camera -z axis.

Be warned that more than one rotation can achieve aim alone. (The issue

is rather complicated, as laid out in Ken Shoemake's article on twist

reduction in Graphics Gems IV.) For example, even if the camera is

already properly aimed you could rotate it around its z axis. The most

direct rotation is given by the non-unit quaternion

q = [(b,-a,0),1-c], to aim -z axis along unit vector (a,b,c)

Normalization is advised, but it will fail for aim vector (0,0,1). In

that case just rotate 180 degrees around the y axis, using

q = [(0,1,0),0]

If the camera is level after rotation by quaternion [(x,y,z),w], the y

component of its rotated x axis will be zero, so

xy+wz = 0

If it is upright, the y component of its rotated y axis will not be

negative, so

ww-xx+yy-zz >= 0

To ensure these two desirable properties, aim with a more sophisticated

non-unit quaternion

[(bs,-at,ab),st], where s = r-c, t = r+1, and r = sqrt(aa+cc).

This can also fail to normalize, in which case normalize instead

[(0,1+c,-b),0]

Unless the aim vector is null, this will succeed. If the aim vector has

not been normalized and its magnitude is m = sqrt(aa+bb+cc), substitute

1->m. That is, use t = r+m and use m+c.

More generally, to rotate unit vector U1 directly to unit vector U2, the

non-unit quaternion will be

q = [U1xU2,1+U1.U2]

Why? If U is a unit vector, then it is normal to a plane through the

origin with equation U.P = 0. Reflection in that plane is given by

reversing the U component of P.

reflect(P,U) = P ^Ư 2(U.P)U

The quaternion product of U1 and U2 is [U1xU2,-U1.U2], so

-2 (U.P) = PU + UP

Noting UU = -1, this gives a quaternion reflection formula.

reflect(P,U) = P + (PU+UP)U

= P ^Ư P + UPU

= UPU

Reflecting with U1 then U2, by U2(U1 P U1)U2, rotates by twice the angle

between the planes, with axis perpendicular to both normals. Noting U1U2

is the conjugate of U2U1, and -q rotates like q, the rotation quaternion

q = -U2U1 = -[U2xU1,-U2.U1] = [U1xU2,U1.U2]

This q fails to aim U1 at U2 by rotating twice as much as needed, but

its square root succeeds. One square root of unit q is 1+q normalized,

geometrically the bisection of the great arc from the identity to q.

There is an inevitable singularity when U2 is the opposite of U1,

because any perpendicular axis gives an equally direct 180 degree

rotation.

[These quaternion methods were provided by Ken Shoemake.]

How do I transform normals?

In 3D, the orientation of a plane in space can be given by a

vector perpendicular to the plane, a "normal vector" or "normal"

for short. Often it is convenient to keep that vector of unit

length, or "normalized"; be careful of the different meanings

of "normality". A smooth surface has a plane tangent to each

point, and by extension a normal to that plane is called a

"surface normal". Graphics code also cheats by associating

artificial normal vectors with the vertices of polygonal models

to simulate the reflection properties of curved surfaces;

these are called "vertex normals".

The "orientation" of a plane has two meanings, both of which

usually apply. Aside from the rotational tilting and turning

meaning, there is also the sense of "side". A closed convex

surface made of polygons has two sides, an inside and an

outside, and normals can be assigned to the polygons in such

a way that they all consistently point outside. This is

often desirable for shading and culling.

When a model is defined in one coordinate system and used in

another, as is commonly done, it may be necessary to transform

normals also. If the change of point coordinates is effected

by means of a rotation plus a translation, one simply rotates

the normals as well (with no translation). Other coordinate

changes need more care.

Although it is possible to use projective transformations to

map a model into world coordinates, ordinarily they are used

only for viewing. It is usually a mistake to apply perspective

to a normal, as shading and culling are best done in world

coordinates for correct results. Also perspective may be

computed using degenerate matrices which are not invertible,

though that need not be the case. For the importance of this,

see below.

The combination of a linear transformation and a translation

is called an affine transformation, and is performed as a

matrix multiplication plus a vector addition:

p' = A(p) = Lp + t.

When the model-to-world point transform is affine, the proper

way to transform normals is with the transpose of the inverse

of L.

n' = (L^{-1})^T n

However that is not enough.

If L includes scaling effects, a unit normal in model space

will usually transform to a non-unit normal, which can cause

problems for shaders and other code. This may need correcting

after the normal is transformed.

If L includes reflection, the inside-outside orientation of

the normal is reversed. This, too, can cause problems, and

may need correcting. The determinant of L will be negative

in this case.

When a complicated distortion is used, it must be approximated

differently at each point in space by a linear transform made

up of partial derivatives, the Jacobian. The matrix for the

Jacobian replaces L in the equation for transforming normals.

Why use the transposed inverse?

Write the dot product of column vectors n and v as a matrix

product n^T v. Write vector v as a difference of points, q-p.

Let p, q, and thus v lie in the desired plane, and let n be

normal to it. Vectors at right angles have zero dot product.

n^T v = 0

The transform v' of v is

v' = (Lq+t)-(Lp+t)

= (Lq-Lp)+(t-t)

= L(q-p)

= Lv

The transform n' of n will remain normal if it satisfies

n'^T v' = n^T v

Let n' = Mn for some M. Then the requirement is

n'^T v' = (Mn)^T (Lv)

= n^T (M^T L) v

= n^T v

This holds if

M^T L = I

where I is the identity. Right multiplying by the inverse

L^{-1} and transposing both sides gives, as claimed,

M = (L^{-1})^T

When L is a rotation, L^{-1} = L^T, so M = L. When L has no

inverse it will still have an "adjoint" to substitute for

for orthogonality purposes, differing only by a scale factor.