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
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
14 * ./output-64 <dense-64.cfm 0.1 50 >t.stl
15 * ./output-125 <best-125.cfm >t.stl 0.1 50
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
26 * We add thickness by adding vertices as follows:
34 * *E *G | :*E *G *E *G *E *G *E
36 * *D ____________|*D ____________ *D ____________ *D ____________ *D __
38 * / \ / :\ / \ / \ / \
39 * \ *A 4/ \: \5 *B / \ *B / \ *B / \
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 __
54 * / \ / :\ / \ / \ / \
55 * \ *A3 4/ \: \5 5*B / \ *B / \ *B / \
58 * *B \ / *B | \ / *A \ / *A \ / *A
59 * \ / : | \ / . ' \ / \ /
60 * _______ *C _____ :_|___ *C' __________ *C ___________ *C _________
68 * _______ *C _____ :_|___ *C ___________ *C ___________ *C _________
71 * *A / \ *A | /4 5\ *B / \ *B / \ *B
72 * / \ 1: / / \ / \ / \
74 * / *B2 2\ / /1 0*A \ / *A \ / *A \ /
75 * \ / \/:/ \ / \ / \ /
76 * *D ____________|*D ____________ *D ___________ *D ___________ *D __
78 * *E *G | :*E *G *E *G *E *G *E
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
88 * The constructions we use are:
90 * A, B: Extend the normal vector of the containing intriangle
91 * from its centroid for thickness.
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.
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.
106 * G: Each F is the mean of those two adjacent Es with the
107 * same angle in their respective Ps.
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.
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.
121 /*---------- declarations and useful subroutines ----------*/
123 #define FOR_SIDE for (side=0; side<2; side++)
130 #define NDEF (NG*2+1)
132 #define OUTPUT_ARRAY_LIST(DO_OUTPUT_ARRAY) \
133 DO_OUTPUT_ARRAY(ovAB) \
134 DO_OUTPUT_ARRAY(ovC) \
135 DO_OUTPUT_ARRAY(ovDEF) \
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 */
145 static double thick; /* in input units */
148 static void outfacet(int rev, const OutVertex *a,
149 const OutVertex *b, const OutVertex *c);
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; }
155 static void normalise_thick(double a[D3]) {
156 /* multiplies a by a scalar so that its magnitude is thick */
158 double multby= thick / magnD(a);
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];
167 K ab[k]= b[k] - a[k];
168 K ac[k]= c[k] - a[k];
172 static OutVertex *invertex2outvertexab(int v0, int e, int side) {
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;
183 if (vchk<0) return 0;
184 int sw= vertices_span_join_p(v0,vref);
185 return &ovAB[vref][ab^sw][side^sw];
188 /*---------- output vertices ----------*/
190 #define Ok(ov, value) ((ov).p[k]= outvertex_coord_check(value))
192 static double outvertex_coord_check(double value) {
193 assert(-10 < value && value < 10);
197 static void compute_outvertices(void) {
198 int v0,k,side,ab,x,y;
200 FOR_VERTEX(v0, INNER) {
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);
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]);
214 FOR_VERTEX(v0, INNER) {
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)
225 OutVertex *ovab= invertex2outvertexab(v0,e,side);
227 assert(!isnan(ovab->p[k]));
228 adjust[k] += ovab->p[k];
232 K adjust[k] -= in[v0][k];
233 normalise_thick(adjust);
234 K Ok(ovC[v0][side], in[v0][k] + adjust[k]);
237 FOR_RIM_VERTEX(y,x,v0, INNER) {
238 double rim[D3], inner[D3], radius_cos[D3], radius_sin[D3];
239 int vback, vfwd, around;
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];
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);
256 assert(vback>=0 && vfwd>=0);
257 K inner[k]= (in[vback][k] + in[vfwd][k]) / 2;
258 K inner[k] -= in[v0][k];
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
263 xprod(radius_cos,rim,inner);
264 xprod(radius_sin,rim,radius_cos);
265 normalise_thick(radius_cos);
266 normalise_thick(radius_sin);
268 int_map *around_map= y ? int_identity_function : defs_aroundmap_swap;
270 for (around=0; around<NDEF; around++) {
271 double angle= around_map(around) * M_PI / (NDEF-1);
272 K Ok(ovDEF[x][!!y][around],
274 cos(angle) * radius_cos[k] +
275 sin(angle) * radius_sin[k]);
278 FOR_RIM_VERTEX(y,x,v0, INNER) {
279 int vfwd= EDGE_END2(v0,0);
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]));
292 /*---------- output facets ----------*/
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 */
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);
307 static OutVertex *invertex2outvertexcd(v0,side) {
308 if (!RIM_VERTEX_P(v0)) return &ovC[v0][side];
310 int around= side ? NDEF-1 : 0;
311 int rimy= !!(v0 & ~XMASK);
312 return &ovDEF[v0 & XMASK][rimy][around];
315 static void outfacets(void) {
316 int v0,e,side,aroung;
318 FOR_VERTEX(v0, INNER) {
319 OutVertex *defs=0, *defs1=0;
321 int_map *defs1aroundmap= 0;
322 if (RIM_VERTEX_P(v0)) {
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;
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 ];
339 outfacets_around(rimy, &gs[aroung], 6,surround);
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;
350 invertex2outvertexcd(v0,side),
351 invertex2outvertexcd(v1,side^vertices_span_join_p(v0,v1)),
352 invertex2outvertexcd(v2,side^vertices_span_join_p(v0,v2)));
358 /*---------- operations on very output vertex ----------*/
360 #define DO_OUTPUT_ARRAY_OUTVERTEX_ARRAY(fn,ovX) \
361 ((fn)(sizeof((ovX))/sizeof(OutVertex), (OutVertex*)(ovX)))
363 static void blank_outvertex_array(int n, OutVertex ovX[n]) {
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)
375 static void transform_outvertex_array(int n, OutVertex ovX[n]) {
377 for (i=0; i<n; i++) {
380 K ovX[i].p[k] *= scale;
383 * double min[D3]= thick;
384 * if (ovX[i].p[k] < min)
386 * for (i=0; i<n; i++) {
387 * K ovX[k].p[k] -= min;
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)
397 /*---------- output file ----------*/
399 static void wr(const void *p, size_t sz) {
400 if (fwrite(p,sz,1,stdout) != 1)
404 #define WR(x) wr((const void*)&(x), sizeof((x)))
406 static void wf(double d) {
407 typedef float ieee754single;
409 assert(sizeof(ieee754single)==4);
412 assert(d >= -1e3 && d <= 1e3);
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
421 # error not little or big endian!
425 static uint32_t noutfacets;
426 static uint32_t noutfacets_counted;
427 static long badfacets;
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; }
437 assert(!isnan(a->p[k]));
438 assert(!isnan(b->p[k]));
439 assert(!isnan(c->p[k]));
442 triangle_normal(normal, a->p, b->p, c->p);
443 double multby= 1/magnD(normal);
451 if (!~noutfacets_counted) return;
453 K normal[k] *= multby;
463 static void write_file(void) {
464 static const char header[80]= "#!/usr/bin/meshlab\n" "binary STL file\n";
466 if (isatty(1)) fail("will not write binary stl to tty!\n");
470 noutfacets_counted=~(uint32_t)0;
475 noutfacets_counted= noutfacets;
478 assert(noutfacets == noutfacets_counted);
480 if (fflush(stdout)) diee("fflush stdout");
483 fprintf(stderr,"%ld degenerate facets!\n",badfacets);
489 /*---------- main program etc. ----------*/
491 int main(int argc, const char *const *argv) {
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 */
498 errno= 0; r= fread(&in,sizeof(in),1,stdin);
499 if (r!=1) diee("fread");
502 compute_outvertices();
503 transform_outvertices();
505 if (fclose(stdout)) diee("fclose stdout");