
This article is Copyright 2002 by Joseph O'Rourke. It may be freely
redistributed in its entirety provided that this copyright notice is
not removed.

Changed items this posting (): 3.14
New items this posting (+): none

History of Changes (approx. last six months):

Changes in 1 Jul 02 posting:
3.14: Correct GIF author info, add URL. [Thanks to Greg Roelofs.]
Changes in 1 May 02 posting:
0.04: Errata for Watt & Watt book added. [Thanks to Jacob Marner.]
5.14: 3D viewing revised by Ken Shoemake.
5.23: Remove (erroneous) 3D medial axis info.
5.25: New article on quaternions by Ken Shoemake.
5.26: New article on camera aiming and quaternions by Ken Shoemake.
6.01: Add (correct) 3D medial axis info. (Thanks to Tamal Dey.)
6.09: Plucker coordinates article revised by Ken Shoemake.
Changes in 15 Apr 02 posting:
3.05: Scaling bitmaps revised by Ken Shoemake.
3.09: Morphing article written by Ken Shoemake.
6.08: Added references on random points on a sphere (Ken Shoemake).
Changes in 1 Apr 02 posting:
1.01: 2D point rotation revised by Ken Shoemake.
1.01: 2D segment intersection revised by Ken Shoemake.
5.01: 3D point rotation revised by Ken Shoemake.
0.07: Greg Ferrar's 3D rendering library no longer available.
Changes in 15 Mar 02 posting:
2.03: Reference Dan Sunday's winding number algorithm.
4.04: More detail on Beziers approximating a circle.
(Thanks to William Gibbons.)
5.22: Added NASA's "Intersect" code for intersecting triangulated
surfaces.
5.23: Updated Cocone software description.
Changes in 15 Feb 02 posting:
5.03: Noted that SutherlandHodgman can clip against any convex polygon.
(Thanks to Ben Landon.)
5.15: More links on simplifying meshes. (Thanks to Stefan Krause.)
Changes in 1 Jan 02 posting:
2.03: Fixed link to Franklin's code. (Thanks to Keith M. Briggs.)
5.13: Update to SWIFT++; add Terdiman's collision lib.
(Thanks to Pierre Terdiman.)
Changes in 1 Nov 01 posting:
6.01,02,03: Update to Qhull 3.1 release (Thanks to Brad Barber.)
Changes in 15 Sep 01 posting:
0.04: "Radiosity: A Programmer's Perspective" out of print.
0.05: CQUANT97 link no longer available; RADBIB info updated.
(Thanks to Ian Ashdown for both.)
2.01: Explained indices in more efficient formula, and restored
Sunday's version. (Thanks to Dan Sunday.)
4.04: Link for approximating a circle via a Bezier curve
(Thanks to John McDonald, Jr.)
5.10: Add in link to Jules Bloomenthal's list of papers for algorithms
that could substitute for the marching cubes algorithm.
5.11: Refer to 5.10. (Thanks to Eric Haines for both.)
Changes in 1 Sep 01 posting:
2.01: Fixed indices in efficient area formula
(Thanks to peter@Glaze.phys.dal.ca.)
2.03: Link to classic "Point in Polygon Strategies" article.
(Thanks to Eric Haines.)
5.09: Additional references for caustics (Thanks to Lars Brinkhoff.)
5.11: New links for marching cubes patent (Thanks to John Stone.)
5.17: Stale link notice.
5.23: New Cocone link for surface reconstruction.

Table of Contents

0. General Information
0.01: Charter of comp.graphics.algorithms
0.02: Are the postings to comp.graphics.algorithms archived?
0.03: How can I get this FAQ?
0.04: What are some musthave books on graphics algorithms?
0.05: Are there any online references?
0.06: Are there other graphics related FAQs?
0.07: Where is all the source?
1. 2D Computations: Points, Segments, Circles, Etc.
1.01: How do I rotate a 2D point?
1.02: How do I find the distance from a point to a line?
1.03: How do I find intersections of 2 2D line segments?
1.04: How do I generate a circle through three points?
1.05: How can the smallest circle enclosing a set of points be found?
1.06: Where can I find graph layout algorithms?
2. 2D Polygon Computations
2.01: How do I find the area of a polygon?
2.02: How can the centroid of a polygon be computed?
2.03: How do I find if a point lies within a polygon?
2.04: How do I find the intersection of two convex polygons?
2.05: How do I do a hidden surface test (backface culling) with 2D points?
2.06: How do I find a single point inside a simple polygon?
2.07: How do I find the orientation of a simple polygon?
2.08: How can I triangulate a simple polygon?
2.09: How can I find the minimum area rectangle enclosing a set of points?
3. 2D Image/Pixel Computations
3.01: How do I rotate a bitmap?
3.02: How do I display a 24 bit image in 8 bits?
3.03: How do I fill the area of an arbitrary shape?
3.04: How do I find the 'edges' in a bitmap?
3.05: How do I enlarge/sharpen/fuzz a bitmap?
3.06: How do I map a texture on to a shape?
3.07: How do I detect a 'corner' in a collection of points?
3.08: Where do I get source to display (raster font format)?
3.09: What is morphing/how is it done?
3.10: How do I quickly draw a filled triangle?
3.11: D Noise functions and turbulence in Solid texturing.
3.12: How do I generate realistic sythetic textures?
3.13: How do I convert between color models (RGB, HLS, CMYK, CIE etc)?
3.14: How is "GIF" pronounced?
4. Curve Computations
4.01: How do I generate a Bezier curve that is parallel to another Bezier?
4.02: How do I split a Bezier at a specific value for t?
4.03: How do I find a t value at a specific point on a Bezier?
4.04: How do I fit a Bezier curve to a circle?
5. 3D computations
5.01: How do I rotate a 3D point?
5.02: What is ARCBALL and where is the source?
5.03: How do I clip a polygon against a rectangle?
5.04: How do I clip a polygon against another polygon?
5.05: How do I find the intersection of a line and a plane?
5.06: How do I determine the intersection between a ray and a triangle?
5.07: How do I determine the intersection between a ray and a sphere?
5.08: How do I find the intersection of a ray and a Bezier surface?
5.09: How do I ray trace caustics?
5.10: What is the marching cubes algorithm?
5.11: What is the status of the patent on the "marching cubes" algorithm?
5.12: How do I do a hidden surface test (backface culling) with 3D points?
5.13: Where can I find algorithms for 3D collision detection?
5.14: How do I perform basic viewing in 3D?
5.15: How do I optimize/simplify a 3D polygon mesh?
5.16: How can I perform volume rendering?
5.17: Where can I get the spline description of the famous teapot etc.?
5.18: How can the distance between two lines in space be computed?
5.19: How can I compute the volume of a polyhedron?
5.20: How can I decompose a polyhedron into convex pieces?
5.21: How can the circumsphere of a tetrahedron be computed?
5.22: How do I determine if two triangles in 3D intersect?
5.23: How can a 3D surface be reconstructed from a collection of points?
5.24: How can I find the smallest sphere enclosing a set of points in 3D?
5.25: What's the big deal with quaternions?
5.26: How can I aim a camera in a specific direction?
6. Geometric Structures and Mathematics
6.01: Where can I get source for Voronoi/Delaunay triangulation?
6.02: Where do I get source for convex hull?
6.03: Where do I get source for halfspace intersection?
6.04: What are barycentric coordinates?
6.05: How do I generate a random point inside a triangle?
6.06: How do I evenly distribute N points on (tesselate) a sphere?
6.07: What are coordinates for the vertices of an icosohedron?
6.08: How do I generate random points on the surface of a sphere?
6.09: What are Plucker coordinates?
7. Contributors
7.01: How can you contribute to this FAQ?
7.02: Contributors. Who made this all possible.
Search e.g. for "Section 6" to find that section.
Search e.g. for "Subject 6.04" to find that item.

Section 0. General Information

Subject 0.01: Charter of comp.graphics.algorithms
comp.graphics.algorithms is an unmoderated newsgroup intended as a forum
for the discussion of the algorithms used in the process of generating
computer graphics. These algorithms may be recently proposed in
published journals or papers, old or previously known algorithms, or
hacks used incidental to the process of computer graphics. The scope of
these algorithms may range from an efficient way to multiply matrices,
all the way to a global illumination method incorporating raytracing,
radiosity, infinite spectrum modeling, and perhaps even mirrored balls
and lime jello.
It is hoped that this group will serve as a forum for programmers and
researchers to exchange ideas and ask questions on recent papers or
current research related to computer graphics.
comp.graphics.algorithms is not:
 for requests for gifs, or other pictures
 for requests for image translator or processing software; see
alt.binaries.pictures* FAQ
alt.binaries.pictures.utilities [now degenerated to pic postings]
alt.graphics.pixutils (image format translation)
comp.sources.misc (image viewing source code)
sci.image.processing
comp.graphics.apps.softimage
fj.comp.image
 for requests for compression software; for these try:
alt.comp.compression
comp.compression
comp.compression.research
 specifically for game development; for this try:
comp.games.development.programming.misc
comp.games.development.programming.algorithms

Subject 0.02: Are the postings to comp.graphics.algorithms archived?
Archives may be found at: http://www.faqs.org/

Subject 0.03: How can I get this FAQ?
The FAQ is posted on the 1st and 15th of every month. The easiest
way to get it is to search back in your news reader for the most
recent posting, with Subject:
comp.graphics.algorithms Frequently Asked Questions
It is posted to comp.graphics.algorithms, and crossposted to
news.answers and comp.answers.
If you can't find it on your newsreader,
you can look at a recent HTML version at the "official" FAQ archive site:
http://www.faqs.org/
The maintainer also keeps a copy of the raw ASCII, always the
latest version, accessible via http://cs.smith.edu/~orourke/FAQ.html .
Finally, you can ftp the FAQ from several sites, including:
ftp://rtfm.mit.edu/pub/faqs/graphics/algorithmsfaq
ftp://mirror.seas.gwu.edu/pub/rtfm/comp/graphics/algorithms
/comp.graphics.algorithms_Frequently_Asked_Questions
The (busy) rtfm.mit.edu site lists many alternative "mirror" sites.
Also can reach the FAQ from http://www.geom.umn.edu/software/cglist/,
which is worth visiting in its own right.

Subject 0.04: What are some musthave books on graphics algorithms?
The keywords in brackets are used to refer to the books in later
questions. They generally refer to the first author except where
it is necessary to resolve ambiguity or in the case of the Gems.
Basic computer graphics, rendering algorithms,

[Foley]
Computer Graphics: Principles and Practice (2nd Ed.),
J.D. Foley, A. van Dam, S.K. Feiner, J.F. Hughes, AddisonWesley
1990, ISBN 0201121107;
Computer Graphics: Principles and Practice, C version
J.D. Foley, A. van Dam, S.K. Feiner, J.F. Hughes, AddisonWesley
ISBN: 0201848406, 1996, 1147 pp.
[Rogers:Procedural]
Procedural Elements for Computer Graphics, Second Edition
David F. Rogers, WCB/McGraw Hill 1998, ISBN 0070535485
[Rogers:Mathematical]
Mathematical Elements for Computer Graphics 2nd Ed.,
David F. Rogers and J. Alan Adams, McGraw Hill 1990, ISBN
0070535302
[Watt:3D]
_3D Computer Graphics, 2nd Edition_,
Alan Watt, AddisonWesley 1993, ISBN 0201631865
[Glassner:RayTracing]
An Introduction to Ray Tracing,
Andrew Glassner (ed.), Academic Press 1989, ISBN 0122861604
[Gems I]
Graphics Gems,
Andrew Glassner (ed.), Academic Press 1990, ISBN 0122861655
http://www.graphicsgems.org/ for all the Gems.
[Gems II]
Graphics Gems II,
James Arvo (ed.), Academic Press 1991, ISBN 012644800
[Gems III]
Graphics Gems III,
David Kirk (ed.), Academic Press 1992, ISBN 0124096700 (with
IBM disk) or 0124096719 (with Mac disk)
See also "AP Professional Graphics CDROM Library,"
Academic Press, ISBN 012059756X, which contains Gems IIII.
[Gems IV]
Graphics Gems IV,
Paul S. Heckbert (ed.), Academic Press 1994, ISBN 0123361559
(with IBM disk) or 0123361567 (with Mac disk)
[Gems V]
Graphic Gems V,
Alan W. Paeth (ed.), Academic Press 1995, ISBN 0125434553
(with IBM disk)
[Watt:Animation]
Advanced Animation and Rendering Techniques,
Alan Watt, Mark Watt, AddisonWesley 1992, ISBN 0201544121
(Unofficial) errata: http://www.rolemaker.dk/other/AART/
[Bartels]
An Introduction to Splines for Use in Computer Graphics and
Geometric Modeling,
Richard H. Bartels, John C. Beatty, Brian A. Barsky, 1987, ISBN
0934613273
[Farin]
Curves and Surfaces for Computer Aided Geometric Design:
A Practical Guide, 4th Edition, Gerald E. Farin, Academic Press
1996. ISBN 0122490541.
[Prusinkiewicz]
The Algorithmic Beauty of Plants,
Przemyslaw W. Prusinkiewicz, Aristid Lindenmayer, SpringerVerlag,
1990, ISBN 0387972978, ISBN 3540972978
[Oliver]
Tricks of the Graphics Gurus,
Dick Oliver, et al. (2) 3.5 PC disks included, $39.95 SAMS Publishing
[Hearn]
Introduction to computer graphics,
Hearn & Baker
[Cohen]
Radiosity and Realistic Imange Sythesis,
Michael F. Cohen, John R. Wallace, Academic Press Professional
1993, ISBN 0121782700 [limited reprint 1999]
[Ashdown]
Radiosity: A Programmer's Perspective
Ian Ashdown, John Wiley & Sons 1994, ISBN 0471304441, 498 pp.
Now (Sep 2001) out of print.
[Sillion]
Radiosity & Global Illumination
Francois X. Sillion snd Claude Puech, Morgan Kaufmann 1994, ISBN
1558602771, 252 pp.
[Ebert]
Texturing and Modeling  A Procedural Approach (2nd Ed.)
David S. Ebert (ed.), F. Kenton Musgrave, Darwyn Peachey, Ken Perlin,
Steven Worley, Academic Press 1998, ISBN 0122287304, Includes CDROM.
[Schroeder]
Visualization Toolkit, 2nd Edition, The: An ObjectOriented Approach to
3D Graphics (Bk/CD) (Professional Description)
William J. Schroeder, Kenneth Martin, and Bill Lorensen,
PrenticeHall 1998, ISBN: 0139546944
See Subject 0.07 for source.
[Anderson]
PC Graphics Unleashed
Scott Anderson. SAMS Publishing, ISBN 0672305704
[Ammeraal]
Computer Graphics for Java Programmers,
Leen Ammeraal, John Wiley 1998, ISBN 0471981427.
Additional information at http://home.wxs.nl/~ammeraal/ .
[Eberly]
3D Game Engine Design: A Practical Approach to RealTime Computer Graphics.
David Eberly, Morgan Kaufmann/Academic Press, 2001.
For image processing,

[Barnsley]
Fractal Image Compression,
Michael F. Barnsley and Lyman P. Hurd, AK Peters, Ltd, 1993 ISBN
1568810008
[Jain]
Fundamentals of Image Processing,
Anil K. Jain, PrenticeHall 1989, ISBN 0133361659
[Castleman]
Digital Image Processing,
Kenneth R. Castleman, PrenticeHall 1996, ISBN(Cloth): 0132114674
(Description and errata at: "http://www.phoenix.net/~castlman")
[Pratt]
Digital Image Processing, Second Edition,
William K. Pratt, WileyInterscience 1991, ISBN 0471857661
[Gonzalez]
Digital Image Processing (3rd Ed.),
Rafael C. Gonzalez, Paul Wintz, AddisonWesley 1992, ISBN
0201508036
[Russ]
The Image Processing Handbook (3rd Ed.),
John C. Russ, CRC Press and IEEE Press 1998, ISBN 0849325323
[Russ & Russ]
The Image Processing Tool Kit v. 3.0
Chris Russ and John Russ, Reindeer Games Inc. 1999, ISBN 192880800X
[Wolberg]
Digital Image Warping,
George Wolberg, IEEE Computer Society Press Monograph 1990, ISBN
0818689447
Computational geometry

[Bowyer]
A Programmer's Geometry,
Adrian Bowyer, John Woodwark, Butterworths 1983,
ISBN 0408012420 Pbk
Out of print, but see:
Introduction to Computing with Geometry,
Adrian Bowyer and John Woodwark, 1993
ISBN 1874728038. Available in PDF:
http://www.inge.com/pubs/index.htm
[Farin & Hansford]
The Geometry Toolbox for Graphics and Modeling
by Gerald E. Farin, Dianne Hansford
A K Peters Ltd; ISBN: 1568810741
[O'Rourke (C)]
Computational Geometry in C (2nd Ed.)
Joseph O'Rourke, Cambridge University Press 1998,
ISBN 0521640105 Pbk, ISBN 0521649765 Hbk
Additional information and code at http://cs.smith.edu/~orourke/ .
[O'Rourke (A)]
Art Gallery Theorems and Algorithms
Joseph O'Rourke, Oxford University Press 1987,
ISBN 0195039653.
[Goodman & O'Rourke]
Handbook of Discrete and Computational Geometry
J. E. Goodman and J. O'Rourke, editors.
CRC Press LLC, July 1997.
ISBN:0849385245
Additional information at http://cs.smith.edu/~orourke/ .
[Samet:Application]
Applications of Spatial Data Structures: Computer Graphics,
Image Processing, and GIS,
Hanan Samet, AddisonWesley, Reading, MA, 1990.
ISBN 0201503000.
[Samet:Design & Analysis]
The Design and Analysis of Spatial Data Structures,
Hanan Samet, AddisonWesley, Reading, MA, 1990.
ISBN 0201502550.
[Mortenson]
Geometric Modeling,
Michael E. Mortenson, Wiley 1985, ISBN 0471882798
[Preparata]
Computational Geometry: An Introduction,
Franco P. Preparata, Michael Ian Shamos, SpringerVerlag 1985,
ISBN 0387961313
[Okabe]
Spatial Tessellations: Concepts and Applications of Voronoi Diagrams,
A. Okabe and B. Boots and K. Sugihara,
John Wiley, Chichester, England, 1992.
[Overmars]
Computational Geometry: Algorithms and Applications
M. de Berg and M. van Kreveld and M. Overmars and O. Schwarzkopf
SpringerVerlag, Berlin, 1997.
[Stolfi]
Oriented Projective Geometry: A Framework for Geometric Computations
Academic Press, 1991.
[Hodge]
Methods of Algebraic Geometry, Volume 1
W.V.D. Hodge and D. Pedoe, Cambridge, 1994.
ISBN 05214690074 Paperback
[Tamassia et al 199?]
Graph Drawing: Algorithms for the Visualization of Graphs
Prentice Hall; ISBN: 0133016153
Algorithms books with chapters on computational geometry

[Cormen et al.]
Introduction to Algorithms,
T. H. Cormen, C. E. Leiserson, R. L. Rivest,
The MIT Press, McGrawHill, 1990.
[Mehlhorn]
Data Structures and Algorithms,
K. Mehlhorn,
SpringerVerlag, 1984.
[Sedgewick]
R. Sedgewick,
Algorithms,
AddisonWesley, 1988.
Solid Modelling

[Mantyla]
Introduction to Solid Modeling
Martti Mantyla, Computer Science Press 1988,
ISBN 0716780151

Subject 0.05: Are there any online references?
The computational geometry community maintains its own
bibliography of publications in or closely related to that
subject. Every four months, additions and corrections are
solicited from users, after which the database is updated and
released anew. As of 7 Nov 200, it contained 13485 bibtex
entries. See Jeff Erickson's page on "Computational Geometry
Bibliographies":
http://compgeom.cs.uiuc.edu/~jeffe/compgeom/biblios.html#geombib
The bibliography can be retrieved from:
ftp://ftp.cs.usask.ca/pub/geometry/geombib.tar.gz  bibliography proper
ftp://ftp.cs.usask.ca/pub/geometry/ocgc19.ps.gz  overview published
in '93 in SIGACT News and the Internat. J. Comput. Geom. Appl.
ftp://ftp.cs.usask.ca/pub/geometry/ftphints  detailed retrieval info
Universitat Politecnica de Catalunya maintains a search engine at:
http://wwwma2.upc.es/~geomc/geombibe.html
The ACM SIGGRAPH Online Bibliography Project, by
Stephen Spencer (biblio@siggraph.org).
The database is available for anonymous FTP from the
ftp://siggraph.org/publications/bibliography directory. Please
download and examine the file READ_ME in that directory for more
specific information concerning the database.
'netlib' is a useful source for algorithms, member inquiries for
SIAM, and bibliographic searches. For information, send mail to
netlib@ornl.gov, with "send index" in the body of the mail
message.
You can also find free sources for numerical computation in C via
ftp://ftp.usc.edu/pub/Cnumanal/ . In particular, grab
numcompfreec.gz in that directory.
Check out Nick Fotis's computer graphics resources FAQ  it's
packed with pointers to all sorts of great computer graphics
stuff. This FAQ is posted biweekly to comp.graphics.
This WWW page contains links to a large number
of computer graphic related pages:
http://www.dataspace.com:84/vlib/compgraphics.html
There's a Computer Science Bibliography Server at:
http://glimpse.cs.arizona.edu:1994/bib/
with Computer Graphics, Vision and Radiosity sections
A comprehensive bibliography of color quantization papers and articles
(CQUANT97) was available at http://www.ledalite.com/library/cgis.htm.
[Link no longer available  replacement? JOR]
Modelling physically based systems for animation:
http://www.cc.gatech.edu/gvu/animation/Animation.html
The University of Manchester NURBS Library:
ftp://unix.hensa.ac.uk/pub/misc/unix/nurbs/
For an implementation of Seidel's algorithm for fast trapezoidation
and triangulation of polygons. You can get the code from:
ftp://ftp.cs.unc.edu/pub/users/narkhede/triangulation.tar.gz
Ray tracing bibliography:
http://www.acm.org/tog/resources/bib/
Quaternions and other comp sci curiosities:
ftp://ftp.netcom.com/pub/hb/hbaker/hakmem/hakmem.html
Directory of Computational Geometry Software,
collected by Nina Amenta (nina@cs.utexas.edu)
Nina Amenta is maintaining a WWW directory to computational
geometry software. The directory lives at The Geometry Center.
It has pointers to lots of convex hull and voronoi diagram programs,
triangulations, collision detection, polygon intersection, smallest
enclosing ball of a point set and other stuff.
http://www.geom.umn.edu/software/cglist/
A compact reference for realtime 3D computer graphics programming:
http://www.math.mcgill.ca/~loisel/
RADBIB is a comprehensive bibliography of radiosity and
related global illumination papers, articles, and
books. It currently includes 1,972 references.
This bibliography is available in BibTex format
(with a release date of 15 Jul 01) from:
http://www.helios32.com/ under "Resources."
The "Electronic Visualization Library" (EVlib) is a domain
secific digital library for Scientific Visualization and
Computer Graphics: http://visinfo.zib.de/
3D Object Intersection: http://www.realtimerendering.com/int/
This page presents information about a wide variety of 3D object/object
intersection tests. Presented in grid form, each axis lists ray, plane,
sphere, triangle, box, frustum, and other objects. For each combination
(e.g. sphere/box), references to articles, books, and online resources
are given.
Ray Tracing News, ed. Eric Haines: http://www.raytracingnews.com .

Subject 0.06: Are there other graphics related FAQs?
BSP Tree FAQ by Bretton Wade
http://reality.sgi.com/bspfaq/
Gamma and Color FAQs by Charles A. Poynton has
ftp://ftp.inforamp.net/pub/users/poynton/doc/colour/
http://www.inforamp.net/~poynton/
The documents are mirrored in Darmstadt, Germany at
ftp://ftp.igd.fhg.de/pub/doc/colour/

Subject 0.07: Where is all the source?
Graphics Gems source code.
http://www.graphicsgems.org
This site is now the offical distribution site for Graphics Gems code.
Master list of Computational Geometry software:
http://www.geom.umn.edu/software/cglist
Described in [Goodman & O'Rourke], Chap. 52.
Jeff Erikson's software list:
http://compgeom.cs.uiuc.edu/~jeffe/compgeom/compgeom.html
Dave Eberly's extensive collection of free geometry, graphics,
and image processing software:
http://www.magicsoftware.com/
General 'stuff'
ftp://wuarchive.wustl.edu/graphics/graphics
There are a number of interesting items in
http://graphics.lcs.mit.edu/~seth including:
 Code for 2D Voronoi, Delaunay, and Convex hull
 Mike Hoymeyer's implementation of Raimund Seidel's
O( d! n ) time linear programming algorithm for
n constraints in d dimensions
 geometric models of UC Berkeley's new computer science
building
Sources to "Computational Geometry in C", by J. O'Rourke
can be found at http://cs.smith.edu/~orourke/books/compgeom.html
or ftp://cs.smith.edu/pub/compgeom .
Greg Ferrar's C++ 3D rendering library seems no longer available
at ftp://ftp.math.ohiostate.edu/pub/users/gregt .
See http://www.flowerfire.com/ADSODA/ instead.
TAGL is a portable and extensible library that provides a subset
of OpenGL functionalities.
ftp://sunsite.unc.edu/pub/packages/programming/graphics/tagl21.tgz
Try ftp://x2ftp.oulu.fi for /pub/msdos/programming/docs/graphpro.lzh by
Michael Abrash. His XSharp package has an implementation of Xiaoulin
Wu's antialiasing algorithm (in C).
Example sources for BSP tree algorithms can be found at
http://reality.sgi.com/bspfaq/, item 24.
Mel Slater (mel@dcs.qmw.ac.uk) also made some implementations of
BSP trees and shadows for static scenes using shadow volumes
code available
http://www.dcs.qmw.ac.uk/~mel/BSP.html
ftp://ftp.dcs.qmw.ac.uk/people/mel/BSP
The Visualization Toolkit (A visualization textbook, C++ library
and Tclbased interpreter) (see [Schroeder]):
http://www.kitware.com/vtk.html
WINGED.ZIP, a C++ implementation of Baumgart's wingededge data structure:
http://www.ledalite.com/library/cgis.htm
CGAL, the Computational Geometry Algorithms Library, is written in C++
and is available at URL http://www.cs.ruu.nl/CGAL/ . It consists of
three parts. The first part is the kernel, which consists of constant
size nonmodifiable geometric primitive objects. The second part is
a collection of basic geometric datastructures and algorithms, which
are parameterized by traits classes that define the interface between
the datastructure or algorithm, and the primitives they use.
The third part consists of nongeometric support facilities.
A C++ NURBS library written by Lavoie Philippe. Version 2.1.
Results may be exported as POVRay, RIB (renderman) or VRML files.
It also offers wrappers to OpenGL:
http://yukon.genie.uottawa.ca/~lavoie/software/nurbs/
Paul Bourke has code for several problems, including isosurface
generation and Delauney triangulation, at:
http://www.swin.edu.au/astronomy/pbourke/geometry/
http://www.swin.edu.au/astronomy/pbourke/modeling/
A nearly comprehensive list of available 3D engines
(most with source code):
http://cg.cs.tuberlin.de/~ki/engines.html
See also 5.17:
Where can I get the spline description of the famous teapot etc.?
Interactive Geometry Software called "Cinderella":
http://www.cinderella.de

Section 1. 2D Computations: Points, Segments, Circles, Etc.

Subject 1.01: How do I rotate a 2D point?
In 2D, you make (X,Y) from (x,y) with a rotation by angle t so:
X = x cos t  y sin t
Y = x sin t + y cos t
As a 2x2 matrix this is very simple. If you want to rotate a
column vector v by t degrees using matrix M, use
M = [cos t sin t]
[sin t cos t]
in the product M v.
If you have a row vector, use the transpose of M (turn rows into
columns and vice versa). If you want to combine rotations, in 2D
you can just add their angles, but in higher dimensions you must
multiply their matrices.

Subject 1.02: How do I find the distance from a point to a line?
Let the point be C (Cx,Cy) and the line be AB (Ax,Ay) to (Bx,By).
Let P be the point of perpendicular projection of C on AB. The parameter
r, which indicates P's position along AB, is computed by the dot product
of AC and AB divided by the square of the length of AB:
(1) AC dot AB
r = 
AB^2
r has the following meaning:
r=0 P = A
r=1 P = B
r<0 P is on the backward extension of AB
r>1 P is on the forward extension of AB
00 C is right of AB
s=0 C is on AB
Compute s as follows:
(AyCy)(BxAx)(AxCx)(ByAy)
s = 
L^2
Then the distance from C to P = s*L.

Subject 1.03: How do I find intersections of 2 2D line segments?
This problem can be extremely easy or extremely difficult; it
depends on your application. If all you want is the intersection
point, the following should work:
Let A,B,C,D be 2space position vectors. Then the directed line
segments AB & CD are given by:
AB=A+r(BA), r in [0,1]
CD=C+s(DC), s in [0,1]
If AB & CD intersect, then
A+r(BA)=C+s(DC), or
Ax+r(BxAx)=Cx+s(DxCx)
Ay+r(ByAy)=Cy+s(DyCy) for some r,s in [0,1]
Solving the above for r and s yields
(AyCy)(DxCx)(AxCx)(DyCy)
r =  (eqn 1)
(BxAx)(DyCy)(ByAy)(DxCx)
(AyCy)(BxAx)(AxCx)(ByAy)
s =  (eqn 2)
(BxAx)(DyCy)(ByAy)(DxCx)
Let P be the position vector of the intersection point, then
P=A+r(BA) or
Px=Ax+r(BxAx)
Py=Ay+r(ByAy)
By examining the values of r & s, you can also determine some
other limiting conditions:
If 0<=r<=1 & 0<=s<=1, intersection exists
r<0 or r>1 or s<0 or s>1 line segments do not intersect
If the denominator in eqn 1 is zero, AB & CD are parallel
If the numerator in eqn 1 is also zero, AB & CD are collinear.
If they are collinear, then the segments may be projected to the x
or yaxis, and overlap of the projected intervals checked.
If the intersection point of the 2 lines are needed (lines in this
context mean infinite lines) regardless whether the two line
segments intersect, then
If r>1, P is located on extension of AB
If r<0, P is located on extension of BA
If s>1, P is located on extension of CD
If s<0, P is located on extension of DC
Also note that the denominators of eqn 1 & 2 are identical.
References:
[O'Rourke (C)] pp. 24951
[Gems III] pp. 199202 "Faster Line Segment Intersection,"

Subject 1.04: How do I generate a circle through three points?
Let the three given points be a, b, c. Use _0 and _1 to represent
x and y coordinates. The coordinates of the center p=(p_0,p_1) of
the circle determined by a, b, and c are:
A = b_0  a_0;
B = b_1  a_1;
C = c_0  a_0;
D = c_1  a_1;
E = A*(a_0 + b_0) + B*(a_1 + b_1);
F = C*(a_0 + c_0) + D*(a_1 + c_1);
G = 2.0*(A*(c_1  b_1)B*(c_0  b_0));
p_0 = (D*E  B*F) / G;
p_1 = (A*F  C*E) / G;
If G is zero then the three points are collinear and no finiteradius
circle through them exists. Otherwise, the radius of the circle is:
r^2 = (a_0  p_0)^2 + (a_1  p_1)^2
Reference:
[O' Rourke (C)] p. 201. Simplified by Jim Ward.

Subject 1.05: How can the smallest circle enclosing a set of points be found?
This circle is often called the minimum spanning circle. It can be
computed in O(n log n) time for n points. The center lies on
the furthestpoint Voronoi diagram. Computing the diagram constrains
the search for the center. Constructing the diagram can be accomplished
by a 3D convex hull algorithm; for that connection, see, e.g.,
[O'Rourke (C), p.195ff]. For direct algorithms, see:
S. Skyum, "A simple algorithm for computing the smallest enclosing circle"
Inform. Process. Lett. 37 (1991) 121125.
J. Rokne, "An Easy Bounding Circle" [Gems II] pp.1416.

Subject 1.06: Where can I find graph layout algorithms?
See the paper by Eades and Tamassia, "Algorithms for Drawing
Graphs: An Annotated Bibliography," Tech Rep CS8909, Dept. of
CS, Brown Univ, Feb. 1989.
A revised version of the annotated bibliography on graph drawing
algorithms by Giuseppe Di Battista, Peter Eades, and Roberto
Tamassia is now available from
ftp://wilma.cs.brown.edu/pub/papers/compgeo/gdbiblio.tex.gz and
ftp://wilma.cs.brown.edu/pub/papers/compgeo/gdbiblio.ps.gz
Commercial software includes the Graph Layout Toolkit from Tom Sawyer
Software http://www.tomsawyer.com and Northwoods Software's GO++
http://www.nwoods.com/go/ .

Section 2. 2D Polygon Computations

Subject 2.01: How do I find the area of a polygon?
The signed area can be computed in linear time by a simple sum.
The key formula is this:
If the coordinates of vertex v_i are x_i and y_i,
twice the signed area of a polygon is given by
2 A( P ) = sum_{i=0}^{n1} (x_i y_{i+1}  y_i x_{i+1}).
Here n is the number of vertices of the polygon, and index
arithmetic is mod n, so that x_n = x_0, etc. A rearrangement
of terms in this equation can save multiplications and operate on
coordinate differences, and so may be both faster and more
accurate:
2 A(P) = sum_{i=0}^{n1} ( x_i (y_{i+1}  y_{i1}) )
Here again modular index arithmetic is implied, with n=0 and 1=n1.
This can be avoided by extending the x[] and y[] arrays up to [n+1]
with x[n]=x[0], y[n]=y[0] and y[n+1]=y[1], and using instead
2 A(P) = sum_{i=1}^{n} ( x_i (y_{i+1}  y_{i1}) )
References: [O' Rourke (C)] Thm. 1.3.3, p. 21; [Gems II] pp. 56:
"The Area of a Simple Polygon", Jon Rokne. Dan Sunday's explanation:
http://GeometryAlgorithms.com/Archive/algorithm_0101/ where
To find the area of a planar polygon not in the xy plane, use:
2 A(P) = abs(N . (sum_{i=0}^{n1} (v_i x v_{i+1})))
where N is a unit vector normal to the plane. The `.' represents the
dot product operator, the `x' represents the cross product operator,
and abs() is the absolute value function.
[Gems II] pp. 170171:
"Area of Planar Polygons and Volume of Polyhedra", Ronald N. Goldman.

Subject 2.02: How can the centroid of a polygon be computed?
The centroid (a.k.a. the center of mass, or center of gravity)
of a polygon can be computed as the weighted sum of the centroids
of a partition of the polygon into triangles. The centroid of a
triangle is simply the average of its three vertices, i.e., it
has coordinates (x1 + x2 + x3)/3 and (y1 + y2 + y3)/3. This
suggests first triangulating the polygon, then forming a sum
of the centroids of each triangle, weighted by the area of
each triangle, the whole sum normalized by the total polygon area.
This indeed works, but there is a simpler method: the triangulation
need not be a partition, but rather can use positively and
negatively oriented triangles (with positive and negative areas),
as is used when computing the area of a polygon. This leads to
a very simple algorithm for computing the centroid, based on a
sum of triangle centroids weighted with their signed area.
The triangles can be taken to be those formed by any fixed point,
e.g., the vertex v0 of the polygon, and the two endpoints of
consecutive edges of the polygon: (v1,v2), (v2,v3), etc. The area
of a triangle with vertices a, b, c is half of this expression:
(b[X]  a[X]) * (c[Y]  a[Y]) 
(c[X]  a[X]) * (b[Y]  a[Y]);
Code available at ftp://cs.smith.edu/pub/code/centroid.c (3K).
Reference: [Gems IV] pp.36; also includes code.

Subject 2.03: How do I find if a point lies within a polygon?
The definitive reference is "Point in Polygon Strategies" by
Eric Haines [Gems IV] pp. 2446. Now also at
http://www.erichaines.com/ptinpoly.
The code in the Sedgewick book Algorithms (2nd Edition, p.354) fails
under certain circumstances. See
http://condor.informatik.UniOldenburg.DE/~stueker/graphic/index.html
for a discussion.
The essence of the raycrossing method is as follows.
Think of standing inside a field with a fence representing the polygon.
Then walk north. If you have to jump the fence you know you are now
outside the poly. If you have to cross again you know you are now
inside again; i.e., if you were inside the field to start with, the total
number of fence jumps you would make will be odd, whereas if you were
ouside the jumps will be even.
The code below is from Wm. Randolph Franklin
(see URL below) with some minor modifications for speed. It returns
1 for strictly interior points, 0 for strictly exterior, and 0 or 1
for points on the boundary. The boundary behavior is complex but
determined; in particular, for a partition of a region into polygons,
each point is "in" exactly one polygon.
(See p.243 of [O'Rourke (C)] for a discussion of boundary behavior.)
int pnpoly(int npol, float *xp, float *yp, float x, float y)
{
int i, j, c = 0;
for (i = 0, j = npol1; i < npol; j = i++) {
if ((((yp[i]<=y) && (y
(xvx,yvy,zvz). 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.

Subject 5.15: 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.
"ReTiling 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. 175184.
"Topological Refinement Procedures for Triangular Finite Element
Procedures", S. A. Cannan, S. N. Muthukrishnan and R. P. Phillips,
Engineering With Computers, No. 12, p. 243255, 1996.
"Progressive Meshes", Hoppe, SIGGRAPH 96,
http://research.microsoft.com/~hoppe/siggraph96pm/paper.htm
Several papers by Michael Garland (quadricbased error metric):
http://graphics.cs.uiuc.edu/~garland/
Demos:
By Stan Melax: http://www.cs.ualberta.ca/~melax/polychop/
By Stefan Krause: http://www.stefankrause.com [Gnu Open Source]
By "klaudius": http://www.klaudius.free.fr

Subject 5.16: How can I perform volume rendering?
Two principal methods can be used:
 Ray casting or fronttoback, 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 backtofront, 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 marchingcubes algorithm, if the extraction of
surfaces from the data set is easy to do (see Subject 5.10).
Here is one algorithm to do fronttoback 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*(1tv)^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. 297321
 IEEE Computer Graphics and application
Vol. 10, Nb. 2, March 1990  pp. 2453
 "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.
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, pointnormal form
const POINT *b, const VECTOR *bdir ) // line B, pointnormal 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 nearparallel 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.0E5F )
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.magicsoftware.com

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 (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 nd convex polytopes,
there is a C implementation by Bueeler and Enge of various
algorithms available at
http://www.Mathpool.UniAugsburg.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:
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
NPhard, a result due to Seidel & Ruppert [Disc. Comput. Geom.
7(3) 227254 (1992).]
o Removing the restriction to diagonals, i.e., permitting
socalled 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: 488507 (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:505526 (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.

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
 
 da^2 [(ba)x(ca)] + ca^2 [(da)x(ba)] + ba^2 [(ca)x(da)] 
 
r= 
 bxax byay bzaz 
2  cxax cyay czaz 
 dxax dyay dzaz 
and
da^2 [(ba)x(ca)] + ca^2 [(da)x(ba)] + ba^2 [(ca)x(da)]
m= a + 
 bxax byay bzaz 
2  cxax cyay czaz 
 dxax dyay dzaz 
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 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 triangletriangle
intersection test," J. Graphics Tools 2(2) 2530 1997. C code at
http://www.acm.org/jgt/papers/Moller97/ . See also
http://www.ce.chalmers.se/staff/tomasm/code/
http://www.magicsoftware.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/surfaceModeling.html#AuxProgs
Based on "Robust and Efficient Cartesian Mesh Generation for
ComponentBased Geometry" AIAA Paper 970196. Michael Aftosmis.

Subject 5.23: 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, 228258.
http://www.acm.org/pubs/citations/journals/tog/1992113/p228meyers/
Code (NUAGES) is available at
http://wwwsop.inria.fr/prisme/logiciel/nuages.html.en
ftp://ftpsop.inria.fr/prisme/NUAGES/Nuages/NUAGES_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, 7178.
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.ohiostate.edu/~tamaldey/cocone.html

Subject 5.24: 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, SpringerVerlag, 1991,
359370]. + 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 UltraSparc II. See also Dave Eberly's direct implementation
of Welzl's algorithm:
http://www.magicsoftware.com/MgcContainment3D.html

Subject 5.25: 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 thickheaded, 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 nonquaternion options are 3x3 special orthogonal matrices (9
numbers with constraints), Euler angles (3 numbers), axisangle (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 axisangle 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,w2w1V2.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 (1t)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(1t))
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
throughthelens 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 nonzero 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/People/Hamilton/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/quatut.ps
Later Shoemake tutorial. More complete than most authors.
Graphics Gems IV, 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/cs348c95fall
/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/quat/quaternion.html
Laura Downs covers the fundamentals.
http://graphics.cs.ucdavis.edu/GraphicsNotes/Quaternions
/Quaternions.htm
Ken Joy covers the fundamentals.
http://www.gg.caltech.edu/STC/rr_sig97.html
Hightech 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 1213 Dec. 1985), pp. 5462.
Subdivision paper in odd place. Author last seen at Pixar.

Subject 5.26: 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 nonunit quaternion
q = [(b,a,0),1c], 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
wwxx+yyzz >= 0
To ensure these two desirable properties, aim with a more sophisticated
nonunit quaternion
[(bs,at,ab),st], where s = rc, 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
nonunit 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
is
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.]

Section 6. Geometric Structures and Mathematics

Subject 6.01: Where can I get source for Voronoi/Delaunay triangulation?
For 2d Delaunay triangulation, try Shewchuk's triangle program. It
includes options for constrained triangulation and quality mesh
generation. It uses exact arithmetic.
The Delaunay triangulation is equivalent to computing the convex hull
of the points lifted to a paraboloid. For nd Delaunay triangulation
try Clarkson's hull program (exact arithmetic) or Barber and Huhdanpaa's
Qhull program (floating point arithmetic). The hull program also
computes Voronoi volumes and alpha shapes. The Qhull program also
computes 2d Voronoi diagrams and nd Voronoi vertices. The output of
both programs may be visualized with Geomview.
There are many other codes for Delaunay triangulation and Voronoi
diagrams. See Amenta's list of computational geometry software.
The Delaunay triangulation satisfies the following property: the
circumcircle of each triangle is empty. The Voronoi diagram is the
closestpoint map, i.e., each Voronoi cell identifies the points that
are closest to an input site. The Voronoi diagram is the dual of
the Delaunay triangulation. Both structures are defined for general
dimension. Delaunay triangulation is an important part of mesh
generation.
There is a FAQ of polyhedral computation explaining how to compute
nd Delaunay triangulation and nd Voronoi diagram using a convex hull
code, and how to use the linear programming technique to
determine the Voronoi cells adjacent to a given Voronoi cell
efficiently for large scale or higher dimensional cases.
Avis' lrs code uses the same file formats as cdd. It
uses exact arithmetic and is useful
for problems with very large output size, since it does not
require storing the output.
On a closely related topic, see
http://www.cis.ohiostate.edu/~tamaldey/medialaxis.htm
for computation of the 3D medial axis from the Voronoi diagram.
References:
Amenta: http://www.geom.umn.edu/software/cglist
Avis: ftp://mutt.cs.mcgill.ca/pub/C/lrs.html
Barber & http://www.geom.umn.edu/locate/qhull
Huhdanpaa ftp://geom.umn.edu/pub/software/qhull.tar.Z
ftp://ftp.geom.umn.edu/pub/software/qhull.zip
Clarkson: http://cm.belllabs.com/netlib/voronoi/hull.html
ftp://cm.belllabs.com/netlib/voronoi/hull.shar.Z
Geomview: http://www.geom.umn.edu/locate/geomview
ftp://geom.umn.edu/pub/software/geomview/
Polyhedral Computation FAQ:
http://www.ifor.math.ethz.ch/ifor/staff/fukuda/fukuda.html
Shewchuk: http://www.cs.cmu.edu/~quake/triangle.html
ftp://cm.belllabs.com/netlib/voronoi/triangle.shar.Z
[O' Rourke (C)] pp. 168204
[Preparata & Shamos] pp. 204225
Chew, L. P. 1987. "Constrained Delaunay Triangulations," Proc. Third
Annual ACM Symposium on Computational Geometry.
Chew, L. P. 1989. "Constrained Delaunay Triangulations," Algorithmica
4:97108. (UPDATED VERSION OF CHEW 1987.)
Fang, TP. and L. A. Piegl. 1994. "Algorithm for Constrained Delaunay
Triangulation," The Visual Computer 10:255265. (RECOMMENDED!)
Frey, W. H. 1987. "Selective Refinement: A New Strategy for Automatic
Node Placement in Graded Triangular Meshes," International Journal for
Numerical Methods in Engineering 24:21832200.
Guibas, L. and J. Stolfi. 1985. "Primitives for the Manipulation of
General Subdivisions and the Computation of Voronoi Diagrams," ACM
Transactions on Graphics 4(2):74123.
Karasick, M., D. Lieber, and L. R. Nackman. 1991. "Efficient Delaunay
Triangulation Using Rational Arithmetic," ACM Transactions on Graphics
10(1):7191.
Lischinski, D. 1994. "Incremental Delaunay Triangulation," Graphics
Gems IV, P. S. Heckbert, Ed. Cambridge, MA: Academic Press Professional.
(INCLUDES C++ SOURCE CODE  THANK YOU, DANI!)
[Okabe]
Schuierer, S. 1989. "Delaunay Triangulation and the Radiosity
Approach," Proc. Eurographics '89, W. Hansmann, F. R. A. Hopgood, and
W. Strasser, Eds. Elsevier Science Publishers, 345353.
Subramanian, G., V. V. S. Raveendra, and M. G. Kamath. 1994. "Robust
Boundary Triangulation and Delaunay Triangulation of Arbitrary
Triangular Domains," International Journal for Numerical Methods in
Engineering 37(10):17791789.
Watson, D. F. and G. M. Philip. 1984. "Survey: Systematic
Triangulations," Computer Vision, Graphics, and Image Processing
26:217223.

Subject 6.02: Where do I get source for convex hull?
For nd convex hulls, try Clarkson's hull program (exact arithmetic)
or Barber and Huhdanpaa's Qhull program (floating point arithmetic).
Qhull 3.1 includes triangulated output and improved
handling of difficult inputs. Qhull computes convex hulls,
Delaunay triangulations, Voronoi diagrams, and halfspace
intersections about a point.
Qhull handles numeric precision problems by merging facets. The output
of both programs may be visualized with Geomview.
In higher dimensions, the number of facets may be much smaller than
the number of lowerdimensional features. When this is the case,
Fukuda's cdd program is much faster than Qhull or hull.
There are many other codes for convex hulls. See Amenta's list of
computational geometry software.
References:
Amenta: http://www.geom.umn.edu/software/cglist/
Barber & http://www.geom.umn.edu/locate/qhull
Huhdanpaa
Clarkson: http://cm.belllabs.com/netlib/voronoi/hull.html
ftp://cm.belllabs.com/netlib/voronoi/hull.shar.Z
Geomview: http://www.geom.umn.edu/locate/geomview
ftp://geom.umn.edu/pub/software/geomview/
Fukuda: http://www.ifor.math.ethz.ch/staff/fukuda/cdd_home/cdd.html
ftp://ftp.ifor.math.ethz.ch/pub/fukuda/cdd/
[O' Rourke (C)] pp. 70167
C code for Graham's algorithm on p.8096.
Threedimensional convex hull discussed in Chapter 4 (p.11367).
C code for the incremental algorithm on p.13060.
[Preparata & Shamos] pp. 95184

Subject 6.03: Where do I get source for halfspace intersection?
For nd halfspace intersection, try Fukuda's cdd program or Barber
and Huhdanpaa's Qhull program. Both use floating point arithmetic.
Fukuda includes code for exact arithmetic. Qhull handles numeric
precision problems by merging facets.
Qhull computes halfspace intersection by computing a convex hull.
The intersection of halfspaces about the origin is equivalent to the
convex hull of the halfspace coefficients divided by their offsets.
See Subject 6.02 for more information.
References:
Barber & http://www.geom.umn.edu/locate/qhull
Huhdanpaa
Fukuda: ftp://ifor13.ethz.ch/pub/fukuda/cdd/
[Preparata & Shamos] pp. 315320

Subject 6.04: What are barycentric coordinates?
Let p1, p2, p3 be the three vertices (corners) of a closed
triangle T (in 2D or 3D or any dimension). Let t1, t2, t3 be real
numbers that sum to 1, with each between 0 and 1: t1 + t2 + t3 = 1,
0 <= ti <= 1. Then the point p = t1*p1 + t2*p2 + t3*p3 lies in
the plane of T and is inside T. The numbers (t1,t2,t3) are called the
barycentric coordinates of p with respect to T. They uniquely identify
p. They can be viewed as masses placed at the vertices whose
center of gravity is p.
For example, let p1=(0,0), p2=(1,0), p3=(3,2). The
barycentric coordinates (1/2,0,1/2) specify the point
p = (0,0)/2 + 0*(1,0) + (3,2)/2 = (3/2,1), the midpoint of the p1p3
edge.
If p is joined to the three vertices, T is partitioned
into three triangles. Call them T1, T2, T3, with Ti not incident
to pi. The areas of these triangles Ti are proportional to the
barycentric coordinates ti of p.
Reference:
[Coxeter, Intro. to Geometry, p.217].

Subject 6.05: How can I generate a random point inside a triangle?
Use barycentric coordinates (see item 53) . Let A, B, C be the
three vertices of your triangle. Any point P inside can be expressed
uniquely as P = aA + bB + cC, where a+b+c=1 and a,b,c are each >= 0.
Knowing a and b permits you to calculate c=1ab. So if you can
generate two random numbers a and b, each in [0,1], such that
their sum <=1, you've got a random point in your triangle.
Generate random a and b independently and uniformly in [0,1]
(just divide the standard C rand() by its max value to get such a
random number.) If a+b>1, replace a by 1a, b by 1b. Let c=1ab.
Then aA + bB + cC is uniformly distributed in triangle ABC: the reflection
step a=1a; b=1b gives a point (a,b) uniformly distributed in the triangle
(0,0)(1,0)(0,1), which is then mapped affinely to ABC.
Now you have barycentric coordinates a,b,c. Compute your point
P = aA + bB + cC.
Reference: [Gems I], Turk, pp. 2428, contains a similar but different
method which requires a square root.

Subject 6.06: How do I evenly distribute N points on (tesselate) a sphere?
"Evenly distributed" doesn't have a single definition. There is a
sense in which only the five Platonic solids achieve regular
tesselations, as they are the only ones whose faces are regular
and equal, with each vertex incident to the same number of faces.
But generally "even distribution" focusses not so much on the
induced tesselation, as it does on the distances and arrangement
of the points/vertices. For example, eight points can be placed
on the sphere further away from one another than is achieved by
the vertices of an inscribed cube: start with an inscribed cube,
twist the top face 45 degrees about the north pole, and then
move the top and bottom points along the surface towards the equator
a bit.
The various definitions of "evenly distributed" lead into moderately
complex mathematics. This topic happens to be a FAQ on sci.math as well
as on comp.graphics.algorithms; a much more extensive and rigorous
discussion written by Dave Rusin can be found at:
http://www.math.niu.edu/~rusin/knownmath/95/sphere.faq
A simple method of tesselating the sphere is repeated subdivision. An
example starts with a unit octahedron. At each stage, every triangle in
the tesselation is divided into 4 smaller triangles by creating 3 new
vertices halfway along the edges. The new vertices are normalized,
"pushing" them out to lie on the sphere. After N steps of subdivision,
there will be 8 * 4^N triangles in the tesselation.
A simple example of tesselation by subdivision is available at
ftp://ftp.cs.unc.edu/pub/users/jon/sphere.c
One frequently used definition of "evenness" is a distribution that
minimizes a 1/r potential energy function of all the points, so that in
a global sense points are as "far away" from each other as possible.
Starting from an arbitrary distribution, we can either use numerical
minimization algorithms or point repulsion, in which all the points
are considered to repel each other according to a 1/r^2 force law and
dynamics are simulated. The algorithm is run until some convergence
criterion is satisfied, and the resulting distribution is our approximation.
Jon Leech developed code to do just this. It is available at
ftp://ftp.cs.unc.edu/pub/users/jon/points.tar.gz.
See his README file for more information. His distribution includes
sample output files for various n < 300, which may be used offtheshelf
if that is all you need.
Another method that is simpler than the above, but attains less
uniformity, is based on spacing the points along a spiral that
encircles the sphere.
Code available from links at http://cs.smith.edu/~orourke/.

Subject 6.07: What are coordinates for the vertices of an icosohedron?
Data on various polyhedra is available at
http://cm.belllabs.com/netlib/polyhedra/index.html, or
http://netlib.belllabs.com/netlib/polyhedra/index.html, or
http://www.netlib.org/polyhedra/index.html
For the icosahedron, the twelve vertices are:
(+ 1, 0, +t), (0, +t, +1), and (+t, +1, 0)
where t = (sqrt(5)1)/2, the golden ratio.
Reference: "Beyond the Third Dimension" by Banchoff, p.168.

Subject 6.08: How do I generate random points on the surface of a sphere?
There are several methods. Perhaps the easiest to understand is
the "rejection method": generate random points in an origin
centered cube with opposite corners (r,r,r) and (r,r,r).
Reject any point p that falls outside of the sphere of radius r.
Scale the vector to lie on the surface. Because the cube to sphere
volume ratio is pi/6, the average number of iterations before an
acceptable vector is found is 6/pi = 1.90986. This essentially
doubles the effort, and makes this method slower than the "trig
method." A timing comparison conducted by Ken Sloan showed that
the trig method runs in about 2/3's the time of the rejection method.
He found that methods based on the use of normal distributions are
twice as slow as the rejection method.
The trig method. This method works only in 3space, but it is
very fast. It depends on the slightly counterintuitive fact (see
proof below) that each of the three coordinates is uniformly
distributed on [1,1] (but the three are not independent,
obviously). Therefore, it suffices to choose one axis (Z, say)
and generate a uniformly distributed value on that axis. This
constrains the chosen point to lie on a circle parallel to the
XY plane, and the obvious trig method may be used to obtain the
remaining coordinates.
(a) Choose z uniformly distributed in [1,1].
(b) Choose t uniformly distributed on [0, 2*pi).
(c) Let r = sqrt(1z^2).
(d) Let x = r * cos(t).
(e) Let y = r * sin(t).
This method uses uniform deviates (faster to generate than normal
deviates), and no set of coordinates is ever rejected.
Here is a proof of correctness for the fact that the zcoordinate is
uniformly distributed. The proof also works for the x and y
coordinates, but note that this works only in 3space.
The area of a surface of revolution in 3space is given by
A = 2 * pi * int_a^b f(x) * sqrt(1 + [f'(x}]^2) dx
Consider a zone of a sphere of radius R. Since we are integrating in
the z direction, we have
f(z) = sqrt(R^2  z^2)
f'(z) = z / sqrt(R^2z^2)
1 + [f'(z)]^2 = r^2 / (R^2z^2)
A = 2 * pi * int_a^b sqrt(R^2z^2) * R/sqrt(R^2z^2) dz
= 2 * pi * R int_a^b dz
= 2 * pi * R * (ba)
= 2 * pi * R * h,
where h = ba is the vertical height of the zone. Notice how the integrand
reduces to a constant. The density is therefore uniform.
Here is simple C code implementing the trig method:
void SpherePoints(int n, double X[], double Y[], double Z[])
{
int i;
double x, y, z, w, t;
for( i=0; i< n; i++ ) {
z = 2.0 * drand48()  1.0;
t = 2.0 * M_PI * drand48();
w = sqrt( 1  z*z );
x = w * cos( t );
y = w * sin( t );
printf("i=%d: x,y,z=\t%10.5lf\t%10.5lf\t%10.5lf\n", i, x,y,z);
X[i] = x; Y[i] = y; Z[i] = z;
}
}
A complete package is available at
ftp://cs.smith.edu/pub/code/sphere.tar.gz (4K),
reachable from http://cs.smith.edu/~orourke/ .
References:
[Knuth, vol. 2]
[Graphics Gems IV] "Uniform Random Rotations"

Subject 6.09: What are Plucker coordinates?
A common convention is to write umlauted u as "ue", so you'll also see
"Pluecker". Lines in 3D can easily be given by listing the coordinates of
two distinct points, for a total of six numbers. Or, they can be given as
the coordinates of two distinct planes, eight numbers. What's wrong with
these? Nothing; but we can do better. Pluecker coordinates are, in a sense,
halfway between these extremes, and can trivially generate either. Neither
extreme is as efficient as Pluecker coordinates for computations.
When Pluecker coordinates generalize to Grassmann coordinates, as laid
out beautifully in [Hodge], Chapter VII, the determinant definition
is clearly the one to use. But 3D lines can use a simpler definition.
Take two distinct points on a line, say
P = (Px,Py,Pz)
Q = (Qx,Qy,Qz)
Think of these as vectors from the origin, if you like. The Pluecker
coordinates for the line are essentially
U = P  Q
V = P x Q
Except for a scale factor, which we ignore, U and V do not depend on the
specific P and Q! Cross products are perpendicular to their factors, so we
always have U.V = 0. In [Stolfi] lines have orientation, so are the same
only if their Pluecker coordinates are related by a positive scale factor.
As determinants of homogeneous coordinates, begin with the 4x2 matrix
[ Px Qx ] row x
[ Py Qy ] row y
[ Pz Qz ] row z
[ 1 1 ] row w
Define Pluecker coordinate Gij as the determinant of rows i and j, in that
order. Notice that Giw = Pi  Qi, which is Ui. Now let (i,j,k) be a cyclic
permutation of (x,y,z), namely (x,y,z) or (y,z,x) or (z,x,y), and notice
that Gij = Vk. Determinants are antisymmetric in the rows, so Gij = Gji.
Thus all possible Pluecker line coordinates are either zero (if i=j) or
components of U or V, perhaps with a sign change. Taking the w component
of a vector as 0, the determinant form will operate just as well on a
point P and vector U as on two points. We can also begin with a 2x4 matrix
whose rows are the coefficients of homogeneous plane equations, E.P=0,
from which come dual coordinates G'ij. Now if (h,i,j,k) is an even
permutation of (x,y,z,w), then Ghi = G'jk. (Just swap U and V!)
Got Pluecker, want points? No problem. At least one component of U is
nonzero, say Ui, which is Giw. Create homogeneous points Pj = Gjw + Gij,
and Qj = Gij. (Don't expect the P and Q that we started with, and don't
expect w=1.) Want plane equations? Let (i,j,k,w) be an even permutation of
(x,y,z,w), so G'jk = Giw. Then create Eh = G'hk, and Fh = G'jh.
Example: Begin with P = (2,4,8) and Q = (2,3,5). Then U = (0,1,3) and
V = (4,6,2). The direct determinant forms are Gxw=0, Gyw=1, Gzw=3,
Gyz=4, Gzx=6, Gxy=2, and the dual forms are G'yz=0, G'zx=1, G'xy=3,
G'xw=4, G'yw=6, G'zw=2. Take Uz = Gzw = G'xy = 3 as a suitable nonzero
element. Then two planes meeting in the line are
(G'xy G'yy G'zy G'wy).P = 0
(G'xx G'xy G'xz G'xw).P = 0
That is, a point P is on the line if it satisfies both these equations:
3 Px + 0 Py + 0 Pz  6 Pw = 0
0 Px + 3 Py  1 Pz  4 Pw = 0
We can also easily determine if two lines meet, or if not, how they pass.
If U1 and V1 are the coordinates of line 1, U2 and V2, of line 2, we look
at the sign of U1.V2 + V1.U2. If it's zero, they meet. The determinant form
reveals even permutations of (x,y,z,w):
G1xw G2yz + G1yw G2zx + G1zw G2xy + G1yz G2xw + G1zx p2yw + G1xy p2zw
Two oriented lines L1 and L2 can interact in three different ways:
L1 might intersect L2, L1 might go clockwise around L2, or L1 might go
counterclockwise around L2. Here are some examples:
 L2  L2  L2
L1  L1  L1 
+> > >
  
V V V
intersect counterclockwise clockwise
 L2  L2  L2
L1  L1  L1 
<+ < <
  
V V V
The first and second rows are just different views of the same lines,
once from the "front" and once from the "back." Here's what they might
look like if you look straight down line L2 (shown here as a dot).
L1 >
o> L1 o L1 o
>
intersect counterclockwise clockwise
The Pluecker coordinates of L1 and L2 give you a quick way to test
which of the three it is.
cw: U1.V2 + V1.U2 < 0
ccw: U1.V2 + V1.U2 > 0
thru: U1.V2 + V1.U2 = 0
So why is this useful? Suppose you want to test if a ray intersects
a triangle in 3space. One way to do this is to represent the ray and
the edges of the triangle with Pluecker coordinates. The ray hits the
triangle if and only if it hits one of the triangle's edges, or it's
"clockwise" from all three edges, or it's "counterclockwise" from all
three edges. For example...
o _
 \ ...in this picture, the ray
 \ is oriented counterclockwise
 \ > from all three edges, so it
 \ must intersect the triangle.
v \
o> o
Using Pluecker coordinates, ray shooting tests like this take only
a few lines of code.
Grassmann coordinates allow similar methods to be used for points,
lines, planes, and so on, and in a space of any dimension. In the
case of lines in 2D, only three coordinates are required:
Gxw = PxQx = G'y
Gyw = PyQy = G'x
Gxy = PxQyPyQx = G'w
Since P and Q are distinct, Giw is nonzero for i = x or y. Then
(Gix,Giy,Giw) is a point P1 on L
(Gxw,Gyw,Gww)+P1 is a point P2 on L
(G'x,G'y,G'w).P = 0 is an equation for L
Two lines in a 2D perspective plane always meet, perhaps in an
ideal point (meaning they're parallel in ordinary 2D). Calling
their homogeneous point of intersection P, the computation from
Grassmann coordinates nicely illustrates the convenience. (See
Subj 1.03, "How do I find intersections of 2 2D line segments?")
For i = x,y,w, set
Pi = G'j H'k ^Ö G'k H'j, where (i,j,k) is even
See [Hodge] for a thorough discussion of the theory, [Stolfi] for
a little theory with a concise implementation for low dimensions
(see Subj. 0.04), and these articles for further discussion:
J. Erickson, Pluecker Coordinates, Ray Tracing News 10(3) 1997,
http://www.acm.org/tog/resources/RTNews/html/rtnv10n3.html.
Ken Shoemake, Pluecker Coordinate Tutorial,
Ray Tracing News 11(1) 1998,
http://www.acm.org/tog/resources/RTNews/html/rtnv11n1.html.

Section 7. Contributors

Subject 7.01: How can you contribute to this FAQ?
Send email to orourke@cs.smith.edu with your suggestions, possible
topics, corrections, or pointers to information.

Subject 7.02: Contributors. Who made this all possible.
Jens Alfke
Nina Amenta
Leen Ammeraal
Scott Anguish
Ian Ashdown
Barak
Brad Barber
James Beech
David Bouman
Paul Bourke
Lars Brinkhoff
Andrew Bromage
Brent Burley
R. Kevin Burton
Gene Caldwell
Ken Clarkson
Robert Day
Tamal Dey
Martin Dillon
Thomas Djafari
Dave Eberly
John Eickemeyer
John E (Edward) Ellis
Jeff Erickson
Ata Etemadi
Hugh Fisher
David N. Fogel
Arne K. Frick
Olexandr Frantchuk
Robert W. Fuentes
Komei Fukuda
William Gibbons
Normand Grégoire
Eric Haines
Jeff Hameluck
Sandy Harris
Luiz Henrique de Figueiredo
Steve Hollasch
Bill Jones
Richard Kinch
Craig Kolb
Stefan Krause
Piyush Kumar
Steve Lamont
Ben Landon
Erik Larsen
Jon Leech
Michael V. Leonov
Sum Lin
Alan J. Livingston
Sebastien Loisel
Fritz Lott
Jacob Marner
Marc Christopher Martin
John McNamara
Samuel Murphy
Alan Murta
S. N. Muthukrishnan
Andrew Myers
David Nixon
Aaron Orenstein
Joseph O'Rourke
Samuel S. Paik
Leonidas Palios
Amitha Perera
Brian Peters
Lavoie Philippe
Christopher Phillips
Tom Plunket
Aaron Quigley
Rudi Bjxrn Rasmussen
Greg Roelofs
Christian von Roques
Dave Seaman
Jonathan R. Shewchuk
Rainer Michael Schmid
Klamer Schutte
Andrzej Serafin
ZhengYu Shan
James Sharman
Ken Shoemake
Jeff Somers
Jon Stone
Dan Sunday
Seth Teller
Saurabh Tendulkar
Yael "YoeL" Touboul
Anson Tsao
Bob van Manen
Remco Veltkamp
Jim Ward
Jason Weiler
Karsten Weiss
Stefan Wolfrum
Previous Editors:
Jon Stone
Anson Tsao
