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

hpp-fcl
Author(s):
autogenerated on Sat Nov 23 2024 03:44:58