ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/gensurf.c
Revision: 1.6
Committed: Fri Mar 2 17:24:02 1990 UTC (34 years, 1 month ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 1.5: +15 -1 lines
Log Message:
Added Bezier cubic function

File Contents

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