output.c

   1 /*
   2  * Generates an STL file
   3  *  usage: ./output+<size> <thickness> <unit-circle-diameter>
   4  *    thickness is in original moebius unit circle units
   5  *    unit circle diameter is in whatever dimensions the STL file uses
   6  *
   7  * We have to:
   8  *   - add some thickness
   9  *   - define a shape for the rim so that it is solid
  10  *   - scale the coordinates
  11  *   - translate the coordinates so they're all positive
  12  */
  13 /*
  14  *  ./output-64 <dense-64.cfm 0.1 50 >t.stl
  15  *  meshlab t.stl
  16  */
  17 /*
  18  * Re STL, see:
  19  *  http://www.ennex.com/~fabbers/StL.asp
  20  *  http://en.wikipedia.org/wiki/STL_%28file_format%29
  21  *    saved here as stl-format-wikipedia.txt
  22  */
  23
  24 /*
  25  * We add thickness by adding vertices as follows:
  26  *
  27  *                     :axis of symmetry
  28  *             sides | :
  29  *               swap| :
  30  *                   | :
  31  *       *F          | : *F              *F              *F              *F
  32  *                   | :
  33  *      *E      *G   | :*E      *G      *E      *G      *E      *G      *E
  34  *                   | :
  35  *    *D ____________|*D ____________ *D ____________ *D ____________ *D __
  36  *    /\          3  |/\   0          /\              /\              /\
  37  *   /  \            / :\            /  \            /  \            /  \
  38  *       \    *A   4/ \: \5   *B    /    \    *B    /    \    *B    /    \
  39  *        \        /   \  \    1   /      \        /      \        /
  40  *         \      /    :\  \      /        \      /        \      /
  41  *    *B    \    /    *B2|  \2  1/   0*A    \    /    *A    \    /    *A
  42  *           \  /      : |   \  / . '        \  /            \  /
  43  *   _______  *C _____ :_|___ *C' __________  *C ___________  *C _________
  44  *            /\       : | 3  /\.  0          /\              /\
  45  *           /  \      :3|   /  \ ` .        /  \            /  \
  46  *    *A    /    \    *A |  /4  5\   5*B    /    \    *B    /    \    *B
  47  *         /      \   1: / /      \        /      \        /      \
  48  *        /        \   :/ /    4   \      /        \      /        \
  49  *       /    *B2  2\  / /1   *A    \    /    *A    \    /    *A    \    /
  50  *   \  /            \/:/     0      \  /            \  /            \  /
  51  *    *C ____________|*C ____________ *C ___________  *C ___________  *C __
  52  *    /\          3  |/\   0          /\              /\              /\
  53  *   /  \            / :\            /  \            /  \            /  \
  54  *       \    *A3  4/ \: \5  5*B    /    \    *B    /    \    *B    /    \
  55  *        \        /   \  \        /      \        /      \        /
  56  *         \      /   4:\  \      /        \      /        \      /
  57  *    *B    \    /    *B |  \    /    *A    \    /    *A    \    /    *A
  58  *           \  /      : |   \  / . '        \  /            \  /
  59  *   _______  *C _____ :_|___ *C' __________  *C ___________  *C _________
  60  *            /\       : | 3  /\   0          /\              /\
  61  *           /  \      : |   /  \            /  \            /  \
  62  *                     :
  63  *         .      .    :   .      .        .      .        .      .
  64  *                     :
  65
  66  *           \  /      : |   \  /            \  /            \  /
  67  *   _______  *C _____ :_|___ *C ___________  *C ___________  *C _________
  68  *            /\       : | 3  /\   0          /\              /\
  69  *           /  \      : |   /  \            /  \            /  \
  70  *    *A    /    \    *A |  /4  5\    *B    /    \    *B    /    \    *B
  71  *         /      \   1: / /      \        /      \        /      \
  72  *        /        \   :/ /        \      /        \      /        \
  73  *       /    *B2  2\  / /1  0*A    \    /    *A    \    /    *A    \    /
  74  *   \  /            \/:/            \  /            \  /            \  /
  75  *    *D ____________|*D ____________ *D ___________  *D ___________  *D __
  76  *                 3 | :   0
  77  *      *E    *G     | :*E    *G        *E    *G        *E    *G        *E
  78  *                   | :
  79  *       *F          | : *F              *F              *F              *F
  80  *                   | :
  81  *
  82  * Each A,B,C,D represents two vertices - one on each side of the
  83  * surface.  Each E,F,G represents a several vertices in an arc around
  84  * the rim.  Digits are `e' values for edges or relative AB
  85  * outvertices.
  86  *
  87  * The constructions we use are:
  88  *
  89  *  A, B:  Extend the normal vector of the containing intriangle
  90  *         from its centroid for thickness.
  91  *
  92  *  C:     Take mean position (centroid) of all surrounding
  93  *         computed A and B, and extend from base point in
  94  *         direction of that centroid for thickness.
  95  *
  96  *  D, E, F:
  97  *         Compute notional rim vector R as mean of the two
  98  *         adjoining rim edges, and consider plane P as
  99  *         that passing through the base point normal to R.
 100  *         Project the centroid of the two adjoining non-rim
 101  *         vertices onto P.  Now D(EF)+ED is the semicircle
 102  *         in P centred on the base point with radius thickness
 103  *         and which is opposite that centroid.
 104  *
 105  *  G:     Each F is the mean of those two adjacent Es with the
 106  *         same angle in their respective Ps.
 107  *
 108  * The outfacets are:
 109  *
 110  *  For each non-rim vertex on each side, the six triangles formed by
 111  *  its C and the surrounding A's and B's.
 112  *
 113  *  For each rim vertex on each side, the two triangles formed by its
 114  *  D and the nearest As and Bs (two As and one B or vice versa).
 115  *
 116  *  For each rim edge on each side, the triangle formed by that edge's
 117  *  ends' Ds and the corresponding A or B.
 118  *
 119  *  For each G, the six triangles formed by that G and the adjacent
 120  *  four Fs (or two Fs and two Ds) and two Es.
 121  */
 122
 123 #include "mgraph.h"
 124 #include "common.h"
 125
 126 /*---------- declarations and useful subroutines ----------*/
 127
 128 #define FOR_SIDE for (side=0; side<2; side++)
 129
 130 typedef struct {
 131   double p[D3];
 132 } OutVertex;
 133
 134 #define NG   10
 135 #define NDEF (NG*2+1)
 136
 137 #define OUTPUT_ARRAY_LIST(DO_OUTPUT_ARRAY)      \
 138   DO_OUTPUT_ARRAY(ovAB)                         \
 139   DO_OUTPUT_ARRAY(ovC)                          \
 140   DO_OUTPUT_ARRAY(ovDEF)                        \
 141   DO_OUTPUT_ARRAY(ovG)
 142
 143 static OutVertex ovAB[N][2][2]; /* indices: vertex to W, A=0 B=1, side */
 144 static OutVertex ovC[N][2];
 145 static OutVertex ovDEF[X][2][NDEF]; /* x, !!y, angle */
 146 static OutVertex ovG[X][2][NG]; /* x, !!y, angle; for G to the East of x */
 147
 148 static Vertices in;
 149
 150 static double thick; /* in input units */
 151 static double scale;
 152
 153 static void outfacet(int rev, const OutVertex *a,
 154                      const OutVertex *b, const OutVertex *c);
 155
 156 typedef int int_map(int);
 157 static int defs_aroundmap_swap(int around) { return NDEF-around; }
 158 static int int_identity_function(int i) { return i; }
 159
 160 static void normalise_thick(double a[D3]) {
 161   /* multiplies a by a scalar so that its magnitude is thick */
 162   int k;
 163   double multby= thick / magnD(a);
 164   K a[k] *= multby;
 165 }
 166
 167 static void triangle_normal(double normal[D3], const double a[D3],
 168                             const double b[D3], const double c[D3]) {
 169   double ab[D3], ac[D3];
 170   int k;
 171
 172   K ab[k]= b[k] - a[k];
 173   K ac[k]= c[k] - a[k];
 174   xprod(normal,ab,ac);
 175 }
 176
 177 static OutVertex *invertex2outvertexab(int v0, int e, int side) {
 178   int vref, vchk, ab;
 179   switch (e) {
 180   case 0:  vref=v0;               vchk=EDGE_END2(v0,1);  ab=0;   break;
 181   case 1:  vref=                  vchk=EDGE_END2(v0,2);  ab=1;   break;
 182   case 2:  vref=EDGE_END2(v0,3);  vchk=EDGE_END2(v0,2);  ab=0;   break;
 183   case 3:  vref=EDGE_END2(v0,3);  vchk=EDGE_END2(v0,4);  ab=1;   break;
 184   case 4:  vref=                  vchk=EDGE_END2(v0,4);  ab=0;   break;
 185   case 5:  vref=v0;               vchk=EDGE_END2(v0,5);  ab=1;   break;
 186   default: abort();
 187   }
 188   if (vchk<0) return 0;
 189   int sw= vertices_span_join_p(v0,vref);
 190   return &ovAB[vref][ab^sw][side^sw];
 191 }
 192
 193 /*---------- output vertices ----------*/
 194
 195 #define Ok(ov, value) ((ov).p[k]= outvertex_coord_check(value))
 196
 197 static double outvertex_coord_check(double value) {
 198   assert(-10 < value && value < 10);
 199   return value;
 200 }
 201
 202 static void compute_outvertices(void) {
 203   int v0,k,side,ab,x,y;
 204
 205   FOR_VERTEX(v0) {
 206     for (ab=0; ab<2; ab++) {
 207       int v1= EDGE_END2(v0, ab?5:0);
 208       int v2= EDGE_END2(v0, ab?0:1);
 209       if (v1<0 || v2<0)
 210         continue;
 211       double normal[D3], centroid[D3];
 212       triangle_normal(normal, in[v0],in[v1],in[v2]);
 213       normalise_thick(normal);
 214       K centroid[k]= (in[v0][k] + in[v1][k] + in[v2][k]) / 3.0;
 215       K Ok(ovAB[v0][ab][0], centroid[k] + normal[k]);
 216       K Ok(ovAB[v0][ab][1], centroid[k] - normal[k]);
 217     }
 218   }
 219   FOR_VERTEX(v0) {
 220     int vw= EDGE_END2(v0,3);
 221     int vnw= EDGE_END2(v0,2);
 222     int vsw= EDGE_END2(v0,4);
 223     if (vnw<0 || vsw<0 || vw<0)
 224       continue;
 225     FOR_SIDE {
 226       double adjust[D3];
 227       int e;
 228       K adjust[k]= 0;
 229       FOR_VPEDGE(e) {
 230         OutVertex *ovab= invertex2outvertexab(v0,e,side);
 231         K {
 232           assert(!isnan(ovab->p[k]));
 233           adjust[k] += ovab->p[k];
 234         }
 235       }
 236       K adjust[k] /= 6;
 237       K adjust[k] -= in[v0][k];
 238       normalise_thick(adjust);
 239       K Ok(ovC[v0][side], in[v0][k] + adjust[k]);
 240     }
 241   }
 242   FOR_RIM_VERTEX(y,x,v0) {
 243     double rim[D3], inner[D3], radius_cos[D3], radius_sin[D3];
 244     int vback, vfwd, around;
 245
 246     /* compute mean rim vector, which is just the vector between
 247      * the two adjacent rim vertex (ignoring the base vertex) */
 248     vback= EDGE_END2(v0,3);
 249     vfwd=  EDGE_END2(v0,0);
 250     assert(vback>=0 && vfwd>=0);
 251     K rim[k]= in[vfwd][k] - in[vback][k];
 252
 253     /* compute the inner centroid */
 254     vback= EDGE_END2(v0,4);
 255     if (vback>=0) { /* North rim */
 256       vfwd= EDGE_END2(v0,5);
 257     } else { /* South rim */
 258       vback= EDGE_END2(v0,2);
 259       vfwd=  EDGE_END2(v0,1);
 260     }
 261     assert(vback>=0 && vfwd>=0);
 262     K inner[k]= (in[vback][k] + in[vfwd][k]) / 2;
 263     K inner[k] -= in[v0][k];
 264
 265     /* we compute the radius cos and sin vectors by cross producting
 266      * the vector to the inner with the rim, and then again, and
 267      * normalising. */
 268     xprod(radius_cos,rim,inner);
 269     xprod(radius_sin,rim,radius_cos);
 270     normalise_thick(radius_cos);
 271     normalise_thick(radius_sin);
 272
 273     int_map *around_map= y ? int_identity_function : defs_aroundmap_swap;
 274
 275     for (around=0; around<NDEF; around++) {
 276       double angle= around_map(around) * M_PI / (NDEF-1);
 277       K Ok(ovDEF[x][!!y][around],
 278            in[v0][k] +
 279            cos(angle) * radius_cos[k] +
 280            sin(angle) * radius_sin[k]);
 281     }
 282   }
 283   FOR_RIM_VERTEX(y,x,v0) {
 284     int vfwd= EDGE_END2(v0,0);
 285     assert(vfwd >= 0);
 286     int around;
 287     for (around=0; around<NG; around++) {
 288       K Ok(ovG[x][!!y][around],
 289            0.5 * (ovDEF[ x          ][!!y][around*2].p[k] +
 290                   ovDEF[vfwd & XMASK][!!y][around*2].p[k]));
 291     }
 292   }
 293 }
 294
 295 /*---------- output facets ----------*/
 296
 297 static void outfacets_around(int reverse, OutVertex *middle,
 298                              int nsurr, OutVertex *surround[nsurr]) {
 299   /* Some entries in surround may be 0, in which case all affected
 300    *  facets will be skipped */
 301   int i;
 302   for (i=0; i<nsurr; i++) {
 303     OutVertex *s0= surround[i];
 304     OutVertex *s1= surround[(i+1) % nsurr];
 305     if (!s0 || !s1) continue;
 306     outfacet(reverse, middle,s0,s1);
 307   }
 308 }
 309
 310 static void outfacets(void) {
 311   int v0,e,side,aroung, k;
 312
 313   FOR_VERTEX(v0) {
 314     OutVertex *defs=0, *defs1=0;
 315     int rimy=-1;
 316     int_map *defs1aroundmap= 0;
 317     if (RIM_VERTEX_P(v0)) {
 318       OutVertex *gs;
 319       rimy= !!(v0 & ~XMASK);
 320       int v1= EDGE_END2(v0,0);  assert(v1>=0);
 321       gs=    ovG  [v0 & XMASK][rimy];
 322       defs=  ovDEF[v0 & XMASK][rimy];
 323       defs1= ovDEF[v1 & XMASK][rimy];
 324       defs1aroundmap= vertices_span_join_p(v0,v1)
 325         ? defs_aroundmap_swap : int_identity_function;
 326
 327       for (aroung=0; aroung<NG; aroung++) {
 328         int around= aroung*2;
 329         OutVertex *surround[6];
 330         for (e=0; e<3; e++) {
 331           surround[e  ]= &defs1[defs1aroundmap(around  +e)];
 332           surround[e+3]= &defs [               around+2-e ];
 333         }
 334         outfacets_around(rimy, &gs[aroung], 6,surround);
 335       }
 336     }
 337
 338     FOR_SIDE {
 339       OutVertex *cd;
 340       if (defs) {
 341         int around= side ? NDEF-1 : 0;
 342         cd= &defs[around];
 343         OutVertex *ab= &ovAB[v0][!rimy][side];
 344         OutVertex *cd1= &defs1[defs1aroundmap(around)];
 345         outfacet(side^rimy,cd,ab,cd1);
 346       } else {
 347         cd= &ovC[v0][side];
 348       }
 349       OutVertex *ab[6];
 350       FOR_VPEDGE(e) {
 351         ab[e]= invertex2outvertexab(v0,e,side);
 352         if (ab[e])
 353           K assert(!isnan(ab[e]->p[k]));
 354       }
 355       outfacets_around(side, cd, 6,ab);
 356     }
 357   }
 358 }
 359
 360 /*---------- operations on very output vertex ----------*/
 361
 362 #define DO_OUTPUT_ARRAY_OUTVERTEX_ARRAY(fn,ovX) \
 363   ((fn)(sizeof((ovX))/sizeof(OutVertex), (OutVertex*)(ovX)))
 364
 365 static void blank_outvertex_array(int n, OutVertex ovX[n]) {
 366   int i, k;
 367   for (i=0; i<n; i++)
 368     K ovX[i].p[k]= NAN;
 369 }
 370
 371 static void blank_outvertices(void) {
 372 #define BLANK_OUTPUT_ARRAY(ovX) \
 373       DO_OUTPUT_ARRAY_OUTVERTEX_ARRAY(blank_outvertex_array, ovX);
 374   OUTPUT_ARRAY_LIST(BLANK_OUTPUT_ARRAY)
 375 }
 376
 377 static void transform_outvertex_array(int n, OutVertex ovX[n]) {
 378   int i, k;
 379   for (i=0; i<n; i++)
 380     K ovX[i].p[k] *= scale;
 381   /*
 382    *   double min[D3]= thick;
 383    *            if (ovX[i].p[k] < min)
 384    *        min= ovX[i].p[k];
 385    *      for (i=0; i<n; i++) {
 386    *  K ovX[k].p[k] -= min;
 387    */
 388 }
 389
 390 static void transform_outvertices(void) {
 391 #define TRANSFORM_OUTPUT_ARRAY(ovX) \
 392       DO_OUTPUT_ARRAY_OUTVERTEX_ARRAY(transform_outvertex_array, ovX);
 393   OUTPUT_ARRAY_LIST(TRANSFORM_OUTPUT_ARRAY)
 394 }
 395
 396 /*---------- output file ----------*/
 397
 398 static void wr(const void *p, size_t sz) {
 399   if (fwrite(p,sz,1,stdout) != 1)
 400     diee("fwrite");
 401 }
 402
 403 #define WR(x) wr((const void*)&(x), sizeof((x)))
 404
 405 static void wf(double d) {
 406   typedef float ieee754single;
 407
 408   assert(sizeof(ieee754single)==4);
 409   assert(!isnan(d));
 410
 411 assert(d >= -1e3 && d <= 1e3);
 412
 413 #if BYTE_ORDER==BIG_ENDIAN
 414   union { uint8_t b[4]; ieee754single f; } value;  value.f= d;
 415   int i;  for (i=3; i>=0; i--) WR(value.b[i]);
 416 #elif BYTE_ORDER==LITTLE_ENDIAN
 417   ieee754single f= d;
 418   WR(f);
 419 #else
 420 # error not little or big endian!
 421 #endif
 422 }
 423
 424 static uint32_t noutfacets;
 425 static uint32_t noutfacets_counted;
 426
 427 static void outfacet(int rev, const OutVertex *a,
 428                      const OutVertex *b, const OutVertex *c) {
 429   if (rev) { outfacet(0, c,b,a); return; }
 430
 431   double normal[D3];
 432   int k;
 433
 434   K {
 435     assert(!isnan(a->p[k]));
 436     assert(!isnan(b->p[k]));
 437     assert(!isnan(c->p[k]));
 438   }
 439
 440   triangle_normal(normal, a->p, b->p, c->p);
 441   double multby= 1/magnD(normal);
 442
 443   if (multby > 1e6)
 444     return;
 445
 446   noutfacets++;
 447   if (!~noutfacets_counted) return;
 448
 449   K normal[k] *= multby;
 450
 451   K wf(normal[k]);
 452   K wf(a->p[k]);
 453   K wf(b->p[k]);
 454   K wf(c->p[k]);
 455   uint16_t attrbc=0;
 456   WR(attrbc);
 457 }
 458
 459 static void write_file(void) {
 460   static const char header[80]= "#!/usr/bin/meshlab\n" "binary STL file\n";
 461
 462   if (isatty(1)) fail("will not write binary stl to tty!\n");
 463
 464   WR(header);
 465
 466   noutfacets_counted=~(uint32_t)0;
 467   noutfacets=0;
 468   outfacets();
 469   WR(noutfacets);
 470
 471   noutfacets_counted= noutfacets;
 472   noutfacets=0;
 473   outfacets();
 474   assert(noutfacets == noutfacets_counted);
 475
 476   if (fflush(stdout)) diee("fflush stdout");
 477 }
 478
 479
 480 /*---------- main program etc. ----------*/
 481
 482 int main(int argc, const char *const *argv) {
 483   int r;
 484
 485   if (argc!=3 || argv[1][0]=='-') { fputs("bad usage\n",stderr); exit(8); }
 486   thick= atof(argv[1]);
 487   scale= atof(argv[2]) * 0.5; /* circle is unit radius but arg is diameter */
 488
 489   errno= 0; r= fread(&in,sizeof(in),1,stdin);
 490   if (r!=1) diee("fread");
 491
 492   blank_outvertices();
 493   compute_outvertices();
 494   transform_outvertices();
 495   write_file();
 496   if (fclose(stdout)) diee("fclose stdout");
 497   return 0;
 498 }