global.c

Quick hull implementation

	  }else {  /* not a procedure because of qh_option (filename, NULL, NULL); */
	    char filename[500], *t= filename;

	    s++;
	    while (*s) {
	      if (t - filename >= sizeof (filename)-2) {
		fprintf (qh ferr, "qhull error: filename for 'TI' too long.\n");
		qh_errexit (qh_ERRinput, NULL, NULL);
	      }
	      if (isspace (*s))
		break;
	      *(t++)= *s++;
	    }
	    *t= '\0';
	    if (!freopen (filename, "r", stdin)) {
	      fprintf (qh ferr, "qhull error: could not open file \"%s\".", filename);
	      qh_errexit (qh_ERRinput, NULL, NULL);
	    }else {
	      qh_option ("TInput-file", NULL, NULL);
	      qh_option (filename, NULL, NULL);
	    }
	  }
	  break;
	case 'O':
	  if (s[0] != ' ' || s[1] == '\"' || isspace (s[1])) {
	    s++;
	    fprintf (qh ferr, "qhull warning: option 'TO' mistyped.\nUse 'TO', one space, file name, and space or end-of-line.\nThe file name may be enclosed in single quotes.\nDo not use double quotes.  Option 'FO' ignored.\n");
	  }else {  /* not a procedure because of qh_option (filename, NULL, NULL); */
	    char filename[500], *t= filename;
	    boolT isquote= False;

	    s++;
	    if (*s == '\'') {
	      isquote= True;
	      s++;
	    }
	    while (*s) {
	      if (t - filename >= sizeof (filename)-2) {
		fprintf (qh ferr, "qhull error: filename for 'TO' too long.\n");
		qh_errexit (qh_ERRinput, NULL, NULL);
	      }
	      if (isquote) {
		if (*s == '\'') {
		  s++;
		  isquote= False;
		  break;
		}
	      }else if (isspace (*s))
		break;
	      *(t++)= *s++;
	    }
	    *t= '\0';
	    if (isquote)
	      fprintf (qh ferr, "qhull error: missing end quote for option 'TO'.  Rest of line ignored.\n");
	    else if (!freopen (filename, "w", stdout)) {
	      fprintf (qh ferr, "qhull error: could not open file \"%s\".", filename);
	      qh_errexit (qh_ERRinput, NULL, NULL);
	    }else {
	      qh_option ("TOutput-file", NULL, NULL);
	      qh_option (filename, NULL, NULL);
	    }
	  }
	  break;
	case 'P':
	  if (!isdigit(*s))
	    fprintf (qh ferr, "qhull warning: missing point id for trace option 'TPn'.  Ignored\n");
	  else {
	    qh TRACEpoint= qh_strtol (s, &s);
            qh_option ("Trace-point", &qh TRACEpoint, NULL);
          }
	  break;
	case 'M':
	  if (!isdigit(*s))
	    fprintf (qh ferr, "qhull warning: missing merge id for trace option 'TMn'.  Ignored\n");
	  else {
	    qh TRACEmerge= qh_strtol (s, &s);
            qh_option ("Trace-merge", &qh TRACEmerge, NULL);
          }
	  break;
	case 'R':
	  if (!isdigit(*s))
	    fprintf (qh ferr, "qhull warning: missing rerun count for trace option 'TRn'.  Ignored\n");
	  else {
	    qh RERUN= qh_strtol (s, &s);
            qh_option ("TRerun", &qh RERUN, NULL);
          }
	  break;
	case 'V':
	  i= qh_strtol (s, &t);
	  if (s == t)
	    fprintf (qh ferr, "qhull warning: missing furthest point id for trace option 'TVn'.  Ignored\n");
	  else if (i < 0) {
	    qh STOPpoint= i - 1;
            qh_option ("TV-stop-before-point", &i, NULL);
	  }else {
	    qh STOPpoint= i + 1;
            qh_option ("TV-stop-after-point", &i, NULL);
          }
          s= t;
	  break;
	case 'W':
	  if (!isdigit(*s))
	    fprintf (qh ferr, "qhull warning: missing max width for trace option 'TWn'.  Ignored\n");
	  else {
 	    qh TRACEdist= (realT) qh_strtod (s, &s);
            qh_option ("TWide-trace", NULL, &qh TRACEdist);
          }
	  break;
	default:
	  s--;
	  fprintf (qh ferr, "qhull warning: unknown 'T' trace option %c, rest ignored\n", (int)s[0]);
	  while (*++s && !isspace(*s));
	  break;
	}
      }
      break;
    default:
      fprintf (qh ferr, "qhull warning: unknown flag %c (%x)\n", (int)s[-1],
	       (int)s[-1]);
      break;
    }
    if (s-1 == prev_s && *s && !isspace(*s)) {
      fprintf (qh ferr, "qhull warning: missing space after flag %c (%x); reserved for menu. Skipped.\n",
	       (int)*prev_s, (int)*prev_s);
      while (*s && !isspace(*s))
	s++;
    }
  }
  if (isgeom && !qh FORCEoutput && qh PRINTout[1])
    fprintf (qh ferr, "qhull warning: additional output formats are not compatible with Geomview\n");
  /* set derived values in qh_initqhull_globals */
} /* initflags */


/*-<a                             href="qh-globa.htm#TOC"
  >-------------------------------</a><a name="initqhull_buffers">-</a>

  qh_initqhull_buffers()
    initialize global memory buffers

  notes:
    must match qh_freebuffers()
*/
void qh_initqhull_buffers (void) {
  int k;

  qh TEMPsize= (qhmem.LASTsize - sizeof (setT))/SETelemsize;
  if (qh TEMPsize <= 0 || qh TEMPsize > qhmem.LASTsize)
    qh TEMPsize= 8;  /* e.g., if qh_NOmem */
  qh other_points= qh_setnew (qh TEMPsize);
  qh del_vertices= qh_setnew (qh TEMPsize);
  qh coplanarset= qh_setnew (qh TEMPsize);
  qh NEARzero= (realT *)qh_memalloc(qh hull_dim * sizeof(realT));
  qh lower_threshold= (realT *)qh_memalloc((qh input_dim+1) * sizeof(realT));
  qh upper_threshold= (realT *)qh_memalloc((qh input_dim+1) * sizeof(realT));
  qh lower_bound= (realT *)qh_memalloc((qh input_dim+1) * sizeof(realT));
  qh upper_bound= (realT *)qh_memalloc((qh input_dim+1) * sizeof(realT));
  for(k= qh input_dim+1; k--; ) {
    qh lower_threshold[k]= -REALmax;
    qh upper_threshold[k]= REALmax;
    qh lower_bound[k]= -REALmax;
    qh upper_bound[k]= REALmax;
  }
  qh gm_matrix= (coordT *)qh_memalloc((qh hull_dim+1) * qh hull_dim * sizeof(coordT));
  qh gm_row= (coordT **)qh_memalloc((qh hull_dim+1) * sizeof(coordT *));
} /* initqhull_buffers */

/*-<a                             href="qh-globa.htm#TOC"
  >-------------------------------</a><a name="initqhull_globals">-</a>

  qh_initqhull_globals( points, numpoints, dim, ismalloc )
    initialize globals
    if ismalloc
      points were malloc'd and qhull should free at end

  returns:
    sets qh.first_point, num_points, input_dim, hull_dim and others
    seeds random number generator (seed=1 if tracing)
    modifies qh.hull_dim if ((qh.DELAUNAY and qh.PROJECTdelaunay) or qh.PROJECTinput)
    adjust user flags as needed
    also checks DIM3 dependencies and constants

  notes:
    do not use qh_point() since an input transformation may move them elsewhere

  see:
    qh_initqhull_start() sets default values for non-zero globals

  design:
    initialize points array from input arguments
    test for qh.ZEROcentrum
      (i.e., use opposite vertex instead of cetrum for convexity testing)
    test for qh.PRINTgood (i.e., only print 'good' facets)
    initialize qh.CENTERtype, qh.normal_size,
      qh.center_size, qh.TRACEpoint/level,
    initialize and test random numbers
    check for conflicting print output options
*/
void qh_initqhull_globals (coordT *points, int numpoints, int dim, boolT ismalloc) {
  int seed, pointsneeded, extra= 0, i, randi, k;
  boolT printgeom= False, printmath= False, printcoplanar= False;
  realT randr;
  realT factorial;

  time_t timedata;

  trace0((qh ferr, "qh_initqhull_globals: for %s | %s\n", qh rbox_command,
      qh qhull_command));
  qh POINTSmalloc= ismalloc;
  qh first_point= points;
  qh num_points= numpoints;
  qh hull_dim= qh input_dim= dim;
  if (!qh NOpremerge && !qh MERGEexact && !qh PREmerge && qh JOGGLEmax > REALmax/2) {
    qh MERGING= True;
    if (qh hull_dim <= 4) {
      qh PREmerge= True;
      qh_option ("_pre-merge", NULL, NULL);
    }else {
      qh MERGEexact= True;
      qh_option ("Qxact_merge", NULL, NULL);
    }
  }else if (qh MERGEexact) 
    qh MERGING= True;
  if (!qh NOpremerge && qh JOGGLEmax > REALmax/2) {
#ifdef qh_NOmerge
    qh JOGGLEmax= 0.0;
#endif
  }
  if (qh TRIangulate && qh JOGGLEmax < REALmax/2 && qh PRINTprecision)
    fprintf(qh ferr, "qhull warning: joggle ('QJ') always produces simplicial output.  Triangulated output ('Qt') does nothing.\n");
  if (qh JOGGLEmax < REALmax/2 && qh DELAUNAY && !qh SCALEinput && !qh SCALElast) {
    qh SCALElast= True;
    qh_option ("Qbbound-last-qj", NULL, NULL);
  }
  if (qh MERGING && !qh POSTmerge && qh premerge_cos > REALmax/2
  && qh premerge_centrum == 0) {
    qh ZEROcentrum= True;
    qh ZEROall_ok= True;
    qh_option ("_zero-centrum", NULL, NULL);
  }
  if (qh JOGGLEmax < REALmax/2 && REALepsilon > 2e-8 && qh PRINTprecision)
    fprintf(qh ferr, "qhull warning: real epsilon, %2.2g, is probably too large for joggle ('QJn')\nRecompile with double precision reals (see user.h).\n",
          REALepsilon);
#ifdef qh_NOmerge
  if (qh MERGING) {
    fprintf (qh ferr, "qhull input error: merging not installed (qh_NOmerge + 'Qx', 'Cn' or 'An')\n");
    qh_errexit (qh_ERRinput, NULL, NULL);
  }
#endif
  if (!(qh PRINTgood || qh PRINTneighbors)) {
    if (qh KEEParea || qh KEEPminArea < REALmax/2 || qh KEEPmerge || qh DELAUNAY
	|| (!qh ONLYgood && (qh GOODvertex || qh GOODpoint))) {
      qh PRINTgood= True;
      qh_option ("Pgood", NULL, NULL);
    }
  }
  if (qh DELAUNAY && qh KEEPcoplanar && !qh KEEPinside) {
    qh KEEPinside= True;
    qh_option ("Qinterior-keep", NULL, NULL);
  }
  if (qh DELAUNAY && qh HALFspace) {
    fprintf (qh ferr, "qhull input error: can not use Delaunay ('d') or Voronoi ('v') with halfspace intersection ('H')\n");
    qh_errexit (qh_ERRinput, NULL, NULL);
  }
  if (!qh DELAUNAY && (qh UPPERdelaunay || qh ATinfinity)) {
    fprintf (qh ferr, "qhull input error: use upper-Delaunay ('Qu') or infinity-point ('Qz') with Delaunay ('d') or Voronoi ('v')\n");
    qh_errexit (qh_ERRinput, NULL, NULL);
  }
  if (qh UPPERdelaunay && qh ATinfinity) {
    fprintf (qh ferr, "qhull input error: can not use infinity-point ('Qz') with upper-Delaunay ('Qu')\n");
    qh_errexit (qh_ERRinput, NULL, NULL);
  }
  if (qh SCALElast && !qh DELAUNAY && qh PRINTprecision)
    fprintf (qh ferr, "qhull input warning: option 'Qbb' (scale-last-coordinate) is normally used with 'd' or 'v'\n");
  qh DOcheckmax= (!qh SKIPcheckmax && qh MERGING );
  qh KEEPnearinside= (qh DOcheckmax && !(qh KEEPinside && qh KEEPcoplanar) 
                          && !qh NOnearinside);
  if (qh MERGING)
    qh CENTERtype= qh_AScentrum;
  else if (qh VORONOI)
    qh CENTERtype= qh_ASvoronoi;
  if (qh TESTvneighbors && !qh MERGING) {
    fprintf(qh ferr, "qhull input error: test vertex neighbors ('Qv') needs a merge option\n");
    qh_errexit (qh_ERRinput, NULL ,NULL);
  }
  if (qh PROJECTinput || (qh DELAUNAY && qh PROJECTdelaunay)) {
    qh hull_dim -= qh PROJECTinput;
    if (qh DELAUNAY) {
      qh hull_dim++;
      extra= 1;
    }
  }
  if (qh hull_dim <= 1) {
    fprintf(qh ferr, "qhull error: dimension %d must be > 1\n", qh hull_dim);
    qh_errexit (qh_ERRinput, NULL, NULL);
  }
  for (k= 2, factorial=1.0; k < qh hull_dim; k++)
    factorial *= k;
  qh AREAfactor= 1.0 / factorial;
  trace2((qh ferr, "qh_initqhull_globals: initialize globals.  dim %d numpoints %d malloc? %d projected %d to hull_dim %d\n",
	dim, numpoints, ismalloc, qh PROJECTinput, qh hull_dim));
  qh normal_size= qh hull_dim * sizeof(coordT);
  qh center_size= qh normal_size - sizeof(coordT);
  pointsneeded= qh hull_dim+1;
  if (qh hull_dim > qh_DIMmergeVertex) {
    qh MERGEvertices= False;
    qh_option ("Q3-no-merge-vertices-dim-high", NULL, NULL);
  }
  if (qh GOODpoint)
    pointsneeded++;
#ifdef qh_NOtrace
  if (qh IStracing) {
    fprintf (qh ferr, "qhull input error: tracing is not installed (qh_NOtrace in user.h)");
    qh_errexit (qh_ERRqhull, NULL, NULL);
  }
#endif
  if (qh RERUN > 1) {
    qh TRACElastrun= qh IStracing; /* qh_build_withrestart duplicates next conditional */
    if (qh IStracing != -1)
      qh IStracing= 0;
  }else if (qh TRACEpoint != -1 || qh TRACEdist < REALmax/2 || qh TRACEmerge) {
    qh TRACElevel= (qh IStracing? qh IStracing : 3);
    qh IStracing= 0;
  }
  if (qh ROTATErandom == 0 || qh ROTATErandom == -1) {
    seed= time (&timedata);
    if (qh ROTATErandom  == -1) {
      seed= -seed;
      qh_option ("QRandom-seed", &seed, NULL );
    }else
      qh_option ("QRotate-random", &seed, NULL);
    qh ROTATErandom= seed;
  }
  seed= qh ROTATErandom;
  if (seed == INT_MIN)    /* default value */
    seed= 1;
  else if (seed < 0)
    seed= -seed;
  qh_RANDOMseed_(seed);
  randr= 0.0;
  for (i= 1000; i--; ) {
    randi= qh_RANDOMint;
    randr += randi;
    if (randi > qh_RANDOMmax) {
      fprintf (qh ferr, "\
qhull configuration error (qh_RANDOMmax in user.h):\n\
   random integer %d > qh_RANDOMmax (%.8g)\n",
	       randi, qh_RANDOMmax);
      qh_errexit (qh_ERRinput, NULL, NULL);
    }
  }
  qh_RANDOMseed_(seed);
  randr = randr/1000;
  if (randr < qh_RANDOMmax/10
  || randr > qh_RANDOMmax * 5)
    fprintf (qh ferr, "\
qhull configuration warning (qh_RANDOMmax in user.h):\n\
   average of 1000 random integers (%.2g) is much different than expected (%.2g).\n\
   Is qh_RANDOMmax (%.2g) wrong?\n",
	     randr, qh_RANDOMmax/2.0, qh_RANDOMmax);
  qh RANDOMa= 2.0 * qh RANDOMfactor/qh_RANDOMmax;
  qh RANDOMb= 1.0 - qh RANDOMfactor;
  if (qh_HASHfactor < 1.1) {
    fprintf(qh ferr, "qhull internal error (qh_initqhull_globals): qh_HASHfactor %d must be at least 1.1.  Qhull uses linear hash probing\n",
      qh_HASHfactor);
    qh_errexit (qh_ERRqhull, NULL, NULL);
  }
  if (numpoints+extra < pointsneeded) {
    fprintf(qh ferr,"qhull input error: not enough points (%d) to construct initial simplex (need %d)\n",
	    numpoints, pointsneeded);
    qh_errexit(qh_ERRinput, NULL, NULL);
  }
  if (qh PRINTtransparent) {
    if (qh hull_dim != 4 || !qh DELAUNAY || qh VORONOI || qh DROPdim >= 0) {
      fprintf(qh ferr,"qhull input error: transparent Delaunay ('Gt') needs 3-d Delaunay ('d') w/o 'GDn'\n");
      qh_errexit(qh_ERRinput, NULL, NULL);
    }
    qh DROPdim = 3;
    qh PRINTridges = True;
  }
  for (i= qh_PRINTEND; i--; ) {
    if (qh PRINTout[i] == qh_PRINTgeom)
      printgeom= True;
    else if (qh PRINTout[i] == qh_PRINTmathematica || qh PRINTout[i] == qh_PRINTmaple)
      printmath= True;
    else if (qh PRINTout[i] == qh_PRINTcoplanars)
      printcoplanar= True;
    else if (qh PRINTout[i] == qh_PRINTpointnearest)
      printcoplanar= True;
    else if (qh PRINTout[i] == qh_PRINTpointintersect && !qh HALFspace) {
      fprintf (qh ferr, "qhull input error: option 'Fp' is only used for \nhalfspace intersection ('Hn,n,n').\n");
      qh_errexit (qh_ERRinput, NULL, NULL);
    }else if (qh PRINTout[i] == qh_PRINTtriangles && (qh HALFspace || qh VORONOI)) {
      fprintf (qh ferr, "qhull input error: option 'Ft' is not available for Voronoi vertices or halfspace intersection\n");
      qh_errexit (qh_ERRinput, NULL, NULL);
    }else if (qh PRINTout[i] == qh_PRINTcentrums && qh VORONOI) {
      fprintf (qh ferr, "qhull input error: option 'FC' is not available for Voronoi vertices ('v')\n");
      qh_errexit (qh_ERRinput, NULL, NULL);
    }else if (qh PRINTout[i] == qh_PRINTvertices) {
      if (qh VORONOI)
        qh_option ("Fvoronoi", NULL, NULL);
      else 
        qh_option ("Fvertices", NULL, NULL);
    }
  }
  if (printcoplanar && qh DELAUNAY && qh JOGGLEmax < REALmax/2) {
    if (qh PRINTprecision) 
      fprintf (qh ferr, "qhull input warning: 'QJ' (joggle) will usually prevent coincident input sites for options 'Fc' and 'FP'\n");
  }