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