ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rsensor.c
Revision: 2.10
Committed: Thu May 19 16:09:38 2011 UTC (12 years, 10 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.9: +14 -4 lines
Log Message:
Bug fix related to zero rows and added check for negative sensitivities

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.10 static const char RCSid[] = "$Id: rsensor.c,v 2.9 2011/05/17 19:34:36 greg Exp $";
3 greg 2.1 #endif
4    
5     /*
6     * Compute sensor signal based on spatial sensitivity.
7     *
8     * Created Feb 2008 for Architectural Energy Corp.
9     */
10    
11     #include "ray.h"
12     #include "source.h"
13     #include "view.h"
14     #include "random.h"
15    
16     #define DEGREE (PI/180.)
17    
18     #define MAXNT 180 /* maximum number of theta divisions */
19     #define MAXNP 360 /* maximum number of phi divisions */
20    
21     extern char *progname; /* global argv[0] */
22     extern int nowarn; /* don't report warnings? */
23    
24     /* current sensor's perspective */
25 greg 2.5 VIEW ourview = {VT_ANG,{0.,0.,0.},{0.,0.,1.},{1.,0.,0.},
26     1.,180.,180.,0.,0.,0.,0.,
27     {0.,0.,0.},{0.,0.,0.},0.,0.};
28 greg 2.1
29     unsigned long nsamps = 10000; /* desired number of initial samples */
30     unsigned long nssamps = 9000; /* number of super-samples */
31 greg 2.2 int ndsamps = 32; /* number of direct samples */
32 greg 2.1 int nprocs = 1; /* number of rendering processes */
33    
34     float *sensor = NULL; /* current sensor data */
35     int sntp[2]; /* number of sensor theta and phi angles */
36     float maxtheta; /* maximum theta value for this sensor */
37 greg 2.9 float tvals[MAXNT+1]; /* theta prob. values (1-D table of 1-cos(t)) */
38     float *pvals = NULL; /* phi prob. values (2-D table in radians) */
39 greg 2.1 int ntheta = 0; /* polar angle divisions */
40     int nphi = 0; /* azimuthal angle divisions */
41     double gscale = 1.; /* global scaling value */
42    
43 greg 2.2 #define s_theta(t) sensor[(t+1)*(sntp[1]+1)]
44     #define s_phi(p) sensor[(p)+1]
45     #define s_val(t,p) sensor[(p)+1+(t+1)*(sntp[1]+1)]
46    
47 greg 2.1 static void comp_sensor(char *sfile);
48    
49     static void
50 greg 2.3 over_options() /* overriding options */
51 greg 2.1 {
52 greg 2.3 directvis = (ndsamps <= 0);
53     do_irrad = 0;
54     }
55    
56     static void
57     print_defaults() /* print out default parameters */
58     {
59     over_options();
60 greg 2.1 printf("-n %-9d\t\t\t# number of processes\n", nprocs);
61     printf("-rd %-9ld\t\t\t# ray directions\n", nsamps);
62     /* printf("-rs %-9ld\t\t\t# ray super-samples\n", nssamps); */
63     printf("-dn %-9d\t\t\t# direct number of samples\n", ndsamps);
64     printf("-vp %f %f %f\t# view point\n",
65     ourview.vp[0], ourview.vp[1], ourview.vp[2]);
66     printf("-vd %f %f %f\t# view direction\n",
67     ourview.vdir[0], ourview.vdir[1], ourview.vdir[2]);
68     printf("-vu %f %f %f\t# view up\n",
69     ourview.vup[0], ourview.vup[1], ourview.vup[2]);
70     printf("-vo %f\t\t\t# view fore clipping distance\n", ourview.vfore);
71     print_rdefaults();
72     }
73    
74 greg 2.7
75     void
76     quit(ec) /* make sure exit is called */
77     int ec;
78     {
79     if (ray_pnprocs > 0) /* close children if any */
80     ray_pclose(0);
81     exit(ec);
82     }
83    
84    
85 greg 2.1 int
86     main(
87     int argc,
88     char *argv[]
89     )
90     {
91     int doheader = 1;
92 greg 2.3 int optwarn = 0;
93 greg 2.1 int i, rval;
94    
95     progname = argv[0];
96     /* set up rendering defaults */
97 greg 2.2 rand_samp = 1;
98 greg 2.6 dstrsrc = 0.65;
99 greg 2.2 srcsizerat = 0.1;
100 greg 2.1 directrelay = 3;
101     ambounce = 1;
102 greg 2.2 maxdepth = -10;
103 greg 2.1 /* get options from command line */
104 greg 2.3 for (i = 1; i < argc; i++) {
105 greg 2.1 while ((rval = expandarg(&argc, &argv, i)) > 0)
106     ;
107     if (rval < 0) {
108     sprintf(errmsg, "cannot expand '%s'", argv[i]);
109     error(SYSTEM, errmsg);
110     }
111 greg 2.3 if (argv[i][0] != '-') {
112     if (i >= argc-1)
113     break; /* final octree argument */
114 greg 2.1 if (!ray_pnprocs) {
115 greg 2.3 over_options();
116 greg 2.1 if (doheader) { /* print header */
117 greg 2.6 newheader("RADIANCE", stdout);
118 greg 2.1 printargs(argc, argv, stdout);
119     fputformat("ascii", stdout);
120     putchar('\n');
121     }
122     /* start process(es) */
123 greg 2.8 if (strcmp(argv[argc-1], "."))
124     ray_pinit(argv[argc-1], nprocs);
125 greg 2.1 }
126 greg 2.3 comp_sensor(argv[i]); /* process a sensor file */
127 greg 2.1 continue;
128     }
129     if (argv[i][1] == 'r') { /* sampling options */
130     if (argv[i][2] == 'd')
131     nsamps = atol(argv[++i]);
132     else if (argv[i][2] == 's')
133     nssamps = atol(argv[++i]);
134     else {
135     sprintf(errmsg, "bad option at '%s'", argv[i]);
136     error(USER, errmsg);
137     }
138     continue;
139     }
140     /* direct component samples */
141     if (argv[i][1] == 'd' && argv[i][2] == 'n') {
142     ndsamps = atoi(argv[++i]);
143     continue;
144     }
145     if (argv[i][1] == 'v') { /* next sensor view */
146     if (argv[i][2] == 'f') {
147     rval = viewfile(argv[++i], &ourview, NULL);
148     if (rval < 0) {
149     sprintf(errmsg,
150     "cannot open view file \"%s\"",
151     argv[i]);
152     error(SYSTEM, errmsg);
153     } else if (rval == 0) {
154     sprintf(errmsg,
155     "bad view file \"%s\"",
156     argv[i]);
157     error(USER, errmsg);
158     }
159     continue;
160     }
161     rval = getviewopt(&ourview, argc-i, argv+i);
162     if (rval >= 0) {
163     i += rval;
164     continue;
165     }
166     sprintf(errmsg, "bad view option at '%s'", argv[i]);
167     error(USER, errmsg);
168     }
169 greg 2.3 if (!strcmp(argv[i], "-w")) { /* toggle warnings */
170     nowarn = !nowarn;
171 greg 2.1 continue;
172     }
173     if (ray_pnprocs) {
174 greg 2.3 if (!optwarn++)
175     error(WARNING,
176 greg 2.1 "rendering options should appear before first sensor");
177     } else if (!strcmp(argv[i], "-defaults")) {
178     print_defaults();
179     return(0);
180     }
181     if (argv[i][1] == 'h') { /* header toggle */
182     doheader = !doheader;
183     continue;
184     }
185     if (!strcmp(argv[i], "-n")) { /* number of processes */
186     nprocs = atoi(argv[++i]);
187     if (nprocs <= 0)
188     error(USER, "illegal number of processes");
189     continue;
190     }
191     rval = getrenderopt(argc-i, argv+i);
192     if (rval < 0) {
193     sprintf(errmsg, "bad render option at '%s'", argv[i]);
194     error(USER, errmsg);
195     }
196     i += rval;
197     }
198 greg 2.8 if (sensor == NULL)
199 greg 2.3 error(USER, i<argc ? "missing sensor file" : "missing octree");
200 greg 2.1 quit(0);
201     }
202    
203     /* Load sensor sensitivities (first row and column are angles) */
204     static float *
205     load_sensor(
206     int ntp[2],
207     char *sfile
208     )
209     {
210 greg 2.10 int warnedneg;
211 greg 2.1 char linebuf[8192];
212     int nelem = 1000;
213     float *sarr = (float *)malloc(sizeof(float)*nelem);
214     FILE *fp;
215     char *cp;
216     int i;
217    
218     fp = frlibopen(sfile);
219     if (fp == NULL) {
220     sprintf(errmsg, "cannot open sensor file '%s'", sfile);
221     error(SYSTEM, errmsg);
222     }
223     fgets(linebuf, sizeof(linebuf), fp);
224     if (!strncmp(linebuf, "Elevation ", 10))
225     fgets(linebuf, sizeof(linebuf), fp);
226     /* get phi values */
227     sarr[0] = .0f;
228     if (strncmp(linebuf, "degrees", 7)) {
229     sprintf(errmsg, "Missing 'degrees' in sensor file '%s'", sfile);
230     error(USER, errmsg);
231     }
232     cp = sskip(linebuf);
233     ntp[1] = 0;
234     for ( ; ; ) {
235     sarr[ntp[1]+1] = atof(cp);
236     cp = fskip(cp);
237     if (cp == NULL)
238     break;
239 greg 2.9 if (ntp[1] > 1 && sarr[ntp[1]+1] <= sarr[ntp[1]]) {
240     sprintf(errmsg,
241     "Phi values not monotinically increasing in sensor file '%s'",
242     sfile);
243     error(USER, errmsg);
244     }
245 greg 2.1 ++ntp[1];
246     }
247 greg 2.10 warnedneg = 0;
248 greg 2.1 ntp[0] = 0; /* get thetas + data */
249     while (fgets(linebuf, sizeof(linebuf), fp) != NULL) {
250     ++ntp[0];
251     if ((ntp[0]+1)*(ntp[1]+1) > nelem) {
252     nelem += (nelem>>2) + ntp[1];
253     sarr = (float *)realloc((void *)sarr,
254     sizeof(float)*nelem);
255     if (sarr == NULL)
256     error(SYSTEM, "out of memory in load_sensor()");
257     }
258     cp = linebuf;
259     i = ntp[0]*(ntp[1]+1);
260     for ( ; ; ) {
261     sarr[i] = atof(cp);
262     cp = fskip(cp);
263     if (cp == NULL)
264     break;
265 greg 2.10 if (i && sarr[i] < .0) {
266     if (!warnedneg++) {
267     sprintf(errmsg,
268     "Negative value(s) in sensor file '%s' (ignored)\n", sfile);
269     error(WARNING, errmsg);
270     }
271     sarr[i] = .0;
272     }
273 greg 2.1 ++i;
274     }
275     if (i == ntp[0]*(ntp[1]+1))
276     break;
277 greg 2.9 if (ntp[0] > 1 && sarr[ntp[0]*(ntp[1]+1)] <=
278     sarr[(ntp[0]-1)*(ntp[1]+1)]) {
279     sprintf(errmsg,
280     "Theta values not monotinically increasing in sensor file '%s'",
281     sfile);
282     error(USER, errmsg);
283     }
284 greg 2.1 if (i != (ntp[0]+1)*(ntp[1]+1)) {
285     sprintf(errmsg,
286     "bad column count near line %d in sensor file '%s'",
287     ntp[0]+1, sfile);
288     error(USER, errmsg);
289     }
290     }
291     nelem = i;
292     fclose(fp);
293     errmsg[0] = '\0'; /* sanity checks */
294     if (ntp[0] <= 0)
295     sprintf(errmsg, "no data in sensor file '%s'", sfile);
296     else if (fabs(sarr[ntp[1]+1]) > FTINY)
297     sprintf(errmsg, "minimum theta must be 0 in sensor file '%s'",
298     sfile);
299     else if (fabs(sarr[1]) > FTINY)
300     sprintf(errmsg, "minimum phi must be 0 in sensor file '%s'",
301     sfile);
302     else if (sarr[ntp[1]] <= FTINY)
303     sprintf(errmsg,
304     "maximum phi must be positive in sensor file '%s'",
305     sfile);
306     else if (sarr[ntp[0]*(ntp[1]+1)] <= FTINY)
307     sprintf(errmsg,
308     "maximum theta must be positive in sensor file '%s'",
309     sfile);
310     if (errmsg[0])
311     error(USER, errmsg);
312     return((float *)realloc((void *)sarr, sizeof(float)*nelem));
313     }
314    
315     /* Initialize probability table */
316     static void
317     init_ptable(
318     char *sfile
319     )
320     {
321     int samptot = nsamps;
322     float *rowp, *rowp1;
323     double rowsum[MAXNT], rowomega[MAXNT];
324     double thdiv[MAXNT+1], phdiv[MAXNP+1];
325     double tsize, psize;
326     double prob, frac, frac1;
327     int i, j, t, p;
328     /* free old table */
329     if (sensor != NULL)
330     free((void *)sensor);
331     if (pvals != NULL)
332     free((void *)pvals);
333     if (sfile == NULL || !*sfile) {
334 greg 2.2 sensor = NULL;
335     sntp[0] = sntp[1] = 0;
336 greg 2.1 pvals = NULL;
337     ntheta = nphi = 0;
338     return;
339     }
340     /* load sensor table */
341     sensor = load_sensor(sntp, sfile);
342     if (sntp[0] > MAXNT) {
343     sprintf(errmsg, "Too many theta rows in sensor file '%s'",
344     sfile);
345     error(INTERNAL, errmsg);
346     }
347     if (sntp[1] > MAXNP) {
348     sprintf(errmsg, "Too many phi columns in sensor file '%s'",
349     sfile);
350     error(INTERNAL, errmsg);
351     }
352     /* compute boundary angles */
353 greg 2.9 maxtheta = DEGREE*(1.5f*s_theta(sntp[0]-1) - 0.5f*s_theta(sntp[0]-2));
354     if (maxtheta > PI)
355     maxtheta = PI;
356 greg 2.1 thdiv[0] = .0;
357     for (t = 1; t < sntp[0]; t++)
358 greg 2.2 thdiv[t] = DEGREE/2.*(s_theta(t-1) + s_theta(t));
359 greg 2.9 thdiv[sntp[0]] = maxtheta;
360     phdiv[0] = DEGREE*(1.5f*s_phi(0) - 0.5f*s_phi(1));
361 greg 2.1 for (p = 1; p < sntp[1]; p++)
362 greg 2.2 phdiv[p] = DEGREE/2.*(s_phi(p-1) + s_phi(p));
363 greg 2.9 phdiv[sntp[1]] = DEGREE*(1.5f*s_phi(sntp[1]-1) - 0.5f*s_phi(sntp[1]-2));
364 greg 2.1 /* size our table */
365 greg 2.9 tsize = 1. - cos(maxtheta);
366     psize = PI*tsize/maxtheta;
367 greg 2.1 if (sntp[0]*sntp[1] < samptot) /* don't overdo resolution */
368     samptot = sntp[0]*sntp[1];
369 greg 2.2 ntheta = (int)(sqrt((double)samptot*tsize/psize) + 0.5);
370 greg 2.1 if (ntheta > MAXNT)
371     ntheta = MAXNT;
372     nphi = samptot/ntheta;
373 greg 2.9 pvals = (float *)malloc(sizeof(float)*(ntheta+1)*(nphi+1));
374 greg 2.1 if (pvals == NULL)
375     error(SYSTEM, "out of memory in init_ptable()");
376     gscale = .0; /* compute our inverse table */
377     for (i = 0; i < sntp[0]; i++) {
378 greg 2.2 rowp = &s_val(i,0);
379 greg 2.9 rowsum[i] = 1e-20;
380 greg 2.1 for (j = 0; j < sntp[1]; j++)
381     rowsum[i] += *rowp++;
382     rowomega[i] = cos(thdiv[i]) - cos(thdiv[i+1]);
383     rowomega[i] *= 2.*PI / (double)sntp[1];
384     gscale += rowsum[i] * rowomega[i];
385     }
386 greg 2.9 if (gscale <= FTINY) {
387     sprintf(errmsg, "Sensor values sum to zero in file '%s'", sfile);
388     error(USER, errmsg);
389     }
390 greg 2.2 for (i = 0; i < ntheta; i++) {
391 greg 2.1 prob = (double)i / (double)ntheta;
392     for (t = 0; t < sntp[0]; t++)
393     if ((prob -= rowsum[t]*rowomega[t]/gscale) <= .0)
394     break;
395     if (t >= sntp[0])
396     error(INTERNAL, "code error 1 in init_ptable()");
397     frac = 1. + prob/(rowsum[t]*rowomega[t]/gscale);
398     tvals[i] = 1. - ( (1.-frac)*cos(thdiv[t]) +
399     frac*cos(thdiv[t+1]) );
400 greg 2.2 /* offset b/c sensor values are centered */
401 greg 2.10 if (!t || (t < sntp[0]-1) & (frac >= 0.5))
402 greg 2.2 frac -= 0.5;
403 greg 2.10 else {
404 greg 2.2 frac += 0.5;
405     --t;
406     }
407 greg 2.9 pvals[i*(nphi+1)] = phdiv[0];
408 greg 2.1 for (j = 1; j < nphi; j++) {
409     prob = (double)j / (double)nphi;
410 greg 2.2 rowp = &s_val(t,0);
411     rowp1 = &s_val(t+1,0);
412 greg 2.9 for (p = 0; p < sntp[1]; p++)
413 greg 2.2 if ((prob -= (1.-frac)*rowp[p]/rowsum[t] +
414     frac*rowp1[p]/rowsum[t+1]) <= .0)
415 greg 2.1 break;
416 greg 2.10 if (p >= sntp[1]) { /* should never happen? */
417 greg 2.9 p = sntp[1] - 1;
418     prob = .5;
419 greg 2.1 }
420 greg 2.9 frac1 = 1. + prob/((1.-frac)*rowp[p]/rowsum[t]
421     + frac*rowp1[p]/rowsum[t+1]);
422     pvals[i*(nphi+1) + j] = (1.-frac1)*phdiv[p] +
423     frac1*phdiv[p+1];
424 greg 2.1 }
425 greg 2.9 pvals[i*(nphi+1) + nphi] = phdiv[sntp[1]];
426 greg 2.1 }
427 greg 2.2 tvals[0] = .0f;
428 greg 2.1 tvals[ntheta] = (float)tsize;
429     }
430    
431     /* Get normalized direction from random variables in [0,1) range */
432     static void
433     get_direc(
434     FVECT dvec,
435     double x,
436     double y
437     )
438     {
439     double xfrac = x*ntheta;
440     int tndx = (int)xfrac;
441     double yfrac = y*nphi;
442     int pndx = (int)yfrac;
443     double rad, phi;
444     FVECT dv;
445     int i;
446    
447     xfrac -= (double)tndx;
448     yfrac -= (double)pndx;
449     pndx += tndx*(nphi+1);
450    
451     dv[2] = 1. - ((1.-xfrac)*tvals[tndx] + xfrac*tvals[tndx+1]);
452     rad = sqrt(1. - dv[2]*dv[2]);
453     phi = (1.-yfrac)*pvals[pndx] + yfrac*pvals[pndx+1];
454     dv[0] = -rad*sin(phi);
455     dv[1] = rad*cos(phi);
456     for (i = 3; i--; )
457     dvec[i] = dv[0]*ourview.hvec[i] +
458     dv[1]*ourview.vvec[i] +
459     dv[2]*ourview.vdir[i] ;
460     }
461    
462     /* Get sensor value in the specified direction (normalized) */
463     static float
464     sens_val(
465     FVECT dvec
466     )
467     {
468     FVECT dv;
469     float theta, phi;
470     int t, p;
471    
472     dv[2] = DOT(dvec, ourview.vdir);
473     theta = (float)((1./DEGREE) * acos(dv[2]));
474     if (theta >= maxtheta)
475     return(.0f);
476     dv[0] = DOT(dvec, ourview.hvec);
477     dv[1] = DOT(dvec, ourview.vvec);
478     phi = (float)((1./DEGREE) * atan2(-dv[0], dv[1]));
479     while (phi < .0f) phi += 360.f;
480     t = (int)(theta/maxtheta * sntp[0]);
481     p = (int)(phi*(1./360.) * sntp[1]);
482     /* hack for non-uniform sensor grid */
483 greg 2.2 while (t+1 < sntp[0] && theta >= s_theta(t+1))
484 greg 2.1 ++t;
485 greg 2.2 while (t-1 >= 0 && theta <= s_theta(t-1))
486 greg 2.1 --t;
487 greg 2.2 while (p+1 < sntp[1] && phi >= s_phi(p+1))
488 greg 2.1 ++p;
489 greg 2.2 while (p-1 >= 0 && phi <= s_phi(p-1))
490 greg 2.1 --p;
491 greg 2.2 return(s_val(t,p));
492 greg 2.1 }
493    
494 greg 2.8 /* Print origin and direction */
495     static void
496     print_ray(
497     FVECT rorg,
498     FVECT rdir
499     )
500     {
501     printf("%.6g %.6g %.6g %.8f %.8f %.8f\n",
502     rorg[0], rorg[1], rorg[2],
503     rdir[0], rdir[1], rdir[2]);
504     }
505    
506 greg 2.1 /* Compute sensor output */
507     static void
508     comp_sensor(
509     char *sfile
510     )
511     {
512     int ndirs = dstrsrc > FTINY ? ndsamps :
513     ndsamps > 0 ? 1 : 0;
514     char *err;
515     int nt, np;
516     COLOR vsum;
517     RAY rr;
518 greg 2.6 double sf;
519 greg 2.1 int i, j;
520     /* set view */
521     ourview.type = VT_ANG;
522     ourview.horiz = ourview.vert = 180.;
523     ourview.hoff = ourview.voff = .0;
524     err = setview(&ourview);
525     if (err != NULL)
526     error(USER, err);
527     /* assign probability table */
528     init_ptable(sfile);
529 greg 2.6 /* stratified MC sampling */
530 greg 2.1 setcolor(vsum, .0f, .0f, .0f);
531     nt = (int)(sqrt((double)nsamps*ntheta/nphi) + .5);
532     np = nsamps/nt;
533 greg 2.6 sf = gscale/nsamps;
534 greg 2.1 for (i = 0; i < nt; i++)
535 greg 2.2 for (j = 0; j < np; j++) {
536 greg 2.5 VCOPY(rr.rorg, ourview.vp);
537 greg 2.2 get_direc(rr.rdir, (i+frandom())/nt, (j+frandom())/np);
538 greg 2.5 if (ourview.vfore > FTINY)
539     VSUM(rr.rorg, rr.rorg, rr.rdir, ourview.vfore);
540 greg 2.8 if (!ray_pnprocs) {
541     print_ray(rr.rorg, rr.rdir);
542     continue;
543     }
544 greg 2.5 rr.rmax = .0;
545 greg 2.1 rayorigin(&rr, PRIMARY, NULL, NULL);
546 greg 2.6 scalecolor(rr.rcoef, sf);
547 greg 2.1 if (ray_pqueue(&rr) == 1)
548     addcolor(vsum, rr.rcol);
549     }
550 greg 2.6 /* remaining rays pure MC */
551     for (i = nsamps - nt*np; i-- > 0; ) {
552     VCOPY(rr.rorg, ourview.vp);
553     get_direc(rr.rdir, frandom(), frandom());
554     if (ourview.vfore > FTINY)
555     VSUM(rr.rorg, rr.rorg, rr.rdir, ourview.vfore);
556 greg 2.8 if (!ray_pnprocs) {
557     print_ray(rr.rorg, rr.rdir);
558     continue;
559     }
560 greg 2.6 rr.rmax = .0;
561     rayorigin(&rr, PRIMARY, NULL, NULL);
562     scalecolor(rr.rcoef, sf);
563     if (ray_pqueue(&rr) == 1)
564     addcolor(vsum, rr.rcol);
565     }
566 greg 2.8 if (!ray_pnprocs) /* just printing rays */
567     return;
568 greg 2.6 /* scale partial result */
569     scalecolor(vsum, sf);
570     /* add direct component */
571 greg 2.1 for (i = ndirs; i-- > 0; ) {
572     SRCINDEX si;
573     initsrcindex(&si);
574     while (srcray(&rr, NULL, &si)) {
575 greg 2.6 sf = sens_val(rr.rdir);
576     if (sf <= FTINY)
577 greg 2.1 continue;
578 greg 2.6 sf *= si.dom/ndirs;
579     scalecolor(rr.rcoef, sf);
580 greg 2.1 if (ray_pqueue(&rr) == 1) {
581     multcolor(rr.rcol, rr.rcoef);
582     addcolor(vsum, rr.rcol);
583     }
584     }
585     }
586 greg 2.6 /* finish our calculation */
587 greg 2.1 while (ray_presult(&rr, 0) > 0) {
588     multcolor(rr.rcol, rr.rcoef);
589     addcolor(vsum, rr.rcol);
590     }
591     /* print our result */
592     printf("%.4e %.4e %.4e\n", colval(vsum,RED),
593     colval(vsum,GRN), colval(vsum,BLU));
594     }