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