chiark / gitweb /
compiles and most output facets are right; needs debugging
[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  * 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 outfacets 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 #define OUTPUT_ARRAY_LIST(DO_OUTPUT_ARRAY)      \
134   DO_OUTPUT_ARRAY(ovAB)                         \
135   DO_OUTPUT_ARRAY(ovC)                          \
136   DO_OUTPUT_ARRAY(ovDEF)                        \
137   DO_OUTPUT_ARRAY(ovG)
138
139 static OutVertex ovAB[N][2][2]; /* indices: vertex to W, A=0 B=1, side */
140 static OutVertex ovC[N][2];
141 static OutVertex ovDEF[X][2][NDEF]; /* x, !!y, angle */
142 static OutVertex ovG[X][2][NG]; /* x, !!y, angle; for G to the East of x */
143
144 static Vertices in;
145
146 static double thick; /* in input units */
147 static double scale;
148
149 static void outfacet(int rev, const OutVertex *a,
150                      const OutVertex *b, const OutVertex *c);
151
152 static void normalise_thick(double a[D3]) {
153   /* multiplies a by a scalar so that its magnitude is thick */
154   int k;
155   double multby= thick / magnD(a);
156   K a[k] *= multby;
157 }
158
159 static void triangle_normal(double normal[D3], const double a[D3],
160                             const double b[D3], const double c[D3]) {
161   double ab[D3], ac[D3];
162   int k;
163   
164   K ab[k]= b[k] - a[k];
165   K ac[k]= c[k] - a[k];
166   xprod(normal,ab,ac);
167 }
168
169 static OutVertex *invertex2outvertexab(int v0, int e, int side) {
170   int vref, vchk, ab;
171   switch (e) {
172   case 0:  vref=v0;               vchk=EDGE_END2(v0,1);  ab=0;   break;
173   case 1:  vref=                  vchk=EDGE_END2(v0,2);  ab=1;   break;
174   case 2:  vref=EDGE_END2(v0,3);  vchk=EDGE_END2(v0,2);  ab=0;   break;
175   case 3:  vref=EDGE_END2(v0,3);  vchk=EDGE_END2(v0,4);  ab=1;   break;
176   case 4:  vref=                  vchk=EDGE_END2(v0,4);  ab=0;   break;
177   case 5:  vref=v0;               vchk=EDGE_END2(v0,5);  ab=1;   break;
178   default: abort();
179   }
180   if (vchk<0) return 0;
181   int sw= vertices_span_join_p(v0,vref);
182   return &ovAB[vref][ab^sw][side^sw];
183 }
184
185 /*---------- output vertices ----------*/
186
187 #define Ok(ov, value) ((ov).p[k]= outvertex_coord_check(value))
188
189 static double outvertex_coord_check(double value) {
190   assert(-10 < value && value < 10);
191   return value;
192 }
193
194 static void compute_outvertices(void) {
195   int v0,k,side,ab,x,y;
196
197   FOR_VERTEX(v0) {
198     for (ab=0; ab<2; ab++) {
199       int v1= EDGE_END2(v0, ab?5:0);
200       int v2= EDGE_END2(v0, ab?0:1);
201       if (v1<0 || v2<0)
202         continue;
203       double normal[D3], centroid[D3];
204       triangle_normal(normal, in[v0],in[v1],in[v2]);
205       normalise_thick(normal);
206       K centroid[k]= (in[v0][k] + in[v1][k] + in[v2][k]) / 3.0;
207       K Ok(ovAB[v0][ab][0], centroid[k] + normal[k]);
208       K Ok(ovAB[v0][ab][1], centroid[k] - normal[k]);
209     }
210   }
211   FOR_VERTEX(v0) {
212     int vw= EDGE_END2(v0,3);
213     int vnw= EDGE_END2(v0,2);
214     int vsw= EDGE_END2(v0,4);
215     if (vnw<0 || vsw<0 || vw<0)
216       continue;
217     FOR_SIDE {
218       double adjust[D3];
219       int e;
220       K adjust[k]= 0;
221       FOR_VPEDGE(e) {
222         OutVertex *ovab= invertex2outvertexab(v0,e,side);
223         K {
224           assert(!isnan(ovab->p[k]));
225           adjust[k] += ovab->p[k];
226         }
227       }
228       K adjust[k] /= 6;
229       K adjust[k] -= in[v0][k];
230       normalise_thick(adjust);
231       K Ok(ovC[v0][side], in[v0][k] + adjust[k]);
232     }
233   }
234   FOR_RIM_VERTEX(y,x,v0) {
235     double rim[D3], inner[D3], radius_cos[D3], radius_sin[D3];
236     int vback, vfwd, around;
237
238     /* compute mean rim vector, which is just the vector between
239      * the two adjacent rim vertex (ignoring the base vertex) */
240     vback= EDGE_END2(v0,3);
241     vfwd=  EDGE_END2(v0,0);
242     assert(vback>=0 && vfwd>=0);
243     K rim[k]= in[vfwd][k] - in[vback][k];
244
245     /* compute the inner centroid */
246     vback= EDGE_END2(v0,4);
247     if (vback>=0) { /* North rim */
248       vfwd= EDGE_END2(v0,5);
249     } else { /* South rim */
250       vback= EDGE_END2(v0,2);
251       vfwd=  EDGE_END2(v0,1);
252     }
253     assert(vback>=0 && vfwd>=0);
254     K inner[k]= (in[vback][k] + in[vfwd][k]) / 2;
255     K inner[k] -= in[v0][k];
256
257     /* we compute the radius cos and sin vectors by cross producting
258      * the vector to the inner with the rim, and then again, and
259      * normalising. */
260     xprod(radius_cos,rim,inner);
261     xprod(radius_sin,rim,radius_cos);
262     normalise_thick(radius_cos);
263     normalise_thick(radius_sin);
264
265     for (around=0; around<NDEF; around++) {
266       double angle= around * M_PI / (NDEF-1);
267       K Ok(ovDEF[x][!!y][around],
268            in[v0][k] +
269            cos(angle) * radius_cos[k] +
270            sin(angle) * radius_sin[k]);
271     }
272   }
273   FOR_RIM_VERTEX(y,x,v0) {
274     int vfwd= EDGE_END2(v0,0);
275     assert(vfwd >= 0);
276     int around;
277     for (around=0; around<NG; around++) {
278       K Ok(ovG[x][!!y][around],
279            0.5 * (ovDEF[ x          ][!!y][around*2].p[k] +
280                   ovDEF[vfwd & XMASK][!!y][around*2].p[k]));
281     }
282   }
283 }
284
285 /*---------- output facets ----------*/
286
287 static void outfacets_around(int reverse, OutVertex *middle,
288                              int nsurr, OutVertex *surround[nsurr]) {
289   /* Some entries in surround may be 0, in which case all affected
290    *  facets will be skipped */
291   int i;
292   for (i=0; i<nsurr; i++) {
293     OutVertex *s0= surround[i];
294     OutVertex *s1= surround[(i+1) % nsurr];
295     if (!s0 || !s1) continue;
296     outfacet(reverse, middle,s0,s1);
297   }
298 }
299
300 static int defs_aroundmap_swap(int around) { return NDEF-around; }
301 static int int_identity_function(int i) { return i; }
302
303 static void outfacets(void) {
304   int v0,e,side,aroung, k;
305   
306   FOR_VERTEX(v0) {
307     OutVertex *defs=0, *defs1=0;
308     int (*defs1aroundmap)(int)=0, rimy=-1;
309     if (RIM_VERTEX_P(v0)) {
310       OutVertex *gs;
311       rimy= !!(v0 & ~XMASK);
312       int v1= EDGE_END2(v0,0);  assert(v1>=0);
313       gs=    ovG  [v0 & XMASK][rimy];
314       defs=  ovDEF[v0 & XMASK][rimy];
315       defs1= ovDEF[v1 & XMASK][rimy];
316       defs1aroundmap= vertices_span_join_p(v0,v1)
317         ? defs_aroundmap_swap : int_identity_function;
318
319       for (aroung=0; aroung<NG; aroung++) {
320         int around= aroung*2;
321         OutVertex *surround[6];
322         for (e=0; e<3; e++) {
323           surround[e  ]= &defs1[defs1aroundmap(around  +e)];
324           surround[e+3]= &defs [               around+2-e ];
325         }
326         outfacets_around(rimy, &gs[aroung], 6,surround);
327       }
328     }
329
330     FOR_SIDE {
331       OutVertex *cd;
332       if (defs) {
333         int around= side ? NDEF-1 : 0;
334         cd= &defs[around];
335         OutVertex *ab= &ovAB[v0][!rimy][side];
336         OutVertex *cd1= &defs1[defs1aroundmap(around)];
337         outfacet(side^rimy,cd,ab,cd1);
338       } else {
339         cd= &ovC[v0][side];
340       }
341       OutVertex *ab[6];
342       FOR_VPEDGE(e) {
343         ab[e]= invertex2outvertexab(v0,e,side);
344         if (ab[e])
345           K assert(!isnan(ab[e]->p[k]));
346       }
347       outfacets_around(side, cd, 6,ab);
348     }
349   }
350 }     
351
352 /*---------- operations on very output vertex ----------*/
353
354 #define DO_OUTPUT_ARRAY_OUTVERTEX_ARRAY(fn,ovX) \
355   ((fn)(sizeof((ovX))/sizeof(OutVertex), (OutVertex*)(ovX)))
356
357 static void blank_outvertex_array(int n, OutVertex ovX[n]) {
358   int i, k;
359   for (i=0; i<n; i++)
360     K ovX[i].p[k]= NAN;
361 }
362
363 static void blank_outvertices(void) {
364 #define BLANK_OUTPUT_ARRAY(ovX) \
365       DO_OUTPUT_ARRAY_OUTVERTEX_ARRAY(blank_outvertex_array, ovX);
366   OUTPUT_ARRAY_LIST(BLANK_OUTPUT_ARRAY)
367 }
368
369 static void transform_outvertex_array(int n, OutVertex ovX[n]) {
370   int i, k;
371   for (i=0; i<n; i++)
372     K ovX[i].p[k] *= scale;
373   /*
374    *   double min[D3]= thick;
375    *            if (ovX[i].p[k] < min)
376    *        min= ovX[i].p[k];
377    *      for (i=0; i<n; i++) {
378    *  K ovX[k].p[k] -= min;
379    */
380 }
381
382 static void transform_outvertices(void) {
383 #define TRANSFORM_OUTPUT_ARRAY(ovX) \
384       DO_OUTPUT_ARRAY_OUTVERTEX_ARRAY(transform_outvertex_array, ovX);
385   OUTPUT_ARRAY_LIST(TRANSFORM_OUTPUT_ARRAY)
386 }
387
388 /*---------- output file ----------*/
389
390 static void wr(const void *p, size_t sz) {
391   if (fwrite(p,sz,1,stdout) != 1)
392     diee("fwrite");
393 }
394
395 #define WR(x) wr((const void*)&(x), sizeof((x)))
396
397 static void wf(double d) {
398   typedef float ieee754single;
399
400   assert(sizeof(ieee754single)==4);
401   assert(!isnan(d));
402
403 assert(d >= -1e3 && d <= 1e3);
404
405 #if BYTE_ORDER==BIG_ENDIAN
406   union { uint8_t b[4]; ieee754single f; } value;  value.f= d;
407   int i;  for (i=3; i>=0; i--) WR(value.b[i]);
408 #elif BYTE_ORDER==LITTLE_ENDIAN
409   ieee754single f= d;
410   WR(f);
411 #else
412 # error not little or big endian!
413 #endif
414 }
415
416 static uint32_t noutfacets;
417 static uint32_t noutfacets_counted;
418
419 static void outfacet(int rev, const OutVertex *a,
420                      const OutVertex *b, const OutVertex *c) {
421   if (rev) { outfacet(0, c,b,a); return; }
422
423   double normal[D3];
424   int k;
425   
426   K {
427     assert(!isnan(a->p[k]));
428     assert(!isnan(b->p[k]));
429     assert(!isnan(c->p[k]));
430   }
431
432   triangle_normal(normal, a->p, b->p, c->p);
433   double multby= 1/magnD(normal);
434   
435   if (multby > 1e6)
436     return;
437
438   noutfacets++;
439   if (!~noutfacets_counted) return;
440   
441   K normal[k] *= multby;
442
443   K wf(normal[k]);
444   K wf(a->p[k]);
445   K wf(b->p[k]);
446   K wf(c->p[k]);
447   uint16_t attrbc=0;
448   WR(attrbc);
449 }
450
451 static void write_file(void) {
452   static const char header[80]= "#!/usr/bin/meshlab\n" "binary STL file\n";
453
454   if (isatty(1)) fail("will not write binary stl to tty!\n");
455
456   WR(header);
457
458   noutfacets_counted=~(uint32_t)0;
459   noutfacets=0;
460   outfacets();
461   WR(noutfacets);
462   
463   noutfacets_counted= noutfacets;
464   noutfacets=0;
465   outfacets();
466   assert(noutfacets == noutfacets_counted);
467
468   if (fflush(stdout)) diee("fflush stdout");
469 }  
470   
471
472 /*---------- main program etc. ----------*/
473
474 int main(int argc, const char *const *argv) {
475   int r;
476
477   if (argc!=3 || argv[1][0]=='-') { fputs("bad usage\n",stderr); exit(8); }
478   thick= atof(argv[1]);
479   scale= atof(argv[2]) * 0.5; /* circle is unit radius but arg is diameter */
480
481   errno= 0; r= fread(&in,sizeof(in),1,stdin);
482   if (r!=1) diee("fread");
483
484   blank_outvertices();
485   compute_outvertices();
486   transform_outvertices();
487   write_file();
488   if (fclose(stdout)) diee("fclose stdout");
489   return 0;
490 }