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