qhull_geom2.cxx

hl2 source code. Do not use it illegal.

    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-c.htm#geom"
  >-------------------------------</a><a name="facetarea">-</a>

  qh_facetarea( 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
  
  design:
    if simplicial
      compute area
    else
      for each ridge
        compute area from centrum to ridge
    negate area if upper Delaunay facet
*/
realT qh_facetarea (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 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 (facet);
    FOREACHridge_(facet->ridges) 
      area += qh_facetarea_simplex (qh hull_dim, centrum, ridge->vertices, 
                 NULL, (ridge->top == facet),  facet->normal, &facet->offset);
    if (qh CENTERtype != qh_AScentrum)
      qh_memfree (centrum, qh normal_size);
  }
  if (facet->upperdelaunay && qh DELAUNAY)
    area= -area;  /* the normal should be [0,...,1] */
  trace4((qh ferr, "qh_facetarea: f%d area %2.2g\n", facet->id, area)); 
  return area;
} /* facetarea */

/*-<a                             href="qh-c.htm#geom"
  >-------------------------------</a><a name="facetarea_simplex">-</a>

  qh_facetarea_simplex( 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 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)
  
  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 (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) {
    ivp_message( "qhull internal error (qh_facetarea_simplex): #points %d != dim %d -1\n", 
               i, dim);
    qh_errexit (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 (rows, dim, &nearzero);
  if (toporient)
    area= -area;
  area *= qh AREAfactor;
  trace4((qh ferr, "qh_facetarea_simplex: area=%2.2g for point p%d, toporient %d, nearzero? %d\n",
	  area, qh_pointid(apex), toporient, nearzero)); 
  return area;
} /* facetarea_simplex */

/*-<a                             href="qh-c.htm#geom"
  >-------------------------------</a><a name="facetcenter">-</a>
  
  qh_facetcenter( 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 (setT *vertices) {
  setT *points= qh_settemp (qh_setsize (vertices));
  vertexT *vertex, **vertexp;
  pointT *center;
  
  FOREACHvertex_(vertices) 
    qh_setappend (&points, vertex->point);
  center= qh_voronoi_center (qh hull_dim-1, points);
  qh_settempfree (&points);
  return center;
} /* facetcenter */

/*-<a                             href="qh-c.htm#geom"
  >-------------------------------</a><a name="findbestsharp">-</a>
  
  qh_findbestsharp( point, bestfacet, bestdist, numpart )
    find best facet on newfacet_list
    skips already visited facets (qh.visit_id) on qh.newfacet_list
    skips upperdelaunay facets unless point is outside

  returns:
    true if could be an acute angle (facets in different quadrants)
    returns bestfacet and distance to facet
    increments numpart by number of distance tests
 
  notes:
    for qh_findbest

  design:
    for all facets on qh.newfacet_list
      if two facets are in different quadrants
        set issharp
      compute distance from point to facet
      update best facet
*/
boolT qh_findbestsharp (pointT *point, facetT **bestfacet, 
           realT *bestdist, int *numpart) {
  facetT *facet;
  realT dist;
  boolT issharp = False;
  int *quadrant, k;
  
  quadrant= (int*)qh_memalloc (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 if (!issharp) {
      for (k= qh hull_dim; k--; ) {
        if (quadrant[ k] != (facet->normal[ k] > 0)) {
          issharp= True;
          break;
        }
      }
    }
    if (facet->visitid != qh visit_id) {
      qh_distplane (point, facet, &dist);
      (*numpart)++;
      if (dist > *bestdist) {
      	if (!facet->upperdelaunay || dist > qh MINoutside) {
 	  *bestdist= dist;
	  *bestfacet= facet;
	}
      }
    }
  }
  qh_memfree( quadrant, qh hull_dim * sizeof(int));
  return issharp;
} /* findbestsharp */
  
/*-<a                             href="qh-c.htm#geom"
  >-------------------------------</a><a name="findgooddist">-</a>
  
  qh_findgooddist( 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 (pointT *point, facetT *facetA, realT *distp, 
               facetT **facetlist) {
  realT bestdist= -REALmax, dist;
  facetT *neighbor, **neighborp, *bestfacet=NULL, *facet;
  boolT goodseen= False;  

  if (facetA->good) {
    zinc_(Zcheckpart);  /* calls from check_bestdist occur after print stats */
    qh_distplane (point, facetA, &bestdist);
    bestfacet= facetA;
    goodseen= True;
  }
  qh_removefacet (facetA);
  qh_appendfacet (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;
      zinc_(Zcheckpart); 
      qh_distplane (point, neighbor, &dist);
      if (dist > 0) {
        qh_removefacet (neighbor);
        qh_appendfacet (neighbor);
        if (neighbor->good) {
          goodseen= True;
          if (dist > bestdist) {
            bestdist= dist;
            bestfacet= neighbor;
          }
        }
      }
    }
  }
  if (bestfacet) {
    *distp= bestdist;
    trace2((qh ferr, "qh_findgooddist: p%d is %2.2g above good facet f%d\n",
      qh_pointid(point), bestdist, bestfacet->id));
    return bestfacet;
  }
  trace4((qh ferr, "qh_findgooddist: no good facet for p%d above f%d\n", 
      qh_pointid(point), facetA->id));
  return NULL;
}  /* findgooddist */
    
/*-<a                             href="qh-c.htm#geom"
  >-------------------------------</a><a name="getarea">-</a>
  
  qh_getarea( facetlist )
    set area of all facets in facetlist
    collect statistics

  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 (facetT *facetlist) {
  realT area;
  realT dist;
  facetT *facet;

  if (qh REPORTfreq)
    ivp_message( "computing area of each facet and volume of the convex hull\n");
  else 
    trace1((qh ferr, "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;
    facet->f.area= area= qh_facetarea (facet);
    facet->isarea= True;
    if (qh DELAUNAY) {
      if (facet->upperdelaunay == qh UPPERdelaunay)
	qh totarea += area;
    }else {
      qh totarea += area;
      qh_distplane (qh interior_point, facet, &dist);
      qh totvol += -dist * area/ qh hull_dim;
    }
    if (qh PRINTstatistics) {
      wadd_(Wareatot, area);
      wmax_(Wareamax, area);
      wmin_(Wareamin, area);
    }
  }
} /* getarea */

/*-<a                             href="qh-c.htm#geom"
  >-------------------------------</a><a name="gram_schmidt">-</a>
  
  qh_gram_schmidt( dim, row )
    implements Gram-Schmidt orthogonalization by rows

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

  notes:
    see Golub & van Loan Algorithm 6.2-2
    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(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++) {