--- /dev/null
+/*
+ * Generates an STL file
+ * usage: ./output+<size> <thickness> <unit-circle-diameter>
+ * thickness is in original moebius unit circle units
+ * unit circle diameter is in whatever dimensions the STL file uses
+ *
+ * We have to:
+ * - add some thickness
+ * - define a shape for the rim so that it is solid
+ * - scale the coordinates
+ * - translate the coordinates so they're all positive
+ */
+/*
+ * Re STL, see:
+ * http://www.ennex.com/~fabbers/StL.asp
+ * http://en.wikipedia.org/wiki/STL_%28file_format%29
+ * saved here as stl-format-wikipedia.txt
+ */
+
+/*
+ * We add thickness by adding vertices as follows:
+ *
+ * :axis of symmetry
+ * sides | :
+ * swap| :
+ * | :
+ * *F | : *F *F *F *F
+ * | :
+ * *E *G | :*E *G *E *G *E *G *E
+ * | :
+ * *D ____________|*D ____________ *D ____________ *D ____________ *D __
+ * /\ 3 |/\ 0 /\ /\ /\
+ * / \ / :\ / \ / \ / \
+ * \ *A 4/ \: \5 *B / \ *B / \ *B / \
+ * \ / \ \ 1 / \ / \ /
+ * \ / :\ \ / \ / \ /
+ * *B \ / *B2| \2 1/ 0*A \ / *A \ / *A
+ * \ / : | \ / . ' \ / \ /
+ * _______ *C _____ :_|___ *C' __________ *C ___________ *C _________
+ * /\ : | 3 /\. 0 /\ /\
+ * / \ :3| / \ ` . / \ / \
+ * *A / \ *A | /4 5\ 5*B / \ *B / \ *B
+ * / \ 1: / / \ / \ / \
+ * / \ :/ / 4 \ / \ / \
+ * / *B2 2\ / /1 *A \ / *A \ / *A \ /
+ * \ / \/:/ 0 \ / \ / \ /
+ * *C ____________|*C ____________ *C ___________ *C ___________ *C __
+ * /\ 3 | : 0 /\ /\ /\
+ * / \ / :\ / \ / \ / \
+ * \ *A3 4/ \: \5 5*B / \ *B / \ *B / \
+ * \ / \ \ / \ / \ /
+ * \ / 4:\ \ / \ / \ /
+ * *B \ / *B | \ / *A \ / *A \ / *A
+ * \ / : | \ / . ' \ / \ /
+ * _______ *C _____ :_|___ *C' __________ *C ___________ *C _________
+ * /\ : | 3 /\. 0 /\ /\
+ * / \ : | / \ ` . / \ / \
+ * :
+ * . . : . . . . . .
+ * :
+
+ * \ / : | \ / . ' \ / \ /
+ * _______ *C _____ :_|___ *C' __________ *C ___________ *C _________
+ * /\ : | 3 /\. 0 /\ /\
+ * / \ : | / \ ` . / \ / \
+ * *A / \ *A | /4 5\ *B / \ *B / \ *B
+ * / \ 1: / / \ / \ / \
+ * / \ :/ / \ / \ / \
+ * / *B2 2\ / /1 0*A \ / *A \ / *A \ /
+ * \ / \/:/ \ / \ / \ /
+ * *D ____________|*D ____________ *D ___________ *D ___________ *D __
+ * 3 | : 0
+ * *E *G | :*E *G *E *G *E *G *E
+ * | :
+ * *F | : *F *F *F *F
+ * | :
+ *
+ * Each A,B,C,D represents two vertices - one on each side of the
+ * surface. Each E,F,G represents a several vertices in an arc around
+ * the rim. Digits are `e' values for edges or relative AB
+ * outvertices.
+ *
+ * The constructions we use are:
+ *
+ * A, B: Extend the normal vector of the containing intriangle
+ * from its centroid for thickness.
+ *
+ * C: Take mean position (centroid) of all surrounding
+ * computed A and B, and extend from base point in
+ * direction of that centroid for thickness.
+ *
+ * D, E, F:
+ * Compute notional rim vector R as mean of the two
+ * adjoining rim edges, and consider plane P as
+ * that passing through the base point normal to R.
+ * Project the centroid of the two adjoining non-rim
+ * vertices onto P. Now D(EF)+ED is the semicircle
+ * in P centred on the base point with radius thickness
+ * and which is opposite that centroid.
+ *
+ * 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).
+ *
+ * For each rim edge on each side, the triangle formed by that edge's
+ * ends' Ds and the corresponding A or B.
+ *
+ * For each G, the six triangles formed by that G and the adjacent
+ * four Fs (or two Fs and two Ds) and two Es.
+ */
+
+#include "mgraph.h"
+#include "common.h"
+
+/*---------- declarations and useful subroutines ----------*/
+
+#define FOR_SIDE for (side=0; side<2; side++)
+
+typedef struct {
+ double p[D3];
+} OutVertex;
+
+#define NG 10
+#define NDEF (NG*2+1)
+
+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 OutVertex ovG[X][2][NG]; /* x, !!y, angle; for G to the East of x */
+
+static Vertices in;
+
+static double thick; /* in input units */
+static double scale;
+
+static void normalise_thick(double a[D3]) {
+ /* multiplies a by a scalar so that its magnitude is thick */
+ int k;
+ double multby= thick / magnD(a);
+ K a[k] *= multby;
+}
+
+static void triangle_normal(double normal[D3], const double a[D3],
+ const double b[D3], const double c[D3]) {
+ double ab[D3], ac[D3];
+ int k;
+
+ K ab[k]= b[k] - a[k];
+ K ac[k]= c[k] - a[k];
+ xprod(normal,ab,ac);
+}
+
+static OutVertex *invertex2outvertexab(int v0, int e, int side) {
+ int vref, 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;
+ default: abort();
+ }
+ if (vref<0) return 0;
+ int sw= VERTICES_SPAN_JOIN_P(v0,vref);
+ return &ovAB[vref][ab^sw][side^sw];
+}
+
+/*---------- output vertices ----------*/
+
+static void compute_outvertices(void) {
+ int v0,k,side,ab,x,y;
+
+ FOR_VERTEX(v0) {
+ 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;
+ 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];
+ }
+ }
+ FOR_VERTEX(v0) {
+ double centroid[D3];
+ 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;
+ continue;
+ }
+ FOR_SIDE {
+ K adjust[k]= 0;
+ FOR_VPEDGE(e) {
+ OutVertex *ovab= invertex2outvertexab(v0,e,side);
+ K adjust[k] += ovab->k[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];
+ }
+ }
+ FOR_RIM_VERTEX(y,x,v0) {
+ 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
+ * the two adjacent rim vertex (ignoring the base vertex) */
+ vback= EDGE_END2(v0,3);
+ vfwd= EDGE_END2(v0,0);
+ assert(vback>=0 && vfwd>=0);
+ K rim[k]= in[fwd][k] - in[vback][k];
+
+ /* compute the inner centroid */
+ vback= EDGE_END2(v0,4);
+ if (vback>=0) { /* North rim */
+ vfwd= EDGE_END2(v0,5);
+ } else { /* South rim */
+ vback= EDGE_END2(v0,2);
+ vfwd= EDGE_END2(v0,1);
+ }
+ assert(vback>=0 && vfwd>=0);
+ K inner[k]= (in[vback][k] + in[vfwd][k]) / 2;
+ K inner[k] -= in[v0][k];
+
+ /* we compute the radius cos and sin vectors by cross producting
+ * the vector to the inner with the rim, and then again, and
+ * normalising. */
+ xprod(radius_cos,rim,inner);
+ xprod(radius_sin,rim,radius_cos);
+ normalise_thick(radius_cos);
+ normalise_thick(radius_sin);
+
+ 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];
+ }
+ }
+ FOR_RIM_VERTEX(y,x,v0) {
+ 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;
+ }
+ }
+}
+
+/*---------- output triangles ----------*/
+
+static void outtriangles_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 */
+ 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);
+ }
+}
+
+static void outtriangles(void) {
+ int v0,e0,e1;
+
+ FOR_VERTEX(v0) {
+ OutVertex *defs=0, *defs1=0;
+ int (*defs1aroundmap)(int)=0, rimy;
+ if (RIM_VERTEX_P(v0)) {
+ 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)
+ ? defs_aroundmap_swap : int_identity_function;
+
+ for (aroung=0; aroung<NG; aroung++) {
+ int around= aroung*2;
+ OutVertex *surround[6];
+ for (e=0; e<3; e++) {
+ surround[e ]= &defs1[defs1aroundmap(around +e)];
+ surround[e+3]= &defs [ around+2-e ];
+ }
+ outtriangles_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];
+ }
+ OutVertex *abs[6];
+ FOR_VPEDGE(e)
+ abs[e0]= invertex2outvertexab(v0,e,side);
+ outtriangles_around(side, cd, 6,abs);
+ }
+ }
+}
+
+/*---------- transformation (scale and perhaps shift) ----------*/
+
+static void scaleshift_outvertex_array(int n, OutVertex ovX[n]) {
+ for (i=0; i<n; i++)
+ K ovX[i].p[k] *= scale;
+ /*
+ * double min[D3]= thick;
+ * if (ovX[i].p[k] < min)
+ * min= ovX[i].p[k];
+ * for (i=0; i<n; i++) {
+ * K ovX[k].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);
+}
+
+/*---------- output file ----------*/
+
+static void wr(const void *p, size_t sz) {
+ if (fwrite(p,sz,1,stdout) != 1)
+ diee("fwrite");
+}
+
+#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;
+ int i; for (i=3; i>=0; i--) WR(value.b[i]);
+#elif defined(LITTLE_ENDIAN)
+ ieee754single f= d;
+ WR(f);
+#else
+# error not little or big endian!
+#endif
+}
+
+static uint32_t nouttriangles;
+static uint32_t nouttriangles_counted;
+
+static void outtriangle(int rev, OutVertex *a, OutVertex *b, OutVertex *c) {
+ if (rev) { outtriangle(0, c,b,a); return; }
+ nouttriangles++;
+ if (!~nouttriangles_counted) return;
+
+ triangle_normal(normal, a.p, b.p, c.p);
+ double multby= 1/magnD(normal);
+ 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;
+ 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!");
+
+ WR(header);
+
+ nouttriangles_counted=~(uint32_t)0;
+ nouttriangles=0;
+ outtriangles();
+ WR(nouttriangles);
+
+ nouttriangles_counted= nouttriangles;
+ nouttriangles=0;
+ outtriangles();
+ assert(nouttriangles == nouttriangles_counted);
+
+ if (fflush(stdout)) diee("fflush stdout");
+}
+
+
+/*---------- main program etc. ----------*/
+
+int main(int argc, const char *const *argv) {
+ int r;
+
+ if (argc!=3 || argv[1][0]=='-') { fputs("bad usage\n",stderr); exit(8); }
+ thick= atof(argv[1]);
+ scale= atof(argv[2]) * 0.5; /* circle is unit radius but arg is diameter */
+
+ errno= 0; r= fread(&in,sizeof(in),1,stdin);
+ if (r!=1) diee("fread");
+
+ compute_outvertices();
+ scaleshift_outvertices();
+ write_file();
+ if (fclose(stdout)) diee("fclose stdout");
+ return 0;
+}
~ian / moebius2.git / blobdiff