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