Quick hull implementation

源代码在线查看: user.c

软件大小: 517 K
上传用户: focuslocus
关键词: implementation Quick hull
下载地址: 免注册下载 普通下载 VIP

相关代码

				/*  -				  >--------------------------------
				
				   user.c 
				   user redefinable functions
				
				   see README.txt  see COPYING.txt for copyright information.
				
				   see qhull.h for data structures, macros, and user-callable functions.
				
				   see user_eg.c, unix.c, and qhull_interface.cpp for examples.
				
				   see user.h for user-definable constants
				
				      use qh_NOmem in mem.h to turn off memory management
				      use qh_NOmerge in user.h to turn off facet merging
				      set qh_KEEPstatistics in user.h to 0 to turn off statistics
				
				   This is unsupported software.  You're welcome to make changes,
				   but you're on your own if something goes wrong.  Use 'Tc' to
				   check frequently.  Usually qhull will report an error if 
				   a data structure becomes inconsistent.  If so, it also reports
				   the last point added to the hull, e.g., 102.  You can then trace
				   the execution of qhull with "T4P102".  
				
				   Please report any errors that you fix to qhull@qhull.org
				
				   call_qhull is a template for calling qhull from within your application
				
				   if you recompile and load this module, then user.o will not be loaded
				   from qhull.a
				
				   you can add additional quick allocation sizes in qh_user_memsizes
				
				   if the other functions here are redefined to not use qh_print...,
				   then io.o will not be loaded from qhull.a.  See user_eg.c for an
				   example.  We recommend keeping io.o for the extra debugging 
				   information it supplies.
				*/
				
				#include "qhull_a.h" 
				
				/*-				  >--------------------------------
				
				  qh_call_qhull( void )
				    template for calling qhull from inside your program
				    remove #if 0, #endif to compile
				
				  returns: 
				    exit code (see qh_ERR... in qhull.h)
				    all memory freed
				
				  notes:
				    This can be called any number of times.  
				
				  see:
				    qh_call_qhull_once()
				    
				*/
				#if 0
				{
				  int dim;	            /* dimension of points */
				  int numpoints;            /* number of points */
				  coordT *points;           /* array of coordinates for each point */
				  boolT ismalloc;           /* True if qhull should free points in qh_freeqhull() or reallocation */
				  char flags[]= "qhull Tv"; /* option flags for qhull, see qh_opt.htm */
				  FILE *outfile= stdout;    /* output from qh_produce_output()
							       use NULL to skip qh_produce_output() */
				  FILE *errfile= stderr;    /* error messages from qhull code */
				  int exitcode;             /* 0 if no error from qhull */
				  facetT *facet;	    /* set by FORALLfacets */
				  int curlong, totlong;	    /* memory remaining after qh_memfreeshort */
				
				  /* initialize dim, numpoints, points[], ismalloc here */
				  exitcode= qh_new_qhull (dim, numpoints, points, ismalloc,
				                      flags, outfile, errfile); 
				  if (!exitcode) {                  /* if no error */
				    /* 'qh facet_list' contains the convex hull */
				    FORALLfacets {
				       /* ... your code ... */
				    }
				  }
				  qh_freeqhull(!qh_ALL);
				  qh_memfreeshort (&curlong, &totlong);
				  if (curlong || totlong) 
				    fprintf (errfile, "qhull internal warning (main): did not free %d bytes of long memory (%d pieces)\n", totlong, curlong);
				}
				#endif
				
				/*-				  >--------------------------------
				
				  qh_new_qhull( dim, numpoints, points, ismalloc, qhull_cmd, outfile, errfile )
				    build new qhull data structure and return exitcode (0 if no errors)
				
				  notes:
				    do not modify points until finished with results.
				      The qhull data structure contains pointers into the points array.
				    do not call qhull functions before qh_new_qhull().
				      The qhull data structure is not initialized until qh_new_qhull().
				
				    outfile may be null
				    qhull_cmd must start with "qhull "
				    projects points to a new point array for Delaunay triangulations ('d' and 'v')
				    transforms points into a new point array for halfspace intersection ('H')
				       
				
				  To allow multiple, concurrent calls to qhull() 
				    - set qh_QHpointer in user.h
				    - use qh_save_qhull and qh_restore_qhull to swap the global data structure between calls.
				    - use qh_freeqhull(qh_ALL) to free intermediate convex hulls
				
				  see:
				    user_eg.c for an example
				*/
				int qh_new_qhull (int dim, int numpoints, coordT *points, boolT ismalloc, 
						char *qhull_cmd, FILE *outfile, FILE *errfile) {
				  int exitcode, hulldim;
				  boolT new_ismalloc;
				  static boolT firstcall = True;
				  coordT *new_points;
				
				  if (firstcall) {
				    qh_meminit (errfile);
				    firstcall= False;
				  }
				  if (strncmp (qhull_cmd,"qhull ", 6)) {
				    fprintf (errfile, "qh_new_qhull: start qhull_cmd argument with \"qhull \"\n");
				    exit(1);
				  }
				  qh_initqhull_start (NULL, outfile, errfile);
				  trace1(( qh ferr, "qh_new_qhull: build new Qhull for %d %d-d points with %s\n", numpoints, dim, qhull_cmd));
				  exitcode = setjmp (qh errexit);
				  if (!exitcode)
				  {
				    qh NOerrexit = False;
				    qh_initflags (qhull_cmd);
				    if (qh DELAUNAY)
				      qh PROJECTdelaunay= True;
				    if (qh HALFspace) {
				      /* points is an array of halfspaces, 
				         the last coordinate of each halfspace is its offset */
				      hulldim= dim-1;
				      qh_setfeasible (hulldim); 
				      new_points= qh_sethalfspace_all (dim, numpoints, points, qh feasible_point);
				      new_ismalloc= True;
				      if (ismalloc)
					free (points);
				    }else {
				      hulldim= dim;
				      new_points= points;
				      new_ismalloc= ismalloc;
				    }
				    qh_init_B (new_points, numpoints, hulldim, new_ismalloc);
				    qh_qhull();
				    qh_check_output();
				    if (outfile)
				      qh_produce_output(); 
				    if (qh VERIFYoutput && !qh STOPpoint && !qh STOPcone)
				      qh_check_points();
				  }
				  qh NOerrexit = True;
				  return exitcode;
				} /* new_qhull */
				
				/*-				  >--------------------------------
				  
				  qh_errexit( exitcode, facet, ridge )
				    report and exit from an error
				    report facet and ridge if non-NULL
				    reports useful information such as last point processed
				    set qh.FORCEoutput to print neighborhood of facet
				
				  see: 
				    qh_errexit2() in qhull.c for printing 2 facets
				
				  design:
				    check for error within error processing
				    compute qh.hulltime
				    print facet and ridge (if any)
				    report commandString, options, qh.furthest_id
				    print summary and statistics (including precision statistics)
				    if qh_ERRsingular
				      print help text for singular data set
				    exit program via long jump (if defined) or exit()      
				*/
				void qh_errexit(int exitcode, facetT *facet, ridgeT *ridge) {
				
				  if (qh ERREXITcalled) {
				    fprintf (qh ferr, "\nqhull error while processing previous error.  Exit program\n");
				    exit(1);
				  }
				  qh ERREXITcalled= True;
				  if (!qh QHULLfinished)
				    qh hulltime= qh_CPUclock - qh hulltime;
				  qh_errprint("ERRONEOUS", facet, NULL, ridge, NULL);
				  fprintf (qh ferr, "\nWhile executing: %s | %s\n", qh rbox_command, qh qhull_command);
				  fprintf(qh ferr, "Options selected for Qhull %s:\n%s\n", qh_version, qh qhull_options);
				  if (qh furthest_id >= 0) {
				    fprintf(qh ferr, "Last point added to hull was p%d.", qh furthest_id);
				    if (zzval_(Ztotmerge))
				      fprintf(qh ferr, "  Last merge was #%d.", zzval_(Ztotmerge));
				    if (qh QHULLfinished)
				      fprintf(qh ferr, "\nQhull has finished constructing the hull.");
				    else if (qh POSTmerging)
				      fprintf(qh ferr, "\nQhull has started post-merging.");
				    fprintf (qh ferr, "\n");
				  }
				  if (qh FORCEoutput && (qh QHULLfinished || (!facet && !ridge)))
				    qh_produce_output();
				  else {
				    if (exitcode != qh_ERRsingular && zzval_(Zsetplane) > qh hull_dim+1) {
				      fprintf (qh ferr, "\nAt error exit:\n");
				      qh_printsummary (qh ferr);
				      if (qh PRINTstatistics) {
					qh_collectstatistics();
					qh_printstatistics(qh ferr, "at error exit");
					qh_memstatistics (qh ferr);
				      }
				    }
				    if (qh PRINTprecision)
				      qh_printstats (qh ferr, qhstat precision, NULL);
				  }
				  if (!exitcode)
				    exitcode= qh_ERRqhull;
				  else if (exitcode == qh_ERRsingular)
				    qh_printhelp_singular(qh ferr);
				  else if (exitcode == qh_ERRprec && !qh PREmerge)
				    qh_printhelp_degenerate (qh ferr);
				  if (qh NOerrexit) {
				    fprintf (qh ferr, "qhull error while ending program.  Exit program\n");
				    exit(1);
				  }
				  qh NOerrexit= True;
				  longjmp(qh errexit, exitcode);
				} /* errexit */
				
				
				/*-				  >--------------------------------
				  
				  qh_errprint( fp, string, atfacet, otherfacet, atridge, atvertex )
				    prints out the information of facets and ridges to fp
				    also prints neighbors and geomview output
				    
				  notes:
				    except for string, any parameter may be NULL
				*/
				void qh_errprint(char *string, facetT *atfacet, facetT *otherfacet, ridgeT *atridge, vertexT *atvertex) {
				  int i;
				
				  if (atfacet) {
				    fprintf(qh ferr, "%s FACET:\n", string);
				    qh_printfacet(qh ferr, atfacet);
				  }
				  if (otherfacet) {
				    fprintf(qh ferr, "%s OTHER FACET:\n", string);
				    qh_printfacet(qh ferr, otherfacet);
				  }
				  if (atridge) {
				    fprintf(qh ferr, "%s RIDGE:\n", string);
				    qh_printridge(qh ferr, atridge);
				    if (atridge->top && atridge->top != atfacet && atridge->top != otherfacet)
				      qh_printfacet(qh ferr, atridge->top);
				    if (atridge->bottom
					&& atridge->bottom != atfacet && atridge->bottom != otherfacet)
				      qh_printfacet(qh ferr, atridge->bottom);
				    if (!atfacet)
				      atfacet= atridge->top;
				    if (!otherfacet)
				      otherfacet= otherfacet_(atridge, atfacet);
				  }
				  if (atvertex) {
				    fprintf(qh ferr, "%s VERTEX:\n", string);
				    qh_printvertex (qh ferr, atvertex);
				  }
				  if (qh fout && qh FORCEoutput && atfacet && !qh QHULLfinished && !qh IStracing) {
				    fprintf(qh ferr, "ERRONEOUS and NEIGHBORING FACETS to output\n");
				    for (i= 0; i < qh_PRINTEND; i++)  /* use fout for geomview output */
				      qh_printneighborhood (qh fout, qh PRINTout[i], atfacet, otherfacet,
							    !qh_ALL);
				  }
				} /* errprint */
				
				
				/*-				  >--------------------------------
				  
				  qh_printfacetlist( fp, facetlist, facets, printall )
				    print all fields for a facet list and/or set of facets to fp
				    if !printall, 
				      only prints good facets
				
				  notes:
				    also prints all vertices
				*/
				void qh_printfacetlist(facetT *facetlist, setT *facets, boolT printall) {
				  facetT *facet, **facetp;
				
				  qh_printbegin (qh ferr, qh_PRINTfacets, facetlist, facets, printall);
				  FORALLfacet_(facetlist)
				    qh_printafacet(qh ferr, qh_PRINTfacets, facet, printall);
				  FOREACHfacet_(facets)
				    qh_printafacet(qh ferr, qh_PRINTfacets, facet, printall);
				  qh_printend (qh ferr, qh_PRINTfacets, facetlist, facets, printall);
				} /* printfacetlist */
				
				
				/*-				  >--------------------------------
				  
				  qh_user_memsizes()
				    allocate up to 10 additional, quick allocation sizes
				
				  notes:
				    increase maximum number of allocations in qh_initqhull_mem()
				*/
				void qh_user_memsizes (void) {
				
				  /* qh_memsize (size); */
				} /* user_memsizes */
				
							

相关资源