triangle.shar

用c语言编写的约束Delauney网格剖分程序源代码

/*  This program may be freely redistributed under the condition that the    */
/*    copyright notices (including this entire header and the copyright      */
/*    notice printed when the `-h' switch is selected) are not removed, and  */
/*    no compensation is received.  Private, research, and institutional     */
/*    use is free.  You may distribute modified versions of this code UNDER  */
/*    THE CONDITION THAT THIS CODE AND ANY MODIFICATIONS MADE TO IT IN THE   */
/*    SAME FILE REMAIN UNDER COPYRIGHT OF THE ORIGINAL AUTHOR, BOTH SOURCE   */
/*    AND OBJECT CODE ARE MADE FREELY AVAILABLE WITHOUT CHARGE, AND CLEAR    */
/*    NOTICE IS GIVEN OF THE MODIFICATIONS.  Distribution of this code as    */
/*    part of a commercial system is permissible ONLY BY DIRECT ARRANGEMENT  */
/*    WITH THE AUTHOR.  (If you are not directly supplying this code to a    */
/*    customer, and you are instead telling them how they can obtain it for  */
/*    free, then you are not required to make any arrangement with me.)      */
/*                                                                           */
/*  Hypertext instructions for Triangle are available on the Web at          */
/*                                                                           */
/*      http://www.cs.cmu.edu/~quake/triangle.html                           */
/*                                                                           */
/*  Some of the references listed below are marked [*].  These are available */
/*    for downloading from the Web page                                      */
/*                                                                           */
/*      http://www.cs.cmu.edu/~quake/triangle.research.html                  */
/*                                                                           */
/*  A paper discussing some aspects of Triangle is available.  See Jonathan  */
/*    Richard Shewchuk, "Triangle:  Engineering a 2D Quality Mesh Generator  */
/*    and Delaunay Triangulator," First Workshop on Applied Computational    */
/*    Geometry, ACM, May 1996.  [*]                                          */
/*                                                                           */
/*  Triangle was created as part of the Archimedes project in the School of  */
/*    Computer Science at Carnegie Mellon University.  Archimedes is a       */
/*    system for compiling parallel finite element solvers.  For further     */
/*    information, see Anja Feldmann, Omar Ghattas, John R. Gilbert, Gary L. */
/*    Miller, David R. O'Hallaron, Eric J. Schwabe, Jonathan R. Shewchuk,    */
/*    and Shang-Hua Teng, "Automated Parallel Solution of Unstructured PDE   */
/*    Problems."  To appear in Communications of the ACM, we hope.           */
/*                                                                           */
/*  The quality mesh generation algorithm is due to Jim Ruppert, "A          */
/*    Delaunay Refinement Algorithm for Quality 2-Dimensional Mesh           */
/*    Generation," Journal of Algorithms 18(3):548-585, May 1995.  [*]       */
/*                                                                           */
/*  My implementation of the divide-and-conquer and incremental Delaunay     */
/*    triangulation algorithms follows closely the presentation of Guibas    */
/*    and Stolfi, even though I use a triangle-based data structure instead  */
/*    of their quad-edge data structure.  (In fact, I originally implemented */
/*    Triangle using the quad-edge data structure, but switching to a        */
/*    triangle-based data structure sped Triangle by a factor of two.)  The  */
/*    mesh manipulation primitives and the two aforementioned Delaunay       */
/*    triangulation algorithms are described by Leonidas J. Guibas and Jorge */
/*    Stolfi, "Primitives for the Manipulation of General Subdivisions and   */
/*    the Computation of Voronoi Diagrams," ACM Transactions on Graphics     */
/*    4(2):74-123, April 1985.                                               */
/*                                                                           */
/*  Their O(n log n) divide-and-conquer algorithm is adapted from Der-Tsai   */
/*    Lee and Bruce J. Schachter, "Two Algorithms for Constructing the       */
/*    Delaunay Triangulation," International Journal of Computer and         */
/*    Information Science 9(3):219-242, 1980.  The idea to improve the       */
/*    divide-and-conquer algorithm by alternating between vertical and       */
/*    horizontal cuts was introduced by Rex A. Dwyer, "A Faster Divide-and-  */
/*    Conquer Algorithm for Constructing Delaunay Triangulations,"           */
/*    Algorithmica 2(2):137-151, 1987.                                       */
/*                                                                           */
/*  The incremental insertion algorithm was first proposed by C. L. Lawson,  */
/*    "Software for C1 Surface Interpolation," in Mathematical Software III, */
/*    John R. Rice, editor, Academic Press, New York, pp. 161-194, 1977.     */
/*    For point location, I use the algorithm of Ernst P. Mucke, Isaac       */
/*    Saias, and Binhai Zhu, "Fast Randomized Point Location Without         */
/*    Preprocessing in Two- and Three-dimensional Delaunay Triangulations,"  */
/*    Proceedings of the Twelfth Annual Symposium on Computational Geometry, */
/*    ACM, May 1996.  [*]  If I were to randomize the order of point         */
/*    insertion (I currently don't bother), their result combined with the   */
/*    result of Leonidas J. Guibas, Donald E. Knuth, and Micha Sharir,       */
/*    "Randomized Incremental Construction of Delaunay and Voronoi           */
/*    Diagrams," Algorithmica 7(4):381-413, 1992, would yield an expected    */
/*    O(n^{4/3}) bound on running time.                                      */
/*                                                                           */
/*  The O(n log n) sweepline Delaunay triangulation algorithm is taken from  */
/*    Steven Fortune, "A Sweepline Algorithm for Voronoi Diagrams",          */
/*    Algorithmica 2(2):153-174, 1987.  A random sample of edges on the      */
/*    boundary of the triangulation are maintained in a splay tree for the   */
/*    purpose of point location.  Splay trees are described by Daniel        */
/*    Dominic Sleator and Robert Endre Tarjan, "Self-Adjusting Binary Search */
/*    Trees," Journal of the ACM 32(3):652-686, July 1985.                   */
/*                                                                           */
/*  The algorithms for exact computation of the signs of determinants are    */
/*    described in Jonathan Richard Shewchuk, "Adaptive Precision Floating-  */
/*    Point Arithmetic and Fast Robust Geometric Predicates," Technical      */
/*    Report CMU-CS-96-140, School of Computer Science, Carnegie Mellon      */
/*    University, Pittsburgh, Pennsylvania, May 1996.  [*]  (Submitted to    */
/*    Discrete & Computational Geometry.)  An abbreviated version appears as */
/*    Jonathan Richard Shewchuk, "Robust Adaptive Floating-Point Geometric   */
/*    Predicates," Proceedings of the Twelfth Annual Symposium on Computa-   */
/*    tional Geometry, ACM, May 1996.  [*]  Many of the ideas for my exact   */
/*    arithmetic routines originate with Douglas M. Priest, "Algorithms for  */
/*    Arbitrary Precision Floating Point Arithmetic," Tenth Symposium on     */
/*    Computer Arithmetic, 132-143, IEEE Computer Society Press, 1991.  [*]  */
/*    Many of the ideas for the correct evaluation of the signs of           */
/*    determinants are taken from Steven Fortune and Christopher J. Van Wyk, */
/*    "Efficient Exact Arithmetic for Computational Geometry," Proceedings   */
/*    of the Ninth Annual Symposium on Computational Geometry, ACM,          */
/*    pp. 163-172, May 1993, and from Steven Fortune, "Numerical Stability   */
/*    of Algorithms for 2D Delaunay Triangulations," International Journal   */
/*    of Computational Geometry & Applications 5(1-2):193-213, March-June    */
/*    1995.                                                                  */
/*                                                                           */
/*  For definitions of and results involving Delaunay triangulations,        */
/*    constrained and conforming versions thereof, and other aspects of      */
/*    triangular mesh generation, see the excellent survey by Marshall Bern  */
/*    and David Eppstein, "Mesh Generation and Optimal Triangulation," in    */
/*    Computing and Euclidean Geometry, Ding-Zhu Du and Frank Hwang,         */
/*    editors, World Scientific, Singapore, pp. 23-90, 1992.                 */
/*                                                                           */
/*  The time for incrementally adding PSLG (planar straight line graph)      */
/*    segments to create a constrained Delaunay triangulation is probably    */
/*    O(n^2) per segment in the worst case and O(n) per edge in the common   */
/*    case, where n is the number of triangles that intersect the segment    */
/*    before it is inserted.  This doesn't count point location, which can   */
/*    be much more expensive.  (This note does not apply to conforming       */
/*    Delaunay triangulations, for which a different method is used to       */
/*    insert segments.)                                                      */
/*                                                                           */
/*  The time for adding segments to a conforming Delaunay triangulation is   */
/*    not clear, but does not depend upon n alone.  In some cases, very      */
/*    small features (like a point lying next to a segment) can cause a      */
/*    single segment to be split an arbitrary number of times.  Of course,   */
/*    floating-point precision is a practical barrier to how much this can   */
/*    happen.                                                                */
/*                                                                           */
/*  The time for deleting a point from a Delaunay triangulation is O(n^2) in */
/*    the worst case and O(n) in the common case, where n is the degree of   */
/*    the point being deleted.  I could improve this to expected O(n) time   */
/*    by "inserting" the neighboring vertices in random order, but n is      */
/*    usually quite small, so it's not worth the bother.  (The O(n) time     */
/*    for random insertion follows from L. Paul Chew, "Building Voronoi      */
/*    Diagrams for Convex Polygons in Linear Expected Time," Technical       */
/*    Report PCS-TR90-147, Department of Mathematics and Computer Science,   */
/*    Dartmouth College, 1990.                                               */
/*                                                                           */
/*  Ruppert's Delaunay refinement algorithm typically generates triangles    */
/*    at a linear rate (constant time per triangle) after the initial        */
/*    triangulation is formed.  There may be pathological cases where more   */
/*    time is required, but these never arise in practice.                   */
/*                                                                           */
/*  The segment intersection formulae are straightforward.  If you want to   */
/*    see them derived, see Franklin Antonio.  "Faster Line Segment          */
/*    Intersection."  In Graphics Gems III (David Kirk, editor), pp. 199-    */
/*    202.  Academic Press, Boston, 1992.                                    */
/*                                                                           */
/*  If you make any improvements to this code, please please please let me   */
/*    know, so that I may obtain the improvements.  Even if you don't change */
/*    the code, I'd still love to hear what it's being used for.             */
/*                                                                           */
/*  Disclaimer:  Neither I nor Carnegie Mellon warrant this code in any way  */
/*    whatsoever.  This code is provided "as-is".  Use at your own risk.     */
/*                                                                           */
/*****************************************************************************/

/* For single precision (which will save some memory and reduce paging),     */
/*   define the symbol SINGLE by using the -DSINGLE compiler switch or by    */
/*   writing "#define SINGLE" below.                                         */
/*                                                                           */
/* For double precision (which will allow you to refine meshes to a smaller  */
/*   edge length), leave SINGLE undefined.                                   */
/*                                                                           */
/* Double precision uses more memory, but improves the resolution of the     */
/*   meshes you can generate with Triangle.  It also reduces the likelihood  */
/*   of a floating exception due to overflow.  Finally, it is much faster    */
/*   than single precision on 64-bit architectures like the DEC Alpha.  I    */
/*   recommend double precision unless you want to generate a mesh for which */
/*   you do not have enough memory.                                          */

/* #define SINGLE */

#ifdef SINGLE
#define REAL float
#else /* not SINGLE */
#define REAL double
#endif /* not SINGLE */

/* If yours is not a Unix system, define the NO_TIMER compiler switch to     */
/*   remove the Unix-specific timing code.                                   */

/* #define NO_TIMER */

/* To insert lots of self-checks for internal errors, define the SELF_CHECK  */
/*   symbol.  This will slow down the program significantly.  It is best to  */
/*   define the symbol using the -DSELF_CHECK compiler switch, but you could */
/*   write "#define SELF_CHECK" below.  If you are modifying this code, I    */
/*   recommend you turn self-checks on.                                      */

/* #define SELF_CHECK */

/* To compile Triangle as a callable object library (triangle.o), define the */
/*   TRILIBRARY symbol.  Read the file triangle.h for details on how to call */
/*   the procedure triangulate() that results.                               */

/* #define TRILIBRARY */

/* It is possible to generate a smaller version of Triangle using one or     */
/*   both of the following symbols.  Define the REDUCED symbol to eliminate  */
/*   all features that are primarily of research interest; specifically, the */
/*   -i, -F, -s, and -C switches.  Define the CDT_ONLY symbol to eliminate   */
/*   all meshing algorithms above and beyond constrained Delaunay            */
/*   triangulation; specifically, the -r, -q, -a, -S, and -s switches.       */
/*   These reductions are most likely to be useful when generating an object */
/*   library (triangle.o) by defining the TRILIBRARY symbol.                 */

/* #define REDUCED */
/* #define CDT_ONLY */

/* On some machines, the exact arithmetic routines might be defeated by the  */
/*   use of internal extended precision floating-point registers.  Sometimes */
/*   this problem can be fixed by defining certain values to be volatile,    */
/*   thus forcing them to be stored to memory and rounded off.  This isn't   */
/*   a great solution, though, as it slows Triangle down.                    */
/*                                                                           */
/* To try this out, write "#define INEXACT volatile" below.  Normally,       */
/*   however, INEXACT should be defined to be nothing.  ("#define INEXACT".) */

#define INEXACT /* Nothing */
/* #define INEXACT volatile */

/* Maximum number of characters in a file name (including the null).         */

#define FILENAMESIZE 512

/* Maximum number of characters in a line read from a file (including the    */
/*   null).                                                                  */

#define INPUTLINESIZE 512

/* For efficiency, a variety of data structures are allocated in bulk.  The  */