X-Git-Url: http://www.chiark.greenend.org.uk/ucgi/~ian/git?p=moebius2.git;a=blobdiff_plain;f=output.c;h=262991a41d978e0b02200a22ca6de92a3ed16830;hp=39a25356b9a9b352a03891c7f2337cb44b0ca3ec;hb=cf1b738d035a22d539a44ca50925fb23cdb7588c;hpb=4f9769b263b5e16c9c94ed3aba3829678a504261

diff --git a/output.c b/output.c
index 39a2535..262991a 100644
--- a/output.c
+++ b/output.c
@@ -10,6 +10,11 @@
  *   - 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
@@ -45,7 +50,7 @@
  *       /    *B2  2\  / /1   *A    \    /    *A    \    /    *A    \    /
  *   \  /            \/:/     0      \  /            \  /            \  /
  *    *C ____________|*C ____________ *C ___________  *C ___________  *C __
- *    /\          3  | :   0          /\              /\              /\
+ *    /\          3  |/\   0          /\              /\              /\
  *   /  \            / :\            /  \            /  \            /  \
  *       \    *A3  4/ \: \5  5*B    /    \    *B    /    \    *B    /    \
  *        \        /   \  \        /      \        /      \        /
@@ -53,16 +58,16 @@
  *    *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: / /      \        /      \        /      \
  *        /        \   :/ /        \      /        \      /        \
@@ -101,16 +106,10 @@
  *  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.
+ * The outfacets are:
  *
- *  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 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.
@@ -130,6 +129,12 @@ typedef struct {
 #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 */
@@ -140,8 +145,12 @@ static Vertices in;
 static double thick; /* in input units */
 static double scale;
 
-static void outtriangle(int rev, const OutVertex *a,
-			const OutVertex *b, const OutVertex *c);
+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 */
@@ -161,23 +170,30 @@ static void triangle_normal(double normal[D3], const double a[D3],
 }
 
 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;
 
@@ -185,38 +201,37 @@ static void compute_outvertices(void) {
     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) {
     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->p[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) {
@@ -250,58 +265,68 @@ static void compute_outvertices(void) {
     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 * M_PI / (NDEF-1);
-      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) {
     int vfwd= EDGE_END2(v0,0);
     assert(vfwd >= 0);
-    int around;
-    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 int defs_aroundmap_swap(int around) { return NDEF-around; }
-static int int_identity_function(int i) { return i; }
+static OutVertex *invertex2outvertexcd(v0,side) {
+  if (!RIM_VERTEX_P(v0)) return &ovC[v0][side];
 
-static void outtriangles(void) {
+  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) {
     OutVertex *defs=0, *defs1=0;
-    int (*defs1aroundmap)(int)=0, rimy=-1;
+    int rimy=-1;
+    int_map *defs1aroundmap= 0;
     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)
+      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++) {
@@ -311,35 +336,49 @@ static void outtriangles(void) {
 	  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 {
-      OutVertex *cd;
-      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 *ab[6];
-      FOR_VPEDGE(e)
-	ab[e]= invertex2outvertexab(v0,e,side);
-      outtriangles_around(side, cd, 6,ab);
     }
   }
 }     
 
-/*---------- 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)
@@ -349,14 +388,10 @@ static void scaleshift_outvertex_array(int n, OutVertex ovX[n]) {
    */
 }
 
-#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 ----------*/
@@ -372,11 +407,14 @@ static void wf(double d) {
   typedef float ieee754single;
 
   assert(sizeof(ieee754single)==4);
+  assert(!isnan(d));
 
-#if defined(BIG_ENDIAN)
+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
@@ -384,20 +422,34 @@ static void wf(double d) {
 #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, const OutVertex *a,
-			const OutVertex *b, const 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; }
 
   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]);
@@ -411,21 +463,26 @@ static void outtriangle(int rev, const OutVertex *a,
 static void write_file(void) {
   static const char header[80]= "#!/usr/bin/meshlab\n" "binary STL file\n";
 
-  if (isatty(1)) fail("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);
+  }
 }  
   
 
@@ -441,8 +498,9 @@ int main(int argc, const char *const *argv) {
   errno= 0; r= fread(&in,sizeof(in),1,stdin);
   if (r!=1) diee("fread");
 
+  blank_outvertices();
   compute_outvertices();
-  scaleshift_outvertices();
+  transform_outvertices();
   write_file();
   if (fclose(stdout)) diee("fclose stdout");
   return 0;