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