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