X-Git-Url: http://www.chiark.greenend.org.uk/ucgi/~ian/git?p=moebius2.git;a=blobdiff_plain;f=output.c;fp=output.c;h=afdc3df5d02bfef56bfdd87d2a730459cc0daa00;hp=0000000000000000000000000000000000000000;hb=6bcc34b96b6674e9ece1c339ada0557f7d2fd1b5;hpb=c620934d56d69abb34aeb475b09471f811d39898

diff --git a/output.c b/output.c
new file mode 100644
index 0000000..afdc3df
--- /dev/null
+++ b/output.c
@@ -0,0 +1,436 @@
+/*
+ * 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;
+}