user_eg_r.c
Go to the documentation of this file.
1 /*<html><pre> -<a href="../libqhull_r/qh-user_r.htm"
2  >-------------------------------</a><a name="TOP">-</a>
3 
4  user_eg_r.c
5  sample code for calling qhull() from an application. Uses reentrant libqhull_r
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_r.c and user_eg2_r.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_r/qhull_ra.h"
50 
51 /*-------------------------------------------------
52 -internal function prototypes
53 */
54 void print_summary(qhT *qh);
55 void makecube(coordT *points, int numpoints, int dim);
56 void makeDelaunay(qhT *qh, coordT *points, int numpoints, int dim, int seed);
57 void findDelaunay(qhT *qh, int dim);
58 void makehalf(coordT *points, int numpoints, int dim);
59 
60 /*-------------------------------------------------
61 -print_summary(qh)
62 */
63 void print_summary(qhT *qh) {
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(qhT *qh, 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_(qh, 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(qhT *qh, 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(qh, dim+1, 1, point);
136  facet= qh_findbestfacet(qh, 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, 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  qhT qh_qh; /* Qhull's data structure. First argument of most calls */
201  qhT *qh= &qh_qh;
202 
203  QHULL_LIB_CHECK
204 
205  qh_zero(qh, errfile);
206 
207  printf("This is the output from user_eg_r.c\n\n\
208 It shows how qhull() may be called from an application using the qhull\n\
209 reentrant library. user_eg is not part of qhull itself. If it appears\n\
210 accidently, please remove user_eg_r.c from your project. If it fails\n\
211 immediately, user_eg_r.c was incorrectly linked to the non-reentrant library.\n\
212 Also try 'user_eg T1 2>&1'\n\n");
213 
214  /*
215  Run 1: convex hull
216  */
217  printf( "\ncompute convex hull of cube after rotating input\n");
218  sprintf(flags, "qhull s Tcv %s", argc >= 2 ? argv[1] : "");
219  numpoints= SIZEcube;
220  makecube(points, numpoints, DIM);
221  for (i=numpoints; i--; )
222  rows[i]= points+dim*i;
223  qh_printmatrix(qh, outfile, "input", rows, numpoints, dim);
224  exitcode= qh_new_qhull(qh, dim, numpoints, points, ismalloc,
225  flags, outfile, errfile);
226  if (!exitcode) { /* if no error */
227  /* 'qh->facet_list' contains the convex hull */
228  print_summary(qh);
229  FORALLfacets {
230  /* ... your code ... */
231  }
232  }
233  qh_freeqhull(qh, !qh_ALL); /* free long memory */
234  qh_memfreeshort(qh, &curlong, &totlong); /* free short memory and memory allocator */
235  if (curlong || totlong)
236  fprintf(errfile, "qhull internal warning (user_eg, #1): did not free %d bytes of long memory (%d pieces)\n", totlong, curlong);
237 
238  /*
239  Run 2: Delaunay triangulation, reusing the previous qh/qh_qh
240  */
241 
242  printf( "\ncompute %d-d Delaunay triangulation\n", dim);
243  sprintf(flags, "qhull s d Tcv %s", argc >= 3 ? argv[2] : "");
244  numpoints= SIZEcube;
245  makeDelaunay(qh, points, numpoints, dim, (int)time(NULL));
246  for (i=numpoints; i--; )
247  rows[i]= points+dim*i;
248  qh_printmatrix(qh, outfile, "input", rows, numpoints, dim);
249  exitcode= qh_new_qhull(qh, dim, numpoints, points, ismalloc,
250  flags, outfile, errfile);
251  if (!exitcode) { /* if no error */
252  /* 'qh->facet_list' contains the convex hull */
253  /* If you want a Voronoi diagram ('v') and do not request output (i.e., outfile=NULL),
254  call qh_setvoronoi_all() after qh_new_qhull(). */
255  print_summary(qh);
256  FORALLfacets {
257  /* ... your code ... */
258  }
259  printf( "\nfind %d-d Delaunay triangle closest to [0.5, 0.5, ...]\n", dim);
260  exitcode= setjmp(qh->errexit);
261  if (!exitcode) {
262  /* Trap Qhull errors in findDelaunay(). Without the setjmp(), Qhull
263  will exit() after reporting an error */
264  qh->NOerrexit= False;
265  findDelaunay(qh, DIM);
266  }
267  qh->NOerrexit= True;
268  }
269  {
270  coordT pointsB[DIM*TOTpoints]; /* array of coordinates for each point */
271 
272  qhT qh_qhB; /* Create a new instance of Qhull (qhB) */
273  qhT *qhB= &qh_qhB;
274  qh_zero(qhB, errfile);
275 
276  printf( "\nCompute a new triangulation as a separate instance of Qhull\n");
277  sprintf(flags, "qhull s d Tcv %s", argc >= 3 ? argv[2] : "");
278  numpoints= SIZEcube;
279  makeDelaunay(qhB, pointsB, numpoints, dim, (int)time(NULL)+1);
280  for (i=numpoints; i--; )
281  rows[i]= pointsB+dim*i;
282  qh_printmatrix(qhB, outfile, "input", rows, numpoints, dim);
283  exitcode= qh_new_qhull(qhB, dim, numpoints, pointsB, ismalloc,
284  flags, outfile, errfile);
285  if (!exitcode)
286  print_summary(qhB);
287  printf( "\nFree memory allocated by the new instance of Qhull, and redisplay the old results.\n");
288  qh_freeqhull(qhB, !qh_ALL); /* free long memory */
289  qh_memfreeshort(qhB, &curlong, &totlong); /* free short memory and memory allocator */
290  if (curlong || totlong)
291  fprintf(errfile, "qhull internal warning (user_eg, #4): did not free %d bytes of long memory (%d pieces)\n", totlong, curlong);
292 
293  printf( "\n\n");
294  print_summary(qh); /* The other instance is unchanged */
295  /* Exiting the block frees qh_qhB */
296  }
297  qh_freeqhull(qh, !qh_ALL); /* free long memory */
298  qh_memfreeshort(qh, &curlong, &totlong); /* free short memory and memory allocator */
299  if (curlong || totlong)
300  fprintf(errfile, "qhull internal warning (user_eg, #2): did not free %d bytes of long memory (%d pieces)\n", totlong, curlong);
301 
302  /*
303  Run 3: halfspace intersection about the origin
304  */
305  printf( "\ncompute halfspace intersection about the origin for a diamond\n");
306  sprintf(flags, "qhull H0 s Tcv %s", argc >= 4 ? argv[3] : "Fp");
307  numpoints= SIZEcube;
308  makehalf(points, numpoints, dim);
309  for (i=numpoints; i--; )
310  rows[i]= points+(dim+1)*i;
311  qh_printmatrix(qh, outfile, "input as halfspace coefficients + offsets", rows, numpoints, dim+1);
312  /* use qh_sethalfspace_all to transform the halfspaces yourself.
313  If so, set 'qh->feasible_point and do not use option 'Hn,...' [it would retransform the halfspaces]
314  */
315  exitcode= qh_new_qhull(qh, dim+1, numpoints, points, ismalloc,
316  flags, outfile, errfile);
317  if (!exitcode)
318  print_summary(qh);
319  qh_freeqhull(qh, !qh_ALL);
320  qh_memfreeshort(qh, &curlong, &totlong);
321  if (curlong || totlong) /* could also check previous runs */
322  fprintf(stderr, "qhull internal warning (user_eg, #3): did not free %d bytes of long memory (%d pieces)\n",
323  totlong, curlong);
324  return exitcode;
325 } /* main */
326 
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_r.c:185
facetT::tricoplanar
flagT tricoplanar
Definition: libqhull.h:314
QHULL_LIB_CHECK
#define QHULL_LIB_CHECK
Definition: libqhull.h:462
coordT
#define coordT
Definition: libqhull.h:80
makecube
void makecube(coordT *points, int numpoints, int dim)
Definition: user_eg_r.c:80
qh_zero
void qh_zero(qhT *qh, FILE *errfile)
Definition: global_r.c:2095
seed
void seed(unsigned int seed_value)
qh_ALL
#define qh_ALL
Definition: libqhull.h:180
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
realT
#define realT
Definition: user.h:154
rows
int rows
vertexT::point
pointT * point
Definition: libqhull.h:399
TOTpoints
#define TOTpoints
Definition: user_eg_r.c:174
qhT
Definition: libqhull.h:465
facetT
Definition: libqhull.h:262
qh_RANDOMseed_
#define qh_RANDOMseed_(seed)
Definition: user.h:284
qh_setdelaunay
void qh_setdelaunay(int dim, int count, pointT *points)
Definition: geom2.c:1806
boolT
#define boolT
Definition: libqhull.h:121
DIM
#define DIM
Definition: user_eg_r.c:171
qh_qh
qhT qh_qh
Definition: global.c:26
FORALLfacets
#define FORALLfacets
Definition: libqhull.h:840
False
#define False
Definition: libqhull.h:128
qh
#define qh
Definition: libqhull.h:457
qh_freeqhull
void qh_freeqhull(boolT allmem)
Definition: global.c:425
qh_errexit
void qh_errexit(int exitcode, facetT *facet, ridgeT *ridge)
Definition: user.c:213
qh_RANDOMint
#define qh_RANDOMint
Definition: user.h:283
qhull_ra.h
makeDelaunay
void makeDelaunay(qhT *qh, coordT *points, int numpoints, int dim, int seed)
Definition: user_eg_r.c:101
qh_printmatrix
void qh_printmatrix(FILE *fp, const char *string, realT **rows, int numrow, int numcol)
Definition: geom2.c:1345
makehalf
void makehalf(coordT *points, int numpoints, int dim)
Definition: user_eg_r.c:155
SIZEcube
#define SIZEcube
Definition: user_eg_r.c:172
qh_RANDOMmax
#define qh_RANDOMmax
Definition: user.h:282
findDelaunay
void findDelaunay(qhT *qh, int dim)
Definition: user_eg_r.c:125
qh_ERRqhull
#define qh_ERRqhull
Definition: libqhull.h:198
print_summary
void print_summary(qhT *qh)
Definition: user_eg_r.c:63
qh_memfreeshort
void qh_memfreeshort(int *curlong, int *totlong)
Definition: mem.c:288
vertexT
Definition: libqhull.h:396
FOREACHvertex_
#define FOREACHvertex_(vertices)
Definition: libqhull.h:950
facetT::normal
coordT * normal
Definition: libqhull.h:274
dim
int dim
facetT::id
unsigned id
Definition: libqhull.h:311
facetT::vertices
setT * vertices
Definition: libqhull.h:295
True
#define True
Definition: libqhull.h:129

hpp-fcl
Author(s):
autogenerated on Fri Aug 2 2024 02:45:15