-
Notifications
You must be signed in to change notification settings - Fork 375
/
Copy pathtriangle.c
16027 lines (14803 loc) · 639 KB
/
triangle.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
/* Preamble:
* (1) This is a slightly modified version of triangle.c that allows it to be
* compiled on both Windows and Linux/OSX/Unix with 32- or 64-bits. The
* problem is the extensive use of "long" which is 4 bytes on Windows and
* 8 bytes on the others. The effects of this problem started to manifest
* (crashes) only on binaries compiled with VS2012 and run on Windows 8.
* Although not tested it is likely that MinGW64 binaries are affected too.
*
* GMT Team, June 2013.
*
* (2) While containing lots of local functions, none are declared static.
* This means they are all exported in the shared library. We have thus
* added the "static" keyword before all the function declarations other
* than triangulate, which is the only one we are using from GMT.
* Likewise we extended the #ifndef CDT_ONLY and #ifndef REDUCED pragmas
* to the internal functions that are used only when those switches are
* actually defined.
*
* GMT Team, March 2016.
*
* Declaration of type LONG:
*/
#if defined(_WIN64) || defined(__MINGW64__)
#define LONG __int64 /* A signed 8-byte integer, under Windows64 */
#define PRIu "I64u" /* printf format for unsigned uint64 */
#else
#define LONG long /* A signed 4 or 8-byte integer */
#define PRIu "lu"
#endif
/*****************************************************************************/
/* */
/* 888888888 ,o, / 888 */
/* 888 88o88o " o8888o 88o8888o o88888o 888 o88888o */
/* 888 888 888 88b 888 888 888 888 888 d888 88b */
/* 888 888 888 o88^o888 888 888 "88888" 888 8888oo888 */
/* 888 888 888 C888 888 888 888 / 888 q888 */
/* 888 888 888 "88o^888 888 888 Cb 888 "88oooo" */
/* "8oo8D */
/* */
/* A Two-Dimensional Quality Mesh Generator and Delaunay Triangulator. */
/* (triangle.c) */
/* */
/* Version 1.6 */
/* July 28, 2005 */
/* */
/* Copyright 1993, 1995, 1997, 1998, 2002, 2005 */
/* Jonathan Richard Shewchuk */
/* 2360 Woolsey #H */
/* Berkeley, California 94705-1927 */
/* [email protected] */
/* */
/* 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 */
/* */
/* Disclaimer: Neither I nor Carnegie Mellon warrant this code in any way */
/* whatsoever. This code is provided "as-is". Use at your own risk. */
/* */
/* Some of the references listed below are marked with an asterisk. [*] */
/* These references are available for downloading from the Web page */
/* */
/* http://www.cs.cmu.edu/~quake/triangle.research.html */
/* */
/* Three papers discussing aspects of Triangle are available. A short */
/* overview appears in "Triangle: Engineering a 2D Quality Mesh */
/* Generator and Delaunay Triangulator," in Applied Computational */
/* Geometry: Towards Geometric Engineering, Ming C. Lin and Dinesh */
/* Manocha, editors, Lecture Notes in Computer Science volume 1148, */
/* pages 203-222, Springer-Verlag, Berlin, May 1996 (from the First ACM */
/* Workshop on Applied Computational Geometry). [*] */
/* */
/* The algorithms are discussed in the greatest detail in "Delaunay */
/* Refinement Algorithms for Triangular Mesh Generation," Computational */
/* Geometry: Theory and Applications 22(1-3):21-74, May 2002. [*] */
/* */
/* More detail about the data structures may be found in my dissertation: */
/* "Delaunay Refinement Mesh Generation," Ph.D. thesis, Technical Report */
/* CMU-CS-97-137, School of Computer Science, Carnegie Mellon University, */
/* Pittsburgh, Pennsylvania, 18 May 1997. [*] */
/* */
/* Triangle was created as part of the Quake Project in the School of */
/* Computer Science at Carnegie Mellon University. For further */
/* information, see Hesheng Bao, Jacobo Bielak, Omar Ghattas, Loukas F. */
/* Kallivokas, David R. O'Hallaron, Jonathan R. Shewchuk, and Jifeng Xu, */
/* "Large-scale Simulation of Elastic Wave Propagation in Heterogeneous */
/* Media on Parallel Computers," Computer Methods in Applied Mechanics */
/* and Engineering 152(1-2):85-102, 22 January 1998. */
/* */
/* Triangle's Delaunay refinement algorithm for quality mesh generation is */
/* a hybrid of one due to Jim Ruppert, "A Delaunay Refinement Algorithm */
/* for Quality 2-Dimensional Mesh Generation," Journal of Algorithms */
/* 18(3):548-585, May 1995 [*], and one due to L. Paul Chew, "Guaranteed- */
/* Quality Mesh Generation for Curved Surfaces," Proceedings of the Ninth */
/* Annual Symposium on Computational Geometry (San Diego, California), */
/* pages 274-280, Association for Computing Machinery, May 1993, */
/* http://portal.acm.org/citation.cfm?id=161150 . */
/* */
/* The Delaunay refinement algorithm has been modified so that it meshes */
/* domains with small input angles well, as described in Gary L. Miller, */
/* Steven E. Pav, and Noel J. Walkington, "When and Why Ruppert's */
/* Algorithm Works," Twelfth International Meshing Roundtable, pages */
/* 91-102, Sandia National Laboratories, September 2003. [*] */
/* */
/* 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 the switch 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, http://portal.acm.org/citation.cfm?id=282923 .*/
/* */
/* 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. Triangle's improvement of 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 vertex */
/* insertion (I currently don't bother), their result combined with the */
/* result of Kenneth L. Clarkson and Peter W. Shor, "Applications of */
/* Random Sampling in Computational Geometry II," Discrete & */
/* Computational Geometry 4(1):387-421, 1989, 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, */
/* http://portal.acm.org/citation.cfm?id=3835 . */
/* */
/* 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," Discrete & */
/* Computational Geometry 18(3):305-363, October 1997. (Also available */
/* as Technical Report CMU-CS-96-140, School of Computer Science, */
/* Carnegie Mellon University, Pittsburgh, Pennsylvania, May 1996.) [*] */
/* An abbreviated version appears as Jonathan Richard Shewchuk, "Robust */
/* Adaptive Floating-Point Geometric Predicates," Proceedings of the */
/* Twelfth Annual Symposium on Computational 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, pp. 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 Computa- */
/* tional 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 Triangu- */
/* lations," International Journal of Computational Geometry & Applica- */
/* tions 5(1-2):193-213, March-June 1995. */
/* */
/* The method of inserting new vertices off-center (not precisely at the */
/* circumcenter of every poor-quality triangle) is from Alper Ungor, */
/* "Off-centers: A New Type of Steiner Points for Computing Size-Optimal */
/* Quality-Guaranteed Delaunay Triangulations," Proceedings of LATIN */
/* 2004 (Buenos Aires, Argentina), April 2004. */
/* */
/* 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(t^2) per segment in the worst case and O(t) per segment in the */
/* common case, where t 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. I could improve this to O(d log d) */
/* time, but d is usually quite small, so it's not worth the bother. */
/* (This note does not apply when the -s switch is used, invoking a */
/* different method is used to insert segments.) */
/* */
/* The time for deleting a vertex from a Delaunay triangulation is O(d^2) */
/* in the worst case and O(d) in the common case, where d is the degree */
/* of the vertex being deleted. I could improve this to O(d log d) time, */
/* but d is usually quite small, so it's not worth the bother. */
/* */
/* 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 */
/* quadratic time is required, but these never arise in practice. */
/* */
/* The geometric predicates (circumcenter calculations, segment */
/* intersection formulae, etc.) appear in my "Lecture Notes on Geometric */
/* Robustness" at http://www.cs.berkeley.edu/~jrs/mesh . */
/* */
/* 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. */
/* */
/*****************************************************************************/
/* 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 until your work is debugged. */
/* #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, -u, -D, -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, my exact arithmetic routines might be defeated by the */
/* use of internal extended precision floating-point registers. The best */
/* way to solve this problem is to set the floating-point registers to use */
/* single or double precision internally. On 80x86 processors, this may */
/* be accomplished by setting the CPU86 symbol for the Microsoft C */
/* compiler, or the LINUX symbol for the gcc compiler running on Linux. */
/* */
/* An inferior solution is to declare certain values as `volatile', thus */
/* forcing them to be stored to memory and rounded off. Unfortunately, */
/* this solution might slow Triangle down quite a bit. To use volatile */
/* values, write "#define INEXACT volatile" below. Normally, however, */
/* INEXACT should be defined to be nothing. ("#define INEXACT".) */
/* */
/* For more discussion, see http://www.cs.cmu.edu/~quake/robust.pc.html . */
/* For yet more discussion, see Section 5 of my paper, "Adaptive Precision */
/* Floating-Point Arithmetic and Fast Robust Geometric Predicates" (also */
/* available as Section 6.6 of my dissertation). */
/* #define CPU86 */
/* #define LINUX */
#define INEXACT /* Nothing */
/* #define INEXACT volatile */
/* Maximum number of characters in a file name (including the null). */
#define FILENAMESIZE 2048
/* Maximum number of characters in a line read from a file (including the */
/* null). */
#define INPUTLINESIZE 1024
/* For efficiency, a variety of data structures are allocated in bulk. The */
/* following constants determine how many of each structure is allocated */
/* at once. */
#define TRIPERBLOCK 4092 /* Number of triangles allocated at once. */
#define SUBSEGPERBLOCK 508 /* Number of subsegments allocated at once. */
#define VERTEXPERBLOCK 4092 /* Number of vertices allocated at once. */
#define VIRUSPERBLOCK 1020 /* Number of virus triangles allocated at once. */
/* Number of encroached subsegments allocated at once. */
#define BADSUBSEGPERBLOCK 252
/* Number of skinny triangles allocated at once. */
#define BADTRIPERBLOCK 4092
/* Number of flipped triangles allocated at once. */
#define FLIPSTACKERPERBLOCK 252
/* Number of splay tree nodes allocated at once. */
#define SPLAYNODEPERBLOCK 508
/* The vertex types. A DEADVERTEX has been deleted entirely. An */
/* UNDEADVERTEX is not part of the mesh, but is written to the output */
/* .node file and affects the node indexing in the other output files. */
#define INPUTVERTEX 0
#define SEGMENTVERTEX 1
#define FREEVERTEX 2
#define DEADVERTEX -32768
#define UNDEADVERTEX -32767
/* The next line is used to outsmart some very stupid compilers. If your */
/* compiler is smarter, feel free to replace the "int" with "void". */
/* Not that it matters. */
#define VOID void
/* Two constants for algorithms based on random sampling. Both constants */
/* have been chosen empirically to optimize their respective algorithms. */
/* Used for the point location scheme of Mucke, Saias, and Zhu, to decide */
/* how large a random sample of triangles to inspect. */
#define SAMPLEFACTOR 11
/* Used in Fortune's sweepline Delaunay algorithm to determine what fraction */
/* of boundary edges should be maintained in the splay tree for point */
/* location on the front. */
#define SAMPLERATE 10
/* A number that speaks for itself, every kissable digit. */
#define PI 3.141592653589793238462643383279502884197169399375105820974944592308
/* Another fave. */
#define SQUAREROOTTWO 1.4142135623730950488016887242096980785696718753769480732
/* And here's one for those of you who are intimidated by math. */
#define ONETHIRD 0.333333333333333333333333333333333333333333333333333333333333
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#ifndef NO_TIMER
#include <sys/time.h>
#endif /* not NO_TIMER */
#include <float.h>
#ifdef LINUX
#include <fpu_control.h>
#endif /* LINUX */
#ifdef TRILIBRARY
#include "triangle.h"
#endif /* TRILIBRARY */
/* A few forward declarations. */
#ifndef TRILIBRARY
char *readline();
char *findfield();
#endif /* not TRILIBRARY */
/* Labels that signify the result of point location. The result of a */
/* search indicates that the point falls in the interior of a triangle, on */
/* an edge, on a vertex, or outside the mesh. */
enum locateresult {INTRIANGLE, ONEDGE, ONVERTEX, OUTSIDE};
/* Labels that signify the result of vertex insertion. The result indicates */
/* that the vertex was inserted with complete success, was inserted but */
/* encroaches upon a subsegment, was not inserted because it lies on a */
/* segment, or was not inserted because another vertex occupies the same */
/* location. */
enum insertvertexresult {SUCCESSFULVERTEX, ENCROACHINGVERTEX, VIOLATINGVERTEX,
DUPLICATEVERTEX};
/* Labels that signify the result of direction finding. The result */
/* indicates that a segment connecting the two query points falls within */
/* the direction triangle, along the left edge of the direction triangle, */
/* or along the right edge of the direction triangle. */
enum finddirectionresult {WITHIN, LEFTCOLLINEAR, RIGHTCOLLINEAR};
/*****************************************************************************/
/* */
/* The basic mesh data structures */
/* */
/* There are three: vertices, triangles, and subsegments (abbreviated */
/* `subseg'). These three data structures, linked by pointers, comprise */
/* the mesh. A vertex simply represents a mesh vertex and its properties. */
/* A triangle is a triangle. A subsegment is a special data structure used */
/* to represent an impenetrable edge of the mesh (perhaps on the outer */
/* boundary, on the boundary of a hole, or part of an internal boundary */
/* separating two triangulated regions). Subsegments represent boundaries, */
/* defined by the user, that triangles may not lie across. */
/* */
/* A triangle consists of a list of three vertices, a list of three */
/* adjoining triangles, a list of three adjoining subsegments (when */
/* segments exist), an arbitrary number of optional user-defined */
/* floating-point attributes, and an optional area constraint. The latter */
/* is an upper bound on the permissible area of each triangle in a region, */
/* used for mesh refinement. */
/* */
/* For a triangle on a boundary of the mesh, some or all of the neighboring */
/* triangles may not be present. For a triangle in the interior of the */
/* mesh, often no neighboring subsegments are present. Such absent */
/* triangles and subsegments are never represented by NULL pointers; they */
/* are represented by two special records: `dummytri', the triangle that */
/* fills "outer space", and `dummysub', the omnipresent subsegment. */
/* `dummytri' and `dummysub' are used for several reasons; for instance, */
/* they can be dereferenced and their contents examined without violating */
/* protected memory. */
/* */
/* However, it is important to understand that a triangle includes other */
/* information as well. The pointers to adjoining vertices, triangles, and */
/* subsegments are ordered in a way that indicates their geometric relation */
/* to each other. Furthermore, each of these pointers contains orientation */
/* information. Each pointer to an adjoining triangle indicates which face */
/* of that triangle is contacted. Similarly, each pointer to an adjoining */
/* subsegment indicates which side of that subsegment is contacted, and how */
/* the subsegment is oriented relative to the triangle. */
/* */
/* The data structure representing a subsegment may be thought to be */
/* abutting the edge of one or two triangle data structures: either */
/* sandwiched between two triangles, or resting against one triangle on an */
/* exterior boundary or hole boundary. */
/* */
/* A subsegment consists of a list of four vertices--the vertices of the */
/* subsegment, and the vertices of the segment it is a part of--a list of */
/* two adjoining subsegments, and a list of two adjoining triangles. One */
/* of the two adjoining triangles may not be present (though there should */
/* always be one), and neighboring subsegments might not be present. */
/* Subsegments also store a user-defined integer "boundary marker". */
/* Typically, this integer is used to indicate what boundary conditions are */
/* to be applied at that location in a finite element simulation. */
/* */
/* Like triangles, subsegments maintain information about the relative */
/* orientation of neighboring objects. */
/* */
/* Vertices are relatively simple. A vertex is a list of floating-point */
/* numbers, starting with the x, and y coordinates, followed by an */
/* arbitrary number of optional user-defined floating-point attributes, */
/* followed by an integer boundary marker. During the segment insertion */
/* phase, there is also a pointer from each vertex to a triangle that may */
/* contain it. Each pointer is not always correct, but when one is, it */
/* speeds up segment insertion. These pointers are assigned values once */
/* at the beginning of the segment insertion phase, and are not used or */
/* updated except during this phase. Edge flipping during segment */
/* insertion will render some of them incorrect. Hence, don't rely upon */
/* them for anything. */
/* */
/* Other than the exception mentioned above, vertices have no information */
/* about what triangles, subfacets, or subsegments they are linked to. */
/* */
/*****************************************************************************/
/*****************************************************************************/
/* */
/* Handles */
/* */
/* The oriented triangle (`otri') and oriented subsegment (`osub') data */
/* structures defined below do not themselves store any part of the mesh. */
/* The mesh itself is made of `triangle's, `subseg's, and `vertex's. */
/* */
/* Oriented triangles and oriented subsegments will usually be referred to */
/* as "handles." A handle is essentially a pointer into the mesh; it */
/* allows you to "hold" one particular part of the mesh. Handles are used */
/* to specify the regions in which one is traversing and modifying the mesh.*/
/* A single `triangle' may be held by many handles, or none at all. (The */
/* latter case is not a memory leak, because the triangle is still */
/* connected to other triangles in the mesh.) */
/* */
/* An `otri' is a handle that holds a triangle. It holds a specific edge */
/* of the triangle. An `osub' is a handle that holds a subsegment. It */
/* holds either the left or right side of the subsegment. */
/* */
/* Navigation about the mesh is accomplished through a set of mesh */
/* manipulation primitives, further below. Many of these primitives take */
/* a handle and produce a new handle that holds the mesh near the first */
/* handle. Other primitives take two handles and glue the corresponding */
/* parts of the mesh together. The orientation of the handles is */
/* important. For instance, when two triangles are glued together by the */
/* bond() primitive, they are glued at the edges on which the handles lie. */
/* */
/* Because vertices have no information about which triangles they are */
/* attached to, I commonly represent a vertex by use of a handle whose */
/* origin is the vertex. A single handle can simultaneously represent a */
/* triangle, an edge, and a vertex. */
/* */
/*****************************************************************************/
/* The triangle data structure. Each triangle contains three pointers to */
/* adjoining triangles, plus three pointers to vertices, plus three */
/* pointers to subsegments (declared below; these pointers are usually */
/* `dummysub'). It may or may not also contain user-defined attributes */
/* and/or a floating-point "area constraint." It may also contain extra */
/* pointers for nodes, when the user asks for high-order elements. */
/* Because the size and structure of a `triangle' is not decided until */
/* runtime, I haven't simply declared the type `triangle' as a struct. */
typedef REAL **triangle; /* Really: typedef triangle *triangle */
/* An oriented triangle: includes a pointer to a triangle and orientation. */
/* The orientation denotes an edge of the triangle. Hence, there are */
/* three possible orientations. By convention, each edge always points */
/* counterclockwise about the corresponding triangle. */
struct otri {
triangle *tri;
int orient; /* Ranges from 0 to 2. */
};
/* The subsegment data structure. Each subsegment contains two pointers to */
/* adjoining subsegments, plus four pointers to vertices, plus two */
/* pointers to adjoining triangles, plus one boundary marker, plus one */
/* segment number. */
typedef REAL **subseg; /* Really: typedef subseg *subseg */
/* An oriented subsegment: includes a pointer to a subsegment and an */
/* orientation. The orientation denotes a side of the edge. Hence, there */
/* are two possible orientations. By convention, the edge is always */
/* directed so that the "side" denoted is the right side of the edge. */
struct osub {
subseg *ss;
int ssorient; /* Ranges from 0 to 1. */
};
/* The vertex data structure. Each vertex is actually an array of REALs. */
/* The number of REALs is unknown until runtime. An integer boundary */
/* marker, and sometimes a pointer to a triangle, is appended after the */
/* REALs. */
typedef REAL *vertex;
/* A queue used to store encroached subsegments. Each subsegment's vertices */
/* are stored so that we can check whether a subsegment is still the same. */
struct badsubseg {
subseg encsubseg; /* An encroached subsegment. */
vertex subsegorg, subsegdest; /* Its two vertices. */
};
/* A queue used to store bad triangles. The key is the square of the cosine */
/* of the smallest angle of the triangle. Each triangle's vertices are */
/* stored so that one can check whether a triangle is still the same. */
struct badtriang {
triangle poortri; /* A skinny or too-large triangle. */
REAL key; /* cos^2 of smallest (apical) angle. */
vertex triangorg, triangdest, triangapex; /* Its three vertices. */
struct badtriang *nexttriang; /* Pointer to next bad triangle. */
};
/* A stack of triangles flipped during the most recent vertex insertion. */
/* The stack is used to undo the vertex insertion if the vertex encroaches */
/* upon a subsegment. */
struct flipstacker {
triangle flippedtri; /* A recently flipped triangle. */
struct flipstacker *prevflip; /* Previous flip in the stack. */
};
/* A node in a heap used to store events for the sweepline Delaunay */
/* algorithm. Nodes do not point directly to their parents or children in */
/* the heap. Instead, each node knows its position in the heap, and can */
/* look up its parent and children in a separate array. The `eventptr' */
/* points either to a `vertex' or to a triangle (in encoded format, so */
/* that an orientation is included). In the latter case, the origin of */
/* the oriented triangle is the apex of a "circle event" of the sweepline */
/* algorithm. To distinguish site events from circle events, all circle */
/* events are given an invalid (smaller than `xmin') x-coordinate `xkey'. */
struct event {
REAL xkey, ykey; /* Coordinates of the event. */
VOID *eventptr; /* Can be a vertex or the location of a circle event. */
int heapposition; /* Marks this event's position in the heap. */
};
/* A node in the splay tree. Each node holds an oriented ghost triangle */
/* that represents a boundary edge of the growing triangulation. When a */
/* circle event covers two boundary edges with a triangle, so that they */
/* are no longer boundary edges, those edges are not immediately deleted */
/* from the tree; rather, they are lazily deleted when they are next */
/* encountered. (Since only a random sample of boundary edges are kept */
/* in the tree, lazy deletion is faster.) `keydest' is used to verify */
/* that a triangle is still the same as when it entered the splay tree; if */
/* it has been rotated (due to a circle event), it no longer represents a */
/* boundary edge and should be deleted. */
struct splaynode {
struct otri keyedge; /* Lprev of an edge on the front. */
vertex keydest; /* Used to verify that splay node is still live. */
struct splaynode *lchild, *rchild; /* Children in splay tree. */
};
/* A type used to allocate memory. firstblock is the first block of items. */
/* nowblock is the block from which items are currently being allocated. */
/* nextitem points to the next slab of free memory for an item. */
/* deaditemstack is the head of a linked list (stack) of deallocated items */
/* that can be recycled. unallocateditems is the number of items that */
/* remain to be allocated from nowblock. */
/* */
/* Traversal is the process of walking through the entire list of items, and */
/* is separate from allocation. Note that a traversal will visit items on */
/* the "deaditemstack" stack as well as live items. pathblock points to */
/* the block currently being traversed. pathitem points to the next item */
/* to be traversed. pathitemsleft is the number of items that remain to */
/* be traversed in pathblock. */
/* */
/* alignbytes determines how new records should be aligned in memory. */
/* itembytes is the length of a record in bytes (after rounding up). */
/* itemsperblock is the number of items allocated at once in a single */
/* block. itemsfirstblock is the number of items in the first block, */
/* which can vary from the others. items is the number of currently */
/* allocated items. maxitems is the maximum number of items that have */
/* been allocated at once; it is the current number of items plus the */
/* number of records kept on deaditemstack. */
struct memorypool {
VOID **firstblock, **nowblock;
VOID *nextitem;
VOID *deaditemstack;
VOID **pathblock;
VOID *pathitem;
size_t alignbytes;
size_t itembytes;
size_t itemsperblock;
size_t itemsfirstblock;
LONG items, maxitems;
size_t unallocateditems;
size_t pathitemsleft;
};
/* Global constants. */
static REAL splitter; /* Used to split REAL factors for exact multiplication. */
static REAL epsilon; /* Floating-point machine epsilon. */
static REAL resulterrbound;
static REAL ccwerrboundA, ccwerrboundB, ccwerrboundC;
static REAL iccerrboundA, iccerrboundB, iccerrboundC;
static REAL o3derrboundA, o3derrboundB, o3derrboundC;
/* Random number seed is not constant, but I've made it global anyway. */
static unsigned LONG randomseed; /* Current random number seed. */
/* Mesh data structure. Triangle operates on only one mesh, but the mesh */
/* structure is used (instead of global variables) to allow reentrancy. */
struct mesh {
/* Variables used to allocate memory for triangles, subsegments, vertices, */
/* viri (triangles being eaten), encroached segments, bad (skinny or too */
/* large) triangles, and splay tree nodes. */
struct memorypool triangles;
struct memorypool subsegs;
struct memorypool vertices;
struct memorypool viri;
struct memorypool badsubsegs;
struct memorypool badtriangles;
struct memorypool flipstackers;
struct memorypool splaynodes;
/* Variables that maintain the bad triangle queues. The queues are */
/* ordered from 4095 (highest priority) to 0 (lowest priority). */
struct badtriang *queuefront[4096];
struct badtriang *queuetail[4096];
int nextnonemptyq[4096];
int firstnonemptyq;
/* Variable that maintains the stack of recently flipped triangles. */
struct flipstacker *lastflip;
/* Other variables. */
REAL xmin, xmax, ymin, ymax; /* x and y bounds. */
REAL xminextreme; /* Nonexistent x value used as a flag in sweepline. */
int invertices; /* Number of input vertices. */
int inelements; /* Number of input triangles. */
int insegments; /* Number of input segments. */
int holes; /* Number of input holes. */
int regions; /* Number of input regions. */
int undeads; /* Number of input vertices that don't appear in the mesh. */
LONG edges; /* Number of output edges. */
int mesh_dim; /* Dimension (ought to be 2). */
int nextras; /* Number of attributes per vertex. */
int eextras; /* Number of attributes per triangle. */
LONG hullsize; /* Number of edges in convex hull. */
int steinerleft; /* Number of Steiner points not yet used. */
int vertexmarkindex; /* Index to find boundary marker of a vertex. */
int vertex2triindex; /* Index to find a triangle adjacent to a vertex. */
int highorderindex; /* Index to find extra nodes for high-order elements. */
int elemattribindex; /* Index to find attributes of a triangle. */
int areaboundindex; /* Index to find area bound of a triangle. */
int checksegments; /* Are there segments in the triangulation yet? */
int checkquality; /* Has quality triangulation begun yet? */
int readnodefile; /* Has a .node file been read? */
LONG samples; /* Number of random samples for point location. */
LONG incirclecount; /* Number of incircle tests performed. */
LONG counterclockcount; /* Number of counterclockwise tests performed. */
LONG orient3dcount; /* Number of 3D orientation tests performed. */
LONG hyperbolacount; /* Number of right-of-hyperbola tests performed. */
LONG circumcentercount; /* Number of circumcenter calculations performed. */
LONG circletopcount; /* Number of circle top calculations performed. */
/* Triangular bounding box vertices. */
vertex infvertex1, infvertex2, infvertex3;
/* Pointer to the `triangle' that occupies all of "outer space." */
triangle *dummytri;
triangle *dummytribase; /* Keep base address so we can free() it later. */
/* Pointer to the omnipresent subsegment. Referenced by any triangle or */
/* subsegment that isn't really connected to a subsegment at that */
/* location. */
subseg *dummysub;
subseg *dummysubbase; /* Keep base address so we can free() it later. */
/* Pointer to a recently visited triangle. Improves point location if */
/* proximate vertices are inserted sequentially. */
struct otri recenttri;
}; /* End of `struct mesh'. */
/* Data structure for command line switches and file names. This structure */
/* is used (instead of global variables) to allow reentrancy. */
struct behavior {
/* Switches for the triangulator. */
/* poly: -p switch. refine: -r switch. */
/* quality: -q switch. */
/* minangle: minimum angle bound, specified after -q switch. */
/* goodangle: cosine squared of minangle. */
/* offconstant: constant used to place off-center Steiner points. */
/* vararea: -a switch without number. */
/* fixedarea: -a switch with number. */
/* maxarea: maximum area bound, specified after -a switch. */
/* usertest: -u switch. */
/* regionattrib: -A switch. convex: -c switch. */
/* weighted: 1 for -w switch, 2 for -W switch. jettison: -j switch */
/* firstnumber: inverse of -z switch. All items are numbered starting */
/* from `firstnumber'. */
/* edgesout: -e switch. voronoi: -v switch. */
/* neighbors: -n switch. geomview: -g switch. */
/* nobound: -B switch. nopolywritten: -P switch. */
/* nonodewritten: -N switch. noelewritten: -E switch. */
/* noiterationnum: -I switch. noholes: -O switch. */
/* noexact: -X switch. */
/* order: element order, specified after -o switch. */
/* nobisect: count of how often -Y switch is selected. */
/* steiner: maximum number of Steiner points, specified after -S switch. */
/* incremental: -i switch. sweepline: -F switch. */
/* dwyer: inverse of -l switch. */
/* splitseg: -s switch. */
/* conformdel: -D switch. docheck: -C switch. */
/* quiet: -Q switch. verbose: count of how often -V switch is selected. */
/* usesegments: -p, -r, -q, or -c switch; determines whether segments are */
/* used at all. */
/* */
/* Read the instructions to find out the meaning of these switches. */
int poly, refine, quality, vararea, fixedarea, usertest;
int regionattrib, convex, weighted, jettison;
int firstnumber;
int edgesout, voronoi, neighbors, geomview;
int nobound, nopolywritten, nonodewritten, noelewritten, noiterationnum;
int noholes, noexact, conformdel;
int incremental, sweepline, dwyer;
int splitseg;
int docheck;
int quiet, verbose;
int usesegments;
int order;
int nobisect;
int steiner;
REAL minangle, goodangle, offconstant;
REAL maxarea;
/* Variables for file names. */
#ifndef TRILIBRARY
char innodefilename[FILENAMESIZE];
char inelefilename[FILENAMESIZE];
char inpolyfilename[FILENAMESIZE];
char areafilename[FILENAMESIZE];
char outnodefilename[FILENAMESIZE];
char outelefilename[FILENAMESIZE];
char outpolyfilename[FILENAMESIZE];
char edgefilename[FILENAMESIZE];
char vnodefilename[FILENAMESIZE];
char vedgefilename[FILENAMESIZE];
char neighborfilename[FILENAMESIZE];
char offfilename[FILENAMESIZE];
#endif /* not TRILIBRARY */
}; /* End of `struct behavior'. */
/*****************************************************************************/
/* */
/* Mesh manipulation primitives. Each triangle contains three pointers to */
/* other triangles, with orientations. Each pointer points not to the */
/* first byte of a triangle, but to one of the first three bytes of a */
/* triangle. It is necessary to extract both the triangle itself and the */
/* orientation. To save memory, I keep both pieces of information in one */
/* pointer. To make this possible, I assume that all triangles are aligned */
/* to four-byte boundaries. The decode() routine below decodes a pointer, */
/* extracting an orientation (in the range 0 to 2) and a pointer to the */
/* beginning of a triangle. The encode() routine compresses a pointer to a */
/* triangle and an orientation into a single pointer. My assumptions that */
/* triangles are four-byte-aligned and that the `unsigned long' type is */
/* long enough to hold a pointer are two of the few kludges in this program.*/
/* */
/* Subsegments are manipulated similarly. A pointer to a subsegment */
/* carries both an address and an orientation in the range 0 to 1. */
/* */
/* The other primitives take an oriented triangle or oriented subsegment, */
/* and return an oriented triangle or oriented subsegment or vertex; or */
/* they change the connections in the data structure. */
/* */
/* Below, triangles and subsegments are denoted by their vertices. The */
/* triangle abc has origin (org) a, destination (dest) b, and apex (apex) */
/* c. These vertices occur in counterclockwise order about the triangle. */
/* The handle abc may simultaneously denote vertex a, edge ab, and triangle */
/* abc. */
/* */
/* Similarly, the subsegment ab has origin (sorg) a and destination (sdest) */
/* b. If ab is thought to be directed upward (with b directly above a), */
/* then the handle ab is thought to grasp the right side of ab, and may */
/* simultaneously denote vertex a and edge ab. */
/* */
/* An asterisk (*) denotes a vertex whose identity is unknown. */
/* */
/* Given this notation, a partial list of mesh manipulation primitives */
/* follows. */
/* */
/* */
/* For triangles: */
/* */
/* sym: Find the abutting triangle; same edge. */
/* sym(abc) -> ba* */
/* */
/* lnext: Find the next edge (counterclockwise) of a triangle. */
/* lnext(abc) -> bca */
/* */
/* lprev: Find the previous edge (clockwise) of a triangle. */
/* lprev(abc) -> cab */
/* */
/* onext: Find the next edge counterclockwise with the same origin. */
/* onext(abc) -> ac* */
/* */
/* oprev: Find the next edge clockwise with the same origin. */
/* oprev(abc) -> a*b */
/* */
/* dnext: Find the next edge counterclockwise with the same destination. */
/* dnext(abc) -> *ba */
/* */
/* dprev: Find the next edge clockwise with the same destination. */
/* dprev(abc) -> cb* */
/* */
/* rnext: Find the next edge (counterclockwise) of the adjacent triangle. */
/* rnext(abc) -> *a* */
/* */
/* rprev: Find the previous edge (clockwise) of the adjacent triangle. */
/* rprev(abc) -> b** */
/* */
/* org: Origin dest: Destination apex: Apex */
/* org(abc) -> a dest(abc) -> b apex(abc) -> c */
/* */
/* bond: Bond two triangles together at the respective handles. */
/* bond(abc, bad) */
/* */
/* */
/* For subsegments: */
/* */
/* ssym: Reverse the orientation of a subsegment. */
/* ssym(ab) -> ba */
/* */
/* spivot: Find adjoining subsegment with the same origin. */
/* spivot(ab) -> a* */
/* */
/* snext: Find next subsegment in sequence. */
/* snext(ab) -> b* */
/* */
/* sorg: Origin sdest: Destination */
/* sorg(ab) -> a sdest(ab) -> b */
/* */
/* sbond: Bond two subsegments together at the respective origins. */
/* sbond(ab, ac) */
/* */
/* */
/* For interacting tetrahedra and subfacets: */
/* */
/* tspivot: Find a subsegment abutting a triangle. */
/* tspivot(abc) -> ba */
/* */
/* stpivot: Find a triangle abutting a subsegment. */
/* stpivot(ab) -> ba* */
/* */
/* tsbond: Bond a triangle to a subsegment. */
/* tsbond(abc, ba) */
/* */
/*****************************************************************************/
/********* Mesh manipulation primitives begin here *********/
/** **/
/** **/
/* Fast lookup arrays to speed some of the mesh manipulation primitives. */
static int plus1mod3[3] = {1, 2, 0};
static int minus1mod3[3] = {2, 0, 1};
/********* Primitives for triangles *********/
/* */
/* */
/* decode() converts a pointer to an oriented triangle. The orientation is */
/* extracted from the two least significant bits of the pointer. */
#define decode(ptr, otri) \
(otri).orient = (int) ((unsigned LONG) (ptr) & (unsigned LONG) 3l); \
(otri).tri = (triangle *) \
((unsigned LONG) (ptr) ^ (unsigned LONG) (otri).orient)
/* encode() compresses an oriented triangle into a single pointer. It */
/* relies on the assumption that all triangles are aligned to four-byte */
/* boundaries, so the two least significant bits of (otri).tri are zero. */
#define encode(otri) \
(triangle) ((unsigned LONG) (otri).tri | (unsigned LONG) (otri).orient)
/* The following handle manipulation primitives are all described by Guibas */
/* and Stolfi. However, Guibas and Stolfi use an edge-based data */
/* structure, whereas I use a triangle-based data structure. */
/* sym() finds the abutting triangle, on the same edge. Note that the edge */
/* direction is necessarily reversed, because the handle specified by an */
/* oriented triangle is directed counterclockwise around the triangle. */
#define sym(otri1, otri2) \
ptr = (otri1).tri[(otri1).orient]; \
decode(ptr, otri2);
#define symself(otri) \
ptr = (otri).tri[(otri).orient]; \
decode(ptr, otri);
/* lnext() finds the next edge (counterclockwise) of a triangle. */
#define lnext(otri1, otri2) \
(otri2).tri = (otri1).tri; \