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