geometry/src/geom2_r.c

/*<html><pre>  -<a                             href="qh-geom_r.htm"
  >-------------------------------</a><a name="TOP">-</a>


   geom2_r.c
   infrequently used geometric routines of qhull

   see qh-geom_r.htm and geom_r.h

   Copyright (c) 1993-2015 The Geometry Center.
   $Id: //main/2015/qhull/src/libqhull_r/geom2_r.c#6 $$Change: 2065 $
   $DateTime: 2016/01/18 13:51:04 $$Author: bbarber $

   frequently used code goes into geom_r.c
*/

#include "qhull_ra.h"

/*================== functions in alphabetic order ============*/

/*-<a                             href="qh-geom_r.htm#TOC"
  >-------------------------------</a><a name="copypoints">-</a>

  qh_copypoints(qh, points, numpoints, dimension)
    return qh_malloc'd copy of points

  notes:
    qh_free the returned points to avoid a memory leak
*/
coordT *qh_copypoints(qhT *qh, coordT *points, int numpoints, int dimension) {
  int size;
  coordT *newpoints;

  size= numpoints * dimension * (int)sizeof(coordT);
  if (!(newpoints= (coordT*)qh_malloc((size_t)size))) {
    qh_fprintf(qh, qh->ferr, 6004, "qhull error: insufficient memory to copy %d points\n",
        numpoints);
    qh_errexit(qh, qh_ERRmem, NULL, NULL);
  }
  memcpy((char *)newpoints, (char *)points, (size_t)size); /* newpoints!=0 by QH6004 */
  return newpoints;
} /* copypoints */

/*-<a                             href="qh-geom_r.htm#TOC"
  >-------------------------------</a><a name="crossproduct">-</a>

  qh_crossproduct( dim, vecA, vecB, vecC )
    crossproduct of 2 dim vectors
    C= A x B

  notes:
    from Glasner, Graphics Gems I, p. 639
    only defined for dim==3
*/
void qh_crossproduct(int dim, realT vecA[3], realT vecB[3], realT vecC[3]){

  if (dim == 3) {
    vecC[0]=   det2_(vecA[1], vecA[2],
                     vecB[1], vecB[2]);
    vecC[1]= - det2_(vecA[0], vecA[2],
                     vecB[0], vecB[2]);
    vecC[2]=   det2_(vecA[0], vecA[1],
                     vecB[0], vecB[1]);
  }
} /* vcross */

/*-<a                             href="qh-geom_r.htm#TOC"
  >-------------------------------</a><a name="determinant">-</a>

  qh_determinant(qh, rows, dim, nearzero )
    compute signed determinant of a square matrix
    uses qh.NEARzero to test for degenerate matrices

  returns:
    determinant
    overwrites rows and the matrix
    if dim == 2 or 3
      nearzero iff determinant < qh->NEARzero[dim-1]
      (!quite correct, not critical)
    if dim >= 4
      nearzero iff diagonal[k] < qh->NEARzero[k]
*/
realT qh_determinant(qhT *qh, realT **rows, int dim, boolT *nearzero) {
  realT det=0;
  int i;
  boolT sign= False;

  *nearzero= False;
  if (dim < 2) {
    qh_fprintf(qh, qh->ferr, 6005, "qhull internal error (qh_determinate): only implemented for dimension >= 2\n");
    qh_errexit(qh, qh_ERRqhull, NULL, NULL);
  }else if (dim == 2) {
    det= det2_(rows[0][0], rows[0][1],
                 rows[1][0], rows[1][1]);
    if (fabs_(det) < 10*qh->NEARzero[1])  /* not really correct, what should this be? */
      *nearzero= True;
  }else if (dim == 3) {
    det= det3_(rows[0][0], rows[0][1], rows[0][2],
                 rows[1][0], rows[1][1], rows[1][2],
                 rows[2][0], rows[2][1], rows[2][2]);
    if (fabs_(det) < 10*qh->NEARzero[2])  /* what should this be?  det 5.5e-12 was flat for qh_maxsimplex of qdelaunay 0,0 27,27 -36,36 -9,63 */
      *nearzero= True;
  }else {
    qh_gausselim(qh, rows, dim, dim, &sign, nearzero);  /* if nearzero, diagonal still ok*/
    det= 1.0;
    for (i=dim; i--; )
      det *= (rows[i])[i];
    if (sign)
      det= -det;
  }
  return det;
} /* determinant */

/*-<a                             href="qh-geom_r.htm#TOC"
  >-------------------------------</a><a name="detjoggle">-</a>

  qh_detjoggle(qh, points, numpoints, dimension )
    determine default max joggle for point array
      as qh_distround * qh_JOGGLEdefault

  returns:
    initial value for JOGGLEmax from points and REALepsilon

  notes:
    computes DISTround since qh_maxmin not called yet
    if qh->SCALElast, last dimension will be scaled later to MAXwidth

    loop duplicated from qh_maxmin
*/
realT qh_detjoggle(qhT *qh, pointT *points, int numpoints, int dimension) {
  realT abscoord, distround, joggle, maxcoord, mincoord;
  pointT *point, *pointtemp;
  realT maxabs= -REALmax;
  realT sumabs= 0;
  realT maxwidth= 0;
  int k;

  for (k=0; k < dimension; k++) {
    if (qh->SCALElast && k == dimension-1)
      abscoord= maxwidth;
    else if (qh->DELAUNAY && k == dimension-1) /* will qh_setdelaunay() */
      abscoord= 2 * maxabs * maxabs;  /* may be low by qh->hull_dim/2 */
    else {
      maxcoord= -REALmax;
      mincoord= REALmax;
      FORALLpoint_(qh, points, numpoints) {
        maximize_(maxcoord, point[k]);
        minimize_(mincoord, point[k]);
      }
      maximize_(maxwidth, maxcoord-mincoord);
      abscoord= fmax_(maxcoord, -mincoord);
    }
    sumabs += abscoord;
    maximize_(maxabs, abscoord);
  } /* for k */
  distround= qh_distround(qh, qh->hull_dim, maxabs, sumabs);
  joggle= distround * qh_JOGGLEdefault;
  maximize_(joggle, REALepsilon * qh_JOGGLEdefault);
  trace2((qh, qh->ferr, 2001, "qh_detjoggle: joggle=%2.2g maxwidth=%2.2g\n", joggle, maxwidth));
  return joggle;
} /* detjoggle */

/*-<a                             href="qh-geom_r.htm#TOC"
  >-------------------------------</a><a name="detroundoff">-</a>

  qh_detroundoff(qh)
    determine maximum roundoff errors from
      REALepsilon, REALmax, REALmin, qh.hull_dim, qh.MAXabs_coord,
      qh.MAXsumcoord, qh.MAXwidth, qh.MINdenom_1

    accounts for qh.SETroundoff, qh.RANDOMdist, qh->MERGEexact
      qh.premerge_cos, qh.postmerge_cos, qh.premerge_centrum,
      qh.postmerge_centrum, qh.MINoutside,
      qh_RATIOnearinside, qh_COPLANARratio, qh_WIDEcoplanar

  returns:
    sets qh.DISTround, etc. (see below)
    appends precision constants to qh.qhull_options

  see:
    qh_maxmin() for qh.NEARzero

  design:
    determine qh.DISTround for distance computations
    determine minimum denominators for qh_divzero
    determine qh.ANGLEround for angle computations
    adjust qh.premerge_cos,... for roundoff error
    determine qh.ONEmerge for maximum error due to a single merge
    determine qh.NEARinside, qh.MAXcoplanar, qh.MINvisible,
      qh.MINoutside, qh.WIDEfacet
    initialize qh.max_vertex and qh.minvertex
*/
void qh_detroundoff(qhT *qh) {

  qh_option(qh, "_max-width", NULL, &qh->MAXwidth);
  if (!qh->SETroundoff) {
    qh->DISTround= qh_distround(qh, qh->hull_dim, qh->MAXabs_coord, qh->MAXsumcoord);
    if (qh->RANDOMdist)
      qh->DISTround += qh->RANDOMfactor * qh->MAXabs_coord;
    qh_option(qh, "Error-roundoff", NULL, &qh->DISTround);
  }
  qh->MINdenom= qh->MINdenom_1 * qh->MAXabs_coord;
  qh->MINdenom_1_2= sqrt(qh->MINdenom_1 * qh->hull_dim) ;  /* if will be normalized */
  qh->MINdenom_2= qh->MINdenom_1_2 * qh->MAXabs_coord;
                                              /* for inner product */
  qh->ANGLEround= 1.01 * qh->hull_dim * REALepsilon;
  if (qh->RANDOMdist)
    qh->ANGLEround += qh->RANDOMfactor;
  if (qh->premerge_cos < REALmax/2) {
    qh->premerge_cos -= qh->ANGLEround;
    if (qh->RANDOMdist)
      qh_option(qh, "Angle-premerge-with-random", NULL, &qh->premerge_cos);
  }
  if (qh->postmerge_cos < REALmax/2) {
    qh->postmerge_cos -= qh->ANGLEround;
    if (qh->RANDOMdist)
      qh_option(qh, "Angle-postmerge-with-random", NULL, &qh->postmerge_cos);
  }
  qh->premerge_centrum += 2 * qh->DISTround;    /*2 for centrum and distplane()*/
  qh->postmerge_centrum += 2 * qh->DISTround;
  if (qh->RANDOMdist && (qh->MERGEexact || qh->PREmerge))
    qh_option(qh, "Centrum-premerge-with-random", NULL, &qh->premerge_centrum);
  if (qh->RANDOMdist && qh->POSTmerge)
    qh_option(qh, "Centrum-postmerge-with-random", NULL, &qh->postmerge_centrum);
  { /* compute ONEmerge, max vertex offset for merging simplicial facets */
    realT maxangle= 1.0, maxrho;

    minimize_(maxangle, qh->premerge_cos);
    minimize_(maxangle, qh->postmerge_cos);
    /* max diameter * sin theta + DISTround for vertex to its hyperplane */
    qh->ONEmerge= sqrt((realT)qh->hull_dim) * qh->MAXwidth *
      sqrt(1.0 - maxangle * maxangle) + qh->DISTround;
    maxrho= qh->hull_dim * qh->premerge_centrum + qh->DISTround;
    maximize_(qh->ONEmerge, maxrho);
    maxrho= qh->hull_dim * qh->postmerge_centrum + qh->DISTround;
    maximize_(qh->ONEmerge, maxrho);
    if (qh->MERGING)
      qh_option(qh, "_one-merge", NULL, &qh->ONEmerge);
  }
  qh->NEARinside= qh->ONEmerge * qh_RATIOnearinside; /* only used if qh->KEEPnearinside */
  if (qh->JOGGLEmax < REALmax/2 && (qh->KEEPcoplanar || qh->KEEPinside)) {
    realT maxdist;             /* adjust qh.NEARinside for joggle */
    qh->KEEPnearinside= True;
    maxdist= sqrt((realT)qh->hull_dim) * qh->JOGGLEmax + qh->DISTround;
    maxdist= 2*maxdist;        /* vertex and coplanar point can joggle in opposite directions */
    maximize_(qh->NEARinside, maxdist);  /* must agree with qh_nearcoplanar() */
  }
  if (qh->KEEPnearinside)
    qh_option(qh, "_near-inside", NULL, &qh->NEARinside);
  if (qh->JOGGLEmax < qh->DISTround) {
    qh_fprintf(qh, qh->ferr, 6006, "qhull error: the joggle for 'QJn', %.2g, is below roundoff for distance computations, %.2g\n",
         qh->JOGGLEmax, qh->DISTround);
    qh_errexit(qh, qh_ERRinput, NULL, NULL);
  }
  if (qh->MINvisible > REALmax/2) {
    if (!qh->MERGING)
      qh->MINvisible= qh->DISTround;
    else if (qh->hull_dim <= 3)
      qh->MINvisible= qh->premerge_centrum;
    else
      qh->MINvisible= qh_COPLANARratio * qh->premerge_centrum;
    if (qh->APPROXhull && qh->MINvisible > qh->MINoutside)
      qh->MINvisible= qh->MINoutside;
    qh_option(qh, "Visible-distance", NULL, &qh->MINvisible);
  }
  if (qh->MAXcoplanar > REALmax/2) {
    qh->MAXcoplanar= qh->MINvisible;
    qh_option(qh, "U-coplanar-distance", NULL, &qh->MAXcoplanar);
  }
  if (!qh->APPROXhull) {             /* user may specify qh->MINoutside */
    qh->MINoutside= 2 * qh->MINvisible;
    if (qh->premerge_cos < REALmax/2)
      maximize_(qh->MINoutside, (1- qh->premerge_cos) * qh->MAXabs_coord);
    qh_option(qh, "Width-outside", NULL, &qh->MINoutside);
  }
  qh->WIDEfacet= qh->MINoutside;
  maximize_(qh->WIDEfacet, qh_WIDEcoplanar * qh->MAXcoplanar);
  maximize_(qh->WIDEfacet, qh_WIDEcoplanar * qh->MINvisible);
  qh_option(qh, "_wide-facet", NULL, &qh->WIDEfacet);
  if (qh->MINvisible > qh->MINoutside + 3 * REALepsilon
  && !qh->BESToutside && !qh->FORCEoutput)
    qh_fprintf(qh, qh->ferr, 7001, "qhull input warning: minimum visibility V%.2g is greater than \nminimum outside W%.2g.  Flipped facets are likely.\n",
             qh->MINvisible, qh->MINoutside);
  qh->max_vertex= qh->DISTround;
  qh->min_vertex= -qh->DISTround;
  /* numeric constants reported in printsummary */
} /* detroundoff */

/*-<a                             href="qh-geom_r.htm#TOC"
  >-------------------------------</a><a name="detsimplex">-</a>

  qh_detsimplex(qh, apex, points, dim, nearzero )
    compute determinant of a simplex with point apex and base points

  returns:
     signed determinant and nearzero from qh_determinant

  notes:
     uses qh.gm_matrix/qh.gm_row (assumes they're big enough)

  design:
    construct qm_matrix by subtracting apex from points
    compute determinate
*/
realT qh_detsimplex(qhT *qh, pointT *apex, setT *points, int dim, boolT *nearzero) {
  pointT *coorda, *coordp, *gmcoord, *point, **pointp;
  coordT **rows;
  int k,  i=0;
  realT det;

  zinc_(Zdetsimplex);
  gmcoord= qh->gm_matrix;
  rows= qh->gm_row;
  FOREACHpoint_(points) {
    if (i == dim)
      break;
    rows[i++]= gmcoord;
    coordp= point;
    coorda= apex;
    for (k=dim; k--; )
      *(gmcoord++)= *coordp++ - *coorda++;
  }
  if (i < dim) {
    qh_fprintf(qh, qh->ferr, 6007, "qhull internal error (qh_detsimplex): #points %d < dimension %d\n",
               i, dim);
    qh_errexit(qh, qh_ERRqhull, NULL, NULL);
  }
  det= qh_determinant(qh, rows, dim, nearzero);
  trace2((qh, qh->ferr, 2002, "qh_detsimplex: det=%2.2g for point p%d, dim %d, nearzero? %d\n",
          det, qh_pointid(qh, apex), dim, *nearzero));
  return det;
} /* detsimplex */

/*-<a                             href="qh-geom_r.htm#TOC"
  >-------------------------------</a><a name="distnorm">-</a>

  qh_distnorm( dim, point, normal, offset )
    return distance from point to hyperplane at normal/offset

  returns:
    dist

  notes:
    dist > 0 if point is outside of hyperplane

  see:
    qh_distplane in geom_r.c
*/
realT qh_distnorm(int dim, pointT *point, pointT *normal, realT *offsetp) {
  coordT *normalp= normal, *coordp= point;
  realT dist;
  int k;

  dist= *offsetp;
  for (k=dim; k--; )
    dist += *(coordp++) * *(normalp++);
  return dist;
} /* distnorm */

/*-<a                             href="qh-geom_r.htm#TOC"
  >-------------------------------</a><a name="distround">-</a>

  qh_distround(qh, dimension, maxabs, maxsumabs )
    compute maximum round-off error for a distance computation
      to a normalized hyperplane
    maxabs is the maximum absolute value of a coordinate
    maxsumabs is the maximum possible sum of absolute coordinate values

  returns:
    max dist round for REALepsilon

  notes:
    calculate roundoff error according to Golub & van Loan, 1983, Lemma 3.2-1, "Rounding Errors"
    use sqrt(dim) since one vector is normalized
      or use maxsumabs since one vector is < 1
*/
realT qh_distround(qhT *qh, int dimension, realT maxabs, realT maxsumabs) {
  realT maxdistsum, maxround;

  maxdistsum= sqrt((realT)dimension) * maxabs;
  minimize_( maxdistsum, maxsumabs);
  maxround= REALepsilon * (dimension * maxdistsum * 1.01 + maxabs);
              /* adds maxabs for offset */
  trace4((qh, qh->ferr, 4008, "qh_distround: %2.2g maxabs %2.2g maxsumabs %2.2g maxdistsum %2.2g\n",
                 maxround, maxabs, maxsumabs, maxdistsum));
  return maxround;
} /* distround */

/*-<a                             href="qh-geom_r.htm#TOC"
  >-------------------------------</a><a name="divzero">-</a>

  qh_divzero( numer, denom, mindenom1, zerodiv )
    divide by a number that's nearly zero
    mindenom1= minimum denominator for dividing into 1.0

  returns:
    quotient
    sets zerodiv and returns 0.0 if it would overflow

  design:
    if numer is nearly zero and abs(numer) < abs(denom)
      return numer/denom
    else if numer is nearly zero
      return 0 and zerodiv
    else if denom/numer non-zero
      return numer/denom
    else
      return 0 and zerodiv
*/
realT qh_divzero(realT numer, realT denom, realT mindenom1, boolT *zerodiv) {
  realT temp, numerx, denomx;


  if (numer < mindenom1 && numer > -mindenom1) {
    numerx= fabs_(numer);
    denomx= fabs_(denom);
    if (numerx < denomx) {
      *zerodiv= False;
      return numer/denom;
    }else {
      *zerodiv= True;
      return 0.0;
    }
  }
  temp= denom/numer;
  if (temp > mindenom1 || temp < -mindenom1) {
    *zerodiv= False;
    return numer/denom;
  }else {
    *zerodiv= True;
    return 0.0;
  }
} /* divzero */


/*-<a                             href="qh-geom_r.htm#TOC"
  >-------------------------------</a><a name="facetarea">-</a>

  qh_facetarea(qh, facet )
    return area for a facet

  notes:
    if non-simplicial,
      uses centrum to triangulate facet and sums the projected areas.
    if (qh->DELAUNAY),
      computes projected area instead for last coordinate
    assumes facet->normal exists
    projecting tricoplanar facets to the hyperplane does not appear to make a difference

  design:
    if simplicial
      compute area
    else
      for each ridge
        compute area from centrum to ridge
    negate area if upper Delaunay facet
*/
realT qh_facetarea(qhT *qh, facetT *facet) {
  vertexT *apex;
  pointT *centrum;
  realT area= 0.0;
  ridgeT *ridge, **ridgep;

  if (facet->simplicial) {
    apex= SETfirstt_(facet->vertices, vertexT);
    area= qh_facetarea_simplex(qh, qh->hull_dim, apex->point, facet->vertices,
                    apex, facet->toporient, facet->normal, &facet->offset);
  }else {
    if (qh->CENTERtype == qh_AScentrum)
      centrum= facet->center;
    else
      centrum= qh_getcentrum(qh, facet);
    FOREACHridge_(facet->ridges)
      area += qh_facetarea_simplex(qh, qh->hull_dim, centrum, ridge->vertices,
                 NULL, (boolT)(ridge->top == facet),  facet->normal, &facet->offset);
    if (qh->CENTERtype != qh_AScentrum)
      qh_memfree(qh, centrum, qh->normal_size);
  }
  if (facet->upperdelaunay && qh->DELAUNAY)
    area= -area;  /* the normal should be [0,...,1] */
  trace4((qh, qh->ferr, 4009, "qh_facetarea: f%d area %2.2g\n", facet->id, area));
  return area;
} /* facetarea */

/*-<a                             href="qh-geom_r.htm#TOC"
  >-------------------------------</a><a name="facetarea_simplex">-</a>

  qh_facetarea_simplex(qh, dim, apex, vertices, notvertex, toporient, normal, offset )
    return area for a simplex defined by
      an apex, a base of vertices, an orientation, and a unit normal
    if simplicial or tricoplanar facet,
      notvertex is defined and it is skipped in vertices

  returns:
    computes area of simplex projected to plane [normal,offset]
    returns 0 if vertex too far below plane (qh->WIDEfacet)
      vertex can't be apex of tricoplanar facet

  notes:
    if (qh->DELAUNAY),
      computes projected area instead for last coordinate
    uses qh->gm_matrix/gm_row and qh->hull_dim
    helper function for qh_facetarea

  design:
    if Notvertex
      translate simplex to apex
    else
      project simplex to normal/offset
      translate simplex to apex
    if Delaunay
      set last row/column to 0 with -1 on diagonal
    else
      set last row to Normal
    compute determinate
    scale and flip sign for area
*/
realT qh_facetarea_simplex(qhT *qh, int dim, coordT *apex, setT *vertices,
        vertexT *notvertex,  boolT toporient, coordT *normal, realT *offset) {
  pointT *coorda, *coordp, *gmcoord;
  coordT **rows, *normalp;
  int k,  i=0;
  realT area, dist;
  vertexT *vertex, **vertexp;
  boolT nearzero;

  gmcoord= qh->gm_matrix;
  rows= qh->gm_row;
  FOREACHvertex_(vertices) {
    if (vertex == notvertex)
      continue;
    rows[i++]= gmcoord;
    coorda= apex;
    coordp= vertex->point;
    normalp= normal;
    if (notvertex) {
      for (k=dim; k--; )
        *(gmcoord++)= *coordp++ - *coorda++;
    }else {
      dist= *offset;
      for (k=dim; k--; )
        dist += *coordp++ * *normalp++;
      if (dist < -qh->WIDEfacet) {
        zinc_(Znoarea);
        return 0.0;
      }
      coordp= vertex->point;
      normalp= normal;
      for (k=dim; k--; )
        *(gmcoord++)= (*coordp++ - dist * *normalp++) - *coorda++;
    }
  }
  if (i != dim-1) {
    qh_fprintf(qh, qh->ferr, 6008, "qhull internal error (qh_facetarea_simplex): #points %d != dim %d -1\n",
               i, dim);
    qh_errexit(qh, qh_ERRqhull, NULL, NULL);
  }
  rows[i]= gmcoord;
  if (qh->DELAUNAY) {
    for (i=0; i < dim-1; i++)
      rows[i][dim-1]= 0.0;
    for (k=dim; k--; )
      *(gmcoord++)= 0.0;
    rows[dim-1][dim-1]= -1.0;
  }else {
    normalp= normal;
    for (k=dim; k--; )
      *(gmcoord++)= *normalp++;
  }
  zinc_(Zdetsimplex);
  area= qh_determinant(qh, rows, dim, &nearzero);
  if (toporient)
    area= -area;
  area *= qh->AREAfactor;
  trace4((qh, qh->ferr, 4010, "qh_facetarea_simplex: area=%2.2g for point p%d, toporient %d, nearzero? %d\n",
          area, qh_pointid(qh, apex), toporient, nearzero));
  return area;
} /* facetarea_simplex */

/*-<a                             href="qh-geom_r.htm#TOC"
  >-------------------------------</a><a name="facetcenter">-</a>

  qh_facetcenter(qh, vertices )
    return Voronoi center (Voronoi vertex) for a facet's vertices

  returns:
    return temporary point equal to the center

  see:
    qh_voronoi_center()
*/
pointT *qh_facetcenter(qhT *qh, setT *vertices) {
  setT *points= qh_settemp(qh, qh_setsize(qh, vertices));
  vertexT *vertex, **vertexp;
  pointT *center;

  FOREACHvertex_(vertices)
    qh_setappend(qh, &points, vertex->point);
  center= qh_voronoi_center(qh, qh->hull_dim-1, points);
  qh_settempfree(qh, &points);
  return center;
} /* facetcenter */

/*-<a                             href="qh-geom_r.htm#TOC"
  >-------------------------------</a><a name="findgooddist">-</a>

  qh_findgooddist(qh, point, facetA, dist, facetlist )
    find best good facet visible for point from facetA
    assumes facetA is visible from point

  returns:
    best facet, i.e., good facet that is furthest from point
      distance to best facet
      NULL if none

    moves good, visible facets (and some other visible facets)
      to end of qh->facet_list

  notes:
    uses qh->visit_id

  design:
    initialize bestfacet if facetA is good
    move facetA to end of facetlist
    for each facet on facetlist
      for each unvisited neighbor of facet
        move visible neighbors to end of facetlist
        update best good neighbor
        if no good neighbors, update best facet
*/
facetT *qh_findgooddist(qhT *qh, pointT *point, facetT *facetA, realT *distp,
               facetT **facetlist) {
  realT bestdist= -REALmax, dist;
  facetT *neighbor, **neighborp, *bestfacet=NULL, *facet;
  boolT goodseen= False;

  if (facetA->good) {
    zzinc_(Zcheckpart);  /* calls from check_bestdist occur after print stats */
    qh_distplane(qh, point, facetA, &bestdist);
    bestfacet= facetA;
    goodseen= True;
  }
  qh_removefacet(qh, facetA);
  qh_appendfacet(qh, facetA);
  *facetlist= facetA;
  facetA->visitid= ++qh->visit_id;
  FORALLfacet_(*facetlist) {
    FOREACHneighbor_(facet) {
      if (neighbor->visitid == qh->visit_id)
        continue;
      neighbor->visitid= qh->visit_id;
      if (goodseen && !neighbor->good)
        continue;
      zzinc_(Zcheckpart);
      qh_distplane(qh, point, neighbor, &dist);
      if (dist > 0) {
        qh_removefacet(qh, neighbor);
        qh_appendfacet(qh, neighbor);
        if (neighbor->good) {
          goodseen= True;
          if (dist > bestdist) {
            bestdist= dist;
            bestfacet= neighbor;
          }
        }
      }
    }
  }
  if (bestfacet) {
    *distp= bestdist;
    trace2((qh, qh->ferr, 2003, "qh_findgooddist: p%d is %2.2g above good facet f%d\n",
      qh_pointid(qh, point), bestdist, bestfacet->id));
    return bestfacet;
  }
  trace4((qh, qh->ferr, 4011, "qh_findgooddist: no good facet for p%d above f%d\n",
      qh_pointid(qh, point), facetA->id));
  return NULL;
}  /* findgooddist */

/*-<a                             href="qh-geom_r.htm#TOC"
  >-------------------------------</a><a name="getarea">-</a>

  qh_getarea(qh, facetlist )
    set area of all facets in facetlist
    collect statistics
    nop if hasAreaVolume

  returns:
    sets qh->totarea/totvol to total area and volume of convex hull
    for Delaunay triangulation, computes projected area of the lower or upper hull
      ignores upper hull if qh->ATinfinity

  notes:
    could compute outer volume by expanding facet area by rays from interior
    the following attempt at perpendicular projection underestimated badly:
      qh.totoutvol += (-dist + facet->maxoutside + qh->DISTround)
                            * area/ qh->hull_dim;
  design:
    for each facet on facetlist
      compute facet->area
      update qh.totarea and qh.totvol
*/
void qh_getarea(qhT *qh, facetT *facetlist) {
  realT area;
  realT dist;
  facetT *facet;

  if (qh->hasAreaVolume)
    return;
  if (qh->REPORTfreq)
    qh_fprintf(qh, qh->ferr, 8020, "computing area of each facet and volume of the convex hull\n");
  else
    trace1((qh, qh->ferr, 1001, "qh_getarea: computing volume and area for each facet\n"));
  qh->totarea= qh->totvol= 0.0;
  FORALLfacet_(facetlist) {
    if (!facet->normal)
      continue;
    if (facet->upperdelaunay && qh->ATinfinity)
      continue;
    if (!facet->isarea) {
      facet->f.area= qh_facetarea(qh, facet);
      facet->isarea= True;
    }
    area= facet->f.area;
    if (qh->DELAUNAY) {
      if (facet->upperdelaunay == qh->UPPERdelaunay)
        qh->totarea += area;
    }else {
      qh->totarea += area;
      qh_distplane(qh, qh->interior_point, facet, &dist);
      qh->totvol += -dist * area/ qh->hull_dim;
    }
    if (qh->PRINTstatistics) {
      wadd_(Wareatot, area);
      wmax_(Wareamax, area);
      wmin_(Wareamin, area);
    }
  }
  qh->hasAreaVolume= True;
} /* getarea */

/*-<a                             href="qh-geom_r.htm#TOC"
  >-------------------------------</a><a name="gram_schmidt">-</a>

  qh_gram_schmidt(qh, dim, row )
    implements Gram-Schmidt orthogonalization by rows

  returns:
    false if zero norm
    overwrites rows[dim][dim]

  notes:
    see Golub & van Loan, 1983, Algorithm 6.2-2, "Modified Gram-Schmidt"
    overflow due to small divisors not handled

  design:
    for each row
      compute norm for row
      if non-zero, normalize row
      for each remaining rowA
        compute inner product of row and rowA
        reduce rowA by row * inner product
*/
boolT qh_gram_schmidt(qhT *qh, int dim, realT **row) {
  realT *rowi, *rowj, norm;
  int i, j, k;

  for (i=0; i < dim; i++) {
    rowi= row[i];
    for (norm= 0.0, k= dim; k--; rowi++)
      norm += *rowi * *rowi;
    norm= sqrt(norm);
    wmin_(Wmindenom, norm);
    if (norm == 0.0)  /* either 0 or overflow due to sqrt */
      return False;
    for (k=dim; k--; )
      *(--rowi) /= norm;
    for (j=i+1; j < dim; j++) {
      rowj= row[j];
      for (norm= 0.0, k=dim; k--; )
        norm += *rowi++ * *rowj++;
      for (k=dim; k--; )
        *(--rowj) -= *(--rowi) * norm;
    }
  }
  return True;
} /* gram_schmidt */


/*-<a                             href="qh-geom_r.htm#TOC"
  >-------------------------------</a><a name="inthresholds">-</a>

  qh_inthresholds(qh, normal, angle )
    return True if normal within qh.lower_/upper_threshold

  returns:
    estimate of angle by summing of threshold diffs
      angle may be NULL
      smaller "angle" is better

  notes:
    invalid if qh.SPLITthresholds

  see:
    qh.lower_threshold in qh_initbuild()
    qh_initthresholds()

  design:
    for each dimension
      test threshold
*/
boolT qh_inthresholds(qhT *qh, coordT *normal, realT *angle) {
  boolT within= True;
  int k;
  realT threshold;

  if (angle)
    *angle= 0.0;
  for (k=0; k < qh->hull_dim; k++) {
    threshold= qh->lower_threshold[k];
    if (threshold > -REALmax/2) {
      if (normal[k] < threshold)
        within= False;
      if (angle) {
        threshold -= normal[k];
        *angle += fabs_(threshold);
      }
    }
    if (qh->upper_threshold[k] < REALmax/2) {
      threshold= qh->upper_threshold[k];
      if (normal[k] > threshold)
        within= False;
      if (angle) {
        threshold -= normal[k];
        *angle += fabs_(threshold);
      }
    }
  }
  return within;
} /* inthresholds */


/*-<a                             href="qh-geom_r.htm#TOC"
  >-------------------------------</a><a name="joggleinput">-</a>

  qh_joggleinput(qh)
    randomly joggle input to Qhull by qh.JOGGLEmax
    initial input is qh.first_point/qh.num_points of qh.hull_dim
      repeated calls use qh.input_points/qh.num_points

  returns:
    joggles points at qh.first_point/qh.num_points
    copies data to qh.input_points/qh.input_malloc if first time
    determines qh.JOGGLEmax if it was zero
    if qh.DELAUNAY
      computes the Delaunay projection of the joggled points

  notes:
    if qh.DELAUNAY, unnecessarily joggles the last coordinate
    the initial 'QJn' may be set larger than qh_JOGGLEmaxincrease

  design:
    if qh.DELAUNAY
      set qh.SCALElast for reduced precision errors
    if first call
      initialize qh.input_points to the original input points
      if qh.JOGGLEmax == 0
        determine default qh.JOGGLEmax
    else
      increase qh.JOGGLEmax according to qh.build_cnt
    joggle the input by adding a random number in [-qh.JOGGLEmax,qh.JOGGLEmax]
    if qh.DELAUNAY
      sets the Delaunay projection
*/
void qh_joggleinput(qhT *qh) {
  int i, seed, size;
  coordT *coordp, *inputp;
  realT randr, randa, randb;

  if (!qh->input_points) { /* first call */
    qh->input_points= qh->first_point;
    qh->input_malloc= qh->POINTSmalloc;
    size= qh->num_points * qh->hull_dim * sizeof(coordT);
    if (!(qh->first_point=(coordT*)qh_malloc((size_t)size))) {
      qh_fprintf(qh, qh->ferr, 6009, "qhull error: insufficient memory to joggle %d points\n",
          qh->num_points);
      qh_errexit(qh, qh_ERRmem, NULL, NULL);
    }
    qh->POINTSmalloc= True;
    if (qh->JOGGLEmax == 0.0) {
      qh->JOGGLEmax= qh_detjoggle(qh, qh->input_points, qh->num_points, qh->hull_dim);
      qh_option(qh, "QJoggle", NULL, &qh->JOGGLEmax);
    }
  }else {                 /* repeated call */
    if (!qh->RERUN && qh->build_cnt > qh_JOGGLEretry) {
      if (((qh->build_cnt-qh_JOGGLEretry-1) % qh_JOGGLEagain) == 0) {
        realT maxjoggle= qh->MAXwidth * qh_JOGGLEmaxincrease;
        if (qh->JOGGLEmax < maxjoggle) {
          qh->JOGGLEmax *= qh_JOGGLEincrease;
          minimize_(qh->JOGGLEmax, maxjoggle);
        }
      }
    }
    qh_option(qh, "QJoggle", NULL, &qh->JOGGLEmax);
  }
  if (qh->build_cnt > 1 && qh->JOGGLEmax > fmax_(qh->MAXwidth/4, 0.1)) {
      qh_fprintf(qh, qh->ferr, 6010, "qhull error: the current joggle for 'QJn', %.2g, is too large for the width\nof the input.  If possible, recompile Qhull with higher-precision reals.\n",
                qh->JOGGLEmax);
      qh_errexit(qh, qh_ERRqhull, NULL, NULL);
  }
  /* for some reason, using qh->ROTATErandom and qh_RANDOMseed does not repeat the run. Use 'TRn' instead */
  seed= qh_RANDOMint;
  qh_option(qh, "_joggle-seed", &seed, NULL);
  trace0((qh, qh->ferr, 6, "qh_joggleinput: joggle input by %2.2g with seed %d\n",
    qh->JOGGLEmax, seed));
  inputp= qh->input_points;
  coordp= qh->first_point;
  randa= 2.0 * qh->JOGGLEmax/qh_RANDOMmax;
  randb= -qh->JOGGLEmax;
  size= qh->num_points * qh->hull_dim;
  for (i=size; i--; ) {
    randr= qh_RANDOMint;
    *(coordp++)= *(inputp++) + (randr * randa + randb);
  }
  if (qh->DELAUNAY) {
    qh->last_low= qh->last_high= qh->last_newhigh= REALmax;
    qh_setdelaunay(qh, qh->hull_dim, qh->num_points, qh->first_point);
  }
} /* joggleinput */

/*-<a                             href="qh-geom_r.htm#TOC"
  >-------------------------------</a><a name="maxabsval">-</a>

  qh_maxabsval( normal, dim )
    return pointer to maximum absolute value of a dim vector
    returns NULL if dim=0
*/
realT *qh_maxabsval(realT *normal, int dim) {
  realT maxval= -REALmax;
  realT *maxp= NULL, *colp, absval;
  int k;

  for (k=dim, colp= normal; k--; colp++) {
    absval= fabs_(*colp);
    if (absval > maxval) {
      maxval= absval;
      maxp= colp;
    }
  }
  return maxp;
} /* maxabsval */


/*-<a                             href="qh-geom_r.htm#TOC"
  >-------------------------------</a><a name="maxmin">-</a>

  qh_maxmin(qh, points, numpoints, dimension )
    return max/min points for each dimension
    determine max and min coordinates

  returns:
    returns a temporary set of max and min points
      may include duplicate points. Does not include qh.GOODpoint
    sets qh.NEARzero, qh.MAXabs_coord, qh.MAXsumcoord, qh.MAXwidth
         qh.MAXlastcoord, qh.MINlastcoord
    initializes qh.max_outside, qh.min_vertex, qh.WAScoplanar, qh.ZEROall_ok

  notes:
    loop duplicated in qh_detjoggle()

  design:
    initialize global precision variables
    checks definition of REAL...
    for each dimension
      for each point
        collect maximum and minimum point
      collect maximum of maximums and minimum of minimums
      determine qh.NEARzero for Gaussian Elimination
*/
setT *qh_maxmin(qhT *qh, pointT *points, int numpoints, int dimension) {
  int k;
  realT maxcoord, temp;
  pointT *minimum, *maximum, *point, *pointtemp;
  setT *set;

  qh->max_outside= 0.0;
  qh->MAXabs_coord= 0.0;
  qh->MAXwidth= -REALmax;
  qh->MAXsumcoord= 0.0;
  qh->min_vertex= 0.0;
  qh->WAScoplanar= False;
  if (qh->ZEROcentrum)
    qh->ZEROall_ok= True;
  if (REALmin < REALepsilon && REALmin < REALmax && REALmin > -REALmax
  && REALmax > 0.0 && -REALmax < 0.0)
    ; /* all ok */
  else {
    qh_fprintf(qh, qh->ferr, 6011, "qhull error: floating point constants in user.h are wrong\n\
REALepsilon %g REALmin %g REALmax %g -REALmax %g\n",
             REALepsilon, REALmin, REALmax, -REALmax);
    qh_errexit(qh, qh_ERRinput, NULL, NULL);
  }
  set= qh_settemp(qh, 2*dimension);
  for (k=0; k < dimension; k++) {
    if (points == qh->GOODpointp)
      minimum= maximum= points + dimension;
    else
      minimum= maximum= points;
    FORALLpoint_(qh, points, numpoints) {
      if (point == qh->GOODpointp)
        continue;
      if (maximum[k] < point[k])
        maximum= point;
      else if (minimum[k] > point[k])
        minimum= point;
    }
    if (k == dimension-1) {
      qh->MINlastcoord= minimum[k];
      qh->MAXlastcoord= maximum[k];
    }
    if (qh->SCALElast && k == dimension-1)
      maxcoord= qh->MAXwidth;
    else {
      maxcoord= fmax_(maximum[k], -minimum[k]);
      if (qh->GOODpointp) {
        temp= fmax_(qh->GOODpointp[k], -qh->GOODpointp[k]);
        maximize_(maxcoord, temp);
      }
      temp= maximum[k] - minimum[k];
      maximize_(qh->MAXwidth, temp);
    }
    maximize_(qh->MAXabs_coord, maxcoord);
    qh->MAXsumcoord += maxcoord;
    qh_setappend(qh, &set, maximum);
    qh_setappend(qh, &set, minimum);
    /* calculation of qh NEARzero is based on Golub & van Loan, 1983,
       Eq. 4.4-13 for "Gaussian elimination with complete pivoting".
       Golub & van Loan say that n^3 can be ignored and 10 be used in
       place of rho */
    qh->NEARzero[k]= 80 * qh->MAXsumcoord * REALepsilon;
  }
  if (qh->IStracing >=1)
    qh_printpoints(qh, qh->ferr, "qh_maxmin: found the max and min points(by dim):", set);
  return(set);
} /* maxmin */

/*-<a                             href="qh-geom_r.htm#TOC"
  >-------------------------------</a><a name="maxouter">-</a>

  qh_maxouter(qh)
    return maximum distance from facet to outer plane
    normally this is qh.max_outside+qh.DISTround
    does not include qh.JOGGLEmax

  see:
    qh_outerinner()

  notes:
    need to add another qh.DISTround if testing actual point with computation

  for joggle:
    qh_setfacetplane() updated qh.max_outer for Wnewvertexmax (max distance to vertex)
    need to use Wnewvertexmax since could have a coplanar point for a high
      facet that is replaced by a low facet
    need to add qh.JOGGLEmax if testing input points
*/
realT qh_maxouter(qhT *qh) {
  realT dist;

  dist= fmax_(qh->max_outside, qh->DISTround);
  dist += qh->DISTround;
  trace4((qh, qh->ferr, 4012, "qh_maxouter: max distance from facet to outer plane is %2.2g max_outside is %2.2g\n", dist, qh->max_outside));
  return dist;
} /* maxouter */

/*-<a                             href="qh-geom_r.htm#TOC"
  >-------------------------------</a><a name="maxsimplex">-</a>

  qh_maxsimplex(qh, dim, maxpoints, points, numpoints, simplex )
    determines maximum simplex for a set of points
    starts from points already in simplex
    skips qh.GOODpointp (assumes that it isn't in maxpoints)

  returns:
    simplex with dim+1 points

  notes:
    assumes at least pointsneeded points in points
    maximizes determinate for x,y,z,w, etc.
    uses maxpoints as long as determinate is clearly non-zero

  design:
    initialize simplex with at least two points
      (find points with max or min x coordinate)
    for each remaining dimension
      add point that maximizes the determinate
        (use points from maxpoints first)
*/
void qh_maxsimplex(qhT *qh, int dim, setT *maxpoints, pointT *points, int numpoints, setT **simplex) {
  pointT *point, **pointp, *pointtemp, *maxpoint, *minx=NULL, *maxx=NULL;
  boolT nearzero, maxnearzero= False;
  int k, sizinit;
  realT maxdet= -REALmax, det, mincoord= REALmax, maxcoord= -REALmax;

  sizinit= qh_setsize(qh, *simplex);
  if (sizinit < 2) {
    if (qh_setsize(qh, maxpoints) >= 2) {
      FOREACHpoint_(maxpoints) {
        if (maxcoord < point[0]) {
          maxcoord= point[0];
          maxx= point;
        }
        if (mincoord > point[0]) {
          mincoord= point[0];
          minx= point;
        }
      }
    }else {
      FORALLpoint_(qh, points, numpoints) {
        if (point == qh->GOODpointp)
          continue;
        if (maxcoord < point[0]) {
          maxcoord= point[0];
          maxx= point;
        }
        if (mincoord > point[0]) {
          mincoord= point[0];
          minx= point;
        }
      }
    }
    qh_setunique(qh, simplex, minx);
    if (qh_setsize(qh, *simplex) < 2)
      qh_setunique(qh, simplex, maxx);
    sizinit= qh_setsize(qh, *simplex);
    if (sizinit < 2) {
      qh_precision(qh, "input has same x coordinate");
      if (zzval_(Zsetplane) > qh->hull_dim+1) {
        qh_fprintf(qh, qh->ferr, 6012, "qhull precision error (qh_maxsimplex for voronoi_center):\n%d points with the same x coordinate.\n",
                 qh_setsize(qh, maxpoints)+numpoints);
        qh_errexit(qh, qh_ERRprec, NULL, NULL);
      }else {
        qh_fprintf(qh, qh->ferr, 6013, "qhull input error: input is less than %d-dimensional since it has the same x coordinate\n", qh->hull_dim);
        qh_errexit(qh, qh_ERRinput, NULL, NULL);
      }
    }
  }
  for (k=sizinit; k < dim+1; k++) {
    maxpoint= NULL;
    maxdet= -REALmax;
    FOREACHpoint_(maxpoints) {
      if (!qh_setin(*simplex, point)) {
        det= qh_detsimplex(qh, point, *simplex, k, &nearzero);
        if ((det= fabs_(det)) > maxdet) {
          maxdet= det;
          maxpoint= point;
          maxnearzero= nearzero;
        }
      }
    }
    if (!maxpoint || maxnearzero) {
      zinc_(Zsearchpoints);
      if (!maxpoint) {
        trace0((qh, qh->ferr, 7, "qh_maxsimplex: searching all points for %d-th initial vertex.\n", k+1));
      }else {
        trace0((qh, qh->ferr, 8, "qh_maxsimplex: searching all points for %d-th initial vertex, better than p%d det %2.2g\n",
                k+1, qh_pointid(qh, maxpoint), maxdet));
      }
      FORALLpoint_(qh, points, numpoints) {
        if (point == qh->GOODpointp)
          continue;
        if (!qh_setin(*simplex, point)) {
          det= qh_detsimplex(qh, point, *simplex, k, &nearzero);
          if ((det= fabs_(det)) > maxdet) {
            maxdet= det;
            maxpoint= point;
            maxnearzero= nearzero;
          }
        }
      }
    } /* !maxpoint */
    if (!maxpoint) {
      qh_fprintf(qh, qh->ferr, 6014, "qhull internal error (qh_maxsimplex): not enough points available\n");
      qh_errexit(qh, qh_ERRqhull, NULL, NULL);
    }
    qh_setappend(qh, simplex, maxpoint);
    trace1((qh, qh->ferr, 1002, "qh_maxsimplex: selected point p%d for %d`th initial vertex, det=%2.2g\n",
            qh_pointid(qh, maxpoint), k+1, maxdet));
  } /* k */
} /* maxsimplex */

/*-<a                             href="qh-geom_r.htm#TOC"
  >-------------------------------</a><a name="minabsval">-</a>

  qh_minabsval( normal, dim )
    return minimum absolute value of a dim vector
*/
realT qh_minabsval(realT *normal, int dim) {
  realT minval= 0;
  realT maxval= 0;
  realT *colp;
  int k;

  for (k=dim, colp=normal; k--; colp++) {
    maximize_(maxval, *colp);
    minimize_(minval, *colp);
  }
  return fmax_(maxval, -minval);
} /* minabsval */


/*-<a                             href="qh-geom_r.htm#TOC"
  >-------------------------------</a><a name="mindiff">-</a>

  qh_mindif(qh, vecA, vecB, dim )
    return index of min abs. difference of two vectors
*/
int qh_mindiff(realT *vecA, realT *vecB, int dim) {
  realT mindiff= REALmax, diff;
  realT *vecAp= vecA, *vecBp= vecB;
  int k, mink= 0;

  for (k=0; k < dim; k++) {
    diff= *vecAp++ - *vecBp++;
    diff= fabs_(diff);
    if (diff < mindiff) {
      mindiff= diff;
      mink= k;
    }
  }
  return mink;
} /* mindiff */


/*-<a                             href="qh-geom_r.htm#TOC"
  >-------------------------------</a><a name="orientoutside">-</a>

  qh_orientoutside(qh, facet  )
    make facet outside oriented via qh.interior_point

  returns:
    True if facet reversed orientation.
*/
boolT qh_orientoutside(qhT *qh, facetT *facet) {
  int k;
  realT dist;

  qh_distplane(qh, qh->interior_point, facet, &dist);
  if (dist > 0) {
    for (k=qh->hull_dim; k--; )
      facet->normal[k]= -facet->normal[k];
    facet->offset= -facet->offset;
    return True;
  }
  return False;
} /* orientoutside */

/*-<a                             href="qh-geom_r.htm#TOC"
  >-------------------------------</a><a name="outerinner">-</a>

  qh_outerinner(qh, facet, outerplane, innerplane  )
    if facet and qh.maxoutdone (i.e., qh_check_maxout)
      returns outer and inner plane for facet
    else
      returns maximum outer and inner plane
    accounts for qh.JOGGLEmax

  see:
    qh_maxouter(qh), qh_check_bestdist(), qh_check_points()

  notes:
    outerplaner or innerplane may be NULL
    facet is const
    Does not error (QhullFacet)

    includes qh.DISTround for actual points
    adds another qh.DISTround if testing with floating point arithmetic
*/
void qh_outerinner(qhT *qh, facetT *facet, realT *outerplane, realT *innerplane) {
  realT dist, mindist;
  vertexT *vertex, **vertexp;

  if (outerplane) {
    if (!qh_MAXoutside || !facet || !qh->maxoutdone) {
      *outerplane= qh_maxouter(qh);       /* includes qh.DISTround */
    }else { /* qh_MAXoutside ... */
#if qh_MAXoutside
      *outerplane= facet->maxoutside + qh->DISTround;
#endif

    }
    if (qh->JOGGLEmax < REALmax/2)
      *outerplane += qh->JOGGLEmax * sqrt((realT)qh->hull_dim);
  }
  if (innerplane) {
    if (facet) {
      mindist= REALmax;
      FOREACHvertex_(facet->vertices) {
        zinc_(Zdistio);
        qh_distplane(qh, vertex->point, facet, &dist);
        minimize_(mindist, dist);
      }
      *innerplane= mindist - qh->DISTround;
    }else
      *innerplane= qh->min_vertex - qh->DISTround;
    if (qh->JOGGLEmax < REALmax/2)
      *innerplane -= qh->JOGGLEmax * sqrt((realT)qh->hull_dim);
  }
} /* outerinner */

/*-<a                             href="qh-geom_r.htm#TOC"
  >-------------------------------</a><a name="pointdist">-</a>

  qh_pointdist( point1, point2, dim )
    return distance between two points

  notes:
    returns distance squared if 'dim' is negative
*/
coordT qh_pointdist(pointT *point1, pointT *point2, int dim) {
  coordT dist, diff;
  int k;

  dist= 0.0;
  for (k= (dim > 0 ? dim : -dim); k--; ) {
    diff= *point1++ - *point2++;
    dist += diff * diff;
  }
  if (dim > 0)
    return(sqrt(dist));
  return dist;
} /* pointdist */


/*-<a                             href="qh-geom_r.htm#TOC"
  >-------------------------------</a><a name="printmatrix">-</a>

  qh_printmatrix(qh, fp, string, rows, numrow, numcol )
    print matrix to fp given by row vectors
    print string as header
    qh may be NULL if fp is defined

  notes:
    print a vector by qh_printmatrix(qh, fp, "", &vect, 1, len)
*/
void qh_printmatrix(qhT *qh, FILE *fp, const char *string, realT **rows, int numrow, int numcol) {
  realT *rowp;
  realT r; /*bug fix*/
  int i,k;

  qh_fprintf(qh, fp, 9001, "%s\n", string);
  for (i=0; i < numrow; i++) {
    rowp= rows[i];
    for (k=0; k < numcol; k++) {
      r= *rowp++;
      qh_fprintf(qh, fp, 9002, "%6.3g ", r);
    }
    qh_fprintf(qh, fp, 9003, "\n");
  }
} /* printmatrix */


/*-<a                             href="qh-geom_r.htm#TOC"
  >-------------------------------</a><a name="printpoints">-</a>

  qh_printpoints(qh, fp, string, points )
    print pointids to fp for a set of points
    if string, prints string and 'p' point ids
*/
void qh_printpoints(qhT *qh, FILE *fp, const char *string, setT *points) {
  pointT *point, **pointp;

  if (string) {
    qh_fprintf(qh, fp, 9004, "%s", string);
    FOREACHpoint_(points)
      qh_fprintf(qh, fp, 9005, " p%d", qh_pointid(qh, point));
    qh_fprintf(qh, fp, 9006, "\n");
  }else {
    FOREACHpoint_(points)
      qh_fprintf(qh, fp, 9007, " %d", qh_pointid(qh, point));
    qh_fprintf(qh, fp, 9008, "\n");
  }
} /* printpoints */


/*-<a                             href="qh-geom_r.htm#TOC"
  >-------------------------------</a><a name="projectinput">-</a>

  qh_projectinput(qh)
    project input points using qh.lower_bound/upper_bound and qh->DELAUNAY
    if qh.lower_bound[k]=qh.upper_bound[k]= 0,
      removes dimension k
    if halfspace intersection
      removes dimension k from qh.feasible_point
    input points in qh->first_point, num_points, input_dim

  returns:
    new point array in qh->first_point of qh->hull_dim coordinates
    sets qh->POINTSmalloc
    if qh->DELAUNAY
      projects points to paraboloid
      lowbound/highbound is also projected
    if qh->ATinfinity
      adds point "at-infinity"
    if qh->POINTSmalloc
      frees old point array

  notes:
    checks that qh.hull_dim agrees with qh.input_dim, PROJECTinput, and DELAUNAY


  design:
    sets project[k] to -1 (delete), 0 (keep), 1 (add for Delaunay)
    determines newdim and newnum for qh->hull_dim and qh->num_points
    projects points to newpoints
    projects qh.lower_bound to itself
    projects qh.upper_bound to itself
    if qh->DELAUNAY
      if qh->ATINFINITY
        projects points to paraboloid
        computes "infinity" point as vertex average and 10% above all points
      else
        uses qh_setdelaunay to project points to paraboloid
*/
void qh_projectinput(qhT *qh) {
  int k,i;
  int newdim= qh->input_dim, newnum= qh->num_points;
  signed char *project;
  int projectsize= (qh->input_dim+1)*sizeof(*project);
  pointT *newpoints, *coord, *infinity;
  realT paraboloid, maxboloid= 0;

  project= (signed char*)qh_memalloc(qh, projectsize);
  memset((char*)project, 0, (size_t)projectsize);
  for (k=0; k < qh->input_dim; k++) {   /* skip Delaunay bound */
    if (qh->lower_bound[k] == 0 && qh->upper_bound[k] == 0) {
      project[k]= -1;
      newdim--;
    }
  }
  if (qh->DELAUNAY) {
    project[k]= 1;
    newdim++;
    if (qh->ATinfinity)
      newnum++;
  }
  if (newdim != qh->hull_dim) {
    qh_memfree(qh, project, projectsize);
    qh_fprintf(qh, qh->ferr, 6015, "qhull internal error (qh_projectinput): dimension after projection %d != hull_dim %d\n", newdim, qh->hull_dim);
    qh_errexit(qh, qh_ERRqhull, NULL, NULL);
  }
  if (!(newpoints= qh->temp_malloc= (coordT*)qh_malloc(newnum*newdim*sizeof(coordT)))){
    qh_memfree(qh, project, projectsize);
    qh_fprintf(qh, qh->ferr, 6016, "qhull error: insufficient memory to project %d points\n",
           qh->num_points);
    qh_errexit(qh, qh_ERRmem, NULL, NULL);
  }
  /* qh_projectpoints throws error if mismatched dimensions */
  qh_projectpoints(qh, project, qh->input_dim+1, qh->first_point,
                    qh->num_points, qh->input_dim, newpoints, newdim);
  trace1((qh, qh->ferr, 1003, "qh_projectinput: updating lower and upper_bound\n"));
  qh_projectpoints(qh, project, qh->input_dim+1, qh->lower_bound,
                    1, qh->input_dim+1, qh->lower_bound, newdim+1);
  qh_projectpoints(qh, project, qh->input_dim+1, qh->upper_bound,
                    1, qh->input_dim+1, qh->upper_bound, newdim+1);
  if (qh->HALFspace) {
    if (!qh->feasible_point) {
      qh_memfree(qh, project, projectsize);
      qh_fprintf(qh, qh->ferr, 6017, "qhull internal error (qh_projectinput): HALFspace defined without qh.feasible_point\n");
      qh_errexit(qh, qh_ERRqhull, NULL, NULL);
    }
    qh_projectpoints(qh, project, qh->input_dim, qh->feasible_point,
                      1, qh->input_dim, qh->feasible_point, newdim);
  }
  qh_memfree(qh, project, projectsize);
  if (qh->POINTSmalloc)
    qh_free(qh->first_point);
  qh->first_point= newpoints;
  qh->POINTSmalloc= True;
  qh->temp_malloc= NULL;
  if (qh->DELAUNAY && qh->ATinfinity) {
    coord= qh->first_point;
    infinity= qh->first_point + qh->hull_dim * qh->num_points;
    for (k=qh->hull_dim-1; k--; )
      infinity[k]= 0.0;
    for (i=qh->num_points; i--; ) {
      paraboloid= 0.0;
      for (k=0; k < qh->hull_dim-1; k++) {
        paraboloid += *coord * *coord;
        infinity[k] += *coord;
        coord++;
      }
      *(coord++)= paraboloid;
      maximize_(maxboloid, paraboloid);
    }
    /* coord == infinity */
    for (k=qh->hull_dim-1; k--; )
      *(coord++) /= qh->num_points;
    *(coord++)= maxboloid * 1.1;
    qh->num_points++;
    trace0((qh, qh->ferr, 9, "qh_projectinput: projected points to paraboloid for Delaunay\n"));
  }else if (qh->DELAUNAY)  /* !qh->ATinfinity */
    qh_setdelaunay(qh, qh->hull_dim, qh->num_points, qh->first_point);
} /* projectinput */


/*-<a                             href="qh-geom_r.htm#TOC"
  >-------------------------------</a><a name="projectpoints">-</a>

  qh_projectpoints(qh, project, n, points, numpoints, dim, newpoints, newdim )
    project points/numpoints/dim to newpoints/newdim
    if project[k] == -1
      delete dimension k
    if project[k] == 1
      add dimension k by duplicating previous column
    n is size of project

  notes:
    newpoints may be points if only adding dimension at end

  design:
    check that 'project' and 'newdim' agree
    for each dimension
      if project == -1
        skip dimension
      else
        determine start of column in newpoints
        determine start of column in points
          if project == +1, duplicate previous column
        copy dimension (column) from points to newpoints
*/
void qh_projectpoints(qhT *qh, signed char *project, int n, realT *points,
        int numpoints, int dim, realT *newpoints, int newdim) {
  int testdim= dim, oldk=0, newk=0, i,j=0,k;
  realT *newp, *oldp;

  for (k=0; k < n; k++)
    testdim += project[k];
  if (testdim != newdim) {
    qh_fprintf(qh, qh->ferr, 6018, "qhull internal error (qh_projectpoints): newdim %d should be %d after projection\n",
      newdim, testdim);
    qh_errexit(qh, qh_ERRqhull, NULL, NULL);
  }
  for (j=0; j<n; j++) {
    if (project[j] == -1)
      oldk++;
    else {
      newp= newpoints+newk++;
      if (project[j] == +1) {
        if (oldk >= dim)
          continue;
        oldp= points+oldk;
      }else
        oldp= points+oldk++;
      for (i=numpoints; i--; ) {
        *newp= *oldp;
        newp += newdim;
        oldp += dim;
      }
    }
    if (oldk >= dim)
      break;
  }
  trace1((qh, qh->ferr, 1004, "qh_projectpoints: projected %d points from dim %d to dim %d\n",
    numpoints, dim, newdim));
} /* projectpoints */


/*-<a                             href="qh-geom_r.htm#TOC"
  >-------------------------------</a><a name="rotateinput">-</a>

  qh_rotateinput(qh, rows )
    rotate input using row matrix
    input points given by qh->first_point, num_points, hull_dim
    assumes rows[dim] is a scratch buffer
    if qh->POINTSmalloc, overwrites input points, else mallocs a new array

  returns:
    rotated input
    sets qh->POINTSmalloc

  design:
    see qh_rotatepoints
*/
void qh_rotateinput(qhT *qh, realT **rows) {

  if (!qh->POINTSmalloc) {
    qh->first_point= qh_copypoints(qh, qh->first_point, qh->num_points, qh->hull_dim);
    qh->POINTSmalloc= True;
  }
  qh_rotatepoints(qh, qh->first_point, qh->num_points, qh->hull_dim, rows);
}  /* rotateinput */

/*-<a                             href="qh-geom_r.htm#TOC"
  >-------------------------------</a><a name="rotatepoints">-</a>

  qh_rotatepoints(qh, points, numpoints, dim, row )
    rotate numpoints points by a d-dim row matrix
    assumes rows[dim] is a scratch buffer

  returns:
    rotated points in place

  design:
    for each point
      for each coordinate
        use row[dim] to compute partial inner product
      for each coordinate
        rotate by partial inner product
*/
void qh_rotatepoints(qhT *qh, realT *points, int numpoints, int dim, realT **row) {
  realT *point, *rowi, *coord= NULL, sum, *newval;
  int i,j,k;

  if (qh->IStracing >= 1)
    qh_printmatrix(qh, qh->ferr, "qh_rotatepoints: rotate points by", row, dim, dim);
  for (point= points, j= numpoints; j--; point += dim) {
    newval= row[dim];
    for (i=0; i < dim; i++) {
      rowi= row[i];
      coord= point;
      for (sum= 0.0, k= dim; k--; )
        sum += *rowi++ * *coord++;
      *(newval++)= sum;
    }
    for (k=dim; k--; )
      *(--coord)= *(--newval);
  }
} /* rotatepoints */


/*-<a                             href="qh-geom_r.htm#TOC"
  >-------------------------------</a><a name="scaleinput">-</a>

  qh_scaleinput(qh)
    scale input points using qh->low_bound/high_bound
    input points given by qh->first_point, num_points, hull_dim
    if qh->POINTSmalloc, overwrites input points, else mallocs a new array

  returns:
    scales coordinates of points to low_bound[k], high_bound[k]
    sets qh->POINTSmalloc

  design:
    see qh_scalepoints
*/
void qh_scaleinput(qhT *qh) {

  if (!qh->POINTSmalloc) {
    qh->first_point= qh_copypoints(qh, qh->first_point, qh->num_points, qh->hull_dim);
    qh->POINTSmalloc= True;
  }
  qh_scalepoints(qh, qh->first_point, qh->num_points, qh->hull_dim,
       qh->lower_bound, qh->upper_bound);
}  /* scaleinput */

/*-<a                             href="qh-geom_r.htm#TOC"
  >-------------------------------</a><a name="scalelast">-</a>

  qh_scalelast(qh, points, numpoints, dim, low, high, newhigh )
    scale last coordinate to [0,m] for Delaunay triangulations
    input points given by points, numpoints, dim

  returns:
    changes scale of last coordinate from [low, high] to [0, newhigh]
    overwrites last coordinate of each point
    saves low/high/newhigh in qh.last_low, etc. for qh_setdelaunay()

  notes:
    when called by qh_setdelaunay, low/high may not match actual data

  design:
    compute scale and shift factors
    apply to last coordinate of each point
*/
void qh_scalelast(qhT *qh, coordT *points, int numpoints, int dim, coordT low,
                   coordT high, coordT newhigh) {
  realT scale, shift;
  coordT *coord;
  int i;
  boolT nearzero= False;

  trace4((qh, qh->ferr, 4013, "qh_scalelast: scale last coordinate from [%2.2g, %2.2g] to [0,%2.2g]\n",
    low, high, newhigh));
  qh->last_low= low;
  qh->last_high= high;
  qh->last_newhigh= newhigh;
  scale= qh_divzero(newhigh, high - low,
                  qh->MINdenom_1, &nearzero);
  if (nearzero) {
    if (qh->DELAUNAY)
      qh_fprintf(qh, qh->ferr, 6019, "qhull input error: can not scale last coordinate.  Input is cocircular\n   or cospherical.   Use option 'Qz' to add a point at infinity.\n");
    else
      qh_fprintf(qh, qh->ferr, 6020, "qhull input error: can not scale last coordinate.  New bounds [0, %2.2g] are too wide for\nexisting bounds [%2.2g, %2.2g] (width %2.2g)\n",
                newhigh, low, high, high-low);
    qh_errexit(qh, qh_ERRinput, NULL, NULL);
  }
  shift= - low * newhigh / (high-low);
  coord= points + dim - 1;
  for (i=numpoints; i--; coord += dim)
    *coord= *coord * scale + shift;
} /* scalelast */

/*-<a                             href="qh-geom_r.htm#TOC"
  >-------------------------------</a><a name="scalepoints">-</a>

  qh_scalepoints(qh, points, numpoints, dim, newlows, newhighs )
    scale points to new lowbound and highbound
    retains old bound when newlow= -REALmax or newhigh= +REALmax

  returns:
    scaled points
    overwrites old points

  design:
    for each coordinate
      compute current low and high bound
      compute scale and shift factors
      scale all points
      enforce new low and high bound for all points
*/
void qh_scalepoints(qhT *qh, pointT *points, int numpoints, int dim,
        realT *newlows, realT *newhighs) {
  int i,k;
  realT shift, scale, *coord, low, high, newlow, newhigh, mincoord, maxcoord;
  boolT nearzero= False;

  for (k=0; k < dim; k++) {
    newhigh= newhighs[k];
    newlow= newlows[k];
    if (newhigh > REALmax/2 && newlow < -REALmax/2)
      continue;
    low= REALmax;
    high= -REALmax;
    for (i=numpoints, coord=points+k; i--; coord += dim) {
      minimize_(low, *coord);
      maximize_(high, *coord);
    }
    if (newhigh > REALmax/2)
      newhigh= high;
    if (newlow < -REALmax/2)
      newlow= low;
    if (qh->DELAUNAY && k == dim-1 && newhigh < newlow) {
      qh_fprintf(qh, qh->ferr, 6021, "qhull input error: 'Qb%d' or 'QB%d' inverts paraboloid since high bound %.2g < low bound %.2g\n",
               k, k, newhigh, newlow);
      qh_errexit(qh, qh_ERRinput, NULL, NULL);
    }
    scale= qh_divzero(newhigh - newlow, high - low,
                  qh->MINdenom_1, &nearzero);
    if (nearzero) {
      qh_fprintf(qh, qh->ferr, 6022, "qhull input error: %d'th dimension's new bounds [%2.2g, %2.2g] too wide for\nexisting bounds [%2.2g, %2.2g]\n",
              k, newlow, newhigh, low, high);
      qh_errexit(qh, qh_ERRinput, NULL, NULL);
    }
    shift= (newlow * high - low * newhigh)/(high-low);
    coord= points+k;
    for (i=numpoints; i--; coord += dim)
      *coord= *coord * scale + shift;
    coord= points+k;
    if (newlow < newhigh) {
      mincoord= newlow;
      maxcoord= newhigh;
    }else {
      mincoord= newhigh;
      maxcoord= newlow;
    }
    for (i=numpoints; i--; coord += dim) {
      minimize_(*coord, maxcoord);  /* because of roundoff error */
      maximize_(*coord, mincoord);
    }
    trace0((qh, qh->ferr, 10, "qh_scalepoints: scaled %d'th coordinate [%2.2g, %2.2g] to [%.2g, %.2g] for %d points by %2.2g and shifted %2.2g\n",
      k, low, high, newlow, newhigh, numpoints, scale, shift));
  }
} /* scalepoints */


/*-<a                             href="qh-geom_r.htm#TOC"
  >-------------------------------</a><a name="setdelaunay">-</a>

  qh_setdelaunay(qh, dim, count, points )
    project count points to dim-d paraboloid for Delaunay triangulation

    dim is one more than the dimension of the input set
    assumes dim is at least 3 (i.e., at least a 2-d Delaunay triangulation)

    points is a dim*count realT array.  The first dim-1 coordinates
    are the coordinates of the first input point.  array[dim] is
    the first coordinate of the second input point.  array[2*dim] is
    the first coordinate of the third input point.

    if qh.last_low defined (i.e., 'Qbb' called qh_scalelast)
      calls qh_scalelast to scale the last coordinate the same as the other points

  returns:
    for each point
      sets point[dim-1] to sum of squares of coordinates
    scale points to 'Qbb' if needed

  notes:
    to project one point, use
      qh_setdelaunay(qh, qh->hull_dim, 1, point)

    Do not use options 'Qbk', 'QBk', or 'QbB' since they scale
    the coordinates after the original projection.

*/
void qh_setdelaunay(qhT *qh, int dim, int count, pointT *points) {
  int i, k;
  coordT *coordp, coord;
  realT paraboloid;

  trace0((qh, qh->ferr, 11, "qh_setdelaunay: project %d points to paraboloid for Delaunay triangulation\n", count));
  coordp= points;
  for (i=0; i < count; i++) {
    coord= *coordp++;
    paraboloid= coord*coord;
    for (k=dim-2; k--; ) {
      coord= *coordp++;
      paraboloid += coord*coord;
    }
    *coordp++ = paraboloid;
  }
  if (qh->last_low < REALmax/2)
    qh_scalelast(qh, points, count, dim, qh->last_low, qh->last_high, qh->last_newhigh);
} /* setdelaunay */


/*-<a                             href="qh-geom_r.htm#TOC"
  >-------------------------------</a><a name="sethalfspace">-</a>

  qh_sethalfspace(qh, dim, coords, nextp, normal, offset, feasible )
    set point to dual of halfspace relative to feasible point
    halfspace is normal coefficients and offset.

  returns:
    false and prints error if feasible point is outside of hull
    overwrites coordinates for point at dim coords
    nextp= next point (coords)
    does not call qh_errexit

  design:
    compute distance from feasible point to halfspace
    divide each normal coefficient by -dist
*/
boolT qh_sethalfspace(qhT *qh, int dim, coordT *coords, coordT **nextp,
         coordT *normal, coordT *offset, coordT *feasible) {
  coordT *normp= normal, *feasiblep= feasible, *coordp= coords;
  realT dist;
  realT r; /*bug fix*/
  int k;
  boolT zerodiv;

  dist= *offset;
  for (k=dim; k--; )
    dist += *(normp++) * *(feasiblep++);
  if (dist > 0)
    goto LABELerroroutside;
  normp= normal;
  if (dist < -qh->MINdenom) {
    for (k=dim; k--; )
      *(coordp++)= *(normp++) / -dist;
  }else {
    for (k=dim; k--; ) {
      *(coordp++)= qh_divzero(*(normp++), -dist, qh->MINdenom_1, &zerodiv);
      if (zerodiv)
        goto LABELerroroutside;
    }
  }
  *nextp= coordp;
  if (qh->IStracing >= 4) {
    qh_fprintf(qh, qh->ferr, 8021, "qh_sethalfspace: halfspace at offset %6.2g to point: ", *offset);
    for (k=dim, coordp=coords; k--; ) {
      r= *coordp++;
      qh_fprintf(qh, qh->ferr, 8022, " %6.2g", r);
    }
    qh_fprintf(qh, qh->ferr, 8023, "\n");
  }
  return True;
LABELerroroutside:
  feasiblep= feasible;
  normp= normal;
  qh_fprintf(qh, qh->ferr, 6023, "qhull input error: feasible point is not clearly inside halfspace\nfeasible point: ");
  for (k=dim; k--; )
    qh_fprintf(qh, qh->ferr, 8024, qh_REAL_1, r=*(feasiblep++));
  qh_fprintf(qh, qh->ferr, 8025, "\n     halfspace: ");
  for (k=dim; k--; )
    qh_fprintf(qh, qh->ferr, 8026, qh_REAL_1, r=*(normp++));
  qh_fprintf(qh, qh->ferr, 8027, "\n     at offset: ");
  qh_fprintf(qh, qh->ferr, 8028, qh_REAL_1, *offset);
  qh_fprintf(qh, qh->ferr, 8029, " and distance: ");
  qh_fprintf(qh, qh->ferr, 8030, qh_REAL_1, dist);
  qh_fprintf(qh, qh->ferr, 8031, "\n");
  return False;
} /* sethalfspace */

/*-<a                             href="qh-geom_r.htm#TOC"
  >-------------------------------</a><a name="sethalfspace_all">-</a>

  qh_sethalfspace_all(qh, dim, count, halfspaces, feasible )
    generate dual for halfspace intersection with feasible point
    array of count halfspaces
      each halfspace is normal coefficients followed by offset
      the origin is inside the halfspace if the offset is negative
    feasible is a point inside all halfspaces (http://www.qhull.org/html/qhalf.htm#notes)

  returns:
    malloc'd array of count X dim-1 points

  notes:
    call before qh_init_B or qh_initqhull_globals
    free memory when done
    unused/untested code: please email bradb@shore.net if this works ok for you
    if using option 'Fp', qh->feasible_point must be set (e.g., to 'feasible')
    qh->feasible_point is a malloc'd array that is freed by qh_freebuffers.

  design:
    see qh_sethalfspace
*/
coordT *qh_sethalfspace_all(qhT *qh, int dim, int count, coordT *halfspaces, pointT *feasible) {
  int i, newdim;
  pointT *newpoints;
  coordT *coordp, *normalp, *offsetp;

  trace0((qh, qh->ferr, 12, "qh_sethalfspace_all: compute dual for halfspace intersection\n"));
  newdim= dim - 1;
  if (!(newpoints=(coordT*)qh_malloc(count*newdim*sizeof(coordT)))){
    qh_fprintf(qh, qh->ferr, 6024, "qhull error: insufficient memory to compute dual of %d halfspaces\n",
          count);
    qh_errexit(qh, qh_ERRmem, NULL, NULL);
  }
  coordp= newpoints;
  normalp= halfspaces;
  for (i=0; i < count; i++) {
    offsetp= normalp + newdim;
    if (!qh_sethalfspace(qh, newdim, coordp, &coordp, normalp, offsetp, feasible)) {
      qh_free(newpoints);  /* feasible is not inside halfspace as reported by qh_sethalfspace */
      qh_fprintf(qh, qh->ferr, 8032, "The halfspace was at index %d\n", i);
      qh_errexit(qh, qh_ERRinput, NULL, NULL);
    }
    normalp= offsetp + 1;
  }
  return newpoints;
} /* sethalfspace_all */


/*-<a                             href="qh-geom_r.htm#TOC"
  >-------------------------------</a><a name="sharpnewfacets">-</a>

  qh_sharpnewfacets(qh)

  returns:
    true if could be an acute angle (facets in different quadrants)

  notes:
    for qh_findbest

  design:
    for all facets on qh.newfacet_list
      if two facets are in different quadrants
        set issharp
*/
boolT qh_sharpnewfacets(qhT *qh) {
  facetT *facet;
  boolT issharp = False;
  int *quadrant, k;

  quadrant= (int*)qh_memalloc(qh, qh->hull_dim * sizeof(int));
  FORALLfacet_(qh->newfacet_list) {
    if (facet == qh->newfacet_list) {
      for (k=qh->hull_dim; k--; )
        quadrant[ k]= (facet->normal[ k] > 0);
    }else {
      for (k=qh->hull_dim; k--; ) {
        if (quadrant[ k] != (facet->normal[ k] > 0)) {
          issharp= True;
          break;
        }
      }
    }
    if (issharp)
      break;
  }
  qh_memfree(qh, quadrant, qh->hull_dim * sizeof(int));
  trace3((qh, qh->ferr, 3001, "qh_sharpnewfacets: %d\n", issharp));
  return issharp;
} /* sharpnewfacets */

/*-<a                             href="qh-geom_r.htm#TOC"
  >-------------------------------</a><a name="voronoi_center">-</a>

  qh_voronoi_center(qh, dim, points )
    return Voronoi center for a set of points
    dim is the orginal dimension of the points
    gh.gm_matrix/qh.gm_row are scratch buffers

  returns:
    center as a temporary point (qh_memalloc)
    if non-simplicial,
      returns center for max simplex of points

  notes:
    only called by qh_facetcenter
    from Bowyer & Woodwark, A Programmer's Geometry, 1983, p. 65

  design:
    if non-simplicial
      determine max simplex for points
    translate point0 of simplex to origin
    compute sum of squares of diagonal
    compute determinate
    compute Voronoi center (see Bowyer & Woodwark)
*/
pointT *qh_voronoi_center(qhT *qh, int dim, setT *points) {
  pointT *point, **pointp, *point0;
  pointT *center= (pointT*)qh_memalloc(qh, qh->center_size);
  setT *simplex;
  int i, j, k, size= qh_setsize(qh, points);
  coordT *gmcoord;
  realT *diffp, sum2, *sum2row, *sum2p, det, factor;
  boolT nearzero, infinite;

  if (size == dim+1)
    simplex= points;
  else if (size < dim+1) {
    qh_memfree(qh, center, qh->center_size);
    qh_fprintf(qh, qh->ferr, 6025, "qhull internal error (qh_voronoi_center):\n  need at least %d points to construct a Voronoi center\n",
             dim+1);
    qh_errexit(qh, qh_ERRqhull, NULL, NULL);
    simplex= points;  /* never executed -- avoids warning */
  }else {
    simplex= qh_settemp(qh, dim+1);
    qh_maxsimplex(qh, dim, points, NULL, 0, &simplex);
  }
  point0= SETfirstt_(simplex, pointT);
  gmcoord= qh->gm_matrix;
  for (k=0; k < dim; k++) {
    qh->gm_row[k]= gmcoord;
    FOREACHpoint_(simplex) {
      if (point != point0)
        *(gmcoord++)= point[k] - point0[k];
    }
  }
  sum2row= gmcoord;
  for (i=0; i < dim; i++) {
    sum2= 0.0;
    for (k=0; k < dim; k++) {
      diffp= qh->gm_row[k] + i;
      sum2 += *diffp * *diffp;
    }
    *(gmcoord++)= sum2;
  }
  det= qh_determinant(qh, qh->gm_row, dim, &nearzero);
  factor= qh_divzero(0.5, det, qh->MINdenom, &infinite);
  if (infinite) {
    for (k=dim; k--; )
      center[k]= qh_INFINITE;
    if (qh->IStracing)
      qh_printpoints(qh, qh->ferr, "qh_voronoi_center: at infinity for ", simplex);
  }else {
    for (i=0; i < dim; i++) {
      gmcoord= qh->gm_matrix;
      sum2p= sum2row;
      for (k=0; k < dim; k++) {
        qh->gm_row[k]= gmcoord;
        if (k == i) {
          for (j=dim; j--; )
            *(gmcoord++)= *sum2p++;
        }else {
          FOREACHpoint_(simplex) {
            if (point != point0)
              *(gmcoord++)= point[k] - point0[k];
          }
        }
      }
      center[i]= qh_determinant(qh, qh->gm_row, dim, &nearzero)*factor + point0[i];
    }
#ifndef qh_NOtrace
    if (qh->IStracing >= 3) {
      qh_fprintf(qh, qh->ferr, 8033, "qh_voronoi_center: det %2.2g factor %2.2g ", det, factor);
      qh_printmatrix(qh, qh->ferr, "center:", &center, 1, dim);
      if (qh->IStracing >= 5) {
        qh_printpoints(qh, qh->ferr, "points", simplex);
        FOREACHpoint_(simplex)
          qh_fprintf(qh, qh->ferr, 8034, "p%d dist %.2g, ", qh_pointid(qh, point),
                   qh_pointdist(point, center, dim));
        qh_fprintf(qh, qh->ferr, 8035, "\n");
      }
    }
#endif
  }
  if (simplex != points)
    qh_settempfree(qh, &simplex);
  return center;
} /* voronoi_center */