This is polymerge.c in view mode; [Download] [Up]
/*
* Squeeze OFF files.
* Merges collinear edges and coplanar faces.
*/
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <stdlib.h> /* for qsort(), malloc() */
extern char *strdup(const char *);
#define VSUB(a, b, dst) (dst)->x = (a)->x - (b)->x, \
(dst)->y = (a)->y - (b)->y, \
(dst)->z = (a)->z - (b)->z
#define VDOT(a, b) (a)->x*(b)->x + (a)->y*(b)->y + (a)->z*(b)->z
#define VCROSS(a, b, dst) (dst)->x = (a)->y * (b)->z - (a)->z * (b)->y, \
(dst)->y = (a)->z * (b)->x - (a)->x * (b)->z, \
(dst)->z = (a)->x * (b)->y - (a)->y * (b)->x
float EPS_LINEAR= 0.0001; /* Max unit-vector cross product for collinearity */
float EPS_NORMAL= 0.0009; /* Max sin^2 theta between normals for coplanarity */
float EPS_POINT= 0.00001;/* Max distance between coincident vertices */
#define F_VERT 0x1
#define F_EDGE 0x2
#define F_FACE 0x4
#define F_EVOLVER 0x8
typedef struct vertex V;
typedef struct face F;
typedef struct faceedge Fe;
typedef struct edge E;
typedef struct point {
float x, y, z, w;
} P;
struct vertex {
P p;
int ref;
int index;
float *more;
};
struct faceedge {
F *face;
Fe *prev, *next; /* double links in face-edges of this face */
E *edge;
int sign; /* direction w.r.t. edge */
P n; /* adopted surface normal at this edge */
Fe *elink;
};
struct face {
Fe *fedges;
char *affix;
};
struct edge {
V *v[2]; /* vertex pointers */
E *link; /* hash table link */
P to; /* normalized edge direction vector */
Fe *feds; /* faceedge's on this edge */
int index; /* for numbering edges for .fe format */
};
int nv, nf, tossedf, nhash;
V *Vs; /* All vertices */
F *Fs; /* All faces */
E **ehash; /* Array of hash-table head pointers */
int nhash;
int flags = F_VERT|F_EDGE|F_FACE;
int debug = 0;
int vdim = 3;
Fe *fedge( F *f, V *v0, V *v1 );
#define New(t) (t *) malloc(sizeof(t))
#define NewN(t, N) (t *) malloc((N)*sizeof(t))
main(argc, argv)
char *argv[];
{
register int i;
int tnv, tne, any;
int binary = 0, vertsize = 0;
char *C = "", *N = "", *four = "";
extern int optind;
extern char *optarg;
while((i = getopt(argc, argv, "bdv:e:f:VEF")) != EOF) {
switch(i) {
case 'b': flags |= F_EVOLVER; break;
case 'd': debug = 1; break;
case 'v': EPS_POINT = atof(optarg); break;
case 'e': EPS_LINEAR = atof(optarg); EPS_LINEAR *= EPS_LINEAR; break;
case 'f': EPS_NORMAL = atof(optarg); EPS_NORMAL *= EPS_NORMAL; break;
case 'V': flags &= ~F_VERT; break;
case 'E': flags &= ~F_EDGE; break;
case 'F': flags &= ~F_FACE; break;
default:
fprintf(stderr, "\
Usage: polymerge [-v vertex_thresh] [-e edge_thresh] [-f face_thresh]\n\
[-V][-E][-F][-d][-b] [infile.off]\n\
Merges coincident vertices, collinear edges, coplanar faces of an OFF object.\n\
-v vertex_thresh: max separation for \"coincident\" vertices (default %g)\n\
-e edge_thresh : max sin(theta) for \"collinear\" edges (default %g)\n\
-f face_thresh : max sin(theta) for \"coplanar\" facet normals (default %.4g)\n\
-V -E -F : don't try to merge vertices/edges/faces\n\
-b : create output in evolver .fe format\n\
-d : debug\n",
EPS_POINT, sqrt(EPS_LINEAR), sqrt(EPS_NORMAL));
exit(1);
}
}
if(optind < argc) {
if(freopen(argv[argc-1], "r", stdin) == NULL) {
fprintf(stderr, "polymerge: can't open input: ");
perror(argv[argc-1]);
exit(1);
}
}
vdim = 3;
for(;;) {
switch( fnextc(stdin, 0) ) {
case 'N': N = "N"; vertsize += 3; goto swallow;
case 'C': C = "C"; vertsize += 4; goto swallow;
case '4':
four = "4";
vdim = 4;
flags &= ~(F_EDGE|F_FACE);
goto swallow;
swallow:
case '{': case '=': fgetc(stdin); break;
case 'a': /* Appearance? Silently swallow keyword */
fexpectstr(stdin, "appearance");
if(fnextc(stdin, 0) == '{') {
int c, brack = 0;
do {
if((c = fgetc(stdin)) == '{') brack++;
else if(c == '}') brack--;
} while(brack > 0 && c != EOF);
}
break;
default: goto done;
}
}
done:
fexpecttoken(stdin, "OFF"); /* swallow optional OFF prefix */
if(fexpecttoken(stdin, "BINARY") == 0) {
binary = 1;
fnextc(stdin, 1);
fgetc(stdin);
}
if(fgetni(stdin, 1, &nv, 0) <= 0 ||
fgetni(stdin, 1, &nf, 0) <= 0 ||
fgetni(stdin, 1, &nhash, 0) <= 0 || nv <= 0 || nf <= 0) {
fprintf(stderr, "polymerge: Input not an OFF file.\n");
exit(1);
}
nhash = (nf + nv + 4) / 2;
if(nhash < 16) nhash = 17;
if(!(nhash & 1)) nhash++;
ehash = NewN(E *, nhash);
bzero((char *)ehash, nhash*sizeof(E *));
Vs = NewN(V, nv);
Fs = NewN(F, nf);
for(i = 0; i < nv; i++) {
V *vp = &Vs[i];
vp->p.w = 1;
if(fgetnf(stdin, vdim, &vp->p, binary) < vdim) {
badvert:
fprintf(stderr, "polymerge: error reading vertex %d/%d\n", i,nv);
exit(1);
}
if(vertsize) {
Vs[i].more = NewN(float, vertsize);
if(fgetnf(stdin, vertsize, Vs[i].more, binary) < vertsize)
goto badvert;
}
Vs[i].ref = 0;
Vs[i].index = i;
}
/*
* Combine vertices
*/
if(flags & F_VERT)
vmerge();
/*
* Load faces
* The fedge() calls here assume we've already merged all appropriate
* vertices.
*/
tossedf = 0;
for(i = 0; i < nf; i++) {
register F *fp = &Fs[i];
register Fe *fe;
int nfv, k, c;
int v, ov, v0;
Fe head;
char *cp;
char aff[512];
if(fgetni(stdin, 1, &nfv, binary) <= 0 || nfv <= 0) {
fprintf(stderr, "polymerge: error reading face %d/%d\n",
i, nf);
exit(1);
}
head.prev = head.next = &head;
fp->fedges = &head;
fgetni(stdin, 1, &v, binary);
if(v < 0 || v >= nv) {
fprintf(stderr, "polymerge: bad vertex %d on face %d\n", v, i);
exit(1);
}
v0 = ov = Vs[v].index; /* Use common vertex if merged */
for(k = 1; k < nfv; k++, ov = v) {
fgetni(stdin, 1, &v, binary);
if(v < 0 || v >= nv) {
fprintf(stderr, "polymerge: bad vertex %d on face %d\n", v, i);
exit(1);
}
v = Vs[v].index; /* Use common vertex if merged */
if(ov == v)
continue;
head.prev->next = fe = fedge(fp, &Vs[ov], &Vs[v]);
fe->prev = head.prev;
fe->next = &head;
head.prev = fe;
}
if(v != v0) {
fe = fedge(fp, &Vs[v], &Vs[v0]);
head.prev->next = fe;
fe->prev = head.prev;
fe->next = &head;
head.prev = fe;
}
head.next->prev = head.prev;
if(head.next == &head) {
/*
* Degenerate face here
*/
fp->fedges = NULL;
tossedf++;
debug && printf("# Face %d degenerate already\n", i);
} else {
head.prev->next = fp->fedges = head.next;
}
if(binary) {
int nfc;
float c;
fgetni(stdin, 1, &nfc, binary);
if(nfc > 4) {
fprintf(stderr, "Bad face color count 0x%x", nfc);
exit(1);
}
for(cp = aff; --nfc >= 0; cp += strlen(cp)) {
fgetnf(stdin, 1, &c, binary);
sprintf(cp, " %g", c);
}
} else {
(void) fnextc(stdin, 1);
for(cp = aff; (c = getc(stdin)) != EOF && c != '\n' && cp < &aff[511]; )
*cp++ = c;
}
*cp = '\0';
fp->affix = (cp > aff) ? strdup(aff) : "";
}
/*
* Compute face normals -- actually corner normals, since we avoid
* assuming faces are planar. As a side effect, we detect & join
* collinear edges.
*/
if(flags & F_EDGE) {
for(i = 0; i < nf; i++) {
normalize(&Fs[i]);
}
}
/*
* Locate edges bounding faces with the same normal direction.
*/
if(flags & F_FACE) {
do {
any = 0;
for(i = 0; i < nhash; i++) {
E *e;
for(e = ehash[i]; e != NULL; e = e->link) {
register Fe *fe, *fee;
for(fe = e->feds; fe != NULL; fe = fe->elink) {
for(fee = fe->elink; fee != NULL; fee = fee->elink) {
register float dn;
dn = VDOT(&fe->n, &fee->n);
if(fabs(1 - dn*dn) < EPS_NORMAL) {
/*
* OK, merge fee into fe
*/
femerge(fe, fee);
any++;
goto another;
}
}
}
another: ;
}
}
debug && printf("# %d faces merged this pass.\n", any);
} while(any);
}
/*
* Scan for unused edges.
*/
for(i = 0; i < nv; i++)
Vs[i].ref = 0;
tne = 0;
for(i = 0; i < nhash; i++) {
E *e;
for(e = ehash[i]; e != NULL; e = e->link) {
if(e->feds != NULL) {
e->v[0]->ref++;
e->v[1]->ref++;
e->index = ++tne;
}
}
}
/*
* Renumber used vertices.
*/
if(flags & 1) {
tnv = 0;
for(i = 0; i < nv; i++)
Vs[i].index = Vs[i].ref ? tnv++ : -i-1;
} else {
tnv = nv;
}
if (flags & F_EVOLVER) /* Produce Brakke's evolver .fe format */
{
int j=0;
if (vdim == 4) printf("space_dimension 4\n");
printf("vertices\n");
for(i = 0; i < nv; i++) {
register V *v = &Vs[i];
int k;
if(v->ref || !(flags & 1)) {
v->index = ++j;
printf("%d\t%g %g %g", v->index, v->p.x, v->p.y, v->p.z);
if(vdim == 4) printf(" %g", v->p.w);
if(vertsize) {
printf(" ");
for(k = 0; k < vertsize; k++)
printf(" %g", v->more[k]);
}
printf("\n");
}
}
printf("\nedges\n");
for(i = 0; i < nhash; i++) {
E *e;
for(e = ehash[i]; e != NULL; e = e->link)
if(e->feds != NULL)
printf("%d\t%d %d\n", e->index,e->v[0]->index,e->v[1]->index);
}
printf("\nfaces\n");
j=0;
for(i=0; i<nf; i++) {
register Fe *fe, *fee;
int k, nfv;
fe = Fs[i].fedges;
if(fe == NULL)
continue;
for(fee = fe, k = 1; (fee = fee->next) != fe; k++)
;
if ((nfv = k)<3)
continue; /* don't print faces of less than 3 sides */
printf("%d", ++j);
for(fee = fe, k = nfv; --k >= 0; fee = fee->next)
printf(" %d", (1-2*fee->sign)*fee->edge->index);
printf("\n");
}
}
else /* Produce OFF format */
{
printf("%s%s%sOFF\n%d %d %d\n", C, N, four, tnv, nf - tossedf, tne);
for(i = 0; i < nv; i++) {
register V *v = &Vs[i];
int k;
if(v->ref || !(flags & F_VERT)) {
printf("%g %g %g", v->p.x, v->p.y, v->p.z);
if(vdim == 4) printf(" %g", v->p.w);
if(vertsize) {
printf(" ");
for(k = 0; k < vertsize; k++)
printf(" %g", v->more[k]);
}
if(debug)
printf("\t# %d [%d] #%d", v->index, v->ref, i);
printf("\n");
}
}
printf("\n");
/* ho hum */
for(i=0; i<nf; i++) {
register Fe *fe, *fee;
int k, nfv;
fe = Fs[i].fedges;
if(fe == NULL)
continue;
for(fee = fe, k = 1; (fee = fee->next) != fe; k++)
;
nfv = k;
printf("%d", nfv);
for(fee = fe, k = nfv; --k >= 0; fee = fee->next)
printf(" %d", fee->edge->v[fee->sign]->index);
printf("\t%s", Fs[i].affix);
if(debug) {
printf(" #");
for(fee = fe, k = nfv; --k >= 0; fee = fee->next)
printf(" %d", fee->edge->v[fee->sign] - Vs);
}
printf("\n");
}
}
exit(1);
}
/*
* Add a new faceedge
*/
Fe *
fedge(F *f, V *v0, V *v1)
{
Fe *fe;
E *e;
int t;
int i0, i1;
float r;
fe = New(Fe);
fe->face = f;
fe->sign = 0;
i0 = v0->index;
i1 = v1->index;
if(i0 > i1) {
V *tv = v0;
v0 = v1;
v1 = tv;
i0 = v0->index;
i1 = v1->index;
fe->sign = 1;
}
t = (unsigned long)(v0->index + v1->index + v0->index*v1->index) % nhash;
/* Symmetric hash function */
for(e = ehash[t]; e != NULL; e = e->link)
if(e->v[0]->index == i0 && e->v[1]->index == i1)
goto gotit;
e = New(E);
e->v[0] = v0;
e->v[1] = v1;
e->feds = NULL;
VSUB(&v0->p, &v1->p, &e->to);
r = VDOT(&e->to, &e->to);
if(r != 0) {
r = 1/sqrt(r);
e->to.x *= r; e->to.y *= r; e->to.z *= r;
} else
debug && printf("# Coincident: %d == %d [%g %g %g]\n", v0->index, v1->index, v0->p.x,v0->p.y,v0->p.z);
e->link = ehash[t];
ehash[t] = e;
gotit:
fe->edge = e;
fe->elink = e->feds;
e->feds = fe;
return fe;
}
/*
* Remove a faceedge from its edge list
*/
unfedge(fe)
register Fe *fe;
{
register Fe **fepp;
for(fepp = &fe->edge->feds; *fepp != NULL; fepp = &(*fepp)->elink) {
if(*fepp == fe) {
*fepp = fe->elink;
break;
}
}
free(fe);
}
/*
* Merge two faces
* We delete these face-edges from both faces
*/
femerge(fe1, fe2)
register Fe *fe1, *fe2;
{
F *f1, *f2;
register Fe *tfe;
if(fe1->face == fe2->face) {
debug && printf("# Merging two edges of face %d -- tossing it.\n",
fe1->face - Fs);
deface(fe1->face);
return;
}
if(fe1->sign == fe2->sign) {
/*
* Messy. To merge these, we need to reverse all the links
* in one of the two faces.
*/
Fe *xfe;
tfe = fe2;
do {
tfe->sign ^= 1;
xfe = tfe->next;
tfe->next = tfe->prev;
tfe->prev = xfe;
tfe = xfe;
} while(tfe != fe2);
}
fe1->prev->next = fe2->next;
fe2->next->prev = fe1->prev;
fe1->next->prev = fe2->prev;
fe2->prev->next = fe1->next;
f1 = fe1->face;
f2 = fe2->face;
if(f1->fedges == fe1)
f1->fedges = fe2->next;
if(debug)
printf("# Merged face %d into %d (vertices %d %d) n1.n2 %g\n",
f2-Fs, f1-Fs,
fe1->edge->v[fe1->sign]->index,
fe1->edge->v[1 - fe1->sign]->index,
VDOT(&fe1->n, &fe2->n));
tfe = fe2->next;
if(tfe == NULL) {
fprintf(stderr, "polymerge: face f2 already deleted?\n");
} else {
do {
tfe->face = f1;
} while((tfe = tfe->next) != fe1->next);
}
f2->fedges = NULL;
tossedf++;
unfedge(fe1);
unfedge(fe2);
/*
* Join collinear edges, recompute normals (might have changed a bit).
*/
normalize(f1);
}
#define PRETTY(x) ((int)(x) - 0x10000000)
fecheck(fe)
Fe *fe;
{
register Fe *fee;
int ne;
int onface = 0;
register F *f;
register E *e;
if(fe == NULL)
return;
f = fe->face;
fprintf(stderr,"0x%x: on face %d (%x); ", fe, f - Fs, f);
fee = fe;
do {
fprintf(stderr," %s%x[%d%s%d] ", (f->fedges == fee) ? "*" : "",
PRETTY(fee),
fee->edge->v[0]->index,
fee->sign ? "<-" : "->",
fee->edge->v[1]->index);
if(fee->face != f)
fprintf(stderr," Fe %x: face %d (%x) != %x\n",
PRETTY(fee), fee->face - Fs, fee->face, f);
if(fee->next->prev != fee)
fprintf(stderr," Fe %x: next %x next->prev %x\n",
PRETTY(fee), PRETTY(fee->next), PRETTY(fee->next->prev));
e = fee->edge;
fee = fee->next;
} while(fee != fe);
fprintf(stderr, "\n");
}
E *
echeck(v0, v1)
register int v0, v1;
{
register E *e;
int t = (v0 + v1 + v0*v1) % nhash;
/* Symmetric hash function */
if(v0 > v1) v0 ^= v1, v1 ^= v0, v0 ^= v1;
for(e = ehash[t]; e != NULL; e = e->link) {
if(e->v[0] == &Vs[v0] && e->v[1] == &Vs[v1]) {
register Fe *fe;
fprintf(stderr, "E 0x%x %d-%d (%d-%d) %x...\n",
e, v0,v1, e->v[0]->index, e->v[1]->index, PRETTY(e->feds));
for(fe = e->feds; fe != NULL; fe = fe->elink) {
fecheck(fe);
}
}
}
return e;
}
normalize(f)
F *f;
{
register Fe *fe;
fe = f->fedges;
if(fe == NULL)
return;
if(fe->prev == fe->next) {
debug && printf("# Face %d already degenerate -- tossing it.\n", f - Fs);
deface(f);
return;
}
do { /* loop over edges on this face */
P n;
float r;
for(;;) {
register Fe *fn, *fee;
register P *pp, *qp;
pp = &fe->edge->to;
fee = fe->next;
qp = &fee->edge->to;
VCROSS(pp, qp, &n);
r = VDOT(&n, &n);
if(r > EPS_LINEAR)
break;
/*
* Join collinear edges; produce a new edge
*/
fn = fedge(f, fe->edge->v[fe->sign], fee->edge->v[1-fee->sign]);
fee->next->prev = fe->prev->next = fn;
fn->next = fee->next;
fn->prev = fe->prev;
if(fe == f->fedges || fee == f->fedges)
f->fedges = fn; /* preserve headiness */
unfedge(fee);
if(fee != fe)
unfedge(fe);
if(fn->prev == fn->next) {
/*
* This face became degenerate -- toss it.
*/
debug && printf("# degenerate face %d\n", f - Fs);
deface(f);
return;
}
fe = fn;
}
r = 1/sqrt(r);
if(n.x < 0 || (n.x == 0 && n.y < 0 || (n.y == 0 && n.z < 0)))
r = -r; /* Canonicalize */
fe->n.x = n.x*r; fe->n.y = n.y*r; fe->n.z = n.z*r;
} while((fe = fe->next) != f->fedges);
}
/*
* Delete a face, erasing all edges.
*/
deface(f)
F *f;
{
register Fe *fe, *fee;
fe = f->fedges;
if(fe == NULL)
return;
do {
fee = fe->next;
unfedge(fe);
fe = fee;
} while(fe != f->fedges);
f->fedges = NULL;
tossedf++;
}
vcmp(p, q)
V **p, **q;
{
register V *vp, *vq;
register float d;
vp = *p;
vq = *q;
d = vp->p.x - vq->p.x;
if(d < 0) return -1;
if(d > 0) return 1;
d = vp->p.y - vq->p.y;
if(d < 0) return -1;
if(d > 0) return 1;
d = vp->p.z - vq->p.z;
if(d < 0) return -1;
if(d > 0) return 1;
if(vp->p.w == vq->p.w) return 0;
return(vp->p.w < vq->p.w ? -1 : 1);
}
vmerge()
{
V **vp;
int i, j;
register V *a, *b;
int nexti;
vp = NewN(V *, nv);
for(i = 0, a = Vs; i < nv; i++) {
a->ref = 0;
vp[i] = a++;
}
qsort(vp, nv, sizeof(V *), vcmp);
/*
* Now all matches will occur within X-runs
*/
for(i = 0; i < nv; i = nexti ? nexti : i+1) {
nexti = 0;
a = vp[i];
if(a->ref)
continue;
for(j = i; ++j < nv; ) {
b = vp[j];
if(b->ref)
continue;
if(b->p.x - a->p.x > EPS_POINT)
break;
if(fabs(a->p.y - b->p.y) < EPS_POINT
&& fabs(a->p.z - b->p.z) < EPS_POINT
&& (vdim == 3 || fabs(a->p.w - b->p.w) < EPS_POINT)) {
debug && printf("# Vtx %d->%d\n", b->index, a->index);
b->index = a->index;
b->ref++;
} else if(!nexti)
nexti = j;
}
}
free(vp);
}
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.