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