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