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