* - scale the coordinates
* - translate the coordinates so they're all positive
*/
+/*
+ * ./output-64 <dense-64.cfm 0.1 50 >t.stl
+ * ./output-125 <best-125.cfm >t.stl 0.1 50
+ * meshlab t.stl
+ */
/*
* Re STL, see:
* http://www.ennex.com/~fabbers/StL.asp
* / *B2 2\ / /1 *A \ / *A \ / *A \ /
* \ / \/:/ 0 \ / \ / \ /
* *C ____________|*C ____________ *C ___________ *C ___________ *C __
- * /\ 3 | : 0 /\ /\ /\
+ * /\ 3 |/\ 0 /\ /\ /\
* / \ / :\ / \ / \ / \
* \ *A3 4/ \: \5 5*B / \ *B / \ *B / \
* \ / \ \ / \ / \ /
* *B \ / *B | \ / *A \ / *A \ / *A
* \ / : | \ / . ' \ / \ /
* _______ *C _____ :_|___ *C' __________ *C ___________ *C _________
- * /\ : | 3 /\. 0 /\ /\
- * / \ : | / \ ` . / \ / \
+ * /\ : | 3 /\ 0 /\ /\
+ * / \ : | / \ / \ / \
* :
* . . : . . . . . .
* :
- * \ / : | \ / . ' \ / \ /
- * _______ *C _____ :_|___ *C' __________ *C ___________ *C _________
- * /\ : | 3 /\. 0 /\ /\
- * / \ : | / \ ` . / \ / \
+ * \ / : | \ / \ / \ /
+ * _______ *C _____ :_|___ *C ___________ *C ___________ *C _________
+ * /\ : | 3 /\ 0 /\ /\
+ * / \ : | / \ / \ / \
* *A / \ *A | /4 5\ *B / \ *B / \ *B
* / \ 1: / / \ / \ / \
* / \ :/ / \ / \ / \
* G: Each F is the mean of those two adjacent Es with the
* same angle in their respective Ps.
*
- * The outtriangles are:
- *
- * For each non-rim vertex on each side, the six triangles formed by
- * its C and the surrounding A's and B's.
- *
- * For each rim vertex on each side, the two triangles formed by its
- * D and the nearest As and Bs (two As and one B or vice versa).
+ * The outfacets are:
*
- * For each rim edge on each side, the triangle formed by that edge's
- * ends' Ds and the corresponding A or B.
+ * For each input vertex on each side, the six (or perhaps only three)
+ * triangles formed by its C or D and the surrounding Cs and/or Ds.
*
* For each G, the six triangles formed by that G and the adjacent
* four Fs (or two Fs and two Ds) and two Es.
#define NG 10
#define NDEF (NG*2+1)
+#define OUTPUT_ARRAY_LIST(DO_OUTPUT_ARRAY) \
+ DO_OUTPUT_ARRAY(ovAB) \
+ DO_OUTPUT_ARRAY(ovC) \
+ DO_OUTPUT_ARRAY(ovDEF) \
+ DO_OUTPUT_ARRAY(ovG)
+
static OutVertex ovAB[N][2][2]; /* indices: vertex to W, A=0 B=1, side */
static OutVertex ovC[N][2];
static OutVertex ovDEF[X][2][NDEF]; /* x, !!y, angle */
static double thick; /* in input units */
static double scale;
+static void outfacet(int rev, const OutVertex *a,
+ const OutVertex *b, const OutVertex *c);
+
+typedef int int_map(int);
+static int defs_aroundmap_swap(int around) { return NDEF-1-around; }
+static int int_identity_function(int i) { return i; }
+
static void normalise_thick(double a[D3]) {
/* multiplies a by a scalar so that its magnitude is thick */
int k;
}
static OutVertex *invertex2outvertexab(int v0, int e, int side) {
- int vref, ab;
+ int vref, vchk, ab;
switch (e) {
- case 0: vref= v0 ; ab=0; break;
- case 1: vref=EDGE_END2(v0,2); ab=1; break;
- case 2: vref=EDGE_END2(v0,3); ab=0; break;
- case 3: vref=EDGE_END2(v0,3); ab=1; break;
- case 4: vref=EDGE_END2(v0,4); ab=0; break;
- case 5: vref= v0 ; ab=1; break;
+ case 0: vref=v0; vchk=EDGE_END2(v0,1); ab=0; break;
+ case 1: vref= vchk=EDGE_END2(v0,2); ab=1; break;
+ case 2: vref=EDGE_END2(v0,3); vchk=EDGE_END2(v0,2); ab=0; break;
+ case 3: vref=EDGE_END2(v0,3); vchk=EDGE_END2(v0,4); ab=1; break;
+ case 4: vref= vchk=EDGE_END2(v0,4); ab=0; break;
+ case 5: vref=v0; vchk=EDGE_END2(v0,5); ab=1; break;
default: abort();
}
- if (vref<0) return 0;
- int sw= VERTICES_SPAN_JOIN_P(v0,vref);
+ if (vchk<0) return 0;
+ int sw= vertices_span_join_p(v0,vref);
return &ovAB[vref][ab^sw][side^sw];
}
/*---------- output vertices ----------*/
+#define Ok(ov, value) ((ov).p[k]= outvertex_coord_check(value))
+
+static double outvertex_coord_check(double value) {
+ assert(-10 < value && value < 10);
+ return value;
+}
+
static void compute_outvertices(void) {
int v0,k,side,ab,x,y;
- FOR_VERTEX(v0) {
+ FOR_VERTEX(v0, INNER) {
for (ab=0; ab<2; ab++) {
int v1= EDGE_END2(v0, ab?5:0);
int v2= EDGE_END2(v0, ab?0:1);
- if (v1<0 || v2<0) {
- K ovAB[v0][ab][0].p[k]= ovAB[v0][ab][1].p[k]= NAN;
+ if (v1<0 || v2<0)
continue;
- }
double normal[D3], centroid[D3];
triangle_normal(normal, in[v0],in[v1],in[v2]);
normalise_thick(normal);
K centroid[k]= (in[v0][k] + in[v1][k] + in[v2][k]) / 3.0;
- K ovAB[v0][ab][0].p[k]= centroid[k] + normal[k];
- K ovAB[v0][ab][1].p[k]= centroid[k] - normal[k];
+ K Ok(ovAB[v0][ab][0], centroid[k] + normal[k]);
+ K Ok(ovAB[v0][ab][1], centroid[k] - normal[k]);
}
}
- FOR_VERTEX(v0) {
- double centroid[D3];
+ FOR_VERTEX(v0, INNER) {
int vw= EDGE_END2(v0,3);
int vnw= EDGE_END2(v0,2);
int vsw= EDGE_END2(v0,4);
- if (vnw<0 || vsw<0 || vw<0) {
- K ovC[v0][0].p[k]= ovC[v0][1].p[k]= NAN;
+ if (vnw<0 || vsw<0 || vw<0)
continue;
- }
FOR_SIDE {
+ double adjust[D3];
+ int e;
K adjust[k]= 0;
FOR_VPEDGE(e) {
OutVertex *ovab= invertex2outvertexab(v0,e,side);
- K adjust[k] += ovab->k[k];
+ K {
+ assert(!isnan(ovab->p[k]));
+ adjust[k] += ovab->p[k];
+ }
}
K adjust[k] /= 6;
K adjust[k] -= in[v0][k];
normalise_thick(adjust);
- K ovC[v0][side].p[k]= in[v0][k] + adjust[k];
+ K Ok(ovC[v0][side], in[v0][k] + adjust[k]);
}
}
- FOR_RIM_VERTEX(y,x,v0) {
- double rim[D3], inner[D3], radius_cos[D3];, radius_sin[D3];
+ FOR_RIM_VERTEX(y,x,v0, INNER) {
+ double rim[D3], inner[D3], radius_cos[D3], radius_sin[D3];
int vback, vfwd, around;
/* compute mean rim vector, which is just the vector between
vback= EDGE_END2(v0,3);
vfwd= EDGE_END2(v0,0);
assert(vback>=0 && vfwd>=0);
- K rim[k]= in[fwd][k] - in[vback][k];
+ K rim[k]= in[vfwd][k] - in[vback][k];
/* compute the inner centroid */
vback= EDGE_END2(v0,4);
normalise_thick(radius_cos);
normalise_thick(radius_sin);
+ int_map *around_map= y ? int_identity_function : defs_aroundmap_swap;
+
for (around=0; around<NDEF; around++) {
- double angle= around * PI / (NDEF-1);
- double result[D3];
- K ovDEF[x][!!y][around].p[k]=
- in[v0][k] +
- cos(angle) * radius_cos[k] +
- sin(angle) * radius_sin[k];
+ double angle= around_map(around) * M_PI / (NDEF-1);
+ K Ok(ovDEF[x][!!y][around],
+ in[v0][k] +
+ cos(angle) * radius_cos[k] +
+ sin(angle) * radius_sin[k]);
}
}
- FOR_RIM_VERTEX(y,x,v0) {
+ FOR_RIM_VERTEX(y,x,v0, INNER) {
int vfwd= EDGE_END2(v0,0);
assert(vfwd >= 0);
- for (around=0; around<NG; around++) {
- K ovG[x][!!y][around].p[k]=
- (ovDEF[ x ][!!y][around*2].p[k] +
- ovDEF[vfwd & XMASK][!!y][around*2].p[k]) / 2;
+ int aroung;
+ int_map *around_map= vertices_span_join_p(v0,vfwd)
+ ? defs_aroundmap_swap : int_identity_function;
+ for (aroung=0; aroung<NG; aroung++) {
+ K Ok(ovG[x][!!y][aroung],
+ 0.5 * (ovDEF[ x ][!! y ][aroung*2+1].p[k] +
+ ovDEF[vfwd & XMASK][!!(vfwd & ~XMASK)][around_map(aroung*2+1)].p[k]));
}
}
}
-/*---------- output triangles ----------*/
+/*---------- output facets ----------*/
-static void outtriangles_around(int reverse, OutVertex *middle,
- int nsurr, OutVertex *surround[nsurr]) {
+static void outfacets_around(int reverse, OutVertex *middle,
+ int nsurr, OutVertex *surround[nsurr]) {
/* Some entries in surround may be 0, in which case all affected
- * triangles will be skipped */
+ * facets will be skipped */
int i;
for (i=0; i<nsurr; i++) {
OutVertex *s0= surround[i];
OutVertex *s1= surround[(i+1) % nsurr];
if (!s0 || !s1) continue;
- outtriangle(reverse, middle,s0,s1);
+ outfacet(reverse, middle,s0,s1);
}
}
-static void outtriangles(void) {
- int v0,e0,e1;
+static OutVertex *invertex2outvertexcd(v0,side) {
+ if (!RIM_VERTEX_P(v0)) return &ovC[v0][side];
+
+ int around= side ? NDEF-1 : 0;
+ int rimy= !!(v0 & ~XMASK);
+ return &ovDEF[v0 & XMASK][rimy][around];
+}
+
+static void outfacets(void) {
+ int v0,e,side,aroung;
- FOR_VERTEX(v0) {
+ FOR_VERTEX(v0, INNER) {
OutVertex *defs=0, *defs1=0;
- int (*defs1aroundmap)(int)=0, rimy;
+ int rimy=-1;
+ int_map *defs1aroundmap= 0;
if (RIM_VERTEX_P(v0)) {
- Outvertex *gs;
+ OutVertex *gs;
rimy= !!(v0 & ~XMASK);
int v1= EDGE_END2(v0,0); assert(v1>=0);
gs= ovG [v0 & XMASK][rimy];
defs= ovDEF[v0 & XMASK][rimy];
- defs1= ovDEF[v1 & XMASK][rimy];
- defs1aroundmap= VERTICES_SPAN_JOIN_P(v0,v1)
+ defs1= ovDEF[v1 & XMASK][!!(v1 & ~XMASK)];
+ defs1aroundmap= vertices_span_join_p(v0,v1)
? defs_aroundmap_swap : int_identity_function;
for (aroung=0; aroung<NG; aroung++) {
surround[e ]= &defs1[defs1aroundmap(around +e)];
surround[e+3]= &defs [ around+2-e ];
}
- outtriangles_around(rimy, &gs[aroung], 6,surround);
+ outfacets_around(rimy, &gs[aroung], 6,surround);
}
}
-
+
FOR_SIDE {
- if (defs) {
- int around= side ? NDEF-1 : 0;
- cd= &defs[around];
- OutVertex *ab= ovAB[v0][!rimy][side];
- OutVertex *cd1= &defs1[defs1aroundmap(around)];
- outtriangle(side^rimy,cd,ab,cd1);
- } else {
- cd= &ovC[v0][side];
+ int ab;
+ for (ab=0; ab<2; ab++) {
+ int v1= EDGE_END2(v0, ab ? 5 : 0);
+ int v2= EDGE_END2(v0, ab ? 0 : 1);
+ if (v1<0 || v2<0) continue;
+ outfacet(side,
+ invertex2outvertexcd(v0,side),
+ invertex2outvertexcd(v1,side^vertices_span_join_p(v0,v1)),
+ invertex2outvertexcd(v2,side^vertices_span_join_p(v0,v2)));
}
- OutVertex *abs[6];
- FOR_VPEDGE(e)
- abs[e0]= invertex2outvertexab(v0,e,side);
- outtriangles_around(side, cd, 6,abs);
}
}
}
-/*---------- transformation (scale and perhaps shift) ----------*/
+/*---------- operations on very output vertex ----------*/
+
+#define DO_OUTPUT_ARRAY_OUTVERTEX_ARRAY(fn,ovX) \
+ ((fn)(sizeof((ovX))/sizeof(OutVertex), (OutVertex*)(ovX)))
-static void scaleshift_outvertex_array(int n, OutVertex ovX[n]) {
+static void blank_outvertex_array(int n, OutVertex ovX[n]) {
+ int i, k;
for (i=0; i<n; i++)
+ K ovX[i].p[k]= NAN;
+}
+
+static void blank_outvertices(void) {
+#define BLANK_OUTPUT_ARRAY(ovX) \
+ DO_OUTPUT_ARRAY_OUTVERTEX_ARRAY(blank_outvertex_array, ovX);
+ OUTPUT_ARRAY_LIST(BLANK_OUTPUT_ARRAY)
+}
+
+static void transform_outvertex_array(int n, OutVertex ovX[n]) {
+ int i, k;
+ for (i=0; i<n; i++) {
+ ovX[i].p[0] *= -1;
+ ovX[i].p[1] *= -1;
K ovX[i].p[k] *= scale;
+ }
/*
* double min[D3]= thick;
* if (ovX[i].p[k] < min)
*/
}
-#define SCALESHIFT_OUTVERTEX_ARRAY(ovX) \
- scaleshift_outvertex_array(sizeof((ovX))/sizeof(OutVertex),(OutVertex*)(ovX))
-
-static void scaleshift_outvertices(void) {
- SCALESHIFT_OUTVERTEX_ARRAY(ovAB);
- SCALESHIFT_OUTVERTEX_ARRAY(ovC);
- SCALESHIFT_OUTVERTEX_ARRAY(ovDEF);
- SCALESHIFT_OUTVERTEX_ARRAY(ovG);
+static void transform_outvertices(void) {
+#define TRANSFORM_OUTPUT_ARRAY(ovX) \
+ DO_OUTPUT_ARRAY_OUTVERTEX_ARRAY(transform_outvertex_array, ovX);
+ OUTPUT_ARRAY_LIST(TRANSFORM_OUTPUT_ARRAY)
}
/*---------- output file ----------*/
#define WR(x) wr((const void*)&(x), sizeof((x)))
static void wf(double d) {
-#if sizeof(float)==4
typedef float ieee754single;
-#endif
-#if defined(BIG_ENDIAN)
- union { Byte b[4]; ieee754single f; } value; value.f= d;
+ assert(sizeof(ieee754single)==4);
+ assert(!isnan(d));
+
+assert(d >= -1e3 && d <= 1e3);
+
+#if BYTE_ORDER==BIG_ENDIAN
+ union { uint8_t b[4]; ieee754single f; } value; value.f= d;
int i; for (i=3; i>=0; i--) WR(value.b[i]);
-#elif defined(LITTLE_ENDIAN)
+#elif BYTE_ORDER==LITTLE_ENDIAN
ieee754single f= d;
WR(f);
#else
#endif
}
-static uint32_t nouttriangles;
-static uint32_t nouttriangles_counted;
+static uint32_t noutfacets;
+static uint32_t noutfacets_counted;
+static long badfacets;
-static void outtriangle(int rev, OutVertex *a, OutVertex *b, OutVertex *c) {
- if (rev) { outtriangle(0, c,b,a); return; }
- nouttriangles++;
- if (!~nouttriangles_counted) return;
+static void outfacet(int rev, const OutVertex *a,
+ const OutVertex *b, const OutVertex *c) {
+ if (rev) { outfacet(0, c,b,a); return; }
- triangle_normal(normal, a.p, b.p, c.p);
+ double normal[D3];
+ int k;
+
+ K {
+ assert(!isnan(a->p[k]));
+ assert(!isnan(b->p[k]));
+ assert(!isnan(c->p[k]));
+ }
+
+ triangle_normal(normal, a->p, b->p, c->p);
double multby= 1/magnD(normal);
+
+ if (multby > 1e6) {
+ badfacets++;
+ return;
+ }
+
+ noutfacets++;
+ if (!~noutfacets_counted) return;
+
K normal[k] *= multby;
K wf(normal[k]);
- K wf(a.p[k]);
- K wf(b.p[k]);
- K wf(c.p[k]);
- uint16_t attrbc;
+ K wf(a->p[k]);
+ K wf(b->p[k]);
+ K wf(c->p[k]);
+ uint16_t attrbc=0;
WR(attrbc);
}
static void write_file(void) {
static const char header[80]= "#!/usr/bin/meshlab\n" "binary STL file\n";
- if (isatty(stdout)) die("will not write binary stl to tty!");
+ if (isatty(1)) fail("will not write binary stl to tty!\n");
WR(header);
- nouttriangles_counted=~(uint32_t)0;
- nouttriangles=0;
- outtriangles();
- WR(nouttriangles);
+ noutfacets_counted=~(uint32_t)0;
+ noutfacets=0;
+ outfacets();
+ WR(noutfacets);
- nouttriangles_counted= nouttriangles;
- nouttriangles=0;
- outtriangles();
- assert(nouttriangles == nouttriangles_counted);
+ noutfacets_counted= noutfacets;
+ noutfacets=0;
+ outfacets();
+ assert(noutfacets == noutfacets_counted);
if (fflush(stdout)) diee("fflush stdout");
+
+ if (badfacets) {
+ fprintf(stderr,"%ld degenerate facets!\n",badfacets);
+ exit(4);
+ }
}
errno= 0; r= fread(&in,sizeof(in),1,stdin);
if (r!=1) diee("fread");
+ mgraph_prepare();
+ blank_outvertices();
compute_outvertices();
- scaleshift_outvertices();
+ transform_outvertices();
write_file();
if (fclose(stdout)) diee("fclose stdout");
return 0;
~ian / moebius2.git / blobdiff