1 /*<html><pre>  -<a                             href="qh-geom.htm"
2   >-------------------------------</a><a name="TOP">-</a>
3 
4 
5    geom2.c
6    infrequently used geometric routines of qhull
7 
8    see qh-geom.htm and geom.h
9 
10    Copyright (c) 1993-2020 The Geometry Center.
11    $Id: //main/2019/qhull/src/libqhull/geom2.c#14 $$Change: 3037 $
12    $DateTime: 2020/09/03 17:28:32 $$Author: bbarber $
13 
14    frequently used code goes into geom.c
15 */
16 
17 #include "qhull_a.h"
18 
19 /*================== functions in alphabetic order ============*/
20 
21 /*-<a                             href="qh-geom.htm#TOC"
22   >-------------------------------</a><a name="copypoints">-</a>
23 
24   qh_copypoints( points, numpoints, dimension )
25     return qh_malloc'd copy of points
26 
27   notes:
28     qh_free the returned points to avoid a memory leak
29 */
qh_copypoints(coordT * points,int numpoints,int dimension)30 coordT *qh_copypoints(coordT *points, int numpoints, int dimension)
31 {
32   int size;
33   coordT *newpoints;
34 
35   size= numpoints * dimension * (int)sizeof(coordT);
36   if (!(newpoints= (coordT *)qh_malloc((size_t)size))) {
37     qh_fprintf(qh ferr, 6004, "qhull error: insufficient memory to copy %d points\n",
38         numpoints);
39     qh_errexit(qh_ERRmem, NULL, NULL);
40   }
41   memcpy((char *)newpoints, (char *)points, (size_t)size); /* newpoints!=0 by QH6004 */
42   return newpoints;
43 } /* copypoints */
44 
45 /*-<a                             href="qh-geom.htm#TOC"
46   >-------------------------------</a><a name="crossproduct">-</a>
47 
48   qh_crossproduct( dim, vecA, vecB, vecC )
49     crossproduct of 2 dim vectors
50     C= A x B
51 
52   notes:
53     from Glasner, Graphics Gems I, p. 639
54     only defined for dim==3
55 */
qh_crossproduct(int dim,realT vecA[3],realT vecB[3],realT vecC[3])56 void qh_crossproduct(int dim, realT vecA[3], realT vecB[3], realT vecC[3]){
57 
58   if (dim == 3) {
59     vecC[0]=   det2_(vecA[1], vecA[2],
60                      vecB[1], vecB[2]);
61     vecC[1]= - det2_(vecA[0], vecA[2],
62                      vecB[0], vecB[2]);
63     vecC[2]=   det2_(vecA[0], vecA[1],
64                      vecB[0], vecB[1]);
65   }
66 } /* vcross */
67 
68 /*-<a                             href="qh-geom.htm#TOC"
69   >-------------------------------</a><a name="determinant">-</a>
70 
71   qh_determinant( rows, dim, nearzero )
72     compute signed determinant of a square matrix
73     uses qh.NEARzero to test for degenerate matrices
74 
75   returns:
76     determinant
77     overwrites rows and the matrix
78     if dim == 2 or 3
79       nearzero iff determinant < qh NEARzero[dim-1]
80       (!quite correct, not critical)
81     if dim >= 4
82       nearzero iff diagonal[k] < qh NEARzero[k]
83 */
qh_determinant(realT ** rows,int dim,boolT * nearzero)84 realT qh_determinant(realT **rows, int dim, boolT *nearzero) {
85   realT det=0;
86   int i;
87   boolT sign= False;
88 
89   *nearzero= False;
90   if (dim < 2) {
91     qh_fprintf(qh ferr, 6005, "qhull internal error (qh_determinate): only implemented for dimension >= 2\n");
92     qh_errexit(qh_ERRqhull, NULL, NULL);
93   }else if (dim == 2) {
94     det= det2_(rows[0][0], rows[0][1],
95                  rows[1][0], rows[1][1]);
96     if (fabs_(det) < 10*qh NEARzero[1])  /* QH11031 FIX: not really correct, what should this be? */
97       *nearzero= True;
98   }else if (dim == 3) {
99     det= det3_(rows[0][0], rows[0][1], rows[0][2],
100                  rows[1][0], rows[1][1], rows[1][2],
101                  rows[2][0], rows[2][1], rows[2][2]);
102     if (fabs_(det) < 10*qh NEARzero[2])  /* QH11031 FIX: what should this be?  det 5.5e-12 was flat for qh_maxsimplex of qdelaunay 0,0 27,27 -36,36 -9,63 */
103       *nearzero= True;
104   }else {
105     qh_gausselim(rows, dim, dim, &sign, nearzero);  /* if nearzero, diagonal still ok */
106     det= 1.0;
107     for (i=dim; i--; )
108       det *= (rows[i])[i];
109     if (sign)
110       det= -det;
111   }
112   return det;
113 } /* determinant */
114 
115 /*-<a                             href="qh-geom.htm#TOC"
116   >-------------------------------</a><a name="detjoggle">-</a>
117 
118   qh_detjoggle( points, numpoints, dimension )
119     determine default max joggle for point array
120       as qh_distround * qh_JOGGLEdefault
121 
122   returns:
123     initial value for JOGGLEmax from points and REALepsilon
124 
125   notes:
126     computes DISTround since qh_maxmin not called yet
127     if qh SCALElast, last dimension will be scaled later to MAXwidth
128 
129     loop duplicated from qh_maxmin
130 */
qh_detjoggle(pointT * points,int numpoints,int dimension)131 realT qh_detjoggle(pointT *points, int numpoints, int dimension) {
132   realT abscoord, distround, joggle, maxcoord, mincoord;
133   pointT *point, *pointtemp;
134   realT maxabs= -REALmax;
135   realT sumabs= 0;
136   realT maxwidth= 0;
137   int k;
138 
139   if (qh SETroundoff)
140     distround= qh DISTround; /* 'En' */
141   else{
142     for (k=0; k < dimension; k++) {
143       if (qh SCALElast && k == dimension-1)
144         abscoord= maxwidth;
145       else if (qh DELAUNAY && k == dimension-1) /* will qh_setdelaunay() */
146         abscoord= 2 * maxabs * maxabs;  /* may be low by qh hull_dim/2 */
147       else {
148         maxcoord= -REALmax;
149         mincoord= REALmax;
150         FORALLpoint_(points, numpoints) {
151           maximize_(maxcoord, point[k]);
152           minimize_(mincoord, point[k]);
153         }
154         maximize_(maxwidth, maxcoord-mincoord);
155         abscoord= fmax_(maxcoord, -mincoord);
156       }
157       sumabs += abscoord;
158       maximize_(maxabs, abscoord);
159     } /* for k */
160     distround= qh_distround(qh hull_dim, maxabs, sumabs);
161   }
162   joggle= distround * qh_JOGGLEdefault;
163   maximize_(joggle, REALepsilon * qh_JOGGLEdefault);
164   trace2((qh ferr, 2001, "qh_detjoggle: joggle=%2.2g maxwidth=%2.2g\n", joggle, maxwidth));
165   return joggle;
166 } /* detjoggle */
167 
168 /*-<a                             href="qh-geom.htm#TOC"
169   >-------------------------------</a><a name="detmaxoutside">-</a>
170 
171   qh_detmaxoutside();
172     determine qh.MAXoutside target for qh_RATIO... tests of distance
173     updates option '_max-outside'
174 
175   notes:
176     called from qh_addpoint and qh_detroundoff
177     accounts for qh.ONEmerge, qh.DISTround, qh.MINoutside ('Wn'), qh.max_outside
178     see qh_maxout for qh.max_outside with qh.DISTround
179 */
180 
qh_detmaxoutside(void)181 void qh_detmaxoutside(void) {
182   realT maxoutside;
183 
184   maxoutside= fmax_(qh max_outside, qh ONEmerge + qh DISTround);
185   maximize_(maxoutside, qh MINoutside);
186   qh MAXoutside= maxoutside;
187   trace3((qh ferr, 3056, "qh_detmaxoutside: MAXoutside %2.2g from qh.max_outside %2.2g, ONEmerge %2.2g, MINoutside %2.2g, DISTround %2.2g\n",
188       qh MAXoutside, qh max_outside, qh ONEmerge, qh MINoutside, qh DISTround));
189 } /* detmaxoutside */
190 
191 /*-<a                             href="qh-geom.htm#TOC"
192   >-------------------------------</a><a name="detroundoff">-</a>
193 
194   qh_detroundoff( )
195     determine maximum roundoff errors from
196       REALepsilon, REALmax, REALmin, qh.hull_dim, qh.MAXabs_coord,
197       qh.MAXsumcoord, qh.MAXwidth, qh.MINdenom_1
198 
199     accounts for qh.SETroundoff, qh.RANDOMdist, qh MERGEexact
200       qh.premerge_cos, qh.postmerge_cos, qh.premerge_centrum,
201       qh.postmerge_centrum, qh.MINoutside,
202       qh_RATIOnearinside, qh_COPLANARratio, qh_WIDEcoplanar
203 
204   returns:
205     sets qh.DISTround, etc. (see below)
206     appends precision constants to qh.qhull_options
207 
208   see:
209     qh_maxmin() for qh.NEARzero
210 
211   design:
212     determine qh.DISTround for distance computations
213     determine minimum denominators for qh_divzero
214     determine qh.ANGLEround for angle computations
215     adjust qh.premerge_cos,... for roundoff error
216     determine qh.ONEmerge for maximum error due to a single merge
217     determine qh.NEARinside, qh.MAXcoplanar, qh.MINvisible,
218       qh.MINoutside, qh.WIDEfacet
219     initialize qh.max_vertex and qh.minvertex
220 */
qh_detroundoff(void)221 void qh_detroundoff(void) {
222 
223   qh_option("_max-width", NULL, &qh MAXwidth);
224   if (!qh SETroundoff) {
225     qh DISTround= qh_distround(qh hull_dim, qh MAXabs_coord, qh MAXsumcoord);
226     qh_option("Error-roundoff", NULL, &qh DISTround);
227   }
228   qh MINdenom= qh MINdenom_1 * qh MAXabs_coord;
229   qh MINdenom_1_2= sqrt(qh MINdenom_1 * qh hull_dim) ;  /* if will be normalized */
230   qh MINdenom_2= qh MINdenom_1_2 * qh MAXabs_coord;
231                                               /* for inner product */
232   qh ANGLEround= 1.01 * qh hull_dim * REALepsilon;
233   if (qh RANDOMdist) {
234     qh ANGLEround += qh RANDOMfactor;
235     trace4((qh ferr, 4096, "qh_detroundoff: increase qh.ANGLEround by option 'R%2.2g'\n", qh RANDOMfactor));
236   }
237   if (qh premerge_cos < REALmax/2) {
238     qh premerge_cos -= qh ANGLEround;
239     if (qh RANDOMdist)
240       qh_option("Angle-premerge-with-random", NULL, &qh premerge_cos);
241   }
242   if (qh postmerge_cos < REALmax/2) {
243     qh postmerge_cos -= qh ANGLEround;
244     if (qh RANDOMdist)
245       qh_option("Angle-postmerge-with-random", NULL, &qh postmerge_cos);
246   }
247   qh premerge_centrum += 2 * qh DISTround;    /*2 for centrum and distplane()*/
248   qh postmerge_centrum += 2 * qh DISTround;
249   if (qh RANDOMdist && (qh MERGEexact || qh PREmerge))
250     qh_option("Centrum-premerge-with-random", NULL, &qh premerge_centrum);
251   if (qh RANDOMdist && qh POSTmerge)
252     qh_option("Centrum-postmerge-with-random", NULL, &qh postmerge_centrum);
253   { /* compute ONEmerge, max vertex offset for merging simplicial facets */
254     realT maxangle= 1.0, maxrho;
255 
256     minimize_(maxangle, qh premerge_cos);
257     minimize_(maxangle, qh postmerge_cos);
258     /* max diameter * sin theta + DISTround for vertex to its hyperplane */
259     qh ONEmerge= sqrt((realT)qh hull_dim) * qh MAXwidth *
260       sqrt(1.0 - maxangle * maxangle) + qh DISTround;
261     maxrho= qh hull_dim * qh premerge_centrum + qh DISTround;
262     maximize_(qh ONEmerge, maxrho);
263     maxrho= qh hull_dim * qh postmerge_centrum + qh DISTround;
264     maximize_(qh ONEmerge, maxrho);
265     if (qh MERGING)
266       qh_option("_one-merge", NULL, &qh ONEmerge);
267   }
268   qh NEARinside= qh ONEmerge * qh_RATIOnearinside; /* only used if qh KEEPnearinside */
269   if (qh JOGGLEmax < REALmax/2 && (qh KEEPcoplanar || qh KEEPinside)) {
270     realT maxdist;             /* adjust qh.NEARinside for joggle */
271     qh KEEPnearinside= True;
272     maxdist= sqrt((realT)qh hull_dim) * qh JOGGLEmax + qh DISTround;
273     maxdist= 2*maxdist;        /* vertex and coplanar point can joggle in opposite directions */
274     maximize_(qh NEARinside, maxdist);  /* must agree with qh_nearcoplanar() */
275   }
276   if (qh KEEPnearinside)
277     qh_option("_near-inside", NULL, &qh NEARinside);
278   if (qh JOGGLEmax < qh DISTround) {
279     qh_fprintf(qh ferr, 6006, "qhull option error: the joggle for 'QJn', %.2g, is below roundoff for distance computations, %.2g\n",
280          qh JOGGLEmax, qh DISTround);
281     qh_errexit(qh_ERRinput, NULL, NULL);
282   }
283   if (qh MINvisible > REALmax/2) {
284     if (!qh MERGING)
285       qh MINvisible= qh DISTround;
286     else if (qh hull_dim <= 3)
287       qh MINvisible= qh premerge_centrum;
288     else
289       qh MINvisible= qh_COPLANARratio * qh premerge_centrum;
290     if (qh APPROXhull && qh MINvisible > qh MINoutside)
291       qh MINvisible= qh MINoutside;
292     qh_option("Visible-distance", NULL, &qh MINvisible);
293   }
294   if (qh MAXcoplanar > REALmax/2) {
295     qh MAXcoplanar= qh MINvisible;
296     qh_option("U-max-coplanar", NULL, &qh MAXcoplanar);
297   }
298   if (!qh APPROXhull) {             /* user may specify qh MINoutside */
299     qh MINoutside= 2 * qh MINvisible;
300     if (qh premerge_cos < REALmax/2)
301       maximize_(qh MINoutside, (1- qh premerge_cos) * qh MAXabs_coord);
302     qh_option("Width-outside", NULL, &qh MINoutside);
303   }
304   qh WIDEfacet= qh MINoutside;
305   maximize_(qh WIDEfacet, qh_WIDEcoplanar * qh MAXcoplanar);
306   maximize_(qh WIDEfacet, qh_WIDEcoplanar * qh MINvisible);
307   qh_option("_wide-facet", NULL, &qh WIDEfacet);
308   if (qh MINvisible > qh MINoutside + 3 * REALepsilon
309   && !qh BESToutside && !qh FORCEoutput)
310     qh_fprintf(qh ferr, 7001, "qhull input warning: minimum visibility V%.2g is greater than \nminimum outside W%.2g.  Flipped facets are likely.\n",
311              qh MINvisible, qh MINoutside);
312   qh max_vertex= qh DISTround;
313   qh min_vertex= -qh DISTround;
314   /* numeric constants reported in printsummary */
315   qh_detmaxoutside();
316 } /* detroundoff */
317 
318 /*-<a                             href="qh-geom.htm#TOC"
319   >-------------------------------</a><a name="detsimplex">-</a>
320 
321   qh_detsimplex( apex, points, dim, nearzero )
322     compute determinant of a simplex with point apex and base points
323 
324   returns:
325      signed determinant and nearzero from qh_determinant
326 
327   notes:
328      called by qh_maxsimplex and qh_initialvertices
329      uses qh.gm_matrix/qh.gm_row (assumes they're big enough)
330 
331   design:
332     construct qm_matrix by subtracting apex from points
333     compute determinate
334 */
qh_detsimplex(pointT * apex,setT * points,int dim,boolT * nearzero)335 realT qh_detsimplex(pointT *apex, setT *points, int dim, boolT *nearzero) {
336   pointT *coorda, *coordp, *gmcoord, *point, **pointp;
337   coordT **rows;
338   int k,  i=0;
339   realT det;
340 
341   zinc_(Zdetsimplex);
342   gmcoord= qh gm_matrix;
343   rows= qh gm_row;
344   FOREACHpoint_(points) {
345     if (i == dim)
346       break;
347     rows[i++]= gmcoord;
348     coordp= point;
349     coorda= apex;
350     for (k=dim; k--; )
351       *(gmcoord++)= *coordp++ - *coorda++;
352   }
353   if (i < dim) {
354     qh_fprintf(qh ferr, 6007, "qhull internal error (qh_detsimplex): #points %d < dimension %d\n",
355                i, dim);
356     qh_errexit(qh_ERRqhull, NULL, NULL);
357   }
358   det= qh_determinant(rows, dim, nearzero);
359   trace2((qh ferr, 2002, "qh_detsimplex: det=%2.2g for point p%d, dim %d, nearzero? %d\n",
360           det, qh_pointid(apex), dim, *nearzero));
361   return det;
362 } /* detsimplex */
363 
364 /*-<a                             href="qh-geom.htm#TOC"
365   >-------------------------------</a><a name="distnorm">-</a>
366 
367   qh_distnorm( dim, point, normal, offset )
368     return distance from point to hyperplane at normal/offset
369 
370   returns:
371     dist
372 
373   notes:
374     dist > 0 if point is outside of hyperplane
375 
376   see:
377     qh_distplane in geom.c
378 */
qh_distnorm(int dim,pointT * point,pointT * normal,realT * offsetp)379 realT qh_distnorm(int dim, pointT *point, pointT *normal, realT *offsetp) {
380   coordT *normalp= normal, *coordp= point;
381   realT dist;
382   int k;
383 
384   dist= *offsetp;
385   for (k=dim; k--; )
386     dist += *(coordp++) * *(normalp++);
387   return dist;
388 } /* distnorm */
389 
390 /*-<a                             href="qh-geom.htm#TOC"
391   >-------------------------------</a><a name="distround">-</a>
392 
393   qh_distround( dimension, maxabs, maxsumabs )
394     compute maximum round-off error for a distance computation
395       to a normalized hyperplane
396     maxabs is the maximum absolute value of a coordinate
397     maxsumabs is the maximum possible sum of absolute coordinate values
398     if qh.RANDOMdist ('Qr'), adjusts qh_distround
399 
400   returns:
401     max dist round for qh.REALepsilon and qh.RANDOMdist
402 
403   notes:
404     calculate roundoff error according to Golub & van Loan, 1983, Lemma 3.2-1, "Rounding Errors"
405     use sqrt(dim) since one vector is normalized
406       or use maxsumabs since one vector is < 1
407 */
qh_distround(int dimension,realT maxabs,realT maxsumabs)408 realT qh_distround(int dimension, realT maxabs, realT maxsumabs) {
409   realT maxdistsum, maxround, delta;
410 
411   maxdistsum= sqrt((realT)dimension) * maxabs;
412   minimize_( maxdistsum, maxsumabs);
413   maxround= REALepsilon * (dimension * maxdistsum * 1.01 + maxabs);
414               /* adds maxabs for offset */
415   if (qh RANDOMdist) {
416     delta= qh RANDOMfactor * maxabs;
417     maxround += delta;
418     trace4((qh ferr, 4092, "qh_distround: increase roundoff by random delta %2.2g for option 'R%2.2g'\n", delta, qh RANDOMfactor));
419   }
420   trace4((qh ferr, 4008, "qh_distround: %2.2g, maxabs %2.2g, maxsumabs %2.2g, maxdistsum %2.2g\n",
421             maxround, maxabs, maxsumabs, maxdistsum));
422   return maxround;
423 } /* distround */
424 
425 /*-<a                             href="qh-geom.htm#TOC"
426   >-------------------------------</a><a name="divzero">-</a>
427 
428   qh_divzero( numer, denom, mindenom1, zerodiv )
429     divide by a number that's nearly zero
430     mindenom1= minimum denominator for dividing into 1.0
431 
432   returns:
433     quotient
434     sets zerodiv and returns 0.0 if it would overflow
435 
436   design:
437     if numer is nearly zero and abs(numer) < abs(denom)
438       return numer/denom
439     else if numer is nearly zero
440       return 0 and zerodiv
441     else if denom/numer non-zero
442       return numer/denom
443     else
444       return 0 and zerodiv
445 */
qh_divzero(realT numer,realT denom,realT mindenom1,boolT * zerodiv)446 realT qh_divzero(realT numer, realT denom, realT mindenom1, boolT *zerodiv) {
447   realT temp, numerx, denomx;
448 
449 
450   if (numer < mindenom1 && numer > -mindenom1) {
451     numerx= fabs_(numer);
452     denomx= fabs_(denom);
453     if (numerx < denomx) {
454       *zerodiv= False;
455       return numer/denom;
456     }else {
457       *zerodiv= True;
458       return 0.0;
459     }
460   }
461   temp= denom/numer;
462   if (temp > mindenom1 || temp < -mindenom1) {
463     *zerodiv= False;
464     return numer/denom;
465   }else {
466     *zerodiv= True;
467     return 0.0;
468   }
469 } /* divzero */
470 
471 
472 /*-<a                             href="qh-geom.htm#TOC"
473   >-------------------------------</a><a name="facetarea">-</a>
474 
475   qh_facetarea( facet )
476     return area for a facet
477 
478   notes:
479     if non-simplicial,
480       uses centrum to triangulate facet and sums the projected areas.
481     if (qh DELAUNAY),
482       computes projected area instead for last coordinate
483     assumes facet->normal exists
484     projecting tricoplanar facets to the hyperplane does not appear to make a difference
485 
486   design:
487     if simplicial
488       compute area
489     else
490       for each ridge
491         compute area from centrum to ridge
492     negate area if upper Delaunay facet
493 */
qh_facetarea(facetT * facet)494 realT qh_facetarea(facetT *facet) {
495   vertexT *apex;
496   pointT *centrum;
497   realT area= 0.0;
498   ridgeT *ridge, **ridgep;
499 
500   if (facet->simplicial) {
501     apex= SETfirstt_(facet->vertices, vertexT);
502     area= qh_facetarea_simplex(qh hull_dim, apex->point, facet->vertices,
503                     apex, facet->toporient, facet->normal, &facet->offset);
504   }else {
505     if (qh CENTERtype == qh_AScentrum)
506       centrum= facet->center;
507     else
508       centrum= qh_getcentrum(facet);
509     FOREACHridge_(facet->ridges)
510       area += qh_facetarea_simplex(qh hull_dim, centrum, ridge->vertices,
511                  NULL, (boolT)(ridge->top == facet),  facet->normal, &facet->offset);
512     if (qh CENTERtype != qh_AScentrum)
513       qh_memfree(centrum, qh normal_size);
514   }
515   if (facet->upperdelaunay && qh DELAUNAY)
516     area= -area;  /* the normal should be [0,...,1] */
517   trace4((qh ferr, 4009, "qh_facetarea: f%d area %2.2g\n", facet->id, area));
518   return area;
519 } /* facetarea */
520 
521 /*-<a                             href="qh-geom.htm#TOC"
522   >-------------------------------</a><a name="facetarea_simplex">-</a>
523 
524   qh_facetarea_simplex( dim, apex, vertices, notvertex, toporient, normal, offset )
525     return area for a simplex defined by
526       an apex, a base of vertices, an orientation, and a unit normal
527     if simplicial or tricoplanar facet,
528       notvertex is defined and it is skipped in vertices
529 
530   returns:
531     computes area of simplex projected to plane [normal,offset]
532     returns 0 if vertex too far below plane (qh WIDEfacet)
533       vertex can't be apex of tricoplanar facet
534 
535   notes:
536     if (qh DELAUNAY),
537       computes projected area instead for last coordinate
538     uses qh gm_matrix/gm_row and qh hull_dim
539     helper function for qh_facetarea
540 
541   design:
542     if Notvertex
543       translate simplex to apex
544     else
545       project simplex to normal/offset
546       translate simplex to apex
547     if Delaunay
548       set last row/column to 0 with -1 on diagonal
549     else
550       set last row to Normal
551     compute determinate
552     scale and flip sign for area
553 */
qh_facetarea_simplex(int dim,coordT * apex,setT * vertices,vertexT * notvertex,boolT toporient,coordT * normal,realT * offset)554 realT qh_facetarea_simplex(int dim, coordT *apex, setT *vertices,
555         vertexT *notvertex,  boolT toporient, coordT *normal, realT *offset) {
556   pointT *coorda, *coordp, *gmcoord;
557   coordT **rows, *normalp;
558   int k,  i=0;
559   realT area, dist;
560   vertexT *vertex, **vertexp;
561   boolT nearzero;
562 
563   gmcoord= qh gm_matrix;
564   rows= qh gm_row;
565   FOREACHvertex_(vertices) {
566     if (vertex == notvertex)
567       continue;
568     rows[i++]= gmcoord;
569     coorda= apex;
570     coordp= vertex->point;
571     normalp= normal;
572     if (notvertex) {
573       for (k=dim; k--; )
574         *(gmcoord++)= *coordp++ - *coorda++;
575     }else {
576       dist= *offset;
577       for (k=dim; k--; )
578         dist += *coordp++ * *normalp++;
579       if (dist < -qh WIDEfacet) {
580         zinc_(Znoarea);
581         return 0.0;
582       }
583       coordp= vertex->point;
584       normalp= normal;
585       for (k=dim; k--; )
586         *(gmcoord++)= (*coordp++ - dist * *normalp++) - *coorda++;
587     }
588   }
589   if (i != dim-1) {
590     qh_fprintf(qh ferr, 6008, "qhull internal error (qh_facetarea_simplex): #points %d != dim %d -1\n",
591                i, dim);
592     qh_errexit(qh_ERRqhull, NULL, NULL);
593   }
594   rows[i]= gmcoord;
595   if (qh DELAUNAY) {
596     for (i=0; i < dim-1; i++)
597       rows[i][dim-1]= 0.0;
598     for (k=dim; k--; )
599       *(gmcoord++)= 0.0;
600     rows[dim-1][dim-1]= -1.0;
601   }else {
602     normalp= normal;
603     for (k=dim; k--; )
604       *(gmcoord++)= *normalp++;
605   }
606   zinc_(Zdetfacetarea);
607   area= qh_determinant(rows, dim, &nearzero);
608   if (toporient)
609     area= -area;
610   area *= qh AREAfactor;
611   trace4((qh ferr, 4010, "qh_facetarea_simplex: area=%2.2g for point p%d, toporient %d, nearzero? %d\n",
612           area, qh_pointid(apex), toporient, nearzero));
613   return area;
614 } /* facetarea_simplex */
615 
616 /*-<a                             href="qh-geom.htm#TOC"
617   >-------------------------------</a><a name="facetcenter">-</a>
618 
619   qh_facetcenter( vertices )
620     return Voronoi center (Voronoi vertex) for a facet's vertices
621 
622   returns:
623     return temporary point equal to the center
624 
625   see:
626     qh_voronoi_center()
627 */
qh_facetcenter(setT * vertices)628 pointT *qh_facetcenter(setT *vertices) {
629   setT *points= qh_settemp(qh_setsize(vertices));
630   vertexT *vertex, **vertexp;
631   pointT *center;
632 
633   FOREACHvertex_(vertices)
634     qh_setappend(&points, vertex->point);
635   center= qh_voronoi_center(qh hull_dim-1, points);
636   qh_settempfree(&points);
637   return center;
638 } /* facetcenter */
639 
640 /*-<a                             href="qh-geom.htm#TOC"
641   >-------------------------------</a><a name="findgooddist">-</a>
642 
643   qh_findgooddist( point, facetA, dist, facetlist )
644     find best good facet visible for point from facetA
645     assumes facetA is visible from point
646 
647   returns:
648     best facet, i.e., good facet that is furthest from point
649       distance to best facet
650       NULL if none
651 
652     moves good, visible facets (and some other visible facets)
653       to end of qh facet_list
654 
655   notes:
656     uses qh visit_id
657 
658   design:
659     initialize bestfacet if facetA is good
660     move facetA to end of facetlist
661     for each facet on facetlist
662       for each unvisited neighbor of facet
663         move visible neighbors to end of facetlist
664         update best good neighbor
665         if no good neighbors, update best facet
666 */
qh_findgooddist(pointT * point,facetT * facetA,realT * distp,facetT ** facetlist)667 facetT *qh_findgooddist(pointT *point, facetT *facetA, realT *distp,
668                facetT **facetlist) {
669   realT bestdist= -REALmax, dist;
670   facetT *neighbor, **neighborp, *bestfacet=NULL, *facet;
671   boolT goodseen= False;
672 
673   if (facetA->good) {
674     zzinc_(Zcheckpart);  /* calls from check_bestdist occur after print stats */
675     qh_distplane(point, facetA, &bestdist);
676     bestfacet= facetA;
677     goodseen= True;
678   }
679   qh_removefacet(facetA);
680   qh_appendfacet(facetA);
681   *facetlist= facetA;
682   facetA->visitid= ++qh visit_id;
683   FORALLfacet_(*facetlist) {
684     FOREACHneighbor_(facet) {
685       if (neighbor->visitid == qh visit_id)
686         continue;
687       neighbor->visitid= qh visit_id;
688       if (goodseen && !neighbor->good)
689         continue;
690       zzinc_(Zcheckpart);
691       qh_distplane(point, neighbor, &dist);
692       if (dist > 0) {
693         qh_removefacet(neighbor);
694         qh_appendfacet(neighbor);
695         if (neighbor->good) {
696           goodseen= True;
697           if (dist > bestdist) {
698             bestdist= dist;
699             bestfacet= neighbor;
700           }
701         }
702       }
703     }
704   }
705   if (bestfacet) {
706     *distp= bestdist;
707     trace2((qh ferr, 2003, "qh_findgooddist: p%d is %2.2g above good facet f%d\n",
708       qh_pointid(point), bestdist, bestfacet->id));
709     return bestfacet;
710   }
711   trace4((qh ferr, 4011, "qh_findgooddist: no good facet for p%d above f%d\n",
712       qh_pointid(point), facetA->id));
713   return NULL;
714 }  /* findgooddist */
715 
716 /*-<a                             href="qh-geom.htm#TOC"
717   >-------------------------------</a><a name="furthestnewvertex">-</a>
718 
719   qh_furthestnewvertex( unvisited, facet, &maxdist )
720     return furthest unvisited, new vertex to a facet
721 
722   return:
723     NULL if no vertex is above facet
724     maxdist to facet
725     updates v.visitid
726 
727   notes:
728     Ignores vertices in facetB
729     Does not change qh.vertex_visit.  Use in conjunction with qh_furthestvertex
730 */
qh_furthestnewvertex(unsigned int unvisited,facetT * facet,realT * maxdistp)731 vertexT *qh_furthestnewvertex(unsigned int unvisited, facetT *facet, realT *maxdistp /* qh.newvertex_list */) {
732   vertexT *maxvertex= NULL, *vertex;
733   coordT dist, maxdist= 0.0;
734 
735   FORALLvertex_(qh newvertex_list) {
736     if (vertex->newfacet && vertex->visitid <= unvisited) {
737       vertex->visitid= qh vertex_visit;
738       qh_distplane(vertex->point, facet, &dist);
739       if (dist > maxdist) {
740         maxdist= dist;
741         maxvertex= vertex;
742       }
743     }
744   }
745   trace4((qh ferr, 4085, "qh_furthestnewvertex: v%d dist %2.2g is furthest new vertex for f%d\n",
746     getid_(maxvertex), maxdist, facet->id));
747   *maxdistp= maxdist;
748   return maxvertex;
749 } /* furthestnewvertex */
750 
751 /*-<a                             href="qh-geom.htm#TOC"
752   >-------------------------------</a><a name="furthestvertex">-</a>
753 
754   qh_furthestvertex( facetA, facetB, &maxdist, &mindist )
755     return furthest vertex in facetA from facetB, or NULL if none
756 
757   return:
758     maxdist and mindist to facetB or 0.0 if none
759     updates qh.vertex_visit
760 
761   notes:
762     Ignores vertices in facetB
763 */
qh_furthestvertex(facetT * facetA,facetT * facetB,realT * maxdistp,realT * mindistp)764 vertexT *qh_furthestvertex(facetT *facetA, facetT *facetB, realT *maxdistp, realT *mindistp) {
765   vertexT *maxvertex= NULL, *vertex, **vertexp;
766   coordT dist, maxdist= -REALmax, mindist= REALmax;
767 
768   qh vertex_visit++;
769   FOREACHvertex_(facetB->vertices)
770     vertex->visitid= qh vertex_visit;
771   FOREACHvertex_(facetA->vertices) {
772     if (vertex->visitid != qh vertex_visit) {
773       vertex->visitid= qh vertex_visit;
774       zzinc_(Zvertextests);
775       qh_distplane(vertex->point, facetB, &dist);
776       if (!maxvertex) {
777         maxdist= dist;
778         mindist= dist;
779         maxvertex= vertex;
780       }else if (dist > maxdist) {
781         maxdist= dist;
782         maxvertex= vertex;
783       }else if (dist < mindist)
784         mindist= dist;
785     }
786   }
787   if (!maxvertex) {
788     trace3((qh ferr, 3067, "qh_furthestvertex: all vertices of f%d are in f%d.  Returning 0.0 for max and mindist\n",
789       facetA->id, facetB->id));
790     maxdist= mindist= 0.0;
791   }else {
792     trace4((qh ferr, 4084, "qh_furthestvertex: v%d dist %2.2g is furthest (mindist %2.2g) of f%d above f%d\n",
793       maxvertex->id, maxdist, mindist, facetA->id, facetB->id));
794   }
795   *maxdistp= maxdist;
796   *mindistp= mindist;
797   return maxvertex;
798 } /* furthestvertex */
799 
800 /*-<a                             href="qh-geom.htm#TOC"
801   >-------------------------------</a><a name="getarea">-</a>
802 
803   qh_getarea( facetlist )
804     set area of all facets in facetlist
805     collect statistics
806     nop if hasAreaVolume
807 
808   returns:
809     sets qh totarea/totvol to total area and volume of convex hull
810     for Delaunay triangulation, computes projected area of the lower or upper hull
811       ignores upper hull if qh ATinfinity
812 
813   notes:
814     could compute outer volume by expanding facet area by rays from interior
815     the following attempt at perpendicular projection underestimated badly:
816       qh.totoutvol += (-dist + facet->maxoutside + qh DISTround)
817                             * area/ qh hull_dim;
818   design:
819     for each facet on facetlist
820       compute facet->area
821       update qh.totarea and qh.totvol
822 */
qh_getarea(facetT * facetlist)823 void qh_getarea(facetT *facetlist) {
824   realT area;
825   realT dist;
826   facetT *facet;
827 
828   if (qh hasAreaVolume)
829     return;
830   if (qh REPORTfreq)
831     qh_fprintf(qh ferr, 8020, "computing area of each facet and volume of the convex hull\n");
832   else
833     trace1((qh ferr, 1001, "qh_getarea: computing area for each facet and its volume to qh.interior_point (dist*area/dim)\n"));
834   qh totarea= qh totvol= 0.0;
835   FORALLfacet_(facetlist) {
836     if (!facet->normal)
837       continue;
838     if (facet->upperdelaunay && qh ATinfinity)
839       continue;
840     if (!facet->isarea) {
841       facet->f.area= qh_facetarea(facet);
842       facet->isarea= True;
843     }
844     area= facet->f.area;
845     if (qh DELAUNAY) {
846       if (facet->upperdelaunay == qh UPPERdelaunay)
847         qh totarea += area;
848     }else {
849       qh totarea += area;
850       qh_distplane(qh interior_point, facet, &dist);
851       qh totvol += -dist * area/ qh hull_dim;
852     }
853     if (qh PRINTstatistics) {
854       wadd_(Wareatot, area);
855       wmax_(Wareamax, area);
856       wmin_(Wareamin, area);
857     }
858   }
859   qh hasAreaVolume= True;
860 } /* getarea */
861 
862 /*-<a                             href="qh-geom.htm#TOC"
863   >-------------------------------</a><a name="gram_schmidt">-</a>
864 
865   qh_gram_schmidt( dim, row )
866     implements Gram-Schmidt orthogonalization by rows
867 
868   returns:
869     false if zero norm
870     overwrites rows[dim][dim]
871 
872   notes:
873     see Golub & van Loan, 1983, Algorithm 6.2-2, "Modified Gram-Schmidt"
874     overflow due to small divisors not handled
875 
876   design:
877     for each row
878       compute norm for row
879       if non-zero, normalize row
880       for each remaining rowA
881         compute inner product of row and rowA
882         reduce rowA by row * inner product
883 */
qh_gram_schmidt(int dim,realT ** row)884 boolT qh_gram_schmidt(int dim, realT **row) {
885   realT *rowi, *rowj, norm;
886   int i, j, k;
887 
888   for (i=0; i < dim; i++) {
889     rowi= row[i];
890     for (norm=0.0, k=dim; k--; rowi++)
891       norm += *rowi * *rowi;
892     norm= sqrt(norm);
893     wmin_(Wmindenom, norm);
894     if (norm == 0.0)  /* either 0 or overflow due to sqrt */
895       return False;
896     for (k=dim; k--; )
897       *(--rowi) /= norm;
898     for (j=i+1; j < dim; j++) {
899       rowj= row[j];
900       for (norm=0.0, k=dim; k--; )
901         norm += *rowi++ * *rowj++;
902       for (k=dim; k--; )
903         *(--rowj) -= *(--rowi) * norm;
904     }
905   }
906   return True;
907 } /* gram_schmidt */
908 
909 
910 /*-<a                             href="qh-geom.htm#TOC"
911   >-------------------------------</a><a name="inthresholds">-</a>
912 
913   qh_inthresholds( normal, angle )
914     return True if normal within qh.lower_/upper_threshold
915 
916   returns:
917     estimate of angle by summing of threshold diffs
918       angle may be NULL
919       smaller "angle" is better
920 
921   notes:
922     invalid if qh.SPLITthresholds
923 
924   see:
925     qh.lower_threshold in qh_initbuild()
926     qh_initthresholds()
927 
928   design:
929     for each dimension
930       test threshold
931 */
qh_inthresholds(coordT * normal,realT * angle)932 boolT qh_inthresholds(coordT *normal, realT *angle) {
933   boolT within= True;
934   int k;
935   realT threshold;
936 
937   if (angle)
938     *angle= 0.0;
939   for (k=0; k < qh hull_dim; k++) {
940     threshold= qh lower_threshold[k];
941     if (threshold > -REALmax/2) {
942       if (normal[k] < threshold)
943         within= False;
944       if (angle) {
945         threshold -= normal[k];
946         *angle += fabs_(threshold);
947       }
948     }
949     if (qh upper_threshold[k] < REALmax/2) {
950       threshold= qh upper_threshold[k];
951       if (normal[k] > threshold)
952         within= False;
953       if (angle) {
954         threshold -= normal[k];
955         *angle += fabs_(threshold);
956       }
957     }
958   }
959   return within;
960 } /* inthresholds */
961 
962 
963 /*-<a                             href="qh-geom.htm#TOC"
964   >-------------------------------</a><a name="joggleinput">-</a>
965 
966   qh_joggleinput( )
967     randomly joggle input to Qhull by qh.JOGGLEmax
968     initial input is qh.first_point/qh.num_points of qh.hull_dim
969       repeated calls use qh.input_points/qh.num_points
970 
971   returns:
972     joggles points at qh.first_point/qh.num_points
973     copies data to qh.input_points/qh.input_malloc if first time
974     determines qh.JOGGLEmax if it was zero
975     if qh.DELAUNAY
976       computes the Delaunay projection of the joggled points
977 
978   notes:
979     if qh.DELAUNAY, unnecessarily joggles the last coordinate
980     the initial 'QJn' may be set larger than qh_JOGGLEmaxincrease
981 
982   design:
983     if qh.DELAUNAY
984       set qh.SCALElast for reduced precision errors
985     if first call
986       initialize qh.input_points to the original input points
987       if qh.JOGGLEmax == 0
988         determine default qh.JOGGLEmax
989     else
990       increase qh.JOGGLEmax according to qh.build_cnt
991     joggle the input by adding a random number in [-qh.JOGGLEmax,qh.JOGGLEmax]
992     if qh.DELAUNAY
993       sets the Delaunay projection
994 */
qh_joggleinput(void)995 void qh_joggleinput(void) {
996   int i, seed, size;
997   coordT *coordp, *inputp;
998   realT randr, randa, randb;
999 
1000   if (!qh input_points) { /* first call */
1001     qh input_points= qh first_point;
1002     qh input_malloc= qh POINTSmalloc;
1003     size= qh num_points * qh hull_dim * (int)sizeof(coordT);
1004     if (!(qh first_point= (coordT *)qh_malloc((size_t)size))) {
1005       qh_fprintf(qh ferr, 6009, "qhull error: insufficient memory to joggle %d points\n",
1006           qh num_points);
1007       qh_errexit(qh_ERRmem, NULL, NULL);
1008     }
1009     qh POINTSmalloc= True;
1010     if (qh JOGGLEmax == 0.0) {
1011       qh JOGGLEmax= qh_detjoggle(qh input_points, qh num_points, qh hull_dim);
1012       qh_option("QJoggle", NULL, &qh JOGGLEmax);
1013     }
1014   }else {                 /* repeated call */
1015     if (!qh RERUN && qh build_cnt > qh_JOGGLEretry) {
1016       if (((qh build_cnt-qh_JOGGLEretry-1) % qh_JOGGLEagain) == 0) {
1017         realT maxjoggle= qh MAXwidth * qh_JOGGLEmaxincrease;
1018         if (qh JOGGLEmax < maxjoggle) {
1019           qh JOGGLEmax *= qh_JOGGLEincrease;
1020           minimize_(qh JOGGLEmax, maxjoggle);
1021         }
1022       }
1023     }
1024     qh_option("QJoggle", NULL, &qh JOGGLEmax);
1025   }
1026   if (qh build_cnt > 1 && qh JOGGLEmax > fmax_(qh MAXwidth/4, 0.1)) {
1027       qh_fprintf(qh ferr, 6010, "qhull input error (qh_joggleinput): the current joggle for 'QJn', %.2g, is too large for the width\nof the input.  If possible, recompile Qhull with higher-precision reals.\n",
1028                 qh JOGGLEmax);
1029       qh_errexit(qh_ERRinput, NULL, NULL);
1030   }
1031   /* for some reason, using qh ROTATErandom and qh_RANDOMseed does not repeat the run. Use 'TRn' instead */
1032   seed= qh_RANDOMint;
1033   qh_option("_joggle-seed", &seed, NULL);
1034   trace0((qh ferr, 6, "qh_joggleinput: joggle input by %4.4g with seed %d\n",
1035     qh JOGGLEmax, seed));
1036   inputp= qh input_points;
1037   coordp= qh first_point;
1038   randa= 2.0 * qh JOGGLEmax/qh_RANDOMmax;
1039   randb= -qh JOGGLEmax;
1040   size= qh num_points * qh hull_dim;
1041   for (i=size; i--; ) {
1042     randr= qh_RANDOMint;
1043     *(coordp++)= *(inputp++) + (randr * randa + randb);
1044   }
1045   if (qh DELAUNAY) {
1046     qh last_low= qh last_high= qh last_newhigh= REALmax;
1047     qh_setdelaunay(qh hull_dim, qh num_points, qh first_point);
1048   }
1049 } /* joggleinput */
1050 
1051 /*-<a                             href="qh-geom.htm#TOC"
1052   >-------------------------------</a><a name="maxabsval">-</a>
1053 
1054   qh_maxabsval( normal, dim )
1055     return pointer to maximum absolute value of a dim vector
1056     returns NULL if dim=0
1057 */
qh_maxabsval(realT * normal,int dim)1058 realT *qh_maxabsval(realT *normal, int dim) {
1059   realT maxval= -REALmax;
1060   realT *maxp= NULL, *colp, absval;
1061   int k;
1062 
1063   for (k=dim, colp= normal; k--; colp++) {
1064     absval= fabs_(*colp);
1065     if (absval > maxval) {
1066       maxval= absval;
1067       maxp= colp;
1068     }
1069   }
1070   return maxp;
1071 } /* maxabsval */
1072 
1073 
1074 /*-<a                             href="qh-geom.htm#TOC"
1075   >-------------------------------</a><a name="maxmin">-</a>
1076 
1077   qh_maxmin( points, numpoints, dimension )
1078     return max/min points for each dimension
1079     determine max and min coordinates
1080 
1081   returns:
1082     returns a temporary set of max and min points
1083       may include duplicate points. Does not include qh.GOODpoint
1084     sets qh.NEARzero, qh.MAXabs_coord, qh.MAXsumcoord, qh.MAXwidth
1085          qh.MAXlastcoord, qh.MINlastcoord
1086     initializes qh.max_outside, qh.min_vertex, qh.WAScoplanar, qh.ZEROall_ok
1087 
1088   notes:
1089     loop duplicated in qh_detjoggle()
1090 
1091   design:
1092     initialize global precision variables
1093     checks definition of REAL...
1094     for each dimension
1095       for each point
1096         collect maximum and minimum point
1097       collect maximum of maximums and minimum of minimums
1098       determine qh.NEARzero for Gaussian Elimination
1099 */
qh_maxmin(pointT * points,int numpoints,int dimension)1100 setT *qh_maxmin(pointT *points, int numpoints, int dimension) {
1101   int k;
1102   realT maxcoord, temp;
1103   pointT *minimum, *maximum, *point, *pointtemp;
1104   setT *set;
1105 
1106   qh max_outside= 0.0;
1107   qh MAXabs_coord= 0.0;
1108   qh MAXwidth= -REALmax;
1109   qh MAXsumcoord= 0.0;
1110   qh min_vertex= 0.0;
1111   qh WAScoplanar= False;
1112   if (qh ZEROcentrum)
1113     qh ZEROall_ok= True;
1114   if (REALmin < REALepsilon && REALmin < REALmax && REALmin > -REALmax
1115   && REALmax > 0.0 && -REALmax < 0.0)
1116     ; /* all ok */
1117   else {
1118     qh_fprintf(qh ferr, 6011, "qhull error: one or more floating point constants in user.h are inconsistent. REALmin %g, -REALmax %g, 0.0, REALepsilon %g, REALmax %g\n",
1119           REALmin, -REALmax, REALepsilon, REALmax);
1120     qh_errexit(qh_ERRinput, NULL, NULL);
1121   }
1122   set= qh_settemp(2*dimension);
1123   trace1((qh ferr, 8082, "qh_maxmin: dim             min             max           width    nearzero  min-point  max-point\n"));
1124   for (k=0; k < dimension; k++) {
1125     if (points == qh GOODpointp)
1126       minimum= maximum= points + dimension;
1127     else
1128       minimum= maximum= points;
1129     FORALLpoint_(points, numpoints) {
1130       if (point == qh GOODpointp)
1131         continue;
1132       if (maximum[k] < point[k])
1133         maximum= point;
1134       else if (minimum[k] > point[k])
1135         minimum= point;
1136     }
1137     if (k == dimension-1) {
1138       qh MINlastcoord= minimum[k];
1139       qh MAXlastcoord= maximum[k];
1140     }
1141     if (qh SCALElast && k == dimension-1)
1142       maxcoord= qh MAXabs_coord;
1143     else {
1144       maxcoord= fmax_(maximum[k], -minimum[k]);
1145       if (qh GOODpointp) {
1146         temp= fmax_(qh GOODpointp[k], -qh GOODpointp[k]);
1147         maximize_(maxcoord, temp);
1148       }
1149       temp= maximum[k] - minimum[k];
1150       maximize_(qh MAXwidth, temp);
1151     }
1152     maximize_(qh MAXabs_coord, maxcoord);
1153     qh MAXsumcoord += maxcoord;
1154     qh_setappend(&set, minimum);
1155     qh_setappend(&set, maximum);
1156     /* calculation of qh NEARzero is based on Golub & van Loan, 1983,
1157        Eq. 4.4-13 for "Gaussian elimination with complete pivoting".
1158        Golub & van Loan say that n^3 can be ignored and 10 be used in
1159        place of rho */
1160     qh NEARzero[k]= 80 * qh MAXsumcoord * REALepsilon;
1161     trace1((qh ferr, 8106, "           %3d % 14.8e % 14.8e % 14.8e  %4.4e  p%-9d p%-d\n",
1162             k, minimum[k], maximum[k], maximum[k]-minimum[k], qh NEARzero[k], qh_pointid(minimum), qh_pointid(maximum)));
1163     if (qh SCALElast && k == dimension-1)
1164       trace1((qh ferr, 8107, "           last coordinate scaled to (%4.4g, %4.4g), width %4.4e for option 'Qbb'\n",
1165             qh MAXabs_coord - qh MAXwidth, qh MAXabs_coord, qh MAXwidth));
1166   }
1167   if (qh IStracing >= 1)
1168     qh_printpoints(qh ferr, "qh_maxmin: found the max and min points (by dim):", set);
1169   return(set);
1170 } /* maxmin */
1171 
1172 /*-<a                             href="qh-geom.htm#TOC"
1173   >-------------------------------</a><a name="maxouter">-</a>
1174 
1175   qh_maxouter( )
1176     return maximum distance from facet to outer plane
1177     normally this is qh.max_outside+qh.DISTround
1178     does not include qh.JOGGLEmax
1179 
1180   see:
1181     qh_outerinner()
1182 
1183   notes:
1184     need to add another qh.DISTround if testing actual point with computation
1185     see qh_detmaxoutside for a qh_RATIO... target
1186 
1187   for joggle:
1188     qh_setfacetplane() updated qh.max_outer for Wnewvertexmax (max distance to vertex)
1189     need to use Wnewvertexmax since could have a coplanar point for a high
1190       facet that is replaced by a low facet
1191     need to add qh.JOGGLEmax if testing input points
1192 */
qh_maxouter(void)1193 realT qh_maxouter(void) {
1194   realT dist;
1195 
1196   dist= fmax_(qh max_outside, qh DISTround);
1197   dist += qh DISTround;
1198   trace4((qh ferr, 4012, "qh_maxouter: max distance from facet to outer plane is %4.4g, qh.max_outside is %4.4g\n", dist, qh max_outside));
1199   return dist;
1200 } /* maxouter */
1201 
1202 /*-<a                             href="qh-geom.htm#TOC"
1203   >-------------------------------</a><a name="maxsimplex">-</a>
1204 
1205   qh_maxsimplex( dim, maxpoints, points, numpoints, simplex )
1206     determines maximum simplex for a set of points
1207     maxpoints is the subset of points with a min or max coordinate
1208     may start with points already in simplex
1209     skips qh.GOODpointp (assumes that it isn't in maxpoints)
1210 
1211   returns:
1212     simplex with dim+1 points
1213 
1214   notes:
1215     called by qh_initialvertices, qh_detvnorm, and qh_voronoi_center
1216     requires qh.MAXwidth to estimate determinate for each vertex
1217     assumes at least needed points in points
1218     maximizes determinate for x,y,z,w, etc.
1219     uses maxpoints as long as determinate is clearly non-zero
1220 
1221   design:
1222     initialize simplex with at least two points
1223       (find points with max or min x coordinate)
1224     create a simplex of dim+1 vertices as follows
1225       add point from maxpoints that maximizes the determinate of the point and the simplex vertices
1226       if last point and maxdet/prevdet < qh_RATIOmaxsimplex (3.0e-2)
1227         flag maybe_falsenarrow
1228       if no maxpoint or maxnearzero or maybe_falsenarrow
1229         search all points for maximum determinate
1230         early exit if maybe_falsenarrow and !maxnearzero and maxdet > prevdet
1231 */
qh_maxsimplex(int dim,setT * maxpoints,pointT * points,int numpoints,setT ** simplex)1232 void qh_maxsimplex(int dim, setT *maxpoints, pointT *points, int numpoints, setT **simplex) {
1233   pointT *point, **pointp, *pointtemp, *maxpoint, *minx=NULL, *maxx=NULL;
1234   boolT nearzero, maxnearzero= False, maybe_falsenarrow;
1235   int i, sizinit;
1236   realT maxdet= -1.0, prevdet= -1.0, det, mincoord= REALmax, maxcoord= -REALmax, mindet, ratio, targetdet;
1237 
1238   if (qh MAXwidth <= 0.0) {
1239     qh_fprintf(qh ferr, 6421, "qhull internal error (qh_maxsimplex): qh.MAXwidth required for qh_maxsimplex.  Used to estimate determinate\n");
1240     qh_errexit(qh_ERRqhull, NULL, NULL);
1241   }
1242   sizinit= qh_setsize(*simplex);
1243   if (sizinit >= 2) {
1244     maxdet= pow(qh MAXwidth, sizinit - 1);
1245   }else {
1246     if (qh_setsize(maxpoints) >= 2) {
1247       FOREACHpoint_(maxpoints) {
1248         if (maxcoord < point[0]) {
1249           maxcoord= point[0];
1250           maxx= point;
1251         }
1252         if (mincoord > point[0]) {
1253           mincoord= point[0];
1254           minx= point;
1255         }
1256       }
1257     }else {
1258       FORALLpoint_(points, numpoints) {
1259         if (point == qh GOODpointp)
1260           continue;
1261         if (maxcoord < point[0]) {
1262           maxcoord= point[0];
1263           maxx= point;
1264         }
1265         if (mincoord > point[0]) {
1266           mincoord= point[0];
1267           minx= point;
1268         }
1269       }
1270     }
1271     maxdet= maxcoord - mincoord;
1272     qh_setunique(simplex, minx);
1273     if (qh_setsize(*simplex) < 2)
1274       qh_setunique(simplex, maxx);
1275     sizinit= qh_setsize(*simplex);
1276     if (sizinit < 2) {
1277       qh_joggle_restart("input has same x coordinate");
1278       if (zzval_(Zsetplane) > qh hull_dim+1) {
1279         qh_fprintf(qh ferr, 6012, "qhull precision error (qh_maxsimplex for voronoi_center): %d points with the same x coordinate %4.4g\n",
1280                  qh_setsize(maxpoints)+numpoints, mincoord);
1281         qh_errexit(qh_ERRprec, NULL, NULL);
1282       }else {
1283         qh_fprintf(qh ferr, 6013, "qhull input error: input is less than %d-dimensional since all points have the same x coordinate %4.4g\n",
1284                  qh hull_dim, mincoord);
1285         qh_errexit(qh_ERRinput, NULL, NULL);
1286       }
1287     }
1288   }
1289   for (i=sizinit; i < dim+1; i++) {
1290     prevdet= maxdet;
1291     maxpoint= NULL;
1292     maxdet= -1.0;
1293     FOREACHpoint_(maxpoints) {
1294       if (!qh_setin(*simplex, point) && point != maxpoint) {
1295         det= qh_detsimplex(point, *simplex, i, &nearzero); /* retests maxpoints if duplicate or multiple iterations */
1296         if ((det= fabs_(det)) > maxdet) {
1297           maxdet= det;
1298           maxpoint= point;
1299           maxnearzero= nearzero;
1300         }
1301       }
1302     }
1303     maybe_falsenarrow= False;
1304     ratio= 1.0;
1305     targetdet= prevdet * qh MAXwidth;
1306     mindet= 10 * qh_RATIOmaxsimplex * targetdet;
1307     if (maxdet > 0.0) {
1308       ratio= maxdet / targetdet;
1309       if (ratio < qh_RATIOmaxsimplex)
1310         maybe_falsenarrow= True;
1311     }
1312     if (!maxpoint || maxnearzero || maybe_falsenarrow) {
1313       zinc_(Zsearchpoints);
1314       if (!maxpoint) {
1315         trace0((qh ferr, 7, "qh_maxsimplex: searching all points for %d-th initial vertex, better than mindet %4.4g, targetdet %4.4g\n",
1316                 i+1, mindet, targetdet));
1317       }else if (qh ALLpoints) {
1318         trace0((qh ferr, 30, "qh_maxsimplex: searching all points ('Qs') for %d-th initial vertex, better than p%d det %4.4g, targetdet %4.4g, ratio %4.4g\n",
1319                 i+1, qh_pointid(maxpoint), maxdet, targetdet, ratio));
1320       }else if (maybe_falsenarrow) {
1321         trace0((qh ferr, 17, "qh_maxsimplex: searching all points for %d-th initial vertex, better than p%d det %4.4g and mindet %4.4g, ratio %4.4g\n",
1322                 i+1, qh_pointid(maxpoint), maxdet, mindet, ratio));
1323       }else {
1324         trace0((qh ferr, 8, "qh_maxsimplex: searching all points for %d-th initial vertex, better than p%d det %2.2g and mindet %4.4g, targetdet %4.4g\n",
1325                 i+1, qh_pointid(maxpoint), maxdet, mindet, targetdet));
1326       }
1327       FORALLpoint_(points, numpoints) {
1328         if (point == qh GOODpointp)
1329           continue;
1330         if (!qh_setin(maxpoints, point) && !qh_setin(*simplex, point)) {
1331           det= qh_detsimplex(point, *simplex, i, &nearzero);
1332           if ((det= fabs_(det)) > maxdet) {
1333             maxdet= det;
1334             maxpoint= point;
1335             maxnearzero= nearzero;
1336             if (!maxnearzero && !qh ALLpoints && maxdet > mindet)
1337               break;
1338           }
1339         }
1340       }
1341     } /* !maxpoint */
1342     if (!maxpoint) {
1343       qh_fprintf(qh ferr, 6014, "qhull internal error (qh_maxsimplex): not enough points available\n");
1344       qh_errexit(qh_ERRqhull, NULL, NULL);
1345     }
1346     qh_setappend(simplex, maxpoint);
1347     trace1((qh ferr, 1002, "qh_maxsimplex: selected point p%d for %d`th initial vertex, det=%4.4g, targetdet=%4.4g, mindet=%4.4g\n",
1348             qh_pointid(maxpoint), i+1, maxdet, prevdet * qh MAXwidth, mindet));
1349   } /* i */
1350 } /* maxsimplex */
1351 
1352 /*-<a                             href="qh-geom.htm#TOC"
1353   >-------------------------------</a><a name="minabsval">-</a>
1354 
1355   qh_minabsval( normal, dim )
1356     return minimum absolute value of a dim vector
1357 */
qh_minabsval(realT * normal,int dim)1358 realT qh_minabsval(realT *normal, int dim) {
1359   realT minval= 0;
1360   realT maxval= 0;
1361   realT *colp;
1362   int k;
1363 
1364   for (k=dim, colp=normal; k--; colp++) {
1365     maximize_(maxval, *colp);
1366     minimize_(minval, *colp);
1367   }
1368   return fmax_(maxval, -minval);
1369 } /* minabsval */
1370 
1371 
1372 /*-<a                             href="qh-geom.htm#TOC"
1373   >-------------------------------</a><a name="mindiff">-</a>
1374 
1375   qh_mindif( vecA, vecB, dim )
1376     return index of min abs. difference of two vectors
1377 */
qh_mindiff(realT * vecA,realT * vecB,int dim)1378 int qh_mindiff(realT *vecA, realT *vecB, int dim) {
1379   realT mindiff= REALmax, diff;
1380   realT *vecAp= vecA, *vecBp= vecB;
1381   int k, mink= 0;
1382 
1383   for (k=0; k < dim; k++) {
1384     diff= *vecAp++ - *vecBp++;
1385     diff= fabs_(diff);
1386     if (diff < mindiff) {
1387       mindiff= diff;
1388       mink= k;
1389     }
1390   }
1391   return mink;
1392 } /* mindiff */
1393 
1394 
1395 
1396 /*-<a                             href="qh-geom.htm#TOC"
1397   >-------------------------------</a><a name="orientoutside">-</a>
1398 
1399   qh_orientoutside( facet  )
1400     make facet outside oriented via qh.interior_point
1401 
1402   returns:
1403     True if facet reversed orientation.
1404 */
qh_orientoutside(facetT * facet)1405 boolT qh_orientoutside(facetT *facet) {
1406   int k;
1407   realT dist;
1408 
1409   qh_distplane(qh interior_point, facet, &dist);
1410   if (dist > 0) {
1411     for (k=qh hull_dim; k--; )
1412       facet->normal[k]= -facet->normal[k];
1413     facet->offset= -facet->offset;
1414     return True;
1415   }
1416   return False;
1417 } /* orientoutside */
1418 
1419 /*-<a                             href="qh-geom.htm#TOC"
1420   >-------------------------------</a><a name="outerinner">-</a>
1421 
1422   qh_outerinner( facet, outerplane, innerplane  )
1423     if facet and qh.maxoutdone (i.e., qh_check_maxout)
1424       returns outer and inner plane for facet
1425     else
1426       returns maximum outer and inner plane
1427     accounts for qh.JOGGLEmax
1428 
1429   see:
1430     qh_maxouter(), qh_check_bestdist(), qh_check_points()
1431 
1432   notes:
1433     outerplaner or innerplane may be NULL
1434     facet is const
1435     Does not error (QhullFacet)
1436 
1437     includes qh.DISTround for actual points
1438     adds another qh.DISTround if testing with floating point arithmetic
1439 */
qh_outerinner(facetT * facet,realT * outerplane,realT * innerplane)1440 void qh_outerinner(facetT *facet, realT *outerplane, realT *innerplane) {
1441   realT dist, mindist;
1442   vertexT *vertex, **vertexp;
1443 
1444   if (outerplane) {
1445     if (!qh_MAXoutside || !facet || !qh maxoutdone) {
1446       *outerplane= qh_maxouter();       /* includes qh.DISTround */
1447     }else { /* qh_MAXoutside ... */
1448 #if qh_MAXoutside
1449       *outerplane= facet->maxoutside + qh DISTround;
1450 #endif
1451 
1452     }
1453     if (qh JOGGLEmax < REALmax/2)
1454       *outerplane += qh JOGGLEmax * sqrt((realT)qh hull_dim);
1455   }
1456   if (innerplane) {
1457     if (facet) {
1458       mindist= REALmax;
1459       FOREACHvertex_(facet->vertices) {
1460         zinc_(Zdistio);
1461         qh_distplane(vertex->point, facet, &dist);
1462         minimize_(mindist, dist);
1463       }
1464       *innerplane= mindist - qh DISTround;
1465     }else
1466       *innerplane= qh min_vertex - qh DISTround;
1467     if (qh JOGGLEmax < REALmax/2)
1468       *innerplane -= qh JOGGLEmax * sqrt((realT)qh hull_dim);
1469   }
1470 } /* outerinner */
1471 
1472 /*-<a                             href="qh-geom.htm#TOC"
1473   >-------------------------------</a><a name="pointdist">-</a>
1474 
1475   qh_pointdist( point1, point2, dim )
1476     return distance between two points
1477 
1478   notes:
1479     returns distance squared if 'dim' is negative
1480 */
qh_pointdist(pointT * point1,pointT * point2,int dim)1481 coordT qh_pointdist(pointT *point1, pointT *point2, int dim) {
1482   coordT dist, diff;
1483   int k;
1484 
1485   dist= 0.0;
1486   for (k= (dim > 0 ? dim : -dim); k--; ) {
1487     diff= *point1++ - *point2++;
1488     dist += diff * diff;
1489   }
1490   if (dim > 0)
1491     return(sqrt(dist));
1492   return dist;
1493 } /* pointdist */
1494 
1495 
1496 /*-<a                             href="qh-geom.htm#TOC"
1497   >-------------------------------</a><a name="printmatrix">-</a>
1498 
1499   qh_printmatrix( fp, string, rows, numrow, numcol )
1500     print matrix to fp given by row vectors
1501     print string as header
1502 
1503   notes:
1504     print a vector by qh_printmatrix(fp, "", &vect, 1, len)
1505 */
qh_printmatrix(FILE * fp,const char * string,realT ** rows,int numrow,int numcol)1506 void qh_printmatrix(FILE *fp, const char *string, realT **rows, int numrow, int numcol) {
1507   realT *rowp;
1508   realT r; /*bug fix*/
1509   int i,k;
1510 
1511   qh_fprintf(fp, 9001, "%s\n", string);
1512   for (i=0; i < numrow; i++) {
1513     rowp= rows[i];
1514     for (k=0; k < numcol; k++) {
1515       r= *rowp++;
1516       qh_fprintf(fp, 9002, "%6.3g ", r);
1517     }
1518     qh_fprintf(fp, 9003, "\n");
1519   }
1520 } /* printmatrix */
1521 
1522 
1523 /*-<a                             href="qh-geom.htm#TOC"
1524   >-------------------------------</a><a name="printpoints">-</a>
1525 
1526   qh_printpoints( fp, string, points )
1527     print pointids to fp for a set of points
1528     if string, prints string and 'p' point ids
1529 */
qh_printpoints(FILE * fp,const char * string,setT * points)1530 void qh_printpoints(FILE *fp, const char *string, setT *points) {
1531   pointT *point, **pointp;
1532 
1533   if (string) {
1534     qh_fprintf(fp, 9004, "%s", string);
1535     FOREACHpoint_(points)
1536       qh_fprintf(fp, 9005, " p%d", qh_pointid(point));
1537     qh_fprintf(fp, 9006, "\n");
1538   }else {
1539     FOREACHpoint_(points)
1540       qh_fprintf(fp, 9007, " %d", qh_pointid(point));
1541     qh_fprintf(fp, 9008, "\n");
1542   }
1543 } /* printpoints */
1544 
1545 
1546 /*-<a                             href="qh-geom.htm#TOC"
1547   >-------------------------------</a><a name="projectinput">-</a>
1548 
1549   qh_projectinput( )
1550     project input points using qh.lower_bound/upper_bound and qh DELAUNAY
1551     if qh.lower_bound[k]=qh.upper_bound[k]= 0,
1552       removes dimension k
1553     if halfspace intersection
1554       removes dimension k from qh.feasible_point
1555     input points in qh first_point, num_points, input_dim
1556 
1557   returns:
1558     new point array in qh first_point of qh hull_dim coordinates
1559     sets qh POINTSmalloc
1560     if qh DELAUNAY
1561       projects points to paraboloid
1562       lowbound/highbound is also projected
1563     if qh ATinfinity
1564       adds point "at-infinity"
1565     if qh POINTSmalloc
1566       frees old point array
1567 
1568   notes:
1569     checks that qh.hull_dim agrees with qh.input_dim, PROJECTinput, and DELAUNAY
1570 
1571 
1572   design:
1573     sets project[k] to -1 (delete), 0 (keep), 1 (add for Delaunay)
1574     determines newdim and newnum for qh hull_dim and qh num_points
1575     projects points to newpoints
1576     projects qh.lower_bound to itself
1577     projects qh.upper_bound to itself
1578     if qh DELAUNAY
1579       if qh ATINFINITY
1580         projects points to paraboloid
1581         computes "infinity" point as vertex average and 10% above all points
1582       else
1583         uses qh_setdelaunay to project points to paraboloid
1584 */
qh_projectinput(void)1585 void qh_projectinput(void) {
1586   int k,i;
1587   int newdim= qh input_dim, newnum= qh num_points;
1588   signed char *project;
1589   int projectsize= (qh input_dim + 1) * (int)sizeof(*project);
1590   pointT *newpoints, *coord, *infinity;
1591   realT paraboloid, maxboloid= 0;
1592 
1593   project= (signed char *)qh_memalloc(projectsize);
1594   memset((char *)project, 0, (size_t)projectsize);
1595   for (k=0; k < qh input_dim; k++) {   /* skip Delaunay bound */
1596     if (qh lower_bound[k] == 0.0 && qh upper_bound[k] == 0.0) {
1597       project[k]= -1;
1598       newdim--;
1599     }
1600   }
1601   if (qh DELAUNAY) {
1602     project[k]= 1;
1603     newdim++;
1604     if (qh ATinfinity)
1605       newnum++;
1606   }
1607   if (newdim != qh hull_dim) {
1608     qh_memfree(project, projectsize);
1609     qh_fprintf(qh ferr, 6015, "qhull internal error (qh_projectinput): dimension after projection %d != hull_dim %d\n", newdim, qh hull_dim);
1610     qh_errexit(qh_ERRqhull, NULL, NULL);
1611   }
1612   if (!(newpoints= qh temp_malloc= (coordT *)qh_malloc((size_t)(newnum * newdim) * sizeof(coordT)))) {
1613     qh_memfree(project, projectsize);
1614     qh_fprintf(qh ferr, 6016, "qhull error: insufficient memory to project %d points\n",
1615            qh num_points);
1616     qh_errexit(qh_ERRmem, NULL, NULL);
1617   }
1618   /* qh_projectpoints throws error if mismatched dimensions */
1619   qh_projectpoints(project, qh input_dim+1, qh first_point,
1620                     qh num_points, qh input_dim, newpoints, newdim);
1621   trace1((qh ferr, 1003, "qh_projectinput: updating lower and upper_bound\n"));
1622   qh_projectpoints(project, qh input_dim+1, qh lower_bound,
1623                     1, qh input_dim+1, qh lower_bound, newdim+1);
1624   qh_projectpoints(project, qh input_dim+1, qh upper_bound,
1625                     1, qh input_dim+1, qh upper_bound, newdim+1);
1626   if (qh HALFspace) {
1627     if (!qh feasible_point) {
1628       qh_memfree(project, projectsize);
1629       qh_fprintf(qh ferr, 6017, "qhull internal error (qh_projectinput): HALFspace defined without qh.feasible_point\n");
1630       qh_errexit(qh_ERRqhull, NULL, NULL);
1631     }
1632     qh_projectpoints(project, qh input_dim, qh feasible_point,
1633                       1, qh input_dim, qh feasible_point, newdim);
1634   }
1635   qh_memfree(project, projectsize);
1636   if (qh POINTSmalloc)
1637     qh_free(qh first_point);
1638   qh first_point= newpoints;
1639   qh POINTSmalloc= True;
1640   qh temp_malloc= NULL;
1641   if (qh DELAUNAY && qh ATinfinity) {
1642     coord= qh first_point;
1643     infinity= qh first_point + qh hull_dim * qh num_points;
1644     for (k=qh hull_dim-1; k--; )
1645       infinity[k]= 0.0;
1646     for (i=qh num_points; i--; ) {
1647       paraboloid= 0.0;
1648       for (k=0; k < qh hull_dim-1; k++) {
1649         paraboloid += *coord * *coord;
1650         infinity[k] += *coord;
1651         coord++;
1652       }
1653       *(coord++)= paraboloid;
1654       maximize_(maxboloid, paraboloid);
1655     }
1656     /* coord == infinity */
1657     for (k=qh hull_dim-1; k--; )
1658       *(coord++) /= qh num_points;
1659     *(coord++)= maxboloid * 1.1;
1660     qh num_points++;
1661     trace0((qh ferr, 9, "qh_projectinput: projected points to paraboloid for Delaunay\n"));
1662   }else if (qh DELAUNAY)  /* !qh ATinfinity */
1663     qh_setdelaunay(qh hull_dim, qh num_points, qh first_point);
1664 } /* projectinput */
1665 
1666 
1667 /*-<a                             href="qh-geom.htm#TOC"
1668   >-------------------------------</a><a name="projectpoints">-</a>
1669 
1670   qh_projectpoints( project, n, points, numpoints, dim, newpoints, newdim )
1671     project points/numpoints/dim to newpoints/newdim
1672     if project[k] == -1
1673       delete dimension k
1674     if project[k] == 1
1675       add dimension k by duplicating previous column
1676     n is size of project
1677 
1678   notes:
1679     newpoints may be points if only adding dimension at end
1680 
1681   design:
1682     check that 'project' and 'newdim' agree
1683     for each dimension
1684       if project == -1
1685         skip dimension
1686       else
1687         determine start of column in newpoints
1688         determine start of column in points
1689           if project == +1, duplicate previous column
1690         copy dimension (column) from points to newpoints
1691 */
qh_projectpoints(signed char * project,int n,realT * points,int numpoints,int dim,realT * newpoints,int newdim)1692 void qh_projectpoints(signed char *project, int n, realT *points,
1693         int numpoints, int dim, realT *newpoints, int newdim) {
1694   int testdim= dim, oldk=0, newk=0, i,j=0,k;
1695   realT *newp, *oldp;
1696 
1697   for (k=0; k < n; k++)
1698     testdim += project[k];
1699   if (testdim != newdim) {
1700     qh_fprintf(qh ferr, 6018, "qhull internal error (qh_projectpoints): newdim %d should be %d after projection\n",
1701       newdim, testdim);
1702     qh_errexit(qh_ERRqhull, NULL, NULL);
1703   }
1704   for (j=0; j<n; j++) {
1705     if (project[j] == -1)
1706       oldk++;
1707     else {
1708       newp= newpoints+newk++;
1709       if (project[j] == +1) {
1710         if (oldk >= dim)
1711           continue;
1712         oldp= points+oldk;
1713       }else
1714         oldp= points+oldk++;
1715       for (i=numpoints; i--; ) {
1716         *newp= *oldp;
1717         newp += newdim;
1718         oldp += dim;
1719       }
1720     }
1721     if (oldk >= dim)
1722       break;
1723   }
1724   trace1((qh ferr, 1004, "qh_projectpoints: projected %d points from dim %d to dim %d\n",
1725     numpoints, dim, newdim));
1726 } /* projectpoints */
1727 
1728 
1729 /*-<a                             href="qh-geom.htm#TOC"
1730   >-------------------------------</a><a name="rotateinput">-</a>
1731 
1732   qh_rotateinput( rows )
1733     rotate input using row matrix
1734     input points given by qh first_point, num_points, hull_dim
1735     assumes rows[dim] is a scratch buffer
1736     if qh POINTSmalloc, overwrites input points, else mallocs a new array
1737 
1738   returns:
1739     rotated input
1740     sets qh POINTSmalloc
1741 
1742   design:
1743     see qh_rotatepoints
1744 */
qh_rotateinput(realT ** rows)1745 void qh_rotateinput(realT **rows) {
1746 
1747   if (!qh POINTSmalloc) {
1748     qh first_point= qh_copypoints(qh first_point, qh num_points, qh hull_dim);
1749     qh POINTSmalloc= True;
1750   }
1751   qh_rotatepoints(qh first_point, qh num_points, qh hull_dim, rows);
1752 }  /* rotateinput */
1753 
1754 /*-<a                             href="qh-geom.htm#TOC"
1755   >-------------------------------</a><a name="rotatepoints">-</a>
1756 
1757   qh_rotatepoints( points, numpoints, dim, row )
1758     rotate numpoints points by a d-dim row matrix
1759     assumes rows[dim] is a scratch buffer
1760 
1761   returns:
1762     rotated points in place
1763 
1764   design:
1765     for each point
1766       for each coordinate
1767         use row[dim] to compute partial inner product
1768       for each coordinate
1769         rotate by partial inner product
1770 */
qh_rotatepoints(realT * points,int numpoints,int dim,realT ** row)1771 void qh_rotatepoints(realT *points, int numpoints, int dim, realT **row) {
1772   realT *point, *rowi, *coord= NULL, sum, *newval;
1773   int i,j,k;
1774 
1775   if (qh IStracing >= 1)
1776     qh_printmatrix(qh ferr, "qh_rotatepoints: rotate points by", row, dim, dim);
1777   for (point=points, j=numpoints; j--; point += dim) {
1778     newval= row[dim];
1779     for (i=0; i < dim; i++) {
1780       rowi= row[i];
1781       coord= point;
1782       for (sum=0.0, k=dim; k--; )
1783         sum += *rowi++ * *coord++;
1784       *(newval++)= sum;
1785     }
1786     for (k=dim; k--; )
1787       *(--coord)= *(--newval);
1788   }
1789 } /* rotatepoints */
1790 
1791 
1792 /*-<a                             href="qh-geom.htm#TOC"
1793   >-------------------------------</a><a name="scaleinput">-</a>
1794 
1795   qh_scaleinput( )
1796     scale input points using qh low_bound/high_bound
1797     input points given by qh first_point, num_points, hull_dim
1798     if qh POINTSmalloc, overwrites input points, else mallocs a new array
1799 
1800   returns:
1801     scales coordinates of points to low_bound[k], high_bound[k]
1802     sets qh POINTSmalloc
1803 
1804   design:
1805     see qh_scalepoints
1806 */
qh_scaleinput(void)1807 void qh_scaleinput(void) {
1808 
1809   if (!qh POINTSmalloc) {
1810     qh first_point= qh_copypoints(qh first_point, qh num_points, qh hull_dim);
1811     qh POINTSmalloc= True;
1812   }
1813   qh_scalepoints(qh first_point, qh num_points, qh hull_dim,
1814        qh lower_bound, qh upper_bound);
1815 }  /* scaleinput */
1816 
1817 /*-<a                             href="qh-geom.htm#TOC"
1818   >-------------------------------</a><a name="scalelast">-</a>
1819 
1820   qh_scalelast( points, numpoints, dim, low, high, newhigh )
1821     scale last coordinate to [0.0, newhigh], for Delaunay triangulation
1822     input points given by points, numpoints, dim
1823 
1824   returns:
1825     changes scale of last coordinate from [low, high] to [0.0, newhigh]
1826     overwrites last coordinate of each point
1827     saves low/high/newhigh in qh.last_low, etc. for qh_setdelaunay()
1828 
1829   notes:
1830     to reduce precision issues, qh_scalelast makes the last coordinate similar to other coordinates
1831       the last coordinate for Delaunay triangulation is the sum of squares of input coordinates
1832       note that the range [0.0, newwidth] is wrong for narrow distributions with large positive coordinates (e.g., [995933.64, 995963.48])
1833 
1834     when called by qh_setdelaunay, low/high may not match the data passed to qh_setdelaunay
1835 
1836   design:
1837     compute scale and shift factors
1838     apply to last coordinate of each point
1839 */
qh_scalelast(coordT * points,int numpoints,int dim,coordT low,coordT high,coordT newhigh)1840 void qh_scalelast(coordT *points, int numpoints, int dim, coordT low,
1841                    coordT high, coordT newhigh) {
1842   realT scale, shift;
1843   coordT *coord, newlow;
1844   int i;
1845   boolT nearzero= False;
1846 
1847   newlow= 0.0;
1848   trace4((qh ferr, 4013, "qh_scalelast: scale last coordinate from [%2.2g, %2.2g] to [%2.2g, %2.2g]\n",
1849     low, high, newlow, newhigh));
1850   qh last_low= low;
1851   qh last_high= high;
1852   qh last_newhigh= newhigh;
1853   scale= qh_divzero(newhigh - newlow, high - low,
1854                   qh MINdenom_1, &nearzero);
1855   if (nearzero) {
1856     if (qh DELAUNAY)
1857       qh_fprintf(qh ferr, 6019, "qhull input error (qh_scalelast): can not scale last coordinate to [%4.4g, %4.4g].  Input is cocircular or cospherical.   Use option 'Qz' to add a point at infinity.\n",
1858              newlow, newhigh);
1859     else
1860       qh_fprintf(qh ferr, 6020, "qhull input error (qh_scalelast): can not scale last coordinate to [%4.4g, %4.4g].  New bounds are too wide for compared to existing bounds [%4.4g, %4.4g] (width %4.4g)\n",
1861              newlow, newhigh, low, high, high-low);
1862     qh_errexit(qh_ERRinput, NULL, NULL);
1863   }
1864   shift= newlow - low * scale;
1865   coord= points + dim - 1;
1866   for (i=numpoints; i--; coord += dim)
1867     *coord= *coord * scale + shift;
1868 } /* scalelast */
1869 
1870 /*-<a                             href="qh-geom.htm#TOC"
1871   >-------------------------------</a><a name="scalepoints">-</a>
1872 
1873   qh_scalepoints( points, numpoints, dim, newlows, newhighs )
1874     scale points to new lowbound and highbound
1875     retains old bound when newlow= -REALmax or newhigh= +REALmax
1876 
1877   returns:
1878     scaled points
1879     overwrites old points
1880 
1881   design:
1882     for each coordinate
1883       compute current low and high bound
1884       compute scale and shift factors
1885       scale all points
1886       enforce new low and high bound for all points
1887 */
qh_scalepoints(pointT * points,int numpoints,int dim,realT * newlows,realT * newhighs)1888 void qh_scalepoints(pointT *points, int numpoints, int dim,
1889         realT *newlows, realT *newhighs) {
1890   int i,k;
1891   realT shift, scale, *coord, low, high, newlow, newhigh, mincoord, maxcoord;
1892   boolT nearzero= False;
1893 
1894   for (k=0; k < dim; k++) {
1895     newhigh= newhighs[k];
1896     newlow= newlows[k];
1897     if (newhigh > REALmax/2 && newlow < -REALmax/2)
1898       continue;
1899     low= REALmax;
1900     high= -REALmax;
1901     for (i=numpoints, coord=points+k; i--; coord += dim) {
1902       minimize_(low, *coord);
1903       maximize_(high, *coord);
1904     }
1905     if (newhigh > REALmax/2)
1906       newhigh= high;
1907     if (newlow < -REALmax/2)
1908       newlow= low;
1909     if (qh DELAUNAY && k == dim-1 && newhigh < newlow) {
1910       qh_fprintf(qh ferr, 6021, "qhull input error: 'Qb%d' or 'QB%d' inverts paraboloid since high bound %.2g < low bound %.2g\n",
1911                k, k, newhigh, newlow);
1912       qh_errexit(qh_ERRinput, NULL, NULL);
1913     }
1914     scale= qh_divzero(newhigh - newlow, high - low,
1915                   qh MINdenom_1, &nearzero);
1916     if (nearzero) {
1917       qh_fprintf(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",
1918               k, newlow, newhigh, low, high);
1919       qh_errexit(qh_ERRinput, NULL, NULL);
1920     }
1921     shift= (newlow * high - low * newhigh)/(high-low);
1922     coord= points+k;
1923     for (i=numpoints; i--; coord += dim)
1924       *coord= *coord * scale + shift;
1925     coord= points+k;
1926     if (newlow < newhigh) {
1927       mincoord= newlow;
1928       maxcoord= newhigh;
1929     }else {
1930       mincoord= newhigh;
1931       maxcoord= newlow;
1932     }
1933     for (i=numpoints; i--; coord += dim) {
1934       minimize_(*coord, maxcoord);  /* because of roundoff error */
1935       maximize_(*coord, mincoord);
1936     }
1937     trace0((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",
1938       k, low, high, newlow, newhigh, numpoints, scale, shift));
1939   }
1940 } /* scalepoints */
1941 
1942 
1943 /*-<a                             href="qh-geom.htm#TOC"
1944   >-------------------------------</a><a name="setdelaunay">-</a>
1945 
1946   qh_setdelaunay( dim, count, points )
1947     project count points to dim-d paraboloid for Delaunay triangulation
1948 
1949     dim is one more than the dimension of the input set
1950     assumes dim is at least 3 (i.e., at least a 2-d Delaunay triangulation)
1951 
1952     points is a dim*count realT array.  The first dim-1 coordinates
1953     are the coordinates of the first input point.  array[dim] is
1954     the first coordinate of the second input point.  array[2*dim] is
1955     the first coordinate of the third input point.
1956 
1957     if qh.last_low defined (i.e., 'Qbb' called qh_scalelast)
1958       calls qh_scalelast to scale the last coordinate the same as the other points
1959 
1960   returns:
1961     for each point
1962       sets point[dim-1] to sum of squares of coordinates
1963     scale points to 'Qbb' if needed
1964 
1965   notes:
1966     to project one point, use
1967       qh_setdelaunay(qh hull_dim, 1, point)
1968 
1969     Do not use options 'Qbk', 'QBk', or 'QbB' since they scale
1970     the coordinates after the original projection.
1971 
1972 */
qh_setdelaunay(int dim,int count,pointT * points)1973 void qh_setdelaunay(int dim, int count, pointT *points) {
1974   int i, k;
1975   coordT *coordp, coord;
1976   realT paraboloid;
1977 
1978   trace0((qh ferr, 11, "qh_setdelaunay: project %d points to paraboloid for Delaunay triangulation\n", count));
1979   coordp= points;
1980   for (i=0; i < count; i++) {
1981     coord= *coordp++;
1982     paraboloid= coord*coord;
1983     for (k=dim-2; k--; ) {
1984       coord= *coordp++;
1985       paraboloid += coord*coord;
1986     }
1987     *coordp++= paraboloid;
1988   }
1989   if (qh last_low < REALmax/2)
1990     qh_scalelast(points, count, dim, qh last_low, qh last_high, qh last_newhigh);
1991 } /* setdelaunay */
1992 
1993 
1994 /*-<a                             href="qh-geom.htm#TOC"
1995   >-------------------------------</a><a name="sethalfspace">-</a>
1996 
1997   qh_sethalfspace( dim, coords, nextp, normal, offset, feasible )
1998     set point to dual of halfspace relative to feasible point
1999     halfspace is normal coefficients and offset.
2000 
2001   returns:
2002     false and prints error if feasible point is outside of hull
2003     overwrites coordinates for point at dim coords
2004     nextp= next point (coords)
2005     does not call qh_errexit
2006 
2007   design:
2008     compute distance from feasible point to halfspace
2009     divide each normal coefficient by -dist
2010 */
qh_sethalfspace(int dim,coordT * coords,coordT ** nextp,coordT * normal,coordT * offset,coordT * feasible)2011 boolT qh_sethalfspace(int dim, coordT *coords, coordT **nextp,
2012          coordT *normal, coordT *offset, coordT *feasible) {
2013   coordT *normp= normal, *feasiblep= feasible, *coordp= coords;
2014   realT dist;
2015   realT r; /*bug fix*/
2016   int k;
2017   boolT zerodiv;
2018 
2019   dist= *offset;
2020   for (k=dim; k--; )
2021     dist += *(normp++) * *(feasiblep++);
2022   if (dist > 0)
2023     goto LABELerroroutside;
2024   normp= normal;
2025   if (dist < -qh MINdenom) {
2026     for (k=dim; k--; )
2027       *(coordp++)= *(normp++) / -dist;
2028   }else {
2029     for (k=dim; k--; ) {
2030       *(coordp++)= qh_divzero(*(normp++), -dist, qh MINdenom_1, &zerodiv);
2031       if (zerodiv)
2032         goto LABELerroroutside;
2033     }
2034   }
2035   *nextp= coordp;
2036 #ifndef qh_NOtrace
2037   if (qh IStracing >= 4) {
2038     qh_fprintf(qh ferr, 8021, "qh_sethalfspace: halfspace at offset %6.2g to point: ", *offset);
2039     for (k=dim, coordp=coords; k--; ) {
2040       r= *coordp++;
2041       qh_fprintf(qh ferr, 8022, " %6.2g", r);
2042     }
2043     qh_fprintf(qh ferr, 8023, "\n");
2044   }
2045 #endif
2046   return True;
2047 LABELerroroutside:
2048   feasiblep= feasible;
2049   normp= normal;
2050   qh_fprintf(qh ferr, 6023, "qhull input error: feasible point is not clearly inside halfspace\nfeasible point: ");
2051   for (k=dim; k--; )
2052     qh_fprintf(qh ferr, 8024, qh_REAL_1, r=*(feasiblep++));
2053   qh_fprintf(qh ferr, 8025, "\n     halfspace: ");
2054   for (k=dim; k--; )
2055     qh_fprintf(qh ferr, 8026, qh_REAL_1, r=*(normp++));
2056   qh_fprintf(qh ferr, 8027, "\n     at offset: ");
2057   qh_fprintf(qh ferr, 8028, qh_REAL_1, *offset);
2058   qh_fprintf(qh ferr, 8029, " and distance: ");
2059   qh_fprintf(qh ferr, 8030, qh_REAL_1, dist);
2060   qh_fprintf(qh ferr, 8031, "\n");
2061   return False;
2062 } /* sethalfspace */
2063 
2064 /*-<a                             href="qh-geom.htm#TOC"
2065   >-------------------------------</a><a name="sethalfspace_all">-</a>
2066 
2067   qh_sethalfspace_all( dim, count, halfspaces, feasible )
2068     generate dual for halfspace intersection with feasible point
2069     array of count halfspaces
2070       each halfspace is normal coefficients followed by offset
2071       the origin is inside the halfspace if the offset is negative
2072     feasible is a point inside all halfspaces (http://www.qhull.org/html/qhalf.htm#notes)
2073 
2074   returns:
2075     malloc'd array of count X dim-1 points
2076 
2077   notes:
2078     call before qh_init_B or qh_initqhull_globals
2079     free memory when done
2080     unused/untested code: please email bradb@shore.net if this works ok for you
2081     if using option 'Fp', qh.feasible_point must be set (e.g., to 'feasible')
2082     qh feasible_point is a malloc'd array that is freed by qh_freebuffers.
2083 
2084   design:
2085     see qh_sethalfspace
2086 */
qh_sethalfspace_all(int dim,int count,coordT * halfspaces,pointT * feasible)2087 coordT *qh_sethalfspace_all(int dim, int count, coordT *halfspaces, pointT *feasible) {
2088   int i, newdim;
2089   pointT *newpoints;
2090   coordT *coordp, *normalp, *offsetp;
2091 
2092   trace0((qh ferr, 12, "qh_sethalfspace_all: compute dual for halfspace intersection\n"));
2093   newdim= dim - 1;
2094   if (!(newpoints= (coordT *)qh_malloc((size_t)(count * newdim) * sizeof(coordT)))){
2095     qh_fprintf(qh ferr, 6024, "qhull error: insufficient memory to compute dual of %d halfspaces\n",
2096           count);
2097     qh_errexit(qh_ERRmem, NULL, NULL);
2098   }
2099   coordp= newpoints;
2100   normalp= halfspaces;
2101   for (i=0; i < count; i++) {
2102     offsetp= normalp + newdim;
2103     if (!qh_sethalfspace(newdim, coordp, &coordp, normalp, offsetp, feasible)) {
2104       qh_free(newpoints);  /* feasible is not inside halfspace as reported by qh_sethalfspace */
2105       qh_fprintf(qh ferr, 8032, "The halfspace was at index %d\n", i);
2106       qh_errexit(qh_ERRinput, NULL, NULL);
2107     }
2108     normalp= offsetp + 1;
2109   }
2110   return newpoints;
2111 } /* sethalfspace_all */
2112 
2113 
2114 /*-<a                             href="qh-geom.htm#TOC"
2115   >-------------------------------</a><a name="sharpnewfacets">-</a>
2116 
2117   qh_sharpnewfacets( )
2118 
2119   returns:
2120     true if could be an acute angle (facets in different quadrants)
2121 
2122   notes:
2123     for qh_findbest
2124 
2125   design:
2126     for all facets on qh.newfacet_list
2127       if two facets are in different quadrants
2128         set issharp
2129 */
qh_sharpnewfacets(void)2130 boolT qh_sharpnewfacets(void) {
2131   facetT *facet;
2132   boolT issharp= False;
2133   int *quadrant, k;
2134 
2135   quadrant= (int *)qh_memalloc(qh hull_dim * (int)sizeof(int));
2136   FORALLfacet_(qh newfacet_list) {
2137     if (facet == qh newfacet_list) {
2138       for (k=qh hull_dim; k--; )
2139         quadrant[ k]= (facet->normal[ k] > 0);
2140     }else {
2141       for (k=qh hull_dim; k--; ) {
2142         if (quadrant[ k] != (facet->normal[ k] > 0)) {
2143           issharp= True;
2144           break;
2145         }
2146       }
2147     }
2148     if (issharp)
2149       break;
2150   }
2151   qh_memfree(quadrant, qh hull_dim * (int)sizeof(int));
2152   trace3((qh ferr, 3001, "qh_sharpnewfacets: %d\n", issharp));
2153   return issharp;
2154 } /* sharpnewfacets */
2155 
2156 /*-<a                             href="qh-geom.htm#TOC"
2157   >-------------------------------</a><a name="vertex_bestdist">-</a>
2158 
2159   qh_vertex_bestdist( vertices )
2160   qh_vertex_bestdist2( vertices, vertexp, vertexp2 )
2161     return nearest distance between vertices
2162     optionally returns vertex and vertex2
2163 
2164   notes:
2165     called by qh_partitioncoplanar, qh_mergefacet, qh_check_maxout, qh_checkpoint
2166 */
qh_vertex_bestdist(setT * vertices)2167 coordT qh_vertex_bestdist(setT *vertices) {
2168   vertexT *vertex, *vertex2;
2169 
2170   return qh_vertex_bestdist2(vertices, &vertex, &vertex2);
2171 } /* vertex_bestdist */
2172 
qh_vertex_bestdist2(setT * vertices,vertexT ** vertexp,vertexT ** vertexp2)2173 coordT qh_vertex_bestdist2(setT *vertices, vertexT **vertexp/*= NULL*/, vertexT **vertexp2/*= NULL*/) {
2174   vertexT *vertex, *vertexA, *bestvertex= NULL, *bestvertex2= NULL;
2175   coordT dist, bestdist= REALmax;
2176   int k, vertex_i, vertex_n;
2177 
2178   FOREACHvertex_i_(vertices) {
2179     for (k= vertex_i+1; k < vertex_n; k++) {
2180       vertexA= SETelemt_(vertices, k, vertexT);
2181       dist= qh_pointdist(vertex->point, vertexA->point, -qh hull_dim);
2182       if (dist < bestdist) {
2183         bestdist= dist;
2184         bestvertex= vertex;
2185         bestvertex2= vertexA;
2186       }
2187     }
2188   }
2189   *vertexp= bestvertex;
2190   *vertexp2= bestvertex2;
2191   return sqrt(bestdist);
2192 } /* vertex_bestdist */
2193 
2194 /*-<a                             href="qh-geom.htm#TOC"
2195   >-------------------------------</a><a name="voronoi_center">-</a>
2196 
2197   qh_voronoi_center( dim, points )
2198     return Voronoi center for a set of points
2199     dim is the orginal dimension of the points
2200     gh.gm_matrix/qh.gm_row are scratch buffers
2201 
2202   returns:
2203     center as a temporary point (qh_memalloc)
2204     if non-simplicial,
2205       returns center for max simplex of points
2206 
2207   notes:
2208     only called by qh_facetcenter
2209     from Bowyer & Woodwark, A Programmer's Geometry, 1983, p. 65
2210 
2211   design:
2212     if non-simplicial
2213       determine max simplex for points
2214     translate point0 of simplex to origin
2215     compute sum of squares of diagonal
2216     compute determinate
2217     compute Voronoi center (see Bowyer & Woodwark)
2218 */
qh_voronoi_center(int dim,setT * points)2219 pointT *qh_voronoi_center(int dim, setT *points) {
2220   pointT *point, **pointp, *point0;
2221   pointT *center= (pointT *)qh_memalloc(qh center_size);
2222   setT *simplex;
2223   int i, j, k, size= qh_setsize(points);
2224   coordT *gmcoord;
2225   realT *diffp, sum2, *sum2row, *sum2p, det, factor;
2226   boolT nearzero, infinite;
2227 
2228   if (size == dim+1)
2229     simplex= points;
2230   else if (size < dim+1) {
2231     qh_memfree(center, qh center_size);
2232     qh_fprintf(qh ferr, 6025, "qhull internal error (qh_voronoi_center):  need at least %d points to construct a Voronoi center\n",
2233              dim+1);
2234     qh_errexit(qh_ERRqhull, NULL, NULL);
2235     simplex= points;  /* never executed -- avoids warning */
2236   }else {
2237     simplex= qh_settemp(dim+1);
2238     qh_maxsimplex(dim, points, NULL, 0, &simplex);
2239   }
2240   point0= SETfirstt_(simplex, pointT);
2241   gmcoord= qh gm_matrix;
2242   for (k=0; k < dim; k++) {
2243     qh gm_row[k]= gmcoord;
2244     FOREACHpoint_(simplex) {
2245       if (point != point0)
2246         *(gmcoord++)= point[k] - point0[k];
2247     }
2248   }
2249   sum2row= gmcoord;
2250   for (i=0; i < dim; i++) {
2251     sum2= 0.0;
2252     for (k=0; k < dim; k++) {
2253       diffp= qh gm_row[k] + i;
2254       sum2 += *diffp * *diffp;
2255     }
2256     *(gmcoord++)= sum2;
2257   }
2258   det= qh_determinant(qh gm_row, dim, &nearzero);
2259   factor= qh_divzero(0.5, det, qh MINdenom, &infinite);
2260   if (infinite) {
2261     for (k=dim; k--; )
2262       center[k]= qh_INFINITE;
2263     if (qh IStracing)
2264       qh_printpoints(qh ferr, "qh_voronoi_center: at infinity for ", simplex);
2265   }else {
2266     for (i=0; i < dim; i++) {
2267       gmcoord= qh gm_matrix;
2268       sum2p= sum2row;
2269       for (k=0; k < dim; k++) {
2270         qh gm_row[k]= gmcoord;
2271         if (k == i) {
2272           for (j=dim; j--; )
2273             *(gmcoord++)= *sum2p++;
2274         }else {
2275           FOREACHpoint_(simplex) {
2276             if (point != point0)
2277               *(gmcoord++)= point[k] - point0[k];
2278           }
2279         }
2280       }
2281       center[i]= qh_determinant(qh gm_row, dim, &nearzero)*factor + point0[i];
2282     }
2283 #ifndef qh_NOtrace
2284     if (qh IStracing >= 3) {
2285       qh_fprintf(qh ferr, 3061, "qh_voronoi_center: det %2.2g factor %2.2g ", det, factor);
2286       qh_printmatrix(qh ferr, "center:", &center, 1, dim);
2287       if (qh IStracing >= 5) {
2288         qh_printpoints(qh ferr, "points", simplex);
2289         FOREACHpoint_(simplex)
2290           qh_fprintf(qh ferr, 8034, "p%d dist %.2g, ", qh_pointid(point),
2291                    qh_pointdist(point, center, dim));
2292         qh_fprintf(qh ferr, 8035, "\n");
2293       }
2294     }
2295 #endif
2296   }
2297   if (simplex != points)
2298     qh_settempfree(&simplex);
2299   return center;
2300 } /* voronoi_center */
2301 
2302