ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/gensurf.c
Revision: 1.8
Committed: Sat May 26 19:44:40 1990 UTC (33 years, 10 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 1.7: +13 -2 lines
Log Message:
eliminated computations outside [0,1] for unsmoothed surfaces

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