user_eg.c
Go to the documentation of this file.
1 /*<html><pre> -<a href="../libqhull/qh-user.htm"
2  >-------------------------------</a><a name="TOP">-</a>
3 
4  user_eg.c
5  sample code for calling qhull() from an application
6 
7  call with:
8 
9  user_eg "cube/diamond options" "delaunay options" "halfspace options"
10 
11  for example:
12 
13  user_eg # return summaries
14 
15  user_eg "n" "o" "Fp" # return normals, OFF, points
16 
17  user_eg "n Qt" "o" "Fp" # triangulated cube
18 
19  user_eg "QR0 p" "QR0 v p" "QR0 Fp" # rotate input and return points
20  # 'v' returns Voronoi
21  # transform is rotated for halfspaces
22 
23  main() makes three runs of qhull.
24 
25  1) compute the convex hull of a cube
26 
27  2a) compute the Delaunay triangulation of random points
28 
29  2b) find the Delaunay triangle closest to a point.
30 
31  3) compute the halfspace intersection of a diamond
32 
33  notes:
34 
35  For another example, see main() in unix.c and user_eg2.c.
36  These examples, call qh_qhull() directly. They allow
37  tighter control on the code loaded with Qhull.
38 
39  For a C++ example, see user_eg3/user_eg3_r.cpp
40 
41  Summaries are sent to stderr if other output formats are used
42 
43  compiled by 'make bin/user_eg'
44 
45  see libqhull.h for data structures, macros, and user-callable functions.
46 */
47 
48 #define qh_QHimport
49 #include "libqhull/qhull_a.h"
50 
51 /*-------------------------------------------------
52 -internal function prototypes
53 */
54 void print_summary(void);
55 void makecube(coordT *points, int numpoints, int dim);
56 void makeDelaunay(coordT *points, int numpoints, int dim, int seed);
57 void findDelaunay(int dim);
58 void makehalf(coordT *points, int numpoints, int dim);
59 
60 /*-------------------------------------------------
61 -print_summary()
62 */
63 void print_summary(void) {
64  facetT *facet;
65  int k;
66 
67  printf("\n%d vertices and %d facets with normals:\n",
68  qh num_vertices, qh num_facets);
69  FORALLfacets {
70  for (k=0; k < qh hull_dim; k++)
71  printf("%6.2g ", facet->normal[k]);
72  printf("\n");
73  }
74 }
75 
76 /*--------------------------------------------------
77 -makecube- set points to vertices of cube
78  points is numpoints X dim
79 */
80 void makecube(coordT *points, int numpoints, int dim) {
81  int j,k;
82  coordT *point;
83 
84  for (j=0; j<numpoints; j++) {
85  point= points + j*dim;
86  for (k=dim; k--; ) {
87  if (j & ( 1 << k))
88  point[k]= 1.0;
89  else
90  point[k]= -1.0;
91  }
92  }
93 } /*.makecube.*/
94 
95 /*--------------------------------------------------
96 -makeDelaunay- set points for dim Delaunay triangulation of random points
97  points is numpoints X dim.
98 notes:
99  makeDelaunay() in user_eg2.c uses qh_setdelaunay() to project points in place.
100 */
101 void makeDelaunay(coordT *points, int numpoints, int dim, int seed) {
102  int j,k;
103  coordT *point, realr;
104 
105  printf("seed: %d\n", seed);
106  qh_RANDOMseed_( seed);
107  for (j=0; j<numpoints; j++) {
108  point= points + j*dim;
109  for (k= 0; k < dim; k++) {
110  realr= qh_RANDOMint;
111  point[k]= 2.0 * realr/(qh_RANDOMmax+1) - 1.0;
112  }
113  }
114 } /*.makeDelaunay.*/
115 
116 /*--------------------------------------------------
117 -findDelaunay- find Delaunay triangle for [0.5,0.5,...]
118  assumes dim < 100
119 notes:
120  calls qh_setdelaunay() to project the point to a parabaloid
121 warning:
122  This is not implemented for tricoplanar facets ('Qt'),
123  See <a href="../html/qh-code.htm#findfacet">locate a facet with qh_findbestfacet()</a>
124 */
125 void findDelaunay(int dim) {
126  int k;
127  coordT point[ 100];
128  boolT isoutside;
129  realT bestdist;
130  facetT *facet;
131  vertexT *vertex, **vertexp;
132 
133  for (k= 0; k < dim; k++)
134  point[k]= 0.5;
135  qh_setdelaunay(dim+1, 1, point);
136  facet= qh_findbestfacet(point, qh_ALL, &bestdist, &isoutside);
137  if (facet->tricoplanar) {
138  fprintf(stderr, "findDelaunay: not implemented for triangulated, non-simplicial Delaunay regions (tricoplanar facet, f%d).\n",
139  facet->id);
140  qh_errexit(qh_ERRqhull, facet, NULL);
141  }
142  FOREACHvertex_(facet->vertices) {
143  for (k=0; k < dim; k++)
144  printf("%5.2f ", vertex->point[k]);
145  printf("\n");
146  }
147 } /*.findDelaunay.*/
148 
149 /*--------------------------------------------------
150 -makehalf- set points to halfspaces for a (dim)-dimensional diamond
151  points is numpoints X dim+1
152 
153  each halfspace consists of dim coefficients followed by an offset
154 */
155 void makehalf(coordT *points, int numpoints, int dim) {
156  int j,k;
157  coordT *point;
158 
159  for (j=0; j<numpoints; j++) {
160  point= points + j*(dim+1);
161  point[dim]= -1.0; /* offset */
162  for (k=dim; k--; ) {
163  if (j & ( 1 << k))
164  point[k]= 1.0;
165  else
166  point[k]= -1.0;
167  }
168  }
169 } /*.makehalf.*/
170 
171 #define DIM 3 /* dimension of points, must be < 31 for SIZEcube */
172 #define SIZEcube (1<<DIM)
173 #define SIZEdiamond (2*DIM)
174 #define TOTpoints (SIZEcube + SIZEdiamond)
175 
176 /*--------------------------------------------------
177 -main- derived from Qhull-template in user.c
178 
179  see program header
180 
181  this contains three runs of Qhull for convex hull, Delaunay
182  triangulation or Voronoi vertices, and halfspace intersection
183 
184 */
185 int main(int argc, char *argv[]) {
186  int dim= DIM; /* dimension of points */
187  int numpoints; /* number of points */
188  coordT points[(DIM+1)*TOTpoints]; /* array of coordinates for each point */
189  coordT *rows[TOTpoints];
190  boolT ismalloc= False; /* True if qhull should free points in qh_freeqhull() or reallocation */
191  char flags[250]; /* option flags for qhull, see qh-quick.htm */
192  FILE *outfile= stdout; /* output from qh_produce_output()
193  use NULL to skip qh_produce_output() */
194  FILE *errfile= stderr; /* error messages from qhull code */
195  int exitcode; /* 0 if no error from qhull */
196  facetT *facet; /* set by FORALLfacets */
197  int curlong, totlong; /* memory remaining after qh_memfreeshort */
198  int i;
199 
200  QHULL_LIB_CHECK
201 
202  printf("This is the output from user_eg.c\n\n\
203 It shows how qhull() may be called from an application using the qhull\n\
204 shared library. user_eg is not part of qhull itself. If it appears\n\
205 accidently, please remove user_eg.c from your project. If it fails\n\
206 immediately, user_eg.c was incorrectly linked to the reentrant library.\n\
207 Also try 'user_eg T1 2>&1'\n\n");
208 
209 #if qh_QHpointer /* see user.h */
210  if (qh_qh){
211  printf("QH6233: Qhull link error. The global variable qh_qh was not initialized\n\
212 to NULL by global.c. Please compile user_eg.c with -Dqh_QHpointer_dllimport\n\
213 as well as -Dqh_QHpointer, or use libqhullstatic, or use a different tool chain.\n\n");
214  return -1;
215  }
216 #endif
217 
218  /*
219  Run 1: convex hull
220  */
221  printf( "\ncompute convex hull of cube after rotating input\n");
222  sprintf(flags, "qhull s Tcv %s", argc >= 2 ? argv[1] : "");
223  numpoints= SIZEcube;
224  makecube(points, numpoints, DIM);
225  for (i=numpoints; i--; )
226  rows[i]= points+dim*i;
227  qh_printmatrix(outfile, "input", rows, numpoints, dim);
228  exitcode= qh_new_qhull(dim, numpoints, points, ismalloc,
229  flags, outfile, errfile);
230  if (!exitcode) { /* if no error */
231  /* 'qh facet_list' contains the convex hull */
232  print_summary();
233  FORALLfacets {
234  /* ... your code ... */
235  }
236  }
237  qh_freeqhull(!qh_ALL); /* free long memory */
238  qh_memfreeshort(&curlong, &totlong); /* free short memory and memory allocator */
239  if (curlong || totlong)
240  fprintf(errfile, "qhull internal warning (user_eg, #1): did not free %d bytes of long memory (%d pieces)\n", totlong, curlong);
241 
242  /*
243  Run 2: Delaunay triangulation
244  */
245 
246  printf( "\ncompute %d-d Delaunay triangulation\n", dim);
247  sprintf(flags, "qhull s d Tcv %s", argc >= 3 ? argv[2] : "");
248  numpoints= SIZEcube;
249  makeDelaunay(points, numpoints, dim, (int)time(NULL));
250  for (i=numpoints; i--; )
251  rows[i]= points+dim*i;
252  qh_printmatrix(outfile, "input", rows, numpoints, dim);
253  exitcode= qh_new_qhull(dim, numpoints, points, ismalloc,
254  flags, outfile, errfile);
255  if (!exitcode) { /* if no error */
256  /* 'qh facet_list' contains the convex hull */
257  /* If you want a Voronoi diagram ('v') and do not request output (i.e., outfile=NULL),
258  call qh_setvoronoi_all() after qh_new_qhull(). */
259  print_summary();
260  FORALLfacets {
261  /* ... your code ... */
262  }
263  printf( "\nfind %d-d Delaunay triangle closest to [0.5, 0.5, ...]\n", dim);
264  exitcode= setjmp(qh errexit);
265  if (!exitcode) {
266  /* Trap Qhull errors in findDelaunay(). Without the setjmp(), Qhull
267  will exit() after reporting an error */
268  qh NOerrexit= False;
269  findDelaunay(DIM);
270  }
271  qh NOerrexit= True;
272  }
273 #if qh_QHpointer /* see user.h */
274  {
275  qhT *oldqhA, *oldqhB;
276  coordT pointsB[DIM*TOTpoints]; /* array of coordinates for each point */
277 
278  printf( "\nsave first triangulation and compute a new triangulation\n");
279  oldqhA= qh_save_qhull();
280  sprintf(flags, "qhull s d Tcv %s", argc >= 3 ? argv[2] : "");
281  numpoints= SIZEcube;
282  makeDelaunay(pointsB, numpoints, dim, (int)time(NULL)+1);
283  for (i=numpoints; i--; )
284  rows[i]= pointsB+dim*i;
285  qh_printmatrix(outfile, "input", rows, numpoints, dim);
286  exitcode= qh_new_qhull(dim, numpoints, pointsB, ismalloc,
287  flags, outfile, errfile);
288  if (!exitcode)
289  print_summary();
290  printf( "\nsave second triangulation and restore first one\n");
291  oldqhB= qh_save_qhull();
292  qh_restore_qhull(&oldqhA);
293  print_summary();
294  printf( "\nfree first triangulation and restore second one.\n");
295  qh_freeqhull(qh_ALL); /* free short and long memory used by first call */
296  /* do not use qh_memfreeshort */
297  qh_restore_qhull(&oldqhB);
298  print_summary();
299  }
300 #endif
301  qh_freeqhull(!qh_ALL); /* free long memory */
302  qh_memfreeshort(&curlong, &totlong); /* free short memory and memory allocator */
303  if (curlong || totlong)
304  fprintf(errfile, "qhull internal warning (user_eg, #2): did not free %d bytes of long memory (%d pieces)\n", totlong, curlong);
305 
306  /*
307  Run 3: halfspace intersection about the origin
308  */
309  printf( "\ncompute halfspace intersection about the origin for a diamond\n");
310  sprintf(flags, "qhull H0 s Tcv %s", argc >= 4 ? argv[3] : "Fp");
311  numpoints= SIZEcube;
312  makehalf(points, numpoints, dim);
313  for (i=numpoints; i--; )
314  rows[i]= points+(dim+1)*i;
315  qh_printmatrix(outfile, "input as halfspace coefficients + offsets", rows, numpoints, dim+1);
316  /* use qh_sethalfspace_all to transform the halfspaces yourself.
317  If so, set 'qh feasible_point and do not use option 'Hn,...' [it would retransform the halfspaces]
318  */
319  exitcode= qh_new_qhull(dim+1, numpoints, points, ismalloc,
320  flags, outfile, errfile);
321  if (!exitcode)
322  print_summary();
323  qh_freeqhull(!qh_ALL);
324  qh_memfreeshort(&curlong, &totlong);
325  if (curlong || totlong) /* could also check previous runs */
326  fprintf(stderr, "qhull internal warning (user_eg, #3): did not free %d bytes of long memory (%d pieces)\n",
327  totlong, curlong);
328  return exitcode;
329 } /* main */
330 
qhT
Definition: libqhull.h:465
qh_findbestfacet
facetT * qh_findbestfacet(pointT *point, boolT bestoutside, realT *bestdist, boolT *isoutside)
Definition: poly2.c:1239
main
int main(int argc, char *argv[])
Definition: user_eg.c:185
qh_ERRqhull
#define qh_ERRqhull
Definition: libqhull.h:198
FORALLfacets
#define FORALLfacets
Definition: libqhull.h:840
boolT
#define boolT
Definition: libqhull.h:121
print_summary
void print_summary(void)
Definition: user_eg.c:63
qh_RANDOMmax
#define qh_RANDOMmax
Definition: user.h:282
FOREACHvertex_
#define FOREACHvertex_(vertices)
Definition: libqhull.h:950
rows
int rows
TOTpoints
#define TOTpoints
Definition: user_eg.c:174
qh_qh
qhT qh_qh
Definition: global.c:26
DIM
#define DIM
Definition: user_eg.c:171
qh_RANDOMint
#define qh_RANDOMint
Definition: user.h:283
dim
int dim
facetT::vertices
setT * vertices
Definition: libqhull.h:295
coordT
#define coordT
Definition: libqhull.h:80
True
#define True
Definition: libqhull.h:129
qh
#define qh
Definition: libqhull.h:457
qh_setdelaunay
void qh_setdelaunay(int dim, int count, pointT *points)
Definition: geom2.c:1806
vertexT::point
pointT * point
Definition: libqhull.h:399
qh_freeqhull
void qh_freeqhull(boolT allmem)
Definition: global.c:425
makecube
void makecube(coordT *points, int numpoints, int dim)
Definition: user_eg.c:80
qh_printmatrix
void qh_printmatrix(FILE *fp, const char *string, realT **rows, int numrow, int numcol)
Definition: geom2.c:1345
findDelaunay
void findDelaunay(int dim)
Definition: user_eg.c:125
qhull_a.h
makeDelaunay
void makeDelaunay(coordT *points, int numpoints, int dim, int seed)
Definition: user_eg.c:101
QHULL_LIB_CHECK
#define QHULL_LIB_CHECK
Definition: libqhull.h:462
SIZEcube
#define SIZEcube
Definition: user_eg.c:172
qh_memfreeshort
void qh_memfreeshort(int *curlong, int *totlong)
Definition: mem.c:288
makehalf
void makehalf(coordT *points, int numpoints, int dim)
Definition: user_eg.c:155
qh_new_qhull
int qh_new_qhull(int dim, int numpoints, coordT *points, boolT ismalloc, char *qhull_cmd, FILE *outfile, FILE *errfile)
Definition: user.c:129
qh_errexit
void qh_errexit(int exitcode, facetT *facet, ridgeT *ridge)
Definition: user.c:213
False
#define False
Definition: libqhull.h:128
qh_ALL
#define qh_ALL
Definition: libqhull.h:180
facetT::id
unsigned id
Definition: libqhull.h:311
qh_RANDOMseed_
#define qh_RANDOMseed_(seed)
Definition: user.h:284
realT
#define realT
Definition: user.h:154
facetT::normal
coordT * normal
Definition: libqhull.h:274
facetT::tricoplanar
flagT tricoplanar
Definition: libqhull.h:314

hpp-fcl
Author(s):
autogenerated on Fri Jun 2 2023 02:39:02