00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015 #include "qhull_a.h"
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 coordT *qh_copypoints (coordT *points, int numpoints, int dimension) {
00026 int size;
00027 coordT *newpoints;
00028
00029 size= numpoints * dimension * sizeof(coordT);
00030 if (!(newpoints=(coordT*)malloc(size))) {
00031 fprintf(qh ferr, "qhull error: insufficient memory to copy %d points\n",
00032 numpoints);
00033 qh_errexit(qh_ERRmem, NULL, NULL);
00034 }
00035 memcpy ((char *)newpoints, (char *)points, size);
00036 return newpoints;
00037 }
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050 void qh_crossproduct (int dim, realT vecA[3], realT vecB[3], realT vecC[3]){
00051
00052 if (dim == 3) {
00053 vecC[0]= det2_(vecA[1], vecA[2],
00054 vecB[1], vecB[2]);
00055 vecC[1]= - det2_(vecA[0], vecA[2],
00056 vecB[0], vecB[2]);
00057 vecC[2]= det2_(vecA[0], vecA[1],
00058 vecB[0], vecB[1]);
00059 }
00060 }
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078 realT qh_determinant (realT **rows, int dim, boolT *nearzero) {
00079 realT det=0;
00080 int i;
00081 boolT sign= False;
00082
00083 *nearzero= False;
00084 if (dim < 2) {
00085 fprintf (qh ferr, "qhull internal error (qh_determinate): only implemented for dimension >= 2\n");
00086 qh_errexit (qh_ERRqhull, NULL, NULL);
00087 }else if (dim == 2) {
00088 det= det2_(rows[0][0], rows[0][1],
00089 rows[1][0], rows[1][1]);
00090 if (fabs_(det) < qh NEARzero[1])
00091 *nearzero= True;
00092 }else if (dim == 3) {
00093 det= det3_(rows[0][0], rows[0][1], rows[0][2],
00094 rows[1][0], rows[1][1], rows[1][2],
00095 rows[2][0], rows[2][1], rows[2][2]);
00096 if (fabs_(det) < qh NEARzero[2])
00097 *nearzero= True;
00098 }else {
00099 qh_gausselim(rows, dim, dim, &sign, nearzero);
00100 det= 1.0;
00101 for (i= dim; i--; )
00102 det *= (rows[i])[i];
00103 if (sign)
00104 det= -det;
00105 }
00106 return det;
00107 }
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125 realT qh_detjoggle (pointT *points, int numpoints, int dimension) {
00126 realT abscoord, distround, joggle, maxcoord, mincoord;
00127 pointT *point, *pointtemp;
00128 realT maxabs= -REALmax;
00129 realT sumabs= 0;
00130 realT maxwidth= 0;
00131 int k;
00132
00133 for (k= 0; k < dimension; k++) {
00134 if (qh SCALElast && k == dimension-1)
00135 abscoord= maxwidth;
00136 else if (qh DELAUNAY && k == dimension-1)
00137 abscoord= 2 * maxabs * maxabs;
00138 else {
00139 maxcoord= -REALmax;
00140 mincoord= REALmax;
00141 FORALLpoint_(points, numpoints) {
00142 maximize_(maxcoord, point[k]);
00143 minimize_(mincoord, point[k]);
00144 }
00145 maximize_(maxwidth, maxcoord-mincoord);
00146 abscoord= fmax_(maxcoord, -mincoord);
00147 }
00148 sumabs += abscoord;
00149 maximize_(maxabs, abscoord);
00150 }
00151 distround= qh_distround (qh hull_dim, maxabs, sumabs);
00152 joggle= distround * qh_JOGGLEdefault;
00153 maximize_(joggle, REALepsilon * qh_JOGGLEdefault);
00154 trace2((qh ferr, "qh_detjoggle: joggle=%2.2g maxwidth=%2.2g\n", joggle, maxwidth));
00155 return joggle;
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 void qh_detroundoff (void) {
00189
00190 qh_option ("_max-width", NULL, &qh MAXwidth);
00191 if (!qh SETroundoff) {
00192 qh DISTround= qh_distround (qh hull_dim, qh MAXabs_coord, qh MAXsumcoord);
00193 if (qh RANDOMdist)
00194 qh DISTround += qh RANDOMfactor * qh MAXabs_coord;
00195 qh_option ("Error-roundoff", NULL, &qh DISTround);
00196 }
00197 qh MINdenom= qh MINdenom_1 * qh MAXabs_coord;
00198 qh MINdenom_1_2= sqrt (qh MINdenom_1 * qh hull_dim) ;
00199 qh MINdenom_2= qh MINdenom_1_2 * qh MAXabs_coord;
00200
00201 qh ANGLEround= 1.01 * qh hull_dim * REALepsilon;
00202 if (qh RANDOMdist)
00203 qh ANGLEround += qh RANDOMfactor;
00204 if (qh premerge_cos < REALmax/2) {
00205 qh premerge_cos -= qh ANGLEround;
00206 if (qh RANDOMdist)
00207 qh_option ("Angle-premerge-with-random", NULL, &qh premerge_cos);
00208 }
00209 if (qh postmerge_cos < REALmax/2) {
00210 qh postmerge_cos -= qh ANGLEround;
00211 if (qh RANDOMdist)
00212 qh_option ("Angle-postmerge-with-random", NULL, &qh postmerge_cos);
00213 }
00214 qh premerge_centrum += 2 * qh DISTround;
00215 qh postmerge_centrum += 2 * qh DISTround;
00216 if (qh RANDOMdist && (qh MERGEexact || qh PREmerge))
00217 qh_option ("Centrum-premerge-with-random", NULL, &qh premerge_centrum);
00218 if (qh RANDOMdist && qh POSTmerge)
00219 qh_option ("Centrum-postmerge-with-random", NULL, &qh postmerge_centrum);
00220 {
00221 realT maxangle= 1.0, maxrho;
00222
00223 minimize_(maxangle, qh premerge_cos);
00224 minimize_(maxangle, qh postmerge_cos);
00225
00226 qh ONEmerge= sqrt (qh hull_dim) * qh MAXwidth *
00227 sqrt (1.0 - maxangle * maxangle) + qh DISTround;
00228 maxrho= qh hull_dim * qh premerge_centrum + qh DISTround;
00229 maximize_(qh ONEmerge, maxrho);
00230 maxrho= qh hull_dim * qh postmerge_centrum + qh DISTround;
00231 maximize_(qh ONEmerge, maxrho);
00232 if (qh MERGING)
00233 qh_option ("_one-merge", NULL, &qh ONEmerge);
00234 }
00235 qh NEARinside= qh ONEmerge * qh_RATIOnearinside;
00236 if (qh JOGGLEmax < REALmax/2 && (qh KEEPcoplanar || qh KEEPinside)) {
00237 realT maxdist;
00238
00239 qh KEEPnearinside= True;
00240 maxdist= sqrt (qh hull_dim) * qh JOGGLEmax + qh DISTround;
00241 maxdist= 2*maxdist;
00242 maximize_(qh NEARinside, maxdist);
00243 }
00244 if (qh KEEPnearinside)
00245 qh_option ("_near-inside", NULL, &qh NEARinside);
00246 if (qh JOGGLEmax < qh DISTround) {
00247 fprintf (qh ferr, "qhull error: the joggle for 'QJn', %.2g, is below roundoff for distance computations, %.2g\n",
00248 qh JOGGLEmax, qh DISTround);
00249 qh_errexit (qh_ERRinput, NULL, NULL);
00250 }
00251 if (qh MINvisible > REALmax/2) {
00252 if (!qh MERGING)
00253 qh MINvisible= qh DISTround;
00254 else if (qh hull_dim <= 3)
00255 qh MINvisible= qh premerge_centrum;
00256 else
00257 qh MINvisible= qh_COPLANARratio * qh premerge_centrum;
00258 if (qh APPROXhull && qh MINvisible > qh MINoutside)
00259 qh MINvisible= qh MINoutside;
00260 qh_option ("Visible-distance", NULL, &qh MINvisible);
00261 }
00262 if (qh MAXcoplanar > REALmax/2) {
00263 qh MAXcoplanar= qh MINvisible;
00264 qh_option ("U-coplanar-distance", NULL, &qh MAXcoplanar);
00265 }
00266 if (!qh APPROXhull) {
00267 qh MINoutside= 2 * qh MINvisible;
00268 if (qh premerge_cos < REALmax/2)
00269 maximize_(qh MINoutside, (1- qh premerge_cos) * qh MAXabs_coord);
00270 qh_option ("Width-outside", NULL, &qh MINoutside);
00271 }
00272 qh WIDEfacet= qh MINoutside;
00273 maximize_(qh WIDEfacet, qh_WIDEcoplanar * qh MAXcoplanar);
00274 maximize_(qh WIDEfacet, qh_WIDEcoplanar * qh MINvisible);
00275 qh_option ("_wide-facet", NULL, &qh WIDEfacet);
00276 if (qh MINvisible > qh MINoutside + 3 * REALepsilon
00277 && !qh BESToutside && !qh FORCEoutput)
00278 fprintf (qh ferr, "qhull input warning: minimum visibility V%.2g is greater than \nminimum outside W%.2g. Flipped facets are likely.\n",
00279 qh MINvisible, qh MINoutside);
00280 qh max_vertex= qh DISTround;
00281 qh min_vertex= -qh DISTround;
00282
00283 }
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301 realT qh_detsimplex(pointT *apex, setT *points, int dim, boolT *nearzero) {
00302 pointT *coorda, *coordp, *gmcoord, *point, **pointp;
00303 coordT **rows;
00304 int k, i=0;
00305 realT det;
00306
00307 zinc_(Zdetsimplex);
00308 gmcoord= qh gm_matrix;
00309 rows= qh gm_row;
00310 FOREACHpoint_(points) {
00311 if (i == dim)
00312 break;
00313 rows[i++]= gmcoord;
00314 coordp= point;
00315 coorda= apex;
00316 for (k= dim; k--; )
00317 *(gmcoord++)= *coordp++ - *coorda++;
00318 }
00319 if (i < dim) {
00320 fprintf (qh ferr, "qhull internal error (qh_detsimplex): #points %d < dimension %d\n",
00321 i, dim);
00322 qh_errexit (qh_ERRqhull, NULL, NULL);
00323 }
00324 det= qh_determinant (rows, dim, nearzero);
00325 trace2((qh ferr, "qh_detsimplex: det=%2.2g for point p%d, dim %d, nearzero? %d\n",
00326 det, qh_pointid(apex), dim, *nearzero));
00327 return det;
00328 }
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345 realT qh_distnorm (int dim, pointT *point, pointT *normal, realT *offsetp) {
00346 coordT *normalp= normal, *coordp= point;
00347 realT dist;
00348 int k;
00349
00350 dist= *offsetp;
00351 for (k= dim; k--; )
00352 dist += *(coordp++) * *(normalp++);
00353 return dist;
00354 }
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374 realT qh_distround (int dimension, realT maxabs, realT maxsumabs) {
00375 realT maxdistsum, maxround;
00376
00377 maxdistsum= sqrt (dimension) * maxabs;
00378 minimize_( maxdistsum, maxsumabs);
00379 maxround= REALepsilon * (dimension * maxdistsum * 1.01 + maxabs);
00380
00381 trace4((qh ferr, "qh_distround: %2.2g maxabs %2.2g maxsumabs %2.2g maxdistsum %2.2g\n",
00382 maxround, maxabs, maxsumabs, maxdistsum));
00383 return maxround;
00384 }
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407 realT qh_divzero (realT numer, realT denom, realT mindenom1, boolT *zerodiv) {
00408 realT temp, numerx, denomx;
00409
00410
00411 if (numer < mindenom1 && numer > -mindenom1) {
00412 numerx= fabs_(numer);
00413 denomx= fabs_(denom);
00414 if (numerx < denomx) {
00415 *zerodiv= False;
00416 return numer/denom;
00417 }else {
00418 *zerodiv= True;
00419 return 0.0;
00420 }
00421 }
00422 temp= denom/numer;
00423 if (temp > mindenom1 || temp < -mindenom1) {
00424 *zerodiv= False;
00425 return numer/denom;
00426 }else {
00427 *zerodiv= True;
00428 return 0.0;
00429 }
00430 }
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445
00446
00447
00448
00449
00450
00451
00452
00453
00454 realT qh_facetarea (facetT *facet) {
00455 vertexT *apex;
00456 pointT *centrum;
00457 realT area= 0.0;
00458 ridgeT *ridge, **ridgep;
00459
00460 if (facet->simplicial) {
00461 apex= SETfirstt_(facet->vertices, vertexT);
00462 area= qh_facetarea_simplex (qh hull_dim, apex->point, facet->vertices,
00463 apex, facet->toporient, facet->normal, &facet->offset);
00464 }else {
00465 if (qh CENTERtype == qh_AScentrum)
00466 centrum= facet->center;
00467 else
00468 centrum= qh_getcentrum (facet);
00469 FOREACHridge_(facet->ridges)
00470 area += qh_facetarea_simplex (qh hull_dim, centrum, ridge->vertices,
00471 NULL, (ridge->top == facet), facet->normal, &facet->offset);
00472 if (qh CENTERtype != qh_AScentrum)
00473 qh_memfree (centrum, qh normal_size);
00474 }
00475 if (facet->upperdelaunay && qh DELAUNAY)
00476 area= -area;
00477 trace4((qh ferr, "qh_facetarea: f%d area %2.2g\n", facet->id, area));
00478 return area;
00479 }
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513 realT qh_facetarea_simplex (int dim, coordT *apex, setT *vertices,
00514 vertexT *notvertex, boolT toporient, coordT *normal, realT *offset) {
00515 pointT *coorda, *coordp, *gmcoord;
00516 coordT **rows, *normalp;
00517 int k, i=0;
00518 realT area, dist;
00519 vertexT *vertex, **vertexp;
00520 boolT nearzero;
00521
00522 gmcoord= qh gm_matrix;
00523 rows= qh gm_row;
00524 FOREACHvertex_(vertices) {
00525 if (vertex == notvertex)
00526 continue;
00527 rows[i++]= gmcoord;
00528 coorda= apex;
00529 coordp= vertex->point;
00530 normalp= normal;
00531 if (notvertex) {
00532 for (k= dim; k--; )
00533 *(gmcoord++)= *coordp++ - *coorda++;
00534 }else {
00535 dist= *offset;
00536 for (k= dim; k--; )
00537 dist += *coordp++ * *normalp++;
00538 if (dist < -qh WIDEfacet) {
00539 zinc_(Znoarea);
00540 return 0.0;
00541 }
00542 coordp= vertex->point;
00543 normalp= normal;
00544 for (k= dim; k--; )
00545 *(gmcoord++)= (*coordp++ - dist * *normalp++) - *coorda++;
00546 }
00547 }
00548 if (i != dim-1) {
00549 fprintf (qh ferr, "qhull internal error (qh_facetarea_simplex): #points %d != dim %d -1\n",
00550 i, dim);
00551 qh_errexit (qh_ERRqhull, NULL, NULL);
00552 }
00553 rows[i]= gmcoord;
00554 if (qh DELAUNAY) {
00555 for (i= 0; i < dim-1; i++)
00556 rows[i][dim-1]= 0.0;
00557 for (k= dim; k--; )
00558 *(gmcoord++)= 0.0;
00559 rows[dim-1][dim-1]= -1.0;
00560 }else {
00561 normalp= normal;
00562 for (k= dim; k--; )
00563 *(gmcoord++)= *normalp++;
00564 }
00565 zinc_(Zdetsimplex);
00566 area= qh_determinant (rows, dim, &nearzero);
00567 if (toporient)
00568 area= -area;
00569 area *= qh AREAfactor;
00570 trace4((qh ferr, "qh_facetarea_simplex: area=%2.2g for point p%d, toporient %d, nearzero? %d\n",
00571 area, qh_pointid(apex), toporient, nearzero));
00572 return area;
00573 }
00574
00575
00576
00577
00578
00579
00580
00581
00582
00583
00584
00585
00586
00587 pointT *qh_facetcenter (setT *vertices) {
00588 setT *points= qh_settemp (qh_setsize (vertices));
00589 vertexT *vertex, **vertexp;
00590 pointT *center;
00591
00592 FOREACHvertex_(vertices)
00593 qh_setappend (&points, vertex->point);
00594 center= qh_voronoi_center (qh hull_dim-1, points);
00595 qh_settempfree (&points);
00596 return center;
00597 }
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620
00621
00622 boolT qh_findbestsharp (pointT *point, facetT **bestfacet,
00623 realT *bestdist, int *numpart) {
00624 facetT *facet;
00625 realT dist;
00626 boolT issharp = False;
00627 int *quadrant, k;
00628
00629 quadrant= (int*)qh_memalloc (qh hull_dim * sizeof(int));
00630 FORALLfacet_(qh newfacet_list) {
00631 if (facet == qh newfacet_list) {
00632 for (k= qh hull_dim; k--; )
00633 quadrant[ k]= (facet->normal[ k] > 0);
00634 }else if (!issharp) {
00635 for (k= qh hull_dim; k--; ) {
00636 if (quadrant[ k] != (facet->normal[ k] > 0)) {
00637 issharp= True;
00638 break;
00639 }
00640 }
00641 }
00642 if (facet->visitid != qh visit_id) {
00643 qh_distplane (point, facet, &dist);
00644 (*numpart)++;
00645 if (dist > *bestdist) {
00646 if (!facet->upperdelaunay || dist > qh MINoutside) {
00647 *bestdist= dist;
00648 *bestfacet= facet;
00649 }
00650 }
00651 }
00652 }
00653 qh_memfree( quadrant, qh hull_dim * sizeof(int));
00654 return issharp;
00655 }
00656
00657
00658
00659
00660
00661
00662
00663
00664
00665
00666
00667
00668
00669
00670
00671
00672
00673
00674
00675
00676
00677
00678
00679
00680
00681
00682
00683
00684 facetT *qh_findgooddist (pointT *point, facetT *facetA, realT *distp,
00685 facetT **facetlist) {
00686 realT bestdist= -REALmax, dist;
00687 facetT *neighbor, **neighborp, *bestfacet=NULL, *facet;
00688 boolT goodseen= False;
00689
00690 if (facetA->good) {
00691 zinc_(Zcheckpart);
00692 qh_distplane (point, facetA, &bestdist);
00693 bestfacet= facetA;
00694 goodseen= True;
00695 }
00696 qh_removefacet (facetA);
00697 qh_appendfacet (facetA);
00698 *facetlist= facetA;
00699 facetA->visitid= ++qh visit_id;
00700 FORALLfacet_(*facetlist) {
00701 FOREACHneighbor_(facet) {
00702 if (neighbor->visitid == qh visit_id)
00703 continue;
00704 neighbor->visitid= qh visit_id;
00705 if (goodseen && !neighbor->good)
00706 continue;
00707 zinc_(Zcheckpart);
00708 qh_distplane (point, neighbor, &dist);
00709 if (dist > 0) {
00710 qh_removefacet (neighbor);
00711 qh_appendfacet (neighbor);
00712 if (neighbor->good) {
00713 goodseen= True;
00714 if (dist > bestdist) {
00715 bestdist= dist;
00716 bestfacet= neighbor;
00717 }
00718 }
00719 }
00720 }
00721 }
00722 if (bestfacet) {
00723 *distp= bestdist;
00724 trace2((qh ferr, "qh_findgooddist: p%d is %2.2g above good facet f%d\n",
00725 qh_pointid(point), bestdist, bestfacet->id));
00726 return bestfacet;
00727 }
00728 trace4((qh ferr, "qh_findgooddist: no good facet for p%d above f%d\n",
00729 qh_pointid(point), facetA->id));
00730 return NULL;
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_getarea (facetT *facetlist) {
00756 realT area;
00757 realT dist;
00758 facetT *facet;
00759
00760 if (qh REPORTfreq)
00761 fprintf (qh ferr, "computing area of each facet and volume of the convex hull\n");
00762 else
00763 trace1((qh ferr, "qh_getarea: computing volume and area for each facet\n"));
00764 qh totarea= qh totvol= 0.0;
00765 FORALLfacet_(facetlist) {
00766 if (!facet->normal)
00767 continue;
00768 if (facet->upperdelaunay && qh ATinfinity)
00769 continue;
00770 facet->f.area= area= qh_facetarea (facet);
00771 facet->isarea= True;
00772 if (qh DELAUNAY) {
00773 if (facet->upperdelaunay == qh UPPERdelaunay)
00774 qh totarea += area;
00775 }else {
00776 qh totarea += area;
00777 qh_distplane (qh interior_point, facet, &dist);
00778 qh totvol += -dist * area/ qh hull_dim;
00779 }
00780 if (qh PRINTstatistics) {
00781 wadd_(Wareatot, area);
00782 wmax_(Wareamax, area);
00783 wmin_(Wareamin, area);
00784 }
00785 }
00786 }
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796
00797
00798
00799
00800
00801
00802
00803
00804
00805
00806
00807
00808
00809
00810 boolT qh_gram_schmidt(int dim, realT **row) {
00811 realT *rowi, *rowj, norm;
00812 int i, j, k;
00813
00814 for(i=0; i < dim; i++) {
00815 rowi= row[i];
00816 for (norm= 0.0, k= dim; k--; rowi++)
00817 norm += *rowi * *rowi;
00818 norm= sqrt(norm);
00819 wmin_(Wmindenom, norm);
00820 if (norm == 0.0)
00821 return False;
00822 for(k= dim; k--; )
00823 *(--rowi) /= norm;
00824 for(j= i+1; j < dim; j++) {
00825 rowj= row[j];
00826 for(norm= 0.0, k=dim; k--; )
00827 norm += *rowi++ * *rowj++;
00828 for(k=dim; k--; )
00829 *(--rowj) -= *(--rowi) * norm;
00830 }
00831 }
00832 return True;
00833 }
00834
00835
00836
00837
00838
00839
00840
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854 boolT qh_inthresholds (coordT *normal, realT *angle) {
00855 boolT within= True;
00856 int k;
00857 realT threshold;
00858
00859 if (angle)
00860 *angle= 0.0;
00861 for(k= 0; k < qh hull_dim; k++) {
00862 threshold= qh lower_threshold[k];
00863 if (threshold > -REALmax/2) {
00864 if (normal[k] < threshold)
00865 within= False;
00866 if (angle) {
00867 threshold -= normal[k];
00868 *angle += fabs_(threshold);
00869 }
00870 }
00871 if (qh upper_threshold[k] < REALmax/2) {
00872 threshold= qh upper_threshold[k];
00873 if (normal[k] > threshold)
00874 within= False;
00875 if (angle) {
00876 threshold -= normal[k];
00877 *angle += fabs_(threshold);
00878 }
00879 }
00880 }
00881 return within;
00882 }
00883
00884
00885
00886
00887
00888
00889
00890
00891
00892
00893
00894
00895
00896
00897
00898
00899
00900
00901
00902
00903
00904
00905
00906
00907
00908
00909
00910
00911
00912
00913
00914
00915
00916
00917 void qh_joggleinput (void) {
00918 int size, i, seed;
00919 coordT *coordp, *inputp;
00920 realT randr, randa, randb;
00921
00922 if (!qh input_points) {
00923 qh input_points= qh first_point;
00924 qh input_malloc= qh POINTSmalloc;
00925 size= qh num_points * qh hull_dim * sizeof(coordT);
00926 if (!(qh first_point=(coordT*)malloc(size))) {
00927 fprintf(qh ferr, "qhull error: insufficient memory to joggle %d points\n",
00928 qh num_points);
00929 qh_errexit(qh_ERRmem, NULL, NULL);
00930 }
00931 qh POINTSmalloc= True;
00932 if (qh JOGGLEmax == 0.0) {
00933 qh JOGGLEmax= qh_detjoggle (qh input_points, qh num_points, qh hull_dim);
00934 qh_option ("QJoggle", NULL, &qh JOGGLEmax);
00935 }
00936 }else {
00937 if (!qh RERUN && qh build_cnt > qh_JOGGLEretry) {
00938 if (((qh build_cnt-qh_JOGGLEretry-1) % qh_JOGGLEagain) == 0) {
00939 realT maxjoggle= qh MAXwidth * qh_JOGGLEmaxincrease;
00940 if (qh JOGGLEmax < maxjoggle) {
00941 qh JOGGLEmax *= qh_JOGGLEincrease;
00942 minimize_(qh JOGGLEmax, maxjoggle);
00943 }
00944 }
00945 }
00946 qh_option ("QJoggle", NULL, &qh JOGGLEmax);
00947 }
00948 if (qh build_cnt > 1 && qh JOGGLEmax > fmax_(qh MAXwidth/4, 0.1)) {
00949 fprintf (qh ferr, "qhull error: the current joggle for 'QJn', %.2g, is too large for the width\nof the input. If possible, recompile Qhull with higher-precision reals.\n",
00950 qh JOGGLEmax);
00951 qh_errexit (qh_ERRqhull, NULL, NULL);
00952 }
00953 seed= qh_RANDOMint;
00954 qh_option ("_joggle-seed", &seed, NULL);
00955 trace0((qh ferr, "qh_joggleinput: joggle input by %2.2g with seed %d\n",
00956 qh JOGGLEmax, seed));
00957 inputp= qh input_points;
00958 coordp= qh first_point;
00959 randa= 2.0 * qh JOGGLEmax/qh_RANDOMmax;
00960 randb= -qh JOGGLEmax;
00961 size= qh num_points * qh hull_dim;
00962 for (i= size; i--; ) {
00963 randr= qh_RANDOMint;
00964 *(coordp++)= *(inputp++) + (randr * randa + randb);
00965 }
00966 if (qh DELAUNAY) {
00967 qh last_low= qh last_high= qh last_newhigh= REALmax;
00968 qh_setdelaunay (qh hull_dim, qh num_points, qh first_point);
00969 }
00970 }
00971
00972
00973
00974
00975
00976
00977
00978
00979 realT *qh_maxabsval (realT *normal, int dim) {
00980 realT maxval= -REALmax;
00981 realT *maxp= NULL, *colp, absval;
00982 int k;
00983
00984 for (k= dim, colp= normal; k--; colp++) {
00985 absval= fabs_(*colp);
00986 if (absval > maxval) {
00987 maxval= absval;
00988 maxp= colp;
00989 }
00990 }
00991 return maxp;
00992 }
00993
00994
00995
00996
00997
00998
00999
01000
01001
01002
01003
01004
01005
01006
01007
01008
01009
01010
01011
01012
01013
01014
01015
01016
01017
01018
01019
01020
01021 setT *qh_maxmin(pointT *points, int numpoints, int dimension) {
01022 int k;
01023 realT maxcoord, temp;
01024 pointT *minimum, *maximum, *point, *pointtemp;
01025 setT *set;
01026
01027 qh max_outside= 0.0;
01028 qh MAXabs_coord= 0.0;
01029 qh MAXwidth= -REALmax;
01030 qh MAXsumcoord= 0.0;
01031 qh min_vertex= 0.0;
01032 qh WAScoplanar= False;
01033 if (qh ZEROcentrum)
01034 qh ZEROall_ok= True;
01035 if (REALmin < REALepsilon && REALmin < REALmax && REALmin > -REALmax
01036 && REALmax > 0.0 && -REALmax < 0.0)
01037 ;
01038 else {
01039 fprintf (qh ferr, "qhull error: floating point constants in user.h are wrong\n\
01040 REALepsilon %g REALmin %g REALmax %g -REALmax %g\n",
01041 REALepsilon, REALmin, REALmax, -REALmax);
01042 qh_errexit (qh_ERRinput, NULL, NULL);
01043 }
01044 set= qh_settemp(2*dimension);
01045 for(k= 0; k < dimension; k++) {
01046 if (points == qh GOODpointp)
01047 minimum= maximum= points + dimension;
01048 else
01049 minimum= maximum= points;
01050 FORALLpoint_(points, numpoints) {
01051 if (point == qh GOODpointp)
01052 continue;
01053 if (maximum[k] < point[k])
01054 maximum= point;
01055 else if (minimum[k] > point[k])
01056 minimum= point;
01057 }
01058 if (k == dimension-1) {
01059 qh MINlastcoord= minimum[k];
01060 qh MAXlastcoord= maximum[k];
01061 }
01062 if (qh SCALElast && k == dimension-1)
01063 maxcoord= qh MAXwidth;
01064 else {
01065 maxcoord= fmax_(maximum[k], -minimum[k]);
01066 if (qh GOODpointp) {
01067 temp= fmax_(qh GOODpointp[k], -qh GOODpointp[k]);
01068 maximize_(maxcoord, temp);
01069 }
01070 temp= maximum[k] - minimum[k];
01071 maximize_(qh MAXwidth, temp);
01072 }
01073 maximize_(qh MAXabs_coord, maxcoord);
01074 qh MAXsumcoord += maxcoord;
01075 qh_setappend (&set, maximum);
01076 qh_setappend (&set, minimum);
01077
01078
01079
01080 qh NEARzero[k]= 80 * qh MAXsumcoord * REALepsilon;
01081 }
01082 if (qh IStracing >=1)
01083 qh_printpoints (qh ferr, "qh_maxmin: found the max and min points (by dim):", set);
01084 return(set);
01085 }
01086
01087
01088
01089
01090
01091
01092
01093
01094
01095
01096
01097
01098
01099
01100
01101
01102
01103
01104
01105
01106
01107 realT qh_maxouter (void) {
01108 realT dist;
01109
01110 dist= fmax_(qh max_outside, qh DISTround);
01111 dist += qh DISTround;
01112 trace4((qh ferr, "qh_maxouter: max distance from facet to outer plane is %2.2g max_outside is %2.2g\n", dist, qh max_outside));
01113 return dist;
01114 }
01115
01116
01117
01118
01119
01120
01121
01122
01123
01124
01125
01126
01127
01128
01129
01130
01131
01132
01133
01134
01135
01136
01137
01138
01139 void qh_maxsimplex (int dim, setT *maxpoints, pointT *points, int numpoints, setT **simplex) {
01140 pointT *point, **pointp, *pointtemp, *maxpoint, *minx=NULL, *maxx=NULL;
01141 boolT nearzero, maxnearzero= False;
01142 int k, sizinit;
01143 realT maxdet= -REALmax, det, mincoord= REALmax, maxcoord= -REALmax;
01144
01145 sizinit= qh_setsize (*simplex);
01146 if (sizinit < 2) {
01147 if (qh_setsize (maxpoints) >= 2) {
01148 FOREACHpoint_(maxpoints) {
01149
01150 if (maxcoord < point[0]) {
01151 maxcoord= point[0];
01152 maxx= point;
01153 }
01154 if (mincoord > point[0]) {
01155 mincoord= point[0];
01156 minx= point;
01157 }
01158 }
01159 }else {
01160 FORALLpoint_(points, numpoints) {
01161 if (point == qh GOODpointp)
01162 continue;
01163 if (maxcoord < point[0]) {
01164 maxcoord= point[0];
01165 maxx= point;
01166 }
01167 if (mincoord > point[0]) {
01168 mincoord= point[0];
01169 minx= point;
01170 }
01171 }
01172 }
01173 qh_setunique (simplex, minx);
01174 if (qh_setsize (*simplex) < 2)
01175 qh_setunique (simplex, maxx);
01176 sizinit= qh_setsize (*simplex);
01177 if (sizinit < 2) {
01178 qh_precision ("input has same x coordinate");
01179 if (zzval_(Zsetplane) > qh hull_dim+1) {
01180 fprintf (qh ferr, "qhull precision error (qh_maxsimplex for voronoi_center):\n%d points with the same x coordinate.\n",
01181 qh_setsize(maxpoints)+numpoints);
01182 qh_errexit (qh_ERRprec, NULL, NULL);
01183 }else {
01184 fprintf (qh ferr, "qhull input error: input is less than %d-dimensional since it has the same x coordinate\n", qh hull_dim);
01185 qh_errexit (qh_ERRinput, NULL, NULL);
01186 }
01187 }
01188 }
01189 for(k= sizinit; k < dim+1; k++) {
01190 maxpoint= NULL;
01191 maxdet= -REALmax;
01192 FOREACHpoint_(maxpoints) {
01193 if (!qh_setin (*simplex, point)) {
01194 det= qh_detsimplex(point, *simplex, k, &nearzero);
01195 if ((det= fabs_(det)) > maxdet) {
01196 maxdet= det;
01197 maxpoint= point;
01198 maxnearzero= nearzero;
01199 }
01200 }
01201 }
01202 if (!maxpoint || maxnearzero) {
01203 zinc_(Zsearchpoints);
01204 if (!maxpoint) {
01205 trace0((qh ferr, "qh_maxsimplex: searching all points for %d-th initial vertex.\n", k+1));
01206 }else {
01207 trace0((qh ferr, "qh_maxsimplex: searching all points for %d-th initial vertex, better than p%d det %2.2g\n",
01208 k+1, qh_pointid(maxpoint), maxdet));
01209 }
01210 FORALLpoint_(points, numpoints) {
01211 if (point == qh GOODpointp)
01212 continue;
01213 if (!qh_setin (*simplex, point)) {
01214 det= qh_detsimplex(point, *simplex, k, &nearzero);
01215 if ((det= fabs_(det)) > maxdet) {
01216 maxdet= det;
01217 maxpoint= point;
01218 maxnearzero= nearzero;
01219 }
01220 }
01221 }
01222 }
01223 if (!maxpoint) {
01224 fprintf (qh ferr, "qhull internal error (qh_maxsimplex): not enough points available\n");
01225 qh_errexit (qh_ERRqhull, NULL, NULL);
01226 }
01227 qh_setappend(simplex, maxpoint);
01228 trace1((qh ferr, "qh_maxsimplex: selected point p%d for %d`th initial vertex, det=%2.2g\n",
01229 qh_pointid(maxpoint), k+1, maxdet));
01230 }
01231 }
01232
01233
01234
01235
01236
01237
01238
01239 realT qh_minabsval (realT *normal, int dim) {
01240 realT minval= 0;
01241 realT maxval= 0;
01242 realT *colp;
01243 int k;
01244
01245 for (k= dim, colp= normal; k--; colp++) {
01246 maximize_(maxval, *colp);
01247 minimize_(minval, *colp);
01248 }
01249 return fmax_(maxval, -minval);
01250 }
01251
01252
01253
01254
01255
01256
01257
01258
01259 int qh_mindiff (realT *vecA, realT *vecB, int dim) {
01260 realT mindiff= REALmax, diff;
01261 realT *vecAp= vecA, *vecBp= vecB;
01262 int k, mink= 0;
01263
01264 for (k= 0; k < dim; k++) {
01265 diff= *vecAp++ - *vecBp++;
01266 diff= fabs_(diff);
01267 if (diff < mindiff) {
01268 mindiff= diff;
01269 mink= k;
01270 }
01271 }
01272 return mink;
01273 }
01274
01275
01276
01277
01278
01279
01280
01281
01282
01283
01284
01285
01286 boolT qh_orientoutside (facetT *facet) {
01287 int k;
01288 realT dist;
01289
01290 qh_distplane (qh interior_point, facet, &dist);
01291 if (dist > 0) {
01292 for (k= qh hull_dim; k--; )
01293 facet->normal[k]= -facet->normal[k];
01294 facet->offset= -facet->offset;
01295 return True;
01296 }
01297 return False;
01298 }
01299
01300
01301
01302
01303
01304
01305
01306
01307
01308
01309
01310
01311
01312
01313
01314
01315
01316
01317
01318
01319
01320 void qh_outerinner (facetT *facet, realT *outerplane, realT *innerplane) {
01321 realT dist, mindist;
01322 vertexT *vertex, **vertexp;
01323
01324 if (outerplane) {
01325 if (!qh_MAXoutside || !facet || !qh maxoutdone) {
01326 *outerplane= qh_maxouter();
01327 }else {
01328 #if qh_MAXoutside
01329 *outerplane= facet->maxoutside + qh DISTround;
01330 #endif
01331
01332 }
01333 if (qh JOGGLEmax < REALmax/2)
01334 *outerplane += qh JOGGLEmax * sqrt (qh hull_dim);
01335 }
01336 if (innerplane) {
01337 if (facet) {
01338 mindist= REALmax;
01339 FOREACHvertex_(facet->vertices) {
01340 zinc_(Zdistio);
01341 qh_distplane (vertex->point, facet, &dist);
01342 minimize_(mindist, dist);
01343 }
01344 *innerplane= mindist - qh DISTround;
01345 }else
01346 *innerplane= qh min_vertex - qh DISTround;
01347 if (qh JOGGLEmax < REALmax/2)
01348 *innerplane -= qh JOGGLEmax * sqrt (qh hull_dim);
01349 }
01350 }
01351
01352
01353
01354
01355
01356
01357
01358
01359
01360
01361 coordT qh_pointdist(pointT *point1, pointT *point2, int dim) {
01362 coordT dist, diff;
01363 int k;
01364
01365 dist= 0.0;
01366 for (k= (dim > 0 ? dim : -dim); k--; ) {
01367 diff= *point1++ - *point2++;
01368 dist += diff * diff;
01369 }
01370 if (dim > 0)
01371 return(sqrt(dist));
01372 return dist;
01373 }
01374
01375
01376
01377
01378
01379
01380
01381
01382
01383
01384
01385
01386 void qh_printmatrix (FILE *fp, char *string, realT **rows, int numrow, int numcol) {
01387 realT *rowp;
01388 realT r;
01389 int i,k;
01390
01391 fprintf (fp, "%s\n", string);
01392 for (i= 0; i < numrow; i++) {
01393 rowp= rows[i];
01394 for (k= 0; k < numcol; k++) {
01395 r= *rowp++;
01396 fprintf (fp, "%6.3g ", r);
01397 }
01398 fprintf (fp, "\n");
01399 }
01400 }
01401
01402
01403
01404
01405
01406
01407
01408
01409
01410 void qh_printpoints (FILE *fp, char *string, setT *points) {
01411 pointT *point, **pointp;
01412
01413 if (string) {
01414 fprintf (fp, "%s", string);
01415 FOREACHpoint_(points)
01416 fprintf (fp, " p%d", qh_pointid(point));
01417 fprintf (fp, "\n");
01418 }else {
01419 FOREACHpoint_(points)
01420 fprintf (fp, " %d", qh_pointid(point));
01421 fprintf (fp, "\n");
01422 }
01423 }
01424
01425
01426
01427
01428
01429
01430
01431
01432
01433
01434
01435
01436
01437
01438
01439
01440
01441
01442
01443
01444
01445
01446
01447
01448
01449
01450
01451
01452
01453
01454
01455
01456
01457
01458
01459
01460
01461
01462
01463
01464
01465 void qh_projectinput (void) {
01466 int k,i;
01467 int newdim= qh input_dim, newnum= qh num_points;
01468 signed char *project;
01469 int size= (qh input_dim+1)*sizeof(*project);
01470 pointT *newpoints, *coord, *infinity;
01471 realT paraboloid, maxboloid= 0;
01472
01473 project= (signed char*)qh_memalloc (size);
01474 memset ((char*)project, 0, size);
01475 for (k= 0; k < qh input_dim; k++) {
01476 if (qh lower_bound[k] == 0 && qh upper_bound[k] == 0) {
01477 project[k]= -1;
01478 newdim--;
01479 }
01480 }
01481 if (qh DELAUNAY) {
01482 project[k]= 1;
01483 newdim++;
01484 if (qh ATinfinity)
01485 newnum++;
01486 }
01487 if (newdim != qh hull_dim) {
01488 fprintf(qh ferr, "qhull internal error (qh_projectinput): dimension after projection %d != hull_dim %d\n", newdim, qh hull_dim);
01489 qh_errexit(qh_ERRqhull, NULL, NULL);
01490 }
01491 if (!(newpoints=(coordT*)malloc(newnum*newdim*sizeof(coordT)))){
01492 fprintf(qh ferr, "qhull error: insufficient memory to project %d points\n",
01493 qh num_points);
01494 qh_errexit(qh_ERRmem, NULL, NULL);
01495 }
01496 qh_projectpoints (project, qh input_dim+1, qh first_point,
01497 qh num_points, qh input_dim, newpoints, newdim);
01498 trace1((qh ferr, "qh_projectinput: updating lower and upper_bound\n"));
01499 qh_projectpoints (project, qh input_dim+1, qh lower_bound,
01500 1, qh input_dim+1, qh lower_bound, newdim+1);
01501 qh_projectpoints (project, qh input_dim+1, qh upper_bound,
01502 1, qh input_dim+1, qh upper_bound, newdim+1);
01503 if (qh HALFspace) {
01504 if (!qh feasible_point) {
01505 fprintf(qh ferr, "qhull internal error (qh_projectinput): HALFspace defined without qh.feasible_point\n");
01506 qh_errexit(qh_ERRqhull, NULL, NULL);
01507 }
01508 qh_projectpoints (project, qh input_dim, qh feasible_point,
01509 1, qh input_dim, qh feasible_point, newdim);
01510 }
01511 qh_memfree(project, ((qh input_dim+1)*sizeof(*project)));
01512 if (qh POINTSmalloc)
01513 free (qh first_point);
01514 qh first_point= newpoints;
01515 qh POINTSmalloc= True;
01516 if (qh DELAUNAY && qh ATinfinity) {
01517 coord= qh first_point;
01518 infinity= qh first_point + qh hull_dim * qh num_points;
01519 for (k=qh hull_dim-1; k--; )
01520 infinity[k]= 0.0;
01521 for (i=qh num_points; i--; ) {
01522 paraboloid= 0.0;
01523 for (k=qh hull_dim-1; k--; ) {
01524 paraboloid += *coord * *coord;
01525 infinity[k] += *coord;
01526 coord++;
01527 }
01528 *(coord++)= paraboloid;
01529 maximize_(maxboloid, paraboloid);
01530 }
01531
01532 for (k=qh hull_dim-1; k--; )
01533 *(coord++) /= qh num_points;
01534 *(coord++)= maxboloid * 1.1;
01535 qh num_points++;
01536 trace0((qh ferr, "qh_projectinput: projected points to paraboloid for Delaunay\n"));
01537 }else if (qh DELAUNAY)
01538 qh_setdelaunay( qh hull_dim, qh num_points, qh first_point);
01539 }
01540
01541
01542
01543
01544
01545
01546
01547
01548
01549
01550
01551
01552
01553
01554
01555
01556
01557
01558
01559
01560
01561
01562
01563
01564
01565
01566
01567 void qh_projectpoints (signed char *project, int n, realT *points,
01568 int numpoints, int dim, realT *newpoints, int newdim) {
01569 int testdim= dim, oldk=0, newk=0, i,j=0,k;
01570 realT *newp, *oldp;
01571
01572 for (k= 0; k < n; k++)
01573 testdim += project[k];
01574 if (testdim != newdim) {
01575 fprintf (qh ferr, "qhull internal error (qh_projectpoints): newdim %d should be %d after projection\n",
01576 newdim, testdim);
01577 qh_errexit (qh_ERRqhull, NULL, NULL);
01578 }
01579 for (j= 0; j<n; j++) {
01580 if (project[j] == -1)
01581 oldk++;
01582 else {
01583 newp= newpoints+newk++;
01584 if (project[j] == +1) {
01585 if (oldk >= dim)
01586 continue;
01587 oldp= points+oldk;
01588 }else
01589 oldp= points+oldk++;
01590 for (i=numpoints; i--; ) {
01591 *newp= *oldp;
01592 newp += newdim;
01593 oldp += dim;
01594 }
01595 }
01596 if (oldk >= dim)
01597 break;
01598 }
01599 trace1((qh ferr, "qh_projectpoints: projected %d points from dim %d to dim %d\n",
01600 numpoints, dim, newdim));
01601 }
01602
01603
01604
01605
01606
01607
01608
01609
01610
01611
01612
01613
01614
01615
01616
01617
01618 int qh_rand_seed= 1;
01619
01620 int qh_rand( void) {
01621 #define qh_rand_a 16807
01622 #define qh_rand_m 2147483647
01623 #define qh_rand_q 127773
01624 #define qh_rand_r 2836
01625 int lo, hi, test;
01626 int seed = qh_rand_seed;
01627
01628 hi = seed / qh_rand_q;
01629 lo = seed % qh_rand_q;
01630 test = qh_rand_a * lo - qh_rand_r * hi;
01631 if (test > 0)
01632 seed= test;
01633 else
01634 seed= test + qh_rand_m;
01635 qh_rand_seed= seed;
01636
01637
01638 return seed;
01639 }
01640
01641 void qh_srand( int seed) {
01642 if (seed < 1)
01643 qh_rand_seed= 1;
01644 else if (seed >= qh_rand_m)
01645 qh_rand_seed= qh_rand_m - 1;
01646 else
01647 qh_rand_seed= seed;
01648 }
01649
01650
01651
01652
01653
01654
01655
01656
01657
01658
01659 realT qh_randomfactor (void) {
01660 realT randr;
01661
01662 randr= qh_RANDOMint;
01663 return randr * qh RANDOMa + qh RANDOMb;
01664 }
01665
01666
01667
01668
01669
01670
01671
01672
01673
01674
01675
01676
01677
01678 void qh_randommatrix (realT *buffer, int dim, realT **rows) {
01679 int i, k;
01680 realT **rowi, *coord, realr;
01681
01682 coord= buffer;
01683 rowi= rows;
01684 for (i=0; i < dim; i++) {
01685 *(rowi++)= coord;
01686 for (k=0; k < dim; k++) {
01687 realr= qh_RANDOMint;
01688 *(coord++)= 2.0 * realr/(qh_RANDOMmax+1) - 1.0;
01689 }
01690 }
01691 *rowi= coord;
01692 }
01693
01694
01695
01696
01697
01698
01699
01700
01701
01702
01703
01704
01705
01706
01707
01708
01709
01710
01711 void qh_rotateinput (realT **rows) {
01712
01713 if (!qh POINTSmalloc) {
01714 qh first_point= qh_copypoints (qh first_point, qh num_points, qh hull_dim);
01715 qh POINTSmalloc= True;
01716 }
01717 qh_rotatepoints (qh first_point, qh num_points, qh hull_dim, rows);
01718 }
01719
01720
01721
01722
01723
01724
01725
01726
01727
01728
01729
01730
01731
01732
01733
01734
01735
01736
01737 void qh_rotatepoints (realT *points, int numpoints, int dim, realT **row) {
01738 realT *point, *rowi, *coord= NULL, sum, *newval;
01739 int i,j,k;
01740
01741 if (qh IStracing >= 1)
01742 qh_printmatrix (qh ferr, "qh_rotatepoints: rotate points by", row, dim, dim);
01743 for (point= points, j= numpoints; j--; point += dim) {
01744 newval= row[dim];
01745 for (i= 0; i < dim; i++) {
01746 rowi= row[i];
01747 coord= point;
01748 for (sum= 0.0, k= dim; k--; )
01749 sum += *rowi++ * *coord++;
01750 *(newval++)= sum;
01751 }
01752 for (k= dim; k--; )
01753 *(--coord)= *(--newval);
01754 }
01755 }
01756
01757
01758
01759
01760
01761
01762
01763
01764
01765
01766
01767
01768
01769
01770
01771
01772
01773 void qh_scaleinput (void) {
01774
01775 if (!qh POINTSmalloc) {
01776 qh first_point= qh_copypoints (qh first_point, qh num_points, qh hull_dim);
01777 qh POINTSmalloc= True;
01778 }
01779 qh_scalepoints (qh first_point, qh num_points, qh hull_dim,
01780 qh lower_bound, qh upper_bound);
01781 }
01782
01783
01784
01785
01786
01787
01788
01789
01790
01791
01792
01793
01794
01795
01796
01797
01798
01799
01800
01801
01802 void qh_scalelast (coordT *points, int numpoints, int dim, coordT low,
01803 coordT high, coordT newhigh) {
01804 realT scale, shift;
01805 coordT *coord;
01806 int i;
01807 boolT nearzero= False;
01808
01809 trace4((qh ferr, "qh_scalelast: scale last coordinate from [%2.2g, %2.2g] to [0,%2.2g]\n",
01810 low, high, newhigh));
01811 qh last_low= low;
01812 qh last_high= high;
01813 qh last_newhigh= newhigh;
01814 scale= qh_divzero (newhigh, high - low,
01815 qh MINdenom_1, &nearzero);
01816 if (nearzero) {
01817 if (qh DELAUNAY)
01818 fprintf (qh ferr, "qhull input error: can not scale last coordinate. Input is cocircular\n or cospherical. Use option 'Qz' to add a point at infinity.\n");
01819 else
01820 fprintf (qh ferr, "qhull input error: can not scale last coordinate. New bounds [0, %2.2g] are too wide for\nexisting bounds [%2.2g, %2.2g] (width %2.2g)\n",
01821 newhigh, low, high, high-low);
01822 qh_errexit (qh_ERRinput, NULL, NULL);
01823 }
01824 shift= - low * newhigh / (high-low);
01825 coord= points + dim - 1;
01826 for (i= numpoints; i--; coord += dim)
01827 *coord= *coord * scale + shift;
01828 }
01829
01830
01831
01832
01833
01834
01835
01836
01837
01838
01839
01840
01841
01842
01843
01844
01845
01846
01847
01848 void qh_scalepoints (pointT *points, int numpoints, int dim,
01849 realT *newlows, realT *newhighs) {
01850 int i,k;
01851 realT shift, scale, *coord, low, high, newlow, newhigh, mincoord, maxcoord;
01852 boolT nearzero= False;
01853
01854 for (k= 0; k < dim; k++) {
01855 newhigh= newhighs[k];
01856 newlow= newlows[k];
01857 if (newhigh > REALmax/2 && newlow < -REALmax/2)
01858 continue;
01859 low= REALmax;
01860 high= -REALmax;
01861 for (i= numpoints, coord= points+k; i--; coord += dim) {
01862 minimize_(low, *coord);
01863 maximize_(high, *coord);
01864 }
01865 if (newhigh > REALmax/2)
01866 newhigh= high;
01867 if (newlow < -REALmax/2)
01868 newlow= low;
01869 if (qh DELAUNAY && k == dim-1 && newhigh < newlow) {
01870 fprintf (qh ferr, "qhull input error: 'Qb%d' or 'QB%d' inverts paraboloid since high bound %.2g < low bound %.2g\n",
01871 k, k, newhigh, newlow);
01872 qh_errexit (qh_ERRinput, NULL, NULL);
01873 }
01874 scale= qh_divzero (newhigh - newlow, high - low,
01875 qh MINdenom_1, &nearzero);
01876 if (nearzero) {
01877 fprintf (qh ferr, "qhull input error: %d'th dimension's new bounds [%2.2g, %2.2g] too wide for\nexisting bounds [%2.2g, %2.2g]\n",
01878 k, newlow, newhigh, low, high);
01879 qh_errexit (qh_ERRinput, NULL, NULL);
01880 }
01881 shift= (newlow * high - low * newhigh)/(high-low);
01882 coord= points+k;
01883 for (i= numpoints; i--; coord += dim)
01884 *coord= *coord * scale + shift;
01885 coord= points+k;
01886 if (newlow < newhigh) {
01887 mincoord= newlow;
01888 maxcoord= newhigh;
01889 }else {
01890 mincoord= newhigh;
01891 maxcoord= newlow;
01892 }
01893 for (i= numpoints; i--; coord += dim) {
01894 minimize_(*coord, maxcoord);
01895 maximize_(*coord, mincoord);
01896 }
01897 trace0((qh ferr, "qh_scalepoints: scaled %d'th coordinate [%2.2g, %2.2g] to [%.2g, %.2g] for %d points by %2.2g and shifted %2.2g\n",
01898 k, low, high, newlow, newhigh, numpoints, scale, shift));
01899 }
01900 }
01901
01902
01903
01904
01905
01906
01907
01908
01909
01910
01911
01912
01913
01914
01915
01916
01917
01918
01919
01920
01921
01922
01923
01924
01925
01926
01927
01928
01929
01930
01931
01932
01933 void qh_setdelaunay (int dim, int count, pointT *points) {
01934 int i, k;
01935 coordT *coordp, coord;
01936 realT paraboloid;
01937
01938 trace0((qh ferr, "qh_setdelaunay: project %d points to paraboloid for Delaunay triangulation\n", count));
01939 coordp= points;
01940 for (i= 0; i < count; i++) {
01941 coord= *coordp++;
01942 paraboloid= coord*coord;
01943 for (k= dim-2; k--; ) {
01944 coord= *coordp++;
01945 paraboloid += coord*coord;
01946 }
01947 *coordp++ = paraboloid;
01948 }
01949 if (qh last_low < REALmax/2)
01950 qh_scalelast (points, count, dim, qh last_low, qh last_high, qh last_newhigh);
01951 }
01952
01953
01954
01955
01956
01957
01958
01959
01960
01961
01962
01963
01964
01965
01966
01967
01968
01969
01970 boolT qh_sethalfspace (int dim, coordT *coords, coordT **nextp,
01971 coordT *normal, coordT *offset, coordT *feasible) {
01972 coordT *normp= normal, *feasiblep= feasible, *coordp= coords;
01973 realT dist;
01974 realT r;
01975 int k;
01976 boolT zerodiv;
01977
01978 dist= *offset;
01979 for (k= dim; k--; )
01980 dist += *(normp++) * *(feasiblep++);
01981 if (dist > 0)
01982 goto LABELerroroutside;
01983 normp= normal;
01984 if (dist < -qh MINdenom) {
01985 for (k= dim; k--; )
01986 *(coordp++)= *(normp++) / -dist;
01987 }else {
01988 for (k= dim; k--; ) {
01989 *(coordp++)= qh_divzero (*(normp++), -dist, qh MINdenom_1, &zerodiv);
01990 if (zerodiv)
01991 goto LABELerroroutside;
01992 }
01993 }
01994 *nextp= coordp;
01995 if (qh IStracing >= 4) {
01996 fprintf (qh ferr, "qh_sethalfspace: halfspace at offset %6.2g to point: ", *offset);
01997 for (k= dim, coordp= coords; k--; ) {
01998 r= *coordp++;
01999 fprintf (qh ferr, " %6.2g", r);
02000 }
02001 fprintf (qh ferr, "\n");
02002 }
02003 return True;
02004 LABELerroroutside:
02005 feasiblep= feasible;
02006 normp= normal;
02007 fprintf(qh ferr, "qhull input error: feasible point is not clearly inside halfspace\nfeasible point: ");
02008 for (k= dim; k--; )
02009 fprintf (qh ferr, qh_REAL_1, r=*(feasiblep++));
02010 fprintf (qh ferr, "\n halfspace: ");
02011 for (k= dim; k--; )
02012 fprintf (qh ferr, qh_REAL_1, r=*(normp++));
02013 fprintf (qh ferr, "\n at offset: ");
02014 fprintf (qh ferr, qh_REAL_1, *offset);
02015 fprintf (qh ferr, " and distance: ");
02016 fprintf (qh ferr, qh_REAL_1, dist);
02017 fprintf (qh ferr, "\n");
02018 return False;
02019 }
02020
02021
02022
02023
02024
02025
02026
02027
02028
02029
02030
02031
02032
02033
02034
02035
02036
02037
02038
02039
02040
02041
02042 coordT *qh_sethalfspace_all (int dim, int count, coordT *halfspaces, pointT *feasible) {
02043 int i, newdim;
02044 pointT *newpoints;
02045 coordT *coordp, *normalp, *offsetp;
02046
02047 trace0((qh ferr, "qh_sethalfspace_all: compute dual for halfspace intersection\n"));
02048 newdim= dim - 1;
02049 if (!(newpoints=(coordT*)malloc(count*newdim*sizeof(coordT)))){
02050 fprintf(qh ferr, "qhull error: insufficient memory to compute dual of %d halfspaces\n",
02051 count);
02052 qh_errexit(qh_ERRmem, NULL, NULL);
02053 }
02054 coordp= newpoints;
02055 normalp= halfspaces;
02056 for (i= 0; i < count; i++) {
02057 offsetp= normalp + newdim;
02058 if (!qh_sethalfspace (newdim, coordp, &coordp, normalp, offsetp, feasible)) {
02059 fprintf (qh ferr, "The halfspace was at index %d\n", i);
02060 qh_errexit (qh_ERRinput, NULL, NULL);
02061 }
02062 normalp= offsetp + 1;
02063 }
02064 return newpoints;
02065 }
02066
02067
02068
02069
02070
02071
02072
02073
02074
02075
02076
02077
02078
02079
02080
02081
02082
02083
02084
02085
02086
02087
02088
02089
02090
02091
02092 pointT *qh_voronoi_center (int dim, setT *points) {
02093 pointT *point, **pointp, *point0;
02094 pointT *center= (pointT*)qh_memalloc (qh center_size);
02095 setT *simplex;
02096 int i, j, k, size= qh_setsize(points);
02097 coordT *gmcoord;
02098 realT *diffp, sum2, *sum2row, *sum2p, det, factor;
02099 boolT nearzero, infinite;
02100
02101 if (size == dim+1)
02102 simplex= points;
02103 else if (size < dim+1) {
02104 fprintf (qh ferr, "qhull internal error (qh_voronoi_center):\n need at least %d points to construct a Voronoi center\n",
02105 dim+1);
02106 qh_errexit (qh_ERRqhull, NULL, NULL);
02107 }else {
02108 simplex= qh_settemp (dim+1);
02109 qh_maxsimplex (dim, points, NULL, 0, &simplex);
02110 }
02111 point0= SETfirstt_(simplex, pointT);
02112 gmcoord= qh gm_matrix;
02113 for (k=0; k < dim; k++) {
02114 qh gm_row[k]= gmcoord;
02115 FOREACHpoint_(simplex) {
02116 if (point != point0)
02117 *(gmcoord++)= point[k] - point0[k];
02118 }
02119 }
02120 sum2row= gmcoord;
02121 for (i=0; i < dim; i++) {
02122 sum2= 0.0;
02123 for (k= 0; k < dim; k++) {
02124 diffp= qh gm_row[k] + i;
02125 sum2 += *diffp * *diffp;
02126 }
02127 *(gmcoord++)= sum2;
02128 }
02129 det= qh_determinant (qh gm_row, dim, &nearzero);
02130 factor= qh_divzero (0.5, det, qh MINdenom, &infinite);
02131 if (infinite) {
02132 for (k=dim; k--; )
02133 center[k]= qh_INFINITE;
02134 if (qh IStracing)
02135 qh_printpoints (qh ferr, "qh_voronoi_center: at infinity for ", simplex);
02136 }else {
02137 for (i=0; i < dim; i++) {
02138 gmcoord= qh gm_matrix;
02139 sum2p= sum2row;
02140 for (k=0; k < dim; k++) {
02141 qh gm_row[k]= gmcoord;
02142 if (k == i) {
02143 for (j= dim; j--; )
02144 *(gmcoord++)= *sum2p++;
02145 }else {
02146 FOREACHpoint_(simplex) {
02147 if (point != point0)
02148 *(gmcoord++)= point[k] - point0[k];
02149 }
02150 }
02151 }
02152 center[i]= qh_determinant (qh gm_row, dim, &nearzero)*factor + point0[i];
02153 }
02154 #ifndef qh_NOtrace
02155 if (qh IStracing >= 3) {
02156 fprintf (qh ferr, "qh_voronoi_center: det %2.2g factor %2.2g ", det, factor);
02157 qh_printmatrix (qh ferr, "center:", ¢er, 1, dim);
02158 if (qh IStracing >= 5) {
02159 qh_printpoints (qh ferr, "points", simplex);
02160 FOREACHpoint_(simplex)
02161 fprintf (qh ferr, "p%d dist %.2g, ", qh_pointid (point),
02162 qh_pointdist (point, center, dim));
02163 fprintf (qh ferr, "\n");
02164 }
02165 }
02166 #endif
02167 }
02168 if (simplex != points)
02169 qh_settempfree (&simplex);
02170 return center;
02171 }
02172