From: orourke@cs.smith.edu (Joseph O'Rourke)Newsgroups: comp.graphics.algorithmsSubject: comp.graphics.algorithms Frequently Asked QuestionsDate: Sat, 15 Feb 2003 16:16:44 +0000 (UTC)Sender: orourke@grendel.csc.smith.eduMessage-ID: Reply-To: orourke@cs.smith.edu(Joseph O'Rourke)Keywords: faq computer graphics algorithms geometryX-Content-Currency: This FAQ changes regularly. See item (0.03) for instructionson how to obtain a new copy via FTP.Posted-By: auto-faq 3.3.1 beta (Perl 5.006)Archive-name: graphics/algorithms-faqPosting-Frequency: bi-weeklyWelcome to the FAQ for comp.graphics.algorithms!Thanks to all who have contributed.  Corrections and contributions(to orourke@cs.smith.edu) always welcome.----------------------------------------------------------------------This article is Copyright 2003 by Joseph O'Rourke.  It may be freelyredistributed in its entirety provided that this copyright notice isnot removed.----------------------------------------------------------------------Changed items this posting (|): 5.04New     items this posting (+): none----------------------------------------------------------------------History of Changes (approx. last six months):----------------------------------------------------------------------Changes in 15 Feb 03 posting:5.04: Fixed broken link in clipping article. [Thanks to Keith Forbes.]Changes in  1 Feb 03 posting:0.04: Ashdown Radiosity back in print. [Thanks to Ian Ashdown.]0.06: Update BSP FAQ links [Thanks to Ken Shoemake.]0.07: Update CGAL links in source article.3.11: Broken link re course based on Perlin's Noise book. [Thanks to Mikkel Gjoel.]6.01: Add CGAL link to Voronoi source article. [Thanks to Andreas Fabri.]7.02: All contributor email addresses removed to protect them from spam.Changes in 15 Jan 03 posting:0.06: Query re reality.sgi.com/bspfaq/ [Thanks to <<a href="mailto:warp@subphase.de" style="color: rgb(51, 102, 153); text-decoration: none;">warp@subphase.de
>.]0.07: Update moved link on WINGED.ZIP. [Thanks to Ben Landon.]0.07: Update Ferrar's ++ 3D rendering library link. [Thanks to F.Iannarilli, Jr.]1.06: Added ref to AT&T Graphviz. [Thanks to Michael Meire.]2.08: Fix sloan ear-clipping link.  [Thanks to logicalink@juno.com.]5.09: Update moved link re caustics. [Thanks to Ben Landon.]5.18: Formula for distance between two 3D lines. [Thanks to Daniel Zwick.]5.27: New article on transforming normals by Ken Shoemake.6.08: Random points on sphere in terms of longitude & latitude [Thanks to Uffe Kousgaard.]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 triangulatedsurfaces.5.23: Updated Cocone software description.Changes in 15 Feb 02 posting:5.03: Noted that Sutherland-Hodgman 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 restoredSunday'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 algorithmsthat 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 Information0.01: Charter of comp.graphics.algorithms0.02: Are the postings to comp.graphics.algorithms archived?0.03: How can I get this FAQ?0.04: What are some must-have 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 Computations2.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 Computations3.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 Computations4.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 computations5.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?5.27: How can I transform normals?6. Geometric Structures and Mathematics6.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. Contributors7.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.algorithmscomp.graphics.algorithms is an unmoderated newsgroup intended as a forumfor the discussion of the algorithms used in the process of generatingcomputer graphics.  These algorithms may be recently proposed inpublished journals or papers, old or previously known algorithms, orhacks used incidental to the process of computer graphics.  The scope ofthese 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 ballsand lime jello.It is hoped that this group will serve as a forum for programmers andresearchers to exchange ideas and ask questions on recent papers orcurrent 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; seealt.binaries.pictures
* FAQalt.binaries.pictures.utilities[now degenerated to pic postings]alt.graphics.pixutils(image format translation)comp.sources.misc (image viewing source code)sci.image.processingcomp.graphics.apps.softimagefj.comp.image- for requests for compression software; for these try:alt.comp.compressioncomp.compressioncomp.compression.research- specifically for game development; for this try:comp.games.development.programming.misccomp.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 easiestway to get it is to search back in your news reader for the mostrecent posting, with Subject:comp.graphics.algorithms Frequently Asked QuestionsIt is posted to comp.graphics.algorithms, and cross-posted tonews.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 thelatest 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/algorithms-faqftp://mirror.seas.gwu.edu/pub/rtfm/comp/graphics/algorithms/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 must-have books on graphics algorithms?The keywords in brackets are used to refer to the books in laterquestions.  They generally refer to the first author except whereit 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, Addison-Wesley1990, ISBN 0-201-12110-7;Computer Graphics: Principles and Practice, C versionJ.D. Foley,  A. van Dam, S.K. Feiner, J.F. Hughes, Addison-WesleyISBN: 0-201-84840-6, 1996, 1147 pp.[Rogers:Procedural]Procedural Elements for Computer Graphics, Second EditionDavid F. Rogers, WCB/McGraw Hill 1998, ISBN 0-07-053548-5[Rogers:Mathematical]Mathematical Elements for Computer Graphics 2nd Ed.,David F. Rogers and J. Alan Adams, McGraw Hill 1990, ISBN0-07-053530-2[Watt:3D]_3D Computer Graphics, 2nd Edition_,Alan Watt, Addison-Wesley 1993, ISBN 0-201-63186-5[Glassner:RayTracing]An Introduction to Ray Tracing,Andrew Glassner (ed.), Academic Press 1989, ISBN 0-12-286160-4[Gems I]Graphics Gems,Andrew Glassner (ed.), Academic Press 1990, ISBN 0-12-286165-5http://www.graphicsgems.org/for all the Gems.[Gems II]Graphics Gems II,James Arvo (ed.), Academic Press 1991, ISBN 0-12-64480-0[Gems III]Graphics Gems III,David Kirk (ed.), Academic Press 1992, ISBN 0-12-409670-0 (withIBM disk) or 0-12-409671-9 (with Mac disk)See also "AP Professional Graphics CD-ROM Library,"Academic Press,  ISBN 0-12-059756-X, which contains Gems I-III.[Gems IV]Graphics Gems IV,Paul S. Heckbert (ed.), Academic Press 1994, ISBN 0-12-336155-9(with IBM disk) or 0-12-336156-7 (with Mac disk)[Gems V]Graphic Gems V,Alan W. Paeth (ed.), Academic Press 1995, ISBN 0-12-543455-3(with IBM disk)[Watt:Animation]Advanced Animation and Rendering Techniques,Alan Watt, Mark Watt, Addison-Wesley 1992, ISBN 0-201-54412-1(Unofficial) errata: http://www.rolemaker.dk/other/AART/[Bartels]An Introduction to Splines for Use in Computer Graphics andGeometric Modeling,Richard H. Bartels, John C. Beatty, Brian A. Barsky, 1987, ISBN0-934613-27-3[Farin]Curves and Surfaces for Computer Aided Geometric Design:A Practical Guide, 4th Edition, Gerald E. Farin, Academic Press1996. ISBN 0122490541.[Prusinkiewicz]The Algorithmic Beauty of Plants,Przemyslaw W. Prusinkiewicz, Aristid Lindenmayer, Springer-Verlag,1990, ISBN 0-387-97297-8, ISBN 3-540-97297-8[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 Professional1993, ISBN 0-12-178270-0 [limited reprint 1999][Ashdown]Radiosity: A Programmer's PerspectiveIan Ashdown, John Wiley & Sons 1994, ISBN 0-471-30444-1, 498 pp.Back in print, Jan 2003.  See www.helios32.com.[sillion]Radiosity & Global IlluminationFrancois X. Sillion snd Claude Puech, Morgan Kaufmann 1994, ISBN1-55860-277-1, 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 0-12-228730-4, Includes CD-ROM.[Schroeder]Visualization Toolkit, 2nd Edition, The: An Object-Oriented Approach to3-D Graphics (Bk/CD) (Professional Description)William J. Schroeder,  Kenneth Martin, and Bill Lorensen,Prentice-Hall 1998, ISBN: 0-13-954694-4See Subject 0.07 for source.[Anderson]PC Graphics UnleashedScott Anderson. SAMS Publishing, ISBN 0-672-30570-4[Ammeraal]Computer Graphics for Java Programmers,Leen Ammeraal, John Wiley 1998, ISBN 0-471-98142-7.Additional information at http://home.wxs.nl/~ammeraal/.[Eberly]3D Game Engine Design: A Practical Approach to Real-Time 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 ISBN1-56881-000-8[Jain]Fundamentals of Image Processing,Anil K. Jain, Prentice-Hall 1989, ISBN 0-13-336165-9[Castleman]Digital Image Processing,Kenneth R. Castleman, Prentice-Hall 1996, ISBN(Cloth): 0-13-211467-4(Description and errata at: "http://www.phoenix.net/~castlman")[Pratt]Digital Image Processing, Second Edition,William K. Pratt, Wiley-Interscience 1991, ISBN 0-471-85766-1[Gonzalez]Digital Image Processing (3rd Ed.),Rafael C. Gonzalez, Paul Wintz, Addison-Wesley 1992, ISBN0-201-50803-6[Russ]The Image Processing Handbook (3rd Ed.),John C. Russ, CRC Press and IEEE Press 1998, ISBN 0-8493-2532-3[Russ & Russ]The Image Processing Tool Kit v. 3.0Chris Russ and John Russ, Reindeer Games Inc. 1999, ISBN 1-928808-00-X[Wolberg]Digital Image Warping,George Wolberg, IEEE Computer Society Press Monograph 1990, ISBN0-8186-8944-7Computational geometry----------------------[Bowyer]A Programmer's Geometry,Adrian Bowyer, John Woodwark, Butterworths 1983,ISBN 0-408-01242-0 PbkOut of print, but see:Introduction to Computing with Geometry,Adrian Bowyer and John Woodwark, 1993ISBN 1-874728-03-8.  Available in PDF:http://www.inge.com/pubs/index.htm[Farin & Hansford]The Geometry Toolbox for Graphics and Modelingby Gerald E. Farin, Dianne HansfordA K Peters Ltd; ISBN: 1568810741[O'Rourke (C)]Computational Geometry in C (2nd Ed.)Joseph O'Rourke, Cambridge University Press 1998,ISBN 0-521-64010-5 Pbk, ISBN 0-521-64976-5 HbkAdditional information and code at http://cs.smith.edu/~orourke/.[O'Rourke (A)]Art Gallery Theorems and AlgorithmsJoseph O'Rourke, Oxford University Press 1987,ISBN 0-19-503965-3.[Goodman & O'Rourke]Handbook of Discrete and Computational GeometryJ. E. Goodman and J. O'Rourke, editors.CRC Press LLC, July 1997.ISBN:0-8493-8524-5Additional information at http://cs.smith.edu/~orourke/ .[Samet:Application]Applications of Spatial Data Structures:  Computer Graphics,Image Processing, and GIS,Hanan Samet, Addison-Wesley, Reading, MA, 1990.ISBN 0-201-50300-0.[Samet:Design & Analysis]The Design and Analysis of Spatial Data Structures,Hanan Samet, Addison-Wesley, Reading, MA, 1990.ISBN 0-201-50255-0.[Mortenson]Geometric Modeling,Michael E. Mortenson, Wiley 1985, ISBN 0-471-88279-8[Preparata]Computational Geometry: An Introduction,Franco P. Preparata, Michael Ian Shamos, Springer-Verlag 1985,ISBN 0-387-96131-3[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 ApplicationsM. de Berg and M. van Kreveld and M. Overmars and O. SchwarzkopfSpringer-Verlag, Berlin, 1997.[Stolfi]Oriented Projective Geometry: A Framework for Geometric ComputationsAcademic Press, 1991.[Hodge]Methods of Algebraic Geometry, Volume 1W.V.D. Hodge and D. Pedoe, Cambridge, 1994.ISBN 0-521-469007-4 Paperback[Tamassia et al 199?]Graph Drawing: Algorithms for the Visualization of GraphsPrentice Hall; ISBN: 0133016153Algorithms books with chapters on computational geometry--------------------------------------------------------[Cormen et al.]Introduction to Algorithms,T. H. Cormen, C. E. Leiserson, R. L. Rivest,The MIT Press, McGraw-Hill, 1990.[Mehlhorn]Data Structures and Algorithms,K. Mehlhorn,Springer-Verlag, 1984.[Sedgewick]R. Sedgewick,Algorithms,Addison-Wesley, 1988.Solid Modelling---------------[Mantyla]Introduction to Solid ModelingMartti Mantyla, Computer Science Press 1988,ISBN 07167-8015-1----------------------------------------------------------------------Subject 0.05: Are there any online references?The computational geometry community maintains its ownbibliography of publications in or closely related to thatsubject.  Every four months, additions and corrections aresolicited from users, after which the database is updated andreleased anew.  As of 7 Nov 200, it contained 13485 bib-texentries.  See Jeff Erickson's page on "Computational GeometryBibliographies":http://compgeom.cs.uiuc.edu/~jeffe/compgeom/biblios.html#geombibThe bibliography can be retrieved from:ftp://ftp.cs.usask.ca/pub/geometry/geombib.tar.gz- bibliography properftp://ftp.cs.usask.ca/pub/geometry/o-cgc19.ps.gz- overview publishedin '93 in SIGACT News and the Internat. J. Comput.  Geom. Appl.ftp://ftp.cs.usask.ca/pub/geometry/ftp-hints- detailed retrieval infoUniversitat Politecnica de Catalunya maintains a search engine at:http://www-ma2.upc.es/~geomc/geombibe.htmlThe ACM SIGGRAPH Online Bibliography Project, byStephen Spencer (biblio@siggraph.org
).The database is available for anonymous FTP from theftp://siggraph.org/publications/ directory.  Pleasedownload and examine the file READ_ME in that directory for morespecific information concerning the database.'netlib' is a useful source for algorithms, member inquiries forSIAM, and bibliographic searches.  For information, send mail tonetlib@ornl.gov
, with "send index" in the body of the mailmessage.You can also find free sources for numerical computation in C viaftp://ftp.usc.edu/pub/C-numanal/. In particular, grabnumcomp-free-c.gz in that directory.Check out Nick Fotis's computer graphics resources FAQ -- it'spacked with pointers to all sorts of great computer graphicsstuff.  This FAQ is posted biweekly to comp.graphics
.This WWW page contains links to a large numberof computer graphic related pages:http://www.dataspace.com:84/vlib/comp-graphics.htmlThere's a Computer Science Bibliography Server at:http://glimpse.cs.arizona.edu:1994/bib/with Computer Graphics, Vision and Radiosity sectionsA 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.htmlThe University of Manchester NURBS Library:ftp://unix.hensa.ac.uk/pub/misc/unix/nurbs/For an implementation of Seidel's algorithm for fast trapezoidationand triangulation of polygons. You can get the code from:ftp://ftp.cs.unc.edu/pub/users/narkhede/triangulation.tar.gzRay tracing bibliography:http://www.acm.org/tog/resources/bib/Quaternions and other comp sci curiosities:ftp://ftp.netcom.com/pub/hb/hbaker/hakmem/Directory of Computational Geometry Software,collected by Nina Amenta (nina@cs.utexas.edu
)Nina Amenta is maintaining a WWW directory to computationalgeometry 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, smallestenclosing ball of a point set and other stuff.http://www.geom.umn.edu/software/cglist/A compact reference for real-time 3D computer graphics programming:http://www.math.mcgill.ca/~loisel/RADBIB is a comprehensive bibliography of radiosity andrelated global illumination papers, articles, andbooks. 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 andComputer Graphics:  http://visinfo.zib.de/3D Object Intersection: http://www.realtimerendering.com/int/This page presents information about a wide variety of 3D object/objectintersection 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 resourcesare 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 Wadehttp://www.andrew.cmu.edu/~rrost/bsp/ftp://ftp.sgi.com/other/bspfaq/faq/bspfaq.htmlftp://ftp.sgi.com/other/bspfaq/faq/bspfaq.txtand seeftp://ftp.sgi.com/other/bspfaq/Gamma and Color FAQs by Charles A. Poynton hasftp://ftp.inforamp.net/pub/users/poynton/doc/colour/http://www.inforamp.net/~poynton/The documents are mirrored in Darmstadt, Germany atftp://ftp.igd.fhg.de/pub/doc/colour/----------------------------------------------------------------------Subject 0.07: Where is all the source?Graphics Gems source code.http://www.graphicsgems.orgThis site is now the offical distribution site for Graphics Gems code.Master list of Computational Geometry software:http://www.geom.umn.edu/software/cglistDescribed in [Goodman & O'Rourke], Chap. 52.Jeff Erikson's software list:http://compgeom.cs.uiuc.edu/~jeffe/compgeom/compgeom.htmlDave Eberly's extensive collection of free geometry, graphics,and image processing software:http://www.magic-software.com/General 'stuff'ftp://wuarchive.wustl.edu/graphics/There are a number of interesting items inhttp://graphics.lcs.mit.edu/~sethincluding:- Code for 2D Voronoi, Delaunay, and Convex hull- Mike Hoymeyer's implementation of Raimund Seidel'sO( d! n ) time linear programming algorithm forn constraints in d dimensions- geometric models of UC Berkeley's new computer sciencebuildingSources to "Computational Geometry in C", by J. O'Rourkecan be found at http://cs.smith.edu/~orourke/books/compgeom.htmlor ftp://cs.smith.edu/pub/compgeom.Greg Ferrar's C++ 3D rendering library is available athttp://www.flowerfire.com/ferrar/Graph3D.htmlTAGL is a portable and extensible library that provides a subsetof Open-GL functionalities.ftp://sunsite.unc.edu/pub/packages/programming/graphics/Try ftp://x2ftp.oulu.fifor /pub/msdos/programming/docs/graphpro.lzh byMichael Abrash. His XSharp package has an implementation of XiaoulinWu's anti-aliasing algorithm (in C).Example sources for BSP tree algorithms can be found athttp://reality.sgi.com/bspfaq/, item 24.Mel Slater (mel@dcs.qmw.ac.uk
) also made some implementations ofBSP trees and shadows for static scenes using shadow volumescode availablehttp://www.dcs.qmw.ac.uk/~mel/BSP.htmlftp://ftp.dcs.qmw.ac.uk/people/mel/BSPThe Visualization Toolkit (A visualization textbook, C++ libraryand Tcl-based interpreter) (see [Schroeder]):http://www.kitware.com/vtk.htmlWINGED.ZIP, a C++ implementation of Baumgart's winged-edge data structure:ftp://ftp.ledalite.com/pub/CGAL, the Computational Geometry Algorithms Library, is written in C++and is available at http://www.cgal.org
.CGAL contains algorithms and data structures for 2D computations(convex hull, Delaunay, constrained Delaunay, Voronoi diagram,regular traingulation, (weighted) Alpha shapes, polytope distance,boolean operations on polygons, decomposition of polygons inmonotone or convex parts, arrangements, etc.), 3D, and arbitrarydimensions.A C++ NURBS library written by Lavoie Philippe. Version 2.1.Results may be exported as POV-Ray, 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 isosurfacegeneration 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.tu-berlin.de/~ki/engines.htmlSee 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 tY = x sin t + y cos tAs a 2x2 matrix this is very simple.  If you want to rotate acolumn vector v by t degrees using matrix M, useM = [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 intocolumns and vice versa).  If you want to combine rotations, in 2Dyou can just add their angles, but in higher dimensions you mustmultiply 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 parameterr, which indicates P's position along AB, is computed by the dot productof AC and AB divided by the square of the length of AB:(1)     AC dot ABr = ---------||AB||^2r has the following meaning:r=0      P = Ar=1      P = Br<0      P is on the backward extension of ABr>1      P is on the forward extension of AB0The length of a line segment in d dimensions, AB is computed by:L = sqrt( (Bx-Ax)^2 + (By-Ay)^2 + ... + (Bd-Ad)^2)so in 2D:L = sqrt( (Bx-Ax)^2 + (By-Ay)^2 )and the dot product of two vectors in d dimensions, U dot V is computed:D = (Ux * Vx) + (Uy * Vy) + ... + (Ud * Vd)so in 2D:D = (Ux * Vx) + (Uy * Vy)So (1) expands to:(Cx-Ax)(Bx-Ax) + (Cy-Ay)(By-Ay)r = -------------------------------L^2The point P can then be found:Px = Ax + r(Bx-Ax)Py = Ay + r(By-Ay)And the distance from A to P = r*L.Use another parameter s to indicate the location along PC, with thefollowing meaning:s<0      C is left of ABs>0      C is right of ABs=0      C is on ABCompute s as follows:(Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay)s = -----------------------------L^2Then 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; itdepends on your application. If all you want is the intersectionpoint, the following should work:Let A,B,C,D be 2-space position vectors.  Then the directed linesegments AB & CD are given by:AB=A+r(B-A), r in [0,1]CD=C+s(D-C), s in [0,1]If AB & CD intersect, thenA+r(B-A)=C+s(D-C), orAx+r(Bx-Ax)=Cx+s(Dx-Cx)Ay+r(By-Ay)=Cy+s(Dy-Cy)  for some r,s in [0,1]Solving the above for r and s yields(Ay-Cy)(Dx-Cx)-(Ax-Cx)(Dy-Cy)r = -----------------------------  (eqn 1)(Bx-Ax)(Dy-Cy)-(By-Ay)(Dx-Cx)(Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay)s = -----------------------------  (eqn 2)(Bx-Ax)(Dy-Cy)-(By-Ay)(Dx-Cx)Let P be the position vector of the intersection point, thenP=A+r(B-A) orPx=Ax+r(Bx-Ax)Py=Ay+r(By-Ay)By examining the values of r & s, you can also determine someother limiting conditions:If 0<=r<=1 & 0<=s<=1, intersection existsr<0 or r>1 or s<0 or s>1 line segments do not intersectIf the denominator in eqn 1 is zero, AB & CD are parallelIf 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 y-axis, and overlap of the projected intervals checked.If the intersection point of the 2 lines are needed (lines in thiscontext mean infinite lines) regardless whether the two linesegments intersect, thenIf r>1, P is located on extension of ABIf r<0, P is located on extension of BAIf s>1, P is located on extension of CDIf s<0, P is located on extension of DCAlso note that the denominators of eqn 1 & 2 are identical.References:[O'Rourke (C)] pp. 249-51[Gems III] pp. 199-202 "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 representx and y coordinates. The coordinates of the center p=(p_0,p_1) ofthe 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 finite-radiuscircle through them exists.  Otherwise, the radius of the circle is:r^2 = (a_0 - p_0)^2 + (a_1 - p_1)^2Reference:[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 becomputed in O(n log n) time for n points.  The center lies onthe furthest-point Voronoi diagram.  Computing the diagram constrainsthe search for the center.  Constructing the diagram can be accomplishedby 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) 121--125.J. Rokne, "An Easy Bounding Circle" [Gems II] pp.14--16.----------------------------------------------------------------------Subject 1.06: Where can I find graph layout algorithms?See the paper by Eades and Tamassia, "Algorithms for DrawingGraphs: An Annotated Bibliography," Tech Rep CS-89-09, Dept.  ofCS, Brown Univ, Feb. 1989.A revised version of the annotated bibliography on graph drawingalgorithms by Giuseppe Di Battista, Peter Eades, and RobertoTamassia is now available fromftp://wilma.cs.brown.edu/pub/papers/compgeo/gdbiblio.tex.gzandftp://wilma.cs.brown.edu/pub/papers/compgeo/gdbiblio.ps.gzCommercial software includes the Graph Layout Toolkit from Tom SawyerSoftware http://www.tomsawyer.comand Northwoods Software's GO++http://www.nwoods.com/go/.Perhaps the best code is the AT&T Research Labs open C source:http://www.research.att.com/sw/tools/graphviz/----------------------------------------------------------------------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 by2 A( P ) = sum_{i=0}^{n-1} (x_i y_{i+1} - y_i x_{i+1}).Here n is the number of vertices of the polygon, and indexarithmetic is mod n, so that x_n = x_0, etc. A rearrangementof terms in this equation can save multiplications and operate oncoordinate differences, and so may be both faster and moreaccurate:2 A(P) = sum_{i=0}^{n-1} ( x_i  (y_{i+1} - y_{i-1}) )Here again modular index arithmetic is implied, with n=0 and -1=n-1.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 instead2 A(P) = sum_{i=1}^{n} ( x_i  (y_{i+1} - y_{i-1}) )References: [O' Rourke (C)] Thm. 1.3.3, p. 21; [Gems II] pp. 5-6:"The Area of a Simple Polygon", Jon Rokne.  Dan Sunday's explanation:http://GeometryAlgorithms.com/Archive/algorithm_0101/whereTo find the area of a planar polygon not in the x-y plane, use:2 A(P) = abs(N . (sum_{i=0}^{n-1} (v_i x v_{i+1})))where N is a unit vector normal to the plane. The `.' represents thedot product operator, the `x' represents the cross product operator,and abs() is the absolute value function.[Gems II] pp. 170-171:"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 centroidsof a partition of the polygon into triangles.  The centroid of atriangle is simply the average of its three vertices, i.e., ithas coordinates (x1 + x2 + x3)/3 and (y1 + y2 + y3)/3.  Thissuggests first triangulating the polygon, then forming a sumof the centroids of each triangle, weighted by the area ofeach triangle, the whole sum normalized by the total polygon area.This indeed works, but there is a simpler method:  the triangulationneed not be a partition, but rather can use positively andnegatively oriented triangles (with positive and negative areas),as is used when computing the area of a polygon.  This leads toa very simple algorithm for computing the centroid, based on asum 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 ofconsecutive edges of the polygon: (v1,v2), (v2,v3), etc.  The areaof 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.3-6; 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" byEric Haines [Gems IV]  pp. 24-46.  Now also athttp://www.erichaines.com/ptinpoly
.The code in the Sedgewick book Algorithms (2nd Edition, p.354) failsunder certain circumstances.  Seehttp://condor.informatik.Uni-Oldenburg.DE/~stueker/graphic/index.htmlfor a discussion.The essence of the ray-crossing 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 nowoutside the poly. If you have to cross again you know you are nowinside again; i.e., if you were inside the field to start with, the totalnumber of fence jumps you would make will be odd, whereas if you wereouside the jumps will be even.The code below is from Wm. Randolph Franklin <<a href="mailto:wrf@ecse.rpi.edu" style="color: rgb(51, 102, 153); text-decoration: none;">wrf@ecse.rpi.edu
>(see URL below) with some minor modifications for speed.  It returns1 for strictly interior points, 0 for strictly exterior, and 0 or 1for points on the boundary.  The boundary behavior is complex butdetermined; 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 = npol-1; i < npol; j = i++) {if ((((yp[i]<=y) && (y((yp[j]<=y) && (y(x < (xp[j] - xp[i]) * (y - yp[i]) / (yp[j] - yp[i]) + xp[i]))c = !c;}return c;}The code may be further accelerated, at some loss in clarity, byavoiding the central computation when the inequality can be deduced,and by replacing the division by a multiplication for those processorswith slow divides.  For code that distinguishes strictly interiorpoints from those on the boundary, see [O'Rourke (C)] pp. 239-245.For a method based on winding number, see Dan Sunday,"Fast Winding Number Test for Point Inclusion in a Polygon,"http://softsurfer.com/algorithms.htm
, March 2001.References:Franklin's code:http://www.ecse.rpi.edu/Homepages/wrf/research/geom/pnpoly.html[Gems IV]  pp. 24-46[O'Rourke (C)] Sec. 7.4.[Glassner:RayTracing]----------------------------------------------------------------------Subject 2.04: How do I find the intersection of two convex polygons?Unlike intersections of general polygons, which might produce aquadratic number of pieces, intersection of convex polygons of nand m vertices always produces a polygon of at most (n+m) vertices.There are a variety of algorithms whose time complexity is proportionalto this size, i.e., linear.  The first, due to Shamos and Hoey, isperhaps the easiest to understand.  Let the two polygons be P andQ.  Draw a horizontal line through every vertex of each.  Thispartitions each into trapezoids (with at most two triangles, oneat the top and one at the bottom).  Now scan down the two polygonsin concert, one trapezoid at a time, and intersect the trapezoidsif they overlap vertically.A more difficult-to-describe algorithm is in [O'Rourke (C)], pp. 252-262.This walks around the boundaries of P and Q in concert, intersectingedges.  An implementation is included in [O'Rourke (C)].----------------------------------------------------------------------Subject 2.05: How do I do a hidden surface test (backface culling) with 2D points?c = (x1-x2)*(y3-y2)-(y1-y2)*(x3-x2)x1,y1, x2,y2, x3,y3 = the first three points of the polygon.If c is positive the polygon is visible. If c is negative thepolygon is invisible (or the other way).----------------------------------------------------------------------Subject 2.06: How do I find a single point inside a simple polygon?Given a simple polygon, find some point inside it.  Here is a methodbased on the proof that there exists an internal diagonal, in[O'Rourke (C), 13-14].  The idea is that the midpoint of a diagonalis interior to the polygon.1. Identify a convex vertex v; let its adjacent vertices be a and b.2. For each other vertex q do:2a. If q is inside avb, compute distance to v (orthogonal to ab).2b. Save point q if distance is a new min.3. If no point is inside, return midpoint of ab, or centroid of avb.4. Else if some point inside, qv is internal: return its midpoint.Code for finding a diagonal is in [O'Rourke (C), 35-39], and is partof many other software packages.  See Subject 0.07: Where is all thesource?----------------------------------------------------------------------Subject 2.07: How do I find the orientation of a simple polygon?Compute the signed area (Subject 2.01).  The orientation iscounter-clockwise if this area is positive.A slightly faster method is based on the observation that it isn'tnecessary to compute the area.  Find the lowest vertex (or, ifthere is more than one vertex with the same lowest coordinate,the rightmost of those vertices) and then take the cross productof the edges fore and aft of it.  Both methods are O(n) for n vertices,but it does seem a waste to add up the total area when a single crossproduct (of just the right edges) suffices.  Code for this isavailable at ftp://cs.smith.edu/pub/code/polyorient.C(2K).The reason that the lowest, rightmost (or any other such extreme) pointworks is that the internal angle at this vertex is necessarily convex,strictly less than pi (even if there are several equally-lowest points).----------------------------------------------------------------------Subject 2.08: How can I triangulate a simple polygon?Triangulation of a polygon partitions its interior into triangleswith disjoint interiors.  Usually one restricts corners of thetriangles to coincide with vertices of the polygon, in which caseevery polygon of n vertices can be triangulated, and all triangulationscontain n-2 triangles, employing n-3 "diagonals" (chords betweenvertices that otherwise do not touch the boundary of the polygon).Triangulations can be constructed by a variety of algorithms,ranging from a naive search for noncrossing diagonals, which isO(n^4), to "ear" clipping, which is O(n^2) and relatively straightforwardto implement [O'Rourke (C), Chap. 1], to partitioning intomonotone polygons, which leads to O(n log n) time complexity[O'Rourke (C), Chap. 2; Overmars, Chap. 3], to several randomizedalgorithms---by Clarkson et al., by Seidel, and by Devillers, whichlead to O(n log* n) complexity---to Chazelle's linear-time algorithm,which has yet to be implemented.There is a tradeoff between code complexity and time complexity.Fortunately, several of the algorithms have been implemented and areavailable:Ear-clipping:http://cs.smith.edu/~orourke/books/compgeom.htmlftp://ftp.cis.uab.edu/pub/sloan/Software/triangulation/src/Seidel's Alg:http://www.cs.unc.edu/~dm/CODE/GEM/chapter.htmlftp://ftp.cs.unc.edu/pub/users/narkhede/triangulation.tar.gzhttp://reality.sgi.com/atul/code/chapter.htmlSee also the collection of triangulation links athttp://www.geom.umn.edu/software/cglist/References:[O'Rourke (C)][Overmars][Gems V]Clarkson, K., Tarjan, R., and VanWyk, C. A fast Las Vegas algorithm fortriangulating a simple polygon. Discrete and Computational Geometry,4(1):423--432, 1989.Clarkson, K., Cole, R., Tarjan, R. Randomized parallel algorithms fortrapezoidal diagrams. Int. J. Comp. Geom. Appl., 117--133, 1992.http://cm.bell-labs.com/cm/cs/who/clarkson/tri.htmlSeidel, R.  (1991), A simple and fast incremental randomized algorithm forcomputing trapezoidal decompositions and for triangulating polygons,Comput. Geom. Theory Appl. 1, 51--64.Devillers, O.  (1992), Randomization yields simple O(n log* n)algorithms for difficult Omega(n) problems,Internat. J. Comput.  Geom. Appl. 2(1), 97--111.Chazelle, B.  (1991), Triangulating a simple polygon in linear time,Discrete Comput. Geom. 6, 485--524.Held, M. (1998) "FIST: Fast Industrial-Strength Triangulation".http://www.cosy.sbg.ac.at/~held/projects/triang/triang.html----------------------------------------------------------------------Subject 2.09: How can I find the minimum area rectangle enclosing a set of points?First take the convex hull of the points.  Let the resulting convexpolygon be P.  It has been known for some time that the minimumarea rectangle enclosing P must have one rectangle side flush with(i.e., collinear with and overlapping) one edge of P.  This geometricfact was used by Godfried Toussaint to develop the "rotating calipers"algorithm in a hard-to-find 1983 paper, "Solving Geometric Problemswith the Rotating Calipers" (Proc. IEEE MELECON).  The algorithmrotates a surrounding rectangle from one flush edge to the next,keeping track of the minimum area for each edge.  It achieves O(n)time (after hull computation).  See the "Rotating Calipers Homepage"http://www.cs.mcgill.ca/~orm/rotcal.frame.html for a descriptionand applet.----------------------------------------------------------------------Section 3. 2D Image/Pixel Computations----------------------------------------------------------------------Subject 3.01: How do I rotate a bitmap?The easiest way, according to the comp.graphics faq, is to takethe rotation transformation and invert it. Then you just iterateover the destination image, apply this inverse transformation andfind which source pixel to copy there.A much nicer way comes from the observation that the rotationmatrix:R(T) = { { cos(T), -sin(T) }, { sin(T), cos(T) } }is formed my multiplying three matrices, namely:R(T) = M1(T) * M2(T) * M3(T)whereM1(T) = { { 1, -tan(T/2) },{ 0, 1         } }M2(T) = { { 1,      0    },{ sin(T), 1    } }M3(T) = { { 1, -tan(T/2) },{ 0,         1 } }Each transformation can be performed in a separate pass, andbecause these transformations are either row-preserving orcolumn-preserving, anti-aliasing is quite easy.Another fast approach is to perform first a column-preserving roation,and then a row-preserving rotation. For an image W pixels wide andH pixels high, this requires W+H BitBlt operations in comparison tothe brute-force rotation, which uses W*H SetPixel operations (and alot of multiplying).Reference:Paeth, A. W., "A Fast Algorithm for General Raster Rotation",Proceedings Graphics Interface '89, Canadian InformationProcessing Society, 1986, 77-81[Note - e-mail copies of this paper are no longer available][Gems I]----------------------------------------------------------------------Subject 3.02: How do I display a 24 bit image in 8 bits?[Gems I] pp. 287-293, "A Simple Method for Color Quantization:Octree Quantization"B. Kurz.  Optimal Color Quantization for Color Displays.Proceedings of the IEEE Conference on Computer Vision and PatternRecognition, 1983, pp. 217-224.[Gems II] pp. 116-125, "Efficient Inverse Color Map Computation"This describes an efficient technique tomap actual colors to a reduced color map,selected by some other technique describedin the other papers.[Gems II] pp. 126-133, "Efficient Statistical Computations forOptimal Color Quantization"Xiaolin Wu.  Color Quantization by Dynamic Programming andPrincipal Analysis.  ACM Transactions on Graphics, Vol. 11, No. 4,October 1992, pp 348-372.----------------------------------------------------------------------Subject 3.03: How do I fill the area of an arbitrary shape?"A Fast Algorithm for the Restoration of Images Based on ChainCodes Description and Its Applications", L.W. Chang & K.L. Leu,Computer Vision, Graphics, and Image Processing, vol.50,pp296-307 (1990)Heckbert, Paul S., Generic Convex Polygon Scan Conversion and Clipping,Graphics Gems, p. 84-86, code: p. 667-680, PolyScan/.Heckbert, Paul S., Concave Polygon Scan Conversion, Graphics Gems, p.87-91, code: p. 681-684, ConcaveScan.c.http://www.acm.org/tog/GraphicsGems/gems/PolyScan/http://www.acm.org/tog/GraphicsGems/gems/ConcaveScan.cFor filling a region of a bitmap, seeHeckbert, Paul S., A Seed Fill Algorithm, Graphics Gems, p. 275-277,code: p. 721-722, SeedFill.c. The code is athttp://www.acm.org/tog/GraphicsGems/gems/SeedFill.c[Gems I][Foley][Hearn]----------------------------------------------------------------------Subject 3.04: How do I find the 'edges' in a bitmap?A simple method is to put the bitmap through the filter:-1    -1    -1-1     8    -1-1    -1    -1This will highlight changes in contrast.  Then any part of thepicture where the absolute filtered value is higher than somethreshold is an "edge".A more appropriate edge detector for noisy images isdescribed by Van Vliet et al. "A nonlinear Laplace operatoras edge detector in noisy images", in Computer Vision,Graphics, and image processing 45, pp. 167-195, 1989.----------------------------------------------------------------------Subject 3.05: How do I enlarge/sharpen/fuzz a bitmap?Sharpening of bitmaps can be done by the following algorithm:I_enh(x,y) = I_fuz(x,y)-k*Laplace(I_fuz(x,y))or in words:  An image can be sharpened by subtracting a positivefraction k of the Laplace from the fuzzy image.One "Laplace" kernel, approximating a Laplacian operator, is:1   1   11  -8   11   1   1The following library implements Fast Gaussian Blurs:MAGIC: An Object-Oriented Library for Image Analysis by David EberlyThe library source code and the documentation (in Latex) are  athttp://www.magic-software.com/The code compiles on Unix systems using g++ and on PCs usingMicrosoft Windows 3.1 and Borland C++.  The fast Gaussian blurringis based on a finite difference method for solving s u_s = s^2 /nabla^2 uwhere s is the standard deviation of the Gaussian (t = s^2/2).  Ittakes advantage of geometrically increasing steps in s (rather thanlinearly increasing steps in t), thus getting to a larger "time" rapidly,but still retaining stability.  Section 4.5 of the  documentation containsthe algorithm description and implementation.A bitmap is a sampled image, a special case of a digital signal,and suffers from two limitations common to all digital signals.First, it cannot provide details at fine enough spacing to exactlyreproduce every continuous image, nor even more detailed sampledimages. And second, each sample approximates the infinitely finevariability of ideal values with a discrete set of ranges encodedin a small number of bits---sometimes just one bit per pixel. Mostbitmaps have another limitation imposed: The values cannot benegative. The resolution limitation is especially important, butsee "How do I display a 24 bit image in 8 bits?" for range issues.The ideal way to enlarge a bitmap is to work from the originalcontinuous image, magnifying and resampling it. The standard wayto do it in practice is to (conceptually) reconstruct a continuousimage from the bitmap, and magnify and resample that instead. Thiswill not give the same results, since details of the original havealready been lost, but it is the best approach possible given analready sampled image. More details are provided below.Both sharpening and fuzzing are examples of filtering. Even morespecifically, they can be both be accomplished with filters whichare linear and shift invariant. A crude way to sharpen along a row(or column) is to set output pixel B[n] to the difference of inputpixels, A[n]-A[n-1]. A similarly crude way to fuzz is to set B[n]to the average of input pixels, 1/2*A[n]+1/2*A[n-1]. In each casethe output is a weighted sum of input pixels, a "convolution". Oneimportant characteristic of such filters is that a sinusoid goingin produces a sinusoid coming out, one of the same frequency. Thusthe Fourier transform, which decomposes a signal into sinusoids ofvarious frequencies, is the key to analysis of these filters. Thesimplest (and most efficient) way to handle the two dimensions ofimages is to operate on first the rows then the columns (or viceversa). Fourier transforms and many filters allow this separation.A filter is linear if it satisfies two simple relations between theinput and output: scaling the input by a factor scales the outputby the same factor, and the sum of two inputs gives the sum of thetwo outputs. A filter is shift invariant if shifting the input up,down, left, or right merely shifts the output the same way. When afilter is both linear and shift invariant, it can be implemented asa convolution, a weighted sum. If you find the output of the filterwhen the input is a single pixel with value one in a sea of zeros,you will know all the weights. This output is the impulse responseof the filter. The Fourier transform of the impulse response givesthe frequency response of the filter. The pattern of weights readoff from the impulse response gives the filter kernel, which willusually be displayed (for image filters) as a 2D stencil array, andit is almost always symmetric around the center. For example, thefollowing filter, approximating a Laplacian (and used for detectingedges), is centered on the negative value.1/6     4/6     1/64/6   -20/6     4/61/6     4/6     1/6The symmetry allows a streamlined implementation. Suppose the inputimage is in A, and the output is to go into B. Then computeB[i][j] = (A[i-1][j-1]+A[i-1][j+1]+A[i+1][j-1]+A[i+1][j+1]+4.0*(A[i-1][j]+A[i][j-1]+A[i][j+1]+A[i+1][j])-20.0*A[i][j])/6.0Ideal blurring is uniform in all directions, in other words it hascircular symmetry. Gaussian blurs are popular, but the obvious codeis slow for wide blurs. A cheap alternative is the following filter(written for rows, but then applied to the columns as well).B[i][j] = ((A[i][j]*2+A[i][j-1]+A[i][j+1])*4+A[i][j-1]+A[i][j+1]-A[i][j-3]-A[i][j+3])/16For sharpening, subtract the results from the original image, whichis equivalent to using the following.B[i][j] = ((A[i][j]*2-A[i][j-1]-A[i][j+1])*4-A[i][j-1]-A[i][j+1]+A[i][j-3]+A[i][j+3])/16Credit for this filter goes to Ken Turkowski and Steve Gabriel.Reconstruction is impossible without some assumptions, and becauseof the importance of sinusoids in filtering it is traditional toassume the continuous image is made of sinusoids mixed together.That makes more sense for sounds, where signal processing began,than it does for images, especially computer images of charactershapes, sharp surface features, and halftoned shading. As pointedout above, often image values cannot be negative, unlike sinusoids.Also, real world images contain noise. The best noise suppressors(and edge detectors) are, ironically, nonlinear filters.The simplest way to double the size of an image is to use each ofthe original pixels twice in its row and in its column. For muchbetter results, try this instead. Put zeros between the originalpixels, then use the blurring filter given a moment ago. But youmight want to divide by 8 instead of 16 (since the zeros will dimthe image otherwise). To instead shrink the image by half (in bothvertical and horizontal), first apply the filter (dividing by 16),then throw away every other pixel. Notice that there are obviousoptimizations involving arithmetic with powers of two, zeros whichare in known locations, and pixels which will be discarded.----------------------------------------------------------------------Subject 3.06: How do I map a texture on to a shape?Paul S. Heckbert, "Survey of Texture Mapping", IEEE ComputerGraphics and Applications V6, #11, Nov. 1986, pp 56-67 revisedfrom Graphics Interface '86 versionEric A. Bier and Kenneth R. Sloan, Jr., "Two-Part TextureMappings", IEEE Computer Graphics and Applications V6 #9, Sept.1986, pp 40-53 (projection parameterizations)----------------------------------------------------------------------Subject 3.07: How do I detect a 'corner' in a collection of points?[Currently empty entry.]----------------------------------------------------------------------Subject 3.08: Where do I get source to display (raster font format)?ftp://oak.oakland.edu/SimTel/msdos/See also James Murray's graphics file formats FAQ:http://www.ora.com/centers/gff/gff-faq/index.htm----------------------------------------------------------------------Subject 3.09: What is morphing/how is it done?Morphing is the name that has come to be applied to the techniqueILM used in the movie "Willow", where one object changes intoanother by changing both its shape and picture detail. It was a2D image manipulation, and has been done in different ways. Thefirst method published was by Thad Beier at PDI. Michael Jacksonfamously used morphing in his music videos, notably "Black orWhite". The word is now used more generally.For more, see [Anderson], Chapter 3, page 59-90, and Beier'shttp://www.hammerhead.com/thad/morph.html----------------------------------------------------------------------Subject 3.10: How do I quickly draw a filled triangle?The easiest way is to render each line separately into an edgebuffer. This buffer is a structure which looks like this (in C):struct { int xmin, xmax; } edgebuffer[YDIM];There is one entry for each scan line on the screen, and eachentry is to be interpreted as a horizontal line to be drawn fromxmin to xmax.Since most people who ask this question are trying to write fastgames on the PC, I'll tell you where to find code.  Look at:ftp::/ftp.uwp.edu/pub/msdos/demos/programming/sourceftp::/ftp.luth.se/pub/msdos/demos (Sweden)ftp::/NCTUCCCA.edu.tw:/PC/uwp/demos/www.wit.com:/mirrors/uwp/pub/msdos/demos
">http://www.wit.com:/mirrors/uwp/pub/msdos/demosftp::/ftp.cdrom.com:/demosSee also Subject 3.03, which describes methods for filling polygons.----------------------------------------------------------------------Subject 3.11: 3D Noise functions and turbulence in Solid texturing.See:ftp://gondwana.ecr.mu.oz.au/pub/ftp://ftp.cis.ohio-state.edu/pub/siggraph92/siggraph92_C23.sharIn it there are implementations of Perlin's noise and turbulencefunctions, (By the man himself) as well as Lewis' sparseconvolution noise function (by D. Peachey) There is also some ofother stuff in there (Musgrave's Earth texture functions, and somestuff on animating gases by Ebert).SPD (Standard Procedural Databases) package:ftp://avalon.chinalake.navy.mil/utils/SPD/ftp://avalon.chinalake.navy.mil/utils/SPD/.Now moved to http://www.viewpoint.com/References:[Ebert]Noise, Hypertexture, Antialiasing and Gesture, (Ken Perlin) inChapter 6, (p.193-), The disk accompanying the book is availablefromftp://archive.cs.umbc.edu/pub/.For more info on this text/code see:http://www.cs.umbc.edu/~ebert/book/book.htmlFor examples from a current course based on this book, see:http://www.seas.gwu.edu/graphics/ProcTexCourse/Linke broken 21Jan03; will remove eventually if not fixed.[Watt:Animation]Three-dimensional Nocie, Chapter 7.2.1Simulating turbulance, Chapter 7.2.2----------------------------------------------------------------------Subject 3.12: How do I generate realistic sythetic textures?For fractal mountains, trees and sea-shells:SPD (Standard Procedural Databases) package:ftp://avalon.chinalake.navy.mil/utils/SPD/ftp://avalon.chinalake.navy.mil/utils/SPD/.Now moved to http://www.viewpoint.com/Reaction-Diffusion Algorithms:For an illustartion of the parameter space of a reactiondiffusion system, check out the Xmorphia page athttp://www.ccsf.caltech.edu/ismap/image.htmlReferences:[Ebert]Entire book devoted to this subject, with RenderMan(TM) and C code.[Watt:Animation]Procedural texture mapping and modelling, Chapter 7"Generating Textures on Arbitrary Surfaces Using Reaction-Diffusion"Greg Turk,   Computer Graphics, Vol. 25, No. 4, pp. 289-298July 1991 (SIGGRAPH '91)http://www.cs.unc.edu:80/~turk/reaction_diffusion/reaction_diffusion.htmlA list of procedural texture synthesis related web pageshttp://www.threedgraphics.com/pixelloom/tex_synth.html----------------------------------------------------------------------Subject 3.13: How do I convert between color models (RGB, HLS, CMYK, CIE etc)?References:[Watt:3D] pp. 313-354[Foley]   pp. 563-603----------------------------------------------------------------------Subject 3.14: How is "GIF" pronounced?"GIF" is an acronymn for "Graphics Interchange Format."  Despite thehard "G" in "Graphics," GIF is pronounced "JIF."  Although we don'thave a direct quote from the official CompuServe specificationreleased June 1987, here is a quote from related CompuServe documentation,for CompuShow, a DOS-based image viewer used shortly thereafter:"The GIF (Graphics Interchange Format), pronounced "JIF", wasdesigned by CompuServe ..."We also have a report that the principal author of the GIF spec,Steve Wilhite, says "it's pronounced JIF (like the peanut butter."See alsohttp://www.60-seconds.com/articles/86.html----------------------------------------------------------------------Section 4. Curve Computations----------------------------------------------------------------------Subject 4.01: How do I generate a Bezier curve that is parallel to another Bezier?You can't.  The only case where this is possible is when theBezier can be represented by a straight line.  And then theparallel 'Bezier' can also be represented by a straight line.The situation is different for the broader class of rationalBezier curves.  For example, these can represent circular arcs,and a parallel offset is just a concentric circular arc, alsorepresentable as a rational Bezier.----------------------------------------------------------------------Subject 4.02: How do I split a Bezier at a specific value for t?A Bezier curve is a parametric polynomial function from theinterval [0..1] to a space, usually 2D or 3D.  Common Beziercurves use cubic polynomials, so have the formf(t) = a3 t^3 + a2 t^2 + a1 t + a0,where the coefficients are points in 3D. Blossoming converts thispolynomial to a more helpful form.  Let s = 1-t, and think oftri-linear interpolation:F([s0,t0],[s1,t1],[s2,t2]) =s0(s1(s2 c30 + t2 c21) + t1(s2 c21 + t2 c12)) +t0(s1(s2 c21 + t2 c12) + t1(s2 c12 + t2 c03))=c30(s0 s1 s2) +c21(s0 s1 t2 + s0 t1 s2 + t0 s1 s2) +c12(s0 t1 t2 + t0 s1 t2 + t0 t1 s2) +c03(t0 t1 t2).The data points c30, c21, c12, and c03 have been used in such away as to make the resulting function give the same value if anytwo arguments, say [s0,t0] and [s2,t2], are swapped. When [1-t,t]is used for all three arguments, the result is the cubic Beziercurve with those data points controlling it:f(t) = F([1-t,t],[1-t,t],[1-t,t])=  (1-t)^3 c30 +3(1-t)^2 t c21 +3(1-t) t^2 c12 +t^3 c03.Notice thatF([1,0],[1,0],[1,0]) = c30,F([1,0],[1,0],[0,1]) = c21,F([1,0],[0,1],[0,1]) = c12, _F([0,1],[0,1],[0,1]) = c03.In other words, cij is obtained by giving F argument t's i ofwhich are 0 and j of which are 1. Since F is a linear polynomialin each argument, we can find f(t) using a series of simple steps.Begin withf000 = c30, f001 = c21, f011 = c12, f111 = c03.Then compute in succession:f00t = s f000 + t f001, f01t = s f001 + t f011,f11t = s f011 + t f111,f0tt = s f00t + t f01t, f1tt = s f01t + t f11t,fttt = s f0tt + t f1tt.This is the de Casteljau algorithm for computing f(t) = fttt.It also has split the curve for the intervals [0..t] and [t..1].The control points for the first interval are f000, f00t, f0tt,fttt, while those for the second interval are fttt, f1tt, f11t,f111.If you evaluate 3 F([1-t,t],[1-t,t],[-1,1]) you will get thederivate of f at t. Similarly, 2*3 F([1-t,t],[-1,1],[-1,1]) givesthe second derivative of f at t, and finally using 1*2*3F([-1,1],[-1,1],[-1,1]) gives the third derivative.This procedure is easily generalized to different degrees,triangular patches, and B-spline curves.----------------------------------------------------------------------Subject 4.03: How do I find a t value at a specific point on a Bezier?In general, you'll need to find t closest to your search point.There are a number of ways you can do this. Look at [Gems I, 607],there's a chapter on finding the nearest point on the Beziercurve.  In my experience, digitizing the Bezier curve is anacceptable method. You can also try recursively subdividing thecurve, see if you point is in the convex hull of the controlpoints, and checking is the control points are close enough to alinear line segment and find the nearest point on the linesegment, using linear interpolation and keeping track of thesubdivision level, you'll be able to find t.----------------------------------------------------------------------Subject 4.04: How do I fit a Bezier curve to a circle?Interestingly enough, Bezier curves can approximate a circle butnot perfectly fit a circle.A common approximation is to use four beziers to model a circle, eachwith control points a distance d=r*4*(sqrt(2)-1)/3 from the end points(where r is the circle radius), and in a direction tangent to thecircle at the end points.  This will ensure the mid-points of theBeziers are on the circle, and that the first derivative is continuous.The radial error in this approximation will be about 0.0273% of thecircle's radius.Michael Goldapp, "Approximation of circular arcs by cubicpolynomials" Computer Aided Geometric Design (#8 1991 pp.227-238)Tor Dokken and Morten Daehlen, "Good Approximations of circles bycurvature-continuous Bezier curves" Computer Aided GeometricDesign (#7 1990 pp. 33-41).See also http://www.whizkidtech.net/bezier/circle/ .----------------------------------------------------------------------Section 5. 3D computations----------------------------------------------------------------------Subject 5.01: How do I rotate a 3D point?Let's assume you want to rotate vectors around the origin of yourcoordinate system. (If you want to rotate around some other point,subtract its coordinates from the point you are rotating, do therotation, and then add back what you subtracted.) In 3D, you neednot only an angle, but also an axis. (In higher dimensions it getsmuch worse, very quickly.)  Actually, you need 3 independentnumbers, and these come in a variety of flavors.  The flavor Irecommend 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 conceptof an axis is unique to 3D. It is a line through the origincontaining all the points which do not move during the rotation.So we know if we are turning forwards or back, we use a vectorpointing out along the line. Suppose you want to use unit vector uas 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 quaternionmultiplication, q v q^-1, or just convert the quaternion to a 3x3matrix M. If the components of q are {(x,y,z),w], then you wantthe matrixM = {{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 basiccomputer graphics texts.  Quaternions are covered briefly in[Foley], and more extensively in several Graphics Gems, and theSIGGRAPH 85 proceedings.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);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);}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);}VECTOR Normalize(VECTOR v){VECTOR u;FPOINT norm = MUL(v.x,v.x)+MUL(v.y,v.y)+MUL(v.z,v.z);FPOINT scl = RECIPROCAL(SQRT(norm));u.x = MUL(v.x,scl); u.y = MUL(v.y,scl); u.z = MUL(v.z,scl);return (u);}QUAT MatrixFromPointsAngle(POINT o, POINT p, ANGLE theta, MATRIX m){QUAT q;VECTOR diff, axis;diff.x = p.x-o.x; diff.y = p.y-o.y; diff.z = p.z-o.z;axis = Normalize(diff);return (MatrixFromAnyAxisAngle(o,axis,theta,m));}----------------------------------------------------------------------Subject 5.02: What is ARCBALL and where is the source?Arcball is a general purpose 3D rotation controller described byKen Shoemake in the Graphics Interface '92 Proceedings.  Itfeatures good behavior, easy implementation, cheap execution, andoptional axis constraints. A Macintosh demo and electronic versionof the original paper (Microsoft Word format) may be obtained fromftp::/ftp.cis.upenn.edu/pub/graphics.Complete source code appears in Graphics Gems IV pp. 175-192. Afairly complete sketch of the code appeared in the originalarticle, in Graphics Interface 92 Proceedings, available fromMorgan Kaufmann.The original arcball code was written for IRIS GL. A translationinto OpenGL/GLUT, and for IRIS Performer, may be found at:http://cs.anu.edu.au/people/Hugh.Fisher/3dstuff/----------------------------------------------------------------------Subject 5.03: How do I clip a polygon against a rectangle?This is the Sutherland-Hodgman algorithm, an old favorite thatshould be covered in any textbook. See the selected list below.According to Vatti (q.v.)  "This[algorithm] produces degenerate edges in certain concave / selfintersecting polygons that need to be removed as a specialextension to the main algorithm" (though this is not a problem ifall you are doing with the results is scan converting them.)It should be noted that the Sutherland-Hodgman algorithmmay 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!)See alsohttp://www.csclub.uwaterloo.ca/u/mpslager/articles/sutherland/wr.html----------------------------------------------------------------------|Subject 5.04: How do I clip a polygon against another polygon?Klamer Schutte, klamer@ph.tn.tudelft.nlhas developed and implementedsome code in C++ to perform clipping of two possibly concave 2Dpolygons. 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.htmlTo 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 extendsthe above to permit holes.Alan Murta released a polygon clipper library (in C) which uses amodified version of the Vatti algorithm:http://www.cs.man.ac.uk/aig/staff/alan/software/index.htmlReferences:Weiler, K. "Polygon Comparison Using a Graph Representation", SIGGRAPH 80pg. 10-18Vatti, Bala R. "A Generic Solution to Polygon Clipping",Communications of the ACM, July 1992, Vol 35, No. 7, pg. 57-63----------------------------------------------------------------------Subject 5.05: How do I find the intersection of a line and a plane?If the plane is defined as:a*x + b*y + c*z + d = 0and the line is defined as:x = x1 + (x2 - x1)*t = x1 + i*ty = y1 + (y2 - y1)*t = y1 + j*tz = z1 + (z2 - z1)*t = z1 + k*tThen just substitute these into the plane equation. You end upwith: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 planeif the numerator is also zero (the point at t=0 satisfies theplane equation), otherwise the line is parallel to the plane.----------------------------------------------------------------------Subject 5.06: How do I determine the intersection between a ray and a triangle?First find the intersection between the ray and the plane in whichthe triangle is situated. Then see if the point of intersection isinside the triangle.Details may be found in [O'Rourke (C)] pp.226-238, whose code is athttp://cs.smith.edu/~orourke/ .Efficient code complete with statistical tests is described in the Mo:ller-Trumbore paper in J. Graphics Tools (C code downloadable from there):http://www.acm.org/jgt/papers/MollerTrumbore97/See also the full paper:http://www.Graphics.Cornell.EDU/pubs/1997/MT97.htmlSee also the "3D Object Intersection" page, described in Subject 0.05.----------------------------------------------------------------------Subject 5.07: How do I determine the intersection between a ray and a sphereGiven a ray defined as:x = x1 + (x2 - x1)*t = x1 + i*ty = y1 + (y2 - y1)*t = y1 + j*tz = z1 + (z2 - z1)*t = z1 + k*tand a sphere defined as:(x - l)**2 + (y - m)**2 + (z - n)**2 = r**2Substituting in gives the quadratic equation:a*t**2 + b*t + c = 0where:a = i**2 + j**2 + k**2b = 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 doesnot intersect the sphere. If it is zero, the line is tangential tothe sphere and if it is greater than zero it intersects at twopoints. Solving the equation and substituting the values of t intothe ray equation will give you the points.Reference:[Glassner:RayTracing]See also the "3D Object Intersection" page, described in Subject 0.05.----------------------------------------------------------------------Subject 5.08: How do I find the intersection of a ray and a Bezier surface?Joy I.K. and Bhetanabhotla M.N., "Ray tracing parametric surfacesutilizing numeric techniques and ray coherence", ComputerGraphics, 16, 1986, 279-286Joy and Bhetanabhotla is only one approach of three major methodclasses: tessellation, direct computation, and a hybrid of thetwo.  Tessellation is relatively easy (you intersect the polygonsmaking up the tessellation) and has no numerical problems, but canbe inaccurate; direct is cheap on memory, but very expensivecomputationally, and can (and usually does) suffer from precisionproblems 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.----------------------------------------------------------------------Subject 5.09: How do I ray trace caustics?See the work of Henrik Wann Jensen athttp://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 forSpherical 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 NaotaInamato and Lining Sun ",title          = "Gemstone Fire: Adaptive Dispersive Ray Tracingof Polyhedrons",year           = "1988",month          = "November",journal        = "The Visual Computer",volume         = "4",number         = "5",pages          = "259--70",keywords       = "caustics",}----------------------------------------------------------------------Subject 5.10: What is the marching cubes algorithm?The marching cubes algorithm is used in volume rendering toconstruct an isosurface from a 3D field of values.The 2D analog would be to take an image, and for each pixel, setit to black if the value is below some threshold, and set it towhite if it's above the threshold.  Then smooth the jagged blackoutlines by skinning them with lines.The marching cubes algorithm tests the corner of each cube (orvoxel) in the scalar field as being either above or below a giventhreshold.  This yields a collection of boxes with classifiedcorners.  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 thatmeets the classification of the cube.  For example, the followingare some 2D examples showing the cubes and their associatedsurface.- ----- +       - ----- -       - ----- +       - ----- +|:::'   |       |:::::::|       |::::   |       |   ':::||:'     |       |:::::::|       |::::   |       |.    ':||       |       |       |       |::::   |       |::.    |+ ----- +       + ----- +       - ----- +       + ----- -The result of the marching cubes algorithm is a smooth surfacethat approximates the isosurface that is constant along a giventhreshold. This is useful for displaying a volume of oil in ageological 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,seehttp://www.unchainedgeometry.com/jbloom/papers/index.htmlunder "Implicit Surface Polygonization."----------------------------------------------------------------------Subject 5.11: What is the status of the patent on the "marching cubes" algorithm?United States Patent Number: 4,710,876Date of Patent: Dec. 1, 1987Inventors: Harvey E. Cline, William E. LorensenAssignee: General Electric CompanyTitle: "System and Method for the Display of Surface Structures ContainedWithin the Interior Region of a Solid Body"Filed: Jun. 5, 1985http://www.delphion.com/Type in "4710876" (w/o commas, w/o quotes) into their search engine.United States Patent Number: 4,885,688Date of Patent: Dec. 5, 1989Inventor: Carl R. CrawfordAssignee: General Electric CompanyTitle: "Minimization of Directed Points Generated in Three-DimensionalDividing Cubes Method"Filed: Nov. 25, 1987Access as above.For alternative, unpatented algorithms, cf. Subj. 5.10.----------------------------------------------------------------------Subject 5.12: How do I do a hidden surface test (backface culling) with 3D points?Just define all points of all polygons in clockwise order.c = (x3*((z1*y2)-(y1*z2))+(y3*((x1*z2)-(z1*x2))+(z3*((y1*x2)-(x1*y2))+x1,y1,z1, x2,y2,z2, x3,y3,z3 = the first three points of thepolygon.If c is positive the polygon is visible. If c is negative thepolygon is invisible (or the other way).----------------------------------------------------------------------Subject 5.13: Where can I find algorithms for 3D collision detection?Check out "proxima", from Purdue, which is a C++ library for 3Dcollision 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 systemhttp://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-dimensionalobjects undergoing rigid motion and deformation. SOLID is designed tobe used in interactive 3D graphics applications, and is especiallysuited for collision detection of objects and worlds described in VRML.Written in standard C++, compiles under GNU g++ version 2.8.1 andVisual C++ 5.0."  See:http://www.win.tue.nl/cs/tt/gino/solid/SWIFT++: a C++ library for collision detection, exact and approximatedistance computation, and contact determination of three-dimensionalpolyhedral objects undergoing rigid motion.Some preliminary results indicate that it is faster than I-COLLIDE andV-CLIP, and more robust than I-COLLIDE.http://www.cs.unc.edu/~geom/SWIFT++ColDet: C++ library for 3D collison detection.  Works on genericpolyhedra, and even polygon soups.   Uses bounding box hierarchiesand 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----------------------------------------------------------------------Subject 5.14: How do I perform basic viewing in 3D?We describe the shape and position of objects using numbers,usually (x,y,z) coordinates. For example, the corners of a cubemight 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 withthese numbers requires mathematical study, but I will try to keepthis simple. At the least, we must understand that we havedesignated some point in space as the origin--coordinates(0,0,0)--and marked off lines in 3 mutually perpendiculardirections 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, 1meter, or 1 parsec; the numbers alone do not tell us.A picture on a screen is two steps removed from the 3D world itdepicts. First, it is a 2D projection; and second, it is a finiteresolution approximation to the infinitely precise projection. Iwill ignore the approximation (sampling) for now. To know what theprojection looks like, we need to know where our viewpoint is, andwhere the plane of the projection is, both in the 3D world. Thinkof it as looking out a window into a scene. As artists discoveredsome 500 years ago, each point in the world appears to be at apoint on the window. If you move your head or look out a differentwindow, everything changes. When the mathematicians understoodwhat the artists were doing, they invented perspective geometry.If your viewpoint is at the origin--(0,0,0)--and the window sitsparallel 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. Distantobjects will have large z values, causing them to shrink in thepicture. That's perspective.The trick is to take an arbitrary viewpoint and plane, andtransform the world so we have the simple viewing situation.There are two steps: move the viewpoint to the origin, then movethe 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 theviewplane, then scale so it is 1 unit away.After all this, we may find ourselves looking out upside- down. Itis traditional to specify some direction in the world or viewplaneas "up", and rotate so the positive y axis points that way (asnearly as possible if it's a world vector). Finally, we have actedso far as if the window was the entire plane instead of a limitedportal. A final shift and scale transforms coordinates in theplane to coordinates on the screen, so that a rectangular regionof interest (our "window") in the plane fills a rectangular regionof the screen (our "canvas" if you like).Details of how to define and perform the rotation of the viewplanehave been left out, but see "How can I aim a camera in a specificdirection?" elsewhere in this FAQ. One simple way to designate aplane is with the point closest to the origin, call it D. Thena point P is on the plane if D.P = D.D; or using d = ||D|| andN = 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 ofthe 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!) Infact ordinarily a clipping box, the "viewing frustum", is usedto 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 4x4homogeneous 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."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 ElementProcedures", S. A. Cannan, S. N. Muthukrishnan and R. P. Phillips,Engineering With Computers, No. 12, p. 243-255, 1996."Progressive Meshes", Hoppe, SIGGRAPH 96,http://research.microsoft.com/~hoppe/siggraph96pm/paper.htmSeveral 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----------------------------------------------------------------------Subject 5.16: How can I perform volume rendering?Two principal methods can be used:- Ray casting or front-to-back, where the volume is behind theprojection plane. A ray is projected from each point in the projectionplane through the volume. The ray accumulates the properties of eachvoxel it passes through.- Object order or back-to-front, where the projection plane is behindthe volume. Each slice of the volume is projected on the projectionplane, from the farest plane to the nearest plane.You can also use the marching-cubes algorithm, if the extraction ofsurfaces 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 enclosesthe volume. From each pixel of the projection plane, cast a raythrough the volume by using a Bresenham 3D algorithm.The ray accumulates properties from each voxel intersected, stoppingwhen the ray exits the volume. The pixel value onthe 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 transparencyC' the new ray color           t' the new ray transparencyCv the voxel color             tv the voxel transparencythen:C' = C + t*Cv       and        t' = t * tvwith initial values: C = 0.0 and t = 1.0An alternate version: instead of C' = C + t*Cv , use :C' = C + t*Cv*(1-tv)^p     with p a float variable.Sometimes this gives the best results.When the ray has accumulated transparency, if it becomes negligible(e.g., t<0.1), the process can stop and proceed to the next ray.References:Bresenham 3D:- http://www.sct.gu.edu.au/~anthony/info/graphics/bresenham.procs- [Gems IV] p. 366Volume rendering:- [Watt:Animation] pp. 297-321- IEEE Computer Graphics and applicationVol. 10, Nb. 2, March 1990 - pp. 24-53- "Volume Visualization"Arie Kaufman - Ed. IEEE Computer Society Press Tutorial- "Efficient Ray Tracing of Volume Data"Marc Levoy - ACM Transactions on Graphics, Vol. 9, Nb 3, July 1990----------------------------------------------------------------------Subject 5.17: Where can I get the spline description of the famous teapot etc.?See the Standard Procedural Databases software, whose officialdistribution site ishttp://www.acm.org/tog/resources/SPD/This database contains much useful 3D code, including spline surfacetessellation, for the teapot.----------------------------------------------------------------------Subject 5.18: How can the distance between two lines in space be computed?Let x_i be points on the respective lines and n_i unit directionvectors 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 thethese points on two 3D lines.  It also classifiesthe lines as parallel, intersecting, or (the generic case) skew.What's listed below shows the main ideas; the full code is athttp://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)intline_line_closest_points3d (POINT *pA, POINT *pB,                   // computed pointsconst POINT *a, const VECTOR *adir,             // line A, point-normal formconst 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 bothvcross ( cdir, adir, bdir );// check for near-parallel linesif ( !vnorm( cdir ) )   {*pA = *a;       // all points are closest*pB = *b;return 0;       // degenerate: lines parallel}// form plane containing line A, parallel to cdirplane_from_two_vectors_and_point ( ac, cdir, adir, a );// form plane containing line B, parallel to cdirplane_from_two_vectors_and_point ( bc, cdir, bdir, b );// closest point on A is line A ^ bcintersect_line_plane ( pA, a, adir, bc );// closest point on B is line B ^ acintersect_line_plane ( pB, b, bdir, ac );// distinguish intersecting, skew linesif ( edist( pA, pB ) < 1.0E-5F )return 1;     // coincident: lines intersectelsereturn 2;     // distinct: lines skew}Also Dave Eberly has code for computing distance between variousgeometric primitives, including MinLineLine(), athttp://www.magic-software.com----------------------------------------------------------------------Subject 5.19: How can I compute the volume of a polyhedron?Assume that the surface is closed, every face is a triangle, andthe vertices of each triangle are oriented ccw from the outside.Let Volume(t,p) be the signed volume of the tetrahedron formedby a point p and a triangle t.  This may be computed by a 4x4determinant, as in [O'Rourke (C), p.26].Choose an arbitrary point p (e.g., the origin), and computethe sum Volume(t_i,p) for every triangle t_i of the surface.  Thatis the volume of the object.  The justification for this claimis nontrivial, but is essentially the same as the justification forthe 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 variousalgorithms available athttp://www.Mathpool.Uni-Augsburg.DE/~enge/Volumen.html .----------------------------------------------------------------------Subject 5.20:  How can I decompose a polyhedron into convex pieces?Usually this problem is interpreted as seeking a collectionof pairwise disjoint convex polyhedra whose set union is theoriginal polyhedron P.  The following facts are known aboutthis difficult problem:o  Not every polyhedron may be partitioned by diagonals intotetrahedra.  A counterexample is due to Scho:nhardt[O'Rourke (A), p.254].o  Determining whether a polyhedron may be so partitioned isNP-hard, a result due to Seidel & Ruppert [Disc. Comput. Geom.7(3) 227-254 (1992).]o  Removing the restriction to diagonals, i.e., permittingso-called Steiner points, there are polyhedra of n verticesthat 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 mostO(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, atftp://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, withax, 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 thecross product of a and b.  Then the center m and radius r of thecircumsphere 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 betweencoordinates.  The advantage is that the relative error incurred in thecomputation of r is also a function of the _differences_ between thevertices, and is not influenced by the _absolute_ coordinates of thevertices.  In most applications, vertices are usually nearer to each otherthan to the origin, so this property is advantageous.Similarly, the formula for m incurs roundoff error proportional to thedifferences between vertices, but not proportional to the absolutecoordinates of the vertices until the final addition.- These expressions are unstable in only one case:  if the denominator isclose to zero.  This instability, which arises if the tetrahedron isnearly degenerate, is unavoidable.  Depending on your application, youmay want to use exact arithmetic to compute the value of the determinant.Seehttp://www.geom.umn.edu/software/cglist/alg.htmlorhttp://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 sideof the plane containing T2, or T2 lies strictly to one side of theplane containing T1, the triangles do not intersect.  Otherwise,compute the line of intersection L between the planes.  Let Ikbe 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-triangleintersection test," J. Graphics Tools 2(2) 25-30 1997.  C code athttp://www.acm.org/jgt/papers/Moller97/.  See alsohttp://www.ce.chalmers.se/staff/tomasm/code/http://www.magic-software.com/MgcIntersection.htmlSee also the "3D Object Intersection" page, described in Subject 0.05.NASA's "Intersect" code will intersect any number of triangulatedsurfaces provided that each of the surfaces is both closed and manifold.http://www.nas.nasa.gov/~aftosmis/cart3d/surfaceModeling.html#AuxProgsBased on "Robust and Efficient Cartesian Mesh Generation forComponent-Based Geometry" AIAA Paper 97-0196. 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 throughthe 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/tog/1992-11-3/p228-meyers/Code (NUAGES) is available athttp://www-sop.inria.fr/prisme/logiciel/nuages.html.enftp://ftp-sop.inria.fr/prisme/NUAGES/Nuages/NUAGES_SRC.tar.gzFor (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 athttp://members.tripod.com/~GeomWiz/www.sites.htmlNew code, Cocone, written with CGAL, based on recent work byN. Amenta, S. Choi, T. K. Dey and N. Leekha:http://www.cis.ohio-state.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 pointsin any dimension can be found by Linear Programming.  This was provedby Emo Welzl in the paper, "Smallest enclosing disks (balls andellipsoids)" [Lecture Notes Comput. Sci., 555, Springer-Verlag, 1991,359--370].  +   Code developed by Bernd Gaertner available (GNUGeneral Public License) at:http://www.inf.ethz.ch/~gaertner/miniball.htmlThis code is remarkably efficient: e.g., 2 seconds for 10^6 points in3D on an Ultra-Sparc II.  See also Dave Eberly's direct implementationof Welzl's algorithm:http://www.magic-software.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 aretheir virtues?"The heat of debate is hard to explain, but it's been happening for manydecades. When Gibbs first deprecated the quaternion product and split itinto a cross product and a dot product, one outraged Victorian calledthe result a "hermaphrodite monster" -- and that before the Internet'sflame wars. Generally, the quaternion advocates seem to feel theopponents are lazy or thick-headed, and that deeper understanding ofquaternions would lead to deeper appreciation. The opponents don'tappreciate that attitude, and seem to feel the advocates are snooty orsheep, and that matrices and such are less abstract and do just fine.(Advocates of Clifford algebra would claim that both sides are mired inthe past.) Passion aside, quaternions have appropriate uses, as do theiralternatives.Someone new to the debate first needs to know what quaternions are, andwhat they're supposed to be good for. Quaternions are a quadruple ofnumbers, 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 squaresof the four components. Some writers confuse this with the magnitude,which is the square root of the norm. Another common misconception isthat only quaternions of unit norm can be used, those with the sum ofthe four squares equal to 1, but that is wrong (though they arepreferred).[U sin a,cos a] rotates by angle 2a around unit vector UPopular non-quaternion options are 3x3 special orthogonal matrices (9numbers with constraints), Euler angles (3 numbers), axis-angle (4numbers), and angular velocity vectors (3 numbers). None of theseoptions actually _are_ rotations, which are physical; they _represent_rotations. The distinction is important because, for example, it iscommon to use an axis-angle with an angle greater than 360 degrees totell an animation system to spin an object more than a full turn,something a matrix cannot say. In mathematics, the usual meaning of arotation would not allow the multiple spin version, which can lead toconfusing debates.q and -q represent the same rotationTwo rotations, the physical things, can be applied one after the other.Assuming the two rotation axes have a least one point in common, theresult will be another rotation. Some rotation representations handlethis gracefully, some don't. For quaternions and matrices, forms ofmultiplication are defined such that the product gives the desiredresult. 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 thetwo 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 quaternionsdo a particularly good job of representing similar rotations withsimilar numbers.similar(q1,q2) = |q1.q2| = | x1x2+y1y2+z1z2+w1w2 |Points in space, the physical things, are normally represented as 3 or 4numbers. The effect of a rotation on a collection of points can becomputed from the representation of the rotation, and here matrices seemfastest, using three dot products. Using their own product twice,quaternions are a bit less efficient. (They are usually converted tomatrices at the last minute.)p2 = q p1 q^(-1)Sequences of rotations can be interpolated, so that the object beingturned is rotated to specific poses at specific times. This motivatedKen Shoemake's early use of quaternions in computer graphics, aspublished in 1985. He used an analog of linear interpolation (sometimescalled "lerp") that he called "Slerp", and also introduced an analog ofa piecewise Bezier curve. A few years later in some course notes hedescribed another curve variation he called "Squad", which still seemsto be popular. Later authors have proposed many alternatives.sin (1-t)A      sin tASlerp(q1,q2;t) = q1 ---------- + q2 ------, cos A = q1.q2sin A         sin ASquad(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 ofcomputations which also depend on rotation representation. Constrainedrotations like a wheel on an axle or the elbow bend of a robot typicallyuse specialized representations, such as an angle alone. In many generalsituations, however, quaternions have proved valuable.2 dq = W q dt, W is the angular velocity vectorUser interfaces for 3D rotation also require a representation. Directmanipulation interfaces typically use angles for jointed figures, butfor freer manipulation may use quaternions, as in Arcball orthrough-the-lens camera control. As Shoemake's _full_ Graphics Gems codefor Arcball demonstrates (with the [CAPS LOCK] key), any rotation can begraphed as an arc on a sphere. (Not to be confused with the quaternionunit sphere in 4D.) Whether quaternions, or any other representation,are helpful for numeric presentation and input seems a matter of tasteand circumstance.q = U2 U1^(-1) = [U1xU2,U1.U2]Because of their metric properties for representing rotations, unitquaternions are most common. Advocates frequently point out that it isfar cheaper to normalize the length of a non-zero quaternion than tobring a matrix back to rotation form. Also Shoemake's later conversioncode cheaply creates a correct rotation matrix from _any_ quaternion(found with his Euler angle code from Graphics Gems, which does the samefor all 24 variations of that representation).Normalize(q) = q/Sqrt(q.q)Comparisons to Euler angles may mention "gimbal lock" (frequentlymisspelled) as a disadvantage quaternions avoid. In the physical worldwhere gyroscopes are mounted on nested pivots, which are called gimbals,locking is a real problem quaternions cannot help. What's usually meantis that because the similarity of rotations organizes them somewhat likea sphere, while similarity of vectors is quite different, an inevitablemisfit plagues Euler angles. This can show up in behavior much likephysical gimbal lock, but also in other ways. The difficulties aretopological, and aiming runs into them as well, even if quaternions areused. Quaternion authors who propose using curves in the vector space ofquaternion logarithms often risk the misfit unawares. Frankly, you mustpick through the literature carefully, whether informal and online orrefereed and printed, because mistakes are tragically common.To explore Graphics Gems code, see "Where is all the source?" in thisFAQ. To read more about quaternions, you have many options. Since theywere discovered in 1843 by Hamilton (the same Irish mathematician andphysicist whose name shows up in quantum mechanics), quaternions havefound many friends, as a web search will reveal. Quaternions can beapproached and applied in numerous different ways, so if you keeplooking it's likely you will find something that suits your taste andyour 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/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.pdfVisualizing 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.htmlMathematical, but not too fast and not too fancy.http://www.cs.berkeley.edu/~laura/cs184/quat/quaternion.htmlLaura Downs covers the fundamentals.http://graphics.cs.ucdavis.edu/GraphicsNotes/Quaternions/Quaternions.htmKen Joy covers the fundamentals.http://www.gg.caltech.edu/STC/rr_sig97.htmlHigh-tech interpolation method. Demanding but rewarding.Duff, Tom: Quaternion Splines for Animating Orientation,Proceedings of the USENIX Association Second Computer GraphicsWorkshop (held in Monterey, CA 12-13 Dec. 1985), pp. 54-62.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 unitvector to line up with another. To aim at an object, you can subtractthe position of the camera from the position of the object to get avector which you then normalize. The vector you want to turn is thecamera forward vector, commonly a unit vector along the camera -z axis.Be warned that more than one rotation can achieve aim alone. (The issueis rather complicated, as laid out in Ken Shoemake's article on twistreduction in Graphics Gems IV.) For example, even if the camera isalready properly aimed you could rotate it around its z axis. The mostdirect rotation is given by the non-unit quaternionq = [(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). Inthat case just rotate 180 degrees around the y axis, usingq = [(0,1,0),0]If the camera is level after rotation by quaternion [(x,y,z),w], the ycomponent of its rotated x axis will be zero, soxy+wz = 0If it is upright, the y component of its rotated y axis will not benegative, soww-xx+yy-zz >= 0To ensure these two desirable properties, aim with a more sophisticatednon-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 hasnot been normalized and its magnitude is m = sqrt(aa+bb+cc), substitute1->m. That is, use t = r+m and use m+c.More generally, to rotate unit vector U1 directly to unit vector U2, thenon-unit quaternion will beq = [U1xU2,1+U1.U2]Why? If U is a unit vector, then it is normal to a plane through theorigin with equation U.P = 0. Reflection in that plane is given byreversing the U component of P.reflect(P,U) = P ^� 2(U.P)UThe quaternion product of U1 and U2 is [U1xU2,-U1.U2], so-2 (U.P) = PU + UPNoting UU = -1, this gives a quaternion reflection formula.reflect(P,U) = P + (PU+UP)U= P ^� P + UPU= UPUReflecting with U1 then U2, by U2(U1 P U1)U2, rotates by twice the anglebetween the planes, with axis perpendicular to both normals. Noting U1U2is the conjugate of U2U1, and -q rotates like q, the rotation quaternionisq = -U2U1 = -[U2xU1,-U2.U1] = [U1xU2,U1.U2]This q fails to aim U1 at U2 by rotating twice as much as needed, butits 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 degreerotation.[These quaternion methods were provided by Ken Shoemake.]----------------------------------------------------------------------Subject 5.27: How do I transform normals?In 3D, the orientation of a plane in space can be given by avector perpendicular to the plane, a "normal vector" or "normal"for short. Often it is convenient to keep that vector of unitlength, or "normalized"; be careful of the different meaningsof "normality". A smooth surface has a plane tangent to eachpoint, and by extension a normal to that plane is called a"surface normal". Graphics code also cheats by associatingartificial normal vectors with the vertices of polygonal modelsto simulate the reflection properties of curved surfaces;these are called "vertex normals".The "orientation" of a plane has two meanings, both of whichusually apply. Aside from the rotational tilting and turningmeaning, there is also the sense of "side". A closed convexsurface made of polygons has two sides, an inside and anoutside, and normals can be assigned to the polygons in sucha way that they all consistently point outside. This isoften desirable for shading and culling.When a model is defined in one coordinate system and used inanother, as is commonly done, it may be necessary to transformnormals also. If the change of point coordinates is effectedby means of a rotation plus a translation, one simply rotatesthe normals as well (with no translation). Other coordinatechanges need more care.Although it is possible to use projective transformations tomap a model into world coordinates, ordinarily they are usedonly for viewing. It is usually a mistake to apply perspectiveto a normal, as shading and culling are best done in worldcoordinates for correct results. Also perspective may becomputed 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 translationis called an affine transformation, and is performed as amatrix multiplication plus a vector addition:p' = A(p) = Lp + t.When the model-to-world point transform is affine, the properway to transform normals is with the transpose of the inverseof L.n' = (L^{-1})^T nHowever that is not enough.If L includes scaling effects, a unit normal in model spacewill usually transform to a non-unit normal, which can causeproblems for shaders and other code. This may need correctingafter the normal is transformed.If L includes reflection, the inside-outside orientation ofthe normal is reversed. This, too, can cause problems, andmay need correcting. The determinant of L will be negativein this case.When a complicated distortion is used, it must be approximateddifferently at each point in space by a linear transform madeup of partial derivatives, the Jacobian. The matrix for theJacobian 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 matrixproduct 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 benormal to it. Vectors at right angles have zero dot product.n^T v = 0The transform v' of v isv' = (Lq+t)-(Lp+t)= (Lq-Lp)+(t-t)= L(q-p)= LvThe transform n' of n will remain normal if it satisfiesn'^T v' = n^T vLet n' = Mn for some M. Then the requirement isn'^T v' = (Mn)^T (Lv)= n^T (M^T L) v= n^T vThis holds ifM^T L = Iwhere I is the identity. Right multiplying by the inverseL^{-1} and transposing both sides gives, as claimed,M = (L^{-1})^TWhen L is a rotation, L^{-1} = L^T, so M = L. When L has noinverse it will still have an "adjoint" to substitute forfor orthogonality purposes, differing only by a scale factor.----------------------------------------------------------------------Section 6. Geometric Structures and Mathematics----------------------------------------------------------------------Subject 6.01: Where can I get source for Voronoi/Delaunay triangulation?For 2-d Delaunay triangulation, try Shewchuk's triangle program.  Itincludes options for constrained triangulation and quality meshgeneration.  It uses exact arithmetic.The Delaunay triangulation is equivalent to computing the convex hullof the points lifted to a paraboloid.  For n-d Delaunay triangulationtry Clarkson's hull program (exact arithmetic) or Barber and Huhdanpaa'sQhull program (floating point arithmetic).  The hull program alsocomputes Voronoi volumes and alpha shapes.  The Qhull program alsocomputes 2-d Voronoi diagrams and n-d Voronoi vertices.  The output ofboth programs may be visualized with Geomview.There are many other codes for Delaunay triangulation and Voronoidiagrams.  See Amenta's list of computational geometry software.Especially noteworthy is the CGAL code: Subject 0.07 andhttp://www.cgal.org.The Delaunay triangulation satisfies the following property: thecircumcircle of each triangle is empty.  The Voronoi diagram is theclosest-point map, i.e., each Voronoi cell identifies the points thatare closest to an input site.  The Voronoi diagram is the dual ofthe Delaunay triangulation.  Both structures are defined for generaldimension.  Delaunay triangulation is an important part of meshgeneration.There is a FAQ of polyhedral computation explaining how to computen-d Delaunay triangulation and n-d Voronoi diagram using a convex hullcode, and how to use the linear programming technique todetermine the Voronoi cells adjacent to a given Voronoi cellefficiently for large scale or higher dimensional cases.Avis' lrs code uses the same file formats as cdd. Ituses exact arithmetic and is usefulfor problems with very large output size, since it does notrequire storing the output.On a closely related topic, seehttp://www.cis.ohio-state.edu/~tamaldey/medialaxis.htmfor computation of the 3D medial axis from the Voronoi diagram.References:Amenta:   http://www.geom.umn.edu/software/cglistAvis:     ftp://mutt.cs.mcgill.ca/pub/C/lrs.htmlBarber &  http://www.geom.umn.edu/locate/qhullHuhdanpaa ftp://geom.umn.edu/pub/software/ftp://ftp.geom.umn.edu/pub/software/Clarkson: http://cm.bell-labs.com/netlib/voronoi/hull.htmlftp://cm.bell-labs.com/netlib/voronoi/hull.shar.ZGeomview: http://www.geom.umn.edu/locate/geomviewftp://geom.umn.edu/pub/software/geomview/Polyhedral Computation FAQ:http://www.ifor.math.ethz.ch/ifor/staff/fukuda/fukuda.htmlShewchuk: http://www.cs.cmu.edu/~quake/triangle.htmlftp://cm.bell-labs.com/netlib/voronoi/triangle.shar.Z[O' Rourke (C)] pp. 168-204[Preparata & Shamos] pp. 204-225Chew, L. P. 1987. "Constrained Delaunay Triangulations," Proc. ThirdAnnual ACM Symposium on Computational Geometry.Chew, L. P. 1989. "Constrained Delaunay Triangulations," Algorithmica4:97-108. (UPDATED VERSION OF CHEW 1987.)Fang, T-P. and L. A. Piegl. 1994. "Algorithm for Constrained DelaunayTriangulation," The Visual Computer 10:255-265. (RECOMMENDED!)Frey, W. H. 1987. "Selective Refinement: A New Strategy for AutomaticNode Placement in Graded Triangular Meshes," International Journal forNumerical Methods in Engineering 24:2183-2200.Guibas, L. and J. Stolfi. 1985. "Primitives for the Manipulation ofGeneral Subdivisions and the Computation of Voronoi Diagrams," ACMTransactions on Graphics 4(2):74-123.Karasick, M., D. Lieber, and L. R. Nackman. 1991. "Efficient DelaunayTriangulation Using Rational Arithmetic," ACM Transactions on Graphics10(1):71-91.Lischinski, D. 1994. "Incremental Delaunay Triangulation," GraphicsGems 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 RadiosityApproach," Proc. Eurographics '89, W. Hansmann, F. R. A. Hopgood, andW. Strasser, Eds. Elsevier Science Publishers, 345-353.Subramanian, G., V. V. S. Raveendra, and M. G. Kamath. 1994. "RobustBoundary Triangulation and Delaunay Triangulation of ArbitraryTriangular Domains," International Journal for Numerical Methods inEngineering 37(10):1779-1789.Watson, D. F. and G. M. Philip. 1984. "Survey: SystematicTriangulations," Computer Vision, Graphics, and Image Processing26:217-223.----------------------------------------------------------------------Subject 6.02: Where do I get source for convex hull?For n-d 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 improvedhandling of difficult inputs.   Qhull computes convex hulls,Delaunay triangulations, Voronoi diagrams, and halfspaceintersections about a point.Qhull handles numeric precision problems by merging facets. The outputof both programs may be visualized with Geomview.In higher dimensions, the number of facets may be much smaller thanthe number of lower-dimensional features.  When this is the case,Fukuda's cdd program is much faster than Qhull or hull.There are many other codes for convex hulls.  See Amenta's list ofcomputational geometry software.References:Amenta:   http://www.geom.umn.edu/software/cglist/Barber &  http://www.geom.umn.edu/locate/qhullHuhdanpaaClarkson: http://cm.bell-labs.com/netlib/voronoi/hull.htmlftp://cm.bell-labs.com/netlib/voronoi/hull.shar.ZGeomview: http://www.geom.umn.edu/locate/geomviewftp://geom.umn.edu/pub/software/geomview/Fukuda:   http://www.ifor.math.ethz.ch/staff/fukuda/cdd_home/cdd.htmlftp://ftp.ifor.math.ethz.ch/pub/fukuda/cdd/[O' Rourke (C)] pp. 70-167C code for Graham's algorithm on p.80-96.Three-dimensional convex hull discussed in Chapter 4 (p.113-67).C code for the incremental algorithm on p.130-60.[Preparata & Shamos] pp. 95-184----------------------------------------------------------------------Subject 6.03: Where do I get source for halfspace intersection?For n-d halfspace intersection, try Fukuda's cdd program or Barberand Huhdanpaa's Qhull program.  Both use floating point arithmetic.Fukuda includes code for exact arithmetic.  Qhull handles numericprecision problems by merging facets.Qhull computes halfspace intersection by computing a convex hull.The intersection of halfspaces about the origin is equivalent to theconvex hull of the halfspace coefficients divided by their offsets.See Subject 6.02 for more information.References:Barber &  http://www.geom.umn.edu/locate/qhullHuhdanpaaFukuda:   ftp://ifor13.ethz.ch/pub/fukuda/cdd/[Preparata & Shamos] pp. 315-320----------------------------------------------------------------------Subject 6.04: What are barycentric coordinates?Let p1, p2, p3 be the three vertices (corners) of a closedtriangle T (in 2D or 3D or any dimension).  Let t1, t2, t3 be realnumbers 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 inthe plane of T and is inside T.  The numbers (t1,t2,t3) are called thebarycentric coordinates of p with respect to T. They uniquely identifyp.  They can be viewed as masses placed at the vertices whosecenter of gravity is p.For example, let p1=(0,0), p2=(1,0), p3=(3,2).  Thebarycentric coordinates (1/2,0,1/2) specify the pointp = (0,0)/2 + 0*(1,0) + (3,2)/2 = (3/2,1), the midpoint of the p1-p3edge.If p is joined to the three vertices, T is partitionedinto three triangles.  Call them T1, T2, T3, with Ti not incidentto pi.  The areas of these triangles Ti are proportional to thebarycentric 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 thethree vertices of your triangle.  Any point P inside can be expresseduniquely as P = aA + bB + cC, where a+b+c=1 and a,b,c are each >= 0.Knowing a and b permits you to calculate c=1-a-b.  So if you cangenerate two random numbers a and b, each in [0,1], such thattheir 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 arandom number.) If a+b>1, replace a by 1-a, b by 1-b.  Let c=1-a-b.Then aA + bB + cC is uniformly distributed in triangle ABC: the reflectionstep a=1-a; b=1-b gives a point (a,b) uniformly distributed in the triangle(0,0)(1,0)(0,1), which is then mapped affinely to ABC.Now you have barycentric coordinates a,b,c.  Compute your pointP = aA + bB + cC.Reference: [Gems I], Turk, pp. 24-28, contains a similar but differentmethod 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 asense in which only the five Platonic solids achieve regulartesselations, as they are the only ones whose faces are regularand equal, with each vertex incident to the same number of faces.But generally "even distribution" focusses not so much on theinduced tesselation, as it does on the distances and arrangementof the points/vertices.  For example, eight points can be placedon the sphere further away from one another than is achieved bythe vertices of an inscribed cube: start with an inscribed cube,twist the top face 45 degrees about the north pole, and thenmove the top and bottom points along the surface towards the equatora bit.The various definitions of "evenly distributed" lead into moderatelycomplex mathematics. This topic happens to be a FAQ on sci.mathas wellas on comp.graphics.algorithms; a much more extensive and rigorousdiscussion written by Dave Rusin can be found at:http://www.math.niu.edu/~rusin/known-math/95/sphere.faqA simple method of tesselating the sphere is repeated subdivision. Anexample starts with a unit octahedron. At each stage, every triangle inthe tesselation is divided into 4 smaller triangles by creating 3 newvertices 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 atftp://ftp.cs.unc.edu/pub/users/jon/sphere.cOne frequently used definition of "evenness" is a distribution thatminimizes a 1/r potential energy function of all the points, so that ina global sense points are as "far away" from each other as possible.Starting from an arbitrary distribution, we can either use numericalminimization algorithms or point repulsion, in which all the pointsare considered to repel each other according to a 1/r^2 force law anddynamics are simulated. The algorithm is run until some convergencecriterion is satisfied, and the resulting distribution is our approximation.Jon Leech developed code to do just this.  It is available atftp://ftp.cs.unc.edu/pub/users/jon/points.tar.gz
.See his README file for more information.  His distribution includessample output files for various n < 300, which may be used off-the-shelfif that is all you need.Another method that is simpler than the above, but attains lessuniformity, is based on spacing the points along a spiral thatencircles 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 athttp://cm.bell-labs.com/netlib/polyhedra/index.html
, orhttp://netlib.bell-labs.com/netlib/polyhedra/index.html
, orhttp://www.netlib.org/polyhedra/index.htmlFor 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 isthe "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 spherevolume ratio is pi/6, the average number of iterations before anacceptable vector is found is 6/pi = 1.90986.  This essentiallydoubles the effort, and makes this method slower than the "trigmethod."  A timing comparison conducted by Ken Sloan showed thatthe 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 aretwice as slow as the rejection method.The trig method.  This method works only in 3-space, but it isvery fast.  It depends on the slightly counterintuitive fact (seeproof below) that each of the three coordinates is uniformlydistributed 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.  Thisconstrains the chosen point to lie on a circle parallel to theX-Y plane, and the obvious trig method may be used to obtain theremaining coordinates.(a) Choose z uniformly distributed in [-1,1].(b) Choose t uniformly distributed on [0, 2*pi).(c) Let r = sqrt(1-z^2).(d) Let x = r * cos(t).(e) Let y = r * sin(t).This method uses uniform deviates (faster to generate than normaldeviates), and no set of coordinates is ever rejected.Here is a proof of correctness for the fact that the z-coordinate isuniformly distributed.  The proof also works for the x- and y-coordinates, but note that this works only in 3-space.The area of a surface of revolution in 3-space is given byA = 2 * pi * int_a^b f(x) * sqrt(1 + [f'(x}]^2) dxConsider a zone of a sphere of radius R.  Since we are integrating inthe z direction, we havef(z) = sqrt(R^2 - z^2)f'(z) = -z / sqrt(R^2-z^2)1 + [f'(z)]^2 = r^2 / (R^2-z^2)A = 2 * pi * int_a^b sqrt(R^2-z^2) * R/sqrt(R^2-z^2) dz= 2 * pi * R int_a^b dz= 2 * pi * R * (b-a)= 2 * pi * R * h,where h = b-a is the vertical height of the zone.  Notice how the integrandreduces 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.5lf/t.5lf/t.5lf/n", i, x,y,z);
X[i] = x; Y[i] = y; Z[i] = z;
}
}A complete package is available atftp://cs.smith.edu/pub/code/sphere.tar.gz(4K),reachable from http://cs.smith.edu/~orourke/ .If one wants to generate the random points in terms of longitudeand latitude in degrees, these equations suffice:Longitude = random * 360 - 180Latitude  = [arccos( random * 2 - 1 )/pi ] * 180 - 90References:[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 oftwo distinct points, for a total of six numbers. Or, they can be given asthe coordinates of two distinct planes, eight numbers. What's wrong withthese? Nothing; but we can do better. Pluecker coordinates are, in a sense,halfway between these extremes, and can trivially generate either. Neitherextreme is as efficient as Pluecker coordinates for computations.When Pluecker coordinates generalize to Grassmann coordinates, as laidout beautifully in [Hodge], Chapter VII, the determinant definitionis clearly the one to use. But 3D lines can use a simpler definition.Take two distinct points on a line, sayP = (Px,Py,Pz)
Q = (Qx,Qy,Qz)Think of these as vectors from the origin, if you like. The Plueckercoordinates for the line are essentiallyU = P - Q
V = P x QExcept for a scale factor, which we ignore, U and V do not depend on thespecific P and Q! Cross products are perpendicular to their factors, so wealways have U.V = 0. In [Stolfi] lines have orientation, so are the sameonly 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 wDefine Pluecker coordinate Gij as the determinant of rows i and j, in thatorder. Notice that Giw = Pi - Qi, which is Ui. Now let (i,j,k) be a cyclicpermutation of (x,y,z), namely (x,y,z) or (y,z,x) or (z,x,y), and noticethat Gij = Vk. Determinants are anti-symmetric in the rows, so Gij = -Gji.Thus all possible Pluecker line coordinates are either zero (if i=j) orcomponents 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 apoint P and vector U as on two points. We can also begin with a 2x4 matrixwhose 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 evenpermutation 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 isnon-zero, say Ui, which is Giw. Create homogeneous points Pj = Gjw + Gij,and Qj = Gij. (Don't expect the P and Q that we started with, and don'texpect 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) andV = (-4,6,-2). The direct determinant forms are Gxw=0, Gyw=1, Gzw=3,Gyz=-4, Gzx=6, Gxy=-2, and the dual forms are G'yz=0, G'zx=1, G'xy=3,G'xw=-4, G'yw=6, G'zw=-2. Take Uz = Gzw = G'xy = 3 as a suitable non-zeroelement. 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 = 0That is, a point P is on the line if it satisfies both these equations:3 Px + 0 Py + 0 Pz - 6 Pw = 00 Px + 3 Py - 1 Pz - 4 Pw = 0We 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 formreveals even permutations of (x,y,z,w):G1xw G2yz + G1yw G2zx + G1zw G2xy + G1yz G2xw + G1zx p2yw + G1xy p2zwTwo 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 mightlook 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 testwhich 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 intersectsa triangle in 3-space. One way to do this is to represent the ray andthe edges of the triangle with Pluecker coordinates. The ray hits thetriangle 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 allthree 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 onlya 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 thecase of lines in 2D, only three coordinates are required:

 Gxw = Px-Qx     = -G'y

 Gyw = Py-Qy     =  G'x

 Gxy = PxQy-PyQx =  G'w

 Since P and Q are distinct, Giw is non-zero 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 anideal point (meaning they're parallel in ordinary 2D). Callingtheir homogeneous point of intersection P, the computation fromGrassmann coordinates nicely illustrates the convenience. (SeeSubj 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] fora 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, possibletopics, corrections, or pointers to information.

转自：http://www.faqs.org/faqs/graphics/algorithms-faq/