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