00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201 #define FILENAMESIZE 2048
00202
00203
00204
00205
00206 #define INPUTLINESIZE 1024
00207
00208
00209
00210
00211
00212 #define TRIPERBLOCK 4092
00213 #define SUBSEGPERBLOCK 508
00214 #define VERTEXPERBLOCK 4092
00215 #define VIRUSPERBLOCK 1020
00216 #define BADSUBSEGPERBLOCK 252
00217 #define BADTRIPERBLOCK 4092
00218 #define FLIPSTACKERPERBLOCK 252
00219 #define SPLAYNODEPERBLOCK 508
00220
00221
00222
00223
00224
00225 #define INPUTVERTEX 0
00226 #define SEGMENTVERTEX 1
00227 #define FREEVERTEX 2
00228 #define DEADVERTEX -32768
00229 #define UNDEADVERTEX -32767
00230
00231
00232
00233
00234
00235
00236
00237 #define SAMPLEFACTOR 11
00238
00239
00240
00241
00242
00243 #define SAMPLERATE 10
00244
00245
00246
00247 #define PI 3.141592653589793238462643383279502884197169399375105820974944592308
00248
00249
00250
00251 #define SQUAREROOTTWO 1.4142135623730950488016887242096980785696718753769480732
00252
00253
00254
00255 #define ONETHIRD 0.333333333333333333333333333333333333333333333333333333333333
00256
00257 #include <stdio.h>
00258 #include <stdlib.h>
00259 #include <string.h>
00260 #include <math.h>
00261
00262 #include "triangle.h"
00263
00264
00265
00266
00267
00268 enum locateresult {INTRIANGLE, ONEDGE, ONVERTEX, OUTSIDE};
00269
00270
00271
00272
00273
00274
00275
00276 enum insertvertexresult {SUCCESSFULVERTEX, ENCROACHINGVERTEX, VIOLATINGVERTEX,
00277 DUPLICATEVERTEX};
00278
00279
00280
00281
00282
00283
00284 enum finddirectionresult {WITHIN, LEFTCOLLINEAR, RIGHTCOLLINEAR};
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403 typedef float **triangle;
00404
00405
00406
00407
00408
00409
00410 struct otri {
00411 triangle *tri;
00412 int orient;
00413 };
00414
00415
00416
00417
00418
00419
00420 typedef float **subseg;
00421
00422
00423
00424
00425
00426
00427 struct osub {
00428 subseg *ss;
00429 int ssorient;
00430 };
00431
00432
00433
00434
00435
00436
00437 typedef float *vertex;
00438
00439
00440
00441
00442 struct badsubseg {
00443 subseg encsubseg;
00444 vertex subsegorg, subsegdest;
00445 };
00446
00447
00448
00449
00450
00451 struct badtriang {
00452 triangle poortri;
00453 float key;
00454 vertex triangorg, triangdest, triangapex;
00455 struct badtriang *nexttriang;
00456 };
00457
00458
00459
00460
00461
00462 struct flipstacker {
00463 triangle flippedtri;
00464 struct flipstacker *prevflip;
00465 };
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477 struct event {
00478 float xkey, ykey;
00479 int *eventptr;
00480 int heapposition;
00481 };
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494 struct splaynode {
00495 struct otri keyedge;
00496 vertex keydest;
00497 struct splaynode *lchild, *rchild;
00498 };
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523 struct memorypool {
00524 int **firstblock, **nowblock;
00525 int *nextitem;
00526 int *deaditemstack;
00527 int **pathblock;
00528 int *pathitem;
00529 int alignbytes;
00530 int itembytes;
00531 int itemsperblock;
00532 int itemsfirstblock;
00533 long items, maxitems;
00534 int unallocateditems;
00535 int pathitemsleft;
00536 };
00537
00538
00539
00540
00541 float splitter;
00542 float epsilon;
00543 float resulterrbound;
00544 float ccwerrboundA, ccwerrboundB, ccwerrboundC;
00545 float iccerrboundA, iccerrboundB, iccerrboundC;
00546 float o3derrboundA, o3derrboundB, o3derrboundC;
00547
00548
00549
00550 unsigned long randomseed;
00551
00552
00553
00554
00555
00556 struct mesh {
00557
00558
00559
00560
00561
00562 struct memorypool triangles;
00563 struct memorypool subsegs;
00564 struct memorypool vertices;
00565 struct memorypool viri;
00566 struct memorypool badsubsegs;
00567 struct memorypool badtriangles;
00568 struct memorypool flipstackers;
00569 struct memorypool splaynodes;
00570
00571
00572
00573
00574 struct badtriang *queuefront[4096];
00575 struct badtriang *queuetail[4096];
00576 int nextnonemptyq[4096];
00577 int firstnonemptyq;
00578
00579
00580
00581 struct flipstacker *lastflip;
00582
00583
00584
00585 float xmin, xmax, ymin, ymax;
00586 float xminextreme;
00587 int invertices;
00588 int inelements;
00589 int insegments;
00590 int holes;
00591 int regions;
00592 int undeads;
00593 long edges;
00594 int mesh_dim;
00595 int nextras;
00596 int eextras;
00597 long hullsize;
00598 int steinerleft;
00599 int vertexmarkindex;
00600 int vertex2triindex;
00601 int highorderindex;
00602 int elemattribindex;
00603 int areaboundindex;
00604 int checksegments;
00605 int checkquality;
00606 int readnodefile;
00607 long samples;
00608
00609 long incirclecount;
00610 long counterclockcount;
00611 long orient3dcount;
00612 long hyperbolacount;
00613 long circumcentercount;
00614 long circletopcount;
00615
00616
00617
00618 vertex infvertex1, infvertex2, infvertex3;
00619
00620
00621
00622 triangle *dummytri;
00623 triangle *dummytribase;
00624
00625
00626
00627
00628
00629 subseg *dummysub;
00630 subseg *dummysubbase;
00631
00632
00633
00634
00635 struct otri recenttri;
00636
00637 };
00638
00639
00640
00641
00642
00643 struct behavior {
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655
00656
00657
00658
00659
00660
00661
00662
00663
00664
00665
00666
00667
00668
00669
00670
00671
00672
00673
00674
00675
00676
00677
00678 int poly, refine, quality, vararea, fixedarea, usertest;
00679 int regionattrib, convex, weighted, jettison;
00680 int firstnumber;
00681 int edgesout, voronoi, neighbors, geomview;
00682 int nobound, nopolywritten, nonodewritten, noelewritten, noiterationnum;
00683 int noholes, noexact, conformdel;
00684 int incremental, sweepline, dwyer;
00685 int splitseg;
00686 int docheck;
00687 int quiet, verbose;
00688 int usesegments;
00689 int order;
00690 int nobisect;
00691 int steiner;
00692 float minangle, goodangle, offconstant;
00693 float maxarea;
00694
00695
00696
00697 };
00698
00699
00700
00701
00702
00703
00704
00705
00706
00707
00708
00709
00710
00711
00712
00713
00714
00715
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729
00730
00731
00732
00733
00734
00735
00736
00737
00738
00739
00740
00741
00742
00743
00744
00745
00746
00747
00748
00749
00750
00751
00752
00753
00754
00755
00756
00757
00758
00759
00760
00761
00762
00763
00764
00765
00766
00767
00768
00769
00770
00771
00772
00773
00774
00775
00776
00777
00778
00779
00780
00781
00782
00783
00784
00785
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796
00797
00798
00799
00800
00801
00802
00803
00804
00805
00806
00810
00811
00812 int plus1mod3[3] = {1, 2, 0};
00813 int minus1mod3[3] = {2, 0, 1};
00814
00815
00816
00817
00818
00819
00820
00821
00822 #define decode(ptr, otri) \
00823 (otri).orient = (int) ((unsigned long) (ptr) & (unsigned long) 3l); \
00824 (otri).tri = (triangle *) \
00825 ((unsigned long) (ptr) ^ (unsigned long) (otri).orient)
00826
00827
00828
00829
00830
00831 #define encode(otri) \
00832 (triangle) ((unsigned long) (otri).tri | (unsigned long) (otri).orient)
00833
00834
00835
00836
00837
00838
00839
00840
00841
00842 #define sym(otri1, otri2) \
00843 ptr = (otri1).tri[(otri1).orient]; \
00844 decode(ptr, otri2);
00845
00846 #define symself(otri) \
00847 ptr = (otri).tri[(otri).orient]; \
00848 decode(ptr, otri);
00849
00850
00851
00852 #define lnext(otri1, otri2) \
00853 (otri2).tri = (otri1).tri; \
00854 (otri2).orient = plus1mod3[(otri1).orient]
00855
00856 #define lnextself(otri) \
00857 (otri).orient = plus1mod3[(otri).orient]
00858
00859
00860
00861 #define lprev(otri1, otri2) \
00862 (otri2).tri = (otri1).tri; \
00863 (otri2).orient = minus1mod3[(otri1).orient]
00864
00865 #define lprevself(otri) \
00866 (otri).orient = minus1mod3[(otri).orient]
00867
00868
00869
00870
00871
00872 #define onext(otri1, otri2) \
00873 lprev(otri1, otri2); \
00874 symself(otri2);
00875
00876 #define onextself(otri) \
00877 lprevself(otri); \
00878 symself(otri);
00879
00880
00881
00882
00883
00884 #define oprev(otri1, otri2) \
00885 sym(otri1, otri2); \
00886 lnextself(otri2);
00887
00888 #define oprevself(otri) \
00889 symself(otri); \
00890 lnextself(otri);
00891
00892
00893
00894
00895
00896 #define dnext(otri1, otri2) \
00897 sym(otri1, otri2); \
00898 lprevself(otri2);
00899
00900 #define dnextself(otri) \
00901 symself(otri); \
00902 lprevself(otri);
00903
00904
00905
00906
00907
00908 #define dprev(otri1, otri2) \
00909 lnext(otri1, otri2); \
00910 symself(otri2);
00911
00912 #define dprevself(otri) \
00913 lnextself(otri); \
00914 symself(otri);
00915
00916
00917
00918
00919
00920 #define rnext(otri1, otri2) \
00921 sym(otri1, otri2); \
00922 lnextself(otri2); \
00923 symself(otri2);
00924
00925 #define rnextself(otri) \
00926 symself(otri); \
00927 lnextself(otri); \
00928 symself(otri);
00929
00930
00931
00932
00933
00934 #define rprev(otri1, otri2) \
00935 sym(otri1, otri2); \
00936 lprevself(otri2); \
00937 symself(otri2);
00938
00939 #define rprevself(otri) \
00940 symself(otri); \
00941 lprevself(otri); \
00942 symself(otri);
00943
00944
00945
00946
00947 #define org(otri, vertexptr) \
00948 vertexptr = (vertex) (otri).tri[plus1mod3[(otri).orient] + 3]
00949
00950 #define dest(otri, vertexptr) \
00951 vertexptr = (vertex) (otri).tri[minus1mod3[(otri).orient] + 3]
00952
00953 #define apex(otri, vertexptr) \
00954 vertexptr = (vertex) (otri).tri[(otri).orient + 3]
00955
00956 #define setorg(otri, vertexptr) \
00957 (otri).tri[plus1mod3[(otri).orient] + 3] = (triangle) vertexptr
00958
00959 #define setdest(otri, vertexptr) \
00960 (otri).tri[minus1mod3[(otri).orient] + 3] = (triangle) vertexptr
00961
00962 #define setapex(otri, vertexptr) \
00963 (otri).tri[(otri).orient + 3] = (triangle) vertexptr
00964
00965
00966
00967 #define bond(otri1, otri2) \
00968 (otri1).tri[(otri1).orient] = encode(otri2); \
00969 (otri2).tri[(otri2).orient] = encode(otri1)
00970
00971
00972
00973
00974
00975
00976 #define dissolve(otri) \
00977 (otri).tri[(otri).orient] = (triangle) m->dummytri
00978
00979
00980
00981 #define otricopy(otri1, otri2) \
00982 (otri2).tri = (otri1).tri; \
00983 (otri2).orient = (otri1).orient
00984
00985
00986
00987 #define otriequal(otri1, otri2) \
00988 (((otri1).tri == (otri2).tri) && \
00989 ((otri1).orient == (otri2).orient))
00990
00991
00992
00993
00994 #define infect(otri) \
00995 (otri).tri[6] = (triangle) \
00996 ((unsigned long) (otri).tri[6] | (unsigned long) 2l)
00997
00998 #define uninfect(otri) \
00999 (otri).tri[6] = (triangle) \
01000 ((unsigned long) (otri).tri[6] & ~ (unsigned long) 2l)
01001
01002
01003
01004 #define infected(otri) \
01005 (((unsigned long) (otri).tri[6] & (unsigned long) 2l) != 0l)
01006
01007
01008
01009 #define elemattribute(otri, attnum) \
01010 ((float *) (otri).tri)[m->elemattribindex + (attnum)]
01011
01012 #define setelemattribute(otri, attnum, value) \
01013 ((float *) (otri).tri)[m->elemattribindex + (attnum)] = value
01014
01015
01016
01017 #define areabound(otri) ((float *) (otri).tri)[m->areaboundindex]
01018
01019 #define setareabound(otri, value) \
01020 ((float *) (otri).tri)[m->areaboundindex] = value
01021
01022
01023
01024
01025
01026
01027 #define deadtri(tria) ((tria)[1] == (triangle) NULL)
01028
01029 #define killtri(tria) \
01030 (tria)[1] = (triangle) NULL; \
01031 (tria)[3] = (triangle) NULL
01032
01033
01034
01035
01036
01037
01038
01039
01040
01041
01042 #define sdecode(sptr, osub) \
01043 (osub).ssorient = (int) ((unsigned long) (sptr) & (unsigned long) 1l); \
01044 (osub).ss = (subseg *) \
01045 ((unsigned long) (sptr) & ~ (unsigned long) 3l)
01046
01047
01048
01049
01050
01051 #define sencode(osub) \
01052 (subseg) ((unsigned long) (osub).ss | (unsigned long) (osub).ssorient)
01053
01054
01055
01056 #define ssym(osub1, osub2) \
01057 (osub2).ss = (osub1).ss; \
01058 (osub2).ssorient = 1 - (osub1).ssorient
01059
01060 #define ssymself(osub) \
01061 (osub).ssorient = 1 - (osub).ssorient
01062
01063
01064
01065
01066 #define spivot(osub1, osub2) \
01067 sptr = (osub1).ss[(osub1).ssorient]; \
01068 sdecode(sptr, osub2)
01069
01070 #define spivotself(osub) \
01071 sptr = (osub).ss[(osub).ssorient]; \
01072 sdecode(sptr, osub)
01073
01074
01075
01076
01077 #define snext(osub1, osub2) \
01078 sptr = (osub1).ss[1 - (osub1).ssorient]; \
01079 sdecode(sptr, osub2)
01080
01081 #define snextself(osub) \
01082 sptr = (osub).ss[1 - (osub).ssorient]; \
01083 sdecode(sptr, osub)
01084
01085
01086
01087
01088 #define sorg(osub, vertexptr) \
01089 vertexptr = (vertex) (osub).ss[2 + (osub).ssorient]
01090
01091 #define sdest(osub, vertexptr) \
01092 vertexptr = (vertex) (osub).ss[3 - (osub).ssorient]
01093
01094 #define setsorg(osub, vertexptr) \
01095 (osub).ss[2 + (osub).ssorient] = (subseg) vertexptr
01096
01097 #define setsdest(osub, vertexptr) \
01098 (osub).ss[3 - (osub).ssorient] = (subseg) vertexptr
01099
01100 #define segorg(osub, vertexptr) \
01101 vertexptr = (vertex) (osub).ss[4 + (osub).ssorient]
01102
01103 #define segdest(osub, vertexptr) \
01104 vertexptr = (vertex) (osub).ss[5 - (osub).ssorient]
01105
01106 #define setsegorg(osub, vertexptr) \
01107 (osub).ss[4 + (osub).ssorient] = (subseg) vertexptr
01108
01109 #define setsegdest(osub, vertexptr) \
01110 (osub).ss[5 - (osub).ssorient] = (subseg) vertexptr
01111
01112
01113
01114
01115
01116 #define mark(osub) (* (int *) ((osub).ss + 8))
01117
01118 #define setmark(osub, value) \
01119 * (int *) ((osub).ss + 8) = value
01120
01121
01122
01123 #define sbond(osub1, osub2) \
01124 (osub1).ss[(osub1).ssorient] = sencode(osub2); \
01125 (osub2).ss[(osub2).ssorient] = sencode(osub1)
01126
01127
01128
01129
01130 #define sdissolve(osub) \
01131 (osub).ss[(osub).ssorient] = (subseg) m->dummysub
01132
01133
01134
01135 #define subsegcopy(osub1, osub2) \
01136 (osub2).ss = (osub1).ss; \
01137 (osub2).ssorient = (osub1).ssorient
01138
01139
01140
01141 #define subsegequal(osub1, osub2) \
01142 (((osub1).ss == (osub2).ss) && \
01143 ((osub1).ssorient == (osub2).ssorient))
01144
01145
01146
01147
01148
01149
01150 #define deadsubseg(sub) ((sub)[1] == (subseg) NULL)
01151
01152 #define killsubseg(sub) \
01153 (sub)[1] = (subseg) NULL; \
01154 (sub)[2] = (subseg) NULL
01155
01156
01157
01158
01159
01160
01161
01162 #define tspivot(otri, osub) \
01163 sptr = (subseg) (otri).tri[6 + (otri).orient]; \
01164 sdecode(sptr, osub)
01165
01166
01167
01168
01169 #define stpivot(osub, otri) \
01170 ptr = (triangle) (osub).ss[6 + (osub).ssorient]; \
01171 decode(ptr, otri)
01172
01173
01174
01175 #define tsbond(otri, osub) \
01176 (otri).tri[6 + (otri).orient] = (triangle) sencode(osub); \
01177 (osub).ss[6 + (osub).ssorient] = (subseg) encode(otri)
01178
01179
01180
01181 #define tsdissolve(otri) \
01182 (otri).tri[6 + (otri).orient] = (triangle) m->dummysub
01183
01184
01185
01186 #define stdissolve(osub) \
01187 (osub).ss[6 + (osub).ssorient] = (subseg) m->dummytri
01188
01189
01190
01191
01192
01193 #define vertexmark(vx) ((int *) (vx))[m->vertexmarkindex]
01194
01195 #define setvertexmark(vx, value) \
01196 ((int *) (vx))[m->vertexmarkindex] = value
01197
01198 #define vertextype(vx) ((int *) (vx))[m->vertexmarkindex + 1]
01199
01200 #define setvertextype(vx, value) \
01201 ((int *) (vx))[m->vertexmarkindex + 1] = value
01202
01203 #define vertex2tri(vx) ((triangle *) (vx))[m->vertex2triindex]
01204
01205 #define setvertex2tri(vx, value) \
01206 ((triangle *) (vx))[m->vertex2triindex] = value
01207
01210
01211
01212
01216 void triexit(int status)
01217 {
01218 exit(status);
01219 }
01220
01221 int *trimalloc(int size)
01222 {
01223 int *memptr;
01224
01225 memptr = (int *) malloc((unsigned int) size);
01226 if (memptr == (int *) NULL) {
01227 printf("Error: Out of memory.\n");
01228 triexit(1);
01229 }
01230 return(memptr);
01231 }
01232
01233 void trifree(int *memptr)
01234 {
01235 free(memptr);
01236 }
01237
01240
01241
01242
01243
01244
01245
01246
01247
01248 void internalerror()
01249 {
01250 printf(" Please report this bug to jrs@cs.berkeley.edu\n");
01251 printf(" Include the message above, your input data set, and the exact\n");
01252 printf(" command line you used to run Triangle.\n");
01253 triexit(1);
01254 }
01255
01256
01257
01258
01259
01260
01261
01262
01263 void parsecommandline(int argc, char **argv, struct behavior *b) {
01264 int i, j;
01265
01266 b->poly = b->refine = b->quality = 0;
01267 b->vararea = b->fixedarea = b->usertest = 0;
01268 b->regionattrib = b->convex = b->weighted = b->jettison = 0;
01269 b->firstnumber = 1;
01270 b->edgesout = b->voronoi = b->neighbors = b->geomview = 0;
01271 b->nobound = b->nopolywritten = b->nonodewritten = b->noelewritten = 0;
01272 b->noiterationnum = 0;
01273 b->noholes = b->noexact = 0;
01274 b->incremental = b->sweepline = 0;
01275 b->dwyer = 1;
01276 b->splitseg = 0;
01277 b->docheck = 0;
01278 b->nobisect = 0;
01279 b->conformdel = 0;
01280 b->steiner = -1;
01281 b->order = 1;
01282 b->minangle = 0.0;
01283 b->maxarea = -1.0;
01284 b->quiet = b->verbose = 0;
01285
01286 for (i = 0; i < argc; i++) {
01287 for (j = 0; argv[i][j] != '\0'; j++) {
01288 if (argv[i][j] == 'p') {
01289 b->poly = 1;
01290 }
01291 if (argv[i][j] == 'A') {
01292 b->regionattrib = 1;
01293 }
01294 if (argv[i][j] == 'c') {
01295 b->convex = 1;
01296 }
01297 if (argv[i][j] == 'w') {
01298 b->weighted = 1;
01299 }
01300 if (argv[i][j] == 'W') {
01301 b->weighted = 2;
01302 }
01303 if (argv[i][j] == 'j') {
01304 b->jettison = 1;
01305 }
01306 if (argv[i][j] == 'z') {
01307 b->firstnumber = 0;
01308 }
01309 if (argv[i][j] == 'e') {
01310 b->edgesout = 1;
01311 }
01312 if (argv[i][j] == 'v') {
01313 b->voronoi = 1;
01314 }
01315 if (argv[i][j] == 'n') {
01316 b->neighbors = 1;
01317 }
01318 if (argv[i][j] == 'g') {
01319 b->geomview = 1;
01320 }
01321 if (argv[i][j] == 'B') {
01322 b->nobound = 1;
01323 }
01324 if (argv[i][j] == 'P') {
01325 b->nopolywritten = 1;
01326 }
01327 if (argv[i][j] == 'N') {
01328 b->nonodewritten = 1;
01329 }
01330 if (argv[i][j] == 'E') {
01331 b->noelewritten = 1;
01332 }
01333 if (argv[i][j] == 'O') {
01334 b->noholes = 1;
01335 }
01336 if (argv[i][j] == 'X') {
01337 b->noexact = 1;
01338 }
01339 if (argv[i][j] == 'o') {
01340 if (argv[i][j + 1] == '2') {
01341 j++;
01342 b->order = 2;
01343 }
01344 }
01345 if (argv[i][j] == 'l') {
01346 b->dwyer = 0;
01347 }
01348 if (argv[i][j] == 'Q') {
01349 b->quiet = 1;
01350 }
01351 if (argv[i][j] == 'V') {
01352 b->verbose++;
01353 }
01354 }
01355 }
01356 b->usesegments = b->poly || b->refine || b->quality || b->convex;
01357 b->goodangle = cos(b->minangle * PI / 180.0);
01358 if (b->goodangle == 1.0) {
01359 b->offconstant = 0.0;
01360 } else {
01361 b->offconstant = 0.475 * sqrt((1.0 + b->goodangle) / (1.0 - b->goodangle));
01362 }
01363 b->goodangle *= b->goodangle;
01364 if (b->refine && b->noiterationnum) {
01365 printf(
01366 "Error: You cannot use the -I switch when refining a triangulation.\n");
01367 triexit(1);
01368 }
01369
01370
01371 if (!b->refine && !b->poly) {
01372 b->vararea = 0;
01373 }
01374
01375
01376 if (b->refine || !b->poly) {
01377 b->regionattrib = 0;
01378 }
01379
01380
01381 if (b->weighted && (b->poly || b->quality)) {
01382 b->weighted = 0;
01383 if (!b->quiet) {
01384 printf("Warning: weighted triangulations (-w, -W) are incompatible\n");
01385 printf(" with PSLGs (-p) and meshing (-q, -a, -u). Weights ignored.\n"
01386 );
01387 }
01388 }
01389 if (b->jettison && b->nonodewritten && !b->quiet) {
01390 printf("Warning: -j and -N switches are somewhat incompatible.\n");
01391 printf(" If any vertices are jettisoned, you will need the output\n");
01392 printf(" .node file to reconstruct the new node indices.");
01393 }
01394 }
01395
01398
01399
01400
01404
01405
01406
01407
01408
01409
01410
01411
01412
01413
01414
01415 void printtriangle(struct mesh *m, struct behavior *b, struct otri *t)
01416 {
01417 struct otri printtri;
01418 struct osub printsh;
01419 vertex printvertex;
01420
01421 printf("triangle x%lx with orientation %d:\n", (unsigned long) t->tri,
01422 t->orient);
01423 decode(t->tri[0], printtri);
01424 if (printtri.tri == m->dummytri) {
01425 printf(" [0] = Outer space\n");
01426 } else {
01427 printf(" [0] = x%lx %d\n", (unsigned long) printtri.tri,
01428 printtri.orient);
01429 }
01430 decode(t->tri[1], printtri);
01431 if (printtri.tri == m->dummytri) {
01432 printf(" [1] = Outer space\n");
01433 } else {
01434 printf(" [1] = x%lx %d\n", (unsigned long) printtri.tri,
01435 printtri.orient);
01436 }
01437 decode(t->tri[2], printtri);
01438 if (printtri.tri == m->dummytri) {
01439 printf(" [2] = Outer space\n");
01440 } else {
01441 printf(" [2] = x%lx %d\n", (unsigned long) printtri.tri,
01442 printtri.orient);
01443 }
01444
01445 org(*t, printvertex);
01446 if (printvertex == (vertex) NULL)
01447 printf(" Origin[%d] = NULL\n", (t->orient + 1) % 3 + 3);
01448 else
01449 printf(" Origin[%d] = x%lx (%.12g, %.12g)\n",
01450 (t->orient + 1) % 3 + 3, (unsigned long) printvertex,
01451 printvertex[0], printvertex[1]);
01452 dest(*t, printvertex);
01453 if (printvertex == (vertex) NULL)
01454 printf(" Dest [%d] = NULL\n", (t->orient + 2) % 3 + 3);
01455 else
01456 printf(" Dest [%d] = x%lx (%.12g, %.12g)\n",
01457 (t->orient + 2) % 3 + 3, (unsigned long) printvertex,
01458 printvertex[0], printvertex[1]);
01459 apex(*t, printvertex);
01460 if (printvertex == (vertex) NULL)
01461 printf(" Apex [%d] = NULL\n", t->orient + 3);
01462 else
01463 printf(" Apex [%d] = x%lx (%.12g, %.12g)\n",
01464 t->orient + 3, (unsigned long) printvertex,
01465 printvertex[0], printvertex[1]);
01466
01467 if (b->usesegments) {
01468 sdecode(t->tri[6], printsh);
01469 if (printsh.ss != m->dummysub) {
01470 printf(" [6] = x%lx %d\n", (unsigned long) printsh.ss,
01471 printsh.ssorient);
01472 }
01473 sdecode(t->tri[7], printsh);
01474 if (printsh.ss != m->dummysub) {
01475 printf(" [7] = x%lx %d\n", (unsigned long) printsh.ss,
01476 printsh.ssorient);
01477 }
01478 sdecode(t->tri[8], printsh);
01479 if (printsh.ss != m->dummysub) {
01480 printf(" [8] = x%lx %d\n", (unsigned long) printsh.ss,
01481 printsh.ssorient);
01482 }
01483 }
01484
01485 if (b->vararea) {
01486 printf(" Area constraint: %.4g\n", areabound(*t));
01487 }
01488 }
01489
01490
01491
01492
01493
01494
01495
01496
01497
01498
01499
01500
01501 void printsubseg(struct mesh *m, struct behavior *b, struct osub *s)
01502 {
01503 struct osub printsh;
01504 struct otri printtri;
01505 vertex printvertex;
01506
01507 printf("subsegment x%lx with orientation %d and mark %d:\n",
01508 (unsigned long) s->ss, s->ssorient, mark(*s));
01509 sdecode(s->ss[0], printsh);
01510 if (printsh.ss == m->dummysub) {
01511 printf(" [0] = No subsegment\n");
01512 } else {
01513 printf(" [0] = x%lx %d\n", (unsigned long) printsh.ss,
01514 printsh.ssorient);
01515 }
01516 sdecode(s->ss[1], printsh);
01517 if (printsh.ss == m->dummysub) {
01518 printf(" [1] = No subsegment\n");
01519 } else {
01520 printf(" [1] = x%lx %d\n", (unsigned long) printsh.ss,
01521 printsh.ssorient);
01522 }
01523
01524 sorg(*s, printvertex);
01525 if (printvertex == (vertex) NULL)
01526 printf(" Origin[%d] = NULL\n", 2 + s->ssorient);
01527 else
01528 printf(" Origin[%d] = x%lx (%.12g, %.12g)\n",
01529 2 + s->ssorient, (unsigned long) printvertex,
01530 printvertex[0], printvertex[1]);
01531 sdest(*s, printvertex);
01532 if (printvertex == (vertex) NULL)
01533 printf(" Dest [%d] = NULL\n", 3 - s->ssorient);
01534 else
01535 printf(" Dest [%d] = x%lx (%.12g, %.12g)\n",
01536 3 - s->ssorient, (unsigned long) printvertex,
01537 printvertex[0], printvertex[1]);
01538
01539 decode(s->ss[6], printtri);
01540 if (printtri.tri == m->dummytri) {
01541 printf(" [6] = Outer space\n");
01542 } else {
01543 printf(" [6] = x%lx %d\n", (unsigned long) printtri.tri,
01544 printtri.orient);
01545 }
01546 decode(s->ss[7], printtri);
01547 if (printtri.tri == m->dummytri) {
01548 printf(" [7] = Outer space\n");
01549 } else {
01550 printf(" [7] = x%lx %d\n", (unsigned long) printtri.tri,
01551 printtri.orient);
01552 }
01553
01554 segorg(*s, printvertex);
01555 if (printvertex == (vertex) NULL)
01556 printf(" Segment origin[%d] = NULL\n", 4 + s->ssorient);
01557 else
01558 printf(" Segment origin[%d] = x%lx (%.12g, %.12g)\n",
01559 4 + s->ssorient, (unsigned long) printvertex,
01560 printvertex[0], printvertex[1]);
01561 segdest(*s, printvertex);
01562 if (printvertex == (vertex) NULL)
01563 printf(" Segment dest [%d] = NULL\n", 5 - s->ssorient);
01564 else
01565 printf(" Segment dest [%d] = x%lx (%.12g, %.12g)\n",
01566 5 - s->ssorient, (unsigned long) printvertex,
01567 printvertex[0], printvertex[1]);
01568 }
01569
01572
01573
01574
01578
01579
01580
01581
01582
01583
01584
01585
01586
01587 void poolzero(struct memorypool *pool)
01588 {
01589 pool->firstblock = (int **) NULL;
01590 pool->nowblock = (int **) NULL;
01591 pool->nextitem = (int *) NULL;
01592 pool->deaditemstack = (int *) NULL;
01593 pool->pathblock = (int **) NULL;
01594 pool->pathitem = (int *) NULL;
01595 pool->alignbytes = 0;
01596 pool->itembytes = 0;
01597 pool->itemsperblock = 0;
01598 pool->itemsfirstblock = 0;
01599 pool->items = 0;
01600 pool->maxitems = 0;
01601 pool->unallocateditems = 0;
01602 pool->pathitemsleft = 0;
01603 }
01604
01605
01606
01607
01608
01609
01610
01611
01612
01613
01614
01615 void poolrestart(struct memorypool *pool)
01616 {
01617 unsigned long alignptr;
01618
01619 pool->items = 0;
01620 pool->maxitems = 0;
01621
01622
01623 pool->nowblock = pool->firstblock;
01624
01625 alignptr = (unsigned long) (pool->nowblock + 1);
01626
01627 pool->nextitem = (int *)
01628 (alignptr + (unsigned long) pool->alignbytes -
01629 (alignptr % (unsigned long) pool->alignbytes));
01630
01631 pool->unallocateditems = pool->itemsfirstblock;
01632
01633 pool->deaditemstack = (int *) NULL;
01634 }
01635
01636
01637
01638
01639
01640
01641
01642
01643
01644
01645
01646
01647
01648
01649
01650
01651
01652
01653
01654
01655 void poolinit(struct memorypool *pool, int bytecount, int itemcount,
01656 int firstitemcount, int alignment)
01657 {
01658
01659
01660
01661
01662 if (alignment > (int)sizeof(int *)) {
01663 pool->alignbytes = alignment;
01664 } else {
01665 pool->alignbytes = sizeof(int *);
01666 }
01667 pool->itembytes = ((bytecount - 1) / pool->alignbytes + 1) *
01668 pool->alignbytes;
01669 pool->itemsperblock = itemcount;
01670 if (firstitemcount == 0) {
01671 pool->itemsfirstblock = itemcount;
01672 } else {
01673 pool->itemsfirstblock = firstitemcount;
01674 }
01675
01676
01677
01678
01679 pool->firstblock = (int **)
01680 trimalloc(pool->itemsfirstblock * pool->itembytes + (int) sizeof(int *) +
01681 pool->alignbytes);
01682
01683 *(pool->firstblock) = (int *) NULL;
01684 poolrestart(pool);
01685 }
01686
01687
01688
01689
01690
01691
01692
01693 void pooldeinit(struct memorypool *pool)
01694 {
01695 while (pool->firstblock != (int **) NULL) {
01696 pool->nowblock = (int **) *(pool->firstblock);
01697 trifree((int *) pool->firstblock);
01698 pool->firstblock = pool->nowblock;
01699 }
01700 }
01701
01702
01703
01704
01705
01706
01707
01708 int *poolalloc(struct memorypool *pool)
01709 {
01710 int *newitem;
01711 int **newblock;
01712 unsigned long alignptr;
01713
01714
01715
01716 if (pool->deaditemstack != (int *) NULL) {
01717 newitem = pool->deaditemstack;
01718 pool->deaditemstack = * (int **) pool->deaditemstack;
01719 } else {
01720
01721 if (pool->unallocateditems == 0) {
01722
01723 if (*(pool->nowblock) == (int *) NULL) {
01724
01725 newblock = (int **) trimalloc(pool->itemsperblock * pool->itembytes +
01726 (int) sizeof(int *) +
01727 pool->alignbytes);
01728 *(pool->nowblock) = (int *) newblock;
01729
01730 *newblock = (int *) NULL;
01731 }
01732
01733
01734 pool->nowblock = (int **) *(pool->nowblock);
01735
01736
01737 alignptr = (unsigned long) (pool->nowblock + 1);
01738
01739 pool->nextitem = (int *)
01740 (alignptr + (unsigned long) pool->alignbytes -
01741 (alignptr % (unsigned long) pool->alignbytes));
01742
01743 pool->unallocateditems = pool->itemsperblock;
01744 }
01745
01746
01747 newitem = pool->nextitem;
01748
01749 pool->nextitem = (int *) ((char *) pool->nextitem + pool->itembytes);
01750 pool->unallocateditems--;
01751 pool->maxitems++;
01752 }
01753 pool->items++;
01754 return newitem;
01755 }
01756
01757
01758
01759
01760
01761
01762
01763
01764
01765 void pooldealloc(struct memorypool *pool, int *dyingitem)
01766 {
01767
01768 *((int **) dyingitem) = pool->deaditemstack;
01769 pool->deaditemstack = dyingitem;
01770 pool->items--;
01771 }
01772
01773
01774
01775
01776
01777
01778
01779
01780
01781 void traversalinit(struct memorypool *pool)
01782 {
01783 unsigned long alignptr;
01784
01785
01786 pool->pathblock = pool->firstblock;
01787
01788 alignptr = (unsigned long) (pool->pathblock + 1);
01789
01790 pool->pathitem = (int *)
01791 (alignptr + (unsigned long) pool->alignbytes -
01792 (alignptr % (unsigned long) pool->alignbytes));
01793
01794 pool->pathitemsleft = pool->itemsfirstblock;
01795 }
01796
01797
01798
01799
01800
01801
01802
01803
01804
01805
01806
01807
01808
01809
01810
01811 int *traverse(struct memorypool *pool)
01812 {
01813 int *newitem;
01814 unsigned long alignptr;
01815
01816
01817 if (pool->pathitem == pool->nextitem) {
01818 return (int *) NULL;
01819 }
01820
01821
01822 if (pool->pathitemsleft == 0) {
01823
01824 pool->pathblock = (int **) *(pool->pathblock);
01825
01826 alignptr = (unsigned long) (pool->pathblock + 1);
01827
01828 pool->pathitem = (int *)
01829 (alignptr + (unsigned long) pool->alignbytes -
01830 (alignptr % (unsigned long) pool->alignbytes));
01831
01832 pool->pathitemsleft = pool->itemsperblock;
01833 }
01834
01835 newitem = pool->pathitem;
01836
01837 pool->pathitem = (int *) ((char *) pool->pathitem + pool->itembytes);
01838 pool->pathitemsleft--;
01839 return newitem;
01840 }
01841
01842
01843
01844
01845
01846
01847
01848
01849
01850
01851
01852
01853
01854
01855
01856
01857
01858
01859
01860
01861
01862
01863
01864
01865
01866
01867
01868
01869
01870 void dummyinit(struct mesh *m, struct behavior *b, int trianglebytes,
01871 int subsegbytes)
01872 {
01873 unsigned long alignptr;
01874
01875
01876 m->dummytribase = (triangle *) trimalloc(trianglebytes +
01877 m->triangles.alignbytes);
01878
01879 alignptr = (unsigned long) m->dummytribase;
01880 m->dummytri = (triangle *)
01881 (alignptr + (unsigned long) m->triangles.alignbytes -
01882 (alignptr % (unsigned long) m->triangles.alignbytes));
01883
01884
01885
01886
01887 m->dummytri[0] = (triangle) m->dummytri;
01888 m->dummytri[1] = (triangle) m->dummytri;
01889 m->dummytri[2] = (triangle) m->dummytri;
01890
01891 m->dummytri[3] = (triangle) NULL;
01892 m->dummytri[4] = (triangle) NULL;
01893 m->dummytri[5] = (triangle) NULL;
01894
01895 if (b->usesegments) {
01896
01897
01898
01899 m->dummysubbase = (subseg *) trimalloc(subsegbytes +
01900 m->subsegs.alignbytes);
01901
01902 alignptr = (unsigned long) m->dummysubbase;
01903 m->dummysub = (subseg *)
01904 (alignptr + (unsigned long) m->subsegs.alignbytes -
01905 (alignptr % (unsigned long) m->subsegs.alignbytes));
01906
01907
01908
01909
01910 m->dummysub[0] = (subseg) m->dummysub;
01911 m->dummysub[1] = (subseg) m->dummysub;
01912
01913 m->dummysub[2] = (subseg) NULL;
01914 m->dummysub[3] = (subseg) NULL;
01915 m->dummysub[4] = (subseg) NULL;
01916 m->dummysub[5] = (subseg) NULL;
01917
01918 m->dummysub[6] = (subseg) m->dummytri;
01919 m->dummysub[7] = (subseg) m->dummytri;
01920
01921 * (int *) (m->dummysub + 8) = 0;
01922
01923
01924
01925 m->dummytri[6] = (triangle) m->dummysub;
01926 m->dummytri[7] = (triangle) m->dummysub;
01927 m->dummytri[8] = (triangle) m->dummysub;
01928 }
01929 }
01930
01931
01932
01933
01934
01935
01936
01937
01938
01939
01940
01941 void initializevertexpool(struct mesh *m, struct behavior *b)
01942 {
01943 int vertexsize;
01944
01945
01946
01947
01948 m->vertexmarkindex = ((m->mesh_dim + m->nextras) * sizeof(float) +
01949 sizeof(int) - 1) /
01950 sizeof(int);
01951 vertexsize = (m->vertexmarkindex + 2) * sizeof(int);
01952 if (b->poly) {
01953
01954
01955 m->vertex2triindex = (vertexsize + sizeof(triangle) - 1) /
01956 sizeof(triangle);
01957 vertexsize = (m->vertex2triindex + 1) * sizeof(triangle);
01958 }
01959
01960
01961 poolinit(&m->vertices, vertexsize, VERTEXPERBLOCK,
01962 m->invertices > VERTEXPERBLOCK ? m->invertices : VERTEXPERBLOCK,
01963 sizeof(float));
01964 }
01965
01966
01967
01968
01969
01970
01971
01972
01973
01974
01975
01976
01977 void initializetrisubpools(struct mesh *m, struct behavior *b)
01978 {
01979 int trisize;
01980
01981
01982
01983
01984
01985 m->highorderindex = 6 + (b->usesegments * 3);
01986
01987 trisize = ((b->order + 1) * (b->order + 2) / 2 + (m->highorderindex - 3)) *
01988 sizeof(triangle);
01989
01990
01991 m->elemattribindex = (trisize + sizeof(float) - 1) / sizeof(float);
01992
01993
01994
01995 m->areaboundindex = m->elemattribindex + m->eextras + b->regionattrib;
01996
01997
01998 if (b->vararea) {
01999 trisize = (m->areaboundindex + 1) * sizeof(float);
02000 } else if (m->eextras + b->regionattrib > 0) {
02001 trisize = m->areaboundindex * sizeof(float);
02002 }
02003
02004
02005
02006
02007 if ((b->voronoi || b->neighbors) &&
02008 (trisize < 6 * (int)sizeof(triangle) + (int)sizeof(int))) {
02009 trisize = 6 * sizeof(triangle) + sizeof(int);
02010 }
02011
02012
02013 poolinit(&m->triangles, trisize, TRIPERBLOCK,
02014 (2 * m->invertices - 2) > TRIPERBLOCK ? (2 * m->invertices - 2) :
02015 TRIPERBLOCK, 4);
02016
02017 if (b->usesegments) {
02018
02019
02020 poolinit(&m->subsegs, 8 * sizeof(triangle) + sizeof(int),
02021 SUBSEGPERBLOCK, SUBSEGPERBLOCK, 4);
02022
02023
02024 dummyinit(m, b, m->triangles.itembytes, m->subsegs.itembytes);
02025 } else {
02026
02027 dummyinit(m, b, m->triangles.itembytes, 0);
02028 }
02029 }
02030
02031
02032
02033
02034
02035
02036
02037 void triangledealloc(struct mesh *m, triangle *dyingtriangle)
02038 {
02039
02040
02041 killtri(dyingtriangle);
02042 pooldealloc(&m->triangles, (int *) dyingtriangle);
02043 }
02044
02045
02046
02047
02048
02049
02050
02051 triangle *triangletraverse(struct mesh *m)
02052 {
02053 triangle *newtriangle;
02054
02055 do {
02056 newtriangle = (triangle *) traverse(&m->triangles);
02057 if (newtriangle == (triangle *) NULL) {
02058 return (triangle *) NULL;
02059 }
02060 } while (deadtri(newtriangle));
02061 return newtriangle;
02062 }
02063
02064
02065
02066
02067
02068
02069
02070 void subsegdealloc(struct mesh *m, subseg *dyingsubseg)
02071 {
02072
02073
02074 killsubseg(dyingsubseg);
02075 pooldealloc(&m->subsegs, (int *) dyingsubseg);
02076 }
02077
02078
02079
02080
02081
02082
02083
02084 subseg *subsegtraverse(struct mesh *m)
02085 {
02086 subseg *newsubseg;
02087
02088 do {
02089 newsubseg = (subseg *) traverse(&m->subsegs);
02090 if (newsubseg == (subseg *) NULL) {
02091 return (subseg *) NULL;
02092 }
02093 } while (deadsubseg(newsubseg));
02094 return newsubseg;
02095 }
02096
02097
02098
02099
02100
02101
02102
02103 void vertexdealloc(struct mesh *m, vertex dyingvertex)
02104 {
02105
02106
02107 setvertextype(dyingvertex, DEADVERTEX);
02108 pooldealloc(&m->vertices, (int *) dyingvertex);
02109 }
02110
02111
02112
02113
02114
02115
02116
02117 vertex vertextraverse(struct mesh *m)
02118 {
02119 vertex newvertex;
02120
02121 do {
02122 newvertex = (vertex) traverse(&m->vertices);
02123 if (newvertex == (vertex) NULL) {
02124 return (vertex) NULL;
02125 }
02126 } while (vertextype(newvertex) == DEADVERTEX);
02127 return newvertex;
02128 }
02129
02130
02131
02132
02133
02134
02135
02136
02137
02138
02139
02140
02141
02142 vertex getvertex(struct mesh *m, struct behavior *b, int number)
02143 {
02144 int **getblock;
02145 char *foundvertex;
02146 unsigned long alignptr;
02147 int current;
02148
02149 getblock = m->vertices.firstblock;
02150 current = b->firstnumber;
02151
02152
02153 if (current + m->vertices.itemsfirstblock <= number) {
02154 getblock = (int **) *getblock;
02155 current += m->vertices.itemsfirstblock;
02156 while (current + m->vertices.itemsperblock <= number) {
02157 getblock = (int **) *getblock;
02158 current += m->vertices.itemsperblock;
02159 }
02160 }
02161
02162
02163 alignptr = (unsigned long) (getblock + 1);
02164 foundvertex = (char *) (alignptr + (unsigned long) m->vertices.alignbytes -
02165 (alignptr % (unsigned long) m->vertices.alignbytes));
02166 return (vertex) (foundvertex + m->vertices.itembytes * (number - current));
02167 }
02168
02169
02170
02171
02172
02173
02174
02175 void triangledeinit(struct mesh *m, struct behavior *b)
02176 {
02177 pooldeinit(&m->triangles);
02178 trifree((int *) m->dummytribase);
02179 if (b->usesegments) {
02180 pooldeinit(&m->subsegs);
02181 trifree((int *) m->dummysubbase);
02182 }
02183 pooldeinit(&m->vertices);
02184 }
02185
02188
02189
02190
02194
02195
02196
02197
02198
02199
02200 void maketriangle(struct mesh *m, struct behavior *b, struct otri *newotri)
02201 {
02202 int i;
02203
02204 newotri->tri = (triangle *) poolalloc(&m->triangles);
02205
02206 newotri->tri[0] = (triangle) m->dummytri;
02207 newotri->tri[1] = (triangle) m->dummytri;
02208 newotri->tri[2] = (triangle) m->dummytri;
02209
02210 newotri->tri[3] = (triangle) NULL;
02211 newotri->tri[4] = (triangle) NULL;
02212 newotri->tri[5] = (triangle) NULL;
02213 if (b->usesegments) {
02214
02215
02216 newotri->tri[6] = (triangle) m->dummysub;
02217 newotri->tri[7] = (triangle) m->dummysub;
02218 newotri->tri[8] = (triangle) m->dummysub;
02219 }
02220 for (i = 0; i < m->eextras; i++) {
02221 setelemattribute(*newotri, i, 0.0);
02222 }
02223 if (b->vararea) {
02224 setareabound(*newotri, -1.0);
02225 }
02226
02227 newotri->orient = 0;
02228 }
02229
02230
02231
02232
02233
02234
02235
02236 void makesubseg(struct mesh *m, struct osub *newsubseg)
02237 {
02238 newsubseg->ss = (subseg *) poolalloc(&m->subsegs);
02239
02240
02241 newsubseg->ss[0] = (subseg) m->dummysub;
02242 newsubseg->ss[1] = (subseg) m->dummysub;
02243
02244 newsubseg->ss[2] = (subseg) NULL;
02245 newsubseg->ss[3] = (subseg) NULL;
02246 newsubseg->ss[4] = (subseg) NULL;
02247 newsubseg->ss[5] = (subseg) NULL;
02248
02249 newsubseg->ss[6] = (subseg) m->dummytri;
02250 newsubseg->ss[7] = (subseg) m->dummytri;
02251
02252 setmark(*newsubseg, 0);
02253
02254 newsubseg->ssorient = 0;
02255 }
02256
02259
02260
02261
02265
02266
02267
02268
02269
02270
02271
02272
02273
02274
02275
02276
02277
02278 #define Absolute(a) ((a) >= 0.0 ? (a) : -(a))
02279
02280
02281
02282
02283
02284
02285
02286
02287
02288
02289
02290
02291
02292
02293
02294 #define Fast_Two_Sum_Tail(a, b, x, y) \
02295 bvirt = x - a; \
02296 y = b - bvirt
02297
02298 #define Fast_Two_Sum(a, b, x, y) \
02299 x = (float) (a + b); \
02300 Fast_Two_Sum_Tail(a, b, x, y)
02301
02302 #define Two_Sum_Tail(a, b, x, y) \
02303 bvirt = (float) (x - a); \
02304 avirt = x - bvirt; \
02305 bround = b - bvirt; \
02306 around = a - avirt; \
02307 y = around + bround
02308
02309 #define Two_Sum(a, b, x, y) \
02310 x = (float) (a + b); \
02311 Two_Sum_Tail(a, b, x, y)
02312
02313 #define Two_Diff_Tail(a, b, x, y) \
02314 bvirt = (float) (a - x); \
02315 avirt = x + bvirt; \
02316 bround = bvirt - b; \
02317 around = a - avirt; \
02318 y = around + bround
02319
02320 #define Two_Diff(a, b, x, y) \
02321 x = (float) (a - b); \
02322 Two_Diff_Tail(a, b, x, y)
02323
02324 #define Split(a, ahi, alo) \
02325 c = (float) (splitter * a); \
02326 abig = (float) (c - a); \
02327 ahi = c - abig; \
02328 alo = a - ahi
02329
02330 #define Two_Product_Tail(a, b, x, y) \
02331 Split(a, ahi, alo); \
02332 Split(b, bhi, blo); \
02333 err1 = x - (ahi * bhi); \
02334 err2 = err1 - (alo * bhi); \
02335 err3 = err2 - (ahi * blo); \
02336 y = (alo * blo) - err3
02337
02338 #define Two_Product(a, b, x, y) \
02339 x = (float) (a * b); \
02340 Two_Product_Tail(a, b, x, y)
02341
02342
02343
02344
02345 #define Two_Product_Presplit(a, b, bhi, blo, x, y) \
02346 x = (float) (a * b); \
02347 Split(a, ahi, alo); \
02348 err1 = x - (ahi * bhi); \
02349 err2 = err1 - (alo * bhi); \
02350 err3 = err2 - (ahi * blo); \
02351 y = (alo * blo) - err3
02352
02353
02354
02355 #define Square_Tail(a, x, y) \
02356 Split(a, ahi, alo); \
02357 err1 = x - (ahi * ahi); \
02358 err3 = err1 - ((ahi + ahi) * alo); \
02359 y = (alo * alo) - err3
02360
02361 #define Square(a, x, y) \
02362 x = (float) (a * a); \
02363 Square_Tail(a, x, y)
02364
02365
02366
02367
02368 #define Two_One_Sum(a1, a0, b, x2, x1, x0) \
02369 Two_Sum(a0, b , _i, x0); \
02370 Two_Sum(a1, _i, x2, x1)
02371
02372 #define Two_One_Diff(a1, a0, b, x2, x1, x0) \
02373 Two_Diff(a0, b , _i, x0); \
02374 Two_Sum( a1, _i, x2, x1)
02375
02376 #define Two_Two_Sum(a1, a0, b1, b0, x3, x2, x1, x0) \
02377 Two_One_Sum(a1, a0, b0, _j, _0, x0); \
02378 Two_One_Sum(_j, _0, b1, x3, x2, x1)
02379
02380 #define Two_Two_Diff(a1, a0, b1, b0, x3, x2, x1, x0) \
02381 Two_One_Diff(a1, a0, b0, _j, _0, x0); \
02382 Two_One_Diff(_j, _0, b1, x3, x2, x1)
02383
02384
02385
02386 #define Two_One_Product(a1, a0, b, x3, x2, x1, x0) \
02387 Split(b, bhi, blo); \
02388 Two_Product_Presplit(a0, b, bhi, blo, _i, x0); \
02389 Two_Product_Presplit(a1, b, bhi, blo, _j, _0); \
02390 Two_Sum(_i, _0, _k, x1); \
02391 Fast_Two_Sum(_j, _k, x3, x2)
02392
02393
02394
02395
02396
02397
02398
02399
02400
02401
02402
02403
02404
02405
02406
02407
02408
02409
02410
02411
02412 void exactinit()
02413 {
02414 float half;
02415 float check, lastcheck;
02416 int every_other;
02417 every_other = 1;
02418 half = 0.5;
02419 epsilon = 1.0;
02420 splitter = 1.0;
02421 check = 1.0;
02422
02423
02424
02425
02426 do {
02427 lastcheck = check;
02428 epsilon *= half;
02429 if (every_other) {
02430 splitter *= 2.0;
02431 }
02432 every_other = !every_other;
02433 check = 1.0 + epsilon;
02434 } while ((check != 1.0) && (check != lastcheck));
02435 splitter += 1.0;
02436
02437 resulterrbound = (3.0 + 8.0 * epsilon) * epsilon;
02438 ccwerrboundA = (3.0 + 16.0 * epsilon) * epsilon;
02439 ccwerrboundB = (2.0 + 12.0 * epsilon) * epsilon;
02440 ccwerrboundC = (9.0 + 64.0 * epsilon) * epsilon * epsilon;
02441 iccerrboundA = (10.0 + 96.0 * epsilon) * epsilon;
02442 iccerrboundB = (4.0 + 48.0 * epsilon) * epsilon;
02443 iccerrboundC = (44.0 + 576.0 * epsilon) * epsilon * epsilon;
02444 o3derrboundA = (7.0 + 56.0 * epsilon) * epsilon;
02445 o3derrboundB = (3.0 + 28.0 * epsilon) * epsilon;
02446 o3derrboundC = (26.0 + 288.0 * epsilon) * epsilon * epsilon;
02447 }
02448
02449
02450
02451
02452
02453
02454
02455
02456
02457
02458
02459
02460
02461
02462
02463 int fast_expansion_sum_zeroelim(int elen, float *e, int flen, float *f, float *h)
02464 {
02465 float Q;
02466 float Qnew;
02467 float hh;
02468 float bvirt;
02469 float avirt, bround, around;
02470 int eindex, findex, hindex;
02471 float enow, fnow;
02472
02473 enow = e[0];
02474 fnow = f[0];
02475 eindex = findex = 0;
02476 if ((fnow > enow) == (fnow > -enow)) {
02477 Q = enow;
02478 enow = e[++eindex];
02479 } else {
02480 Q = fnow;
02481 fnow = f[++findex];
02482 }
02483 hindex = 0;
02484 if ((eindex < elen) && (findex < flen)) {
02485 if ((fnow > enow) == (fnow > -enow)) {
02486 Fast_Two_Sum(enow, Q, Qnew, hh);
02487 enow = e[++eindex];
02488 } else {
02489 Fast_Two_Sum(fnow, Q, Qnew, hh);
02490 fnow = f[++findex];
02491 }
02492 Q = Qnew;
02493 if (hh != 0.0) {
02494 h[hindex++] = hh;
02495 }
02496 while ((eindex < elen) && (findex < flen)) {
02497 if ((fnow > enow) == (fnow > -enow)) {
02498 Two_Sum(Q, enow, Qnew, hh);
02499 enow = e[++eindex];
02500 } else {
02501 Two_Sum(Q, fnow, Qnew, hh);
02502 fnow = f[++findex];
02503 }
02504 Q = Qnew;
02505 if (hh != 0.0) {
02506 h[hindex++] = hh;
02507 }
02508 }
02509 }
02510 while (eindex < elen) {
02511 Two_Sum(Q, enow, Qnew, hh);
02512 enow = e[++eindex];
02513 Q = Qnew;
02514 if (hh != 0.0) {
02515 h[hindex++] = hh;
02516 }
02517 }
02518 while (findex < flen) {
02519 Two_Sum(Q, fnow, Qnew, hh);
02520 fnow = f[++findex];
02521 Q = Qnew;
02522 if (hh != 0.0) {
02523 h[hindex++] = hh;
02524 }
02525 }
02526 if ((Q != 0.0) || (hindex == 0)) {
02527 h[hindex++] = Q;
02528 }
02529 return hindex;
02530 }
02531
02532
02533
02534
02535
02536
02537
02538
02539
02540
02541
02542
02543
02544
02545
02546
02547 int scale_expansion_zeroelim(int elen, float *e, float b, float *h)
02548 {
02549 float Q, sum;
02550 float hh;
02551 float product1;
02552 float product0;
02553 int eindex, hindex;
02554 float enow;
02555 float bvirt;
02556 float avirt, bround, around;
02557 float c;
02558 float abig;
02559 float ahi, alo, bhi, blo;
02560 float err1, err2, err3;
02561
02562 Split(b, bhi, blo);
02563 Two_Product_Presplit(e[0], b, bhi, blo, Q, hh);
02564 hindex = 0;
02565 if (hh != 0) {
02566 h[hindex++] = hh;
02567 }
02568 for (eindex = 1; eindex < elen; eindex++) {
02569 enow = e[eindex];
02570 Two_Product_Presplit(enow, b, bhi, blo, product1, product0);
02571 Two_Sum(Q, product0, sum, hh);
02572 if (hh != 0) {
02573 h[hindex++] = hh;
02574 }
02575 Fast_Two_Sum(product1, sum, Q, hh);
02576 if (hh != 0) {
02577 h[hindex++] = hh;
02578 }
02579 }
02580 if ((Q != 0.0) || (hindex == 0)) {
02581 h[hindex++] = Q;
02582 }
02583 return hindex;
02584 }
02585
02586
02587
02588
02589
02590
02591
02592
02593
02594 float estimate(int elen, float *e)
02595 {
02596 float Q;
02597 int eindex;
02598 Q = e[0];
02599 for (eindex = 1; eindex < elen; eindex++) {
02600 Q += e[eindex];
02601 }
02602 return Q;
02603 }
02604
02605
02606
02607
02608
02609
02610
02611
02612
02613
02614
02615
02616
02617
02618
02619
02620
02621
02622
02623
02624
02625 float counterclockwiseadapt(vertex pa, vertex pb, vertex pc, float detsum)
02626 {
02627 float acx, acy, bcx, bcy;
02628 float acxtail, acytail, bcxtail, bcytail;
02629 float detleft, detright;
02630 float detlefttail, detrighttail;
02631 float det, errbound;
02632 float B[4], C1[8], C2[12], D[16];
02633 float B3;
02634 int C1length, C2length, Dlength;
02635 float u[4];
02636 float u3;
02637 float s1, t1;
02638 float s0, t0;
02639
02640 float bvirt;
02641 float avirt, bround, around;
02642 float c;
02643 float abig;
02644 float ahi, alo, bhi, blo;
02645 float err1, err2, err3;
02646 float _i, _j;
02647 float _0;
02648
02649 acx = (float) (pa[0] - pc[0]);
02650 bcx = (float) (pb[0] - pc[0]);
02651 acy = (float) (pa[1] - pc[1]);
02652 bcy = (float) (pb[1] - pc[1]);
02653
02654 Two_Product(acx, bcy, detleft, detlefttail);
02655 Two_Product(acy, bcx, detright, detrighttail);
02656
02657 Two_Two_Diff(detleft, detlefttail, detright, detrighttail,
02658 B3, B[2], B[1], B[0]);
02659 B[3] = B3;
02660
02661 det = estimate(4, B);
02662 errbound = ccwerrboundB * detsum;
02663 if ((det >= errbound) || (-det >= errbound)) {
02664 return det;
02665 }
02666
02667 Two_Diff_Tail(pa[0], pc[0], acx, acxtail);
02668 Two_Diff_Tail(pb[0], pc[0], bcx, bcxtail);
02669 Two_Diff_Tail(pa[1], pc[1], acy, acytail);
02670 Two_Diff_Tail(pb[1], pc[1], bcy, bcytail);
02671
02672 if ((acxtail == 0.0) && (acytail == 0.0)
02673 && (bcxtail == 0.0) && (bcytail == 0.0)) {
02674 return det;
02675 }
02676
02677 errbound = ccwerrboundC * detsum + resulterrbound * Absolute(det);
02678 det += (acx * bcytail + bcy * acxtail)
02679 - (acy * bcxtail + bcx * acytail);
02680 if ((det >= errbound) || (-det >= errbound)) {
02681 return det;
02682 }
02683
02684 Two_Product(acxtail, bcy, s1, s0);
02685 Two_Product(acytail, bcx, t1, t0);
02686 Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]);
02687 u[3] = u3;
02688 C1length = fast_expansion_sum_zeroelim(4, B, 4, u, C1);
02689
02690 Two_Product(acx, bcytail, s1, s0);
02691 Two_Product(acy, bcxtail, t1, t0);
02692 Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]);
02693 u[3] = u3;
02694 C2length = fast_expansion_sum_zeroelim(C1length, C1, 4, u, C2);
02695
02696 Two_Product(acxtail, bcytail, s1, s0);
02697 Two_Product(acytail, bcxtail, t1, t0);
02698 Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]);
02699 u[3] = u3;
02700 Dlength = fast_expansion_sum_zeroelim(C2length, C2, 4, u, D);
02701
02702 return(D[Dlength - 1]);
02703 }
02704
02705 float counterclockwise(struct mesh *m, struct behavior *b,
02706 vertex pa, vertex pb, vertex pc)
02707 {
02708 float detleft, detright, det;
02709 float detsum, errbound;
02710
02711 m->counterclockcount++;
02712
02713 detleft = (pa[0] - pc[0]) * (pb[1] - pc[1]);
02714 detright = (pa[1] - pc[1]) * (pb[0] - pc[0]);
02715 det = detleft - detright;
02716
02717 if (b->noexact) {
02718 return det;
02719 }
02720
02721 if (detleft > 0.0) {
02722 if (detright <= 0.0) {
02723 return det;
02724 } else {
02725 detsum = detleft + detright;
02726 }
02727 } else if (detleft < 0.0) {
02728 if (detright >= 0.0) {
02729 return det;
02730 } else {
02731 detsum = -detleft - detright;
02732 }
02733 } else {
02734 return det;
02735 }
02736
02737 errbound = ccwerrboundA * detsum;
02738 if ((det >= errbound) || (-det >= errbound)) {
02739 return det;
02740 }
02741
02742 return counterclockwiseadapt(pa, pb, pc, detsum);
02743 }
02744
02745
02746
02747
02748
02749
02750
02751
02752
02753
02754
02755
02756
02757
02758
02759
02760
02761
02762
02763
02764 float incircleadapt(vertex pa, vertex pb, vertex pc, vertex pd, float permanent)
02765 {
02766 float adx, bdx, cdx, ady, bdy, cdy;
02767 float det, errbound;
02768
02769 float bdxcdy1, cdxbdy1, cdxady1, adxcdy1, adxbdy1, bdxady1;
02770 float bdxcdy0, cdxbdy0, cdxady0, adxcdy0, adxbdy0, bdxady0;
02771 float bc[4], ca[4], ab[4];
02772 float bc3, ca3, ab3;
02773 float axbc[8], axxbc[16], aybc[8], ayybc[16], adet[32];
02774 int axbclen, axxbclen, aybclen, ayybclen, alen;
02775 float bxca[8], bxxca[16], byca[8], byyca[16], bdet[32];
02776 int bxcalen, bxxcalen, bycalen, byycalen, blen;
02777 float cxab[8], cxxab[16], cyab[8], cyyab[16], cdet[32];
02778 int cxablen, cxxablen, cyablen, cyyablen, clen;
02779 float abdet[64];
02780 int ablen;
02781 float fin1[1152], fin2[1152];
02782 float *finnow, *finother, *finswap;
02783 int finlength;
02784
02785 float adxtail, bdxtail, cdxtail, adytail, bdytail, cdytail;
02786 float adxadx1, adyady1, bdxbdx1, bdybdy1, cdxcdx1, cdycdy1;
02787 float adxadx0, adyady0, bdxbdx0, bdybdy0, cdxcdx0, cdycdy0;
02788 float aa[4], bb[4], cc[4];
02789 float aa3, bb3, cc3;
02790 float ti1, tj1;
02791 float ti0, tj0;
02792 float u[4], v[4];
02793 float u3, v3;
02794 float temp8[8], temp16a[16], temp16b[16], temp16c[16];
02795 float temp32a[32], temp32b[32], temp48[48], temp64[64];
02796 int temp8len, temp16alen, temp16blen, temp16clen;
02797 int temp32alen, temp32blen, temp48len, temp64len;
02798 float axtbb[8], axtcc[8], aytbb[8], aytcc[8];
02799 int axtbblen, axtcclen, aytbblen, aytcclen;
02800 float bxtaa[8], bxtcc[8], bytaa[8], bytcc[8];
02801 int bxtaalen, bxtcclen, bytaalen, bytcclen;
02802 float cxtaa[8], cxtbb[8], cytaa[8], cytbb[8];
02803 int cxtaalen, cxtbblen, cytaalen, cytbblen;
02804 float axtbc[8], aytbc[8], bxtca[8], bytca[8], cxtab[8], cytab[8];
02805 int axtbclen, aytbclen, bxtcalen, bytcalen, cxtablen, cytablen;
02806 float axtbct[16], aytbct[16], bxtcat[16], bytcat[16], cxtabt[16], cytabt[16];
02807 int axtbctlen, aytbctlen, bxtcatlen, bytcatlen, cxtabtlen, cytabtlen;
02808 float axtbctt[8], aytbctt[8], bxtcatt[8];
02809 float bytcatt[8], cxtabtt[8], cytabtt[8];
02810 int axtbcttlen, aytbcttlen, bxtcattlen, bytcattlen, cxtabttlen, cytabttlen;
02811 float abt[8], bct[8], cat[8];
02812 int abtlen, bctlen, catlen;
02813 float abtt[4], bctt[4], catt[4];
02814 int abttlen, bcttlen, cattlen;
02815 float abtt3, bctt3, catt3;
02816 float negate;
02817
02818 float bvirt;
02819 float avirt, bround, around;
02820 float c;
02821 float abig;
02822 float ahi, alo, bhi, blo;
02823 float err1, err2, err3;
02824 float _i, _j;
02825 float _0;
02826
02827 adx = (float) (pa[0] - pd[0]);
02828 bdx = (float) (pb[0] - pd[0]);
02829 cdx = (float) (pc[0] - pd[0]);
02830 ady = (float) (pa[1] - pd[1]);
02831 bdy = (float) (pb[1] - pd[1]);
02832 cdy = (float) (pc[1] - pd[1]);
02833
02834 Two_Product(bdx, cdy, bdxcdy1, bdxcdy0);
02835 Two_Product(cdx, bdy, cdxbdy1, cdxbdy0);
02836 Two_Two_Diff(bdxcdy1, bdxcdy0, cdxbdy1, cdxbdy0, bc3, bc[2], bc[1], bc[0]);
02837 bc[3] = bc3;
02838 axbclen = scale_expansion_zeroelim(4, bc, adx, axbc);
02839 axxbclen = scale_expansion_zeroelim(axbclen, axbc, adx, axxbc);
02840 aybclen = scale_expansion_zeroelim(4, bc, ady, aybc);
02841 ayybclen = scale_expansion_zeroelim(aybclen, aybc, ady, ayybc);
02842 alen = fast_expansion_sum_zeroelim(axxbclen, axxbc, ayybclen, ayybc, adet);
02843
02844 Two_Product(cdx, ady, cdxady1, cdxady0);
02845 Two_Product(adx, cdy, adxcdy1, adxcdy0);
02846 Two_Two_Diff(cdxady1, cdxady0, adxcdy1, adxcdy0, ca3, ca[2], ca[1], ca[0]);
02847 ca[3] = ca3;
02848 bxcalen = scale_expansion_zeroelim(4, ca, bdx, bxca);
02849 bxxcalen = scale_expansion_zeroelim(bxcalen, bxca, bdx, bxxca);
02850 bycalen = scale_expansion_zeroelim(4, ca, bdy, byca);
02851 byycalen = scale_expansion_zeroelim(bycalen, byca, bdy, byyca);
02852 blen = fast_expansion_sum_zeroelim(bxxcalen, bxxca, byycalen, byyca, bdet);
02853
02854 Two_Product(adx, bdy, adxbdy1, adxbdy0);
02855 Two_Product(bdx, ady, bdxady1, bdxady0);
02856 Two_Two_Diff(adxbdy1, adxbdy0, bdxady1, bdxady0, ab3, ab[2], ab[1], ab[0]);
02857 ab[3] = ab3;
02858 cxablen = scale_expansion_zeroelim(4, ab, cdx, cxab);
02859 cxxablen = scale_expansion_zeroelim(cxablen, cxab, cdx, cxxab);
02860 cyablen = scale_expansion_zeroelim(4, ab, cdy, cyab);
02861 cyyablen = scale_expansion_zeroelim(cyablen, cyab, cdy, cyyab);
02862 clen = fast_expansion_sum_zeroelim(cxxablen, cxxab, cyyablen, cyyab, cdet);
02863
02864 ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
02865 finlength = fast_expansion_sum_zeroelim(ablen, abdet, clen, cdet, fin1);
02866
02867 det = estimate(finlength, fin1);
02868 errbound = iccerrboundB * permanent;
02869 if ((det >= errbound) || (-det >= errbound)) {
02870 return det;
02871 }
02872
02873 Two_Diff_Tail(pa[0], pd[0], adx, adxtail);
02874 Two_Diff_Tail(pa[1], pd[1], ady, adytail);
02875 Two_Diff_Tail(pb[0], pd[0], bdx, bdxtail);
02876 Two_Diff_Tail(pb[1], pd[1], bdy, bdytail);
02877 Two_Diff_Tail(pc[0], pd[0], cdx, cdxtail);
02878 Two_Diff_Tail(pc[1], pd[1], cdy, cdytail);
02879 if ((adxtail == 0.0) && (bdxtail == 0.0) && (cdxtail == 0.0)
02880 && (adytail == 0.0) && (bdytail == 0.0) && (cdytail == 0.0)) {
02881 return det;
02882 }
02883
02884 errbound = iccerrboundC * permanent + resulterrbound * Absolute(det);
02885 det += ((adx * adx + ady * ady) * ((bdx * cdytail + cdy * bdxtail)
02886 - (bdy * cdxtail + cdx * bdytail))
02887 + 2.0 * (adx * adxtail + ady * adytail) * (bdx * cdy - bdy * cdx))
02888 + ((bdx * bdx + bdy * bdy) * ((cdx * adytail + ady * cdxtail)
02889 - (cdy * adxtail + adx * cdytail))
02890 + 2.0 * (bdx * bdxtail + bdy * bdytail) * (cdx * ady - cdy * adx))
02891 + ((cdx * cdx + cdy * cdy) * ((adx * bdytail + bdy * adxtail)
02892 - (ady * bdxtail + bdx * adytail))
02893 + 2.0 * (cdx * cdxtail + cdy * cdytail) * (adx * bdy - ady * bdx));
02894 if ((det >= errbound) || (-det >= errbound)) {
02895 return det;
02896 }
02897
02898 finnow = fin1;
02899 finother = fin2;
02900
02901 if ((bdxtail != 0.0) || (bdytail != 0.0)
02902 || (cdxtail != 0.0) || (cdytail != 0.0)) {
02903 Square(adx, adxadx1, adxadx0);
02904 Square(ady, adyady1, adyady0);
02905 Two_Two_Sum(adxadx1, adxadx0, adyady1, adyady0, aa3, aa[2], aa[1], aa[0]);
02906 aa[3] = aa3;
02907 }
02908 if ((cdxtail != 0.0) || (cdytail != 0.0)
02909 || (adxtail != 0.0) || (adytail != 0.0)) {
02910 Square(bdx, bdxbdx1, bdxbdx0);
02911 Square(bdy, bdybdy1, bdybdy0);
02912 Two_Two_Sum(bdxbdx1, bdxbdx0, bdybdy1, bdybdy0, bb3, bb[2], bb[1], bb[0]);
02913 bb[3] = bb3;
02914 }
02915 if ((adxtail != 0.0) || (adytail != 0.0)
02916 || (bdxtail != 0.0) || (bdytail != 0.0)) {
02917 Square(cdx, cdxcdx1, cdxcdx0);
02918 Square(cdy, cdycdy1, cdycdy0);
02919 Two_Two_Sum(cdxcdx1, cdxcdx0, cdycdy1, cdycdy0, cc3, cc[2], cc[1], cc[0]);
02920 cc[3] = cc3;
02921 }
02922
02923 if (adxtail != 0.0) {
02924 axtbclen = scale_expansion_zeroelim(4, bc, adxtail, axtbc);
02925 temp16alen = scale_expansion_zeroelim(axtbclen, axtbc, 2.0 * adx,
02926 temp16a);
02927
02928 axtcclen = scale_expansion_zeroelim(4, cc, adxtail, axtcc);
02929 temp16blen = scale_expansion_zeroelim(axtcclen, axtcc, bdy, temp16b);
02930
02931 axtbblen = scale_expansion_zeroelim(4, bb, adxtail, axtbb);
02932 temp16clen = scale_expansion_zeroelim(axtbblen, axtbb, -cdy, temp16c);
02933
02934 temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
02935 temp16blen, temp16b, temp32a);
02936 temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
02937 temp32alen, temp32a, temp48);
02938 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
02939 temp48, finother);
02940 finswap = finnow; finnow = finother; finother = finswap;
02941 }
02942 if (adytail != 0.0) {
02943 aytbclen = scale_expansion_zeroelim(4, bc, adytail, aytbc);
02944 temp16alen = scale_expansion_zeroelim(aytbclen, aytbc, 2.0 * ady,
02945 temp16a);
02946
02947 aytbblen = scale_expansion_zeroelim(4, bb, adytail, aytbb);
02948 temp16blen = scale_expansion_zeroelim(aytbblen, aytbb, cdx, temp16b);
02949
02950 aytcclen = scale_expansion_zeroelim(4, cc, adytail, aytcc);
02951 temp16clen = scale_expansion_zeroelim(aytcclen, aytcc, -bdx, temp16c);
02952
02953 temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
02954 temp16blen, temp16b, temp32a);
02955 temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
02956 temp32alen, temp32a, temp48);
02957 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
02958 temp48, finother);
02959 finswap = finnow; finnow = finother; finother = finswap;
02960 }
02961 if (bdxtail != 0.0) {
02962 bxtcalen = scale_expansion_zeroelim(4, ca, bdxtail, bxtca);
02963 temp16alen = scale_expansion_zeroelim(bxtcalen, bxtca, 2.0 * bdx,
02964 temp16a);
02965
02966 bxtaalen = scale_expansion_zeroelim(4, aa, bdxtail, bxtaa);
02967 temp16blen = scale_expansion_zeroelim(bxtaalen, bxtaa, cdy, temp16b);
02968
02969 bxtcclen = scale_expansion_zeroelim(4, cc, bdxtail, bxtcc);
02970 temp16clen = scale_expansion_zeroelim(bxtcclen, bxtcc, -ady, temp16c);
02971
02972 temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
02973 temp16blen, temp16b, temp32a);
02974 temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
02975 temp32alen, temp32a, temp48);
02976 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
02977 temp48, finother);
02978 finswap = finnow; finnow = finother; finother = finswap;
02979 }
02980 if (bdytail != 0.0) {
02981 bytcalen = scale_expansion_zeroelim(4, ca, bdytail, bytca);
02982 temp16alen = scale_expansion_zeroelim(bytcalen, bytca, 2.0 * bdy,
02983 temp16a);
02984
02985 bytcclen = scale_expansion_zeroelim(4, cc, bdytail, bytcc);
02986 temp16blen = scale_expansion_zeroelim(bytcclen, bytcc, adx, temp16b);
02987
02988 bytaalen = scale_expansion_zeroelim(4, aa, bdytail, bytaa);
02989 temp16clen = scale_expansion_zeroelim(bytaalen, bytaa, -cdx, temp16c);
02990
02991 temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
02992 temp16blen, temp16b, temp32a);
02993 temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
02994 temp32alen, temp32a, temp48);
02995 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
02996 temp48, finother);
02997 finswap = finnow; finnow = finother; finother = finswap;
02998 }
02999 if (cdxtail != 0.0) {
03000 cxtablen = scale_expansion_zeroelim(4, ab, cdxtail, cxtab);
03001 temp16alen = scale_expansion_zeroelim(cxtablen, cxtab, 2.0 * cdx,
03002 temp16a);
03003
03004 cxtbblen = scale_expansion_zeroelim(4, bb, cdxtail, cxtbb);
03005 temp16blen = scale_expansion_zeroelim(cxtbblen, cxtbb, ady, temp16b);
03006
03007 cxtaalen = scale_expansion_zeroelim(4, aa, cdxtail, cxtaa);
03008 temp16clen = scale_expansion_zeroelim(cxtaalen, cxtaa, -bdy, temp16c);
03009
03010 temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
03011 temp16blen, temp16b, temp32a);
03012 temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
03013 temp32alen, temp32a, temp48);
03014 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
03015 temp48, finother);
03016 finswap = finnow; finnow = finother; finother = finswap;
03017 }
03018 if (cdytail != 0.0) {
03019 cytablen = scale_expansion_zeroelim(4, ab, cdytail, cytab);
03020 temp16alen = scale_expansion_zeroelim(cytablen, cytab, 2.0 * cdy,
03021 temp16a);
03022
03023 cytaalen = scale_expansion_zeroelim(4, aa, cdytail, cytaa);
03024 temp16blen = scale_expansion_zeroelim(cytaalen, cytaa, bdx, temp16b);
03025
03026 cytbblen = scale_expansion_zeroelim(4, bb, cdytail, cytbb);
03027 temp16clen = scale_expansion_zeroelim(cytbblen, cytbb, -adx, temp16c);
03028
03029 temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
03030 temp16blen, temp16b, temp32a);
03031 temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
03032 temp32alen, temp32a, temp48);
03033 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
03034 temp48, finother);
03035 finswap = finnow; finnow = finother; finother = finswap;
03036 }
03037
03038 if ((adxtail != 0.0) || (adytail != 0.0)) {
03039 if ((bdxtail != 0.0) || (bdytail != 0.0)
03040 || (cdxtail != 0.0) || (cdytail != 0.0)) {
03041 Two_Product(bdxtail, cdy, ti1, ti0);
03042 Two_Product(bdx, cdytail, tj1, tj0);
03043 Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
03044 u[3] = u3;
03045 negate = -bdy;
03046 Two_Product(cdxtail, negate, ti1, ti0);
03047 negate = -bdytail;
03048 Two_Product(cdx, negate, tj1, tj0);
03049 Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
03050 v[3] = v3;
03051 bctlen = fast_expansion_sum_zeroelim(4, u, 4, v, bct);
03052
03053 Two_Product(bdxtail, cdytail, ti1, ti0);
03054 Two_Product(cdxtail, bdytail, tj1, tj0);
03055 Two_Two_Diff(ti1, ti0, tj1, tj0, bctt3, bctt[2], bctt[1], bctt[0]);
03056 bctt[3] = bctt3;
03057 bcttlen = 4;
03058 } else {
03059 bct[0] = 0.0;
03060 bctlen = 1;
03061 bctt[0] = 0.0;
03062 bcttlen = 1;
03063 }
03064
03065 if (adxtail != 0.0) {
03066 temp16alen = scale_expansion_zeroelim(axtbclen, axtbc, adxtail, temp16a);
03067 axtbctlen = scale_expansion_zeroelim(bctlen, bct, adxtail, axtbct);
03068 temp32alen = scale_expansion_zeroelim(axtbctlen, axtbct, 2.0 * adx,
03069 temp32a);
03070 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
03071 temp32alen, temp32a, temp48);
03072 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
03073 temp48, finother);
03074 finswap = finnow; finnow = finother; finother = finswap;
03075 if (bdytail != 0.0) {
03076 temp8len = scale_expansion_zeroelim(4, cc, adxtail, temp8);
03077 temp16alen = scale_expansion_zeroelim(temp8len, temp8, bdytail,
03078 temp16a);
03079 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
03080 temp16a, finother);
03081 finswap = finnow; finnow = finother; finother = finswap;
03082 }
03083 if (cdytail != 0.0) {
03084 temp8len = scale_expansion_zeroelim(4, bb, -adxtail, temp8);
03085 temp16alen = scale_expansion_zeroelim(temp8len, temp8, cdytail,
03086 temp16a);
03087 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
03088 temp16a, finother);
03089 finswap = finnow; finnow = finother; finother = finswap;
03090 }
03091
03092 temp32alen = scale_expansion_zeroelim(axtbctlen, axtbct, adxtail,
03093 temp32a);
03094 axtbcttlen = scale_expansion_zeroelim(bcttlen, bctt, adxtail, axtbctt);
03095 temp16alen = scale_expansion_zeroelim(axtbcttlen, axtbctt, 2.0 * adx,
03096 temp16a);
03097 temp16blen = scale_expansion_zeroelim(axtbcttlen, axtbctt, adxtail,
03098 temp16b);
03099 temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
03100 temp16blen, temp16b, temp32b);
03101 temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
03102 temp32blen, temp32b, temp64);
03103 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
03104 temp64, finother);
03105 finswap = finnow; finnow = finother; finother = finswap;
03106 }
03107 if (adytail != 0.0) {
03108 temp16alen = scale_expansion_zeroelim(aytbclen, aytbc, adytail, temp16a);
03109 aytbctlen = scale_expansion_zeroelim(bctlen, bct, adytail, aytbct);
03110 temp32alen = scale_expansion_zeroelim(aytbctlen, aytbct, 2.0 * ady,
03111 temp32a);
03112 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
03113 temp32alen, temp32a, temp48);
03114 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
03115 temp48, finother);
03116 finswap = finnow; finnow = finother; finother = finswap;
03117
03118
03119 temp32alen = scale_expansion_zeroelim(aytbctlen, aytbct, adytail,
03120 temp32a);
03121 aytbcttlen = scale_expansion_zeroelim(bcttlen, bctt, adytail, aytbctt);
03122 temp16alen = scale_expansion_zeroelim(aytbcttlen, aytbctt, 2.0 * ady,
03123 temp16a);
03124 temp16blen = scale_expansion_zeroelim(aytbcttlen, aytbctt, adytail,
03125 temp16b);
03126 temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
03127 temp16blen, temp16b, temp32b);
03128 temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
03129 temp32blen, temp32b, temp64);
03130 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
03131 temp64, finother);
03132 finswap = finnow; finnow = finother; finother = finswap;
03133 }
03134 }
03135 if ((bdxtail != 0.0) || (bdytail != 0.0)) {
03136 if ((cdxtail != 0.0) || (cdytail != 0.0)
03137 || (adxtail != 0.0) || (adytail != 0.0)) {
03138 Two_Product(cdxtail, ady, ti1, ti0);
03139 Two_Product(cdx, adytail, tj1, tj0);
03140 Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
03141 u[3] = u3;
03142 negate = -cdy;
03143 Two_Product(adxtail, negate, ti1, ti0);
03144 negate = -cdytail;
03145 Two_Product(adx, negate, tj1, tj0);
03146 Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
03147 v[3] = v3;
03148 catlen = fast_expansion_sum_zeroelim(4, u, 4, v, cat);
03149
03150 Two_Product(cdxtail, adytail, ti1, ti0);
03151 Two_Product(adxtail, cdytail, tj1, tj0);
03152 Two_Two_Diff(ti1, ti0, tj1, tj0, catt3, catt[2], catt[1], catt[0]);
03153 catt[3] = catt3;
03154 cattlen = 4;
03155 } else {
03156 cat[0] = 0.0;
03157 catlen = 1;
03158 catt[0] = 0.0;
03159 cattlen = 1;
03160 }
03161
03162 if (bdxtail != 0.0) {
03163 temp16alen = scale_expansion_zeroelim(bxtcalen, bxtca, bdxtail, temp16a);
03164 bxtcatlen = scale_expansion_zeroelim(catlen, cat, bdxtail, bxtcat);
03165 temp32alen = scale_expansion_zeroelim(bxtcatlen, bxtcat, 2.0 * bdx,
03166 temp32a);
03167 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
03168 temp32alen, temp32a, temp48);
03169 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
03170 temp48, finother);
03171 finswap = finnow; finnow = finother; finother = finswap;
03172 if (cdytail != 0.0) {
03173 temp8len = scale_expansion_zeroelim(4, aa, bdxtail, temp8);
03174 temp16alen = scale_expansion_zeroelim(temp8len, temp8, cdytail,
03175 temp16a);
03176 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
03177 temp16a, finother);
03178 finswap = finnow; finnow = finother; finother = finswap;
03179 }
03180 if (adytail != 0.0) {
03181 temp8len = scale_expansion_zeroelim(4, cc, -bdxtail, temp8);
03182 temp16alen = scale_expansion_zeroelim(temp8len, temp8, adytail,
03183 temp16a);
03184 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
03185 temp16a, finother);
03186 finswap = finnow; finnow = finother; finother = finswap;
03187 }
03188
03189 temp32alen = scale_expansion_zeroelim(bxtcatlen, bxtcat, bdxtail,
03190 temp32a);
03191 bxtcattlen = scale_expansion_zeroelim(cattlen, catt, bdxtail, bxtcatt);
03192 temp16alen = scale_expansion_zeroelim(bxtcattlen, bxtcatt, 2.0 * bdx,
03193 temp16a);
03194 temp16blen = scale_expansion_zeroelim(bxtcattlen, bxtcatt, bdxtail,
03195 temp16b);
03196 temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
03197 temp16blen, temp16b, temp32b);
03198 temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
03199 temp32blen, temp32b, temp64);
03200 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
03201 temp64, finother);
03202 finswap = finnow; finnow = finother; finother = finswap;
03203 }
03204 if (bdytail != 0.0) {
03205 temp16alen = scale_expansion_zeroelim(bytcalen, bytca, bdytail, temp16a);
03206 bytcatlen = scale_expansion_zeroelim(catlen, cat, bdytail, bytcat);
03207 temp32alen = scale_expansion_zeroelim(bytcatlen, bytcat, 2.0 * bdy,
03208 temp32a);
03209 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
03210 temp32alen, temp32a, temp48);
03211 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
03212 temp48, finother);
03213 finswap = finnow; finnow = finother; finother = finswap;
03214
03215
03216 temp32alen = scale_expansion_zeroelim(bytcatlen, bytcat, bdytail,
03217 temp32a);
03218 bytcattlen = scale_expansion_zeroelim(cattlen, catt, bdytail, bytcatt);
03219 temp16alen = scale_expansion_zeroelim(bytcattlen, bytcatt, 2.0 * bdy,
03220 temp16a);
03221 temp16blen = scale_expansion_zeroelim(bytcattlen, bytcatt, bdytail,
03222 temp16b);
03223 temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
03224 temp16blen, temp16b, temp32b);
03225 temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
03226 temp32blen, temp32b, temp64);
03227 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
03228 temp64, finother);
03229 finswap = finnow; finnow = finother; finother = finswap;
03230 }
03231 }
03232 if ((cdxtail != 0.0) || (cdytail != 0.0)) {
03233 if ((adxtail != 0.0) || (adytail != 0.0)
03234 || (bdxtail != 0.0) || (bdytail != 0.0)) {
03235 Two_Product(adxtail, bdy, ti1, ti0);
03236 Two_Product(adx, bdytail, tj1, tj0);
03237 Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
03238 u[3] = u3;
03239 negate = -ady;
03240 Two_Product(bdxtail, negate, ti1, ti0);
03241 negate = -adytail;
03242 Two_Product(bdx, negate, tj1, tj0);
03243 Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
03244 v[3] = v3;
03245 abtlen = fast_expansion_sum_zeroelim(4, u, 4, v, abt);
03246
03247 Two_Product(adxtail, bdytail, ti1, ti0);
03248 Two_Product(bdxtail, adytail, tj1, tj0);
03249 Two_Two_Diff(ti1, ti0, tj1, tj0, abtt3, abtt[2], abtt[1], abtt[0]);
03250 abtt[3] = abtt3;
03251 abttlen = 4;
03252 } else {
03253 abt[0] = 0.0;
03254 abtlen = 1;
03255 abtt[0] = 0.0;
03256 abttlen = 1;
03257 }
03258
03259 if (cdxtail != 0.0) {
03260 temp16alen = scale_expansion_zeroelim(cxtablen, cxtab, cdxtail, temp16a);
03261 cxtabtlen = scale_expansion_zeroelim(abtlen, abt, cdxtail, cxtabt);
03262 temp32alen = scale_expansion_zeroelim(cxtabtlen, cxtabt, 2.0 * cdx,
03263 temp32a);
03264 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
03265 temp32alen, temp32a, temp48);
03266 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
03267 temp48, finother);
03268 finswap = finnow; finnow = finother; finother = finswap;
03269 if (adytail != 0.0) {
03270 temp8len = scale_expansion_zeroelim(4, bb, cdxtail, temp8);
03271 temp16alen = scale_expansion_zeroelim(temp8len, temp8, adytail,
03272 temp16a);
03273 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
03274 temp16a, finother);
03275 finswap = finnow; finnow = finother; finother = finswap;
03276 }
03277 if (bdytail != 0.0) {
03278 temp8len = scale_expansion_zeroelim(4, aa, -cdxtail, temp8);
03279 temp16alen = scale_expansion_zeroelim(temp8len, temp8, bdytail,
03280 temp16a);
03281 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
03282 temp16a, finother);
03283 finswap = finnow; finnow = finother; finother = finswap;
03284 }
03285
03286 temp32alen = scale_expansion_zeroelim(cxtabtlen, cxtabt, cdxtail,
03287 temp32a);
03288 cxtabttlen = scale_expansion_zeroelim(abttlen, abtt, cdxtail, cxtabtt);
03289 temp16alen = scale_expansion_zeroelim(cxtabttlen, cxtabtt, 2.0 * cdx,
03290 temp16a);
03291 temp16blen = scale_expansion_zeroelim(cxtabttlen, cxtabtt, cdxtail,
03292 temp16b);
03293 temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
03294 temp16blen, temp16b, temp32b);
03295 temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
03296 temp32blen, temp32b, temp64);
03297 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
03298 temp64, finother);
03299 finswap = finnow; finnow = finother; finother = finswap;
03300 }
03301 if (cdytail != 0.0) {
03302 temp16alen = scale_expansion_zeroelim(cytablen, cytab, cdytail, temp16a);
03303 cytabtlen = scale_expansion_zeroelim(abtlen, abt, cdytail, cytabt);
03304 temp32alen = scale_expansion_zeroelim(cytabtlen, cytabt, 2.0 * cdy,
03305 temp32a);
03306 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
03307 temp32alen, temp32a, temp48);
03308 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
03309 temp48, finother);
03310 finswap = finnow; finnow = finother; finother = finswap;
03311
03312
03313 temp32alen = scale_expansion_zeroelim(cytabtlen, cytabt, cdytail,
03314 temp32a);
03315 cytabttlen = scale_expansion_zeroelim(abttlen, abtt, cdytail, cytabtt);
03316 temp16alen = scale_expansion_zeroelim(cytabttlen, cytabtt, 2.0 * cdy,
03317 temp16a);
03318 temp16blen = scale_expansion_zeroelim(cytabttlen, cytabtt, cdytail,
03319 temp16b);
03320 temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
03321 temp16blen, temp16b, temp32b);
03322 temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
03323 temp32blen, temp32b, temp64);
03324 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
03325 temp64, finother);
03326 finswap = finnow; finnow = finother; finother = finswap;
03327 }
03328 }
03329
03330 return finnow[finlength - 1];
03331 }
03332
03333 float incircle(struct mesh *m, struct behavior *b,
03334 vertex pa, vertex pb, vertex pc, vertex pd)
03335 {
03336 float adx, bdx, cdx, ady, bdy, cdy;
03337 float bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady;
03338 float alift, blift, clift;
03339 float det;
03340 float permanent, errbound;
03341
03342 m->incirclecount++;
03343
03344 adx = pa[0] - pd[0];
03345 bdx = pb[0] - pd[0];
03346 cdx = pc[0] - pd[0];
03347 ady = pa[1] - pd[1];
03348 bdy = pb[1] - pd[1];
03349 cdy = pc[1] - pd[1];
03350
03351 bdxcdy = bdx * cdy;
03352 cdxbdy = cdx * bdy;
03353 alift = adx * adx + ady * ady;
03354
03355 cdxady = cdx * ady;
03356 adxcdy = adx * cdy;
03357 blift = bdx * bdx + bdy * bdy;
03358
03359 adxbdy = adx * bdy;
03360 bdxady = bdx * ady;
03361 clift = cdx * cdx + cdy * cdy;
03362
03363 det = alift * (bdxcdy - cdxbdy)
03364 + blift * (cdxady - adxcdy)
03365 + clift * (adxbdy - bdxady);
03366
03367 if (b->noexact) {
03368 return det;
03369 }
03370
03371 permanent = (Absolute(bdxcdy) + Absolute(cdxbdy)) * alift
03372 + (Absolute(cdxady) + Absolute(adxcdy)) * blift
03373 + (Absolute(adxbdy) + Absolute(bdxady)) * clift;
03374 errbound = iccerrboundA * permanent;
03375 if ((det > errbound) || (-det > errbound)) {
03376 return det;
03377 }
03378
03379 return incircleadapt(pa, pb, pc, pd, permanent);
03380 }
03381
03382
03383
03384
03385
03386
03387
03388
03389
03390
03391
03392
03393
03394
03395
03396
03397
03398
03399
03400
03401
03402
03403
03404 float orient3dadapt(vertex pa, vertex pb, vertex pc, vertex pd,
03405 float aheight, float bheight, float cheight, float dheight,
03406 float permanent)
03407 {
03408 float adx, bdx, cdx, ady, bdy, cdy, adheight, bdheight, cdheight;
03409 float det, errbound;
03410
03411 float bdxcdy1, cdxbdy1, cdxady1, adxcdy1, adxbdy1, bdxady1;
03412 float bdxcdy0, cdxbdy0, cdxady0, adxcdy0, adxbdy0, bdxady0;
03413 float bc[4], ca[4], ab[4];
03414 float bc3, ca3, ab3;
03415 float adet[8], bdet[8], cdet[8];
03416 int alen, blen, clen;
03417 float abdet[16];
03418 int ablen;
03419 float *finnow, *finother, *finswap;
03420 float fin1[192], fin2[192];
03421 int finlength;
03422
03423 float adxtail, bdxtail, cdxtail;
03424 float adytail, bdytail, cdytail;
03425 float adheighttail, bdheighttail, cdheighttail;
03426 float at_blarge, at_clarge;
03427 float bt_clarge, bt_alarge;
03428 float ct_alarge, ct_blarge;
03429 float at_b[4], at_c[4], bt_c[4], bt_a[4], ct_a[4], ct_b[4];
03430 int at_blen, at_clen, bt_clen, bt_alen, ct_alen, ct_blen;
03431 float bdxt_cdy1, cdxt_bdy1, cdxt_ady1;
03432 float adxt_cdy1, adxt_bdy1, bdxt_ady1;
03433 float bdxt_cdy0, cdxt_bdy0, cdxt_ady0;
03434 float adxt_cdy0, adxt_bdy0, bdxt_ady0;
03435 float bdyt_cdx1, cdyt_bdx1, cdyt_adx1;
03436 float adyt_cdx1, adyt_bdx1, bdyt_adx1;
03437 float bdyt_cdx0, cdyt_bdx0, cdyt_adx0;
03438 float adyt_cdx0, adyt_bdx0, bdyt_adx0;
03439 float bct[8], cat[8], abt[8];
03440 int bctlen, catlen, abtlen;
03441 float bdxt_cdyt1, cdxt_bdyt1, cdxt_adyt1;
03442 float adxt_cdyt1, adxt_bdyt1, bdxt_adyt1;
03443 float bdxt_cdyt0, cdxt_bdyt0, cdxt_adyt0;
03444 float adxt_cdyt0, adxt_bdyt0, bdxt_adyt0;
03445 float u[4], v[12], w[16];
03446 float u3;
03447 int vlength, wlength;
03448 float negate;
03449
03450 float bvirt;
03451 float avirt, bround, around;
03452 float c;
03453 float abig;
03454 float ahi, alo, bhi, blo;
03455 float err1, err2, err3;
03456 float _i, _j, _k;
03457 float _0;
03458
03459 adx = (float) (pa[0] - pd[0]);
03460 bdx = (float) (pb[0] - pd[0]);
03461 cdx = (float) (pc[0] - pd[0]);
03462 ady = (float) (pa[1] - pd[1]);
03463 bdy = (float) (pb[1] - pd[1]);
03464 cdy = (float) (pc[1] - pd[1]);
03465 adheight = (float) (aheight - dheight);
03466 bdheight = (float) (bheight - dheight);
03467 cdheight = (float) (cheight - dheight);
03468
03469 Two_Product(bdx, cdy, bdxcdy1, bdxcdy0);
03470 Two_Product(cdx, bdy, cdxbdy1, cdxbdy0);
03471 Two_Two_Diff(bdxcdy1, bdxcdy0, cdxbdy1, cdxbdy0, bc3, bc[2], bc[1], bc[0]);
03472 bc[3] = bc3;
03473 alen = scale_expansion_zeroelim(4, bc, adheight, adet);
03474
03475 Two_Product(cdx, ady, cdxady1, cdxady0);
03476 Two_Product(adx, cdy, adxcdy1, adxcdy0);
03477 Two_Two_Diff(cdxady1, cdxady0, adxcdy1, adxcdy0, ca3, ca[2], ca[1], ca[0]);
03478 ca[3] = ca3;
03479 blen = scale_expansion_zeroelim(4, ca, bdheight, bdet);
03480
03481 Two_Product(adx, bdy, adxbdy1, adxbdy0);
03482 Two_Product(bdx, ady, bdxady1, bdxady0);
03483 Two_Two_Diff(adxbdy1, adxbdy0, bdxady1, bdxady0, ab3, ab[2], ab[1], ab[0]);
03484 ab[3] = ab3;
03485 clen = scale_expansion_zeroelim(4, ab, cdheight, cdet);
03486
03487 ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
03488 finlength = fast_expansion_sum_zeroelim(ablen, abdet, clen, cdet, fin1);
03489
03490 det = estimate(finlength, fin1);
03491 errbound = o3derrboundB * permanent;
03492 if ((det >= errbound) || (-det >= errbound)) {
03493 return det;
03494 }
03495
03496 Two_Diff_Tail(pa[0], pd[0], adx, adxtail);
03497 Two_Diff_Tail(pb[0], pd[0], bdx, bdxtail);
03498 Two_Diff_Tail(pc[0], pd[0], cdx, cdxtail);
03499 Two_Diff_Tail(pa[1], pd[1], ady, adytail);
03500 Two_Diff_Tail(pb[1], pd[1], bdy, bdytail);
03501 Two_Diff_Tail(pc[1], pd[1], cdy, cdytail);
03502 Two_Diff_Tail(aheight, dheight, adheight, adheighttail);
03503 Two_Diff_Tail(bheight, dheight, bdheight, bdheighttail);
03504 Two_Diff_Tail(cheight, dheight, cdheight, cdheighttail);
03505
03506 if ((adxtail == 0.0) && (bdxtail == 0.0) && (cdxtail == 0.0) &&
03507 (adytail == 0.0) && (bdytail == 0.0) && (cdytail == 0.0) &&
03508 (adheighttail == 0.0) &&
03509 (bdheighttail == 0.0) &&
03510 (cdheighttail == 0.0)) {
03511 return det;
03512 }
03513
03514 errbound = o3derrboundC * permanent + resulterrbound * Absolute(det);
03515 det += (adheight * ((bdx * cdytail + cdy * bdxtail) -
03516 (bdy * cdxtail + cdx * bdytail)) +
03517 adheighttail * (bdx * cdy - bdy * cdx)) +
03518 (bdheight * ((cdx * adytail + ady * cdxtail) -
03519 (cdy * adxtail + adx * cdytail)) +
03520 bdheighttail * (cdx * ady - cdy * adx)) +
03521 (cdheight * ((adx * bdytail + bdy * adxtail) -
03522 (ady * bdxtail + bdx * adytail)) +
03523 cdheighttail * (adx * bdy - ady * bdx));
03524 if ((det >= errbound) || (-det >= errbound)) {
03525 return det;
03526 }
03527
03528 finnow = fin1;
03529 finother = fin2;
03530
03531 if (adxtail == 0.0) {
03532 if (adytail == 0.0) {
03533 at_b[0] = 0.0;
03534 at_blen = 1;
03535 at_c[0] = 0.0;
03536 at_clen = 1;
03537 } else {
03538 negate = -adytail;
03539 Two_Product(negate, bdx, at_blarge, at_b[0]);
03540 at_b[1] = at_blarge;
03541 at_blen = 2;
03542 Two_Product(adytail, cdx, at_clarge, at_c[0]);
03543 at_c[1] = at_clarge;
03544 at_clen = 2;
03545 }
03546 } else {
03547 if (adytail == 0.0) {
03548 Two_Product(adxtail, bdy, at_blarge, at_b[0]);
03549 at_b[1] = at_blarge;
03550 at_blen = 2;
03551 negate = -adxtail;
03552 Two_Product(negate, cdy, at_clarge, at_c[0]);
03553 at_c[1] = at_clarge;
03554 at_clen = 2;
03555 } else {
03556 Two_Product(adxtail, bdy, adxt_bdy1, adxt_bdy0);
03557 Two_Product(adytail, bdx, adyt_bdx1, adyt_bdx0);
03558 Two_Two_Diff(adxt_bdy1, adxt_bdy0, adyt_bdx1, adyt_bdx0,
03559 at_blarge, at_b[2], at_b[1], at_b[0]);
03560 at_b[3] = at_blarge;
03561 at_blen = 4;
03562 Two_Product(adytail, cdx, adyt_cdx1, adyt_cdx0);
03563 Two_Product(adxtail, cdy, adxt_cdy1, adxt_cdy0);
03564 Two_Two_Diff(adyt_cdx1, adyt_cdx0, adxt_cdy1, adxt_cdy0,
03565 at_clarge, at_c[2], at_c[1], at_c[0]);
03566 at_c[3] = at_clarge;
03567 at_clen = 4;
03568 }
03569 }
03570 if (bdxtail == 0.0) {
03571 if (bdytail == 0.0) {
03572 bt_c[0] = 0.0;
03573 bt_clen = 1;
03574 bt_a[0] = 0.0;
03575 bt_alen = 1;
03576 } else {
03577 negate = -bdytail;
03578 Two_Product(negate, cdx, bt_clarge, bt_c[0]);
03579 bt_c[1] = bt_clarge;
03580 bt_clen = 2;
03581 Two_Product(bdytail, adx, bt_alarge, bt_a[0]);
03582 bt_a[1] = bt_alarge;
03583 bt_alen = 2;
03584 }
03585 } else {
03586 if (bdytail == 0.0) {
03587 Two_Product(bdxtail, cdy, bt_clarge, bt_c[0]);
03588 bt_c[1] = bt_clarge;
03589 bt_clen = 2;
03590 negate = -bdxtail;
03591 Two_Product(negate, ady, bt_alarge, bt_a[0]);
03592 bt_a[1] = bt_alarge;
03593 bt_alen = 2;
03594 } else {
03595 Two_Product(bdxtail, cdy, bdxt_cdy1, bdxt_cdy0);
03596 Two_Product(bdytail, cdx, bdyt_cdx1, bdyt_cdx0);
03597 Two_Two_Diff(bdxt_cdy1, bdxt_cdy0, bdyt_cdx1, bdyt_cdx0,
03598 bt_clarge, bt_c[2], bt_c[1], bt_c[0]);
03599 bt_c[3] = bt_clarge;
03600 bt_clen = 4;
03601 Two_Product(bdytail, adx, bdyt_adx1, bdyt_adx0);
03602 Two_Product(bdxtail, ady, bdxt_ady1, bdxt_ady0);
03603 Two_Two_Diff(bdyt_adx1, bdyt_adx0, bdxt_ady1, bdxt_ady0,
03604 bt_alarge, bt_a[2], bt_a[1], bt_a[0]);
03605 bt_a[3] = bt_alarge;
03606 bt_alen = 4;
03607 }
03608 }
03609 if (cdxtail == 0.0) {
03610 if (cdytail == 0.0) {
03611 ct_a[0] = 0.0;
03612 ct_alen = 1;
03613 ct_b[0] = 0.0;
03614 ct_blen = 1;
03615 } else {
03616 negate = -cdytail;
03617 Two_Product(negate, adx, ct_alarge, ct_a[0]);
03618 ct_a[1] = ct_alarge;
03619 ct_alen = 2;
03620 Two_Product(cdytail, bdx, ct_blarge, ct_b[0]);
03621 ct_b[1] = ct_blarge;
03622 ct_blen = 2;
03623 }
03624 } else {
03625 if (cdytail == 0.0) {
03626 Two_Product(cdxtail, ady, ct_alarge, ct_a[0]);
03627 ct_a[1] = ct_alarge;
03628 ct_alen = 2;
03629 negate = -cdxtail;
03630 Two_Product(negate, bdy, ct_blarge, ct_b[0]);
03631 ct_b[1] = ct_blarge;
03632 ct_blen = 2;
03633 } else {
03634 Two_Product(cdxtail, ady, cdxt_ady1, cdxt_ady0);
03635 Two_Product(cdytail, adx, cdyt_adx1, cdyt_adx0);
03636 Two_Two_Diff(cdxt_ady1, cdxt_ady0, cdyt_adx1, cdyt_adx0,
03637 ct_alarge, ct_a[2], ct_a[1], ct_a[0]);
03638 ct_a[3] = ct_alarge;
03639 ct_alen = 4;
03640 Two_Product(cdytail, bdx, cdyt_bdx1, cdyt_bdx0);
03641 Two_Product(cdxtail, bdy, cdxt_bdy1, cdxt_bdy0);
03642 Two_Two_Diff(cdyt_bdx1, cdyt_bdx0, cdxt_bdy1, cdxt_bdy0,
03643 ct_blarge, ct_b[2], ct_b[1], ct_b[0]);
03644 ct_b[3] = ct_blarge;
03645 ct_blen = 4;
03646 }
03647 }
03648
03649 bctlen = fast_expansion_sum_zeroelim(bt_clen, bt_c, ct_blen, ct_b, bct);
03650 wlength = scale_expansion_zeroelim(bctlen, bct, adheight, w);
03651 finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
03652 finother);
03653 finswap = finnow; finnow = finother; finother = finswap;
03654
03655 catlen = fast_expansion_sum_zeroelim(ct_alen, ct_a, at_clen, at_c, cat);
03656 wlength = scale_expansion_zeroelim(catlen, cat, bdheight, w);
03657 finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
03658 finother);
03659 finswap = finnow; finnow = finother; finother = finswap;
03660
03661 abtlen = fast_expansion_sum_zeroelim(at_blen, at_b, bt_alen, bt_a, abt);
03662 wlength = scale_expansion_zeroelim(abtlen, abt, cdheight, w);
03663 finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
03664 finother);
03665 finswap = finnow; finnow = finother; finother = finswap;
03666
03667 if (adheighttail != 0.0) {
03668 vlength = scale_expansion_zeroelim(4, bc, adheighttail, v);
03669 finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v,
03670 finother);
03671 finswap = finnow; finnow = finother; finother = finswap;
03672 }
03673 if (bdheighttail != 0.0) {
03674 vlength = scale_expansion_zeroelim(4, ca, bdheighttail, v);
03675 finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v,
03676 finother);
03677 finswap = finnow; finnow = finother; finother = finswap;
03678 }
03679 if (cdheighttail != 0.0) {
03680 vlength = scale_expansion_zeroelim(4, ab, cdheighttail, v);
03681 finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v,
03682 finother);
03683 finswap = finnow; finnow = finother; finother = finswap;
03684 }
03685
03686 if (adxtail != 0.0) {
03687 if (bdytail != 0.0) {
03688 Two_Product(adxtail, bdytail, adxt_bdyt1, adxt_bdyt0);
03689 Two_One_Product(adxt_bdyt1, adxt_bdyt0, cdheight, u3, u[2], u[1], u[0]);
03690 u[3] = u3;
03691 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
03692 finother);
03693 finswap = finnow; finnow = finother; finother = finswap;
03694 if (cdheighttail != 0.0) {
03695 Two_One_Product(adxt_bdyt1, adxt_bdyt0, cdheighttail,
03696 u3, u[2], u[1], u[0]);
03697 u[3] = u3;
03698 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
03699 finother);
03700 finswap = finnow; finnow = finother; finother = finswap;
03701 }
03702 }
03703 if (cdytail != 0.0) {
03704 negate = -adxtail;
03705 Two_Product(negate, cdytail, adxt_cdyt1, adxt_cdyt0);
03706 Two_One_Product(adxt_cdyt1, adxt_cdyt0, bdheight, u3, u[2], u[1], u[0]);
03707 u[3] = u3;
03708 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
03709 finother);
03710 finswap = finnow; finnow = finother; finother = finswap;
03711 if (bdheighttail != 0.0) {
03712 Two_One_Product(adxt_cdyt1, adxt_cdyt0, bdheighttail,
03713 u3, u[2], u[1], u[0]);
03714 u[3] = u3;
03715 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
03716 finother);
03717 finswap = finnow; finnow = finother; finother = finswap;
03718 }
03719 }
03720 }
03721 if (bdxtail != 0.0) {
03722 if (cdytail != 0.0) {
03723 Two_Product(bdxtail, cdytail, bdxt_cdyt1, bdxt_cdyt0);
03724 Two_One_Product(bdxt_cdyt1, bdxt_cdyt0, adheight, u3, u[2], u[1], u[0]);
03725 u[3] = u3;
03726 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
03727 finother);
03728 finswap = finnow; finnow = finother; finother = finswap;
03729 if (adheighttail != 0.0) {
03730 Two_One_Product(bdxt_cdyt1, bdxt_cdyt0, adheighttail,
03731 u3, u[2], u[1], u[0]);
03732 u[3] = u3;
03733 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
03734 finother);
03735 finswap = finnow; finnow = finother; finother = finswap;
03736 }
03737 }
03738 if (adytail != 0.0) {
03739 negate = -bdxtail;
03740 Two_Product(negate, adytail, bdxt_adyt1, bdxt_adyt0);
03741 Two_One_Product(bdxt_adyt1, bdxt_adyt0, cdheight, u3, u[2], u[1], u[0]);
03742 u[3] = u3;
03743 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
03744 finother);
03745 finswap = finnow; finnow = finother; finother = finswap;
03746 if (cdheighttail != 0.0) {
03747 Two_One_Product(bdxt_adyt1, bdxt_adyt0, cdheighttail,
03748 u3, u[2], u[1], u[0]);
03749 u[3] = u3;
03750 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
03751 finother);
03752 finswap = finnow; finnow = finother; finother = finswap;
03753 }
03754 }
03755 }
03756 if (cdxtail != 0.0) {
03757 if (adytail != 0.0) {
03758 Two_Product(cdxtail, adytail, cdxt_adyt1, cdxt_adyt0);
03759 Two_One_Product(cdxt_adyt1, cdxt_adyt0, bdheight, u3, u[2], u[1], u[0]);
03760 u[3] = u3;
03761 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
03762 finother);
03763 finswap = finnow; finnow = finother; finother = finswap;
03764 if (bdheighttail != 0.0) {
03765 Two_One_Product(cdxt_adyt1, cdxt_adyt0, bdheighttail,
03766 u3, u[2], u[1], u[0]);
03767 u[3] = u3;
03768 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
03769 finother);
03770 finswap = finnow; finnow = finother; finother = finswap;
03771 }
03772 }
03773 if (bdytail != 0.0) {
03774 negate = -cdxtail;
03775 Two_Product(negate, bdytail, cdxt_bdyt1, cdxt_bdyt0);
03776 Two_One_Product(cdxt_bdyt1, cdxt_bdyt0, adheight, u3, u[2], u[1], u[0]);
03777 u[3] = u3;
03778 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
03779 finother);
03780 finswap = finnow; finnow = finother; finother = finswap;
03781 if (adheighttail != 0.0) {
03782 Two_One_Product(cdxt_bdyt1, cdxt_bdyt0, adheighttail,
03783 u3, u[2], u[1], u[0]);
03784 u[3] = u3;
03785 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
03786 finother);
03787 finswap = finnow; finnow = finother; finother = finswap;
03788 }
03789 }
03790 }
03791
03792 if (adheighttail != 0.0) {
03793 wlength = scale_expansion_zeroelim(bctlen, bct, adheighttail, w);
03794 finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
03795 finother);
03796 finswap = finnow; finnow = finother; finother = finswap;
03797 }
03798 if (bdheighttail != 0.0) {
03799 wlength = scale_expansion_zeroelim(catlen, cat, bdheighttail, w);
03800 finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
03801 finother);
03802 finswap = finnow; finnow = finother; finother = finswap;
03803 }
03804 if (cdheighttail != 0.0) {
03805 wlength = scale_expansion_zeroelim(abtlen, abt, cdheighttail, w);
03806 finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
03807 finother);
03808 finswap = finnow; finnow = finother; finother = finswap;
03809 }
03810
03811 return finnow[finlength - 1];
03812 }
03813
03814 float orient3d(struct mesh *m, struct behavior *b,
03815 vertex pa, vertex pb, vertex pc, vertex pd,
03816 float aheight, float bheight, float cheight, float dheight)
03817 {
03818 float adx, bdx, cdx, ady, bdy, cdy, adheight, bdheight, cdheight;
03819 float bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady;
03820 float det;
03821 float permanent, errbound;
03822
03823 m->orient3dcount++;
03824
03825 adx = pa[0] - pd[0];
03826 bdx = pb[0] - pd[0];
03827 cdx = pc[0] - pd[0];
03828 ady = pa[1] - pd[1];
03829 bdy = pb[1] - pd[1];
03830 cdy = pc[1] - pd[1];
03831 adheight = aheight - dheight;
03832 bdheight = bheight - dheight;
03833 cdheight = cheight - dheight;
03834
03835 bdxcdy = bdx * cdy;
03836 cdxbdy = cdx * bdy;
03837
03838 cdxady = cdx * ady;
03839 adxcdy = adx * cdy;
03840
03841 adxbdy = adx * bdy;
03842 bdxady = bdx * ady;
03843
03844 det = adheight * (bdxcdy - cdxbdy)
03845 + bdheight * (cdxady - adxcdy)
03846 + cdheight * (adxbdy - bdxady);
03847
03848 if (b->noexact) {
03849 return det;
03850 }
03851
03852 permanent = (Absolute(bdxcdy) + Absolute(cdxbdy)) * Absolute(adheight)
03853 + (Absolute(cdxady) + Absolute(adxcdy)) * Absolute(bdheight)
03854 + (Absolute(adxbdy) + Absolute(bdxady)) * Absolute(cdheight);
03855 errbound = o3derrboundA * permanent;
03856 if ((det > errbound) || (-det > errbound)) {
03857 return det;
03858 }
03859
03860 return orient3dadapt(pa, pb, pc, pd, aheight, bheight, cheight, dheight,
03861 permanent);
03862 }
03863
03864
03865
03866
03867
03868
03869
03870
03871
03872
03873
03874
03875
03876
03877
03878
03879
03880
03881
03882 float nonregular(struct mesh *m, struct behavior *b,
03883 vertex pa, vertex pb, vertex pc, vertex pd)
03884 {
03885 if (b->weighted == 0) {
03886 return incircle(m, b, pa, pb, pc, pd);
03887 } else if (b->weighted == 1) {
03888 return orient3d(m, b, pa, pb, pc, pd,
03889 pa[0] * pa[0] + pa[1] * pa[1] - pa[2],
03890 pb[0] * pb[0] + pb[1] * pb[1] - pb[2],
03891 pc[0] * pc[0] + pc[1] * pc[1] - pc[2],
03892 pd[0] * pd[0] + pd[1] * pd[1] - pd[2]);
03893 } else {
03894 return orient3d(m, b, pa, pb, pc, pd, pa[2], pb[2], pc[2], pd[2]);
03895 }
03896 }
03897
03898
03899
03900
03901
03902
03903
03904
03905
03906
03907
03908
03909
03910
03911
03912 void findcircumcenter(struct mesh *m, struct behavior *b,
03913 vertex torg, vertex tdest, vertex tapex,
03914 vertex circumcenter, float *xi, float *eta, int offcenter)
03915 {
03916 float xdo, ydo, xao, yao;
03917 float dodist, aodist, dadist;
03918 float denominator;
03919 float dx, dy, dxoff, dyoff;
03920
03921 m->circumcentercount++;
03922
03923
03924 xdo = tdest[0] - torg[0];
03925 ydo = tdest[1] - torg[1];
03926 xao = tapex[0] - torg[0];
03927 yao = tapex[1] - torg[1];
03928 dodist = xdo * xdo + ydo * ydo;
03929 aodist = xao * xao + yao * yao;
03930 dadist = (tdest[0] - tapex[0]) * (tdest[0] - tapex[0]) +
03931 (tdest[1] - tapex[1]) * (tdest[1] - tapex[1]);
03932 if (b->noexact) {
03933 denominator = 0.5 / (xdo * yao - xao * ydo);
03934 } else {
03935
03936
03937
03938 denominator = 0.5 / counterclockwise(m, b, tdest, tapex, torg);
03939
03940 m->counterclockcount--;
03941 }
03942 dx = (yao * dodist - ydo * aodist) * denominator;
03943 dy = (xdo * aodist - xao * dodist) * denominator;
03944
03945
03946
03947
03948
03949
03950 if ((dodist < aodist) && (dodist < dadist)) {
03951 if (offcenter && (b->offconstant > 0.0)) {
03952
03953 dxoff = 0.5 * xdo - b->offconstant * ydo;
03954 dyoff = 0.5 * ydo + b->offconstant * xdo;
03955
03956
03957 if (dxoff * dxoff + dyoff * dyoff < dx * dx + dy * dy) {
03958 dx = dxoff;
03959 dy = dyoff;
03960 }
03961 }
03962 } else if (aodist < dadist) {
03963 if (offcenter && (b->offconstant > 0.0)) {
03964 dxoff = 0.5 * xao + b->offconstant * yao;
03965 dyoff = 0.5 * yao - b->offconstant * xao;
03966
03967
03968 if (dxoff * dxoff + dyoff * dyoff < dx * dx + dy * dy) {
03969 dx = dxoff;
03970 dy = dyoff;
03971 }
03972 }
03973 } else {
03974 if (offcenter && (b->offconstant > 0.0)) {
03975 dxoff = 0.5 * (tapex[0] - tdest[0]) -
03976 b->offconstant * (tapex[1] - tdest[1]);
03977 dyoff = 0.5 * (tapex[1] - tdest[1]) +
03978 b->offconstant * (tapex[0] - tdest[0]);
03979
03980
03981 if (dxoff * dxoff + dyoff * dyoff <
03982 (dx - xdo) * (dx - xdo) + (dy - ydo) * (dy - ydo)) {
03983 dx = xdo + dxoff;
03984 dy = ydo + dyoff;
03985 }
03986 }
03987 }
03988
03989 circumcenter[0] = torg[0] + dx;
03990 circumcenter[1] = torg[1] + dy;
03991
03992
03993
03994
03995
03996
03997 *xi = (yao * dx - xao * dy) * (2.0 * denominator);
03998 *eta = (xdo * dy - ydo * dx) * (2.0 * denominator);
03999 }
04000
04003
04004
04005
04006
04007
04008
04009
04010
04011 void triangleinit(struct mesh *m)
04012 {
04013 poolzero(&m->vertices);
04014 poolzero(&m->triangles);
04015 poolzero(&m->subsegs);
04016 poolzero(&m->viri);
04017 poolzero(&m->badsubsegs);
04018 poolzero(&m->badtriangles);
04019 poolzero(&m->flipstackers);
04020 poolzero(&m->splaynodes);
04021
04022 m->recenttri.tri = (triangle *) NULL;
04023 m->undeads = 0;
04024 m->samples = 1;
04025 m->checksegments = 0;
04026 m->checkquality = 0;
04027 m->incirclecount = m->counterclockcount = m->orient3dcount = 0;
04028 m->hyperbolacount = m->circletopcount = m->circumcentercount = 0;
04029 randomseed = 1;
04030
04031 exactinit();
04032 }
04033
04034
04035
04036
04037
04038
04039
04040
04041
04042
04043
04044 unsigned long randomnation(unsigned int choices)
04045 {
04046 randomseed = (randomseed * 1366l + 150889l) % 714025l;
04047 return randomseed / (714025l / choices + 1);
04048 }
04049
04050
04054
04055
04056
04057
04058
04059
04060
04061
04062
04063
04064
04065
04066
04067
04068 void makevertexmap(struct mesh *m, struct behavior *b)
04069 {
04070 struct otri triangleloop;
04071 vertex triorg;
04072
04073 if (b->verbose) {
04074 printf(" Constructing mapping from vertices to triangles.\n");
04075 }
04076 traversalinit(&m->triangles);
04077 triangleloop.tri = triangletraverse(m);
04078 while (triangleloop.tri != (triangle *) NULL) {
04079
04080 for (triangleloop.orient = 0; triangleloop.orient < 3;
04081 triangleloop.orient++) {
04082 org(triangleloop, triorg);
04083 setvertex2tri(triorg, encode(triangleloop));
04084 }
04085 triangleloop.tri = triangletraverse(m);
04086 }
04087 }
04088
04089
04090
04091
04092
04093
04094
04095
04096
04097
04098
04099
04100
04101
04102
04103
04104
04105
04106
04107
04108
04109
04110
04111
04112
04113
04114
04115
04116
04117
04118
04119
04120
04121
04122
04123
04124
04125
04126
04127
04128
04129
04130
04131
04132
04133
04134
04135
04136
04137
04138
04139
04140
04141
04142
04143
04144
04145
04146
04147
04148
04149
04150
04151
04152
04153
04154
04155
04156 enum locateresult preciselocate(struct mesh *m, struct behavior *b,
04157 vertex searchpoint, struct otri *searchtri,
04158 int stopatsubsegment)
04159 {
04160 struct otri backtracktri;
04161 struct osub checkedge;
04162 vertex forg, fdest, fapex;
04163 float orgorient, destorient;
04164 int moveleft;
04165 triangle ptr;
04166 subseg sptr;
04167
04168 if (b->verbose > 2) {
04169 printf(" Searching for point (%.12g, %.12g).\n",
04170 searchpoint[0], searchpoint[1]);
04171 }
04172
04173 org(*searchtri, forg);
04174 dest(*searchtri, fdest);
04175 apex(*searchtri, fapex);
04176 while (1) {
04177 if (b->verbose > 2) {
04178 printf(" At (%.12g, %.12g) (%.12g, %.12g) (%.12g, %.12g)\n",
04179 forg[0], forg[1], fdest[0], fdest[1], fapex[0], fapex[1]);
04180 }
04181
04182 if ((fapex[0] == searchpoint[0]) && (fapex[1] == searchpoint[1])) {
04183 lprevself(*searchtri);
04184 return ONVERTEX;
04185 }
04186
04187
04188 destorient = counterclockwise(m, b, forg, fapex, searchpoint);
04189
04190
04191 orgorient = counterclockwise(m, b, fapex, fdest, searchpoint);
04192 if (destorient > 0.0) {
04193 if (orgorient > 0.0) {
04194
04195
04196
04197
04198
04199 moveleft = (fapex[0] - searchpoint[0]) * (fdest[0] - forg[0]) +
04200 (fapex[1] - searchpoint[1]) * (fdest[1] - forg[1]) > 0.0;
04201 } else {
04202 moveleft = 1;
04203 }
04204 } else {
04205 if (orgorient > 0.0) {
04206 moveleft = 0;
04207 } else {
04208
04209
04210 if (destorient == 0.0) {
04211 lprevself(*searchtri);
04212 return ONEDGE;
04213 }
04214 if (orgorient == 0.0) {
04215 lnextself(*searchtri);
04216 return ONEDGE;
04217 }
04218 return INTRIANGLE;
04219 }
04220 }
04221
04222
04223
04224
04225 if (moveleft) {
04226 lprev(*searchtri, backtracktri);
04227 fdest = fapex;
04228 } else {
04229 lnext(*searchtri, backtracktri);
04230 forg = fapex;
04231 }
04232 sym(backtracktri, *searchtri);
04233
04234 if (m->checksegments && stopatsubsegment) {
04235
04236 tspivot(backtracktri, checkedge);
04237 if (checkedge.ss != m->dummysub) {
04238
04239 otricopy(backtracktri, *searchtri);
04240 return OUTSIDE;
04241 }
04242 }
04243
04244 if (searchtri->tri == m->dummytri) {
04245
04246 otricopy(backtracktri, *searchtri);
04247 return OUTSIDE;
04248 }
04249
04250 apex(*searchtri, fapex);
04251 }
04252 }
04253
04254
04255
04256
04257
04258
04259
04260
04261
04262
04263
04264
04265
04266
04267
04268
04269
04270
04271
04272
04273
04274
04275
04276
04277
04278
04279
04280
04281
04282
04283
04284
04285
04286
04287
04288
04289
04290 enum locateresult locate(struct mesh *m, struct behavior *b,
04291 vertex searchpoint, struct otri *searchtri)
04292 {
04293 int **sampleblock;
04294 char *firsttri;
04295 struct otri sampletri;
04296 vertex torg, tdest;
04297 unsigned long alignptr;
04298 float searchdist, dist;
04299 float ahead;
04300 long samplesperblock, totalsamplesleft, samplesleft;
04301 long population, totalpopulation;
04302 triangle ptr;
04303
04304 if (b->verbose > 2) {
04305 printf(" Randomly sampling for a triangle near point (%.12g, %.12g).\n",
04306 searchpoint[0], searchpoint[1]);
04307 }
04308
04309
04310 org(*searchtri, torg);
04311 searchdist = (searchpoint[0] - torg[0]) * (searchpoint[0] - torg[0]) +
04312 (searchpoint[1] - torg[1]) * (searchpoint[1] - torg[1]);
04313 if (b->verbose > 2) {
04314 printf(" Boundary triangle has origin (%.12g, %.12g).\n",
04315 torg[0], torg[1]);
04316 }
04317
04318
04319
04320 if (m->recenttri.tri != (triangle *) NULL) {
04321 if (!deadtri(m->recenttri.tri)) {
04322 org(m->recenttri, torg);
04323 if ((torg[0] == searchpoint[0]) && (torg[1] == searchpoint[1])) {
04324 otricopy(m->recenttri, *searchtri);
04325 return ONVERTEX;
04326 }
04327 dist = (searchpoint[0] - torg[0]) * (searchpoint[0] - torg[0]) +
04328 (searchpoint[1] - torg[1]) * (searchpoint[1] - torg[1]);
04329 if (dist < searchdist) {
04330 otricopy(m->recenttri, *searchtri);
04331 searchdist = dist;
04332 if (b->verbose > 2) {
04333 printf(" Choosing recent triangle with origin (%.12g, %.12g).\n",
04334 torg[0], torg[1]);
04335 }
04336 }
04337 }
04338 }
04339
04340
04341
04342
04343
04344 while (SAMPLEFACTOR * m->samples * m->samples * m->samples <
04345 m->triangles.items) {
04346 m->samples++;
04347 }
04348
04349
04350
04351
04352
04353 samplesperblock = (m->samples * TRIPERBLOCK - 1) / m->triangles.maxitems + 1;
04354
04355
04356 samplesleft = (m->samples * m->triangles.itemsfirstblock - 1) /
04357 m->triangles.maxitems + 1;
04358 totalsamplesleft = m->samples;
04359 population = m->triangles.itemsfirstblock;
04360 totalpopulation = m->triangles.maxitems;
04361 sampleblock = m->triangles.firstblock;
04362 sampletri.orient = 0;
04363 while (totalsamplesleft > 0) {
04364
04365 if (population > totalpopulation) {
04366 population = totalpopulation;
04367 }
04368
04369 alignptr = (unsigned long) (sampleblock + 1);
04370 firsttri = (char *) (alignptr +
04371 (unsigned long) m->triangles.alignbytes -
04372 (alignptr %
04373 (unsigned long) m->triangles.alignbytes));
04374
04375
04376 do {
04377 sampletri.tri = (triangle *) (firsttri +
04378 (randomnation((unsigned int) population) *
04379 m->triangles.itembytes));
04380 if (!deadtri(sampletri.tri)) {
04381 org(sampletri, torg);
04382 dist = (searchpoint[0] - torg[0]) * (searchpoint[0] - torg[0]) +
04383 (searchpoint[1] - torg[1]) * (searchpoint[1] - torg[1]);
04384 if (dist < searchdist) {
04385 otricopy(sampletri, *searchtri);
04386 searchdist = dist;
04387 if (b->verbose > 2) {
04388 printf(" Choosing triangle with origin (%.12g, %.12g).\n",
04389 torg[0], torg[1]);
04390 }
04391 }
04392 }
04393
04394 samplesleft--;
04395 totalsamplesleft--;
04396 } while ((samplesleft > 0) && (totalsamplesleft > 0));
04397
04398 if (totalsamplesleft > 0) {
04399 sampleblock = (int **) *sampleblock;
04400 samplesleft = samplesperblock;
04401 totalpopulation -= population;
04402 population = TRIPERBLOCK;
04403 }
04404 }
04405
04406
04407 org(*searchtri, torg);
04408 dest(*searchtri, tdest);
04409
04410 if ((torg[0] == searchpoint[0]) && (torg[1] == searchpoint[1])) {
04411 return ONVERTEX;
04412 }
04413 if ((tdest[0] == searchpoint[0]) && (tdest[1] == searchpoint[1])) {
04414 lnextself(*searchtri);
04415 return ONVERTEX;
04416 }
04417
04418 ahead = counterclockwise(m, b, torg, tdest, searchpoint);
04419 if (ahead < 0.0) {
04420
04421
04422 symself(*searchtri);
04423 } else if (ahead == 0.0) {
04424
04425 if (((torg[0] < searchpoint[0]) == (searchpoint[0] < tdest[0])) &&
04426 ((torg[1] < searchpoint[1]) == (searchpoint[1] < tdest[1]))) {
04427 return ONEDGE;
04428 }
04429 }
04430 return preciselocate(m, b, searchpoint, searchtri, 0);
04431 }
04432
04435
04436
04437
04441
04442
04443
04444
04445
04446
04447
04448
04449
04450
04451
04452 void insertsubseg(struct mesh *m, struct behavior *b, struct otri *tri,
04453 int subsegmark)
04454 {
04455 struct otri oppotri;
04456 struct osub newsubseg;
04457 vertex triorg, tridest;
04458 triangle ptr;
04459 subseg sptr;
04460
04461 org(*tri, triorg);
04462 dest(*tri, tridest);
04463
04464 if (vertexmark(triorg) == 0) {
04465 setvertexmark(triorg, subsegmark);
04466 }
04467 if (vertexmark(tridest) == 0) {
04468 setvertexmark(tridest, subsegmark);
04469 }
04470
04471 tspivot(*tri, newsubseg);
04472 if (newsubseg.ss == m->dummysub) {
04473
04474 makesubseg(m, &newsubseg);
04475 setsorg(newsubseg, tridest);
04476 setsdest(newsubseg, triorg);
04477 setsegorg(newsubseg, tridest);
04478 setsegdest(newsubseg, triorg);
04479
04480
04481
04482
04483 tsbond(*tri, newsubseg);
04484 sym(*tri, oppotri);
04485 ssymself(newsubseg);
04486 tsbond(oppotri, newsubseg);
04487 setmark(newsubseg, subsegmark);
04488 if (b->verbose > 2) {
04489 printf(" Inserting new ");
04490 printsubseg(m, b, &newsubseg);
04491 }
04492 } else {
04493 if (mark(newsubseg) == 0) {
04494 setmark(newsubseg, subsegmark);
04495 }
04496 }
04497 }
04498
04499
04500
04501
04502
04503
04504
04505
04506
04507
04508
04509
04510
04511
04512
04513
04514
04515
04516
04517
04518
04519
04520
04521
04522
04523
04524
04525
04526
04527
04528
04529
04530
04531
04532
04533
04534
04535
04536
04537
04538
04539
04540
04541
04542
04543
04544
04545
04546
04547 void flip(struct mesh *m, struct behavior *b, struct otri *flipedge)
04548 {
04549 struct otri botleft, botright;
04550 struct otri topleft, topright;
04551 struct otri top;
04552 struct otri botlcasing, botrcasing;
04553 struct otri toplcasing, toprcasing;
04554 struct osub botlsubseg, botrsubseg;
04555 struct osub toplsubseg, toprsubseg;
04556 vertex leftvertex, rightvertex, botvertex;
04557 vertex farvertex;
04558 triangle ptr;
04559 subseg sptr;
04560
04561
04562 org(*flipedge, rightvertex);
04563 dest(*flipedge, leftvertex);
04564 apex(*flipedge, botvertex);
04565 sym(*flipedge, top);
04566 apex(top, farvertex);
04567
04568
04569 lprev(top, topleft);
04570 sym(topleft, toplcasing);
04571 lnext(top, topright);
04572 sym(topright, toprcasing);
04573 lnext(*flipedge, botleft);
04574 sym(botleft, botlcasing);
04575 lprev(*flipedge, botright);
04576 sym(botright, botrcasing);
04577
04578 bond(topleft, botlcasing);
04579 bond(botleft, botrcasing);
04580 bond(botright, toprcasing);
04581 bond(topright, toplcasing);
04582
04583 if (m->checksegments) {
04584
04585 tspivot(topleft, toplsubseg);
04586 tspivot(botleft, botlsubseg);
04587 tspivot(botright, botrsubseg);
04588 tspivot(topright, toprsubseg);
04589 if (toplsubseg.ss == m->dummysub) {
04590 tsdissolve(topright);
04591 } else {
04592 tsbond(topright, toplsubseg);
04593 }
04594 if (botlsubseg.ss == m->dummysub) {
04595 tsdissolve(topleft);
04596 } else {
04597 tsbond(topleft, botlsubseg);
04598 }
04599 if (botrsubseg.ss == m->dummysub) {
04600 tsdissolve(botleft);
04601 } else {
04602 tsbond(botleft, botrsubseg);
04603 }
04604 if (toprsubseg.ss == m->dummysub) {
04605 tsdissolve(botright);
04606 } else {
04607 tsbond(botright, toprsubseg);
04608 }
04609 }
04610
04611
04612 setorg(*flipedge, farvertex);
04613 setdest(*flipedge, botvertex);
04614 setapex(*flipedge, rightvertex);
04615 setorg(top, botvertex);
04616 setdest(top, farvertex);
04617 setapex(top, leftvertex);
04618 if (b->verbose > 2) {
04619 printf(" Edge flip results in left ");
04620 printtriangle(m, b, &top);
04621 printf(" and right ");
04622 printtriangle(m, b, flipedge);
04623 }
04624 }
04625
04626
04627
04628
04629
04630
04631
04632
04633
04634
04635
04636
04637
04638
04639
04640
04641
04642
04643
04644
04645
04646
04647
04648
04649
04650
04651
04652
04653
04654
04655
04656
04657
04658
04659 void unflip(struct mesh *m, struct behavior *b, struct otri *flipedge)
04660 {
04661 struct otri botleft, botright;
04662 struct otri topleft, topright;
04663 struct otri top;
04664 struct otri botlcasing, botrcasing;
04665 struct otri toplcasing, toprcasing;
04666 struct osub botlsubseg, botrsubseg;
04667 struct osub toplsubseg, toprsubseg;
04668 vertex leftvertex, rightvertex, botvertex;
04669 vertex farvertex;
04670 triangle ptr;
04671 subseg sptr;
04672
04673
04674 org(*flipedge, rightvertex);
04675 dest(*flipedge, leftvertex);
04676 apex(*flipedge, botvertex);
04677 sym(*flipedge, top);
04678 apex(top, farvertex);
04679
04680
04681 lprev(top, topleft);
04682 sym(topleft, toplcasing);
04683 lnext(top, topright);
04684 sym(topright, toprcasing);
04685 lnext(*flipedge, botleft);
04686 sym(botleft, botlcasing);
04687 lprev(*flipedge, botright);
04688 sym(botright, botrcasing);
04689
04690 bond(topleft, toprcasing);
04691 bond(botleft, toplcasing);
04692 bond(botright, botlcasing);
04693 bond(topright, botrcasing);
04694
04695 if (m->checksegments) {
04696
04697 tspivot(topleft, toplsubseg);
04698 tspivot(botleft, botlsubseg);
04699 tspivot(botright, botrsubseg);
04700 tspivot(topright, toprsubseg);
04701 if (toplsubseg.ss == m->dummysub) {
04702 tsdissolve(botleft);
04703 } else {
04704 tsbond(botleft, toplsubseg);
04705 }
04706 if (botlsubseg.ss == m->dummysub) {
04707 tsdissolve(botright);
04708 } else {
04709 tsbond(botright, botlsubseg);
04710 }
04711 if (botrsubseg.ss == m->dummysub) {
04712 tsdissolve(topright);
04713 } else {
04714 tsbond(topright, botrsubseg);
04715 }
04716 if (toprsubseg.ss == m->dummysub) {
04717 tsdissolve(topleft);
04718 } else {
04719 tsbond(topleft, toprsubseg);
04720 }
04721 }
04722
04723
04724 setorg(*flipedge, botvertex);
04725 setdest(*flipedge, farvertex);
04726 setapex(*flipedge, leftvertex);
04727 setorg(top, farvertex);
04728 setdest(top, botvertex);
04729 setapex(top, rightvertex);
04730 if (b->verbose > 2) {
04731 printf(" Edge unflip results in left ");
04732 printtriangle(m, b, flipedge);
04733 printf(" and right ");
04734 printtriangle(m, b, &top);
04735 }
04736 }
04737
04738
04739
04740
04741
04742
04743
04744
04745
04746
04747
04748
04749
04750
04751
04752
04753
04754
04755
04756
04757
04758
04759
04760
04761
04762
04763
04764
04765
04766
04767
04768
04769
04770
04771
04772
04773
04774
04775
04776
04777
04778
04779
04780
04781
04782
04783
04784
04785 enum insertvertexresult insertvertex(struct mesh *m, struct behavior *b,
04786 vertex newvertex, struct otri *searchtri,
04787 struct osub *splitseg,
04788 int segmentflaws, int triflaws)
04789 {
04790 struct otri horiz;
04791 struct otri top;
04792 struct otri botleft, botright;
04793 struct otri topleft, topright;
04794 struct otri newbotleft, newbotright;
04795 struct otri newtopright;
04796 struct otri botlcasing, botrcasing;
04797 struct otri toplcasing, toprcasing;
04798 struct otri testtri;
04799 struct osub botlsubseg, botrsubseg;
04800 struct osub toplsubseg, toprsubseg;
04801 struct osub brokensubseg;
04802 struct osub checksubseg;
04803 struct osub rightsubseg;
04804 struct osub newsubseg;
04805 struct badsubseg *encroached;
04806 struct flipstacker *newflip;
04807 vertex first;
04808 vertex leftvertex, rightvertex, botvertex, topvertex, farvertex;
04809 vertex segmentorg, segmentdest;
04810 float attrib;
04811 float area;
04812 enum insertvertexresult success;
04813 enum locateresult intersect;
04814 int doflip;
04815 int mirrorflag;
04816 int enq;
04817 int i;
04818 triangle ptr;
04819 subseg sptr;
04820
04821 if (b->verbose > 1) {
04822 printf(" Inserting (%.12g, %.12g).\n", newvertex[0], newvertex[1]);
04823 }
04824
04825 if (splitseg == (struct osub *) NULL) {
04826
04827
04828 if (searchtri->tri == m->dummytri) {
04829
04830 horiz.tri = m->dummytri;
04831 horiz.orient = 0;
04832 symself(horiz);
04833
04834 intersect = locate(m, b, newvertex, &horiz);
04835 } else {
04836
04837 otricopy(*searchtri, horiz);
04838 intersect = preciselocate(m, b, newvertex, &horiz, 1);
04839 }
04840 } else {
04841
04842
04843 otricopy(*searchtri, horiz);
04844 intersect = ONEDGE;
04845 }
04846
04847 if (intersect == ONVERTEX) {
04848
04849
04850 otricopy(horiz, *searchtri);
04851 otricopy(horiz, m->recenttri);
04852 return DUPLICATEVERTEX;
04853 }
04854 if ((intersect == ONEDGE) || (intersect == OUTSIDE)) {
04855
04856 if (m->checksegments && (splitseg == (struct osub *) NULL)) {
04857
04858 tspivot(horiz, brokensubseg);
04859 if (brokensubseg.ss != m->dummysub) {
04860
04861 if (segmentflaws) {
04862 enq = b->nobisect != 2;
04863 if (enq && (b->nobisect == 1)) {
04864
04865
04866 sym(horiz, testtri);
04867 enq = testtri.tri != m->dummytri;
04868 }
04869 if (enq) {
04870
04871 encroached = (struct badsubseg *) poolalloc(&m->badsubsegs);
04872 encroached->encsubseg = sencode(brokensubseg);
04873 sorg(brokensubseg, encroached->subsegorg);
04874 sdest(brokensubseg, encroached->subsegdest);
04875 if (b->verbose > 2) {
04876 printf(
04877 " Queueing encroached subsegment (%.12g, %.12g) (%.12g, %.12g).\n",
04878 encroached->subsegorg[0], encroached->subsegorg[1],
04879 encroached->subsegdest[0], encroached->subsegdest[1]);
04880 }
04881 }
04882 }
04883
04884
04885 otricopy(horiz, *searchtri);
04886 otricopy(horiz, m->recenttri);
04887 return VIOLATINGVERTEX;
04888 }
04889 }
04890
04891
04892
04893 lprev(horiz, botright);
04894 sym(botright, botrcasing);
04895 sym(horiz, topright);
04896
04897 mirrorflag = topright.tri != m->dummytri;
04898 if (mirrorflag) {
04899 lnextself(topright);
04900 sym(topright, toprcasing);
04901 maketriangle(m, b, &newtopright);
04902 } else {
04903
04904 m->hullsize++;
04905 }
04906 maketriangle(m, b, &newbotright);
04907
04908
04909 org(horiz, rightvertex);
04910 dest(horiz, leftvertex);
04911 apex(horiz, botvertex);
04912 setorg(newbotright, botvertex);
04913 setdest(newbotright, rightvertex);
04914 setapex(newbotright, newvertex);
04915 setorg(horiz, newvertex);
04916 for (i = 0; i < m->eextras; i++) {
04917
04918 setelemattribute(newbotright, i, elemattribute(botright, i));
04919 }
04920 if (b->vararea) {
04921
04922 setareabound(newbotright, areabound(botright));
04923 }
04924 if (mirrorflag) {
04925 dest(topright, topvertex);
04926 setorg(newtopright, rightvertex);
04927 setdest(newtopright, topvertex);
04928 setapex(newtopright, newvertex);
04929 setorg(topright, newvertex);
04930 for (i = 0; i < m->eextras; i++) {
04931
04932 setelemattribute(newtopright, i, elemattribute(topright, i));
04933 }
04934 if (b->vararea) {
04935
04936 setareabound(newtopright, areabound(topright));
04937 }
04938 }
04939
04940
04941
04942 if (m->checksegments) {
04943 tspivot(botright, botrsubseg);
04944 if (botrsubseg.ss != m->dummysub) {
04945 tsdissolve(botright);
04946 tsbond(newbotright, botrsubseg);
04947 }
04948 if (mirrorflag) {
04949 tspivot(topright, toprsubseg);
04950 if (toprsubseg.ss != m->dummysub) {
04951 tsdissolve(topright);
04952 tsbond(newtopright, toprsubseg);
04953 }
04954 }
04955 }
04956
04957
04958 bond(newbotright, botrcasing);
04959 lprevself(newbotright);
04960 bond(newbotright, botright);
04961 lprevself(newbotright);
04962 if (mirrorflag) {
04963 bond(newtopright, toprcasing);
04964 lnextself(newtopright);
04965 bond(newtopright, topright);
04966 lnextself(newtopright);
04967 bond(newtopright, newbotright);
04968 }
04969
04970 if (splitseg != (struct osub *) NULL) {
04971
04972 setsdest(*splitseg, newvertex);
04973 segorg(*splitseg, segmentorg);
04974 segdest(*splitseg, segmentdest);
04975 ssymself(*splitseg);
04976 spivot(*splitseg, rightsubseg);
04977 insertsubseg(m, b, &newbotright, mark(*splitseg));
04978 tspivot(newbotright, newsubseg);
04979 setsegorg(newsubseg, segmentorg);
04980 setsegdest(newsubseg, segmentdest);
04981 sbond(*splitseg, newsubseg);
04982 ssymself(newsubseg);
04983 sbond(newsubseg, rightsubseg);
04984 ssymself(*splitseg);
04985
04986
04987 if (vertexmark(newvertex) == 0) {
04988 setvertexmark(newvertex, mark(*splitseg));
04989 }
04990 }
04991
04992 if (m->checkquality) {
04993 poolrestart(&m->flipstackers);
04994 m->lastflip = (struct flipstacker *) poolalloc(&m->flipstackers);
04995 m->lastflip->flippedtri = encode(horiz);
04996 m->lastflip->prevflip = (struct flipstacker *) &insertvertex;
04997 }
04998 if (b->verbose > 2) {
04999 printf(" Updating bottom left ");
05000 printtriangle(m, b, &botright);
05001 if (mirrorflag) {
05002 printf(" Updating top left ");
05003 printtriangle(m, b, &topright);
05004 printf(" Creating top right ");
05005 printtriangle(m, b, &newtopright);
05006 }
05007 printf(" Creating bottom right ");
05008 printtriangle(m, b, &newbotright);
05009 }
05010
05011
05012
05013 lnextself(horiz);
05014 } else {
05015
05016 lnext(horiz, botleft);
05017 lprev(horiz, botright);
05018 sym(botleft, botlcasing);
05019 sym(botright, botrcasing);
05020 maketriangle(m, b, &newbotleft);
05021 maketriangle(m, b, &newbotright);
05022
05023
05024 org(horiz, rightvertex);
05025 dest(horiz, leftvertex);
05026 apex(horiz, botvertex);
05027 setorg(newbotleft, leftvertex);
05028 setdest(newbotleft, botvertex);
05029 setapex(newbotleft, newvertex);
05030 setorg(newbotright, botvertex);
05031 setdest(newbotright, rightvertex);
05032 setapex(newbotright, newvertex);
05033 setapex(horiz, newvertex);
05034 for (i = 0; i < m->eextras; i++) {
05035
05036 attrib = elemattribute(horiz, i);
05037 setelemattribute(newbotleft, i, attrib);
05038 setelemattribute(newbotright, i, attrib);
05039 }
05040 if (b->vararea) {
05041
05042 area = areabound(horiz);
05043 setareabound(newbotleft, area);
05044 setareabound(newbotright, area);
05045 }
05046
05047
05048
05049 if (m->checksegments) {
05050 tspivot(botleft, botlsubseg);
05051 if (botlsubseg.ss != m->dummysub) {
05052 tsdissolve(botleft);
05053 tsbond(newbotleft, botlsubseg);
05054 }
05055 tspivot(botright, botrsubseg);
05056 if (botrsubseg.ss != m->dummysub) {
05057 tsdissolve(botright);
05058 tsbond(newbotright, botrsubseg);
05059 }
05060 }
05061
05062
05063 bond(newbotleft, botlcasing);
05064 bond(newbotright, botrcasing);
05065 lnextself(newbotleft);
05066 lprevself(newbotright);
05067 bond(newbotleft, newbotright);
05068 lnextself(newbotleft);
05069 bond(botleft, newbotleft);
05070 lprevself(newbotright);
05071 bond(botright, newbotright);
05072
05073 if (m->checkquality) {
05074 poolrestart(&m->flipstackers);
05075 m->lastflip = (struct flipstacker *) poolalloc(&m->flipstackers);
05076 m->lastflip->flippedtri = encode(horiz);
05077 m->lastflip->prevflip = (struct flipstacker *) NULL;
05078 }
05079 if (b->verbose > 2) {
05080 printf(" Updating top ");
05081 printtriangle(m, b, &horiz);
05082 printf(" Creating left ");
05083 printtriangle(m, b, &newbotleft);
05084 printf(" Creating right ");
05085 printtriangle(m, b, &newbotright);
05086 }
05087 }
05088
05089
05090
05091 success = SUCCESSFULVERTEX;
05092
05093
05094
05095
05096 org(horiz, first);
05097 rightvertex = first;
05098 dest(horiz, leftvertex);
05099
05100 while (1) {
05101
05102 doflip = 1;
05103
05104 if (m->checksegments) {
05105
05106 tspivot(horiz, checksubseg);
05107 if (checksubseg.ss != m->dummysub) {
05108
05109 doflip = 0;
05110 }
05111 }
05112
05113 if (doflip) {
05114
05115 sym(horiz, top);
05116 if (top.tri == m->dummytri) {
05117
05118 doflip = 0;
05119 } else {
05120
05121 apex(top, farvertex);
05122
05123
05124
05125
05126
05127 if ((leftvertex == m->infvertex1) || (leftvertex == m->infvertex2) ||
05128 (leftvertex == m->infvertex3)) {
05129
05130
05131
05132
05133 doflip = counterclockwise(m, b, newvertex, rightvertex, farvertex)
05134 > 0.0;
05135 } else if ((rightvertex == m->infvertex1) ||
05136 (rightvertex == m->infvertex2) ||
05137 (rightvertex == m->infvertex3)) {
05138
05139
05140
05141
05142 doflip = counterclockwise(m, b, farvertex, leftvertex, newvertex)
05143 > 0.0;
05144 } else if ((farvertex == m->infvertex1) ||
05145 (farvertex == m->infvertex2) ||
05146 (farvertex == m->infvertex3)) {
05147
05148
05149 doflip = 0;
05150 } else {
05151
05152 doflip = incircle(m, b, leftvertex, newvertex, rightvertex,
05153 farvertex) > 0.0;
05154 }
05155 if (doflip) {
05156
05157
05158
05159 lprev(top, topleft);
05160 sym(topleft, toplcasing);
05161 lnext(top, topright);
05162 sym(topright, toprcasing);
05163 lnext(horiz, botleft);
05164 sym(botleft, botlcasing);
05165 lprev(horiz, botright);
05166 sym(botright, botrcasing);
05167
05168 bond(topleft, botlcasing);
05169 bond(botleft, botrcasing);
05170 bond(botright, toprcasing);
05171 bond(topright, toplcasing);
05172 if (m->checksegments) {
05173
05174 tspivot(topleft, toplsubseg);
05175 tspivot(botleft, botlsubseg);
05176 tspivot(botright, botrsubseg);
05177 tspivot(topright, toprsubseg);
05178 if (toplsubseg.ss == m->dummysub) {
05179 tsdissolve(topright);
05180 } else {
05181 tsbond(topright, toplsubseg);
05182 }
05183 if (botlsubseg.ss == m->dummysub) {
05184 tsdissolve(topleft);
05185 } else {
05186 tsbond(topleft, botlsubseg);
05187 }
05188 if (botrsubseg.ss == m->dummysub) {
05189 tsdissolve(botleft);
05190 } else {
05191 tsbond(botleft, botrsubseg);
05192 }
05193 if (toprsubseg.ss == m->dummysub) {
05194 tsdissolve(botright);
05195 } else {
05196 tsbond(botright, toprsubseg);
05197 }
05198 }
05199
05200 setorg(horiz, farvertex);
05201 setdest(horiz, newvertex);
05202 setapex(horiz, rightvertex);
05203 setorg(top, newvertex);
05204 setdest(top, farvertex);
05205 setapex(top, leftvertex);
05206 for (i = 0; i < m->eextras; i++) {
05207
05208 attrib = 0.5 * (elemattribute(top, i) + elemattribute(horiz, i));
05209 setelemattribute(top, i, attrib);
05210 setelemattribute(horiz, i, attrib);
05211 }
05212 if (b->vararea) {
05213 if ((areabound(top) <= 0.0) || (areabound(horiz) <= 0.0)) {
05214 area = -1.0;
05215 } else {
05216
05217
05218
05219 area = 0.5 * (areabound(top) + areabound(horiz));
05220 }
05221 setareabound(top, area);
05222 setareabound(horiz, area);
05223 }
05224
05225 if (m->checkquality) {
05226 newflip = (struct flipstacker *) poolalloc(&m->flipstackers);
05227 newflip->flippedtri = encode(horiz);
05228 newflip->prevflip = m->lastflip;
05229 m->lastflip = newflip;
05230 }
05231 if (b->verbose > 2) {
05232 printf(" Edge flip results in left ");
05233 lnextself(topleft);
05234 printtriangle(m, b, &topleft);
05235 printf(" and right ");
05236 printtriangle(m, b, &horiz);
05237 }
05238
05239
05240
05241 lprevself(horiz);
05242 leftvertex = farvertex;
05243 }
05244 }
05245 }
05246 if (!doflip) {
05247
05248
05249 lnextself(horiz);
05250 sym(horiz, testtri);
05251
05252
05253
05254 if ((leftvertex == first) || (testtri.tri == m->dummytri)) {
05255
05256 lnext(horiz, *searchtri);
05257 lnext(horiz, m->recenttri);
05258 return success;
05259 }
05260
05261 lnext(testtri, horiz);
05262 rightvertex = leftvertex;
05263 dest(horiz, leftvertex);
05264 }
05265 }
05266 }
05267
05268
05269
05270
05271
05272
05273
05274
05275
05276
05277
05278
05279
05280
05281
05282
05283
05284
05285
05286
05287
05288
05289
05290
05291
05292
05293
05294
05295
05296
05297
05298
05299
05300
05301
05302
05303
05304
05305
05306
05307
05308
05309
05310
05311
05312
05313
05314
05315
05316
05317
05318
05319
05320
05321
05322
05323
05324
05325
05326
05327
05328
05329
05330
05331
05332 void triangulatepolygon(struct mesh *m, struct behavior *b,
05333 struct otri *firstedge, struct otri *lastedge,
05334 int edgecount, int doflip, int triflaws)
05335 {
05336 struct otri testtri;
05337 struct otri besttri;
05338 struct otri tempedge;
05339 vertex leftbasevertex, rightbasevertex;
05340 vertex testvertex;
05341 vertex bestvertex;
05342 int bestnumber;
05343 int i;
05344 triangle ptr;
05345
05346
05347 apex(*lastedge, leftbasevertex);
05348 dest(*firstedge, rightbasevertex);
05349 if (b->verbose > 2) {
05350 printf(" Triangulating interior polygon at edge\n");
05351 printf(" (%.12g, %.12g) (%.12g, %.12g)\n", leftbasevertex[0],
05352 leftbasevertex[1], rightbasevertex[0], rightbasevertex[1]);
05353 }
05354
05355 onext(*firstedge, besttri);
05356 dest(besttri, bestvertex);
05357 otricopy(besttri, testtri);
05358 bestnumber = 1;
05359 for (i = 2; i <= edgecount - 2; i++) {
05360 onextself(testtri);
05361 dest(testtri, testvertex);
05362
05363 if (incircle(m, b, leftbasevertex, rightbasevertex, bestvertex,
05364 testvertex) > 0.0) {
05365 otricopy(testtri, besttri);
05366 bestvertex = testvertex;
05367 bestnumber = i;
05368 }
05369 }
05370 if (b->verbose > 2) {
05371 printf(" Connecting edge to (%.12g, %.12g)\n", bestvertex[0],
05372 bestvertex[1]);
05373 }
05374 if (bestnumber > 1) {
05375
05376 oprev(besttri, tempedge);
05377 triangulatepolygon(m, b, firstedge, &tempedge, bestnumber + 1, 1,
05378 triflaws);
05379 }
05380 if (bestnumber < edgecount - 2) {
05381
05382 sym(besttri, tempedge);
05383 triangulatepolygon(m, b, &besttri, lastedge, edgecount - bestnumber, 1,
05384 triflaws);
05385
05386 sym(tempedge, besttri);
05387 }
05388 if (doflip) {
05389
05390 flip(m, b, &besttri);
05391 }
05392
05393 otricopy(besttri, *lastedge);
05394 }
05395
05398
05399
05400
05404
05405
05406
05407
05408
05409
05410
05411
05412
05413
05414
05415
05416
05417
05418
05419
05420
05421
05422
05423
05424
05425
05426
05427
05428
05429
05430
05431
05432
05433
05434
05435
05436
05437
05438
05439
05440
05441
05442
05443
05444
05445 void vertexsort(vertex *sortarray, int arraysize)
05446 {
05447 int left, right;
05448 int pivot;
05449 float pivotx, pivoty;
05450 vertex temp;
05451
05452 if (arraysize == 2) {
05453
05454 if ((sortarray[0][0] > sortarray[1][0]) ||
05455 ((sortarray[0][0] == sortarray[1][0]) &&
05456 (sortarray[0][1] > sortarray[1][1]))) {
05457 temp = sortarray[1];
05458 sortarray[1] = sortarray[0];
05459 sortarray[0] = temp;
05460 }
05461 return;
05462 }
05463
05464 pivot = (int) randomnation((unsigned int) arraysize);
05465 pivotx = sortarray[pivot][0];
05466 pivoty = sortarray[pivot][1];
05467
05468 left = -1;
05469 right = arraysize;
05470 while (left < right) {
05471
05472 do {
05473 left++;
05474 } while ((left <= right) && ((sortarray[left][0] < pivotx) ||
05475 ((sortarray[left][0] == pivotx) &&
05476 (sortarray[left][1] < pivoty))));
05477
05478 do {
05479 right--;
05480 } while ((left <= right) && ((sortarray[right][0] > pivotx) ||
05481 ((sortarray[right][0] == pivotx) &&
05482 (sortarray[right][1] > pivoty))));
05483 if (left < right) {
05484
05485 temp = sortarray[left];
05486 sortarray[left] = sortarray[right];
05487 sortarray[right] = temp;
05488 }
05489 }
05490 if (left > 1) {
05491
05492 vertexsort(sortarray, left);
05493 }
05494 if (right < arraysize - 2) {
05495
05496 vertexsort(&sortarray[right + 1], arraysize - right - 1);
05497 }
05498 }
05499
05500
05501
05502
05503
05504
05505
05506
05507
05508
05509
05510
05511
05512 void vertexmedian(vertex *sortarray, int arraysize, int median, int axis)
05513 {
05514 int left, right;
05515 int pivot;
05516 float pivot1, pivot2;
05517 vertex temp;
05518
05519 if (arraysize == 2) {
05520
05521 if ((sortarray[0][axis] > sortarray[1][axis]) ||
05522 ((sortarray[0][axis] == sortarray[1][axis]) &&
05523 (sortarray[0][1 - axis] > sortarray[1][1 - axis]))) {
05524 temp = sortarray[1];
05525 sortarray[1] = sortarray[0];
05526 sortarray[0] = temp;
05527 }
05528 return;
05529 }
05530
05531 pivot = (int) randomnation((unsigned int) arraysize);
05532 pivot1 = sortarray[pivot][axis];
05533 pivot2 = sortarray[pivot][1 - axis];
05534
05535 left = -1;
05536 right = arraysize;
05537 while (left < right) {
05538
05539 do {
05540 left++;
05541 } while ((left <= right) && ((sortarray[left][axis] < pivot1) ||
05542 ((sortarray[left][axis] == pivot1) &&
05543 (sortarray[left][1 - axis] < pivot2))));
05544
05545 do {
05546 right--;
05547 } while ((left <= right) && ((sortarray[right][axis] > pivot1) ||
05548 ((sortarray[right][axis] == pivot1) &&
05549 (sortarray[right][1 - axis] > pivot2))));
05550 if (left < right) {
05551
05552 temp = sortarray[left];
05553 sortarray[left] = sortarray[right];
05554 sortarray[right] = temp;
05555 }
05556 }
05557
05558
05559 if (left > median) {
05560
05561 vertexmedian(sortarray, left, median, axis);
05562 }
05563 if (right < median - 1) {
05564
05565 vertexmedian(&sortarray[right + 1], arraysize - right - 1,
05566 median - right - 1, axis);
05567 }
05568 }
05569
05570
05571
05572
05573
05574
05575
05576
05577
05578
05579
05580
05581 void alternateaxes(vertex *sortarray, int arraysize, int axis)
05582 {
05583 int divider;
05584
05585 divider = arraysize >> 1;
05586 if (arraysize <= 3) {
05587
05588
05589 axis = 0;
05590 }
05591
05592 vertexmedian(sortarray, arraysize, divider, axis);
05593
05594 if (arraysize - divider >= 2) {
05595 if (divider >= 2) {
05596 alternateaxes(sortarray, divider, 1 - axis);
05597 }
05598 alternateaxes(&sortarray[divider], arraysize - divider, 1 - axis);
05599 }
05600 }
05601
05602
05603
05604
05605
05606
05607
05608
05609
05610
05611
05612
05613
05614
05615
05616
05617
05618
05619
05620
05621
05622
05623
05624
05625
05626
05627
05628
05629
05630
05631
05632
05633
05634
05635
05636
05637 void mergehulls(struct mesh *m, struct behavior *b, struct otri *farleft,
05638 struct otri *innerleft, struct otri *innerright,
05639 struct otri *farright, int axis)
05640 {
05641 struct otri leftcand, rightcand;
05642 struct otri baseedge;
05643 struct otri nextedge;
05644 struct otri sidecasing, topcasing, outercasing;
05645 struct otri checkedge;
05646 vertex innerleftdest;
05647 vertex innerrightorg;
05648 vertex innerleftapex, innerrightapex;
05649 vertex farleftpt, farrightpt;
05650 vertex farleftapex, farrightapex;
05651 vertex lowerleft, lowerright;
05652 vertex upperleft, upperright;
05653 vertex nextapex;
05654 vertex checkvertex;
05655 int changemade;
05656 int badedge;
05657 int leftfinished, rightfinished;
05658 triangle ptr;
05659
05660 dest(*innerleft, innerleftdest);
05661 apex(*innerleft, innerleftapex);
05662 org(*innerright, innerrightorg);
05663 apex(*innerright, innerrightapex);
05664
05665 if (b->dwyer && (axis == 1)) {
05666 org(*farleft, farleftpt);
05667 apex(*farleft, farleftapex);
05668 dest(*farright, farrightpt);
05669 apex(*farright, farrightapex);
05670
05671
05672
05673 while (farleftapex[1] < farleftpt[1]) {
05674 lnextself(*farleft);
05675 symself(*farleft);
05676 farleftpt = farleftapex;
05677 apex(*farleft, farleftapex);
05678 }
05679 sym(*innerleft, checkedge);
05680 apex(checkedge, checkvertex);
05681 while (checkvertex[1] > innerleftdest[1]) {
05682 lnext(checkedge, *innerleft);
05683 innerleftapex = innerleftdest;
05684 innerleftdest = checkvertex;
05685 sym(*innerleft, checkedge);
05686 apex(checkedge, checkvertex);
05687 }
05688 while (innerrightapex[1] < innerrightorg[1]) {
05689 lnextself(*innerright);
05690 symself(*innerright);
05691 innerrightorg = innerrightapex;
05692 apex(*innerright, innerrightapex);
05693 }
05694 sym(*farright, checkedge);
05695 apex(checkedge, checkvertex);
05696 while (checkvertex[1] > farrightpt[1]) {
05697 lnext(checkedge, *farright);
05698 farrightapex = farrightpt;
05699 farrightpt = checkvertex;
05700 sym(*farright, checkedge);
05701 apex(checkedge, checkvertex);
05702 }
05703 }
05704
05705 do {
05706 changemade = 0;
05707
05708 if (counterclockwise(m, b, innerleftdest, innerleftapex, innerrightorg) >
05709 0.0) {
05710 lprevself(*innerleft);
05711 symself(*innerleft);
05712 innerleftdest = innerleftapex;
05713 apex(*innerleft, innerleftapex);
05714 changemade = 1;
05715 }
05716
05717 if (counterclockwise(m, b, innerrightapex, innerrightorg, innerleftdest) >
05718 0.0) {
05719 lnextself(*innerright);
05720 symself(*innerright);
05721 innerrightorg = innerrightapex;
05722 apex(*innerright, innerrightapex);
05723 changemade = 1;
05724 }
05725 } while (changemade);
05726
05727 sym(*innerleft, leftcand);
05728 sym(*innerright, rightcand);
05729
05730 maketriangle(m, b, &baseedge);
05731
05732 bond(baseedge, *innerleft);
05733 lnextself(baseedge);
05734 bond(baseedge, *innerright);
05735 lnextself(baseedge);
05736 setorg(baseedge, innerrightorg);
05737 setdest(baseedge, innerleftdest);
05738
05739 if (b->verbose > 2) {
05740 printf(" Creating base bounding ");
05741 printtriangle(m, b, &baseedge);
05742 }
05743
05744 org(*farleft, farleftpt);
05745 if (innerleftdest == farleftpt) {
05746 lnext(baseedge, *farleft);
05747 }
05748 dest(*farright, farrightpt);
05749 if (innerrightorg == farrightpt) {
05750 lprev(baseedge, *farright);
05751 }
05752
05753 lowerleft = innerleftdest;
05754 lowerright = innerrightorg;
05755
05756 apex(leftcand, upperleft);
05757 apex(rightcand, upperright);
05758
05759 while (1) {
05760
05761
05762
05763
05764 leftfinished = counterclockwise(m, b, upperleft, lowerleft, lowerright) <=
05765 0.0;
05766 rightfinished = counterclockwise(m, b, upperright, lowerleft, lowerright)
05767 <= 0.0;
05768 if (leftfinished && rightfinished) {
05769
05770 maketriangle(m, b, &nextedge);
05771 setorg(nextedge, lowerleft);
05772 setdest(nextedge, lowerright);
05773
05774
05775 bond(nextedge, baseedge);
05776 lnextself(nextedge);
05777 bond(nextedge, rightcand);
05778 lnextself(nextedge);
05779 bond(nextedge, leftcand);
05780 if (b->verbose > 2) {
05781 printf(" Creating top bounding ");
05782 printtriangle(m, b, &nextedge);
05783 }
05784
05785 if (b->dwyer && (axis == 1)) {
05786 org(*farleft, farleftpt);
05787 apex(*farleft, farleftapex);
05788 dest(*farright, farrightpt);
05789 apex(*farright, farrightapex);
05790 sym(*farleft, checkedge);
05791 apex(checkedge, checkvertex);
05792
05793
05794
05795 while (checkvertex[0] < farleftpt[0]) {
05796 lprev(checkedge, *farleft);
05797 farleftapex = farleftpt;
05798 farleftpt = checkvertex;
05799 sym(*farleft, checkedge);
05800 apex(checkedge, checkvertex);
05801 }
05802 while (farrightapex[0] > farrightpt[0]) {
05803 lprevself(*farright);
05804 symself(*farright);
05805 farrightpt = farrightapex;
05806 apex(*farright, farrightapex);
05807 }
05808 }
05809 return;
05810 }
05811
05812 if (!leftfinished) {
05813
05814 lprev(leftcand, nextedge);
05815 symself(nextedge);
05816 apex(nextedge, nextapex);
05817
05818
05819 if (nextapex != (vertex) NULL) {
05820
05821 badedge = incircle(m, b, lowerleft, lowerright, upperleft, nextapex) >
05822 0.0;
05823 while (badedge) {
05824
05825
05826 lnextself(nextedge);
05827 sym(nextedge, topcasing);
05828 lnextself(nextedge);
05829 sym(nextedge, sidecasing);
05830 bond(nextedge, topcasing);
05831 bond(leftcand, sidecasing);
05832 lnextself(leftcand);
05833 sym(leftcand, outercasing);
05834 lprevself(nextedge);
05835 bond(nextedge, outercasing);
05836
05837 setorg(leftcand, lowerleft);
05838 setdest(leftcand, NULL);
05839 setapex(leftcand, nextapex);
05840 setorg(nextedge, NULL);
05841 setdest(nextedge, upperleft);
05842 setapex(nextedge, nextapex);
05843
05844 upperleft = nextapex;
05845
05846 otricopy(sidecasing, nextedge);
05847 apex(nextedge, nextapex);
05848 if (nextapex != (vertex) NULL) {
05849
05850 badedge = incircle(m, b, lowerleft, lowerright, upperleft,
05851 nextapex) > 0.0;
05852 } else {
05853
05854 badedge = 0;
05855 }
05856 }
05857 }
05858 }
05859
05860 if (!rightfinished) {
05861
05862 lnext(rightcand, nextedge);
05863 symself(nextedge);
05864 apex(nextedge, nextapex);
05865
05866
05867 if (nextapex != (vertex) NULL) {
05868
05869 badedge = incircle(m, b, lowerleft, lowerright, upperright, nextapex) >
05870 0.0;
05871 while (badedge) {
05872
05873
05874 lprevself(nextedge);
05875 sym(nextedge, topcasing);
05876 lprevself(nextedge);
05877 sym(nextedge, sidecasing);
05878 bond(nextedge, topcasing);
05879 bond(rightcand, sidecasing);
05880 lprevself(rightcand);
05881 sym(rightcand, outercasing);
05882 lnextself(nextedge);
05883 bond(nextedge, outercasing);
05884
05885 setorg(rightcand, NULL);
05886 setdest(rightcand, lowerright);
05887 setapex(rightcand, nextapex);
05888 setorg(nextedge, upperright);
05889 setdest(nextedge, NULL);
05890 setapex(nextedge, nextapex);
05891
05892 upperright = nextapex;
05893
05894 otricopy(sidecasing, nextedge);
05895 apex(nextedge, nextapex);
05896 if (nextapex != (vertex) NULL) {
05897
05898 badedge = incircle(m, b, lowerleft, lowerright, upperright,
05899 nextapex) > 0.0;
05900 } else {
05901
05902 badedge = 0;
05903 }
05904 }
05905 }
05906 }
05907 if (leftfinished || (!rightfinished &&
05908 (incircle(m, b, upperleft, lowerleft, lowerright, upperright) >
05909 0.0))) {
05910
05911
05912 bond(baseedge, rightcand);
05913 lprev(rightcand, baseedge);
05914 setdest(baseedge, lowerleft);
05915 lowerright = upperright;
05916 sym(baseedge, rightcand);
05917 apex(rightcand, upperright);
05918 } else {
05919
05920
05921 bond(baseedge, leftcand);
05922 lnext(leftcand, baseedge);
05923 setorg(baseedge, lowerright);
05924 lowerleft = upperleft;
05925 sym(baseedge, leftcand);
05926 apex(leftcand, upperleft);
05927 }
05928 if (b->verbose > 2) {
05929 printf(" Connecting ");
05930 printtriangle(m, b, &baseedge);
05931 }
05932 }
05933 }
05934
05935
05936
05937
05938
05939
05940
05941
05942
05943
05944
05945
05946
05947
05948
05949
05950
05951
05952 void divconqrecurse(struct mesh *m, struct behavior *b, vertex *sortarray,
05953 int vertices, int axis,
05954 struct otri *farleft, struct otri *farright)
05955 {
05956 struct otri midtri, tri1, tri2, tri3;
05957 struct otri innerleft, innerright;
05958 float area;
05959 int divider;
05960
05961 if (b->verbose > 2) {
05962 printf(" Triangulating %d vertices.\n", vertices);
05963 }
05964 if (vertices == 2) {
05965
05966
05967 maketriangle(m, b, farleft);
05968 setorg(*farleft, sortarray[0]);
05969 setdest(*farleft, sortarray[1]);
05970
05971 maketriangle(m, b, farright);
05972 setorg(*farright, sortarray[1]);
05973 setdest(*farright, sortarray[0]);
05974
05975 bond(*farleft, *farright);
05976 lprevself(*farleft);
05977 lnextself(*farright);
05978 bond(*farleft, *farright);
05979 lprevself(*farleft);
05980 lnextself(*farright);
05981 bond(*farleft, *farright);
05982 if (b->verbose > 2) {
05983 printf(" Creating ");
05984 printtriangle(m, b, farleft);
05985 printf(" Creating ");
05986 printtriangle(m, b, farright);
05987 }
05988
05989 lprev(*farright, *farleft);
05990 return;
05991 } else if (vertices == 3) {
05992
05993
05994
05995 maketriangle(m, b, &midtri);
05996 maketriangle(m, b, &tri1);
05997 maketriangle(m, b, &tri2);
05998 maketriangle(m, b, &tri3);
05999 area = counterclockwise(m, b, sortarray[0], sortarray[1], sortarray[2]);
06000 if (area == 0.0) {
06001
06002 setorg(midtri, sortarray[0]);
06003 setdest(midtri, sortarray[1]);
06004 setorg(tri1, sortarray[1]);
06005 setdest(tri1, sortarray[0]);
06006 setorg(tri2, sortarray[2]);
06007 setdest(tri2, sortarray[1]);
06008 setorg(tri3, sortarray[1]);
06009 setdest(tri3, sortarray[2]);
06010
06011 bond(midtri, tri1);
06012 bond(tri2, tri3);
06013 lnextself(midtri);
06014 lprevself(tri1);
06015 lnextself(tri2);
06016 lprevself(tri3);
06017 bond(midtri, tri3);
06018 bond(tri1, tri2);
06019 lnextself(midtri);
06020 lprevself(tri1);
06021 lnextself(tri2);
06022 lprevself(tri3);
06023 bond(midtri, tri1);
06024 bond(tri2, tri3);
06025
06026 otricopy(tri1, *farleft);
06027
06028 otricopy(tri2, *farright);
06029 } else {
06030
06031
06032 setorg(midtri, sortarray[0]);
06033 setdest(tri1, sortarray[0]);
06034 setorg(tri3, sortarray[0]);
06035
06036 if (area > 0.0) {
06037
06038 setdest(midtri, sortarray[1]);
06039 setorg(tri1, sortarray[1]);
06040 setdest(tri2, sortarray[1]);
06041 setapex(midtri, sortarray[2]);
06042 setorg(tri2, sortarray[2]);
06043 setdest(tri3, sortarray[2]);
06044 } else {
06045
06046 setdest(midtri, sortarray[2]);
06047 setorg(tri1, sortarray[2]);
06048 setdest(tri2, sortarray[2]);
06049 setapex(midtri, sortarray[1]);
06050 setorg(tri2, sortarray[1]);
06051 setdest(tri3, sortarray[1]);
06052 }
06053
06054 bond(midtri, tri1);
06055 lnextself(midtri);
06056 bond(midtri, tri2);
06057 lnextself(midtri);
06058 bond(midtri, tri3);
06059 lprevself(tri1);
06060 lnextself(tri2);
06061 bond(tri1, tri2);
06062 lprevself(tri1);
06063 lprevself(tri3);
06064 bond(tri1, tri3);
06065 lnextself(tri2);
06066 lprevself(tri3);
06067 bond(tri2, tri3);
06068
06069 otricopy(tri1, *farleft);
06070
06071 if (area > 0.0) {
06072 otricopy(tri2, *farright);
06073 } else {
06074 lnext(*farleft, *farright);
06075 }
06076 }
06077 if (b->verbose > 2) {
06078 printf(" Creating ");
06079 printtriangle(m, b, &midtri);
06080 printf(" Creating ");
06081 printtriangle(m, b, &tri1);
06082 printf(" Creating ");
06083 printtriangle(m, b, &tri2);
06084 printf(" Creating ");
06085 printtriangle(m, b, &tri3);
06086 }
06087 return;
06088 } else {
06089
06090 divider = vertices >> 1;
06091
06092 divconqrecurse(m, b, sortarray, divider, 1 - axis, farleft, &innerleft);
06093 divconqrecurse(m, b, &sortarray[divider], vertices - divider, 1 - axis,
06094 &innerright, farright);
06095 if (b->verbose > 1) {
06096 printf(" Joining triangulations with %d and %d vertices.\n", divider,
06097 vertices - divider);
06098 }
06099
06100 mergehulls(m, b, farleft, &innerleft, &innerright, farright, axis);
06101 }
06102 }
06103
06104 long removeghosts(struct mesh *m, struct behavior *b, struct otri *startghost)
06105 {
06106 struct otri searchedge;
06107 struct otri dissolveedge;
06108 struct otri deadtriangle;
06109 vertex markorg;
06110 long hullsize;
06111 triangle ptr;
06112
06113 if (b->verbose) {
06114 printf(" Removing ghost triangles.\n");
06115 }
06116
06117 lprev(*startghost, searchedge);
06118 symself(searchedge);
06119 m->dummytri[0] = encode(searchedge);
06120
06121 otricopy(*startghost, dissolveedge);
06122 hullsize = 0;
06123 do {
06124 hullsize++;
06125 lnext(dissolveedge, deadtriangle);
06126 lprevself(dissolveedge);
06127 symself(dissolveedge);
06128
06129
06130 if (!b->poly) {
06131
06132 if (dissolveedge.tri != m->dummytri) {
06133 org(dissolveedge, markorg);
06134 if (vertexmark(markorg) == 0) {
06135 setvertexmark(markorg, 1);
06136 }
06137 }
06138 }
06139
06140 dissolve(dissolveedge);
06141
06142 sym(deadtriangle, dissolveedge);
06143
06144 triangledealloc(m, deadtriangle.tri);
06145 } while (!otriequal(dissolveedge, *startghost));
06146 return hullsize;
06147 }
06148
06149
06150
06151
06152
06153
06154
06155
06156
06157
06158
06159 long divconqdelaunay(struct mesh *m, struct behavior *b)
06160 {
06161 vertex *sortarray;
06162 struct otri hullleft, hullright;
06163 int divider;
06164 int i, j;
06165
06166 if (b->verbose) {
06167 printf(" Sorting vertices.\n");
06168 }
06169
06170
06171 sortarray = (vertex *) trimalloc(m->invertices * (int) sizeof(vertex));
06172 traversalinit(&m->vertices);
06173 for (i = 0; i < m->invertices; i++) {
06174 sortarray[i] = vertextraverse(m);
06175 }
06176
06177 vertexsort(sortarray, m->invertices);
06178
06179 i = 0;
06180 for (j = 1; j < m->invertices; j++) {
06181 if ((sortarray[i][0] == sortarray[j][0])
06182 && (sortarray[i][1] == sortarray[j][1])) {
06183 if (!b->quiet) {
06184 printf(
06185 "Warning: A duplicate vertex at (%.12g, %.12g) appeared and was ignored.\n",
06186 sortarray[j][0], sortarray[j][1]);
06187 }
06188 setvertextype(sortarray[j], UNDEADVERTEX);
06189 m->undeads++;
06190 } else {
06191 i++;
06192 sortarray[i] = sortarray[j];
06193 }
06194 }
06195 i++;
06196 if (b->dwyer) {
06197
06198 divider = i >> 1;
06199 if (i - divider >= 2) {
06200 if (divider >= 2) {
06201 alternateaxes(sortarray, divider, 1);
06202 }
06203 alternateaxes(&sortarray[divider], i - divider, 1);
06204 }
06205 }
06206
06207 if (b->verbose) {
06208 printf(" Forming triangulation.\n");
06209 }
06210
06211
06212 divconqrecurse(m, b, sortarray, i, 0, &hullleft, &hullright);
06213 trifree((int *) sortarray);
06214
06215 return removeghosts(m, b, &hullleft);
06216 }
06217
06220
06221
06222
06226
06227
06228
06229
06230
06231
06232 long delaunay(struct mesh *m, struct behavior *b)
06233 {
06234 long hulledges;
06235
06236 m->eextras = 0;
06237 initializetrisubpools(m, b);
06238
06239 if (!b->quiet) {
06240 printf(
06241 "Constructing Delaunay triangulation by divide-and-conquer method.\n");
06242 }
06243 hulledges = divconqdelaunay(m, b);
06244
06245 if (m->triangles.items == 0) {
06246
06247 return 0l;
06248 } else {
06249 return hulledges;
06250 }
06251 }
06252
06255
06256
06257
06261
06262
06263
06264
06265
06266
06267
06268
06269
06270
06271
06272
06273
06274
06275
06276
06277
06278 enum finddirectionresult finddirection(struct mesh *m, struct behavior *b,
06279 struct otri *searchtri,
06280 vertex searchpoint)
06281 {
06282 struct otri checktri;
06283 vertex startvertex;
06284 vertex leftvertex, rightvertex;
06285 float leftccw, rightccw;
06286 int leftflag, rightflag;
06287 triangle ptr;
06288
06289 org(*searchtri, startvertex);
06290 dest(*searchtri, rightvertex);
06291 apex(*searchtri, leftvertex);
06292
06293 leftccw = counterclockwise(m, b, searchpoint, startvertex, leftvertex);
06294 leftflag = leftccw > 0.0;
06295
06296 rightccw = counterclockwise(m, b, startvertex, searchpoint, rightvertex);
06297 rightflag = rightccw > 0.0;
06298 if (leftflag && rightflag) {
06299
06300
06301 onext(*searchtri, checktri);
06302 if (checktri.tri == m->dummytri) {
06303 leftflag = 0;
06304 } else {
06305 rightflag = 0;
06306 }
06307 }
06308 while (leftflag) {
06309
06310 onextself(*searchtri);
06311 if (searchtri->tri == m->dummytri) {
06312 printf("Internal error in finddirection(): Unable to find a\n");
06313 printf(" triangle leading from (%.12g, %.12g) to", startvertex[0],
06314 startvertex[1]);
06315 printf(" (%.12g, %.12g).\n", searchpoint[0], searchpoint[1]);
06316 internalerror();
06317 }
06318 apex(*searchtri, leftvertex);
06319 rightccw = leftccw;
06320 leftccw = counterclockwise(m, b, searchpoint, startvertex, leftvertex);
06321 leftflag = leftccw > 0.0;
06322 }
06323 while (rightflag) {
06324
06325 oprevself(*searchtri);
06326 if (searchtri->tri == m->dummytri) {
06327 printf("Internal error in finddirection(): Unable to find a\n");
06328 printf(" triangle leading from (%.12g, %.12g) to", startvertex[0],
06329 startvertex[1]);
06330 printf(" (%.12g, %.12g).\n", searchpoint[0], searchpoint[1]);
06331 internalerror();
06332 }
06333 dest(*searchtri, rightvertex);
06334 leftccw = rightccw;
06335 rightccw = counterclockwise(m, b, startvertex, searchpoint, rightvertex);
06336 rightflag = rightccw > 0.0;
06337 }
06338 if (leftccw == 0.0) {
06339 return LEFTCOLLINEAR;
06340 } else if (rightccw == 0.0) {
06341 return RIGHTCOLLINEAR;
06342 } else {
06343 return WITHIN;
06344 }
06345 }
06346
06347
06348
06349
06350
06351
06352
06353
06354
06355
06356
06357
06358
06359
06360
06361
06362
06363
06364 void segmentintersection(struct mesh *m, struct behavior *b,
06365 struct otri *splittri, struct osub *splitsubseg,
06366 vertex endpoint2)
06367 {
06368 struct osub opposubseg;
06369 vertex endpoint1;
06370 vertex torg, tdest;
06371 vertex leftvertex, rightvertex;
06372 vertex newvertex;
06373 enum insertvertexresult success;
06374 enum finddirectionresult collinear;
06375 float ex, ey;
06376 float tx, ty;
06377 float etx, ety;
06378 float split, denom;
06379 int i;
06380 triangle ptr;
06381 subseg sptr;
06382
06383
06384 apex(*splittri, endpoint1);
06385 org(*splittri, torg);
06386 dest(*splittri, tdest);
06387
06388 tx = tdest[0] - torg[0];
06389 ty = tdest[1] - torg[1];
06390 ex = endpoint2[0] - endpoint1[0];
06391 ey = endpoint2[1] - endpoint1[1];
06392 etx = torg[0] - endpoint2[0];
06393 ety = torg[1] - endpoint2[1];
06394 denom = ty * ex - tx * ey;
06395 if (denom == 0.0) {
06396 printf("Internal error in segmentintersection():");
06397 printf(" Attempt to find intersection of parallel segments.\n");
06398 internalerror();
06399 }
06400 split = (ey * etx - ex * ety) / denom;
06401
06402 newvertex = (vertex) poolalloc(&m->vertices);
06403
06404 for (i = 0; i < 2 + m->nextras; i++) {
06405 newvertex[i] = torg[i] + split * (tdest[i] - torg[i]);
06406 }
06407 setvertexmark(newvertex, mark(*splitsubseg));
06408 setvertextype(newvertex, INPUTVERTEX);
06409 if (b->verbose > 1) {
06410 printf(
06411 " Splitting subsegment (%.12g, %.12g) (%.12g, %.12g) at (%.12g, %.12g).\n",
06412 torg[0], torg[1], tdest[0], tdest[1], newvertex[0], newvertex[1]);
06413 }
06414
06415 success = insertvertex(m, b, newvertex, splittri, splitsubseg, 0, 0);
06416 if (success != SUCCESSFULVERTEX) {
06417 printf("Internal error in segmentintersection():\n");
06418 printf(" Failure to split a segment.\n");
06419 internalerror();
06420 }
06421
06422 setvertex2tri(newvertex, encode(*splittri));
06423 if (m->steinerleft > 0) {
06424 m->steinerleft--;
06425 }
06426
06427
06428 ssymself(*splitsubseg);
06429 spivot(*splitsubseg, opposubseg);
06430 sdissolve(*splitsubseg);
06431 sdissolve(opposubseg);
06432 do {
06433 setsegorg(*splitsubseg, newvertex);
06434 snextself(*splitsubseg);
06435 } while (splitsubseg->ss != m->dummysub);
06436 do {
06437 setsegorg(opposubseg, newvertex);
06438 snextself(opposubseg);
06439 } while (opposubseg.ss != m->dummysub);
06440
06441
06442
06443 collinear = finddirection(m, b, splittri, endpoint1);
06444 dest(*splittri, rightvertex);
06445 apex(*splittri, leftvertex);
06446 if ((leftvertex[0] == endpoint1[0]) && (leftvertex[1] == endpoint1[1])) {
06447 onextself(*splittri);
06448 } else if ((rightvertex[0] != endpoint1[0]) ||
06449 (rightvertex[1] != endpoint1[1])) {
06450 printf("Internal error in segmentintersection():\n");
06451 printf(" Topological inconsistency after splitting a segment.\n");
06452 internalerror();
06453 }
06454
06455 }
06456
06457
06458
06459
06460
06461
06462
06463
06464
06465
06466
06467
06468
06469
06470
06471
06472
06473
06474
06475
06476
06477
06478
06479
06480
06481 int scoutsegment(struct mesh *m, struct behavior *b, struct otri *searchtri,
06482 vertex endpoint2, int newmark)
06483 {
06484 struct otri crosstri;
06485 struct osub crosssubseg;
06486 vertex leftvertex, rightvertex;
06487 enum finddirectionresult collinear;
06488 subseg sptr;
06489
06490 collinear = finddirection(m, b, searchtri, endpoint2);
06491 dest(*searchtri, rightvertex);
06492 apex(*searchtri, leftvertex);
06493 if (((leftvertex[0] == endpoint2[0]) && (leftvertex[1] == endpoint2[1])) ||
06494 ((rightvertex[0] == endpoint2[0]) && (rightvertex[1] == endpoint2[1]))) {
06495
06496 if ((leftvertex[0] == endpoint2[0]) && (leftvertex[1] == endpoint2[1])) {
06497 lprevself(*searchtri);
06498 }
06499
06500 insertsubseg(m, b, searchtri, newmark);
06501 return 1;
06502 } else if (collinear == LEFTCOLLINEAR) {
06503
06504
06505 lprevself(*searchtri);
06506 insertsubseg(m, b, searchtri, newmark);
06507
06508 return scoutsegment(m, b, searchtri, endpoint2, newmark);
06509 } else if (collinear == RIGHTCOLLINEAR) {
06510
06511 insertsubseg(m, b, searchtri, newmark);
06512
06513 lnextself(*searchtri);
06514
06515 return scoutsegment(m, b, searchtri, endpoint2, newmark);
06516 } else {
06517 lnext(*searchtri, crosstri);
06518 tspivot(crosstri, crosssubseg);
06519
06520 if (crosssubseg.ss == m->dummysub) {
06521 return 0;
06522 } else {
06523
06524 segmentintersection(m, b, &crosstri, &crosssubseg, endpoint2);
06525 otricopy(crosstri, *searchtri);
06526 insertsubseg(m, b, searchtri, newmark);
06527
06528 return scoutsegment(m, b, searchtri, endpoint2, newmark);
06529 }
06530 }
06531 }
06532
06533
06534
06535
06536
06537
06538
06539
06540
06541
06542
06543
06544
06545
06546
06547
06548
06549
06550
06551
06552
06553
06554
06555
06556
06557
06558
06559
06560
06561
06562
06563
06564
06565
06566
06567
06568
06569
06570
06571 void delaunayfixup(struct mesh *m, struct behavior *b,
06572 struct otri *fixuptri, int leftside)
06573 {
06574 struct otri neartri;
06575 struct otri fartri;
06576 struct osub faredge;
06577 vertex nearvertex, leftvertex, rightvertex, farvertex;
06578 triangle ptr;
06579 subseg sptr;
06580
06581 lnext(*fixuptri, neartri);
06582 sym(neartri, fartri);
06583
06584 if (fartri.tri == m->dummytri) {
06585 return;
06586 }
06587 tspivot(neartri, faredge);
06588 if (faredge.ss != m->dummysub) {
06589 return;
06590 }
06591
06592 apex(neartri, nearvertex);
06593 org(neartri, leftvertex);
06594 dest(neartri, rightvertex);
06595 apex(fartri, farvertex);
06596
06597 if (leftside) {
06598 if (counterclockwise(m, b, nearvertex, leftvertex, farvertex) <= 0.0) {
06599
06600
06601 return;
06602 }
06603 } else {
06604 if (counterclockwise(m, b, farvertex, rightvertex, nearvertex) <= 0.0) {
06605
06606
06607 return;
06608 }
06609 }
06610 if (counterclockwise(m, b, rightvertex, leftvertex, farvertex) > 0.0) {
06611
06612
06613
06614
06615 if (incircle(m, b, leftvertex, farvertex, rightvertex, nearvertex) <=
06616 0.0) {
06617 return;
06618 }
06619
06620 }
06621 flip(m, b, &neartri);
06622 lprevself(*fixuptri);
06623
06624 delaunayfixup(m, b, fixuptri, leftside);
06625 delaunayfixup(m, b, &fartri, leftside);
06626 }
06627
06628
06629
06630
06631
06632
06633
06634
06635
06636
06637
06638
06639
06640
06641
06642
06643
06644
06645
06646
06647
06648
06649
06650
06651
06652
06653
06654
06655
06656
06657
06658
06659
06660
06661
06662
06663
06664
06665
06666
06667
06668
06669
06670
06671
06672
06673
06674
06675
06676
06677
06678
06679
06680
06681
06682 void constrainededge(struct mesh *m, struct behavior *b,
06683 struct otri *starttri, vertex endpoint2, int newmark)
06684 {
06685 struct otri fixuptri, fixuptri2;
06686 struct osub crosssubseg;
06687 vertex endpoint1;
06688 vertex farvertex;
06689 float area;
06690 int collision;
06691 int done;
06692 triangle ptr;
06693 subseg sptr;
06694
06695 org(*starttri, endpoint1);
06696 lnext(*starttri, fixuptri);
06697 flip(m, b, &fixuptri);
06698
06699
06700 collision = 0;
06701 done = 0;
06702 do {
06703 org(fixuptri, farvertex);
06704
06705
06706 if ((farvertex[0] == endpoint2[0]) && (farvertex[1] == endpoint2[1])) {
06707 oprev(fixuptri, fixuptri2);
06708
06709 delaunayfixup(m, b, &fixuptri, 0);
06710 delaunayfixup(m, b, &fixuptri2, 1);
06711 done = 1;
06712 } else {
06713
06714
06715
06716 area = counterclockwise(m, b, endpoint1, endpoint2, farvertex);
06717 if (area == 0.0) {
06718
06719 collision = 1;
06720 oprev(fixuptri, fixuptri2);
06721
06722 delaunayfixup(m, b, &fixuptri, 0);
06723 delaunayfixup(m, b, &fixuptri2, 1);
06724 done = 1;
06725 } else {
06726 if (area > 0.0) {
06727 oprev(fixuptri, fixuptri2);
06728
06729
06730 delaunayfixup(m, b, &fixuptri2, 1);
06731
06732
06733
06734 lprevself(fixuptri);
06735 } else {
06736 delaunayfixup(m, b, &fixuptri, 0);
06737
06738
06739
06740 oprevself(fixuptri);
06741 }
06742
06743 tspivot(fixuptri, crosssubseg);
06744 if (crosssubseg.ss == m->dummysub) {
06745 flip(m, b, &fixuptri);
06746 } else {
06747
06748 collision = 1;
06749
06750 segmentintersection(m, b, &fixuptri, &crosssubseg, endpoint2);
06751 done = 1;
06752 }
06753 }
06754 }
06755 } while (!done);
06756
06757 insertsubseg(m, b, &fixuptri, newmark);
06758
06759
06760 if (collision) {
06761
06762 if (!scoutsegment(m, b, &fixuptri, endpoint2, newmark)) {
06763 constrainededge(m, b, &fixuptri, endpoint2, newmark);
06764 }
06765 }
06766 }
06767
06768
06769
06770
06771
06772
06773
06774 void insertsegment(struct mesh *m, struct behavior *b,
06775 vertex endpoint1, vertex endpoint2, int newmark)
06776 {
06777 struct otri searchtri1, searchtri2;
06778 triangle encodedtri;
06779 vertex checkvertex;
06780 triangle ptr;
06781
06782 if (b->verbose > 1) {
06783 printf(" Connecting (%.12g, %.12g) to (%.12g, %.12g).\n",
06784 endpoint1[0], endpoint1[1], endpoint2[0], endpoint2[1]);
06785 }
06786
06787
06788 checkvertex = (vertex) NULL;
06789 encodedtri = vertex2tri(endpoint1);
06790 if (encodedtri != (triangle) NULL) {
06791 decode(encodedtri, searchtri1);
06792 org(searchtri1, checkvertex);
06793 }
06794 if (checkvertex != endpoint1) {
06795
06796 searchtri1.tri = m->dummytri;
06797 searchtri1.orient = 0;
06798 symself(searchtri1);
06799
06800 if (locate(m, b, endpoint1, &searchtri1) != ONVERTEX) {
06801 printf(
06802 "Internal error in insertsegment(): Unable to locate PSLG vertex\n");
06803 printf(" (%.12g, %.12g) in triangulation.\n",
06804 endpoint1[0], endpoint1[1]);
06805 internalerror();
06806 }
06807 }
06808
06809 otricopy(searchtri1, m->recenttri);
06810
06811
06812 if (scoutsegment(m, b, &searchtri1, endpoint2, newmark)) {
06813
06814 return;
06815 }
06816
06817
06818 org(searchtri1, endpoint1);
06819
06820
06821 checkvertex = (vertex) NULL;
06822 encodedtri = vertex2tri(endpoint2);
06823 if (encodedtri != (triangle) NULL) {
06824 decode(encodedtri, searchtri2);
06825 org(searchtri2, checkvertex);
06826 }
06827 if (checkvertex != endpoint2) {
06828
06829 searchtri2.tri = m->dummytri;
06830 searchtri2.orient = 0;
06831 symself(searchtri2);
06832
06833 if (locate(m, b, endpoint2, &searchtri2) != ONVERTEX) {
06834 printf(
06835 "Internal error in insertsegment(): Unable to locate PSLG vertex\n");
06836 printf(" (%.12g, %.12g) in triangulation.\n",
06837 endpoint2[0], endpoint2[1]);
06838 internalerror();
06839 }
06840 }
06841
06842 otricopy(searchtri2, m->recenttri);
06843
06844
06845 if (scoutsegment(m, b, &searchtri2, endpoint1, newmark)) {
06846
06847 return;
06848 }
06849
06850
06851 org(searchtri2, endpoint2);
06852
06853
06854 constrainededge(m, b, &searchtri1, endpoint2, newmark);
06855 }
06856
06857
06858
06859
06860
06861
06862
06863 void markhull(struct mesh *m, struct behavior *b)
06864 {
06865 struct otri hulltri;
06866 struct otri nexttri;
06867 struct otri starttri;
06868 triangle ptr;
06869
06870
06871 hulltri.tri = m->dummytri;
06872 hulltri.orient = 0;
06873 symself(hulltri);
06874
06875 otricopy(hulltri, starttri);
06876
06877 do {
06878
06879 insertsubseg(m, b, &hulltri, 1);
06880
06881 lnextself(hulltri);
06882 oprev(hulltri, nexttri);
06883 while (nexttri.tri != m->dummytri) {
06884 otricopy(nexttri, hulltri);
06885 oprev(hulltri, nexttri);
06886 }
06887 } while (!otriequal(hulltri, starttri));
06888 }
06889
06890
06891
06892
06893
06894
06895
06896
06897
06898
06899
06900 void formskeleton(struct mesh *m, struct behavior *b, int *segmentlist,
06901 int *segmentmarkerlist, int numberofsegments)
06902 {
06903 char polyfilename[6];
06904 int index;
06905 vertex endpoint1, endpoint2;
06906 int segmentmarkers;
06907 int end1, end2;
06908 int boundmarker;
06909 int i;
06910
06911 if (b->poly) {
06912 if (!b->quiet) {
06913 printf("Recovering segments in Delaunay triangulation.\n");
06914 }
06915 strcpy(polyfilename, "input");
06916 m->insegments = numberofsegments;
06917 segmentmarkers = segmentmarkerlist != (int *) NULL;
06918 index = 0;
06919
06920
06921 if (m->triangles.items == 0) {
06922 return;
06923 }
06924
06925
06926
06927 if (m->insegments > 0) {
06928 makevertexmap(m, b);
06929 if (b->verbose) {
06930 printf(" Recovering PSLG segments.\n");
06931 }
06932 }
06933
06934 boundmarker = 0;
06935
06936 for (i = 0; i < m->insegments; i++) {
06937 end1 = segmentlist[index++];
06938 end2 = segmentlist[index++];
06939 if (segmentmarkers) {
06940 boundmarker = segmentmarkerlist[i];
06941 }
06942 if ((end1 < b->firstnumber) ||
06943 (end1 >= b->firstnumber + m->invertices)) {
06944 if (!b->quiet) {
06945 printf("Warning: Invalid first endpoint of segment %d in %s.\n",
06946 b->firstnumber + i, polyfilename);
06947 }
06948 } else if ((end2 < b->firstnumber) ||
06949 (end2 >= b->firstnumber + m->invertices)) {
06950 if (!b->quiet) {
06951 printf("Warning: Invalid second endpoint of segment %d in %s.\n",
06952 b->firstnumber + i, polyfilename);
06953 }
06954 } else {
06955
06956 endpoint1 = getvertex(m, b, end1);
06957 endpoint2 = getvertex(m, b, end2);
06958 if ((endpoint1[0] == endpoint2[0]) && (endpoint1[1] == endpoint2[1])) {
06959 if (!b->quiet) {
06960 printf("Warning: Endpoints of segment %d are coincident in %s.\n",
06961 b->firstnumber + i, polyfilename);
06962 }
06963 } else {
06964 insertsegment(m, b, endpoint1, endpoint2, boundmarker);
06965 }
06966 }
06967 }
06968 } else {
06969 m->insegments = 0;
06970 }
06971 if (b->convex || !b->poly) {
06972
06973 if (b->verbose) {
06974 printf(" Enclosing convex hull with segments.\n");
06975 }
06976 markhull(m, b);
06977 }
06978 }
06979
06982
06983
06984
06988
06989
06990
06991
06992
06993
06994
06995
06996 void infecthull(struct mesh *m, struct behavior *b)
06997 {
06998 struct otri hulltri;
06999 struct otri nexttri;
07000 struct otri starttri;
07001 struct osub hullsubseg;
07002 triangle **deadtriangle;
07003 vertex horg, hdest;
07004 triangle ptr;
07005 subseg sptr;
07006
07007 if (b->verbose) {
07008 printf(" Marking concavities (external triangles) for elimination.\n");
07009 }
07010
07011 hulltri.tri = m->dummytri;
07012 hulltri.orient = 0;
07013 symself(hulltri);
07014
07015 otricopy(hulltri, starttri);
07016
07017 do {
07018
07019 if (!infected(hulltri)) {
07020
07021 tspivot(hulltri, hullsubseg);
07022 if (hullsubseg.ss == m->dummysub) {
07023
07024 if (!infected(hulltri)) {
07025 infect(hulltri);
07026 deadtriangle = (triangle **) poolalloc(&m->viri);
07027 *deadtriangle = hulltri.tri;
07028 }
07029 } else {
07030
07031 if (mark(hullsubseg) == 0) {
07032 setmark(hullsubseg, 1);
07033 org(hulltri, horg);
07034 dest(hulltri, hdest);
07035 if (vertexmark(horg) == 0) {
07036 setvertexmark(horg, 1);
07037 }
07038 if (vertexmark(hdest) == 0) {
07039 setvertexmark(hdest, 1);
07040 }
07041 }
07042 }
07043 }
07044
07045 lnextself(hulltri);
07046 oprev(hulltri, nexttri);
07047 while (nexttri.tri != m->dummytri) {
07048 otricopy(nexttri, hulltri);
07049 oprev(hulltri, nexttri);
07050 }
07051 } while (!otriequal(hulltri, starttri));
07052 }
07053
07054
07055
07056
07057
07058
07059
07060
07061
07062
07063
07064
07065
07066
07067
07068
07069
07070
07071 void plague(struct mesh *m, struct behavior *b)
07072 {
07073 struct otri testtri;
07074 struct otri neighbor;
07075 triangle **virusloop;
07076 triangle **deadtriangle;
07077 struct osub neighborsubseg;
07078 vertex testvertex;
07079 vertex norg, ndest;
07080 vertex deadorg, deaddest, deadapex;
07081 int killorg;
07082 triangle ptr;
07083 subseg sptr;
07084
07085 if (b->verbose) {
07086 printf(" Marking neighbors of marked triangles.\n");
07087 }
07088
07089
07090 traversalinit(&m->viri);
07091 virusloop = (triangle **) traverse(&m->viri);
07092 while (virusloop != (triangle **) NULL) {
07093 testtri.tri = *virusloop;
07094
07095
07096
07097
07098 uninfect(testtri);
07099 if (b->verbose > 2) {
07100
07101
07102 testtri.orient = 0;
07103 org(testtri, deadorg);
07104 dest(testtri, deaddest);
07105 apex(testtri, deadapex);
07106 printf(" Checking (%.12g, %.12g) (%.12g, %.12g) (%.12g, %.12g)\n",
07107 deadorg[0], deadorg[1], deaddest[0], deaddest[1],
07108 deadapex[0], deadapex[1]);
07109 }
07110
07111 for (testtri.orient = 0; testtri.orient < 3; testtri.orient++) {
07112
07113 sym(testtri, neighbor);
07114
07115 tspivot(testtri, neighborsubseg);
07116
07117 if ((neighbor.tri == m->dummytri) || infected(neighbor)) {
07118 if (neighborsubseg.ss != m->dummysub) {
07119
07120
07121
07122 subsegdealloc(m, neighborsubseg.ss);
07123 if (neighbor.tri != m->dummytri) {
07124
07125
07126 uninfect(neighbor);
07127 tsdissolve(neighbor);
07128 infect(neighbor);
07129 }
07130 }
07131 } else {
07132 if (neighborsubseg.ss == m->dummysub) {
07133
07134
07135 if (b->verbose > 2) {
07136 org(neighbor, deadorg);
07137 dest(neighbor, deaddest);
07138 apex(neighbor, deadapex);
07139 printf(
07140 " Marking (%.12g, %.12g) (%.12g, %.12g) (%.12g, %.12g)\n",
07141 deadorg[0], deadorg[1], deaddest[0], deaddest[1],
07142 deadapex[0], deadapex[1]);
07143 }
07144 infect(neighbor);
07145
07146 deadtriangle = (triangle **) poolalloc(&m->viri);
07147 *deadtriangle = neighbor.tri;
07148 } else {
07149
07150 stdissolve(neighborsubseg);
07151
07152 if (mark(neighborsubseg) == 0) {
07153 setmark(neighborsubseg, 1);
07154 }
07155 org(neighbor, norg);
07156 dest(neighbor, ndest);
07157 if (vertexmark(norg) == 0) {
07158 setvertexmark(norg, 1);
07159 }
07160 if (vertexmark(ndest) == 0) {
07161 setvertexmark(ndest, 1);
07162 }
07163 }
07164 }
07165 }
07166
07167
07168 infect(testtri);
07169 virusloop = (triangle **) traverse(&m->viri);
07170 }
07171
07172 if (b->verbose) {
07173 printf(" Deleting marked triangles.\n");
07174 }
07175
07176 traversalinit(&m->viri);
07177 virusloop = (triangle **) traverse(&m->viri);
07178 while (virusloop != (triangle **) NULL) {
07179 testtri.tri = *virusloop;
07180
07181
07182
07183
07184 for (testtri.orient = 0; testtri.orient < 3; testtri.orient++) {
07185 org(testtri, testvertex);
07186
07187 if (testvertex != (vertex) NULL) {
07188 killorg = 1;
07189
07190 setorg(testtri, NULL);
07191
07192 onext(testtri, neighbor);
07193
07194 while ((neighbor.tri != m->dummytri) &&
07195 (!otriequal(neighbor, testtri))) {
07196 if (infected(neighbor)) {
07197
07198 setorg(neighbor, NULL);
07199 } else {
07200
07201 killorg = 0;
07202 }
07203
07204 onextself(neighbor);
07205 }
07206
07207 if (neighbor.tri == m->dummytri) {
07208
07209 oprev(testtri, neighbor);
07210
07211 while (neighbor.tri != m->dummytri) {
07212 if (infected(neighbor)) {
07213
07214 setorg(neighbor, NULL);
07215 } else {
07216
07217 killorg = 0;
07218 }
07219
07220 oprevself(neighbor);
07221 }
07222 }
07223 if (killorg) {
07224 if (b->verbose > 1) {
07225 printf(" Deleting vertex (%.12g, %.12g)\n",
07226 testvertex[0], testvertex[1]);
07227 }
07228 setvertextype(testvertex, UNDEADVERTEX);
07229 m->undeads++;
07230 }
07231 }
07232 }
07233
07234
07235
07236 for (testtri.orient = 0; testtri.orient < 3; testtri.orient++) {
07237 sym(testtri, neighbor);
07238 if (neighbor.tri == m->dummytri) {
07239
07240
07241
07242 m->hullsize--;
07243 } else {
07244
07245 dissolve(neighbor);
07246
07247
07248 m->hullsize++;
07249 }
07250 }
07251
07252 triangledealloc(m, testtri.tri);
07253 virusloop = (triangle **) traverse(&m->viri);
07254 }
07255
07256 poolrestart(&m->viri);
07257 }
07258
07259
07260
07261
07262
07263
07264
07265
07266
07267
07268
07269
07270
07271
07272
07273
07274 void regionplague(struct mesh *m, struct behavior *b,
07275 float attribute, float area)
07276 {
07277 struct otri testtri;
07278 struct otri neighbor;
07279 triangle **virusloop;
07280 triangle **regiontri;
07281 struct osub neighborsubseg;
07282 vertex regionorg, regiondest, regionapex;
07283 triangle ptr;
07284 subseg sptr;
07285
07286 if (b->verbose > 1) {
07287 printf(" Marking neighbors of marked triangles.\n");
07288 }
07289
07290
07291
07292 traversalinit(&m->viri);
07293 virusloop = (triangle **) traverse(&m->viri);
07294 while (virusloop != (triangle **) NULL) {
07295 testtri.tri = *virusloop;
07296
07297
07298
07299
07300 uninfect(testtri);
07301 if (b->regionattrib) {
07302
07303 setelemattribute(testtri, m->eextras, attribute);
07304 }
07305 if (b->vararea) {
07306
07307 setareabound(testtri, area);
07308 }
07309 if (b->verbose > 2) {
07310
07311
07312 testtri.orient = 0;
07313 org(testtri, regionorg);
07314 dest(testtri, regiondest);
07315 apex(testtri, regionapex);
07316 printf(" Checking (%.12g, %.12g) (%.12g, %.12g) (%.12g, %.12g)\n",
07317 regionorg[0], regionorg[1], regiondest[0], regiondest[1],
07318 regionapex[0], regionapex[1]);
07319 }
07320
07321 for (testtri.orient = 0; testtri.orient < 3; testtri.orient++) {
07322
07323 sym(testtri, neighbor);
07324
07325 tspivot(testtri, neighborsubseg);
07326
07327
07328 if ((neighbor.tri != m->dummytri) && !infected(neighbor)
07329 && (neighborsubseg.ss == m->dummysub)) {
07330 if (b->verbose > 2) {
07331 org(neighbor, regionorg);
07332 dest(neighbor, regiondest);
07333 apex(neighbor, regionapex);
07334 printf(" Marking (%.12g, %.12g) (%.12g, %.12g) (%.12g, %.12g)\n",
07335 regionorg[0], regionorg[1], regiondest[0], regiondest[1],
07336 regionapex[0], regionapex[1]);
07337 }
07338
07339 infect(neighbor);
07340
07341 regiontri = (triangle **) poolalloc(&m->viri);
07342 *regiontri = neighbor.tri;
07343 }
07344 }
07345
07346
07347 infect(testtri);
07348 virusloop = (triangle **) traverse(&m->viri);
07349 }
07350
07351
07352 if (b->verbose > 1) {
07353 printf(" Unmarking marked triangles.\n");
07354 }
07355 traversalinit(&m->viri);
07356 virusloop = (triangle **) traverse(&m->viri);
07357 while (virusloop != (triangle **) NULL) {
07358 testtri.tri = *virusloop;
07359 uninfect(testtri);
07360 virusloop = (triangle **) traverse(&m->viri);
07361 }
07362
07363 poolrestart(&m->viri);
07364 }
07365
07366
07367
07368
07369
07370
07371
07372
07373
07374
07375
07376
07377
07378 void carveholes(struct mesh *m, struct behavior *b, float *holelist, int holes,
07379 float *regionlist, int regions)
07380 {
07381 struct otri searchtri;
07382 struct otri triangleloop;
07383 struct otri *regiontris;
07384 triangle **holetri;
07385 triangle **regiontri;
07386 vertex searchorg, searchdest;
07387 enum locateresult intersect;
07388 int i;
07389 triangle ptr;
07390
07391 if (!(b->quiet || (b->noholes && b->convex))) {
07392 printf("Removing unwanted triangles.\n");
07393 if (b->verbose && (holes > 0)) {
07394 printf(" Marking holes for elimination.\n");
07395 }
07396 }
07397
07398 if (regions > 0) {
07399
07400 regiontris = (struct otri *) trimalloc(regions *
07401 (int) sizeof(struct otri));
07402 } else {
07403 regiontris = (struct otri *) NULL;
07404 }
07405
07406 if (((holes > 0) && !b->noholes) || !b->convex || (regions > 0)) {
07407
07408
07409 poolinit(&m->viri, sizeof(triangle *), VIRUSPERBLOCK, VIRUSPERBLOCK, 0);
07410 }
07411
07412 if (!b->convex) {
07413
07414
07415 infecthull(m, b);
07416 }
07417
07418 if ((holes > 0) && !b->noholes) {
07419
07420 for (i = 0; i < 2 * holes; i += 2) {
07421
07422 if ((holelist[i] >= m->xmin) && (holelist[i] <= m->xmax)
07423 && (holelist[i + 1] >= m->ymin) && (holelist[i + 1] <= m->ymax)) {
07424
07425 searchtri.tri = m->dummytri;
07426 searchtri.orient = 0;
07427 symself(searchtri);
07428
07429
07430
07431 org(searchtri, searchorg);
07432 dest(searchtri, searchdest);
07433 if (counterclockwise(m, b, searchorg, searchdest, &holelist[i]) >
07434 0.0) {
07435
07436 intersect = locate(m, b, &holelist[i], &searchtri);
07437 if ((intersect != OUTSIDE) && (!infected(searchtri))) {
07438
07439
07440 infect(searchtri);
07441 holetri = (triangle **) poolalloc(&m->viri);
07442 *holetri = searchtri.tri;
07443 }
07444 }
07445 }
07446 }
07447 }
07448
07449
07450
07451
07452
07453
07454
07455 if (regions > 0) {
07456
07457 for (i = 0; i < regions; i++) {
07458 regiontris[i].tri = m->dummytri;
07459
07460 if ((regionlist[4 * i] >= m->xmin) && (regionlist[4 * i] <= m->xmax) &&
07461 (regionlist[4 * i + 1] >= m->ymin) &&
07462 (regionlist[4 * i + 1] <= m->ymax)) {
07463
07464 searchtri.tri = m->dummytri;
07465 searchtri.orient = 0;
07466 symself(searchtri);
07467
07468
07469
07470 org(searchtri, searchorg);
07471 dest(searchtri, searchdest);
07472 if (counterclockwise(m, b, searchorg, searchdest, ®ionlist[4 * i]) >
07473 0.0) {
07474
07475 intersect = locate(m, b, ®ionlist[4 * i], &searchtri);
07476 if ((intersect != OUTSIDE) && (!infected(searchtri))) {
07477
07478
07479 otricopy(searchtri, regiontris[i]);
07480 }
07481 }
07482 }
07483 }
07484 }
07485
07486 if (m->viri.items > 0) {
07487
07488 plague(m, b);
07489 }
07490
07491
07492 if (regions > 0) {
07493 if (!b->quiet) {
07494 if (b->regionattrib) {
07495 if (b->vararea) {
07496 printf("Spreading regional attributes and area constraints.\n");
07497 } else {
07498 printf("Spreading regional attributes.\n");
07499 }
07500 } else {
07501 printf("Spreading regional area constraints.\n");
07502 }
07503 }
07504 if (b->regionattrib && !b->refine) {
07505
07506 traversalinit(&m->triangles);
07507 triangleloop.orient = 0;
07508 triangleloop.tri = triangletraverse(m);
07509 while (triangleloop.tri != (triangle *) NULL) {
07510 setelemattribute(triangleloop, m->eextras, 0.0);
07511 triangleloop.tri = triangletraverse(m);
07512 }
07513 }
07514 for (i = 0; i < regions; i++) {
07515 if (regiontris[i].tri != m->dummytri) {
07516
07517
07518 if (!deadtri(regiontris[i].tri)) {
07519
07520 infect(regiontris[i]);
07521 regiontri = (triangle **) poolalloc(&m->viri);
07522 *regiontri = regiontris[i].tri;
07523
07524 regionplague(m, b, regionlist[4 * i + 2], regionlist[4 * i + 3]);
07525
07526 }
07527 }
07528 }
07529 if (b->regionattrib && !b->refine) {
07530
07531 m->eextras++;
07532 }
07533 }
07534
07535
07536 if (((holes > 0) && !b->noholes) || !b->convex || (regions > 0)) {
07537 pooldeinit(&m->viri);
07538 }
07539 if (regions > 0) {
07540 trifree((int *) regiontris);
07541 }
07542 }
07543
07546
07547
07548
07549
07550
07551
07552
07553
07554 void highorder(struct mesh *m, struct behavior *b)
07555 {
07556 struct otri triangleloop, trisym;
07557 struct osub checkmark;
07558 vertex newvertex;
07559 vertex torg, tdest;
07560 int i;
07561 triangle ptr;
07562 subseg sptr;
07563
07564 if (!b->quiet) {
07565 printf("Adding vertices for second-order triangles.\n");
07566 }
07567
07568
07569
07570
07571
07572 m->vertices.deaditemstack = (int *) NULL;
07573
07574 traversalinit(&m->triangles);
07575 triangleloop.tri = triangletraverse(m);
07576
07577
07578
07579
07580
07581
07582 while (triangleloop.tri != (triangle *) NULL) {
07583 for (triangleloop.orient = 0; triangleloop.orient < 3;
07584 triangleloop.orient++) {
07585 sym(triangleloop, trisym);
07586 if ((triangleloop.tri < trisym.tri) || (trisym.tri == m->dummytri)) {
07587 org(triangleloop, torg);
07588 dest(triangleloop, tdest);
07589
07590
07591 newvertex = (vertex) poolalloc(&m->vertices);
07592 for (i = 0; i < 2 + m->nextras; i++) {
07593 newvertex[i] = 0.5 * (torg[i] + tdest[i]);
07594 }
07595
07596
07597 setvertexmark(newvertex, trisym.tri == m->dummytri);
07598 setvertextype(newvertex,
07599 trisym.tri == m->dummytri ? FREEVERTEX : SEGMENTVERTEX);
07600 if (b->usesegments) {
07601 tspivot(triangleloop, checkmark);
07602
07603 if (checkmark.ss != m->dummysub) {
07604 setvertexmark(newvertex, mark(checkmark));
07605 setvertextype(newvertex, SEGMENTVERTEX);
07606 }
07607 }
07608 if (b->verbose > 1) {
07609 printf(" Creating (%.12g, %.12g).\n", newvertex[0], newvertex[1]);
07610 }
07611
07612 triangleloop.tri[m->highorderindex + triangleloop.orient] =
07613 (triangle) newvertex;
07614 if (trisym.tri != m->dummytri) {
07615 trisym.tri[m->highorderindex + trisym.orient] = (triangle) newvertex;
07616 }
07617 }
07618 }
07619 triangleloop.tri = triangletraverse(m);
07620 }
07621 }
07622
07623
07627
07628
07629
07630
07631
07632
07633 void transfernodes(struct mesh *m, struct behavior *b, float *pointlist,
07634 float *pointattriblist, int *pointmarkerlist,
07635 int numberofpoints, int numberofpointattribs)
07636 {
07637 vertex vertexloop;
07638 float x, y;
07639 int i, j;
07640 int coordindex;
07641 int attribindex;
07642
07643 m->invertices = numberofpoints;
07644 m->mesh_dim = 2;
07645 m->nextras = numberofpointattribs;
07646 m->readnodefile = 0;
07647 if (m->invertices < 3) {
07648 printf("Error: Input must have at least three input vertices.\n");
07649 triexit(1);
07650 }
07651 if (m->nextras == 0) {
07652 b->weighted = 0;
07653 }
07654
07655 initializevertexpool(m, b);
07656
07657
07658 coordindex = 0;
07659 attribindex = 0;
07660 for (i = 0; i < m->invertices; i++) {
07661 vertexloop = (vertex) poolalloc(&m->vertices);
07662
07663 x = vertexloop[0] = pointlist[coordindex++];
07664 y = vertexloop[1] = pointlist[coordindex++];
07665
07666 for (j = 0; j < numberofpointattribs; j++) {
07667 vertexloop[2 + j] = pointattriblist[attribindex++];
07668 }
07669 if (pointmarkerlist != (int *) NULL) {
07670
07671 setvertexmark(vertexloop, pointmarkerlist[i]);
07672 } else {
07673
07674 setvertexmark(vertexloop, 0);
07675 }
07676 setvertextype(vertexloop, INPUTVERTEX);
07677
07678 if (i == 0) {
07679 m->xmin = m->xmax = x;
07680 m->ymin = m->ymax = y;
07681 } else {
07682 m->xmin = (x < m->xmin) ? x : m->xmin;
07683 m->xmax = (x > m->xmax) ? x : m->xmax;
07684 m->ymin = (y < m->ymin) ? y : m->ymin;
07685 m->ymax = (y > m->ymax) ? y : m->ymax;
07686 }
07687 }
07688
07689
07690
07691 m->xminextreme = 10 * m->xmin - 9 * m->xmax;
07692 }
07693
07694
07695
07696
07697
07698
07699
07700
07701
07702
07703 void writenodes(struct mesh *m, struct behavior *b, float **pointlist,
07704 float **pointattriblist, int **pointmarkerlist)
07705 {
07706 float *plist;
07707 float *palist;
07708 int *pmlist;
07709 int coordindex;
07710 int attribindex;
07711 vertex vertexloop;
07712 long outvertices;
07713 int vertexnumber;
07714 int i;
07715
07716 if (b->jettison) {
07717 outvertices = m->vertices.items - m->undeads;
07718 } else {
07719 outvertices = m->vertices.items;
07720 }
07721
07722 if (!b->quiet) {
07723 printf("Writing vertices.\n");
07724 }
07725
07726 if (*pointlist == (float *) NULL) {
07727 *pointlist = (float *) trimalloc((int) (outvertices * 2 * sizeof(float)));
07728 }
07729
07730 if ((m->nextras > 0) && (*pointattriblist == (float *) NULL)) {
07731 *pointattriblist = (float *) trimalloc((int) (outvertices * m->nextras *
07732 sizeof(float)));
07733 }
07734
07735 if (!b->nobound && (*pointmarkerlist == (int *) NULL)) {
07736 *pointmarkerlist = (int *) trimalloc((int) (outvertices * sizeof(int)));
07737 }
07738 plist = *pointlist;
07739 palist = *pointattriblist;
07740 pmlist = *pointmarkerlist;
07741 coordindex = 0;
07742 attribindex = 0;
07743 traversalinit(&m->vertices);
07744 vertexnumber = b->firstnumber;
07745 vertexloop = vertextraverse(m);
07746 while (vertexloop != (vertex) NULL) {
07747 if (!b->jettison || (vertextype(vertexloop) != UNDEADVERTEX)) {
07748
07749 plist[coordindex++] = vertexloop[0];
07750 plist[coordindex++] = vertexloop[1];
07751
07752 for (i = 0; i < m->nextras; i++) {
07753 palist[attribindex++] = vertexloop[2 + i];
07754 }
07755 if (!b->nobound) {
07756
07757 pmlist[vertexnumber - b->firstnumber] = vertexmark(vertexloop);
07758 }
07759 setvertexmark(vertexloop, vertexnumber);
07760 vertexnumber++;
07761 }
07762 vertexloop = vertextraverse(m);
07763 }
07764 }
07765
07766
07767
07768
07769
07770
07771
07772
07773
07774
07775
07776 void numbernodes(struct mesh *m, struct behavior *b)
07777 {
07778 vertex vertexloop;
07779 int vertexnumber;
07780
07781 traversalinit(&m->vertices);
07782 vertexnumber = b->firstnumber;
07783 vertexloop = vertextraverse(m);
07784 while (vertexloop != (vertex) NULL) {
07785 setvertexmark(vertexloop, vertexnumber);
07786 if (!b->jettison || (vertextype(vertexloop) != UNDEADVERTEX)) {
07787 vertexnumber++;
07788 }
07789 vertexloop = vertextraverse(m);
07790 }
07791 }
07792
07793
07794
07795
07796
07797
07798
07799 void writeelements(struct mesh *m, struct behavior *b,
07800 int **trianglelist, float **triangleattriblist)
07801 {
07802 int *tlist;
07803 float *talist;
07804 int vertexindex;
07805 int attribindex;
07806 struct otri triangleloop;
07807 vertex p1, p2, p3;
07808 vertex mid1, mid2, mid3;
07809 long elementnumber;
07810 int i;
07811
07812 if (!b->quiet) {
07813 printf("Writing triangles.\n");
07814 }
07815
07816 if (*trianglelist == (int *) NULL) {
07817 *trianglelist = (int *) trimalloc((int) (m->triangles.items *
07818 ((b->order + 1) * (b->order + 2) /
07819 2) * sizeof(int)));
07820 }
07821
07822 if ((m->eextras > 0) && (*triangleattriblist == (float *) NULL)) {
07823 *triangleattriblist = (float *) trimalloc((int) (m->triangles.items *
07824 m->eextras *
07825 sizeof(float)));
07826 }
07827 tlist = *trianglelist;
07828 talist = *triangleattriblist;
07829 vertexindex = 0;
07830 attribindex = 0;
07831 traversalinit(&m->triangles);
07832 triangleloop.tri = triangletraverse(m);
07833 triangleloop.orient = 0;
07834 elementnumber = b->firstnumber;
07835 while (triangleloop.tri != (triangle *) NULL) {
07836 org(triangleloop, p1);
07837 dest(triangleloop, p2);
07838 apex(triangleloop, p3);
07839 if (b->order == 1) {
07840 tlist[vertexindex++] = vertexmark(p1);
07841 tlist[vertexindex++] = vertexmark(p2);
07842 tlist[vertexindex++] = vertexmark(p3);
07843 } else {
07844 mid1 = (vertex) triangleloop.tri[m->highorderindex + 1];
07845 mid2 = (vertex) triangleloop.tri[m->highorderindex + 2];
07846 mid3 = (vertex) triangleloop.tri[m->highorderindex];
07847 tlist[vertexindex++] = vertexmark(p1);
07848 tlist[vertexindex++] = vertexmark(p2);
07849 tlist[vertexindex++] = vertexmark(p3);
07850 tlist[vertexindex++] = vertexmark(mid1);
07851 tlist[vertexindex++] = vertexmark(mid2);
07852 tlist[vertexindex++] = vertexmark(mid3);
07853 }
07854
07855 for (i = 0; i < m->eextras; i++) {
07856 talist[attribindex++] = elemattribute(triangleloop, i);
07857 }
07858 triangleloop.tri = triangletraverse(m);
07859 elementnumber++;
07860 }
07861 }
07862
07863
07864
07865
07866
07867
07868
07869 void writepoly(struct mesh *m, struct behavior *b,
07870 int **segmentlist, int **segmentmarkerlist)
07871 {
07872 int *slist;
07873 int *smlist;
07874 int index;
07875 struct osub subsegloop;
07876 vertex endpoint1, endpoint2;
07877 long subsegnumber;
07878
07879 if (!b->quiet) {
07880 printf("Writing segments.\n");
07881 }
07882
07883 if (*segmentlist == (int *) NULL) {
07884 *segmentlist = (int *) trimalloc((int) (m->subsegs.items * 2 *
07885 sizeof(int)));
07886 }
07887
07888 if (!b->nobound && (*segmentmarkerlist == (int *) NULL)) {
07889 *segmentmarkerlist = (int *) trimalloc((int) (m->subsegs.items *
07890 sizeof(int)));
07891 }
07892 slist = *segmentlist;
07893 smlist = *segmentmarkerlist;
07894 index = 0;
07895
07896 traversalinit(&m->subsegs);
07897 subsegloop.ss = subsegtraverse(m);
07898 subsegloop.ssorient = 0;
07899 subsegnumber = b->firstnumber;
07900 while (subsegloop.ss != (subseg *) NULL) {
07901 sorg(subsegloop, endpoint1);
07902 sdest(subsegloop, endpoint2);
07903
07904 slist[index++] = vertexmark(endpoint1);
07905 slist[index++] = vertexmark(endpoint2);
07906 if (!b->nobound) {
07907
07908 smlist[subsegnumber - b->firstnumber] = mark(subsegloop);
07909 }
07910 subsegloop.ss = subsegtraverse(m);
07911 subsegnumber++;
07912 }
07913 }
07914
07915
07916
07917
07918
07919
07920
07921 void writeedges(struct mesh *m, struct behavior *b,
07922 int **edgelist, int **edgemarkerlist)
07923 {
07924 int *elist;
07925 int *emlist;
07926 int index;
07927 struct otri triangleloop, trisym;
07928 struct osub checkmark;
07929 vertex p1, p2;
07930 long edgenumber;
07931 triangle ptr;
07932 subseg sptr;
07933
07934 if (!b->quiet) {
07935 printf("Writing edges.\n");
07936 }
07937
07938 if (*edgelist == (int *) NULL) {
07939 *edgelist = (int *) trimalloc((int) (m->edges * 2 * sizeof(int)));
07940 }
07941
07942 if (!b->nobound && (*edgemarkerlist == (int *) NULL)) {
07943 *edgemarkerlist = (int *) trimalloc((int) (m->edges * sizeof(int)));
07944 }
07945 elist = *edgelist;
07946 emlist = *edgemarkerlist;
07947 index = 0;
07948
07949 traversalinit(&m->triangles);
07950 triangleloop.tri = triangletraverse(m);
07951 edgenumber = b->firstnumber;
07952
07953
07954
07955
07956
07957
07958 while (triangleloop.tri != (triangle *) NULL) {
07959 for (triangleloop.orient = 0; triangleloop.orient < 3;
07960 triangleloop.orient++) {
07961 sym(triangleloop, trisym);
07962 if ((triangleloop.tri < trisym.tri) || (trisym.tri == m->dummytri)) {
07963 org(triangleloop, p1);
07964 dest(triangleloop, p2);
07965 elist[index++] = vertexmark(p1);
07966 elist[index++] = vertexmark(p2);
07967 if (b->nobound) {
07968 } else {
07969
07970
07971 if (b->usesegments) {
07972 tspivot(triangleloop, checkmark);
07973 if (checkmark.ss == m->dummysub) {
07974 emlist[edgenumber - b->firstnumber] = 0;
07975 } else {
07976 emlist[edgenumber - b->firstnumber] = mark(checkmark);
07977 }
07978 } else {
07979 emlist[edgenumber - b->firstnumber] = trisym.tri == m->dummytri;
07980 }
07981 }
07982 edgenumber++;
07983 }
07984 }
07985 triangleloop.tri = triangletraverse(m);
07986 }
07987 }
07988
07989
07990
07991
07992
07993
07994
07995
07996
07997
07998
07999
08000
08001
08002
08003
08004
08005 void writevoronoi(struct mesh *m, struct behavior *b, float **vpointlist,
08006 float **vpointattriblist, int **vpointmarkerlist,
08007 int **vedgelist, int **vedgemarkerlist, float **vnormlist)
08008 {
08009 float *plist;
08010 float *palist;
08011 int *elist;
08012 float *normlist;
08013 int coordindex;
08014 int attribindex;
08015 struct otri triangleloop, trisym;
08016 vertex torg, tdest, tapex;
08017 float circumcenter[2];
08018 float xi, eta;
08019 long vnodenumber, vedgenumber;
08020 int p1, p2;
08021 int i;
08022 triangle ptr;
08023
08024 if (!b->quiet) {
08025 printf("Writing Voronoi vertices.\n");
08026 }
08027
08028 if (*vpointlist == (float *) NULL) {
08029 *vpointlist = (float *) trimalloc((int) (m->triangles.items * 2 *
08030 sizeof(float)));
08031 }
08032
08033 if (*vpointattriblist == (float *) NULL) {
08034 *vpointattriblist = (float *) trimalloc((int) (m->triangles.items *
08035 m->nextras * sizeof(float)));
08036 }
08037 *vpointmarkerlist = (int *) NULL;
08038 plist = *vpointlist;
08039 palist = *vpointattriblist;
08040 coordindex = 0;
08041 attribindex = 0;
08042
08043 traversalinit(&m->triangles);
08044 triangleloop.tri = triangletraverse(m);
08045 triangleloop.orient = 0;
08046 vnodenumber = b->firstnumber;
08047 while (triangleloop.tri != (triangle *) NULL) {
08048 org(triangleloop, torg);
08049 dest(triangleloop, tdest);
08050 apex(triangleloop, tapex);
08051 findcircumcenter(m, b, torg, tdest, tapex, circumcenter, &xi, &eta, 0);
08052
08053
08054 plist[coordindex++] = circumcenter[0];
08055 plist[coordindex++] = circumcenter[1];
08056 for (i = 2; i < 2 + m->nextras; i++) {
08057
08058 palist[attribindex++] = torg[i] + xi * (tdest[i] - torg[i])
08059 + eta * (tapex[i] - torg[i]);
08060 }
08061
08062 * (int *) (triangleloop.tri + 6) = (int) vnodenumber;
08063 triangleloop.tri = triangletraverse(m);
08064 vnodenumber++;
08065 }
08066
08067 if (!b->quiet) {
08068 printf("Writing Voronoi edges.\n");
08069 }
08070
08071 if (*vedgelist == (int *) NULL) {
08072 *vedgelist = (int *) trimalloc((int) (m->edges * 2 * sizeof(int)));
08073 }
08074 *vedgemarkerlist = (int *) NULL;
08075
08076 if (*vnormlist == (float *) NULL) {
08077 *vnormlist = (float *) trimalloc((int) (m->edges * 2 * sizeof(float)));
08078 }
08079 elist = *vedgelist;
08080 normlist = *vnormlist;
08081 coordindex = 0;
08082
08083 traversalinit(&m->triangles);
08084 triangleloop.tri = triangletraverse(m);
08085 vedgenumber = b->firstnumber;
08086
08087
08088
08089
08090
08091
08092 while (triangleloop.tri != (triangle *) NULL) {
08093 for (triangleloop.orient = 0; triangleloop.orient < 3;
08094 triangleloop.orient++) {
08095 sym(triangleloop, trisym);
08096 if ((triangleloop.tri < trisym.tri) || (trisym.tri == m->dummytri)) {
08097
08098 p1 = * (int *) (triangleloop.tri + 6);
08099 if (trisym.tri == m->dummytri) {
08100 org(triangleloop, torg);
08101 dest(triangleloop, tdest);
08102
08103 elist[coordindex] = p1;
08104 normlist[coordindex++] = tdest[1] - torg[1];
08105 elist[coordindex] = -1;
08106 normlist[coordindex++] = torg[0] - tdest[0];
08107 } else {
08108
08109 p2 = * (int *) (trisym.tri + 6);
08110
08111 elist[coordindex] = p1;
08112 normlist[coordindex++] = 0.0;
08113 elist[coordindex] = p2;
08114 normlist[coordindex++] = 0.0;
08115 }
08116 vedgenumber++;
08117 }
08118 }
08119 triangleloop.tri = triangletraverse(m);
08120 }
08121 }
08122
08123
08124 void writeneighbors(struct mesh *m, struct behavior *b, int **neighborlist)
08125 {
08126 int *nlist;
08127 int index;
08128 struct otri triangleloop, trisym;
08129 long elementnumber;
08130 int neighbor1, neighbor2, neighbor3;
08131 triangle ptr;
08132
08133 if (!b->quiet) {
08134 printf("Writing neighbors.\n");
08135 }
08136
08137 if (*neighborlist == (int *) NULL) {
08138 *neighborlist = (int *) trimalloc((int) (m->triangles.items * 3 *
08139 sizeof(int)));
08140 }
08141 nlist = *neighborlist;
08142 index = 0;
08143
08144 traversalinit(&m->triangles);
08145 triangleloop.tri = triangletraverse(m);
08146 triangleloop.orient = 0;
08147 elementnumber = b->firstnumber;
08148 while (triangleloop.tri != (triangle *) NULL) {
08149 * (int *) (triangleloop.tri + 6) = (int) elementnumber;
08150 triangleloop.tri = triangletraverse(m);
08151 elementnumber++;
08152 }
08153 * (int *) (m->dummytri + 6) = -1;
08154
08155 traversalinit(&m->triangles);
08156 triangleloop.tri = triangletraverse(m);
08157 elementnumber = b->firstnumber;
08158 while (triangleloop.tri != (triangle *) NULL) {
08159 triangleloop.orient = 1;
08160 sym(triangleloop, trisym);
08161 neighbor1 = * (int *) (trisym.tri + 6);
08162 triangleloop.orient = 2;
08163 sym(triangleloop, trisym);
08164 neighbor2 = * (int *) (trisym.tri + 6);
08165 triangleloop.orient = 0;
08166 sym(triangleloop, trisym);
08167 neighbor3 = * (int *) (trisym.tri + 6);
08168 nlist[index++] = neighbor1;
08169 nlist[index++] = neighbor2;
08170 nlist[index++] = neighbor3;
08171
08172 triangleloop.tri = triangletraverse(m);
08173 elementnumber++;
08174 }
08175 }
08176
08179
08180
08181
08182
08183
08184
08185
08186
08187 void quality_statistics(struct mesh *m, struct behavior *b)
08188 {
08189 struct otri triangleloop;
08190 vertex p[3];
08191 float cossquaretable[8];
08192 float ratiotable[16];
08193 float dx[3], dy[3];
08194 float edgelength[3];
08195 float dotproduct;
08196 float cossquare;
08197 float triarea;
08198 float shortest, longest;
08199 float trilongest2;
08200 float smallestarea, biggestarea;
08201 float triminaltitude2;
08202 float minaltitude;
08203 float triaspect2;
08204 float worstaspect;
08205 float smallestangle, biggestangle;
08206 float radconst, degconst;
08207 int angletable[18];
08208 int aspecttable[16];
08209 int aspectindex;
08210 int tendegree;
08211 int acutebiggest;
08212 int i, ii, j, k;
08213
08214 printf("Mesh quality statistics:\n\n");
08215 radconst = PI / 18.0;
08216 degconst = 180.0 / PI;
08217 for (i = 0; i < 8; i++) {
08218 cossquaretable[i] = cos(radconst * (float) (i + 1));
08219 cossquaretable[i] = cossquaretable[i] * cossquaretable[i];
08220 }
08221 for (i = 0; i < 18; i++) {
08222 angletable[i] = 0;
08223 }
08224
08225 ratiotable[0] = 1.5; ratiotable[1] = 2.0;
08226 ratiotable[2] = 2.5; ratiotable[3] = 3.0;
08227 ratiotable[4] = 4.0; ratiotable[5] = 6.0;
08228 ratiotable[6] = 10.0; ratiotable[7] = 15.0;
08229 ratiotable[8] = 25.0; ratiotable[9] = 50.0;
08230 ratiotable[10] = 100.0; ratiotable[11] = 300.0;
08231 ratiotable[12] = 1000.0; ratiotable[13] = 10000.0;
08232 ratiotable[14] = 100000.0; ratiotable[15] = 0.0;
08233 for (i = 0; i < 16; i++) {
08234 aspecttable[i] = 0;
08235 }
08236
08237 worstaspect = 0.0;
08238 minaltitude = m->xmax - m->xmin + m->ymax - m->ymin;
08239 minaltitude = minaltitude * minaltitude;
08240 shortest = minaltitude;
08241 longest = 0.0;
08242 smallestarea = minaltitude;
08243 biggestarea = 0.0;
08244 worstaspect = 0.0;
08245 smallestangle = 0.0;
08246 biggestangle = 2.0;
08247 acutebiggest = 1;
08248
08249 traversalinit(&m->triangles);
08250 triangleloop.tri = triangletraverse(m);
08251 triangleloop.orient = 0;
08252 while (triangleloop.tri != (triangle *) NULL) {
08253 org(triangleloop, p[0]);
08254 dest(triangleloop, p[1]);
08255 apex(triangleloop, p[2]);
08256 trilongest2 = 0.0;
08257
08258 for (i = 0; i < 3; i++) {
08259 j = plus1mod3[i];
08260 k = minus1mod3[i];
08261 dx[i] = p[j][0] - p[k][0];
08262 dy[i] = p[j][1] - p[k][1];
08263 edgelength[i] = dx[i] * dx[i] + dy[i] * dy[i];
08264 if (edgelength[i] > trilongest2) {
08265 trilongest2 = edgelength[i];
08266 }
08267 if (edgelength[i] > longest) {
08268 longest = edgelength[i];
08269 }
08270 if (edgelength[i] < shortest) {
08271 shortest = edgelength[i];
08272 }
08273 }
08274
08275 triarea = counterclockwise(m, b, p[0], p[1], p[2]);
08276 if (triarea < smallestarea) {
08277 smallestarea = triarea;
08278 }
08279 if (triarea > biggestarea) {
08280 biggestarea = triarea;
08281 }
08282 triminaltitude2 = triarea * triarea / trilongest2;
08283 if (triminaltitude2 < minaltitude) {
08284 minaltitude = triminaltitude2;
08285 }
08286 triaspect2 = trilongest2 / triminaltitude2;
08287 if (triaspect2 > worstaspect) {
08288 worstaspect = triaspect2;
08289 }
08290 aspectindex = 0;
08291 while ((triaspect2 > ratiotable[aspectindex] * ratiotable[aspectindex])
08292 && (aspectindex < 15)) {
08293 aspectindex++;
08294 }
08295 aspecttable[aspectindex]++;
08296
08297 for (i = 0; i < 3; i++) {
08298 j = plus1mod3[i];
08299 k = minus1mod3[i];
08300 dotproduct = dx[j] * dx[k] + dy[j] * dy[k];
08301 cossquare = dotproduct * dotproduct / (edgelength[j] * edgelength[k]);
08302 tendegree = 8;
08303 for (ii = 7; ii >= 0; ii--) {
08304 if (cossquare > cossquaretable[ii]) {
08305 tendegree = ii;
08306 }
08307 }
08308 if (dotproduct <= 0.0) {
08309 angletable[tendegree]++;
08310 if (cossquare > smallestangle) {
08311 smallestangle = cossquare;
08312 }
08313 if (acutebiggest && (cossquare < biggestangle)) {
08314 biggestangle = cossquare;
08315 }
08316 } else {
08317 angletable[17 - tendegree]++;
08318 if (acutebiggest || (cossquare > biggestangle)) {
08319 biggestangle = cossquare;
08320 acutebiggest = 0;
08321 }
08322 }
08323 }
08324 triangleloop.tri = triangletraverse(m);
08325 }
08326
08327 shortest = sqrt(shortest);
08328 longest = sqrt(longest);
08329 minaltitude = sqrt(minaltitude);
08330 worstaspect = sqrt(worstaspect);
08331 smallestarea *= 0.5;
08332 biggestarea *= 0.5;
08333 if (smallestangle >= 1.0) {
08334 smallestangle = 0.0;
08335 } else {
08336 smallestangle = degconst * acos(sqrt(smallestangle));
08337 }
08338 if (biggestangle >= 1.0) {
08339 biggestangle = 180.0;
08340 } else {
08341 if (acutebiggest) {
08342 biggestangle = degconst * acos(sqrt(biggestangle));
08343 } else {
08344 biggestangle = 180.0 - degconst * acos(sqrt(biggestangle));
08345 }
08346 }
08347
08348 printf(" Smallest area: %16.5g | Largest area: %16.5g\n",
08349 smallestarea, biggestarea);
08350 printf(" Shortest edge: %16.5g | Longest edge: %16.5g\n",
08351 shortest, longest);
08352 printf(" Shortest altitude: %12.5g | Largest aspect ratio: %8.5g\n\n",
08353 minaltitude, worstaspect);
08354
08355 printf(" Triangle aspect ratio histogram:\n");
08356 printf(" 1.1547 - %-6.6g : %8d | %6.6g - %-6.6g : %8d\n",
08357 ratiotable[0], aspecttable[0], ratiotable[7], ratiotable[8],
08358 aspecttable[8]);
08359 for (i = 1; i < 7; i++) {
08360 printf(" %6.6g - %-6.6g : %8d | %6.6g - %-6.6g : %8d\n",
08361 ratiotable[i - 1], ratiotable[i], aspecttable[i],
08362 ratiotable[i + 7], ratiotable[i + 8], aspecttable[i + 8]);
08363 }
08364 printf(" %6.6g - %-6.6g : %8d | %6.6g - : %8d\n",
08365 ratiotable[6], ratiotable[7], aspecttable[7], ratiotable[14],
08366 aspecttable[15]);
08367 printf(" (Aspect ratio is longest edge divided by shortest altitude)\n\n");
08368
08369 printf(" Smallest angle: %15.5g | Largest angle: %15.5g\n\n",
08370 smallestangle, biggestangle);
08371
08372 printf(" Angle histogram:\n");
08373 for (i = 0; i < 9; i++) {
08374 printf(" %3d - %3d degrees: %8d | %3d - %3d degrees: %8d\n",
08375 i * 10, i * 10 + 10, angletable[i],
08376 i * 10 + 90, i * 10 + 100, angletable[i + 9]);
08377 }
08378 printf("\n");
08379 }
08380
08381
08382
08383
08384
08385
08386
08387 void statistics(struct mesh *m, struct behavior *b)
08388 {
08389 printf("\nStatistics:\n\n");
08390 printf(" Input vertices: %d\n", m->invertices);
08391 if (b->refine) {
08392 printf(" Input triangles: %d\n", m->inelements);
08393 }
08394 if (b->poly) {
08395 printf(" Input segments: %d\n", m->insegments);
08396 if (!b->refine) {
08397 printf(" Input holes: %d\n", m->holes);
08398 }
08399 }
08400
08401 printf("\n Mesh vertices: %ld\n", m->vertices.items - m->undeads);
08402 printf(" Mesh triangles: %ld\n", m->triangles.items);
08403 printf(" Mesh edges: %ld\n", m->edges);
08404 printf(" Mesh exterior boundary edges: %ld\n", m->hullsize);
08405 if (b->poly || b->refine) {
08406 printf(" Mesh interior boundary edges: %ld\n",
08407 m->subsegs.items - m->hullsize);
08408 printf(" Mesh subsegments (constrained edges): %ld\n",
08409 m->subsegs.items);
08410 }
08411 printf("\n");
08412
08413 if (b->verbose) {
08414 quality_statistics(m, b);
08415 printf("Memory allocation statistics:\n\n");
08416 printf(" Maximum number of vertices: %ld\n", m->vertices.maxitems);
08417 printf(" Maximum number of triangles: %ld\n", m->triangles.maxitems);
08418 if (m->subsegs.maxitems > 0) {
08419 printf(" Maximum number of subsegments: %ld\n", m->subsegs.maxitems);
08420 }
08421 if (m->viri.maxitems > 0) {
08422 printf(" Maximum number of viri: %ld\n", m->viri.maxitems);
08423 }
08424 if (m->badsubsegs.maxitems > 0) {
08425 printf(" Maximum number of encroached subsegments: %ld\n",
08426 m->badsubsegs.maxitems);
08427 }
08428 if (m->badtriangles.maxitems > 0) {
08429 printf(" Maximum number of bad triangles: %ld\n",
08430 m->badtriangles.maxitems);
08431 }
08432 if (m->flipstackers.maxitems > 0) {
08433 printf(" Maximum number of stacked triangle flips: %ld\n",
08434 m->flipstackers.maxitems);
08435 }
08436 if (m->splaynodes.maxitems > 0) {
08437 printf(" Maximum number of splay tree nodes: %ld\n",
08438 m->splaynodes.maxitems);
08439 }
08440 printf(" Approximate heap memory use (bytes): %ld\n\n",
08441 m->vertices.maxitems * m->vertices.itembytes +
08442 m->triangles.maxitems * m->triangles.itembytes +
08443 m->subsegs.maxitems * m->subsegs.itembytes +
08444 m->viri.maxitems * m->viri.itembytes +
08445 m->badsubsegs.maxitems * m->badsubsegs.itembytes +
08446 m->badtriangles.maxitems * m->badtriangles.itembytes +
08447 m->flipstackers.maxitems * m->flipstackers.itembytes +
08448 m->splaynodes.maxitems * m->splaynodes.itembytes);
08449
08450 printf("Algorithmic statistics:\n\n");
08451 if (!b->weighted) {
08452 printf(" Number of incircle tests: %ld\n", m->incirclecount);
08453 } else {
08454 printf(" Number of 3D orientation tests: %ld\n", m->orient3dcount);
08455 }
08456 printf(" Number of 2D orientation tests: %ld\n", m->counterclockcount);
08457 if (m->hyperbolacount > 0) {
08458 printf(" Number of right-of-hyperbola tests: %ld\n",
08459 m->hyperbolacount);
08460 }
08461 if (m->circletopcount > 0) {
08462 printf(" Number of circle top computations: %ld\n",
08463 m->circletopcount);
08464 }
08465 if (m->circumcentercount > 0) {
08466 printf(" Number of triangle circumcenter computations: %ld\n",
08467 m->circumcentercount);
08468 }
08469 printf("\n");
08470 }
08471 }
08472
08473
08474
08475
08476
08477
08478
08479
08480
08481
08482
08483
08484
08485
08486
08487
08488
08489
08490
08491
08492
08493
08494
08495
08496
08497
08498 void triangulate(char *triswitches, struct triangulateio *in,
08499 struct triangulateio *out, struct triangulateio *vorout)
08500 {
08501 struct mesh m;
08502 struct behavior b;
08503 float *holearray;
08504 float *regionarray;
08505
08506 triangleinit(&m);
08507 parsecommandline(1, &triswitches, &b);
08508 m.steinerleft = b.steiner;
08509
08510 transfernodes(&m, &b, in->pointlist, in->pointattributelist,
08511 in->pointmarkerlist, in->numberofpoints,
08512 in->numberofpointattributes);
08513
08514 m.hullsize = delaunay(&m, &b);
08515
08516
08517 m.infvertex1 = (vertex) NULL;
08518 m.infvertex2 = (vertex) NULL;
08519 m.infvertex3 = (vertex) NULL;
08520
08521 if (b.usesegments) {
08522 m.checksegments = 1;
08523 if (!b.refine) {
08524
08525 formskeleton(&m, &b, in->segmentlist,
08526 in->segmentmarkerlist, in->numberofsegments);
08527 }
08528 }
08529
08530 if (b.poly && (m.triangles.items > 0)) {
08531 holearray = in->holelist;
08532 m.holes = in->numberofholes;
08533 regionarray = in->regionlist;
08534 m.regions = in->numberofregions;
08535 if (!b.refine) {
08536
08537 carveholes(&m, &b, holearray, m.holes, regionarray, m.regions);
08538 }
08539 } else {
08540
08541
08542
08543 m.holes = 0;
08544 m.regions = 0;
08545 }
08546
08547
08548 m.edges = (3l * m.triangles.items + m.hullsize) / 2l;
08549
08550 if (b.order > 1) {
08551 highorder(&m, &b);
08552 }
08553 if (!b.quiet) {
08554 printf("\n");
08555 }
08556
08557 if (b.jettison) {
08558 out->numberofpoints = m.vertices.items - m.undeads;
08559 } else {
08560 out->numberofpoints = m.vertices.items;
08561 }
08562 out->numberofpointattributes = m.nextras;
08563 out->numberoftriangles = m.triangles.items;
08564 out->numberofcorners = (b.order + 1) * (b.order + 2) / 2;
08565 out->numberoftriangleattributes = m.eextras;
08566 out->numberofedges = m.edges;
08567 if (b.usesegments) {
08568 out->numberofsegments = m.subsegs.items;
08569 } else {
08570 out->numberofsegments = m.hullsize;
08571 }
08572 if (vorout != (struct triangulateio *) NULL) {
08573 vorout->numberofpoints = m.triangles.items;
08574 vorout->numberofpointattributes = m.nextras;
08575 vorout->numberofedges = m.edges;
08576 }
08577
08578
08579 if (b.nonodewritten || (b.noiterationnum && m.readnodefile)) {
08580 if (!b.quiet) {
08581 printf("NOT writing vertices.\n");
08582 }
08583 numbernodes(&m, &b);
08584 } else {
08585
08586 writenodes(&m, &b, &out->pointlist, &out->pointattributelist,
08587 &out->pointmarkerlist);
08588 }
08589 if (b.noelewritten) {
08590 if (!b.quiet) {
08591 printf("NOT writing triangles.\n");
08592 }
08593 } else {
08594 writeelements(&m, &b, &out->trianglelist, &out->triangleattributelist);
08595 }
08596
08597
08598 if (b.poly || b.convex) {
08599
08600 if (b.nopolywritten || b.noiterationnum) {
08601 if (!b.quiet) {
08602 printf("NOT writing segments.\n");
08603 }
08604 } else {
08605 writepoly(&m, &b, &out->segmentlist, &out->segmentmarkerlist);
08606 out->numberofholes = m.holes;
08607 out->numberofregions = m.regions;
08608 if (b.poly) {
08609 out->holelist = in->holelist;
08610 out->regionlist = in->regionlist;
08611 } else {
08612 out->holelist = (float *) NULL;
08613 out->regionlist = (float *) NULL;
08614 }
08615 }
08616 }
08617 if (b.edgesout) {
08618 writeedges(&m, &b, &out->edgelist, &out->edgemarkerlist);
08619 }
08620 if (b.voronoi) {
08621 writevoronoi(&m, &b, &vorout->pointlist, &vorout->pointattributelist,
08622 &vorout->pointmarkerlist, &vorout->edgelist,
08623 &vorout->edgemarkerlist, &vorout->normlist);
08624 }
08625 if (b.neighbors) {
08626 writeneighbors(&m, &b, &out->neighborlist);
08627 }
08628
08629 if (!b.quiet) {
08630 statistics(&m, &b);
08631 }
08632
08633 triangledeinit(&m, &b);
08634 }