Context Navigation

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

Visit:

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

Line
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
60	void 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	*/
168	boolT 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	*/
289	void 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	*/
363	void 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	*/
434	void 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\
451	At %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\
475	At %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	*/
520	void 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	*/
555	void 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\
629	QhullPoint 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	*/
661	pointT 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	*/
776	void 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	*/
895	void 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	*/
1004	void 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	*/
1116	void 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	*/
1179	void 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	*/
1203	void 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\
1241	Furthest-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\
1244	Voronoi 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\
1261	Furthest-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\
1264	Delaunay 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\
1280	Halfspace 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\
1300	Convex 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.

Download in other formats:

Update cookies preferences