Free cookie consent management tool by TermsFeed Policy Generator

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

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

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