34 size= numpoints * dimension * (int)
sizeof(
coordT);
36 qh_fprintf(
qh,
qh->ferr, 6004,
"qhull error: insufficient memory to copy %d points\n",
40 memcpy((
char *)newpoints, (
char *)points, (
size_t)size);
58 vecC[0]=
det2_(vecA[1], vecA[2],
60 vecC[1]= -
det2_(vecA[0], vecA[2],
62 vecC[2]=
det2_(vecA[0], vecA[1],
90 qh_fprintf(
qh,
qh->ferr, 6005,
"qhull internal error (qh_determinate): only implemented for dimension >= 2\n");
95 if (
fabs_(det) < 10*
qh->NEARzero[1])
101 if (
fabs_(det) < 10*
qh->NEARzero[2])
131 realT abscoord, distround, joggle, maxcoord, mincoord;
132 pointT *point, *pointtemp;
138 for (k=0; k < dimension; k++) {
139 if (
qh->SCALElast && k == dimension-1)
141 else if (
qh->DELAUNAY && k == dimension-1)
142 abscoord= 2 * maxabs * maxabs;
151 abscoord=
fmax_(maxcoord, -mincoord);
159 trace2((
qh,
qh->ferr, 2001,
"qh_detjoggle: joggle=%2.2g maxwidth=%2.2g\n", joggle, maxwidth));
196 if (!
qh->SETroundoff) {
199 qh->DISTround +=
qh->RANDOMfactor *
qh->MAXabs_coord;
202 qh->MINdenom=
qh->MINdenom_1 *
qh->MAXabs_coord;
203 qh->MINdenom_1_2= sqrt(
qh->MINdenom_1 *
qh->hull_dim) ;
204 qh->MINdenom_2=
qh->MINdenom_1_2 *
qh->MAXabs_coord;
208 qh->ANGLEround +=
qh->RANDOMfactor;
210 qh->premerge_cos -=
qh->ANGLEround;
212 qh_option(
qh,
"Angle-premerge-with-random", NULL, &
qh->premerge_cos);
215 qh->postmerge_cos -=
qh->ANGLEround;
217 qh_option(
qh,
"Angle-postmerge-with-random", NULL, &
qh->postmerge_cos);
219 qh->premerge_centrum += 2 *
qh->DISTround;
220 qh->postmerge_centrum += 2 *
qh->DISTround;
221 if (
qh->RANDOMdist && (
qh->MERGEexact ||
qh->PREmerge))
222 qh_option(
qh,
"Centrum-premerge-with-random", NULL, &
qh->premerge_centrum);
223 if (
qh->RANDOMdist &&
qh->POSTmerge)
224 qh_option(
qh,
"Centrum-postmerge-with-random", NULL, &
qh->postmerge_centrum);
226 realT maxangle= 1.0, maxrho;
231 qh->ONEmerge= sqrt((
realT)
qh->hull_dim) *
qh->MAXwidth *
232 sqrt(1.0 - maxangle * maxangle) +
qh->DISTround;
233 maxrho=
qh->hull_dim *
qh->premerge_centrum +
qh->DISTround;
235 maxrho=
qh->hull_dim *
qh->postmerge_centrum +
qh->DISTround;
241 if (
qh->JOGGLEmax <
REALmax/2 && (
qh->KEEPcoplanar ||
qh->KEEPinside)) {
244 maxdist= sqrt((
realT)
qh->hull_dim) *
qh->JOGGLEmax +
qh->DISTround;
248 if (
qh->KEEPnearinside)
250 if (
qh->JOGGLEmax <
qh->DISTround) {
251 qh_fprintf(
qh,
qh->ferr, 6006,
"qhull error: the joggle for 'QJn', %.2g, is below roundoff for distance computations, %.2g\n",
252 qh->JOGGLEmax,
qh->DISTround);
257 qh->MINvisible=
qh->DISTround;
258 else if (
qh->hull_dim <= 3)
259 qh->MINvisible=
qh->premerge_centrum;
262 if (
qh->APPROXhull &&
qh->MINvisible >
qh->MINoutside)
263 qh->MINvisible=
qh->MINoutside;
267 qh->MAXcoplanar=
qh->MINvisible;
268 qh_option(
qh,
"U-coplanar-distance", NULL, &
qh->MAXcoplanar);
270 if (!
qh->APPROXhull) {
271 qh->MINoutside= 2 *
qh->MINvisible;
276 qh->WIDEfacet=
qh->MINoutside;
281 && !
qh->BESToutside && !
qh->FORCEoutput)
282 qh_fprintf(
qh,
qh->ferr, 7001,
"qhull input warning: minimum visibility V%.2g is greater than \nminimum outside W%.2g. Flipped facets are likely.\n",
283 qh->MINvisible,
qh->MINoutside);
284 qh->max_vertex=
qh->DISTround;
285 qh->min_vertex= -
qh->DISTround;
306 pointT *coorda, *coordp, *gmcoord, *point, **pointp;
312 gmcoord=
qh->gm_matrix;
321 *(gmcoord++)= *coordp++ - *coorda++;
324 qh_fprintf(
qh,
qh->ferr, 6007,
"qhull internal error (qh_detsimplex): #points %d < dimension %d\n",
329 trace2((
qh,
qh->ferr, 2002,
"qh_detsimplex: det=%2.2g for point p%d, dim %d, nearzero? %d\n",
350 coordT *normalp= normal, *coordp= point;
356 dist += *(coordp++) * *(normalp++);
378 realT maxdistsum, maxround;
380 maxdistsum= sqrt((
realT)dimension) * maxabs;
382 maxround=
REALepsilon * (dimension * maxdistsum * 1.01 + maxabs);
384 trace4((
qh,
qh->ferr, 4008,
"qh_distround: %2.2g maxabs %2.2g maxsumabs %2.2g maxdistsum %2.2g\n",
385 maxround, maxabs, maxsumabs, maxdistsum));
411 realT temp, numerx, denomx;
414 if (numer < mindenom1 && numer > -mindenom1) {
415 numerx=
fabs_(numer);
416 denomx=
fabs_(denom);
417 if (numerx < denomx) {
426 if (temp > mindenom1 || temp < -mindenom1) {
481 trace4((
qh,
qh->ferr, 4009,
"qh_facetarea: f%d area %2.2g\n", facet->
id, area));
520 pointT *coorda, *coordp, *gmcoord;
527 gmcoord=
qh->gm_matrix;
530 if (vertex == notvertex)
534 coordp= vertex->
point;
538 *(gmcoord++)= *coordp++ - *coorda++;
542 dist += *coordp++ * *normalp++;
543 if (dist < -qh->WIDEfacet) {
547 coordp= vertex->
point;
550 *(gmcoord++)= (*coordp++ - dist * *normalp++) - *coorda++;
554 qh_fprintf(
qh,
qh->ferr, 6008,
"qhull internal error (qh_facetarea_simplex): #points %d != dim %d -1\n",
560 for (i=0; i <
dim-1; i++)
568 *(gmcoord++)= *normalp++;
574 area *=
qh->AREAfactor;
575 trace4((
qh,
qh->ferr, 4010,
"qh_facetarea_simplex: area=%2.2g for point p%d, toporient %d, nearzero? %d\n",
634 facetT *neighbor, **neighborp, *bestfacet=NULL, *facet;
652 if (goodseen && !neighbor->
good)
659 if (neighbor->
good) {
661 if (dist > bestdist) {
671 trace2((
qh,
qh->ferr, 2003,
"qh_findgooddist: p%d is %2.2g above good facet f%d\n",
675 trace4((
qh,
qh->ferr, 4011,
"qh_findgooddist: no good facet for p%d above f%d\n",
708 if (
qh->hasAreaVolume)
711 qh_fprintf(
qh,
qh->ferr, 8020,
"computing area of each facet and volume of the convex hull\n");
713 trace1((
qh,
qh->ferr, 1001,
"qh_getarea: computing volume and area for each facet\n"));
714 qh->totarea=
qh->totvol= 0.0;
731 qh->totvol += -dist * area/
qh->hull_dim;
733 if (
qh->PRINTstatistics) {
765 realT *rowi, *rowj, norm;
768 for (i=0; i <
dim; i++) {
770 for (norm= 0.0, k=
dim; k--; rowi++)
771 norm += *rowi * *rowi;
778 for (j=i+1; j <
dim; j++) {
780 for (norm= 0.0, k=
dim; k--; )
781 norm += *rowi++ * *rowj++;
783 *(--rowj) -= *(--rowi) * norm;
819 for (k=0; k <
qh->hull_dim; k++) {
820 threshold=
qh->lower_threshold[k];
822 if (normal[k] < threshold)
825 threshold -= normal[k];
826 *angle +=
fabs_(threshold);
829 if (
qh->upper_threshold[k] <
REALmax/2) {
830 threshold=
qh->upper_threshold[k];
831 if (normal[k] > threshold)
834 threshold -= normal[k];
835 *angle +=
fabs_(threshold);
878 realT randr, randa, randb;
880 if (!
qh->input_points) {
881 qh->input_points=
qh->first_point;
882 qh->input_malloc=
qh->POINTSmalloc;
883 size=
qh->num_points *
qh->hull_dim *
sizeof(
coordT);
885 qh_fprintf(
qh,
qh->ferr, 6009,
"qhull error: insufficient memory to joggle %d points\n",
890 if (
qh->JOGGLEmax == 0.0) {
898 if (
qh->JOGGLEmax < maxjoggle) {
906 if (
qh->build_cnt > 1 &&
qh->JOGGLEmax >
fmax_(
qh->MAXwidth/4, 0.1)) {
907 qh_fprintf(
qh,
qh->ferr, 6010,
"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",
914 trace0((
qh,
qh->ferr, 6,
"qh_joggleinput: joggle input by %2.2g with seed %d\n",
916 inputp=
qh->input_points;
917 coordp=
qh->first_point;
919 randb= -
qh->JOGGLEmax;
920 size=
qh->num_points *
qh->hull_dim;
921 for (i=size; i--; ) {
923 *(coordp++)= *(inputp++) + (randr * randa + randb);
940 realT *maxp= NULL, *colp, absval;
943 for (k=
dim, colp= normal; k--; colp++) {
944 absval=
fabs_(*colp);
945 if (absval > maxval) {
982 realT maxcoord, temp;
983 pointT *minimum, *maximum, *point, *pointtemp;
986 qh->max_outside= 0.0;
987 qh->MAXabs_coord= 0.0;
989 qh->MAXsumcoord= 0.0;
998 qh_fprintf(
qh,
qh->ferr, 6011,
"qhull error: floating point constants in user.h are wrong\n\
999 REALepsilon %g REALmin %g REALmax %g -REALmax %g\n",
1004 for (k=0; k < dimension; k++) {
1005 if (points ==
qh->GOODpointp)
1006 minimum= maximum= points + dimension;
1008 minimum= maximum= points;
1010 if (point ==
qh->GOODpointp)
1012 if (maximum[k] < point[k])
1014 else if (minimum[k] > point[k])
1017 if (k == dimension-1) {
1018 qh->MINlastcoord= minimum[k];
1019 qh->MAXlastcoord= maximum[k];
1021 if (
qh->SCALElast && k == dimension-1)
1022 maxcoord=
qh->MAXwidth;
1024 maxcoord=
fmax_(maximum[k], -minimum[k]);
1025 if (
qh->GOODpointp) {
1026 temp=
fmax_(
qh->GOODpointp[k], -
qh->GOODpointp[k]);
1029 temp= maximum[k] - minimum[k];
1033 qh->MAXsumcoord += maxcoord;
1042 if (
qh->IStracing >=1)
1043 qh_printpoints(
qh,
qh->ferr,
"qh_maxmin: found the max and min points(by dim):", set);
1070 dist=
fmax_(
qh->max_outside,
qh->DISTround);
1071 dist +=
qh->DISTround;
1072 trace4((
qh,
qh->ferr, 4012,
"qh_maxouter: max distance from facet to outer plane is %2.2g max_outside is %2.2g\n", dist,
qh->max_outside));
1100 pointT *point, **pointp, *pointtemp, *maxpoint, *minx=NULL, *maxx=NULL;
1109 if (maxcoord < point[0]) {
1113 if (mincoord > point[0]) {
1120 if (point ==
qh->GOODpointp)
1122 if (maxcoord < point[0]) {
1126 if (mincoord > point[0]) {
1139 qh_fprintf(
qh,
qh->ferr, 6012,
"qhull precision error (qh_maxsimplex for voronoi_center):\n%d points with the same x coordinate.\n",
1143 qh_fprintf(
qh,
qh->ferr, 6013,
"qhull input error: input is less than %d-dimensional since it has the same x coordinate\n",
qh->hull_dim);
1148 for (k=sizinit; k <
dim+1; k++) {
1154 if ((det=
fabs_(det)) > maxdet) {
1157 maxnearzero= nearzero;
1161 if (!maxpoint || maxnearzero) {
1164 trace0((
qh,
qh->ferr, 7,
"qh_maxsimplex: searching all points for %d-th initial vertex.\n", k+1));
1166 trace0((
qh,
qh->ferr, 8,
"qh_maxsimplex: searching all points for %d-th initial vertex, better than p%d det %2.2g\n",
1170 if (point ==
qh->GOODpointp)
1174 if ((det=
fabs_(det)) > maxdet) {
1177 maxnearzero= nearzero;
1183 qh_fprintf(
qh,
qh->ferr, 6014,
"qhull internal error (qh_maxsimplex): not enough points available\n");
1187 trace1((
qh,
qh->ferr, 1002,
"qh_maxsimplex: selected point p%d for %d`th initial vertex, det=%2.2g\n",
1204 for (k=
dim, colp=normal; k--; colp++) {
1208 return fmax_(maxval, -minval);
1220 realT *vecAp= vecA, *vecBp= vecB;
1223 for (k=0; k <
dim; k++) {
1224 diff= *vecAp++ - *vecBp++;
1226 if (diff < mindiff) {
1251 for (k=
qh->hull_dim; k--; )
1281 realT dist, mindist;
1294 *outerplane +=
qh->JOGGLEmax * sqrt((
realT)
qh->hull_dim);
1304 *innerplane= mindist -
qh->DISTround;
1306 *innerplane=
qh->min_vertex -
qh->DISTround;
1308 *innerplane -=
qh->JOGGLEmax * sqrt((
realT)
qh->hull_dim);
1326 for (k= (
dim > 0 ?
dim : -
dim); k--; ) {
1327 diff= *point1++ - *point2++;
1328 dist += diff * diff;
1353 for (i=0; i < numrow; i++) {
1355 for (k=0; k < numcol; k++) {
1428 int newdim=
qh->input_dim, newnum=
qh->num_points;
1429 signed char *project;
1430 int projectsize= (
qh->input_dim+1)*
sizeof(*project);
1431 pointT *newpoints, *coord, *infinity;
1432 realT paraboloid, maxboloid= 0;
1435 memset((
char*)project, 0, (
size_t)projectsize);
1436 for (k=0; k <
qh->input_dim; k++) {
1437 if (
qh->lower_bound[k] == 0 &&
qh->upper_bound[k] == 0) {
1448 if (newdim !=
qh->hull_dim) {
1450 qh_fprintf(
qh,
qh->ferr, 6015,
"qhull internal error (qh_projectinput): dimension after projection %d != hull_dim %d\n", newdim,
qh->hull_dim);
1455 qh_fprintf(
qh,
qh->ferr, 6016,
"qhull error: insufficient memory to project %d points\n",
1461 qh->num_points,
qh->input_dim, newpoints, newdim);
1462 trace1((
qh,
qh->ferr, 1003,
"qh_projectinput: updating lower and upper_bound\n"));
1464 1,
qh->input_dim+1,
qh->lower_bound, newdim+1);
1466 1,
qh->input_dim+1,
qh->upper_bound, newdim+1);
1467 if (
qh->HALFspace) {
1468 if (!
qh->feasible_point) {
1470 qh_fprintf(
qh,
qh->ferr, 6017,
"qhull internal error (qh_projectinput): HALFspace defined without qh.feasible_point\n");
1474 1,
qh->input_dim,
qh->feasible_point, newdim);
1477 if (
qh->POINTSmalloc)
1479 qh->first_point= newpoints;
1481 qh->temp_malloc= NULL;
1482 if (
qh->DELAUNAY &&
qh->ATinfinity) {
1483 coord=
qh->first_point;
1484 infinity=
qh->first_point +
qh->hull_dim *
qh->num_points;
1485 for (k=
qh->hull_dim-1; k--; )
1487 for (i=
qh->num_points; i--; ) {
1489 for (k=0; k <
qh->hull_dim-1; k++) {
1490 paraboloid += *coord * *coord;
1491 infinity[k] += *coord;
1494 *(coord++)= paraboloid;
1498 for (k=
qh->hull_dim-1; k--; )
1499 *(coord++) /=
qh->num_points;
1500 *(coord++)= maxboloid * 1.1;
1502 trace0((
qh,
qh->ferr, 9,
"qh_projectinput: projected points to paraboloid for Delaunay\n"));
1503 }
else if (
qh->DELAUNAY)
1534 int numpoints,
int dim,
realT *newpoints,
int newdim) {
1535 int testdim=
dim, oldk=0, newk=0, i,j=0,k;
1538 for (k=0; k < n; k++)
1539 testdim += project[k];
1540 if (testdim != newdim) {
1541 qh_fprintf(
qh,
qh->ferr, 6018,
"qhull internal error (qh_projectpoints): newdim %d should be %d after projection\n",
1545 for (j=0; j<n; j++) {
1546 if (project[j] == -1)
1549 newp= newpoints+newk++;
1550 if (project[j] == +1) {
1555 oldp= points+oldk++;
1556 for (i=numpoints; i--; ) {
1565 trace1((
qh,
qh->ferr, 1004,
"qh_projectpoints: projected %d points from dim %d to dim %d\n",
1566 numpoints,
dim, newdim));
1588 if (!
qh->POINTSmalloc) {
1613 realT *point, *rowi, *coord= NULL, sum, *newval;
1616 if (
qh->IStracing >= 1)
1618 for (point= points, j= numpoints; j--; point +=
dim) {
1620 for (i=0; i <
dim; i++) {
1623 for (sum= 0.0, k=
dim; k--; )
1624 sum += *rowi++ * *coord++;
1628 *(--coord)= *(--newval);
1650 if (!
qh->POINTSmalloc) {
1655 qh->lower_bound,
qh->upper_bound);
1684 trace4((
qh,
qh->ferr, 4013,
"qh_scalelast: scale last coordinate from [%2.2g, %2.2g] to [0,%2.2g]\n",
1685 low, high, newhigh));
1687 qh->last_high= high;
1688 qh->last_newhigh= newhigh;
1690 qh->MINdenom_1, &nearzero);
1693 qh_fprintf(
qh,
qh->ferr, 6019,
"qhull input error: can not scale last coordinate. Input is cocircular\n or cospherical. Use option 'Qz' to add a point at infinity.\n");
1695 qh_fprintf(
qh,
qh->ferr, 6020,
"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",
1696 newhigh, low, high, high-low);
1699 shift= - low * newhigh / (high-low);
1700 coord= points +
dim - 1;
1701 for (i=numpoints; i--; coord +=
dim)
1702 *coord= *coord * scale + shift;
1726 realT shift, scale, *coord, low, high, newlow, newhigh, mincoord, maxcoord;
1729 for (k=0; k <
dim; k++) {
1730 newhigh= newhighs[k];
1736 for (i=numpoints, coord=points+k; i--; coord +=
dim) {
1744 if (
qh->DELAUNAY && k ==
dim-1 && newhigh < newlow) {
1745 qh_fprintf(
qh,
qh->ferr, 6021,
"qhull input error: 'Qb%d' or 'QB%d' inverts paraboloid since high bound %.2g < low bound %.2g\n",
1746 k, k, newhigh, newlow);
1749 scale=
qh_divzero(newhigh - newlow, high - low,
1750 qh->MINdenom_1, &nearzero);
1752 qh_fprintf(
qh,
qh->ferr, 6022,
"qhull input error: %d'th dimension's new bounds [%2.2g, %2.2g] too wide for\nexisting bounds [%2.2g, %2.2g]\n",
1753 k, newlow, newhigh, low, high);
1756 shift= (newlow * high - low * newhigh)/(high-low);
1758 for (i=numpoints; i--; coord +=
dim)
1759 *coord= *coord * scale + shift;
1761 if (newlow < newhigh) {
1768 for (i=numpoints; i--; coord +=
dim) {
1772 trace0((
qh,
qh->ferr, 10,
"qh_scalepoints: scaled %d'th coordinate [%2.2g, %2.2g] to [%.2g, %.2g] for %d points by %2.2g and shifted %2.2g\n",
1773 k, low, high, newlow, newhigh, numpoints, scale, shift));
1813 trace0((
qh,
qh->ferr, 11,
"qh_setdelaunay: project %d points to paraboloid for Delaunay triangulation\n", count));
1815 for (i=0; i < count; i++) {
1817 paraboloid= coord*coord;
1818 for (k=
dim-2; k--; ) {
1820 paraboloid += coord*coord;
1822 *coordp++ = paraboloid;
1848 coordT *normp= normal, *feasiblep= feasible, *coordp= coords;
1856 dist += *(normp++) * *(feasiblep++);
1858 goto LABELerroroutside;
1860 if (dist < -qh->MINdenom) {
1862 *(coordp++)= *(normp++) / -dist;
1864 for (k=
dim; k--; ) {
1865 *(coordp++)=
qh_divzero(*(normp++), -dist,
qh->MINdenom_1, &zerodiv);
1867 goto LABELerroroutside;
1871 if (
qh->IStracing >= 4) {
1872 qh_fprintf(
qh,
qh->ferr, 8021,
"qh_sethalfspace: halfspace at offset %6.2g to point: ", *offset);
1873 for (k=
dim, coordp=coords; k--; ) {
1881 feasiblep= feasible;
1883 qh_fprintf(
qh,
qh->ferr, 6023,
"qhull input error: feasible point is not clearly inside halfspace\nfeasible point: ");
1923 coordT *coordp, *normalp, *offsetp;
1925 trace0((
qh,
qh->ferr, 12,
"qh_sethalfspace_all: compute dual for halfspace intersection\n"));
1928 qh_fprintf(
qh,
qh->ferr, 6024,
"qhull error: insufficient memory to compute dual of %d halfspaces\n",
1933 normalp= halfspaces;
1934 for (i=0; i < count; i++) {
1935 offsetp= normalp + newdim;
1936 if (!
qh_sethalfspace(
qh, newdim, coordp, &coordp, normalp, offsetp, feasible)) {
1938 qh_fprintf(
qh,
qh->ferr, 8032,
"The halfspace was at index %d\n", i);
1941 normalp= offsetp + 1;
1970 if (facet ==
qh->newfacet_list) {
1971 for (k=
qh->hull_dim; k--; )
1972 quadrant[ k]= (facet->
normal[ k] > 0);
1974 for (k=
qh->hull_dim; k--; ) {
1975 if (quadrant[ k] != (facet->
normal[ k] > 0)) {
1985 trace3((
qh,
qh->ferr, 3001,
"qh_sharpnewfacets: %d\n", issharp));
2015 pointT *point, **pointp, *point0;
2020 realT *diffp, sum2, *sum2row, *sum2p, det, factor;
2021 boolT nearzero, infinite;
2025 else if (size <
dim+1) {
2027 qh_fprintf(
qh,
qh->ferr, 6025,
"qhull internal error (qh_voronoi_center):\n need at least %d points to construct a Voronoi center\n",
2036 gmcoord=
qh->gm_matrix;
2037 for (k=0; k <
dim; k++) {
2038 qh->gm_row[k]= gmcoord;
2040 if (point != point0)
2041 *(gmcoord++)= point[k] - point0[k];
2045 for (i=0; i <
dim; i++) {
2047 for (k=0; k <
dim; k++) {
2048 diffp=
qh->gm_row[k] + i;
2049 sum2 += *diffp * *diffp;
2061 for (i=0; i <
dim; i++) {
2062 gmcoord=
qh->gm_matrix;
2064 for (k=0; k <
dim; k++) {
2065 qh->gm_row[k]= gmcoord;
2068 *(gmcoord++)= *sum2p++;
2071 if (point != point0)
2072 *(gmcoord++)= point[k] - point0[k];
2079 if (
qh->IStracing >= 3) {
2080 qh_fprintf(
qh,
qh->ferr, 8033,
"qh_voronoi_center: det %2.2g factor %2.2g ", det, factor);
2082 if (
qh->IStracing >= 5) {
2092 if (simplex != points)