Line data Source code
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_r.h, unix_r.c
10 :
11 : see qhull_ra.h for internal functions
12 :
13 : Copyright (c) 1993-2020 The Geometry Center.
14 : $Id: //main/2019/qhull/src/libqhull_r/libqhull_r.c#17 $$Change: 2953 $
15 : $DateTime: 2020/05/21 22:05:32 $$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 :
60 4 : void qh_qhull(qhT *qh) {
61 : int numoutside;
62 :
63 4 : qh->hulltime= qh_CPUclock;
64 4 : if (qh->RERUN || qh->JOGGLEmax < REALmax/2)
65 0 : qh_build_withrestart(qh);
66 : else {
67 4 : qh_initbuild(qh);
68 3 : qh_buildhull(qh);
69 : }
70 3 : if (!qh->STOPadd && !qh->STOPcone && !qh->STOPpoint) {
71 3 : if (qh->ZEROall_ok && !qh->TESTvneighbors && qh->MERGEexact)
72 0 : qh_checkzero(qh, qh_ALL);
73 3 : if (qh->ZEROall_ok && !qh->TESTvneighbors && !qh->WAScoplanar) {
74 0 : 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 0 : qh->DOcheckmax= False;
76 : }else {
77 3 : qh_initmergesets(qh /* qh.facet_mergeset,degen_mergeset,vertex_mergeset */);
78 3 : if (qh->MERGEexact || (qh->hull_dim > qh_DIMreduceBuild && qh->PREmerge))
79 0 : qh_postmerge(qh, "First post-merge", qh->premerge_centrum, qh->premerge_cos,
80 0 : (qh->POSTmerge ? False : qh->TESTvneighbors)); /* calls qh_reducevertices */
81 3 : else if (!qh->POSTmerge && qh->TESTvneighbors)
82 0 : qh_postmerge(qh, "For testing vertex neighbors", qh->premerge_centrum,
83 : qh->premerge_cos, True); /* calls qh_test_vneighbors */
84 3 : if (qh->POSTmerge)
85 0 : qh_postmerge(qh, "For post-merging", qh->postmerge_centrum,
86 : qh->postmerge_cos, qh->TESTvneighbors);
87 3 : if (qh->visible_list == qh->facet_list) { /* qh_postmerge was called */
88 0 : qh->findbestnew= True;
89 0 : qh_partitionvisible(qh, !qh_ALL, &numoutside /* qh.visible_list */);
90 0 : qh->findbestnew= False;
91 0 : qh_deletevisible(qh /* qh.visible_list */); /* stops at first !f.visible */
92 0 : qh_resetlists(qh, False, qh_RESETvisible /* qh.visible_list newvertex_list qh.newfacet_list */);
93 : }
94 3 : qh_all_vertexmerges(qh, -1, NULL, NULL);
95 3 : qh_freemergesets(qh);
96 : }
97 3 : if (qh->TRACEpoint == qh_IDunknown && qh->TRACElevel > qh->IStracing) {
98 0 : qh->IStracing= qh->TRACElevel;
99 0 : qh_fprintf(qh, qh->ferr, 2112, "qh_qhull: finished qh_buildhull and qh_postmerge, start tracing (TP-1)\n");
100 : }
101 3 : if (qh->DOcheckmax){
102 3 : if (qh->REPORTfreq) {
103 0 : qh_buildtracing(qh, NULL, NULL);
104 0 : qh_fprintf(qh, qh->ferr, 8115, "\nTesting all coplanar points.\n");
105 : }
106 3 : qh_check_maxout(qh);
107 : }
108 3 : if (qh->KEEPnearinside && !qh->maxoutdone)
109 0 : qh_nearcoplanar(qh);
110 : }
111 3 : if (qh_setsize(qh, qh->qhmem.tempstack) != 0) {
112 0 : qh_fprintf(qh, qh->ferr, 6164, "qhull internal error (qh_qhull): temporary sets not empty(%d) at end of Qhull\n",
113 : qh_setsize(qh, qh->qhmem.tempstack));
114 0 : qh_errexit(qh, qh_ERRqhull, NULL, NULL);
115 : }
116 3 : qh->hulltime= qh_CPUclock - qh->hulltime;
117 3 : qh->QHULLfinished= True;
118 3 : trace1((qh, qh->ferr, 1036, "Qhull: algorithm completed\n"));
119 3 : } /* qhull */
120 :
121 : /*-<a href="qh-qhull_r.htm#TOC"
122 : >-------------------------------</a><a name="addpoint">-</a>
123 :
124 : qh_addpoint(qh, furthest, facet, checkdist )
125 : add point (usually furthest point) above facet to hull
126 : if checkdist,
127 : check that point is above facet.
128 : if point is not outside of the hull, uses qh_partitioncoplanar()
129 : assumes that facet is defined by qh_findbestfacet()
130 : else if facet specified,
131 : assumes that point is above facet (major damage if below)
132 : for Delaunay triangulations,
133 : Use qh_setdelaunay() to lift point to paraboloid and scale by 'Qbb' if needed
134 : Do not use options 'Qbk', 'QBk', or 'QbB' since they scale the coordinates.
135 :
136 : returns:
137 : returns False if user requested an early termination
138 : qh.visible_list, newfacet_list, delvertex_list, NEWfacets may be defined
139 : updates qh.facet_list, qh.num_facets, qh.vertex_list, qh.num_vertices
140 : clear qh.maxoutdone (will need to call qh_check_maxout() for facet->maxoutside)
141 : if unknown point, adds a pointer to qh.other_points
142 : do not deallocate the point's coordinates
143 :
144 : notes:
145 : called from qh_initbuild, qh_buildhull, and qh_addpoint
146 : tail recursive call if merged a pinchedvertex due to a duplicated ridge
147 : no more than qh.num_vertices calls (QH6296)
148 : assumes point is near its best facet and not at a local minimum of a lens
149 : distributions. Use qh_findbestfacet to avoid this case.
150 : uses qh.visible_list, qh.newfacet_list, qh.delvertex_list, qh.NEWfacets
151 : if called from a user application after qh_qhull and 'QJ' (joggle),
152 : facet merging for precision problems is disabled by default
153 :
154 : design:
155 : exit if qh.STOPadd vertices 'TAn'
156 : add point to other_points if needed
157 : if checkdist
158 : if point not above facet
159 : partition coplanar point
160 : exit
161 : exit if pre STOPpoint requested
162 : find horizon and visible facets for point
163 : build cone of new facets to the horizon
164 : exit if build cone fails due to qh.ONLYgood
165 : tail recursive call if build cone fails due to pinched vertices
166 : exit if STOPcone requested
167 : merge non-convex new facets
168 : if merge found, many merges, or 'Qf'
169 : use qh_findbestnew() instead of qh_findbest()
170 : partition outside points from visible facets
171 : delete visible facets
172 : check polyhedron if requested
173 : exit if post STOPpoint requested
174 : reset working lists of facets and vertices
175 : */
176 14660 : boolT qh_addpoint(qhT *qh, pointT *furthest, facetT *facet, boolT checkdist) {
177 : realT dist, pbalance;
178 : facetT *replacefacet, *newfacet;
179 : vertexT *apex;
180 14660 : boolT isoutside= False;
181 : int numpart, numpoints, goodvisible, goodhorizon, apexpointid;
182 :
183 14660 : qh->maxoutdone= False;
184 14660 : if (qh_pointid(qh, furthest) == qh_IDunknown)
185 0 : qh_setappend(qh, &qh->other_points, furthest);
186 14660 : if (!facet) {
187 0 : qh_fprintf(qh, qh->ferr, 6213, "qhull internal error (qh_addpoint): NULL facet. Need to call qh_findbestfacet first\n");
188 0 : qh_errexit(qh, qh_ERRqhull, NULL, NULL);
189 : }
190 14660 : qh_detmaxoutside(qh);
191 14660 : if (checkdist) {
192 0 : facet= qh_findbest(qh, furthest, facet, !qh_ALL, !qh_ISnewfacets, !qh_NOupper,
193 : &dist, &isoutside, &numpart);
194 0 : zzadd_(Zpartition, numpart);
195 0 : if (!isoutside) {
196 0 : zinc_(Znotmax); /* last point of outsideset is no longer furthest. */
197 0 : facet->notfurthest= True;
198 0 : qh_partitioncoplanar(qh, furthest, facet, &dist, qh->findbestnew);
199 0 : return True;
200 : }
201 : }
202 14660 : qh_buildtracing(qh, furthest, facet);
203 14660 : if (qh->STOPpoint < 0 && qh->furthest_id == -qh->STOPpoint-1) {
204 0 : facet->notfurthest= True;
205 0 : return False;
206 : }
207 14660 : qh_findhorizon(qh, furthest, facet, &goodvisible, &goodhorizon);
208 14660 : if (qh->ONLYgood && !qh->GOODclosest && !(goodvisible+goodhorizon)) {
209 0 : zinc_(Znotgood);
210 0 : facet->notfurthest= True;
211 : /* last point of outsideset is no longer furthest. This is ok
212 : since all points of the outside are likely to be bad */
213 0 : qh_resetlists(qh, False, qh_RESETvisible /* qh.visible_list newvertex_list qh.newfacet_list */);
214 0 : return True;
215 : }
216 14660 : apex= qh_buildcone(qh, furthest, facet, goodhorizon, &replacefacet);
217 : /* qh.newfacet_list, visible_list, newvertex_list */
218 14660 : if (!apex) {
219 0 : if (qh->ONLYgood)
220 0 : return True; /* ignore this furthest point, a good new facet was not found */
221 0 : if (replacefacet) {
222 0 : if (qh->retry_addpoint++ >= qh->num_vertices) {
223 0 : qh_fprintf(qh, qh->ferr, 6296, "qhull internal error (qh_addpoint): infinite loop (%d retries) of merging pinched vertices due to dupridge for point p%d, facet f%d, and %d vertices\n",
224 : qh->retry_addpoint, qh_pointid(qh, furthest), facet->id, qh->num_vertices);
225 0 : qh_errexit(qh, qh_ERRqhull, facet, NULL);
226 : }
227 : /* retry qh_addpoint after resolving a dupridge via qh_merge_pinchedvertices */
228 0 : return qh_addpoint(qh, furthest, replacefacet, True /* checkdisk */);
229 : }
230 0 : qh->retry_addpoint= 0;
231 0 : return True; /* ignore this furthest point, resolved a dupridge by making furthest a coplanar point */
232 : }
233 14660 : if (qh->retry_addpoint) {
234 0 : zinc_(Zretryadd);
235 0 : zadd_(Zretryaddtot, qh->retry_addpoint);
236 0 : zmax_(Zretryaddmax, qh->retry_addpoint);
237 0 : qh->retry_addpoint= 0;
238 : }
239 14660 : apexpointid= qh_pointid(qh, apex->point);
240 14660 : zzinc_(Zprocessed);
241 14660 : if (qh->STOPcone && qh->furthest_id == qh->STOPcone-1) {
242 0 : facet->notfurthest= True;
243 0 : return False; /* visible_list etc. still defined */
244 : }
245 14660 : qh->findbestnew= False;
246 14660 : if (qh->PREmerge || qh->MERGEexact) {
247 14660 : qh_initmergesets(qh /* qh.facet_mergeset,degen_mergeset,vertex_mergeset */);
248 14660 : qh_premerge(qh, apexpointid, qh->premerge_centrum, qh->premerge_cos /* qh.newfacet_list */);
249 14660 : if (qh_USEfindbestnew)
250 14603 : qh->findbestnew= True;
251 : else {
252 237 : FORALLnew_facets {
253 210 : if (!newfacet->simplicial) {
254 30 : qh->findbestnew= True; /* use qh_findbestnew instead of qh_findbest*/
255 30 : break;
256 : }
257 : }
258 : }
259 0 : }else if (qh->BESToutside)
260 0 : qh->findbestnew= True;
261 14660 : if (qh->IStracing >= 4)
262 0 : qh_checkpolygon(qh, qh->visible_list);
263 14660 : qh_partitionvisible(qh, !qh_ALL, &numpoints /* qh.visible_list */);
264 14660 : qh->findbestnew= False;
265 14660 : qh->findbest_notsharp= False;
266 14660 : zinc_(Zpbalance);
267 14660 : pbalance= numpoints - (realT) qh->hull_dim /* assumes all points extreme */
268 14660 : * (qh->num_points - qh->num_vertices)/qh->num_vertices;
269 14660 : wadd_(Wpbalance, pbalance);
270 14660 : wadd_(Wpbalance2, pbalance * pbalance);
271 14660 : qh_deletevisible(qh /* qh.visible_list */);
272 14660 : zmax_(Zmaxvertex, qh->num_vertices);
273 14660 : qh->NEWfacets= False;
274 14660 : if (qh->IStracing >= 4) {
275 0 : if (qh->num_facets < 200)
276 0 : qh_printlists(qh);
277 0 : qh_printfacetlist(qh, qh->newfacet_list, NULL, True);
278 0 : qh_checkpolygon(qh, qh->facet_list);
279 14660 : }else if (qh->CHECKfrequently) {
280 0 : if (qh->num_facets < 1000)
281 0 : qh_checkpolygon(qh, qh->facet_list);
282 : else
283 0 : qh_checkpolygon(qh, qh->newfacet_list);
284 : }
285 14660 : if (qh->STOPpoint > 0 && qh->furthest_id == qh->STOPpoint-1 && qh_setsize(qh, qh->vertex_mergeset) > 0)
286 0 : return False;
287 14660 : qh_resetlists(qh, True, qh_RESETvisible /* qh.visible_list newvertex_list qh.newfacet_list */);
288 14660 : if (qh->facet_mergeset) {
289 : /* vertex merges occur after facet merges (qh_premerge) and qh_resetlists */
290 14660 : qh_all_vertexmerges(qh, apexpointid, NULL, NULL);
291 14660 : qh_freemergesets(qh);
292 : }
293 : /* qh_triangulate(qh); to test qh.TRInormals */
294 14660 : if (qh->STOPpoint > 0 && qh->furthest_id == qh->STOPpoint-1)
295 0 : return False;
296 14660 : trace2((qh, qh->ferr, 2056, "qh_addpoint: added p%d to convex hull with point balance %2.2g\n",
297 : qh_pointid(qh, furthest), pbalance));
298 14660 : return True;
299 : } /* addpoint */
300 :
301 : /*-<a href="qh-qhull_r.htm#TOC"
302 : >-------------------------------</a><a name="build_withrestart">-</a>
303 :
304 : qh_build_withrestart(qh)
305 : allow restarts due to qh.JOGGLEmax while calling qh_buildhull()
306 : qh_errexit always undoes qh_build_withrestart()
307 : qh.FIRSTpoint/qh.NUMpoints is point array
308 : it may be moved by qh_joggleinput
309 : */
310 0 : void qh_build_withrestart(qhT *qh) {
311 : int restart;
312 : vertexT *vertex, **vertexp;
313 :
314 0 : qh->ALLOWrestart= True;
315 : while (True) {
316 0 : restart= setjmp(qh->restartexit); /* simple statement for CRAY J916 */
317 0 : if (restart) { /* only from qh_joggle_restart() */
318 0 : qh->last_errcode= qh_ERRnone;
319 0 : zzinc_(Zretry);
320 0 : wmax_(Wretrymax, qh->JOGGLEmax);
321 : /* QH7078 warns about using 'TCn' with 'QJn' */
322 0 : qh->STOPcone= qh_IDunknown; /* if break from joggle, prevents normal output */
323 0 : FOREACHvertex_(qh->del_vertices) {
324 0 : if (vertex->point && !vertex->partitioned)
325 0 : vertex->partitioned= True; /* avoid error in qh_freebuild -> qh_delvertex */
326 : }
327 : }
328 0 : if (!qh->RERUN && qh->JOGGLEmax < REALmax/2) {
329 0 : if (qh->build_cnt > qh_JOGGLEmaxretry) {
330 0 : qh_fprintf(qh, qh->ferr, 6229, "qhull input error: %d attempts to construct a convex hull with joggled input. Increase joggle above 'QJ%2.2g' or modify qh_JOGGLE... parameters in user_r.h\n",
331 : qh->build_cnt, qh->JOGGLEmax);
332 0 : qh_errexit(qh, qh_ERRinput, NULL, NULL);
333 : }
334 0 : if (qh->build_cnt && !restart)
335 0 : break;
336 0 : }else if (qh->build_cnt && qh->build_cnt >= qh->RERUN)
337 0 : break;
338 0 : qh->STOPcone= 0;
339 0 : qh_freebuild(qh, True); /* first call is a nop */
340 0 : qh->build_cnt++;
341 0 : if (!qh->qhull_optionsiz)
342 0 : qh->qhull_optionsiz= (int)strlen(qh->qhull_options); /* WARN64 */
343 : else {
344 0 : qh->qhull_options[qh->qhull_optionsiz]= '\0';
345 0 : qh->qhull_optionlen= qh_OPTIONline; /* starts a new line */
346 : }
347 0 : qh_option(qh, "_run", &qh->build_cnt, NULL);
348 0 : if (qh->build_cnt == qh->RERUN) {
349 0 : qh->IStracing= qh->TRACElastrun; /* duplicated from qh_initqhull_globals */
350 0 : if (qh->TRACEpoint != qh_IDnone || qh->TRACEdist < REALmax/2 || qh->TRACEmerge) {
351 0 : qh->TRACElevel= (qh->IStracing? qh->IStracing : 3);
352 0 : qh->IStracing= 0;
353 : }
354 0 : qh->qhmem.IStracing= qh->IStracing;
355 : }
356 0 : if (qh->JOGGLEmax < REALmax/2)
357 0 : qh_joggleinput(qh);
358 0 : qh_initbuild(qh);
359 0 : qh_buildhull(qh);
360 0 : if (qh->JOGGLEmax < REALmax/2 && !qh->MERGING)
361 0 : qh_checkconvex(qh, qh->facet_list, qh_ALGORITHMfault);
362 : }
363 0 : qh->ALLOWrestart= False;
364 0 : } /* qh_build_withrestart */
365 :
366 : /*-<a href="qh-qhull_r.htm#TOC"
367 : >-------------------------------</a><a name="buildcone">-</a>
368 :
369 : qh_buildcone(qh, furthest, facet, goodhorizon, &replacefacet )
370 : build cone of new facets from furthest to the horizon
371 : goodhorizon is count of good, horizon facets from qh_find_horizon
372 :
373 : returns:
374 : returns apex of cone with qh.newfacet_list and qh.first_newfacet (f.id)
375 : returns NULL if qh.ONLYgood and no good facets
376 : returns NULL and retryfacet if merging pinched vertices will resolve a dupridge
377 : a horizon vertex was nearly adjacent to another vertex
378 : will retry qh_addpoint
379 : returns NULL if resolve a dupridge by making furthest a coplanar point
380 : furthest was nearly adjacent to an existing vertex
381 : updates qh.degen_mergeset (MRGridge) if resolve a dupridge by merging facets
382 : updates qh.newfacet_list, visible_list, newvertex_list
383 : updates qh.facet_list, vertex_list, num_facets, num_vertices
384 :
385 : notes:
386 : called by qh_addpoint
387 : see qh_triangulate, it triangulates non-simplicial facets in post-processing
388 :
389 : design:
390 : make new facets for point to horizon
391 : compute balance statistics
392 : make hyperplanes for point
393 : exit if qh.ONLYgood and not good (qh_buildcone_onlygood)
394 : match neighboring new facets
395 : if dupridges
396 : exit if !qh.IGNOREpinched and dupridge resolved by coplanar furthest
397 : retry qh_buildcone if !qh.IGNOREpinched and dupridge resolved by qh_buildcone_mergepinched
398 : otherwise dupridges resolved by merging facets
399 : update vertex neighbors and delete interior vertices
400 : */
401 14660 : vertexT *qh_buildcone(qhT *qh, pointT *furthest, facetT *facet, int goodhorizon, facetT **retryfacet) {
402 : vertexT *apex;
403 : realT newbalance;
404 : int numnew;
405 :
406 14660 : *retryfacet= NULL;
407 14660 : qh->first_newfacet= qh->facet_id;
408 14660 : qh->NEWtentative= (qh->MERGEpinched || qh->ONLYgood); /* cleared by qh_attachnewfacets or qh_resetlists */
409 14660 : apex= qh_makenewfacets(qh, furthest /* qh.newfacet_list visible_list, attaches new facets if !qh.NEWtentative */);
410 14660 : numnew= (int)(qh->facet_id - qh->first_newfacet);
411 14660 : newbalance= numnew - (realT)(qh->num_facets - qh->num_visible) * qh->hull_dim / qh->num_vertices;
412 : /* newbalance statistics updated below if the new facets are accepted */
413 14660 : if (qh->ONLYgood) { /* qh.MERGEpinched is false by QH6362 */
414 0 : if (!qh_buildcone_onlygood(qh, apex, goodhorizon /* qh.newfacet_list */)) {
415 0 : facet->notfurthest= True;
416 0 : return NULL;
417 : }
418 14660 : }else if(qh->MERGEpinched) {
419 : #ifndef qh_NOmerge
420 0 : if (qh_buildcone_mergepinched(qh, apex, facet, retryfacet /* qh.newfacet_list */))
421 0 : return NULL;
422 : #else
423 : qh_fprintf(qh, qh->ferr, 6375, "qhull option error (qh_buildcone): option 'Q14' (qh.MERGEpinched) is not available due to qh_NOmerge\n");
424 : qh_errexit(qh, qh_ERRinput, NULL, NULL);
425 : #endif
426 : }else {
427 : /* qh_makenewfacets attached new facets to the horizon */
428 14660 : qh_matchnewfacets(qh); /* ignore returned value. qh_forcedmerges will merge dupridges if any */
429 14660 : qh_makenewplanes(qh /* qh.newfacet_list */);
430 14660 : qh_update_vertexneighbors_cone(qh);
431 : }
432 14660 : wadd_(Wnewbalance, newbalance);
433 14660 : wadd_(Wnewbalance2, newbalance * newbalance);
434 14660 : trace2((qh, qh->ferr, 2067, "qh_buildcone: created %d newfacets for p%d(v%d) new facet balance %2.2g\n",
435 : numnew, qh_pointid(qh, furthest), apex->id, newbalance));
436 14660 : return apex;
437 : } /* buildcone */
438 :
439 : #ifndef qh_NOmerge
440 : /*-<a href="qh-qhull_r.htm#TOC"
441 : >-------------------------------</a><a name="buildcone_mergepinched">-</a>
442 :
443 : qh_buildcone_mergepinched(qh, apex, facet, maxdupdist, &retryfacet )
444 : build cone of new facets from furthest to the horizon
445 : maxdupdist>0.0 for merging dupridges (qh_matchdupridge)
446 :
447 : returns:
448 : returns True if merged a pinched vertex and deleted the cone of new facets
449 : if retryfacet is set
450 : a dupridge was resolved by qh_merge_pinchedvertices
451 : retry qh_addpoint
452 : otherwise
453 : apex/furthest was partitioned as a coplanar point
454 : ignore this furthest point
455 : returns False if no dupridges or if dupridges will be resolved by MRGridge
456 : updates qh.facet_list, qh.num_facets, qh.vertex_list, qh.num_vertices
457 :
458 : notes:
459 : only called from qh_buildcone with qh.MERGEpinched
460 :
461 : design:
462 : match neighboring new facets
463 : if matching detected dupridges with a wide merge (qh_RATIOtrypinched)
464 : if pinched vertices (i.e., nearly adjacent)
465 : delete the cone of new facets
466 : delete the apex and reset the facet lists
467 : if coplanar, pinched apex
468 : partition the apex as a coplanar point
469 : else
470 : repeatedly merge the nearest pair of pinched vertices and subsequent facet merges
471 : return True
472 : otherwise
473 : MRGridge are better than vertex merge, but may report an error
474 : attach new facets
475 : make hyperplanes for point
476 : update vertex neighbors and delete interior vertices
477 : */
478 0 : boolT qh_buildcone_mergepinched(qhT *qh, vertexT *apex, facetT *facet, facetT **retryfacet) {
479 : facetT *newfacet, *nextfacet;
480 : pointT *apexpoint;
481 : coordT maxdupdist;
482 : int apexpointid;
483 : boolT iscoplanar;
484 :
485 0 : *retryfacet= NULL;
486 0 : maxdupdist= qh_matchnewfacets(qh);
487 0 : if (maxdupdist > qh_RATIOtrypinched * qh->ONEmerge) { /* one or more dupridges with a wide merge */
488 0 : if (qh->IStracing >= 4 && qh->num_facets < 1000)
489 0 : qh_printlists(qh);
490 0 : qh_initmergesets(qh /* qh.facet_mergeset,degen_mergeset,vertex_mergeset */);
491 0 : if (qh_getpinchedmerges(qh, apex, maxdupdist, &iscoplanar /* qh.newfacet_list, qh.vertex_mergeset */)) {
492 0 : for (newfacet=qh->newfacet_list; newfacet && newfacet->next; newfacet= nextfacet) {
493 0 : nextfacet= newfacet->next;
494 0 : qh_delfacet(qh, newfacet);
495 : }
496 0 : apexpoint= apex->point;
497 0 : apexpointid= qh_pointid(qh, apexpoint);
498 0 : qh_delvertex(qh, apex);
499 0 : qh_resetlists(qh, False, qh_RESETvisible /* qh.visible_list newvertex_list qh.newfacet_list */);
500 0 : if (iscoplanar) {
501 0 : zinc_(Zpinchedapex);
502 0 : facet->notfurthest= True;
503 0 : qh_partitioncoplanar(qh, apexpoint, facet, NULL, qh->findbestnew);
504 : }else {
505 0 : qh_all_vertexmerges(qh, apexpointid, facet, retryfacet);
506 : }
507 0 : qh_freemergesets(qh); /* errors if not empty */
508 0 : return True;
509 : }
510 : /* MRGridge are better than vertex merge, but may report an error */
511 0 : qh_freemergesets(qh);
512 : }
513 0 : qh_attachnewfacets(qh /* qh.visible_list */);
514 0 : qh_makenewplanes(qh /* qh.newfacet_list */);
515 0 : qh_update_vertexneighbors_cone(qh);
516 0 : return False;
517 : } /* buildcone_mergepinched */
518 : #endif /* !qh_NOmerge */
519 :
520 : /*-<a href="qh-qhull_r.htm#TOC"
521 : >-------------------------------</a><a name="buildcone_onlygood">-</a>
522 :
523 : qh_buildcone_onlygood(qh, apex, goodhorizon )
524 : build cone of good, new facets from apex and its qh.newfacet_list to the horizon
525 : goodhorizon is count of good, horizon facets from qh_find_horizon
526 :
527 : returns:
528 : False if a f.good facet or a qh.GOODclosest facet is not found
529 : updates qh.facet_list, qh.num_facets, qh.vertex_list, qh.num_vertices
530 :
531 : notes:
532 : called from qh_buildcone
533 : QH11030 FIX: Review effect of qh.GOODclosest on qh_buildcone_onlygood ('Qg'). qh_findgood preserves old value if didn't find a good facet. See qh_findgood_all for disabling
534 :
535 : design:
536 : make hyperplanes for point
537 : if qh_findgood fails to find a f.good facet or a qh.GOODclosest facet
538 : delete cone of new facets
539 : return NULL (ignores apex)
540 : else
541 : attach cone to horizon
542 : match neighboring new facets
543 : */
544 0 : boolT qh_buildcone_onlygood(qhT *qh, vertexT *apex, int goodhorizon) {
545 : facetT *newfacet, *nextfacet;
546 :
547 0 : qh_makenewplanes(qh /* qh.newfacet_list */);
548 0 : if(qh_findgood(qh, qh->newfacet_list, goodhorizon) == 0) {
549 0 : if (!qh->GOODclosest) {
550 0 : for (newfacet=qh->newfacet_list; newfacet && newfacet->next; newfacet= nextfacet) {
551 0 : nextfacet= newfacet->next;
552 0 : qh_delfacet(qh, newfacet);
553 : }
554 0 : qh_delvertex(qh, apex);
555 0 : qh_resetlists(qh, False /*no stats*/, qh_RESETvisible /* qh.visible_list newvertex_list qh.newfacet_list */);
556 0 : zinc_(Znotgoodnew);
557 : /* !good outside points dropped from hull */
558 0 : return False;
559 : }
560 : }
561 0 : qh_attachnewfacets(qh /* qh.visible_list */);
562 0 : qh_matchnewfacets(qh); /* ignore returned value. qh_forcedmerges will merge dupridges if any */
563 0 : qh_update_vertexneighbors_cone(qh);
564 0 : return True;
565 : } /* buildcone_onlygood */
566 :
567 : /*-<a href="qh-qhull_r.htm#TOC"
568 : >-------------------------------</a><a name="buildhull">-</a>
569 :
570 : qh_buildhull(qh)
571 : construct a convex hull by adding outside points one at a time
572 :
573 : returns:
574 :
575 : notes:
576 : may be called multiple times
577 : checks facet and vertex lists for incorrect flags
578 : to recover from STOPcone, call qh_deletevisible and qh_resetlists
579 :
580 : design:
581 : check visible facet and newfacet flags
582 : check newfacet vertex flags and qh.STOPcone/STOPpoint
583 : for each facet with a furthest outside point
584 : add point to facet
585 : exit if qh.STOPcone or qh.STOPpoint requested
586 : if qh.NARROWhull for initial simplex
587 : partition remaining outside points to coplanar sets
588 : */
589 3 : void qh_buildhull(qhT *qh) {
590 : facetT *facet;
591 : pointT *furthest;
592 : vertexT *vertex;
593 : int id;
594 :
595 3 : trace1((qh, qh->ferr, 1037, "qh_buildhull: start build hull\n"));
596 15 : FORALLfacets {
597 12 : if (facet->visible || facet->newfacet) {
598 0 : qh_fprintf(qh, qh->ferr, 6165, "qhull internal error (qh_buildhull): visible or new facet f%d in facet list\n",
599 0 : facet->id);
600 0 : qh_errexit(qh, qh_ERRqhull, facet, NULL);
601 : }
602 : }
603 15 : FORALLvertices {
604 12 : if (vertex->newfacet) {
605 0 : qh_fprintf(qh, qh->ferr, 6166, "qhull internal error (qh_buildhull): new vertex f%d in vertex list\n",
606 : vertex->id);
607 0 : qh_errprint(qh, "ERRONEOUS", NULL, NULL, NULL, vertex);
608 0 : qh_errexit(qh, qh_ERRqhull, NULL, NULL);
609 : }
610 12 : id= qh_pointid(qh, vertex->point);
611 12 : if ((qh->STOPpoint>0 && id == qh->STOPpoint-1) ||
612 12 : (qh->STOPpoint<0 && id == -qh->STOPpoint-1) ||
613 12 : (qh->STOPcone>0 && id == qh->STOPcone-1)) {
614 0 : trace1((qh, qh->ferr, 1038,"qh_buildhull: stop point or cone P%d in initial hull\n", id));
615 0 : return;
616 : }
617 : }
618 3 : qh->facet_next= qh->facet_list; /* advance facet when processed */
619 14663 : while ((furthest= qh_nextfurthest(qh, &facet))) {
620 14660 : qh->num_outside--; /* if ONLYmax, furthest may not be outside */
621 14660 : if (qh->STOPadd>0 && (qh->num_vertices - qh->hull_dim - 1 >= qh->STOPadd - 1)) {
622 0 : trace1((qh, qh->ferr, 1059, "qh_buildhull: stop after adding %d vertices\n", qh->STOPadd-1));
623 0 : return;
624 : }
625 14660 : if (!qh_addpoint(qh, furthest, facet, qh->ONLYmax))
626 0 : break;
627 : }
628 3 : if (qh->NARROWhull) /* move points from outsideset to coplanarset */
629 0 : qh_outcoplanar(qh /* facet_list */ );
630 3 : if (qh->num_outside && !furthest) {
631 0 : qh_fprintf(qh, qh->ferr, 6167, "qhull internal error (qh_buildhull): %d outside points were never processed.\n", qh->num_outside);
632 0 : qh_errexit(qh, qh_ERRqhull, NULL, NULL);
633 : }
634 3 : trace1((qh, qh->ferr, 1039, "qh_buildhull: completed the hull construction\n"));
635 : } /* buildhull */
636 :
637 :
638 : /*-<a href="qh-qhull_r.htm#TOC"
639 : >-------------------------------</a><a name="buildtracing">-</a>
640 :
641 : qh_buildtracing(qh, furthest, facet )
642 : trace an iteration of qh_buildhull() for furthest point and facet
643 : if !furthest, prints progress message
644 :
645 : returns:
646 : tracks progress with qh.lastreport, lastcpu, lastfacets, lastmerges, lastplanes, lastdist
647 : updates qh.furthest_id (-3 if furthest is NULL)
648 : also resets visit_id, vertext_visit on wrap around
649 :
650 : see:
651 : qh_tracemerging()
652 :
653 : design:
654 : if !furthest
655 : print progress message
656 : exit
657 : if 'TFn' iteration
658 : print progress message
659 : else if tracing
660 : trace furthest point and facet
661 : reset qh.visit_id and qh.vertex_visit if overflow may occur
662 : set qh.furthest_id for tracing
663 : */
664 14660 : void qh_buildtracing(qhT *qh, pointT *furthest, facetT *facet) {
665 14660 : realT dist= 0;
666 : double cpu;
667 : int total, furthestid;
668 : time_t timedata;
669 : struct tm *tp;
670 : vertexT *vertex;
671 :
672 14660 : qh->old_randomdist= qh->RANDOMdist;
673 14660 : qh->RANDOMdist= False;
674 14660 : if (!furthest) {
675 0 : time(&timedata);
676 0 : tp= localtime(&timedata);
677 0 : cpu= (double)qh_CPUclock - (double)qh->hulltime;
678 0 : cpu /= (double)qh_SECticks;
679 0 : total= zzval_(Ztotmerge) - zzval_(Zcyclehorizon) + zzval_(Zcyclefacettot);
680 0 : qh_fprintf(qh, qh->ferr, 8118, "\n\
681 : At %02d:%02d:%02d & %2.5g CPU secs, qhull has created %d facets and merged %d.\n\
682 : The current hull contains %d facets and %d vertices. Last point was p%d\n",
683 0 : tp->tm_hour, tp->tm_min, tp->tm_sec, cpu, qh->facet_id -1,
684 : total, qh->num_facets, qh->num_vertices, qh->furthest_id);
685 0 : return;
686 : }
687 14660 : furthestid= qh_pointid(qh, furthest);
688 : #ifndef qh_NOtrace
689 14660 : if (qh->TRACEpoint == furthestid) {
690 0 : trace1((qh, qh->ferr, 1053, "qh_buildtracing: start trace T%d for point TP%d above facet f%d\n", qh->TRACElevel, furthestid, facet->id));
691 0 : qh->IStracing= qh->TRACElevel;
692 0 : qh->qhmem.IStracing= qh->TRACElevel;
693 14660 : }else if (qh->TRACEpoint != qh_IDnone && qh->TRACEdist < REALmax/2) {
694 0 : qh->IStracing= 0;
695 0 : qh->qhmem.IStracing= 0;
696 : }
697 : #endif
698 14660 : if (qh->REPORTfreq && (qh->facet_id-1 > qh->lastreport + (unsigned int)qh->REPORTfreq)) {
699 0 : qh->lastreport= qh->facet_id-1;
700 0 : time(&timedata);
701 0 : tp= localtime(&timedata);
702 0 : cpu= (double)qh_CPUclock - (double)qh->hulltime;
703 0 : cpu /= (double)qh_SECticks;
704 0 : total= zzval_(Ztotmerge) - zzval_(Zcyclehorizon) + zzval_(Zcyclefacettot);
705 0 : zinc_(Zdistio);
706 0 : qh_distplane(qh, furthest, facet, &dist);
707 0 : qh_fprintf(qh, qh->ferr, 8119, "\n\
708 : At %02d:%02d:%02d & %2.5g CPU secs, qhull has created %d facets and merged %d.\n\
709 : The current hull contains %d facets and %d vertices. There are %d\n\
710 : outside points. Next is point p%d(v%d), %2.2g above f%d.\n",
711 0 : tp->tm_hour, tp->tm_min, tp->tm_sec, cpu, qh->facet_id -1,
712 0 : total, qh->num_facets, qh->num_vertices, qh->num_outside+1,
713 0 : furthestid, qh->vertex_id, dist, getid_(facet));
714 14660 : }else if (qh->IStracing >=1) {
715 0 : cpu= (double)qh_CPUclock - (double)qh->hulltime;
716 0 : cpu /= (double)qh_SECticks;
717 0 : qh_distplane(qh, furthest, facet, &dist);
718 0 : qh_fprintf(qh, qh->ferr, 1049, "qh_addpoint: add p%d(v%d) %2.2g above f%d to hull of %d facets, %d merges, %d outside at %4.4g CPU secs. Previous p%d(v%d) delta %4.4g CPU, %d facets, %d merges, %d hyperplanes, %d distplanes, %d retries\n",
719 0 : furthestid, qh->vertex_id, dist, getid_(facet), qh->num_facets, zzval_(Ztotmerge), qh->num_outside+1, cpu, qh->furthest_id, qh->vertex_id - 1,
720 0 : cpu - qh->lastcpu, qh->num_facets - qh->lastfacets, zzval_(Ztotmerge) - qh->lastmerges, zzval_(Zsetplane) - qh->lastplanes, zzval_(Zdistplane) - qh->lastdist, qh->retry_addpoint);
721 0 : qh->lastcpu= cpu;
722 0 : qh->lastfacets= qh->num_facets;
723 0 : qh->lastmerges= zzval_(Ztotmerge);
724 0 : qh->lastplanes= zzval_(Zsetplane);
725 0 : qh->lastdist= zzval_(Zdistplane);
726 : }
727 14660 : zmax_(Zvisit2max, (int)qh->visit_id/2);
728 14660 : if (qh->visit_id > (unsigned int) INT_MAX) { /* 31 bits */
729 0 : zinc_(Zvisit);
730 0 : if (!qh_checklists(qh, qh->facet_list)) {
731 0 : qh_fprintf(qh, qh->ferr, 6370, "qhull internal error: qh_checklists failed on reset of qh.visit_id %u\n", qh->visit_id);
732 0 : qh_errexit(qh, qh_ERRqhull, NULL, NULL);
733 : }
734 0 : qh->visit_id= 0;
735 0 : FORALLfacets
736 0 : facet->visitid= 0;
737 : }
738 14660 : zmax_(Zvvisit2max, (int)qh->vertex_visit/2);
739 14660 : if (qh->vertex_visit > (unsigned int) INT_MAX) { /* 31 bits */
740 0 : zinc_(Zvvisit);
741 0 : if (qh->visit_id && !qh_checklists(qh, qh->facet_list)) {
742 0 : qh_fprintf(qh, qh->ferr, 6371, "qhull internal error: qh_checklists failed on reset of qh.vertex_visit %u\n", qh->vertex_visit);
743 0 : qh_errexit(qh, qh_ERRqhull, NULL, NULL);
744 : }
745 0 : qh->vertex_visit= 0;
746 0 : FORALLvertices
747 0 : vertex->visitid= 0;
748 : }
749 14660 : qh->furthest_id= furthestid;
750 14660 : qh->RANDOMdist= qh->old_randomdist;
751 : } /* buildtracing */
752 :
753 : /*-<a href="qh-qhull_r.htm#TOC"
754 : >-------------------------------</a><a name="errexit2">-</a>
755 :
756 : qh_errexit2(qh, exitcode, facet, otherfacet )
757 : return exitcode to system after an error
758 : report two facets
759 :
760 : returns:
761 : assumes exitcode non-zero
762 :
763 : see:
764 : normally use qh_errexit() in user_r.c(reports a facet and a ridge)
765 : */
766 0 : void qh_errexit2(qhT *qh, int exitcode, facetT *facet, facetT *otherfacet) {
767 0 : qh->tracefacet= NULL; /* avoid infinite recursion through qh_fprintf */
768 0 : qh->traceridge= NULL;
769 0 : qh->tracevertex= NULL;
770 0 : qh_errprint(qh, "ERRONEOUS", facet, otherfacet, NULL, NULL);
771 0 : qh_errexit(qh, exitcode, NULL, NULL);
772 0 : } /* errexit2 */
773 :
774 :
775 : /*-<a href="qh-qhull_r.htm#TOC"
776 : >-------------------------------</a><a name="findhorizon">-</a>
777 :
778 : qh_findhorizon(qh, point, facet, goodvisible, goodhorizon )
779 : given a visible facet, find the point's horizon and visible facets
780 : for all facets, !facet-visible
781 :
782 : returns:
783 : returns qh.visible_list/num_visible with all visible facets
784 : marks visible facets with ->visible
785 : updates count of good visible and good horizon facets
786 : updates qh.max_outside, qh.max_vertex, facet->maxoutside
787 :
788 : see:
789 : similar to qh_delpoint()
790 :
791 : design:
792 : move facet to qh.visible_list at end of qh.facet_list
793 : for all visible facets
794 : for each unvisited neighbor of a visible facet
795 : compute distance of point to neighbor
796 : if point above neighbor
797 : move neighbor to end of qh.visible_list
798 : else if point is coplanar with neighbor
799 : update qh.max_outside, qh.max_vertex, neighbor->maxoutside
800 : mark neighbor coplanar (will create a samecycle later)
801 : update horizon statistics
802 : */
803 14660 : void qh_findhorizon(qhT *qh, pointT *point, facetT *facet, int *goodvisible, int *goodhorizon) {
804 : facetT *neighbor, **neighborp, *visible;
805 14660 : int numhorizon= 0, coplanar= 0;
806 : realT dist;
807 :
808 14660 : trace1((qh, qh->ferr, 1040, "qh_findhorizon: find horizon for point p%d facet f%d\n",qh_pointid(qh, point),facet->id));
809 14660 : *goodvisible= *goodhorizon= 0;
810 14660 : zinc_(Ztotvisible);
811 14660 : qh_removefacet(qh, facet); /* visible_list at end of qh->facet_list */
812 14660 : qh_appendfacet(qh, facet);
813 14660 : qh->num_visible= 1;
814 14660 : if (facet->good)
815 14660 : (*goodvisible)++;
816 14660 : qh->visible_list= facet;
817 14660 : facet->visible= True;
818 14660 : facet->f.replace= NULL;
819 14660 : if (qh->IStracing >=4)
820 0 : qh_errprint(qh, "visible", facet, NULL, NULL, NULL);
821 14660 : qh->visit_id++;
822 53569 : FORALLvisible_facets {
823 38909 : if (visible->tricoplanar && !qh->TRInormals) {
824 0 : qh_fprintf(qh, qh->ferr, 6230, "qhull internal error (qh_findhorizon): does not work for tricoplanar facets. Use option 'Q11'\n");
825 0 : qh_errexit(qh, qh_ERRqhull, visible, NULL);
826 : }
827 38909 : if (qh_setsize(qh, visible->neighbors) == 0) {
828 0 : qh_fprintf(qh, qh->ferr, 6295, "qhull internal error (qh_findhorizon): visible facet f%d does not have neighbors\n", visible->id);
829 0 : qh_errexit(qh, qh_ERRqhull, visible, NULL);
830 : }
831 38909 : visible->visitid= qh->visit_id;
832 160811 : FOREACHneighbor_(visible) {
833 121902 : if (neighbor->visitid == qh->visit_id)
834 24313 : continue;
835 97589 : neighbor->visitid= qh->visit_id;
836 97589 : zzinc_(Znumvisibility);
837 97589 : qh_distplane(qh, point, neighbor, &dist);
838 97589 : if (dist > qh->MINvisible) {
839 24249 : zinc_(Ztotvisible);
840 24249 : qh_removefacet(qh, neighbor); /* append to end of qh->visible_list */
841 24249 : qh_appendfacet(qh, neighbor);
842 24249 : neighbor->visible= True;
843 24249 : neighbor->f.replace= NULL;
844 24249 : qh->num_visible++;
845 24249 : if (neighbor->good)
846 24249 : (*goodvisible)++;
847 24249 : if (qh->IStracing >=4)
848 0 : qh_errprint(qh, "visible", neighbor, NULL, NULL, NULL);
849 : }else {
850 73340 : if (dist >= -qh->MAXcoplanar) {
851 20047 : neighbor->coplanarhorizon= True;
852 20047 : zzinc_(Zcoplanarhorizon);
853 20047 : qh_joggle_restart(qh, "coplanar horizon");
854 20047 : coplanar++;
855 20047 : if (qh->MERGING) {
856 20047 : if (dist > 0) {
857 9224 : maximize_(qh->max_outside, dist);
858 9224 : maximize_(qh->max_vertex, dist);
859 : #if qh_MAXoutside
860 9224 : maximize_(neighbor->maxoutside, dist);
861 : #endif
862 : }else
863 10823 : minimize_(qh->min_vertex, dist); /* due to merge later */
864 : }
865 20047 : trace2((qh, qh->ferr, 2057, "qh_findhorizon: point p%d is coplanar to horizon f%d, dist=%2.7g < qh->MINvisible(%2.7g)\n",
866 : qh_pointid(qh, point), neighbor->id, dist, qh->MINvisible));
867 : }else
868 53293 : neighbor->coplanarhorizon= False;
869 73340 : zinc_(Ztothorizon);
870 73340 : numhorizon++;
871 73340 : if (neighbor->good)
872 73340 : (*goodhorizon)++;
873 73340 : if (qh->IStracing >=4)
874 0 : qh_errprint(qh, "horizon", neighbor, NULL, NULL, NULL);
875 : }
876 : }
877 : }
878 14660 : if (!numhorizon) {
879 0 : qh_joggle_restart(qh, "empty horizon");
880 0 : qh_fprintf(qh, qh->ferr, 6168, "qhull topology error (qh_findhorizon): empty horizon for p%d. It was above all facets.\n", qh_pointid(qh, point));
881 0 : if (qh->num_facets < 100) {
882 0 : qh_printfacetlist(qh, qh->facet_list, NULL, True);
883 : }
884 0 : qh_errexit(qh, qh_ERRtopology, NULL, NULL);
885 : }
886 14660 : trace1((qh, qh->ferr, 1041, "qh_findhorizon: %d horizon facets(good %d), %d visible(good %d), %d coplanar\n",
887 : numhorizon, *goodhorizon, qh->num_visible, *goodvisible, coplanar));
888 14660 : if (qh->IStracing >= 4 && qh->num_facets < 100)
889 0 : qh_printlists(qh);
890 14660 : } /* findhorizon */
891 :
892 : /*-<a href="qh-qhull_r.htm#TOC"
893 : >-------------------------------</a><a name="joggle_restart">-</a>
894 :
895 : qh_joggle_restart(qh, reason )
896 : if joggle ('QJn') and not merging, restart on precision and topology errors
897 : */
898 20049 : void qh_joggle_restart(qhT *qh, const char *reason) {
899 :
900 20049 : if (qh->JOGGLEmax < REALmax/2) {
901 0 : if (qh->ALLOWrestart && !qh->PREmerge && !qh->MERGEexact) {
902 0 : trace0((qh, qh->ferr, 26, "qh_joggle_restart: qhull restart because of %s\n", reason));
903 : /* May be called repeatedly if qh->ALLOWrestart */
904 0 : longjmp(qh->restartexit, qh_ERRprec);
905 : }
906 : }
907 20049 : } /* qh_joggle_restart */
908 :
909 : /*-<a href="qh-qhull_r.htm#TOC"
910 : >-------------------------------</a><a name="nextfurthest">-</a>
911 :
912 : qh_nextfurthest(qh, visible )
913 : returns next furthest point and visible facet for qh_addpoint()
914 : starts search at qh.facet_next
915 :
916 : returns:
917 : removes furthest point from outside set
918 : NULL if none available
919 : advances qh.facet_next over facets with empty outside sets
920 :
921 : design:
922 : for each facet from qh.facet_next
923 : if empty outside set
924 : advance qh.facet_next
925 : else if qh.NARROWhull
926 : determine furthest outside point
927 : if furthest point is not outside
928 : advance qh.facet_next(point will be coplanar)
929 : remove furthest point from outside set
930 : */
931 14663 : pointT *qh_nextfurthest(qhT *qh, facetT **visible) {
932 : facetT *facet;
933 14663 : int size, idx, loopcount= 0;
934 : realT randr, dist;
935 : pointT *furthest;
936 :
937 42296 : while ((facet= qh->facet_next) != qh->facet_tail) {
938 42293 : if (!facet || loopcount++ > qh->num_facets) {
939 0 : qh_fprintf(qh, qh->ferr, 6406, "qhull internal error (qh_nextfurthest): null facet or infinite loop detected for qh.facet_next f%d facet_tail f%d\n",
940 0 : getid_(facet), getid_(qh->facet_tail));
941 0 : qh_printlists(qh);
942 0 : qh_errexit2(qh, qh_ERRqhull, facet, qh->facet_tail);
943 : }
944 42293 : if (!facet->outsideset) {
945 27633 : qh->facet_next= facet->next;
946 27633 : continue;
947 : }
948 14660 : SETreturnsize_(facet->outsideset, size);
949 14660 : if (!size) {
950 0 : qh_setfree(qh, &facet->outsideset);
951 0 : qh->facet_next= facet->next;
952 0 : continue;
953 : }
954 14660 : if (qh->NARROWhull) {
955 0 : if (facet->notfurthest)
956 0 : qh_furthestout(qh, facet);
957 0 : furthest= (pointT *)qh_setlast(facet->outsideset);
958 : #if qh_COMPUTEfurthest
959 : qh_distplane(qh, furthest, facet, &dist);
960 : zinc_(Zcomputefurthest);
961 : #else
962 0 : dist= facet->furthestdist;
963 : #endif
964 0 : if (dist < qh->MINoutside) { /* remainder of outside set is coplanar for qh_outcoplanar */
965 0 : qh->facet_next= facet->next;
966 0 : continue;
967 : }
968 : }
969 14660 : if (!qh->RANDOMoutside && !qh->VIRTUALmemory) {
970 14660 : if (qh->PICKfurthest) {
971 0 : qh_furthestnext(qh /* qh.facet_list */);
972 0 : facet= qh->facet_next;
973 : }
974 14660 : *visible= facet;
975 14660 : return ((pointT *)qh_setdellast(facet->outsideset));
976 : }
977 0 : if (qh->RANDOMoutside) {
978 0 : int outcoplanar= 0;
979 0 : if (qh->NARROWhull) {
980 0 : FORALLfacets {
981 0 : if (facet == qh->facet_next)
982 0 : break;
983 0 : if (facet->outsideset)
984 0 : outcoplanar += qh_setsize(qh, facet->outsideset);
985 : }
986 : }
987 0 : randr= qh_RANDOMint;
988 0 : randr= randr/(qh_RANDOMmax+1);
989 0 : randr= floor((qh->num_outside - outcoplanar) * randr);
990 0 : idx= (int)randr;
991 0 : FORALLfacet_(qh->facet_next) {
992 0 : if (facet->outsideset) {
993 0 : SETreturnsize_(facet->outsideset, size);
994 0 : if (!size)
995 0 : qh_setfree(qh, &facet->outsideset);
996 0 : else if (size > idx) {
997 0 : *visible= facet;
998 0 : return ((pointT *)qh_setdelnth(qh, facet->outsideset, idx));
999 : }else
1000 0 : idx -= size;
1001 : }
1002 : }
1003 0 : 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",
1004 : qh->num_outside, idx+1, randr);
1005 0 : qh_errexit(qh, qh_ERRqhull, NULL, NULL);
1006 : }else { /* VIRTUALmemory */
1007 0 : facet= qh->facet_tail->previous;
1008 0 : if (!(furthest= (pointT *)qh_setdellast(facet->outsideset))) {
1009 0 : if (facet->outsideset)
1010 0 : qh_setfree(qh, &facet->outsideset);
1011 0 : qh_removefacet(qh, facet);
1012 0 : qh_prependfacet(qh, facet, &qh->facet_list);
1013 0 : continue;
1014 : }
1015 0 : *visible= facet;
1016 0 : return furthest;
1017 : }
1018 : }
1019 3 : return NULL;
1020 : } /* nextfurthest */
1021 :
1022 : /*-<a href="qh-qhull_r.htm#TOC"
1023 : >-------------------------------</a><a name="partitionall">-</a>
1024 :
1025 : qh_partitionall(qh, vertices, points, numpoints )
1026 : partitions all points in points/numpoints to the outsidesets of facets
1027 : vertices= vertices in qh.facet_list(!partitioned)
1028 :
1029 : returns:
1030 : builds facet->outsideset
1031 : does not partition qh.GOODpoint
1032 : if qh.ONLYgood && !qh.MERGING,
1033 : does not partition qh.GOODvertex
1034 :
1035 : notes:
1036 : faster if qh.facet_list sorted by anticipated size of outside set
1037 :
1038 : design:
1039 : initialize pointset with all points
1040 : remove vertices from pointset
1041 : remove qh.GOODpointp from pointset (unless it's qh.STOPcone or qh.STOPpoint)
1042 : for all facets
1043 : for all remaining points in pointset
1044 : compute distance from point to facet
1045 : if point is outside facet
1046 : remove point from pointset (by not reappending)
1047 : update bestpoint
1048 : append point or old bestpoint to facet's outside set
1049 : append bestpoint to facet's outside set (furthest)
1050 : for all points remaining in pointset
1051 : partition point into facets' outside sets and coplanar sets
1052 : */
1053 3 : void qh_partitionall(qhT *qh, setT *vertices, pointT *points, int numpoints){
1054 : setT *pointset;
1055 : vertexT *vertex, **vertexp;
1056 : pointT *point, **pointp, *bestpoint;
1057 : int size, point_i, point_n, point_end, remaining, i, id;
1058 : facetT *facet;
1059 3 : realT bestdist= -REALmax, dist, distoutside;
1060 :
1061 3 : trace1((qh, qh->ferr, 1042, "qh_partitionall: partition all points into outside sets\n"));
1062 3 : pointset= qh_settemp(qh, numpoints);
1063 3 : qh->num_outside= 0;
1064 3 : pointp= SETaddr_(pointset, pointT);
1065 14676 : for (i=numpoints, point= points; i--; point += qh->hull_dim)
1066 14673 : *(pointp++)= point;
1067 3 : qh_settruncate(qh, pointset, numpoints);
1068 15 : FOREACHvertex_(vertices) {
1069 12 : if ((id= qh_pointid(qh, vertex->point)) >= 0)
1070 12 : SETelem_(pointset, id)= NULL;
1071 : }
1072 3 : id= qh_pointid(qh, qh->GOODpointp);
1073 3 : if (id >=0 && qh->STOPcone-1 != id && -qh->STOPpoint-1 != id)
1074 0 : SETelem_(pointset, id)= NULL;
1075 3 : if (qh->GOODvertexp && qh->ONLYgood && !qh->MERGING) { /* matches qhull()*/
1076 0 : if ((id= qh_pointid(qh, qh->GOODvertexp)) >= 0)
1077 0 : SETelem_(pointset, id)= NULL;
1078 : }
1079 3 : if (!qh->BESToutside) { /* matches conditional for qh_partitionpoint below */
1080 3 : distoutside= qh_DISToutside; /* multiple of qh.MINoutside & qh.max_outside, see user_r.h */
1081 3 : zval_(Ztotpartition)= qh->num_points - qh->hull_dim - 1; /*misses GOOD... */
1082 3 : remaining= qh->num_facets;
1083 3 : point_end= numpoints;
1084 15 : FORALLfacets {
1085 12 : size= point_end/(remaining--) + 100;
1086 12 : facet->outsideset= qh_setnew(qh, size);
1087 12 : bestpoint= NULL;
1088 12 : point_end= 0;
1089 14699 : FOREACHpoint_i_(qh, pointset) {
1090 14687 : if (point) {
1091 14675 : zzinc_(Zpartitionall);
1092 14675 : qh_distplane(qh, point, facet, &dist);
1093 14675 : if (dist < distoutside)
1094 14 : SETelem_(pointset, point_end++)= point;
1095 : else {
1096 14661 : qh->num_outside++;
1097 14661 : if (!bestpoint) {
1098 6 : bestpoint= point;
1099 6 : bestdist= dist;
1100 14655 : }else if (dist > bestdist) {
1101 533 : qh_setappend(qh, &facet->outsideset, bestpoint);
1102 533 : bestpoint= point;
1103 533 : bestdist= dist;
1104 : }else
1105 14122 : qh_setappend(qh, &facet->outsideset, point);
1106 : }
1107 : }
1108 : }
1109 12 : if (bestpoint) {
1110 6 : qh_setappend(qh, &facet->outsideset, bestpoint);
1111 : #if !qh_COMPUTEfurthest
1112 6 : facet->furthestdist= bestdist;
1113 : #endif
1114 : }else
1115 6 : qh_setfree(qh, &facet->outsideset);
1116 12 : qh_settruncate(qh, pointset, point_end);
1117 : }
1118 : }
1119 : /* if !qh->BESToutside, pointset contains points not assigned to outsideset */
1120 3 : if (qh->BESToutside || qh->MERGING || qh->KEEPcoplanar || qh->KEEPinside || qh->KEEPnearinside) {
1121 3 : qh->findbestnew= True;
1122 3 : FOREACHpoint_i_(qh, pointset) {
1123 0 : if (point)
1124 0 : qh_partitionpoint(qh, point, qh->facet_list);
1125 : }
1126 3 : qh->findbestnew= False;
1127 : }
1128 3 : zzadd_(Zpartitionall, zzval_(Zpartition));
1129 3 : zzval_(Zpartition)= 0;
1130 3 : qh_settempfree(qh, &pointset);
1131 3 : if (qh->IStracing >= 4)
1132 0 : qh_printfacetlist(qh, qh->facet_list, NULL, True);
1133 3 : } /* partitionall */
1134 :
1135 :
1136 : /*-<a href="qh-qhull_r.htm#TOC"
1137 : >-------------------------------</a><a name="partitioncoplanar">-</a>
1138 :
1139 : qh_partitioncoplanar(qh, point, facet, dist, allnew )
1140 : partition coplanar point to a facet
1141 : dist is distance from point to facet
1142 : if dist NULL,
1143 : searches for bestfacet and does nothing if inside
1144 : if allnew (qh.findbestnew)
1145 : searches new facets instead of using qh_findbest()
1146 :
1147 : returns:
1148 : qh.max_ouside updated
1149 : if qh.KEEPcoplanar or qh.KEEPinside
1150 : point assigned to best coplanarset
1151 : qh.repart_facetid==0 (for detecting infinite recursion via qh_partitionpoint)
1152 :
1153 : notes:
1154 : facet->maxoutside is updated at end by qh_check_maxout
1155 :
1156 : design:
1157 : if dist undefined
1158 : find best facet for point
1159 : if point sufficiently below facet (depends on qh.NEARinside and qh.KEEPinside)
1160 : exit
1161 : if keeping coplanar/nearinside/inside points
1162 : if point is above furthest coplanar point
1163 : append point to coplanar set (it is the new furthest)
1164 : update qh.max_outside
1165 : else
1166 : append point one before end of coplanar set
1167 : else if point is clearly outside of qh.max_outside and bestfacet->coplanarset
1168 : and bestfacet is more than perpendicular to facet
1169 : repartition the point using qh_findbest() -- it may be put on an outsideset
1170 : else
1171 : update qh.max_outside
1172 : */
1173 4 : void qh_partitioncoplanar(qhT *qh, pointT *point, facetT *facet, realT *dist, boolT allnew) {
1174 : facetT *bestfacet;
1175 : pointT *oldfurthest;
1176 4 : realT bestdist, angle, nearest, dist2= 0.0;
1177 4 : int numpart= 0;
1178 4 : boolT isoutside, oldfindbest, repartition= False;
1179 :
1180 4 : trace4((qh, qh->ferr, 4090, "qh_partitioncoplanar: partition coplanar point p%d starting with f%d dist? %2.2g, allnew? %d, gh.repart_facetid f%d\n",
1181 : qh_pointid(qh, point), facet->id, (dist ? *dist : 0.0), allnew, qh->repart_facetid));
1182 4 : qh->WAScoplanar= True;
1183 4 : if (!dist) {
1184 3 : if (allnew)
1185 1 : bestfacet= qh_findbestnew(qh, point, facet, &bestdist, qh_ALL, &isoutside, &numpart);
1186 : else
1187 2 : bestfacet= qh_findbest(qh, point, facet, qh_ALL, !qh_ISnewfacets, qh->DELAUNAY,
1188 : &bestdist, &isoutside, &numpart);
1189 3 : zinc_(Ztotpartcoplanar);
1190 3 : zzadd_(Zpartcoplanar, numpart);
1191 3 : if (!qh->DELAUNAY && !qh->KEEPinside) { /* for 'd', bestdist skips upperDelaunay facets */
1192 0 : if (qh->KEEPnearinside) {
1193 0 : if (bestdist < -qh->NEARinside) {
1194 0 : zinc_(Zcoplanarinside);
1195 0 : trace4((qh, qh->ferr, 4062, "qh_partitioncoplanar: point p%d is more than near-inside facet f%d dist %2.2g allnew? %d\n",
1196 : qh_pointid(qh, point), bestfacet->id, bestdist, allnew));
1197 0 : qh->repart_facetid= 0;
1198 0 : return;
1199 : }
1200 0 : }else if (bestdist < -qh->MAXcoplanar) {
1201 0 : trace4((qh, qh->ferr, 4063, "qh_partitioncoplanar: point p%d is inside facet f%d dist %2.2g allnew? %d\n",
1202 : qh_pointid(qh, point), bestfacet->id, bestdist, allnew));
1203 0 : zinc_(Zcoplanarinside);
1204 0 : qh->repart_facetid= 0;
1205 0 : return;
1206 : }
1207 : }
1208 : }else {
1209 1 : bestfacet= facet;
1210 1 : bestdist= *dist;
1211 : }
1212 4 : if(bestfacet->visible){
1213 0 : qh_fprintf(qh, qh->ferr, 6405, "qhull internal error (qh_partitioncoplanar): cannot partition coplanar p%d of f%d into visible facet f%d\n",
1214 : qh_pointid(qh, point), facet->id, bestfacet->id);
1215 0 : qh_errexit2(qh, qh_ERRqhull, facet, bestfacet);
1216 : }
1217 4 : if (bestdist > qh->max_outside) {
1218 1 : if (!dist && facet != bestfacet) { /* can't be recursive from qh_partitionpoint since facet != bestfacet */
1219 0 : zinc_(Zpartangle);
1220 0 : angle= qh_getangle(qh, facet->normal, bestfacet->normal);
1221 0 : if (angle < 0) {
1222 : /* nearest= qh_vertex_bestdist(qh, bestfacet->vertices) */;
1223 : /* typically due to deleted vertex and coplanar facets, e.g.,
1224 : RBOX 1000 s Z1 G1e-13 t1001185205 | QHULL Tv */
1225 0 : zinc_(Zpartcorner);
1226 0 : trace2((qh, qh->ferr, 2058, "qh_partitioncoplanar: repartition coplanar point p%d from f%d as an outside point above corner facet f%d dist %2.2g with angle %2.2g\n",
1227 : qh_pointid(qh, point), facet->id, bestfacet->id, bestdist, angle));
1228 0 : repartition= True;
1229 : }
1230 : }
1231 1 : if (!repartition) {
1232 1 : if (bestdist > qh->MAXoutside * qh_RATIOcoplanaroutside) {
1233 0 : nearest= qh_vertex_bestdist(qh, bestfacet->vertices);
1234 0 : if (facet->id == bestfacet->id) {
1235 0 : if (facet->id == qh->repart_facetid) {
1236 0 : qh_fprintf(qh, qh->ferr, 6404, "Qhull internal error (qh_partitioncoplanar): infinite loop due to recursive call to qh_partitionpoint. Repartition point p%d from f%d as a outside point dist %2.2g nearest vertices %2.2g\n",
1237 : qh_pointid(qh, point), facet->id, bestdist, nearest);
1238 0 : qh_errexit(qh, qh_ERRqhull, facet, NULL);
1239 : }
1240 0 : qh->repart_facetid= facet->id; /* reset after call to qh_partitionpoint */
1241 : }
1242 0 : if (point == qh->coplanar_apex) {
1243 : /* otherwise may loop indefinitely, the point is well above a facet, yet near a vertex */
1244 0 : qh_fprintf(qh, qh->ferr, 6425, "Qhull topology error (qh_partitioncoplanar): can not repartition coplanar point p%d from f%d as outside point above f%d. It previously failed to form a cone of facets, dist %2.2g, nearest vertices %2.2g\n",
1245 : qh_pointid(qh, point), facet->id, bestfacet->id, bestdist, nearest);
1246 0 : qh_errexit(qh, qh_ERRtopology, facet, NULL);
1247 : }
1248 0 : if (nearest < 2 * qh->MAXoutside * qh_RATIOcoplanaroutside) {
1249 0 : zinc_(Zparttwisted);
1250 0 : qh_fprintf(qh, qh->ferr, 7085, "Qhull precision warning: repartition coplanar point p%d from f%d as an outside point above twisted facet f%d dist %2.2g nearest vertices %2.2g\n",
1251 : qh_pointid(qh, point), facet->id, bestfacet->id, bestdist, nearest);
1252 : }else {
1253 0 : zinc_(Zparthidden);
1254 0 : qh_fprintf(qh, qh->ferr, 7086, "Qhull precision warning: repartition coplanar point p%d from f%d as an outside point above hidden facet f%d dist %2.2g nearest vertices %2.2g\n",
1255 : qh_pointid(qh, point), facet->id, bestfacet->id, bestdist, nearest);
1256 : }
1257 0 : repartition= True;
1258 : }
1259 : }
1260 1 : if (repartition) {
1261 0 : oldfindbest= qh->findbestnew;
1262 0 : qh->findbestnew= False;
1263 0 : qh_partitionpoint(qh, point, bestfacet);
1264 0 : qh->findbestnew= oldfindbest;
1265 0 : qh->repart_facetid= 0;
1266 0 : return;
1267 : }
1268 1 : qh->repart_facetid= 0;
1269 1 : qh->max_outside= bestdist;
1270 1 : if (bestdist > qh->TRACEdist || qh->IStracing >= 3) {
1271 0 : qh_fprintf(qh, qh->ferr, 3041, "qh_partitioncoplanar: == p%d from f%d increases qh.max_outside to %2.2g of f%d last p%d\n",
1272 : qh_pointid(qh, point), facet->id, bestdist, bestfacet->id, qh->furthest_id);
1273 0 : qh_errprint(qh, "DISTANT", facet, bestfacet, NULL, NULL);
1274 : }
1275 : }
1276 4 : if (qh->KEEPcoplanar + qh->KEEPinside + qh->KEEPnearinside) {
1277 4 : oldfurthest= (pointT *)qh_setlast(bestfacet->coplanarset);
1278 4 : if (oldfurthest) {
1279 0 : zinc_(Zcomputefurthest);
1280 0 : qh_distplane(qh, oldfurthest, bestfacet, &dist2);
1281 : }
1282 4 : if (!oldfurthest || dist2 < bestdist)
1283 4 : qh_setappend(qh, &bestfacet->coplanarset, point);
1284 : else
1285 0 : qh_setappend2ndlast(qh, &bestfacet->coplanarset, point);
1286 : }
1287 4 : trace4((qh, qh->ferr, 4064, "qh_partitioncoplanar: point p%d is coplanar with facet f%d (or inside) dist %2.2g\n",
1288 : qh_pointid(qh, point), bestfacet->id, bestdist));
1289 : } /* partitioncoplanar */
1290 :
1291 : /*-<a href="qh-qhull_r.htm#TOC"
1292 : >-------------------------------</a><a name="partitionpoint">-</a>
1293 :
1294 : qh_partitionpoint(qh, point, facet )
1295 : assigns point to an outside set, coplanar set, or inside set (i.e., dropt)
1296 : if qh.findbestnew
1297 : uses qh_findbestnew() to search all new facets
1298 : else
1299 : uses qh_findbest()
1300 :
1301 : notes:
1302 : after qh_distplane(), this and qh_findbest() are most expensive in 3-d
1303 :
1304 : design:
1305 : find best facet for point
1306 : (either exhaustive search of new facets or directed search from facet)
1307 : if qh.NARROWhull
1308 : retain coplanar and nearinside points as outside points
1309 : if point is outside bestfacet
1310 : if point above furthest point for bestfacet
1311 : append point to outside set (it becomes the new furthest)
1312 : if outside set was empty
1313 : move bestfacet to end of qh.facet_list (i.e., after qh.facet_next)
1314 : update bestfacet->furthestdist
1315 : else
1316 : append point one before end of outside set
1317 : else if point is coplanar to bestfacet
1318 : if keeping coplanar points or need to update qh.max_outside
1319 : partition coplanar point into bestfacet
1320 : else if near-inside point
1321 : partition as coplanar point into bestfacet
1322 : else is an inside point
1323 : if keeping inside points
1324 : partition as coplanar point into bestfacet
1325 : */
1326 208091 : void qh_partitionpoint(qhT *qh, pointT *point, facetT *facet) {
1327 : realT bestdist, previousdist;
1328 208091 : boolT isoutside, isnewoutside= False;
1329 : facetT *bestfacet;
1330 : int numpart;
1331 :
1332 208091 : if (qh->findbestnew)
1333 184782 : bestfacet= qh_findbestnew(qh, point, facet, &bestdist, qh->BESToutside, &isoutside, &numpart);
1334 : else
1335 23309 : bestfacet= qh_findbest(qh, point, facet, qh->BESToutside, qh_ISnewfacets, !qh_NOupper,
1336 : &bestdist, &isoutside, &numpart);
1337 208091 : zinc_(Ztotpartition);
1338 208091 : zzadd_(Zpartition, numpart);
1339 208091 : if(bestfacet->visible){
1340 0 : qh_fprintf(qh, qh->ferr, 6293, "qhull internal error (qh_partitionpoint): cannot partition p%d of f%d into visible facet f%d\n",
1341 : qh_pointid(qh, point), facet->id, bestfacet->id);
1342 0 : qh_errexit2(qh, qh_ERRqhull, facet, bestfacet);
1343 : }
1344 208091 : if (qh->NARROWhull) {
1345 0 : if (qh->DELAUNAY && !isoutside && bestdist >= -qh->MAXcoplanar)
1346 0 : qh_joggle_restart(qh, "nearly incident point (narrow hull)");
1347 0 : if (qh->KEEPnearinside) {
1348 0 : if (bestdist >= -qh->NEARinside)
1349 0 : isoutside= True;
1350 0 : }else if (bestdist >= -qh->MAXcoplanar)
1351 0 : isoutside= True;
1352 : }
1353 :
1354 208091 : if (isoutside) {
1355 208090 : if (!bestfacet->outsideset
1356 181577 : || !qh_setlast(bestfacet->outsideset)) { /* empty outside set */
1357 26513 : qh_setappend(qh, &(bestfacet->outsideset), point);
1358 26513 : if (!qh->NARROWhull || bestdist > qh->MINoutside)
1359 26513 : isnewoutside= True;
1360 : #if !qh_COMPUTEfurthest
1361 26513 : bestfacet->furthestdist= bestdist;
1362 : #endif
1363 : }else {
1364 : #if qh_COMPUTEfurthest
1365 : zinc_(Zcomputefurthest);
1366 : qh_distplane(qh, oldfurthest, bestfacet, &previousdist);
1367 : if (previousdist < bestdist)
1368 : qh_setappend(qh, &(bestfacet->outsideset), point);
1369 : else
1370 : qh_setappend2ndlast(qh, &(bestfacet->outsideset), point);
1371 : #else
1372 181577 : previousdist= bestfacet->furthestdist;
1373 181577 : if (previousdist < bestdist) {
1374 31526 : qh_setappend(qh, &(bestfacet->outsideset), point);
1375 31526 : bestfacet->furthestdist= bestdist;
1376 31526 : if (qh->NARROWhull && previousdist < qh->MINoutside && bestdist >= qh->MINoutside)
1377 0 : isnewoutside= True;
1378 : }else
1379 150051 : qh_setappend2ndlast(qh, &(bestfacet->outsideset), point);
1380 : #endif
1381 : }
1382 208090 : if (isnewoutside && qh->facet_next != bestfacet) {
1383 26513 : if (bestfacet->newfacet) {
1384 26513 : if (qh->facet_next->newfacet)
1385 0 : qh->facet_next= qh->newfacet_list; /* make sure it's after qh.facet_next */
1386 : }else {
1387 0 : qh_removefacet(qh, bestfacet); /* make sure it's after qh.facet_next */
1388 0 : qh_appendfacet(qh, bestfacet);
1389 0 : if(qh->newfacet_list){
1390 0 : bestfacet->newfacet= True;
1391 : }
1392 : }
1393 : }
1394 208090 : qh->num_outside++;
1395 208090 : trace4((qh, qh->ferr, 4065, "qh_partitionpoint: point p%d is outside facet f%d newfacet? %d, newoutside? %d (or narrowhull)\n",
1396 : qh_pointid(qh, point), bestfacet->id, bestfacet->newfacet, isnewoutside));
1397 1 : }else if (qh->DELAUNAY || bestdist >= -qh->MAXcoplanar) { /* for 'd', bestdist skips upperDelaunay facets */
1398 1 : if (qh->DELAUNAY)
1399 1 : qh_joggle_restart(qh, "nearly incident point");
1400 : /* allow coplanar points with joggle, may be interior */
1401 1 : zzinc_(Zcoplanarpart);
1402 1 : if ((qh->KEEPcoplanar + qh->KEEPnearinside) || bestdist > qh->max_outside)
1403 1 : qh_partitioncoplanar(qh, point, bestfacet, &bestdist, qh->findbestnew);
1404 : else {
1405 0 : trace4((qh, qh->ferr, 4066, "qh_partitionpoint: point p%d is coplanar to facet f%d (dropped)\n",
1406 : qh_pointid(qh, point), bestfacet->id));
1407 : }
1408 0 : }else if (qh->KEEPnearinside && bestdist >= -qh->NEARinside) {
1409 0 : zinc_(Zpartnear);
1410 0 : qh_partitioncoplanar(qh, point, bestfacet, &bestdist, qh->findbestnew);
1411 : }else {
1412 0 : zinc_(Zpartinside);
1413 0 : trace4((qh, qh->ferr, 4067, "qh_partitionpoint: point p%d is inside all facets, closest to f%d dist %2.2g\n",
1414 : qh_pointid(qh, point), bestfacet->id, bestdist));
1415 0 : if (qh->KEEPinside)
1416 0 : qh_partitioncoplanar(qh, point, bestfacet, &bestdist, qh->findbestnew);
1417 : }
1418 208091 : } /* partitionpoint */
1419 :
1420 : /*-<a href="qh-qhull_r.htm#TOC"
1421 : >-------------------------------</a><a name="partitionvisible">-</a>
1422 :
1423 : qh_partitionvisible(qh, allpoints, numoutside )
1424 : partitions outside points in visible facets (qh.visible_list) to qh.newfacet_list
1425 : if keeping coplanar/near-inside/inside points
1426 : partitions coplanar points; repartitions if 'allpoints' (not used)
1427 : 1st neighbor (if any) of visible facets points to a horizon facet or a new facet
1428 :
1429 : returns:
1430 : updates outside sets and coplanar sets of qh.newfacet_list
1431 : updates qh.num_outside (count of outside points)
1432 : does not truncate f.outsideset, f.coplanarset, or qh.del_vertices (see qh_deletevisible)
1433 :
1434 : notes:
1435 : called by qh_qhull, qh_addpoint, and qh_all_vertexmerges
1436 : qh.findbest_notsharp should be clear (extra work if set)
1437 :
1438 : design:
1439 : for all visible facets with outside set or coplanar set
1440 : select a newfacet for visible facet
1441 : if outside set
1442 : partition outside set into new facets
1443 : if coplanar set and keeping coplanar/near-inside/inside points
1444 : if allpoints
1445 : partition coplanar set into new facets, may be assigned outside
1446 : else
1447 : partition coplanar set into coplanar sets of new facets
1448 : for each deleted vertex
1449 : if allpoints
1450 : partition vertex into new facets, may be assigned outside
1451 : else
1452 : partition vertex into coplanar sets of new facets
1453 : */
1454 14660 : void qh_partitionvisible(qhT *qh, boolT allpoints, int *numoutside /* qh.visible_list */) {
1455 : facetT *visible, *newfacet;
1456 : pointT *point, **pointp;
1457 14660 : int delsize, coplanar=0, size;
1458 : vertexT *vertex, **vertexp;
1459 :
1460 14660 : trace3((qh, qh->ferr, 3042, "qh_partitionvisible: partition outside and coplanar points of visible and merged facets f%d into new facets f%d\n",
1461 : qh->visible_list->id, qh->newfacet_list->id));
1462 14660 : if (qh->ONLYmax)
1463 0 : maximize_(qh->MINoutside, qh->max_vertex);
1464 14660 : *numoutside= 0;
1465 73617 : FORALLvisible_facets {
1466 58957 : if (!visible->outsideset && !visible->coplanarset)
1467 32435 : continue;
1468 26522 : newfacet= qh_getreplacement(qh, visible);
1469 26522 : if (!newfacet)
1470 637 : newfacet= qh->newfacet_list;
1471 26522 : if (!newfacet->next) {
1472 0 : qh_fprintf(qh, qh->ferr, 6170, "qhull topology error (qh_partitionvisible): all new facets deleted as\n degenerate facets. Can not continue.\n");
1473 0 : qh_errexit(qh, qh_ERRtopology, NULL, NULL);
1474 : }
1475 26522 : if (visible->outsideset) {
1476 26519 : size= qh_setsize(qh, visible->outsideset);
1477 26519 : *numoutside += size;
1478 26519 : qh->num_outside -= size;
1479 234610 : FOREACHpoint_(visible->outsideset)
1480 208091 : qh_partitionpoint(qh, point, newfacet);
1481 : }
1482 26522 : if (visible->coplanarset && (qh->KEEPcoplanar + qh->KEEPinside + qh->KEEPnearinside)) {
1483 3 : size= qh_setsize(qh, visible->coplanarset);
1484 3 : coplanar += size;
1485 6 : FOREACHpoint_(visible->coplanarset) {
1486 3 : if (allpoints) /* not used */
1487 0 : qh_partitionpoint(qh, point, newfacet);
1488 : else
1489 3 : qh_partitioncoplanar(qh, point, newfacet, NULL, qh->findbestnew);
1490 : }
1491 : }
1492 : }
1493 14660 : delsize= qh_setsize(qh, qh->del_vertices);
1494 14660 : if (delsize > 0) {
1495 0 : trace3((qh, qh->ferr, 3049, "qh_partitionvisible: partition %d deleted vertices as coplanar? %d points into new facets f%d\n",
1496 : delsize, !allpoints, qh->newfacet_list->id));
1497 0 : FOREACHvertex_(qh->del_vertices) {
1498 0 : if (vertex->point && !vertex->partitioned) {
1499 0 : if (!qh->newfacet_list || qh->newfacet_list == qh->facet_tail) {
1500 0 : qh_fprintf(qh, qh->ferr, 6284, "qhull internal error (qh_partitionvisible): all new facets deleted or none defined. Can not partition deleted v%d.\n", vertex->id);
1501 0 : qh_errexit(qh, qh_ERRqhull, NULL, NULL);
1502 : }
1503 0 : if (allpoints) /* not used */
1504 : /* [apr'2019] infinite loop if vertex recreates the same facets from the same horizon
1505 : e.g., qh_partitionpoint if qh.DELAUNAY with qh.MERGEindependent for all mergetype, ../eg/qtest.sh t427764 '1000 s W1e-13 D3' 'd' */
1506 0 : qh_partitionpoint(qh, vertex->point, qh->newfacet_list);
1507 : else
1508 0 : qh_partitioncoplanar(qh, vertex->point, qh->newfacet_list, NULL, qh_ALL); /* search all new facets */
1509 0 : vertex->partitioned= True;
1510 : }
1511 : }
1512 : }
1513 14660 : trace1((qh, qh->ferr, 1043,"qh_partitionvisible: partitioned %d points from outsidesets, %d points from coplanarsets, and %d deleted vertices\n", *numoutside, coplanar, delsize));
1514 14660 : } /* partitionvisible */
1515 :
1516 : /*-<a href="qh-qhull_r.htm#TOC"
1517 : >-------------------------------</a><a name="printsummary">-</a>
1518 :
1519 : qh_printsummary(qh, fp )
1520 : prints summary to fp
1521 :
1522 : notes:
1523 : not in io_r.c so that user_eg.c can prevent io_r.c from loading
1524 : qh_printsummary and qh_countfacets must match counts
1525 : updates qh.facet_visit to detect infinite loop
1526 :
1527 : design:
1528 : determine number of points, vertices, and coplanar points
1529 : print summary
1530 : */
1531 0 : void qh_printsummary(qhT *qh, FILE *fp) {
1532 : realT ratio, outerplane, innerplane;
1533 : double cpu;
1534 0 : int size, id, nummerged, numpinched, numvertices, numcoplanars= 0, nonsimplicial=0, numdelaunay= 0;
1535 : facetT *facet;
1536 : const char *s;
1537 0 : int numdel= zzval_(Zdelvertextot);
1538 0 : int numtricoplanars= 0;
1539 : boolT goodused;
1540 :
1541 0 : size= qh->num_points + qh_setsize(qh, qh->other_points);
1542 0 : numvertices= qh->num_vertices - qh_setsize(qh, qh->del_vertices);
1543 0 : id= qh_pointid(qh, qh->GOODpointp);
1544 0 : if (!qh_checklists(qh, qh->facet_list) && !qh->ERREXITcalled) {
1545 0 : qh_fprintf(qh, fp, 6372, "qhull internal error: qh_checklists failed at qh_printsummary\n");
1546 0 : if (qh->num_facets < 4000)
1547 0 : qh_printlists(qh);
1548 0 : qh_errexit(qh, qh_ERRqhull, NULL, NULL);
1549 : }
1550 0 : if (qh->DELAUNAY && qh->ERREXITcalled) {
1551 : /* update f.good and determine qh.num_good as in qh_findgood_all */
1552 0 : FORALLfacets {
1553 0 : if (facet->visible)
1554 0 : facet->good= False; /* will be deleted */
1555 0 : else if (facet->good) {
1556 0 : if (facet->normal && !qh_inthresholds(qh, facet->normal, NULL))
1557 0 : facet->good= False;
1558 : else
1559 0 : numdelaunay++;
1560 : }
1561 : }
1562 0 : qh->num_good= numdelaunay;
1563 : }
1564 0 : FORALLfacets {
1565 0 : if (facet->coplanarset)
1566 0 : numcoplanars += qh_setsize(qh, facet->coplanarset);
1567 0 : if (facet->good) {
1568 0 : if (facet->simplicial) {
1569 0 : if (facet->keepcentrum && facet->tricoplanar)
1570 0 : numtricoplanars++;
1571 0 : }else if (qh_setsize(qh, facet->vertices) != qh->hull_dim)
1572 0 : nonsimplicial++;
1573 : }
1574 : }
1575 0 : if (id >=0 && qh->STOPcone-1 != id && -qh->STOPpoint-1 != id)
1576 0 : size--;
1577 0 : if (qh->STOPadd || qh->STOPcone || qh->STOPpoint)
1578 0 : qh_fprintf(qh, fp, 9288, "\nEarly exit due to 'TAn', 'TVn', 'TCn', 'TRn', or precision error with 'QJn'.");
1579 0 : goodused= False;
1580 0 : if (qh->ERREXITcalled)
1581 : ; /* qh_findgood_all not called */
1582 0 : else if (qh->UPPERdelaunay) {
1583 0 : if (qh->GOODvertex || qh->GOODpoint || qh->SPLITthresholds)
1584 0 : goodused= True;
1585 0 : }else if (qh->DELAUNAY) {
1586 0 : if (qh->GOODvertex || qh->GOODpoint || qh->GOODthreshold)
1587 0 : goodused= True;
1588 0 : }else if (qh->num_good > 0 || qh->GOODthreshold)
1589 0 : goodused= True;
1590 0 : nummerged= zzval_(Ztotmerge) - zzval_(Zcyclehorizon) + zzval_(Zcyclefacettot);
1591 0 : if (qh->VORONOI) {
1592 0 : if (qh->UPPERdelaunay)
1593 0 : qh_fprintf(qh, fp, 9289, "\n\
1594 : Furthest-site Voronoi vertices by the convex hull of %d points in %d-d:\n\n", size, qh->hull_dim);
1595 : else
1596 0 : qh_fprintf(qh, fp, 9290, "\n\
1597 : Voronoi diagram by the convex hull of %d points in %d-d:\n\n", size, qh->hull_dim);
1598 0 : qh_fprintf(qh, fp, 9291, " Number of Voronoi regions%s: %d\n",
1599 0 : qh->ATinfinity ? " and at-infinity" : "", numvertices);
1600 0 : if (numdel)
1601 0 : qh_fprintf(qh, fp, 9292, " Total number of deleted points due to merging: %d\n", numdel);
1602 0 : if (numcoplanars - numdel > 0)
1603 0 : qh_fprintf(qh, fp, 9293, " Number of nearly incident points: %d\n", numcoplanars - numdel);
1604 0 : else if (size - numvertices - numdel > 0)
1605 0 : qh_fprintf(qh, fp, 9294, " Total number of nearly incident points: %d\n", size - numvertices - numdel);
1606 0 : qh_fprintf(qh, fp, 9295, " Number of%s Voronoi vertices: %d\n",
1607 : goodused ? " 'good'" : "", qh->num_good);
1608 0 : if (nonsimplicial)
1609 0 : qh_fprintf(qh, fp, 9296, " Number of%s non-simplicial Voronoi vertices: %d\n",
1610 : goodused ? " 'good'" : "", nonsimplicial);
1611 0 : }else if (qh->DELAUNAY) {
1612 0 : if (qh->UPPERdelaunay)
1613 0 : qh_fprintf(qh, fp, 9297, "\n\
1614 : Furthest-site Delaunay triangulation by the convex hull of %d points in %d-d:\n\n", size, qh->hull_dim);
1615 : else
1616 0 : qh_fprintf(qh, fp, 9298, "\n\
1617 : Delaunay triangulation by the convex hull of %d points in %d-d:\n\n", size, qh->hull_dim);
1618 0 : qh_fprintf(qh, fp, 9299, " Number of input sites%s: %d\n",
1619 0 : qh->ATinfinity ? " and at-infinity" : "", numvertices);
1620 0 : if (numdel)
1621 0 : qh_fprintf(qh, fp, 9300, " Total number of deleted points due to merging: %d\n", numdel);
1622 0 : if (numcoplanars - numdel > 0)
1623 0 : qh_fprintf(qh, fp, 9301, " Number of nearly incident points: %d\n", numcoplanars - numdel);
1624 0 : else if (size - numvertices - numdel > 0)
1625 0 : qh_fprintf(qh, fp, 9302, " Total number of nearly incident points: %d\n", size - numvertices - numdel);
1626 0 : qh_fprintf(qh, fp, 9303, " Number of%s Delaunay regions: %d\n",
1627 : goodused ? " 'good'" : "", qh->num_good);
1628 0 : if (nonsimplicial)
1629 0 : qh_fprintf(qh, fp, 9304, " Number of%s non-simplicial Delaunay regions: %d\n",
1630 : goodused ? " 'good'" : "", nonsimplicial);
1631 0 : }else if (qh->HALFspace) {
1632 0 : qh_fprintf(qh, fp, 9305, "\n\
1633 : Halfspace intersection by the convex hull of %d points in %d-d:\n\n", size, qh->hull_dim);
1634 0 : qh_fprintf(qh, fp, 9306, " Number of halfspaces: %d\n", size);
1635 0 : qh_fprintf(qh, fp, 9307, " Number of non-redundant halfspaces: %d\n", numvertices);
1636 0 : if (numcoplanars) {
1637 0 : if (qh->KEEPinside && qh->KEEPcoplanar)
1638 0 : s= "similar and redundant";
1639 0 : else if (qh->KEEPinside)
1640 0 : s= "redundant";
1641 : else
1642 0 : s= "similar";
1643 0 : qh_fprintf(qh, fp, 9308, " Number of %s halfspaces: %d\n", s, numcoplanars);
1644 : }
1645 0 : qh_fprintf(qh, fp, 9309, " Number of intersection points: %d\n", qh->num_facets - qh->num_visible);
1646 0 : if (goodused)
1647 0 : qh_fprintf(qh, fp, 9310, " Number of 'good' intersection points: %d\n", qh->num_good);
1648 0 : if (nonsimplicial)
1649 0 : qh_fprintf(qh, fp, 9311, " Number of%s non-simplicial intersection points: %d\n",
1650 : goodused ? " 'good'" : "", nonsimplicial);
1651 : }else {
1652 0 : qh_fprintf(qh, fp, 9312, "\n\
1653 : Convex hull of %d points in %d-d:\n\n", size, qh->hull_dim);
1654 0 : qh_fprintf(qh, fp, 9313, " Number of vertices: %d\n", numvertices);
1655 0 : if (numcoplanars) {
1656 0 : if (qh->KEEPinside && qh->KEEPcoplanar)
1657 0 : s= "coplanar and interior";
1658 0 : else if (qh->KEEPinside)
1659 0 : s= "interior";
1660 : else
1661 0 : s= "coplanar";
1662 0 : qh_fprintf(qh, fp, 9314, " Number of %s points: %d\n", s, numcoplanars);
1663 : }
1664 0 : qh_fprintf(qh, fp, 9315, " Number of facets: %d\n", qh->num_facets - qh->num_visible);
1665 0 : if (goodused)
1666 0 : qh_fprintf(qh, fp, 9316, " Number of 'good' facets: %d\n", qh->num_good);
1667 0 : if (nonsimplicial)
1668 0 : qh_fprintf(qh, fp, 9317, " Number of%s non-simplicial facets: %d\n",
1669 : goodused ? " 'good'" : "", nonsimplicial);
1670 : }
1671 0 : if (numtricoplanars)
1672 0 : qh_fprintf(qh, fp, 9318, " Number of triangulated facets: %d\n", numtricoplanars);
1673 0 : qh_fprintf(qh, fp, 9319, "\nStatistics for: %s | %s",
1674 0 : qh->rbox_command, qh->qhull_command);
1675 0 : if (qh->ROTATErandom != INT_MIN)
1676 0 : qh_fprintf(qh, fp, 9320, " QR%d\n\n", qh->ROTATErandom);
1677 : else
1678 0 : qh_fprintf(qh, fp, 9321, "\n\n");
1679 0 : qh_fprintf(qh, fp, 9322, " Number of points processed: %d\n", zzval_(Zprocessed));
1680 0 : qh_fprintf(qh, fp, 9323, " Number of hyperplanes created: %d\n", zzval_(Zsetplane));
1681 0 : if (qh->DELAUNAY)
1682 0 : qh_fprintf(qh, fp, 9324, " Number of facets in hull: %d\n", qh->num_facets - qh->num_visible);
1683 0 : qh_fprintf(qh, fp, 9325, " Number of distance tests for qhull: %d\n", zzval_(Zpartition)+
1684 0 : zzval_(Zpartitionall)+zzval_(Znumvisibility)+zzval_(Zpartcoplanar));
1685 : #if 0 /* NOTE: must print before printstatistics() */
1686 : {realT stddev, ave;
1687 : qh_fprintf(qh, fp, 9326, " average new facet balance: %2.2g\n",
1688 : wval_(Wnewbalance)/zval_(Zprocessed));
1689 : stddev= qh_stddev(zval_(Zprocessed), wval_(Wnewbalance),
1690 : wval_(Wnewbalance2), &ave);
1691 : qh_fprintf(qh, fp, 9327, " new facet standard deviation: %2.2g\n", stddev);
1692 : qh_fprintf(qh, fp, 9328, " average partition balance: %2.2g\n",
1693 : wval_(Wpbalance)/zval_(Zpbalance));
1694 : stddev= qh_stddev(zval_(Zpbalance), wval_(Wpbalance),
1695 : wval_(Wpbalance2), &ave);
1696 : qh_fprintf(qh, fp, 9329, " partition standard deviation: %2.2g\n", stddev);
1697 : }
1698 : #endif
1699 0 : if (nummerged) {
1700 0 : qh_fprintf(qh, fp, 9330," Number of distance tests for merging: %d\n",zzval_(Zbestdist)+
1701 0 : zzval_(Zcentrumtests)+zzval_(Zvertextests)+zzval_(Zdistcheck)+zzval_(Zdistzero));
1702 0 : qh_fprintf(qh, fp, 9331," Number of distance tests for checking: %d\n",zzval_(Zcheckpart)+zzval_(Zdistconvex));
1703 0 : qh_fprintf(qh, fp, 9332," Number of merged facets: %d\n", nummerged);
1704 : }
1705 0 : numpinched= zzval_(Zpinchduplicate) + zzval_(Zpinchedvertex);
1706 0 : if (numpinched)
1707 0 : qh_fprintf(qh, fp, 9375," Number of merged pinched vertices: %d\n", numpinched);
1708 0 : if (!qh->RANDOMoutside && qh->QHULLfinished) {
1709 0 : cpu= (double)qh->hulltime;
1710 0 : cpu /= (double)qh_SECticks;
1711 0 : wval_(Wcpu)= cpu;
1712 0 : qh_fprintf(qh, fp, 9333, " CPU seconds to compute hull (after input): %2.4g\n", cpu);
1713 : }
1714 0 : if (qh->RERUN) {
1715 0 : if (!qh->PREmerge && !qh->MERGEexact)
1716 0 : qh_fprintf(qh, fp, 9334, " Percentage of runs with precision errors: %4.1f\n",
1717 0 : zzval_(Zretry)*100.0/qh->build_cnt); /* careful of order */
1718 0 : }else if (qh->JOGGLEmax < REALmax/2) {
1719 0 : if (zzval_(Zretry))
1720 0 : qh_fprintf(qh, fp, 9335, " After %d retries, input joggled by: %2.2g\n",
1721 : zzval_(Zretry), qh->JOGGLEmax);
1722 : else
1723 0 : qh_fprintf(qh, fp, 9336, " Input joggled by: %2.2g\n", qh->JOGGLEmax);
1724 : }
1725 0 : if (qh->totarea != 0.0)
1726 0 : qh_fprintf(qh, fp, 9337, " %s facet area: %2.8g\n",
1727 0 : zzval_(Ztotmerge) ? "Approximate" : "Total", qh->totarea);
1728 0 : if (qh->totvol != 0.0)
1729 0 : qh_fprintf(qh, fp, 9338, " %s volume: %2.8g\n",
1730 0 : zzval_(Ztotmerge) ? "Approximate" : "Total", qh->totvol);
1731 0 : if (qh->MERGING) {
1732 0 : qh_outerinner(qh, NULL, &outerplane, &innerplane);
1733 0 : if (outerplane > 2 * qh->DISTround) {
1734 0 : qh_fprintf(qh, fp, 9339, " Maximum distance of point above facet: %2.2g", outerplane);
1735 0 : ratio= outerplane/(qh->ONEmerge + qh->DISTround);
1736 : /* don't report ratio if MINoutside is large */
1737 0 : if (ratio > 0.05 && 2* qh->ONEmerge > qh->MINoutside && qh->JOGGLEmax > REALmax/2)
1738 0 : qh_fprintf(qh, fp, 9340, " (%.1fx)\n", ratio);
1739 : else
1740 0 : qh_fprintf(qh, fp, 9341, "\n");
1741 : }
1742 0 : if (innerplane < -2 * qh->DISTround) {
1743 0 : qh_fprintf(qh, fp, 9342, " Maximum distance of vertex below facet: %2.2g", innerplane);
1744 0 : ratio= -innerplane/(qh->ONEmerge+qh->DISTround);
1745 0 : if (ratio > 0.05 && qh->JOGGLEmax > REALmax/2)
1746 0 : qh_fprintf(qh, fp, 9343, " (%.1fx)\n", ratio);
1747 : else
1748 0 : qh_fprintf(qh, fp, 9344, "\n");
1749 : }
1750 : }
1751 0 : qh_fprintf(qh, fp, 9345, "\n");
1752 0 : } /* printsummary */
1753 :
1754 :
|
