1 /*<html><pre> -<a href="qh-user.htm"
2 >-------------------------------</a><a name="TOP">-</a>
3
4 user.c
5 user redefinable functions
6
7 see user2.c for qh_fprintf, qh_malloc, qh_free
8
9 see README.txt see COPYING.txt for copyright information.
10
11 see libqhull.h for data structures, macros, and user-callable functions.
12
13 see user_eg.c, unix.c, and qhull_interface.cpp for examples.
14
15 see user.h for user-definable constants
16
17 use qh_NOmem in mem.h to turn off memory management
18 use qh_NOmerge in user.h to turn off facet merging
19 set qh_KEEPstatistics in user.h to 0 to turn off statistics
20
21 This is unsupported software. You're welcome to make changes,
22 but you're on your own if something goes wrong. Use 'Tc' to
23 check frequently. Usually qhull will report an error if
24 a data structure becomes inconsistent. If so, it also reports
25 the last point added to the hull, e.g., 102. You can then trace
26 the execution of qhull with "T4P102".
27
28 Please report any errors that you fix to qhull@qhull.org
29
30 call_qhull is a template for calling qhull from within your application
31
32 if you recompile and load this module, then user.o will not be loaded
33 from qhull.a
34
35 you can add additional quick allocation sizes in qh_user_memsizes
36
37 if the other functions here are redefined to not use qh_print...,
38 then io.o will not be loaded from qhull.a. See user_eg.c for an
39 example. We recommend keeping io.o for the extra debugging
40 information it supplies.
41 */
42
43 #include "qhull_a.h"
44
45 #include <stdarg.h>
46
47 /*-<a href="qh-user.htm#TOC"
48 >-------------------------------</a><a name="call_qhull">-</a>
49
50 qh_call_qhull( void )
51 template for calling qhull from inside your program
52 remove #if 0, #endif to compile
53
54 returns:
55 exit code(see qh_ERR... in libqhull.h)
56 all memory freed
57
58 notes:
59 This can be called any number of times.
60
61 see:
62 qh_call_qhull_once()
63
64 */
65 #if 0
66 {
67 int dim; /* dimension of points */
68 int numpoints; /* number of points */
69 coordT *points; /* array of coordinates for each point */
70 boolT ismalloc; /* True if qhull should free points in qh_freeqhull() or reallocation */
71 char flags[]= "qhull Tv"; /* option flags for qhull, see qh_opt.htm */
72 FILE *outfile= stdout; /* output from qh_produce_output()
73 use NULL to skip qh_produce_output() */
74 FILE *errfile= stderr; /* error messages from qhull code */
75 int exitcode; /* 0 if no error from qhull */
76 facetT *facet; /* set by FORALLfacets */
77 int curlong, totlong; /* memory remaining after qh_memfreeshort */
78
79 /* initialize dim, numpoints, points[], ismalloc here */
80 exitcode= qh_new_qhull(dim, numpoints, points, ismalloc,
81 flags, outfile, errfile);
82 if (!exitcode) { /* if no error */
83 /* 'qh facet_list' contains the convex hull */
84 FORALLfacets {
85 /* ... your code ... */
86 }
87 }
88 qh_freeqhull(!qh_ALL);
89 qh_memfreeshort(&curlong, &totlong);
90 if (curlong || totlong)
91 qh_fprintf(errfile, 7068, "qhull internal warning (main): did not free %d bytes of long memory(%d pieces)\n", totlong, curlong);
92 }
93 #endif
94
95 /*-<a href="qh-user.htm#TOC"
96 >-------------------------------</a><a name="new_qhull">-</a>
97
98 qh_new_qhull( dim, numpoints, points, ismalloc, qhull_cmd, outfile, errfile )
99 build new qhull data structure and return exitcode (0 if no errors)
100
101 notes:
102 do not modify points until finished with results.
103 The qhull data structure contains pointers into the points array.
104 do not call qhull functions before qh_new_qhull().
105 The qhull data structure is not initialized until qh_new_qhull().
106
107 outfile may be null
108 qhull_cmd must start with "qhull "
109 projects points to a new point array for Delaunay triangulations ('d' and 'v')
110 transforms points into a new point array for halfspace intersection ('H')
111
112
113 To allow multiple, concurrent calls to qhull()
114 - set qh_QHpointer in user.h
115 - use qh_save_qhull and qh_restore_qhull to swap the global data structure between calls.
116 - use qh_freeqhull(qh_ALL) to free intermediate convex hulls
117
118 see:
119 user_eg.c for an example
120 */
qh_new_qhull(int dim,int numpoints,coordT * points,boolT ismalloc,char * qhull_cmd,FILE * outfile,FILE * errfile)121 int qh_new_qhull(int dim, int numpoints, coordT *points, boolT ismalloc,
122 char *qhull_cmd, FILE *outfile, FILE *errfile) {
123 int exitcode, hulldim;
124 boolT new_ismalloc;
125 static boolT firstcall = True;
126 coordT *new_points;
127
128 if (firstcall) {
129 qh_meminit(errfile);
130 firstcall= False;
131 }
132 if (strncmp(qhull_cmd,"qhull ", (size_t)6)) {
133 qh_fprintf(errfile, 6186, "qhull error (qh_new_qhull): start qhull_cmd argument with \"qhull \"\n");
134 qh_exit(qh_ERRinput);
135 }
136 qh_initqhull_start(NULL, outfile, errfile);
137 trace1((qh ferr, 1044, "qh_new_qhull: build new Qhull for %d %d-d points with %s\n", numpoints, dim, qhull_cmd));
138 exitcode = setjmp(qh errexit);
139 if (!exitcode)
140 {
141 qh NOerrexit = False;
142 qh_initflags(qhull_cmd);
143 if (qh DELAUNAY)
144 qh PROJECTdelaunay= True;
145 if (qh HALFspace) {
146 /* points is an array of halfspaces,
147 the last coordinate of each halfspace is its offset */
148 hulldim= dim-1;
149 qh_setfeasible(hulldim);
150 new_points= qh_sethalfspace_all(dim, numpoints, points, qh feasible_point);
151 new_ismalloc= True;
152 if (ismalloc)
153 qh_free(points);
154 }else {
155 hulldim= dim;
156 new_points= points;
157 new_ismalloc= ismalloc;
158 }
159 qh_init_B(new_points, numpoints, hulldim, new_ismalloc);
160 qh_qhull();
161 qh_check_output();
162 if (outfile)
163 qh_produce_output();
164 if (qh VERIFYoutput && !qh STOPpoint && !qh STOPcone)
165 qh_check_points();
166 }
167 qh NOerrexit = True;
168 return exitcode;
169 } /* new_qhull */
170
171 /*-<a href="qh-user.htm#TOC"
172 >-------------------------------</a><a name="errexit">-</a>
173
174 qh_errexit( exitcode, facet, ridge )
175 report and exit from an error
176 report facet and ridge if non-NULL
177 reports useful information such as last point processed
178 set qh.FORCEoutput to print neighborhood of facet
179
180 see:
181 qh_errexit2() in libqhull.c for printing 2 facets
182
183 design:
184 check for error within error processing
185 compute qh.hulltime
186 print facet and ridge (if any)
187 report commandString, options, qh.furthest_id
188 print summary and statistics (including precision statistics)
189 if qh_ERRsingular
190 print help text for singular data set
191 exit program via long jump (if defined) or exit()
192 */
qh_errexit(int exitcode,facetT * facet,ridgeT * ridge)193 void qh_errexit(int exitcode, facetT *facet, ridgeT *ridge) {
194
195 if (qh ERREXITcalled) {
196 qh_fprintf(qh ferr, 8126, "\nqhull error while processing previous error. Exit program\n");
197 qh_exit(qh_ERRqhull);
198 }
199 qh ERREXITcalled= True;
200 if (!qh QHULLfinished)
201 qh hulltime= qh_CPUclock - qh hulltime;
202 qh_errprint("ERRONEOUS", facet, NULL, ridge, NULL);
203 qh_fprintf(qh ferr, 8127, "\nWhile executing: %s | %s\n", qh rbox_command, qh qhull_command);
204 qh_fprintf(qh ferr, 8128, "Options selected for Qhull %s:\n%s\n", qh_version, qh qhull_options);
205 if (qh furthest_id >= 0) {
206 qh_fprintf(qh ferr, 8129, "Last point added to hull was p%d.", qh furthest_id);
207 if (zzval_(Ztotmerge))
208 qh_fprintf(qh ferr, 8130, " Last merge was #%d.", zzval_(Ztotmerge));
209 if (qh QHULLfinished)
210 qh_fprintf(qh ferr, 8131, "\nQhull has finished constructing the hull.");
211 else if (qh POSTmerging)
212 qh_fprintf(qh ferr, 8132, "\nQhull has started post-merging.");
213 qh_fprintf(qh ferr, 8133, "\n");
214 }
215 if (qh FORCEoutput && (qh QHULLfinished || (!facet && !ridge)))
216 qh_produce_output();
217 else {
218 if (exitcode != qh_ERRsingular && zzval_(Zsetplane) > qh hull_dim+1) {
219 qh_fprintf(qh ferr, 8134, "\nAt error exit:\n");
220 qh_printsummary(qh ferr);
221 if (qh PRINTstatistics) {
222 qh_collectstatistics();
223 qh_printstatistics(qh ferr, "at error exit");
224 qh_memstatistics(qh ferr);
225 }
226 }
227 if (qh PRINTprecision)
228 qh_printstats(qh ferr, qhstat precision, NULL);
229 }
230 if (!exitcode)
231 exitcode= qh_ERRqhull;
232 else if (exitcode == qh_ERRsingular)
233 qh_printhelp_singular(qh ferr);
234 else if (exitcode == qh_ERRprec && !qh PREmerge)
235 qh_printhelp_degenerate(qh ferr);
236 if (qh NOerrexit) {
237 qh_fprintf(qh ferr, 6187, "qhull error while ending program. Exit program\n");
238 qh_exit(qh_ERRqhull);
239 }
240 qh ERREXITcalled= False;
241 qh NOerrexit= True;
242 longjmp(qh errexit, exitcode);
243 } /* errexit */
244
245
246 /*-<a href="qh-user.htm#TOC"
247 >-------------------------------</a><a name="errprint">-</a>
248
249 qh_errprint( fp, string, atfacet, otherfacet, atridge, atvertex )
250 prints out the information of facets and ridges to fp
251 also prints neighbors and geomview output
252
253 notes:
254 except for string, any parameter may be NULL
255 */
qh_errprint(const char * string,facetT * atfacet,facetT * otherfacet,ridgeT * atridge,vertexT * atvertex)256 void qh_errprint(const char *string, facetT *atfacet, facetT *otherfacet, ridgeT *atridge, vertexT *atvertex) {
257 int i;
258
259 if (atfacet) {
260 qh_fprintf(qh ferr, 8135, "%s FACET:\n", string);
261 qh_printfacet(qh ferr, atfacet);
262 }
263 if (otherfacet) {
264 qh_fprintf(qh ferr, 8136, "%s OTHER FACET:\n", string);
265 qh_printfacet(qh ferr, otherfacet);
266 }
267 if (atridge) {
268 qh_fprintf(qh ferr, 8137, "%s RIDGE:\n", string);
269 qh_printridge(qh ferr, atridge);
270 if (atridge->top && atridge->top != atfacet && atridge->top != otherfacet)
271 qh_printfacet(qh ferr, atridge->top);
272 if (atridge->bottom
273 && atridge->bottom != atfacet && atridge->bottom != otherfacet)
274 qh_printfacet(qh ferr, atridge->bottom);
275 if (!atfacet)
276 atfacet= atridge->top;
277 if (!otherfacet)
278 otherfacet= otherfacet_(atridge, atfacet);
279 }
280 if (atvertex) {
281 qh_fprintf(qh ferr, 8138, "%s VERTEX:\n", string);
282 qh_printvertex(qh ferr, atvertex);
283 }
284 if (qh fout && qh FORCEoutput && atfacet && !qh QHULLfinished && !qh IStracing) {
285 qh_fprintf(qh ferr, 8139, "ERRONEOUS and NEIGHBORING FACETS to output\n");
286 for (i=0; i < qh_PRINTEND; i++) /* use fout for geomview output */
287 qh_printneighborhood(qh fout, qh PRINTout[i], atfacet, otherfacet,
288 !qh_ALL);
289 }
290 } /* errprint */
291
292
293 /*-<a href="qh-user.htm#TOC"
294 >-------------------------------</a><a name="printfacetlist">-</a>
295
296 qh_printfacetlist( fp, facetlist, facets, printall )
297 print all fields for a facet list and/or set of facets to fp
298 if !printall,
299 only prints good facets
300
301 notes:
302 also prints all vertices
303 */
qh_printfacetlist(facetT * facetlist,setT * facets,boolT printall)304 void qh_printfacetlist(facetT *facetlist, setT *facets, boolT printall) {
305 facetT *facet, **facetp;
306
307 qh_printbegin(qh ferr, qh_PRINTfacets, facetlist, facets, printall);
308 FORALLfacet_(facetlist)
309 qh_printafacet(qh ferr, qh_PRINTfacets, facet, printall);
310 FOREACHfacet_(facets)
311 qh_printafacet(qh ferr, qh_PRINTfacets, facet, printall);
312 qh_printend(qh ferr, qh_PRINTfacets, facetlist, facets, printall);
313 } /* printfacetlist */
314
315
316 /*-<a href="qh-io.htm#TOC"
317 >-------------------------------</a><a name="printhelp_degenerate">-</a>
318
319 qh_printhelp_degenerate( fp )
320 prints descriptive message for precision error
321
322 notes:
323 no message if qh_QUICKhelp
324 */
qh_printhelp_degenerate(FILE * fp)325 void qh_printhelp_degenerate(FILE *fp) {
326
327 if (qh MERGEexact || qh PREmerge || qh JOGGLEmax < REALmax/2)
328 qh_fprintf(fp, 9368, "\n\
329 A Qhull error has occurred. Qhull should have corrected the above\n\
330 precision error. Please send the input and all of the output to\n\
331 qhull_bug@qhull.org\n");
332 else if (!qh_QUICKhelp) {
333 qh_fprintf(fp, 9369, "\n\
334 Precision problems were detected during construction of the convex hull.\n\
335 This occurs because convex hull algorithms assume that calculations are\n\
336 exact, but floating-point arithmetic has roundoff errors.\n\
337 \n\
338 To correct for precision problems, do not use 'Q0'. By default, Qhull\n\
339 selects 'C-0' or 'Qx' and merges non-convex facets. With option 'QJ',\n\
340 Qhull joggles the input to prevent precision problems. See \"Imprecision\n\
341 in Qhull\" (qh-impre.htm).\n\
342 \n\
343 If you use 'Q0', the output may include\n\
344 coplanar ridges, concave ridges, and flipped facets. In 4-d and higher,\n\
345 Qhull may produce a ridge with four neighbors or two facets with the same \n\
346 vertices. Qhull reports these events when they occur. It stops when a\n\
347 concave ridge, flipped facet, or duplicate facet occurs.\n");
348 #if REALfloat
349 qh_fprintf(fp, 9370, "\
350 \n\
351 Qhull is currently using single precision arithmetic. The following\n\
352 will probably remove the precision problems:\n\
353 - recompile qhull for realT precision(#define REALfloat 0 in user.h).\n");
354 #endif
355 if (qh DELAUNAY && !qh SCALElast && qh MAXabs_coord > 1e4)
356 qh_fprintf(fp, 9371, "\
357 \n\
358 When computing the Delaunay triangulation of coordinates > 1.0,\n\
359 - use 'Qbb' to scale the last coordinate to [0,m] (max previous coordinate)\n");
360 if (qh DELAUNAY && !qh ATinfinity)
361 qh_fprintf(fp, 9372, "\
362 When computing the Delaunay triangulation:\n\
363 - use 'Qz' to add a point at-infinity. This reduces precision problems.\n");
364
365 qh_fprintf(fp, 9373, "\
366 \n\
367 If you need triangular output:\n\
368 - use option 'Qt' to triangulate the output\n\
369 - use option 'QJ' to joggle the input points and remove precision errors\n\
370 - use option 'Ft'. It triangulates non-simplicial facets with added points.\n\
371 \n\
372 If you must use 'Q0',\n\
373 try one or more of the following options. They can not guarantee an output.\n\
374 - use 'QbB' to scale the input to a cube.\n\
375 - use 'Po' to produce output and prevent partitioning for flipped facets\n\
376 - use 'V0' to set min. distance to visible facet as 0 instead of roundoff\n\
377 - use 'En' to specify a maximum roundoff error less than %2.2g.\n\
378 - options 'Qf', 'Qbb', and 'QR0' may also help\n",
379 qh DISTround);
380 qh_fprintf(fp, 9374, "\
381 \n\
382 To guarantee simplicial output:\n\
383 - use option 'Qt' to triangulate the output\n\
384 - use option 'QJ' to joggle the input points and remove precision errors\n\
385 - use option 'Ft' to triangulate the output by adding points\n\
386 - use exact arithmetic (see \"Imprecision in Qhull\", qh-impre.htm)\n\
387 ");
388 }
389 } /* printhelp_degenerate */
390
391
392 /*-<a href="qh-globa.htm#TOC"
393 >-------------------------------</a><a name="printhelp_narrowhull">-</a>
394
395 qh_printhelp_narrowhull( minangle )
396 Warn about a narrow hull
397
398 notes:
399 Alternatively, reduce qh_WARNnarrow in user.h
400
401 */
qh_printhelp_narrowhull(FILE * fp,realT minangle)402 void qh_printhelp_narrowhull(FILE *fp, realT minangle) {
403
404 qh_fprintf(fp, 9375, "qhull precision warning: \n\
405 The initial hull is narrow (cosine of min. angle is %.16f).\n\
406 A coplanar point may lead to a wide facet. Options 'QbB' (scale to unit box)\n\
407 or 'Qbb' (scale last coordinate) may remove this warning. Use 'Pp' to skip\n\
408 this warning. See 'Limitations' in qh-impre.htm.\n",
409 -minangle); /* convert from angle between normals to angle between facets */
410 } /* printhelp_narrowhull */
411
412 /*-<a href="qh-io.htm#TOC"
413 >-------------------------------</a><a name="printhelp_singular">-</a>
414
415 qh_printhelp_singular( fp )
416 prints descriptive message for singular input
417 */
qh_printhelp_singular(FILE * fp)418 void qh_printhelp_singular(FILE *fp) {
419 facetT *facet;
420 vertexT *vertex, **vertexp;
421 realT min, max, *coord, dist;
422 int i,k;
423
424 qh_fprintf(fp, 9376, "\n\
425 The input to qhull appears to be less than %d dimensional, or a\n\
426 computation has overflowed.\n\n\
427 Qhull could not construct a clearly convex simplex from points:\n",
428 qh hull_dim);
429 qh_printvertexlist(fp, "", qh facet_list, NULL, qh_ALL);
430 if (!qh_QUICKhelp)
431 qh_fprintf(fp, 9377, "\n\
432 The center point is coplanar with a facet, or a vertex is coplanar\n\
433 with a neighboring facet. The maximum round off error for\n\
434 computing distances is %2.2g. The center point, facets and distances\n\
435 to the center point are as follows:\n\n", qh DISTround);
436 qh_printpointid(fp, "center point", qh hull_dim, qh interior_point, -1);
437 qh_fprintf(fp, 9378, "\n");
438 FORALLfacets {
439 qh_fprintf(fp, 9379, "facet");
440 FOREACHvertex_(facet->vertices)
441 qh_fprintf(fp, 9380, " p%d", qh_pointid(vertex->point));
442 zinc_(Zdistio);
443 qh_distplane(qh interior_point, facet, &dist);
444 qh_fprintf(fp, 9381, " distance= %4.2g\n", dist);
445 }
446 if (!qh_QUICKhelp) {
447 if (qh HALFspace)
448 qh_fprintf(fp, 9382, "\n\
449 These points are the dual of the given halfspaces. They indicate that\n\
450 the intersection is degenerate.\n");
451 qh_fprintf(fp, 9383,"\n\
452 These points either have a maximum or minimum x-coordinate, or\n\
453 they maximize the determinant for k coordinates. Trial points\n\
454 are first selected from points that maximize a coordinate.\n");
455 if (qh hull_dim >= qh_INITIALmax)
456 qh_fprintf(fp, 9384, "\n\
457 Because of the high dimension, the min x-coordinate and max-coordinate\n\
458 points are used if the determinant is non-zero. Option 'Qs' will\n\
459 do a better, though much slower, job. Instead of 'Qs', you can change\n\
460 the points by randomly rotating the input with 'QR0'.\n");
461 }
462 qh_fprintf(fp, 9385, "\nThe min and max coordinates for each dimension are:\n");
463 for (k=0; k < qh hull_dim; k++) {
464 min= REALmax;
465 max= -REALmin;
466 for (i=qh num_points, coord= qh first_point+k; i--; coord += qh hull_dim) {
467 maximize_(max, *coord);
468 minimize_(min, *coord);
469 }
470 qh_fprintf(fp, 9386, " %d: %8.4g %8.4g difference= %4.4g\n", k, min, max, max-min);
471 }
472 if (!qh_QUICKhelp) {
473 qh_fprintf(fp, 9387, "\n\
474 If the input should be full dimensional, you have several options that\n\
475 may determine an initial simplex:\n\
476 - use 'QJ' to joggle the input and make it full dimensional\n\
477 - use 'QbB' to scale the points to the unit cube\n\
478 - use 'QR0' to randomly rotate the input for different maximum points\n\
479 - use 'Qs' to search all points for the initial simplex\n\
480 - use 'En' to specify a maximum roundoff error less than %2.2g.\n\
481 - trace execution with 'T3' to see the determinant for each point.\n",
482 qh DISTround);
483 #if REALfloat
484 qh_fprintf(fp, 9388, "\
485 - recompile qhull for realT precision(#define REALfloat 0 in libqhull.h).\n");
486 #endif
487 qh_fprintf(fp, 9389, "\n\
488 If the input is lower dimensional:\n\
489 - use 'QJ' to joggle the input and make it full dimensional\n\
490 - use 'Qbk:0Bk:0' to delete coordinate k from the input. You should\n\
491 pick the coordinate with the least range. The hull will have the\n\
492 correct topology.\n\
493 - determine the flat containing the points, rotate the points\n\
494 into a coordinate plane, and delete the other coordinates.\n\
495 - add one or more points to make the input full dimensional.\n\
496 ");
497 if (qh DELAUNAY && !qh ATinfinity)
498 qh_fprintf(fp, 9390, "\n\n\
499 This is a Delaunay triangulation and the input is co-circular or co-spherical:\n\
500 - use 'Qz' to add a point \"at infinity\" (i.e., above the paraboloid)\n\
501 - or use 'QJ' to joggle the input and avoid co-circular data\n");
502 }
503 } /* printhelp_singular */
504
505 /*-<a href="qh-globa.htm#TOC"
506 >-------------------------------</a><a name="user_memsizes">-</a>
507
508 qh_user_memsizes()
509 allocate up to 10 additional, quick allocation sizes
510
511 notes:
512 increase maximum number of allocations in qh_initqhull_mem()
513 */
qh_user_memsizes(void)514 void qh_user_memsizes(void) {
515
516 /* qh_memsize(size); */
517 } /* user_memsizes */
518
519
520