Reverted incorrect merge (missing files)
[blender.git] / extern / qhull / src / qhull.c
1 /*<html><pre>  -<a                             href="qh-qhull.htm"
2   >-------------------------------</a><a name="TOP">-</a>
3
4    qhull.c
5    Quickhull algorithm for convex hulls
6
7    qhull() and top-level routines
8
9    see qh-qhull.htm, qhull.h, unix.c
10
11    see qhull_a.h for internal functions
12
13    copyright (c) 1993-2002 The Geometry Center        
14 */
15
16 #include "qhull_a.h" 
17
18 /*============= functions in alphabetic order after qhull() =======*/
19
20 /*-<a                             href="qh-qhull.htm#TOC"
21   >-------------------------------</a><a name="qhull">-</a>
22   
23   qh_qhull()
24     compute DIM3 convex hull of qh.num_points starting at qh.first_point
25     qh contains all global options and variables
26
27   returns:
28     returns polyhedron
29       qh.facet_list, qh.num_facets, qh.vertex_list, qh.num_vertices,
30     
31     returns global variables
32       qh.hulltime, qh.max_outside, qh.interior_point, qh.max_vertex, qh.min_vertex
33     
34     returns precision constants
35       qh.ANGLEround, centrum_radius, cos_max, DISTround, MAXabs_coord, ONEmerge
36
37   notes:
38     unless needed for output
39       qh.max_vertex and qh.min_vertex are max/min due to merges
40       
41   see:
42     to add individual points to either qh.num_points
43       use qh_addpoint()
44       
45     if qh.GETarea
46       qh_produceoutput() returns qh.totarea and qh.totvol via qh_getarea()
47
48   design:
49     record starting time
50     initialize hull and partition points
51     build convex hull
52     unless early termination
53       update facet->maxoutside for vertices, coplanar, and near-inside points
54     error if temporary sets exist
55     record end time
56 */
57 void qh_qhull (void) {
58   int numoutside;
59
60   qh hulltime= qh_CPUclock;
61   if (qh RERUN || qh JOGGLEmax < REALmax/2) 
62     qh_build_withrestart();
63   else {
64     qh_initbuild();
65     qh_buildhull();
66   }
67   if (!qh STOPpoint && !qh STOPcone) {
68     if (qh ZEROall_ok && !qh TESTvneighbors && qh MERGEexact)
69       qh_checkzero( qh_ALL);
70     if (qh ZEROall_ok && !qh TESTvneighbors && !qh WAScoplanar) {
71       trace2((qh ferr, "qh_qhull: all facets are clearly convex and no coplanar points.  Post-merging and check of maxout not needed.\n"));
72       qh DOcheckmax= False;
73     }else {
74       if (qh MERGEexact || (qh hull_dim > qh_DIMreduceBuild && qh PREmerge))
75         qh_postmerge ("First post-merge", qh premerge_centrum, qh premerge_cos, 
76              (qh POSTmerge ? False : qh TESTvneighbors));
77       else if (!qh POSTmerge && qh TESTvneighbors) 
78         qh_postmerge ("For testing vertex neighbors", qh premerge_centrum,
79              qh premerge_cos, True); 
80       if (qh POSTmerge)
81         qh_postmerge ("For post-merging", qh postmerge_centrum, 
82              qh postmerge_cos, qh TESTvneighbors);
83       if (qh visible_list == qh facet_list) { /* i.e., merging done */
84         qh findbestnew= True;
85         qh_partitionvisible (/*visible_list, newfacet_list*/ !qh_ALL, &numoutside);
86         qh findbestnew= False;
87         qh_deletevisible (/*qh visible_list*/);
88         qh_resetlists (False, qh_RESETvisible /*qh visible_list newvertex_list newfacet_list */);
89       }
90     }
91     if (qh DOcheckmax){
92       if (qh REPORTfreq) {
93         qh_buildtracing (NULL, NULL); 
94         fprintf (qh ferr, "\nTesting all coplanar points.\n");
95       }
96       qh_check_maxout();
97     }
98     if (qh KEEPnearinside && !qh maxoutdone)  
99       qh_nearcoplanar();
100   }
101   if (qh_setsize ((setT*)qhmem.tempstack) != 0) {
102     fprintf (qh ferr, "qhull internal error (qh_qhull): temporary sets not empty (%d)\n",
103              qh_setsize ((setT*)qhmem.tempstack));
104     qh_errexit (qh_ERRqhull, NULL, NULL);
105   }
106   qh hulltime= qh_CPUclock - qh hulltime;
107   qh QHULLfinished= True;
108   trace1((qh ferr, "qh_qhull: algorithm completed\n"));
109 } /* qhull */
110
111 /*-<a                             href="qh-qhull.htm#TOC"
112   >-------------------------------</a><a name="addpoint">-</a>
113   
114   qh_addpoint( furthest, facet, checkdist )
115     add point (usually furthest point) above facet to hull 
116     if checkdist, 
117       check that point is above facet.
118       if point is not outside of the hull, uses qh_partitioncoplanar()
119       assumes that facet is defined by qh_findbestfacet()
120     else if facet specified,
121       assumes that point is above facet (major damage if below)
122     for Delaunay triangulations, 
123       Use qh_setdelaunay() to lift point to paraboloid and scale by 'Qbb' if needed
124       Do not use options 'Qbk', 'QBk', or 'QbB' since they scale the coordinates. 
125
126   returns:
127     returns False if user requested an early termination
128      qh.visible_list, newfacet_list, delvertex_list, NEWfacets may be defined
129     updates qh.facet_list, qh.num_facets, qh.vertex_list, qh.num_vertices
130     clear qh.maxoutdone (will need to call qh_check_maxout() for facet->maxoutside)
131     if unknown point, adds a pointer to qh.other_points
132       do not deallocate the point's coordinates
133   
134   notes:
135     assumes point is near its best facet and not at a local minimum of a lens
136       distributions.  Use qh_findbestfacet to avoid this case.
137     uses qh.visible_list, qh.newfacet_list, qh.delvertex_list, qh.NEWfacets
138
139   see also:
140     qh_triangulate() -- triangulate non-simplicial facets
141
142   design:
143     check point in qh.first_point/.num_points
144     if checkdist
145       if point not above facet
146         partition coplanar point 
147         exit
148     exit if pre STOPpoint requested
149     find horizon and visible facets for point
150     make new facets for point to horizon
151     make hyperplanes for point
152     compute balance statistics
153     match neighboring new facets
154     update vertex neighbors and delete interior vertices
155     exit if STOPcone requested
156     merge non-convex new facets
157     if merge found, many merges, or 'Qf'
158        use qh_findbestnew() instead of qh_findbest()
159     partition outside points from visible facets
160     delete visible facets
161     check polyhedron if requested
162     exit if post STOPpoint requested
163     reset working lists of facets and vertices
164 */
165 boolT qh_addpoint (pointT *furthest, facetT *facet, boolT checkdist) {
166   int goodvisible, goodhorizon;
167   vertexT *vertex;
168   facetT *newfacet;
169   realT dist, newbalance, pbalance;
170   boolT isoutside= False;
171   int numpart, numpoints, numnew, firstnew;
172
173   qh maxoutdone= False;
174   if (qh_pointid (furthest) == -1)
175     qh_setappend (&qh other_points, furthest);
176   if (!facet) {
177     fprintf (qh ferr, "qh_addpoint: NULL facet.  Need to call qh_findbestfacet first\n");
178     qh_errexit (qh_ERRqhull, NULL, NULL);
179   }
180   if (checkdist) {
181     facet= qh_findbest (furthest, facet, !qh_ALL, !qh_ISnewfacets, !qh_NOupper,
182                         &dist, &isoutside, &numpart);
183     zzadd_(Zpartition, numpart);
184     if (!isoutside) {
185       zinc_(Znotmax);  /* last point of outsideset is no longer furthest. */
186       facet->notfurthest= True;
187       qh_partitioncoplanar (furthest, facet, &dist);
188       return True;
189     }
190   }
191   qh_buildtracing (furthest, facet);
192   if (qh STOPpoint < 0 && qh furthest_id == -qh STOPpoint-1) {
193     facet->notfurthest= True;
194     return False;
195   }
196   qh_findhorizon (furthest, facet, &goodvisible, &goodhorizon); 
197   if (qh ONLYgood && !(goodvisible+goodhorizon) && !qh GOODclosest) {
198     zinc_(Znotgood);  
199     facet->notfurthest= True;
200     /* last point of outsideset is no longer furthest.  This is ok
201        since all points of the outside are likely to be bad */
202     qh_resetlists (False, qh_RESETvisible /*qh visible_list newvertex_list newfacet_list */);
203     return True;
204   }
205   zzinc_(Zprocessed);
206   firstnew= qh facet_id;
207   vertex= qh_makenewfacets (furthest /*visible_list, attaches if !ONLYgood */);
208   qh_makenewplanes (/* newfacet_list */);
209   numnew= qh facet_id - firstnew;
210   newbalance= numnew - (realT) (qh num_facets-qh num_visible)
211                          * qh hull_dim/qh num_vertices;
212   wadd_(Wnewbalance, newbalance);
213   wadd_(Wnewbalance2, newbalance * newbalance);
214   if (qh ONLYgood 
215   && !qh_findgood (qh newfacet_list, goodhorizon) && !qh GOODclosest) {
216     FORALLnew_facets 
217       qh_delfacet (newfacet);
218     qh_delvertex (vertex);
219     qh_resetlists (True, qh_RESETvisible /*qh visible_list newvertex_list newfacet_list */);
220     zinc_(Znotgoodnew);
221     facet->notfurthest= True;
222     return True;
223   }
224   if (qh ONLYgood)
225     qh_attachnewfacets(/*visible_list*/);
226   qh_matchnewfacets();
227   qh_updatevertices();
228   if (qh STOPcone && qh furthest_id == qh STOPcone-1) {
229     facet->notfurthest= True;
230     return False;  /* visible_list etc. still defined */
231   }
232   qh findbestnew= False;
233   if (qh PREmerge || qh MERGEexact) {
234     qh_premerge (vertex, qh premerge_centrum, qh premerge_cos);
235     if (qh_USEfindbestnew)
236       qh findbestnew= True;
237     else {
238       FORALLnew_facets {
239         if (!newfacet->simplicial) {
240           qh findbestnew= True;  /* use qh_findbestnew instead of qh_findbest*/
241           break;
242         }
243       }
244     }
245   }else if (qh BESToutside)
246     qh findbestnew= True;
247   qh_partitionvisible (/*visible_list, newfacet_list*/ !qh_ALL, &numpoints);
248   qh findbestnew= False;
249   qh findbest_notsharp= False;
250   zinc_(Zpbalance);
251   pbalance= numpoints - (realT) qh hull_dim /* assumes all points extreme */
252                 * (qh num_points - qh num_vertices)/qh num_vertices;
253   wadd_(Wpbalance, pbalance);
254   wadd_(Wpbalance2, pbalance * pbalance);
255   qh_deletevisible (/*qh visible_list*/);
256   zmax_(Zmaxvertex, qh num_vertices);
257   qh NEWfacets= False;
258   if (qh IStracing >= 4) {
259     if (qh num_facets < 2000)
260       qh_printlists();
261     qh_printfacetlist (qh newfacet_list, NULL, True);
262     qh_checkpolygon (qh facet_list);
263   }else if (qh CHECKfrequently) {
264     if (qh num_facets < 50)
265       qh_checkpolygon (qh facet_list);
266     else
267       qh_checkpolygon (qh newfacet_list);
268   }
269   if (qh STOPpoint > 0 && qh furthest_id == qh STOPpoint-1) 
270     return False; 
271   qh_resetlists (True, qh_RESETvisible /*qh visible_list newvertex_list newfacet_list */);
272   /* qh_triangulate(); to test qh.TRInormals */
273   trace2((qh ferr, "qh_addpoint: added p%d new facets %d new balance %2.2g point balance %2.2g\n",
274     qh_pointid (furthest), numnew, newbalance, pbalance));
275   return True;
276 } /* addpoint */
277
278 /*-<a                             href="qh-qhull.htm#TOC"
279   >-------------------------------</a><a name="build_withrestart">-</a>
280   
281   qh_build_withrestart()
282     allow restarts due to qh.JOGGLEmax while calling qh_buildhull()
283     qh.FIRSTpoint/qh.NUMpoints is point array
284         it may be moved by qh_joggleinput()
285 */
286 void qh_build_withrestart (void) {
287   int restart;
288
289   qh ALLOWrestart= True;
290   while (True) {
291     restart= setjmp (qh restartexit); /* simple statement for CRAY J916 */
292     if (restart) {       /* only from qh_precision() */
293       zzinc_(Zretry);
294       wmax_(Wretrymax, qh JOGGLEmax);
295       qh ERREXITcalled= False;
296       qh STOPcone= True; /* if break, prevents normal output */
297     }
298     if (!qh RERUN && qh JOGGLEmax < REALmax/2) {
299       if (qh build_cnt > qh_JOGGLEmaxretry) {
300         fprintf(qh ferr, "\n\
301 qhull precision error: %d attempts to construct a convex hull\n\
302         with joggled input.  Increase joggle above 'QJ%2.2g'\n\
303         or modify qh_JOGGLE... parameters in user.h\n",
304            qh build_cnt, qh JOGGLEmax);
305         qh_errexit (qh_ERRqhull, NULL, NULL);
306       }
307       if (qh build_cnt && !restart)
308         break;
309     }else if (qh build_cnt && qh build_cnt >= qh RERUN)
310       break;
311     qh STOPcone= False;
312     qh_freebuild (True);  /* first call is a nop */
313     qh build_cnt++;
314     if (!qh qhull_optionsiz)
315       qh qhull_optionsiz= strlen (qh qhull_options);
316     else { 
317       qh qhull_options [qh qhull_optionsiz]= '\0';
318       qh qhull_optionlen= 80;
319     }
320     qh_option("_run", &qh build_cnt, NULL);
321     if (qh build_cnt == qh RERUN) {
322       qh IStracing= qh TRACElastrun;  /* duplicated from qh_initqhull_globals */
323       if (qh TRACEpoint != -1 || qh TRACEdist < REALmax/2 || qh TRACEmerge) {
324         qh TRACElevel= (qh IStracing? qh IStracing : 3);
325         qh IStracing= 0;
326       }
327       qhmem.IStracing= qh IStracing;
328     }
329     if (qh JOGGLEmax < REALmax/2)
330       qh_joggleinput();
331     qh_initbuild();
332     qh_buildhull();
333     if (qh JOGGLEmax < REALmax/2 && !qh MERGING)
334       qh_checkconvex (qh facet_list, qh_ALGORITHMfault);
335   }
336   qh ALLOWrestart= False;
337 } /* qh_build_withrestart */
338
339 /*-<a                             href="qh-qhull.htm#TOC"
340   >-------------------------------</a><a name="buildhull">-</a>
341   
342   qh_buildhull()
343     construct a convex hull by adding outside points one at a time
344
345   returns:
346   
347   notes:
348     may be called multiple times
349     checks facet and vertex lists for incorrect flags
350     to recover from STOPcone, call qh_deletevisible and qh_resetlists
351
352   design:
353     check visible facet and newfacet flags
354     check newlist vertex flags and qh.STOPcone/STOPpoint
355     for each facet with a furthest outside point
356       add point to facet
357       exit if qh.STOPcone or qh.STOPpoint requested
358     if qh.NARROWhull for initial simplex
359       partition remaining outside points to coplanar sets
360 */
361 void qh_buildhull(void) {
362   facetT *facet;
363   pointT *furthest;
364   vertexT *vertex;
365   int id;
366   
367   trace1((qh ferr, "qh_buildhull: start build hull\n"));
368   FORALLfacets {
369     if (facet->visible || facet->newfacet) {
370       fprintf (qh ferr, "qhull internal error (qh_buildhull): visible or new facet f%d in facet list\n",
371                    facet->id);    
372       qh_errexit (qh_ERRqhull, facet, NULL);
373     }
374   }
375   FORALLvertices {
376     if (vertex->newlist) {
377       fprintf (qh ferr, "qhull internal error (qh_buildhull): new vertex f%d in vertex list\n",
378                    vertex->id);
379       qh_errprint ("ERRONEOUS", NULL, NULL, NULL, vertex);
380       qh_errexit (qh_ERRqhull, NULL, NULL);
381     }
382     id= qh_pointid (vertex->point);
383     if ((qh STOPpoint>0 && id == qh STOPpoint-1) ||
384         (qh STOPpoint<0 && id == -qh STOPpoint-1) ||
385         (qh STOPcone>0 && id == qh STOPcone-1)) {
386       trace1((qh ferr,"qh_buildhull: stop point or cone P%d in initial hull\n", id));
387       return;
388     }
389   }
390   qh facet_next= qh facet_list;      /* advance facet when processed */
391   while ((furthest= qh_nextfurthest (&facet))) {
392     qh num_outside--;  /* if ONLYmax, furthest may not be outside */
393     if (!qh_addpoint (furthest, facet, qh ONLYmax))
394       break;
395   }
396   if (qh NARROWhull) /* move points from outsideset to coplanarset */
397     qh_outcoplanar( /* facet_list */ );
398   if (qh num_outside && !furthest) {
399     fprintf (qh ferr, "qhull internal error (qh_buildhull): %d outside points were never processed.\n", qh num_outside);
400     qh_errexit (qh_ERRqhull, NULL, NULL);
401   }
402   trace1((qh ferr, "qh_buildhull: completed the hull construction\n"));
403 } /* buildhull */
404   
405
406 /*-<a                             href="qh-qhull.htm#TOC"
407   >-------------------------------</a><a name="buildtracing">-</a>
408   
409   qh_buildtracing( furthest, facet )
410     trace an iteration of qh_buildhull() for furthest point and facet
411     if !furthest, prints progress message
412
413   returns:
414     tracks progress with qh.lastreport
415     updates qh.furthest_id (-3 if furthest is NULL)
416     also resets visit_id, vertext_visit on wrap around
417
418   see:
419     qh_tracemerging()
420
421   design:
422     if !furthest
423       print progress message
424       exit
425     if 'TFn' iteration
426       print progress message
427     else if tracing
428       trace furthest point and facet
429     reset qh.visit_id and qh.vertex_visit if overflow may occur
430     set qh.furthest_id for tracing
431 */
432 void qh_buildtracing (pointT *furthest, facetT *facet) {
433   realT dist= 0;
434   float cpu;
435   int total, furthestid;
436   time_t timedata;
437   struct tm *tp;
438   vertexT *vertex;
439
440   qh old_randomdist= qh RANDOMdist;
441   qh RANDOMdist= False;
442   if (!furthest) {
443     time (&timedata);
444     tp= localtime (&timedata);
445     cpu= qh_CPUclock - qh hulltime;
446     cpu /= qh_SECticks;
447     total= zzval_(Ztotmerge) - zzval_(Zcyclehorizon) + zzval_(Zcyclefacettot);
448     fprintf (qh ferr, "\n\
449 At %02d:%02d:%02d & %2.5g CPU secs, qhull has created %d facets and merged %d.\n\
450  The current hull contains %d facets and %d vertices.  Last point was p%d\n",
451       tp->tm_hour, tp->tm_min, tp->tm_sec, cpu, qh facet_id -1,
452       total, qh num_facets, qh num_vertices, qh furthest_id);
453     return;
454   }
455   furthestid= qh_pointid (furthest);
456   if (qh TRACEpoint == furthestid) {
457     qh IStracing= qh TRACElevel;
458     qhmem.IStracing= qh TRACElevel;
459   }else if (qh TRACEpoint != -1 && qh TRACEdist < REALmax/2) {
460     qh IStracing= 0;
461     qhmem.IStracing= 0;
462   }
463   if (qh REPORTfreq && (qh facet_id-1 > qh lastreport+qh REPORTfreq)) {
464     qh lastreport= qh facet_id-1;
465     time (&timedata);
466     tp= localtime (&timedata);
467     cpu= qh_CPUclock - qh hulltime;
468     cpu /= qh_SECticks;
469     total= zzval_(Ztotmerge) - zzval_(Zcyclehorizon) + zzval_(Zcyclefacettot);
470     zinc_(Zdistio);
471     qh_distplane (furthest, facet, &dist);
472     fprintf (qh ferr, "\n\
473 At %02d:%02d:%02d & %2.5g CPU secs, qhull has created %d facets and merged %d.\n\
474  The current hull contains %d facets and %d vertices.  There are %d\n\
475  outside points.  Next is point p%d (v%d), %2.2g above f%d.\n",
476       tp->tm_hour, tp->tm_min, tp->tm_sec, cpu, qh facet_id -1,
477       total, qh num_facets, qh num_vertices, qh num_outside+1,
478       furthestid, qh vertex_id, dist, getid_(facet));
479   }else if (qh IStracing >=1) {
480     cpu= qh_CPUclock - qh hulltime;
481     cpu /= qh_SECticks;
482     qh_distplane (furthest, facet, &dist);
483     fprintf (qh ferr, "qh_addpoint: add p%d (v%d) to hull of %d facets (%2.2g above f%d) and %d outside at %4.4g CPU secs.  Previous was p%d.\n",
484       furthestid, qh vertex_id, qh num_facets, dist,
485       getid_(facet), qh num_outside+1, cpu, qh furthest_id);
486   }
487   if (qh visit_id > (unsigned) INT_MAX) {
488     qh visit_id= 0;
489     FORALLfacets
490       facet->visitid= qh visit_id;
491   }
492   if (qh vertex_visit > (unsigned) INT_MAX) {
493     qh vertex_visit= 0;
494     FORALLvertices
495       vertex->visitid= qh vertex_visit;
496   }
497   qh furthest_id= furthestid;
498   qh RANDOMdist= qh old_randomdist;
499 } /* buildtracing */
500
501 /*-<a                             href="qh-qhull.htm#TOC"
502   >-------------------------------</a><a name="errexit2">-</a>
503   
504   qh_errexit2( exitcode, facet, otherfacet )
505     return exitcode to system after an error
506     report two facets
507
508   returns:
509     assumes exitcode non-zero
510
511   see:
512     normally use qh_errexit() in user.c (reports a facet and a ridge)
513 */
514 void qh_errexit2(int exitcode, facetT *facet, facetT *otherfacet) {
515   
516   qh_errprint("ERRONEOUS", facet, otherfacet, NULL, NULL);
517   qh_errexit (exitcode, NULL, NULL);
518 } /* errexit2 */
519
520
521 /*-<a                             href="qh-qhull.htm#TOC"
522   >-------------------------------</a><a name="findhorizon">-</a>
523   
524   qh_findhorizon( point, facet, goodvisible, goodhorizon )
525     given a visible facet, find the point's horizon and visible facets
526     for all facets, !facet-visible
527
528   returns:
529     returns qh.visible_list/num_visible with all visible facets 
530       marks visible facets with ->visible 
531     updates count of good visible and good horizon facets
532     updates qh.max_outside, qh.max_vertex, facet->maxoutside
533
534   see:
535     similar to qh_delpoint()
536
537   design:
538     move facet to qh.visible_list at end of qh.facet_list
539     for all visible facets
540      for each unvisited neighbor of a visible facet
541        compute distance of point to neighbor
542        if point above neighbor
543          move neighbor to end of qh.visible_list
544        else if point is coplanar with neighbor
545          update qh.max_outside, qh.max_vertex, neighbor->maxoutside
546          mark neighbor coplanar (will create a samecycle later)
547          update horizon statistics         
548 */
549 void qh_findhorizon(pointT *point, facetT *facet, int *goodvisible, int *goodhorizon) {
550   facetT *neighbor, **neighborp, *visible;
551   int numhorizon= 0, coplanar= 0;
552   realT dist;
553   
554   trace1((qh ferr,"qh_findhorizon: find horizon for point p%d facet f%d\n",qh_pointid(point),facet->id));
555   *goodvisible= *goodhorizon= 0;
556   zinc_(Ztotvisible);
557   qh_removefacet(facet);  /* visible_list at end of qh facet_list */
558   qh_appendfacet(facet);
559   qh num_visible= 1;
560   if (facet->good)
561     (*goodvisible)++;
562   qh visible_list= facet;
563   facet->visible= True;
564   facet->f.replace= NULL;
565   if (qh IStracing >=4)
566     qh_errprint ("visible", facet, NULL, NULL, NULL);
567   qh visit_id++;
568   FORALLvisible_facets {
569     if (visible->tricoplanar && !qh TRInormals) {
570       fprintf (qh ferr, "qh_findhorizon: does not work for tricoplanar facets.  Use option 'Q11'\n");
571       qh_errexit (qh_ERRqhull, visible, NULL);
572     }
573     visible->visitid= qh visit_id;
574     FOREACHneighbor_(visible) {
575       if (neighbor->visitid == qh visit_id) 
576         continue;
577       neighbor->visitid= qh visit_id;
578       zzinc_(Znumvisibility);
579       qh_distplane(point, neighbor, &dist);
580       if (dist > qh MINvisible) {
581         zinc_(Ztotvisible);
582         qh_removefacet(neighbor);  /* append to end of qh visible_list */
583         qh_appendfacet(neighbor);
584         neighbor->visible= True;
585         neighbor->f.replace= NULL;
586         qh num_visible++;
587         if (neighbor->good)
588           (*goodvisible)++;
589         if (qh IStracing >=4)
590           qh_errprint ("visible", neighbor, NULL, NULL, NULL);
591       }else {
592         if (dist > - qh MAXcoplanar) {
593           neighbor->coplanar= True;
594           zzinc_(Zcoplanarhorizon);
595           qh_precision ("coplanar horizon");
596           coplanar++;
597           if (qh MERGING) {
598             if (dist > 0) {
599               maximize_(qh max_outside, dist);
600               maximize_(qh max_vertex, dist);
601 #if qh_MAXoutside
602               maximize_(neighbor->maxoutside, dist);
603 #endif
604             }else
605               minimize_(qh min_vertex, dist);  /* due to merge later */
606           }
607           trace2((qh ferr, "qh_findhorizon: point p%d is coplanar to horizon f%d, dist=%2.7g < qh MINvisible (%2.7g)\n",
608               qh_pointid(point), neighbor->id, dist, qh MINvisible));
609         }else
610           neighbor->coplanar= False;
611         zinc_(Ztothorizon);
612         numhorizon++;
613         if (neighbor->good)
614           (*goodhorizon)++;
615         if (qh IStracing >=4)
616           qh_errprint ("horizon", neighbor, NULL, NULL, NULL);
617       }
618     }
619   }
620   if (!numhorizon) {
621     qh_precision ("empty horizon");
622     fprintf(qh ferr, "qhull precision error (qh_findhorizon): empty horizon\n\
623 Point p%d was above all facets.\n", qh_pointid(point));
624     qh_printfacetlist (qh facet_list, NULL, True);
625     qh_errexit(qh_ERRprec, NULL, NULL);
626   }
627   trace1((qh ferr, "qh_findhorizon: %d horizon facets (good %d), %d visible (good %d), %d coplanar\n", 
628        numhorizon, *goodhorizon, qh num_visible, *goodvisible, coplanar));
629   if (qh IStracing >= 4 && qh num_facets < 50) 
630     qh_printlists ();
631 } /* findhorizon */
632
633 /*-<a                             href="qh-qhull.htm#TOC"
634   >-------------------------------</a><a name="nextfurthest">-</a>
635   
636   qh_nextfurthest( visible )
637     returns next furthest point and visible facet for qh_addpoint()
638     starts search at qh.facet_next
639
640   returns:
641     removes furthest point from outside set
642     NULL if none available
643     advances qh.facet_next over facets with empty outside sets  
644
645   design:
646     for each facet from qh.facet_next
647       if empty outside set
648         advance qh.facet_next
649       else if qh.NARROWhull
650         determine furthest outside point
651         if furthest point is not outside
652           advance qh.facet_next (point will be coplanar)
653     remove furthest point from outside set
654 */
655 pointT *qh_nextfurthest (facetT **visible) {
656   facetT *facet;
657   int size, index;
658   realT randr, dist;
659   pointT *furthest;
660
661   while ((facet= qh facet_next) != qh facet_tail) {
662     if (!facet->outsideset) {
663       qh facet_next= facet->next;
664       continue;
665     }
666     SETreturnsize_(facet->outsideset, size);
667     if (!size) {
668       qh_setfree (&facet->outsideset);
669       qh facet_next= facet->next;
670       continue;
671     }
672     if (qh NARROWhull) {
673       if (facet->notfurthest) 
674         qh_furthestout (facet);
675       furthest= (pointT*)qh_setlast (facet->outsideset);
676 #if qh_COMPUTEfurthest
677       qh_distplane (furthest, facet, &dist);
678       zinc_(Zcomputefurthest);
679 #else
680       dist= facet->furthestdist;
681 #endif
682       if (dist < qh MINoutside) { /* remainder of outside set is coplanar for qh_outcoplanar */
683         qh facet_next= facet->next;
684         continue;
685       }
686     }
687     if (!qh RANDOMoutside && !qh VIRTUALmemory) {
688       if (qh PICKfurthest) {
689         qh_furthestnext (/* qh facet_list */);
690         facet= qh facet_next;
691       }
692       *visible= facet;
693       return ((pointT*)qh_setdellast (facet->outsideset));
694     }
695     if (qh RANDOMoutside) {
696       int outcoplanar = 0;
697       if (qh NARROWhull) {
698         FORALLfacets {
699           if (facet == qh facet_next)
700             break;
701           if (facet->outsideset)
702             outcoplanar += qh_setsize( facet->outsideset);
703         }
704       }
705       randr= qh_RANDOMint;
706       randr= randr/(qh_RANDOMmax+1);
707       index= (int)floor((qh num_outside - outcoplanar) * randr);
708       FORALLfacet_(qh facet_next) {
709         if (facet->outsideset) {
710           SETreturnsize_(facet->outsideset, size);
711           if (!size)
712             qh_setfree (&facet->outsideset);
713           else if (size > index) {
714             *visible= facet;
715             return ((pointT*)qh_setdelnth (facet->outsideset, index));
716           }else
717             index -= size;
718         }
719       }
720       fprintf (qh ferr, "qhull internal error (qh_nextfurthest): num_outside %d is too low\nby at least %d, or a random real %g >= 1.0\n",
721               qh num_outside, index+1, randr);
722       qh_errexit (qh_ERRqhull, NULL, NULL);
723     }else { /* VIRTUALmemory */
724       facet= qh facet_tail->previous;
725       if (!(furthest= (pointT*)qh_setdellast(facet->outsideset))) {
726         if (facet->outsideset)
727           qh_setfree (&facet->outsideset);
728         qh_removefacet (facet);
729         qh_prependfacet (facet, &qh facet_list);
730         continue;
731       }
732       *visible= facet;
733       return furthest;
734     }
735   }
736   return NULL;
737 } /* nextfurthest */
738
739 /*-<a                             href="qh-qhull.htm#TOC"
740   >-------------------------------</a><a name="partitionall">-</a>
741   
742   qh_partitionall( vertices, points, numpoints )
743     partitions all points in points/numpoints to the outsidesets of facets
744     vertices= vertices in qh.facet_list (not partitioned)
745
746   returns:
747     builds facet->outsideset
748     does not partition qh.GOODpoint
749     if qh.ONLYgood && !qh.MERGING, 
750       does not partition qh.GOODvertex
751
752   notes:
753     faster if qh.facet_list sorted by anticipated size of outside set
754
755   design:
756     initialize pointset with all points
757     remove vertices from pointset
758     remove qh.GOODpointp from pointset (unless it's qh.STOPcone or qh.STOPpoint)
759     for all facets
760       for all remaining points in pointset
761         compute distance from point to facet
762         if point is outside facet
763           remove point from pointset (by not reappending)
764           update bestpoint
765           append point or old bestpoint to facet's outside set
766       append bestpoint to facet's outside set (furthest)
767     for all points remaining in pointset
768       partition point into facets' outside sets and coplanar sets
769 */
770 void qh_partitionall(setT *vertices, pointT *points, int numpoints){
771   setT *pointset;
772   vertexT *vertex, **vertexp;
773   pointT *point, **pointp, *bestpoint;
774   int size, point_i, point_n, point_end, remaining, i, id;
775   facetT *facet;
776   realT bestdist= -REALmax, dist, distoutside;
777     
778   trace1((qh ferr, "qh_partitionall: partition all points into outside sets\n"));
779   pointset= qh_settemp (numpoints);
780   qh num_outside= 0;
781   pointp= SETaddr_(pointset, pointT);
782   for (i=numpoints, point= points; i--; point += qh hull_dim)
783     *(pointp++)= point;
784   qh_settruncate (pointset, numpoints);
785   FOREACHvertex_(vertices) {
786     if ((id= qh_pointid(vertex->point)) >= 0)
787       SETelem_(pointset, id)= NULL;
788   }
789   id= qh_pointid (qh GOODpointp);
790   if (id >=0 && qh STOPcone-1 != id && -qh STOPpoint-1 != id)
791     SETelem_(pointset, id)= NULL;
792   if (qh GOODvertexp && qh ONLYgood && !qh MERGING) { /* matches qhull()*/
793     if ((id= qh_pointid(qh GOODvertexp)) >= 0)
794       SETelem_(pointset, id)= NULL;
795   }
796   if (!qh BESToutside) {  /* matches conditional for qh_partitionpoint below */
797     distoutside= qh_DISToutside; /* multiple of qh.MINoutside & qh.max_outside, see user.h */
798     zval_(Ztotpartition)= qh num_points - qh hull_dim - 1; /*misses GOOD... */
799     remaining= qh num_facets;
800     point_end= numpoints;
801     FORALLfacets {
802       size= point_end/(remaining--) + 100;
803       facet->outsideset= qh_setnew (size);
804       bestpoint= NULL;
805       point_end= 0;
806       FOREACHpoint_i_(pointset) {
807         if (point) {
808           zzinc_(Zpartitionall);
809           qh_distplane (point, facet, &dist);
810           if (dist < distoutside)
811             SETelem_(pointset, point_end++)= point;
812           else {
813             qh num_outside++;
814             if (!bestpoint) {
815               bestpoint= point;
816               bestdist= dist;
817             }else if (dist > bestdist) {
818               qh_setappend (&facet->outsideset, bestpoint);
819               bestpoint= point;
820               bestdist= dist;
821             }else 
822               qh_setappend (&facet->outsideset, point);
823           }
824         }
825       }
826       if (bestpoint) {
827         qh_setappend (&facet->outsideset, bestpoint);
828 #if !qh_COMPUTEfurthest
829         facet->furthestdist= bestdist;
830 #endif
831       }else
832         qh_setfree (&facet->outsideset);
833       qh_settruncate (pointset, point_end);
834     }
835   }
836   /* if !qh BESToutside, pointset contains points not assigned to outsideset */
837   if (qh BESToutside || qh MERGING || qh KEEPcoplanar || qh KEEPinside) {
838     qh findbestnew= True;
839     FOREACHpoint_i_(pointset) { 
840       if (point)
841         qh_partitionpoint(point, qh facet_list);
842     }
843     qh findbestnew= False;
844   }
845   zzadd_(Zpartitionall, zzval_(Zpartition));
846   zzval_(Zpartition)= 0;
847   qh_settempfree(&pointset);
848   if (qh IStracing >= 4)
849     qh_printfacetlist (qh facet_list, NULL, True);
850 } /* partitionall */
851
852
853 /*-<a                             href="qh-qhull.htm#TOC"
854   >-------------------------------</a><a name="partitioncoplanar">-</a>
855   
856   qh_partitioncoplanar( point, facet, dist )
857     partition coplanar point to a facet
858     dist is distance from point to facet
859     if dist NULL, 
860       searches for bestfacet and does nothing if inside
861     if qh.findbestnew set, 
862       searches new facets instead of using qh_findbest()
863
864   returns:
865     qh.max_ouside updated
866     if qh.KEEPcoplanar or qh.KEEPinside
867       point assigned to best coplanarset
868   
869   notes:
870     facet->maxoutside is updated at end by qh_check_maxout
871
872   design:
873     if dist undefined
874       find best facet for point
875       if point sufficiently below facet (depends on qh.NEARinside and qh.KEEPinside)
876         exit
877     if keeping coplanar/nearinside/inside points
878       if point is above furthest coplanar point
879         append point to coplanar set (it is the new furthest)
880         update qh.max_outside
881       else
882         append point one before end of coplanar set
883     else if point is clearly outside of qh.max_outside and bestfacet->coplanarset
884     and bestfacet is more than perpendicular to facet
885       repartition the point using qh_findbest() -- it may be put on an outsideset
886     else
887       update qh.max_outside
888 */
889 void qh_partitioncoplanar (pointT *point, facetT *facet, realT *dist) {
890   facetT *bestfacet;
891   pointT *oldfurthest;
892   realT bestdist, dist2, angle;
893   int numpart= 0, oldfindbest;
894   boolT isoutside;
895
896   qh WAScoplanar= True;
897   if (!dist) {
898     if (qh findbestnew)
899       bestfacet= qh_findbestnew (point, facet, &bestdist, qh_ALL, &isoutside, &numpart);
900     else
901       bestfacet= qh_findbest (point, facet, qh_ALL, !qh_ISnewfacets, qh DELAUNAY, 
902                           &bestdist, &isoutside, &numpart);
903     zinc_(Ztotpartcoplanar);
904     zzadd_(Zpartcoplanar, numpart);
905     if (!qh DELAUNAY && !qh KEEPinside) { /*  for 'd', bestdist skips upperDelaunay facets */
906       if (qh KEEPnearinside) {
907         if (bestdist < -qh NEARinside) { 
908           zinc_(Zcoplanarinside);
909           trace4((qh ferr, "qh_partitioncoplanar: point p%d is more than near-inside facet f%d dist %2.2g findbestnew %d\n",
910                   qh_pointid(point), bestfacet->id, bestdist, qh findbestnew));
911           return;
912         }
913       }else if (bestdist < -qh MAXcoplanar) {
914           trace4((qh ferr, "qh_partitioncoplanar: point p%d is inside facet f%d dist %2.2g findbestnew %d\n",
915                   qh_pointid(point), bestfacet->id, bestdist, qh findbestnew));
916         zinc_(Zcoplanarinside);
917         return;
918       }
919     }
920   }else {
921     bestfacet= facet;
922     bestdist= *dist;
923   }
924   if (bestdist > qh max_outside) {
925     if (!dist && facet != bestfacet) { 
926       zinc_(Zpartangle);
927       angle= qh_getangle(facet->normal, bestfacet->normal);
928       if (angle < 0) {
929         /* typically due to deleted vertex and coplanar facets, e.g.,
930              RBOX 1000 s Z1 G1e-13 t1001185205 | QHULL Tv */
931         zinc_(Zpartflip);
932         trace2((qh ferr, "qh_partitioncoplanar: repartition point p%d from f%d.  It is above flipped facet f%d dist %2.2g\n",
933                 qh_pointid(point), facet->id, bestfacet->id, bestdist));
934         oldfindbest= qh findbestnew;
935         qh findbestnew= False;
936         qh_partitionpoint(point, bestfacet);
937         qh findbestnew= oldfindbest;
938         return;
939       }
940     }
941     qh max_outside= bestdist;
942     if (bestdist > qh TRACEdist) {
943       fprintf (qh ferr, "qh_partitioncoplanar: ====== p%d from f%d increases max_outside to %2.2g of f%d last p%d\n",
944                      qh_pointid(point), facet->id, bestdist, bestfacet->id, qh furthest_id);
945       qh_errprint ("DISTANT", facet, bestfacet, NULL, NULL);
946     }
947   }
948   if (qh KEEPcoplanar + qh KEEPinside + qh KEEPnearinside) {
949     oldfurthest= (pointT*)qh_setlast (bestfacet->coplanarset);
950     if (oldfurthest) {
951       zinc_(Zcomputefurthest);
952       qh_distplane (oldfurthest, bestfacet, &dist2);
953     }
954     if (!oldfurthest || dist2 < bestdist) 
955       qh_setappend(&bestfacet->coplanarset, point);
956     else 
957       qh_setappend2ndlast(&bestfacet->coplanarset, point);
958   }
959   trace4((qh ferr, "qh_partitioncoplanar: point p%d is coplanar with facet f%d (or inside) dist %2.2g\n",
960           qh_pointid(point), bestfacet->id, bestdist));
961 } /* partitioncoplanar */
962
963 /*-<a                             href="qh-qhull.htm#TOC"
964   >-------------------------------</a><a name="partitionpoint">-</a>
965   
966   qh_partitionpoint( point, facet )
967     assigns point to an outside set, coplanar set, or inside set (i.e., dropt)
968     if qh.findbestnew
969       uses qh_findbestnew() to search all new facets
970     else
971       uses qh_findbest()
972   
973   notes:
974     after qh_distplane(), this and qh_findbest() are most expensive in 3-d
975
976   design:
977     find best facet for point 
978       (either exhaustive search of new facets or directed search from facet)
979     if qh.NARROWhull
980       retain coplanar and nearinside points as outside points
981     if point is outside bestfacet
982       if point above furthest point for bestfacet
983         append point to outside set (it becomes the new furthest)
984         if outside set was empty
985           move bestfacet to end of qh.facet_list (i.e., after qh.facet_next)
986         update bestfacet->furthestdist
987       else
988         append point one before end of outside set
989     else if point is coplanar to bestfacet
990       if keeping coplanar points or need to update qh.max_outside
991         partition coplanar point into bestfacet
992     else if near-inside point        
993       partition as coplanar point into bestfacet
994     else is an inside point
995       if keeping inside points 
996         partition as coplanar point into bestfacet
997 */
998 void qh_partitionpoint (pointT *point, facetT *facet) {
999   realT bestdist;
1000   boolT isoutside;
1001   facetT *bestfacet;
1002   int numpart;
1003 #if qh_COMPUTEfurthest
1004   realT dist;
1005 #endif
1006
1007   if (qh findbestnew)
1008     bestfacet= qh_findbestnew (point, facet, &bestdist, qh BESToutside, &isoutside, &numpart);
1009   else
1010     bestfacet= qh_findbest (point, facet, qh BESToutside, qh_ISnewfacets, !qh_NOupper,
1011                           &bestdist, &isoutside, &numpart);
1012   zinc_(Ztotpartition);
1013   zzadd_(Zpartition, numpart);
1014   if (qh NARROWhull) {
1015     if (qh DELAUNAY && !isoutside && bestdist >= -qh MAXcoplanar)
1016       qh_precision ("nearly incident point (narrow hull)");
1017     if (qh KEEPnearinside) {
1018       if (bestdist >= -qh NEARinside)
1019         isoutside= True;
1020     }else if (bestdist >= -qh MAXcoplanar)
1021       isoutside= True;
1022   }
1023
1024   if (isoutside) {
1025     if (!bestfacet->outsideset 
1026     || !qh_setlast (bestfacet->outsideset)) {
1027       qh_setappend(&(bestfacet->outsideset), point);
1028       if (!bestfacet->newfacet) {
1029         qh_removefacet (bestfacet);  /* make sure it's after qh facet_next */
1030         qh_appendfacet (bestfacet);
1031       }
1032 #if !qh_COMPUTEfurthest
1033       bestfacet->furthestdist= bestdist;
1034 #endif
1035     }else {
1036 #if qh_COMPUTEfurthest
1037       zinc_(Zcomputefurthest);
1038       qh_distplane (oldfurthest, bestfacet, &dist);
1039       if (dist < bestdist) 
1040         qh_setappend(&(bestfacet->outsideset), point);
1041       else
1042         qh_setappend2ndlast(&(bestfacet->outsideset), point);
1043 #else
1044       if (bestfacet->furthestdist < bestdist) {
1045         qh_setappend(&(bestfacet->outsideset), point);
1046         bestfacet->furthestdist= bestdist;
1047       }else
1048         qh_setappend2ndlast(&(bestfacet->outsideset), point);
1049 #endif
1050     }
1051     qh num_outside++;
1052     trace4((qh ferr, "qh_partitionpoint: point p%d is outside facet f%d new? %d(or narrowhull)\n",
1053           qh_pointid(point), bestfacet->id, bestfacet->newfacet));
1054   }else if (qh DELAUNAY || bestdist >= -qh MAXcoplanar) { /* for 'd', bestdist skips upperDelaunay facets */
1055     zzinc_(Zcoplanarpart);
1056     if (qh DELAUNAY)
1057       qh_precision ("nearly incident point");
1058     if ((qh KEEPcoplanar + qh KEEPnearinside) || bestdist > qh max_outside) 
1059       qh_partitioncoplanar (point, bestfacet, &bestdist);
1060     else {
1061       trace4((qh ferr, "qh_partitionpoint: point p%d is coplanar to facet f%d (dropped)\n",
1062           qh_pointid(point), bestfacet->id));
1063     }
1064   }else if (qh KEEPnearinside && bestdist > -qh NEARinside) {
1065     zinc_(Zpartnear);
1066     qh_partitioncoplanar (point, bestfacet, &bestdist);
1067   }else {
1068     zinc_(Zpartinside);
1069     trace4((qh ferr, "qh_partitionpoint: point p%d is inside all facets, closest to f%d dist %2.2g\n",
1070           qh_pointid(point), bestfacet->id, bestdist));
1071     if (qh KEEPinside)
1072       qh_partitioncoplanar (point, bestfacet, &bestdist);
1073   }
1074 } /* partitionpoint */
1075
1076 /*-<a                             href="qh-qhull.htm#TOC"
1077   >-------------------------------</a><a name="partitionvisible">-</a>
1078   
1079   qh_partitionvisible( allpoints, numoutside )
1080     partitions points in visible facets to qh.newfacet_list
1081     qh.visible_list= visible facets
1082     for visible facets
1083       1st neighbor (if any) points to a horizon facet or a new facet
1084     if allpoints (not used),
1085       repartitions coplanar points
1086
1087   returns:
1088     updates outside sets and coplanar sets of qh.newfacet_list
1089     updates qh.num_outside (count of outside points)
1090   
1091   notes:
1092     qh.findbest_notsharp should be clear (extra work if set)
1093
1094   design:
1095     for all visible facets with outside set or coplanar set
1096       select a newfacet for visible facet
1097       if outside set
1098         partition outside set into new facets
1099       if coplanar set and keeping coplanar/near-inside/inside points
1100         if allpoints
1101           partition coplanar set into new facets, may be assigned outside
1102         else
1103           partition coplanar set into coplanar sets of new facets
1104     for each deleted vertex
1105       if allpoints
1106         partition vertex into new facets, may be assigned outside
1107       else
1108         partition vertex into coplanar sets of new facets
1109 */
1110 void qh_partitionvisible(/*visible_list*/ boolT allpoints, int *numoutside) {
1111   facetT *visible, *newfacet;
1112   pointT *point, **pointp;
1113   int coplanar=0, size;
1114   unsigned count;
1115   vertexT *vertex, **vertexp;
1116   
1117   if (qh ONLYmax)
1118     maximize_(qh MINoutside, qh max_vertex);
1119   *numoutside= 0;
1120   FORALLvisible_facets {
1121     if (!visible->outsideset && !visible->coplanarset)
1122       continue;
1123     newfacet= visible->f.replace;
1124     count= 0;
1125     while (newfacet && newfacet->visible) {
1126       newfacet= newfacet->f.replace;
1127       if (count++ > qh facet_id)
1128         qh_infiniteloop (visible);
1129     }
1130     if (!newfacet)
1131       newfacet= qh newfacet_list;
1132     if (newfacet == qh facet_tail) {
1133       fprintf (qh ferr, "qhull precision error (qh_partitionvisible): all new facets deleted as\n        degenerate facets. Can not continue.\n");
1134       qh_errexit (qh_ERRprec, NULL, NULL);
1135     }
1136     if (visible->outsideset) {
1137       size= qh_setsize (visible->outsideset);
1138       *numoutside += size;
1139       qh num_outside -= size;
1140       FOREACHpoint_(visible->outsideset) 
1141         qh_partitionpoint (point, newfacet);
1142     }
1143     if (visible->coplanarset && (qh KEEPcoplanar + qh KEEPinside + qh KEEPnearinside)) {
1144       size= qh_setsize (visible->coplanarset);
1145       coplanar += size;
1146       FOREACHpoint_(visible->coplanarset) {
1147         if (allpoints) /* not used */
1148           qh_partitionpoint (point, newfacet);
1149         else
1150           qh_partitioncoplanar (point, newfacet, NULL);
1151       }
1152     }
1153   }
1154   FOREACHvertex_(qh del_vertices) {
1155     if (vertex->point) {
1156       if (allpoints) /* not used */
1157         qh_partitionpoint (vertex->point, qh newfacet_list);
1158       else
1159         qh_partitioncoplanar (vertex->point, qh newfacet_list, NULL);
1160     }
1161   }
1162   trace1((qh ferr,"qh_partitionvisible: partitioned %d points from outsidesets and %d points from coplanarsets\n", *numoutside, coplanar));
1163 } /* partitionvisible */
1164
1165
1166
1167 /*-<a                             href="qh-qhull.htm#TOC"
1168   >-------------------------------</a><a name="precision">-</a>
1169   
1170   qh_precision( reason )
1171     restart on precision errors if not merging and if 'QJn'
1172 */
1173 void qh_precision (char *reason) {
1174
1175   if (qh ALLOWrestart && !qh PREmerge && !qh MERGEexact) {
1176     if (qh JOGGLEmax < REALmax/2) {
1177       trace0((qh ferr, "qh_precision: qhull restart because of %s\n", reason));
1178       longjmp(qh restartexit, qh_ERRprec);
1179     }
1180   }
1181 } /* qh_precision */
1182
1183 /*-<a                             href="qh-qhull.htm#TOC"
1184   >-------------------------------</a><a name="printsummary">-</a>
1185   
1186   qh_printsummary( fp )
1187     prints summary to fp
1188
1189   notes:
1190     not in io.c so that user_eg.c can prevent io.c from loading
1191     qh_printsummary and qh_countfacets must match counts
1192
1193   design:
1194     determine number of points, vertices, and coplanar points
1195     print summary
1196 */
1197 void qh_printsummary(FILE *fp) {
1198   realT ratio, outerplane, innerplane;
1199   float cpu;
1200   int size, id, nummerged, numvertices, numcoplanars= 0, nonsimplicial=0;
1201   int goodused;
1202   facetT *facet;
1203   char *s;
1204   int numdel= zzval_(Zdelvertextot);
1205   int numtricoplanars= 0;
1206
1207   size= qh num_points + qh_setsize (qh other_points);
1208   numvertices= qh num_vertices - qh_setsize (qh del_vertices);
1209   id= qh_pointid (qh GOODpointp);
1210   FORALLfacets {
1211     if (facet->coplanarset)
1212       numcoplanars += qh_setsize( facet->coplanarset);
1213     if (facet->good) {
1214       if (facet->simplicial) {
1215         if (facet->keepcentrum && facet->tricoplanar)
1216           numtricoplanars++;
1217       }else if (qh_setsize(facet->vertices) != qh hull_dim)
1218         nonsimplicial++;
1219     }
1220   }
1221   if (id >=0 && qh STOPcone-1 != id && -qh STOPpoint-1 != id)
1222     size--;
1223   if (qh STOPcone || qh STOPpoint)
1224       fprintf (fp, "\nAt a premature exit due to 'TVn', 'TCn', 'TRn', or precision error.");
1225   if (qh UPPERdelaunay)
1226     goodused= qh GOODvertex + qh GOODpoint + qh SPLITthresholds;
1227   else if (qh DELAUNAY)
1228     goodused= qh GOODvertex + qh GOODpoint + qh GOODthreshold;
1229   else
1230     goodused= qh num_good;
1231   nummerged= zzval_(Ztotmerge) - zzval_(Zcyclehorizon) + zzval_(Zcyclefacettot);
1232   if (qh VORONOI) {
1233     if (qh UPPERdelaunay)
1234       fprintf (fp, "\n\
1235 Furthest-site Voronoi vertices by the convex hull of %d points in %d-d:\n\n", size, qh hull_dim);
1236     else
1237       fprintf (fp, "\n\
1238 Voronoi diagram by the convex hull of %d points in %d-d:\n\n", size, qh hull_dim);
1239     fprintf(fp, "  Number of Voronoi regions%s: %d\n",
1240               qh ATinfinity ? " and at-infinity" : "", numvertices);
1241     if (numdel)
1242       fprintf(fp, "  Total number of deleted points due to merging: %d\n", numdel); 
1243     if (numcoplanars - numdel > 0) 
1244       fprintf(fp, "  Number of nearly incident points: %d\n", numcoplanars - numdel); 
1245     else if (size - numvertices - numdel > 0) 
1246       fprintf(fp, "  Total number of nearly incident points: %d\n", size - numvertices - numdel); 
1247     fprintf(fp, "  Number of%s Voronoi vertices: %d\n", 
1248               goodused ? " 'good'" : "", qh num_good);
1249     if (nonsimplicial) 
1250       fprintf(fp, "  Number of%s non-simplicial Voronoi vertices: %d\n", 
1251               goodused ? " 'good'" : "", nonsimplicial);
1252   }else if (qh DELAUNAY) {
1253     if (qh UPPERdelaunay)
1254       fprintf (fp, "\n\
1255 Furthest-site Delaunay triangulation by the convex hull of %d points in %d-d:\n\n", size, qh hull_dim);
1256     else
1257       fprintf (fp, "\n\
1258 Delaunay triangulation by the convex hull of %d points in %d-d:\n\n", size, qh hull_dim);
1259     fprintf(fp, "  Number of input sites%s: %d\n", 
1260               qh ATinfinity ? " and at-infinity" : "", numvertices);
1261     if (numdel)
1262       fprintf(fp, "  Total number of deleted points due to merging: %d\n", numdel); 
1263     if (numcoplanars - numdel > 0) 
1264       fprintf(fp, "  Number of nearly incident points: %d\n", numcoplanars - numdel); 
1265     else if (size - numvertices - numdel > 0)
1266       fprintf(fp, "  Total number of nearly incident points: %d\n", size - numvertices - numdel); 
1267     fprintf(fp, "  Number of%s Delaunay regions: %d\n", 
1268               goodused ? " 'good'" : "", qh num_good);
1269     if (nonsimplicial) 
1270       fprintf(fp, "  Number of%s non-simplicial Delaunay regions: %d\n", 
1271               goodused ? " 'good'" : "", nonsimplicial);
1272   }else if (qh HALFspace) {
1273     fprintf (fp, "\n\
1274 Halfspace intersection by the convex hull of %d points in %d-d:\n\n", size, qh hull_dim);
1275     fprintf(fp, "  Number of halfspaces: %d\n", size);
1276     fprintf(fp, "  Number of non-redundant halfspaces: %d\n", numvertices);
1277     if (numcoplanars) {
1278       if (qh KEEPinside && qh KEEPcoplanar)
1279         s= "similar and redundant";
1280       else if (qh KEEPinside)
1281         s= "redundant";
1282       else
1283         s= "similar"; 
1284       fprintf(fp, "  Number of %s halfspaces: %d\n", s, numcoplanars);
1285     } 
1286     fprintf(fp, "  Number of intersection points: %d\n", qh num_facets - qh num_visible);
1287     if (goodused)
1288       fprintf(fp, "  Number of 'good' intersection points: %d\n", qh num_good);
1289     if (nonsimplicial) 
1290       fprintf(fp, "  Number of%s non-simplicial intersection points: %d\n", 
1291               goodused ? " 'good'" : "", nonsimplicial);
1292   }else {
1293     fprintf (fp, "\n\
1294 Convex hull of %d points in %d-d:\n\n", size, qh hull_dim);
1295     fprintf(fp, "  Number of vertices: %d\n", numvertices);
1296     if (numcoplanars) {
1297       if (qh KEEPinside && qh KEEPcoplanar)
1298         s= "coplanar and interior";
1299       else if (qh KEEPinside)
1300         s= "interior";
1301       else
1302         s= "coplanar"; 
1303       fprintf(fp, "  Number of %s points: %d\n", s, numcoplanars);
1304     } 
1305     fprintf(fp, "  Number of facets: %d\n", qh num_facets - qh num_visible);
1306     if (goodused)
1307       fprintf(fp, "  Number of 'good' facets: %d\n", qh num_good);
1308     if (nonsimplicial) 
1309       fprintf(fp, "  Number of%s non-simplicial facets: %d\n", 
1310               goodused ? " 'good'" : "", nonsimplicial);
1311   }
1312   if (numtricoplanars)
1313       fprintf(fp, "  Number of triangulated facets: %d\n", numtricoplanars);
1314   fprintf(fp, "\nStatistics for: %s | %s", 
1315                       qh rbox_command, qh qhull_command);
1316   if (qh ROTATErandom != INT_MIN)
1317     fprintf(fp, " QR%d\n\n", qh ROTATErandom);
1318   else
1319     fprintf(fp, "\n\n");
1320   fprintf(fp, "  Number of points processed: %d\n", zzval_(Zprocessed));
1321   fprintf(fp, "  Number of hyperplanes created: %d\n", zzval_(Zsetplane));
1322   if (qh DELAUNAY)
1323     fprintf(fp, "  Number of facets in hull: %d\n", qh num_facets - qh num_visible);
1324   fprintf(fp, "  Number of distance tests for qhull: %d\n", zzval_(Zpartition)+
1325       zzval_(Zpartitionall)+zzval_(Znumvisibility)+zzval_(Zpartcoplanar));
1326 #if 0  /* NOTE: must print before printstatistics() */
1327   {realT stddev, ave;
1328   fprintf(fp, "  average new facet balance: %2.2g\n",
1329           wval_(Wnewbalance)/zval_(Zprocessed));
1330   stddev= qh_stddev (zval_(Zprocessed), wval_(Wnewbalance), 
1331                                  wval_(Wnewbalance2), &ave);
1332   fprintf(fp, "  new facet standard deviation: %2.2g\n", stddev);
1333   fprintf(fp, "  average partition balance: %2.2g\n",
1334           wval_(Wpbalance)/zval_(Zpbalance));
1335   stddev= qh_stddev (zval_(Zpbalance), wval_(Wpbalance), 
1336                                  wval_(Wpbalance2), &ave);
1337   fprintf(fp, "  partition standard deviation: %2.2g\n", stddev);
1338   }
1339 #endif
1340   if (nummerged) {
1341     fprintf(fp,"  Number of distance tests for merging: %d\n",zzval_(Zbestdist)+
1342           zzval_(Zcentrumtests)+zzval_(Zdistconvex)+zzval_(Zdistcheck)+
1343           zzval_(Zdistzero));
1344     fprintf(fp,"  Number of distance tests for checking: %d\n",zzval_(Zcheckpart));
1345     fprintf(fp,"  Number of merged facets: %d\n", nummerged);
1346   }
1347   if (!qh RANDOMoutside && qh QHULLfinished) {
1348     cpu= qh hulltime;
1349     cpu /= qh_SECticks;
1350     wval_(Wcpu)= cpu;
1351     fprintf (fp, "  CPU seconds to compute hull (after input): %2.4g\n", cpu);
1352   }
1353   if (qh RERUN) {
1354     if (!qh PREmerge && !qh MERGEexact)
1355       fprintf(fp, "  Percentage of runs with precision errors: %4.1f\n",
1356            zzval_(Zretry)*100.0/qh build_cnt);  /* careful of order */
1357   }else if (qh JOGGLEmax < REALmax/2) {
1358     if (zzval_(Zretry))
1359       fprintf(fp, "  After %d retries, input joggled by: %2.2g\n",
1360          zzval_(Zretry), qh JOGGLEmax);
1361     else
1362       fprintf(fp, "  Input joggled by: %2.2g\n", qh JOGGLEmax);
1363   }
1364   if (qh totarea != 0.0) 
1365     fprintf(fp, "  %s facet area:   %2.8g\n", 
1366             zzval_(Ztotmerge) ? "Approximate" : "Total", qh totarea);
1367   if (qh totvol != 0.0) 
1368     fprintf(fp, "  %s volume:       %2.8g\n", 
1369             zzval_(Ztotmerge) ? "Approximate" : "Total", qh totvol);
1370   if (qh MERGING) {
1371     qh_outerinner (NULL, &outerplane, &innerplane);
1372     if (outerplane > 2 * qh DISTround) {
1373       fprintf(fp, "  Maximum distance of %spoint above facet: %2.2g", 
1374             (qh QHULLfinished ? "" : "merged "), outerplane);
1375       ratio= outerplane/(qh ONEmerge + qh DISTround);
1376       /* don't report ratio if MINoutside is large */
1377       if (ratio > 0.05 && 2* qh ONEmerge > qh MINoutside && qh JOGGLEmax > REALmax/2)
1378         fprintf (fp, " (%.1fx)\n", ratio);
1379       else
1380         fprintf (fp, "\n");
1381     }
1382     if (innerplane < -2 * qh DISTround) {
1383       fprintf(fp, "  Maximum distance of %svertex below facet: %2.2g", 
1384             (qh QHULLfinished ? "" : "merged "), innerplane);
1385       ratio= -innerplane/(qh ONEmerge+qh DISTround);
1386       if (ratio > 0.05 && qh JOGGLEmax > REALmax/2)
1387         fprintf (fp, " (%.1fx)\n", ratio);
1388       else
1389         fprintf (fp, "\n");
1390     }
1391   }
1392   fprintf(fp, "\n");
1393 } /* printsummary */
1394
1395