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
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
21 * We add thickness by adding vertices as follows:
29 * *E *G | :*E *G *E *G *E *G *E
31 * *D ____________|*D ____________ *D ____________ *D ____________ *D __
33 * / \ / :\ / \ / \ / \
34 * \ *A 4/ \: \5 *B / \ *B / \ *B / \
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 __
49 * / \ / :\ / \ / \ / \
50 * \ *A3 4/ \: \5 5*B / \ *B / \ *B / \
53 * *B \ / *B | \ / *A \ / *A \ / *A
54 * \ / : | \ / . ' \ / \ /
55 * _______ *C _____ :_|___ *C' __________ *C ___________ *C _________
56 * /\ : | 3 /\. 0 /\ /\
57 * / \ : | / \ ` . / \ / \
62 * \ / : | \ / . ' \ / \ /
63 * _______ *C _____ :_|___ *C' __________ *C ___________ *C _________
64 * /\ : | 3 /\. 0 /\ /\
65 * / \ : | / \ ` . / \ / \
66 * *A / \ *A | /4 5\ *B / \ *B / \ *B
67 * / \ 1: / / \ / \ / \
69 * / *B2 2\ / /1 0*A \ / *A \ / *A \ /
70 * \ / \/:/ \ / \ / \ /
71 * *D ____________|*D ____________ *D ___________ *D ___________ *D __
73 * *E *G | :*E *G *E *G *E *G *E
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
83 * The constructions we use are:
85 * A, B: Extend the normal vector of the containing intriangle
86 * from its centroid for thickness.
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.
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.
101 * G: Each F is the mean of those two adjacent Es with the
102 * same angle in their respective Ps.
104 * The outtriangles are:
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.
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).
112 * For each rim edge on each side, the triangle formed by that edge's
113 * ends' Ds and the corresponding A or B.
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.
122 /*---------- declarations and useful subroutines ----------*/
124 #define FOR_SIDE for (side=0; side<2; side++)
131 #define NDEF (NG*2+1)
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 */
140 static double thick; /* in input units */
143 static void outtriangle(int rev, const OutVertex *a,
144 const OutVertex *b, const OutVertex *c);
146 static void normalise_thick(double a[D3]) {
147 /* multiplies a by a scalar so that its magnitude is thick */
149 double multby= thick / magnD(a);
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];
158 K ab[k]= b[k] - a[k];
159 K ac[k]= c[k] - a[k];
163 static OutVertex *invertex2outvertexab(int v0, int e, int side) {
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;
174 if (vref<0) return 0;
175 int sw= VERTICES_SPAN_JOIN_P(v0,vref);
176 return &ovAB[vref][ab^sw][side^sw];
179 /*---------- output vertices ----------*/
181 static void compute_outvertices(void) {
182 int v0,k,side,ab,x,y;
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);
189 K ovAB[v0][ab][0].p[k]= ovAB[v0][ab][1].p[k]= NAN;
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];
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;
213 OutVertex *ovab= invertex2outvertexab(v0,e,side);
214 K adjust[k] += ovab->p[k];
217 K adjust[k] -= in[v0][k];
218 normalise_thick(adjust);
219 K ovC[v0][side].p[k]= in[v0][k] + adjust[k];
222 FOR_RIM_VERTEX(y,x,v0) {
223 double rim[D3], inner[D3], radius_cos[D3], radius_sin[D3];
224 int vback, vfwd, around;
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];
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);
241 assert(vback>=0 && vfwd>=0);
242 K inner[k]= (in[vback][k] + in[vfwd][k]) / 2;
243 K inner[k] -= in[v0][k];
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
248 xprod(radius_cos,rim,inner);
249 xprod(radius_sin,rim,radius_cos);
250 normalise_thick(radius_cos);
251 normalise_thick(radius_sin);
253 for (around=0; around<NDEF; around++) {
254 double angle= around * M_PI / (NDEF-1);
255 K ovDEF[x][!!y][around].p[k]=
257 cos(angle) * radius_cos[k] +
258 sin(angle) * radius_sin[k];
261 FOR_RIM_VERTEX(y,x,v0) {
262 int vfwd= EDGE_END2(v0,0);
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;
273 /*---------- output triangles ----------*/
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 */
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);
288 static int defs_aroundmap_swap(int around) { return NDEF-around; }
289 static int int_identity_function(int i) { return i; }
291 static void outtriangles(void) {
292 int v0,e,side,aroung;
295 OutVertex *defs=0, *defs1=0;
296 int (*defs1aroundmap)(int)=0, rimy=-1;
297 if (RIM_VERTEX_P(v0)) {
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;
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 ];
314 outtriangles_around(rimy, &gs[aroung], 6,surround);
321 int around= side ? NDEF-1 : 0;
323 OutVertex *ab= &ovAB[v0][!rimy][side];
324 OutVertex *cd1= &defs1[defs1aroundmap(around)];
325 outtriangle(side^rimy,cd,ab,cd1);
331 ab[e]= invertex2outvertexab(v0,e,side);
332 outtriangles_around(side, cd, 6,ab);
337 /*---------- transformation (scale and perhaps shift) ----------*/
339 static void scaleshift_outvertex_array(int n, OutVertex ovX[n]) {
342 K ovX[i].p[k] *= scale;
344 * double min[D3]= thick;
345 * if (ovX[i].p[k] < min)
347 * for (i=0; i<n; i++) {
348 * K ovX[k].p[k] -= min;
352 #define SCALESHIFT_OUTVERTEX_ARRAY(ovX) \
353 scaleshift_outvertex_array(sizeof((ovX))/sizeof(OutVertex),(OutVertex*)(ovX))
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);
362 /*---------- output file ----------*/
364 static void wr(const void *p, size_t sz) {
365 if (fwrite(p,sz,1,stdout) != 1)
369 #define WR(x) wr((const void*)&(x), sizeof((x)))
371 static void wf(double d) {
372 typedef float ieee754single;
374 assert(sizeof(ieee754single)==4);
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)
383 # error not little or big endian!
387 static uint32_t nouttriangles;
388 static uint32_t nouttriangles_counted;
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; }
394 if (!~nouttriangles_counted) return;
399 triangle_normal(normal, a->p, b->p, c->p);
400 double multby= 1/magnD(normal);
401 K normal[k] *= multby;
411 static void write_file(void) {
412 static const char header[80]= "#!/usr/bin/meshlab\n" "binary STL file\n";
414 if (isatty(1)) fail("will not write binary stl to tty!");
418 nouttriangles_counted=~(uint32_t)0;
423 nouttriangles_counted= nouttriangles;
426 assert(nouttriangles == nouttriangles_counted);
428 if (fflush(stdout)) diee("fflush stdout");
432 /*---------- main program etc. ----------*/
434 int main(int argc, const char *const *argv) {
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 */
441 errno= 0; r= fread(&in,sizeof(in),1,stdin);
442 if (r!=1) diee("fread");
444 compute_outvertices();
445 scaleshift_outvertices();
447 if (fclose(stdout)) diee("fclose stdout");