ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rsensor.c
Revision: 2.25
Committed: Tue Jun 3 21:31:51 2025 UTC (5 days, 18 hours ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 2.24: +3 -4 lines
Log Message:
refactor: More consistent use of global char * progname and fixargv0()

File Contents

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