home *** CD-ROM | disk | FTP | other *** search
/ ftp.pasteur.org/FAQ/ / ftp-pasteur-org-FAQ.zip / FAQ / graphics / algorithms-faq next >
Encoding:
Internet Message Format  |  2003-02-16  |  153.7 KB

  1. From: orourke@cs.smith.edu (Joseph O'Rourke)
  2. Newsgroups: comp.graphics.algorithms,news.answers,comp.answers
  3. Subject: comp.graphics.algorithms Frequently Asked Questions
  4. Supersedes: <graphics/algorithms-faq-1-1045285201@cs.smith.edu>
  5. Followup-To: comp.graphics.algorithms
  6. Date: Sat, 15 Feb 2003 16:16:44 +0000 (UTC)
  7. Lines: 3587
  8. Sender: orourke@grendel.csc.smith.edu
  9. Approved: news-answers-request@MIT.EDU
  10. Distribution: world
  11. Expires: 03/02/03 11:16:44 EDT
  12. Message-ID: <graphics/algorithms-faq-1-1045325804@cs.smith.edu>
  13. Reply-To: orourke@cs.smith.edu (Joseph O'Rourke)
  14. Keywords: faq computer graphics algorithms geometry
  15. X-Content-Currency: This FAQ changes regularly. See item (0.03) for instructions
  16.    on how to obtain a new copy via FTP.
  17. NNTP-Posting-Host: grendel.csc.smith.edu
  18. Organization: Smith College, Northampton Mass, USA
  19. Path: senator-bedfellow.mit.edu!bloom-beacon.mit.edu!aanews.merit.edu!elk.ncren.net!news.umass.edu!news.smith.edu!grendel.csc.smith.edu!not-for-mail
  20. Xref: senator-bedfellow.mit.edu comp.graphics.algorithms:134327 news.answers:246121 comp.answers:52786
  21.  
  22. Posted-By: auto-faq 3.3.1 beta (Perl 5.006)
  23. Archive-name: graphics/algorithms-faq
  24. Posting-Frequency: bi-weekly
  25.  
  26. Welcome to the FAQ for comp.graphics.algorithms!
  27.  
  28. Thanks to all who have contributed.  Corrections and contributions
  29. (to orourke@cs.smith.edu) always welcome.  
  30.  
  31. ----------------------------------------------------------------------
  32. This article is Copyright 2003 by Joseph O'Rourke.  It may be freely 
  33. redistributed in its entirety provided that this copyright notice is 
  34. not removed.
  35. ----------------------------------------------------------------------
  36.  
  37.   Changed items this posting (|): 5.04
  38.   New     items this posting (+): none
  39.  
  40. ----------------------------------------------------------------------
  41. History of Changes (approx. last six months):
  42. ----------------------------------------------------------------------
  43. Changes in 15 Feb 03 posting:
  44.   5.04: Fixed broken link in clipping article. [Thanks to Keith Forbes.]
  45. Changes in  1 Feb 03 posting:
  46.   0.04: Ashdown Radiosity back in print. [Thanks to Ian Ashdown.]
  47.   0.06: Update BSP FAQ links [Thanks to Ken Shoemake.]
  48.   0.07: Update CGAL links in source article.
  49.   3.11: Broken link re course based on Perlin's Noise book. [Thanks to Mikkel Gjoel.]
  50.   6.01: Add CGAL link to Voronoi source article. [Thanks to Andreas Fabri.]
  51.   7.02: All contributor email addresses removed to protect them from spam.
  52. Changes in 15 Jan 03 posting:
  53.   0.06: Query re reality.sgi.com/bspfaq/ [Thanks to <warp@subphase.de>.]
  54.   0.07: Update moved link on WINGED.ZIP. [Thanks to Ben Landon.]
  55.   0.07: Update Ferrar's ++ 3D rendering library link. [Thanks to F.Iannarilli, Jr.]
  56.   1.06: Added ref to AT&T Graphviz. [Thanks to Michael Meire.]
  57.   2.08: Fix sloan ear-clipping link.  [Thanks to logicalink@juno.com.]
  58.   5.09: Update moved link re caustics. [Thanks to Ben Landon.]
  59.   5.18: Formula for distance between two 3D lines. [Thanks to Daniel Zwick.]
  60.   5.27: New article on transforming normals by Ken Shoemake.
  61.   6.08: Random points on sphere in terms of longitude & latitude [Thanks to Uffe Kousgaard.]
  62. Changes in  1 Jul 02 posting:
  63.   3.14: Correct GIF author info, add URL. [Thanks to Greg Roelofs.]
  64. Changes in  1 May 02 posting:
  65.   0.04: Errata for Watt & Watt book added. [Thanks to Jacob Marner.]
  66.   5.14: 3D viewing revised by Ken Shoemake.
  67.   5.23: Remove (erroneous) 3D medial axis info.
  68.   5.25: New article on quaternions by Ken Shoemake.
  69.   5.26: New article on camera aiming and quaternions by Ken Shoemake.
  70.   6.01: Add (correct) 3D medial axis info.  (Thanks to Tamal Dey.)
  71.   6.09: Plucker coordinates article revised by Ken Shoemake.
  72. Changes in 15 Apr 02 posting:
  73.   3.05: Scaling bitmaps revised by Ken Shoemake.
  74.   3.09: Morphing article written by Ken Shoemake.
  75.   6.08: Added references on random points on a sphere (Ken Shoemake).
  76. Changes in  1 Apr 02 posting:
  77.   1.01: 2D point rotation revised by Ken Shoemake.
  78.   1.01: 2D segment intersection revised by Ken Shoemake.
  79.   5.01: 3D point rotation revised by Ken Shoemake.
  80.   0.07: Greg Ferrar's 3D rendering library no longer available.
  81. Changes in 15 Mar 02 posting:
  82.   2.03: Reference Dan Sunday's winding number algorithm.
  83.   4.04: More detail on Beziers approximating a circle.
  84.         (Thanks to William Gibbons.)
  85.   5.22: Added NASA's "Intersect" code for intersecting triangulated
  86.         surfaces.
  87.   5.23: Updated Cocone software description.
  88. Changes in 15 Feb 02 posting:
  89.   5.03: Noted that Sutherland-Hodgman can clip against any convex polygon.
  90.         (Thanks to Ben Landon.)
  91.   5.15: More links on simplifying meshes. (Thanks to Stefan Krause.)
  92. Changes in  1 Jan 02 posting:
  93.   2.03: Fixed link to Franklin's code. (Thanks to Keith M. Briggs.)
  94.   5.13: Update to SWIFT++; add Terdiman's collision lib.
  95.         (Thanks to Pierre Terdiman.)
  96. Changes in  1 Nov 01 posting:
  97.   6.01,02,03: Update to Qhull 3.1 release (Thanks to Brad Barber.)
  98. Changes in 15 Sep 01 posting:
  99.   0.04: "Radiosity: A Programmer's Perspective" out of print.
  100.   0.05: CQUANT97 link no longer available; RADBIB info updated.
  101.         (Thanks to Ian Ashdown for both.)
  102.   2.01: Explained indices in more efficient formula, and restored
  103.         Sunday's version. (Thanks to Dan Sunday.)
  104.   4.04: Link for approximating a circle via a Bezier curve
  105.         (Thanks to John McDonald, Jr.)
  106.   5.10: Add in link to Jules Bloomenthal's list of papers for algorithms
  107.         that could substitute for the marching cubes algorithm.
  108.   5.11: Refer to 5.10. (Thanks to Eric Haines for both.)
  109. Changes in  1 Sep 01 posting:
  110.   2.01: Fixed indices in efficient area formula 
  111.     (Thanks to peter@Glaze.phys.dal.ca.)
  112.   2.03: Link to classic "Point in Polygon Strategies" article.
  113.         (Thanks to Eric Haines.)
  114.   5.09: Additional references for caustics (Thanks to Lars Brinkhoff.)
  115.   5.11: New links for marching cubes patent (Thanks to John Stone.)
  116.   5.17: Stale link notice.
  117.   5.23: New Cocone link for surface reconstruction.
  118. ----------------------------------------------------------------------
  119. Table of Contents
  120. ----------------------------------------------------------------------
  121.  
  122. 0. General Information
  123.    0.01: Charter of comp.graphics.algorithms
  124.    0.02: Are the postings to comp.graphics.algorithms archived?
  125.    0.03: How can I get this FAQ?
  126.    0.04: What are some must-have books on graphics algorithms?
  127.    0.05: Are there any online references?
  128.    0.06: Are there other graphics related FAQs?
  129.    0.07: Where is all the source?
  130.  
  131. 1. 2D Computations: Points, Segments, Circles, Etc.
  132.    1.01: How do I rotate a 2D point?
  133.    1.02: How do I find the distance from a point to a line?
  134.    1.03: How do I find intersections of 2 2D line segments?
  135.    1.04: How do I generate a circle through three points?
  136.    1.05: How can the smallest circle enclosing a set of points be found?
  137.    1.06: Where can I find graph layout algorithms?
  138.  
  139. 2. 2D Polygon Computations
  140.    2.01: How do I find the area of a polygon?
  141.    2.02: How can the centroid of a polygon be computed?
  142.    2.03: How do I find if a point lies within a polygon?
  143.    2.04: How do I find the intersection of two convex polygons?
  144.    2.05: How do I do a hidden surface test (backface culling) with 2D points?
  145.    2.06: How do I find a single point inside a simple polygon?
  146.    2.07: How do I find the orientation of a simple polygon?
  147.    2.08: How can I triangulate a simple polygon?
  148.    2.09: How can I find the minimum area rectangle enclosing a set of points?
  149.  
  150. 3. 2D Image/Pixel Computations
  151.    3.01: How do I rotate a bitmap?
  152.    3.02: How do I display a 24 bit image in 8 bits?
  153.    3.03: How do I fill the area of an arbitrary shape?
  154.    3.04: How do I find the 'edges' in a bitmap?
  155.    3.05: How do I enlarge/sharpen/fuzz a bitmap?
  156.    3.06: How do I map a texture on to a shape?
  157.    3.07: How do I detect a 'corner' in a collection of points?
  158.    3.08: Where do I get source to display (raster font format)?
  159.    3.09: What is morphing/how is it done?
  160.    3.10: How do I quickly draw a filled triangle?
  161.    3.11: D Noise functions and turbulence in Solid texturing.
  162.    3.12: How do I generate realistic sythetic textures?
  163.    3.13: How do I convert between color models (RGB, HLS, CMYK, CIE etc)?
  164.    3.14: How is "GIF" pronounced?
  165.  
  166. 4. Curve Computations
  167.    4.01: How do I generate a Bezier curve that is parallel to another Bezier?
  168.    4.02: How do I split a Bezier at a specific value for t?
  169.    4.03: How do I find a t value at a specific point on a Bezier?
  170.    4.04: How do I fit a Bezier curve to a circle?
  171.  
  172. 5. 3D computations
  173.    5.01: How do I rotate a 3D point?
  174.    5.02: What is ARCBALL and where is the source?
  175.    5.03: How do I clip a polygon against a rectangle?
  176.    5.04: How do I clip a polygon against another polygon?
  177.    5.05: How do I find the intersection of a line and a plane?
  178.    5.06: How do I determine the intersection between a ray and a triangle?
  179.    5.07: How do I determine the intersection between a ray and a sphere?
  180.    5.08: How do I find the intersection of a ray and a Bezier surface?
  181.    5.09: How do I ray trace caustics?
  182.    5.10: What is the marching cubes algorithm?
  183.    5.11: What is the status of the patent on the "marching cubes" algorithm?
  184.    5.12: How do I do a hidden surface test (backface culling) with 3D points?
  185.    5.13: Where can I find algorithms for 3D collision detection?
  186.    5.14: How do I perform basic viewing in 3D?
  187.    5.15: How do I optimize/simplify a 3D polygon mesh?
  188.    5.16: How can I perform volume rendering?
  189.    5.17: Where can I get the spline description of the famous teapot etc.?
  190.    5.18: How can the distance between two lines in space be computed?
  191.    5.19: How can I compute the volume of a polyhedron?
  192.    5.20: How can I decompose a polyhedron into convex pieces?
  193.    5.21: How can the circumsphere of a tetrahedron be computed?
  194.    5.22: How do I determine if two triangles in 3D intersect?
  195.    5.23: How can a 3D surface be reconstructed from a collection of points?
  196.    5.24: How can I find the smallest sphere enclosing a set of points in 3D? 
  197.    5.25: What's the big deal with quaternions?
  198.    5.26: How can I aim a camera in a specific direction?
  199.    5.27: How can I transform normals?
  200.  
  201.  
  202. 6. Geometric Structures and Mathematics
  203.    6.01: Where can I get source for Voronoi/Delaunay triangulation?
  204.    6.02: Where do I get source for convex hull?
  205.    6.03: Where do I get source for halfspace intersection?
  206.    6.04: What are barycentric coordinates?
  207.    6.05: How do I generate a random point inside a triangle?
  208.    6.06: How do I evenly distribute N points on (tesselate) a sphere?
  209.    6.07: What are coordinates for the vertices of an icosohedron?
  210.    6.08: How do I generate random points on the surface of a sphere?
  211.    6.09: What are Plucker coordinates?
  212.  
  213. 7. Contributors
  214.    7.01: How can you contribute to this FAQ?
  215.    7.02: Contributors.  Who made this all possible.
  216.  
  217. Search e.g. for "Section 6" to find that section.
  218. Search e.g. for "Subject 6.04" to find that item.
  219. ----------------------------------------------------------------------
  220. Section 0. General Information
  221. ----------------------------------------------------------------------
  222. Subject 0.01: Charter of comp.graphics.algorithms
  223.  
  224.     comp.graphics.algorithms is an unmoderated newsgroup intended as a forum
  225.     for the discussion of the algorithms used in the process of generating
  226.     computer graphics.  These algorithms may be recently proposed in
  227.     published journals or papers, old or previously known algorithms, or
  228.     hacks used incidental to the process of computer graphics.  The scope of
  229.     these algorithms may range from an efficient way to multiply matrices,
  230.     all the way to a global illumination method incorporating raytracing,
  231.     radiosity, infinite spectrum modeling, and perhaps even mirrored balls
  232.     and lime jello.
  233.  
  234.     It is hoped that this group will serve as a forum for programmers and
  235.     researchers to exchange ideas and ask questions on recent papers or
  236.     current research related to computer graphics.
  237.  
  238.     comp.graphics.algorithms is not:
  239.  
  240.      - for requests for gifs, or other pictures
  241.      - for requests for image translator or processing software; see
  242.             alt.binaries.pictures* FAQ  
  243.             alt.binaries.pictures.utilities [now degenerated to pic postings]
  244.             alt.graphics.pixutils (image format translation)
  245.             comp.sources.misc (image viewing source code)
  246.             sci.image.processing
  247.             comp.graphics.apps.softimage
  248.             fj.comp.image
  249.      - for requests for compression software; for these try:
  250.             alt.comp.compression
  251.             comp.compression
  252.             comp.compression.research
  253.      - specifically for game development; for this try:
  254.             comp.games.development.programming.misc
  255.             comp.games.development.programming.algorithms
  256.  
  257.  
  258. ----------------------------------------------------------------------
  259. Subject 0.02: Are the postings to comp.graphics.algorithms archived?
  260.  
  261.     Archives may be found at: http://www.faqs.org/
  262.  
  263. ----------------------------------------------------------------------
  264. Subject 0.03: How can I get this FAQ?
  265.  
  266.     The FAQ is posted on the 1st and 15th of every month.  The easiest
  267.     way to get it is to search back in your news reader for the most
  268.     recent posting, with Subject: 
  269.           comp.graphics.algorithms Frequently Asked Questions
  270.     It is posted to comp.graphics.algorithms, and cross-posted to
  271.     news.answers and comp.answers.  
  272.  
  273.     If you can't find it on your newsreader,
  274.     you can look at a recent HTML version at the "official" FAQ archive site:
  275.       http://www.faqs.org/
  276.     The maintainer also keeps a copy of the raw ASCII, always the
  277.     latest version, accessible via http://cs.smith.edu/~orourke/FAQ.html .
  278.  
  279.     Finally, you can ftp the FAQ from several sites, including:
  280.  
  281.       ftp://rtfm.mit.edu/pub/faqs/graphics/algorithms-faq
  282.       ftp://mirror.seas.gwu.edu/pub/rtfm/comp/graphics/algorithms/comp.graphics.algorithms_Frequently_Asked_Questions
  283.  
  284.     The (busy) rtfm.mit.edu site lists many alternative "mirror" sites.
  285.     Also can reach the FAQ from http://www.geom.umn.edu/software/cglist/,
  286.     which is worth visiting in its own right.
  287.  
  288. ----------------------------------------------------------------------
  289. Subject 0.04: What are some must-have books on graphics algorithms?
  290.  
  291.     The keywords in brackets are used to refer to the books in later
  292.     questions.  They generally refer to the first author except where
  293.     it is necessary to resolve ambiguity or in the case of the Gems.
  294.  
  295.  
  296.     Basic computer graphics, rendering algorithms,
  297.     ----------------------------------------------
  298.  
  299.     [Foley]
  300.     Computer Graphics: Principles and Practice (2nd Ed.),
  301.     J.D. Foley, A. van Dam, S.K. Feiner, J.F. Hughes, Addison-Wesley
  302.     1990, ISBN 0-201-12110-7;
  303.     Computer Graphics: Principles and Practice, C version
  304.     J.D. Foley,  A. van Dam, S.K. Feiner, J.F. Hughes, Addison-Wesley 
  305.     ISBN: 0-201-84840-6, 1996, 1147 pp.
  306.  
  307.     [Rogers:Procedural]
  308.     Procedural Elements for Computer Graphics, Second Edition
  309.     David F. Rogers, WCB/McGraw Hill 1998, ISBN 0-07-053548-5
  310.  
  311.     [Rogers:Mathematical]
  312.     Mathematical Elements for Computer Graphics 2nd Ed.,
  313.     David F. Rogers and J. Alan Adams, McGraw Hill 1990, ISBN
  314.     0-07-053530-2
  315.  
  316.     [Watt:3D]
  317.     _3D Computer Graphics, 2nd Edition_,
  318.     Alan Watt, Addison-Wesley 1993, ISBN 0-201-63186-5
  319.  
  320.     [Glassner:RayTracing]
  321.     An Introduction to Ray Tracing,
  322.     Andrew Glassner (ed.), Academic Press 1989, ISBN 0-12-286160-4
  323.  
  324.     [Gems I]
  325.     Graphics Gems,
  326.     Andrew Glassner (ed.), Academic Press 1990, ISBN 0-12-286165-5
  327.     http://www.graphicsgems.org/ for all the Gems.
  328.  
  329.     [Gems II]
  330.     Graphics Gems II,
  331.     James Arvo (ed.), Academic Press 1991, ISBN 0-12-64480-0
  332.  
  333.     [Gems III]
  334.     Graphics Gems III,
  335.     David Kirk (ed.), Academic Press 1992, ISBN 0-12-409670-0 (with
  336.     IBM disk) or 0-12-409671-9 (with Mac disk)
  337.     See also "AP Professional Graphics CD-ROM Library,"
  338.     Academic Press,  ISBN 0-12-059756-X, which contains Gems I-III.
  339.  
  340.     [Gems IV]
  341.     Graphics Gems IV,
  342.     Paul S. Heckbert (ed.), Academic Press 1994, ISBN 0-12-336155-9
  343.     (with IBM disk) or 0-12-336156-7 (with Mac disk)
  344.  
  345.     [Gems V]
  346.     Graphic Gems V,
  347.     Alan W. Paeth (ed.), Academic Press 1995, ISBN 0-12-543455-3
  348.     (with IBM disk)
  349.  
  350.     [Watt:Animation]
  351.     Advanced Animation and Rendering Techniques,
  352.     Alan Watt, Mark Watt, Addison-Wesley 1992, ISBN 0-201-54412-1
  353.     (Unofficial) errata: http://www.rolemaker.dk/other/AART/
  354.  
  355.     [Bartels]
  356.     An Introduction to Splines for Use in Computer Graphics and
  357.         Geometric Modeling,
  358.     Richard H. Bartels, John C. Beatty, Brian A. Barsky, 1987, ISBN
  359.     0-934613-27-3
  360.  
  361.     [Farin]
  362.     Curves and Surfaces for Computer Aided Geometric Design:
  363.     A Practical Guide, 4th Edition, Gerald E. Farin, Academic Press
  364.     1996. ISBN 0122490541.
  365.  
  366.     [Prusinkiewicz]
  367.     The Algorithmic Beauty of Plants,
  368.     Przemyslaw W. Prusinkiewicz, Aristid Lindenmayer, Springer-Verlag,
  369.     1990, ISBN 0-387-97297-8, ISBN 3-540-97297-8
  370.  
  371.     [Oliver]
  372.     Tricks of the Graphics Gurus,
  373.     Dick Oliver, et al. (2) 3.5 PC disks included, $39.95 SAMS Publishing
  374.  
  375.     [Hearn]
  376.     Introduction to computer graphics,
  377.     Hearn & Baker
  378.  
  379.     [Cohen]
  380.     Radiosity and Realistic Imange Sythesis,
  381.     Michael F. Cohen, John R. Wallace, Academic Press Professional
  382.     1993, ISBN 0-12-178270-0 [limited reprint 1999]
  383.  
  384.     [Ashdown]
  385.     Radiosity: A Programmer's Perspective
  386.     Ian Ashdown, John Wiley & Sons 1994, ISBN 0-471-30444-1, 498 pp.
  387.     Back in print, Jan 2003.  See www.helios32.com.
  388.  
  389.     [sillion]
  390.     Radiosity & Global Illumination
  391.     Francois X. Sillion snd Claude Puech, Morgan Kaufmann 1994, ISBN
  392.     1-55860-277-1, 252 pp.
  393.  
  394.     [Ebert]
  395.     Texturing and Modeling - A Procedural Approach (2nd Ed.)
  396.     David S. Ebert (ed.), F. Kenton Musgrave, Darwyn Peachey, Ken Perlin,
  397.     Steven Worley, Academic Press 1998, ISBN 0-12-228730-4, Includes CD-ROM.
  398.  
  399.     [Schroeder]
  400.     Visualization Toolkit, 2nd Edition, The: An Object-Oriented Approach to
  401.     3-D Graphics (Bk/CD) (Professional Description)
  402.     William J. Schroeder,  Kenneth Martin, and Bill Lorensen,
  403.     Prentice-Hall 1998, ISBN: 0-13-954694-4
  404.     See Subject 0.07 for source.
  405.  
  406.     [Anderson]
  407.     PC Graphics Unleashed
  408.     Scott Anderson. SAMS Publishing, ISBN 0-672-30570-4
  409.  
  410.     [Ammeraal]
  411.     Computer Graphics for Java Programmers,
  412.     Leen Ammeraal, John Wiley 1998, ISBN 0-471-98142-7.
  413.     Additional information at http://home.wxs.nl/~ammeraal/ .
  414.  
  415.     [Eberly]
  416.     3D Game Engine Design: A Practical Approach to Real-Time Computer Graphics.
  417.     David Eberly, Morgan Kaufmann/Academic Press, 2001.
  418.  
  419.     For image processing,
  420.     ---------------------
  421.  
  422.     [Barnsley]
  423.     Fractal Image Compression,
  424.     Michael F. Barnsley and Lyman P. Hurd, AK Peters, Ltd, 1993 ISBN
  425.     1-56881-000-8
  426.  
  427.     [Jain]
  428.     Fundamentals of Image Processing,
  429.     Anil K. Jain, Prentice-Hall 1989, ISBN 0-13-336165-9
  430.  
  431.     [Castleman]
  432.     Digital Image Processing,
  433.     Kenneth R. Castleman, Prentice-Hall 1996, ISBN(Cloth): 0-13-211467-4
  434.     (Description and errata at: "http://www.phoenix.net/~castlman")
  435.  
  436.     [Pratt]
  437.     Digital Image Processing, Second Edition,
  438.     William K. Pratt, Wiley-Interscience 1991, ISBN 0-471-85766-1
  439.  
  440.     [Gonzalez]
  441.     Digital Image Processing (3rd Ed.),
  442.     Rafael C. Gonzalez, Paul Wintz, Addison-Wesley 1992, ISBN
  443.     0-201-50803-6
  444.  
  445.     [Russ]
  446.     The Image Processing Handbook (3rd Ed.),
  447.     John C. Russ, CRC Press and IEEE Press 1998, ISBN 0-8493-2532-3
  448.     [Russ & Russ]
  449.     The Image Processing Tool Kit v. 3.0
  450.     Chris Russ and John Russ, Reindeer Games Inc. 1999, ISBN 1-928808-00-X
  451.  
  452.  
  453.     [Wolberg]
  454.     Digital Image Warping,
  455.     George Wolberg, IEEE Computer Society Press Monograph 1990, ISBN
  456.     0-8186-8944-7
  457.  
  458.  
  459.     Computational geometry
  460.     ----------------------
  461.  
  462.     [Bowyer]
  463.     A Programmer's Geometry,
  464.     Adrian Bowyer, John Woodwark, Butterworths 1983, 
  465.     ISBN 0-408-01242-0 Pbk
  466.     Out of print, but see:
  467.     Introduction to Computing with Geometry,
  468.     Adrian Bowyer and John Woodwark, 1993
  469.     ISBN 1-874728-03-8.  Available in PDF:
  470.        http://www.inge.com/pubs/index.htm
  471.  
  472.     [Farin & Hansford]
  473.     The Geometry Toolbox for Graphics and Modeling
  474.     by Gerald E. Farin, Dianne Hansford
  475.     A K Peters Ltd; ISBN: 1568810741 
  476.  
  477.     [O'Rourke (C)]
  478.     Computational Geometry in C (2nd Ed.)
  479.     Joseph O'Rourke, Cambridge University Press 1998, 
  480.     ISBN 0-521-64010-5 Pbk, ISBN 0-521-64976-5 Hbk
  481.     Additional information and code at http://cs.smith.edu/~orourke/ .
  482.  
  483.     [O'Rourke (A)]
  484.     Art Gallery Theorems and Algorithms
  485.     Joseph O'Rourke, Oxford University Press 1987,
  486.     ISBN 0-19-503965-3.
  487.  
  488.     [Goodman & O'Rourke]
  489.     Handbook of Discrete and Computational Geometry
  490.     J. E. Goodman and J. O'Rourke, editors.
  491.     CRC Press LLC, July 1997.
  492.     ISBN:0-8493-8524-5
  493.     Additional information at http://cs.smith.edu/~orourke/ .
  494.  
  495.     [Samet:Application]
  496.     Applications of Spatial Data Structures:  Computer Graphics, 
  497.     Image Processing, and GIS, 
  498.     Hanan Samet, Addison-Wesley, Reading, MA, 1990.
  499.     ISBN 0-201-50300-0.
  500.  
  501.     [Samet:Design & Analysis]
  502.     The Design and Analysis of Spatial Data Structures,
  503.     Hanan Samet, Addison-Wesley, Reading, MA, 1990.
  504.     ISBN 0-201-50255-0.
  505.  
  506.     [Mortenson]
  507.     Geometric Modeling,
  508.     Michael E. Mortenson, Wiley 1985, ISBN 0-471-88279-8
  509.  
  510.     [Preparata]
  511.     Computational Geometry: An Introduction,
  512.     Franco P. Preparata, Michael Ian Shamos, Springer-Verlag 1985,
  513.     ISBN 0-387-96131-3
  514.  
  515.     [Okabe]
  516.     Spatial Tessellations: Concepts and Applications of Voronoi Diagrams,
  517.     A. Okabe and B. Boots and K. Sugihara,
  518.     John Wiley, Chichester, England, 1992.
  519.  
  520.     [Overmars]
  521.     Computational Geometry: Algorithms and Applications
  522.     M. de Berg and M. van Kreveld and M. Overmars and O. Schwarzkopf
  523.     Springer-Verlag, Berlin, 1997.
  524.  
  525.     [Stolfi]
  526.     Oriented Projective Geometry: A Framework for Geometric Computations
  527.     Academic Press, 1991.
  528.  
  529.     [Hodge]
  530.     Methods of Algebraic Geometry, Volume 1
  531.     W.V.D. Hodge and D. Pedoe, Cambridge, 1994.
  532.     ISBN 0-521-469007-4 Paperback
  533.  
  534.     [Tamassia et al 199?]
  535.     Graph Drawing: Algorithms for the Visualization of Graphs
  536.     Prentice Hall; ISBN: 0133016153
  537.  
  538.     Algorithms books with chapters on computational geometry
  539.     --------------------------------------------------------
  540.  
  541.     [Cormen et al.]
  542.     Introduction to Algorithms,
  543.     T. H. Cormen, C. E. Leiserson, R. L. Rivest,
  544.     The MIT Press, McGraw-Hill, 1990.
  545.  
  546.     [Mehlhorn]
  547.     Data Structures and Algorithms,
  548.     K. Mehlhorn,
  549.     Springer-Verlag, 1984.
  550.  
  551.     [Sedgewick]
  552.     R. Sedgewick,
  553.     Algorithms,
  554.     Addison-Wesley, 1988.
  555.  
  556.     Solid Modelling
  557.     ---------------
  558.  
  559.     [Mantyla]
  560.     Introduction to Solid Modeling
  561.     Martti Mantyla, Computer Science Press 1988,
  562.     ISBN 07167-8015-1
  563.  
  564. ----------------------------------------------------------------------
  565. Subject 0.05: Are there any online references?
  566.  
  567.     The computational geometry community maintains its own
  568.     bibliography of publications in or closely related to that
  569.     subject.  Every four months, additions and corrections are
  570.     solicited from users, after which the database is updated and
  571.     released anew.  As of 7 Nov 200, it contained 13485 bib-tex
  572.     entries.  See Jeff Erickson's page on "Computational Geometry
  573.     Bibliographies": 
  574.       http://compgeom.cs.uiuc.edu/~jeffe/compgeom/biblios.html#geombib
  575.     The bibliography can be retrieved from:
  576.  
  577.     ftp://ftp.cs.usask.ca/pub/geometry/geombib.tar.gz - bibliography proper
  578.     ftp://ftp.cs.usask.ca/pub/geometry/o-cgc19.ps.gz  - overview published
  579.         in '93 in SIGACT News and the Internat. J. Comput.  Geom. Appl.
  580.     ftp://ftp.cs.usask.ca/pub/geometry/ftp-hints      - detailed retrieval info
  581.  
  582.     Universitat Politecnica de Catalunya maintains a search engine at:
  583.        http://www-ma2.upc.es/~geomc/geombibe.html
  584.  
  585.     The ACM SIGGRAPH Online Bibliography Project, by
  586.     Stephen Spencer (biblio@siggraph.org).
  587.     The database is available for anonymous FTP from the
  588.     ftp://siggraph.org/publications/bibliography directory.  Please
  589.     download and examine the file READ_ME in that directory for more
  590.     specific information concerning the database.
  591.  
  592.     'netlib' is a useful source for algorithms, member inquiries for
  593.     SIAM, and bibliographic searches.  For information, send mail to
  594.     netlib@ornl.gov, with "send index" in the body of the mail
  595.     message.
  596.  
  597.     You can also find free sources for numerical computation in C via
  598.     ftp://ftp.usc.edu/pub/C-numanal/ . In particular, grab
  599.     numcomp-free-c.gz in that directory.
  600.  
  601.     Check out Nick Fotis's computer graphics resources FAQ -- it's
  602.     packed with pointers to all sorts of great computer graphics
  603.     stuff.  This FAQ is posted biweekly to comp.graphics.
  604.  
  605.     This WWW page contains links to a large number
  606.     of computer graphic related pages:
  607.     http://www.dataspace.com:84/vlib/comp-graphics.html
  608.  
  609.     There's a Computer Science Bibliography Server at:
  610.     http://glimpse.cs.arizona.edu:1994/bib/
  611.     with Computer Graphics, Vision and Radiosity sections
  612.  
  613.     A comprehensive bibliography of color quantization papers and articles 
  614.     (CQUANT97) was available at http://www.ledalite.com/library-/cgis.htm.
  615.     [Link no longer available -- replacement? --JOR]
  616.  
  617.     Modelling physically based systems for animation:
  618.     http://www.cc.gatech.edu/gvu/animation/Animation.html
  619.  
  620.     The University of Manchester NURBS Library:
  621.     ftp://unix.hensa.ac.uk/pub/misc/unix/nurbs/
  622.  
  623.     For an implementation of Seidel's algorithm for fast trapezoidation
  624.     and triangulation of polygons. You can get the code from:
  625.     ftp://ftp.cs.unc.edu/pub/users/narkhede/triangulation.tar.gz
  626.  
  627.     Ray tracing bibliography:
  628.     http://www.acm.org/tog/resources/bib/
  629.  
  630.     Quaternions and other comp sci curiosities:
  631.     ftp://ftp.netcom.com/pub/hb/hbaker/hakmem/hakmem.html
  632.  
  633.     Directory of Computational Geometry Software,
  634.     collected by Nina Amenta (nina@cs.utexas.edu)
  635.     Nina Amenta is maintaining a WWW directory to computational 
  636.     geometry software. The directory lives at The Geometry Center. 
  637.     It has pointers to lots of convex hull and voronoi diagram programs, 
  638.     triangulations, collision detection, polygon intersection, smallest 
  639.     enclosing ball of a point set and other stuff.
  640.     http://www.geom.umn.edu/software/cglist/
  641.  
  642.     A compact reference for real-time 3D computer graphics programming:
  643.     http://www.math.mcgill.ca/~loisel/
  644.  
  645.     RADBIB is a comprehensive bibliography of radiosity and
  646.     related global illumination papers, articles, and
  647.     books. It currently includes 1,972 references.
  648.     This bibliography is available in BibTex format 
  649.     (with a release date of 15 Jul 01) from:
  650.       http://www.helios32.com/   under "Resources."
  651.  
  652.     The "Electronic Visualization Library" (EVlib) is a domain-
  653.     secific digital library for Scientific Visualization and 
  654.     Computer Graphics:  http://visinfo.zib.de/
  655.  
  656.     3D Object Intersection: http://www.realtimerendering.com/int/ 
  657.     This page presents information about a wide variety of 3D object/object
  658.     intersection tests. Presented in grid form, each axis lists ray, plane,
  659.     sphere, triangle, box, frustum, and other objects. For each combination
  660.     (e.g. sphere/box), references to articles, books, and online resources
  661.     are given.
  662.  
  663.     Ray Tracing News, ed. Eric Haines:  http://www.raytracingnews.com .
  664.  
  665. ----------------------------------------------------------------------
  666. Subject 0.06: Are there other graphics related FAQs?
  667.  
  668.     BSP Tree FAQ by Bretton Wade
  669.      http://www.andrew.cmu.edu/~rrost/bsp/
  670.      ftp://ftp.sgi.com/other/bspfaq/faq/bspfaq.html
  671.      ftp://ftp.sgi.com/other/bspfaq/faq/bspfaq.txt
  672.      and see
  673.      ftp://ftp.sgi.com/other/bspfaq/
  674.  
  675.     Gamma and Color FAQs by Charles A. Poynton has
  676.         ftp://ftp.inforamp.net/pub/users/poynton/doc/colour/
  677.         http://www.inforamp.net/~poynton/
  678.  
  679.     The documents are mirrored in Darmstadt, Germany at
  680.         ftp://ftp.igd.fhg.de/pub/doc/colour/
  681.  
  682. ----------------------------------------------------------------------
  683. Subject 0.07: Where is all the source?
  684.  
  685.     Graphics Gems source code.
  686.     http://www.graphicsgems.org
  687.     This site is now the offical distribution site for Graphics Gems code.
  688.  
  689.     Master list of Computational Geometry software:
  690.     http://www.geom.umn.edu/software/cglist
  691.     Described in [Goodman & O'Rourke], Chap. 52.
  692.  
  693.     Jeff Erikson's software list:
  694.     http://compgeom.cs.uiuc.edu/~jeffe/compgeom/compgeom.html
  695.  
  696.     Dave Eberly's extensive collection of free geometry, graphics,
  697.     and image processing software:
  698.     http://www.magic-software.com/
  699.  
  700.     General 'stuff'
  701.     ftp://wuarchive.wustl.edu/graphics/graphics
  702.  
  703.     There are a number of interesting items in
  704.     http://graphics.lcs.mit.edu/~seth including:
  705.     - Code for 2D Voronoi, Delaunay, and Convex hull
  706.     - Mike Hoymeyer's implementation of Raimund Seidel's
  707.       O( d! n ) time linear programming algorithm for
  708.       n constraints in d dimensions
  709.     - geometric models of UC Berkeley's new computer science
  710.       building
  711.  
  712.     Sources to "Computational Geometry in C", by J. O'Rourke
  713.     can be found at http://cs.smith.edu/~orourke/books/compgeom.html
  714.     or ftp://cs.smith.edu/pub/compgeom .
  715.  
  716.     Greg Ferrar's C++ 3D rendering library is available at
  717.     http://www.flowerfire.com/ferrar/Graph3D.html  
  718.  
  719.     TAGL is a portable and extensible library that provides a subset
  720.     of Open-GL functionalities.
  721.     ftp://sunsite.unc.edu/pub/packages/programming/graphics/tagl21.tgz
  722.  
  723.     Try ftp://x2ftp.oulu.fi  for /pub/msdos/programming/docs/graphpro.lzh by
  724.     Michael Abrash. His XSharp package has an implementation of Xiaoulin
  725.     Wu's anti-aliasing algorithm (in C).
  726.  
  727.     Example sources for BSP tree algorithms can be found at
  728.     http://reality.sgi.com/bspfaq/, item 24.
  729.  
  730.     Mel Slater (mel@dcs.qmw.ac.uk) also made some implementations of
  731.     BSP trees and shadows for static scenes using shadow volumes
  732.     code available
  733.     http://www.dcs.qmw.ac.uk/~mel/BSP.html
  734.     ftp://ftp.dcs.qmw.ac.uk/people/mel/BSP
  735.  
  736.     The Visualization Toolkit (A visualization textbook, C++ library 
  737.     and Tcl-based interpreter) (see [Schroeder]): 
  738.     http://www.kitware.com/vtk.html
  739.  
  740.     WINGED.ZIP, a C++ implementation of Baumgart's winged-edge data structure:
  741.     ftp://ftp.ledalite.com/pub/winged.zip   
  742.  
  743.     CGAL, the Computational Geometry Algorithms Library, is written in C++ 
  744.     and is available at http://www.cgal.org.  
  745.     CGAL contains algorithms and data structures for 2D computations
  746.     (convex hull, Delaunay, constrained Delaunay, Voronoi diagram, 
  747.     regular traingulation, (weighted) Alpha shapes, polytope distance, 
  748.     boolean operations on polygons, decomposition of polygons in
  749.     monotone or convex parts, arrangements, etc.), 3D, and arbitrary
  750.     dimensions.
  751.  
  752.     A C++ NURBS library written by Lavoie Philippe. Version 2.1. 
  753.     Results may be exported as POV-Ray, RIB (renderman) or VRML files.
  754.     It also offers wrappers to OpenGL:
  755.        http://yukon.genie.uottawa.ca/~lavoie/software/nurbs/
  756.  
  757.     Paul Bourke has code for several problems, including isosurface 
  758.     generation and Delauney triangulation, at:
  759.        http://www.swin.edu.au/astronomy/pbourke/geometry/
  760.        http://www.swin.edu.au/astronomy/pbourke/modeling/
  761.  
  762.     A nearly comprehensive list of available 3D engines 
  763.     (most with source code):
  764.        http://cg.cs.tu-berlin.de/~ki/engines.html
  765.  
  766.     See also 5.17: 
  767.         Where can I get the spline description of the famous teapot etc.?
  768.  
  769.     Interactive Geometry Software called "Cinderella":
  770.         http://www.cinderella.de
  771.  
  772. ----------------------------------------------------------------------
  773. Section 1. 2D Computations: Points, Segments, Circles, Etc.
  774. ----------------------------------------------------------------------
  775. Subject 1.01: How do I rotate a 2D point?
  776.  
  777.     In 2D, you make (X,Y) from (x,y) with a rotation by angle t so:
  778.         X = x cos t - y sin t
  779.         Y = x sin t + y cos t
  780.     As a 2x2 matrix this is very simple.  If you want to rotate a
  781.     column vector v by t degrees using matrix M, use
  782.         M = [cos t  -sin t]
  783.             [sin t   cos t]
  784.     in the product M v.
  785.  
  786.     If you have a row vector, use the transpose of M (turn rows into
  787.     columns and vice versa).  If you want to combine rotations, in 2D
  788.     you can just add their angles, but in higher dimensions you must
  789.     multiply their matrices.
  790.  
  791. ----------------------------------------------------------------------
  792. Subject 1.02: How do I find the distance from a point to a line?
  793.  
  794.  
  795.     Let the point be C (Cx,Cy) and the line be AB (Ax,Ay) to (Bx,By).
  796.     Let P be the point of perpendicular projection of C on AB.  The parameter
  797.     r, which indicates P's position along AB, is computed by the dot product 
  798.     of AC and AB divided by the square of the length of AB:
  799.     
  800.     (1)     AC dot AB
  801.         r = ---------  
  802.             ||AB||^2
  803.     
  804.     r has the following meaning:
  805.     
  806.         r=0      P = A
  807.         r=1      P = B
  808.         r<0      P is on the backward extension of AB
  809.         r>1      P is on the forward extension of AB
  810.         0<r<1    P is interior to AB
  811.     
  812.     The length of a line segment in d dimensions, AB is computed by:
  813.     
  814.         L = sqrt( (Bx-Ax)^2 + (By-Ay)^2 + ... + (Bd-Ad)^2)
  815.  
  816.     so in 2D:   
  817.     
  818.         L = sqrt( (Bx-Ax)^2 + (By-Ay)^2 )
  819.     
  820.     and the dot product of two vectors in d dimensions, U dot V is computed:
  821.     
  822.         D = (Ux * Vx) + (Uy * Vy) + ... + (Ud * Vd)
  823.     
  824.     so in 2D:   
  825.     
  826.         D = (Ux * Vx) + (Uy * Vy) 
  827.     
  828.     So (1) expands to:
  829.     
  830.             (Cx-Ax)(Bx-Ax) + (Cy-Ay)(By-Ay)
  831.         r = -------------------------------
  832.                           L^2
  833.  
  834.     The point P can then be found:
  835.  
  836.         Px = Ax + r(Bx-Ax)
  837.         Py = Ay + r(By-Ay)
  838.  
  839.     And the distance from A to P = r*L.
  840.  
  841.     Use another parameter s to indicate the location along PC, with the 
  842.     following meaning:
  843.            s<0      C is left of AB
  844.            s>0      C is right of AB
  845.            s=0      C is on AB
  846.  
  847.     Compute s as follows:
  848.  
  849.             (Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay)
  850.         s = -----------------------------
  851.                         L^2
  852.  
  853.  
  854.     Then the distance from C to P = |s|*L.
  855.  
  856.  
  857. ----------------------------------------------------------------------
  858. Subject 1.03: How do I find intersections of 2 2D line segments?
  859.  
  860.     This problem can be extremely easy or extremely difficult; it
  861.     depends on your application. If all you want is the intersection
  862.     point, the following should work:
  863.  
  864.     Let A,B,C,D be 2-space position vectors.  Then the directed line
  865.     segments AB & CD are given by:
  866.  
  867.         AB=A+r(B-A), r in [0,1]
  868.         CD=C+s(D-C), s in [0,1]
  869.  
  870.     If AB & CD intersect, then
  871.  
  872.         A+r(B-A)=C+s(D-C), or
  873.  
  874.         Ax+r(Bx-Ax)=Cx+s(Dx-Cx)
  875.         Ay+r(By-Ay)=Cy+s(Dy-Cy)  for some r,s in [0,1]
  876.  
  877.     Solving the above for r and s yields
  878.  
  879.             (Ay-Cy)(Dx-Cx)-(Ax-Cx)(Dy-Cy)
  880.         r = -----------------------------  (eqn 1)
  881.             (Bx-Ax)(Dy-Cy)-(By-Ay)(Dx-Cx)
  882.  
  883.             (Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay)
  884.         s = -----------------------------  (eqn 2)
  885.             (Bx-Ax)(Dy-Cy)-(By-Ay)(Dx-Cx)
  886.  
  887.     Let P be the position vector of the intersection point, then
  888.  
  889.         P=A+r(B-A) or
  890.  
  891.         Px=Ax+r(Bx-Ax)
  892.         Py=Ay+r(By-Ay)
  893.  
  894.     By examining the values of r & s, you can also determine some
  895.     other limiting conditions:
  896.  
  897.         If 0<=r<=1 & 0<=s<=1, intersection exists
  898.             r<0 or r>1 or s<0 or s>1 line segments do not intersect
  899.  
  900.         If the denominator in eqn 1 is zero, AB & CD are parallel
  901.         If the numerator in eqn 1 is also zero, AB & CD are collinear.
  902.  
  903.     If they are collinear, then the segments may be projected to the x- 
  904.     or y-axis, and overlap of the projected intervals checked.
  905.  
  906.     If the intersection point of the 2 lines are needed (lines in this
  907.     context mean infinite lines) regardless whether the two line
  908.     segments intersect, then
  909.  
  910.         If r>1, P is located on extension of AB
  911.         If r<0, P is located on extension of BA
  912.         If s>1, P is located on extension of CD
  913.         If s<0, P is located on extension of DC
  914.  
  915.     Also note that the denominators of eqn 1 & 2 are identical.
  916.  
  917.     References:
  918.  
  919.     [O'Rourke (C)] pp. 249-51
  920.     [Gems III] pp. 199-202 "Faster Line Segment Intersection,"
  921.  
  922. ----------------------------------------------------------------------
  923. Subject 1.04: How do I generate a circle through three points?
  924.  
  925.     Let the three given points be a, b, c.  Use _0 and _1 to represent
  926.     x and y coordinates. The coordinates of the center p=(p_0,p_1) of
  927.     the circle determined by a, b, and c are:
  928.  
  929.         A = b_0 - a_0;
  930.         B = b_1 - a_1;
  931.         C = c_0 - a_0;
  932.         D = c_1 - a_1;
  933.     
  934.         E = A*(a_0 + b_0) + B*(a_1 + b_1);
  935.         F = C*(a_0 + c_0) + D*(a_1 + c_1);
  936.     
  937.         G = 2.0*(A*(c_1 - b_1)-B*(c_0 - b_0));
  938.     
  939.         p_0 = (D*E - B*F) / G;
  940.         p_1 = (A*F - C*E) / G;
  941.  
  942.     If G is zero then the three points are collinear and no finite-radius
  943.     circle through them exists.  Otherwise, the radius of the circle is:
  944.  
  945.             r^2 = (a_0 - p_0)^2 + (a_1 - p_1)^2
  946.  
  947.     Reference:
  948.  
  949.     [O' Rourke (C)] p. 201. Simplified by Jim Ward.
  950.  
  951. ----------------------------------------------------------------------
  952. Subject 1.05: How can the smallest circle enclosing a set of points be found?
  953.  
  954.     This circle is often called the minimum spanning circle.  It can be
  955.     computed in O(n log n) time for n points.  The center lies on
  956.     the furthest-point Voronoi diagram.  Computing the diagram constrains
  957.     the search for the center.  Constructing the diagram can be accomplished
  958.     by a 3D convex hull algorithm; for that connection, see, e.g.,
  959.     [O'Rourke (C), p.195ff].  For direct algorithms, see:
  960.       S. Skyum, "A simple algorithm for computing the smallest enclosing circle"
  961.       Inform. Process. Lett. 37 (1991) 121--125.
  962.       J. Rokne, "An Easy Bounding Circle" [Gems II] pp.14--16.
  963.  
  964. ----------------------------------------------------------------------
  965. Subject 1.06: Where can I find graph layout algorithms?
  966.  
  967.     See the paper by Eades and Tamassia, "Algorithms for Drawing
  968.     Graphs: An Annotated Bibliography," Tech Rep CS-89-09, Dept.  of
  969.     CS, Brown Univ, Feb. 1989.
  970.  
  971.     A revised version of the annotated bibliography on graph drawing
  972.     algorithms by Giuseppe Di Battista, Peter Eades, and Roberto
  973.     Tamassia is now available from
  974.  
  975.     ftp://wilma.cs.brown.edu/pub/papers/compgeo/gdbiblio.tex.gz and
  976.     ftp://wilma.cs.brown.edu/pub/papers/compgeo/gdbiblio.ps.gz
  977.  
  978.     Commercial software includes the Graph Layout Toolkit from Tom Sawyer
  979.     Software http://www.tomsawyer.com and Northwoods Software's GO++
  980.     http://www.nwoods.com/go/ .
  981.  
  982.     Perhaps the best code is the AT&T Research Labs open C source:
  983.     http://www.research.att.com/sw/tools/graphviz/
  984.  
  985. ----------------------------------------------------------------------
  986. Section 2. 2D Polygon Computations
  987. ----------------------------------------------------------------------
  988. Subject 2.01: How do I find the area of a polygon?
  989.  
  990.     The signed area can be computed in linear time by a simple sum.
  991.     The key formula is this:
  992.  
  993.         If the coordinates of vertex v_i are x_i and y_i,
  994.         twice the signed area of a polygon is given by
  995.  
  996.         2 A( P ) = sum_{i=0}^{n-1} (x_i y_{i+1} - y_i x_{i+1}).
  997.  
  998.     Here n is the number of vertices of the polygon, and index
  999.     arithmetic is mod n, so that x_n = x_0, etc. A rearrangement
  1000.     of terms in this equation can save multiplications and operate on
  1001.     coordinate differences, and so may be both faster and more
  1002.     accurate:
  1003.  
  1004.        2 A(P) = sum_{i=0}^{n-1} ( x_i  (y_{i+1} - y_{i-1}) )
  1005.  
  1006.     Here again modular index arithmetic is implied, with n=0 and -1=n-1.
  1007.     This can be avoided by extending the x[] and y[] arrays up to [n+1]
  1008.     with x[n]=x[0], y[n]=y[0] and y[n+1]=y[1], and using instead
  1009.  
  1010.        2 A(P) = sum_{i=1}^{n} ( x_i  (y_{i+1} - y_{i-1}) )
  1011.  
  1012.  
  1013.     References: [O' Rourke (C)] Thm. 1.3.3, p. 21; [Gems II] pp. 5-6:
  1014.     "The Area of a Simple Polygon", Jon Rokne.  Dan Sunday's explanation:
  1015.        http://GeometryAlgorithms.com/Archive/algorithm_0101/  where
  1016.  
  1017.     To find the area of a planar polygon not in the x-y plane, use:
  1018.     
  1019.        2 A(P) = abs(N . (sum_{i=0}^{n-1} (v_i x v_{i+1})))
  1020.     
  1021.        where N is a unit vector normal to the plane. The `.' represents the
  1022.        dot product operator, the `x' represents the cross product operator,
  1023.        and abs() is the absolute value function.
  1024.     
  1025.     [Gems II] pp. 170-171:
  1026.     "Area of Planar Polygons and Volume of Polyhedra", Ronald N. Goldman.
  1027.  
  1028. ----------------------------------------------------------------------
  1029. Subject 2.02: How can the centroid of a polygon be computed?
  1030.  
  1031.     The centroid (a.k.a. the center of mass, or center of gravity)
  1032.     of a polygon can be computed as the weighted sum of the centroids
  1033.     of a partition of the polygon into triangles.  The centroid of a
  1034.     triangle is simply the average of its three vertices, i.e., it
  1035.     has coordinates (x1 + x2 + x3)/3 and (y1 + y2 + y3)/3.  This 
  1036.     suggests first triangulating the polygon, then forming a sum
  1037.     of the centroids of each triangle, weighted by the area of
  1038.     each triangle, the whole sum normalized by the total polygon area.
  1039.     This indeed works, but there is a simpler method:  the triangulation
  1040.     need not be a partition, but rather can use positively and
  1041.     negatively oriented triangles (with positive and negative areas),
  1042.     as is used when computing the area of a polygon.  This leads to
  1043.     a very simple algorithm for computing the centroid, based on a
  1044.     sum of triangle centroids weighted with their signed area.
  1045.     The triangles can be taken to be those formed by any fixed point,
  1046.     e.g., the vertex v0 of the polygon, and the two endpoints of 
  1047.     consecutive edges of the polygon: (v1,v2), (v2,v3), etc.  The area 
  1048.     of a triangle with vertices a, b, c is half of this expression:
  1049.                 (b[X] - a[X]) * (c[Y] - a[Y]) -
  1050.                 (c[X] - a[X]) * (b[Y] - a[Y]);
  1051.  
  1052.     Code available at ftp://cs.smith.edu/pub/code/centroid.c (3K).
  1053.     Reference: [Gems IV] pp.3-6; also includes code.
  1054.  
  1055. ----------------------------------------------------------------------
  1056. Subject 2.03: How do I find if a point lies within a polygon?
  1057.  
  1058.     The definitive reference is "Point in Polygon Strategies" by
  1059.     Eric Haines [Gems IV]  pp. 24-46.  Now also at 
  1060.        http://www.erichaines.com/ptinpoly.
  1061.     The code in the Sedgewick book Algorithms (2nd Edition, p.354) fails
  1062.     under certain circumstances.  See 
  1063.        http://condor.informatik.Uni-Oldenburg.DE/~stueker/graphic/index.html
  1064.     for a discussion.
  1065.  
  1066.     The essence of the ray-crossing method is as follows.
  1067.     Think of standing inside a field with a fence representing the polygon. 
  1068.     Then walk north. If you have to jump the fence you know you are now 
  1069.     outside the poly. If you have to cross again you know you are now 
  1070.     inside again; i.e., if you were inside the field to start with, the total 
  1071.     number of fence jumps you would make will be odd, whereas if you were 
  1072.     ouside the jumps will be even.
  1073.  
  1074.     The code below is from Wm. Randolph Franklin <wrf@ecse.rpi.edu>
  1075.     (see URL below) with some minor modifications for speed.  It returns 
  1076.     1 for strictly interior points, 0 for strictly exterior, and 0 or 1 
  1077.     for points on the boundary.  The boundary behavior is complex but 
  1078.     determined; in particular, for a partition of a region into polygons, 
  1079.     each point is "in" exactly one polygon.  
  1080.     (See p.243 of [O'Rourke (C)] for a discussion of boundary behavior.)
  1081.  
  1082.     int pnpoly(int npol, float *xp, float *yp, float x, float y)
  1083.     {
  1084.       int i, j, c = 0;
  1085.       for (i = 0, j = npol-1; i < npol; j = i++) {
  1086.         if ((((yp[i]<=y) && (y<yp[j])) ||
  1087.              ((yp[j]<=y) && (y<yp[i]))) &&
  1088.             (x < (xp[j] - xp[i]) * (y - yp[i]) / (yp[j] - yp[i]) + xp[i]))
  1089.  
  1090.           c = !c;
  1091.       }
  1092.       return c;
  1093.     }
  1094.  
  1095.     The code may be further accelerated, at some loss in clarity, by
  1096.     avoiding the central computation when the inequality can be deduced,
  1097.     and by replacing the division by a multiplication for those processors
  1098.     with slow divides.  For code that distinguishes strictly interior
  1099.     points from those on the boundary, see [O'Rourke (C)] pp. 239-245.
  1100.     For a method based on winding number, see Dan Sunday,
  1101.     "Fast Winding Number Test for Point Inclusion in a Polygon,"
  1102.     http://softsurfer.com/algorithms.htm, March 2001.
  1103.  
  1104.     References:
  1105.     Franklin's code: 
  1106.        http://www.ecse.rpi.edu/Homepages/wrf/research/geom/pnpoly.html
  1107.     [Gems IV]  pp. 24-46
  1108.     [O'Rourke (C)] Sec. 7.4.
  1109.     [Glassner:RayTracing]
  1110.  
  1111. ----------------------------------------------------------------------
  1112. Subject 2.04: How do I find the intersection of two convex polygons?
  1113.  
  1114.     Unlike intersections of general polygons, which might produce a
  1115.     quadratic number of pieces, intersection of convex polygons of n
  1116.     and m vertices always produces a polygon of at most (n+m) vertices.
  1117.     There are a variety of algorithms whose time complexity is proportional
  1118.     to this size, i.e., linear.  The first, due to Shamos and Hoey, is
  1119.     perhaps the easiest to understand.  Let the two polygons be P and
  1120.     Q.  Draw a horizontal line through every vertex of each.  This
  1121.     partitions each into trapezoids (with at most two triangles, one
  1122.     at the top and one at the bottom).  Now scan down the two polygons
  1123.     in concert, one trapezoid at a time, and intersect the trapezoids
  1124.     if they overlap vertically.
  1125.        A more difficult-to-describe algorithm is in [O'Rourke (C)], pp. 252-262.
  1126.     This walks around the boundaries of P and Q in concert, intersecting
  1127.     edges.  An implementation is included in [O'Rourke (C)].
  1128.  
  1129.  
  1130. ----------------------------------------------------------------------
  1131. Subject 2.05: How do I do a hidden surface test (backface culling) with 2D points?
  1132.  
  1133.     c = (x1-x2)*(y3-y2)-(y1-y2)*(x3-x2)
  1134.  
  1135.     x1,y1, x2,y2, x3,y3 = the first three points of the polygon.
  1136.  
  1137.     If c is positive the polygon is visible. If c is negative the
  1138.     polygon is invisible (or the other way).
  1139.  
  1140.  
  1141. ----------------------------------------------------------------------
  1142. Subject 2.06: How do I find a single point inside a simple polygon?
  1143.  
  1144.    Given a simple polygon, find some point inside it.  Here is a method
  1145.    based on the proof that there exists an internal diagonal, in
  1146.    [O'Rourke (C), 13-14].  The idea is that the midpoint of a diagonal
  1147.    is interior to the polygon.
  1148.  
  1149.    1. Identify a convex vertex v; let its adjacent vertices be a and b.
  1150.    2. For each other vertex q do:
  1151.        2a. If q is inside avb, compute distance to v (orthogonal to ab).
  1152.        2b. Save point q if distance is a new min.
  1153.    3. If no point is inside, return midpoint of ab, or centroid of avb.
  1154.    4. Else if some point inside, qv is internal: return its midpoint.
  1155.  
  1156.    Code for finding a diagonal is in [O'Rourke (C), 35-39], and is part
  1157.    of many other software packages.  See Subject 0.07: Where is all the 
  1158.    source?
  1159.  
  1160. ----------------------------------------------------------------------
  1161. Subject 2.07: How do I find the orientation of a simple polygon?
  1162.  
  1163.     Compute the signed area (Subject 2.01).  The orientation is 
  1164.     counter-clockwise if this area is positive.  
  1165.  
  1166.     A slightly faster method is based on the observation that it isn't
  1167.     necessary to compute the area.  Find the lowest vertex (or, if 
  1168.     there is more than one vertex with the same lowest coordinate, 
  1169.     the rightmost of those vertices) and then take the cross product 
  1170.     of the edges fore and aft of it.  Both methods are O(n) for n vertices, 
  1171.     but it does seem a waste to add up the total area when a single cross
  1172.     product (of just the right edges) suffices.  Code for this is
  1173.     available at ftp://cs.smith.edu/pub/code/polyorient.C (2K).
  1174.  
  1175.     The reason that the lowest, rightmost (or any other such extreme) point
  1176.     works is that the internal angle at this vertex is necessarily convex, 
  1177.     strictly less than pi (even if there are several equally-lowest points).
  1178.  
  1179. ----------------------------------------------------------------------
  1180. Subject 2.08: How can I triangulate a simple polygon?
  1181.  
  1182.     Triangulation of a polygon partitions its interior into triangles
  1183.     with disjoint interiors.  Usually one restricts corners of the
  1184.     triangles to coincide with vertices of the polygon, in which case
  1185.     every polygon of n vertices can be triangulated, and all triangulations
  1186.     contain n-2 triangles, employing n-3 "diagonals" (chords between
  1187.     vertices that otherwise do not touch the boundary of the polygon).
  1188.     
  1189.     Triangulations can be constructed by a variety of algorithms,
  1190.     ranging from a naive search for noncrossing diagonals, which is
  1191.     O(n^4), to "ear" clipping, which is O(n^2) and relatively straightforward
  1192.     to implement [O'Rourke (C), Chap. 1], to partitioning into
  1193.     monotone polygons, which leads to O(n log n) time complexity
  1194.     [O'Rourke (C), Chap. 2; Overmars, Chap. 3], to several randomized 
  1195.     algorithms---by Clarkson et al., by Seidel, and by Devillers, which 
  1196.     lead to O(n log* n) complexity---to Chazelle's linear-time algorithm, 
  1197.     which has yet to be implemented.
  1198.  
  1199.     There is a tradeoff between code complexity and time complexity.
  1200.     Fortunately, several of the algorithms have been implemented and are
  1201.     available:
  1202.          Ear-clipping: 
  1203.                 http://cs.smith.edu/~orourke/books/compgeom.html 
  1204.                 ftp://ftp.cis.uab.edu/pub/sloan/Software/triangulation/src/
  1205.          Seidel's Alg: 
  1206.                 http://www.cs.unc.edu/~dm/CODE/GEM/chapter.html
  1207.                 ftp://ftp.cs.unc.edu/pub/users/narkhede/triangulation.tar.gz
  1208.                 http://reality.sgi.com/atul/code/chapter.html
  1209.     See also the collection of triangulation links at
  1210.                 http://www.geom.umn.edu/software/cglist/
  1211.     
  1212.     References:
  1213.  
  1214.     [O'Rourke (C)]
  1215.     [Overmars]
  1216.     [Gems V]
  1217.     Clarkson, K., Tarjan, R., and VanWyk, C. A fast Las Vegas algorithm for
  1218.       triangulating a simple polygon. Discrete and Computational Geometry,
  1219.       4(1):423--432, 1989.
  1220.     Clarkson, K., Cole, R., Tarjan, R. Randomized parallel algorithms for
  1221.       trapezoidal diagrams. Int. J. Comp. Geom. Appl., 117--133, 1992.
  1222.       http://cm.bell-labs.com/cm/cs/who/clarkson/tri.html
  1223.     Seidel, R.  (1991), A simple and fast incremental randomized algorithm for
  1224.       computing trapezoidal decompositions and for triangulating polygons,
  1225.       Comput. Geom. Theory Appl. 1, 51--64.
  1226.     Devillers, O.  (1992), Randomization yields simple O(n log* n)
  1227.       algorithms for difficult Omega(n) problems, 
  1228.       Internat. J. Comput.  Geom. Appl. 2(1), 97--111.
  1229.     Chazelle, B.  (1991), Triangulating a simple polygon in linear time,
  1230.       Discrete Comput. Geom. 6, 485--524.
  1231.     Held, M. (1998) "FIST: Fast Industrial-Strength Triangulation". 
  1232.       http://www.cosy.sbg.ac.at/~held/projects/triang/triang.html
  1233.  
  1234. ----------------------------------------------------------------------
  1235. Subject 2.09: How can I find the minimum area rectangle enclosing a set of points?
  1236.     First take the convex hull of the points.  Let the resulting convex
  1237.     polygon be P.  It has been known for some time that the minimum
  1238.     area rectangle enclosing P must have one rectangle side flush with
  1239.     (i.e., collinear with and overlapping) one edge of P.  This geometric
  1240.     fact was used by Godfried Toussaint to develop the "rotating calipers"
  1241.     algorithm in a hard-to-find 1983 paper, "Solving Geometric Problems 
  1242.     with the Rotating Calipers" (Proc. IEEE MELECON).  The algorithm
  1243.     rotates a surrounding rectangle from one flush edge to the next,
  1244.     keeping track of the minimum area for each edge.  It achieves O(n)
  1245.     time (after hull computation).  See the "Rotating Calipers Homepage"
  1246.     http://www.cs.mcgill.ca/~orm/rotcal.frame.html for a description
  1247.     and applet.
  1248.  
  1249. ----------------------------------------------------------------------
  1250. Section 3. 2D Image/Pixel Computations
  1251. ----------------------------------------------------------------------
  1252. Subject 3.01: How do I rotate a bitmap?
  1253.  
  1254.     The easiest way, according to the comp.graphics faq, is to take
  1255.     the rotation transformation and invert it. Then you just iterate
  1256.     over the destination image, apply this inverse transformation and
  1257.     find which source pixel to copy there.
  1258.  
  1259.     A much nicer way comes from the observation that the rotation
  1260.     matrix:
  1261.  
  1262.         R(T) = { { cos(T), -sin(T) }, { sin(T), cos(T) } }
  1263.  
  1264.     is formed my multiplying three matrices, namely:
  1265.  
  1266.         R(T) = M1(T) * M2(T) * M3(T)
  1267.  
  1268.     where
  1269.  
  1270.         M1(T) = { { 1, -tan(T/2) },
  1271.                   { 0, 1         } }
  1272.         M2(T) = { { 1,      0    },
  1273.                   { sin(T), 1    } }
  1274.         M3(T) = { { 1, -tan(T/2) },
  1275.                   { 0,         1 } }
  1276.  
  1277.     Each transformation can be performed in a separate pass, and
  1278.     because these transformations are either row-preserving or
  1279.     column-preserving, anti-aliasing is quite easy.
  1280.  
  1281.     Another fast approach is to perform first a column-preserving roation, 
  1282.     and then a row-preserving rotation. For an image W pixels wide and 
  1283.     H pixels high, this requires W+H BitBlt operations in comparison to 
  1284.     the brute-force rotation, which uses W*H SetPixel operations (and a 
  1285.     lot of multiplying).
  1286.  
  1287.     Reference:
  1288.  
  1289.     Paeth, A. W., "A Fast Algorithm for General Raster Rotation",
  1290.     Proceedings Graphics Interface '89, Canadian Information
  1291.     Processing Society, 1986, 77-81
  1292.     [Note - e-mail copies of this paper are no longer available]
  1293.  
  1294.     [Gems I]
  1295.  
  1296. ----------------------------------------------------------------------
  1297. Subject 3.02: How do I display a 24 bit image in 8 bits?
  1298.  
  1299.     [Gems I] pp. 287-293, "A Simple Method for Color Quantization:
  1300.     Octree Quantization"
  1301.  
  1302.     B. Kurz.  Optimal Color Quantization for Color Displays.
  1303.     Proceedings of the IEEE Conference on Computer Vision and Pattern
  1304.     Recognition, 1983, pp. 217-224.
  1305.  
  1306.     [Gems II] pp. 116-125, "Efficient Inverse Color Map Computation"
  1307.  
  1308.         This describes an efficient technique to
  1309.         map actual colors to a reduced color map,
  1310.         selected by some other technique described
  1311.         in the other papers.
  1312.  
  1313.     [Gems II] pp. 126-133, "Efficient Statistical Computations for
  1314.     Optimal Color Quantization"
  1315.  
  1316.     Xiaolin Wu.  Color Quantization by Dynamic Programming and
  1317.     Principal Analysis.  ACM Transactions on Graphics, Vol. 11, No. 4,
  1318.     October 1992, pp 348-372.
  1319.  
  1320.  
  1321. ----------------------------------------------------------------------
  1322. Subject 3.03: How do I fill the area of an arbitrary shape?
  1323.  
  1324.  
  1325.     "A Fast Algorithm for the Restoration of Images Based on Chain
  1326.         Codes Description and Its Applications", L.W. Chang & K.L. Leu,
  1327.         Computer Vision, Graphics, and Image Processing, vol.50,
  1328.         pp296-307 (1990)
  1329.  
  1330.      Heckbert, Paul S., Generic Convex Polygon Scan Conversion and Clipping,
  1331.         Graphics Gems, p. 84-86, code: p. 667-680, PolyScan/.
  1332.      Heckbert, Paul S., Concave Polygon Scan Conversion, Graphics Gems, p.
  1333.         87-91, code: p. 681-684, ConcaveScan.c.
  1334.      http://www.acm.org/tog/GraphicsGems/gems/PolyScan/
  1335.      http://www.acm.org/tog/GraphicsGems/gems/ConcaveScan.c
  1336.  
  1337.      For filling a region of a bitmap, see
  1338.      Heckbert, Paul S., A Seed Fill Algorithm, Graphics Gems, p. 275-277,
  1339.         code: p. 721-722, SeedFill.c. The code is at
  1340.         http://www.acm.org/tog/GraphicsGems/gems/SeedFill.c
  1341.  
  1342.  
  1343.     [Gems I]
  1344.     [Foley]
  1345.     [Hearn]
  1346.  
  1347. ----------------------------------------------------------------------
  1348. Subject 3.04: How do I find the 'edges' in a bitmap?
  1349.  
  1350.     A simple method is to put the bitmap through the filter:
  1351.  
  1352.         -1    -1    -1
  1353.         -1     8    -1
  1354.         -1    -1    -1
  1355.  
  1356.     This will highlight changes in contrast.  Then any part of the
  1357.     picture where the absolute filtered value is higher than some
  1358.     threshold is an "edge".
  1359.  
  1360.     A more appropriate edge detector for noisy images is
  1361.     described by Van Vliet et al. "A nonlinear Laplace operator
  1362.     as edge detector in noisy images", in Computer Vision,
  1363.     Graphics, and image processing 45, pp. 167-195, 1989.
  1364.  
  1365.  
  1366. ----------------------------------------------------------------------
  1367. Subject 3.05: How do I enlarge/sharpen/fuzz a bitmap?
  1368.  
  1369.     Sharpening of bitmaps can be done by the following algorithm:
  1370.  
  1371.     I_enh(x,y) = I_fuz(x,y)-k*Laplace(I_fuz(x,y))
  1372.  
  1373.     or in words:  An image can be sharpened by subtracting a positive
  1374.     fraction k of the Laplace from the fuzzy image.
  1375.  
  1376.     One "Laplace" kernel, approximating a Laplacian operator, is:
  1377.         1   1   1
  1378.         1  -8   1
  1379.         1   1   1
  1380.  
  1381.  
  1382.     The following library implements Fast Gaussian Blurs:
  1383.  
  1384.     MAGIC: An Object-Oriented Library for Image Analysis by David Eberly
  1385.  
  1386.     The library source code and the documentation (in Latex) are  at
  1387.         http://www.magic-software.com/
  1388.     The code compiles on Unix systems using g++ and on PCs using
  1389.     Microsoft Windows 3.1 and Borland C++.  The fast Gaussian blurring
  1390.     is based on a finite difference method for solving s u_s = s^2 \nabla^2 u
  1391.     where s is the standard deviation of the Gaussian (t = s^2/2).  It
  1392.     takes advantage of geometrically increasing steps in s (rather than
  1393.     linearly increasing steps in t), thus getting to a larger "time" rapidly,
  1394.     but still retaining stability.  Section 4.5 of the  documentation contains
  1395.     the algorithm description and implementation.
  1396.  
  1397.      A bitmap is a sampled image, a special case of a digital signal,
  1398.      and suffers from two limitations common to all digital signals.
  1399.      First, it cannot provide details at fine enough spacing to exactly
  1400.      reproduce every continuous image, nor even more detailed sampled
  1401.      images. And second, each sample approximates the infinitely fine
  1402.      variability of ideal values with a discrete set of ranges encoded
  1403.      in a small number of bits---sometimes just one bit per pixel. Most
  1404.      bitmaps have another limitation imposed: The values cannot be
  1405.      negative. The resolution limitation is especially important, but
  1406.      see "How do I display a 24 bit image in 8 bits?" for range issues.
  1407.  
  1408.      The ideal way to enlarge a bitmap is to work from the original
  1409.      continuous image, magnifying and resampling it. The standard way
  1410.      to do it in practice is to (conceptually) reconstruct a continuous
  1411.      image from the bitmap, and magnify and resample that instead. This
  1412.      will not give the same results, since details of the original have
  1413.      already been lost, but it is the best approach possible given an
  1414.      already sampled image. More details are provided below.
  1415.  
  1416.      Both sharpening and fuzzing are examples of filtering. Even more
  1417.      specifically, they can be both be accomplished with filters which
  1418.      are linear and shift invariant. A crude way to sharpen along a row
  1419.      (or column) is to set output pixel B[n] to the difference of input
  1420.      pixels, A[n]-A[n-1]. A similarly crude way to fuzz is to set B[n]
  1421.      to the average of input pixels, 1/2*A[n]+1/2*A[n-1]. In each case
  1422.      the output is a weighted sum of input pixels, a "convolution". One
  1423.      important characteristic of such filters is that a sinusoid going
  1424.      in produces a sinusoid coming out, one of the same frequency. Thus
  1425.      the Fourier transform, which decomposes a signal into sinusoids of
  1426.      various frequencies, is the key to analysis of these filters. The
  1427.      simplest (and most efficient) way to handle the two dimensions of
  1428.      images is to operate on first the rows then the columns (or vice
  1429.      versa). Fourier transforms and many filters allow this separation.
  1430.  
  1431.      A filter is linear if it satisfies two simple relations between the
  1432.      input and output: scaling the input by a factor scales the output
  1433.      by the same factor, and the sum of two inputs gives the sum of the
  1434.      two outputs. A filter is shift invariant if shifting the input up,
  1435.      down, left, or right merely shifts the output the same way. When a
  1436.      filter is both linear and shift invariant, it can be implemented as
  1437.      a convolution, a weighted sum. If you find the output of the filter
  1438.      when the input is a single pixel with value one in a sea of zeros,
  1439.      you will know all the weights. This output is the impulse response
  1440.      of the filter. The Fourier transform of the impulse response gives
  1441.      the frequency response of the filter. The pattern of weights read
  1442.      off from the impulse response gives the filter kernel, which will
  1443.      usually be displayed (for image filters) as a 2D stencil array, and
  1444.      it is almost always symmetric around the center. For example, the
  1445.      following filter, approximating a Laplacian (and used for detecting
  1446.      edges), is centered on the negative value.
  1447.              1/6     4/6     1/6
  1448.              4/6   -20/6     4/6
  1449.              1/6     4/6     1/6
  1450.      The symmetry allows a streamlined implementation. Suppose the input
  1451.      image is in A, and the output is to go into B. Then compute
  1452.          B[i][j] = (A[i-1][j-1]+A[i-1][j+1]+A[i+1][j-1]+A[i+1][j+1]
  1453.                     +4.0*(A[i-1][j]+A[i][j-1]+A[i][j+1]+A[i+1][j])
  1454.                     -20.0*A[i][j])/6.0
  1455.  
  1456.      Ideal blurring is uniform in all directions, in other words it has
  1457.      circular symmetry. Gaussian blurs are popular, but the obvious code
  1458.      is slow for wide blurs. A cheap alternative is the following filter
  1459.      (written for rows, but then applied to the columns as well).
  1460.          B[i][j] = ((A[i][j]*2+A[i][j-1]+A[i][j+1])*4
  1461.                      +A[i][j-1]+A[i][j+1]-A[i][j-3]-A[i][j+3])/16
  1462.      For sharpening, subtract the results from the original image, which
  1463.      is equivalent to using the following.
  1464.          B[i][j] = ((A[i][j]*2-A[i][j-1]-A[i][j+1])*4
  1465.                      -A[i][j-1]-A[i][j+1]+A[i][j-3]+A[i][j+3])/16
  1466.      Credit for this filter goes to Ken Turkowski and Steve Gabriel.
  1467.  
  1468.      Reconstruction is impossible without some assumptions, and because
  1469.      of the importance of sinusoids in filtering it is traditional to
  1470.      assume the continuous image is made of sinusoids mixed together.
  1471.      That makes more sense for sounds, where signal processing began,
  1472.      than it does for images, especially computer images of character
  1473.      shapes, sharp surface features, and halftoned shading. As pointed
  1474.      out above, often image values cannot be negative, unlike sinusoids.
  1475.      Also, real world images contain noise. The best noise suppressors
  1476.      (and edge detectors) are, ironically, nonlinear filters.
  1477.  
  1478.      The simplest way to double the size of an image is to use each of
  1479.      the original pixels twice in its row and in its column. For much
  1480.      better results, try this instead. Put zeros between the original
  1481.      pixels, then use the blurring filter given a moment ago. But you
  1482.      might want to divide by 8 instead of 16 (since the zeros will dim
  1483.      the image otherwise). To instead shrink the image by half (in both
  1484.      vertical and horizontal), first apply the filter (dividing by 16),
  1485.      then throw away every other pixel. Notice that there are obvious
  1486.      optimizations involving arithmetic with powers of two, zeros which
  1487.      are in known locations, and pixels which will be discarded.
  1488.  
  1489. ----------------------------------------------------------------------
  1490. Subject 3.06: How do I map a texture on to a shape?
  1491.  
  1492.     Paul S. Heckbert, "Survey of Texture Mapping", IEEE Computer
  1493.     Graphics and Applications V6, #11, Nov. 1986, pp 56-67 revised
  1494.     from Graphics Interface '86 version
  1495.  
  1496.     Eric A. Bier and Kenneth R. Sloan, Jr., "Two-Part Texture
  1497.     Mappings", IEEE Computer Graphics and Applications V6 #9, Sept.
  1498.     1986, pp 40-53 (projection parameterizations)
  1499.  
  1500.  
  1501. ----------------------------------------------------------------------
  1502. Subject 3.07: How do I detect a 'corner' in a collection of points?
  1503.  
  1504.     [Currently empty entry.]
  1505.  
  1506. ----------------------------------------------------------------------
  1507. Subject 3.08: Where do I get source to display (raster font format)?
  1508.  
  1509.     ftp://oak.oakland.edu/SimTel/msdos/graphics
  1510.     See also James Murray's graphics file formats FAQ:
  1511.         http://www.ora.com/centers/gff/gff-faq/index.htm
  1512.  
  1513. ----------------------------------------------------------------------
  1514. Subject 3.09: What is morphing/how is it done?
  1515.  
  1516.     Morphing is the name that has come to be applied to the technique
  1517.     ILM used in the movie "Willow", where one object changes into
  1518.     another by changing both its shape and picture detail. It was a
  1519.     2D image manipulation, and has been done in different ways. The
  1520.     first method published was by Thad Beier at PDI. Michael Jackson
  1521.     famously used morphing in his music videos, notably "Black or
  1522.     White". The word is now used more generally.
  1523.  
  1524.     For more, see [Anderson], Chapter 3, page 59-90, and Beier's
  1525.     http://www.hammerhead.com/thad/morph.html
  1526.  
  1527. ----------------------------------------------------------------------
  1528. Subject 3.10: How do I quickly draw a filled triangle?
  1529.  
  1530.     The easiest way is to render each line separately into an edge
  1531.     buffer. This buffer is a structure which looks like this (in C):
  1532.  
  1533.         struct { int xmin, xmax; } edgebuffer[YDIM];
  1534.  
  1535.     There is one entry for each scan line on the screen, and each
  1536.     entry is to be interpreted as a horizontal line to be drawn from
  1537.     xmin to xmax.
  1538.  
  1539.     Since most people who ask this question are trying to write fast
  1540.     games on the PC, I'll tell you where to find code.  Look at:
  1541.  
  1542.         ftp::/ftp.uwp.edu/pub/msdos/demos/programming/source
  1543.         ftp::/ftp.luth.se/pub/msdos/demos (Sweden)
  1544.         ftp::/NCTUCCCA.edu.tw:/PC/uwp/demos
  1545.         http://www.wit.com:/mirrors/uwp/pub/msdos/demos
  1546.         ftp::/ftp.cdrom.com:/demos
  1547.  
  1548.     See also Subject 3.03, which describes methods for filling polygons.
  1549.  
  1550. ----------------------------------------------------------------------
  1551. Subject 3.11: 3D Noise functions and turbulence in Solid texturing.
  1552.  
  1553.     See:
  1554.         ftp://gondwana.ecr.mu.oz.au/pub/siggraph92_C23.shar.gz
  1555.         ftp://ftp.cis.ohio-state.edu/pub/siggraph92/siggraph92_C23.shar
  1556.  
  1557.     In it there are implementations of Perlin's noise and turbulence
  1558.     functions, (By the man himself) as well as Lewis' sparse
  1559.     convolution noise function (by D. Peachey) There is also some of
  1560.     other stuff in there (Musgrave's Earth texture functions, and some
  1561.     stuff on animating gases by Ebert).
  1562.  
  1563.     SPD (Standard Procedural Databases) package:
  1564.         ftp://avalon.chinalake.navy.mil/utils/SPD/SPD33f4.tar.Z
  1565.         ftp://avalon.chinalake.navy.mil/utils/SPD/spd33f4.zip.
  1566.         Now moved to http://www.viewpoint.com/
  1567.  
  1568.  
  1569.     References:
  1570.  
  1571.     [Ebert]
  1572.     Noise, Hypertexture, Antialiasing and Gesture, (Ken Perlin) in
  1573.     Chapter 6, (p.193-), The disk accompanying the book is available
  1574.     from
  1575.         ftp://archive.cs.umbc.edu/pub/texture.
  1576.  
  1577.     For more info on this text/code see:
  1578.         http://www.cs.umbc.edu/~ebert/book/book.html
  1579.  
  1580.     For examples from a current course based on this book, see:
  1581.         http://www.seas.gwu.edu/graphics/ProcTexCourse/
  1582.     Linke broken 21Jan03; will remove eventually if not fixed.
  1583.  
  1584.  
  1585.     [Watt:Animation]
  1586.     Three-dimensional Nocie, Chapter 7.2.1
  1587.     Simulating turbulance, Chapter 7.2.2
  1588.  
  1589. ----------------------------------------------------------------------
  1590. Subject 3.12: How do I generate realistic sythetic textures?
  1591.  
  1592.     For fractal mountains, trees and sea-shells:
  1593.  
  1594.       SPD (Standard Procedural Databases) package:
  1595.         ftp://avalon.chinalake.navy.mil/utils/SPD/SPD33f4.tar.Z
  1596.         ftp://avalon.chinalake.navy.mil/utils/SPD/spd33f4.zip.
  1597.         Now moved to http://www.viewpoint.com/
  1598.  
  1599.  
  1600.  
  1601.     Reaction-Diffusion Algorithms:
  1602.       For an illustartion of the parameter space of a reaction
  1603.       diffusion system, check out the Xmorphia page at
  1604.         http://www.ccsf.caltech.edu/ismap/image.html
  1605.  
  1606.  
  1607.     References:
  1608.  
  1609.     [Ebert]
  1610.     Entire book devoted to this subject, with RenderMan(TM) and C code.
  1611.  
  1612.     [Watt:Animation]
  1613.     Procedural texture mapping and modelling, Chapter 7
  1614.  
  1615.     "Generating Textures on Arbitrary Surfaces Using Reaction-Diffusion"
  1616.      Greg Turk,   Computer Graphics, Vol. 25, No. 4, pp. 289-298
  1617.      July 1991 (SIGGRAPH '91)
  1618.         http://www.cs.unc.edu:80/~turk/reaction_diffusion/reaction_diffusion.html
  1619.  
  1620.     A list of procedural texture synthesis related web pages
  1621.          http://www.threedgraphics.com/pixelloom/tex_synth.html
  1622.  
  1623. ----------------------------------------------------------------------
  1624. Subject 3.13: How do I convert between color models (RGB, HLS, CMYK, CIE etc)?
  1625.  
  1626.     References:
  1627.     [Watt:3D] pp. 313-354
  1628.     [Foley]   pp. 563-603
  1629.  
  1630. ----------------------------------------------------------------------
  1631. Subject 3.14: How is "GIF" pronounced?
  1632.  
  1633.    "GIF" is an acronymn for "Graphics Interchange Format."  Despite the
  1634.    hard "G" in "Graphics," GIF is pronounced "JIF."  Although we don't
  1635.    have a direct quote from the official CompuServe specification
  1636.    released June 1987, here is a quote from related CompuServe documentation,
  1637.    for CompuShow, a DOS-based image viewer used shortly thereafter: 
  1638.      "The GIF (Graphics Interchange Format), pronounced "JIF", was 
  1639.       designed by CompuServe ..."
  1640.    We also have a report that the principal author of the GIF spec,
  1641.    Steve Wilhite, says "it's pronounced JIF (like the peanut butter."
  1642.  
  1643.    See also  
  1644.      http://www.60-seconds.com/articles/86.html
  1645.  
  1646. ----------------------------------------------------------------------
  1647. Section 4. Curve Computations
  1648. ----------------------------------------------------------------------
  1649. Subject 4.01: How do I generate a Bezier curve that is parallel to another Bezier?
  1650.  
  1651.     You can't.  The only case where this is possible is when the
  1652.     Bezier can be represented by a straight line.  And then the
  1653.     parallel 'Bezier' can also be represented by a straight line.
  1654.  
  1655.     The situation is different for the broader class of rational
  1656.     Bezier curves.  For example, these can represent circular arcs,
  1657.     and a parallel offset is just a concentric circular arc, also
  1658.     representable as a rational Bezier.
  1659.  
  1660.  
  1661. ----------------------------------------------------------------------
  1662. Subject 4.02: How do I split a Bezier at a specific value for t?
  1663.  
  1664.     A Bezier curve is a parametric polynomial function from the
  1665.     interval [0..1] to a space, usually 2D or 3D.  Common Bezier
  1666.     curves use cubic polynomials, so have the form
  1667.  
  1668.         f(t) = a3 t^3 + a2 t^2 + a1 t + a0,
  1669.  
  1670.     where the coefficients are points in 3D. Blossoming converts this
  1671.     polynomial to a more helpful form.  Let s = 1-t, and think of
  1672.     tri-linear interpolation:
  1673.  
  1674.         F([s0,t0],[s1,t1],[s2,t2]) =
  1675.             s0(s1(s2 c30 + t2 c21) + t1(s2 c21 + t2 c12)) +
  1676.             t0(s1(s2 c21 + t2 c12) + t1(s2 c12 + t2 c03))
  1677.             =
  1678.             c30(s0 s1 s2) +
  1679.             c21(s0 s1 t2 + s0 t1 s2 + t0 s1 s2) +
  1680.             c12(s0 t1 t2 + t0 s1 t2 + t0 t1 s2) +
  1681.             c03(t0 t1 t2).
  1682.  
  1683.     The data points c30, c21, c12, and c03 have been used in such a
  1684.     way as to make the resulting function give the same value if any
  1685.     two arguments, say [s0,t0] and [s2,t2], are swapped. When [1-t,t]
  1686.     is used for all three arguments, the result is the cubic Bezier
  1687.     curve with those data points controlling it:
  1688.  
  1689.           f(t) = F([1-t,t],[1-t,t],[1-t,t])
  1690.                =  (1-t)^3 c30 +
  1691.                  3(1-t)^2 t c21 +
  1692.                  3(1-t) t^2 c12 +
  1693.                  t^3 c03.
  1694.  
  1695.     Notice that
  1696.            F([1,0],[1,0],[1,0]) = c30,
  1697.            F([1,0],[1,0],[0,1]) = c21,
  1698.            F([1,0],[0,1],[0,1]) = c12, _
  1699.            F([0,1],[0,1],[0,1]) = c03.
  1700.  
  1701.     In other words, cij is obtained by giving F argument t's i of
  1702.     which are 0 and j of which are 1. Since F is a linear polynomial
  1703.     in each argument, we can find f(t) using a series of simple steps.
  1704.     Begin with
  1705.  
  1706.         f000 = c30, f001 = c21, f011 = c12, f111 = c03.
  1707.  
  1708.     Then compute in succession:
  1709.  
  1710.         f00t = s f000 + t f001, f01t = s f001 + t f011,
  1711.         f11t = s f011 + t f111,
  1712.         f0tt = s f00t + t f01t, f1tt = s f01t + t f11t,
  1713.         fttt = s f0tt + t f1tt.
  1714.  
  1715.     This is the de Casteljau algorithm for computing f(t) = fttt.
  1716.  
  1717.     It also has split the curve for the intervals [0..t] and [t..1].
  1718.     The control points for the first interval are f000, f00t, f0tt,
  1719.     fttt, while those for the second interval are fttt, f1tt, f11t,
  1720.     f111.
  1721.  
  1722.     If you evaluate 3 F([1-t,t],[1-t,t],[-1,1]) you will get the
  1723.     derivate of f at t. Similarly, 2*3 F([1-t,t],[-1,1],[-1,1]) gives
  1724.     the second derivative of f at t, and finally using 1*2*3
  1725.     F([-1,1],[-1,1],[-1,1]) gives the third derivative.
  1726.  
  1727.     This procedure is easily generalized to different degrees,
  1728.     triangular patches, and B-spline curves.
  1729.  
  1730.  
  1731. ----------------------------------------------------------------------
  1732. Subject 4.03: How do I find a t value at a specific point on a Bezier?
  1733.  
  1734.     In general, you'll need to find t closest to your search point.
  1735.     There are a number of ways you can do this. Look at [Gems I, 607],
  1736.     there's a chapter on finding the nearest point on the Bezier
  1737.     curve.  In my experience, digitizing the Bezier curve is an
  1738.     acceptable method. You can also try recursively subdividing the
  1739.     curve, see if you point is in the convex hull of the control
  1740.     points, and checking is the control points are close enough to a
  1741.     linear line segment and find the nearest point on the line
  1742.     segment, using linear interpolation and keeping track of the
  1743.     subdivision level, you'll be able to find t.
  1744.  
  1745.  
  1746. ----------------------------------------------------------------------
  1747. Subject 4.04: How do I fit a Bezier curve to a circle?
  1748.  
  1749.     Interestingly enough, Bezier curves can approximate a circle but
  1750.     not perfectly fit a circle.
  1751.     A common approximation is to use four beziers to model a circle, each
  1752.     with control points a distance d=r*4*(sqrt(2)-1)/3 from the end points
  1753.     (where r is the circle radius), and in a direction tangent to the 
  1754.     circle at the end points.  This will ensure the mid-points of the 
  1755.     Beziers are on the circle, and that the first derivative is continuous.  
  1756.     The radial error in this approximation will be about 0.0273% of the 
  1757.     circle's radius.
  1758.  
  1759.     Michael Goldapp, "Approximation of circular arcs by cubic
  1760.     polynomials" Computer Aided Geometric Design (#8 1991 pp.227-238)
  1761.  
  1762.     Tor Dokken and Morten Daehlen, "Good Approximations of circles by
  1763.     curvature-continuous Bezier curves" Computer Aided Geometric
  1764.     Design (#7 1990 pp. 33-41).
  1765.  
  1766.     See also http://www.whizkidtech.net/bezier/circle/ .
  1767.  
  1768.  
  1769. ----------------------------------------------------------------------
  1770. Section 5. 3D computations
  1771. ----------------------------------------------------------------------
  1772. Subject 5.01: How do I rotate a 3D point?
  1773.  
  1774.     Let's assume you want to rotate vectors around the origin of your
  1775.     coordinate system. (If you want to rotate around some other point,
  1776.     subtract its coordinates from the point you are rotating, do the
  1777.     rotation, and then add back what you subtracted.) In 3D, you need
  1778.     not only an angle, but also an axis. (In higher dimensions it gets
  1779.     much worse, very quickly.)  Actually, you need 3 independent
  1780.     numbers, and these come in a variety of flavors.  The flavor I
  1781.     recommend is unit quaternions: 4 numbers that square and add up to
  1782.     +1.  You can write these as [(x,y,z),w], with 4 real numbers, or
  1783.     [v,w], with v, a 3D vector pointing along the axis. The concept
  1784.     of an axis is unique to 3D. It is a line through the origin
  1785.     containing all the points which do not move during the rotation.
  1786.     So we know if we are turning forwards or back, we use a vector
  1787.     pointing out along the line. Suppose you want to use unit vector u
  1788.     as your axis, and rotate by 2t degrees.  (Yes, that's twice t.)
  1789.     Make a quaternion [u sin t, cos t]. You can use the quaternion --
  1790.     call it q -- directly on a vector v with quaternion
  1791.     multiplication, q v q^-1, or just convert the quaternion to a 3x3
  1792.     matrix M. If the components of q are {(x,y,z),w], then you want
  1793.     the matrix
  1794.  
  1795.         M = {{1-2(yy+zz),  2(xy-wz),  2(xz+wy)},
  1796.              {  2(xy+wz),1-2(xx+zz),  2(yz-wx)},
  1797.              {  2(xz-wy),  2(yz+wx),1-2(xx+yy)}}.
  1798.  
  1799.     Rotations, translations, and much more are explained in all basic
  1800.     computer graphics texts.  Quaternions are covered briefly in
  1801.     [Foley], and more extensively in several Graphics Gems, and the
  1802.     SIGGRAPH 85 proceedings.
  1803.  
  1804.     /* The following routine converts an angle and a unit axis vector
  1805.         * to a matrix, returning the corresponding unit quaternion at no
  1806.         * extra cost. It is written in such a way as to allow both fixed
  1807.         * point and floating point versions to be created by appropriate
  1808.         * definitions of FPOINT, ANGLE, VECTOR, QUAT, MATRIX, MUL, HALF,
  1809.         * TWICE, COS, SIN, ONE, and ZERO.
  1810.         * The following is an example of floating point definitions.
  1811.         #define FPOINT double
  1812.         #define ANGLE FPOINT
  1813.         #define VECTOR QUAT
  1814.         typedef struct {FPOINT x,y,z,w;} QUAT;
  1815.         enum Indices {X,Y,Z,W};
  1816.         typedef FPOINT MATRIX[4][4];
  1817.         #define MUL(a,b) ((a)*(b))
  1818.         #define HALF(a) ((a)*0.5)
  1819.         #define TWICE(a) ((a)*2.0)
  1820.         #define COS cos
  1821.         #define SIN sin
  1822.         #define ONE 1.0
  1823.         #define ZERO 0.0
  1824.     */
  1825.     QUAT MatrixFromAxisAngle(VECTOR axis, ANGLE theta, MATRIX m)
  1826.     {
  1827.         QUAT q;
  1828.         ANGLE halfTheta = HALF(theta);
  1829.         FPOINT cosHalfTheta = COS(halfTheta);
  1830.         FPOINT sinHalfTheta = SIN(halfTheta);
  1831.         FPOINT xs, ys, zs, wx, wy, wz, xx, xy, xz, yy, yz, zz;
  1832.         q.x = MUL(axis.x,sinHalfTheta);
  1833.         q.y = MUL(axis.y,sinHalfTheta);
  1834.         q.z = MUL(axis.z,sinHalfTheta);
  1835.         q.w = cosHalfTheta;
  1836.         xs = TWICE(q.x);  ys = TWICE(q.y);  zs = TWICE(q.z);
  1837.         wx = MUL(q.w,xs); wy = MUL(q.w,ys); wz = MUL(q.w,zs);
  1838.         xx = MUL(q.x,xs); xy = MUL(q.x,ys); xz = MUL(q.x,zs);
  1839.         yy = MUL(q.y,ys); yz = MUL(q.y,zs); zz = MUL(q.z,zs);
  1840.         m[X][X] = ONE - (yy + zz); m[X][Y] = xy - wz; m[X][Z] = xz + wy;
  1841.         m[Y][X] = xy + wz; m[Y][Y] = ONE - (xx + zz); m[Y][Z] = yz - wx;
  1842.         m[Z][X] = xz - wy; m[Z][Y] = yz + wx; m[Z][Z] = ONE - (xx + yy);
  1843.         /* Fill in remainder of 4x4 homogeneous transform matrix. */
  1844.         m[W][X] = m[W][Y] = m[W][Z] = m[X][W] = m[Y][W] = m[Z][W] = ZERO;
  1845.         m[W][W] = ONE;
  1846.         return (q);
  1847.     }
  1848.     /* The routine just given, MatrixFromAxisAngle, performs rotation about
  1849.         * an axis passing through the origin, so only a unit vector was needed
  1850.         * in addition to the angle. To rotate about an axis not containing the
  1851.         * origin, a point on the axis is also needed, as in the following. For
  1852.         * mathematical purity, the type POINT is used, but may be defined as:
  1853.         #define POINT VECTOR
  1854.     */
  1855.     QUAT MatrixFromAnyAxisAngle(POINT o, VECTOR axis, ANGLE theta, MATRIX m)
  1856.     {
  1857.         QUAT q;
  1858.         q = MatrixFromAxisAngle(axis,theta,m);
  1859.         m[X][W] = o.x-(MUL(m[X][X],o.x)+MUL(m[X][Y],o.y)+MUL(m[X][Z],o.z));
  1860.         m[Y][W] = o.y-(MUL(m[Y][X],o.x)+MUL(m[Y][Y],o.y)+MUL(m[Y][Z],o.z));
  1861.         m[Z][W] = o.x-(MUL(m[Z][X],o.x)+MUL(m[Z][Y],o.y)+MUL(m[Z][Z],o.z));
  1862.         return (q);
  1863.     }
  1864.     /* An axis can also be specified by a pair of points, with the direction
  1865.         * along the line assumed from the ordering of the points. Although both
  1866.         * the previous routines assumed the axis vector was unit length without
  1867.         * checking, this routine must cope with a more delicate possibility. If
  1868.         * the two points are identical, or even nearly so, the axis is unknown.
  1869.         * For now the auxiliary routine which makes a unit vector ignores that.
  1870.         * It needs definitions like the following for floating point.
  1871.         #define SQRT sqrt
  1872.         #define RECIPROCAL(a) (1.0/(a))
  1873.     */
  1874.     VECTOR Normalize(VECTOR v)
  1875.     {
  1876.         VECTOR u;
  1877.         FPOINT norm = MUL(v.x,v.x)+MUL(v.y,v.y)+MUL(v.z,v.z);
  1878.         /* Better to test for (near-)zero norm before taking reciprocal. */
  1879.         FPOINT scl = RECIPROCAL(SQRT(norm));
  1880.         u.x = MUL(v.x,scl); u.y = MUL(v.y,scl); u.z = MUL(v.z,scl);
  1881.         return (u);
  1882.     }
  1883.     QUAT MatrixFromPointsAngle(POINT o, POINT p, ANGLE theta, MATRIX m)
  1884.     {
  1885.         QUAT q;
  1886.         VECTOR diff, axis;
  1887.         diff.x = p.x-o.x; diff.y = p.y-o.y; diff.z = p.z-o.z;
  1888.         axis = Normalize(diff);
  1889.         return (MatrixFromAnyAxisAngle(o,axis,theta,m));
  1890.     }
  1891.  
  1892. ----------------------------------------------------------------------
  1893. Subject 5.02: What is ARCBALL and where is the source?
  1894.  
  1895.     Arcball is a general purpose 3D rotation controller described by
  1896.     Ken Shoemake in the Graphics Interface '92 Proceedings.  It
  1897.     features good behavior, easy implementation, cheap execution, and
  1898.     optional axis constraints. A Macintosh demo and electronic version
  1899.     of the original paper (Microsoft Word format) may be obtained from
  1900.     ftp::/ftp.cis.upenn.edu/pub/graphics.
  1901.  
  1902.     Complete source code appears in Graphics Gems IV pp. 175-192. A
  1903.     fairly complete sketch of the code appeared in the original
  1904.     article, in Graphics Interface 92 Proceedings, available from
  1905.     Morgan Kaufmann.
  1906.  
  1907.     The original arcball code was written for IRIS GL. A translation
  1908.     into OpenGL/GLUT, and for IRIS Performer, may be found at:
  1909.         http://cs.anu.edu.au/people/Hugh.Fisher/3dstuff/
  1910.  
  1911.  
  1912. ----------------------------------------------------------------------
  1913. Subject 5.03: How do I clip a polygon against a rectangle?
  1914.  
  1915.     This is the Sutherland-Hodgman algorithm, an old favorite that
  1916.     should be covered in any textbook. See the selected list below.
  1917.     According to Vatti (q.v.)  "This
  1918.     [algorithm] produces degenerate edges in certain concave / self
  1919.     intersecting polygons that need to be removed as a special
  1920.     extension to the main algorithm" (though this is not a problem if
  1921.     all you are doing with the results is scan converting them.)
  1922.  
  1923.     It should be noted that the Sutherland-Hodgman algorithm
  1924.     may be used to clip a polygon against any convex polygon.
  1925.     Cf. also Subject 5.04.
  1926.  
  1927.     [Foley, van Dam]: Section 3.14.1 (pp 124 - 126)
  1928.     [Hearn]: Section 6-8, pp 237 - 242 (with actual C code!)
  1929.     See also
  1930.       http://www.csclub.uwaterloo.ca/u/mpslager/articles/sutherland/wr.html 
  1931.  
  1932. ----------------------------------------------------------------------
  1933. |Subject 5.04: How do I clip a polygon against another polygon?
  1934.  
  1935.     Klamer Schutte, klamer@ph.tn.tudelft.nl has developed and implemented
  1936.     some code in C++ to perform clipping of two possibly concave 2D
  1937.     polygons. A description can be found at:
  1938.     http://www.ph.tn.tudelft.nl:/People/klamer/clippoly_entry.html
  1939.     To compile the source code you will need a C++ compiler with templates,
  1940.     such as g++. The source code is available at:
  1941.     ftp://ftp.ph.tn.tudelft.nl/pub/klamer/clippoly.tar.gz
  1942. |   See also http://home.attbi.com/~msleonov/pbcomp.html, which extends
  1943.     the above to permit holes.
  1944.  
  1945.     Alan Murta released a polygon clipper library (in C) which uses a 
  1946.     modified version of the Vatti algorithm:
  1947.       http://www.cs.man.ac.uk/aig/staff/alan/software/index.html
  1948.  
  1949.     References:
  1950.  
  1951.     Weiler, K. "Polygon Comparison Using a Graph Representation", SIGGRAPH 80
  1952.     pg. 10-18
  1953.  
  1954.     Vatti, Bala R. "A Generic Solution to Polygon Clipping",
  1955.     Communications of the ACM, July 1992, Vol 35, No. 7, pg. 57-63
  1956.  
  1957. ----------------------------------------------------------------------
  1958. Subject 5.05: How do I find the intersection of a line and a plane?
  1959.  
  1960.     If the plane is defined as:
  1961.  
  1962.         a*x + b*y + c*z + d = 0
  1963.  
  1964.     and the line is defined as:
  1965.  
  1966.         x = x1 + (x2 - x1)*t = x1 + i*t
  1967.         y = y1 + (y2 - y1)*t = y1 + j*t
  1968.         z = z1 + (z2 - z1)*t = z1 + k*t
  1969.  
  1970.     Then just substitute these into the plane equation. You end up
  1971.     with:
  1972.  
  1973.         t = - (a*x1 + b*y1 + c*z1 + d)/(a*i + b*j + c*k)
  1974.  
  1975.     When the denominator is zero, the line is contained in the plane 
  1976.     if the numerator is also zero (the point at t=0 satisfies the
  1977.     plane equation), otherwise the line is parallel to the plane.
  1978.  
  1979.  
  1980. ----------------------------------------------------------------------
  1981. Subject 5.06: How do I determine the intersection between a ray and a triangle?
  1982.  
  1983.     First find the intersection between the ray and the plane in which
  1984.     the triangle is situated. Then see if the point of intersection is
  1985.     inside the triangle.
  1986.     Details may be found in [O'Rourke (C)] pp.226-238, whose code is at
  1987.        http://cs.smith.edu/~orourke/ .
  1988.     Efficient code complete with statistical tests is described in the Mo:ller-
  1989.     Trumbore paper in J. Graphics Tools (C code downloadable from there):
  1990.        http://www.acm.org/jgt/papers/MollerTrumbore97/
  1991.     See also the full paper:
  1992.        http://www.Graphics.Cornell.EDU/pubs/1997/MT97.html
  1993.     See also the "3D Object Intersection" page, described in Subject 0.05.
  1994.  
  1995.  
  1996. ----------------------------------------------------------------------
  1997. Subject 5.07: How do I determine the intersection between a ray and a sphere
  1998.  
  1999.     Given a ray defined as:
  2000.  
  2001.         x = x1 + (x2 - x1)*t = x1 + i*t
  2002.         y = y1 + (y2 - y1)*t = y1 + j*t
  2003.         z = z1 + (z2 - z1)*t = z1 + k*t
  2004.  
  2005.     and a sphere defined as:
  2006.  
  2007.         (x - l)**2 + (y - m)**2 + (z - n)**2 = r**2
  2008.  
  2009.     Substituting in gives the quadratic equation:
  2010.  
  2011.         a*t**2 + b*t + c = 0
  2012.  
  2013.     where:
  2014.  
  2015.         a = i**2 + j**2 + k**2
  2016.         b = 2*i*(x1 - l) + 2*j*(y1 - m) + 2*k*(z1 - n)
  2017.         c = (x1-l)**2 + (y1-m)**2 + (z1-n)**2 - r**2;
  2018.  
  2019.     If the discriminant of this equation is less than 0, the line does
  2020.     not intersect the sphere. If it is zero, the line is tangential to
  2021.     the sphere and if it is greater than zero it intersects at two
  2022.     points. Solving the equation and substituting the values of t into
  2023.     the ray equation will give you the points.
  2024.  
  2025.     Reference:
  2026.     [Glassner:RayTracing]
  2027.  
  2028.     See also the "3D Object Intersection" page, described in Subject 0.05.
  2029.  
  2030. ----------------------------------------------------------------------
  2031. Subject 5.08: How do I find the intersection of a ray and a Bezier surface?
  2032.  
  2033.     Joy I.K. and Bhetanabhotla M.N., "Ray tracing parametric surfaces
  2034.     utilizing numeric techniques and ray coherence", Computer
  2035.     Graphics, 16, 1986, 279-286
  2036.  
  2037.     Joy and Bhetanabhotla is only one approach of three major method
  2038.     classes: tessellation, direct computation, and a hybrid of the
  2039.     two.  Tessellation is relatively easy (you intersect the polygons
  2040.     making up the tessellation) and has no numerical problems, but can
  2041.     be inaccurate; direct is cheap on memory, but very expensive
  2042.     computationally, and can (and usually does) suffer from precision
  2043.     problems within the root solver; hybrids try to blend the two.
  2044.  
  2045.     Reference:
  2046.     [Glassner:RayTracing]
  2047.  
  2048.     See also the "3D Object Intersection" page, described in Subject 0.05.
  2049.  
  2050. ----------------------------------------------------------------------
  2051. Subject 5.09: How do I ray trace caustics?
  2052.  
  2053.     See the work of Henrik Wann Jensen at 
  2054.        http://graphics.ucsd.edu/~henrik/ 
  2055.  
  2056.     @inproceedings{j-rcnls-96
  2057.     , author =      "Henrik Wann Jensen"
  2058.     , title =       "Rendering Caustics on Non-Lambertian Surfaces"
  2059.     , booktitle =   "Proc. Graphics Interface '96"
  2060.     , pages =       "116--121"
  2061.     , location =    "Toronto"
  2062.     , year =         1996
  2063.     }
  2064.  
  2065.     Metropolis Light Transport handles this phenomenon well:
  2066.       http://www-graphics.stanford.edu/papers/metro/
  2067.  
  2068.     Bidirectional path tracing also handles caustics.
  2069.       http://graphics.stanford.EDU/papers/veach_thesis/ (Chapter 10)
  2070.       http://www.graphics.cornell.edu/~eric/thesis/
  2071.  
  2072.  
  2073.     Some older references:
  2074.  
  2075.     An expensive answer:
  2076.     @InProceedings{mitchell-1992-illumination,
  2077.       author         = "Don P. Mitchell and Pat Hanrahan",
  2078.       title          = "Illumination From Curved Reflectors",
  2079.       year           = "1992",
  2080.       month          = "July",
  2081.       volume         = "26",
  2082.       booktitle      = "Computer Graphics (SIGGRAPH '92 Proceedings)",
  2083.       pages          = "283--291",
  2084.       keywords       = "caustics, interval arithmetic, ray tracing",
  2085.       editor         = "Edwin E. Catmull",
  2086.     }
  2087.  
  2088.     A cheat:
  2089.     @Article{inakage-1986-caustics,
  2090.       author         = "Masa Inakage",
  2091.       title          = "Caustics and Specular Reflection Models for
  2092.                         Spherical Objects and Lenses ",
  2093.       pages          = "379--383",
  2094.       journal        = "The Visual Computer",
  2095.       volume         = "2",
  2096.       number         = "6",
  2097.       year           = "1986",
  2098.       keywords       = "ray tracing effects",
  2099.     }
  2100.  
  2101.     Very specialized:
  2102.     @Article{yuan-1988-gemstone,
  2103.       author         = "Ying Yuan and Tosiyasu L. Kunii and Naota
  2104.                         Inamato and Lining Sun ",
  2105.       title          = "Gemstone Fire: Adaptive Dispersive Ray Tracing
  2106.                         of Polyhedrons",
  2107.       year           = "1988",
  2108.       month          = "November",
  2109.       journal        = "The Visual Computer",
  2110.       volume         = "4",
  2111.       number         = "5",
  2112.       pages          = "259--70",
  2113.       keywords       = "caustics",
  2114.     }
  2115.  
  2116.  
  2117. ----------------------------------------------------------------------
  2118. Subject 5.10: What is the marching cubes algorithm?
  2119.  
  2120.     The marching cubes algorithm is used in volume rendering to
  2121.     construct an isosurface from a 3D field of values.
  2122.  
  2123.     The 2D analog would be to take an image, and for each pixel, set
  2124.     it to black if the value is below some threshold, and set it to
  2125.     white if it's above the threshold.  Then smooth the jagged black
  2126.     outlines by skinning them with lines.
  2127.  
  2128.     The marching cubes algorithm tests the corner of each cube (or
  2129.     voxel) in the scalar field as being either above or below a given
  2130.     threshold.  This yields a collection of boxes with classified
  2131.     corners.  Since there are eight corners with one of two states,
  2132.     there are 256 different possible combinations for each cube.
  2133.     Then, for each cube, you replace the cube with a surface that
  2134.     meets the classification of the cube.  For example, the following
  2135.     are some 2D examples showing the cubes and their associated
  2136.     surface.
  2137.  
  2138.         - ----- +       - ----- -       - ----- +       - ----- +
  2139.         |:::'   |       |:::::::|       |::::   |       |   ':::|
  2140.         |:'     |       |:::::::|       |::::   |       |.    ':|
  2141.         |       |       |       |       |::::   |       |::.    |
  2142.         + ----- +       + ----- +       - ----- +       + ----- -
  2143.  
  2144.     The result of the marching cubes algorithm is a smooth surface
  2145.     that approximates the isosurface that is constant along a given
  2146.     threshold. This is useful for displaying a volume of oil in a
  2147.     geological volume, for example.
  2148.  
  2149.     References:
  2150.     "Marching Cubes: A High Resolution 3D Surface  Construction Algorithm",
  2151.     William E. Lorensen and Harvey E. Cline,
  2152.     Computer Graphics (Proceedings of  SIGGRAPH '87), Vol. 21, No. 4, pp. 163-169.
  2153.  
  2154.     [Watt:Animation] pp. 302-305 and 313-321
  2155.     [Schroeder]
  2156.  
  2157.     For alternatives to the (patented; cf. Subj. 5.11) marching cubes algorithm, 
  2158.     see
  2159.        http://www.unchainedgeometry.com/jbloom/papers/index.html
  2160.     under "Implicit Surface Polygonization."
  2161.  
  2162. ----------------------------------------------------------------------
  2163. Subject 5.11: What is the status of the patent on the "marching cubes" algorithm?
  2164.  
  2165.     United States Patent Number: 4,710,876
  2166.     Date of Patent: Dec. 1, 1987
  2167.     Inventors: Harvey E. Cline, William E. Lorensen
  2168.     Assignee: General Electric Company
  2169.     Title: "System and Method for the Display of Surface Structures Contained
  2170.             Within the Interior Region of a Solid Body"
  2171.     Filed: Jun. 5, 1985
  2172.       http://www.delphion.com/
  2173.       Type in "4710876" (w/o commas, w/o quotes) into their search engine.
  2174.  
  2175.     United States Patent Number: 4,885,688
  2176.     Date of Patent: Dec. 5, 1989
  2177.     Inventor: Carl R. Crawford
  2178.     Assignee: General Electric Company
  2179.     Title: "Minimization of Directed Points Generated in Three-Dimensional
  2180.             Dividing Cubes Method"
  2181.     Filed: Nov. 25, 1987
  2182.       Access as above.
  2183.  
  2184.     For alternative, unpatented algorithms, cf. Subj. 5.10.
  2185.  
  2186. ----------------------------------------------------------------------
  2187. Subject 5.12: How do I do a hidden surface test (backface culling) with 3D points?
  2188.  
  2189.     Just define all points of all polygons in clockwise order.
  2190.  
  2191.             c = (x3*((z1*y2)-(y1*z2))+
  2192.                 (y3*((x1*z2)-(z1*x2))+
  2193.                 (z3*((y1*x2)-(x1*y2))+
  2194.  
  2195.     x1,y1,z1, x2,y2,z2, x3,y3,z3 = the first three points of the
  2196.     polygon.
  2197.  
  2198.     If c is positive the polygon is visible. If c is negative the
  2199.     polygon is invisible (or the other way).
  2200.  
  2201.  
  2202. ----------------------------------------------------------------------
  2203. Subject 5.13: Where can I find algorithms for 3D collision detection?
  2204.  
  2205.     Check out "proxima", from Purdue, which is a C++ library for 3D
  2206.     collision detection for arbitrary polyhedra.  It's a nice system;
  2207.     the algorithms are sophisticated, but the code is of modest size.
  2208.  
  2209.         ftp://ftp.cs.purdue.edu/pub/pse/PROX/prox9.1.tar.Z
  2210.  
  2211.     For information about the I_COLLIDE 3D collision detection system
  2212.         http://www.cs.unc.edu/~geom/I_COLLIDE.html
  2213.  
  2214.     "Fast Collision Detection of Moving Convex Polyhedra", Rich Rabbitz,
  2215.     Graphics Gems IV, pages 83-109, includes source in C.
  2216.  
  2217.     SOLID: "a library for collision detection of three-dimensional
  2218.     objects undergoing rigid motion and deformation. SOLID is designed to 
  2219.     be used in interactive 3D graphics applications, and is especially 
  2220.     suited for collision detection of objects and worlds described in VRML.
  2221.     Written in standard C++, compiles under GNU g++ version 2.8.1 and 
  2222.     Visual C++ 5.0."  See:
  2223.         http://www.win.tue.nl/cs/tt/gino/solid/
  2224.  
  2225.     SWIFT++: a C++ library for collision detection, exact and approximate
  2226.     distance computation, and contact determination of three-dimensional
  2227.     polyhedral objects undergoing rigid motion.
  2228.     Some preliminary results indicate that it is faster than I-COLLIDE and 
  2229.     V-CLIP, and more robust than I-COLLIDE.
  2230.        http://www.cs.unc.edu/~geom/SWIFT++
  2231.  
  2232.     ColDet: C++ library for 3D collison detection.  Works on generic 
  2233.     polyhedra, and even polygon soups.   Uses bounding box hierarchies 
  2234.     and triangle intersection tests.   Released as open source under LGPL.
  2235.     Tested on Windows, MacOS, and Linux. http://photoneffect.com/coldet/ .
  2236.  
  2237.     Terdiman's lib, which might need less RAM than the above:
  2238.        http://www.codercorner.com/Opcode.htm
  2239.  
  2240.  
  2241. ----------------------------------------------------------------------
  2242. Subject 5.14: How do I perform basic viewing in 3D?
  2243.  
  2244.     We describe the shape and position of objects using numbers,
  2245.     usually (x,y,z) coordinates. For example, the corners of a cube
  2246.     might be {(0,0,0), (1,0,0), (1,1,0), (0,1,0), (0,0,1), (1,0,1),
  2247.     (1,1,1), (0,1,1)}. A deep understanding of what we are saying with
  2248.     these numbers requires mathematical study, but I will try to keep
  2249.     this simple. At the least, we must understand that we have
  2250.     designated some point in space as the origin--coordinates
  2251.     (0,0,0)--and marked off lines in 3 mutually perpendicular
  2252.     directions using equally spaced units to give us (x,y,z) values.
  2253.     It might be helpful to know if we are using 1 to mean 1 foot, 1
  2254.     meter, or 1 parsec; the numbers alone do not tell us.
  2255.  
  2256.     A picture on a screen is two steps removed from the 3D world it
  2257.     depicts. First, it is a 2D projection; and second, it is a finite
  2258.     resolution approximation to the infinitely precise projection. I
  2259.     will ignore the approximation (sampling) for now. To know what the
  2260.     projection looks like, we need to know where our viewpoint is, and
  2261.     where the plane of the projection is, both in the 3D world. Think
  2262.     of it as looking out a window into a scene. As artists discovered
  2263.     some 500 years ago, each point in the world appears to be at a
  2264.     point on the window. If you move your head or look out a different
  2265.     window, everything changes. When the mathematicians understood
  2266.     what the artists were doing, they invented perspective geometry.
  2267.  
  2268.     If your viewpoint is at the origin--(0,0,0)--and the window sits
  2269.     parallel to the x-y plane but at z=1, projection is no more than
  2270.     (x,y,z) in the world appears at (x/z,y/z,1) on the plane. Distant
  2271.     objects will have large z values, causing them to shrink in the
  2272.     picture. That's perspective.
  2273.  
  2274.     The trick is to take an arbitrary viewpoint and plane, and
  2275.     transform the world so we have the simple viewing situation.
  2276.     There are two steps: move the viewpoint to the origin, then move
  2277.     the viewplane to the z=1 plane. If the viewpoint is at (vx,vy,vz),
  2278.     transform every point by the translation (x,y,z) -->
  2279.     (x-vx,y-vy,z-vz). This includes the viewpoint and the viewplane.
  2280.     Now we need to rotate so that the z axis points straight at the
  2281.     viewplane, then scale so it is 1 unit away.
  2282.  
  2283.     After all this, we may find ourselves looking out upside- down. It
  2284.     is traditional to specify some direction in the world or viewplane
  2285.     as "up", and rotate so the positive y axis points that way (as
  2286.     nearly as possible if it's a world vector). Finally, we have acted
  2287.     so far as if the window was the entire plane instead of a limited
  2288.     portal. A final shift and scale transforms coordinates in the
  2289.     plane to coordinates on the screen, so that a rectangular region
  2290.     of interest (our "window") in the plane fills a rectangular region
  2291.     of the screen (our "canvas" if you like).
  2292.  
  2293.     Details of how to define and perform the rotation of the viewplane
  2294.     have been left out, but see "How can I aim a camera in a specific
  2295.     direction?" elsewhere in this FAQ. One simple way to designate a
  2296.     plane is with the point closest to the origin, call it D. Then
  2297.     a point P is on the plane if D.P = D.D; or using d = ||D|| and
  2298.     N = D/d, if N.P = d. Aim the camera with N, and scale with d.
  2299.  
  2300.     A further practical difficulty is the need to clip away parts of
  2301.     the world behind us, so -(x,y,z) doesn't pop up at (x/z,y/z,1).
  2302.     (Notice the mathematics of projection alone would allow that!) In
  2303.     fact ordinarily a clipping box, the "viewing frustum", is used
  2304.     to eliminate parts of the scene outside the window left or right,
  2305.     top or bottom, and too close or too far.
  2306.  
  2307.     All the viewing transformations can be done using translation,
  2308.     rotation, scale, and a final perspective divide. If a 4x4
  2309.     homogeneous matrix is used, it can represent everything needed,
  2310.     which saves a lot of work. 
  2311.  
  2312. ----------------------------------------------------------------------
  2313. Subject 5.15: How do I optimize/simplify a 3D polygon mesh?
  2314.  
  2315.     References:
  2316.     "Mesh Optimization" Hoppe, DeRose Duchamp, McDonald, Stuetzle,
  2317.      ACM COMPUTER GRAPHICS Proceedings, Annual Conference Series, 1993.
  2318.  
  2319.     "Re-Tiling Polygonal Surfaces",
  2320.      Greg Turk, ACM Computer Graphics, 26, 2, July 1992
  2321.  
  2322.     "Decimation of Triangle Meshes", Schroeder, Zarge, Lorensen,
  2323.     ACM Computer Graphics, 26, 2 July 1992
  2324.  
  2325.     "Simplification of Objects Rendered by Polygonal Approximations",
  2326.     Michael J. DeHaemer, Jr. and Michael J. Zyda, Computer & Graphics,
  2327.     Vol. 15, No. 2, 1991, Great Britain: Pergamon Press, pp. 175-184.
  2328.  
  2329.     "Topological Refinement Procedures for Triangular Finite Element
  2330.     Procedures", S. A. Cannan, S. N. Muthukrishnan and R. P. Phillips,
  2331.     Engineering With Computers, No. 12, p. 243-255, 1996.
  2332.  
  2333.     "Progressive Meshes", Hoppe, SIGGRAPH 96,
  2334.       http://research.microsoft.com/~hoppe/siggraph96pm/paper.htm
  2335.  
  2336.     Several papers by Michael Garland (quadric-based error metric):
  2337.       http://graphics.cs.uiuc.edu/~garland/
  2338.  
  2339.     Demos:
  2340.       By Stan Melax: http://www.cs.ualberta.ca/~melax/polychop/
  2341.       By Stefan Krause: http://www.stefan-krause.com  [Gnu Open Source]
  2342.       By "klaudius": http://www.klaudius.free.fr
  2343.  
  2344.  
  2345. ----------------------------------------------------------------------
  2346. Subject 5.16: How can I perform volume rendering?
  2347.  
  2348.     Two principal methods can be used:
  2349.     - Ray casting or front-to-back, where the volume is behind the
  2350.       projection plane. A ray is projected from each point in the projection
  2351.       plane through the volume. The ray accumulates the properties of each
  2352.       voxel it passes through.
  2353.     - Object order or back-to-front, where the projection plane is behind
  2354.       the volume. Each slice of the volume is projected on the projection
  2355.       plane, from the farest plane to the nearest plane.
  2356.  
  2357.     You can also use the marching-cubes algorithm, if the extraction of
  2358.     surfaces from the data set is easy to do (see Subject 5.10).
  2359.  
  2360.     Here is one algorithm to do front-to-back volume rendering:
  2361.  
  2362.     Set up a projection plane as a plane tangent to a sphere that encloses 
  2363.     the volume. From each pixel of the projection plane, cast a ray 
  2364.     through the volume by using a Bresenham 3D algorithm.
  2365.     The ray accumulates properties from each voxel intersected, stopping
  2366.     when the ray exits the volume. The pixel value on
  2367.     the projection plane determines the final color of the ray.
  2368.  
  2369.     For unshaded images (i.e., without gradient and light computations),
  2370.     if     C  is the ray color            t  the ray transparency
  2371.         C' the new ray color           t' the new ray transparency
  2372.         Cv the voxel color             tv the voxel transparency
  2373.     then:
  2374.         C' = C + t*Cv       and        t' = t * tv
  2375.     with initial values: C = 0.0 and t = 1.0
  2376.  
  2377.     An alternate version: instead of C' = C + t*Cv , use :
  2378.         C' = C + t*Cv*(1-tv)^p     with p a float variable.
  2379.     Sometimes this gives the best results.
  2380.     When the ray has accumulated transparency, if it becomes negligible
  2381.     (e.g., t<0.1), the process can stop and proceed to the next ray.
  2382.  
  2383.     References:
  2384.  
  2385.     Bresenham 3D:
  2386.       - http://www.sct.gu.edu.au/~anthony/info/graphics/bresenham.procs
  2387.       - [Gems IV] p. 366
  2388.     Volume rendering:
  2389.       - [Watt:Animation] pp. 297-321
  2390.       - IEEE Computer Graphics and application
  2391.         Vol. 10, Nb. 2, March 1990 - pp. 24-53
  2392.       - "Volume Visualization"
  2393.         Arie Kaufman - Ed. IEEE Computer Society Press Tutorial
  2394.       - "Efficient Ray Tracing of Volume Data"
  2395.         Marc Levoy - ACM Transactions on Graphics, Vol. 9, Nb 3, July 1990
  2396.  
  2397. ----------------------------------------------------------------------
  2398. Subject 5.17: Where can I get the spline description of the famous teapot etc.?
  2399.  
  2400.     See the Standard Procedural Databases software, whose official
  2401.     distribution site is 
  2402.         http://www.acm.org/tog/resources/SPD/
  2403.     This database contains much useful 3D code, including spline surface
  2404.     tessellation, for the teapot.
  2405.  
  2406. ----------------------------------------------------------------------
  2407. Subject 5.18: How can the distance between two lines in space be computed?
  2408.  
  2409.     Let x_i be points on the respective lines and n_i unit direction         
  2410.     vectors along the lines.  Then the distance is
  2411.        | (x_1 - x_0)^T (n_1 X n_0) | / || n_1 X n_0 ||.
  2412.     Often one wants the points of closest approach as well as the distance.
  2413.     The following is robust C code from Seth Teller that computes the 
  2414.     these points on two 3D lines.  It also classifies 
  2415.     the lines as parallel, intersecting, or (the generic case) skew.
  2416.     What's listed below shows the main ideas; the full code is at
  2417.       http://graphics.lcs.mit.edu/~seth/geomlib/linelinecp.c
  2418.  
  2419.     // computes pB ON line B closest to line A
  2420.     // computes pA ON line A closest to line B
  2421.     // return: 0 if parallel; 1 if coincident; 2 if generic (i.e., skew)
  2422.     int
  2423.     line_line_closest_points3d (
  2424.         POINT *pA, POINT *pB,            // computed points
  2425.         const POINT *a, const VECTOR *adir,        // line A, point-normal form
  2426.         const POINT *b, const VECTOR *bdir )    // line B, point-normal form
  2427.     {
  2428.     static VECTOR Cdir, *cdir = &Cdir;
  2429.     static PLANE Ac, *ac = &Ac, Bc, *bc = &Bc;
  2430.     
  2431.         // connecting line is perpendicular to both
  2432.         vcross ( cdir, adir, bdir );
  2433.     
  2434.         // check for near-parallel lines
  2435.         if ( !vnorm( cdir ) )   {
  2436.             *pA = *a;    // all points are closest
  2437.             *pB = *b;
  2438.             return 0;    // degenerate: lines parallel
  2439.             }
  2440.     
  2441.         // form plane containing line A, parallel to cdir
  2442.         plane_from_two_vectors_and_point ( ac, cdir, adir, a );
  2443.     
  2444.         // form plane containing line B, parallel to cdir
  2445.         plane_from_two_vectors_and_point ( bc, cdir, bdir, b );
  2446.     
  2447.         // closest point on A is line A ^ bc
  2448.         intersect_line_plane ( pA, a, adir, bc );
  2449.     
  2450.         // closest point on B is line B ^ ac
  2451.         intersect_line_plane ( pB, b, bdir, ac );
  2452.     
  2453.         // distinguish intersecting, skew lines
  2454.         if ( edist( pA, pB ) < 1.0E-5F )
  2455.             return 1;     // coincident: lines intersect
  2456.         else
  2457.             return 2;     // distinct: lines skew
  2458.     }
  2459.  
  2460.     Also Dave Eberly has code for computing distance between various 
  2461.     geometric primitives, including MinLineLine(), at
  2462.       http://www.magic-software.com
  2463.  
  2464.  
  2465. ----------------------------------------------------------------------
  2466. Subject 5.19: How can I compute the volume of a polyhedron?
  2467.  
  2468.     Assume that the surface is closed, every face is a triangle, and
  2469.     the vertices of each triangle are oriented ccw from the outside.
  2470.     Let Volume(t,p) be the signed volume of the tetrahedron formed
  2471.     by a point p and a triangle t.  This may be computed by a 4x4
  2472.     determinant, as in [O'Rourke (C), p.26].
  2473.             Choose an arbitrary point p (e.g., the origin), and compute
  2474.     the sum Volume(t_i,p) for every triangle t_i of the surface.  That
  2475.     is the volume of the object.  The justification for this claim
  2476.     is nontrivial, but is essentially the same as the justification for
  2477.     the computation of the area of a polygon (Subject 2.01).
  2478.  
  2479.     C Code available at http://cs.smith.edu/~orourke/ 
  2480.     and http://www.acm.org/jgt/papers/Mirtich96/ .
  2481.  
  2482.     For computing the volumes of n-d convex polytopes, 
  2483.     there is a C implementation by Bueeler and Enge of various
  2484.     algorithms available at
  2485.  
  2486.     http://www.Mathpool.Uni-Augsburg.DE/~enge/Volumen.html .
  2487.  
  2488. ----------------------------------------------------------------------
  2489. Subject 5.20:  How can I decompose a polyhedron into convex pieces?
  2490.  
  2491.     Usually this problem is interpreted as seeking a collection
  2492.     of pairwise disjoint convex polyhedra whose set union is the
  2493.     original polyhedron P.  The following facts are known about
  2494.     this difficult problem:
  2495.  
  2496.        o  Not every polyhedron may be partitioned by diagonals into
  2497.           tetrahedra.  A counterexample is due to Scho:nhardt
  2498.           [O'Rourke (A), p.254].
  2499.  
  2500.        o  Determining whether a polyhedron may be so partitioned is
  2501.           NP-hard, a result due to Seidel & Ruppert [Disc. Comput. Geom. 
  2502.           7(3) 227-254 (1992).]
  2503.  
  2504.        o  Removing the restriction to diagonals, i.e., permitting
  2505.           so-called Steiner points, there are polyhedra of n vertices
  2506.           that require n^2 convex pieces in any decomposition.
  2507.           This was established by Chazelle [SIAM J. Comput. 
  2508.           13: 488-507 (1984)]. See also [O'Rourke (A), p.256]
  2509.  
  2510.        o  An algorithm of Palios & Chazelle guarantees at most
  2511.           O(n+r^2) pieces, where r is the number of reflex edges
  2512.           (i.e., grooves). [Disc. Comput. Geom. 5:505-526 (1990).]
  2513.  
  2514.        o  Barry Joe's geompack package, at 
  2515.              ftp://ftp.cs.ualberta.ca/pub/geompack,
  2516.           includes a 3D convex decomposition Fortran program. 
  2517.  
  2518.        o  There seems to be no other publicly available code.
  2519.  
  2520. ----------------------------------------------------------------------
  2521. Subject 5.21: How can the circumsphere of a tetrahedron be computed?
  2522.  
  2523.  
  2524.     Let a, b, c, and d be the corners of the tetrahedron, with
  2525.     ax, ay, and az the components of a, and likewise for b, c, and d.  
  2526.     Let |a| denote the Euclidean norm of a, and let a x b denote the 
  2527.     cross product of a and b.  Then the center m and radius r of the
  2528.     circumsphere are given by
  2529.     
  2530.        |                                                                       |
  2531.        | |d-a|^2 [(b-a)x(c-a)] + |c-a|^2 [(d-a)x(b-a)] + |b-a|^2 [(c-a)x(d-a)] |
  2532.        |                                                                       |
  2533.     r= -------------------------------------------------------------------------
  2534.                                  | bx-ax  by-ay  bz-az |
  2535.                                2 | cx-ax  cy-ay  cz-az |
  2536.                                  | dx-ax  dy-ay  dz-az |
  2537.    
  2538.     and
  2539.    
  2540.            |d-a|^2 [(b-a)x(c-a)] + |c-a|^2 [(d-a)x(b-a)] + |b-a|^2 [(c-a)x(d-a)]
  2541.     m= a + ---------------------------------------------------------------------
  2542.                                    | bx-ax  by-ay  bz-az |
  2543.                                  2 | cx-ax  cy-ay  cz-az |
  2544.                                    | dx-ax  dy-ay  dz-az |
  2545.     
  2546.     Some notes on stability:
  2547.     
  2548.     - Note that the expression for r is purely a function of differences between
  2549.       coordinates.  The advantage is that the relative error incurred in the
  2550.       computation of r is also a function of the _differences_ between the
  2551.       vertices, and is not influenced by the _absolute_ coordinates of the
  2552.       vertices.  In most applications, vertices are usually nearer to each other
  2553.       than to the origin, so this property is advantageous.
  2554.     
  2555.       Similarly, the formula for m incurs roundoff error proportional to the
  2556.       differences between vertices, but not proportional to the absolute 
  2557.       coordinates of the vertices until the final addition.
  2558.     
  2559.     - These expressions are unstable in only one case:  if the denominator is
  2560.       close to zero.  This instability, which arises if the tetrahedron is 
  2561.       nearly degenerate, is unavoidable.  Depending on your application, you 
  2562.       may want to use exact arithmetic to compute the value of the determinant.
  2563.       See
  2564.         http://www.geom.umn.edu/software/cglist/alg.html
  2565.       or
  2566.         http://www.cs.cmu.edu/~quake/robust.html
  2567.  
  2568. ----------------------------------------------------------------------
  2569. Subject 5.22: How do I determine if two triangles in 3D intersect?
  2570.  
  2571.     Let the two triangles be T1 and T2.  If T1 lies strictly to one side 
  2572.     of the plane containing T2, or T2 lies strictly to one side of the 
  2573.     plane containing T1, the triangles do not intersect.  Otherwise,
  2574.     compute the line of intersection L between the planes.  Let Ik
  2575.     be the interval (Tk inter L), k=1,2.  Either interval may be empty.
  2576.     T1 and T2 intersect iff I1 and I2 overlap.
  2577.  
  2578.     This method is decribed in Tomas Mo:ller, "A fast triangle-triangle
  2579.     intersection test," J. Graphics Tools 2(2) 25-30 1997.  C code at
  2580.     http://www.acm.org/jgt/papers/Moller97/ .  See also
  2581.     http://www.ce.chalmers.se/staff/tomasm/code/ 
  2582.     http://www.magic-software.com/MgcIntersection.html
  2583.  
  2584.     See also the "3D Object Intersection" page, described in Subject 0.05.
  2585.  
  2586.     NASA's "Intersect" code will intersect any number of triangulated
  2587.     surfaces provided that each of the surfaces is both closed and manifold.
  2588.       http://www.nas.nasa.gov/~aftosmis/cart3d/surfaceModeling.html#AuxProgs
  2589.     Based on "Robust and Efficient Cartesian Mesh Generation for 
  2590.     Component-Based Geometry" AIAA Paper 97-0196. Michael Aftosmis.
  2591.  
  2592. ----------------------------------------------------------------------
  2593. Subject 5.23:  How can a 3D surface be reconstructed from a collection of points?
  2594.  
  2595.     This is a difficult problem.  There are two main variants: 
  2596.     (1) when the points are organized into parallel slices through 
  2597.     the object; (2) when the points are unorganized.  
  2598.     For (1), see this survey:
  2599.        D. Meyers, S. Skinner, K. Sloan. "Surfaces from Contours"
  2600.        ACM Trans. Graph. 11(3) Jul 1992, 228--258.
  2601.        http://www.acm.org/pubs/citations/journals/tog/1992-11-3/p228-meyers/
  2602.     Code (NUAGES) is available at 
  2603.        http://www-sop.inria.fr/prisme/logiciel/nuages.html.en
  2604.        ftp://ftp-sop.inria.fr/prisme/NUAGES/Nuages/NUAGES_SRC.tar.gz
  2605.     For (2), see this paper:
  2606.        H. Hoppe, T. DeRose, T. Duchamp, J. McDonald, W. Stuetzle
  2607.        "Surface reconstruction from unorganized points"
  2608.        Proc. SIGGRAPH '92, 71--78.
  2609.     and P. Kumar's collection of links at
  2610.        http://members.tripod.com/~GeomWiz/www.sites.html
  2611.     New code, Cocone, written with CGAL, based on recent work by 
  2612.     N. Amenta, S. Choi, T. K. Dey and N. Leekha:
  2613.        http://www.cis.ohio-state.edu/~tamaldey/cocone.html
  2614.  
  2615. ----------------------------------------------------------------------
  2616. Subject 5.24: How can I find the smallest sphere enclosing a set of points in 3D? 
  2617.  
  2618.     Although not obvious, the smallest sphere enclosing a set of points
  2619.     in any dimension can be found by Linear Programming.  This was proved
  2620.     by Emo Welzl in the paper, "Smallest enclosing disks (balls and 
  2621.     ellipsoids)" [Lecture Notes Comput. Sci., 555, Springer-Verlag, 1991, 
  2622.     359--370].  +   Code developed by Bernd Gaertner available (GNU 
  2623.     General Public License) at:
  2624.         http://www.inf.ethz.ch/~gaertner/miniball.html
  2625.     This code is remarkably efficient: e.g., 2 seconds for 10^6 points in 
  2626.     3D on an Ultra-Sparc II.  See also Dave Eberly's direct implementation
  2627.     of Welzl's algorithm:
  2628.         http://www.magic-software.com/MgcContainment3D.html
  2629.  
  2630. ----------------------------------------------------------------------
  2631. Subject 5.25: What's the big deal with quaternions?
  2632.  
  2633.     This could mean "Why do they evoke such heated debate?" or "What are
  2634.     their virtues?"
  2635.  
  2636.     The heat of debate is hard to explain, but it's been happening for many
  2637.     decades. When Gibbs first deprecated the quaternion product and split it
  2638.     into a cross product and a dot product, one outraged Victorian called
  2639.     the result a "hermaphrodite monster" -- and that before the Internet's
  2640.     flame wars. Generally, the quaternion advocates seem to feel the
  2641.     opponents are lazy or thick-headed, and that deeper understanding of
  2642.     quaternions would lead to deeper appreciation. The opponents don't
  2643.     appreciate that attitude, and seem to feel the advocates are snooty or
  2644.     sheep, and that matrices and such are less abstract and do just fine.
  2645.     (Advocates of Clifford algebra would claim that both sides are mired in
  2646.     the past.) Passion aside, quaternions have appropriate uses, as do their
  2647.     alternatives.
  2648.  
  2649.     Someone new to the debate first needs to know what quaternions are, and
  2650.     what they're supposed to be good for. Quaternions are a quadruple of
  2651.     numbers, used to represent 3D rotations.
  2652.         q = [x,y,z,w] = [(x,y,z),w] = [V,w]
  2653.  
  2654.     The "norm" of a quaternion N(q) is conventionally the sum of the squares
  2655.     of the four components. Some writers confuse this with the magnitude,
  2656.     which is the square root of the norm. Another common misconception is
  2657.     that only quaternions of unit norm can be used, those with the sum of
  2658.     the four squares equal to 1, but that is wrong (though they are
  2659.     preferred).
  2660.         [U sin a,cos a] rotates by angle 2a around unit vector U
  2661.  
  2662.     Popular non-quaternion options are 3x3 special orthogonal matrices (9
  2663.     numbers with constraints), Euler angles (3 numbers), axis-angle (4
  2664.     numbers), and angular velocity vectors (3 numbers). None of these
  2665.     options actually _are_ rotations, which are physical; they _represent_
  2666.     rotations. The distinction is important because, for example, it is
  2667.     common to use an axis-angle with an angle greater than 360 degrees to
  2668.     tell an animation system to spin an object more than a full turn,
  2669.     something a matrix cannot say. In mathematics, the usual meaning of a
  2670.     rotation would not allow the multiple spin version, which can lead to
  2671.     confusing debates.
  2672.         q and -q represent the same rotation
  2673.  
  2674.     Two rotations, the physical things, can be applied one after the other.
  2675.     Assuming the two rotation axes have a least one point in common, the
  2676.     result will be another rotation. Some rotation representations handle
  2677.     this gracefully, some don't. For quaternions and matrices, forms of
  2678.     multiplication are defined such that the product gives the desired
  2679.     result. For Euler angles especially there is no simple computation.
  2680.         q = q2 q1 = [V2,w2][V1,w1] = [V2xV1+w2V1+w1V2,w2w1-V2.V1]
  2681.  
  2682.     Every rotation has a reverse rotation, such that the combination of the
  2683.     two leaves an object as it was. (Rotations are an algebraic "group".)
  2684.     Euler angles make reversals difficult to compute. Other representations,
  2685.     including quaternions, make them simple.
  2686.         reverse([V,w]) = [V,-w]
  2687.         q^(-1) = [-V,w]/(V.V+ww)
  2688.  
  2689.     Two physical rotations are also more or less similar. Unit quaternions
  2690.     do a particularly good job of representing similar rotations with
  2691.     similar numbers.
  2692.         similar(q1,q2) = |q1.q2| = | x1x2+y1y2+z1z2+w1w2 |
  2693.  
  2694.     Points in space, the physical things, are normally represented as 3 or 4
  2695.     numbers. The effect of a rotation on a collection of points can be
  2696.     computed from the representation of the rotation, and here matrices seem
  2697.     fastest, using three dot products. Using their own product twice,
  2698.     quaternions are a bit less efficient. (They are usually converted to
  2699.     matrices at the last minute.)
  2700.         p2 = q p1 q^(-1)
  2701.  
  2702.     Sequences of rotations can be interpolated, so that the object being
  2703.     turned is rotated to specific poses at specific times. This motivated
  2704.     Ken Shoemake's early use of quaternions in computer graphics, as
  2705.     published in 1985. He used an analog of linear interpolation (sometimes
  2706.     called "lerp") that he called "Slerp", and also introduced an analog of
  2707.     a piecewise Bezier curve. A few years later in some course notes he
  2708.     described another curve variation he called "Squad", which still seems
  2709.     to be popular. Later authors have proposed many alternatives.
  2710.                             sin (1-t)A      sin tA
  2711.         Slerp(q1,q2;t) = q1 ---------- + q2 ------, cos A = q1.q2
  2712.                                sin A         sin A
  2713.  
  2714.         Squad(q1,a1,b2,q2;t) = Slerp(Slerp(q1,q2;t),
  2715.                                      Slerp(a1,b2;t);
  2716.                                      2t(1-t))
  2717.  
  2718.     Physics simulation, aerospace control, and robotics are examples of
  2719.     computations which also depend on rotation representation. Constrained
  2720.     rotations like a wheel on an axle or the elbow bend of a robot typically
  2721.     use specialized representations, such as an angle alone. In many general
  2722.     situations, however, quaternions have proved valuable.
  2723.         2 dq = W q dt, W is the angular velocity vector
  2724.  
  2725.     User interfaces for 3D rotation also require a representation. Direct
  2726.     manipulation interfaces typically use angles for jointed figures, but
  2727.     for freer manipulation may use quaternions, as in Arcball or
  2728.     through-the-lens camera control. As Shoemake's _full_ Graphics Gems code
  2729.     for Arcball demonstrates (with the [CAPS LOCK] key), any rotation can be
  2730.     graphed as an arc on a sphere. (Not to be confused with the quaternion
  2731.     unit sphere in 4D.) Whether quaternions, or any other representation,
  2732.     are helpful for numeric presentation and input seems a matter of taste
  2733.     and circumstance.
  2734.         q = U2 U1^(-1) = [U1xU2,U1.U2]
  2735.  
  2736.     Because of their metric properties for representing rotations, unit
  2737.     quaternions are most common. Advocates frequently point out that it is
  2738.     far cheaper to normalize the length of a non-zero quaternion than to
  2739.     bring a matrix back to rotation form. Also Shoemake's later conversion
  2740.     code cheaply creates a correct rotation matrix from _any_ quaternion
  2741.     (found with his Euler angle code from Graphics Gems, which does the same
  2742.     for all 24 variations of that representation).
  2743.         Normalize(q) = q/Sqrt(q.q)
  2744.  
  2745.     Comparisons to Euler angles may mention "gimbal lock" (frequently
  2746.     misspelled) as a disadvantage quaternions avoid. In the physical world
  2747.     where gyroscopes are mounted on nested pivots, which are called gimbals,
  2748.     locking is a real problem quaternions cannot help. What's usually meant
  2749.     is that because the similarity of rotations organizes them somewhat like
  2750.     a sphere, while similarity of vectors is quite different, an inevitable
  2751.     misfit plagues Euler angles. This can show up in behavior much like
  2752.     physical gimbal lock, but also in other ways. The difficulties are
  2753.     topological, and aiming runs into them as well, even if quaternions are
  2754.     used. Quaternion authors who propose using curves in the vector space of
  2755.     quaternion logarithms often risk the misfit unawares. Frankly, you must
  2756.     pick through the literature carefully, whether informal and online or
  2757.     refereed and printed, because mistakes are tragically common.
  2758.  
  2759.     To explore Graphics Gems code, see "Where is all the source?" in this
  2760.     FAQ. To read more about quaternions, you have many options. Since they
  2761.     were discovered in 1843 by Hamilton (the same Irish mathematician and
  2762.     physicist whose name shows up in quantum mechanics), quaternions have
  2763.     found many friends, as a web search will reveal. Quaternions can be
  2764.     approached and applied in numerous different ways, so if you keep
  2765.     looking it's likely you will find something that suits your taste and
  2766.     your needs.
  2767.  
  2768.       (Subject 0.04) [Eberly], Ch. 2.
  2769.       http://www.maths.tcd.ie/pub/HistMath/People/Hamilton/OnQuat/
  2770.         Hamilton's original paper. Not easy for today's readers.
  2771.       K. Shoemake. Animating Rotation with Quaternion Curves.
  2772.       Proceedings of Siggraph 85.
  2773.         Original animation method. Covers all the basics.
  2774.       ftp://ftp.cis.upenn.edu/pub/graphics/shoemake/quatut.ps
  2775.         Later Shoemake tutorial. More complete than most authors.
  2776.       Graphics Gems I-V, various authors, various articles.
  2777.         As usual, a helpful source of code and discussion.
  2778.       ftp://ftp.netcom.com/pub/hb/hbaker/quaternion/
  2779.         Henry Baker collects good quaternion stuff. Access spotty.
  2780.       http://linux.rice.edu/~rahul/hbaker/quaternion/
  2781.         Henry Baker collection with more reliable access.
  2782.       http://www.eecs.wsu.edu/~hart/papers/vqr.pdf
  2783.         Visualizing quaternion rotation. May help build intuition.
  2784.       http://graphics.stanford.edu/courses/cs348c-95-fall
  2785.       /software/quatdemo/
  2786.         The GL code implementing above Hart et al. paper.
  2787.       http://www.diku.dk/students/myth/quat.html
  2788.         Mathematical, but not too fast and not too fancy.
  2789.       http://www.cs.berkeley.edu/~laura/cs184/quat/quaternion.html
  2790.         Laura Downs covers the fundamentals.
  2791.       http://graphics.cs.ucdavis.edu/GraphicsNotes/Quaternions
  2792.       /Quaternions.htm
  2793.         Ken Joy covers the fundamentals.
  2794.       http://www.gg.caltech.edu/STC/rr_sig97.html
  2795.         High-tech interpolation method. Demanding but rewarding.
  2796.       Duff, Tom: Quaternion Splines for Animating Orientation,
  2797.       Proceedings of the USENIX Association Second Computer Graphics
  2798.       Workshop (held in Monterey, CA 12-13 Dec. 1985), pp. 54-62.
  2799.         Subdivision paper in odd place. Author last seen at Pixar.
  2800.  
  2801. ----------------------------------------------------------------------
  2802. Subject 5.26:  How can I aim a camera in a specific direction?
  2803.  
  2804.     What's needed is a method for creating a rotation that turns one unit
  2805.     vector to line up with another. To aim at an object, you can subtract
  2806.     the position of the camera from the position of the object to get a
  2807.     vector which you then normalize. The vector you want to turn is the
  2808.     camera forward vector, commonly a unit vector along the camera -z axis.
  2809.     Be warned that more than one rotation can achieve aim alone. (The issue
  2810.     is rather complicated, as laid out in Ken Shoemake's article on twist
  2811.     reduction in Graphics Gems IV.) For example, even if the camera is
  2812.     already properly aimed you could rotate it around its z axis. The most
  2813.     direct rotation is given by the non-unit quaternion
  2814.         q = [(b,-a,0),1-c], to aim -z axis along unit vector (a,b,c)
  2815.     Normalization is advised, but it will fail for aim vector (0,0,1). In
  2816.     that case just rotate 180 degrees around the y axis, using
  2817.         q = [(0,1,0),0]
  2818.  
  2819.     If the camera is level after rotation by quaternion [(x,y,z),w], the y
  2820.     component of its rotated x axis will be zero, so
  2821.         xy+wz = 0
  2822.     If it is upright, the y component of its rotated y axis will not be
  2823.     negative, so
  2824.         ww-xx+yy-zz >= 0
  2825.  
  2826.     To ensure these two desirable properties, aim with a more sophisticated
  2827.     non-unit quaternion
  2828.         [(bs,-at,ab),st], where s = r-c, t = r+1, and r = sqrt(aa+cc).
  2829.     This can also fail to normalize, in which case normalize instead
  2830.         [(0,1+c,-b),0]
  2831.     Unless the aim vector is null, this will succeed. If the aim vector has
  2832.     not been normalized and its magnitude is m = sqrt(aa+bb+cc), substitute
  2833.     1->m. That is, use t = r+m and use m+c.
  2834.  
  2835.     More generally, to rotate unit vector U1 directly to unit vector U2, the
  2836.     non-unit quaternion will be
  2837.         q = [U1xU2,1+U1.U2]
  2838.  
  2839.     Why? If U is a unit vector, then it is normal to a plane through the
  2840.     origin with equation U.P = 0. Reflection in that plane is given by
  2841.     reversing the U component of P.
  2842.         reflect(P,U) = P ^╓ 2(U.P)U
  2843.     The quaternion product of U1 and U2 is [U1xU2,-U1.U2], so
  2844.         -2 (U.P) = PU + UP
  2845.     Noting UU = -1, this gives a quaternion reflection formula.
  2846.         reflect(P,U) = P + (PU+UP)U
  2847.                      = P ^╓ P + UPU
  2848.                      = UPU
  2849.     Reflecting with U1 then U2, by U2(U1 P U1)U2, rotates by twice the angle
  2850.     between the planes, with axis perpendicular to both normals. Noting U1U2
  2851.     is the conjugate of U2U1, and -q rotates like q, the rotation quaternion
  2852.     is
  2853.         q = -U2U1 = -[U2xU1,-U2.U1] = [U1xU2,U1.U2]
  2854.     This q fails to aim U1 at U2 by rotating twice as much as needed, but
  2855.     its square root succeeds. One square root of unit q is 1+q normalized,
  2856.     geometrically the bisection of the great arc from the identity to q.
  2857.     There is an inevitable singularity when U2 is the opposite of U1,
  2858.     because any perpendicular axis gives an equally direct 180 degree
  2859.     rotation.
  2860.  
  2861.     [These quaternion methods were provided by Ken Shoemake.]
  2862.  
  2863. ----------------------------------------------------------------------
  2864. Subject 5.27: How do I transform normals?
  2865.  
  2866.     In 3D, the orientation of a plane in space can be given by a
  2867.     vector perpendicular to the plane, a "normal vector" or "normal"
  2868.     for short. Often it is convenient to keep that vector of unit
  2869.     length, or "normalized"; be careful of the different meanings
  2870.     of "normality". A smooth surface has a plane tangent to each
  2871.     point, and by extension a normal to that plane is called a
  2872.     "surface normal". Graphics code also cheats by associating
  2873.     artificial normal vectors with the vertices of polygonal models
  2874.     to simulate the reflection properties of curved surfaces;
  2875.     these are called "vertex normals".
  2876.     
  2877.     The "orientation" of a plane has two meanings, both of which
  2878.     usually apply. Aside from the rotational tilting and turning
  2879.     meaning, there is also the sense of "side". A closed convex
  2880.     surface made of polygons has two sides, an inside and an
  2881.     outside, and normals can be assigned to the polygons in such
  2882.     a way that they all consistently point outside. This is
  2883.     often desirable for shading and culling.
  2884.     
  2885.     When a model is defined in one coordinate system and used in
  2886.     another, as is commonly done, it may be necessary to transform
  2887.     normals also. If the change of point coordinates is effected
  2888.     by means of a rotation plus a translation, one simply rotates
  2889.     the normals as well (with no translation). Other coordinate
  2890.     changes need more care.
  2891.     
  2892.     Although it is possible to use projective transformations to
  2893.     map a model into world coordinates, ordinarily they are used
  2894.     only for viewing. It is usually a mistake to apply perspective
  2895.     to a normal, as shading and culling are best done in world
  2896.     coordinates for correct results. Also perspective may be
  2897.     computed using degenerate matrices which are not invertible,
  2898.     though that need not be the case. For the importance of this,
  2899.     see below.
  2900.     
  2901.     The combination of a linear transformation and a translation
  2902.     is called an affine transformation, and is performed as a
  2903.     matrix multiplication plus a vector addition:
  2904.       p' = A(p) = Lp + t.
  2905.     When the model-to-world point transform is affine, the proper
  2906.     way to transform normals is with the transpose of the inverse
  2907.     of L.
  2908.       n' = (L^{-1})^T n
  2909.     However that is not enough.
  2910.     
  2911.     If L includes scaling effects, a unit normal in model space
  2912.     will usually transform to a non-unit normal, which can cause
  2913.     problems for shaders and other code. This may need correcting
  2914.     after the normal is transformed.
  2915.     
  2916.     If L includes reflection, the inside-outside orientation of
  2917.     the normal is reversed. This, too, can cause problems, and
  2918.     may need correcting. The determinant of L will be negative
  2919.     in this case.
  2920.     
  2921.     When a complicated distortion is used, it must be approximated
  2922.     differently at each point in space by a linear transform made
  2923.     up of partial derivatives, the Jacobian. The matrix for the
  2924.     Jacobian replaces L in the equation for transforming normals.
  2925.     
  2926.     Why use the transposed inverse?
  2927.     Write the dot product of column vectors n and v as a matrix
  2928.     product n^T v. Write vector v as a difference of points, q-p.
  2929.     Let p, q, and thus v lie in the desired plane, and let n be
  2930.     normal to it. Vectors at right angles have zero dot product.
  2931.       n^T v = 0
  2932.     The transform v' of v is
  2933.       v' = (Lq+t)-(Lp+t)
  2934.          = (Lq-Lp)+(t-t)
  2935.          = L(q-p)
  2936.          = Lv
  2937.     The transform n' of n will remain normal if it satisfies
  2938.       n'^T v' = n^T v
  2939.     Let n' = Mn for some M. Then the requirement is
  2940.       n'^T v' = (Mn)^T (Lv)
  2941.               = n^T (M^T L) v
  2942.               = n^T v
  2943.     This holds if
  2944.       M^T L = I
  2945.     where I is the identity. Right multiplying by the inverse
  2946.     L^{-1} and transposing both sides gives, as claimed,
  2947.       M = (L^{-1})^T
  2948.     When L is a rotation, L^{-1} = L^T, so M = L. When L has no
  2949.     inverse it will still have an "adjoint" to substitute for
  2950.     for orthogonality purposes, differing only by a scale factor.
  2951.  
  2952. ----------------------------------------------------------------------
  2953. Section 6. Geometric Structures and Mathematics
  2954. ----------------------------------------------------------------------
  2955. Subject 6.01: Where can I get source for Voronoi/Delaunay triangulation?
  2956.  
  2957.     For 2-d Delaunay triangulation, try Shewchuk's triangle program.  It
  2958.     includes options for constrained triangulation and quality mesh
  2959.     generation.  It uses exact arithmetic.
  2960.  
  2961.     The Delaunay triangulation is equivalent to computing the convex hull
  2962.     of the points lifted to a paraboloid.  For n-d Delaunay triangulation
  2963.     try Clarkson's hull program (exact arithmetic) or Barber and Huhdanpaa's
  2964.     Qhull program (floating point arithmetic).  The hull program also
  2965.     computes Voronoi volumes and alpha shapes.  The Qhull program also
  2966.     computes 2-d Voronoi diagrams and n-d Voronoi vertices.  The output of
  2967.     both programs may be visualized with Geomview.
  2968.  
  2969.     There are many other codes for Delaunay triangulation and Voronoi
  2970.     diagrams.  See Amenta's list of computational geometry software.
  2971.     Especially noteworthy is the CGAL code: Subject 0.07 and
  2972.     http://www.cgal.org.
  2973.  
  2974.     The Delaunay triangulation satisfies the following property: the
  2975.     circumcircle of each triangle is empty.  The Voronoi diagram is the
  2976.     closest-point map, i.e., each Voronoi cell identifies the points that
  2977.     are closest to an input site.  The Voronoi diagram is the dual of
  2978.     the Delaunay triangulation.  Both structures are defined for general
  2979.     dimension.  Delaunay triangulation is an important part of mesh
  2980.     generation.
  2981.  
  2982.     There is a FAQ of polyhedral computation explaining how to compute
  2983.     n-d Delaunay triangulation and n-d Voronoi diagram using a convex hull
  2984.     code, and how to use the linear programming technique to 
  2985.     determine the Voronoi cells adjacent to a given Voronoi cell
  2986.     efficiently for large scale or higher dimensional cases.
  2987.  
  2988.     Avis' lrs code uses the same file formats as cdd. It
  2989.     uses exact arithmetic and is useful
  2990.     for problems with very large output size, since it does not
  2991.     require storing the output.
  2992.  
  2993.     On a closely related topic, see
  2994.       http://www.cis.ohio-state.edu/~tamaldey/medialaxis.htm
  2995.     for computation of the 3D medial axis from the Voronoi diagram.
  2996.  
  2997.  
  2998.     References:
  2999.  
  3000.     Amenta:   http://www.geom.umn.edu/software/cglist
  3001.  
  3002.     Avis:     ftp://mutt.cs.mcgill.ca/pub/C/lrs.html
  3003.  
  3004.     Barber &  http://www.geom.umn.edu/locate/qhull
  3005.     Huhdanpaa ftp://geom.umn.edu/pub/software/qhull.tar.Z
  3006.               ftp://ftp.geom.umn.edu/pub/software/qhull.zip
  3007.  
  3008.     Clarkson: http://cm.bell-labs.com/netlib/voronoi/hull.html
  3009.               ftp://cm.bell-labs.com/netlib/voronoi/hull.shar.Z
  3010.  
  3011.     Geomview: http://www.geom.umn.edu/locate/geomview
  3012.               ftp://geom.umn.edu/pub/software/geomview/
  3013.  
  3014.     Polyhedral Computation FAQ:
  3015.               http://www.ifor.math.ethz.ch/ifor/staff/fukuda/fukuda.html
  3016.  
  3017.     Shewchuk: http://www.cs.cmu.edu/~quake/triangle.html
  3018.               ftp://cm.bell-labs.com/netlib/voronoi/triangle.shar.Z
  3019.  
  3020.  
  3021.     [O' Rourke (C)] pp. 168-204
  3022.  
  3023.     [Preparata & Shamos] pp. 204-225
  3024.  
  3025.     Chew, L. P. 1987. "Constrained Delaunay Triangulations," Proc. Third
  3026.     Annual ACM Symposium on Computational Geometry.
  3027.  
  3028.     Chew, L. P. 1989. "Constrained Delaunay Triangulations," Algorithmica
  3029.     4:97-108. (UPDATED VERSION OF CHEW 1987.)
  3030.  
  3031.     Fang, T-P. and L. A. Piegl. 1994. "Algorithm for Constrained Delaunay
  3032.     Triangulation," The Visual Computer 10:255-265. (RECOMMENDED!)
  3033.  
  3034.     Frey, W. H. 1987. "Selective Refinement: A New Strategy for Automatic
  3035.     Node Placement in Graded Triangular Meshes," International Journal for
  3036.     Numerical Methods in Engineering 24:2183-2200.
  3037.  
  3038.     Guibas, L. and J. Stolfi. 1985. "Primitives for the Manipulation of
  3039.     General Subdivisions and the Computation of Voronoi Diagrams," ACM
  3040.     Transactions on Graphics 4(2):74-123.
  3041.  
  3042.     Karasick, M., D. Lieber, and L. R. Nackman. 1991. "Efficient Delaunay
  3043.     Triangulation Using Rational Arithmetic," ACM Transactions on Graphics
  3044.     10(1):71-91.
  3045.  
  3046.     Lischinski, D. 1994. "Incremental Delaunay Triangulation," Graphics
  3047.     Gems IV, P. S. Heckbert, Ed. Cambridge, MA: Academic Press Professional.
  3048.     (INCLUDES C++ SOURCE CODE -- THANK YOU, DANI!)
  3049.  
  3050.     [Okabe]
  3051.  
  3052.     Schuierer, S. 1989. "Delaunay Triangulation and the Radiosity
  3053.     Approach," Proc. Eurographics '89, W. Hansmann, F. R. A. Hopgood, and
  3054.     W. Strasser, Eds. Elsevier Science Publishers, 345-353.
  3055.  
  3056.     Subramanian, G., V. V. S. Raveendra, and M. G. Kamath. 1994. "Robust
  3057.     Boundary Triangulation and Delaunay Triangulation of Arbitrary
  3058.     Triangular Domains," International Journal for Numerical Methods in
  3059.     Engineering 37(10):1779-1789.
  3060.  
  3061.     Watson, D. F. and G. M. Philip. 1984. "Survey: Systematic
  3062.     Triangulations," Computer Vision, Graphics, and Image Processing
  3063.     26:217-223.
  3064.  
  3065. ----------------------------------------------------------------------
  3066. Subject 6.02: Where do I get source for convex hull?
  3067.  
  3068.     For n-d convex hulls, try Clarkson's hull program (exact arithmetic)
  3069.     or Barber and Huhdanpaa's Qhull program (floating point arithmetic).
  3070.     Qhull 3.1 includes triangulated output and improved
  3071.     handling of difficult inputs.   Qhull computes convex hulls,
  3072.     Delaunay triangulations, Voronoi diagrams, and halfspace
  3073.     intersections about a point.
  3074.  
  3075.     Qhull handles numeric precision problems by merging facets. The output
  3076.     of both programs may be visualized with Geomview.
  3077.  
  3078.     In higher dimensions, the number of facets may be much smaller than
  3079.     the number of lower-dimensional features.  When this is the case,
  3080.     Fukuda's cdd program is much faster than Qhull or hull.
  3081.  
  3082.     There are many other codes for convex hulls.  See Amenta's list of
  3083.     computational geometry software.
  3084.  
  3085.     References:
  3086.  
  3087.     Amenta:   http://www.geom.umn.edu/software/cglist/
  3088.  
  3089.     Barber &  http://www.geom.umn.edu/locate/qhull
  3090.     Huhdanpaa 
  3091.  
  3092.     Clarkson: http://cm.bell-labs.com/netlib/voronoi/hull.html
  3093.               ftp://cm.bell-labs.com/netlib/voronoi/hull.shar.Z
  3094.  
  3095.     Geomview: http://www.geom.umn.edu/locate/geomview
  3096.               ftp://geom.umn.edu/pub/software/geomview/
  3097.  
  3098.     Fukuda:   http://www.ifor.math.ethz.ch/staff/fukuda/cdd_home/cdd.html
  3099.               ftp://ftp.ifor.math.ethz.ch/pub/fukuda/cdd/
  3100.  
  3101.  
  3102.     [O' Rourke (C)] pp. 70-167
  3103.     C code for Graham's algorithm on p.80-96.
  3104.     Three-dimensional convex hull discussed in Chapter 4 (p.113-67).
  3105.     C code for the incremental algorithm on p.130-60.
  3106.  
  3107.     [Preparata & Shamos] pp. 95-184
  3108.  
  3109.  
  3110. ----------------------------------------------------------------------
  3111. Subject 6.03: Where do I get source for halfspace intersection?
  3112.  
  3113.     For n-d halfspace intersection, try Fukuda's cdd program or Barber
  3114.     and Huhdanpaa's Qhull program.  Both use floating point arithmetic.
  3115.     Fukuda includes code for exact arithmetic.  Qhull handles numeric
  3116.     precision problems by merging facets.
  3117.  
  3118.     Qhull computes halfspace intersection by computing a convex hull.
  3119.     The intersection of halfspaces about the origin is equivalent to the
  3120.     convex hull of the halfspace coefficients divided by their offsets.
  3121.     See Subject 6.02 for more information.
  3122.  
  3123.     References:
  3124.  
  3125.     Barber &  http://www.geom.umn.edu/locate/qhull
  3126.     Huhdanpaa 
  3127.  
  3128.     Fukuda:   ftp://ifor13.ethz.ch/pub/fukuda/cdd/
  3129.  
  3130.     [Preparata & Shamos] pp. 315-320
  3131.  
  3132.  
  3133.  
  3134. ----------------------------------------------------------------------
  3135. Subject 6.04: What are barycentric coordinates?
  3136.  
  3137.         Let p1, p2, p3 be the three vertices (corners) of a closed 
  3138.     triangle T (in 2D or 3D or any dimension).  Let t1, t2, t3 be real 
  3139.     numbers that sum to 1, with each between 0 and 1:  t1 + t2 + t3 = 1,
  3140.     0 <= ti <= 1.  Then the point p = t1*p1 + t2*p2 + t3*p3 lies in
  3141.     the plane of T and is inside T.  The numbers (t1,t2,t3) are called the
  3142.     barycentric coordinates of p with respect to T. They uniquely identify
  3143.     p.  They can be viewed as masses placed at the vertices whose
  3144.     center of gravity is p.
  3145.         For example, let p1=(0,0), p2=(1,0), p3=(3,2).  The
  3146.     barycentric coordinates (1/2,0,1/2) specify the point
  3147.     p = (0,0)/2 + 0*(1,0) + (3,2)/2 = (3/2,1), the midpoint of the p1-p3 
  3148.     edge.
  3149.         If p is joined to the three vertices, T is partitioned
  3150.     into three triangles.  Call them T1, T2, T3, with Ti not incident
  3151.     to pi.  The areas of these triangles Ti are proportional to the
  3152.     barycentric coordinates ti of p.  
  3153.     
  3154.     Reference:
  3155.     [Coxeter, Intro. to Geometry, p.217].
  3156.  
  3157. ----------------------------------------------------------------------
  3158. Subject 6.05: How can I generate a random point inside a triangle?
  3159.  
  3160.         Use barycentric coordinates (see item 53) .  Let A, B, C be the
  3161.     three vertices of your triangle.  Any point P inside can be expressed
  3162.     uniquely as P = aA + bB + cC, where a+b+c=1 and a,b,c are each >= 0.
  3163.     Knowing a and b permits you to calculate c=1-a-b.  So if you can
  3164.     generate two random numbers a and b, each in [0,1], such that
  3165.     their sum <=1, you've got a random point in your triangle.
  3166.         Generate random a and b independently and uniformly in [0,1] 
  3167.     (just divide the standard C rand() by its max value to get such a 
  3168.     random number.) If a+b>1, replace a by 1-a, b by 1-b.  Let c=1-a-b.  
  3169.     Then aA + bB + cC is uniformly distributed in triangle ABC: the reflection
  3170.     step a=1-a; b=1-b gives a point (a,b) uniformly distributed in the triangle
  3171.     (0,0)(1,0)(0,1), which is then mapped affinely to ABC.
  3172.     Now you have barycentric coordinates a,b,c.  Compute your point 
  3173.     P = aA + bB + cC.
  3174.  
  3175.     Reference: [Gems I], Turk, pp. 24-28, contains a similar but different
  3176.                method which requires a square root.
  3177.  
  3178. ----------------------------------------------------------------------
  3179. Subject 6.06: How do I evenly distribute N points on (tesselate) a sphere?
  3180.     
  3181.     "Evenly distributed" doesn't have a single definition.  There is a
  3182.     sense in which only the five Platonic solids achieve regular
  3183.     tesselations, as they are the only ones whose faces are regular
  3184.     and equal, with each vertex incident to the same number of faces.
  3185.     But generally "even distribution" focusses not so much on the
  3186.     induced tesselation, as it does on the distances and arrangement
  3187.     of the points/vertices.  For example, eight points can be placed
  3188.     on the sphere further away from one another than is achieved by
  3189.     the vertices of an inscribed cube: start with an inscribed cube,
  3190.     twist the top face 45 degrees about the north pole, and then
  3191.     move the top and bottom points along the surface towards the equator 
  3192.     a bit.
  3193.  
  3194.     The various definitions of "evenly distributed" lead into moderately 
  3195.     complex mathematics. This topic happens to be a FAQ on sci.math as well 
  3196.     as on comp.graphics.algorithms; a much more extensive and rigorous 
  3197.     discussion written by Dave Rusin can be found at:
  3198.     http://www.math.niu.edu/~rusin/known-math/95/sphere.faq
  3199.  
  3200.     A simple method of tesselating the sphere is repeated subdivision. An
  3201.     example starts with a unit octahedron. At each stage, every triangle in
  3202.     the tesselation is divided into 4 smaller triangles by creating 3 new
  3203.     vertices halfway along the edges. The new vertices are normalized,
  3204.     "pushing" them out to lie on the sphere. After N steps of subdivision,
  3205.     there will be 8 * 4^N triangles in the tesselation.
  3206.  
  3207.     A simple example of tesselation by subdivision is available at
  3208.        ftp://ftp.cs.unc.edu/pub/users/jon/sphere.c
  3209.  
  3210.     One frequently used definition of "evenness" is a distribution that 
  3211.     minimizes a 1/r potential energy function of all the points, so that in 
  3212.     a global sense points are as "far away" from each other as possible. 
  3213.     Starting from an arbitrary distribution, we can either use numerical 
  3214.     minimization algorithms or point repulsion, in which all the points 
  3215.     are considered to repel each other according to a 1/r^2 force law and 
  3216.     dynamics are simulated. The algorithm is run until some convergence 
  3217.     criterion is satisfied, and the resulting distribution is our approximation.
  3218.     
  3219.     Jon Leech developed code to do just this.  It is available at
  3220.     ftp://ftp.cs.unc.edu/pub/users/jon/points.tar.gz.
  3221.     See his README file for more information.  His distribution includes
  3222.     sample output files for various n < 300, which may be used off-the-shelf
  3223.     if that is all you need.
  3224.  
  3225.     Another method that is simpler than the above, but attains less
  3226.     uniformity, is based on spacing the points along a spiral that
  3227.     encircles the sphere.
  3228.     Code available from links at http://cs.smith.edu/~orourke/.
  3229.  
  3230. ----------------------------------------------------------------------
  3231. Subject 6.07: What are coordinates for the vertices of an icosohedron?
  3232.  
  3233.     Data on various polyhedra is available at
  3234.     http://cm.bell-labs.com/netlib/polyhedra/index.html, or
  3235.     http://netlib.bell-labs.com/netlib/polyhedra/index.html, or
  3236.     http://www.netlib.org/polyhedra/index.html
  3237.     For the icosahedron, the twelve vertices are:
  3238.  
  3239.       (+- 1, 0, +-t), (0, +-t, +-1), and (+-t, +-1, 0)
  3240.  
  3241.     where t = (sqrt(5)-1)/2, the golden ratio.
  3242.     Reference: "Beyond the Third Dimension" by Banchoff, p.168.
  3243.     
  3244.     
  3245. ----------------------------------------------------------------------
  3246. Subject 6.08: How do I generate random points on the surface of a sphere?
  3247.  
  3248.     There are several methods.  Perhaps the easiest to understand is
  3249.     the "rejection method":  generate random points in an origin-
  3250.     centered cube with opposite corners (r,r,r) and (-r,-r,-r).
  3251.     Reject any point p that falls outside of the sphere of radius r.
  3252.     Scale the vector to lie on the surface.  Because the cube to sphere
  3253.     volume ratio is pi/6, the average number of iterations before an
  3254.     acceptable vector is found is 6/pi = 1.90986.  This essentially
  3255.     doubles the effort, and makes this method slower than the "trig
  3256.     method."  A timing comparison conducted by Ken Sloan showed that
  3257.     the trig method runs in about 2/3's the time of the rejection method.
  3258.     He found that methods based on the use of normal distributions are
  3259.     twice as slow as the rejection method.
  3260.  
  3261.     The trig method.  This method works only in 3-space, but it is
  3262.     very fast.  It depends on the slightly counterintuitive fact (see
  3263.     proof below) that each of the three coordinates is uniformly
  3264.     distributed on [-1,1] (but the three are not independent,
  3265.     obviously).  Therefore, it suffices to choose one axis (Z, say)
  3266.     and generate a uniformly distributed value on that axis.  This
  3267.     constrains the chosen point to lie on a circle parallel to the
  3268.     X-Y plane, and the obvious trig method may be used to obtain the
  3269.     remaining coordinates.
  3270.     
  3271.         (a) Choose z uniformly distributed in [-1,1].
  3272.         (b) Choose t uniformly distributed on [0, 2*pi).
  3273.         (c) Let r = sqrt(1-z^2).
  3274.         (d) Let x = r * cos(t).
  3275.         (e) Let y = r * sin(t).
  3276.     
  3277.     This method uses uniform deviates (faster to generate than normal
  3278.     deviates), and no set of coordinates is ever rejected.  
  3279.     
  3280.     Here is a proof of correctness for the fact that the z-coordinate is
  3281.     uniformly distributed.  The proof also works for the x- and y-
  3282.     coordinates, but note that this works only in 3-space.
  3283.     
  3284.     The area of a surface of revolution in 3-space is given by
  3285.     
  3286.         A = 2 * pi * int_a^b f(x) * sqrt(1 + [f'(x}]^2) dx
  3287.     
  3288.     Consider a zone of a sphere of radius R.  Since we are integrating in
  3289.     the z direction, we have
  3290.     
  3291.         f(z) = sqrt(R^2 - z^2)
  3292.         f'(z) = -z / sqrt(R^2-z^2)
  3293.     
  3294.         1 + [f'(z)]^2 = r^2 / (R^2-z^2)
  3295.     
  3296.         A = 2 * pi * int_a^b sqrt(R^2-z^2) * R/sqrt(R^2-z^2) dz
  3297.           = 2 * pi * R int_a^b dz
  3298.           = 2 * pi * R * (b-a)
  3299.           = 2 * pi * R * h, 
  3300.           
  3301.     where h = b-a is the vertical height of the zone.  Notice how the integrand
  3302.     reduces to a constant.  The density is therefore uniform.
  3303.  
  3304.     Here is simple C code implementing the trig method:
  3305.     
  3306.     void SpherePoints(int n, double X[], double Y[], double Z[])
  3307.     {
  3308.       int i;
  3309.       double x, y, z, w, t;
  3310.     
  3311.       for( i=0; i< n; i++ ) {
  3312.         z = 2.0 * drand48() - 1.0;
  3313.         t = 2.0 * M_PI * drand48();
  3314.         w = sqrt( 1 - z*z );
  3315.         x = w * cos( t );
  3316.         y = w * sin( t );
  3317.         printf("i=%d:  x,y,z=\t%10.5lf\t%10.5lf\t%10.5lf\n", i, x,y,z);
  3318.         X[i] = x; Y[i] = y; Z[i] = z;
  3319.       }
  3320.     }
  3321.    
  3322.     A complete package is available at
  3323.     ftp://cs.smith.edu/pub/code/sphere.tar.gz (4K),
  3324.     reachable from http://cs.smith.edu/~orourke/ .
  3325.  
  3326.     If one wants to generate the random points in terms of longitude
  3327.     and latitude in degrees, these equations suffice:
  3328.       Longitude = random * 360 - 180
  3329.       Latitude  = [arccos( random * 2 - 1 )/pi ] * 180 - 90
  3330.  
  3331.  
  3332.     References:
  3333.     [Knuth, vol. 2]
  3334.     [Graphics Gems IV] "Uniform Random Rotations"
  3335.  
  3336. ----------------------------------------------------------------------
  3337. Subject 6.09: What are Plucker coordinates?
  3338.  
  3339.     A common convention is to write umlauted u as "ue", so you'll also see
  3340.     "Pluecker". Lines in 3D can easily be given by listing the coordinates of
  3341.     two distinct points, for a total of six numbers. Or, they can be given as
  3342.     the coordinates of two distinct planes, eight numbers. What's wrong with
  3343.     these? Nothing; but we can do better. Pluecker coordinates are, in a sense,
  3344.     halfway between these extremes, and can trivially generate either. Neither
  3345.     extreme is as efficient as Pluecker coordinates for computations.
  3346.  
  3347.     When Pluecker coordinates generalize to Grassmann coordinates, as laid
  3348.     out beautifully in [Hodge], Chapter VII, the determinant definition
  3349.     is clearly the one to use. But 3D lines can use a simpler definition.
  3350.     Take two distinct points on a line, say
  3351.  
  3352.         P = (Px,Py,Pz)
  3353.         Q = (Qx,Qy,Qz)
  3354.  
  3355.     Think of these as vectors from the origin, if you like. The Pluecker
  3356.     coordinates for the line are essentially
  3357.  
  3358.         U = P - Q
  3359.         V = P x Q
  3360.  
  3361.     Except for a scale factor, which we ignore, U and V do not depend on the
  3362.     specific P and Q! Cross products are perpendicular to their factors, so we
  3363.     always have U.V = 0. In [Stolfi] lines have orientation, so are the same
  3364.     only if their Pluecker coordinates are related by a positive scale factor.
  3365.  
  3366.     As determinants of homogeneous coordinates, begin with the 4x2 matrix
  3367.  
  3368.         [ Px  Qx ] row x
  3369.         [ Py  Qy ] row y
  3370.         [ Pz  Qz ] row z
  3371.         [  1   1 ] row w
  3372.  
  3373.     Define Pluecker coordinate Gij as the determinant of rows i and j, in that
  3374.     order. Notice that Giw = Pi - Qi, which is Ui. Now let (i,j,k) be a cyclic
  3375.     permutation of (x,y,z), namely (x,y,z) or (y,z,x) or (z,x,y), and notice
  3376.     that Gij = Vk. Determinants are anti-symmetric in the rows, so Gij = -Gji.
  3377.     Thus all possible Pluecker line coordinates are either zero (if i=j) or
  3378.     components of U or V, perhaps with a sign change. Taking the w component
  3379.     of a vector as 0, the determinant form will operate just as well on a
  3380.     point P and vector U as on two points. We can also begin with a 2x4 matrix
  3381.     whose rows are the coefficients of homogeneous plane equations, E.P=0,
  3382.     from which come dual coordinates G'ij. Now if (h,i,j,k) is an even
  3383.     permutation of (x,y,z,w), then Ghi = G'jk. (Just swap U and V!)
  3384.  
  3385.     Got Pluecker, want points? No problem. At least one component of U is
  3386.     non-zero, say Ui, which is Giw. Create homogeneous points Pj = Gjw + Gij,
  3387.     and Qj = Gij. (Don't expect the P and Q that we started with, and don't
  3388.     expect w=1.) Want plane equations? Let (i,j,k,w) be an even permutation of
  3389.     (x,y,z,w), so G'jk = Giw. Then create Eh = G'hk, and Fh = G'jh.
  3390.  
  3391.     Example: Begin with P = (2,4,8) and Q = (2,3,5). Then U = (0,1,3) and
  3392.     V = (-4,6,-2). The direct determinant forms are Gxw=0, Gyw=1, Gzw=3,
  3393.     Gyz=-4, Gzx=6, Gxy=-2, and the dual forms are G'yz=0, G'zx=1, G'xy=3,
  3394.     G'xw=-4, G'yw=6, G'zw=-2. Take Uz = Gzw = G'xy = 3 as a suitable non-zero
  3395.     element. Then two planes meeting in the line are
  3396.  
  3397.         (G'xy  G'yy  G'zy  G'wy).P = 0
  3398.         (G'xx  G'xy  G'xz  G'xw).P = 0
  3399.  
  3400.     That is, a point P is on the line if it satisfies both these equations:
  3401.  
  3402.         3 Px + 0 Py + 0 Pz - 6 Pw = 0
  3403.         0 Px + 3 Py - 1 Pz - 4 Pw = 0
  3404.  
  3405.     We can also easily determine if two lines meet, or if not, how they pass.
  3406.     If U1 and V1 are the coordinates of line 1, U2 and V2, of line 2, we look
  3407.     at the sign of U1.V2 + V1.U2. If it's zero, they meet. The determinant form
  3408.     reveals even permutations of (x,y,z,w):
  3409.         G1xw G2yz + G1yw G2zx + G1zw G2xy + G1yz G2xw + G1zx p2yw + G1xy p2zw
  3410.  
  3411.     Two oriented lines L1 and L2 can interact in three different ways:
  3412.     L1 might intersect L2, L1 might go clockwise around L2, or L1 might go
  3413.     counterclockwise around L2. Here are some examples:
  3414.  
  3415.              | L2            | L2            | L2
  3416.         L1   |          L1   |          L1   |
  3417.         -----+----->    ----------->    -----|----->
  3418.              |               |               |
  3419.              V               V               V
  3420.          intersect    counterclockwise   clockwise
  3421.              | L2            | L2            | L2
  3422.         L1   |          L1   |          L1   |
  3423.         <----+-----     <----|------    <-----------
  3424.              |               |               |
  3425.              V               V               V
  3426.  
  3427.     The first and second rows are just different views of the same lines,
  3428.     once from the "front" and once from the "back." Here's what they might
  3429.     look like if you look straight down line L2 (shown here as a dot).
  3430.  
  3431.         L1                               ---------->
  3432.         -----o---->     L1   o           L1   o
  3433.                         ---------->
  3434.          intersect    counterclockwise    clockwise
  3435.  
  3436.  
  3437.     The Pluecker coordinates of L1 and L2 give you a quick way to test
  3438.     which of the three it is.
  3439.  
  3440.         cw:   U1.V2 + V1.U2 < 0
  3441.         ccw:  U1.V2 + V1.U2 > 0
  3442.         thru: U1.V2 + V1.U2 = 0
  3443.  
  3444.     So why is this useful? Suppose you want to test if a ray intersects
  3445.     a triangle in 3-space. One way to do this is to represent the ray and
  3446.     the edges of the triangle with Pluecker coordinates. The ray hits the
  3447.     triangle if and only if it hits one of the triangle's edges, or it's
  3448.     "clockwise" from all three edges, or it's "counterclockwise" from all
  3449.     three edges. For example...
  3450.  
  3451.           o  _
  3452.           | |\              ...in this picture, the ray
  3453.           |   \            is oriented counterclockwise
  3454.         ------ \ -->        from all three edges, so it
  3455.           |     \         must intersect the triangle.
  3456.           v      \
  3457.           o-----> o
  3458.  
  3459.     Using Pluecker coordinates, ray shooting tests like this take only
  3460.     a few lines of code.
  3461.  
  3462.     Grassmann coordinates allow similar methods to be used for points,
  3463.     lines, planes, and so on, and in a space of any dimension. In the
  3464.     case of lines in 2D, only three coordinates are required:
  3465.  
  3466.         Gxw = Px-Qx     = -G'y
  3467.         Gyw = Py-Qy     =  G'x
  3468.         Gxy = PxQy-PyQx =  G'w
  3469.  
  3470.     Since P and Q are distinct, Giw is non-zero for i = x or y. Then
  3471.  
  3472.         (Gix,Giy,Giw) is a point P1 on L
  3473.         (Gxw,Gyw,Gww)+P1 is a point P2 on L
  3474.         (G'x,G'y,G'w).P = 0 is an equation for L
  3475.  
  3476.     Two lines in a 2D perspective plane always meet, perhaps in an
  3477.     ideal point (meaning they're parallel in ordinary 2D). Calling
  3478.     their homogeneous point of intersection P, the computation from
  3479.     Grassmann coordinates nicely illustrates the convenience. (See
  3480.     Subj 1.03, "How do I find intersections of 2 2D line segments?")
  3481.     For i = x,y,w, set
  3482.  
  3483.         Pi = G'j H'k ^╓ G'k H'j, where (i,j,k) is even
  3484.  
  3485.     See [Hodge] for a thorough discussion of the theory, [Stolfi] for
  3486.     a little theory with a concise implementation for low dimensions
  3487.     (see Subj. 0.04), and these articles for further discussion:
  3488.       J. Erickson, Pluecker Coordinates, Ray Tracing News 10(3) 1997,
  3489.       http://www.acm.org/tog/resources/RTNews/html/rtnv10n3.html.
  3490.  
  3491.       Ken Shoemake, Pluecker Coordinate Tutorial,
  3492.       Ray Tracing News 11(1) 1998,
  3493.       http://www.acm.org/tog/resources/RTNews/html/rtnv11n1.html.
  3494.  
  3495.  
  3496. ----------------------------------------------------------------------
  3497. Section 7. Contributors
  3498. ----------------------------------------------------------------------
  3499. Subject 7.01: How can you contribute to this FAQ?
  3500.  
  3501.     Send email to orourke@cs.smith.edu with your suggestions, possible
  3502.     topics, corrections, or pointers to information.  
  3503.  
  3504. ----------------------------------------------------------------------
  3505. Subject 7.02: Contributors.  Who made this all possible.
  3506.  
  3507.     [All email addresses now removed to protect the authors from spam.]
  3508.     Jens Alfke 
  3509.     Nina Amenta 
  3510.     Leen Ammeraal 
  3511.     Scott Anguish 
  3512.     Ian Ashdown 
  3513.     Barak 
  3514.     Brad Barber 
  3515.     James Beech 
  3516.     David Bouman 
  3517.     Paul Bourke 
  3518.     Lars Brinkhoff 
  3519.     Andrew Bromage 
  3520.     Brent Burley 
  3521.     R. Kevin Burton 
  3522.     Gene Caldwell 
  3523.     Ken Clarkson 
  3524.     Robert Day 
  3525.     Tamal Dey 
  3526.     Martin Dillon 
  3527.     Thomas Djafari 
  3528.     Dave Eberly 
  3529.     John Eickemeyer 
  3530.     John E (Edward) Ellis 
  3531.     Jeff Erickson 
  3532.     Ata Etemadi 
  3533.     Hugh Fisher 
  3534.     David N. Fogel 
  3535.     Arne K. Frick 
  3536.     Olexandr Frantchuk 
  3537.     Robert W. Fuentes 
  3538.     Komei Fukuda 
  3539.     William Gibbons 
  3540.     Normand GrΘgoire 
  3541.     Eric Haines 
  3542.     Jeff Hameluck 
  3543.     Sandy Harris 
  3544.     Luiz Henrique de Figueiredo 
  3545.     Steve Hollasch 
  3546.     Bill Jones 
  3547.     Richard Kinch 
  3548.     Craig Kolb 
  3549.     Uffe Kousgaard 
  3550.     Stefan Krause 
  3551.     Piyush Kumar 
  3552.     Steve Lamont 
  3553.     Ben Landon 
  3554.     Erik Larsen 
  3555.     Jon Leech 
  3556.     Michael V. Leonov 
  3557.     Sum Lin 
  3558.     Alan J. Livingston 
  3559.     Sebastien Loisel 
  3560.     Fritz Lott 
  3561.     Jacob Marner 
  3562.     Marc Christopher Martin 
  3563.     John McNamara 
  3564.     Samuel Murphy 
  3565.     Alan Murta 
  3566.     S. N. Muthukrishnan 
  3567.     Andrew Myers 
  3568.     David Nixon 
  3569.     Aaron Orenstein 
  3570.     Joseph O'Rourke 
  3571.     Samuel S. Paik 
  3572.     Leonidas Palios 
  3573.     Amitha Perera 
  3574.     Brian Peters 
  3575.     Lavoie Philippe 
  3576.     Christopher Phillips 
  3577.     Tom Plunket 
  3578.     Aaron Quigley 
  3579.     Rudi Bjxrn Rasmussen 
  3580.     Greg Roelofs 
  3581.     Christian von Roques 
  3582.     Dave Seaman 
  3583.     Jonathan R. Shewchuk 
  3584.     Rainer Michael Schmid 
  3585.     Klamer Schutte 
  3586.     Andrzej Serafin 
  3587.     ZhengYu Shan 
  3588.     James Sharman 
  3589.     Ken Shoemake 
  3590.     Jeff Somers 
  3591.     Jon Stone 
  3592.     Dan Sunday 
  3593.     Seth Teller 
  3594.     Saurabh Tendulkar 
  3595.     Yael "YoeL" Touboul 
  3596.     Anson Tsao 
  3597.     Bob van Manen 
  3598.     Remco Veltkamp 
  3599.     Jim Ward 
  3600.     Jason Weiler 
  3601.     Karsten Weiss 
  3602.     Stefan Wolfrum 
  3603.     Daniel S. Zwick 
  3604.  
  3605.     Previous Editors:
  3606.     Jon Stone 
  3607.     Anson Tsao 
  3608.  
  3609.