1 /*<html><pre> -<a href="qh-geom_r.htm"
2 >-------------------------------</a><a name="TOP">-</a>
3
4 geom_r.c
5 geometric routines of qhull
6
7 see qh-geom_r.htm and geom_r.h
8
9 Copyright (c) 1993-2019 The Geometry Center.
10 $Id: //main/2019/qhull/src/libqhull_r/geom_r.c#4 $$Change: 2712 $
11 $DateTime: 2019/06/28 12:57:00 $$Author: bbarber $
12
13 infrequent code goes into geom2_r.c
14 */
15
16 #include "qhull_ra.h"
17
18 /*-<a href="qh-geom_r.htm#TOC"
19 >-------------------------------</a><a name="distplane">-</a>
20
21 qh_distplane(qh, point, facet, dist )
22 return distance from point to facet
23
24 returns:
25 dist
26 if qh.RANDOMdist, joggles result
27
28 notes:
29 dist > 0 if point is above facet (i.e., outside)
30 does not error (for qh_sortfacets, qh_outerinner)
31 for nearly coplanar points, the returned values may be duplicates
32 for example pairs of nearly incident points, rbox 175 C1,2e-13 t1538759579 | qhull d T4
33 622 qh_distplane: e-014 # count of two or more duplicate values for unique calls
34 258 qh_distplane: e-015
35 38 qh_distplane: e-016
36 40 qh_distplane: e-017
37 6 qh_distplane: e-018
38 5 qh_distplane: -e-018
39 33 qh_distplane: -e-017
40 3153 qh_distplane: -2.775557561562891e-017 # duplicated value for 3153 unique calls
41 42 qh_distplane: -e-016
42 307 qh_distplane: -e-015
43 1271 qh_distplane: -e-014
44 13 qh_distplane: -e-013
45
46 see:
47 qh_distnorm in geom2_r.c
48 qh_distplane [geom_r.c], QhullFacet::distance, and QhullHyperplane::distance are copies
49 */
qh_distplane(qhT * qh,pointT * point,facetT * facet,realT * dist)50 void qh_distplane(qhT *qh, pointT *point, facetT *facet, realT *dist) {
51 coordT *normal= facet->normal, *coordp, randr;
52 int k;
53
54 switch (qh->hull_dim){
55 case 2:
56 *dist= facet->offset + point[0] * normal[0] + point[1] * normal[1];
57 break;
58 case 3:
59 *dist= facet->offset + point[0] * normal[0] + point[1] * normal[1] + point[2] * normal[2];
60 break;
61 case 4:
62 *dist= facet->offset+point[0]*normal[0]+point[1]*normal[1]+point[2]*normal[2]+point[3]*normal[3];
63 break;
64 case 5:
65 *dist= facet->offset+point[0]*normal[0]+point[1]*normal[1]+point[2]*normal[2]+point[3]*normal[3]+point[4]*normal[4];
66 break;
67 case 6:
68 *dist= facet->offset+point[0]*normal[0]+point[1]*normal[1]+point[2]*normal[2]+point[3]*normal[3]+point[4]*normal[4]+point[5]*normal[5];
69 break;
70 case 7:
71 *dist= facet->offset+point[0]*normal[0]+point[1]*normal[1]+point[2]*normal[2]+point[3]*normal[3]+point[4]*normal[4]+point[5]*normal[5]+point[6]*normal[6];
72 break;
73 case 8:
74 *dist= facet->offset+point[0]*normal[0]+point[1]*normal[1]+point[2]*normal[2]+point[3]*normal[3]+point[4]*normal[4]+point[5]*normal[5]+point[6]*normal[6]+point[7]*normal[7];
75 break;
76 default:
77 *dist= facet->offset;
78 coordp= point;
79 for (k=qh->hull_dim; k--; )
80 *dist += *coordp++ * *normal++;
81 break;
82 }
83 zzinc_(Zdistplane);
84 if (!qh->RANDOMdist && qh->IStracing < 4)
85 return;
86 if (qh->RANDOMdist) {
87 randr= qh_RANDOMint;
88 *dist += (2.0 * randr / qh_RANDOMmax - 1.0) *
89 qh->RANDOMfactor * qh->MAXabs_coord;
90 }
91 #ifndef qh_NOtrace
92 if (qh->IStracing >= 4) {
93 qh_fprintf(qh, qh->ferr, 8001, "qh_distplane: ");
94 qh_fprintf(qh, qh->ferr, 8002, qh_REAL_1, *dist);
95 qh_fprintf(qh, qh->ferr, 8003, "from p%d to f%d\n", qh_pointid(qh, point), facet->id);
96 }
97 #endif
98 return;
99 } /* distplane */
100
101
102 /*-<a href="qh-geom_r.htm#TOC"
103 >-------------------------------</a><a name="findbest">-</a>
104
105 qh_findbest(qh, point, startfacet, bestoutside, qh_ISnewfacets, qh_NOupper, dist, isoutside, numpart )
106 find facet that is furthest below a point
107 for upperDelaunay facets
108 returns facet only if !qh_NOupper and clearly above
109
110 input:
111 starts search at 'startfacet' (can not be flipped)
112 if !bestoutside(qh_ALL), stops at qh.MINoutside
113
114 returns:
115 best facet (reports error if NULL)
116 early out if isoutside defined and bestdist > qh.MINoutside
117 dist is distance to facet
118 isoutside is true if point is outside of facet
119 numpart counts the number of distance tests
120
121 see also:
122 qh_findbestnew()
123
124 notes:
125 If merging (testhorizon), searches horizon facets of coplanar best facets because
126 after qh_distplane, this and qh_partitionpoint are the most expensive in 3-d
127 avoid calls to distplane, function calls, and real number operations.
128 caller traces result
129 Optimized for outside points. Tried recording a search set for qh_findhorizon.
130 Made code more complicated.
131
132 when called by qh_partitionvisible():
133 indicated by qh_ISnewfacets
134 qh.newfacet_list is list of simplicial, new facets
135 qh_findbestnew set if qh_sharpnewfacets returns True (to use qh_findbestnew)
136 qh.bestfacet_notsharp set if qh_sharpnewfacets returns False
137
138 when called by qh_findfacet(), qh_partitionpoint(), qh_partitioncoplanar(),
139 qh_check_bestdist(), qh_addpoint()
140 indicated by !qh_ISnewfacets
141 returns best facet in neighborhood of given facet
142 this is best facet overall if dist >= -qh.MAXcoplanar
143 or hull has at least a "spherical" curvature
144
145 design:
146 initialize and test for early exit
147 repeat while there are better facets
148 for each neighbor of facet
149 exit if outside facet found
150 test for better facet
151 if point is inside and partitioning
152 test for new facets with a "sharp" intersection
153 if so, future calls go to qh_findbestnew()
154 test horizon facets
155 */
qh_findbest(qhT * qh,pointT * point,facetT * startfacet,boolT bestoutside,boolT isnewfacets,boolT noupper,realT * dist,boolT * isoutside,int * numpart)156 facetT *qh_findbest(qhT *qh, pointT *point, facetT *startfacet,
157 boolT bestoutside, boolT isnewfacets, boolT noupper,
158 realT *dist, boolT *isoutside, int *numpart) {
159 realT bestdist= -REALmax/2 /* avoid underflow */;
160 facetT *facet, *neighbor, **neighborp;
161 facetT *bestfacet= NULL, *lastfacet= NULL;
162 int oldtrace= qh->IStracing;
163 unsigned int visitid= ++qh->visit_id;
164 int numpartnew=0;
165 boolT testhorizon= True; /* needed if precise, e.g., rbox c D6 | qhull Q0 Tv */
166
167 zinc_(Zfindbest);
168 #ifndef qh_NOtrace
169 if (qh->IStracing >= 4 || (qh->TRACElevel && qh->TRACEpoint >= 0 && qh->TRACEpoint == qh_pointid(qh, point))) {
170 if (qh->TRACElevel > qh->IStracing)
171 qh->IStracing= qh->TRACElevel;
172 qh_fprintf(qh, qh->ferr, 8004, "qh_findbest: point p%d starting at f%d isnewfacets? %d, unless %d exit if > %2.2g,",
173 qh_pointid(qh, point), startfacet->id, isnewfacets, bestoutside, qh->MINoutside);
174 qh_fprintf(qh, qh->ferr, 8005, " testhorizon? %d, noupper? %d,", testhorizon, noupper);
175 qh_fprintf(qh, qh->ferr, 8006, " Last qh_addpoint p%d,", qh->furthest_id);
176 qh_fprintf(qh, qh->ferr, 8007, " Last merge #%d, max_outside %2.2g\n", zzval_(Ztotmerge), qh->max_outside);
177 }
178 #endif
179 if (isoutside)
180 *isoutside= True;
181 if (!startfacet->flipped) { /* test startfacet before testing its neighbors */
182 *numpart= 1;
183 qh_distplane(qh, point, startfacet, dist); /* this code is duplicated below */
184 if (!bestoutside && *dist >= qh->MINoutside
185 && (!startfacet->upperdelaunay || !noupper)) {
186 bestfacet= startfacet;
187 goto LABELreturn_best;
188 }
189 bestdist= *dist;
190 if (!startfacet->upperdelaunay) {
191 bestfacet= startfacet;
192 }
193 }else
194 *numpart= 0;
195 startfacet->visitid= visitid;
196 facet= startfacet;
197 while (facet) {
198 trace4((qh, qh->ferr, 4001, "qh_findbest: neighbors of f%d, bestdist %2.2g f%d\n",
199 facet->id, bestdist, getid_(bestfacet)));
200 lastfacet= facet;
201 FOREACHneighbor_(facet) {
202 if (!neighbor->newfacet && isnewfacets)
203 continue;
204 if (neighbor->visitid == visitid)
205 continue;
206 neighbor->visitid= visitid;
207 if (!neighbor->flipped) { /* code duplicated above */
208 (*numpart)++;
209 qh_distplane(qh, point, neighbor, dist);
210 if (*dist > bestdist) {
211 if (!bestoutside && *dist >= qh->MINoutside
212 && (!neighbor->upperdelaunay || !noupper)) {
213 bestfacet= neighbor;
214 goto LABELreturn_best;
215 }
216 if (!neighbor->upperdelaunay) {
217 bestfacet= neighbor;
218 bestdist= *dist;
219 break; /* switch to neighbor */
220 }else if (!bestfacet) {
221 bestdist= *dist;
222 break; /* switch to neighbor */
223 }
224 } /* end of *dist>bestdist */
225 } /* end of !flipped */
226 } /* end of FOREACHneighbor */
227 facet= neighbor; /* non-NULL only if *dist>bestdist */
228 } /* end of while facet (directed search) */
229 if (isnewfacets) {
230 if (!bestfacet) { /* startfacet is upperdelaunay (or flipped) w/o !flipped newfacet neighbors */
231 bestdist= -REALmax/2;
232 bestfacet= qh_findbestnew(qh, point, qh->newfacet_list, &bestdist, bestoutside, isoutside, &numpartnew);
233 testhorizon= False; /* qh_findbestnew calls qh_findbesthorizon */
234 }else if (!qh->findbest_notsharp && bestdist < -qh->DISTround) {
235 if (qh_sharpnewfacets(qh)) {
236 /* seldom used, qh_findbestnew will retest all facets */
237 zinc_(Zfindnewsharp);
238 bestfacet= qh_findbestnew(qh, point, bestfacet, &bestdist, bestoutside, isoutside, &numpartnew);
239 testhorizon= False; /* qh_findbestnew calls qh_findbesthorizon */
240 qh->findbestnew= True;
241 }else
242 qh->findbest_notsharp= True;
243 }
244 }
245 if (!bestfacet)
246 bestfacet= qh_findbestlower(qh, lastfacet, point, &bestdist, numpart); /* lastfacet is non-NULL because startfacet is non-NULL */
247 if (testhorizon) /* qh_findbestnew not called */
248 bestfacet= qh_findbesthorizon(qh, !qh_IScheckmax, point, bestfacet, noupper, &bestdist, &numpartnew);
249 *dist= bestdist;
250 if (isoutside && bestdist < qh->MINoutside)
251 *isoutside= False;
252 LABELreturn_best:
253 zadd_(Zfindbesttot, *numpart);
254 zmax_(Zfindbestmax, *numpart);
255 (*numpart) += numpartnew;
256 qh->IStracing= oldtrace;
257 return bestfacet;
258 } /* findbest */
259
260
261 /*-<a href="qh-geom_r.htm#TOC"
262 >-------------------------------</a><a name="findbesthorizon">-</a>
263
264 qh_findbesthorizon(qh, qh_IScheckmax, point, startfacet, qh_NOupper, &bestdist, &numpart )
265 search coplanar and better horizon facets from startfacet/bestdist
266 ischeckmax turns off statistics and minsearch update
267 all arguments must be initialized, including *bestdist and *numpart
268 qh.coplanarfacetset used to maintain current search set, reset whenever best facet is substantially better
269 returns(ischeckmax):
270 best facet
271 updates f.maxoutside for neighbors of searched facets (if qh_MAXoutside)
272 returns(!ischeckmax):
273 best facet that is not upperdelaunay or newfacet (qh.first_newfacet)
274 allows upperdelaunay that is clearly outside
275 returns:
276 bestdist is distance to bestfacet
277 numpart -- updates number of distance tests
278
279 notes:
280 called by qh_findbest if point is not outside a facet (directed search)
281 called by qh_findbestnew if point is not outside a new facet
282 called by qh_check_maxout for each point in hull
283 called by qh_check_bestdist for each point in hull (rarely used)
284
285 no early out -- use qh_findbest() or qh_findbestnew()
286 Searches coplanar or better horizon facets
287
288 when called by qh_check_maxout() (qh_IScheckmax)
289 startfacet must be closest to the point
290 Otherwise, if point is beyond and below startfacet, startfacet may be a local minimum
291 even though other facets are below the point.
292 updates facet->maxoutside for good, visited facets
293 may return NULL
294
295 searchdist is qh.max_outside + 2 * DISTround
296 + max( MINvisible('Vn'), MAXcoplanar('Un'));
297 This setting is a guess. It must be at least max_outside + 2*DISTround
298 because a facet may have a geometric neighbor across a vertex
299
300 design:
301 for each horizon facet of coplanar best facets
302 continue if clearly inside
303 unless upperdelaunay or clearly outside
304 update best facet
305 */
qh_findbesthorizon(qhT * qh,boolT ischeckmax,pointT * point,facetT * startfacet,boolT noupper,realT * bestdist,int * numpart)306 facetT *qh_findbesthorizon(qhT *qh, boolT ischeckmax, pointT* point, facetT *startfacet, boolT noupper, realT *bestdist, int *numpart) {
307 facetT *bestfacet= startfacet;
308 realT dist;
309 facetT *neighbor, **neighborp, *facet;
310 facetT *nextfacet= NULL; /* optimize last facet of coplanarfacetset */
311 int numpartinit= *numpart, coplanarfacetset_size, numcoplanar= 0, numfacet= 0;
312 unsigned int visitid= ++qh->visit_id;
313 boolT newbest= False; /* for tracing */
314 realT minsearch, searchdist; /* skip facets that are too far from point */
315 boolT is_5x_minsearch;
316
317 if (!ischeckmax) {
318 zinc_(Zfindhorizon);
319 }else {
320 #if qh_MAXoutside
321 if ((!qh->ONLYgood || startfacet->good) && *bestdist > startfacet->maxoutside)
322 startfacet->maxoutside= *bestdist;
323 #endif
324 }
325 searchdist= qh_SEARCHdist; /* an expression, a multiple of qh.max_outside and precision constants */
326 minsearch= *bestdist - searchdist;
327 if (ischeckmax) {
328 /* Always check coplanar facets. Needed for RBOX 1000 s Z1 G1e-13 t996564279 | QHULL Tv */
329 minimize_(minsearch, -searchdist);
330 }
331 coplanarfacetset_size= 0;
332 startfacet->visitid= visitid;
333 facet= startfacet;
334 while (True) {
335 numfacet++;
336 is_5x_minsearch= (ischeckmax && facet->nummerge > 10 && qh_setsize(qh, facet->neighbors) > 100); /* QH11033 FIX: qh_findbesthorizon: many tests for facets with many merges and neighbors. Can hide coplanar facets, e.g., 'rbox 1000 s Z1 G1e-13' with 4400+ neighbors */
337 trace4((qh, qh->ferr, 4002, "qh_findbesthorizon: test neighbors of f%d bestdist %2.2g f%d ischeckmax? %d noupper? %d minsearch %2.2g is_5x? %d searchdist %2.2g\n",
338 facet->id, *bestdist, getid_(bestfacet), ischeckmax, noupper,
339 minsearch, is_5x_minsearch, searchdist));
340 FOREACHneighbor_(facet) {
341 if (neighbor->visitid == visitid)
342 continue;
343 neighbor->visitid= visitid;
344 if (!neighbor->flipped) { /* neighbors of flipped facets always searched via nextfacet */
345 qh_distplane(qh, point, neighbor, &dist); /* duplicate qh_distpane for new facets, they may be coplanar */
346 (*numpart)++;
347 if (dist > *bestdist) {
348 if (!neighbor->upperdelaunay || ischeckmax || (!noupper && dist >= qh->MINoutside)) {
349 if (!ischeckmax) {
350 minsearch= dist - searchdist;
351 if (dist > *bestdist + searchdist) {
352 zinc_(Zfindjump); /* everything in qh.coplanarfacetset at least searchdist below */
353 coplanarfacetset_size= 0;
354 }
355 }
356 bestfacet= neighbor;
357 *bestdist= dist;
358 newbest= True;
359 }
360 }else if (is_5x_minsearch) {
361 if (dist < 5 * minsearch)
362 continue; /* skip this neighbor, do not set nextfacet. dist is negative */
363 }else if (dist < minsearch)
364 continue; /* skip this neighbor, do not set nextfacet. If ischeckmax, dist can't be positive */
365 #if qh_MAXoutside
366 if (ischeckmax && dist > neighbor->maxoutside)
367 neighbor->maxoutside= dist;
368 #endif
369 } /* end of !flipped, need to search neighbor */
370 if (nextfacet) {
371 numcoplanar++;
372 if (!coplanarfacetset_size++) {
373 SETfirst_(qh->coplanarfacetset)= nextfacet;
374 SETtruncate_(qh->coplanarfacetset, 1);
375 }else
376 qh_setappend(qh, &qh->coplanarfacetset, nextfacet); /* Was needed for RBOX 1000 s W1e-13 P0 t996547055 | QHULL d Qbb Qc Tv
377 and RBOX 1000 s Z1 G1e-13 t996564279 | qhull Tv */
378 }
379 nextfacet= neighbor;
380 } /* end of EACHneighbor */
381 facet= nextfacet;
382 if (facet)
383 nextfacet= NULL;
384 else if (!coplanarfacetset_size)
385 break;
386 else if (!--coplanarfacetset_size) {
387 facet= SETfirstt_(qh->coplanarfacetset, facetT);
388 SETtruncate_(qh->coplanarfacetset, 0);
389 }else
390 facet= (facetT *)qh_setdellast(qh->coplanarfacetset);
391 } /* while True, i.e., "for each facet in qh.coplanarfacetset" */
392 if (!ischeckmax) {
393 zadd_(Zfindhorizontot, *numpart - numpartinit);
394 zmax_(Zfindhorizonmax, *numpart - numpartinit);
395 if (newbest)
396 zinc_(Znewbesthorizon);
397 }
398 trace4((qh, qh->ferr, 4003, "qh_findbesthorizon: p%d, newbest? %d, bestfacet f%d, bestdist %2.2g, numfacet %d, coplanarfacets %d, numdist %d\n",
399 qh_pointid(qh, point), newbest, getid_(bestfacet), *bestdist, numfacet, numcoplanar, *numpart - numpartinit));
400 return bestfacet;
401 } /* findbesthorizon */
402
403 /*-<a href="qh-geom_r.htm#TOC"
404 >-------------------------------</a><a name="findbestnew">-</a>
405
406 qh_findbestnew(qh, point, startfacet, dist, isoutside, numpart )
407 find best newfacet for point
408 searches all of qh.newfacet_list starting at startfacet
409 searches horizon facets of coplanar best newfacets
410 searches all facets if startfacet == qh.facet_list
411 returns:
412 best new or horizon facet that is not upperdelaunay
413 early out if isoutside and not 'Qf'
414 dist is distance to facet
415 isoutside is true if point is outside of facet
416 numpart is number of distance tests
417
418 notes:
419 Always used for merged new facets (see qh_USEfindbestnew)
420 Avoids upperdelaunay facet unless (isoutside and outside)
421
422 Uses qh.visit_id, qh.coplanarfacetset.
423 If share visit_id with qh_findbest, coplanarfacetset is incorrect.
424
425 If merging (testhorizon), searches horizon facets of coplanar best facets because
426 a point maybe coplanar to the bestfacet, below its horizon facet,
427 and above a horizon facet of a coplanar newfacet. For example,
428 rbox 1000 s Z1 G1e-13 | qhull
429 rbox 1000 s W1e-13 P0 t992110337 | QHULL d Qbb Qc
430
431 qh_findbestnew() used if
432 qh_sharpnewfacets -- newfacets contains a sharp angle
433 if many merges, qh_premerge found a merge, or 'Qf' (qh.findbestnew)
434
435 see also:
436 qh_partitionall() and qh_findbest()
437
438 design:
439 for each new facet starting from startfacet
440 test distance from point to facet
441 return facet if clearly outside
442 unless upperdelaunay and a lowerdelaunay exists
443 update best facet
444 test horizon facets
445 */
qh_findbestnew(qhT * qh,pointT * point,facetT * startfacet,realT * dist,boolT bestoutside,boolT * isoutside,int * numpart)446 facetT *qh_findbestnew(qhT *qh, pointT *point, facetT *startfacet,
447 realT *dist, boolT bestoutside, boolT *isoutside, int *numpart) {
448 realT bestdist= -REALmax/2;
449 facetT *bestfacet= NULL, *facet;
450 int oldtrace= qh->IStracing, i;
451 unsigned int visitid= ++qh->visit_id;
452 realT distoutside= 0.0;
453 boolT isdistoutside; /* True if distoutside is defined */
454 boolT testhorizon= True; /* needed if precise, e.g., rbox c D6 | qhull Q0 Tv */
455
456 if (!startfacet || !startfacet->next) {
457 if (qh->MERGING) {
458 qh_fprintf(qh, qh->ferr, 6001, "qhull topology error (qh_findbestnew): merging has formed and deleted a cone of new facets. Can not continue.\n");
459 qh_errexit(qh, qh_ERRtopology, NULL, NULL);
460 }else {
461 qh_fprintf(qh, qh->ferr, 6002, "qhull internal error (qh_findbestnew): no new facets for point p%d\n",
462 qh->furthest_id);
463 qh_errexit(qh, qh_ERRqhull, NULL, NULL);
464 }
465 }
466 zinc_(Zfindnew);
467 if (qh->BESToutside || bestoutside)
468 isdistoutside= False;
469 else {
470 isdistoutside= True;
471 distoutside= qh_DISToutside; /* multiple of qh.MINoutside & qh.max_outside, see user_r.h */
472 }
473 if (isoutside)
474 *isoutside= True;
475 *numpart= 0;
476 #ifndef qh_NOtrace
477 if (qh->IStracing >= 4 || (qh->TRACElevel && qh->TRACEpoint >= 0 && qh->TRACEpoint == qh_pointid(qh, point))) {
478 if (qh->TRACElevel > qh->IStracing)
479 qh->IStracing= qh->TRACElevel;
480 qh_fprintf(qh, qh->ferr, 8008, "qh_findbestnew: point p%d facet f%d. Stop? %d if dist > %2.2g,",
481 qh_pointid(qh, point), startfacet->id, isdistoutside, distoutside);
482 qh_fprintf(qh, qh->ferr, 8009, " Last qh_addpoint p%d, qh.visit_id %d, vertex_visit %d,", qh->furthest_id, visitid, qh->vertex_visit);
483 qh_fprintf(qh, qh->ferr, 8010, " Last merge #%d\n", zzval_(Ztotmerge));
484 }
485 #endif
486 /* visit all new facets starting with startfacet, maybe qh->facet_list */
487 for (i=0, facet=startfacet; i < 2; i++, facet= qh->newfacet_list) {
488 FORALLfacet_(facet) {
489 if (facet == startfacet && i)
490 break;
491 facet->visitid= visitid;
492 if (!facet->flipped) {
493 qh_distplane(qh, point, facet, dist);
494 (*numpart)++;
495 if (*dist > bestdist) {
496 if (!facet->upperdelaunay || *dist >= qh->MINoutside) {
497 bestfacet= facet;
498 if (isdistoutside && *dist >= distoutside)
499 goto LABELreturn_bestnew;
500 bestdist= *dist;
501 }
502 }
503 } /* end of !flipped */
504 } /* FORALLfacet from startfacet or qh->newfacet_list */
505 }
506 if (testhorizon || !bestfacet) /* testhorizon is always True. Keep the same code as qh_findbest */
507 bestfacet= qh_findbesthorizon(qh, !qh_IScheckmax, point, bestfacet ? bestfacet : startfacet,
508 !qh_NOupper, &bestdist, numpart);
509 *dist= bestdist;
510 if (isoutside && *dist < qh->MINoutside)
511 *isoutside= False;
512 LABELreturn_bestnew:
513 zadd_(Zfindnewtot, *numpart);
514 zmax_(Zfindnewmax, *numpart);
515 trace4((qh, qh->ferr, 4004, "qh_findbestnew: bestfacet f%d bestdist %2.2g for p%d f%d bestoutside? %d \n",
516 getid_(bestfacet), *dist, qh_pointid(qh, point), startfacet->id, bestoutside));
517 qh->IStracing= oldtrace;
518 return bestfacet;
519 } /* findbestnew */
520
521 /* ============ hyperplane functions -- keep code together [?] ============ */
522
523 /*-<a href="qh-geom_r.htm#TOC"
524 >-------------------------------</a><a name="backnormal">-</a>
525
526 qh_backnormal(qh, rows, numrow, numcol, sign, normal, nearzero )
527 given an upper-triangular rows array and a sign,
528 solve for normal equation x using back substitution over rows U
529
530 returns:
531 normal= x
532
533 if will not be able to divzero() when normalized(qh.MINdenom_2 and qh.MINdenom_1_2),
534 if fails on last row
535 this means that the hyperplane intersects [0,..,1]
536 sets last coordinate of normal to sign
537 otherwise
538 sets tail of normal to [...,sign,0,...], i.e., solves for b= [0...0]
539 sets nearzero
540
541 notes:
542 assumes numrow == numcol-1
543
544 see Golub & van Loan, 1983, Eq. 4.4-9 for "Gaussian elimination with complete pivoting"
545
546 solves Ux=b where Ax=b and PA=LU
547 b= [0,...,0,sign or 0] (sign is either -1 or +1)
548 last row of A= [0,...,0,1]
549
550 1) Ly=Pb == y=b since P only permutes the 0's of b
551
552 design:
553 for each row from end
554 perform back substitution
555 if near zero
556 use qh_divzero for division
557 if zero divide and not last row
558 set tail of normal to 0
559 */
qh_backnormal(qhT * qh,realT ** rows,int numrow,int numcol,boolT sign,coordT * normal,boolT * nearzero)560 void qh_backnormal(qhT *qh, realT **rows, int numrow, int numcol, boolT sign,
561 coordT *normal, boolT *nearzero) {
562 int i, j;
563 coordT *normalp, *normal_tail, *ai, *ak;
564 realT diagonal;
565 boolT waszero;
566 int zerocol= -1;
567
568 normalp= normal + numcol - 1;
569 *normalp--= (sign ? -1.0 : 1.0);
570 for (i=numrow; i--; ) {
571 *normalp= 0.0;
572 ai= rows[i] + i + 1;
573 ak= normalp+1;
574 for (j=i+1; j < numcol; j++)
575 *normalp -= *ai++ * *ak++;
576 diagonal= (rows[i])[i];
577 if (fabs_(diagonal) > qh->MINdenom_2)
578 *(normalp--) /= diagonal;
579 else {
580 waszero= False;
581 *normalp= qh_divzero(*normalp, diagonal, qh->MINdenom_1_2, &waszero);
582 if (waszero) {
583 zerocol= i;
584 *(normalp--)= (sign ? -1.0 : 1.0);
585 for (normal_tail= normalp+2; normal_tail < normal + numcol; normal_tail++)
586 *normal_tail= 0.0;
587 }else
588 normalp--;
589 }
590 }
591 if (zerocol != -1) {
592 *nearzero= True;
593 trace4((qh, qh->ferr, 4005, "qh_backnormal: zero diagonal at column %d.\n", i));
594 zzinc_(Zback0);
595 qh_joggle_restart(qh, "zero diagonal on back substitution");
596 }
597 } /* backnormal */
598
599 /*-<a href="qh-geom_r.htm#TOC"
600 >-------------------------------</a><a name="gausselim">-</a>
601
602 qh_gausselim(qh, rows, numrow, numcol, sign )
603 Gaussian elimination with partial pivoting
604
605 returns:
606 rows is upper triangular (includes row exchanges)
607 flips sign for each row exchange
608 sets nearzero if pivot[k] < qh.NEARzero[k], else clears it
609
610 notes:
611 if nearzero, the determinant's sign may be incorrect.
612 assumes numrow <= numcol
613
614 design:
615 for each row
616 determine pivot and exchange rows if necessary
617 test for near zero
618 perform gaussian elimination step
619 */
qh_gausselim(qhT * qh,realT ** rows,int numrow,int numcol,boolT * sign,boolT * nearzero)620 void qh_gausselim(qhT *qh, realT **rows, int numrow, int numcol, boolT *sign, boolT *nearzero) {
621 realT *ai, *ak, *rowp, *pivotrow;
622 realT n, pivot, pivot_abs= 0.0, temp;
623 int i, j, k, pivoti, flip=0;
624
625 *nearzero= False;
626 for (k=0; k < numrow; k++) {
627 pivot_abs= fabs_((rows[k])[k]);
628 pivoti= k;
629 for (i=k+1; i < numrow; i++) {
630 if ((temp= fabs_((rows[i])[k])) > pivot_abs) {
631 pivot_abs= temp;
632 pivoti= i;
633 }
634 }
635 if (pivoti != k) {
636 rowp= rows[pivoti];
637 rows[pivoti]= rows[k];
638 rows[k]= rowp;
639 *sign ^= 1;
640 flip ^= 1;
641 }
642 if (pivot_abs <= qh->NEARzero[k]) {
643 *nearzero= True;
644 if (pivot_abs == 0.0) { /* remainder of column == 0 */
645 #ifndef qh_NOtrace
646 if (qh->IStracing >= 4) {
647 qh_fprintf(qh, qh->ferr, 8011, "qh_gausselim: 0 pivot at column %d. (%2.2g < %2.2g)\n", k, pivot_abs, qh->DISTround);
648 qh_printmatrix(qh, qh->ferr, "Matrix:", rows, numrow, numcol);
649 }
650 #endif
651 zzinc_(Zgauss0);
652 qh_joggle_restart(qh, "zero pivot for Gaussian elimination");
653 goto LABELnextcol;
654 }
655 }
656 pivotrow= rows[k] + k;
657 pivot= *pivotrow++; /* signed value of pivot, and remainder of row */
658 for (i=k+1; i < numrow; i++) {
659 ai= rows[i] + k;
660 ak= pivotrow;
661 n= (*ai++)/pivot; /* divzero() not needed since |pivot| >= |*ai| */
662 for (j= numcol - (k+1); j--; )
663 *ai++ -= n * *ak++;
664 }
665 LABELnextcol:
666 ;
667 }
668 wmin_(Wmindenom, pivot_abs); /* last pivot element */
669 if (qh->IStracing >= 5)
670 qh_printmatrix(qh, qh->ferr, "qh_gausselem: result", rows, numrow, numcol);
671 } /* gausselim */
672
673
674 /*-<a href="qh-geom_r.htm#TOC"
675 >-------------------------------</a><a name="getangle">-</a>
676
677 qh_getangle(qh, vect1, vect2 )
678 returns the dot product of two vectors
679 if qh.RANDOMdist, joggles result
680
681 notes:
682 the angle may be > 1.0 or < -1.0 because of roundoff errors
683
684 */
qh_getangle(qhT * qh,pointT * vect1,pointT * vect2)685 realT qh_getangle(qhT *qh, pointT *vect1, pointT *vect2) {
686 realT angle= 0, randr;
687 int k;
688
689 for (k=qh->hull_dim; k--; )
690 angle += *vect1++ * *vect2++;
691 if (qh->RANDOMdist) {
692 randr= qh_RANDOMint;
693 angle += (2.0 * randr / qh_RANDOMmax - 1.0) *
694 qh->RANDOMfactor;
695 }
696 trace4((qh, qh->ferr, 4006, "qh_getangle: %4.4g\n", angle));
697 return(angle);
698 } /* getangle */
699
700
701 /*-<a href="qh-geom_r.htm#TOC"
702 >-------------------------------</a><a name="getcenter">-</a>
703
704 qh_getcenter(qh, vertices )
705 returns arithmetic center of a set of vertices as a new point
706
707 notes:
708 allocates point array for center
709 */
qh_getcenter(qhT * qh,setT * vertices)710 pointT *qh_getcenter(qhT *qh, setT *vertices) {
711 int k;
712 pointT *center, *coord;
713 vertexT *vertex, **vertexp;
714 int count= qh_setsize(qh, vertices);
715
716 if (count < 2) {
717 qh_fprintf(qh, qh->ferr, 6003, "qhull internal error (qh_getcenter): not defined for %d points\n", count);
718 qh_errexit(qh, qh_ERRqhull, NULL, NULL);
719 }
720 center= (pointT *)qh_memalloc(qh, qh->normal_size);
721 for (k=0; k < qh->hull_dim; k++) {
722 coord= center+k;
723 *coord= 0.0;
724 FOREACHvertex_(vertices)
725 *coord += vertex->point[k];
726 *coord /= count; /* count>=2 by QH6003 */
727 }
728 return(center);
729 } /* getcenter */
730
731
732 /*-<a href="qh-geom_r.htm#TOC"
733 >-------------------------------</a><a name="getcentrum">-</a>
734
735 qh_getcentrum(qh, facet )
736 returns the centrum for a facet as a new point
737
738 notes:
739 allocates the centrum
740 */
qh_getcentrum(qhT * qh,facetT * facet)741 pointT *qh_getcentrum(qhT *qh, facetT *facet) {
742 realT dist;
743 pointT *centrum, *point;
744
745 point= qh_getcenter(qh, facet->vertices);
746 zzinc_(Zcentrumtests);
747 qh_distplane(qh, point, facet, &dist);
748 centrum= qh_projectpoint(qh, point, facet, dist);
749 qh_memfree(qh, point, qh->normal_size);
750 trace4((qh, qh->ferr, 4007, "qh_getcentrum: for f%d, %d vertices dist= %2.2g\n",
751 facet->id, qh_setsize(qh, facet->vertices), dist));
752 return centrum;
753 } /* getcentrum */
754
755
756 /*-<a href="qh-geom_r.htm#TOC"
757 >-------------------------------</a><a name="getdistance">-</a>
758
759 qh_getdistance(qh, facet, neighbor, mindist, maxdist )
760 returns the min and max distance to neighbor of non-neighbor vertices in facet
761
762 returns:
763 the max absolute value
764
765 design:
766 for each vertex of facet that is not in neighbor
767 test the distance from vertex to neighbor
768 */
qh_getdistance(qhT * qh,facetT * facet,facetT * neighbor,coordT * mindist,coordT * maxdist)769 coordT qh_getdistance(qhT *qh, facetT *facet, facetT *neighbor, coordT *mindist, coordT *maxdist) {
770 vertexT *vertex, **vertexp;
771 coordT dist, maxd, mind;
772
773 FOREACHvertex_(facet->vertices)
774 vertex->seen= False;
775 FOREACHvertex_(neighbor->vertices)
776 vertex->seen= True;
777 mind= 0.0;
778 maxd= 0.0;
779 FOREACHvertex_(facet->vertices) {
780 if (!vertex->seen) {
781 zzinc_(Zbestdist);
782 qh_distplane(qh, vertex->point, neighbor, &dist);
783 if (dist < mind)
784 mind= dist;
785 else if (dist > maxd)
786 maxd= dist;
787 }
788 }
789 *mindist= mind;
790 *maxdist= maxd;
791 mind= -mind;
792 if (maxd > mind)
793 return maxd;
794 else
795 return mind;
796 } /* getdistance */
797
798
799 /*-<a href="qh-geom_r.htm#TOC"
800 >-------------------------------</a><a name="normalize">-</a>
801
802 qh_normalize(qh, normal, dim, toporient )
803 normalize a vector and report if too small
804 does not use min norm
805
806 see:
807 qh_normalize2
808 */
qh_normalize(qhT * qh,coordT * normal,int dim,boolT toporient)809 void qh_normalize(qhT *qh, coordT *normal, int dim, boolT toporient) {
810 qh_normalize2(qh, normal, dim, toporient, NULL, NULL);
811 } /* normalize */
812
813 /*-<a href="qh-geom_r.htm#TOC"
814 >-------------------------------</a><a name="normalize2">-</a>
815
816 qh_normalize2(qh, normal, dim, toporient, minnorm, ismin )
817 normalize a vector and report if too small
818 qh.MINdenom/MINdenom1 are the upper limits for divide overflow
819
820 returns:
821 normalized vector
822 flips sign if !toporient
823 if minnorm non-NULL,
824 sets ismin if normal < minnorm
825
826 notes:
827 if zero norm
828 sets all elements to sqrt(1.0/dim)
829 if divide by zero (divzero())
830 sets largest element to +/-1
831 bumps Znearlysingular
832
833 design:
834 computes norm
835 test for minnorm
836 if not near zero
837 normalizes normal
838 else if zero norm
839 sets normal to standard value
840 else
841 uses qh_divzero to normalize
842 if nearzero
843 sets norm to direction of maximum value
844 */
qh_normalize2(qhT * qh,coordT * normal,int dim,boolT toporient,realT * minnorm,boolT * ismin)845 void qh_normalize2(qhT *qh, coordT *normal, int dim, boolT toporient,
846 realT *minnorm, boolT *ismin) {
847 int k;
848 realT *colp, *maxp, norm= 0, temp, *norm1, *norm2, *norm3;
849 boolT zerodiv;
850
851 norm1= normal+1;
852 norm2= normal+2;
853 norm3= normal+3;
854 if (dim == 2)
855 norm= sqrt((*normal)*(*normal) + (*norm1)*(*norm1));
856 else if (dim == 3)
857 norm= sqrt((*normal)*(*normal) + (*norm1)*(*norm1) + (*norm2)*(*norm2));
858 else if (dim == 4) {
859 norm= sqrt((*normal)*(*normal) + (*norm1)*(*norm1) + (*norm2)*(*norm2)
860 + (*norm3)*(*norm3));
861 }else if (dim > 4) {
862 norm= (*normal)*(*normal) + (*norm1)*(*norm1) + (*norm2)*(*norm2)
863 + (*norm3)*(*norm3);
864 for (k=dim-4, colp=normal+4; k--; colp++)
865 norm += (*colp) * (*colp);
866 norm= sqrt(norm);
867 }
868 if (minnorm) {
869 if (norm < *minnorm)
870 *ismin= True;
871 else
872 *ismin= False;
873 }
874 wmin_(Wmindenom, norm);
875 if (norm > qh->MINdenom) {
876 if (!toporient)
877 norm= -norm;
878 *normal /= norm;
879 *norm1 /= norm;
880 if (dim == 2)
881 ; /* all done */
882 else if (dim == 3)
883 *norm2 /= norm;
884 else if (dim == 4) {
885 *norm2 /= norm;
886 *norm3 /= norm;
887 }else if (dim >4) {
888 *norm2 /= norm;
889 *norm3 /= norm;
890 for (k=dim-4, colp=normal+4; k--; )
891 *colp++ /= norm;
892 }
893 }else if (norm == 0.0) {
894 temp= sqrt(1.0/dim);
895 for (k=dim, colp=normal; k--; )
896 *colp++= temp;
897 }else {
898 if (!toporient)
899 norm= -norm;
900 for (k=dim, colp=normal; k--; colp++) { /* k used below */
901 temp= qh_divzero(*colp, norm, qh->MINdenom_1, &zerodiv);
902 if (!zerodiv)
903 *colp= temp;
904 else {
905 maxp= qh_maxabsval(normal, dim);
906 temp= ((*maxp * norm >= 0.0) ? 1.0 : -1.0);
907 for (k=dim, colp=normal; k--; colp++)
908 *colp= 0.0;
909 *maxp= temp;
910 zzinc_(Znearlysingular);
911 /* qh_joggle_restart ignored for Znearlysingular, normal part of qh_sethyperplane_gauss */
912 trace0((qh, qh->ferr, 1, "qh_normalize: norm=%2.2g too small during p%d\n",
913 norm, qh->furthest_id));
914 return;
915 }
916 }
917 }
918 } /* normalize */
919
920
921 /*-<a href="qh-geom_r.htm#TOC"
922 >-------------------------------</a><a name="projectpoint">-</a>
923
924 qh_projectpoint(qh, point, facet, dist )
925 project point onto a facet by dist
926
927 returns:
928 returns a new point
929
930 notes:
931 if dist= distplane(point,facet)
932 this projects point to hyperplane
933 assumes qh_memfree_() is valid for normal_size
934 */
qh_projectpoint(qhT * qh,pointT * point,facetT * facet,realT dist)935 pointT *qh_projectpoint(qhT *qh, pointT *point, facetT *facet, realT dist) {
936 pointT *newpoint, *np, *normal;
937 int normsize= qh->normal_size;
938 int k;
939 void **freelistp; /* used if !qh_NOmem by qh_memalloc_() */
940
941 qh_memalloc_(qh, normsize, freelistp, newpoint, pointT);
942 np= newpoint;
943 normal= facet->normal;
944 for (k=qh->hull_dim; k--; )
945 *(np++)= *point++ - dist * *normal++;
946 return(newpoint);
947 } /* projectpoint */
948
949
950 /*-<a href="qh-geom_r.htm#TOC"
951 >-------------------------------</a><a name="setfacetplane">-</a>
952
953 qh_setfacetplane(qh, facet )
954 sets the hyperplane for a facet
955 if qh.RANDOMdist, joggles hyperplane
956
957 notes:
958 uses global buffers qh.gm_matrix and qh.gm_row
959 overwrites facet->normal if already defined
960 updates Wnewvertex if PRINTstatistics
961 sets facet->upperdelaunay if upper envelope of Delaunay triangulation
962
963 design:
964 copy vertex coordinates to qh.gm_matrix/gm_row
965 compute determinate
966 if nearzero
967 recompute determinate with gaussian elimination
968 if nearzero
969 force outside orientation by testing interior point
970 */
qh_setfacetplane(qhT * qh,facetT * facet)971 void qh_setfacetplane(qhT *qh, facetT *facet) {
972 pointT *point;
973 vertexT *vertex, **vertexp;
974 int normsize= qh->normal_size;
975 int k,i, oldtrace= 0;
976 realT dist;
977 void **freelistp; /* used if !qh_NOmem by qh_memalloc_() */
978 coordT *coord, *gmcoord;
979 pointT *point0= SETfirstt_(facet->vertices, vertexT)->point;
980 boolT nearzero= False;
981
982 zzinc_(Zsetplane);
983 if (!facet->normal)
984 qh_memalloc_(qh, normsize, freelistp, facet->normal, coordT);
985 #ifndef qh_NOtrace
986 if (facet == qh->tracefacet) {
987 oldtrace= qh->IStracing;
988 qh->IStracing= 5;
989 qh_fprintf(qh, qh->ferr, 8012, "qh_setfacetplane: facet f%d created.\n", facet->id);
990 qh_fprintf(qh, qh->ferr, 8013, " Last point added to hull was p%d.", qh->furthest_id);
991 if (zzval_(Ztotmerge))
992 qh_fprintf(qh, qh->ferr, 8014, " Last merge was #%d.", zzval_(Ztotmerge));
993 qh_fprintf(qh, qh->ferr, 8015, "\n\nCurrent summary is:\n");
994 qh_printsummary(qh, qh->ferr);
995 }
996 #endif
997 if (qh->hull_dim <= 4) {
998 i= 0;
999 if (qh->RANDOMdist) {
1000 gmcoord= qh->gm_matrix;
1001 FOREACHvertex_(facet->vertices) {
1002 qh->gm_row[i++]= gmcoord;
1003 coord= vertex->point;
1004 for (k=qh->hull_dim; k--; )
1005 *(gmcoord++)= *coord++ * qh_randomfactor(qh, qh->RANDOMa, qh->RANDOMb);
1006 }
1007 }else {
1008 FOREACHvertex_(facet->vertices)
1009 qh->gm_row[i++]= vertex->point;
1010 }
1011 qh_sethyperplane_det(qh, qh->hull_dim, qh->gm_row, point0, facet->toporient,
1012 facet->normal, &facet->offset, &nearzero);
1013 }
1014 if (qh->hull_dim > 4 || nearzero) {
1015 i= 0;
1016 gmcoord= qh->gm_matrix;
1017 FOREACHvertex_(facet->vertices) {
1018 if (vertex->point != point0) {
1019 qh->gm_row[i++]= gmcoord;
1020 coord= vertex->point;
1021 point= point0;
1022 for (k=qh->hull_dim; k--; )
1023 *(gmcoord++)= *coord++ - *point++;
1024 }
1025 }
1026 qh->gm_row[i]= gmcoord; /* for areasimplex */
1027 if (qh->RANDOMdist) {
1028 gmcoord= qh->gm_matrix;
1029 for (i=qh->hull_dim-1; i--; ) {
1030 for (k=qh->hull_dim; k--; )
1031 *(gmcoord++) *= qh_randomfactor(qh, qh->RANDOMa, qh->RANDOMb);
1032 }
1033 }
1034 qh_sethyperplane_gauss(qh, qh->hull_dim, qh->gm_row, point0, facet->toporient,
1035 facet->normal, &facet->offset, &nearzero);
1036 if (nearzero) {
1037 if (qh_orientoutside(qh, facet)) {
1038 trace0((qh, qh->ferr, 2, "qh_setfacetplane: flipped orientation due to nearzero gauss and interior_point test. During p%d\n", qh->furthest_id));
1039 /* this is part of using Gaussian Elimination. For example in 5-d
1040 1 1 1 1 0
1041 1 1 1 1 1
1042 0 0 0 1 0
1043 0 1 0 0 0
1044 1 0 0 0 0
1045 norm= 0.38 0.38 -0.76 0.38 0
1046 has a determinate of 1, but g.e. after subtracting pt. 0 has
1047 0's in the diagonal, even with full pivoting. It does work
1048 if you subtract pt. 4 instead. */
1049 }
1050 }
1051 }
1052 facet->upperdelaunay= False;
1053 if (qh->DELAUNAY) {
1054 if (qh->UPPERdelaunay) { /* matches qh_triangulate_facet and qh.lower_threshold in qh_initbuild */
1055 if (facet->normal[qh->hull_dim -1] >= qh->ANGLEround * qh_ZEROdelaunay)
1056 facet->upperdelaunay= True;
1057 }else {
1058 if (facet->normal[qh->hull_dim -1] > -qh->ANGLEround * qh_ZEROdelaunay)
1059 facet->upperdelaunay= True;
1060 }
1061 }
1062 if (qh->PRINTstatistics || qh->IStracing || qh->TRACElevel || qh->JOGGLEmax < REALmax) {
1063 qh->old_randomdist= qh->RANDOMdist;
1064 qh->RANDOMdist= False;
1065 FOREACHvertex_(facet->vertices) {
1066 if (vertex->point != point0) {
1067 boolT istrace= False;
1068 zinc_(Zdiststat);
1069 qh_distplane(qh, vertex->point, facet, &dist);
1070 dist= fabs_(dist);
1071 zinc_(Znewvertex);
1072 wadd_(Wnewvertex, dist);
1073 if (dist > wwval_(Wnewvertexmax)) {
1074 wwval_(Wnewvertexmax)= dist;
1075 if (dist > qh->max_outside) {
1076 qh->max_outside= dist; /* used by qh_maxouter(qh) */
1077 if (dist > qh->TRACEdist)
1078 istrace= True;
1079 }
1080 }else if (-dist > qh->TRACEdist)
1081 istrace= True;
1082 if (istrace) {
1083 qh_fprintf(qh, qh->ferr, 3060, "qh_setfacetplane: ====== vertex p%d(v%d) increases max_outside to %2.2g for new facet f%d last p%d\n",
1084 qh_pointid(qh, vertex->point), vertex->id, dist, facet->id, qh->furthest_id);
1085 qh_errprint(qh, "DISTANT", facet, NULL, NULL, NULL);
1086 }
1087 }
1088 }
1089 qh->RANDOMdist= qh->old_randomdist;
1090 }
1091 #ifndef qh_NOtrace
1092 if (qh->IStracing >= 4) {
1093 qh_fprintf(qh, qh->ferr, 8017, "qh_setfacetplane: f%d offset %2.2g normal: ",
1094 facet->id, facet->offset);
1095 for (k=0; k < qh->hull_dim; k++)
1096 qh_fprintf(qh, qh->ferr, 8018, "%2.2g ", facet->normal[k]);
1097 qh_fprintf(qh, qh->ferr, 8019, "\n");
1098 }
1099 #endif
1100 qh_checkflipped(qh, facet, NULL, qh_ALL);
1101 if (facet == qh->tracefacet) {
1102 qh->IStracing= oldtrace;
1103 qh_printfacet(qh, qh->ferr, facet);
1104 }
1105 } /* setfacetplane */
1106
1107
1108 /*-<a href="qh-geom_r.htm#TOC"
1109 >-------------------------------</a><a name="sethyperplane_det">-</a>
1110
1111 qh_sethyperplane_det(qh, dim, rows, point0, toporient, normal, offset, nearzero )
1112 given dim X dim array indexed by rows[], one row per point,
1113 toporient(flips all signs),
1114 and point0 (any row)
1115 set normalized hyperplane equation from oriented simplex
1116
1117 returns:
1118 normal (normalized)
1119 offset (places point0 on the hyperplane)
1120 sets nearzero if hyperplane not through points
1121
1122 notes:
1123 only defined for dim == 2..4
1124 rows[] is not modified
1125 solves det(P-V_0, V_n-V_0, ..., V_1-V_0)=0, i.e. every point is on hyperplane
1126 see Bower & Woodworth, A programmer's geometry, Butterworths 1983.
1127
1128 derivation of 3-d minnorm
1129 Goal: all vertices V_i within qh.one_merge of hyperplane
1130 Plan: exactly translate the facet so that V_0 is the origin
1131 exactly rotate the facet so that V_1 is on the x-axis and y_2=0.
1132 exactly rotate the effective perturbation to only effect n_0
1133 this introduces a factor of sqrt(3)
1134 n_0 = ((y_2-y_0)*(z_1-z_0) - (z_2-z_0)*(y_1-y_0)) / norm
1135 Let M_d be the max coordinate difference
1136 Let M_a be the greater of M_d and the max abs. coordinate
1137 Let u be machine roundoff and distround be max error for distance computation
1138 The max error for n_0 is sqrt(3) u M_a M_d / norm. n_1 is approx. 1 and n_2 is approx. 0
1139 The max error for distance of V_1 is sqrt(3) u M_a M_d M_d / norm. Offset=0 at origin
1140 Then minnorm = 1.8 u M_a M_d M_d / qh.ONEmerge
1141 Note that qh.one_merge is approx. 45.5 u M_a and norm is usually about M_d M_d
1142
1143 derivation of 4-d minnorm
1144 same as above except rotate the facet so that V_1 on x-axis and w_2, y_3, w_3=0
1145 [if two vertices fixed on x-axis, can rotate the other two in yzw.]
1146 n_0 = det3_(...) = y_2 det2_(z_1, w_1, z_3, w_3) = - y_2 w_1 z_3
1147 [all other terms contain at least two factors nearly zero.]
1148 The max error for n_0 is sqrt(4) u M_a M_d M_d / norm
1149 Then minnorm = 2 u M_a M_d M_d M_d / qh.ONEmerge
1150 Note that qh.one_merge is approx. 82 u M_a and norm is usually about M_d M_d M_d
1151 */
qh_sethyperplane_det(qhT * qh,int dim,coordT ** rows,coordT * point0,boolT toporient,coordT * normal,realT * offset,boolT * nearzero)1152 void qh_sethyperplane_det(qhT *qh, int dim, coordT **rows, coordT *point0,
1153 boolT toporient, coordT *normal, realT *offset, boolT *nearzero) {
1154 realT maxround, dist;
1155 int i;
1156 pointT *point;
1157
1158
1159 if (dim == 2) {
1160 normal[0]= dY(1,0);
1161 normal[1]= dX(0,1);
1162 qh_normalize2(qh, normal, dim, toporient, NULL, NULL);
1163 *offset= -(point0[0]*normal[0]+point0[1]*normal[1]);
1164 *nearzero= False; /* since nearzero norm => incident points */
1165 }else if (dim == 3) {
1166 normal[0]= det2_(dY(2,0), dZ(2,0),
1167 dY(1,0), dZ(1,0));
1168 normal[1]= det2_(dX(1,0), dZ(1,0),
1169 dX(2,0), dZ(2,0));
1170 normal[2]= det2_(dX(2,0), dY(2,0),
1171 dX(1,0), dY(1,0));
1172 qh_normalize2(qh, normal, dim, toporient, NULL, NULL);
1173 *offset= -(point0[0]*normal[0] + point0[1]*normal[1]
1174 + point0[2]*normal[2]);
1175 maxround= qh->DISTround;
1176 for (i=dim; i--; ) {
1177 point= rows[i];
1178 if (point != point0) {
1179 dist= *offset + (point[0]*normal[0] + point[1]*normal[1]
1180 + point[2]*normal[2]);
1181 if (dist > maxround || dist < -maxround) {
1182 *nearzero= True;
1183 break;
1184 }
1185 }
1186 }
1187 }else if (dim == 4) {
1188 normal[0]= - det3_(dY(2,0), dZ(2,0), dW(2,0),
1189 dY(1,0), dZ(1,0), dW(1,0),
1190 dY(3,0), dZ(3,0), dW(3,0));
1191 normal[1]= det3_(dX(2,0), dZ(2,0), dW(2,0),
1192 dX(1,0), dZ(1,0), dW(1,0),
1193 dX(3,0), dZ(3,0), dW(3,0));
1194 normal[2]= - det3_(dX(2,0), dY(2,0), dW(2,0),
1195 dX(1,0), dY(1,0), dW(1,0),
1196 dX(3,0), dY(3,0), dW(3,0));
1197 normal[3]= det3_(dX(2,0), dY(2,0), dZ(2,0),
1198 dX(1,0), dY(1,0), dZ(1,0),
1199 dX(3,0), dY(3,0), dZ(3,0));
1200 qh_normalize2(qh, normal, dim, toporient, NULL, NULL);
1201 *offset= -(point0[0]*normal[0] + point0[1]*normal[1]
1202 + point0[2]*normal[2] + point0[3]*normal[3]);
1203 maxround= qh->DISTround;
1204 for (i=dim; i--; ) {
1205 point= rows[i];
1206 if (point != point0) {
1207 dist= *offset + (point[0]*normal[0] + point[1]*normal[1]
1208 + point[2]*normal[2] + point[3]*normal[3]);
1209 if (dist > maxround || dist < -maxround) {
1210 *nearzero= True;
1211 break;
1212 }
1213 }
1214 }
1215 }
1216 if (*nearzero) {
1217 zzinc_(Zminnorm);
1218 /* qh_joggle_restart not needed, will call qh_sethyperplane_gauss instead */
1219 trace0((qh, qh->ferr, 3, "qh_sethyperplane_det: degenerate norm during p%d, use qh_sethyperplane_gauss instead.\n", qh->furthest_id));
1220 }
1221 } /* sethyperplane_det */
1222
1223
1224 /*-<a href="qh-geom_r.htm#TOC"
1225 >-------------------------------</a><a name="sethyperplane_gauss">-</a>
1226
1227 qh_sethyperplane_gauss(qh, dim, rows, point0, toporient, normal, offset, nearzero )
1228 given(dim-1) X dim array of rows[i]= V_{i+1} - V_0 (point0)
1229 set normalized hyperplane equation from oriented simplex
1230
1231 returns:
1232 normal (normalized)
1233 offset (places point0 on the hyperplane)
1234
1235 notes:
1236 if nearzero
1237 orientation may be incorrect because of incorrect sign flips in gausselim
1238 solves [V_n-V_0,...,V_1-V_0, 0 .. 0 1] * N == [0 .. 0 1]
1239 or [V_n-V_0,...,V_1-V_0, 0 .. 0 1] * N == [0]
1240 i.e., N is normal to the hyperplane, and the unnormalized
1241 distance to [0 .. 1] is either 1 or 0
1242
1243 design:
1244 perform gaussian elimination
1245 flip sign for negative values
1246 perform back substitution
1247 normalize result
1248 compute offset
1249 */
qh_sethyperplane_gauss(qhT * qh,int dim,coordT ** rows,pointT * point0,boolT toporient,coordT * normal,coordT * offset,boolT * nearzero)1250 void qh_sethyperplane_gauss(qhT *qh, int dim, coordT **rows, pointT *point0,
1251 boolT toporient, coordT *normal, coordT *offset, boolT *nearzero) {
1252 coordT *pointcoord, *normalcoef;
1253 int k;
1254 boolT sign= toporient, nearzero2= False;
1255
1256 qh_gausselim(qh, rows, dim-1, dim, &sign, nearzero);
1257 for (k=dim-1; k--; ) {
1258 if ((rows[k])[k] < 0)
1259 sign ^= 1;
1260 }
1261 if (*nearzero) {
1262 zzinc_(Znearlysingular);
1263 /* qh_joggle_restart ignored for Znearlysingular, normal part of qh_sethyperplane_gauss */
1264 trace0((qh, qh->ferr, 4, "qh_sethyperplane_gauss: nearly singular or axis parallel hyperplane during p%d.\n", qh->furthest_id));
1265 qh_backnormal(qh, rows, dim-1, dim, sign, normal, &nearzero2);
1266 }else {
1267 qh_backnormal(qh, rows, dim-1, dim, sign, normal, &nearzero2);
1268 if (nearzero2) {
1269 zzinc_(Znearlysingular);
1270 trace0((qh, qh->ferr, 5, "qh_sethyperplane_gauss: singular or axis parallel hyperplane at normalization during p%d.\n", qh->furthest_id));
1271 }
1272 }
1273 if (nearzero2)
1274 *nearzero= True;
1275 qh_normalize2(qh, normal, dim, True, NULL, NULL);
1276 pointcoord= point0;
1277 normalcoef= normal;
1278 *offset= -(*pointcoord++ * *normalcoef++);
1279 for (k=dim-1; k--; )
1280 *offset -= *pointcoord++ * *normalcoef++;
1281 } /* sethyperplane_gauss */
1282
1283
1284
1285