ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/gensurf.c
Revision: 1.4
Committed: Wed Oct 18 18:49:09 1989 UTC (34 years, 5 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 1.3: +51 -79 lines
Log Message:
improved smoothing at edges, improved stability of invmat()

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