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

# Content
1 /* Copyright (c) 1989 Regents of the University of California */
2
3 #ifndef lint
4 static char SCCSid[] = "$SunId$ LBL";
5 #endif
6
7 /*
8 * 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 "standard.h"
20
21 #define XNAME "X_" /* x function name */
22 #define YNAME "Y_" /* y function name */
23 #define ZNAME "Z_" /* z function name */
24
25 #define ABS(x) ((x)>=0 ? (x) : -(x))
26
27 #define pvect(p) printf(vformat, (p)[0], (p)[1], (p)[2])
28
29 char vformat[] = "%15.9g %15.9g %15.9g\n";
30 char tsargs[] = "4 surf_dx surf_dy surf_dz surf.cal\n";
31 char texname[] = "Phong";
32
33 int smooth = 0; /* apply smoothing? */
34
35 char *modname, *surfname;
36
37 double funvalue(), l_hermite(), l_bezier(), argument();
38
39 typedef struct {
40 FVECT p; /* vertex position */
41 FVECT n; /* average normal */
42 } POINT;
43
44
45 main(argc, argv)
46 int argc;
47 char *argv[];
48 {
49 POINT *row0, *row1, *row2, *rp;
50 int i, j, m, n;
51 char stmp[256];
52
53 varset("PI", PI);
54 funset("hermite", 5, l_hermite);
55 funset("bezier", 5, l_bezier);
56
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 else if (!strcmp(argv[i], "-s"))
66 smooth++;
67 else
68 goto userror;
69
70 modname = argv[1];
71 surfname = argv[2];
72 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 row0 = (POINT *)malloc((n+3)*sizeof(POINT));
84 row1 = (POINT *)malloc((n+3)*sizeof(POINT));
85 row2 = (POINT *)malloc((n+3)*sizeof(POINT));
86 if (row0 == NULL || row1 == NULL || row2 == NULL) {
87 fprintf(stderr, "%s: out of memory\n", argv[0]);
88 quit(1);
89 }
90 row0++; row1++; row2++;
91 /* print header */
92 printhead(argc, argv);
93 /* initialize */
94 comprow(-1.0/m, row0, n);
95 comprow(0.0, row1, n);
96 comprow(1.0/m, row2, n);
97 compnorms(row0, row1, row2, n);
98 /* for each row */
99 for (i = 0; i < m; i++) {
100 /* compute next row */
101 rp = row0;
102 row0 = row1;
103 row1 = row2;
104 row2 = rp;
105 comprow((double)(i+2)/m, row2, n);
106 compnorms(row0, row1, row2, n);
107
108 for (j = 0; j < n; j++) {
109 /* 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 }
117 }
118
119 quit(0);
120
121 userror:
122 fprintf(stderr, "Usage: %s material name ", argv[0]);
123 fprintf(stderr, "x(s,t) y(s,t) z(s,t) m n [-s][-e expr][-f file]\n");
124 quit(1);
125 }
126
127
128 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 comprow(s, row, siz) /* compute row of values */
216 double s;
217 register POINT *row;
218 int siz;
219 {
220 double st[2];
221 register int i;
222 /* compute one past each end */
223 st[0] = s;
224 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 }
230 }
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 register int i;
239
240 if (!smooth) /* not needed if no smoothing */
241 return;
242 /* compute middle points */
243 while (siz-- >= 0) {
244 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 double eqnmat[4][4];
272 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 ax = ABS(v1[0]) > ABS(v1[1]) ? 0 : 1;
283 ax = ABS(v1[ax]) > ABS(v1[2]) ? ax : 2;
284 /* 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 if (!invmat(eqnmat, eqnmat))
303 return(-1); /* no solution */
304 /* compute result matrix */
305 for (j = 0; j < 4; j++)
306 for (i = 0; i < 3; i++)
307 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 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 double m4tmp[4][4];
330 register int i,j,k;
331 register double temp;
332
333 bcopy((char *)mat, (char *)m4tmp, sizeof(m4tmp));
334 /* 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
339 for(i = 0; i < 4; i++) {
340 /* 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
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
372 #undef SWAP
373 }
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 }
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 }