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

hpp-fcl
Author(s):
autogenerated on Fri Feb 14 2025 03:45:50