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