1 /*<html><pre> -<a href="qh-geom.htm"
2 >-------------------------------</a><a name="TOP">-</a>
3
4
5 geom2.c
6 infrequently used geometric routines of qhull
7
8 see qh-geom.htm and geom.h
9
10 copyright (c) 1993-2003 The Geometry Center
11
12 frequently used code goes into geom.c
13 */
14
15 #include "qhull_a.h"
16
17 /*================== functions in alphabetic order ============*/
18
19 /*-<a href="qh-geom.htm#TOC"
20 >-------------------------------</a><a name="copypoints">-</a>
21
22 qh_copypoints( points, numpoints, dimension)
23 return malloc'd copy of points
24 */
qh_copypoints(coordT * points,int numpoints,int dimension)25 coordT *qh_copypoints (coordT *points, int numpoints, int dimension) {
26 int size;
27 coordT *newpoints;
28
29 size= numpoints * dimension * sizeof(coordT);
30 if (!(newpoints=(coordT*)malloc(size))) {
31 fprintf(qh ferr, "qhull error: insufficient memory to copy %d points\n",
32 numpoints);
33 qh_errexit(qh_ERRmem, NULL, NULL);
34 }
35 memcpy ((char *)newpoints, (char *)points, size);
36 return newpoints;
37 } /* copypoints */
38
39 /*-<a href="qh-geom.htm#TOC"
40 >-------------------------------</a><a name="crossproduct">-</a>
41
42 qh_crossproduct( dim, vecA, vecB, vecC )
43 crossproduct of 2 dim vectors
44 C= A x B
45
46 notes:
47 from Glasner, Graphics Gems I, p. 639
48 only defined for dim==3
49 */
qh_crossproduct(int dim,realT vecA[3],realT vecB[3],realT vecC[3])50 void qh_crossproduct (int dim, realT vecA[3], realT vecB[3], realT vecC[3]){
51
52 if (dim == 3) {
53 vecC[0]= det2_(vecA[1], vecA[2],
54 vecB[1], vecB[2]);
55 vecC[1]= - det2_(vecA[0], vecA[2],
56 vecB[0], vecB[2]);
57 vecC[2]= det2_(vecA[0], vecA[1],
58 vecB[0], vecB[1]);
59 }
60 } /* vcross */
61
62 /*-<a href="qh-geom.htm#TOC"
63 >-------------------------------</a><a name="determinant">-</a>
64
65 qh_determinant( rows, dim, nearzero )
66 compute signed determinant of a square matrix
67 uses qh.NEARzero to test for degenerate matrices
68
69 returns:
70 determinant
71 overwrites rows and the matrix
72 if dim == 2 or 3
73 nearzero iff determinant < qh NEARzero[dim-1]
74 (not quite correct, not critical)
75 if dim >= 4
76 nearzero iff diagonal[k] < qh NEARzero[k]
77 */
qh_determinant(realT ** rows,int dim,boolT * nearzero)78 realT qh_determinant (realT **rows, int dim, boolT *nearzero) {
79 realT det=0;
80 int i;
81 boolT sign= False;
82
83 *nearzero= False;
84 if (dim < 2) {
85 fprintf (qh ferr, "qhull internal error (qh_determinate): only implemented for dimension >= 2\n");
86 qh_errexit (qh_ERRqhull, NULL, NULL);
87 }else if (dim == 2) {
88 det= det2_(rows[0][0], rows[0][1],
89 rows[1][0], rows[1][1]);
90 if (fabs_(det) < qh NEARzero[1]) /* not really correct, what should this be? */
91 *nearzero= True;
92 }else if (dim == 3) {
93 det= det3_(rows[0][0], rows[0][1], rows[0][2],
94 rows[1][0], rows[1][1], rows[1][2],
95 rows[2][0], rows[2][1], rows[2][2]);
96 if (fabs_(det) < qh NEARzero[2]) /* not really correct, what should this be? */
97 *nearzero= True;
98 }else {
99 qh_gausselim(rows, dim, dim, &sign, nearzero); /* if nearzero, diagonal still ok*/
100 det= 1.0;
101 for (i= dim; i--; )
102 det *= (rows[i])[i];
103 if (sign)
104 det= -det;
105 }
106 return det;
107 } /* determinant */
108
109 /*-<a href="qh-geom.htm#TOC"
110 >-------------------------------</a><a name="detjoggle">-</a>
111
112 qh_detjoggle( points, numpoints, dimension )
113 determine default max joggle for point array
114 as qh_distround * qh_JOGGLEdefault
115
116 returns:
117 initial value for JOGGLEmax from points and REALepsilon
118
119 notes:
120 computes DISTround since qh_maxmin not called yet
121 if qh SCALElast, last dimension will be scaled later to MAXwidth
122
123 loop duplicated from qh_maxmin
124 */
qh_detjoggle(pointT * points,int numpoints,int dimension)125 realT qh_detjoggle (pointT *points, int numpoints, int dimension) {
126 realT abscoord, distround, joggle, maxcoord, mincoord;
127 pointT *point, *pointtemp;
128 realT maxabs= -REALmax;
129 realT sumabs= 0;
130 realT maxwidth= 0;
131 int k;
132
133 for (k= 0; k < dimension; k++) {
134 if (qh SCALElast && k == dimension-1)
135 abscoord= maxwidth;
136 else if (qh DELAUNAY && k == dimension-1) /* will qh_setdelaunay() */
137 abscoord= 2 * maxabs * maxabs; /* may be low by qh hull_dim/2 */
138 else {
139 maxcoord= -REALmax;
140 mincoord= REALmax;
141 FORALLpoint_(points, numpoints) {
142 maximize_(maxcoord, point[k]);
143 minimize_(mincoord, point[k]);
144 }
145 maximize_(maxwidth, maxcoord-mincoord);
146 abscoord= fmax_(maxcoord, -mincoord);
147 }
148 sumabs += abscoord;
149 maximize_(maxabs, abscoord);
150 } /* for k */
151 distround= qh_distround (qh hull_dim, maxabs, sumabs);
152 joggle= distround * qh_JOGGLEdefault;
153 maximize_(joggle, REALepsilon * qh_JOGGLEdefault);
154 trace2((qh ferr, "qh_detjoggle: joggle=%2.2g maxwidth=%2.2g\n", joggle, maxwidth));
155 return joggle;
156 } /* detjoggle */
157
158 /*-<a href="qh-geom.htm#TOC"
159 >-------------------------------</a><a name="detroundoff">-</a>
160
161 qh_detroundoff()
162 determine maximum roundoff errors from
163 REALepsilon, REALmax, REALmin, qh.hull_dim, qh.MAXabs_coord,
164 qh.MAXsumcoord, qh.MAXwidth, qh.MINdenom_1
165
166 accounts for qh.SETroundoff, qh.RANDOMdist, qh MERGEexact
167 qh.premerge_cos, qh.postmerge_cos, qh.premerge_centrum,
168 qh.postmerge_centrum, qh.MINoutside,
169 qh_RATIOnearinside, qh_COPLANARratio, qh_WIDEcoplanar
170
171 returns:
172 sets qh.DISTround, etc. (see below)
173 appends precision constants to qh.qhull_options
174
175 see:
176 qh_maxmin() for qh.NEARzero
177
178 design:
179 determine qh.DISTround for distance computations
180 determine minimum denominators for qh_divzero
181 determine qh.ANGLEround for angle computations
182 adjust qh.premerge_cos,... for roundoff error
183 determine qh.ONEmerge for maximum error due to a single merge
184 determine qh.NEARinside, qh.MAXcoplanar, qh.MINvisible,
185 qh.MINoutside, qh.WIDEfacet
186 initialize qh.max_vertex and qh.minvertex
187 */
qh_detroundoff(void)188 void qh_detroundoff (void) {
189
190 qh_option ("_max-width", NULL, &qh MAXwidth);
191 if (!qh SETroundoff) {
192 qh DISTround= qh_distround (qh hull_dim, qh MAXabs_coord, qh MAXsumcoord);
193 if (qh RANDOMdist)
194 qh DISTround += qh RANDOMfactor * qh MAXabs_coord;
195 qh_option ("Error-roundoff", NULL, &qh DISTround);
196 }
197 qh MINdenom= qh MINdenom_1 * qh MAXabs_coord;
198 qh MINdenom_1_2= sqrt (qh MINdenom_1 * qh hull_dim) ; /* if will be normalized */
199 qh MINdenom_2= qh MINdenom_1_2 * qh MAXabs_coord;
200 /* for inner product */
201 qh ANGLEround= 1.01 * qh hull_dim * REALepsilon;
202 if (qh RANDOMdist)
203 qh ANGLEround += qh RANDOMfactor;
204 if (qh premerge_cos < REALmax/2) {
205 qh premerge_cos -= qh ANGLEround;
206 if (qh RANDOMdist)
207 qh_option ("Angle-premerge-with-random", NULL, &qh premerge_cos);
208 }
209 if (qh postmerge_cos < REALmax/2) {
210 qh postmerge_cos -= qh ANGLEround;
211 if (qh RANDOMdist)
212 qh_option ("Angle-postmerge-with-random", NULL, &qh postmerge_cos);
213 }
214 qh premerge_centrum += 2 * qh DISTround; /*2 for centrum and distplane()*/
215 qh postmerge_centrum += 2 * qh DISTround;
216 if (qh RANDOMdist && (qh MERGEexact || qh PREmerge))
217 qh_option ("Centrum-premerge-with-random", NULL, &qh premerge_centrum);
218 if (qh RANDOMdist && qh POSTmerge)
219 qh_option ("Centrum-postmerge-with-random", NULL, &qh postmerge_centrum);
220 { /* compute ONEmerge, max vertex offset for merging simplicial facets */
221 realT maxangle= 1.0, maxrho;
222
223 minimize_(maxangle, qh premerge_cos);
224 minimize_(maxangle, qh postmerge_cos);
225 /* max diameter * sin theta + DISTround for vertex to its hyperplane */
226 qh ONEmerge= sqrt (qh hull_dim) * qh MAXwidth *
227 sqrt (1.0 - maxangle * maxangle) + qh DISTround;
228 maxrho= qh hull_dim * qh premerge_centrum + qh DISTround;
229 maximize_(qh ONEmerge, maxrho);
230 maxrho= qh hull_dim * qh postmerge_centrum + qh DISTround;
231 maximize_(qh ONEmerge, maxrho);
232 if (qh MERGING)
233 qh_option ("_one-merge", NULL, &qh ONEmerge);
234 }
235 qh NEARinside= qh ONEmerge * qh_RATIOnearinside; /* only used if qh KEEPnearinside */
236 if (qh JOGGLEmax < REALmax/2 && (qh KEEPcoplanar || qh KEEPinside)) {
237 realT maxdist; /* adjust qh.NEARinside for joggle */
238 qh KEEPnearinside= True;
239 maxdist= sqrt (qh hull_dim) * qh JOGGLEmax + qh DISTround;
240 maxdist= 2*maxdist; /* vertex and coplanar point can joggle in opposite directions */
241 maximize_(qh NEARinside, maxdist); /* must agree with qh_nearcoplanar() */
242 }
243 if (qh KEEPnearinside)
244 qh_option ("_near-inside", NULL, &qh NEARinside);
245 if (qh JOGGLEmax < qh DISTround) {
246 fprintf (qh ferr, "qhull error: the joggle for 'QJn', %.2g, is below roundoff for distance computations, %.2g\n",
247 qh JOGGLEmax, qh DISTround);
248 qh_errexit (qh_ERRinput, NULL, NULL);
249 }
250 if (qh MINvisible > REALmax/2) {
251 if (!qh MERGING)
252 qh MINvisible= qh DISTround;
253 else if (qh hull_dim <= 3)
254 qh MINvisible= qh premerge_centrum;
255 else
256 qh MINvisible= qh_COPLANARratio * qh premerge_centrum;
257 if (qh APPROXhull && qh MINvisible > qh MINoutside)
258 qh MINvisible= qh MINoutside;
259 qh_option ("Visible-distance", NULL, &qh MINvisible);
260 }
261 if (qh MAXcoplanar > REALmax/2) {
262 qh MAXcoplanar= qh MINvisible;
263 qh_option ("U-coplanar-distance", NULL, &qh MAXcoplanar);
264 }
265 if (!qh APPROXhull) { /* user may specify qh MINoutside */
266 qh MINoutside= 2 * qh MINvisible;
267 if (qh premerge_cos < REALmax/2)
268 maximize_(qh MINoutside, (1- qh premerge_cos) * qh MAXabs_coord);
269 qh_option ("Width-outside", NULL, &qh MINoutside);
270 }
271 qh WIDEfacet= qh MINoutside;
272 maximize_(qh WIDEfacet, qh_WIDEcoplanar * qh MAXcoplanar);
273 maximize_(qh WIDEfacet, qh_WIDEcoplanar * qh MINvisible);
274 qh_option ("_wide-facet", NULL, &qh WIDEfacet);
275 if (qh MINvisible > qh MINoutside + 3 * REALepsilon
276 && !qh BESToutside && !qh FORCEoutput)
277 fprintf (qh ferr, "qhull input warning: minimum visibility V%.2g is greater than \nminimum outside W%.2g. Flipped facets are likely.\n",
278 qh MINvisible, qh MINoutside);
279 qh max_vertex= qh DISTround;
280 qh min_vertex= -qh DISTround;
281 /* numeric constants reported in printsummary */
282 } /* detroundoff */
283
284 /*-<a href="qh-geom.htm#TOC"
285 >-------------------------------</a><a name="detsimplex">-</a>
286
287 qh_detsimplex( apex, points, dim, nearzero )
288 compute determinant of a simplex with point apex and base points
289
290 returns:
291 signed determinant and nearzero from qh_determinant
292
293 notes:
294 uses qh.gm_matrix/qh.gm_row (assumes they're big enough)
295
296 design:
297 construct qm_matrix by subtracting apex from points
298 compute determinate
299 */
qh_detsimplex(pointT * apex,setT * points,int dim,boolT * nearzero)300 realT qh_detsimplex(pointT *apex, setT *points, int dim, boolT *nearzero) {
301 pointT *coorda, *coordp, *gmcoord, *point, **pointp;
302 coordT **rows;
303 int k, i=0;
304 realT det;
305
306 zinc_(Zdetsimplex);
307 gmcoord= qh gm_matrix;
308 rows= qh gm_row;
309 FOREACHpoint_(points) {
310 if (i == dim)
311 break;
312 rows[i++]= gmcoord;
313 coordp= point;
314 coorda= apex;
315 for (k= dim; k--; )
316 *(gmcoord++)= *coordp++ - *coorda++;
317 }
318 if (i < dim) {
319 fprintf (qh ferr, "qhull internal error (qh_detsimplex): #points %d < dimension %d\n",
320 i, dim);
321 qh_errexit (qh_ERRqhull, NULL, NULL);
322 }
323 det= qh_determinant (rows, dim, nearzero);
324 trace2((qh ferr, "qh_detsimplex: det=%2.2g for point p%d, dim %d, nearzero? %d\n",
325 det, qh_pointid(apex), dim, *nearzero));
326 return det;
327 } /* detsimplex */
328
329 /*-<a href="qh-geom.htm#TOC"
330 >-------------------------------</a><a name="distnorm">-</a>
331
332 qh_distnorm( dim, point, normal, offset )
333 return distance from point to hyperplane at normal/offset
334
335 returns:
336 dist
337
338 notes:
339 dist > 0 if point is outside of hyperplane
340
341 see:
342 qh_distplane in geom.c
343 */
qh_distnorm(int dim,pointT * point,pointT * normal,realT * offsetp)344 realT qh_distnorm (int dim, pointT *point, pointT *normal, realT *offsetp) {
345 coordT *normalp= normal, *coordp= point;
346 realT dist;
347 int k;
348
349 dist= *offsetp;
350 for (k= dim; k--; )
351 dist += *(coordp++) * *(normalp++);
352 return dist;
353 } /* distnorm */
354
355 /*-<a href="qh-geom.htm#TOC"
356 >-------------------------------</a><a name="distround">-</a>
357
358 qh_distround ( dimension, maxabs, maxsumabs )
359 compute maximum round-off error for a distance computation
360 to a normalized hyperplane
361 maxabs is the maximum absolute value of a coordinate
362 maxsumabs is the maximum possible sum of absolute coordinate values
363
364 returns:
365 max dist round for REALepsilon
366
367 notes:
368 calculate roundoff error according to
369 Lemma 3.2-1 of Golub and van Loan "Matrix Computation"
370 use sqrt(dim) since one vector is normalized
371 or use maxsumabs since one vector is < 1
372 */
qh_distround(int dimension,realT maxabs,realT maxsumabs)373 realT qh_distround (int dimension, realT maxabs, realT maxsumabs) {
374 realT maxdistsum, maxround;
375
376 maxdistsum= sqrt (dimension) * maxabs;
377 minimize_( maxdistsum, maxsumabs);
378 maxround= REALepsilon * (dimension * maxdistsum * 1.01 + maxabs);
379 /* adds maxabs for offset */
380 trace4((qh ferr, "qh_distround: %2.2g maxabs %2.2g maxsumabs %2.2g maxdistsum %2.2g\n",
381 maxround, maxabs, maxsumabs, maxdistsum));
382 return maxround;
383 } /* distround */
384
385 /*-<a href="qh-geom.htm#TOC"
386 >-------------------------------</a><a name="divzero">-</a>
387
388 qh_divzero( numer, denom, mindenom1, zerodiv )
389 divide by a number that's nearly zero
390 mindenom1= minimum denominator for dividing into 1.0
391
392 returns:
393 quotient
394 sets zerodiv and returns 0.0 if it would overflow
395
396 design:
397 if numer is nearly zero and abs(numer) < abs(denom)
398 return numer/denom
399 else if numer is nearly zero
400 return 0 and zerodiv
401 else if denom/numer non-zero
402 return numer/denom
403 else
404 return 0 and zerodiv
405 */
qh_divzero(realT numer,realT denom,realT mindenom1,boolT * zerodiv)406 realT qh_divzero (realT numer, realT denom, realT mindenom1, boolT *zerodiv) {
407 realT temp, numerx, denomx;
408
409
410 if (numer < mindenom1 && numer > -mindenom1) {
411 numerx= fabs_(numer);
412 denomx= fabs_(denom);
413 if (numerx < denomx) {
414 *zerodiv= False;
415 return numer/denom;
416 }else {
417 *zerodiv= True;
418 return 0.0;
419 }
420 }
421 temp= denom/numer;
422 if (temp > mindenom1 || temp < -mindenom1) {
423 *zerodiv= False;
424 return numer/denom;
425 }else {
426 *zerodiv= True;
427 return 0.0;
428 }
429 } /* divzero */
430
431
432 /*-<a href="qh-geom.htm#TOC"
433 >-------------------------------</a><a name="facetarea">-</a>
434
435 qh_facetarea( facet )
436 return area for a facet
437
438 notes:
439 if non-simplicial,
440 uses centrum to triangulate facet and sums the projected areas.
441 if (qh DELAUNAY),
442 computes projected area instead for last coordinate
443 assumes facet->normal exists
444 projecting tricoplanar facets to the hyperplane does not appear to make a difference
445
446 design:
447 if simplicial
448 compute area
449 else
450 for each ridge
451 compute area from centrum to ridge
452 negate area if upper Delaunay facet
453 */
qh_facetarea(facetT * facet)454 realT qh_facetarea (facetT *facet) {
455 vertexT *apex;
456 pointT *centrum;
457 realT area= 0.0;
458 ridgeT *ridge, **ridgep;
459
460 if (facet->simplicial) {
461 apex= SETfirstt_(facet->vertices, vertexT);
462 area= qh_facetarea_simplex (qh hull_dim, apex->point, facet->vertices,
463 apex, facet->toporient, facet->normal, &facet->offset);
464 }else {
465 if (qh CENTERtype == qh_AScentrum)
466 centrum= facet->center;
467 else
468 centrum= qh_getcentrum (facet);
469 FOREACHridge_(facet->ridges)
470 area += qh_facetarea_simplex (qh hull_dim, centrum, ridge->vertices,
471 NULL, (ridge->top == facet), facet->normal, &facet->offset);
472 if (qh CENTERtype != qh_AScentrum)
473 qh_memfree (centrum, qh normal_size);
474 }
475 if (facet->upperdelaunay && qh DELAUNAY)
476 area= -area; /* the normal should be [0,...,1] */
477 trace4((qh ferr, "qh_facetarea: f%d area %2.2g\n", facet->id, area));
478 return area;
479 } /* facetarea */
480
481 /*-<a href="qh-geom.htm#TOC"
482 >-------------------------------</a><a name="facetarea_simplex">-</a>
483
484 qh_facetarea_simplex( dim, apex, vertices, notvertex, toporient, normal, offset )
485 return area for a simplex defined by
486 an apex, a base of vertices, an orientation, and a unit normal
487 if simplicial or tricoplanar facet,
488 notvertex is defined and it is skipped in vertices
489
490 returns:
491 computes area of simplex projected to plane [normal,offset]
492 returns 0 if vertex too far below plane (qh WIDEfacet)
493 vertex can't be apex of tricoplanar facet
494
495 notes:
496 if (qh DELAUNAY),
497 computes projected area instead for last coordinate
498 uses qh gm_matrix/gm_row and qh hull_dim
499 helper function for qh_facetarea
500
501 design:
502 if Notvertex
503 translate simplex to apex
504 else
505 project simplex to normal/offset
506 translate simplex to apex
507 if Delaunay
508 set last row/column to 0 with -1 on diagonal
509 else
510 set last row to Normal
511 compute determinate
512 scale and flip sign for area
513 */
qh_facetarea_simplex(int dim,coordT * apex,setT * vertices,vertexT * notvertex,boolT toporient,coordT * normal,realT * offset)514 realT qh_facetarea_simplex (int dim, coordT *apex, setT *vertices,
515 vertexT *notvertex, boolT toporient, coordT *normal, realT *offset) {
516 pointT *coorda, *coordp, *gmcoord;
517 coordT **rows, *normalp;
518 int k, i=0;
519 realT area, dist;
520 vertexT *vertex, **vertexp;
521 boolT nearzero;
522
523 gmcoord= qh gm_matrix;
524 rows= qh gm_row;
525 FOREACHvertex_(vertices) {
526 if (vertex == notvertex)
527 continue;
528 rows[i++]= gmcoord;
529 coorda= apex;
530 coordp= vertex->point;
531 normalp= normal;
532 if (notvertex) {
533 for (k= dim; k--; )
534 *(gmcoord++)= *coordp++ - *coorda++;
535 }else {
536 dist= *offset;
537 for (k= dim; k--; )
538 dist += *coordp++ * *normalp++;
539 if (dist < -qh WIDEfacet) {
540 zinc_(Znoarea);
541 return 0.0;
542 }
543 coordp= vertex->point;
544 normalp= normal;
545 for (k= dim; k--; )
546 *(gmcoord++)= (*coordp++ - dist * *normalp++) - *coorda++;
547 }
548 }
549 if (i != dim-1) {
550 fprintf (qh ferr, "qhull internal error (qh_facetarea_simplex): #points %d != dim %d -1\n",
551 i, dim);
552 qh_errexit (qh_ERRqhull, NULL, NULL);
553 }
554 rows[i]= gmcoord;
555 if (qh DELAUNAY) {
556 for (i= 0; i < dim-1; i++)
557 rows[i][dim-1]= 0.0;
558 for (k= dim; k--; )
559 *(gmcoord++)= 0.0;
560 rows[dim-1][dim-1]= -1.0;
561 }else {
562 normalp= normal;
563 for (k= dim; k--; )
564 *(gmcoord++)= *normalp++;
565 }
566 zinc_(Zdetsimplex);
567 area= qh_determinant (rows, dim, &nearzero);
568 if (toporient)
569 area= -area;
570 area *= qh AREAfactor;
571 trace4((qh ferr, "qh_facetarea_simplex: area=%2.2g for point p%d, toporient %d, nearzero? %d\n",
572 area, qh_pointid(apex), toporient, nearzero));
573 return area;
574 } /* facetarea_simplex */
575
576 /*-<a href="qh-geom.htm#TOC"
577 >-------------------------------</a><a name="facetcenter">-</a>
578
579 qh_facetcenter( vertices )
580 return Voronoi center (Voronoi vertex) for a facet's vertices
581
582 returns:
583 return temporary point equal to the center
584
585 see:
586 qh_voronoi_center()
587 */
qh_facetcenter(setT * vertices)588 pointT *qh_facetcenter (setT *vertices) {
589 setT *points= qh_settemp (qh_setsize (vertices));
590 vertexT *vertex, **vertexp;
591 pointT *center;
592
593 FOREACHvertex_(vertices)
594 qh_setappend (&points, vertex->point);
595 center= qh_voronoi_center (qh hull_dim-1, points);
596 qh_settempfree (&points);
597 return center;
598 } /* facetcenter */
599
600 /*-<a href="qh-geom.htm#TOC"
601 >-------------------------------</a><a name="findgooddist">-</a>
602
603 qh_findgooddist( point, facetA, dist, facetlist )
604 find best good facet visible for point from facetA
605 assumes facetA is visible from point
606
607 returns:
608 best facet, i.e., good facet that is furthest from point
609 distance to best facet
610 NULL if none
611
612 moves good, visible facets (and some other visible facets)
613 to end of qh facet_list
614
615 notes:
616 uses qh visit_id
617
618 design:
619 initialize bestfacet if facetA is good
620 move facetA to end of facetlist
621 for each facet on facetlist
622 for each unvisited neighbor of facet
623 move visible neighbors to end of facetlist
624 update best good neighbor
625 if no good neighbors, update best facet
626 */
qh_findgooddist(pointT * point,facetT * facetA,realT * distp,facetT ** facetlist)627 facetT *qh_findgooddist (pointT *point, facetT *facetA, realT *distp,
628 facetT **facetlist) {
629 realT bestdist= -REALmax, dist;
630 facetT *neighbor, **neighborp, *bestfacet=NULL, *facet;
631 boolT goodseen= False;
632
633 if (facetA->good) {
634 zinc_(Zcheckpart); /* calls from check_bestdist occur after print stats */
635 qh_distplane (point, facetA, &bestdist);
636 bestfacet= facetA;
637 goodseen= True;
638 }
639 qh_removefacet (facetA);
640 qh_appendfacet (facetA);
641 *facetlist= facetA;
642 facetA->visitid= ++qh visit_id;
643 FORALLfacet_(*facetlist) {
644 FOREACHneighbor_(facet) {
645 if (neighbor->visitid == qh visit_id)
646 continue;
647 neighbor->visitid= qh visit_id;
648 if (goodseen && !neighbor->good)
649 continue;
650 zinc_(Zcheckpart);
651 qh_distplane (point, neighbor, &dist);
652 if (dist > 0) {
653 qh_removefacet (neighbor);
654 qh_appendfacet (neighbor);
655 if (neighbor->good) {
656 goodseen= True;
657 if (dist > bestdist) {
658 bestdist= dist;
659 bestfacet= neighbor;
660 }
661 }
662 }
663 }
664 }
665 if (bestfacet) {
666 *distp= bestdist;
667 trace2((qh ferr, "qh_findgooddist: p%d is %2.2g above good facet f%d\n",
668 qh_pointid(point), bestdist, bestfacet->id));
669 return bestfacet;
670 }
671 trace4((qh ferr, "qh_findgooddist: no good facet for p%d above f%d\n",
672 qh_pointid(point), facetA->id));
673 return NULL;
674 } /* findgooddist */
675
676 /*-<a href="qh-geom.htm#TOC"
677 >-------------------------------</a><a name="getarea">-</a>
678
679 qh_getarea( facetlist )
680 set area of all facets in facetlist
681 collect statistics
682
683 returns:
684 sets qh totarea/totvol to total area and volume of convex hull
685 for Delaunay triangulation, computes projected area of the lower or upper hull
686 ignores upper hull if qh ATinfinity
687
688 notes:
689 could compute outer volume by expanding facet area by rays from interior
690 the following attempt at perpendicular projection underestimated badly:
691 qh.totoutvol += (-dist + facet->maxoutside + qh DISTround)
692 * area/ qh hull_dim;
693 design:
694 for each facet on facetlist
695 compute facet->area
696 update qh.totarea and qh.totvol
697 */
qh_getarea(facetT * facetlist)698 void qh_getarea (facetT *facetlist) {
699 realT area;
700 realT dist;
701 facetT *facet;
702
703 if (qh REPORTfreq)
704 fprintf (qh ferr, "computing area of each facet and volume of the convex hull\n");
705 else
706 trace1((qh ferr, "qh_getarea: computing volume and area for each facet\n"));
707 qh totarea= qh totvol= 0.0;
708 FORALLfacet_(facetlist) {
709 if (!facet->normal)
710 continue;
711 if (facet->upperdelaunay && qh ATinfinity)
712 continue;
713 facet->f.area= area= qh_facetarea (facet);
714 facet->isarea= True;
715 if (qh DELAUNAY) {
716 if (facet->upperdelaunay == qh UPPERdelaunay)
717 qh totarea += area;
718 }else {
719 qh totarea += area;
720 qh_distplane (qh interior_point, facet, &dist);
721 qh totvol += -dist * area/ qh hull_dim;
722 }
723 if (qh PRINTstatistics) {
724 wadd_(Wareatot, area);
725 wmax_(Wareamax, area);
726 wmin_(Wareamin, area);
727 }
728 }
729 } /* getarea */
730
731 /*-<a href="qh-geom.htm#TOC"
732 >-------------------------------</a><a name="gram_schmidt">-</a>
733
734 qh_gram_schmidt( dim, row )
735 implements Gram-Schmidt orthogonalization by rows
736
737 returns:
738 false if zero norm
739 overwrites rows[dim][dim]
740
741 notes:
742 see Golub & van Loan Algorithm 6.2-2
743 overflow due to small divisors not handled
744
745 design:
746 for each row
747 compute norm for row
748 if non-zero, normalize row
749 for each remaining rowA
750 compute inner product of row and rowA
751 reduce rowA by row * inner product
752 */
qh_gram_schmidt(int dim,realT ** row)753 boolT qh_gram_schmidt(int dim, realT **row) {
754 realT *rowi, *rowj, norm;
755 int i, j, k;
756
757 for(i=0; i < dim; i++) {
758 rowi= row[i];
759 for (norm= 0.0, k= dim; k--; rowi++)
760 norm += *rowi * *rowi;
761 norm= sqrt(norm);
762 wmin_(Wmindenom, norm);
763 if (norm == 0.0) /* either 0 or overflow due to sqrt */
764 return False;
765 for(k= dim; k--; )
766 *(--rowi) /= norm;
767 for(j= i+1; j < dim; j++) {
768 rowj= row[j];
769 for(norm= 0.0, k=dim; k--; )
770 norm += *rowi++ * *rowj++;
771 for(k=dim; k--; )
772 *(--rowj) -= *(--rowi) * norm;
773 }
774 }
775 return True;
776 } /* gram_schmidt */
777
778
779 /*-<a href="qh-geom.htm#TOC"
780 >-------------------------------</a><a name="inthresholds">-</a>
781
782 qh_inthresholds( normal, angle )
783 return True if normal within qh.lower_/upper_threshold
784
785 returns:
786 estimate of angle by summing of threshold diffs
787 angle may be NULL
788 smaller "angle" is better
789
790 notes:
791 invalid if qh.SPLITthresholds
792
793 see:
794 qh.lower_threshold in qh_initbuild()
795 qh_initthresholds()
796
797 design:
798 for each dimension
799 test threshold
800 */
qh_inthresholds(coordT * normal,realT * angle)801 boolT qh_inthresholds (coordT *normal, realT *angle) {
802 boolT within= True;
803 int k;
804 realT threshold;
805
806 if (angle)
807 *angle= 0.0;
808 for(k= 0; k < qh hull_dim; k++) {
809 threshold= qh lower_threshold[k];
810 if (threshold > -REALmax/2) {
811 if (normal[k] < threshold)
812 within= False;
813 if (angle) {
814 threshold -= normal[k];
815 *angle += fabs_(threshold);
816 }
817 }
818 if (qh upper_threshold[k] < REALmax/2) {
819 threshold= qh upper_threshold[k];
820 if (normal[k] > threshold)
821 within= False;
822 if (angle) {
823 threshold -= normal[k];
824 *angle += fabs_(threshold);
825 }
826 }
827 }
828 return within;
829 } /* inthresholds */
830
831
832 /*-<a href="qh-geom.htm#TOC"
833 >-------------------------------</a><a name="joggleinput">-</a>
834
835 qh_joggleinput()
836 randomly joggle input to Qhull by qh.JOGGLEmax
837 initial input is qh.first_point/qh.num_points of qh.hull_dim
838 repeated calls use qh.input_points/qh.num_points
839
840 returns:
841 joggles points at qh.first_point/qh.num_points
842 copies data to qh.input_points/qh.input_malloc if first time
843 determines qh.JOGGLEmax if it was zero
844 if qh.DELAUNAY
845 computes the Delaunay projection of the joggled points
846
847 notes:
848 if qh.DELAUNAY, unnecessarily joggles the last coordinate
849 the initial 'QJn' may be set larger than qh_JOGGLEmaxincrease
850
851 design:
852 if qh.DELAUNAY
853 set qh.SCALElast for reduced precision errors
854 if first call
855 initialize qh.input_points to the original input points
856 if qh.JOGGLEmax == 0
857 determine default qh.JOGGLEmax
858 else
859 increase qh.JOGGLEmax according to qh.build_cnt
860 joggle the input by adding a random number in [-qh.JOGGLEmax,qh.JOGGLEmax]
861 if qh.DELAUNAY
862 sets the Delaunay projection
863 */
qh_joggleinput(void)864 void qh_joggleinput (void) {
865 int size, i, seed;
866 coordT *coordp, *inputp;
867 realT randr, randa, randb;
868
869 if (!qh input_points) { /* first call */
870 qh input_points= qh first_point;
871 qh input_malloc= qh POINTSmalloc;
872 size= qh num_points * qh hull_dim * sizeof(coordT);
873 if (!(qh first_point=(coordT*)malloc(size))) {
874 fprintf(qh ferr, "qhull error: insufficient memory to joggle %d points\n",
875 qh num_points);
876 qh_errexit(qh_ERRmem, NULL, NULL);
877 }
878 qh POINTSmalloc= True;
879 if (qh JOGGLEmax == 0.0) {
880 qh JOGGLEmax= qh_detjoggle (qh input_points, qh num_points, qh hull_dim);
881 qh_option ("QJoggle", NULL, &qh JOGGLEmax);
882 }
883 }else { /* repeated call */
884 if (!qh RERUN && qh build_cnt > qh_JOGGLEretry) {
885 if (((qh build_cnt-qh_JOGGLEretry-1) % qh_JOGGLEagain) == 0) {
886 realT maxjoggle= qh MAXwidth * qh_JOGGLEmaxincrease;
887 if (qh JOGGLEmax < maxjoggle) {
888 qh JOGGLEmax *= qh_JOGGLEincrease;
889 minimize_(qh JOGGLEmax, maxjoggle);
890 }
891 }
892 }
893 qh_option ("QJoggle", NULL, &qh JOGGLEmax);
894 }
895 if (qh build_cnt > 1 && qh JOGGLEmax > fmax_(qh MAXwidth/4, 0.1)) {
896 fprintf (qh ferr, "qhull error: the current joggle for 'QJn', %.2g, is too large for the width\nof the input. If possible, recompile Qhull with higher-precision reals.\n",
897 qh JOGGLEmax);
898 qh_errexit (qh_ERRqhull, NULL, NULL);
899 }
900 /* for some reason, using qh ROTATErandom and qh_RANDOMseed does not repeat the run. Use 'TRn' instead */
901 seed= qh_RANDOMint;
902 qh_option ("_joggle-seed", &seed, NULL);
903 trace0((qh ferr, "qh_joggleinput: joggle input by %2.2g with seed %d\n",
904 qh JOGGLEmax, seed));
905 inputp= qh input_points;
906 coordp= qh first_point;
907 randa= 2.0 * qh JOGGLEmax/qh_RANDOMmax;
908 randb= -qh JOGGLEmax;
909 size= qh num_points * qh hull_dim;
910 for (i= size; i--; ) {
911 randr= qh_RANDOMint;
912 *(coordp++)= *(inputp++) + (randr * randa + randb);
913 }
914 if (qh DELAUNAY) {
915 qh last_low= qh last_high= qh last_newhigh= REALmax;
916 qh_setdelaunay (qh hull_dim, qh num_points, qh first_point);
917 }
918 } /* joggleinput */
919
920 /*-<a href="qh-geom.htm#TOC"
921 >-------------------------------</a><a name="maxabsval">-</a>
922
923 qh_maxabsval( normal, dim )
924 return pointer to maximum absolute value of a dim vector
925 returns NULL if dim=0
926 */
qh_maxabsval(realT * normal,int dim)927 realT *qh_maxabsval (realT *normal, int dim) {
928 realT maxval= -REALmax;
929 realT *maxp= NULL, *colp, absval;
930 int k;
931
932 for (k= dim, colp= normal; k--; colp++) {
933 absval= fabs_(*colp);
934 if (absval > maxval) {
935 maxval= absval;
936 maxp= colp;
937 }
938 }
939 return maxp;
940 } /* maxabsval */
941
942
943 /*-<a href="qh-geom.htm#TOC"
944 >-------------------------------</a><a name="maxmin">-</a>
945
946 qh_maxmin( points, numpoints, dimension )
947 return max/min points for each dimension
948 determine max and min coordinates
949
950 returns:
951 returns a temporary set of max and min points
952 may include duplicate points. Does not include qh.GOODpoint
953 sets qh.NEARzero, qh.MAXabs_coord, qh.MAXsumcoord, qh.MAXwidth
954 qh.MAXlastcoord, qh.MINlastcoord
955 initializes qh.max_outside, qh.min_vertex, qh.WAScoplanar, qh.ZEROall_ok
956
957 notes:
958 loop duplicated in qh_detjoggle()
959
960 design:
961 initialize global precision variables
962 checks definition of REAL...
963 for each dimension
964 for each point
965 collect maximum and minimum point
966 collect maximum of maximums and minimum of minimums
967 determine qh.NEARzero for Gaussian Elimination
968 */
qh_maxmin(pointT * points,int numpoints,int dimension)969 setT *qh_maxmin(pointT *points, int numpoints, int dimension) {
970 int k;
971 realT maxcoord, temp;
972 pointT *minimum, *maximum, *point, *pointtemp;
973 setT *set;
974
975 qh max_outside= 0.0;
976 qh MAXabs_coord= 0.0;
977 qh MAXwidth= -REALmax;
978 qh MAXsumcoord= 0.0;
979 qh min_vertex= 0.0;
980 qh WAScoplanar= False;
981 if (qh ZEROcentrum)
982 qh ZEROall_ok= True;
983 if (REALmin < REALepsilon && REALmin < REALmax && REALmin > -REALmax
984 && REALmax > 0.0 && -REALmax < 0.0)
985 ; /* all ok */
986 else {
987 fprintf (qh ferr, "qhull error: floating point constants in user.h are wrong\n\
988 REALepsilon %g REALmin %g REALmax %g -REALmax %g\n",
989 REALepsilon, REALmin, REALmax, -REALmax);
990 qh_errexit (qh_ERRinput, NULL, NULL);
991 }
992 set= qh_settemp(2*dimension);
993 for(k= 0; k < dimension; k++) {
994 if (points == qh GOODpointp)
995 minimum= maximum= points + dimension;
996 else
997 minimum= maximum= points;
998 FORALLpoint_(points, numpoints) {
999 if (point == qh GOODpointp)
1000 continue;
1001 if (maximum[k] < point[k])
1002 maximum= point;
1003 else if (minimum[k] > point[k])
1004 minimum= point;
1005 }
1006 if (k == dimension-1) {
1007 qh MINlastcoord= minimum[k];
1008 qh MAXlastcoord= maximum[k];
1009 }
1010 if (qh SCALElast && k == dimension-1)
1011 maxcoord= qh MAXwidth;
1012 else {
1013 maxcoord= fmax_(maximum[k], -minimum[k]);
1014 if (qh GOODpointp) {
1015 temp= fmax_(qh GOODpointp[k], -qh GOODpointp[k]);
1016 maximize_(maxcoord, temp);
1017 }
1018 temp= maximum[k] - minimum[k];
1019 maximize_(qh MAXwidth, temp);
1020 }
1021 maximize_(qh MAXabs_coord, maxcoord);
1022 qh MAXsumcoord += maxcoord;
1023 qh_setappend (&set, maximum);
1024 qh_setappend (&set, minimum);
1025 /* calculation of qh NEARzero is based on error formula 4.4-13 of
1026 Golub & van Loan, authors say n^3 can be ignored and 10 be used in
1027 place of rho */
1028 qh NEARzero[k]= 80 * qh MAXsumcoord * REALepsilon;
1029 }
1030 if (qh IStracing >=1)
1031 qh_printpoints (qh ferr, "qh_maxmin: found the max and min points (by dim):", set);
1032 return(set);
1033 } /* maxmin */
1034
1035 /*-<a href="qh-geom.htm#TOC"
1036 >-------------------------------</a><a name="maxouter">-</a>
1037
1038 qh_maxouter()
1039 return maximum distance from facet to outer plane
1040 normally this is qh.max_outside+qh.DISTround
1041 does not include qh.JOGGLEmax
1042
1043 see:
1044 qh_outerinner()
1045
1046 notes:
1047 need to add another qh.DISTround if testing actual point with computation
1048
1049 for joggle:
1050 qh_setfacetplane() updated qh.max_outer for Wnewvertexmax (max distance to vertex)
1051 need to use Wnewvertexmax since could have a coplanar point for a high
1052 facet that is replaced by a low facet
1053 need to add qh.JOGGLEmax if testing input points
1054 */
qh_maxouter(void)1055 realT qh_maxouter (void) {
1056 realT dist;
1057
1058 dist= fmax_(qh max_outside, qh DISTround);
1059 dist += qh DISTround;
1060 trace4((qh ferr, "qh_maxouter: max distance from facet to outer plane is %2.2g max_outside is %2.2g\n", dist, qh max_outside));
1061 return dist;
1062 } /* maxouter */
1063
1064 /*-<a href="qh-geom.htm#TOC"
1065 >-------------------------------</a><a name="maxsimplex">-</a>
1066
1067 qh_maxsimplex( dim, maxpoints, points, numpoints, simplex )
1068 determines maximum simplex for a set of points
1069 starts from points already in simplex
1070 skips qh.GOODpointp (assumes that it isn't in maxpoints)
1071
1072 returns:
1073 simplex with dim+1 points
1074
1075 notes:
1076 assumes at least pointsneeded points in points
1077 maximizes determinate for x,y,z,w, etc.
1078 uses maxpoints as long as determinate is clearly non-zero
1079
1080 design:
1081 initialize simplex with at least two points
1082 (find points with max or min x coordinate)
1083 for each remaining dimension
1084 add point that maximizes the determinate
1085 (use points from maxpoints first)
1086 */
qh_maxsimplex(int dim,setT * maxpoints,pointT * points,int numpoints,setT ** simplex)1087 void qh_maxsimplex (int dim, setT *maxpoints, pointT *points, int numpoints, setT **simplex) {
1088 pointT *point, **pointp, *pointtemp, *maxpoint, *minx=NULL, *maxx=NULL;
1089 boolT nearzero, maxnearzero= False;
1090 int k, sizinit;
1091 realT maxdet= -REALmax, det, mincoord= REALmax, maxcoord= -REALmax;
1092
1093 sizinit= qh_setsize (*simplex);
1094 if (sizinit < 2) {
1095 if (qh_setsize (maxpoints) >= 2) {
1096 FOREACHpoint_(maxpoints) {
1097 if (maxcoord < point[0]) {
1098 maxcoord= point[0];
1099 maxx= point;
1100 }
1101 if (mincoord > point[0]) {
1102 mincoord= point[0];
1103 minx= point;
1104 }
1105 }
1106 }else {
1107 FORALLpoint_(points, numpoints) {
1108 if (point == qh GOODpointp)
1109 continue;
1110 if (maxcoord < point[0]) {
1111 maxcoord= point[0];
1112 maxx= point;
1113 }
1114 if (mincoord > point[0]) {
1115 mincoord= point[0];
1116 minx= point;
1117 }
1118 }
1119 }
1120 qh_setunique (simplex, minx);
1121 if (qh_setsize (*simplex) < 2)
1122 qh_setunique (simplex, maxx);
1123 sizinit= qh_setsize (*simplex);
1124 if (sizinit < 2) {
1125 qh_precision ("input has same x coordinate");
1126 if (zzval_(Zsetplane) > qh hull_dim+1) {
1127 fprintf (qh ferr, "qhull precision error (qh_maxsimplex for voronoi_center):\n%d points with the same x coordinate.\n",
1128 qh_setsize(maxpoints)+numpoints);
1129 qh_errexit (qh_ERRprec, NULL, NULL);
1130 }else {
1131 fprintf (qh ferr, "qhull input error: input is less than %d-dimensional since it has the same x coordinate\n", qh hull_dim);
1132 qh_errexit (qh_ERRinput, NULL, NULL);
1133 }
1134 }
1135 }
1136 for(k= sizinit; k < dim+1; k++) {
1137 maxpoint= NULL;
1138 maxdet= -REALmax;
1139 FOREACHpoint_(maxpoints) {
1140 if (!qh_setin (*simplex, point)) {
1141 det= qh_detsimplex(point, *simplex, k, &nearzero);
1142 if ((det= fabs_(det)) > maxdet) {
1143 maxdet= det;
1144 maxpoint= point;
1145 maxnearzero= nearzero;
1146 }
1147 }
1148 }
1149 if (!maxpoint || maxnearzero) {
1150 zinc_(Zsearchpoints);
1151 if (!maxpoint) {
1152 trace0((qh ferr, "qh_maxsimplex: searching all points for %d-th initial vertex.\n", k+1));
1153 }else {
1154 trace0((qh ferr, "qh_maxsimplex: searching all points for %d-th initial vertex, better than p%d det %2.2g\n",
1155 k+1, qh_pointid(maxpoint), maxdet));
1156 }
1157 FORALLpoint_(points, numpoints) {
1158 if (point == qh GOODpointp)
1159 continue;
1160 if (!qh_setin (*simplex, point)) {
1161 det= qh_detsimplex(point, *simplex, k, &nearzero);
1162 if ((det= fabs_(det)) > maxdet) {
1163 maxdet= det;
1164 maxpoint= point;
1165 maxnearzero= nearzero;
1166 }
1167 }
1168 }
1169 } /* !maxpoint */
1170 if (!maxpoint) {
1171 fprintf (qh ferr, "qhull internal error (qh_maxsimplex): not enough points available\n");
1172 qh_errexit (qh_ERRqhull, NULL, NULL);
1173 }
1174 qh_setappend(simplex, maxpoint);
1175 trace1((qh ferr, "qh_maxsimplex: selected point p%d for %d`th initial vertex, det=%2.2g\n",
1176 qh_pointid(maxpoint), k+1, maxdet));
1177 } /* k */
1178 } /* maxsimplex */
1179
1180 /*-<a href="qh-geom.htm#TOC"
1181 >-------------------------------</a><a name="minabsval">-</a>
1182
1183 qh_minabsval( normal, dim )
1184 return minimum absolute value of a dim vector
1185 */
qh_minabsval(realT * normal,int dim)1186 realT qh_minabsval (realT *normal, int dim) {
1187 realT minval= 0;
1188 realT maxval= 0;
1189 realT *colp;
1190 int k;
1191
1192 for (k= dim, colp= normal; k--; colp++) {
1193 maximize_(maxval, *colp);
1194 minimize_(minval, *colp);
1195 }
1196 return fmax_(maxval, -minval);
1197 } /* minabsval */
1198
1199
1200 /*-<a href="qh-geom.htm#TOC"
1201 >-------------------------------</a><a name="mindiff">-</a>
1202
1203 qh_mindif( vecA, vecB, dim )
1204 return index of min abs. difference of two vectors
1205 */
qh_mindiff(realT * vecA,realT * vecB,int dim)1206 int qh_mindiff (realT *vecA, realT *vecB, int dim) {
1207 realT mindiff= REALmax, diff;
1208 realT *vecAp= vecA, *vecBp= vecB;
1209 int k, mink= 0;
1210
1211 for (k= 0; k < dim; k++) {
1212 diff= *vecAp++ - *vecBp++;
1213 diff= fabs_(diff);
1214 if (diff < mindiff) {
1215 mindiff= diff;
1216 mink= k;
1217 }
1218 }
1219 return mink;
1220 } /* mindiff */
1221
1222
1223
1224 /*-<a href="qh-geom.htm#TOC"
1225 >-------------------------------</a><a name="orientoutside">-</a>
1226
1227 qh_orientoutside( facet )
1228 make facet outside oriented via qh.interior_point
1229
1230 returns:
1231 True if facet reversed orientation.
1232 */
qh_orientoutside(facetT * facet)1233 boolT qh_orientoutside (facetT *facet) {
1234 int k;
1235 realT dist;
1236
1237 qh_distplane (qh interior_point, facet, &dist);
1238 if (dist > 0) {
1239 for (k= qh hull_dim; k--; )
1240 facet->normal[k]= -facet->normal[k];
1241 facet->offset= -facet->offset;
1242 return True;
1243 }
1244 return False;
1245 } /* orientoutside */
1246
1247 /*-<a href="qh-geom.htm#TOC"
1248 >-------------------------------</a><a name="outerinner">-</a>
1249
1250 qh_outerinner( facet, outerplane, innerplane )
1251 if facet and qh.maxoutdone (i.e., qh_check_maxout)
1252 returns outer and inner plane for facet
1253 else
1254 returns maximum outer and inner plane
1255 accounts for qh.JOGGLEmax
1256
1257 see:
1258 qh_maxouter(), qh_check_bestdist(), qh_check_points()
1259
1260 notes:
1261 outerplaner or innerplane may be NULL
1262
1263 includes qh.DISTround for actual points
1264 adds another qh.DISTround if testing with floating point arithmetic
1265 */
qh_outerinner(facetT * facet,realT * outerplane,realT * innerplane)1266 void qh_outerinner (facetT *facet, realT *outerplane, realT *innerplane) {
1267 realT dist, mindist;
1268 vertexT *vertex, **vertexp;
1269
1270 if (outerplane) {
1271 if (!qh_MAXoutside || !facet || !qh maxoutdone) {
1272 *outerplane= qh_maxouter(); /* includes qh.DISTround */
1273 }else { /* qh_MAXoutside ... */
1274 #if qh_MAXoutside
1275 *outerplane= facet->maxoutside + qh DISTround;
1276 #endif
1277
1278 }
1279 if (qh JOGGLEmax < REALmax/2)
1280 *outerplane += qh JOGGLEmax * sqrt (qh hull_dim);
1281 }
1282 if (innerplane) {
1283 if (facet) {
1284 mindist= REALmax;
1285 FOREACHvertex_(facet->vertices) {
1286 zinc_(Zdistio);
1287 qh_distplane (vertex->point, facet, &dist);
1288 minimize_(mindist, dist);
1289 }
1290 *innerplane= mindist - qh DISTround;
1291 }else
1292 *innerplane= qh min_vertex - qh DISTround;
1293 if (qh JOGGLEmax < REALmax/2)
1294 *innerplane -= qh JOGGLEmax * sqrt (qh hull_dim);
1295 }
1296 } /* outerinner */
1297
1298 /*-<a href="qh-geom.htm#TOC"
1299 >-------------------------------</a><a name="pointdist">-</a>
1300
1301 qh_pointdist( point1, point2, dim )
1302 return distance between two points
1303
1304 notes:
1305 returns distance squared if 'dim' is negative
1306 */
qh_pointdist(pointT * point1,pointT * point2,int dim)1307 coordT qh_pointdist(pointT *point1, pointT *point2, int dim) {
1308 coordT dist, diff;
1309 int k;
1310
1311 dist= 0.0;
1312 for (k= (dim > 0 ? dim : -dim); k--; ) {
1313 diff= *point1++ - *point2++;
1314 dist += diff * diff;
1315 }
1316 if (dim > 0)
1317 return(sqrt(dist));
1318 return dist;
1319 } /* pointdist */
1320
1321
1322 /*-<a href="qh-geom.htm#TOC"
1323 >-------------------------------</a><a name="printmatrix">-</a>
1324
1325 qh_printmatrix( fp, string, rows, numrow, numcol )
1326 print matrix to fp given by row vectors
1327 print string as header
1328
1329 notes:
1330 print a vector by qh_printmatrix(fp, "", &vect, 1, len)
1331 */
qh_printmatrix(FILE * fp,char * string,realT ** rows,int numrow,int numcol)1332 void qh_printmatrix (FILE *fp, char *string, realT **rows, int numrow, int numcol) {
1333 realT *rowp;
1334 realT r; /*bug fix*/
1335 int i,k;
1336
1337 fprintf (fp, "%s\n", string);
1338 for (i= 0; i < numrow; i++) {
1339 rowp= rows[i];
1340 for (k= 0; k < numcol; k++) {
1341 r= *rowp++;
1342 fprintf (fp, "%6.3g ", r);
1343 }
1344 fprintf (fp, "\n");
1345 }
1346 } /* printmatrix */
1347
1348
1349 /*-<a href="qh-geom.htm#TOC"
1350 >-------------------------------</a><a name="printpoints">-</a>
1351
1352 qh_printpoints( fp, string, points )
1353 print pointids to fp for a set of points
1354 if string, prints string and 'p' point ids
1355 */
qh_printpoints(FILE * fp,char * string,setT * points)1356 void qh_printpoints (FILE *fp, char *string, setT *points) {
1357 pointT *point, **pointp;
1358
1359 if (string) {
1360 fprintf (fp, "%s", string);
1361 FOREACHpoint_(points)
1362 fprintf (fp, " p%d", qh_pointid(point));
1363 fprintf (fp, "\n");
1364 }else {
1365 FOREACHpoint_(points)
1366 fprintf (fp, " %d", qh_pointid(point));
1367 fprintf (fp, "\n");
1368 }
1369 } /* printpoints */
1370
1371
1372 /*-<a href="qh-geom.htm#TOC"
1373 >-------------------------------</a><a name="projectinput">-</a>
1374
1375 qh_projectinput()
1376 project input points using qh.lower_bound/upper_bound and qh DELAUNAY
1377 if qh.lower_bound[k]=qh.upper_bound[k]= 0,
1378 removes dimension k
1379 if halfspace intersection
1380 removes dimension k from qh.feasible_point
1381 input points in qh first_point, num_points, input_dim
1382
1383 returns:
1384 new point array in qh first_point of qh hull_dim coordinates
1385 sets qh POINTSmalloc
1386 if qh DELAUNAY
1387 projects points to paraboloid
1388 lowbound/highbound is also projected
1389 if qh ATinfinity
1390 adds point "at-infinity"
1391 if qh POINTSmalloc
1392 frees old point array
1393
1394 notes:
1395 checks that qh.hull_dim agrees with qh.input_dim, PROJECTinput, and DELAUNAY
1396
1397
1398 design:
1399 sets project[k] to -1 (delete), 0 (keep), 1 (add for Delaunay)
1400 determines newdim and newnum for qh hull_dim and qh num_points
1401 projects points to newpoints
1402 projects qh.lower_bound to itself
1403 projects qh.upper_bound to itself
1404 if qh DELAUNAY
1405 if qh ATINFINITY
1406 projects points to paraboloid
1407 computes "infinity" point as vertex average and 10% above all points
1408 else
1409 uses qh_setdelaunay to project points to paraboloid
1410 */
qh_projectinput(void)1411 void qh_projectinput (void) {
1412 int k,i;
1413 int newdim= qh input_dim, newnum= qh num_points;
1414 signed char *project;
1415 int size= (qh input_dim+1)*sizeof(*project);
1416 pointT *newpoints, *coord, *infinity;
1417 realT paraboloid, maxboloid= 0;
1418
1419 project= (signed char*)qh_memalloc (size);
1420 memset ((char*)project, 0, size);
1421 for (k= 0; k < qh input_dim; k++) { /* skip Delaunay bound */
1422 if (qh lower_bound[k] == 0 && qh upper_bound[k] == 0) {
1423 project[k]= -1;
1424 newdim--;
1425 }
1426 }
1427 if (qh DELAUNAY) {
1428 project[k]= 1;
1429 newdim++;
1430 if (qh ATinfinity)
1431 newnum++;
1432 }
1433 if (newdim != qh hull_dim) {
1434 fprintf(qh ferr, "qhull internal error (qh_projectinput): dimension after projection %d != hull_dim %d\n", newdim, qh hull_dim);
1435 qh_errexit(qh_ERRqhull, NULL, NULL);
1436 }
1437 if (!(newpoints=(coordT*)malloc(newnum*newdim*sizeof(coordT)))){
1438 fprintf(qh ferr, "qhull error: insufficient memory to project %d points\n",
1439 qh num_points);
1440 qh_errexit(qh_ERRmem, NULL, NULL);
1441 }
1442 qh_projectpoints (project, qh input_dim+1, qh first_point,
1443 qh num_points, qh input_dim, newpoints, newdim);
1444 trace1((qh ferr, "qh_projectinput: updating lower and upper_bound\n"));
1445 qh_projectpoints (project, qh input_dim+1, qh lower_bound,
1446 1, qh input_dim+1, qh lower_bound, newdim+1);
1447 qh_projectpoints (project, qh input_dim+1, qh upper_bound,
1448 1, qh input_dim+1, qh upper_bound, newdim+1);
1449 if (qh HALFspace) {
1450 if (!qh feasible_point) {
1451 fprintf(qh ferr, "qhull internal error (qh_projectinput): HALFspace defined without qh.feasible_point\n");
1452 qh_errexit(qh_ERRqhull, NULL, NULL);
1453 }
1454 qh_projectpoints (project, qh input_dim, qh feasible_point,
1455 1, qh input_dim, qh feasible_point, newdim);
1456 }
1457 qh_memfree(project, ((qh input_dim+1)*sizeof(*project)));
1458 if (qh POINTSmalloc)
1459 free (qh first_point);
1460 qh first_point= newpoints;
1461 qh POINTSmalloc= True;
1462 if (qh DELAUNAY && qh ATinfinity) {
1463 coord= qh first_point;
1464 infinity= qh first_point + qh hull_dim * qh num_points;
1465 for (k=qh hull_dim-1; k--; )
1466 infinity[k]= 0.0;
1467 for (i=qh num_points; i--; ) {
1468 paraboloid= 0.0;
1469 for (k=0; k < qh hull_dim-1; k++) {
1470 paraboloid += *coord * *coord;
1471 infinity[k] += *coord;
1472 coord++;
1473 }
1474 *(coord++)= paraboloid;
1475 maximize_(maxboloid, paraboloid);
1476 }
1477 /* coord == infinity */
1478 for (k=qh hull_dim-1; k--; )
1479 *(coord++) /= qh num_points;
1480 *(coord++)= maxboloid * 1.1;
1481 qh num_points++;
1482 trace0((qh ferr, "qh_projectinput: projected points to paraboloid for Delaunay\n"));
1483 }else if (qh DELAUNAY) /* !qh ATinfinity */
1484 qh_setdelaunay( qh hull_dim, qh num_points, qh first_point);
1485 } /* projectinput */
1486
1487
1488 /*-<a href="qh-geom.htm#TOC"
1489 >-------------------------------</a><a name="projectpoints">-</a>
1490
1491 qh_projectpoints( project, n, points, numpoints, dim, newpoints, newdim )
1492 project points/numpoints/dim to newpoints/newdim
1493 if project[k] == -1
1494 delete dimension k
1495 if project[k] == 1
1496 add dimension k by duplicating previous column
1497 n is size of project
1498
1499 notes:
1500 newpoints may be points if only adding dimension at end
1501
1502 design:
1503 check that 'project' and 'newdim' agree
1504 for each dimension
1505 if project == -1
1506 skip dimension
1507 else
1508 determine start of column in newpoints
1509 determine start of column in points
1510 if project == +1, duplicate previous column
1511 copy dimension (column) from points to newpoints
1512 */
qh_projectpoints(signed char * project,int n,realT * points,int numpoints,int dim,realT * newpoints,int newdim)1513 void qh_projectpoints (signed char *project, int n, realT *points,
1514 int numpoints, int dim, realT *newpoints, int newdim) {
1515 int testdim= dim, oldk=0, newk=0, i,j=0,k;
1516 realT *newp, *oldp;
1517
1518 for (k= 0; k < n; k++)
1519 testdim += project[k];
1520 if (testdim != newdim) {
1521 fprintf (qh ferr, "qhull internal error (qh_projectpoints): newdim %d should be %d after projection\n",
1522 newdim, testdim);
1523 qh_errexit (qh_ERRqhull, NULL, NULL);
1524 }
1525 for (j= 0; j<n; j++) {
1526 if (project[j] == -1)
1527 oldk++;
1528 else {
1529 newp= newpoints+newk++;
1530 if (project[j] == +1) {
1531 if (oldk >= dim)
1532 continue;
1533 oldp= points+oldk;
1534 }else
1535 oldp= points+oldk++;
1536 for (i=numpoints; i--; ) {
1537 *newp= *oldp;
1538 newp += newdim;
1539 oldp += dim;
1540 }
1541 }
1542 if (oldk >= dim)
1543 break;
1544 }
1545 trace1((qh ferr, "qh_projectpoints: projected %d points from dim %d to dim %d\n",
1546 numpoints, dim, newdim));
1547 } /* projectpoints */
1548
1549
1550 /*-<a href="qh-geom.htm#TOC"
1551 >-------------------------------</a><a name="rand">-</a>
1552
1553 qh_rand()
1554 qh_srand( seed )
1555 generate pseudo-random number between 1 and 2^31 -2
1556
1557 notes:
1558 from Park & Miller's minimimal standard random number generator
1559 Communications of the ACM, 31:1192-1201, 1988.
1560 does not use 0 or 2^31 -1
1561 this is silently enforced by qh_srand()
1562 can make 'Rn' much faster by moving qh_rand to qh_distplane
1563 */
1564 int qh_rand_seed= 1; /* define as global variable instead of using qh */
1565
qh_rand(void)1566 int qh_rand( void) {
1567 #define qh_rand_a 16807
1568 #define qh_rand_m 2147483647
1569 #define qh_rand_q 127773 /* m div a */
1570 #define qh_rand_r 2836 /* m mod a */
1571 int lo, hi, test;
1572 int seed = qh_rand_seed;
1573
1574 hi = seed / qh_rand_q; /* seed div q */
1575 lo = seed % qh_rand_q; /* seed mod q */
1576 test = qh_rand_a * lo - qh_rand_r * hi;
1577 if (test > 0)
1578 seed= test;
1579 else
1580 seed= test + qh_rand_m;
1581 qh_rand_seed= seed;
1582 /* seed = seed < qh_RANDOMmax/2 ? 0 : qh_RANDOMmax; for testing */
1583 /* seed = qh_RANDOMmax; for testing */
1584 return seed;
1585 } /* rand */
1586
qh_srand(int seed)1587 void qh_srand( int seed) {
1588 if (seed < 1)
1589 qh_rand_seed= 1;
1590 else if (seed >= qh_rand_m)
1591 qh_rand_seed= qh_rand_m - 1;
1592 else
1593 qh_rand_seed= seed;
1594 } /* qh_srand */
1595
1596 /*-<a href="qh-geom.htm#TOC"
1597 >-------------------------------</a><a name="randomfactor">-</a>
1598
1599 qh_randomfactor()
1600 return a random factor within qh.RANDOMmax of 1.0
1601
1602 notes:
1603 qh.RANDOMa/b are defined in global.c
1604 */
qh_randomfactor(void)1605 realT qh_randomfactor (void) {
1606 realT randr;
1607
1608 randr= qh_RANDOMint;
1609 return randr * qh RANDOMa + qh RANDOMb;
1610 } /* randomfactor */
1611
1612 /*-<a href="qh-geom.htm#TOC"
1613 >-------------------------------</a><a name="randommatrix">-</a>
1614
1615 qh_randommatrix( buffer, dim, rows )
1616 generate a random dim X dim matrix in range [-1,1]
1617 assumes buffer is [dim+1, dim]
1618
1619 returns:
1620 sets buffer to random numbers
1621 sets rows to rows of buffer
1622 sets row[dim] as scratch row
1623 */
qh_randommatrix(realT * buffer,int dim,realT ** rows)1624 void qh_randommatrix (realT *buffer, int dim, realT **rows) {
1625 int i, k;
1626 realT **rowi, *coord, realr;
1627
1628 coord= buffer;
1629 rowi= rows;
1630 for (i=0; i < dim; i++) {
1631 *(rowi++)= coord;
1632 for (k=0; k < dim; k++) {
1633 realr= qh_RANDOMint;
1634 *(coord++)= 2.0 * realr/(qh_RANDOMmax+1) - 1.0;
1635 }
1636 }
1637 *rowi= coord;
1638 } /* randommatrix */
1639
1640
1641 /*-<a href="qh-geom.htm#TOC"
1642 >-------------------------------</a><a name="rotateinput">-</a>
1643
1644 qh_rotateinput( rows )
1645 rotate input using row matrix
1646 input points given by qh first_point, num_points, hull_dim
1647 assumes rows[dim] is a scratch buffer
1648 if qh POINTSmalloc, overwrites input points, else mallocs a new array
1649
1650 returns:
1651 rotated input
1652 sets qh POINTSmalloc
1653
1654 design:
1655 see qh_rotatepoints
1656 */
qh_rotateinput(realT ** rows)1657 void qh_rotateinput (realT **rows) {
1658
1659 if (!qh POINTSmalloc) {
1660 qh first_point= qh_copypoints (qh first_point, qh num_points, qh hull_dim);
1661 qh POINTSmalloc= True;
1662 }
1663 qh_rotatepoints (qh first_point, qh num_points, qh hull_dim, rows);
1664 } /* rotateinput */
1665
1666 /*-<a href="qh-geom.htm#TOC"
1667 >-------------------------------</a><a name="rotatepoints">-</a>
1668
1669 qh_rotatepoints( points, numpoints, dim, row )
1670 rotate numpoints points by a d-dim row matrix
1671 assumes rows[dim] is a scratch buffer
1672
1673 returns:
1674 rotated points in place
1675
1676 design:
1677 for each point
1678 for each coordinate
1679 use row[dim] to compute partial inner product
1680 for each coordinate
1681 rotate by partial inner product
1682 */
qh_rotatepoints(realT * points,int numpoints,int dim,realT ** row)1683 void qh_rotatepoints (realT *points, int numpoints, int dim, realT **row) {
1684 realT *point, *rowi, *coord= NULL, sum, *newval;
1685 int i,j,k;
1686
1687 if (qh IStracing >= 1)
1688 qh_printmatrix (qh ferr, "qh_rotatepoints: rotate points by", row, dim, dim);
1689 for (point= points, j= numpoints; j--; point += dim) {
1690 newval= row[dim];
1691 for (i= 0; i < dim; i++) {
1692 rowi= row[i];
1693 coord= point;
1694 for (sum= 0.0, k= dim; k--; )
1695 sum += *rowi++ * *coord++;
1696 *(newval++)= sum;
1697 }
1698 for (k= dim; k--; )
1699 *(--coord)= *(--newval);
1700 }
1701 } /* rotatepoints */
1702
1703
1704 /*-<a href="qh-geom.htm#TOC"
1705 >-------------------------------</a><a name="scaleinput">-</a>
1706
1707 qh_scaleinput()
1708 scale input points using qh low_bound/high_bound
1709 input points given by qh first_point, num_points, hull_dim
1710 if qh POINTSmalloc, overwrites input points, else mallocs a new array
1711
1712 returns:
1713 scales coordinates of points to low_bound[k], high_bound[k]
1714 sets qh POINTSmalloc
1715
1716 design:
1717 see qh_scalepoints
1718 */
qh_scaleinput(void)1719 void qh_scaleinput (void) {
1720
1721 if (!qh POINTSmalloc) {
1722 qh first_point= qh_copypoints (qh first_point, qh num_points, qh hull_dim);
1723 qh POINTSmalloc= True;
1724 }
1725 qh_scalepoints (qh first_point, qh num_points, qh hull_dim,
1726 qh lower_bound, qh upper_bound);
1727 } /* scaleinput */
1728
1729 /*-<a href="qh-geom.htm#TOC"
1730 >-------------------------------</a><a name="scalelast">-</a>
1731
1732 qh_scalelast( points, numpoints, dim, low, high, newhigh )
1733 scale last coordinate to [0,m] for Delaunay triangulations
1734 input points given by points, numpoints, dim
1735
1736 returns:
1737 changes scale of last coordinate from [low, high] to [0, newhigh]
1738 overwrites last coordinate of each point
1739 saves low/high/newhigh in qh.last_low, etc. for qh_setdelaunay()
1740
1741 notes:
1742 when called by qh_setdelaunay, low/high may not match actual data
1743
1744 design:
1745 compute scale and shift factors
1746 apply to last coordinate of each point
1747 */
qh_scalelast(coordT * points,int numpoints,int dim,coordT low,coordT high,coordT newhigh)1748 void qh_scalelast (coordT *points, int numpoints, int dim, coordT low,
1749 coordT high, coordT newhigh) {
1750 realT scale, shift;
1751 coordT *coord;
1752 int i;
1753 boolT nearzero= False;
1754
1755 trace4((qh ferr, "qh_scalelast: scale last coordinate from [%2.2g, %2.2g] to [0,%2.2g]\n",
1756 low, high, newhigh));
1757 qh last_low= low;
1758 qh last_high= high;
1759 qh last_newhigh= newhigh;
1760 scale= qh_divzero (newhigh, high - low,
1761 qh MINdenom_1, &nearzero);
1762 if (nearzero) {
1763 if (qh DELAUNAY)
1764 fprintf (qh ferr, "qhull input error: can not scale last coordinate. Input is cocircular\n or cospherical. Use option 'Qz' to add a point at infinity.\n");
1765 else
1766 fprintf (qh ferr, "qhull input error: can not scale last coordinate. New bounds [0, %2.2g] are too wide for\nexisting bounds [%2.2g, %2.2g] (width %2.2g)\n",
1767 newhigh, low, high, high-low);
1768 qh_errexit (qh_ERRinput, NULL, NULL);
1769 }
1770 shift= - low * newhigh / (high-low);
1771 coord= points + dim - 1;
1772 for (i= numpoints; i--; coord += dim)
1773 *coord= *coord * scale + shift;
1774 } /* scalelast */
1775
1776 /*-<a href="qh-geom.htm#TOC"
1777 >-------------------------------</a><a name="scalepoints">-</a>
1778
1779 qh_scalepoints( points, numpoints, dim, newlows, newhighs )
1780 scale points to new lowbound and highbound
1781 retains old bound when newlow= -REALmax or newhigh= +REALmax
1782
1783 returns:
1784 scaled points
1785 overwrites old points
1786
1787 design:
1788 for each coordinate
1789 compute current low and high bound
1790 compute scale and shift factors
1791 scale all points
1792 enforce new low and high bound for all points
1793 */
qh_scalepoints(pointT * points,int numpoints,int dim,realT * newlows,realT * newhighs)1794 void qh_scalepoints (pointT *points, int numpoints, int dim,
1795 realT *newlows, realT *newhighs) {
1796 int i,k;
1797 realT shift, scale, *coord, low, high, newlow, newhigh, mincoord, maxcoord;
1798 boolT nearzero= False;
1799
1800 for (k= 0; k < dim; k++) {
1801 newhigh= newhighs[k];
1802 newlow= newlows[k];
1803 if (newhigh > REALmax/2 && newlow < -REALmax/2)
1804 continue;
1805 low= REALmax;
1806 high= -REALmax;
1807 for (i= numpoints, coord= points+k; i--; coord += dim) {
1808 minimize_(low, *coord);
1809 maximize_(high, *coord);
1810 }
1811 if (newhigh > REALmax/2)
1812 newhigh= high;
1813 if (newlow < -REALmax/2)
1814 newlow= low;
1815 if (qh DELAUNAY && k == dim-1 && newhigh < newlow) {
1816 fprintf (qh ferr, "qhull input error: 'Qb%d' or 'QB%d' inverts paraboloid since high bound %.2g < low bound %.2g\n",
1817 k, k, newhigh, newlow);
1818 qh_errexit (qh_ERRinput, NULL, NULL);
1819 }
1820 scale= qh_divzero (newhigh - newlow, high - low,
1821 qh MINdenom_1, &nearzero);
1822 if (nearzero) {
1823 fprintf (qh ferr, "qhull input error: %d'th dimension's new bounds [%2.2g, %2.2g] too wide for\nexisting bounds [%2.2g, %2.2g]\n",
1824 k, newlow, newhigh, low, high);
1825 qh_errexit (qh_ERRinput, NULL, NULL);
1826 }
1827 shift= (newlow * high - low * newhigh)/(high-low);
1828 coord= points+k;
1829 for (i= numpoints; i--; coord += dim)
1830 *coord= *coord * scale + shift;
1831 coord= points+k;
1832 if (newlow < newhigh) {
1833 mincoord= newlow;
1834 maxcoord= newhigh;
1835 }else {
1836 mincoord= newhigh;
1837 maxcoord= newlow;
1838 }
1839 for (i= numpoints; i--; coord += dim) {
1840 minimize_(*coord, maxcoord); /* because of roundoff error */
1841 maximize_(*coord, mincoord);
1842 }
1843 trace0((qh ferr, "qh_scalepoints: scaled %d'th coordinate [%2.2g, %2.2g] to [%.2g, %.2g] for %d points by %2.2g and shifted %2.2g\n",
1844 k, low, high, newlow, newhigh, numpoints, scale, shift));
1845 }
1846 } /* scalepoints */
1847
1848
1849 /*-<a href="qh-geom.htm#TOC"
1850 >-------------------------------</a><a name="setdelaunay">-</a>
1851
1852 qh_setdelaunay( dim, count, points )
1853 project count points to dim-d paraboloid for Delaunay triangulation
1854
1855 dim is one more than the dimension of the input set
1856 assumes dim is at least 3 (i.e., at least a 2-d Delaunay triangulation)
1857
1858 points is a dim*count realT array. The first dim-1 coordinates
1859 are the coordinates of the first input point. array[dim] is
1860 the first coordinate of the second input point. array[2*dim] is
1861 the first coordinate of the third input point.
1862
1863 if qh.last_low defined (i.e., 'Qbb' called qh_scalelast)
1864 calls qh_scalelast to scale the last coordinate the same as the other points
1865
1866 returns:
1867 for each point
1868 sets point[dim-1] to sum of squares of coordinates
1869 scale points to 'Qbb' if needed
1870
1871 notes:
1872 to project one point, use
1873 qh_setdelaunay (qh hull_dim, 1, point)
1874
1875 Do not use options 'Qbk', 'QBk', or 'QbB' since they scale
1876 the coordinates after the original projection.
1877
1878 */
qh_setdelaunay(int dim,int count,pointT * points)1879 void qh_setdelaunay (int dim, int count, pointT *points) {
1880 int i, k;
1881 coordT *coordp, coord;
1882 realT paraboloid;
1883
1884 trace0((qh ferr, "qh_setdelaunay: project %d points to paraboloid for Delaunay triangulation\n", count));
1885 coordp= points;
1886 for (i= 0; i < count; i++) {
1887 coord= *coordp++;
1888 paraboloid= coord*coord;
1889 for (k= dim-2; k--; ) {
1890 coord= *coordp++;
1891 paraboloid += coord*coord;
1892 }
1893 *coordp++ = paraboloid;
1894 }
1895 if (qh last_low < REALmax/2)
1896 qh_scalelast (points, count, dim, qh last_low, qh last_high, qh last_newhigh);
1897 } /* setdelaunay */
1898
1899
1900 /*-<a href="qh-geom.htm#TOC"
1901 >-------------------------------</a><a name="sethalfspace">-</a>
1902
1903 qh_sethalfspace( dim, coords, nextp, normal, offset, feasible )
1904 set point to dual of halfspace relative to feasible point
1905 halfspace is normal coefficients and offset.
1906
1907 returns:
1908 false if feasible point is outside of hull (error message already reported)
1909 overwrites coordinates for point at dim coords
1910 nextp= next point (coords)
1911
1912 design:
1913 compute distance from feasible point to halfspace
1914 divide each normal coefficient by -dist
1915 */
qh_sethalfspace(int dim,coordT * coords,coordT ** nextp,coordT * normal,coordT * offset,coordT * feasible)1916 boolT qh_sethalfspace (int dim, coordT *coords, coordT **nextp,
1917 coordT *normal, coordT *offset, coordT *feasible) {
1918 coordT *normp= normal, *feasiblep= feasible, *coordp= coords;
1919 realT dist;
1920 realT r; /*bug fix*/
1921 int k;
1922 boolT zerodiv;
1923
1924 dist= *offset;
1925 for (k= dim; k--; )
1926 dist += *(normp++) * *(feasiblep++);
1927 if (dist > 0)
1928 goto LABELerroroutside;
1929 normp= normal;
1930 if (dist < -qh MINdenom) {
1931 for (k= dim; k--; )
1932 *(coordp++)= *(normp++) / -dist;
1933 }else {
1934 for (k= dim; k--; ) {
1935 *(coordp++)= qh_divzero (*(normp++), -dist, qh MINdenom_1, &zerodiv);
1936 if (zerodiv)
1937 goto LABELerroroutside;
1938 }
1939 }
1940 *nextp= coordp;
1941 if (qh IStracing >= 4) {
1942 fprintf (qh ferr, "qh_sethalfspace: halfspace at offset %6.2g to point: ", *offset);
1943 for (k= dim, coordp= coords; k--; ) {
1944 r= *coordp++;
1945 fprintf (qh ferr, " %6.2g", r);
1946 }
1947 fprintf (qh ferr, "\n");
1948 }
1949 return True;
1950 LABELerroroutside:
1951 feasiblep= feasible;
1952 normp= normal;
1953 fprintf(qh ferr, "qhull input error: feasible point is not clearly inside halfspace\nfeasible point: ");
1954 for (k= dim; k--; )
1955 fprintf (qh ferr, qh_REAL_1, r=*(feasiblep++));
1956 fprintf (qh ferr, "\n halfspace: ");
1957 for (k= dim; k--; )
1958 fprintf (qh ferr, qh_REAL_1, r=*(normp++));
1959 fprintf (qh ferr, "\n at offset: ");
1960 fprintf (qh ferr, qh_REAL_1, *offset);
1961 fprintf (qh ferr, " and distance: ");
1962 fprintf (qh ferr, qh_REAL_1, dist);
1963 fprintf (qh ferr, "\n");
1964 return False;
1965 } /* sethalfspace */
1966
1967 /*-<a href="qh-geom.htm#TOC"
1968 >-------------------------------</a><a name="sethalfspace_all">-</a>
1969
1970 qh_sethalfspace_all( dim, count, halfspaces, feasible )
1971 generate dual for halfspace intersection with feasible point
1972 array of count halfspaces
1973 each halfspace is normal coefficients followed by offset
1974 the origin is inside the halfspace if the offset is negative
1975
1976 returns:
1977 malloc'd array of count X dim-1 points
1978
1979 notes:
1980 call before qh_init_B or qh_initqhull_globals
1981 unused/untested code: please email bradb@shore.net if this works ok for you
1982 If using option 'Fp', also set qh feasible_point. It is a malloc'd array
1983 that is freed by qh_freebuffers.
1984
1985 design:
1986 see qh_sethalfspace
1987 */
qh_sethalfspace_all(int dim,int count,coordT * halfspaces,pointT * feasible)1988 coordT *qh_sethalfspace_all (int dim, int count, coordT *halfspaces, pointT *feasible) {
1989 int i, newdim;
1990 pointT *newpoints;
1991 coordT *coordp, *normalp, *offsetp;
1992
1993 trace0((qh ferr, "qh_sethalfspace_all: compute dual for halfspace intersection\n"));
1994 newdim= dim - 1;
1995 if (!(newpoints=(coordT*)malloc(count*newdim*sizeof(coordT)))){
1996 fprintf(qh ferr, "qhull error: insufficient memory to compute dual of %d halfspaces\n",
1997 count);
1998 qh_errexit(qh_ERRmem, NULL, NULL);
1999 }
2000 coordp= newpoints;
2001 normalp= halfspaces;
2002 for (i= 0; i < count; i++) {
2003 offsetp= normalp + newdim;
2004 if (!qh_sethalfspace (newdim, coordp, &coordp, normalp, offsetp, feasible)) {
2005 fprintf (qh ferr, "The halfspace was at index %d\n", i);
2006 qh_errexit (qh_ERRinput, NULL, NULL);
2007 }
2008 normalp= offsetp + 1;
2009 }
2010 return newpoints;
2011 } /* sethalfspace_all */
2012
2013
2014 /*-<a href="qh-geom.htm#TOC"
2015 >-------------------------------</a><a name="sharpnewfacets">-</a>
2016
2017 qh_sharpnewfacets()
2018
2019 returns:
2020 true if could be an acute angle (facets in different quadrants)
2021
2022 notes:
2023 for qh_findbest
2024
2025 design:
2026 for all facets on qh.newfacet_list
2027 if two facets are in different quadrants
2028 set issharp
2029 */
qh_sharpnewfacets()2030 boolT qh_sharpnewfacets () {
2031 facetT *facet;
2032 boolT issharp = False;
2033 int *quadrant, k;
2034
2035 quadrant= (int*)qh_memalloc (qh hull_dim * sizeof(int));
2036 FORALLfacet_(qh newfacet_list) {
2037 if (facet == qh newfacet_list) {
2038 for (k= qh hull_dim; k--; )
2039 quadrant[ k]= (facet->normal[ k] > 0);
2040 }else {
2041 for (k= qh hull_dim; k--; ) {
2042 if (quadrant[ k] != (facet->normal[ k] > 0)) {
2043 issharp= True;
2044 break;
2045 }
2046 }
2047 }
2048 if (issharp)
2049 break;
2050 }
2051 qh_memfree( quadrant, qh hull_dim * sizeof(int));
2052 trace3((qh ferr, "qh_sharpnewfacets: %d\n", issharp));
2053 return issharp;
2054 } /* sharpnewfacets */
2055
2056 /*-<a href="qh-geom.htm#TOC"
2057 >-------------------------------</a><a name="voronoi_center">-</a>
2058
2059 qh_voronoi_center( dim, points )
2060 return Voronoi center for a set of points
2061 dim is the orginal dimension of the points
2062 gh.gm_matrix/qh.gm_row are scratch buffers
2063
2064 returns:
2065 center as a temporary point
2066 if non-simplicial,
2067 returns center for max simplex of points
2068
2069 notes:
2070 from Bowyer & Woodwark, A Programmer's Geometry, 1983, p. 65
2071
2072 design:
2073 if non-simplicial
2074 determine max simplex for points
2075 translate point0 of simplex to origin
2076 compute sum of squares of diagonal
2077 compute determinate
2078 compute Voronoi center (see Bowyer & Woodwark)
2079 */
qh_voronoi_center(int dim,setT * points)2080 pointT *qh_voronoi_center (int dim, setT *points) {
2081 pointT *point, **pointp, *point0;
2082 pointT *center= (pointT*)qh_memalloc (qh center_size);
2083 setT *simplex;
2084 int i, j, k, size= qh_setsize(points);
2085 coordT *gmcoord;
2086 realT *diffp, sum2, *sum2row, *sum2p, det, factor;
2087 boolT nearzero, infinite;
2088
2089 if (size == dim+1)
2090 simplex= points;
2091 else if (size < dim+1) {
2092 fprintf (qh ferr, "qhull internal error (qh_voronoi_center):\n need at least %d points to construct a Voronoi center\n",
2093 dim+1);
2094 qh_errexit (qh_ERRqhull, NULL, NULL);
2095 }else {
2096 simplex= qh_settemp (dim+1);
2097 qh_maxsimplex (dim, points, NULL, 0, &simplex);
2098 }
2099 point0= SETfirstt_(simplex, pointT);
2100 gmcoord= qh gm_matrix;
2101 for (k=0; k < dim; k++) {
2102 qh gm_row[k]= gmcoord;
2103 FOREACHpoint_(simplex) {
2104 if (point != point0)
2105 *(gmcoord++)= point[k] - point0[k];
2106 }
2107 }
2108 sum2row= gmcoord;
2109 for (i=0; i < dim; i++) {
2110 sum2= 0.0;
2111 for (k= 0; k < dim; k++) {
2112 diffp= qh gm_row[k] + i;
2113 sum2 += *diffp * *diffp;
2114 }
2115 *(gmcoord++)= sum2;
2116 }
2117 det= qh_determinant (qh gm_row, dim, &nearzero);
2118 factor= qh_divzero (0.5, det, qh MINdenom, &infinite);
2119 if (infinite) {
2120 for (k=dim; k--; )
2121 center[k]= qh_INFINITE;
2122 if (qh IStracing)
2123 qh_printpoints (qh ferr, "qh_voronoi_center: at infinity for ", simplex);
2124 }else {
2125 for (i=0; i < dim; i++) {
2126 gmcoord= qh gm_matrix;
2127 sum2p= sum2row;
2128 for (k=0; k < dim; k++) {
2129 qh gm_row[k]= gmcoord;
2130 if (k == i) {
2131 for (j= dim; j--; )
2132 *(gmcoord++)= *sum2p++;
2133 }else {
2134 FOREACHpoint_(simplex) {
2135 if (point != point0)
2136 *(gmcoord++)= point[k] - point0[k];
2137 }
2138 }
2139 }
2140 center[i]= qh_determinant (qh gm_row, dim, &nearzero)*factor + point0[i];
2141 }
2142 #ifndef qh_NOtrace
2143 if (qh IStracing >= 3) {
2144 fprintf (qh ferr, "qh_voronoi_center: det %2.2g factor %2.2g ", det, factor);
2145 qh_printmatrix (qh ferr, "center:", ¢er, 1, dim);
2146 if (qh IStracing >= 5) {
2147 qh_printpoints (qh ferr, "points", simplex);
2148 FOREACHpoint_(simplex)
2149 fprintf (qh ferr, "p%d dist %.2g, ", qh_pointid (point),
2150 qh_pointdist (point, center, dim));
2151 fprintf (qh ferr, "\n");
2152 }
2153 }
2154 #endif
2155 }
2156 if (simplex != points)
2157 qh_settempfree (&simplex);
2158 return center;
2159 } /* voronoi_center */
2160
2161