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