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-2015 The Geometry Center.
11 $Id: //main/2015/qhull/src/libqhull/geom2.c#6 $$Change: 2065 $
12 $DateTime: 2016/01/18 13:51:04 $$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 notes:
27 qh_free the returned points to avoid a memory leak
28 */
qh_copypoints(coordT * points,int numpoints,int dimension)29 coordT *qh_copypoints(coordT *points, int numpoints, int dimension) {
30 int size;
31 coordT *newpoints;
32
33 size= numpoints * dimension * (int)sizeof(coordT);
34 if (!(newpoints=(coordT*)qh_malloc((size_t)size))) {
35 qh_fprintf(qh ferr, 6004, "qhull error: insufficient memory to copy %d points\n",
36 numpoints);
37 qh_errexit(qh_ERRmem, NULL, NULL);
38 }
39 memcpy((char *)newpoints, (char *)points, (size_t)size); /* newpoints!=0 by QH6004 */
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) < 10*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) < 10*qh NEARzero[2]) /* what should this be? det 5.5e-12 was flat for qh_maxsimplex of qdelaunay 0,0 27,27 -36,36 -9,63 */
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 Golub & van Loan, 1983, Lemma 3.2-1, "Rounding Errors"
373 use sqrt(dim) since one vector is normalized
374 or use maxsumabs since one vector is < 1
375 */
qh_distround(int dimension,realT maxabs,realT maxsumabs)376 realT qh_distround(int dimension, realT maxabs, realT maxsumabs) {
377 realT maxdistsum, maxround;
378
379 maxdistsum= sqrt((realT)dimension) * maxabs;
380 minimize_( maxdistsum, maxsumabs);
381 maxround= REALepsilon * (dimension * maxdistsum * 1.01 + maxabs);
382 /* adds maxabs for offset */
383 trace4((qh ferr, 4008, "qh_distround: %2.2g maxabs %2.2g maxsumabs %2.2g maxdistsum %2.2g\n",
384 maxround, maxabs, maxsumabs, maxdistsum));
385 return maxround;
386 } /* distround */
387
388 /*-<a href="qh-geom.htm#TOC"
389 >-------------------------------</a><a name="divzero">-</a>
390
391 qh_divzero( numer, denom, mindenom1, zerodiv )
392 divide by a number that's nearly zero
393 mindenom1= minimum denominator for dividing into 1.0
394
395 returns:
396 quotient
397 sets zerodiv and returns 0.0 if it would overflow
398
399 design:
400 if numer is nearly zero and abs(numer) < abs(denom)
401 return numer/denom
402 else if numer is nearly zero
403 return 0 and zerodiv
404 else if denom/numer non-zero
405 return numer/denom
406 else
407 return 0 and zerodiv
408 */
qh_divzero(realT numer,realT denom,realT mindenom1,boolT * zerodiv)409 realT qh_divzero(realT numer, realT denom, realT mindenom1, boolT *zerodiv) {
410 realT temp, numerx, denomx;
411
412
413 if (numer < mindenom1 && numer > -mindenom1) {
414 numerx= fabs_(numer);
415 denomx= fabs_(denom);
416 if (numerx < denomx) {
417 *zerodiv= False;
418 return numer/denom;
419 }else {
420 *zerodiv= True;
421 return 0.0;
422 }
423 }
424 temp= denom/numer;
425 if (temp > mindenom1 || temp < -mindenom1) {
426 *zerodiv= False;
427 return numer/denom;
428 }else {
429 *zerodiv= True;
430 return 0.0;
431 }
432 } /* divzero */
433
434
435 /*-<a href="qh-geom.htm#TOC"
436 >-------------------------------</a><a name="facetarea">-</a>
437
438 qh_facetarea( facet )
439 return area for a facet
440
441 notes:
442 if non-simplicial,
443 uses centrum to triangulate facet and sums the projected areas.
444 if (qh DELAUNAY),
445 computes projected area instead for last coordinate
446 assumes facet->normal exists
447 projecting tricoplanar facets to the hyperplane does not appear to make a difference
448
449 design:
450 if simplicial
451 compute area
452 else
453 for each ridge
454 compute area from centrum to ridge
455 negate area if upper Delaunay facet
456 */
qh_facetarea(facetT * facet)457 realT qh_facetarea(facetT *facet) {
458 vertexT *apex;
459 pointT *centrum;
460 realT area= 0.0;
461 ridgeT *ridge, **ridgep;
462
463 if (facet->simplicial) {
464 apex= SETfirstt_(facet->vertices, vertexT);
465 area= qh_facetarea_simplex(qh hull_dim, apex->point, facet->vertices,
466 apex, facet->toporient, facet->normal, &facet->offset);
467 }else {
468 if (qh CENTERtype == qh_AScentrum)
469 centrum= facet->center;
470 else
471 centrum= qh_getcentrum(facet);
472 FOREACHridge_(facet->ridges)
473 area += qh_facetarea_simplex(qh hull_dim, centrum, ridge->vertices,
474 NULL, (boolT)(ridge->top == facet), facet->normal, &facet->offset);
475 if (qh CENTERtype != qh_AScentrum)
476 qh_memfree(centrum, qh normal_size);
477 }
478 if (facet->upperdelaunay && qh DELAUNAY)
479 area= -area; /* the normal should be [0,...,1] */
480 trace4((qh ferr, 4009, "qh_facetarea: f%d area %2.2g\n", facet->id, area));
481 return area;
482 } /* facetarea */
483
484 /*-<a href="qh-geom.htm#TOC"
485 >-------------------------------</a><a name="facetarea_simplex">-</a>
486
487 qh_facetarea_simplex( dim, apex, vertices, notvertex, toporient, normal, offset )
488 return area for a simplex defined by
489 an apex, a base of vertices, an orientation, and a unit normal
490 if simplicial or tricoplanar facet,
491 notvertex is defined and it is skipped in vertices
492
493 returns:
494 computes area of simplex projected to plane [normal,offset]
495 returns 0 if vertex too far below plane (qh WIDEfacet)
496 vertex can't be apex of tricoplanar facet
497
498 notes:
499 if (qh DELAUNAY),
500 computes projected area instead for last coordinate
501 uses qh gm_matrix/gm_row and qh hull_dim
502 helper function for qh_facetarea
503
504 design:
505 if Notvertex
506 translate simplex to apex
507 else
508 project simplex to normal/offset
509 translate simplex to apex
510 if Delaunay
511 set last row/column to 0 with -1 on diagonal
512 else
513 set last row to Normal
514 compute determinate
515 scale and flip sign for area
516 */
qh_facetarea_simplex(int dim,coordT * apex,setT * vertices,vertexT * notvertex,boolT toporient,coordT * normal,realT * offset)517 realT qh_facetarea_simplex(int dim, coordT *apex, setT *vertices,
518 vertexT *notvertex, boolT toporient, coordT *normal, realT *offset) {
519 pointT *coorda, *coordp, *gmcoord;
520 coordT **rows, *normalp;
521 int k, i=0;
522 realT area, dist;
523 vertexT *vertex, **vertexp;
524 boolT nearzero;
525
526 gmcoord= qh gm_matrix;
527 rows= qh gm_row;
528 FOREACHvertex_(vertices) {
529 if (vertex == notvertex)
530 continue;
531 rows[i++]= gmcoord;
532 coorda= apex;
533 coordp= vertex->point;
534 normalp= normal;
535 if (notvertex) {
536 for (k=dim; k--; )
537 *(gmcoord++)= *coordp++ - *coorda++;
538 }else {
539 dist= *offset;
540 for (k=dim; k--; )
541 dist += *coordp++ * *normalp++;
542 if (dist < -qh WIDEfacet) {
543 zinc_(Znoarea);
544 return 0.0;
545 }
546 coordp= vertex->point;
547 normalp= normal;
548 for (k=dim; k--; )
549 *(gmcoord++)= (*coordp++ - dist * *normalp++) - *coorda++;
550 }
551 }
552 if (i != dim-1) {
553 qh_fprintf(qh ferr, 6008, "qhull internal error (qh_facetarea_simplex): #points %d != dim %d -1\n",
554 i, dim);
555 qh_errexit(qh_ERRqhull, NULL, NULL);
556 }
557 rows[i]= gmcoord;
558 if (qh DELAUNAY) {
559 for (i=0; i < dim-1; i++)
560 rows[i][dim-1]= 0.0;
561 for (k=dim; k--; )
562 *(gmcoord++)= 0.0;
563 rows[dim-1][dim-1]= -1.0;
564 }else {
565 normalp= normal;
566 for (k=dim; k--; )
567 *(gmcoord++)= *normalp++;
568 }
569 zinc_(Zdetsimplex);
570 area= qh_determinant(rows, dim, &nearzero);
571 if (toporient)
572 area= -area;
573 area *= qh AREAfactor;
574 trace4((qh ferr, 4010, "qh_facetarea_simplex: area=%2.2g for point p%d, toporient %d, nearzero? %d\n",
575 area, qh_pointid(apex), toporient, nearzero));
576 return area;
577 } /* facetarea_simplex */
578
579 /*-<a href="qh-geom.htm#TOC"
580 >-------------------------------</a><a name="facetcenter">-</a>
581
582 qh_facetcenter( vertices )
583 return Voronoi center (Voronoi vertex) for a facet's vertices
584
585 returns:
586 return temporary point equal to the center
587
588 see:
589 qh_voronoi_center()
590 */
qh_facetcenter(setT * vertices)591 pointT *qh_facetcenter(setT *vertices) {
592 setT *points= qh_settemp(qh_setsize(vertices));
593 vertexT *vertex, **vertexp;
594 pointT *center;
595
596 FOREACHvertex_(vertices)
597 qh_setappend(&points, vertex->point);
598 center= qh_voronoi_center(qh hull_dim-1, points);
599 qh_settempfree(&points);
600 return center;
601 } /* facetcenter */
602
603 /*-<a href="qh-geom.htm#TOC"
604 >-------------------------------</a><a name="findgooddist">-</a>
605
606 qh_findgooddist( point, facetA, dist, facetlist )
607 find best good facet visible for point from facetA
608 assumes facetA is visible from point
609
610 returns:
611 best facet, i.e., good facet that is furthest from point
612 distance to best facet
613 NULL if none
614
615 moves good, visible facets (and some other visible facets)
616 to end of qh facet_list
617
618 notes:
619 uses qh visit_id
620
621 design:
622 initialize bestfacet if facetA is good
623 move facetA to end of facetlist
624 for each facet on facetlist
625 for each unvisited neighbor of facet
626 move visible neighbors to end of facetlist
627 update best good neighbor
628 if no good neighbors, update best facet
629 */
qh_findgooddist(pointT * point,facetT * facetA,realT * distp,facetT ** facetlist)630 facetT *qh_findgooddist(pointT *point, facetT *facetA, realT *distp,
631 facetT **facetlist) {
632 realT bestdist= -REALmax, dist;
633 facetT *neighbor, **neighborp, *bestfacet=NULL, *facet;
634 boolT goodseen= False;
635
636 if (facetA->good) {
637 zzinc_(Zcheckpart); /* calls from check_bestdist occur after print stats */
638 qh_distplane(point, facetA, &bestdist);
639 bestfacet= facetA;
640 goodseen= True;
641 }
642 qh_removefacet(facetA);
643 qh_appendfacet(facetA);
644 *facetlist= facetA;
645 facetA->visitid= ++qh visit_id;
646 FORALLfacet_(*facetlist) {
647 FOREACHneighbor_(facet) {
648 if (neighbor->visitid == qh visit_id)
649 continue;
650 neighbor->visitid= qh visit_id;
651 if (goodseen && !neighbor->good)
652 continue;
653 zzinc_(Zcheckpart);
654 qh_distplane(point, neighbor, &dist);
655 if (dist > 0) {
656 qh_removefacet(neighbor);
657 qh_appendfacet(neighbor);
658 if (neighbor->good) {
659 goodseen= True;
660 if (dist > bestdist) {
661 bestdist= dist;
662 bestfacet= neighbor;
663 }
664 }
665 }
666 }
667 }
668 if (bestfacet) {
669 *distp= bestdist;
670 trace2((qh ferr, 2003, "qh_findgooddist: p%d is %2.2g above good facet f%d\n",
671 qh_pointid(point), bestdist, bestfacet->id));
672 return bestfacet;
673 }
674 trace4((qh ferr, 4011, "qh_findgooddist: no good facet for p%d above f%d\n",
675 qh_pointid(point), facetA->id));
676 return NULL;
677 } /* findgooddist */
678
679 /*-<a href="qh-geom.htm#TOC"
680 >-------------------------------</a><a name="getarea">-</a>
681
682 qh_getarea( facetlist )
683 set area of all facets in facetlist
684 collect statistics
685 nop if hasAreaVolume
686
687 returns:
688 sets qh totarea/totvol to total area and volume of convex hull
689 for Delaunay triangulation, computes projected area of the lower or upper hull
690 ignores upper hull if qh ATinfinity
691
692 notes:
693 could compute outer volume by expanding facet area by rays from interior
694 the following attempt at perpendicular projection underestimated badly:
695 qh.totoutvol += (-dist + facet->maxoutside + qh DISTround)
696 * area/ qh hull_dim;
697 design:
698 for each facet on facetlist
699 compute facet->area
700 update qh.totarea and qh.totvol
701 */
qh_getarea(facetT * facetlist)702 void qh_getarea(facetT *facetlist) {
703 realT area;
704 realT dist;
705 facetT *facet;
706
707 if (qh hasAreaVolume)
708 return;
709 if (qh REPORTfreq)
710 qh_fprintf(qh ferr, 8020, "computing area of each facet and volume of the convex hull\n");
711 else
712 trace1((qh ferr, 1001, "qh_getarea: computing volume and area for each facet\n"));
713 qh totarea= qh totvol= 0.0;
714 FORALLfacet_(facetlist) {
715 if (!facet->normal)
716 continue;
717 if (facet->upperdelaunay && qh ATinfinity)
718 continue;
719 if (!facet->isarea) {
720 facet->f.area= qh_facetarea(facet);
721 facet->isarea= True;
722 }
723 area= facet->f.area;
724 if (qh DELAUNAY) {
725 if (facet->upperdelaunay == qh UPPERdelaunay)
726 qh totarea += area;
727 }else {
728 qh totarea += area;
729 qh_distplane(qh interior_point, facet, &dist);
730 qh totvol += -dist * area/ qh hull_dim;
731 }
732 if (qh PRINTstatistics) {
733 wadd_(Wareatot, area);
734 wmax_(Wareamax, area);
735 wmin_(Wareamin, area);
736 }
737 }
738 qh hasAreaVolume= True;
739 } /* getarea */
740
741 /*-<a href="qh-geom.htm#TOC"
742 >-------------------------------</a><a name="gram_schmidt">-</a>
743
744 qh_gram_schmidt( dim, row )
745 implements Gram-Schmidt orthogonalization by rows
746
747 returns:
748 false if zero norm
749 overwrites rows[dim][dim]
750
751 notes:
752 see Golub & van Loan, 1983, Algorithm 6.2-2, "Modified Gram-Schmidt"
753 overflow due to small divisors not handled
754
755 design:
756 for each row
757 compute norm for row
758 if non-zero, normalize row
759 for each remaining rowA
760 compute inner product of row and rowA
761 reduce rowA by row * inner product
762 */
qh_gram_schmidt(int dim,realT ** row)763 boolT qh_gram_schmidt(int dim, realT **row) {
764 realT *rowi, *rowj, norm;
765 int i, j, k;
766
767 for (i=0; i < dim; i++) {
768 rowi= row[i];
769 for (norm= 0.0, k= dim; k--; rowi++)
770 norm += *rowi * *rowi;
771 norm= sqrt(norm);
772 wmin_(Wmindenom, norm);
773 if (norm == 0.0) /* either 0 or overflow due to sqrt */
774 return False;
775 for (k=dim; k--; )
776 *(--rowi) /= norm;
777 for (j=i+1; j < dim; j++) {
778 rowj= row[j];
779 for (norm= 0.0, k=dim; k--; )
780 norm += *rowi++ * *rowj++;
781 for (k=dim; k--; )
782 *(--rowj) -= *(--rowi) * norm;
783 }
784 }
785 return True;
786 } /* gram_schmidt */
787
788
789 /*-<a href="qh-geom.htm#TOC"
790 >-------------------------------</a><a name="inthresholds">-</a>
791
792 qh_inthresholds( normal, angle )
793 return True if normal within qh.lower_/upper_threshold
794
795 returns:
796 estimate of angle by summing of threshold diffs
797 angle may be NULL
798 smaller "angle" is better
799
800 notes:
801 invalid if qh.SPLITthresholds
802
803 see:
804 qh.lower_threshold in qh_initbuild()
805 qh_initthresholds()
806
807 design:
808 for each dimension
809 test threshold
810 */
qh_inthresholds(coordT * normal,realT * angle)811 boolT qh_inthresholds(coordT *normal, realT *angle) {
812 boolT within= True;
813 int k;
814 realT threshold;
815
816 if (angle)
817 *angle= 0.0;
818 for (k=0; k < qh hull_dim; k++) {
819 threshold= qh lower_threshold[k];
820 if (threshold > -REALmax/2) {
821 if (normal[k] < threshold)
822 within= False;
823 if (angle) {
824 threshold -= normal[k];
825 *angle += fabs_(threshold);
826 }
827 }
828 if (qh upper_threshold[k] < REALmax/2) {
829 threshold= qh upper_threshold[k];
830 if (normal[k] > threshold)
831 within= False;
832 if (angle) {
833 threshold -= normal[k];
834 *angle += fabs_(threshold);
835 }
836 }
837 }
838 return within;
839 } /* inthresholds */
840
841
842 /*-<a href="qh-geom.htm#TOC"
843 >-------------------------------</a><a name="joggleinput">-</a>
844
845 qh_joggleinput()
846 randomly joggle input to Qhull by qh.JOGGLEmax
847 initial input is qh.first_point/qh.num_points of qh.hull_dim
848 repeated calls use qh.input_points/qh.num_points
849
850 returns:
851 joggles points at qh.first_point/qh.num_points
852 copies data to qh.input_points/qh.input_malloc if first time
853 determines qh.JOGGLEmax if it was zero
854 if qh.DELAUNAY
855 computes the Delaunay projection of the joggled points
856
857 notes:
858 if qh.DELAUNAY, unnecessarily joggles the last coordinate
859 the initial 'QJn' may be set larger than qh_JOGGLEmaxincrease
860
861 design:
862 if qh.DELAUNAY
863 set qh.SCALElast for reduced precision errors
864 if first call
865 initialize qh.input_points to the original input points
866 if qh.JOGGLEmax == 0
867 determine default qh.JOGGLEmax
868 else
869 increase qh.JOGGLEmax according to qh.build_cnt
870 joggle the input by adding a random number in [-qh.JOGGLEmax,qh.JOGGLEmax]
871 if qh.DELAUNAY
872 sets the Delaunay projection
873 */
qh_joggleinput(void)874 void qh_joggleinput(void) {
875 int i, seed, size;
876 coordT *coordp, *inputp;
877 realT randr, randa, randb;
878
879 if (!qh input_points) { /* first call */
880 qh input_points= qh first_point;
881 qh input_malloc= qh POINTSmalloc;
882 size= qh num_points * qh hull_dim * sizeof(coordT);
883 if (!(qh first_point=(coordT*)qh_malloc((size_t)size))) {
884 qh_fprintf(qh ferr, 6009, "qhull error: insufficient memory to joggle %d points\n",
885 qh num_points);
886 qh_errexit(qh_ERRmem, NULL, NULL);
887 }
888 qh POINTSmalloc= True;
889 if (qh JOGGLEmax == 0.0) {
890 qh JOGGLEmax= qh_detjoggle(qh input_points, qh num_points, qh hull_dim);
891 qh_option("QJoggle", NULL, &qh JOGGLEmax);
892 }
893 }else { /* repeated call */
894 if (!qh RERUN && qh build_cnt > qh_JOGGLEretry) {
895 if (((qh build_cnt-qh_JOGGLEretry-1) % qh_JOGGLEagain) == 0) {
896 realT maxjoggle= qh MAXwidth * qh_JOGGLEmaxincrease;
897 if (qh JOGGLEmax < maxjoggle) {
898 qh JOGGLEmax *= qh_JOGGLEincrease;
899 minimize_(qh JOGGLEmax, maxjoggle);
900 }
901 }
902 }
903 qh_option("QJoggle", NULL, &qh JOGGLEmax);
904 }
905 if (qh build_cnt > 1 && qh JOGGLEmax > fmax_(qh MAXwidth/4, 0.1)) {
906 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",
907 qh JOGGLEmax);
908 qh_errexit(qh_ERRqhull, NULL, NULL);
909 }
910 /* for some reason, using qh ROTATErandom and qh_RANDOMseed does not repeat the run. Use 'TRn' instead */
911 seed= qh_RANDOMint;
912 qh_option("_joggle-seed", &seed, NULL);
913 trace0((qh ferr, 6, "qh_joggleinput: joggle input by %2.2g with seed %d\n",
914 qh JOGGLEmax, seed));
915 inputp= qh input_points;
916 coordp= qh first_point;
917 randa= 2.0 * qh JOGGLEmax/qh_RANDOMmax;
918 randb= -qh JOGGLEmax;
919 size= qh num_points * qh hull_dim;
920 for (i=size; i--; ) {
921 randr= qh_RANDOMint;
922 *(coordp++)= *(inputp++) + (randr * randa + randb);
923 }
924 if (qh DELAUNAY) {
925 qh last_low= qh last_high= qh last_newhigh= REALmax;
926 qh_setdelaunay(qh hull_dim, qh num_points, qh first_point);
927 }
928 } /* joggleinput */
929
930 /*-<a href="qh-geom.htm#TOC"
931 >-------------------------------</a><a name="maxabsval">-</a>
932
933 qh_maxabsval( normal, dim )
934 return pointer to maximum absolute value of a dim vector
935 returns NULL if dim=0
936 */
qh_maxabsval(realT * normal,int dim)937 realT *qh_maxabsval(realT *normal, int dim) {
938 realT maxval= -REALmax;
939 realT *maxp= NULL, *colp, absval;
940 int k;
941
942 for (k=dim, colp= normal; k--; colp++) {
943 absval= fabs_(*colp);
944 if (absval > maxval) {
945 maxval= absval;
946 maxp= colp;
947 }
948 }
949 return maxp;
950 } /* maxabsval */
951
952
953 /*-<a href="qh-geom.htm#TOC"
954 >-------------------------------</a><a name="maxmin">-</a>
955
956 qh_maxmin( points, numpoints, dimension )
957 return max/min points for each dimension
958 determine max and min coordinates
959
960 returns:
961 returns a temporary set of max and min points
962 may include duplicate points. Does not include qh.GOODpoint
963 sets qh.NEARzero, qh.MAXabs_coord, qh.MAXsumcoord, qh.MAXwidth
964 qh.MAXlastcoord, qh.MINlastcoord
965 initializes qh.max_outside, qh.min_vertex, qh.WAScoplanar, qh.ZEROall_ok
966
967 notes:
968 loop duplicated in qh_detjoggle()
969
970 design:
971 initialize global precision variables
972 checks definition of REAL...
973 for each dimension
974 for each point
975 collect maximum and minimum point
976 collect maximum of maximums and minimum of minimums
977 determine qh.NEARzero for Gaussian Elimination
978 */
qh_maxmin(pointT * points,int numpoints,int dimension)979 setT *qh_maxmin(pointT *points, int numpoints, int dimension) {
980 int k;
981 realT maxcoord, temp;
982 pointT *minimum, *maximum, *point, *pointtemp;
983 setT *set;
984
985 qh max_outside= 0.0;
986 qh MAXabs_coord= 0.0;
987 qh MAXwidth= -REALmax;
988 qh MAXsumcoord= 0.0;
989 qh min_vertex= 0.0;
990 qh WAScoplanar= False;
991 if (qh ZEROcentrum)
992 qh ZEROall_ok= True;
993 if (REALmin < REALepsilon && REALmin < REALmax && REALmin > -REALmax
994 && REALmax > 0.0 && -REALmax < 0.0)
995 ; /* all ok */
996 else {
997 qh_fprintf(qh ferr, 6011, "qhull error: floating point constants in user.h are wrong\n\
998 REALepsilon %g REALmin %g REALmax %g -REALmax %g\n",
999 REALepsilon, REALmin, REALmax, -REALmax);
1000 qh_errexit(qh_ERRinput, NULL, NULL);
1001 }
1002 set= qh_settemp(2*dimension);
1003 for (k=0; k < dimension; k++) {
1004 if (points == qh GOODpointp)
1005 minimum= maximum= points + dimension;
1006 else
1007 minimum= maximum= points;
1008 FORALLpoint_(points, numpoints) {
1009 if (point == qh GOODpointp)
1010 continue;
1011 if (maximum[k] < point[k])
1012 maximum= point;
1013 else if (minimum[k] > point[k])
1014 minimum= point;
1015 }
1016 if (k == dimension-1) {
1017 qh MINlastcoord= minimum[k];
1018 qh MAXlastcoord= maximum[k];
1019 }
1020 if (qh SCALElast && k == dimension-1)
1021 maxcoord= qh MAXwidth;
1022 else {
1023 maxcoord= fmax_(maximum[k], -minimum[k]);
1024 if (qh GOODpointp) {
1025 temp= fmax_(qh GOODpointp[k], -qh GOODpointp[k]);
1026 maximize_(maxcoord, temp);
1027 }
1028 temp= maximum[k] - minimum[k];
1029 maximize_(qh MAXwidth, temp);
1030 }
1031 maximize_(qh MAXabs_coord, maxcoord);
1032 qh MAXsumcoord += maxcoord;
1033 qh_setappend(&set, maximum);
1034 qh_setappend(&set, minimum);
1035 /* calculation of qh NEARzero is based on Golub & van Loan, 1983,
1036 Eq. 4.4-13 for "Gaussian elimination with complete pivoting".
1037 Golub & van Loan say that 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 for (k= (dim > 0 ? dim : -dim); k--; ) {
1326 diff= *point1++ - *point2++;
1327 dist += diff * diff;
1328 }
1329 if (dim > 0)
1330 return(sqrt(dist));
1331 return dist;
1332 } /* pointdist */
1333
1334
1335 /*-<a href="qh-geom.htm#TOC"
1336 >-------------------------------</a><a name="printmatrix">-</a>
1337
1338 qh_printmatrix( fp, string, rows, numrow, numcol )
1339 print matrix to fp given by row vectors
1340 print string as header
1341
1342 notes:
1343 print a vector by qh_printmatrix(fp, "", &vect, 1, len)
1344 */
qh_printmatrix(FILE * fp,const char * string,realT ** rows,int numrow,int numcol)1345 void qh_printmatrix(FILE *fp, const char *string, realT **rows, int numrow, int numcol) {
1346 realT *rowp;
1347 realT r; /*bug fix*/
1348 int i,k;
1349
1350 qh_fprintf(fp, 9001, "%s\n", string);
1351 for (i=0; i < numrow; i++) {
1352 rowp= rows[i];
1353 for (k=0; k < numcol; k++) {
1354 r= *rowp++;
1355 qh_fprintf(fp, 9002, "%6.3g ", r);
1356 }
1357 qh_fprintf(fp, 9003, "\n");
1358 }
1359 } /* printmatrix */
1360
1361
1362 /*-<a href="qh-geom.htm#TOC"
1363 >-------------------------------</a><a name="printpoints">-</a>
1364
1365 qh_printpoints( fp, string, points )
1366 print pointids to fp for a set of points
1367 if string, prints string and 'p' point ids
1368 */
qh_printpoints(FILE * fp,const char * string,setT * points)1369 void qh_printpoints(FILE *fp, const char *string, setT *points) {
1370 pointT *point, **pointp;
1371
1372 if (string) {
1373 qh_fprintf(fp, 9004, "%s", string);
1374 FOREACHpoint_(points)
1375 qh_fprintf(fp, 9005, " p%d", qh_pointid(point));
1376 qh_fprintf(fp, 9006, "\n");
1377 }else {
1378 FOREACHpoint_(points)
1379 qh_fprintf(fp, 9007, " %d", qh_pointid(point));
1380 qh_fprintf(fp, 9008, "\n");
1381 }
1382 } /* printpoints */
1383
1384
1385 /*-<a href="qh-geom.htm#TOC"
1386 >-------------------------------</a><a name="projectinput">-</a>
1387
1388 qh_projectinput()
1389 project input points using qh.lower_bound/upper_bound and qh DELAUNAY
1390 if qh.lower_bound[k]=qh.upper_bound[k]= 0,
1391 removes dimension k
1392 if halfspace intersection
1393 removes dimension k from qh.feasible_point
1394 input points in qh first_point, num_points, input_dim
1395
1396 returns:
1397 new point array in qh first_point of qh hull_dim coordinates
1398 sets qh POINTSmalloc
1399 if qh DELAUNAY
1400 projects points to paraboloid
1401 lowbound/highbound is also projected
1402 if qh ATinfinity
1403 adds point "at-infinity"
1404 if qh POINTSmalloc
1405 frees old point array
1406
1407 notes:
1408 checks that qh.hull_dim agrees with qh.input_dim, PROJECTinput, and DELAUNAY
1409
1410
1411 design:
1412 sets project[k] to -1 (delete), 0 (keep), 1 (add for Delaunay)
1413 determines newdim and newnum for qh hull_dim and qh num_points
1414 projects points to newpoints
1415 projects qh.lower_bound to itself
1416 projects qh.upper_bound to itself
1417 if qh DELAUNAY
1418 if qh ATINFINITY
1419 projects points to paraboloid
1420 computes "infinity" point as vertex average and 10% above all points
1421 else
1422 uses qh_setdelaunay to project points to paraboloid
1423 */
qh_projectinput(void)1424 void qh_projectinput(void) {
1425 int k,i;
1426 int newdim= qh input_dim, newnum= qh num_points;
1427 signed char *project;
1428 int projectsize= (qh input_dim+1)*sizeof(*project);
1429 pointT *newpoints, *coord, *infinity;
1430 realT paraboloid, maxboloid= 0;
1431
1432 project= (signed char*)qh_memalloc(projectsize);
1433 memset((char*)project, 0, (size_t)projectsize);
1434 for (k=0; k < qh input_dim; k++) { /* skip Delaunay bound */
1435 if (qh lower_bound[k] == 0 && qh upper_bound[k] == 0) {
1436 project[k]= -1;
1437 newdim--;
1438 }
1439 }
1440 if (qh DELAUNAY) {
1441 project[k]= 1;
1442 newdim++;
1443 if (qh ATinfinity)
1444 newnum++;
1445 }
1446 if (newdim != qh hull_dim) {
1447 qh_memfree(project, projectsize);
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= qh temp_malloc= (coordT*)qh_malloc(newnum*newdim*sizeof(coordT)))){
1452 qh_memfree(project, projectsize);
1453 qh_fprintf(qh ferr, 6016, "qhull error: insufficient memory to project %d points\n",
1454 qh num_points);
1455 qh_errexit(qh_ERRmem, NULL, NULL);
1456 }
1457 /* qh_projectpoints throws error if mismatched dimensions */
1458 qh_projectpoints(project, qh input_dim+1, qh first_point,
1459 qh num_points, qh input_dim, newpoints, newdim);
1460 trace1((qh ferr, 1003, "qh_projectinput: updating lower and upper_bound\n"));
1461 qh_projectpoints(project, qh input_dim+1, qh lower_bound,
1462 1, qh input_dim+1, qh lower_bound, newdim+1);
1463 qh_projectpoints(project, qh input_dim+1, qh upper_bound,
1464 1, qh input_dim+1, qh upper_bound, newdim+1);
1465 if (qh HALFspace) {
1466 if (!qh feasible_point) {
1467 qh_memfree(project, projectsize);
1468 qh_fprintf(qh ferr, 6017, "qhull internal error (qh_projectinput): HALFspace defined without qh.feasible_point\n");
1469 qh_errexit(qh_ERRqhull, NULL, NULL);
1470 }
1471 qh_projectpoints(project, qh input_dim, qh feasible_point,
1472 1, qh input_dim, qh feasible_point, newdim);
1473 }
1474 qh_memfree(project, projectsize);
1475 if (qh POINTSmalloc)
1476 qh_free(qh first_point);
1477 qh first_point= newpoints;
1478 qh POINTSmalloc= True;
1479 qh temp_malloc= NULL;
1480 if (qh DELAUNAY && qh ATinfinity) {
1481 coord= qh first_point;
1482 infinity= qh first_point + qh hull_dim * qh num_points;
1483 for (k=qh hull_dim-1; k--; )
1484 infinity[k]= 0.0;
1485 for (i=qh num_points; i--; ) {
1486 paraboloid= 0.0;
1487 for (k=0; k < qh hull_dim-1; k++) {
1488 paraboloid += *coord * *coord;
1489 infinity[k] += *coord;
1490 coord++;
1491 }
1492 *(coord++)= paraboloid;
1493 maximize_(maxboloid, paraboloid);
1494 }
1495 /* coord == infinity */
1496 for (k=qh hull_dim-1; k--; )
1497 *(coord++) /= qh num_points;
1498 *(coord++)= maxboloid * 1.1;
1499 qh num_points++;
1500 trace0((qh ferr, 9, "qh_projectinput: projected points to paraboloid for Delaunay\n"));
1501 }else if (qh DELAUNAY) /* !qh ATinfinity */
1502 qh_setdelaunay( qh hull_dim, qh num_points, qh first_point);
1503 } /* projectinput */
1504
1505
1506 /*-<a href="qh-geom.htm#TOC"
1507 >-------------------------------</a><a name="projectpoints">-</a>
1508
1509 qh_projectpoints( project, n, points, numpoints, dim, newpoints, newdim )
1510 project points/numpoints/dim to newpoints/newdim
1511 if project[k] == -1
1512 delete dimension k
1513 if project[k] == 1
1514 add dimension k by duplicating previous column
1515 n is size of project
1516
1517 notes:
1518 newpoints may be points if only adding dimension at end
1519
1520 design:
1521 check that 'project' and 'newdim' agree
1522 for each dimension
1523 if project == -1
1524 skip dimension
1525 else
1526 determine start of column in newpoints
1527 determine start of column in points
1528 if project == +1, duplicate previous column
1529 copy dimension (column) from points to newpoints
1530 */
qh_projectpoints(signed char * project,int n,realT * points,int numpoints,int dim,realT * newpoints,int newdim)1531 void qh_projectpoints(signed char *project, int n, realT *points,
1532 int numpoints, int dim, realT *newpoints, int newdim) {
1533 int testdim= dim, oldk=0, newk=0, i,j=0,k;
1534 realT *newp, *oldp;
1535
1536 for (k=0; k < n; k++)
1537 testdim += project[k];
1538 if (testdim != newdim) {
1539 qh_fprintf(qh ferr, 6018, "qhull internal error (qh_projectpoints): newdim %d should be %d after projection\n",
1540 newdim, testdim);
1541 qh_errexit(qh_ERRqhull, NULL, NULL);
1542 }
1543 for (j=0; j<n; j++) {
1544 if (project[j] == -1)
1545 oldk++;
1546 else {
1547 newp= newpoints+newk++;
1548 if (project[j] == +1) {
1549 if (oldk >= dim)
1550 continue;
1551 oldp= points+oldk;
1552 }else
1553 oldp= points+oldk++;
1554 for (i=numpoints; i--; ) {
1555 *newp= *oldp;
1556 newp += newdim;
1557 oldp += dim;
1558 }
1559 }
1560 if (oldk >= dim)
1561 break;
1562 }
1563 trace1((qh ferr, 1004, "qh_projectpoints: projected %d points from dim %d to dim %d\n",
1564 numpoints, dim, newdim));
1565 } /* projectpoints */
1566
1567
1568 /*-<a href="qh-geom.htm#TOC"
1569 >-------------------------------</a><a name="rotateinput">-</a>
1570
1571 qh_rotateinput( rows )
1572 rotate input using row matrix
1573 input points given by qh first_point, num_points, hull_dim
1574 assumes rows[dim] is a scratch buffer
1575 if qh POINTSmalloc, overwrites input points, else mallocs a new array
1576
1577 returns:
1578 rotated input
1579 sets qh POINTSmalloc
1580
1581 design:
1582 see qh_rotatepoints
1583 */
qh_rotateinput(realT ** rows)1584 void qh_rotateinput(realT **rows) {
1585
1586 if (!qh POINTSmalloc) {
1587 qh first_point= qh_copypoints(qh first_point, qh num_points, qh hull_dim);
1588 qh POINTSmalloc= True;
1589 }
1590 qh_rotatepoints(qh first_point, qh num_points, qh hull_dim, rows);
1591 } /* rotateinput */
1592
1593 /*-<a href="qh-geom.htm#TOC"
1594 >-------------------------------</a><a name="rotatepoints">-</a>
1595
1596 qh_rotatepoints( points, numpoints, dim, row )
1597 rotate numpoints points by a d-dim row matrix
1598 assumes rows[dim] is a scratch buffer
1599
1600 returns:
1601 rotated points in place
1602
1603 design:
1604 for each point
1605 for each coordinate
1606 use row[dim] to compute partial inner product
1607 for each coordinate
1608 rotate by partial inner product
1609 */
qh_rotatepoints(realT * points,int numpoints,int dim,realT ** row)1610 void qh_rotatepoints(realT *points, int numpoints, int dim, realT **row) {
1611 realT *point, *rowi, *coord= NULL, sum, *newval;
1612 int i,j,k;
1613
1614 if (qh IStracing >= 1)
1615 qh_printmatrix(qh ferr, "qh_rotatepoints: rotate points by", row, dim, dim);
1616 for (point= points, j= numpoints; j--; point += dim) {
1617 newval= row[dim];
1618 for (i=0; i < dim; i++) {
1619 rowi= row[i];
1620 coord= point;
1621 for (sum= 0.0, k= dim; k--; )
1622 sum += *rowi++ * *coord++;
1623 *(newval++)= sum;
1624 }
1625 for (k=dim; k--; )
1626 *(--coord)= *(--newval);
1627 }
1628 } /* rotatepoints */
1629
1630
1631 /*-<a href="qh-geom.htm#TOC"
1632 >-------------------------------</a><a name="scaleinput">-</a>
1633
1634 qh_scaleinput()
1635 scale input points using qh low_bound/high_bound
1636 input points given by qh first_point, num_points, hull_dim
1637 if qh POINTSmalloc, overwrites input points, else mallocs a new array
1638
1639 returns:
1640 scales coordinates of points to low_bound[k], high_bound[k]
1641 sets qh POINTSmalloc
1642
1643 design:
1644 see qh_scalepoints
1645 */
qh_scaleinput(void)1646 void qh_scaleinput(void) {
1647
1648 if (!qh POINTSmalloc) {
1649 qh first_point= qh_copypoints(qh first_point, qh num_points, qh hull_dim);
1650 qh POINTSmalloc= True;
1651 }
1652 qh_scalepoints(qh first_point, qh num_points, qh hull_dim,
1653 qh lower_bound, qh upper_bound);
1654 } /* scaleinput */
1655
1656 /*-<a href="qh-geom.htm#TOC"
1657 >-------------------------------</a><a name="scalelast">-</a>
1658
1659 qh_scalelast( points, numpoints, dim, low, high, newhigh )
1660 scale last coordinate to [0,m] for Delaunay triangulations
1661 input points given by points, numpoints, dim
1662
1663 returns:
1664 changes scale of last coordinate from [low, high] to [0, newhigh]
1665 overwrites last coordinate of each point
1666 saves low/high/newhigh in qh.last_low, etc. for qh_setdelaunay()
1667
1668 notes:
1669 when called by qh_setdelaunay, low/high may not match actual data
1670
1671 design:
1672 compute scale and shift factors
1673 apply to last coordinate of each point
1674 */
qh_scalelast(coordT * points,int numpoints,int dim,coordT low,coordT high,coordT newhigh)1675 void qh_scalelast(coordT *points, int numpoints, int dim, coordT low,
1676 coordT high, coordT newhigh) {
1677 realT scale, shift;
1678 coordT *coord;
1679 int i;
1680 boolT nearzero= False;
1681
1682 trace4((qh ferr, 4013, "qh_scalelast: scale last coordinate from [%2.2g, %2.2g] to [0,%2.2g]\n",
1683 low, high, newhigh));
1684 qh last_low= low;
1685 qh last_high= high;
1686 qh last_newhigh= newhigh;
1687 scale= qh_divzero(newhigh, high - low,
1688 qh MINdenom_1, &nearzero);
1689 if (nearzero) {
1690 if (qh DELAUNAY)
1691 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");
1692 else
1693 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",
1694 newhigh, low, high, high-low);
1695 qh_errexit(qh_ERRinput, NULL, NULL);
1696 }
1697 shift= - low * newhigh / (high-low);
1698 coord= points + dim - 1;
1699 for (i=numpoints; i--; coord += dim)
1700 *coord= *coord * scale + shift;
1701 } /* scalelast */
1702
1703 /*-<a href="qh-geom.htm#TOC"
1704 >-------------------------------</a><a name="scalepoints">-</a>
1705
1706 qh_scalepoints( points, numpoints, dim, newlows, newhighs )
1707 scale points to new lowbound and highbound
1708 retains old bound when newlow= -REALmax or newhigh= +REALmax
1709
1710 returns:
1711 scaled points
1712 overwrites old points
1713
1714 design:
1715 for each coordinate
1716 compute current low and high bound
1717 compute scale and shift factors
1718 scale all points
1719 enforce new low and high bound for all points
1720 */
qh_scalepoints(pointT * points,int numpoints,int dim,realT * newlows,realT * newhighs)1721 void qh_scalepoints(pointT *points, int numpoints, int dim,
1722 realT *newlows, realT *newhighs) {
1723 int i,k;
1724 realT shift, scale, *coord, low, high, newlow, newhigh, mincoord, maxcoord;
1725 boolT nearzero= False;
1726
1727 for (k=0; k < dim; k++) {
1728 newhigh= newhighs[k];
1729 newlow= newlows[k];
1730 if (newhigh > REALmax/2 && newlow < -REALmax/2)
1731 continue;
1732 low= REALmax;
1733 high= -REALmax;
1734 for (i=numpoints, coord=points+k; i--; coord += dim) {
1735 minimize_(low, *coord);
1736 maximize_(high, *coord);
1737 }
1738 if (newhigh > REALmax/2)
1739 newhigh= high;
1740 if (newlow < -REALmax/2)
1741 newlow= low;
1742 if (qh DELAUNAY && k == dim-1 && newhigh < newlow) {
1743 qh_fprintf(qh ferr, 6021, "qhull input error: 'Qb%d' or 'QB%d' inverts paraboloid since high bound %.2g < low bound %.2g\n",
1744 k, k, newhigh, newlow);
1745 qh_errexit(qh_ERRinput, NULL, NULL);
1746 }
1747 scale= qh_divzero(newhigh - newlow, high - low,
1748 qh MINdenom_1, &nearzero);
1749 if (nearzero) {
1750 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",
1751 k, newlow, newhigh, low, high);
1752 qh_errexit(qh_ERRinput, NULL, NULL);
1753 }
1754 shift= (newlow * high - low * newhigh)/(high-low);
1755 coord= points+k;
1756 for (i=numpoints; i--; coord += dim)
1757 *coord= *coord * scale + shift;
1758 coord= points+k;
1759 if (newlow < newhigh) {
1760 mincoord= newlow;
1761 maxcoord= newhigh;
1762 }else {
1763 mincoord= newhigh;
1764 maxcoord= newlow;
1765 }
1766 for (i=numpoints; i--; coord += dim) {
1767 minimize_(*coord, maxcoord); /* because of roundoff error */
1768 maximize_(*coord, mincoord);
1769 }
1770 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",
1771 k, low, high, newlow, newhigh, numpoints, scale, shift));
1772 }
1773 } /* scalepoints */
1774
1775
1776 /*-<a href="qh-geom.htm#TOC"
1777 >-------------------------------</a><a name="setdelaunay">-</a>
1778
1779 qh_setdelaunay( dim, count, points )
1780 project count points to dim-d paraboloid for Delaunay triangulation
1781
1782 dim is one more than the dimension of the input set
1783 assumes dim is at least 3 (i.e., at least a 2-d Delaunay triangulation)
1784
1785 points is a dim*count realT array. The first dim-1 coordinates
1786 are the coordinates of the first input point. array[dim] is
1787 the first coordinate of the second input point. array[2*dim] is
1788 the first coordinate of the third input point.
1789
1790 if qh.last_low defined (i.e., 'Qbb' called qh_scalelast)
1791 calls qh_scalelast to scale the last coordinate the same as the other points
1792
1793 returns:
1794 for each point
1795 sets point[dim-1] to sum of squares of coordinates
1796 scale points to 'Qbb' if needed
1797
1798 notes:
1799 to project one point, use
1800 qh_setdelaunay(qh hull_dim, 1, point)
1801
1802 Do not use options 'Qbk', 'QBk', or 'QbB' since they scale
1803 the coordinates after the original projection.
1804
1805 */
qh_setdelaunay(int dim,int count,pointT * points)1806 void qh_setdelaunay(int dim, int count, pointT *points) {
1807 int i, k;
1808 coordT *coordp, coord;
1809 realT paraboloid;
1810
1811 trace0((qh ferr, 11, "qh_setdelaunay: project %d points to paraboloid for Delaunay triangulation\n", count));
1812 coordp= points;
1813 for (i=0; i < count; i++) {
1814 coord= *coordp++;
1815 paraboloid= coord*coord;
1816 for (k=dim-2; k--; ) {
1817 coord= *coordp++;
1818 paraboloid += coord*coord;
1819 }
1820 *coordp++ = paraboloid;
1821 }
1822 if (qh last_low < REALmax/2)
1823 qh_scalelast(points, count, dim, qh last_low, qh last_high, qh last_newhigh);
1824 } /* setdelaunay */
1825
1826
1827 /*-<a href="qh-geom.htm#TOC"
1828 >-------------------------------</a><a name="sethalfspace">-</a>
1829
1830 qh_sethalfspace( dim, coords, nextp, normal, offset, feasible )
1831 set point to dual of halfspace relative to feasible point
1832 halfspace is normal coefficients and offset.
1833
1834 returns:
1835 false and prints error if feasible point is outside of hull
1836 overwrites coordinates for point at dim coords
1837 nextp= next point (coords)
1838 does not call qh_errexit
1839
1840 design:
1841 compute distance from feasible point to halfspace
1842 divide each normal coefficient by -dist
1843 */
qh_sethalfspace(int dim,coordT * coords,coordT ** nextp,coordT * normal,coordT * offset,coordT * feasible)1844 boolT qh_sethalfspace(int dim, coordT *coords, coordT **nextp,
1845 coordT *normal, coordT *offset, coordT *feasible) {
1846 coordT *normp= normal, *feasiblep= feasible, *coordp= coords;
1847 realT dist;
1848 realT r; /*bug fix*/
1849 int k;
1850 boolT zerodiv;
1851
1852 dist= *offset;
1853 for (k=dim; k--; )
1854 dist += *(normp++) * *(feasiblep++);
1855 if (dist > 0)
1856 goto LABELerroroutside;
1857 normp= normal;
1858 if (dist < -qh MINdenom) {
1859 for (k=dim; k--; )
1860 *(coordp++)= *(normp++) / -dist;
1861 }else {
1862 for (k=dim; k--; ) {
1863 *(coordp++)= qh_divzero(*(normp++), -dist, qh MINdenom_1, &zerodiv);
1864 if (zerodiv)
1865 goto LABELerroroutside;
1866 }
1867 }
1868 *nextp= coordp;
1869 if (qh IStracing >= 4) {
1870 qh_fprintf(qh ferr, 8021, "qh_sethalfspace: halfspace at offset %6.2g to point: ", *offset);
1871 for (k=dim, coordp=coords; k--; ) {
1872 r= *coordp++;
1873 qh_fprintf(qh ferr, 8022, " %6.2g", r);
1874 }
1875 qh_fprintf(qh ferr, 8023, "\n");
1876 }
1877 return True;
1878 LABELerroroutside:
1879 feasiblep= feasible;
1880 normp= normal;
1881 qh_fprintf(qh ferr, 6023, "qhull input error: feasible point is not clearly inside halfspace\nfeasible point: ");
1882 for (k=dim; k--; )
1883 qh_fprintf(qh ferr, 8024, qh_REAL_1, r=*(feasiblep++));
1884 qh_fprintf(qh ferr, 8025, "\n halfspace: ");
1885 for (k=dim; k--; )
1886 qh_fprintf(qh ferr, 8026, qh_REAL_1, r=*(normp++));
1887 qh_fprintf(qh ferr, 8027, "\n at offset: ");
1888 qh_fprintf(qh ferr, 8028, qh_REAL_1, *offset);
1889 qh_fprintf(qh ferr, 8029, " and distance: ");
1890 qh_fprintf(qh ferr, 8030, qh_REAL_1, dist);
1891 qh_fprintf(qh ferr, 8031, "\n");
1892 return False;
1893 } /* sethalfspace */
1894
1895 /*-<a href="qh-geom.htm#TOC"
1896 >-------------------------------</a><a name="sethalfspace_all">-</a>
1897
1898 qh_sethalfspace_all( dim, count, halfspaces, feasible )
1899 generate dual for halfspace intersection with feasible point
1900 array of count halfspaces
1901 each halfspace is normal coefficients followed by offset
1902 the origin is inside the halfspace if the offset is negative
1903 feasible is a point inside all halfspaces (http://www.qhull.org/html/qhalf.htm#notes)
1904
1905 returns:
1906 malloc'd array of count X dim-1 points
1907
1908 notes:
1909 call before qh_init_B or qh_initqhull_globals
1910 free memory when done
1911 unused/untested code: please email bradb@shore.net if this works ok for you
1912 if using option 'Fp', qh->feasible_point must be set (e.g., to 'feasible')
1913 qh->feasible_point is a malloc'd array that is freed by qh_freebuffers.
1914
1915 design:
1916 see qh_sethalfspace
1917 */
qh_sethalfspace_all(int dim,int count,coordT * halfspaces,pointT * feasible)1918 coordT *qh_sethalfspace_all(int dim, int count, coordT *halfspaces, pointT *feasible) {
1919 int i, newdim;
1920 pointT *newpoints;
1921 coordT *coordp, *normalp, *offsetp;
1922
1923 trace0((qh ferr, 12, "qh_sethalfspace_all: compute dual for halfspace intersection\n"));
1924 newdim= dim - 1;
1925 if (!(newpoints=(coordT*)qh_malloc(count*newdim*sizeof(coordT)))){
1926 qh_fprintf(qh ferr, 6024, "qhull error: insufficient memory to compute dual of %d halfspaces\n",
1927 count);
1928 qh_errexit(qh_ERRmem, NULL, NULL);
1929 }
1930 coordp= newpoints;
1931 normalp= halfspaces;
1932 for (i=0; i < count; i++) {
1933 offsetp= normalp + newdim;
1934 if (!qh_sethalfspace(newdim, coordp, &coordp, normalp, offsetp, feasible)) {
1935 qh_free(newpoints); /* feasible is not inside halfspace as reported by qh_sethalfspace */
1936 qh_fprintf(qh ferr, 8032, "The halfspace was at index %d\n", i);
1937 qh_errexit(qh_ERRinput, NULL, NULL);
1938 }
1939 normalp= offsetp + 1;
1940 }
1941 return newpoints;
1942 } /* sethalfspace_all */
1943
1944
1945 /*-<a href="qh-geom.htm#TOC"
1946 >-------------------------------</a><a name="sharpnewfacets">-</a>
1947
1948 qh_sharpnewfacets()
1949
1950 returns:
1951 true if could be an acute angle (facets in different quadrants)
1952
1953 notes:
1954 for qh_findbest
1955
1956 design:
1957 for all facets on qh.newfacet_list
1958 if two facets are in different quadrants
1959 set issharp
1960 */
qh_sharpnewfacets(void)1961 boolT qh_sharpnewfacets(void) {
1962 facetT *facet;
1963 boolT issharp = False;
1964 int *quadrant, k;
1965
1966 quadrant= (int*)qh_memalloc(qh hull_dim * sizeof(int));
1967 FORALLfacet_(qh newfacet_list) {
1968 if (facet == qh newfacet_list) {
1969 for (k=qh hull_dim; k--; )
1970 quadrant[ k]= (facet->normal[ k] > 0);
1971 }else {
1972 for (k=qh hull_dim; k--; ) {
1973 if (quadrant[ k] != (facet->normal[ k] > 0)) {
1974 issharp= True;
1975 break;
1976 }
1977 }
1978 }
1979 if (issharp)
1980 break;
1981 }
1982 qh_memfree( quadrant, qh hull_dim * sizeof(int));
1983 trace3((qh ferr, 3001, "qh_sharpnewfacets: %d\n", issharp));
1984 return issharp;
1985 } /* sharpnewfacets */
1986
1987 /*-<a href="qh-geom.htm#TOC"
1988 >-------------------------------</a><a name="voronoi_center">-</a>
1989
1990 qh_voronoi_center( dim, points )
1991 return Voronoi center for a set of points
1992 dim is the orginal dimension of the points
1993 gh.gm_matrix/qh.gm_row are scratch buffers
1994
1995 returns:
1996 center as a temporary point (qh_memalloc)
1997 if non-simplicial,
1998 returns center for max simplex of points
1999
2000 notes:
2001 only called by qh_facetcenter
2002 from Bowyer & Woodwark, A Programmer's Geometry, 1983, p. 65
2003
2004 design:
2005 if non-simplicial
2006 determine max simplex for points
2007 translate point0 of simplex to origin
2008 compute sum of squares of diagonal
2009 compute determinate
2010 compute Voronoi center (see Bowyer & Woodwark)
2011 */
qh_voronoi_center(int dim,setT * points)2012 pointT *qh_voronoi_center(int dim, setT *points) {
2013 pointT *point, **pointp, *point0;
2014 pointT *center= (pointT*)qh_memalloc(qh center_size);
2015 setT *simplex;
2016 int i, j, k, size= qh_setsize(points);
2017 coordT *gmcoord;
2018 realT *diffp, sum2, *sum2row, *sum2p, det, factor;
2019 boolT nearzero, infinite;
2020
2021 if (size == dim+1)
2022 simplex= points;
2023 else if (size < dim+1) {
2024 qh_memfree(center, qh center_size);
2025 qh_fprintf(qh ferr, 6025, "qhull internal error (qh_voronoi_center):\n need at least %d points to construct a Voronoi center\n",
2026 dim+1);
2027 qh_errexit(qh_ERRqhull, NULL, NULL);
2028 simplex= points; /* never executed -- avoids warning */
2029 }else {
2030 simplex= qh_settemp(dim+1);
2031 qh_maxsimplex(dim, points, NULL, 0, &simplex);
2032 }
2033 point0= SETfirstt_(simplex, pointT);
2034 gmcoord= qh gm_matrix;
2035 for (k=0; k < dim; k++) {
2036 qh gm_row[k]= gmcoord;
2037 FOREACHpoint_(simplex) {
2038 if (point != point0)
2039 *(gmcoord++)= point[k] - point0[k];
2040 }
2041 }
2042 sum2row= gmcoord;
2043 for (i=0; i < dim; i++) {
2044 sum2= 0.0;
2045 for (k=0; k < dim; k++) {
2046 diffp= qh gm_row[k] + i;
2047 sum2 += *diffp * *diffp;
2048 }
2049 *(gmcoord++)= sum2;
2050 }
2051 det= qh_determinant(qh gm_row, dim, &nearzero);
2052 factor= qh_divzero(0.5, det, qh MINdenom, &infinite);
2053 if (infinite) {
2054 for (k=dim; k--; )
2055 center[k]= qh_INFINITE;
2056 if (qh IStracing)
2057 qh_printpoints(qh ferr, "qh_voronoi_center: at infinity for ", simplex);
2058 }else {
2059 for (i=0; i < dim; i++) {
2060 gmcoord= qh gm_matrix;
2061 sum2p= sum2row;
2062 for (k=0; k < dim; k++) {
2063 qh gm_row[k]= gmcoord;
2064 if (k == i) {
2065 for (j=dim; j--; )
2066 *(gmcoord++)= *sum2p++;
2067 }else {
2068 FOREACHpoint_(simplex) {
2069 if (point != point0)
2070 *(gmcoord++)= point[k] - point0[k];
2071 }
2072 }
2073 }
2074 center[i]= qh_determinant(qh gm_row, dim, &nearzero)*factor + point0[i];
2075 }
2076 #ifndef qh_NOtrace
2077 if (qh IStracing >= 3) {
2078 qh_fprintf(qh ferr, 8033, "qh_voronoi_center: det %2.2g factor %2.2g ", det, factor);
2079 qh_printmatrix(qh ferr, "center:", ¢er, 1, dim);
2080 if (qh IStracing >= 5) {
2081 qh_printpoints(qh ferr, "points", simplex);
2082 FOREACHpoint_(simplex)
2083 qh_fprintf(qh ferr, 8034, "p%d dist %.2g, ", qh_pointid(point),
2084 qh_pointdist(point, center, dim));
2085 qh_fprintf(qh ferr, 8035, "\n");
2086 }
2087 }
2088 #endif
2089 }
2090 if (simplex != points)
2091 qh_settempfree(&simplex);
2092 return center;
2093 } /* voronoi_center */
2094
2095