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