ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/gensurf.c
Revision: 1.13
Committed: Tue Apr 23 13:04:58 1991 UTC (32 years, 11 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 1.12: +1 -1 lines
Log Message:
changed arguments to varset()

File Contents

# User Rev Content
1 greg 1.1 #ifndef lint
2     static char SCCSid[] = "$SunId$ LBL";
3     #endif
4 greg 1.2
5 greg 1.7 /* Copyright (c) 1989 Regents of the University of California */
6    
7 greg 1.2 /*
8 greg 1.1 * 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 greg 1.5 #include "standard.h"
20 greg 1.1
21     #define XNAME "X_" /* x function name */
22     #define YNAME "Y_" /* y function name */
23     #define ZNAME "Z_" /* z function name */
24    
25 greg 1.4 #define ABS(x) ((x)>=0 ? (x) : -(x))
26    
27 greg 1.3 #define pvect(p) printf(vformat, (p)[0], (p)[1], (p)[2])
28 greg 1.1
29     char vformat[] = "%15.9g %15.9g %15.9g\n";
30 greg 1.3 char tsargs[] = "4 surf_dx surf_dy surf_dz surf.cal\n";
31     char texname[] = "Phong";
32 greg 1.1
33 greg 1.3 int smooth = 0; /* apply smoothing? */
34 greg 1.1
35 greg 1.3 char *modname, *surfname;
36 greg 1.1
37 greg 1.7 double funvalue(), l_hermite(), l_bezier(), l_bspline(), argument();
38 greg 1.3
39     typedef struct {
40     FVECT p; /* vertex position */
41     FVECT n; /* average normal */
42     } POINT;
43    
44    
45 greg 1.1 main(argc, argv)
46     int argc;
47     char *argv[];
48     {
49 greg 1.9 extern long eclock;
50 greg 1.3 POINT *row0, *row1, *row2, *rp;
51 greg 1.1 int i, j, m, n;
52     char stmp[256];
53    
54 greg 1.13 varset("PI", ':', PI);
55 greg 1.1 funset("hermite", 5, l_hermite);
56 greg 1.6 funset("bezier", 5, l_bezier);
57 greg 1.7 funset("bspline", 5, l_bspline);
58 greg 1.1
59     if (argc < 8)
60     goto userror;
61    
62     for (i = 8; i < argc; i++)
63     if (!strcmp(argv[i], "-e"))
64 greg 1.10 scompile(argv[++i], NULL, 0);
65 greg 1.1 else if (!strcmp(argv[i], "-f"))
66     fcompile(argv[++i]);
67 greg 1.3 else if (!strcmp(argv[i], "-s"))
68     smooth++;
69 greg 1.1 else
70     goto userror;
71    
72 greg 1.3 modname = argv[1];
73     surfname = argv[2];
74 greg 1.1 sprintf(stmp, "%s(s,t)=%s;", XNAME, argv[3]);
75 greg 1.10 scompile(stmp, NULL, 0);
76 greg 1.1 sprintf(stmp, "%s(s,t)=%s;", YNAME, argv[4]);
77 greg 1.10 scompile(stmp, NULL, 0);
78 greg 1.1 sprintf(stmp, "%s(s,t)=%s;", ZNAME, argv[5]);
79 greg 1.10 scompile(stmp, NULL, 0);
80 greg 1.1 m = atoi(argv[6]);
81     n = atoi(argv[7]);
82     if (m <= 0 || n <= 0)
83     goto userror;
84    
85 greg 1.4 row0 = (POINT *)malloc((n+3)*sizeof(POINT));
86     row1 = (POINT *)malloc((n+3)*sizeof(POINT));
87     row2 = (POINT *)malloc((n+3)*sizeof(POINT));
88 greg 1.3 if (row0 == NULL || row1 == NULL || row2 == NULL) {
89 greg 1.1 fprintf(stderr, "%s: out of memory\n", argv[0]);
90     quit(1);
91     }
92 greg 1.4 row0++; row1++; row2++;
93 greg 1.3 /* print header */
94 greg 1.1 printhead(argc, argv);
95 greg 1.9 eclock = 0;
96 greg 1.4 /* initialize */
97     comprow(-1.0/m, row0, n);
98 greg 1.3 comprow(0.0, row1, n);
99     comprow(1.0/m, row2, n);
100 greg 1.4 compnorms(row0, row1, row2, n);
101 greg 1.3 /* for each row */
102 greg 1.1 for (i = 0; i < m; i++) {
103     /* compute next row */
104 greg 1.3 rp = row0;
105 greg 1.1 row0 = row1;
106 greg 1.3 row1 = row2;
107     row2 = rp;
108 greg 1.4 comprow((double)(i+2)/m, row2, n);
109     compnorms(row0, row1, row2, n);
110 greg 1.1
111     for (j = 0; j < n; j++) {
112 greg 1.3 /* 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 greg 1.1 }
120     }
121    
122     quit(0);
123    
124     userror:
125     fprintf(stderr, "Usage: %s material name ", argv[0]);
126 greg 1.3 fprintf(stderr, "x(s,t) y(s,t) z(s,t) m n [-s][-e expr][-f file]\n");
127 greg 1.1 quit(1);
128     }
129    
130    
131 greg 1.3 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 greg 1.1 comprow(s, row, siz) /* compute row of values */
219     double s;
220 greg 1.3 register POINT *row;
221 greg 1.1 int siz;
222     {
223 greg 1.4 double st[2];
224 greg 1.8 int end;
225 greg 1.4 register int i;
226 greg 1.8
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 greg 1.1 st[0] = s;
237 greg 1.8 while (i <= end) {
238 greg 1.4 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 greg 1.8 i++;
243 greg 1.1 }
244 greg 1.3 }
245    
246    
247     compnorms(r0, r1, r2, siz) /* compute row of averaged normals */
248     register POINT *r0, *r1, *r2;
249     int siz;
250     {
251 greg 1.11 FVECT v1, v2;
252 greg 1.4 register int i;
253 greg 1.3
254     if (!smooth) /* not needed if no smoothing */
255     return;
256     /* compute middle points */
257 greg 1.4 while (siz-- >= 0) {
258 greg 1.11 fvsum(v1, r2[0].p, r0[0].p, -1.0);
259     fvsum(v2, r1[1].p, r1[-1].p, -1.0);
260 greg 1.3 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 greg 1.12 MAT4 eqnmat;
277 greg 1.3 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 greg 1.4 ax = ABS(v1[0]) > ABS(v1[1]) ? 0 : 1;
288     ax = ABS(v1[ax]) > ABS(v1[2]) ? ax : 2;
289 greg 1.3 /* 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 greg 1.4 if (!invmat(eqnmat, eqnmat))
308 greg 1.3 return(-1); /* no solution */
309     /* compute result matrix */
310     for (j = 0; j < 4; j++)
311     for (i = 0; i < 3; i++)
312 greg 1.4 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 greg 1.3 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 greg 1.12 MAT4 inverse, mat;
331 greg 1.3 {
332     #define SWAP(a,b,t) (t=a,a=b,b=t)
333    
334 greg 1.12 MAT4 m4tmp;
335 greg 1.3 register int i,j,k;
336     register double temp;
337    
338 greg 1.12 copymat4(m4tmp, mat);
339 greg 1.4 /* 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 greg 1.3
344     for(i = 0; i < 4; i++) {
345 greg 1.11 /* Look for row with largest pivot and swap rows */
346 greg 1.4 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 greg 1.11 if(j == -1) /* No replacing row -> no inverse */
353 greg 1.4 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 greg 1.3
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 greg 1.4
377 greg 1.3 #undef SWAP
378 greg 1.1 }
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 greg 1.6 }
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 greg 1.7 }
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 greg 1.1 }