00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014 #include "qhull_a.h"
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 void qh_backnormal (realT **rows, int numrow, int numcol, boolT sign,
00054 coordT *normal, boolT *nearzero) {
00055 int i, j;
00056 coordT *normalp, *normal_tail, *ai, *ak;
00057 realT diagonal;
00058 boolT waszero;
00059 int zerocol= -1;
00060
00061 normalp= normal + numcol - 1;
00062 *normalp--= (sign ? -1.0 : 1.0);
00063 for(i= numrow; i--; ) {
00064 *normalp= 0.0;
00065 ai= rows[i] + i + 1;
00066 ak= normalp+1;
00067 for(j= i+1; j < numcol; j++)
00068 *normalp -= *ai++ * *ak++;
00069 diagonal= (rows[i])[i];
00070 if (fabs_(diagonal) > qh MINdenom_2)
00071 *(normalp--) /= diagonal;
00072 else {
00073 waszero= False;
00074 *normalp= qh_divzero (*normalp, diagonal, qh MINdenom_1_2, &waszero);
00075 if (waszero) {
00076 zerocol= i;
00077 *(normalp--)= (sign ? -1.0 : 1.0);
00078 for (normal_tail= normalp+2; normal_tail < normal + numcol; normal_tail++)
00079 *normal_tail= 0.0;
00080 }else
00081 normalp--;
00082 }
00083 }
00084 if (zerocol != -1) {
00085 zzinc_(Zback0);
00086 *nearzero= True;
00087 trace4((qh ferr, "qh_backnormal: zero diagonal at column %d.\n", i));
00088 qh_precision ("zero diagonal on back substitution");
00089 }
00090 }
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109 void qh_distplane (pointT *point, facetT *facet, realT *dist) {
00110 coordT *normal= facet->normal, *coordp, randr;
00111 int k;
00112
00113 switch(qh hull_dim){
00114 case 2:
00115 *dist= facet->offset + point[0] * normal[0] + point[1] * normal[1];
00116 break;
00117 case 3:
00118 *dist= facet->offset + point[0] * normal[0] + point[1] * normal[1] + point[2] * normal[2];
00119 break;
00120 case 4:
00121 *dist= facet->offset+point[0]*normal[0]+point[1]*normal[1]+point[2]*normal[2]+point[3]*normal[3];
00122 break;
00123 case 5:
00124 *dist= facet->offset+point[0]*normal[0]+point[1]*normal[1]+point[2]*normal[2]+point[3]*normal[3]+point[4]*normal[4];
00125 break;
00126 case 6:
00127 *dist= facet->offset+point[0]*normal[0]+point[1]*normal[1]+point[2]*normal[2]+point[3]*normal[3]+point[4]*normal[4]+point[5]*normal[5];
00128 break;
00129 case 7:
00130 *dist= facet->offset+point[0]*normal[0]+point[1]*normal[1]+point[2]*normal[2]+point[3]*normal[3]+point[4]*normal[4]+point[5]*normal[5]+point[6]*normal[6];
00131 break;
00132 case 8:
00133 *dist= facet->offset+point[0]*normal[0]+point[1]*normal[1]+point[2]*normal[2]+point[3]*normal[3]+point[4]*normal[4]+point[5]*normal[5]+point[6]*normal[6]+point[7]*normal[7];
00134 break;
00135 default:
00136 *dist= facet->offset;
00137 coordp= point;
00138 for (k= qh hull_dim; k--; )
00139 *dist += *coordp++ * *normal++;
00140 break;
00141 }
00142 zinc_(Zdistplane);
00143 if (!qh RANDOMdist && qh IStracing < 4)
00144 return;
00145 if (qh RANDOMdist) {
00146 randr= qh_RANDOMint;
00147 *dist += (2.0 * randr / qh_RANDOMmax - 1.0) *
00148 qh RANDOMfactor * qh MAXabs_coord;
00149 }
00150 if (qh IStracing >= 4) {
00151 fprintf (qh ferr, "qh_distplane: ");
00152 fprintf (qh ferr, qh_REAL_1, *dist);
00153 fprintf (qh ferr, "from p%d to f%d\n", qh_pointid(point), facet->id);
00154 }
00155 return;
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
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227 facetT *qh_findbest (pointT *point, facetT *startfacet,
00228 boolT bestoutside, boolT newfacets, boolT noupper,
00229 realT *dist, boolT *isoutside, int *numpart) {
00230 realT bestdist= -REALmax/2 , searchdist;
00231 realT cutoff, mincutoff;
00232 facetT *facet, *neighbor, **neighborp, *bestfacet= NULL;
00233 int oldtrace= qh IStracing;
00234 int searchsize= 0;
00235 boolT newbest;
00236 boolT ischeckmax= bestoutside && !newfacets && !isoutside;
00237 boolT ispartition= newfacets && isoutside;
00238 boolT isfindfacet= !newfacets && isoutside;
00239 boolT testhorizon = ispartition && (bestoutside || qh APPROXhull || qh MERGING);
00240
00241 if (!ischeckmax && !ispartition && !isfindfacet) {
00242 fprintf (qh ferr, "qhull internal error (qh_findbest): unknown combination of arguments\n");
00243 qh_errexit (qh_ERRqhull, startfacet, NULL);
00244 }
00245 if (qh TRACElevel && qh TRACEpoint >= 0 && qh TRACEpoint == qh_pointid (point)) {
00246 qh IStracing= qh TRACElevel;
00247 fprintf (qh ferr, "qh_findbest: point p%d starting at f%d bestoutside? %d newfacets %d\n",
00248 qh TRACEpoint, startfacet->id, bestoutside, newfacets);
00249 fprintf(qh ferr, " ischeckmax %d ispartition %d isfindfacet %d testhorizon %d\n",
00250 ischeckmax, ispartition, isfindfacet, testhorizon);
00251 fprintf (qh ferr, " Last point added to hull was p%d.", qh furthest_id);
00252 fprintf(qh ferr, " Last merge was #%d.\n", zzval_(Ztotmerge));
00253 }
00254 if (isoutside)
00255 *isoutside= True;
00256
00257 if (!startfacet->flipped) {
00258 *numpart= 1;
00259 qh_distplane (point, startfacet, dist);
00260 if (!startfacet->upperdelaunay || (!noupper && *dist >= qh MINoutside)) {
00261 bestdist= *dist;
00262 bestfacet= startfacet;
00263 if (!bestoutside && *dist >= qh MINoutside)
00264 goto LABELreturn_best;
00265 }
00266 #if qh_MAXoutside
00267 if (ischeckmax && (!qh ONLYgood || startfacet->good) && *dist > startfacet->maxoutside)
00268 startfacet->maxoutside= *dist;
00269 #endif
00270 }
00271 if (ispartition)
00272 searchdist= 2 * qh DISTround;
00273 else
00274 searchdist= qh max_outside + 2 * qh DISTround
00275 + fmax_( qh MINvisible, qh MAXcoplanar);
00276 cutoff= bestdist - searchdist;
00277 mincutoff= 0;
00278 if (ischeckmax) {
00279 mincutoff= -(qh DISTround - fmax_(qh MINvisible, qh MAXcoplanar));
00280 minimize_(cutoff, mincutoff);
00281 }
00282 startfacet->visitid= ++qh visit_id;
00283 facet= startfacet;
00284 do {
00285 if (True) {
00286 LABELrestart:
00287 newbest= False;
00288 trace4((qh ferr, "qh_findbest: neighbors of f%d, bestdist %2.2g cutoff %2.2g searchdist %2.2g\n",
00289 facet->id, bestdist, cutoff, searchdist));
00290 FOREACHneighbor_(facet) {
00291 if (ispartition && !neighbor->newfacet)
00292 continue;
00293 if (!neighbor->flipped) {
00294 if (neighbor->visitid == qh visit_id)
00295 continue;
00296 neighbor->visitid= qh visit_id;
00297 (*numpart)++;
00298 qh_distplane (point, neighbor, dist);
00299 if (!bestoutside && *dist >= qh MINoutside
00300 && (!noupper || !facet->upperdelaunay)) {
00301 bestfacet= neighbor;
00302 goto LABELreturn_best;
00303 }
00304 #if qh_MAXoutside
00305 if (ischeckmax) {
00306 if ((!qh ONLYgood || neighbor->good)
00307 && *dist > neighbor->maxoutside)
00308 neighbor->maxoutside= *dist;
00309 else if (bestfacet && *dist < cutoff)
00310 continue;
00311 }else
00312 #endif
00313 if (bestfacet && *dist < cutoff)
00314 continue;
00315 if (*dist > bestdist) {
00316 if (!neighbor->upperdelaunay
00317 || (bestoutside && !noupper && *dist >= qh MINoutside)) {
00318 if (ischeckmax && qh_MAXoutside) {
00319 bestdist= *dist;
00320 bestfacet= neighbor;
00321 cutoff= bestdist - searchdist;
00322 minimize_(cutoff, mincutoff);
00323 }else if (*dist > bestdist + searchdist) {
00324 bestdist= *dist;
00325 bestfacet= neighbor;
00326 cutoff= bestdist - searchdist;
00327 searchsize= 0;
00328 facet= neighbor;
00329 if (newbest)
00330 facet->visitid= ++qh visit_id;
00331 goto LABELrestart;
00332 }else {
00333 bestdist= *dist;
00334 bestfacet= neighbor;
00335 cutoff= bestdist - searchdist;
00336 }
00337 newbest= True;
00338 }
00339 }
00340 }
00341 if (!searchsize++) {
00342 SETfirst_(qh searchset) = neighbor;
00343 qh_settruncate (qh searchset, 1);
00344 }else
00345 qh_setappend (&qh searchset, neighbor);
00346 }
00347 }
00348 }while
00349 (searchsize && (facet= (facetT*)qh_setdellast (qh searchset)));
00350 if (!ischeckmax) {
00351 if (!bestfacet) {
00352 fprintf (qh ferr, "qh_findbest: point p%d starting at f%d bestoutside? %d newfacets %d\n",
00353 qh TRACEpoint, startfacet->id, bestoutside, newfacets);
00354 fprintf(qh ferr, "\n\
00355 qh_findbest: all neighbors of facet %d are flipped or upper Delaunay.\n\
00356 Please report this error to qhull_bug@geom.umn.edu with the input and all of the output.\n",
00357 startfacet->id);
00358 qh FORCEoutput= True;
00359 qh_errexit (qh_ERRqhull, startfacet, NULL);
00360 }
00361 if (ispartition && !qh findbest_notsharp && bestdist < - qh DISTround) {
00362 if (qh_findbestsharp ( point, &bestfacet, &bestdist, numpart))
00363 qh findbestnew= True;
00364 else
00365 qh findbest_notsharp= True;
00366 }
00367 if (testhorizon) {
00368 facet= SETfirstt_(bestfacet->neighbors, facetT);
00369 trace4((qh ferr, "qh_findbest: horizon facet f%d\n", facet->id));
00370 (*numpart)++;
00371 qh_distplane (point, facet, dist);
00372 if (*dist > bestdist
00373 && (!facet->upperdelaunay || (!noupper && *dist >= qh MINoutside))) {
00374 bestdist= *dist;
00375 bestfacet= facet;
00376 }
00377 }
00378 }
00379 *dist= bestdist;
00380 if (isoutside && bestdist < qh MINoutside)
00381 *isoutside= False;
00382 LABELreturn_best:
00383 qh IStracing= oldtrace;
00384 return bestfacet;
00385 }
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423 facetT *qh_findbestnew (pointT *point, facetT *startfacet,
00424 realT *dist, boolT *isoutside, int *numpart) {
00425 realT bestdist= -REALmax, bestdist2= -REALmax;
00426 facetT *neighbor, **neighborp, *bestfacet= NULL, *newfacet, *facet;
00427 facetT *bestfacet2= NULL;
00428 int oldtrace= qh IStracing, i;
00429 realT distoutside;
00430
00431 if (!startfacet) {
00432 if (qh MERGING)
00433 fprintf(qh ferr, "qhull precision error (qh_findbestnew): merging has formed and deleted an independent cycle of facets. Can not continue.\n");
00434 else
00435 fprintf(qh ferr, "qhull internal error (qh_findbestnew): no new facets for point p%d\n",
00436 qh furthest_id);
00437 qh_errexit (qh_ERRqhull, NULL, NULL);
00438 }
00439 if (qh BESToutside || !isoutside)
00440 distoutside= REALmax;
00441 else if (qh MERGING)
00442 distoutside= qh_DISToutside;
00443 else
00444 distoutside= qh MINoutside;
00445 if (qh TRACElevel && qh TRACEpoint >= 0 && qh TRACEpoint == qh_pointid (point)) {
00446 qh IStracing= qh TRACElevel;
00447 fprintf(qh ferr, "qh_findbestnew: point p%d facet f%d. Stop if dist > %2.2g\n",
00448 qh TRACEpoint, startfacet->id, distoutside);
00449 fprintf(qh ferr, " Last point added to hull was p%d.", qh furthest_id);
00450 fprintf(qh ferr, " Last merge was #%d.\n", zzval_(Ztotmerge));
00451 }
00452 if (isoutside)
00453 *isoutside= True;
00454 *numpart= 0;
00455
00456
00457 for (i= 0, facet= startfacet; i < 2; i++, facet= qh newfacet_list) {
00458 FORALLfacet_(facet) {
00459 if (facet == startfacet && i)
00460 break;
00461 qh_distplane (point, facet, dist);
00462 (*numpart)++;
00463 if (facet->upperdelaunay) {
00464 if (*dist > bestdist2) {
00465 bestdist2= *dist;
00466 bestfacet2= facet;
00467 if (*dist >= distoutside) {
00468 bestfacet= facet;
00469 goto LABELreturn_bestnew;
00470 }
00471 }
00472 }else if (*dist > bestdist) {
00473 bestdist= *dist;
00474 bestfacet= facet;
00475 if (*dist >= distoutside)
00476 goto LABELreturn_bestnew;
00477 }
00478 }
00479 }
00480 newfacet= bestfacet ? bestfacet : bestfacet2;
00481
00482 FOREACHneighbor_(newfacet) {
00483 if (!neighbor->newfacet) {
00484 qh_distplane (point, neighbor, dist);
00485 (*numpart)++;
00486 if (neighbor->upperdelaunay) {
00487 if (*dist > bestdist2) {
00488 bestdist2= *dist;
00489 bestfacet2= neighbor;
00490 }
00491 }else if (*dist > bestdist) {
00492 bestdist= *dist;
00493 bestfacet= neighbor;
00494 }
00495 }
00496 }
00497 if (!bestfacet
00498 || (isoutside && bestdist2 >= qh MINoutside && bestdist2 > bestdist)) {
00499 *dist= bestdist2;
00500 bestfacet= bestfacet2;
00501 }else
00502 *dist= bestdist;
00503 if (isoutside && *dist < qh MINoutside)
00504 *isoutside= False;
00505 LABELreturn_bestnew:
00506 qh IStracing= oldtrace;
00507 return bestfacet;
00508 }
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530
00531
00532 void qh_gausselim(realT **rows, int numrow, int numcol, boolT *sign, boolT *nearzero) {
00533 realT *ai, *ak, *rowp, *pivotrow;
00534 realT n, pivot, pivot_abs= 0.0, temp;
00535 int i, j, k, pivoti, flip=0;
00536
00537 *nearzero= False;
00538 for(k= 0; k < numrow; k++) {
00539 pivot_abs= fabs_((rows[k])[k]);
00540 pivoti= k;
00541 for(i= k+1; i < numrow; i++) {
00542 if ((temp= fabs_((rows[i])[k])) > pivot_abs) {
00543 pivot_abs= temp;
00544 pivoti= i;
00545 }
00546 }
00547 if (pivoti != k) {
00548 rowp= rows[pivoti];
00549 rows[pivoti]= rows[k];
00550 rows[k]= rowp;
00551 *sign ^= 1;
00552 flip ^= 1;
00553 }
00554 if (pivot_abs <= qh NEARzero[k]) {
00555 *nearzero= True;
00556 if (pivot_abs == 0.0) {
00557 if (qh IStracing >= 4) {
00558 fprintf (qh ferr, "qh_gausselim: 0 pivot at column %d. (%2.2g < %2.2g)\n", k, pivot_abs, qh DISTround);
00559 qh_printmatrix (qh ferr, "Matrix:", rows, numrow, numcol);
00560 }
00561 zzinc_(Zgauss0);
00562 qh_precision ("zero pivot for Gaussian elimination");
00563 goto LABELnextcol;
00564 }
00565 }
00566 pivotrow= rows[k] + k;
00567 pivot= *pivotrow++;
00568 for(i= k+1; i < numrow; i++) {
00569 ai= rows[i] + k;
00570 ak= pivotrow;
00571 n= (*ai++)/pivot;
00572 for(j= numcol - (k+1); j--; )
00573 *ai++ -= n * *ak++;
00574 }
00575 LABELnextcol:
00576 ;
00577 }
00578 wmin_(Wmindenom, pivot_abs);
00579 if (qh IStracing >= 5)
00580 qh_printmatrix (qh ferr, "qh_gausselem: result", rows, numrow, numcol);
00581 }
00582
00583
00584
00585
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595 realT qh_getangle(pointT *vect1, pointT *vect2) {
00596 realT angle= 0, randr;
00597 int k;
00598
00599 for(k= qh hull_dim; k--; )
00600 angle += *vect1++ * *vect2++;
00601 if (qh RANDOMdist) {
00602 randr= qh_RANDOMint;
00603 angle += (2.0 * randr / qh_RANDOMmax - 1.0) *
00604 qh RANDOMfactor;
00605 }
00606 trace4((qh ferr, "qh_getangle: %2.2g\n", angle));
00607 return(angle);
00608 }
00609
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620 pointT *qh_getcenter(setT *vertices) {
00621 int k;
00622 pointT *center, *coord;
00623 vertexT *vertex, **vertexp;
00624 int count= qh_setsize(vertices);
00625
00626 if (count < 2) {
00627 fprintf (qh ferr, "qhull internal error (qh_getcenter): not defined for %d points\n", count);
00628 qh_errexit (qh_ERRqhull, NULL, NULL);
00629 }
00630 center= (pointT *)qh_memalloc(qh normal_size);
00631 for (k=0; k < qh hull_dim; k++) {
00632 coord= center+k;
00633 *coord= 0.0;
00634 FOREACHvertex_(vertices)
00635 *coord += vertex->point[k];
00636 *coord /= count;
00637 }
00638 return(center);
00639 }
00640
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651 pointT *qh_getcentrum(facetT *facet) {
00652 realT dist;
00653 pointT *centrum, *point;
00654
00655 point= qh_getcenter(facet->vertices);
00656 zzinc_(Zcentrumtests);
00657 qh_distplane (point, facet, &dist);
00658 centrum= qh_projectpoint(point, facet, dist);
00659 qh_memfree(point, qh normal_size);
00660 trace4((qh ferr, "qh_getcentrum: for f%d, %d vertices dist= %2.2g\n",
00661 facet->id, qh_setsize(facet->vertices), dist));
00662 return centrum;
00663 }
00664
00665
00666
00667
00668
00669
00670
00671
00672
00673
00674
00675
00676
00677
00678
00679 realT qh_getdistance(facetT *facet, facetT *neighbor, realT *mindist, realT *maxdist) {
00680 vertexT *vertex, **vertexp;
00681 realT dist, maxd, mind;
00682
00683 FOREACHvertex_(facet->vertices)
00684 vertex->seen= False;
00685 FOREACHvertex_(neighbor->vertices)
00686 vertex->seen= True;
00687 mind= 0.0;
00688 maxd= 0.0;
00689 FOREACHvertex_(facet->vertices) {
00690 if (!vertex->seen) {
00691 zzinc_(Zbestdist);
00692 qh_distplane(vertex->point, neighbor, &dist);
00693 if (dist < mind)
00694 mind= dist;
00695 else if (dist > maxd)
00696 maxd= dist;
00697 }
00698 }
00699 *mindist= mind;
00700 *maxdist= maxd;
00701 mind= -mind;
00702 if (maxd > mind)
00703 return maxd;
00704 else
00705 return mind;
00706 }
00707
00708
00709
00710
00711
00712
00713
00714
00715
00716
00717
00718
00719 void qh_normalize (coordT *normal, int dim, boolT toporient) {
00720 qh_normalize2( normal, dim, toporient, NULL, NULL);
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 void qh_normalize2 (coordT *normal, int dim, boolT toporient,
00756 realT *minnorm, boolT *ismin) {
00757 int k;
00758 realT *colp, *maxp, norm= 0, temp, *norm1, *norm2, *norm3;
00759 boolT zerodiv;
00760
00761 norm1= normal+1;
00762 norm2= normal+2;
00763 norm3= normal+3;
00764 if (dim == 2)
00765 norm= sqrt((*normal)*(*normal) + (*norm1)*(*norm1));
00766 else if (dim == 3)
00767 norm= sqrt((*normal)*(*normal) + (*norm1)*(*norm1) + (*norm2)*(*norm2));
00768 else if (dim == 4) {
00769 norm= sqrt((*normal)*(*normal) + (*norm1)*(*norm1) + (*norm2)*(*norm2)
00770 + (*norm3)*(*norm3));
00771 }else if (dim > 4) {
00772 norm= (*normal)*(*normal) + (*norm1)*(*norm1) + (*norm2)*(*norm2)
00773 + (*norm3)*(*norm3);
00774 for (k= dim-4, colp= normal+4; k--; colp++)
00775 norm += (*colp) * (*colp);
00776 norm= sqrt(norm);
00777 }
00778 if (minnorm) {
00779 if (norm < *minnorm)
00780 *ismin= True;
00781 else
00782 *ismin= False;
00783 }
00784 wmin_(Wmindenom, norm);
00785 if (norm > qh MINdenom) {
00786 if (!toporient)
00787 norm= -norm;
00788 *normal /= norm;
00789 *norm1 /= norm;
00790 if (dim == 2)
00791 ;
00792 else if (dim == 3)
00793 *norm2 /= norm;
00794 else if (dim == 4) {
00795 *norm2 /= norm;
00796 *norm3 /= norm;
00797 }else if (dim >4) {
00798 *norm2 /= norm;
00799 *norm3 /= norm;
00800 for (k= dim-4, colp= normal+4; k--; )
00801 *colp++ /= norm;
00802 }
00803 }else if (norm == 0.0) {
00804 temp= sqrt (1.0/dim);
00805 for (k= dim, colp= normal; k--; )
00806 *colp++ = temp;
00807 }else {
00808 if (!toporient)
00809 norm= -norm;
00810 for (k= dim, colp= normal; k--; colp++) {
00811 temp= qh_divzero (*colp, norm, qh MINdenom_1, &zerodiv);
00812 if (!zerodiv)
00813 *colp= temp;
00814 else {
00815 maxp= qh_maxabsval(normal, dim);
00816 temp= ((*maxp * norm >= 0.0) ? 1.0 : -1.0);
00817 for (k= dim, colp= normal; k--; colp++)
00818 *colp= 0.0;
00819 *maxp= temp;
00820 zzinc_(Znearlysingular);
00821 trace0((qh ferr, "qh_normalize: norm=%2.2g too small during p%d\n",
00822 norm, qh furthest_id));
00823 return;
00824 }
00825 }
00826 }
00827 }
00828
00829
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839
00840
00841
00842
00843
00844 pointT *qh_projectpoint(pointT *point, facetT *facet, realT dist) {
00845 pointT *newpoint, *np, *normal;
00846 int normsize= qh normal_size,k;
00847 void **freelistp;
00848
00849 qh_memalloc_(normsize, freelistp, newpoint, pointT);
00850 np= newpoint;
00851 normal= facet->normal;
00852 for(k= qh hull_dim; k--; )
00853 *(np++)= *point++ - dist * *normal++;
00854 return(newpoint);
00855 }
00856
00857
00858
00859
00860
00861
00862
00863
00864
00865
00866
00867
00868
00869
00870
00871
00872
00873
00874
00875
00876
00877
00878
00879 void qh_setfacetplane(facetT *facet) {
00880 pointT *point;
00881 vertexT *vertex, **vertexp;
00882 int k,i, normsize= qh normal_size, oldtrace= 0;
00883 realT dist;
00884 void **freelistp;
00885 coordT *coord, *gmcoord;
00886 pointT *point0= SETfirstt_(facet->vertices, vertexT)->point;
00887 boolT nearzero= False;
00888
00889 zzinc_(Zsetplane);
00890 if (!facet->normal)
00891 qh_memalloc_(normsize, freelistp, facet->normal, coordT);
00892 if (facet == qh tracefacet) {
00893 oldtrace= qh IStracing;
00894 qh IStracing= 5;
00895 fprintf (qh ferr, "qh_setfacetplane: facet f%d created.\n", facet->id);
00896 fprintf (qh ferr, " Last point added to hull was p%d.", qh furthest_id);
00897 if (zzval_(Ztotmerge))
00898 fprintf(qh ferr, " Last merge was #%d.", zzval_(Ztotmerge));
00899 fprintf (qh ferr, "\n\nCurrent summary is:\n");
00900 qh_printsummary (qh ferr);
00901 }
00902 if (qh hull_dim <= 4) {
00903 i= 0;
00904 if (qh RANDOMdist) {
00905 gmcoord= qh gm_matrix;
00906 FOREACHvertex_(facet->vertices) {
00907 qh gm_row[i++]= gmcoord;
00908 coord= vertex->point;
00909 for (k= qh hull_dim; k--; )
00910 *(gmcoord++)= *coord++ * qh_randomfactor();
00911 }
00912 }else {
00913 FOREACHvertex_(facet->vertices)
00914 qh gm_row[i++]= vertex->point;
00915 }
00916 qh_sethyperplane_det(qh hull_dim, qh gm_row, point0, facet->toporient,
00917 facet->normal, &facet->offset, &nearzero);
00918 }
00919 if (qh hull_dim > 4 || nearzero) {
00920 i= 0;
00921 gmcoord= qh gm_matrix;
00922 FOREACHvertex_(facet->vertices) {
00923 if (vertex->point != point0) {
00924 qh gm_row[i++]= gmcoord;
00925 coord= vertex->point;
00926 point= point0;
00927 for(k= qh hull_dim; k--; )
00928 *(gmcoord++)= *coord++ - *point++;
00929 }
00930 }
00931 qh gm_row[i]= gmcoord;
00932 if (qh RANDOMdist) {
00933 gmcoord= qh gm_matrix;
00934 for (i= qh hull_dim-1; i--; ) {
00935 for (k= qh hull_dim; k--; )
00936 *(gmcoord++) *= qh_randomfactor();
00937 }
00938 }
00939 qh_sethyperplane_gauss(qh hull_dim, qh gm_row, point0, facet->toporient,
00940 facet->normal, &facet->offset, &nearzero);
00941 if (nearzero) {
00942 if (qh_orientoutside (facet)) {
00943 trace0((qh ferr, "qh_setfacetplane: flipped orientation after testing interior_point during p%d\n", qh furthest_id));
00944
00945
00946
00947
00948
00949
00950
00951
00952
00953
00954 }
00955 }
00956 }
00957 facet->upperdelaunay= False;
00958 if (qh DELAUNAY) {
00959 if (qh UPPERdelaunay) {
00960 if (facet->normal[qh hull_dim -1] >= qh ANGLEround * qh_ZEROdelaunay)
00961 facet->upperdelaunay= True;
00962 }else {
00963 if (facet->normal[qh hull_dim -1] > -qh ANGLEround * qh_ZEROdelaunay)
00964 facet->upperdelaunay= True;
00965 }
00966 }
00967 if (qh PRINTstatistics || qh IStracing || qh TRACElevel || qh JOGGLEmax < REALmax) {
00968 qh old_randomdist= qh RANDOMdist;
00969 qh RANDOMdist= False;
00970 FOREACHvertex_(facet->vertices) {
00971 if (vertex->point != point0) {
00972 boolT istrace= False;
00973 zinc_(Zdiststat);
00974 qh_distplane(vertex->point, facet, &dist);
00975 dist= fabs_(dist);
00976 zinc_(Znewvertex);
00977 wadd_(Wnewvertex, dist);
00978 if (dist > wwval_(Wnewvertexmax)) {
00979 wwval_(Wnewvertexmax)= dist;
00980 if (dist > qh max_outside) {
00981 qh max_outside= dist;
00982 if (dist > qh TRACEdist)
00983 istrace= True;
00984 }
00985 }else if (-dist > qh TRACEdist)
00986 istrace= True;
00987 if (istrace) {
00988 fprintf (qh ferr, "qh_setfacetplane: ====== vertex p%d (v%d) increases max_outside to %2.2g for new facet f%d last p%d\n",
00989 qh_pointid(vertex->point), vertex->id, dist, facet->id, qh furthest_id);
00990 qh_errprint ("DISTANT", facet, NULL, NULL, NULL);
00991 }
00992 }
00993 }
00994 qh RANDOMdist= qh old_randomdist;
00995 }
00996 if (qh IStracing >= 3) {
00997 fprintf (qh ferr, "qh_setfacetplane: f%d offset %2.2g normal: ",
00998 facet->id, facet->offset);
00999 for (k=0; k < qh hull_dim; k++)
01000 fprintf (qh ferr, "%2.2g ", facet->normal[k]);
01001 fprintf (qh ferr, "\n");
01002 }
01003 if (facet == qh tracefacet)
01004 qh IStracing= oldtrace;
01005 }
01006
01007
01008
01009
01010
01011
01012
01013
01014
01015
01016
01017
01018
01019
01020
01021
01022
01023
01024
01025
01026
01027
01028
01029
01030
01031
01032
01033
01034
01035
01036
01037
01038
01039
01040
01041
01042
01043
01044
01045
01046
01047
01048
01049
01050
01051
01052 void qh_sethyperplane_det (int dim, coordT **rows, coordT *point0,
01053 boolT toporient, coordT *normal, realT *offset, boolT *nearzero) {
01054 realT maxround, dist;
01055 int i;
01056 pointT *point;
01057
01058
01059 if (dim == 2) {
01060 normal[0]= dY(1,0);
01061 normal[1]= dX(0,1);
01062 qh_normalize2 (normal, dim, toporient, NULL, NULL);
01063 *offset= -(point0[0]*normal[0]+point0[1]*normal[1]);
01064 *nearzero= False;
01065 }else if (dim == 3) {
01066 normal[0]= det2_(dY(2,0), dZ(2,0),
01067 dY(1,0), dZ(1,0));
01068 normal[1]= det2_(dX(1,0), dZ(1,0),
01069 dX(2,0), dZ(2,0));
01070 normal[2]= det2_(dX(2,0), dY(2,0),
01071 dX(1,0), dY(1,0));
01072 qh_normalize2 (normal, dim, toporient, NULL, NULL);
01073 *offset= -(point0[0]*normal[0] + point0[1]*normal[1]
01074 + point0[2]*normal[2]);
01075 maxround= qh DISTround;
01076 for (i=dim; i--; ) {
01077 point= rows[i];
01078 if (point != point0) {
01079 dist= *offset + (point[0]*normal[0] + point[1]*normal[1]
01080 + point[2]*normal[2]);
01081 if (dist > maxround || dist < -maxround) {
01082 *nearzero= True;
01083 break;
01084 }
01085 }
01086 }
01087 }else if (dim == 4) {
01088 normal[0]= - det3_(dY(2,0), dZ(2,0), dW(2,0),
01089 dY(1,0), dZ(1,0), dW(1,0),
01090 dY(3,0), dZ(3,0), dW(3,0));
01091 normal[1]= det3_(dX(2,0), dZ(2,0), dW(2,0),
01092 dX(1,0), dZ(1,0), dW(1,0),
01093 dX(3,0), dZ(3,0), dW(3,0));
01094 normal[2]= - det3_(dX(2,0), dY(2,0), dW(2,0),
01095 dX(1,0), dY(1,0), dW(1,0),
01096 dX(3,0), dY(3,0), dW(3,0));
01097 normal[3]= det3_(dX(2,0), dY(2,0), dZ(2,0),
01098 dX(1,0), dY(1,0), dZ(1,0),
01099 dX(3,0), dY(3,0), dZ(3,0));
01100 qh_normalize2 (normal, dim, toporient, NULL, NULL);
01101 *offset= -(point0[0]*normal[0] + point0[1]*normal[1]
01102 + point0[2]*normal[2] + point0[3]*normal[3]);
01103 maxround= qh DISTround;
01104 for (i=dim; i--; ) {
01105 point= rows[i];
01106 if (point != point0) {
01107 dist= *offset + (point[0]*normal[0] + point[1]*normal[1]
01108 + point[2]*normal[2] + point[3]*normal[3]);
01109 if (dist > maxround || dist < -maxround) {
01110 *nearzero= True;
01111 break;
01112 }
01113 }
01114 }
01115 }
01116 if (*nearzero) {
01117 zzinc_(Zminnorm);
01118 trace0((qh ferr, "qh_sethyperplane_det: degenerate norm during p%d.\n", qh furthest_id));
01119 zzinc_(Znearlysingular);
01120 }
01121 }
01122
01123
01124
01125
01126
01127
01128
01129
01130
01131
01132
01133
01134
01135
01136
01137
01138
01139
01140
01141
01142
01143
01144
01145
01146
01147
01148
01149
01150 void qh_sethyperplane_gauss (int dim, coordT **rows, pointT *point0,
01151 boolT toporient, coordT *normal, coordT *offset, boolT *nearzero) {
01152 coordT *pointcoord, *normalcoef;
01153 int k;
01154 boolT sign= toporient, nearzero2= False;
01155
01156 qh_gausselim(rows, dim-1, dim, &sign, nearzero);
01157 for(k= dim-1; k--; ) {
01158 if ((rows[k])[k] < 0)
01159 sign ^= 1;
01160 }
01161 if (*nearzero) {
01162 zzinc_(Znearlysingular);
01163 trace0((qh ferr, "qh_sethyperplane_gauss: nearly singular or axis parallel hyperplane during p%d.\n", qh furthest_id));
01164 qh_backnormal(rows, dim-1, dim, sign, normal, &nearzero2);
01165 }else {
01166 qh_backnormal(rows, dim-1, dim, sign, normal, &nearzero2);
01167 if (nearzero2) {
01168 zzinc_(Znearlysingular);
01169 trace0((qh ferr, "qh_sethyperplane_gauss: singular or axis parallel hyperplane at normalization during p%d.\n", qh furthest_id));
01170 }
01171 }
01172 if (nearzero2)
01173 *nearzero= True;
01174 qh_normalize2(normal, dim, True, NULL, NULL);
01175 pointcoord= point0;
01176 normalcoef= normal;
01177 *offset= -(*pointcoord++ * *normalcoef++);
01178 for(k= dim-1; k--; )
01179 *offset -= *pointcoord++ * *normalcoef++;
01180 }
01181
01182