Line data Source code
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-2020 The Geometry Center.
10 : $Id: //main/2019/qhull/src/libqhull_r/geom_r.c#5 $$Change: 2953 $
11 : $DateTime: 2020/05/21 22:05:32 $$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 : */
50 1324710 : void qh_distplane(qhT *qh, pointT *point, facetT *facet, realT *dist) {
51 1324710 : coordT *normal= facet->normal, *coordp, randr;
52 : int k;
53 :
54 1324710 : switch (qh->hull_dim){
55 0 : case 2:
56 0 : *dist= facet->offset + point[0] * normal[0] + point[1] * normal[1];
57 0 : break;
58 1324710 : case 3:
59 1324710 : *dist= facet->offset + point[0] * normal[0] + point[1] * normal[1] + point[2] * normal[2];
60 1324710 : break;
61 0 : case 4:
62 0 : *dist= facet->offset+point[0]*normal[0]+point[1]*normal[1]+point[2]*normal[2]+point[3]*normal[3];
63 0 : break;
64 0 : case 5:
65 0 : *dist= facet->offset+point[0]*normal[0]+point[1]*normal[1]+point[2]*normal[2]+point[3]*normal[3]+point[4]*normal[4];
66 0 : break;
67 0 : case 6:
68 0 : *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 0 : break;
70 0 : case 7:
71 0 : *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 0 : break;
73 0 : case 8:
74 0 : *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 0 : break;
76 0 : default:
77 0 : *dist= facet->offset;
78 0 : coordp= point;
79 0 : for (k=qh->hull_dim; k--; )
80 0 : *dist += *coordp++ * *normal++;
81 0 : break;
82 : }
83 1324710 : zzinc_(Zdistplane);
84 1324710 : if (!qh->RANDOMdist && qh->IStracing < 4)
85 1324710 : return;
86 0 : if (qh->RANDOMdist) {
87 0 : randr= qh_RANDOMint;
88 0 : *dist += (2.0 * randr / qh_RANDOMmax - 1.0) *
89 0 : qh->RANDOMfactor * qh->MAXabs_coord;
90 : }
91 : #ifndef qh_NOtrace
92 0 : if (qh->IStracing >= 4) {
93 0 : qh_fprintf(qh, qh->ferr, 8001, "qh_distplane: ");
94 0 : qh_fprintf(qh, qh->ferr, 8002, qh_REAL_1, *dist);
95 0 : qh_fprintf(qh, qh->ferr, 8003, "from p%d to f%d\n", qh_pointid(qh, point), facet->id);
96 : }
97 : #endif
98 0 : 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 : */
156 23311 : facetT *qh_findbest(qhT *qh, pointT *point, facetT *startfacet,
157 : boolT bestoutside, boolT isnewfacets, boolT noupper,
158 : realT *dist, boolT *isoutside, int *numpart) {
159 23311 : realT bestdist= -REALmax/2 /* avoid underflow */;
160 : facetT *facet, *neighbor, **neighborp;
161 23311 : facetT *bestfacet= NULL, *lastfacet= NULL;
162 23311 : int oldtrace= qh->IStracing;
163 23311 : unsigned int visitid= ++qh->visit_id;
164 23311 : int numpartnew=0;
165 23311 : boolT testhorizon= True; /* needed if precise, e.g., rbox c D6 | qhull Q0 Tv */
166 :
167 23311 : zinc_(Zfindbest);
168 : #ifndef qh_NOtrace
169 23311 : if (qh->IStracing >= 4 || (qh->TRACElevel && qh->TRACEpoint >= 0 && qh->TRACEpoint == qh_pointid(qh, point))) {
170 0 : if (qh->TRACElevel > qh->IStracing)
171 0 : qh->IStracing= qh->TRACElevel;
172 0 : 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 0 : qh_fprintf(qh, qh->ferr, 8005, " testhorizon? %d, noupper? %d,", testhorizon, noupper);
175 0 : qh_fprintf(qh, qh->ferr, 8006, " Last qh_addpoint p%d,", qh->furthest_id);
176 0 : qh_fprintf(qh, qh->ferr, 8007, " Last merge #%d, max_outside %2.2g\n", zzval_(Ztotmerge), qh->max_outside);
177 : }
178 : #endif
179 23311 : if (isoutside)
180 23311 : *isoutside= True;
181 23311 : if (!startfacet->flipped) { /* test startfacet before testing its neighbors */
182 23311 : *numpart= 1;
183 23311 : qh_distplane(qh, point, startfacet, dist); /* this code is duplicated below */
184 23311 : if (!bestoutside && *dist >= qh->MINoutside
185 10386 : && (!startfacet->upperdelaunay || !noupper)) {
186 10386 : bestfacet= startfacet;
187 10386 : goto LABELreturn_best;
188 : }
189 12925 : bestdist= *dist;
190 12925 : if (!startfacet->upperdelaunay) {
191 12920 : bestfacet= startfacet;
192 : }
193 : }else
194 0 : *numpart= 0;
195 12925 : startfacet->visitid= visitid;
196 12925 : facet= startfacet;
197 16772 : while (facet) {
198 16769 : trace4((qh, qh->ferr, 4001, "qh_findbest: neighbors of f%d, bestdist %2.2g f%d\n",
199 : facet->id, bestdist, getid_(bestfacet)));
200 16769 : lastfacet= facet;
201 38719 : FOREACHneighbor_(facet) {
202 38716 : if (!neighbor->newfacet && isnewfacets)
203 16767 : continue;
204 21949 : if (neighbor->visitid == visitid)
205 1173 : continue;
206 20776 : neighbor->visitid= visitid;
207 20776 : if (!neighbor->flipped) { /* code duplicated above */
208 20776 : (*numpart)++;
209 20776 : qh_distplane(qh, point, neighbor, dist);
210 20776 : if (*dist > bestdist) {
211 16825 : if (!bestoutside && *dist >= qh->MINoutside
212 12922 : && (!neighbor->upperdelaunay || !noupper)) {
213 12922 : bestfacet= neighbor;
214 12922 : goto LABELreturn_best;
215 : }
216 3903 : if (!neighbor->upperdelaunay) {
217 3843 : bestfacet= neighbor;
218 3843 : bestdist= *dist;
219 3843 : break; /* switch to neighbor */
220 60 : }else if (!bestfacet) {
221 1 : bestdist= *dist;
222 1 : break; /* switch to neighbor */
223 : }
224 : } /* end of *dist>bestdist */
225 : } /* end of !flipped */
226 : } /* end of FOREACHneighbor */
227 3847 : facet= neighbor; /* non-NULL only if *dist>bestdist */
228 : } /* end of while facet (directed search) */
229 3 : if (isnewfacets) {
230 1 : if (!bestfacet) { /* startfacet is upperdelaunay (or flipped) w/o !flipped newfacet neighbors */
231 0 : bestdist= -REALmax/2;
232 0 : bestfacet= qh_findbestnew(qh, point, qh->newfacet_list, &bestdist, bestoutside, isoutside, &numpartnew);
233 0 : testhorizon= False; /* qh_findbestnew calls qh_findbesthorizon */
234 1 : }else if (!qh->findbest_notsharp && bestdist < -qh->DISTround) {
235 0 : if (qh_sharpnewfacets(qh)) {
236 : /* seldom used, qh_findbestnew will retest all facets */
237 0 : zinc_(Zfindnewsharp);
238 0 : bestfacet= qh_findbestnew(qh, point, bestfacet, &bestdist, bestoutside, isoutside, &numpartnew);
239 0 : testhorizon= False; /* qh_findbestnew calls qh_findbesthorizon */
240 0 : qh->findbestnew= True;
241 : }else
242 0 : qh->findbest_notsharp= True;
243 : }
244 : }
245 3 : if (!bestfacet)
246 0 : bestfacet= qh_findbestlower(qh, lastfacet, point, &bestdist, numpart); /* lastfacet is non-NULL because startfacet is non-NULL */
247 3 : if (testhorizon) /* qh_findbestnew not called */
248 3 : bestfacet= qh_findbesthorizon(qh, !qh_IScheckmax, point, bestfacet, noupper, &bestdist, &numpartnew);
249 3 : *dist= bestdist;
250 3 : if (isoutside && bestdist < qh->MINoutside)
251 3 : *isoutside= False;
252 0 : LABELreturn_best:
253 23311 : zadd_(Zfindbesttot, *numpart);
254 23311 : zmax_(Zfindbestmax, *numpart);
255 23311 : (*numpart) += numpartnew;
256 23311 : qh->IStracing= oldtrace;
257 23311 : 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 : */
306 14677 : facetT *qh_findbesthorizon(qhT *qh, boolT ischeckmax, pointT* point, facetT *startfacet, boolT noupper, realT *bestdist, int *numpart) {
307 14677 : facetT *bestfacet= startfacet;
308 : realT dist;
309 : facetT *neighbor, **neighborp, *facet;
310 14677 : facetT *nextfacet= NULL; /* optimize last facet of coplanarfacetset */
311 14677 : int numpartinit= *numpart, coplanarfacetset_size, numcoplanar= 0, numfacet= 0;
312 14677 : unsigned int visitid= ++qh->visit_id;
313 14677 : boolT newbest= False; /* for tracing */
314 : realT minsearch, searchdist; /* skip facets that are too far from point */
315 : boolT is_5x_minsearch;
316 :
317 14677 : if (!ischeckmax) {
318 4 : zinc_(Zfindhorizon);
319 : }else {
320 : #if qh_MAXoutside
321 14673 : if ((!qh->ONLYgood || startfacet->good) && *bestdist > startfacet->maxoutside)
322 0 : startfacet->maxoutside= *bestdist;
323 : #endif
324 : }
325 14677 : searchdist= qh_SEARCHdist; /* an expression, a multiple of qh.max_outside and precision constants */
326 14677 : minsearch= *bestdist - searchdist;
327 14677 : if (ischeckmax) {
328 : /* Always check coplanar facets. Needed for RBOX 1000 s Z1 G1e-13 t996564279 | QHULL Tv */
329 14673 : minimize_(minsearch, -searchdist);
330 : }
331 14677 : coplanarfacetset_size= 0;
332 14677 : startfacet->visitid= visitid;
333 14677 : facet= startfacet;
334 : while (True) {
335 58266 : numfacet++;
336 58266 : 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 58266 : 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 347782 : FOREACHneighbor_(facet) {
341 289516 : if (neighbor->visitid == visitid)
342 74378 : continue;
343 215138 : neighbor->visitid= visitid;
344 215138 : if (!neighbor->flipped) { /* neighbors of flipped facets always searched via nextfacet */
345 215138 : qh_distplane(qh, point, neighbor, &dist); /* duplicate qh_distpane for new facets, they may be coplanar */
346 215138 : (*numpart)++;
347 215138 : if (dist > *bestdist) {
348 15948 : if (!neighbor->upperdelaunay || ischeckmax || (!noupper && dist >= qh->MINoutside)) {
349 15946 : if (!ischeckmax) {
350 0 : minsearch= dist - searchdist;
351 0 : if (dist > *bestdist + searchdist) {
352 0 : zinc_(Zfindjump); /* everything in qh.coplanarfacetset at least searchdist below */
353 0 : coplanarfacetset_size= 0;
354 : }
355 : }
356 15946 : bestfacet= neighbor;
357 15946 : *bestdist= dist;
358 15946 : newbest= True;
359 : }
360 199190 : }else if (is_5x_minsearch) {
361 56807 : if (dist < 5 * minsearch)
362 56783 : continue; /* skip this neighbor, do not set nextfacet. dist is negative */
363 142383 : }else if (dist < minsearch)
364 114766 : continue; /* skip this neighbor, do not set nextfacet. If ischeckmax, dist can't be positive */
365 : #if qh_MAXoutside
366 43589 : if (ischeckmax && dist > neighbor->maxoutside)
367 0 : neighbor->maxoutside= dist;
368 : #endif
369 : } /* end of !flipped, need to search neighbor */
370 43589 : if (nextfacet) {
371 14678 : numcoplanar++;
372 14678 : if (!coplanarfacetset_size++) {
373 14677 : SETfirst_(qh->coplanarfacetset)= nextfacet;
374 14677 : SETtruncate_(qh->coplanarfacetset, 1);
375 : }else
376 1 : 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 43589 : nextfacet= neighbor;
380 : } /* end of EACHneighbor */
381 58266 : facet= nextfacet;
382 58266 : if (facet)
383 28911 : nextfacet= NULL;
384 29355 : else if (!coplanarfacetset_size)
385 14677 : break;
386 14678 : else if (!--coplanarfacetset_size) {
387 14677 : facet= SETfirstt_(qh->coplanarfacetset, facetT);
388 14677 : SETtruncate_(qh->coplanarfacetset, 0);
389 : }else
390 1 : facet= (facetT *)qh_setdellast(qh->coplanarfacetset);
391 : } /* while True, i.e., "for each facet in qh.coplanarfacetset" */
392 14677 : if (!ischeckmax) {
393 4 : zadd_(Zfindhorizontot, *numpart - numpartinit);
394 4 : zmax_(Zfindhorizonmax, *numpart - numpartinit);
395 4 : if (newbest)
396 0 : zinc_(Znewbesthorizon);
397 : }
398 14677 : 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 14677 : 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 : */
446 184783 : facetT *qh_findbestnew(qhT *qh, pointT *point, facetT *startfacet,
447 : realT *dist, boolT bestoutside, boolT *isoutside, int *numpart) {
448 184783 : realT bestdist= -REALmax/2;
449 184783 : facetT *bestfacet= NULL, *facet;
450 184783 : int oldtrace= qh->IStracing, i;
451 184783 : unsigned int visitid= ++qh->visit_id;
452 184783 : 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 184783 : if (!startfacet || !startfacet->next) {
457 0 : if (qh->MERGING) {
458 0 : 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 0 : qh_errexit(qh, qh_ERRtopology, NULL, NULL);
460 : }else {
461 0 : qh_fprintf(qh, qh->ferr, 6002, "qhull internal error (qh_findbestnew): no new facets for point p%d\n",
462 : qh->furthest_id);
463 0 : qh_errexit(qh, qh_ERRqhull, NULL, NULL);
464 : }
465 : }
466 184783 : zinc_(Zfindnew);
467 184783 : if (qh->BESToutside || bestoutside)
468 1 : isdistoutside= False;
469 : else {
470 184782 : isdistoutside= True;
471 184782 : distoutside= qh_DISToutside; /* multiple of qh.MINoutside & qh.max_outside, see user_r.h */
472 : }
473 184783 : if (isoutside)
474 184783 : *isoutside= True;
475 184783 : *numpart= 0;
476 : #ifndef qh_NOtrace
477 184783 : if (qh->IStracing >= 4 || (qh->TRACElevel && qh->TRACEpoint >= 0 && qh->TRACEpoint == qh_pointid(qh, point))) {
478 0 : if (qh->TRACElevel > qh->IStracing)
479 0 : qh->IStracing= qh->TRACElevel;
480 0 : 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 0 : 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 0 : 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 214807 : for (i=0, facet=startfacet; i < 2; i++, facet= qh->newfacet_list) {
488 413958 : FORALLfacet_(facet) {
489 383935 : if (facet == startfacet && i)
490 1 : break;
491 383934 : facet->visitid= visitid;
492 383934 : if (!facet->flipped) {
493 383934 : qh_distplane(qh, point, facet, dist);
494 383934 : (*numpart)++;
495 383934 : if (*dist > bestdist) {
496 321428 : if (!facet->upperdelaunay || *dist >= qh->MINoutside) {
497 316692 : bestfacet= facet;
498 316692 : if (isdistoutside && *dist >= distoutside)
499 184782 : goto LABELreturn_bestnew;
500 131910 : 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 1 : bestfacet= qh_findbesthorizon(qh, !qh_IScheckmax, point, bestfacet ? bestfacet : startfacet,
508 : !qh_NOupper, &bestdist, numpart);
509 1 : *dist= bestdist;
510 1 : if (isoutside && *dist < qh->MINoutside)
511 1 : *isoutside= False;
512 0 : LABELreturn_bestnew:
513 184783 : zadd_(Zfindnewtot, *numpart);
514 184783 : zmax_(Zfindnewmax, *numpart);
515 184783 : 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 184783 : qh->IStracing= oldtrace;
518 184783 : 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 : */
560 0 : 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 0 : int zerocol= -1;
567 :
568 0 : normalp= normal + numcol - 1;
569 0 : *normalp--= (sign ? -1.0 : 1.0);
570 0 : for (i=numrow; i--; ) {
571 0 : *normalp= 0.0;
572 0 : ai= rows[i] + i + 1;
573 0 : ak= normalp+1;
574 0 : for (j=i+1; j < numcol; j++)
575 0 : *normalp -= *ai++ * *ak++;
576 0 : diagonal= (rows[i])[i];
577 0 : if (fabs_(diagonal) > qh->MINdenom_2)
578 0 : *(normalp--) /= diagonal;
579 : else {
580 0 : waszero= False;
581 0 : *normalp= qh_divzero(*normalp, diagonal, qh->MINdenom_1_2, &waszero);
582 0 : if (waszero) {
583 0 : zerocol= i;
584 0 : *(normalp--)= (sign ? -1.0 : 1.0);
585 0 : for (normal_tail= normalp+2; normal_tail < normal + numcol; normal_tail++)
586 0 : *normal_tail= 0.0;
587 : }else
588 0 : normalp--;
589 : }
590 : }
591 0 : if (zerocol != -1) {
592 0 : *nearzero= True;
593 0 : trace4((qh, qh->ferr, 4005, "qh_backnormal: zero diagonal at column %d.\n", i));
594 0 : zzinc_(Zback0);
595 0 : qh_joggle_restart(qh, "zero diagonal on back substitution");
596 : }
597 0 : } /* 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 : */
620 0 : void qh_gausselim(qhT *qh, realT **rows, int numrow, int numcol, boolT *sign, boolT *nearzero) {
621 : realT *ai, *ak, *rowp, *pivotrow;
622 0 : realT n, pivot, pivot_abs= 0.0, temp;
623 0 : int i, j, k, pivoti, flip=0;
624 :
625 0 : *nearzero= False;
626 0 : for (k=0; k < numrow; k++) {
627 0 : pivot_abs= fabs_((rows[k])[k]);
628 0 : pivoti= k;
629 0 : for (i=k+1; i < numrow; i++) {
630 0 : if ((temp= fabs_((rows[i])[k])) > pivot_abs) {
631 0 : pivot_abs= temp;
632 0 : pivoti= i;
633 : }
634 : }
635 0 : if (pivoti != k) {
636 0 : rowp= rows[pivoti];
637 0 : rows[pivoti]= rows[k];
638 0 : rows[k]= rowp;
639 0 : *sign ^= 1;
640 0 : flip ^= 1;
641 : }
642 0 : if (pivot_abs <= qh->NEARzero[k]) {
643 0 : *nearzero= True;
644 0 : if (pivot_abs == 0.0) { /* remainder of column == 0 */
645 : #ifndef qh_NOtrace
646 0 : if (qh->IStracing >= 4) {
647 0 : qh_fprintf(qh, qh->ferr, 8011, "qh_gausselim: 0 pivot at column %d. (%2.2g < %2.2g)\n", k, pivot_abs, qh->DISTround);
648 0 : qh_printmatrix(qh, qh->ferr, "Matrix:", rows, numrow, numcol);
649 : }
650 : #endif
651 0 : zzinc_(Zgauss0);
652 0 : qh_joggle_restart(qh, "zero pivot for Gaussian elimination");
653 0 : goto LABELnextcol;
654 : }
655 : }
656 0 : pivotrow= rows[k] + k;
657 0 : pivot= *pivotrow++; /* signed value of pivot, and remainder of row */
658 0 : for (i=k+1; i < numrow; i++) {
659 0 : ai= rows[i] + k;
660 0 : ak= pivotrow;
661 0 : n= (*ai++)/pivot; /* divzero() not needed since |pivot| >= |*ai| */
662 0 : for (j= numcol - (k+1); j--; )
663 0 : *ai++ -= n * *ak++;
664 : }
665 0 : LABELnextcol:
666 : ;
667 : }
668 0 : wmin_(Wmindenom, pivot_abs); /* last pivot element */
669 0 : if (qh->IStracing >= 5)
670 0 : qh_printmatrix(qh, qh->ferr, "qh_gausselem: result", rows, numrow, numcol);
671 0 : } /* 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 : */
685 36 : realT qh_getangle(qhT *qh, pointT *vect1, pointT *vect2) {
686 36 : realT angle= 0, randr;
687 : int k;
688 :
689 144 : for (k=qh->hull_dim; k--; )
690 108 : angle += *vect1++ * *vect2++;
691 36 : if (qh->RANDOMdist) {
692 0 : randr= qh_RANDOMint;
693 0 : angle += (2.0 * randr / qh_RANDOMmax - 1.0) *
694 0 : qh->RANDOMfactor;
695 : }
696 36 : trace4((qh, qh->ferr, 4006, "qh_getangle: %4.4g\n", angle));
697 36 : 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 : */
710 65920 : pointT *qh_getcenter(qhT *qh, setT *vertices) {
711 : int k;
712 : pointT *center, *coord;
713 : vertexT *vertex, **vertexp;
714 65920 : int count= qh_setsize(qh, vertices);
715 :
716 65920 : if (count < 2) {
717 0 : qh_fprintf(qh, qh->ferr, 6003, "qhull internal error (qh_getcenter): not defined for %d points\n", count);
718 0 : qh_errexit(qh, qh_ERRqhull, NULL, NULL);
719 : }
720 65920 : center= (pointT *)qh_memalloc(qh, qh->normal_size);
721 263680 : for (k=0; k < qh->hull_dim; k++) {
722 197760 : coord= center+k;
723 197760 : *coord= 0.0;
724 851778 : FOREACHvertex_(vertices)
725 654018 : *coord += vertex->point[k];
726 197760 : *coord /= count; /* count>=2 by QH6003 */
727 : }
728 65920 : 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 : */
741 65916 : pointT *qh_getcentrum(qhT *qh, facetT *facet) {
742 : realT dist;
743 : pointT *centrum, *point;
744 :
745 65916 : point= qh_getcenter(qh, facet->vertices);
746 65916 : zzinc_(Zcentrumtests);
747 65916 : qh_distplane(qh, point, facet, &dist);
748 65916 : centrum= qh_projectpoint(qh, point, facet, dist);
749 65916 : qh_memfree(qh, point, qh->normal_size);
750 65916 : 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 65916 : 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 : */
769 6 : coordT qh_getdistance(qhT *qh, facetT *facet, facetT *neighbor, coordT *mindist, coordT *maxdist) {
770 : vertexT *vertex, **vertexp;
771 : coordT dist, maxd, mind;
772 :
773 24 : FOREACHvertex_(facet->vertices)
774 18 : vertex->seen= False;
775 24 : FOREACHvertex_(neighbor->vertices)
776 18 : vertex->seen= True;
777 6 : mind= 0.0;
778 6 : maxd= 0.0;
779 24 : FOREACHvertex_(facet->vertices) {
780 18 : if (!vertex->seen) {
781 6 : zzinc_(Zbestdist);
782 6 : qh_distplane(qh, vertex->point, neighbor, &dist);
783 6 : if (dist < mind)
784 6 : mind= dist;
785 0 : else if (dist > maxd)
786 0 : maxd= dist;
787 : }
788 : }
789 6 : *mindist= mind;
790 6 : *maxdist= maxd;
791 6 : mind= -mind;
792 6 : if (maxd > mind)
793 0 : return maxd;
794 : else
795 6 : 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 : */
809 0 : void qh_normalize(qhT *qh, coordT *normal, int dim, boolT toporient) {
810 0 : qh_normalize2(qh, normal, dim, toporient, NULL, NULL);
811 0 : } /* 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 : */
845 53373 : void qh_normalize2(qhT *qh, coordT *normal, int dim, boolT toporient,
846 : realT *minnorm, boolT *ismin) {
847 : int k;
848 53373 : realT *colp, *maxp, norm= 0, temp, *norm1, *norm2, *norm3;
849 : boolT zerodiv;
850 :
851 53373 : norm1= normal+1;
852 53373 : norm2= normal+2;
853 53373 : norm3= normal+3;
854 53373 : if (dim == 2)
855 0 : norm= sqrt((*normal)*(*normal) + (*norm1)*(*norm1));
856 53373 : else if (dim == 3)
857 53373 : norm= sqrt((*normal)*(*normal) + (*norm1)*(*norm1) + (*norm2)*(*norm2));
858 0 : else if (dim == 4) {
859 0 : norm= sqrt((*normal)*(*normal) + (*norm1)*(*norm1) + (*norm2)*(*norm2)
860 0 : + (*norm3)*(*norm3));
861 0 : }else if (dim > 4) {
862 0 : norm= (*normal)*(*normal) + (*norm1)*(*norm1) + (*norm2)*(*norm2)
863 0 : + (*norm3)*(*norm3);
864 0 : for (k=dim-4, colp=normal+4; k--; colp++)
865 0 : norm += (*colp) * (*colp);
866 0 : norm= sqrt(norm);
867 : }
868 53373 : if (minnorm) {
869 0 : if (norm < *minnorm)
870 0 : *ismin= True;
871 : else
872 0 : *ismin= False;
873 : }
874 53373 : wmin_(Wmindenom, norm);
875 53373 : if (norm > qh->MINdenom) {
876 53373 : if (!toporient)
877 26649 : norm= -norm;
878 53373 : *normal /= norm;
879 53373 : *norm1 /= norm;
880 53373 : if (dim == 2)
881 : ; /* all done */
882 53373 : else if (dim == 3)
883 53373 : *norm2 /= norm;
884 0 : else if (dim == 4) {
885 0 : *norm2 /= norm;
886 0 : *norm3 /= norm;
887 0 : }else if (dim >4) {
888 0 : *norm2 /= norm;
889 0 : *norm3 /= norm;
890 0 : for (k=dim-4, colp=normal+4; k--; )
891 0 : *colp++ /= norm;
892 : }
893 0 : }else if (norm == 0.0) {
894 0 : temp= sqrt(1.0/dim);
895 0 : for (k=dim, colp=normal; k--; )
896 0 : *colp++= temp;
897 : }else {
898 0 : if (!toporient)
899 0 : norm= -norm;
900 0 : for (k=dim, colp=normal; k--; colp++) { /* k used below */
901 0 : temp= qh_divzero(*colp, norm, qh->MINdenom_1, &zerodiv);
902 0 : if (!zerodiv)
903 0 : *colp= temp;
904 : else {
905 0 : maxp= qh_maxabsval(normal, dim);
906 0 : temp= ((*maxp * norm >= 0.0) ? 1.0 : -1.0);
907 0 : for (k=dim, colp=normal; k--; colp++)
908 0 : *colp= 0.0;
909 0 : *maxp= temp;
910 0 : zzinc_(Znearlysingular);
911 : /* qh_joggle_restart ignored for Znearlysingular, normal part of qh_sethyperplane_gauss */
912 0 : trace0((qh, qh->ferr, 1, "qh_normalize: norm=%2.2g too small during p%d\n",
913 : norm, qh->furthest_id));
914 0 : 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 : */
935 65916 : pointT *qh_projectpoint(qhT *qh, pointT *point, facetT *facet, realT dist) {
936 : pointT *newpoint, *np, *normal;
937 65916 : int normsize= qh->normal_size;
938 : int k;
939 : void **freelistp; /* used if !qh_NOmem by qh_memalloc_() */
940 :
941 65916 : qh_memalloc_(qh, normsize, freelistp, newpoint, pointT);
942 65916 : np= newpoint;
943 65916 : normal= facet->normal;
944 263664 : for (k=qh->hull_dim; k--; )
945 197748 : *(np++)= *point++ - dist * *normal++;
946 65916 : 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 : */
971 53373 : void qh_setfacetplane(qhT *qh, facetT *facet) {
972 : pointT *point;
973 : vertexT *vertex, **vertexp;
974 53373 : int normsize= qh->normal_size;
975 53373 : int k,i, oldtrace= 0;
976 : realT dist;
977 : void **freelistp; /* used if !qh_NOmem by qh_memalloc_() */
978 : coordT *coord, *gmcoord;
979 53373 : pointT *point0= SETfirstt_(facet->vertices, vertexT)->point;
980 53373 : boolT nearzero= False;
981 :
982 53373 : zzinc_(Zsetplane);
983 53373 : if (!facet->normal)
984 53373 : qh_memalloc_(qh, normsize, freelistp, facet->normal, coordT);
985 : #ifndef qh_NOtrace
986 53373 : if (facet == qh->tracefacet) {
987 0 : oldtrace= qh->IStracing;
988 0 : qh->IStracing= 5;
989 0 : qh_fprintf(qh, qh->ferr, 8012, "qh_setfacetplane: facet f%d created.\n", facet->id);
990 0 : qh_fprintf(qh, qh->ferr, 8013, " Last point added to hull was p%d.", qh->furthest_id);
991 0 : if (zzval_(Ztotmerge))
992 0 : qh_fprintf(qh, qh->ferr, 8014, " Last merge was #%d.", zzval_(Ztotmerge));
993 0 : qh_fprintf(qh, qh->ferr, 8015, "\n\nCurrent summary is:\n");
994 0 : qh_printsummary(qh, qh->ferr);
995 : }
996 : #endif
997 53373 : if (qh->hull_dim <= 4) {
998 53373 : i= 0;
999 53373 : if (qh->RANDOMdist) {
1000 0 : gmcoord= qh->gm_matrix;
1001 0 : FOREACHvertex_(facet->vertices) {
1002 0 : qh->gm_row[i++]= gmcoord;
1003 0 : coord= vertex->point;
1004 0 : for (k=qh->hull_dim; k--; )
1005 0 : *(gmcoord++)= *coord++ * qh_randomfactor(qh, qh->RANDOMa, qh->RANDOMb);
1006 : }
1007 : }else {
1008 213492 : FOREACHvertex_(facet->vertices)
1009 160119 : qh->gm_row[i++]= vertex->point;
1010 : }
1011 53373 : qh_sethyperplane_det(qh, qh->hull_dim, qh->gm_row, point0, facet->toporient,
1012 : facet->normal, &facet->offset, &nearzero);
1013 : }
1014 53373 : if (qh->hull_dim > 4 || nearzero) {
1015 0 : i= 0;
1016 0 : gmcoord= qh->gm_matrix;
1017 0 : FOREACHvertex_(facet->vertices) {
1018 0 : if (vertex->point != point0) {
1019 0 : qh->gm_row[i++]= gmcoord;
1020 0 : coord= vertex->point;
1021 0 : point= point0;
1022 0 : for (k=qh->hull_dim; k--; )
1023 0 : *(gmcoord++)= *coord++ - *point++;
1024 : }
1025 : }
1026 0 : qh->gm_row[i]= gmcoord; /* for areasimplex */
1027 0 : if (qh->RANDOMdist) {
1028 0 : gmcoord= qh->gm_matrix;
1029 0 : for (i=qh->hull_dim-1; i--; ) {
1030 0 : for (k=qh->hull_dim; k--; )
1031 0 : *(gmcoord++) *= qh_randomfactor(qh, qh->RANDOMa, qh->RANDOMb);
1032 : }
1033 : }
1034 0 : qh_sethyperplane_gauss(qh, qh->hull_dim, qh->gm_row, point0, facet->toporient,
1035 : facet->normal, &facet->offset, &nearzero);
1036 0 : if (nearzero) {
1037 0 : if (qh_orientoutside(qh, facet)) {
1038 0 : 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 53373 : facet->upperdelaunay= False;
1053 53373 : if (qh->DELAUNAY) {
1054 53373 : if (qh->UPPERdelaunay) { /* matches qh_triangulate_facet and qh.lower_threshold in qh_initbuild */
1055 0 : if (facet->normal[qh->hull_dim -1] >= qh->ANGLEround * qh_ZEROdelaunay)
1056 0 : facet->upperdelaunay= True;
1057 : }else {
1058 53373 : if (facet->normal[qh->hull_dim -1] > -qh->ANGLEround * qh_ZEROdelaunay)
1059 52 : facet->upperdelaunay= True;
1060 : }
1061 : }
1062 53373 : if (qh->PRINTstatistics || qh->IStracing || qh->TRACElevel || qh->JOGGLEmax < REALmax) {
1063 0 : qh->old_randomdist= qh->RANDOMdist;
1064 0 : qh->RANDOMdist= False;
1065 0 : FOREACHvertex_(facet->vertices) {
1066 0 : if (vertex->point != point0) {
1067 0 : boolT istrace= False;
1068 0 : zinc_(Zdiststat);
1069 0 : qh_distplane(qh, vertex->point, facet, &dist);
1070 0 : dist= fabs_(dist);
1071 0 : zinc_(Znewvertex);
1072 0 : wadd_(Wnewvertex, dist);
1073 0 : if (dist > wwval_(Wnewvertexmax)) {
1074 0 : wwval_(Wnewvertexmax)= dist;
1075 0 : if (dist > qh->max_outside) {
1076 0 : qh->max_outside= dist; /* used by qh_maxouter(qh) */
1077 0 : if (dist > qh->TRACEdist)
1078 0 : istrace= True;
1079 : }
1080 0 : }else if (-dist > qh->TRACEdist)
1081 0 : istrace= True;
1082 0 : if (istrace) {
1083 0 : 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 0 : qh_errprint(qh, "DISTANT", facet, NULL, NULL, NULL);
1086 : }
1087 : }
1088 : }
1089 0 : qh->RANDOMdist= qh->old_randomdist;
1090 : }
1091 : #ifndef qh_NOtrace
1092 53373 : if (qh->IStracing >= 4) {
1093 0 : qh_fprintf(qh, qh->ferr, 8017, "qh_setfacetplane: f%d offset %2.2g normal: ",
1094 : facet->id, facet->offset);
1095 0 : for (k=0; k < qh->hull_dim; k++)
1096 0 : qh_fprintf(qh, qh->ferr, 8018, "%2.2g ", facet->normal[k]);
1097 0 : qh_fprintf(qh, qh->ferr, 8019, "\n");
1098 : }
1099 : #endif
1100 53373 : qh_checkflipped(qh, facet, NULL, qh_ALL);
1101 53373 : if (facet == qh->tracefacet) {
1102 0 : qh->IStracing= oldtrace;
1103 0 : qh_printfacet(qh, qh->ferr, facet);
1104 : }
1105 53373 : } /* 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 : */
1152 53373 : 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 53373 : if (dim == 2) {
1160 0 : normal[0]= dY(1,0);
1161 0 : normal[1]= dX(0,1);
1162 0 : qh_normalize2(qh, normal, dim, toporient, NULL, NULL);
1163 0 : *offset= -(point0[0]*normal[0]+point0[1]*normal[1]);
1164 0 : *nearzero= False; /* since nearzero norm => incident points */
1165 53373 : }else if (dim == 3) {
1166 53373 : normal[0]= det2_(dY(2,0), dZ(2,0),
1167 : dY(1,0), dZ(1,0));
1168 53373 : normal[1]= det2_(dX(1,0), dZ(1,0),
1169 : dX(2,0), dZ(2,0));
1170 53373 : normal[2]= det2_(dX(2,0), dY(2,0),
1171 : dX(1,0), dY(1,0));
1172 53373 : qh_normalize2(qh, normal, dim, toporient, NULL, NULL);
1173 53373 : *offset= -(point0[0]*normal[0] + point0[1]*normal[1]
1174 53373 : + point0[2]*normal[2]);
1175 53373 : maxround= qh->DISTround;
1176 213492 : for (i=dim; i--; ) {
1177 160119 : point= rows[i];
1178 160119 : if (point != point0) {
1179 106746 : dist= *offset + (point[0]*normal[0] + point[1]*normal[1]
1180 106746 : + point[2]*normal[2]);
1181 106746 : if (dist > maxround || dist < -maxround) {
1182 0 : *nearzero= True;
1183 0 : break;
1184 : }
1185 : }
1186 : }
1187 0 : }else if (dim == 4) {
1188 0 : 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 0 : 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 0 : 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 0 : 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 0 : qh_normalize2(qh, normal, dim, toporient, NULL, NULL);
1201 0 : *offset= -(point0[0]*normal[0] + point0[1]*normal[1]
1202 0 : + point0[2]*normal[2] + point0[3]*normal[3]);
1203 0 : maxround= qh->DISTround;
1204 0 : for (i=dim; i--; ) {
1205 0 : point= rows[i];
1206 0 : if (point != point0) {
1207 0 : dist= *offset + (point[0]*normal[0] + point[1]*normal[1]
1208 0 : + point[2]*normal[2] + point[3]*normal[3]);
1209 0 : if (dist > maxround || dist < -maxround) {
1210 0 : *nearzero= True;
1211 0 : break;
1212 : }
1213 : }
1214 : }
1215 : }
1216 53373 : if (*nearzero) {
1217 0 : zzinc_(Zminnorm);
1218 : /* qh_joggle_restart not needed, will call qh_sethyperplane_gauss instead */
1219 0 : trace0((qh, qh->ferr, 3, "qh_sethyperplane_det: degenerate norm during p%d, use qh_sethyperplane_gauss instead.\n", qh->furthest_id));
1220 : }
1221 53373 : } /* 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 : */
1250 0 : 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 0 : boolT sign= toporient, nearzero2= False;
1255 :
1256 0 : qh_gausselim(qh, rows, dim-1, dim, &sign, nearzero);
1257 0 : for (k=dim-1; k--; ) {
1258 0 : if ((rows[k])[k] < 0)
1259 0 : sign ^= 1;
1260 : }
1261 0 : if (*nearzero) {
1262 0 : zzinc_(Znearlysingular);
1263 : /* qh_joggle_restart ignored for Znearlysingular, normal part of qh_sethyperplane_gauss */
1264 0 : trace0((qh, qh->ferr, 4, "qh_sethyperplane_gauss: nearly singular or axis parallel hyperplane during p%d.\n", qh->furthest_id));
1265 0 : qh_backnormal(qh, rows, dim-1, dim, sign, normal, &nearzero2);
1266 : }else {
1267 0 : qh_backnormal(qh, rows, dim-1, dim, sign, normal, &nearzero2);
1268 0 : if (nearzero2) {
1269 0 : zzinc_(Znearlysingular);
1270 0 : 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 0 : if (nearzero2)
1274 0 : *nearzero= True;
1275 0 : qh_normalize2(qh, normal, dim, True, NULL, NULL);
1276 0 : pointcoord= point0;
1277 0 : normalcoef= normal;
1278 0 : *offset= -(*pointcoord++ * *normalcoef++);
1279 0 : for (k=dim-1; k--; )
1280 0 : *offset -= *pointcoord++ * *normalcoef++;
1281 0 : } /* sethyperplane_gauss */
1282 :
1283 :
1284 :
|
