ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/gensurf.c
Revision: 1.10
Committed: Fri Jul 20 14:35:59 1990 UTC (33 years, 8 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 1.9: +4 -4 lines
Log Message:
changed arguments to scompile()

File Contents

# Content
1 #ifndef lint
2 static char SCCSid[] = "$SunId$ LBL";
3 #endif
4
5 /* Copyright (c) 1989 Regents of the University of California */
6
7 /*
8 * gensurf.c - program to generate functional surfaces
9 *
10 * Parametric functions x(s,t), y(s,t) and z(s,t)
11 * specify the surface, which is tesselated into an m by n
12 * array of paired triangles.
13 * The surface normal is defined by the right hand
14 * rule applied to (s,t).
15 *
16 * 4/3/87
17 */
18
19 #include "standard.h"
20
21 #define XNAME "X_" /* x function name */
22 #define YNAME "Y_" /* y function name */
23 #define ZNAME "Z_" /* z function name */
24
25 #define ABS(x) ((x)>=0 ? (x) : -(x))
26
27 #define pvect(p) printf(vformat, (p)[0], (p)[1], (p)[2])
28
29 char vformat[] = "%15.9g %15.9g %15.9g\n";
30 char tsargs[] = "4 surf_dx surf_dy surf_dz surf.cal\n";
31 char texname[] = "Phong";
32
33 int smooth = 0; /* apply smoothing? */
34
35 char *modname, *surfname;
36
37 double funvalue(), l_hermite(), l_bezier(), l_bspline(), argument();
38
39 typedef struct {
40 FVECT p; /* vertex position */
41 FVECT n; /* average normal */
42 } POINT;
43
44
45 main(argc, argv)
46 int argc;
47 char *argv[];
48 {
49 extern long eclock;
50 POINT *row0, *row1, *row2, *rp;
51 int i, j, m, n;
52 char stmp[256];
53
54 varset("PI", PI);
55 funset("hermite", 5, l_hermite);
56 funset("bezier", 5, l_bezier);
57 funset("bspline", 5, l_bspline);
58
59 if (argc < 8)
60 goto userror;
61
62 for (i = 8; i < argc; i++)
63 if (!strcmp(argv[i], "-e"))
64 scompile(argv[++i], NULL, 0);
65 else if (!strcmp(argv[i], "-f"))
66 fcompile(argv[++i]);
67 else if (!strcmp(argv[i], "-s"))
68 smooth++;
69 else
70 goto userror;
71
72 modname = argv[1];
73 surfname = argv[2];
74 sprintf(stmp, "%s(s,t)=%s;", XNAME, argv[3]);
75 scompile(stmp, NULL, 0);
76 sprintf(stmp, "%s(s,t)=%s;", YNAME, argv[4]);
77 scompile(stmp, NULL, 0);
78 sprintf(stmp, "%s(s,t)=%s;", ZNAME, argv[5]);
79 scompile(stmp, NULL, 0);
80 m = atoi(argv[6]);
81 n = atoi(argv[7]);
82 if (m <= 0 || n <= 0)
83 goto userror;
84
85 row0 = (POINT *)malloc((n+3)*sizeof(POINT));
86 row1 = (POINT *)malloc((n+3)*sizeof(POINT));
87 row2 = (POINT *)malloc((n+3)*sizeof(POINT));
88 if (row0 == NULL || row1 == NULL || row2 == NULL) {
89 fprintf(stderr, "%s: out of memory\n", argv[0]);
90 quit(1);
91 }
92 row0++; row1++; row2++;
93 /* print header */
94 printhead(argc, argv);
95 eclock = 0;
96 /* initialize */
97 comprow(-1.0/m, row0, n);
98 comprow(0.0, row1, n);
99 comprow(1.0/m, row2, n);
100 compnorms(row0, row1, row2, n);
101 /* for each row */
102 for (i = 0; i < m; i++) {
103 /* compute next row */
104 rp = row0;
105 row0 = row1;
106 row1 = row2;
107 row2 = rp;
108 comprow((double)(i+2)/m, row2, n);
109 compnorms(row0, row1, row2, n);
110
111 for (j = 0; j < n; j++) {
112 /* put polygons */
113 if ((i+j) & 1)
114 putsquare(&row0[j], &row1[j],
115 &row0[j+1], &row1[j+1]);
116 else
117 putsquare(&row1[j], &row1[j+1],
118 &row0[j], &row0[j+1]);
119 }
120 }
121
122 quit(0);
123
124 userror:
125 fprintf(stderr, "Usage: %s material name ", argv[0]);
126 fprintf(stderr, "x(s,t) y(s,t) z(s,t) m n [-s][-e expr][-f file]\n");
127 quit(1);
128 }
129
130
131 putsquare(p0, p1, p2, p3) /* put out a square */
132 POINT *p0, *p1, *p2, *p3;
133 {
134 static int nout = 0;
135 FVECT norm[4];
136 int axis;
137 FVECT v1, v2, vc1, vc2;
138 int ok1, ok2;
139 /* compute exact normals */
140 fvsum(v1, p1->p, p0->p, -1.0);
141 fvsum(v2, p2->p, p0->p, -1.0);
142 fcross(vc1, v1, v2);
143 ok1 = normalize(vc1) != 0.0;
144 fvsum(v1, p2->p, p3->p, -1.0);
145 fvsum(v2, p1->p, p3->p, -1.0);
146 fcross(vc2, v1, v2);
147 ok2 = normalize(vc2) != 0.0;
148 if (!(ok1 | ok2))
149 return;
150 /* compute normal interpolation */
151 axis = norminterp(norm, p0, p1, p2, p3);
152
153 /* put out quadrilateral? */
154 if (ok1 & ok2 && fdot(vc1,vc2) >= 1.0-FTINY*FTINY) {
155 printf("\n%s ", modname);
156 if (axis != -1) {
157 printf("texfunc %s\n", texname);
158 printf(tsargs);
159 printf("0\n13\t%d\n", axis);
160 pvect(norm[0]);
161 pvect(norm[1]);
162 pvect(norm[2]);
163 fvsum(v1, norm[3], vc1, -0.5);
164 fvsum(v1, v1, vc2, -0.5);
165 pvect(v1);
166 printf("\n%s ", texname);
167 }
168 printf("polygon %s.%d\n", surfname, ++nout);
169 printf("0\n0\n12\n");
170 pvect(p0->p);
171 pvect(p1->p);
172 pvect(p3->p);
173 pvect(p2->p);
174 return;
175 }
176 /* put out triangles? */
177 if (ok1) {
178 printf("\n%s ", modname);
179 if (axis != -1) {
180 printf("texfunc %s\n", texname);
181 printf(tsargs);
182 printf("0\n13\t%d\n", axis);
183 pvect(norm[0]);
184 pvect(norm[1]);
185 pvect(norm[2]);
186 fvsum(v1, norm[3], vc1, -1.0);
187 pvect(v1);
188 printf("\n%s ", texname);
189 }
190 printf("polygon %s.%d\n", surfname, ++nout);
191 printf("0\n0\n9\n");
192 pvect(p0->p);
193 pvect(p1->p);
194 pvect(p2->p);
195 }
196 if (ok2) {
197 printf("\n%s ", modname);
198 if (axis != -1) {
199 printf("texfunc %s\n", texname);
200 printf(tsargs);
201 printf("0\n13\t%d\n", axis);
202 pvect(norm[0]);
203 pvect(norm[1]);
204 pvect(norm[2]);
205 fvsum(v2, norm[3], vc2, -1.0);
206 pvect(v2);
207 printf("\n%s ", texname);
208 }
209 printf("polygon %s.%d\n", surfname, ++nout);
210 printf("0\n0\n9\n");
211 pvect(p2->p);
212 pvect(p1->p);
213 pvect(p3->p);
214 }
215 }
216
217
218 comprow(s, row, siz) /* compute row of values */
219 double s;
220 register POINT *row;
221 int siz;
222 {
223 double st[2];
224 int end;
225 register int i;
226
227 if (smooth) {
228 i = -1; /* compute one past each end */
229 end = siz+1;
230 } else {
231 if (s < -FTINY || s > 1.0+FTINY)
232 return;
233 i = 0;
234 end = siz;
235 }
236 st[0] = s;
237 while (i <= end) {
238 st[1] = (double)i/siz;
239 row[i].p[0] = funvalue(XNAME, 2, st);
240 row[i].p[1] = funvalue(YNAME, 2, st);
241 row[i].p[2] = funvalue(ZNAME, 2, st);
242 i++;
243 }
244 }
245
246
247 compnorms(r0, r1, r2, siz) /* compute row of averaged normals */
248 register POINT *r0, *r1, *r2;
249 int siz;
250 {
251 FVECT v1, v2, vc;
252 register int i;
253
254 if (!smooth) /* not needed if no smoothing */
255 return;
256 /* compute middle points */
257 while (siz-- >= 0) {
258 fvsum(v1, r2[0].p, r1[0].p, -1.0);
259 fvsum(v2, r1[1].p, r1[0].p, -1.0);
260 fcross(r1[0].n, v1, v2);
261 fvsum(v1, r0[0].p, r1[0].p, -1.0);
262 fcross(vc, v2, v1);
263 fvsum(r1[0].n, r1[0].n, vc, 1.0);
264 fvsum(v2, r1[-1].p, r1[0].p, -1.0);
265 fcross(vc, v1, v2);
266 fvsum(r1[0].n, r1[0].n, vc, 1.0);
267 fvsum(v1, r2[0].p, r1[0].p, -1.0);
268 fcross(vc, v2, v1);
269 fvsum(r1[0].n, r1[0].n, vc, 1.0);
270 normalize(r1[0].n);
271 r0++; r1++; r2++;
272 }
273 }
274
275
276 int
277 norminterp(resmat, p0, p1, p2, p3) /* compute normal interpolation */
278 register FVECT resmat[4];
279 POINT *p0, *p1, *p2, *p3;
280 {
281 #define u ((ax+1)%3)
282 #define v ((ax+2)%3)
283
284 register int ax;
285 double eqnmat[4][4];
286 FVECT v1;
287 register int i, j;
288
289 if (!smooth) /* no interpolation if no smoothing */
290 return(-1);
291 /* find dominant axis */
292 VCOPY(v1, p0->n);
293 fvsum(v1, v1, p1->n, 1.0);
294 fvsum(v1, v1, p2->n, 1.0);
295 fvsum(v1, v1, p3->n, 1.0);
296 ax = ABS(v1[0]) > ABS(v1[1]) ? 0 : 1;
297 ax = ABS(v1[ax]) > ABS(v1[2]) ? ax : 2;
298 /* assign equation matrix */
299 eqnmat[0][0] = p0->p[u]*p0->p[v];
300 eqnmat[0][1] = p0->p[u];
301 eqnmat[0][2] = p0->p[v];
302 eqnmat[0][3] = 1.0;
303 eqnmat[1][0] = p1->p[u]*p1->p[v];
304 eqnmat[1][1] = p1->p[u];
305 eqnmat[1][2] = p1->p[v];
306 eqnmat[1][3] = 1.0;
307 eqnmat[2][0] = p2->p[u]*p2->p[v];
308 eqnmat[2][1] = p2->p[u];
309 eqnmat[2][2] = p2->p[v];
310 eqnmat[2][3] = 1.0;
311 eqnmat[3][0] = p3->p[u]*p3->p[v];
312 eqnmat[3][1] = p3->p[u];
313 eqnmat[3][2] = p3->p[v];
314 eqnmat[3][3] = 1.0;
315 /* invert matrix (solve system) */
316 if (!invmat(eqnmat, eqnmat))
317 return(-1); /* no solution */
318 /* compute result matrix */
319 for (j = 0; j < 4; j++)
320 for (i = 0; i < 3; i++)
321 resmat[j][i] = eqnmat[j][0]*p0->n[i] +
322 eqnmat[j][1]*p1->n[i] +
323 eqnmat[j][2]*p2->n[i] +
324 eqnmat[j][3]*p3->n[i];
325 return(ax);
326
327 #undef u
328 #undef v
329 }
330
331
332 /*
333 * invmat - computes the inverse of mat into inverse. Returns 1
334 * if there exists an inverse, 0 otherwise. It uses Gaussian Elimination
335 * method.
336 */
337
338 invmat(inverse,mat)
339 double mat[4][4],inverse[4][4];
340 {
341 #define SWAP(a,b,t) (t=a,a=b,b=t)
342
343 double m4tmp[4][4];
344 register int i,j,k;
345 register double temp;
346
347 bcopy((char *)mat, (char *)m4tmp, sizeof(m4tmp));
348 /* set inverse to identity */
349 for (i = 0; i < 4; i++)
350 for (j = 0; j < 4; j++)
351 inverse[i][j] = i==j ? 1.0 : 0.0;
352
353 for(i = 0; i < 4; i++) {
354 /* Look for raw with largest pivot and swap raws */
355 temp = FTINY; j = -1;
356 for(k = i; k < 4; k++)
357 if(ABS(m4tmp[k][i]) > temp) {
358 temp = ABS(m4tmp[k][i]);
359 j = k;
360 }
361 if(j == -1) /* No replacing raw -> no inverse */
362 return(0);
363 if (j != i)
364 for(k = 0; k < 4; k++) {
365 SWAP(m4tmp[i][k],m4tmp[j][k],temp);
366 SWAP(inverse[i][k],inverse[j][k],temp);
367 }
368
369 temp = m4tmp[i][i];
370 for(k = 0; k < 4; k++) {
371 m4tmp[i][k] /= temp;
372 inverse[i][k] /= temp;
373 }
374 for(j = 0; j < 4; j++) {
375 if(j != i) {
376 temp = m4tmp[j][i];
377 for(k = 0; k < 4; k++) {
378 m4tmp[j][k] -= m4tmp[i][k]*temp;
379 inverse[j][k] -= inverse[i][k]*temp;
380 }
381 }
382 }
383 }
384 return(1);
385
386 #undef SWAP
387 }
388
389
390 eputs(msg)
391 char *msg;
392 {
393 fputs(msg, stderr);
394 }
395
396
397 wputs(msg)
398 char *msg;
399 {
400 eputs(msg);
401 }
402
403
404 quit(code)
405 {
406 exit(code);
407 }
408
409
410 printhead(ac, av) /* print command header */
411 register int ac;
412 register char **av;
413 {
414 putchar('#');
415 while (ac--) {
416 putchar(' ');
417 fputs(*av++, stdout);
418 }
419 putchar('\n');
420 }
421
422
423 double
424 l_hermite()
425 {
426 double t;
427
428 t = argument(5);
429 return( argument(1)*((2.0*t-3.0)*t*t+1.0) +
430 argument(2)*(-2.0*t+3.0)*t*t +
431 argument(3)*((t-2.0)*t+1.0)*t +
432 argument(4)*(t-1.0)*t*t );
433 }
434
435
436 double
437 l_bezier()
438 {
439 double t;
440
441 t = argument(5);
442 return( argument(1) * (1.+t*(-3.+t*(3.-t))) +
443 argument(2) * 3.*t*(1.+t*(-2.+t)) +
444 argument(3) * 3.*t*t*(1.-t) +
445 argument(4) * t*t*t );
446 }
447
448
449 double
450 l_bspline()
451 {
452 double t;
453
454 t = argument(5);
455 return( argument(1) * (1./6.+t*(-1./2.+t*(1./2.-1./6.*t))) +
456 argument(2) * (2./3.+t*t*(-1.+1./2.*t)) +
457 argument(3) * (1./6.+t*(1./2.+t*(1./2.-1./2.*t))) +
458 argument(4) * (1./6.*t*t*t) );
459 }