Free cookie consent management tool by TermsFeed Policy Generator

source: branches/HeuristicLab.Analysis.AlgorithmBehavior/qhull-2012.1/src/libqhull/geom2.c @ 10635

Last change on this file since 10635 was 10207, checked in by ascheibe, 11 years ago

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