poly2.c

Quick hull implementation

      if (testouter) {
#if qh_MAXoutside
	maxoutside= facet->maxoutside + 2* qh DISTround;
	/* one DISTround to actual point and another to computed point */
#endif
      }
      FORALLpoints {
	if (point != qh GOODpointp)
	  qh_check_point (point, facet, &maxoutside, &maxdist, &errfacet1, &errfacet2);
      }
      FOREACHpoint_(qh other_points) {
	if (point != qh GOODpointp)
	  qh_check_point (point, facet, &maxoutside, &maxdist, &errfacet1, &errfacet2);
      }
    }
    if (maxdist > qh outside_err) {
      fprintf( qh ferr, "qhull precision error (qh_check_points): a coplanar point is %6.2g from convex hull.  The maximum value (qh.outside_err) is %6.2g\n",
                maxdist, qh outside_err );
      qh_errexit2( qh_ERRprec, errfacet1, errfacet2 );
    }else if (errfacet1 && qh outside_err > REALmax/2)
        qh_errexit2( qh_ERRprec, errfacet1, errfacet2 );
    else if (errfacet1)
        ;  /* the error was logged to qh.ferr but does not effect the output */
    trace0((qh ferr, "qh_check_points: max distance outside %2.2g\n", maxdist));
  }
} /* check_points */


/*-<a                             href="qh-poly.htm#TOC"
  >-------------------------------</a><a name="checkconvex">-</a>
  
  qh_checkconvex( facetlist, fault )
    check that each ridge in facetlist is convex
    fault = qh_DATAfault if reporting errors
          = qh_ALGORITHMfault otherwise

  returns:
    counts Zconcaveridges and Zcoplanarridges
    errors if concaveridge or if merging an coplanar ridge

  note:
    if not merging, 
      tests vertices for neighboring simplicial facets
    else if ZEROcentrum, 
      tests vertices for neighboring simplicial   facets
    else 
      tests centrums of neighboring facets

  design:
    for all facets
      report flipped facets
      if ZEROcentrum and simplicial neighbors
        test vertices for neighboring simplicial facets
      else
        test centrum against all neighbors 
*/
void qh_checkconvex(facetT *facetlist, int fault) {
  facetT *facet, *neighbor, **neighborp, *errfacet1=NULL, *errfacet2=NULL;
  vertexT *vertex;
  realT dist;
  pointT *centrum;
  boolT waserror= False, centrum_warning= False, tempcentrum= False, allsimplicial;
  int neighbor_i;

  trace1((qh ferr, "qh_checkconvex: check all ridges are convex\n"));
  if (!qh RERUN) {
    zzval_(Zconcaveridges)= 0;
    zzval_(Zcoplanarridges)= 0;
  }
  FORALLfacet_(facetlist) {
    if (facet->flipped) {
      qh_precision ("flipped facet");
      fprintf (qh ferr, "qhull precision error: f%d is flipped (interior point is outside)\n",
	       facet->id);
      errfacet1= facet;
      waserror= True;
      continue;
    }
    if (qh MERGING && (!qh ZEROcentrum || !facet->simplicial || facet->tricoplanar))
      allsimplicial= False;
    else {
      allsimplicial= True;
      neighbor_i= 0;
      FOREACHneighbor_(facet) {
        vertex= SETelemt_(facet->vertices, neighbor_i++, vertexT);
	if (!neighbor->simplicial || neighbor->tricoplanar) {
	  allsimplicial= False;
	  continue;
	}
        qh_distplane (vertex->point, neighbor, &dist);
        if (dist > -qh DISTround) {
	  if (fault == qh_DATAfault) {
            qh_precision ("coplanar or concave ridge");
	    fprintf (qh ferr, "qhull precision error: initial simplex is not convex. Distance=%.2g\n", dist);
	    qh_errexit(qh_ERRsingular, NULL, NULL);
	  }
          if (dist > qh DISTround) {
            zzinc_(Zconcaveridges);
            qh_precision ("concave ridge");
            fprintf (qh ferr, "qhull precision error: f%d is concave to f%d, since p%d (v%d) is %6.4g above\n",
              facet->id, neighbor->id, qh_pointid(vertex->point), vertex->id, dist);
            errfacet1= facet;
            errfacet2= neighbor;
            waserror= True;
          }else if (qh ZEROcentrum) {
            if (dist > 0) {     /* qh_checkzero checks that dist < - qh DISTround */
              zzinc_(Zcoplanarridges); 
              qh_precision ("coplanar ridge");
              fprintf (qh ferr, "qhull precision error: f%d is clearly not convex to f%d, since p%d (v%d) is %6.4g above\n",
                facet->id, neighbor->id, qh_pointid(vertex->point), vertex->id, dist);
              errfacet1= facet;
              errfacet2= neighbor;
              waserror= True;
	    }
	  }else {              
            zzinc_(Zcoplanarridges);
            qh_precision ("coplanar ridge");
            trace0((qh ferr, "qhull precision error: f%d may be coplanar to f%d, since p%d (v%d) is within %6.4g during p%d\n",
              facet->id, neighbor->id, qh_pointid(vertex->point), vertex->id, dist, qh furthest_id));
          }
        }
      }
    }
    if (!allsimplicial) {
      if (qh CENTERtype == qh_AScentrum) {
        if (!facet->center)
          facet->center= qh_getcentrum (facet);
        centrum= facet->center;
      }else {
	if (!centrum_warning && (!facet->simplicial || facet->tricoplanar)) {
	   centrum_warning= True;
	   fprintf (qh ferr, "qhull note: recomputing centrums for convexity test.  This may lead to false, precision errors.\n");
	}
        centrum= qh_getcentrum(facet);
        tempcentrum= True;
      }
      FOREACHneighbor_(facet) {
	if (qh ZEROcentrum && facet->simplicial && neighbor->simplicial)
	  continue;
	if (facet->tricoplanar || neighbor->tricoplanar)
	  continue;
        zzinc_(Zdistconvex);
        qh_distplane (centrum, neighbor, &dist);
        if (dist > qh DISTround) {
          zzinc_(Zconcaveridges);
          qh_precision ("concave ridge");
          fprintf (qh ferr, "qhull precision error: f%d is concave to f%d.  Centrum of f%d is %6.4g above f%d\n",
            facet->id, neighbor->id, facet->id, dist, neighbor->id);
          errfacet1= facet;
          errfacet2= neighbor;
          waserror= True;
	}else if (dist >= 0.0) {   /* if arithmetic always rounds the same,
				     can test against centrum radius instead */
          zzinc_(Zcoplanarridges);
          qh_precision ("coplanar ridge");
          fprintf (qh ferr, "qhull precision error: f%d is coplanar or concave to f%d.  Centrum of f%d is %6.4g above f%d\n",
            facet->id, neighbor->id, facet->id, dist, neighbor->id);
	  errfacet1= facet;
	  errfacet2= neighbor;
	  waserror= True;
        }
      }
      if (tempcentrum)
        qh_memfree(centrum, qh normal_size);
    }
  }
  if (waserror && !qh FORCEoutput)
    qh_errexit2 (qh_ERRprec, errfacet1, errfacet2);
} /* checkconvex */


/*-<a                             href="qh-poly.htm#TOC"
  >-------------------------------</a><a name="checkfacet">-</a>
  
  qh_checkfacet( facet, newmerge, waserror )
    checks for consistency errors in facet
    newmerge set if from merge.c

  returns:
    sets waserror if any error occurs

  checks:
    vertex ids are inverse sorted
    unless newmerge, at least hull_dim neighbors and vertices (exactly if simplicial)
    if non-simplicial, at least as many ridges as neighbors
    neighbors are not duplicated
    ridges are not duplicated
    in 3-d, ridges=verticies
    (qh.hull_dim-1) ridge vertices
    neighbors are reciprocated
    ridge neighbors are facet neighbors and a ridge for every neighbor
    simplicial neighbors match facetintersect
    vertex intersection matches vertices of common ridges 
    vertex neighbors and facet vertices agree
    all ridges have distinct vertex sets

  notes:  
    uses neighbor->seen

  design:
    check sets
    check vertices
    check sizes of neighbors and vertices
    check for qh_MERGEridge and qh_DUPLICATEridge flags
    check neighbor set
    check ridge set
    check ridges, neighbors, and vertices
*/
void qh_checkfacet(facetT *facet, boolT newmerge, boolT *waserrorp) {
  facetT *neighbor, **neighborp, *errother=NULL;
  ridgeT *ridge, **ridgep, *errridge= NULL, *ridge2;
  vertexT *vertex, **vertexp;
  unsigned previousid= INT_MAX;
  int numneighbors, numvertices, numridges=0, numRvertices=0;
  boolT waserror= False;
  int skipA, skipB, ridge_i, ridge_n, i;
  setT *intersection;

  if (facet->visible) {
    fprintf (qh ferr, "qhull internal error (qh_checkfacet): facet f%d is on the visible_list\n",
      facet->id);
    qh_errexit (qh_ERRqhull, facet, NULL);
  }
  if (!facet->normal) {
    fprintf (qh ferr, "qhull internal error (qh_checkfacet): facet f%d does not have  a normal\n",
      facet->id);
    waserror= True;
  }
  qh_setcheck (facet->vertices, "vertices for f", facet->id);
  qh_setcheck (facet->ridges, "ridges for f", facet->id);
  qh_setcheck (facet->outsideset, "outsideset for f", facet->id);
  qh_setcheck (facet->coplanarset, "coplanarset for f", facet->id);
  qh_setcheck (facet->neighbors, "neighbors for f", facet->id);
  FOREACHvertex_(facet->vertices) {
    if (vertex->deleted) {
      fprintf(qh ferr, "qhull internal error (qh_checkfacet): deleted vertex v%d in f%d\n", vertex->id, facet->id);
      qh_errprint ("ERRONEOUS", NULL, NULL, NULL, vertex);
      waserror= True;
    }
    if (vertex->id >= previousid) {
      fprintf(qh ferr, "qhull internal error (qh_checkfacet): vertices of f%d are not in descending id order at v%d\n", facet->id, vertex->id);
      waserror= True;
      break;
    }
    previousid= vertex->id;
  }
  numneighbors= qh_setsize(facet->neighbors);
  numvertices= qh_setsize(facet->vertices);
  numridges= qh_setsize(facet->ridges);
  if (facet->simplicial) {
    if (numvertices+numneighbors != 2*qh hull_dim 
    && !facet->degenerate && !facet->redundant) {
      fprintf(qh ferr, "qhull internal error (qh_checkfacet): for simplicial facet f%d, #vertices %d + #neighbors %d != 2*qh hull_dim\n", 
                facet->id, numvertices, numneighbors);
      qh_setprint (qh ferr, "", facet->neighbors);
      waserror= True;
    }
  }else { /* non-simplicial */
    if (!newmerge 
    &&(numvertices < qh hull_dim || numneighbors < qh hull_dim)
    && !facet->degenerate && !facet->redundant) {
      fprintf(qh ferr, "qhull internal error (qh_checkfacet): for facet f%d, #vertices %d or #neighbors %d < qh hull_dim\n",
         facet->id, numvertices, numneighbors);
       waserror= True;
    }
    /* in 3-d, can get a vertex twice in an edge list, e.g., RBOX 1000 s W1e-13 t995849315 D2 | QHULL d Tc Tv TP624 TW1e-13 T4 */
    if (numridges < numneighbors
    ||(qh hull_dim == 3 && numvertices > numridges && !qh NEWfacets)
    ||(qh hull_dim == 2 && numridges + numvertices + numneighbors != 6)) {
      if (!facet->degenerate && !facet->redundant) {
	fprintf(qh ferr, "qhull internal error (qh_checkfacet): for facet f%d, #ridges %d < #neighbors %d or (3-d) > #vertices %d or (2-d) not all 2\n",
	    facet->id, numridges, numneighbors, numvertices);
	waserror= True;
      }
    }
  }
  FOREACHneighbor_(facet) {
    if (neighbor == qh_MERGEridge || neighbor == qh_DUPLICATEridge) {
      fprintf(qh ferr, "qhull internal error (qh_checkfacet): facet f%d still has a MERGE or DUP neighbor\n", facet->id);
      qh_errexit (qh_ERRqhull, facet, NULL);
    }
    neighbor->seen= True;
  }
  FOREACHneighbor_(facet) {
    if (!qh_setin(neighbor->neighbors, facet)) {
      fprintf(qh ferr, "qhull internal error (qh_checkfacet): facet f%d has neighbor f%d, but f%d does not have neighbor f%d\n",
	      facet->id, neighbor->id, neighbor->id, facet->id);
      errother= neighbor;
      waserror= True;
    }
    if (!neighbor->seen) {
      fprintf(qh ferr, "qhull internal error (qh_checkfacet): facet f%d has a duplicate neighbor f%d\n",
	      facet->id, neighbor->id);
      errother= neighbor;
      waserror= True;
    }    
    neighbor->seen= False;
  }
  FOREACHridge_(facet->ridges) {
    qh_setcheck (ridge->vertices, "vertices for r", ridge->id);
    ridge->seen= False;
  }
  FOREACHridge_(facet->ridges) {
    if (ridge->seen) {
      fprintf(qh ferr, "qhull internal error (qh_checkfacet): facet f%d has a duplicate ridge r%d\n",
	      facet->id, ridge->id);
      errridge= ridge;
      waserror= True;
    }    
    ridge->seen= True;
    numRvertices= qh_setsize(ridge->vertices);
    if (numRvertices != qh hull_dim - 1) {
      fprintf(qh ferr, "qhull internal error (qh_checkfacet): ridge between f%d and f%d has %d vertices\n", 
                ridge->top->id, ridge->bottom->id, numRvertices);
      errridge= ridge;
      waserror= True;
    }
    neighbor= otherfacet_(ridge, facet);
    neighbor->seen= True;
    if (!qh_setin(facet->neighbors, neighbor)) {
      fprintf(qh ferr, "qhull internal error (qh_checkfacet): for facet f%d, neighbor f%d of ridge r%d not in facet\n",
           facet->id, neighbor->id, ridge->id);
      errridge= ridge;
      waserror= True;
    }
  }
  if (!facet->simplicial) {
    FOREACHneighbor_(facet) {
      if (!neighbor->seen) {
        fprintf(qh ferr, "qhull internal error (qh_checkfacet): facet f%d does not have a ridge for neighbor f%d\n",
	      facet->id, neighbor->id);
	errother= neighbor;
        waserror= True;
      }
      intersection= qh_vertexintersect_new(facet->vertices, neighbor->vertices);
      qh_settemppush (intersection);
      FOREACHvertex_(facet->vertices) {
	vertex->seen= False;
	vertex->seen2= False;
      }
      FOREACHvertex_(intersection)
	vertex->seen= True;
      FOREACHridge_(facet->ridges) {
	if (neighbor != otherfacet_(ridge, facet))
	    continue;
	FOREACHvertex_(ridge->vertices) {
	  if (!vertex->seen) {
	    fprintf (qh ferr, "qhull internal error (qh_checkfacet): vertex v%d in r%d not in f%d intersect f%d\n",
  	          vertex->id, ridge->id, facet->id, neighbor->id);
	    qh_errexit (qh_ERRqhull, facet, ridge);
	  }
	  vertex->seen2= True;
	}
      }
      if (!newmerge) {
	FOREACHvertex_(intersection) {
	  if (!vertex->seen2) {
	    if (qh IStracing >=3 || !qh MERGING) {
	      fprintf (qh ferr, "qhull precision error (qh_checkfacet): vertex v%d in f%d intersect f%d but\n\
 not in a ridge.  This is ok under merging.  Last point was p%d\n",
		     vertex->id, facet->id, neighbor->id, qh furthest_id);
	      if (!qh FORCEoutput && !qh MERGING) {
		qh_errprint ("ERRONEOUS", facet, neighbor, NULL, vertex);
		if (!qh MERGING)
		  qh_errexit (qh_ERRqhull, NULL, NULL);
	      }
	    }
	  }
	}
      }