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