PageRenderTime 71ms CodeModel.GetById 25ms RepoModel.GetById 0ms app.codeStats 2ms

/contrib/gtkgensurf/triangle.c

https://gitlab.com/illwieckz/netradiant
C | 12631 lines | 8857 code | 874 blank | 2900 comment | 1792 complexity | 8758741fd6acfa480a2e2328d8de068f MD5 | raw file
Possible License(s): GPL-2.0, LGPL-2.1, LGPL-2.0

Large files files are truncated, but you can click here to view the full file

    

  1. #define ANSI_DECLARATORS
    2. 

  2. /*****************************************************************************/
    3. 

  3. /*                                                                           */
    4. 

  4. /*      888888888        ,o,                          / 888                  */
    5. 

  5. /*         888    88o88o  "    o8888o  88o8888o o88888o 888  o88888o         */
    6. 

  6. /*         888    888    888       88b 888  888 888 888 888 d888  88b        */
    7. 

  7. /*         888    888    888  o88^o888 888  888 "88888" 888 8888oo888        */
    8. 

  8. /*         888    888    888 C888  888 888  888  /      888 q888             */
    9. 

  9. /*         888    888    888  "88o^888 888  888 Cb      888  "88oooo"        */
    10. 

  10. /*                                              "8oo8D                       */
    11. 

  11. /*                                                                           */
    12. 

  12. /*  A Two-Dimensional Quality Mesh Generator and Delaunay Triangulator.      */
    13. 

  13. /*  (triangle.c)                                                             */
    14. 

  14. /*                                                                           */
    15. 

  15. /*  Version 1.3                                                              */
    16. 

  16. /*  July 19, 1996                                                            */
    17. 

  17. /*                                                                           */
    18. 

  18. /*  Copyright 1996                                                           */
    19. 

  19. /*  Jonathan Richard Shewchuk                                                */
    20. 

  20. /*  School of Computer Science                                               */
    21. 

  21. /*  Carnegie Mellon University                                               */
    22. 

  22. /*  5000 Forbes Avenue                                                       */
    23. 

  23. /*  Pittsburgh, Pennsylvania  15213-3891                                     */
    24. 

  24. /*  jrs@cs.cmu.edu                                                           */
    25. 

  25. /*                                                                           */
    26. 

  26. /*  This program may be freely redistributed under the condition that the    */
    27. 

  27. /*    copyright notices (including this entire header and the copyright      */
    28. 

  28. /*    notice printed when the `-h' switch is selected) are not removed, and  */
    29. 

  29. /*    no compensation is received.  Private, research, and institutional     */
    30. 

  30. /*    use is free.  You may distribute modified versions of this code UNDER  */
    31. 

  31. /*    THE CONDITION THAT THIS CODE AND ANY MODIFICATIONS MADE TO IT IN THE   */
    32. 

  32. /*    SAME FILE REMAIN UNDER COPYRIGHT OF THE ORIGINAL AUTHOR, BOTH SOURCE   */
    33. 

  33. /*    AND OBJECT CODE ARE MADE FREELY AVAILABLE WITHOUT CHARGE, AND CLEAR    */
    34. 

  34. /*    NOTICE IS GIVEN OF THE MODIFICATIONS.  Distribution of this code as    */
    35. 

  35. /*    part of a commercial system is permissible ONLY BY DIRECT ARRANGEMENT  */
    36. 

  36. /*    WITH THE AUTHOR.  (If you are not directly supplying this code to a    */
    37. 

  37. /*    customer, and you are instead telling them how they can obtain it for  */
    38. 

  38. /*    free, then you are not required to make any arrangement with me.)      */
    39. 

  39. /*                                                                           */
    40. 

  40. /*  Hypertext instructions for Triangle are available on the Web at          */
    41. 

  41. /*                                                                           */
    42. 

  42. /*      http://www.cs.cmu.edu/~quake/triangle.html                           */
    43. 

  43. /*                                                                           */
    44. 

  44. /*  Some of the references listed below are marked [*].  These are available */
    45. 

  45. /*    for downloading from the Web page                                      */
    46. 

  46. /*                                                                           */
    47. 

  47. /*      http://www.cs.cmu.edu/~quake/triangle.research.html                  */
    48. 

  48. /*                                                                           */
    49. 

  49. /*  A paper discussing some aspects of Triangle is available.  See Jonathan  */
    50. 

  50. /*    Richard Shewchuk, "Triangle:  Engineering a 2D Quality Mesh Generator  */
    51. 

  51. /*    and Delaunay Triangulator," First Workshop on Applied Computational    */
    52. 

  52. /*    Geometry, ACM, May 1996.  [*]                                          */
    53. 

  53. /*                                                                           */
    54. 

  54. /*  Triangle was created as part of the Archimedes project in the School of  */
    55. 

  55. /*    Computer Science at Carnegie Mellon University.  Archimedes is a       */
    56. 

  56. /*    system for compiling parallel finite element solvers.  For further     */
    57. 

  57. /*    information, see Anja Feldmann, Omar Ghattas, John R. Gilbert, Gary L. */
    58. 

  58. /*    Miller, David R. O'Hallaron, Eric J. Schwabe, Jonathan R. Shewchuk,    */
    59. 

  59. /*    and Shang-Hua Teng, "Automated Parallel Solution of Unstructured PDE   */
    60. 

  60. /*    Problems."  To appear in Communications of the ACM, we hope.           */
    61. 

  61. /*                                                                           */
    62. 

  62. /*  The quality mesh generation algorithm is due to Jim Ruppert, "A          */
    63. 

  63. /*    Delaunay Refinement Algorithm for Quality 2-Dimensional Mesh           */
    64. 

  64. /*    Generation," Journal of Algorithms 18(3):548-585, May 1995.  [*]       */
    65. 

  65. /*                                                                           */
    66. 

  66. /*  My implementation of the divide-and-conquer and incremental Delaunay     */
    67. 

  67. /*    triangulation algorithms follows closely the presentation of Guibas    */
    68. 

  68. /*    and Stolfi, even though I use a triangle-based data structure instead  */
    69. 

  69. /*    of their quad-edge data structure.  (In fact, I originally implemented */
    70. 

  70. /*    Triangle using the quad-edge data structure, but switching to a        */
    71. 

  71. /*    triangle-based data structure sped Triangle by a factor of two.)  The  */
    72. 

  72. /*    mesh manipulation primitives and the two aforementioned Delaunay       */
    73. 

  73. /*    triangulation algorithms are described by Leonidas J. Guibas and Jorge */
    74. 

  74. /*    Stolfi, "Primitives for the Manipulation of General Subdivisions and   */
    75. 

  75. /*    the Computation of Voronoi Diagrams," ACM Transactions on Graphics     */
    76. 

  76. /*    4(2):74-123, April 1985.                                               */
    77. 

  77. /*                                                                           */
    78. 

  78. /*  Their O(n log n) divide-and-conquer algorithm is adapted from Der-Tsai   */
    79. 

  79. /*    Lee and Bruce J. Schachter, "Two Algorithms for Constructing the       */
    80. 

  80. /*    Delaunay Triangulation," International Journal of Computer and         */
    81. 

  81. /*    Information Science 9(3):219-242, 1980.  The idea to improve the       */
    82. 

  82. /*    divide-and-conquer algorithm by alternating between vertical and       */
    83. 

  83. /*    horizontal cuts was introduced by Rex A. Dwyer, "A Faster Divide-and-  */
    84. 

  84. /*    Conquer Algorithm for Constructing Delaunay Triangulations,"           */
    85. 

  85. /*    Algorithmica 2(2):137-151, 1987.                                       */
    86. 

  86. /*                                                                           */
    87. 

  87. /*  The incremental insertion algorithm was first proposed by C. L. Lawson,  */
    88. 

  88. /*    "Software for C1 Surface Interpolation," in Mathematical Software III, */
    89. 

  89. /*    John R. Rice, editor, Academic Press, New York, pp. 161-194, 1977.     */
    90. 

  90. /*    For point location, I use the algorithm of Ernst P. Mucke, Isaac       */
    91. 

  91. /*    Saias, and Binhai Zhu, "Fast Randomized Point Location Without         */
    92. 

  92. /*    Preprocessing in Two- and Three-dimensional Delaunay Triangulations,"  */
    93. 

  93. /*    Proceedings of the Twelfth Annual Symposium on Computational Geometry, */
    94. 

  94. /*    ACM, May 1996.  [*]  If I were to randomize the order of point         */
    95. 

  95. /*    insertion (I currently don't bother), their result combined with the   */
    96. 

  96. /*    result of Leonidas J. Guibas, Donald E. Knuth, and Micha Sharir,       */
    97. 

  97. /*    "Randomized Incremental Construction of Delaunay and Voronoi           */
    98. 

  98. /*    Diagrams," Algorithmica 7(4):381-413, 1992, would yield an expected    */
    99. 

  99. /*    O(n^{4/3}) bound on running time.                                      */
    100. 

  100. /*                                                                           */
    101. 

  101. /*  The O(n log n) sweepline Delaunay triangulation algorithm is taken from  */
    102. 

  102. /*    Steven Fortune, "A Sweepline Algorithm for Voronoi Diagrams",          */
    103. 

  103. /*    Algorithmica 2(2):153-174, 1987.  A random sample of edges on the      */
    104. 

  104. /*    boundary of the triangulation are maintained in a splay tree for the   */
    105. 

  105. /*    purpose of point location.  Splay trees are described by Daniel        */
    106. 

  106. /*    Dominic Sleator and Robert Endre Tarjan, "Self-Adjusting Binary Search */
    107. 

  107. /*    Trees," Journal of the ACM 32(3):652-686, July 1985.                   */
    108. 

  108. /*                                                                           */
    109. 

  109. /*  The algorithms for exact computation of the signs of determinants are    */
    110. 

  110. /*    described in Jonathan Richard Shewchuk, "Adaptive Precision Floating-  */
    111. 

  111. /*    Point Arithmetic and Fast Robust Geometric Predicates," Technical      */
    112. 

  112. /*    Report CMU-CS-96-140, School of Computer Science, Carnegie Mellon      */
    113. 

  113. /*    University, Pittsburgh, Pennsylvania, May 1996.  [*]  (Submitted to    */
    114. 

  114. /*    Discrete & Computational Geometry.)  An abbreviated version appears as */
    115. 

  115. /*    Jonathan Richard Shewchuk, "Robust Adaptive Floating-Point Geometric   */
    116. 

  116. /*    Predicates," Proceedings of the Twelfth Annual Symposium on Computa-   */
    117. 

  117. /*    tional Geometry, ACM, May 1996.  [*]  Many of the ideas for my exact   */
    118. 

  118. /*    arithmetic routines originate with Douglas M. Priest, "Algorithms for  */
    119. 

  119. /*    Arbitrary Precision Floating Point Arithmetic," Tenth Symposium on     */
    120. 

  120. /*    Computer Arithmetic, 132-143, IEEE Computer Society Press, 1991.  [*]  */
    121. 

  121. /*    Many of the ideas for the correct evaluation of the signs of           */
    122. 

  122. /*    determinants are taken from Steven Fortune and Christopher J. Van Wyk, */
    123. 

  123. /*    "Efficient Exact Arithmetic for Computational Geometry," Proceedings   */
    124. 

  124. /*    of the Ninth Annual Symposium on Computational Geometry, ACM,          */
    125. 

  125. /*    pp. 163-172, May 1993, and from Steven Fortune, "Numerical Stability   */
    126. 

  126. /*    of Algorithms for 2D Delaunay Triangulations," International Journal   */
    127. 

  127. /*    of Computational Geometry & Applications 5(1-2):193-213, March-June    */
    128. 

  128. /*    1995.                                                                  */
    129. 

  129. /*                                                                           */
    130. 

  130. /*  For definitions of and results involving Delaunay triangulations,        */
    131. 

  131. /*    constrained and conforming versions thereof, and other aspects of      */
    132. 

  132. /*    triangular mesh generation, see the excellent survey by Marshall Bern  */
    133. 

  133. /*    and David Eppstein, "Mesh Generation and Optimal Triangulation," in    */
    134. 

  134. /*    Computing and Euclidean Geometry, Ding-Zhu Du and Frank Hwang,         */
    135. 

  135. /*    editors, World Scientific, Singapore, pp. 23-90, 1992.                 */
    136. 

  136. /*                                                                           */
    137. 

  137. /*  The time for incrementally adding PSLG (planar straight line graph)      */
    138. 

  138. /*    segments to create a constrained Delaunay triangulation is probably    */
    139. 

  139. /*    O(n^2) per segment in the worst case and O(n) per edge in the common   */
    140. 

  140. /*    case, where n is the number of triangles that intersect the segment    */
    141. 

  141. /*    before it is inserted.  This doesn't count point location, which can   */
    142. 

  142. /*    be much more expensive.  (This note does not apply to conforming       */
    143. 

  143. /*    Delaunay triangulations, for which a different method is used to       */
    144. 

  144. /*    insert segments.)                                                      */
    145. 

  145. /*                                                                           */
    146. 

  146. /*  The time for adding segments to a conforming Delaunay triangulation is   */
    147. 

  147. /*    not clear, but does not depend upon n alone.  In some cases, very      */
    148. 

  148. /*    small features (like a point lying next to a segment) can cause a      */
    149. 

  149. /*    single segment to be split an arbitrary number of times.  Of course,   */
    150. 

  150. /*    floating-point precision is a practical barrier to how much this can   */
    151. 

  151. /*    happen.                                                                */
    152. 

  152. /*                                                                           */
    153. 

  153. /*  The time for deleting a point from a Delaunay triangulation is O(n^2) in */
    154. 

  154. /*    the worst case and O(n) in the common case, where n is the degree of   */
    155. 

  155. /*    the point being deleted.  I could improve this to expected O(n) time   */
    156. 

  156. /*    by "inserting" the neighboring vertices in random order, but n is      */
    157. 

  157. /*    usually quite small, so it's not worth the bother.  (The O(n) time     */
    158. 

  158. /*    for random insertion follows from L. Paul Chew, "Building Voronoi      */
    159. 

  159. /*    Diagrams for Convex Polygons in Linear Expected Time," Technical       */
    160. 

  160. /*    Report PCS-TR90-147, Department of Mathematics and Computer Science,   */
    161. 

  161. /*    Dartmouth College, 1990.                                               */
    162. 

  162. /*                                                                           */
    163. 

  163. /*  Ruppert's Delaunay refinement algorithm typically generates triangles    */
    164. 

  164. /*    at a linear rate (constant time per triangle) after the initial        */
    165. 

  165. /*    triangulation is formed.  There may be pathological cases where more   */
    166. 

  166. /*    time is required, but these never arise in practice.                   */
    167. 

  167. /*                                                                           */
    168. 

  168. /*  The segment intersection formulae are straightforward.  If you want to   */
    169. 

  169. /*    see them derived, see Franklin Antonio.  "Faster Line Segment          */
    170. 

  170. /*    Intersection."  In Graphics Gems III (David Kirk, editor), pp. 199-    */
    171. 

  171. /*    202.  Academic Press, Boston, 1992.                                    */
    172. 

  172. /*                                                                           */
    173. 

  173. /*  If you make any improvements to this code, please please please let me   */
    174. 

  174. /*    know, so that I may obtain the improvements.  Even if you don't change */
    175. 

  175. /*    the code, I'd still love to hear what it's being used for.             */
    176. 

  176. /*                                                                           */
    177. 

  177. /*  Disclaimer:  Neither I nor Carnegie Mellon warrant this code in any way  */
    178. 

  178. /*    whatsoever.  This code is provided "as-is".  Use at your own risk.     */
    179. 

  179. /*                                                                           */
    180. 

  180. /*****************************************************************************/
    181. 

  181. 

  182. /* For single precision (which will save some memory and reduce paging),     */
    183. 

  183. /*   define the symbol SINGLE by using the -DSINGLE compiler switch or by    */
    184. 

  184. /*   writing "#define SINGLE" below.                                         */
    185. 

  185. /*                                                                           */
    186. 

  186. /* For double precision (which will allow you to refine meshes to a smaller  */
    187. 

  187. /*   edge length), leave SINGLE undefined.                                   */
    188. 

  188. /*                                                                           */
    189. 

  189. /* Double precision uses more memory, but improves the resolution of the     */
    190. 

  190. /*   meshes you can generate with Triangle.  It also reduces the likelihood  */
    191. 

  191. /*   of a floating exception due to overflow.  Finally, it is much faster    */
    192. 

  192. /*   than single precision on 64-bit architectures like the DEC Alpha.  I    */
    193. 

  193. /*   recommend double precision unless you want to generate a mesh for which */
    194. 

  194. /*   you do not have enough memory.                                          */
    195. 

  195. 

  196. #define SINGLE
    197. 

  197. 

  198. #ifdef SINGLE
    199. 

  199. #define REAL float
    200. 

  200. #else /* not SINGLE */
    201. 

  201. #define REAL double
    202. 

  202. #endif /* not SINGLE */
    203. 

  203. 

  204. /* If yours is not a Unix system, define the NO_TIMER compiler switch to     */
    205. 

  205. /*   remove the Unix-specific timing code.                                   */
    206. 

  206. 

  207. #define NO_TIMER
    208. 

  208. 

  209. /* To insert lots of self-checks for internal errors, define the SELF_CHECK  */
    210. 

  210. /*   symbol.  This will slow down the program significantly.  It is best to  */
    211. 

  211. /*   define the symbol using the -DSELF_CHECK compiler switch, but you could */
    212. 

  212. /*   write "#define SELF_CHECK" below.  If you are modifying this code, I    */
    213. 

  213. /*   recommend you turn self-checks on.                                      */
    214. 

  214. 

  215. /* #define SELF_CHECK */
    216. 

  216. 

  217. /* To compile Triangle as a callable object library (triangle.o), define the */
    218. 

  218. /*   TRILIBRARY symbol.  Read the file triangle.h for details on how to call */
    219. 

  219. /*   the procedure triangulate() that results.                               */
    220. 

  220. 

  221. #define TRILIBRARY
    222. 

  222. 

  223. /* It is possible to generate a smaller version of Triangle using one or     */
    224. 

  224. /*   both of the following symbols.  Define the REDUCED symbol to eliminate  */
    225. 

  225. /*   all features that are primarily of research interest; specifically, the */
    226. 

  226. /*   -i, -F, -s, and -C switches.  Define the CDT_ONLY symbol to eliminate   */
    227. 

  227. /*   all meshing algorithms above and beyond constrained Delaunay            */
    228. 

  228. /*   triangulation; specifically, the -r, -q, -a, -S, and -s switches.       */
    229. 

  229. /*   These reductions are most likely to be useful when generating an object */
    230. 

  230. /*   library (triangle.o) by defining the TRILIBRARY symbol.                 */
    231. 

  231. 

  232. #define REDUCED
    233. 

  233. #define CDT_ONLY
    234. 

  234. 

  235. /* On some machines, the exact arithmetic routines might be defeated by the  */
    236. 

  236. /*   use of internal extended precision floating-point registers.  Sometimes */
    237. 

  237. /*   this problem can be fixed by defining certain values to be volatile,    */
    238. 

  238. /*   thus forcing them to be stored to memory and rounded off.  This isn't   */
    239. 

  239. /*   a great solution, though, as it slows Triangle down.                    */
    240. 

  240. /*                                                                           */
    241. 

  241. /* To try this out, write "#define INEXACT volatile" below.  Normally,       */
    242. 

  242. /*   however, INEXACT should be defined to be nothing.  ("#define INEXACT".) */
    243. 

  243. 

  244. #define INEXACT /* Nothing */
    245. 

  245. /* #define INEXACT volatile */
    246. 

  246. 

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

  248. 

  249. #define FILENAMESIZE 512
    250. 

  250. 

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

  252. /*   null).                                                                  */
    253. 

  253. 

  254. #define INPUTLINESIZE 512
    255. 

  255. 

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

  257. /*   following constants determine how many of each structure is allocated   */
    258. 

  258. /*   at once.                                                                */
    259. 

  259. 

  260. #define TRIPERBLOCK 4092           /* Number of triangles allocated at once. */
    261. 

  261. #define SHELLEPERBLOCK 508       /* Number of shell edges allocated at once. */
    262. 

  262. #define POINTPERBLOCK 4092            /* Number of points allocated at once. */
    263. 

  263. #define VIRUSPERBLOCK 1020   /* Number of virus triangles allocated at once. */
    264. 

  264. /* Number of encroached segments allocated at once. */
    265. 

  265. #define BADSEGMENTPERBLOCK 252
    266. 

  266. /* Number of skinny triangles allocated at once. */
    267. 

  267. #define BADTRIPERBLOCK 4092
    268. 

  268. /* Number of splay tree nodes allocated at once. */
    269. 

  269. #define SPLAYNODEPERBLOCK 508
    270. 

  270. 

  271. /* The point marker DEADPOINT is an arbitrary number chosen large enough to  */
    272. 

  272. /*   (hopefully) not conflict with user boundary markers.  Make sure that it */
    273. 

  273. /*   is small enough to fit into your machine's integer size.                */
    274. 

  274. 

  275. #define DEADPOINT -1073741824
    276. 

  276. 

  277. /* The next line is used to outsmart some very stupid compilers.  If your    */
    278. 

  278. /*   compiler is smarter, feel free to replace the "int" with "void".        */
    279. 

  279. /*   Not that it matters.                                                    */
    280. 

  280. 

  281. #define VOID int
    282. 

  282. 

  283. /* Two constants for algorithms based on random sampling.  Both constants    */
    284. 

  284. /*   have been chosen empirically to optimize their respective algorithms.   */
    285. 

  285. 

  286. /* Used for the point location scheme of Mucke, Saias, and Zhu, to decide    */
    287. 

  287. /*   how large a random sample of triangles to inspect.                      */
    288. 

  288. #define SAMPLEFACTOR 11
    289. 

  289. /* Used in Fortune's sweepline Delaunay algorithm to determine what fraction */
    290. 

  290. /*   of boundary edges should be maintained in the splay tree for point      */
    291. 

  291. /*   location on the front.                                                  */
    292. 

  292. #define SAMPLERATE 10
    293. 

  293. 

  294. /* A number that speaks for itself, every kissable digit.                    */
    295. 

  295. 

  296. #define PI 3.141592653589793238462643383279502884197169399375105820974944592308
    297. 

  297. 

  298. /* Another fave.                                                             */
    299. 

  299. 

  300. #define SQUAREROOTTWO 1.4142135623730950488016887242096980785696718753769480732
    301. 

  301. 

  302. /* And here's one for those of you who are intimidated by math.              */
    303. 

  303. 

  304. #define ONETHIRD 0.333333333333333333333333333333333333333333333333333333333333
    305. 

  305. 

  306. #include <stdio.h>
    307. 

  307. #include <string.h>
    308. 

  308. #include <math.h>
    309. 

  309. #ifndef NO_TIMER
    310. 

  310. #include <sys/time.h>
    311. 

  311. #endif /* NO_TIMER */
    312. 

  312. #ifdef TRILIBRARY
    313. 

  313. #include "triangle.h"
    314. 

  314. #endif /* TRILIBRARY */
    315. 

  315. 

  316. /* The following obscenity seems to be necessary to ensure that this program */
    317. 

  317. /* will port to Dec Alphas running OSF/1, because their stdio.h file commits */
    318. 

  318. /* the unpardonable sin of including stdlib.h.  Hence, malloc(), free(), and */
    319. 

  319. /* exit() may or may not already be defined at this point.  I declare these  */
    320. 

  320. /* functions explicitly because some non-ANSI C compilers lack stdlib.h.     */
    321. 

  321. 

  322. #ifndef _STDLIB_H_
    323. 

  323. extern void *malloc();
    324. 

  324. extern void free();
    325. 

  325. extern void exit();
    326. 

  326. extern double strtod();
    327. 

  327. extern long strtol();
    328. 

  328. #endif /* _STDLIB_H_ */
    329. 

  329. 

  330. /* A few forward declarations.                                               */
    331. 

  331. 

  332. void poolrestart();
    333. 

  333. #ifndef TRILIBRARY
    334. 

  334. char *readline();
    335. 

  335. char *findfield();
    336. 

  336. #endif /* not TRILIBRARY */
    337. 

  337. 

  338. /* Labels that signify whether a record consists primarily of pointers or of */
    339. 

  339. /*   floating-point words.  Used to make decisions about data alignment.     */
    340. 

  340. 

  341. enum wordtype {POINTER, FLOATINGPOINT};
    342. 

  342. 

  343. /* Labels that signify the result of point location.  The result of a        */
    344. 

  344. /*   search indicates that the point falls in the interior of a triangle, on */
    345. 

  345. /*   an edge, on a vertex, or outside the mesh.                              */
    346. 

  346. 

  347. enum locateresult {INTRIANGLE, ONEDGE, ONVERTEX, OUTSIDE};
    348. 

  348. 

  349. /* Labels that signify the result of site insertion.  The result indicates   */
    350. 

  350. /*   that the point was inserted with complete success, was inserted but     */
    351. 

  351. /*   encroaches on a segment, was not inserted because it lies on a segment, */
    352. 

  352. /*   or was not inserted because another point occupies the same location.   */
    353. 

  353. 

  354. enum insertsiteresult {SUCCESSFULPOINT, ENCROACHINGPOINT, VIOLATINGPOINT,
    355. 

  355. 					   DUPLICATEPOINT};
    356. 

  356. 

  357. /* Labels that signify the result of direction finding.  The result          */
    358. 

  358. /*   indicates that a segment connecting the two query points falls within   */
    359. 

  359. /*   the direction triangle, along the left edge of the direction triangle,  */
    360. 

  360. /*   or along the right edge of the direction triangle.                      */
    361. 

  361. 

  362. enum finddirectionresult {WITHIN, LEFTCOLLINEAR, RIGHTCOLLINEAR};
    363. 

  363. 

  364. /* Labels that signify the result of the circumcenter computation routine.   */
    365. 

  365. /*   The return value indicates which edge of the triangle is shortest.      */
    366. 

  366. 

  367. enum circumcenterresult {OPPOSITEORG, OPPOSITEDEST, OPPOSITEAPEX};
    368. 

  368. 

  369. /*****************************************************************************/
    370. 

  370. /*                                                                           */
    371. 

  371. /*  The basic mesh data structures                                           */
    372. 

  372. /*                                                                           */
    373. 

  373. /*  There are three:  points, triangles, and shell edges (abbreviated        */
    374. 

  374. /*  `shelle').  These three data structures, linked by pointers, comprise    */
    375. 

  375. /*  the mesh.  A point simply represents a point in space and its properties.*/
    376. 

  376. /*  A triangle is a triangle.  A shell edge is a special data structure used */
    377. 

  377. /*  to represent impenetrable segments in the mesh (including the outer      */
    378. 

  378. /*  boundary, boundaries of holes, and internal boundaries separating two    */
    379. 

  379. /*  triangulated regions).  Shell edges represent boundaries defined by the  */
    380. 

  380. /*  user that triangles may not lie across.                                  */
    381. 

  381. /*                                                                           */
    382. 

  382. /*  A triangle consists of a list of three vertices, a list of three         */
    383. 

  383. /*  adjoining triangles, a list of three adjoining shell edges (when shell   */
    384. 

  384. /*  edges are used), an arbitrary number of optional user-defined floating-  */
    385. 

  385. /*  point attributes, and an optional area constraint.  The latter is an     */
    386. 

  386. /*  upper bound on the permissible area of each triangle in a region, used   */
    387. 

  387. /*  for mesh refinement.                                                     */
    388. 

  388. /*                                                                           */
    389. 

  389. /*  For a triangle on a boundary of the mesh, some or all of the neighboring */
    390. 

  390. /*  triangles may not be present.  For a triangle in the interior of the     */
    391. 

  391. /*  mesh, often no neighboring shell edges are present.  Such absent         */
    392. 

  392. /*  triangles and shell edges are never represented by NULL pointers; they   */
    393. 

  393. /*  are represented by two special records:  `dummytri', the triangle that   */
    394. 

  394. /*  fills "outer space", and `dummysh', the omnipresent shell edge.          */
    395. 

  395. /*  `dummytri' and `dummysh' are used for several reasons; for instance,     */
    396. 

  396. /*  they can be dereferenced and their contents examined without causing the */
    397. 

  397. /*  memory protection exception that would occur if NULL were dereferenced.  */
    398. 

  398. /*                                                                           */
    399. 

  399. /*  However, it is important to understand that a triangle includes other    */
    400. 

  400. /*  information as well.  The pointers to adjoining vertices, triangles, and */
    401. 

  401. /*  shell edges are ordered in a way that indicates their geometric relation */
    402. 

  402. /*  to each other.  Furthermore, each of these pointers contains orientation */
    403. 

  403. /*  information.  Each pointer to an adjoining triangle indicates which face */
    404. 

  404. /*  of that triangle is contacted.  Similarly, each pointer to an adjoining  */
    405. 

  405. /*  shell edge indicates which side of that shell edge is contacted, and how */
    406. 

  406. /*  the shell edge is oriented relative to the triangle.                     */
    407. 

  407. /*                                                                           */
    408. 

  408. /*  Shell edges are found abutting edges of triangles; either sandwiched     */
    409. 

  409. /*  between two triangles, or resting against one triangle on an exterior    */
    410. 

  410. /*  boundary or hole boundary.                                               */
    411. 

  411. /*                                                                           */
    412. 

  412. /*  A shell edge consists of a list of two vertices, a list of two           */
    413. 

  413. /*  adjoining shell edges, and a list of two adjoining triangles.  One of    */
    414. 

  414. /*  the two adjoining triangles may not be present (though there should      */
    415. 

  415. /*  always be one), and neighboring shell edges might not be present.        */
    416. 

  416. /*  Shell edges also store a user-defined integer "boundary marker".         */
    417. 

  417. /*  Typically, this integer is used to indicate what sort of boundary        */
    418. 

  418. /*  conditions are to be applied at that location in a finite element        */
    419. 

  419. /*  simulation.                                                              */
    420. 

  420. /*                                                                           */
    421. 

  421. /*  Like triangles, shell edges maintain information about the relative      */
    422. 

  422. /*  orientation of neighboring objects.                                      */
    423. 

  423. /*                                                                           */
    424. 

  424. /*  Points are relatively simple.  A point is a list of floating point       */
    425. 

  425. /*  numbers, starting with the x, and y coordinates, followed by an          */
    426. 

  426. /*  arbitrary number of optional user-defined floating-point attributes,     */
    427. 

  427. /*  followed by an integer boundary marker.  During the segment insertion    */
    428. 

  428. /*  phase, there is also a pointer from each point to a triangle that may    */
    429. 

  429. /*  contain it.  Each pointer is not always correct, but when one is, it     */
    430. 

  430. /*  speeds up segment insertion.  These pointers are assigned values once    */
    431. 

  431. /*  at the beginning of the segment insertion phase, and are not used or     */
    432. 

  432. /*  updated at any other time.  Edge swapping during segment insertion will  */
    433. 

  433. /*  render some of them incorrect.  Hence, don't rely upon them for          */
    434. 

  434. /*  anything.  For the most part, points do not have any information about   */
    435. 

  435. /*  what triangles or shell edges they are linked to.                        */
    436. 

  436. /*                                                                           */
    437. 

  437. /*****************************************************************************/
    438. 

  438. 

  439. /*****************************************************************************/
    440. 

  440. /*                                                                           */
    441. 

  441. /*  Handles                                                                  */
    442. 

  442. /*                                                                           */
    443. 

  443. /*  The oriented triangle (`triedge') and oriented shell edge (`edge') data  */
    444. 

  444. /*  structures defined below do not themselves store any part of the mesh.   */
    445. 

  445. /*  The mesh itself is made of `triangle's, `shelle's, and `point's.         */
    446. 

  446. /*                                                                           */
    447. 

  447. /*  Oriented triangles and oriented shell edges will usually be referred to  */
    448. 

  448. /*  as "handles".  A handle is essentially a pointer into the mesh; it       */
    449. 

  449. /*  allows you to "hold" one particular part of the mesh.  Handles are used  */
    450. 

  450. /*  to specify the regions in which one is traversing and modifying the mesh.*/
    451. 

  451. /*  A single `triangle' may be held by many handles, or none at all.  (The   */
    452. 

  452. /*  latter case is not a memory leak, because the triangle is still          */
    453. 

  453. /*  connected to other triangles in the mesh.)                               */
    454. 

  454. /*                                                                           */
    455. 

  455. /*  A `triedge' is a handle that holds a triangle.  It holds a specific side */
    456. 

  456. /*  of the triangle.  An `edge' is a handle that holds a shell edge.  It     */
    457. 

  457. /*  holds either the left or right side of the edge.                         */
    458. 

  458. /*                                                                           */
    459. 

  459. /*  Navigation about the mesh is accomplished through a set of mesh          */
    460. 

  460. /*  manipulation primitives, further below.  Many of these primitives take   */
    461. 

  461. /*  a handle and produce a new handle that holds the mesh near the first     */
    462. 

  462. /*  handle.  Other primitives take two handles and glue the corresponding    */
    463. 

  463. /*  parts of the mesh together.  The exact position of the handles is        */
    464. 

  464. /*  important.  For instance, when two triangles are glued together by the   */
    465. 

  465. /*  bond() primitive, they are glued by the sides on which the handles lie.  */
    466. 

  466. /*                                                                           */
    467. 

  467. /*  Because points have no information about which triangles they are        */
    468. 

  468. /*  attached to, I commonly represent a point by use of a handle whose       */
    469. 

  469. /*  origin is the point.  A single handle can simultaneously represent a     */
    470. 

  470. /*  triangle, an edge, and a point.                                          */
    471. 

  471. /*                                                                           */
    472. 

  472. /*****************************************************************************/
    473. 

  473. 

  474. /* The triangle data structure.  Each triangle contains three pointers to    */
    475. 

  475. /*   adjoining triangles, plus three pointers to vertex points, plus three   */
    476. 

  476. /*   pointers to shell edges (defined below; these pointers are usually      */
    477. 

  477. /*   `dummysh').  It may or may not also contain user-defined attributes     */
    478. 

  478. /*   and/or a floating-point "area constraint".  It may also contain extra   */
    479. 

  479. /*   pointers for nodes, when the user asks for high-order elements.         */
    480. 

  480. /*   Because the size and structure of a `triangle' is not decided until     */
    481. 

  481. /*   runtime, I haven't simply defined the type `triangle' to be a struct.   */
    482. 

  482. 

  483. typedef REAL **triangle;            /* Really:  typedef triangle *triangle   */
    484. 

  484. 

  485. /* An oriented triangle:  includes a pointer to a triangle and orientation.  */
    486. 

  486. /*   The orientation denotes an edge of the triangle.  Hence, there are      */
    487. 

  487. /*   three possible orientations.  By convention, each edge is always        */
    488. 

  488. /*   directed to point counterclockwise about the corresponding triangle.    */
    489. 

  489. 

  490. struct triedge {
    491. 

  491. 	triangle *tri;
    492. 

  492. 	int orient;                                       /* Ranges from 0 to 2. */
    493. 

  493. };
    494. 

  494. 

  495. /* The shell data structure.  Each shell edge contains two pointers to       */
    496. 

  496. /*   adjoining shell edges, plus two pointers to vertex points, plus two     */
    497. 

  497. /*   pointers to adjoining triangles, plus one shell marker.                 */
    498. 

  498. 

  499. typedef REAL **shelle;                  /* Really:  typedef shelle *shelle   */
    500. 

  500. 

  501. /* An oriented shell edge:  includes a pointer to a shell edge and an        */
    502. 

  502. /*   orientation.  The orientation denotes a side of the edge.  Hence, there */
    503. 

  503. /*   are two possible orientations.  By convention, the edge is always       */
    504. 

  504. /*   directed so that the "side" denoted is the right side of the edge.      */
    505. 

  505. 

  506. struct edge {
    507. 

  507. 	shelle *sh;
    508. 

  508. 	int shorient;                                     /* Ranges from 0 to 1. */
    509. 

  509. };
    510. 

  510. 

  511. /* The point data structure.  Each point is actually an array of REALs.      */
    512. 

  512. /*   The number of REALs is unknown until runtime.  An integer boundary      */
    513. 

  513. /*   marker, and sometimes a pointer to a triangle, is appended after the    */
    514. 

  514. /*   REALs.                                                                  */
    515. 

  515. 

  516. typedef REAL *point;
    517. 

  517. 

  518. /* A queue used to store encroached segments.  Each segment's vertices are   */
    519. 

  519. /*   stored so that one can check whether a segment is still the same.       */
    520. 

  520. 

  521. struct badsegment {
    522. 

  522. 	struct edge encsegment;                        /* An encroached segment. */
    523. 

  523. 	point segorg, segdest;                              /* The two vertices. */
    524. 

  524. 	struct badsegment *nextsegment;   /* Pointer to next encroached segment. */
    525. 

  525. };
    526. 

  526. 

  527. /* A queue used to store bad triangles.  The key is the square of the cosine */
    528. 

  528. /*   of the smallest angle of the triangle.  Each triangle's vertices are    */
    529. 

  529. /*   stored so that one can check whether a triangle is still the same.      */
    530. 

  530. 

  531. struct badface {
    532. 

  532. 	struct triedge badfacetri;                            /* A bad triangle. */
    533. 

  533. 	REAL key;                           /* cos^2 of smallest (apical) angle. */
    534. 

  534. 	point faceorg, facedest, faceapex;                /* The three vertices. */
    535. 

  535. 	struct badface *nextface;               /* Pointer to next bad triangle. */
    536. 

  536. };
    537. 

  537. 

  538. /* A node in a heap used to store events for the sweepline Delaunay          */
    539. 

  539. /*   algorithm.  Nodes do not point directly to their parents or children in */
    540. 

  540. /*   the heap.  Instead, each node knows its position in the heap, and can   */
    541. 

  541. /*   look up its parent and children in a separate array.  The `eventptr'    */
    542. 

  542. /*   points either to a `point' or to a triangle (in encoded format, so that */
    543. 

  543. /*   an orientation is included).  In the latter case, the origin of the     */
    544. 

  544. /*   oriented triangle is the apex of a "circle event" of the sweepline      */
    545. 

  545. /*   algorithm.  To distinguish site events from circle events, all circle   */
    546. 

  546. /*   events are given an invalid (smaller than `xmin') x-coordinate `xkey'.  */
    547. 

  547. 

  548. struct event {
    549. 

  549. 	REAL xkey, ykey;                            /* Coordinates of the event. */
    550. 

  550. 	VOID *eventptr;     /* Can be a point or the location of a circle event. */
    551. 

  551. 	int heapposition;            /* Marks this event's position in the heap. */
    552. 

  552. };
    553. 

  553. 

  554. /* A node in the splay tree.  Each node holds an oriented ghost triangle     */
    555. 

  555. /*   that represents a boundary edge of the growing triangulation.  When a   */
    556. 

  556. /*   circle event covers two boundary edges with a triangle, so that they    */
    557. 

  557. /*   are no longer boundary edges, those edges are not immediately deleted   */
    558. 

  558. /*   from the tree; rather, they are lazily deleted when they are next       */
    559. 

  559. /*   encountered.  (Since only a random sample of boundary edges are kept    */
    560. 

  560. /*   in the tree, lazy deletion is faster.)  `keydest' is used to verify     */
    561. 

  561. /*   that a triangle is still the same as when it entered the splay tree; if */
    562. 

  562. /*   it has been rotated (due to a circle event), it no longer represents a  */
    563. 

  563. /*   boundary edge and should be deleted.                                    */
    564. 

  564. 

  565. struct splaynode {
    566. 

  566. 	struct triedge keyedge;                /* Lprev of an edge on the front. */
    567. 

  567. 	point keydest;          /* Used to verify that splay node is still live. */
    568. 

  568. 	struct splaynode *lchild, *rchild;            /* Children in splay tree. */
    569. 

  569. };
    570. 

  570. 

  571. /* A type used to allocate memory.  firstblock is the first block of items.  */
    572. 

  572. /*   nowblock is the block from which items are currently being allocated.   */
    573. 

  573. /*   nextitem points to the next slab of free memory for an item.            */
    574. 

  574. /*   deaditemstack is the head of a linked list (stack) of deallocated items */
    575. 

  575. /*   that can be recycled.  unallocateditems is the number of items that     */
    576. 

  576. /*   remain to be allocated from nowblock.                                   */
    577. 

  577. /*                                                                           */
    578. 

  578. /* Traversal is the process of walking through the entire list of items, and */
    579. 

  579. /*   is separate from allocation.  Note that a traversal will visit items on */
    580. 

  580. /*   the "deaditemstack" stack as well as live items.  pathblock points to   */
    581. 

  581. /*   the block currently being traversed.  pathitem points to the next item  */
    582. 

  582. /*   to be traversed.  pathitemsleft is the number of items that remain to   */
    583. 

  583. /*   be traversed in pathblock.                                              */
    584. 

  584. /*                                                                           */
    585. 

  585. /* itemwordtype is set to POINTER or FLOATINGPOINT, and is used to suggest   */
    586. 

  586. /*   what sort of word the record is primarily made up of.  alignbytes       */
    587. 

  587. /*   determines how new records should be aligned in memory.  itembytes and  */
    588. 

  588. /*   itemwords are the length of a record in bytes (after rounding up) and   */
    589. 

  589. /*   words.  itemsperblock is the number of items allocated at once in a     */
    590. 

  590. /*   single block.  items is the number of currently allocated items.        */
    591. 

  591. /*   maxitems is the maximum number of items that have been allocated at     */
    592. 

  592. /*   once; it is the current number of items plus the number of records kept */
    593. 

  593. /*   on deaditemstack.                                                       */
    594. 

  594. 

  595. struct memorypool {
    596. 

  596. 	VOID **firstblock, **nowblock;
    597. 

  597. 	VOID *nextitem;
    598. 

  598. 	VOID *deaditemstack;
    599. 

  599. 	VOID **pathblock;
    600. 

  600. 	VOID *pathitem;
    601. 

  601. 	enum wordtype itemwordtype;
    602. 

  602. 	int alignbytes;
    603. 

  603. 	int itembytes, itemwords;
    604. 

  604. 	int itemsperblock;
    605. 

  605. 	long items, maxitems;
    606. 

  606. 	int unallocateditems;
    607. 

  607. 	int pathitemsleft;
    608. 

  608. };
    609. 

  609. 

  610. /* Variables used to allocate memory for triangles, shell edges, points,     */
    611. 

  611. /*   viri (triangles being eaten), bad (encroached) segments, bad (skinny    */
    612. 

  612. /*   or too large) triangles, and splay tree nodes.                          */
    613. 

  613. 

  614. static struct memorypool triangles;
    615. 

  615. static struct memorypool shelles;
    616. 

  616. static struct memorypool points;
    617. 

  617. static struct memorypool viri;
    618. 

  618. static struct memorypool badsegments;
    619. 

  619. static struct memorypool badtriangles;
    620. 

  620. static struct memorypool splaynodes;
    621. 

  621. 

  622. /* Variables that maintain the bad triangle queues.  The tails are pointers  */
    623. 

  623. /*   to the pointers that have to be filled in to enqueue an item.           */
    624. 

  624. 

  625. static struct badface *queuefront[64];
    626. 

  626. static struct badface **queuetail[64];
    627. 

  627. 

  628. static REAL xmin, xmax, ymin, ymax;                              /* x and y bounds. */
    629. 

  629. static REAL xminextreme;        /* Nonexistent x value used as a flag in sweepline. */
    630. 

  630. static int inpoints;                                     /* Number of input points. */
    631. 

  631. static int inelements;                                /* Number of input triangles. */
    632. 

  632. static int insegments;                                 /* Number of input segments. */
    633. 

  633. static int holes;                                         /* Number of input holes. */
    634. 

  634. static int regions;                                     /* Number of input regions. */
    635. 

  635. static long edges;                                       /* Number of output edges. */
    636. 

  636. static int mesh_dim;                                  /* Dimension (ought to be 2). */
    637. 

  637. static int nextras;                              /* Number of attributes per point. */
    638. 

  638. static int eextras;                           /* Number of attributes per triangle. */
    639. 

  639. static long hullsize;                            /* Number of edges of convex hull. */
    640. 

  640. static int triwords;                                   /* Total words per triangle. */
    641. 

  641. static int shwords;                                  /* Total words per shell edge. */
    642. 

  642. static int pointmarkindex;             /* Index to find boundary marker of a point. */
    643. 

  643. static int point2triindex;         /* Index to find a triangle adjacent to a point. */
    644. 

  644. static int highorderindex;    /* Index to find extra nodes for high-order elements. */
    645. 

  645. static int elemattribindex;              /* Index to find attributes of a triangle. */
    646. 

  646. static int areaboundindex;               /* Index to find area bound of a triangle. */
    647. 

  647. static int checksegments;           /* Are there segments in the triangulation yet? */
    648. 

  648. static int readnodefile;                             /* Has a .node file been read? */
    649. 

  649. static long samples;                /* Number of random samples for point location. */
    650. 

  650. static unsigned long randomseed;                     /* Current random number seed. */
    651. 

  651. 

  652. static REAL splitter;       /* Used to split REAL factors for exact multiplication. */
    653. 

  653. static REAL epsilon;                             /* Floating-point machine epsilon. */
    654. 

  654. static REAL resulterrbound;
    655. 

  655. static REAL ccwerrboundA, ccwerrboundB, ccwerrboundC;
    656. 

  656. static REAL iccerrboundA, iccerrboundB, iccerrboundC;
    657. 

  657. 

  658. static long incirclecount;                   /* Number of incircle tests performed. */
    659. 

  659. static long counterclockcount;       /* Number of counterclockwise tests performed. */
    660. 

  660. static long hyperbolacount;        /* Number of right-of-hyperbola tests performed. */
    661. 

  661. static long circumcentercount;    /* Number of circumcenter calculations performed. */
    662. 

  662. static long circletopcount;         /* Number of circle top calculations performed. */
    663. 

  663. 

  664. /* Switches for the triangulator.                                            */
    665. 

  665. /*   poly: -p switch.  refine: -r switch.                                    */
    666. 

  666. /*   quality: -q switch.                                                     */
    667. 

  667. /*     minangle: minimum angle bound, specified after -q switch.             */
    668. 

  668. /*     goodangle: cosine squared of minangle.                                */
    669. 

  669. /*   vararea: -a switch without number.                                      */
    670. 

  670. /*   fixedarea: -a switch with number.                                       */
    671. 

  671. /*     maxarea: maximum area bound, specified after -a switch.               */
    672. 

  672. /*   regionattrib: -A switch.  convex: -c switch.                            */
    673. 

  673. /*   firstnumber: inverse of -z switch.  All items are numbered starting     */
    674. 

  674. /*     from firstnumber.                                                     */
    675. 

  675. /*   edgesout: -e switch.  voronoi: -v switch.                               */
    676. 

  676. /*   neighbors: -n switch.  geomview: -g switch.                             */
    677. 

  677. /*   nobound: -B switch.  nopolywritten: -P switch.                          */
    678. 

  678. /*   nonodewritten: -N switch.  noelewritten: -E switch.                     */
    679. 

  679. /*   noiterationnum: -I switch.  noholes: -O switch.                         */
    680. 

  680. /*   noexact: -X switch.                                                     */
    681. 

  681. /*   order: element order, specified after -o switch.                        */
    682. 

  682. /*   nobisect: count of how often -Y switch is selected.                     */
    683. 

  683. /*   steiner: maximum number of Steiner points, specified after -S switch.   */
    684. 

  684. /*     steinerleft: number of Steiner points not yet used.                   */
    685. 

  685. /*   incremental: -i switch.  sweepline: -F switch.                          */
    686. 

  686. /*   dwyer: inverse of -l switch.                                            */
    687. 

  687. /*   splitseg: -s switch.                                                    */
    688. 

  688. /*   docheck: -C switch.                                                     */
    689. 

  689. /*   quiet: -Q switch.  verbose: count of how often -V switch is selected.   */
    690. 

  690. /*   useshelles: -p, -r, -q, or -c switch; determines whether shell edges    */
    691. 

  691. /*     are used at all.                                                      */
    692. 

  692. /*                                                                           */
    693. 

  693. /* Read the instructions to find out the meaning of these switches.          */
    694. 

  694. 

  695. static int poly, refine, quality, vararea, fixedarea, regionattrib, convex;
    696. 

  696. static int firstnumber;
    697. 

  697. static int edgesout, voronoi, neighbors, geomview;
    698. 

  698. static int nobound, nopolywritten, nonodewritten, noelewritten, noiterationnum;
    699. 

  699. static int noholes, noexact;
    700. 

  700. static int incremental, sweepline, dwyer;
    701. 

  701. static int splitseg;
    702. 

  702. static int docheck;
    703. 

  703. static int quiet, verbose;
    704. 

  704. static int useshelles;
    705. 

  705. static int order;
    706. 

  706. static int nobisect;
    707. 

  707. static int steiner, steinerleft;
    708. 

  708. static REAL minangle, goodangle;
    709. 

  709. static REAL maxarea;
    710. 

  710. 

  711. /* Variables for file names.                                                 */
    712. 

  712. 

  713. #ifndef TRILIBRARY
    714. 

  714. char innodefilename[FILENAMESIZE];
    715. 

  715. char inelefilename[FILENAMESIZE];
    716. 

  716. char inpolyfilename[FILENAMESIZE];
    717. 

  717. char areafilename[FILENAMESIZE];
    718. 

  718. char outnodefilename[FILENAMESIZE];
    719. 

  719. char outelefilename[FILENAMESIZE];
    720. 

  720. char outpolyfilename[FILENAMESIZE];
    721. 

  721. char edgefilename[FILENAMESIZE];
    722. 

  722. char vnodefilename[FILENAMESIZE];
    723. 

  723. char vedgefilename[FILENAMESIZE];
    724. 

  724. char neighborfilename[FILENAMESIZE];
    725. 

  725. char offfilename[FILENAMESIZE];
    726. 

  726. #endif /* not TRILIBRARY */
    727. 

  727. 

  728. /* Triangular bounding box points.                                           */
    729. 

  729. 

  730. static point infpoint1, infpoint2, infpoint3;
    731. 

  731. 

  732. /* Pointer to the `triangle' that occupies all of "outer space".             */
    733. 

  733. 

  734. static triangle *dummytri;
    735. 

  735. static triangle *dummytribase;      /* Keep base address so we can free() it later. */
    736. 

  736. 

  737. /* Pointer to the omnipresent shell edge.  Referenced by any triangle or     */
    738. 

  738. /*   shell edge that isn't really connected to a shell edge at that          */
    739. 

  739. /*   location.                                                               */
    740. 

  740. 

  741. static shelle *dummysh;
    742. 

  742. static shelle *dummyshbase;         /* Keep base address so we can free() it later. */
    743. 

  743. 

  744. /* Pointer to a recently visited triangle.  Improves point location if       */
    745. 

  745. /*   proximate points are inserted sequentially.                             */
    746. 

  746. 

  747. static struct triedge recenttri;
    748. 

  748. 

  749. /*****************************************************************************/
    750. 

  750. /*                                                                           */
    751. 

  751. /*  Mesh manipulation primitives.  Each triangle contains three pointers to  */
    752. 

  752. /*  other triangles, with orientations.  Each pointer points not to the      */
    753. 

  753. /*  first byte of a triangle, but to one of the first three bytes of a       */
    754. 

  754. /*  triangle.  It is necessary to extract both the triangle itself and the   */
    755. 

  755. /*  orientation.  To save memory, I keep both pieces of information in one   */
    756. 

  756. /*  pointer.  To make this possible, I assume that all triangles are aligned */
    757. 

  757. /*  to four-byte boundaries.  The `decode' routine below decodes a pointer,  */
    758. 

  758. /*  extracting an orientation (in the range 0 to 2) and a pointer to the     */
    759. 

  759. /*  beginning of a triangle.  The `encode' routine compresses a pointer to a */
    760. 

  760. /*  triangle and an orientation into a single pointer.  My assumptions that  */
    761. 

  761. /*  triangles are four-byte-aligned and that the `unsigned long' type is     */
    762. 

  762. /*  long enough to hold a pointer are two of the few kludges in this program.*/
    763. 

  763. /*                                                                           */
    764. 

  764. /*  Shell edges are manipulated similarly.  A pointer to a shell edge        */
    765. 

  765. /*  carries both an address and an orientation in the range 0 to 1.          */
    766. 

  766. /*                                                                           */
    767. 

  767. /*  The other primitives take an oriented triangle or oriented shell edge,   */
    768. 

  768. /*  and return an oriented triangle or oriented shell edge or point; or they */
    769. 

  769. /*  change the connections in the data structure.                            */
    770. 

  770. /*                                                                           */
    771. 

  771. /*****************************************************************************/
    772. 

  772. 

  773. /********* Mesh manipulation primitives begin here                   *********/
    774. 

  774. /**                                                                         **/
    775. 

  775. /**                                                                         **/
    776. 

  776. 

  777. /* Fast lookup arrays to speed some of the mesh manipulation primitives.     */
    778. 

  778. 

  779. int plus1mod3[3] = {1, 2, 0};
    780. 

  780. int minus1mod3[3] = {2, 0, 1};
    781. 

  781. 

  782. /********* Primitives for triangles                                  *********/
    783. 

  783. /*                                                                           */
    784. 

  784. /*                                                                           */
    785. 

  785. 

  786. /* decode() converts a pointer to an oriented triangle.  The orientation is  */
    787. 

  787. /*   extracted from the two least significant bits of the pointer.           */
    788. 

  788. 

  789. #define decode( ptr, triedge )													\
    790. 

  790. 	( triedge ).orient = (int) ( (unsigned long) ( ptr ) & (unsigned long) 3l );	  \
    791. 

  791. 	( triedge ).tri = (triangle *)												  \
    792. 

  792. 					  ( (unsigned long) ( ptr ) ^ (unsigned long) ( triedge ).orient )
    793. 

  793. 

  794. /* encode() compresses an oriented triangle into a single pointer.  It       */
    795. 

  795. /*   relies on the assumption that all triangles are aligned to four-byte    */
    796. 

  796. /*   boundaries, so the two least significant bits of (triedge).tri are zero.*/
    797. 

  797. 

  798. #define encode( triedge )														\
    799. 

  799. 	(triangle) ( (unsigned long) ( triedge ).tri | (unsigned long) ( triedge ).orient )
    800. 

  800. 

  801. /* The following edge manipulation primitives are all described by Guibas    */
    802. 

  802. /*   and Stolfi.  However, they use an edge-based data structure, whereas I  */
    803. 

  803. /*   am using a triangle-based data structure.                               */
    804. 

  804. 

  805. /* sym() finds the abutting triangle, on the same edge.  Note that the       */
    806. 

  806. /*   edge direction is necessarily reversed…
    807.

Large files files are truncated, but you can click here to view the full file