ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/gensurf.c
Revision: 2.21
Committed: Sun Mar 11 22:27:17 2012 UTC (12 years ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.20: +17 -23 lines
Log Message:
Fixed bug in .OBJ output for degenerate normals

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: gensurf.c,v 2.20 2011/08/15 19:48:06 greg Exp $";
3 #endif
4 /*
5 * gensurf.c - program to generate functional surfaces
6 *
7 * Parametric functions x(s,t), y(s,t) and z(s,t)
8 * specify the surface, which is tesselated into an m by n
9 * array of paired triangles.
10 * The surface normal is defined by the right hand
11 * rule applied to (s,t).
12 *
13 * 4/3/87
14 *
15 * 4/16/02 Added conditional vertex output
16 */
17
18 #include "standard.h"
19
20 #include "paths.h"
21 #include "resolu.h"
22 #include "rterror.h"
23 #include "calcomp.h"
24
25 char XNAME[] = "X`SYS"; /* x function name */
26 char YNAME[] = "Y`SYS"; /* y function name */
27 char ZNAME[] = "Z`SYS"; /* z function name */
28
29 char VNAME[] = "valid"; /* valid vertex name */
30
31 #define ABS(x) ((x)>=0 ? (x) : -(x))
32
33 #define ZEROVECT(v) (DOT(v,v) <= FTINY*FTINY)
34
35 #define pvect(p) printf(vformat, (p)[0], (p)[1], (p)[2])
36
37 char vformat[] = "%18.12g %18.12g %18.12g\n";
38 char tsargs[] = "4 surf_dx surf_dy surf_dz surf.cal";
39 char texname[] = "Phong";
40
41 int smooth = 0; /* apply smoothing? */
42 int objout = 0; /* output .OBJ format? */
43
44 char *modname, *surfname;
45
46 /* recorded data flags */
47 #define HASBORDER 01
48 #define TRIPLETS 02
49 /* a data structure */
50 struct {
51 int flags; /* data type */
52 short m, n; /* number of s and t values */
53 RREAL *data; /* the data itself, s major sort */
54 } datarec; /* our recorded data */
55
56 /* XXX this is redundant with rt/noise3.c, should go to a library */
57 double l_hermite(), l_bezier(), l_bspline(), l_dataval();
58
59 typedef struct {
60 int valid; /* point is valid (vertex number) */
61 int nvalid; /* normal is valid */
62 FVECT p; /* vertex position */
63 FVECT n; /* average normal */
64 RREAL uv[2]; /* (u,v) position */
65 } POINT;
66
67
68 void loaddata(char *file, int m, int n, int pointsize);
69 double l_dataval(char *nam);
70 void putobjrow(POINT *rp, int n);
71 void putsquare(POINT *p0, POINT *p1, POINT *p2, POINT *p3);
72 void comprow(double s, POINT *row, int siz);
73 void compnorms(POINT *r0, POINT *r1, POINT *r2, int siz);
74 int norminterp(FVECT resmat[4], POINT *p0, POINT *p1, POINT *p2, POINT *p3);
75
76
77 int
78 main(argc, argv)
79 int argc;
80 char *argv[];
81 {
82 POINT *row0, *row1, *row2, *rp;
83 int i, j, m, n;
84 char stmp[256];
85
86 varset("PI", ':', PI);
87 funset("hermite", 5, ':', l_hermite);
88 funset("bezier", 5, ':', l_bezier);
89 funset("bspline", 5, ':', l_bspline);
90
91 if (argc < 8)
92 goto userror;
93
94 for (i = 8; i < argc; i++)
95 if (!strcmp(argv[i], "-e"))
96 scompile(argv[++i], NULL, 0);
97 else if (!strcmp(argv[i], "-f"))
98 fcompile(argv[++i]);
99 else if (!strcmp(argv[i], "-s"))
100 smooth++;
101 else if (!strcmp(argv[i], "-o"))
102 objout++;
103 else
104 goto userror;
105
106 modname = argv[1];
107 surfname = argv[2];
108 m = atoi(argv[6]);
109 n = atoi(argv[7]);
110 if (m <= 0 || n <= 0)
111 goto userror;
112 if (!strcmp(argv[5], "-") || access(argv[5], 4) == 0) { /* file? */
113 funset(ZNAME, 2, ':', l_dataval);
114 if (!strcmp(argv[5],argv[3]) && !strcmp(argv[5],argv[4])) {
115 loaddata(argv[5], m, n, 3);
116 funset(XNAME, 2, ':', l_dataval);
117 funset(YNAME, 2, ':', l_dataval);
118 } else {
119 loaddata(argv[5], m, n, 1);
120 sprintf(stmp, "%s(s,t)=%s;", XNAME, argv[3]);
121 scompile(stmp, NULL, 0);
122 sprintf(stmp, "%s(s,t)=%s;", YNAME, argv[4]);
123 scompile(stmp, NULL, 0);
124 }
125 } else {
126 sprintf(stmp, "%s(s,t)=%s;", XNAME, argv[3]);
127 scompile(stmp, NULL, 0);
128 sprintf(stmp, "%s(s,t)=%s;", YNAME, argv[4]);
129 scompile(stmp, NULL, 0);
130 sprintf(stmp, "%s(s,t)=%s;", ZNAME, argv[5]);
131 scompile(stmp, NULL, 0);
132 }
133 row0 = (POINT *)malloc((n+3)*sizeof(POINT));
134 row1 = (POINT *)malloc((n+3)*sizeof(POINT));
135 row2 = (POINT *)malloc((n+3)*sizeof(POINT));
136 if (row0 == NULL || row1 == NULL || row2 == NULL) {
137 fprintf(stderr, "%s: out of memory\n", argv[0]);
138 quit(1);
139 }
140 row0++; row1++; row2++;
141 /* print header */
142 fputs("# ", stdout);
143 printargs(argc, argv, stdout);
144 eclock = 0;
145 /* initialize */
146 comprow(-1.0/m, row0, n);
147 comprow(0.0, row1, n);
148 comprow(1.0/m, row2, n);
149 compnorms(row0, row1, row2, n);
150 if (objout) {
151 printf("\nusemtl %s\n\n", modname);
152 putobjrow(row1, n);
153 }
154 /* for each row */
155 for (i = 0; i < m; i++) {
156 /* compute next row */
157 rp = row0;
158 row0 = row1;
159 row1 = row2;
160 row2 = rp;
161 comprow((double)(i+2)/m, row2, n);
162 compnorms(row0, row1, row2, n);
163 if (objout)
164 putobjrow(row1, n);
165
166 for (j = 0; j < n; j++) {
167 int orient = (j & 1);
168 /* put polygons */
169 if (!(row0[j].valid && row1[j+1].valid))
170 orient = 1;
171 else if (!(row1[j].valid && row0[j+1].valid))
172 orient = 0;
173 if (orient)
174 putsquare(&row0[j], &row1[j],
175 &row0[j+1], &row1[j+1]);
176 else
177 putsquare(&row1[j], &row1[j+1],
178 &row0[j], &row0[j+1]);
179 }
180 }
181
182 return 0;
183
184 userror:
185 fprintf(stderr, "Usage: %s material name ", argv[0]);
186 fprintf(stderr, "x(s,t) y(s,t) z(s,t) m n [-s][-o][-e expr][-f file]\n");
187 return 1;
188 }
189
190
191 void
192 loaddata( /* load point data from file */
193 char *file,
194 int m,
195 int n,
196 int pointsize
197 )
198 {
199 FILE *fp;
200 char word[64];
201 register int size;
202 register RREAL *dp;
203
204 datarec.flags = HASBORDER; /* assume border values */
205 datarec.m = m+1;
206 datarec.n = n+1;
207 size = datarec.m*datarec.n*pointsize;
208 if (pointsize == 3)
209 datarec.flags |= TRIPLETS;
210 dp = (RREAL *)malloc(size*sizeof(RREAL));
211 if ((datarec.data = dp) == NULL) {
212 fputs("Out of memory\n", stderr);
213 exit(1);
214 }
215 if (!strcmp(file, "-")) {
216 file = "<stdin>";
217 fp = stdin;
218 } else if ((fp = fopen(file, "r")) == NULL) {
219 fputs(file, stderr);
220 fputs(": cannot open\n", stderr);
221 exit(1);
222 }
223 while (size > 0 && fgetword(word, sizeof(word), fp) != NULL) {
224 if (!isflt(word)) {
225 fprintf(stderr, "%s: garbled data value: %s\n",
226 file, word);
227 exit(1);
228 }
229 *dp++ = atof(word);
230 size--;
231 }
232 if (size == (m+n+1)*pointsize) { /* no border after all */
233 dp = (RREAL *)realloc(datarec.data,
234 m*n*pointsize*sizeof(RREAL));
235 if (dp != NULL)
236 datarec.data = dp;
237 datarec.flags &= ~HASBORDER;
238 datarec.m = m;
239 datarec.n = n;
240 size = 0;
241 }
242 if (datarec.m < 2 || datarec.n < 2 || size != 0 ||
243 fgetword(word, sizeof(word), fp) != NULL) {
244 fputs(file, stderr);
245 fputs(": bad number of data points\n", stderr);
246 exit(1);
247 }
248 fclose(fp);
249 }
250
251
252 double
253 l_dataval( /* return recorded data value */
254 char *nam
255 )
256 {
257 double u, v;
258 register int i, j;
259 register RREAL *dp;
260 double d00, d01, d10, d11;
261 /* compute coordinates */
262 u = argument(1); v = argument(2);
263 if (datarec.flags & HASBORDER) {
264 i = u *= datarec.m-1;
265 j = v *= datarec.n-1;
266 } else {
267 i = u = u*datarec.m - .5;
268 j = v = v*datarec.n - .5;
269 }
270 if (i < 0) i = 0;
271 else if (i > datarec.m-2) i = datarec.m-2;
272 if (j < 0) j = 0;
273 else if (j > datarec.n-2) j = datarec.n-2;
274 /* compute value */
275 if (datarec.flags & TRIPLETS) {
276 dp = datarec.data + 3*(j*datarec.m + i);
277 if (nam == ZNAME)
278 dp += 2;
279 else if (nam == YNAME)
280 dp++;
281 d00 = dp[0]; d01 = dp[3];
282 dp += 3*datarec.m;
283 d10 = dp[0]; d11 = dp[3];
284 } else {
285 dp = datarec.data + j*datarec.m + i;
286 d00 = dp[0]; d01 = dp[1];
287 dp += datarec.m;
288 d10 = dp[0]; d11 = dp[1];
289 }
290 /* bilinear interpolation */
291 return((j+1-v)*((i+1-u)*d00+(u-i)*d01)+(v-j)*((i+1-u)*d10+(u-i)*d11));
292 }
293
294
295 void
296 putobjrow( /* output vertex row to .OBJ */
297 register POINT *rp,
298 int n
299 )
300 {
301 static int nverts = 0;
302 static int nnorms = 0;
303
304 for ( ; n-- >= 0; rp++) {
305 if (!rp->valid)
306 continue;
307 fputs("v ", stdout);
308 pvect(rp->p);
309 if (smooth && !ZEROVECT(rp->n)) {
310 printf("\tvn %.9g %.9g %.9g\n",
311 rp->n[0], rp->n[1], rp->n[2]);
312 rp->nvalid = ++nnorms;
313 } else
314 rp->nvalid = 0;
315 printf("\tvt %.9g %.9g\n", rp->uv[0], rp->uv[1]);
316 rp->valid = ++nverts;
317 }
318 }
319
320
321 void
322 putsquare( /* put out a square */
323 POINT *p0,
324 POINT *p1,
325 POINT *p2,
326 POINT *p3
327 )
328 {
329 static int nout = 0;
330 FVECT norm[4];
331 int axis;
332 FVECT v1, v2, vc1, vc2;
333 int ok1, ok2;
334 /* compute exact normals */
335 ok1 = (p0->valid && p1->valid && p2->valid);
336 if (ok1) {
337 VSUB(v1, p1->p, p0->p);
338 VSUB(v2, p2->p, p0->p);
339 fcross(vc1, v1, v2);
340 ok1 = (normalize(vc1) != 0.0);
341 }
342 ok2 = (p1->valid && p2->valid && p3->valid);
343 if (ok2) {
344 VSUB(v1, p2->p, p3->p);
345 VSUB(v2, p1->p, p3->p);
346 fcross(vc2, v1, v2);
347 ok2 = (normalize(vc2) != 0.0);
348 }
349 if (!(ok1 | ok2))
350 return;
351 if (objout) { /* output .OBJ faces */
352 if (ok1 & ok2 && fdot(vc1,vc2) >= 1.0-FTINY*FTINY) {
353 printf("f %d/%d/%d %d/%d/%d %d/%d/%d %d/%d/%d\n",
354 p0->valid, p0->valid, p0->nvalid,
355 p1->valid, p1->valid, p1->nvalid,
356 p3->valid, p3->valid, p3->nvalid,
357 p2->valid, p2->valid, p2->nvalid);
358 return;
359 }
360 if (ok1)
361 printf("f %d/%d/%d %d/%d/%d %d/%d/%d\n",
362 p0->valid, p0->valid, p0->nvalid,
363 p1->valid, p1->valid, p1->nvalid,
364 p2->valid, p2->valid, p2->nvalid);
365 if (ok2)
366 printf("f %d/%d/%d %d/%d/%d %d/%d/%d\n",
367 p2->valid, p2->valid, p2->nvalid,
368 p1->valid, p1->valid, p1->nvalid,
369 p3->valid, p3->valid, p3->nvalid);
370 return;
371 }
372 /* compute normal interpolation */
373 axis = norminterp(norm, p0, p1, p2, p3);
374
375 /* put out quadrilateral? */
376 if (ok1 & ok2 && fdot(vc1,vc2) >= 1.0-FTINY*FTINY) {
377 printf("\n%s ", modname);
378 if (axis != -1) {
379 printf("texfunc %s\n%s\n", texname, tsargs);
380 printf("0\n13\t%d\n", axis);
381 pvect(norm[0]);
382 pvect(norm[1]);
383 pvect(norm[2]);
384 fvsum(v1, norm[3], vc1, -0.5);
385 fvsum(v1, v1, vc2, -0.5);
386 pvect(v1);
387 printf("\n%s ", texname);
388 }
389 printf("polygon %s.%d\n", surfname, ++nout);
390 printf("0\n0\n12\n");
391 pvect(p0->p);
392 pvect(p1->p);
393 pvect(p3->p);
394 pvect(p2->p);
395 return;
396 }
397 /* put out triangles? */
398 if (ok1) {
399 printf("\n%s ", modname);
400 if (axis != -1) {
401 printf("texfunc %s\n%s\n", texname, tsargs);
402 printf("0\n13\t%d\n", axis);
403 pvect(norm[0]);
404 pvect(norm[1]);
405 pvect(norm[2]);
406 fvsum(v1, norm[3], vc1, -1.0);
407 pvect(v1);
408 printf("\n%s ", texname);
409 }
410 printf("polygon %s.%d\n", surfname, ++nout);
411 printf("0\n0\n9\n");
412 pvect(p0->p);
413 pvect(p1->p);
414 pvect(p2->p);
415 }
416 if (ok2) {
417 printf("\n%s ", modname);
418 if (axis != -1) {
419 printf("texfunc %s\n%s\n", texname, tsargs);
420 printf("0\n13\t%d\n", axis);
421 pvect(norm[0]);
422 pvect(norm[1]);
423 pvect(norm[2]);
424 fvsum(v2, norm[3], vc2, -1.0);
425 pvect(v2);
426 printf("\n%s ", texname);
427 }
428 printf("polygon %s.%d\n", surfname, ++nout);
429 printf("0\n0\n9\n");
430 pvect(p2->p);
431 pvect(p1->p);
432 pvect(p3->p);
433 }
434 }
435
436
437 void
438 comprow( /* compute row of values */
439 double s,
440 register POINT *row,
441 int siz
442 )
443 {
444 double st[2];
445 int end;
446 int checkvalid;
447 register int i;
448
449 if (smooth) {
450 i = -1; /* compute one past each end */
451 end = siz+1;
452 } else {
453 if (s < -FTINY || s > 1.0+FTINY)
454 return;
455 i = 0;
456 end = siz;
457 }
458 st[0] = s;
459 checkvalid = (fundefined(VNAME) == 2);
460 while (i <= end) {
461 st[1] = (double)i/siz;
462 if (checkvalid && funvalue(VNAME, 2, st) <= 0.0) {
463 row[i].valid = 0;
464 row[i].p[0] = row[i].p[1] = row[i].p[2] = 0.0;
465 row[i].uv[0] = row[i].uv[1] = 0.0;
466 } else {
467 row[i].valid = 1;
468 row[i].p[0] = funvalue(XNAME, 2, st);
469 row[i].p[1] = funvalue(YNAME, 2, st);
470 row[i].p[2] = funvalue(ZNAME, 2, st);
471 row[i].uv[0] = st[0];
472 row[i].uv[1] = st[1];
473 }
474 i++;
475 }
476 }
477
478
479 void
480 compnorms( /* compute row of averaged normals */
481 register POINT *r0,
482 register POINT *r1,
483 register POINT *r2,
484 int siz
485 )
486 {
487 FVECT v1, v2;
488
489 if (!smooth) /* not needed if no smoothing */
490 return;
491 /* compute row 1 normals */
492 while (siz-- >= 0) {
493 if (!r1[0].valid)
494 continue;
495 if (!r0[0].valid) {
496 if (!r2[0].valid) {
497 r1[0].n[0] = r1[0].n[1] = r1[0].n[2] = 0.0;
498 continue;
499 }
500 fvsum(v1, r2[0].p, r1[0].p, -1.0);
501 } else if (!r2[0].valid)
502 fvsum(v1, r1[0].p, r0[0].p, -1.0);
503 else
504 fvsum(v1, r2[0].p, r0[0].p, -1.0);
505 if (!r1[-1].valid) {
506 if (!r1[1].valid) {
507 r1[0].n[0] = r1[0].n[1] = r1[0].n[2] = 0.0;
508 continue;
509 }
510 fvsum(v2, r1[1].p, r1[0].p, -1.0);
511 } else if (!r1[1].valid)
512 fvsum(v2, r1[0].p, r1[-1].p, -1.0);
513 else
514 fvsum(v2, r1[1].p, r1[-1].p, -1.0);
515 fcross(r1[0].n, v1, v2);
516 normalize(r1[0].n);
517 r0++; r1++; r2++;
518 }
519 }
520
521
522 int
523 norminterp( /* compute normal interpolation */
524 register FVECT resmat[4],
525 POINT *p0,
526 POINT *p1,
527 POINT *p2,
528 POINT *p3
529 )
530 {
531 #define u ((ax+1)%3)
532 #define v ((ax+2)%3)
533
534 register int ax;
535 MAT4 eqnmat;
536 FVECT v1;
537 register int i, j;
538
539 if (!smooth) /* no interpolation if no smoothing */
540 return(-1);
541 /* find dominant axis */
542 VCOPY(v1, p0->n);
543 fvsum(v1, v1, p1->n, 1.0);
544 fvsum(v1, v1, p2->n, 1.0);
545 fvsum(v1, v1, p3->n, 1.0);
546 ax = ABS(v1[0]) > ABS(v1[1]) ? 0 : 1;
547 ax = ABS(v1[ax]) > ABS(v1[2]) ? ax : 2;
548 /* assign equation matrix */
549 eqnmat[0][0] = p0->p[u]*p0->p[v];
550 eqnmat[0][1] = p0->p[u];
551 eqnmat[0][2] = p0->p[v];
552 eqnmat[0][3] = 1.0;
553 eqnmat[1][0] = p1->p[u]*p1->p[v];
554 eqnmat[1][1] = p1->p[u];
555 eqnmat[1][2] = p1->p[v];
556 eqnmat[1][3] = 1.0;
557 eqnmat[2][0] = p2->p[u]*p2->p[v];
558 eqnmat[2][1] = p2->p[u];
559 eqnmat[2][2] = p2->p[v];
560 eqnmat[2][3] = 1.0;
561 eqnmat[3][0] = p3->p[u]*p3->p[v];
562 eqnmat[3][1] = p3->p[u];
563 eqnmat[3][2] = p3->p[v];
564 eqnmat[3][3] = 1.0;
565 /* invert matrix (solve system) */
566 if (!invmat4(eqnmat, eqnmat))
567 return(-1); /* no solution */
568 /* compute result matrix */
569 for (j = 0; j < 4; j++)
570 for (i = 0; i < 3; i++)
571 resmat[j][i] = eqnmat[j][0]*p0->n[i] +
572 eqnmat[j][1]*p1->n[i] +
573 eqnmat[j][2]*p2->n[i] +
574 eqnmat[j][3]*p3->n[i];
575 return(ax);
576
577 #undef u
578 #undef v
579 }
580
581
582 double
583 l_hermite(char *nm)
584 {
585 double t;
586
587 t = argument(5);
588 return( argument(1)*((2.0*t-3.0)*t*t+1.0) +
589 argument(2)*(-2.0*t+3.0)*t*t +
590 argument(3)*((t-2.0)*t+1.0)*t +
591 argument(4)*(t-1.0)*t*t );
592 }
593
594
595 double
596 l_bezier(char *nm)
597 {
598 double t;
599
600 t = argument(5);
601 return( argument(1) * (1.+t*(-3.+t*(3.-t))) +
602 argument(2) * 3.*t*(1.+t*(-2.+t)) +
603 argument(3) * 3.*t*t*(1.-t) +
604 argument(4) * t*t*t );
605 }
606
607
608 double
609 l_bspline(char *nm)
610 {
611 double t;
612
613 t = argument(5);
614 return( argument(1) * (1./6.+t*(-1./2.+t*(1./2.-1./6.*t))) +
615 argument(2) * (2./3.+t*t*(-1.+1./2.*t)) +
616 argument(3) * (1./6.+t*(1./2.+t*(1./2.-1./2.*t))) +
617 argument(4) * (1./6.*t*t*t) );
618 }