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