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

hpp-fcl
Author(s):
autogenerated on Fri Aug 2 2024 02:45:13