Free cookie consent management tool by TermsFeed Policy Generator

source: branches/HeuristicLab.Analysis.AlgorithmBehavior/qhull-2012.1/src/libqhull/geom.c @ 11303

Last change on this file since 11303 was 10207, checked in by ascheibe, 11 years ago

#1886 added a unit test for volume calculation and the qhull library
File size: 42.6 KB
Line 
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-2012 The Geometry Center.
10   $Id: //main/2011/qhull/src/libqhull/geom.c#3 $$Change: 1464 $
11   $DateTime: 2012/01/25 22:58:41 $$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*/
36void 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*/
140facetT *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;
234LABELreturn_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*/
281facetT *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*/
412facetT *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)
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;
474LABELreturn_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 4.4-9 for back substitution
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*/
521void 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*/
581void 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*/
644realT 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*/
669pointT *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;
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*/
700pointT *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*/
728realT 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*/
768void 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*/
804void 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*/
893pointT *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 !qh_NOmem */
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*/
929void 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 !qh_NOmem */
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*/
1103void 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*/
1201void 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
Note: See TracBrowser for help on using the repository browser.