io.c

这是一个新的用于图像处理的工具箱

  /*
    assume precise calculations in io.c with roundoff covered by qh_GEOMepsilon
    mindist may be off by qh DISTround.  Maxoutside is calculated elsewhere
    so a DISTround error may have occured.
  */
  if (qh MERGING) {
#if qh_MAXoutside
    maxdist= facet->maxoutside+qh PRINTradius+ qh DISTround;
#else
    maxdist= qh max_outside+qh PRINTradius+ qh DISTround;
#endif
    if (qh PRINTcoplanar || qh PRINTspheres)
      maxdist += qh maxmaxcoord * qh_GEOMepsilon; 
  }else 
    maxdist= 0;
  if (qh MERGING)
    mindist= -qh PRINTradius - qh DISTround;
  else
    mindist= 0;
  if (qh PRINTouter || (!qh PRINTnoplanes && !qh PRINTinner))
    qh_printfacet3geom_points(fp, points, facet, maxdist, color);
  if (qh PRINTinner || (!qh PRINTnoplanes && !qh PRINTouter &&
              maxdist - mindist > 2 * qh maxmaxcoord * qh_GEOMepsilon)) {
    for (k= 3; k--; )
      color[k]= 1.0 - color[k];
    qh_printfacet3geom_points(fp, points, facet, mindist, color);
  }
  qh_settempfree(&points);
  qh_settempfree(&vertices);
  if ((qh DOintersections || qh PRINTridges)
  && (!facet->visible || !qh NEWfacets)) {
    facet->visitid= qh visit_id;
    FOREACHneighbor_(facet) {
      if (neighbor->visitid != qh visit_id) {
	vertices= qh_setnew_delnthsorted (facet->vertices, qh hull_dim,
	                  SETindex_(facet->neighbors, neighbor), 0);
        if (qh DOintersections)
	   qh_printhyperplaneintersection(fp, facet, neighbor, vertices, black); 
        if (qh PRINTridges) {
          vertexA= SETfirst_(vertices);
          vertexB= SETsecond_(vertices);
          qh_printline3geom (fp, vertexA->point, vertexB->point, green);
        }
	qh_setfree(&vertices);
      }
    }
  }
} /* printfacet3geom_simplicial */

/*-----------------------------------------
-printfacet3vertex- print vertex->point for a 3-d facet (may be non-simplicial)
  prints number of vertices first if format == qh_PRINToff
*/
void qh_printfacet3vertex(FILE *fp, facetT *facet, int format) {
  vertexT *vertex, **vertexp;
  setT *vertices;

  vertices= qh_facet3vertex (facet);
  if (format == qh_PRINToff)
    fprintf (fp, "%d ", qh_setsize (vertices));
  FOREACHvertex_(vertices) 
    fprintf (fp, "%d ", qh_pointid(vertex->point));
  fprintf (fp, "\n");
  qh_settempfree(&vertices);
} /* printfacet3vertex */


/*----------------------------------------
-printfacet3math- print 3-d Mathematica output for a facet
  may be non-simplicial
*/
void qh_printfacet3math (FILE *fp, facetT *facet, int notfirst) {
  vertexT *vertex, **vertexp;
  setT *points, *vertices;
  pointT *point, **pointp;
  boolT firstpoint= True;
  realT dist;
  
  if (notfirst)
    fprintf(fp, ",\n");
  vertices= qh_facet3vertex (facet);
  points= qh_settemp (qh_setsize (vertices));
  FOREACHvertex_(vertices) {
    zinc_(Zdistio);
    qh_distplane(vertex->point, facet, &dist);
    point= qh_projectpoint(vertex->point, facet, dist);
    qh_setappend (&points, point);
  }
  fprintf(fp, "Polygon[{");
  FOREACHpoint_(points) {
    if (firstpoint)
      firstpoint= False;
    else
      fprintf(fp, ",\n");
    fprintf(fp, "{%10.8g, %10.8g, %10.8g}", point[0], point[1], point[2]);
  }
  FOREACHpoint_(points)
    qh_memfree (point, qh normal_size);
  qh_settempfree(&points);
  qh_settempfree(&vertices);
  fprintf(fp, "}]");
} /* printfacet3math */


/*-----------------------------------------
-printfacet4geom_nonsimplicial- print Geomview 4OFF file for a 4d nonsimplicial facet
  prints all ridges to unvisited neighbors (qh visit_id)
    must agree with printend4geom()
  prints in OFF format if DROPdim
*/
void qh_printfacet4geom_nonsimplicial(FILE *fp, facetT *facet, realT color[3]) {
  facetT *neighbor;
  ridgeT *ridge, **ridgep;
  vertexT *vertex, **vertexp;
  pointT *point;
  int k;
  realT dist;
  
  facet->visitid= qh visit_id;
  if (qh PRINTnoplanes || (facet->visible && qh NEWfacets))
    return;
  FOREACHridge_(facet->ridges) {
    neighbor= otherfacet_(ridge, facet);
    if (neighbor->visitid == qh visit_id) 
      continue;
    if (qh DOintersections)
      qh_printhyperplaneintersection(fp, facet, neighbor, ridge->vertices, color);
    else {
      if (qh DROPdim >= 0) 
	fprintf(fp, "OFF 3 1 1 # f%d\n", facet->id);
      else {
	qh printoutvar++;
	fprintf (fp, "# r%d between f%d f%d\n", ridge->id, facet->id, neighbor->id);
      }
      FOREACHvertex_(ridge->vertices) {
	zinc_(Zdistio);
	qh_distplane(vertex->point,facet, &dist);
	point=qh_projectpoint(vertex->point,facet, dist);
	for(k= 0; k < qh hull_dim; k++) {
	  if (k != qh DROPdim)
  	    fprintf(fp, "%8.4g ", point[k]);
  	}
	fprintf (fp, "\n");
	qh_memfree (point, qh normal_size);
      }
      if (qh DROPdim >= 0)
        fprintf(fp, "3 0 1 2 %8.4g %8.4g %8.4g\n", color[0], color[1], color[2]);
    }
  }
} /* printfacet4geom_nonsimplicial */


/*-----------------------------------------
-printfacet4geom_simplicial- print Geomview 4OFF file for a 4d simplicial facet
  prints triangles for unvisited neighbors (qh visit_id)
    must agree with printend4geom()
*/
void qh_printfacet4geom_simplicial(FILE *fp, facetT *facet, realT color[3]) {
  setT *vertices;
  facetT *neighbor, **neighborp;
  vertexT *vertex, **vertexp;
  int k;
  
  facet->visitid= qh visit_id;
  if (qh PRINTnoplanes || (facet->visible && qh NEWfacets))
    return;
  FOREACHneighbor_(facet) {
    if (neighbor->visitid == qh visit_id)
      continue;
    vertices= qh_setnew_delnthsorted (facet->vertices, qh hull_dim,
	                  SETindex_(facet->neighbors, neighbor), 0);
    if (qh DOintersections)
      qh_printhyperplaneintersection(fp, facet, neighbor, vertices, color);
    else {
      if (qh DROPdim >= 0) 
	fprintf(fp, "OFF 3 1 1 # ridge between f%d f%d\n",
		facet->id, neighbor->id);
      else {
	qh printoutvar++;
	fprintf (fp, "# ridge between f%d f%d\n", facet->id, neighbor->id);
      }
      FOREACHvertex_(vertices) {
	for(k= 0; k < qh hull_dim; k++) {
	  if (k != qh DROPdim)
  	    fprintf(fp, "%8.4g ", vertex->point[k]);
  	}
	fprintf (fp, "\n");
      }
      if (qh DROPdim >= 0) 
        fprintf(fp, "3 0 1 2 %8.4g %8.4g %8.4g\n", color[0], color[1], color[2]);
    }
    qh_setfree(&vertices);
  }
} /* printfacet4geom_simplicial */


/*-----------------------------------------
-printfacetNvertex_nonsimplicial- print vertices for an N-d non-simplicial facet
  triangulates each ridge to the id
*/
void qh_printfacetNvertex_nonsimplicial(FILE *fp, facetT *facet, int id) {
  vertexT *vertex, **vertexp;
  ridgeT *ridge, **ridgep;

  if (facet->visible && qh NEWfacets)
    return;
  FOREACHridge_(facet->ridges) {
    fprintf(fp, "%d ", id);
    if ((ridge->top == facet) ^ qh_ORIENTclock) {
      FOREACHvertex_(ridge->vertices)
        fprintf(fp, "%d ", qh_pointid(vertex->point));
    }else {
      FOREACHvertexreverse12_(ridge->vertices)
        fprintf(fp, "%d ", qh_pointid(vertex->point));
    }
    fprintf(fp, "\n");
  }
} /* printfacetNvertex_nonsimplicial */


/*-----------------------------------------
-printfacetNvertex_simplicial- print vertices for an N-d simplicial facet
  also prints PRINToff format for non-simplicial facets
*/
void qh_printfacetNvertex_simplicial(FILE *fp, facetT *facet, int format) {
  vertexT *vertex, **vertexp;

  if (format == qh_PRINToff)
    fprintf (fp, "%d ", qh_setsize (facet->vertices));
  if ((facet->toporient ^ qh_ORIENTclock) || !facet->simplicial) {
    FOREACHvertex_(facet->vertices)
      fprintf(fp, "%d ", qh_pointid(vertex->point));
  }else {
    FOREACHvertexreverse12_(facet->vertices)
      fprintf(fp, "%d ", qh_pointid(vertex->point));
  }
  fprintf(fp, "\n");
} /* printfacetNvertex_simplicial */


/*-----------------------------------------
-printfacetheader- prints header fields of a facet to fp
*/
void qh_printfacetheader(FILE *fp, facetT *facet) {
  pointT *point, **pointp, *furthest;
  facetT *neighbor, **neighborp;
  realT dist;

  if (facet == qh_MERGEridge) {
    fprintf (fp, " MERGEridge\n");
    return;
  }else if (facet == qh_DUPLICATEridge) {
    fprintf (fp, " DUPLICATEridge\n");
    return;
  }else if (!facet) {
    fprintf (fp, " NULLfacet\n");
    return;
  }
  fprintf(fp, "- f%d\n", facet->id);
  fprintf(fp, "    - flags:");
  if (facet->toporient) 
    fprintf(fp, " top");
  else
    fprintf(fp, " bottom");
  if (facet->simplicial)
    fprintf(fp, " simplicial");
  if (facet->upperdelaunay)
    fprintf(fp, " upperDelaunay");
  if (facet->visible)
    fprintf(fp, " visible");
  if (facet->newfacet)
    fprintf(fp, " new");
  if (facet->tested)
    fprintf(fp, " tested");
  if (!facet->good)
    fprintf(fp, " notG");
  if (facet->seen)
    fprintf(fp, " seen");
  if (facet->coplanar)
    fprintf(fp, " coplanar");
  if (facet->mergehorizon)
    fprintf(fp, " mergehorizon");
  if (facet->keepcentrum)
    fprintf(fp, " keepcentrum");
  if (facet->dupridge)
    fprintf(fp, " dupridge");
  if (facet->mergeridge && !facet->mergeridge2)
    fprintf(fp, " mergeridge1");
  if (facet->mergeridge2)
    fprintf(fp, " mergeridge2");
  if (facet->newmerge)
    fprintf(fp, " newmerge");
  if (facet->flipped) 
    fprintf(fp, " flipped");
  if (facet->degenerate)
    fprintf(fp, " degenerate");
  if (facet->redundant)
    fprintf(fp, " redundant");
  fprintf(fp, "\n");
  if (facet->isarea)
    fprintf(fp, "    - area: %2.2g\n", facet->f.area);
  else if (qh NEWfacets && facet->visible && facet->f.replace)
    fprintf(fp, "    - replacement: f%d\n", facet->f.replace->id);
  else if (facet->newfacet) {
    if (facet->f.samecycle && facet->f.samecycle != facet)
      fprintf(fp, "    - shares same visible/horizon as f%d\n", facet->f.samecycle->id);
  }else if (facet->f.newcycle)
    fprintf(fp, "    - was horizon to f%d\n", facet->f.newcycle->id);
  if (facet->nummerge)
    fprintf(fp, "    - merges: %d\n", facet->nummerge);
  qh_printpoint(fp, "    - normal: ", facet->normal);
  fprintf(fp, "    - offset: %10.7g\n", facet->offset);
  if (qh CENTERtype == qh_ASvoronoi || facet->center)
    qh_printcenter (fp, qh_PRINTfacets, "    - center: ", facet);
#if qh_MAXoutside
  if (facet->maxoutside > qh DISTround)
    fprintf(fp, "    - maxoutside: %10.7g\n", facet->maxoutside);
#endif
  if (!SETempty_(facet->outsideset)) {
    furthest= (pointT*)qh_setlast(facet->outsideset);
    if (qh_setsize (facet->outsideset) < 6) {
      fprintf(fp, "    - outside set (furthest p%d):\n", qh_pointid(furthest));
      FOREACHpoint_(facet->outsideset)
	qh_printpoint(fp, "     ", point);
    }else if (qh_setsize (facet->outsideset) < 21) {
      qh_printpoints(fp, "    - outside set:", facet->outsideset);
    }else {
      fprintf(fp, "    - outside set:  %d points.", qh_setsize(facet->outsideset));
      qh_printpoint(fp, "  Furthest", furthest);
    }
#if !qh_COMPUTEfurthest
    fprintf(fp, "    - furthest distance= %2.2g\n", facet->furthestdist);
#endif
  }
  if (!SETempty_(facet->coplanarset)) {
    furthest= (pointT*)qh_setlast(facet->coplanarset);
    if (qh_setsize (facet->coplanarset) < 6) {
      fprintf(fp, "    - coplanar set (furthest p%d):\n", qh_pointid(furthest));
      FOREACHpoint_(facet->coplanarset)
	qh_printpoint(fp, "     ", point);
    }else if (qh_setsize (facet->coplanarset) < 21) {
      qh_printpoints(fp, "    - coplanar set:", facet->coplanarset);
    }else {
      fprintf(fp, "    - coplanar set:  %d points.", qh_setsize(facet->coplanarset));
      qh_printpoint(fp, "  Furthest", furthest);
    }
    zinc_(Zdistio);
    qh_distplane (furthest, facet, &dist);
    fprintf(fp, "      furthest distance= %2.2g\n", dist);
  }
  qh_printvertices (fp, "    - vertices:", facet->vertices);
  fprintf(fp, "    - neighboring facets: ");
  FOREACHneighbor_(facet) {
    if (neighbor == qh_MERGEridge)
      fprintf(fp, " MERGE");
    else if (neighbor == qh_DUPLICATEridge)
      fprintf(fp, " DUP");
    else
      fprintf(fp, " f%d", neighbor->id);
  }
  fprintf(fp, "\n");
} /* printfacetheader */


/*-----------------------------------------
-printfacetridges- prints ridges of a facet to fp
  ridges printed in neighbor order
  assumes the ridges exist
*/
void qh_printfacetridges(FILE *fp, facetT *facet) {
  facetT *neighbor, **neighborp;
  ridgeT *ridge, **ridgep, *firstridge;
  int numridges= 0;


  if (facet->visible && qh NEWfacets) {
    fprintf(fp, "    - ridges (ids may be garbage):");
    FOREACHridge_(facet->ridges)
      fprintf(fp, " r%d", ridge->id);
    fprintf(fp, "\n");
  }else {
    fprintf(fp, "    - ridges:\n");
    FOREACHridge_(facet->ridges)
      ridge->seen= False;
    if (qh hull_dim == 3) {
      ridge= firstridge= SETfirst_(facet->ridges);
      while (ridge && !ridge->seen) {
	ridge->seen= True;
	qh_printridge(fp, ridge);
	numridges++;