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