📄 user_eg2.c
字号:
/*<html><pre> -<a href="qh-c.htm#user"
>-------------------------------</a><a name="TOP">-</a>
user_eg2.c
sample code for calling qhull() from an application.
Everything here can be done more simply with qh_new_qhull() [see
user_eg.c]. The format used here gives the caller more control
over Qhull. See unix.c for another example.
call with:
user_eg2 "cube/diamond options" "delaunay options" "halfspace options"
for example:
user_eg2 # return summaries
user_eg2 "n" "o" "Fp" # return normals, OFF, points
user_eg2 "QR0 p" "QR0 v p" "QR0 Fp" # rotate input and return points
# 'v' returns Voronoi
# transform is rotated for halfspaces
main() makes three runs of qhull.
1) compute the convex hull of a cube, and incrementally add a diamond
2a) compute the Delaunay triangulation of random points, and add points.
2b) find the Delaunay triangle closest to a point.
3) compute the halfspace intersection of a diamond, and add a cube
notes:
summaries are sent to stderr if other output formats are used
derived from unix.c and compiled by 'make user_eg2'
see qhull.h for data structures, macros, and user-callable functions.
If you want to control all output to stdio and input to stdin,
set the #if below to "1" and delete all lines that contain "io.c".
This prevents the loading of io.o. Qhull will
still write to 'qh ferr' (stderr) for error reporting and tracing.
Defining #if 1, also prevents user.o from being loaded.
*/
#include "qhull_a.h"
char qh_version[] = "user_eg2 98/7/31"; /* used for error messages */
/*-------------------------------------------------
-internal function prototypes
*/
void print_summary (void);
void makecube (coordT *points, int numpoints, int dim);
void adddiamond (coordT *points, int numpoints, int numnew, int dim);
void makeDelaunay (coordT *points, int numpoints, int dim);
void addDelaunay (coordT *points, int numpoints, int numnew, int dim);
void findDelaunay (int dim);
void makehalf (coordT *points, int numpoints, int dim);
void addhalf (coordT *points, int numpoints, int numnew, int dim, coordT *feasible);
/*-------------------------------------------------
-print_summary()
*/
void print_summary (void) {
facetT *facet;
int k;
printf ("\n%d vertices and %d facets with normals:\n",
qh num_vertices, qh num_facets);
FORALLfacets {
for (k=0; k < qh hull_dim; k++)
printf ("%6.2g ", facet->normal[k]);
printf ("\n");
}
}
/*--------------------------------------------------
-makecube- set points to vertices of cube
points is numpoints X dim
*/
void makecube (coordT *points, int numpoints, int dim) {
int j,k;
coordT *point;
for (j=0; j<numpoints; j++) {
point= points + j*dim;
for (k=dim; k--; ) {
if (j & ( 1 << k))
point[k]= 1.0;
else
point[k]= -1.0;
}
}
} /*.makecube.*/
/*--------------------------------------------------
-adddiamond- add diamond to convex hull
points is numpoints+numnew X dim.
notes:
qh_addpoint() does not make a copy of the point coordinates.
For inside points and some outside points, qh_findbestfacet performs
an exhaustive search for a visible facet. Algorithms that retain
previously constructed hulls should be faster for on-line construction
of the convex hull.
*/
void adddiamond (coordT *points, int numpoints, int numnew, int dim) {
int j,k;
coordT *point;
facetT *facet;
boolT isoutside;
realT bestdist;
for (j= 0; j < numnew ; j++) {
point= points + (numpoints+j)*dim;
if (points == qh first_point) /* in case of 'QRn' */
qh num_points= numpoints+j+1;
/* qh num_points sets the size of the points array. You may
allocate the points elsewhere. If so, qh_addpoint records
the point's address in qh other_points
*/
for (k=dim; k--; ) {
if (j/2 == k)
point[k]= (j & 1) ? 2.0 : -2.0;
else
point[k]= 0.0;
}
facet= qh_findbestfacet (point, !qh_ALL, &bestdist, &isoutside);
if (isoutside) {
if (!qh_addpoint (point, facet, False))
break; /* user requested an early exit with 'TVn' or 'TCn' */
}
printf ("%d vertices and %d facets\n",
qh num_vertices, qh num_facets);
/* qh_produce_output(); */
}
if (qh DOcheckmax)
qh_check_maxout();
else if (qh KEEPnearinside)
qh_nearcoplanar();
} /*.adddiamond.*/
/*--------------------------------------------------
-makeDelaunay- set points for dim-1 Delaunay triangulation of random points
points is numpoints X dim. Each point is projected to a paraboloid.
*/
void makeDelaunay (coordT *points, int numpoints, int dim) {
int j,k, seed;
coordT *point, realr;
seed= time(NULL);
printf ("seed: %d\n", seed);
qh_RANDOMseed_( seed);
for (j=0; j<numpoints; j++) {
point= points + j*dim;
for (k= 0; k < dim-1; k++) {
realr= qh_RANDOMint;
point[k]= 2.0 * realr/(qh_RANDOMmax+1) - 1.0;
}
}
qh_setdelaunay (dim, numpoints, points);
} /*.makeDelaunay.*/
/*--------------------------------------------------
-addDelaunay- add points to dim-1 Delaunay triangulation
points is numpoints+numnew X dim. Each point is projected to a paraboloid.
notes:
qh_addpoint() does not make a copy of the point coordinates.
Since qh_addpoint() is not given a visible facet, it performs a directed
search of all facets. Algorithms that retain previously
constructed hulls may be faster.
*/
void addDelaunay (coordT *points, int numpoints, int numnew, int dim) {
int j,k;
coordT *point, realr;
facetT *facet;
realT bestdist;
boolT isoutside;
for (j= 0; j < numnew ; j++) {
point= points + (numpoints+j)*dim;
if (points == qh first_point) /* in case of 'QRn' */
qh num_points= numpoints+j+1;
/* qh num_points sets the size of the points array. You may
allocate the point elsewhere. If so, qh_addpoint records
the point's address in qh other_points
*/
for (k= 0; k < dim-1; k++) {
realr= qh_RANDOMint;
point[k]= 2.0 * realr/(qh_RANDOMmax+1) - 1.0;
}
qh_setdelaunay (dim, 1, point);
facet= qh_findbestfacet (point, !qh_ALL, &bestdist, &isoutside);
if (isoutside) {
if (!qh_addpoint (point, facet, False))
break; /* user requested an early exit with 'TVn' or 'TCn' */
}
qh_printpoint (stdout, "added point", point);
printf ("%d points, %d extra points, %d vertices, and %d facets in total\n",
qh num_points, qh_setsize (qh other_points),
qh num_vertices, qh num_facets);
/* qh_produce_output(); */
}
if (qh DOcheckmax)
qh_check_maxout();
else if (qh KEEPnearinside)
qh_nearcoplanar();
} /*.addDelaunay.*/
/*--------------------------------------------------
-findDelaunay- find Delaunay triangle for [0.5,0.5,...]
assumes dim < 100
*/
void findDelaunay (int dim) {
int k;
coordT point[ 100];
boolT isoutside;
realT bestdist;
facetT *facet;
vertexT *vertex, **vertexp;
for (k= 0; k < dim-1; k++)
point[k]= 0.5;
qh_setdelaunay (dim, 1, point);
facet= qh_findbestfacet (point, qh_ALL, &bestdist, &isoutside);
FOREACHvertex_(facet->vertices) {
for (k=0; k < dim-1; k++)
printf ("%5.2f ", vertex->point[k]);
printf ("\n");
}
} /*.findDelaunay.*/
/*--------------------------------------------------
-makehalf- set points to halfspaces for a (dim)-d diamond
points is numpoints X dim+1
each halfspace consists of dim coefficients followed by an offset
*/
void makehalf (coordT *points, int numpoints, int dim) {
int j,k;
coordT *point;
for (j=0; j<numpoints; j++) {
point= points + j*(dim+1);
point[dim]= -1.0; /* offset */
for (k=dim; k--; ) {
if (j & ( 1 << k))
point[k]= 1.0;
else
point[k]= -1.0;
}
}
} /*.makehalf.*/
/*--------------------------------------------------
-addhalf- add halfspaces for a (dim)-d cube to the intersection
points is numpoints+numnew X dim+1
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -