ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/gensurf.c
Revision: 1.9
Committed: Mon Jul 9 09:55:47 1990 UTC (33 years, 8 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 1.8: +2 -0 lines
Log Message:
set eclock variable for more efficient computations

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     varset("PI", PI);
55     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     scompile(NULL, argv[++i]);
65     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     scompile(NULL, stmp);
76     sprintf(stmp, "%s(s,t)=%s;", YNAME, argv[4]);
77     scompile(NULL, stmp);
78     sprintf(stmp, "%s(s,t)=%s;", ZNAME, argv[5]);
79     scompile(NULL, stmp);
80     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     FVECT v1, v2, vc;
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.3 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 greg 1.4 double eqnmat[4][4];
286 greg 1.3 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 greg 1.4 ax = ABS(v1[0]) > ABS(v1[1]) ? 0 : 1;
297     ax = ABS(v1[ax]) > ABS(v1[2]) ? ax : 2;
298 greg 1.3 /* 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 greg 1.4 if (!invmat(eqnmat, eqnmat))
317 greg 1.3 return(-1); /* no solution */
318     /* compute result matrix */
319     for (j = 0; j < 4; j++)
320     for (i = 0; i < 3; i++)
321 greg 1.4 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 greg 1.3 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 greg 1.4 double m4tmp[4][4];
344 greg 1.3 register int i,j,k;
345     register double temp;
346    
347 greg 1.5 bcopy((char *)mat, (char *)m4tmp, sizeof(m4tmp));
348 greg 1.4 /* 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 greg 1.3
353     for(i = 0; i < 4; i++) {
354 greg 1.4 /* 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 greg 1.3
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 greg 1.4
386 greg 1.3 #undef SWAP
387 greg 1.1 }
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 greg 1.6 }
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 greg 1.7 }
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 greg 1.1 }