ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/gensurf.c
Revision: 1.14
Committed: Tue Apr 23 15:51:11 1991 UTC (32 years, 11 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 1.13: +3 -3 lines
Log Message:
changed parameters to funset() call

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;
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, r0[0].p, -1.0);
259 fvsum(v2, r1[1].p, r1[-1].p, -1.0);
260 fcross(r1[0].n, v1, v2);
261 normalize(r1[0].n);
262 r0++; r1++; r2++;
263 }
264 }
265
266
267 int
268 norminterp(resmat, p0, p1, p2, p3) /* compute normal interpolation */
269 register FVECT resmat[4];
270 POINT *p0, *p1, *p2, *p3;
271 {
272 #define u ((ax+1)%3)
273 #define v ((ax+2)%3)
274
275 register int ax;
276 MAT4 eqnmat;
277 FVECT v1;
278 register int i, j;
279
280 if (!smooth) /* no interpolation if no smoothing */
281 return(-1);
282 /* find dominant axis */
283 VCOPY(v1, p0->n);
284 fvsum(v1, v1, p1->n, 1.0);
285 fvsum(v1, v1, p2->n, 1.0);
286 fvsum(v1, v1, p3->n, 1.0);
287 ax = ABS(v1[0]) > ABS(v1[1]) ? 0 : 1;
288 ax = ABS(v1[ax]) > ABS(v1[2]) ? ax : 2;
289 /* assign equation matrix */
290 eqnmat[0][0] = p0->p[u]*p0->p[v];
291 eqnmat[0][1] = p0->p[u];
292 eqnmat[0][2] = p0->p[v];
293 eqnmat[0][3] = 1.0;
294 eqnmat[1][0] = p1->p[u]*p1->p[v];
295 eqnmat[1][1] = p1->p[u];
296 eqnmat[1][2] = p1->p[v];
297 eqnmat[1][3] = 1.0;
298 eqnmat[2][0] = p2->p[u]*p2->p[v];
299 eqnmat[2][1] = p2->p[u];
300 eqnmat[2][2] = p2->p[v];
301 eqnmat[2][3] = 1.0;
302 eqnmat[3][0] = p3->p[u]*p3->p[v];
303 eqnmat[3][1] = p3->p[u];
304 eqnmat[3][2] = p3->p[v];
305 eqnmat[3][3] = 1.0;
306 /* invert matrix (solve system) */
307 if (!invmat(eqnmat, eqnmat))
308 return(-1); /* no solution */
309 /* compute result matrix */
310 for (j = 0; j < 4; j++)
311 for (i = 0; i < 3; i++)
312 resmat[j][i] = eqnmat[j][0]*p0->n[i] +
313 eqnmat[j][1]*p1->n[i] +
314 eqnmat[j][2]*p2->n[i] +
315 eqnmat[j][3]*p3->n[i];
316 return(ax);
317
318 #undef u
319 #undef v
320 }
321
322
323 /*
324 * invmat - computes the inverse of mat into inverse. Returns 1
325 * if there exists an inverse, 0 otherwise. It uses Gaussian Elimination
326 * method.
327 */
328
329 invmat(inverse,mat)
330 MAT4 inverse, mat;
331 {
332 #define SWAP(a,b,t) (t=a,a=b,b=t)
333
334 MAT4 m4tmp;
335 register int i,j,k;
336 register double temp;
337
338 copymat4(m4tmp, mat);
339 /* set inverse to identity */
340 for (i = 0; i < 4; i++)
341 for (j = 0; j < 4; j++)
342 inverse[i][j] = i==j ? 1.0 : 0.0;
343
344 for(i = 0; i < 4; i++) {
345 /* Look for row with largest pivot and swap rows */
346 temp = FTINY; j = -1;
347 for(k = i; k < 4; k++)
348 if(ABS(m4tmp[k][i]) > temp) {
349 temp = ABS(m4tmp[k][i]);
350 j = k;
351 }
352 if(j == -1) /* No replacing row -> no inverse */
353 return(0);
354 if (j != i)
355 for(k = 0; k < 4; k++) {
356 SWAP(m4tmp[i][k],m4tmp[j][k],temp);
357 SWAP(inverse[i][k],inverse[j][k],temp);
358 }
359
360 temp = m4tmp[i][i];
361 for(k = 0; k < 4; k++) {
362 m4tmp[i][k] /= temp;
363 inverse[i][k] /= temp;
364 }
365 for(j = 0; j < 4; j++) {
366 if(j != i) {
367 temp = m4tmp[j][i];
368 for(k = 0; k < 4; k++) {
369 m4tmp[j][k] -= m4tmp[i][k]*temp;
370 inverse[j][k] -= inverse[i][k]*temp;
371 }
372 }
373 }
374 }
375 return(1);
376
377 #undef SWAP
378 }
379
380
381 eputs(msg)
382 char *msg;
383 {
384 fputs(msg, stderr);
385 }
386
387
388 wputs(msg)
389 char *msg;
390 {
391 eputs(msg);
392 }
393
394
395 quit(code)
396 {
397 exit(code);
398 }
399
400
401 printhead(ac, av) /* print command header */
402 register int ac;
403 register char **av;
404 {
405 putchar('#');
406 while (ac--) {
407 putchar(' ');
408 fputs(*av++, stdout);
409 }
410 putchar('\n');
411 }
412
413
414 double
415 l_hermite()
416 {
417 double t;
418
419 t = argument(5);
420 return( argument(1)*((2.0*t-3.0)*t*t+1.0) +
421 argument(2)*(-2.0*t+3.0)*t*t +
422 argument(3)*((t-2.0)*t+1.0)*t +
423 argument(4)*(t-1.0)*t*t );
424 }
425
426
427 double
428 l_bezier()
429 {
430 double t;
431
432 t = argument(5);
433 return( argument(1) * (1.+t*(-3.+t*(3.-t))) +
434 argument(2) * 3.*t*(1.+t*(-2.+t)) +
435 argument(3) * 3.*t*t*(1.-t) +
436 argument(4) * t*t*t );
437 }
438
439
440 double
441 l_bspline()
442 {
443 double t;
444
445 t = argument(5);
446 return( argument(1) * (1./6.+t*(-1./2.+t*(1./2.-1./6.*t))) +
447 argument(2) * (2./3.+t*t*(-1.+1./2.*t)) +
448 argument(3) * (1./6.+t*(1./2.+t*(1./2.-1./2.*t))) +
449 argument(4) * (1./6.*t*t*t) );
450 }