ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/ranimove2.c
Revision: 3.8
Committed: Tue Apr 19 01:15:07 2005 UTC (18 years, 11 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 3.7: +2 -2 lines
Log Message:
Extensive changes to enable rtrace -oTW option for tracking ray contributions

File Contents

# User Rev Content
1 greg 3.1 #ifndef lint
2 greg 3.8 static const char RCSid[] = "$Id: ranimove2.c,v 3.7 2004/06/16 09:50:33 greg Exp $";
3 greg 3.1 #endif
4     /*
5     * ranimove2.c
6     *
7     * Frame refinement routines for ranimate(1).
8     *
9     * Created by Gregory Ward on Wed Jan 08 2003.
10     */
11    
12 greg 3.2 #include "copyright.h"
13 greg 3.1
14 schorsch 3.3 #include <string.h>
15    
16 greg 3.1 #include "ranimove.h"
17     #include "random.h"
18    
19    
20     #define HL_ERR 0.32 /* highlight error threshold */
21    
22     int cerrzero; /* is cerrmap all zeroes? */
23    
24 schorsch 3.6 static int ppri_cmp(const void *pp1, const void *pp2);
25     static int ray_refine(int n);
26     static long refine_rays(long nrays);
27 greg 3.1
28 schorsch 3.6
29     extern int
30     refine_first(void) /* initial refinement pass */
31 greg 3.1 {
32     int *esamp = (int *)zprev; /* OK to reuse */
33     int hl_erri = errori(HL_ERR);
34     int nextra = 0;
35     int x, y, xp, yp;
36     int neigh;
37     register int n, np;
38    
39     if (sizeof(int) < sizeof(*zprev))
40     error(CONSISTENCY, "code error in refine_first");
41     if (!silent) {
42     printf("\tFirst refinement pass...");
43     fflush(stdout);
44     }
45 schorsch 3.3 memset((void *)esamp, '\0', sizeof(int)*hres*vres);
46 greg 3.1 /*
47     * In our initial pass, we look for lower error pixels from
48     * the same objects in the previous frame, and copy them here.
49     */
50     for (y = vres; y--; )
51     for (x = hres; x--; ) {
52     n = fndx(x, y);
53     if (obuffer[n] == OVOID)
54     continue;
55     if (xmbuffer[n] == MO_UNK)
56     continue;
57     xp = x + xmbuffer[n];
58 schorsch 3.5 if ((xp < 0) | (xp >= hres))
59 greg 3.1 continue;
60     yp = y + ymbuffer[n];
61 schorsch 3.5 if ((yp < 0) | (yp >= vres))
62 greg 3.1 continue;
63     np = fndx(xp, yp);
64     /* make sure we hit same object */
65     if (oprev[np] != obuffer[n])
66     continue;
67     /* is previous frame error lower? */
68     if (aprev[np] < AMIN + ATIDIFF)
69     continue;
70     if (aprev[np] <= abuffer[n] + ATIDIFF)
71     continue;
72     /* shadow & highlight detection */
73     if (abuffer[n] > hl_erri &&
74     getclosest(&neigh, 1, x, y) &&
75     bigdiff(cbuffer[neigh], cprev[np],
76     HL_ERR*(.9+.2*frandom())))
77     continue;
78     abuffer[n] = aprev[np] - ATIDIFF;
79     copycolor(cbuffer[n], cprev[np]);
80     esamp[n] = 1; /* record extrapolated sample */
81     nextra++;
82     }
83     for (n = hres*vres; n--; ) /* update sample counts */
84     if (esamp[n])
85     sbuffer[n] = 1;
86     if (!silent)
87     printf("extrapolated %d pixels\n", nextra);
88     return(1);
89     }
90    
91    
92     /*
93     * We use a recursive computation of the conspicuity
94     * map to avoid associated memory costs and simplify
95     * coding. We create a virtual image pyramid, pooling
96     * variance calculations, etc. The top of the pyramid
97     * corresponds to the foveal resolution, as there should
98     * not be any interesting mechanisms above this level.
99     */
100    
101     #define CSF_C0 1.14
102     #define CSF_C1 0.67
103     #define CSF_C2 1.7
104     #define CSF_S1 6.1
105     #define CSF_S2 7.3
106     #define CSF_P1 45.9
107     #define CSF_PC (30./45.9*CSF_P1)
108     #define CSF_VR0 0.15
109     #define CSF_VRC 80.
110    
111     struct ConspSum {
112     COLOR vsum; /* value sum */
113     COLOR v2sum; /* value^2 sum */
114     long nsamp; /* number of samples */
115     long xmsum; /* x-motion sum */
116     long ymsum; /* y-motion sum */
117     int npix; /* number of pixels */
118     double hls; /* high-level saliency */
119     };
120    
121     static double pixel_deg; /* base pixel frequency */
122     static int fhsiz, fvsiz; /* foveal subimage size */
123    
124 schorsch 3.6 static void clr_consp(struct ConspSum *cs);
125     static void sum_consp(struct ConspSum *cdest, struct ConspSum *cs);
126     static void est_consp(int x0, int y0, int x1, int y1, struct ConspSum *cs);
127     static void subconspicuity(int x0, int y0, int x1, int y1, struct ConspSum *cs);
128    
129 greg 3.1 static void
130 schorsch 3.6 clr_consp( /* initialize a conspicuity sum */
131     register struct ConspSum *cs
132     )
133 greg 3.1 {
134     if (cs == NULL)
135     return;
136     setcolor(cs->vsum, 0., 0., 0.);
137     setcolor(cs->v2sum, 0., 0., 0.);
138     cs->nsamp = 0;
139     cs->xmsum = cs->ymsum = 0;
140     cs->npix = 0;
141     cs->hls = 0;
142     }
143    
144     static void
145 schorsch 3.6 sum_consp( /* sum in conspicuity result */
146     register struct ConspSum *cdest,
147     register struct ConspSum *cs
148     )
149 greg 3.1 {
150 schorsch 3.5 if ((cdest == NULL) | (cs == NULL))
151 greg 3.1 return;
152     addcolor(cdest->vsum, cs->vsum);
153     addcolor(cdest->v2sum, cs->v2sum);
154     cdest->nsamp += cs->nsamp;
155     cdest->xmsum += cs->xmsum;
156     cdest->ymsum += cs->ymsum;
157     cdest->npix += cs->npix;
158     if (cs->hls > cdest->hls)
159     cdest->hls = cs->hls;
160     }
161    
162     static void
163 schorsch 3.6 est_consp( /* estimate error conspicuity & update */
164     int x0,
165     int y0,
166     int x1,
167     int y1,
168     register struct ConspSum *cs
169     )
170 greg 3.1 {
171     double rad2, mtn2, cpd, vm, vr, csf, eest;
172     /* do we care? */
173     if (cs->hls <= FTINY)
174     return;
175     /* get relative error */
176     if (cs->nsamp < NSAMPOK) {
177     int neigh[NSAMPOK]; /* gather neighbors */
178     eest = comperr(neigh,
179     getclosest(neigh, NSAMPOK, (x0+x1)>>1, (y0+y1)>>1),
180     cs->nsamp);
181     } else
182     eest = estimaterr(cs->vsum, cs->v2sum, cs->nsamp, cs->nsamp);
183    
184 schorsch 3.5 if ((x0 == x1-1) & (y0 == y1-1)) { /* update pixel error */
185 greg 3.1 int n = fndx(x0, y0);
186     int ai;
187     int ne;
188     if (sbuffer[n] >= 255) {
189     abuffer[n] = ADISTANT;
190     } else {
191     ai = errori(eest);
192     if (ai < AMIN) ai = AMIN;
193     else if (ai >= ADISTANT) ai = ADISTANT-1;
194     abuffer[n] = ai;
195     /* can't improve on closest */
196     if (!cs->nsamp && getclosest(&ne, 1, x0, y0) &&
197     abuffer[ne] < ai &&
198     abuffer[ne] >= AMIN)
199     abuffer[n] = abuffer[ne];
200     }
201     }
202     /* compute radius^2 */
203     rad2 = 0.125*((x1-x0)*(x1-x0) + (y1-y0)*(y1-y0));
204    
205     /* average motion^2 */
206     mtn2 = (double)cs->xmsum*cs->xmsum + (double)cs->ymsum*cs->ymsum;
207     mtn2 /= (double)(cs->npix*cs->npix);
208     /* motion blur hides us? */
209     if (mblur*mblur*mtn2 >= 4.*rad2)
210     return;
211     /* too small to see? */
212     cpd = pixel_deg * pixel_deg / rad2;
213     if (cpd > CSF_PC*CSF_PC)
214     return;
215     cpd = sqrt(cpd);
216     /* compute CSF [Daley98] */
217     vm = rate * sqrt(mtn2) / pixel_deg;
218     vr = cs->hls/hlsmax*vm + CSF_VR0; /* use hls tracking eff. */
219     if (vr > CSF_VRC) vr = CSF_VRC;
220     vr = vm - vr;
221     if (vr < 0) vr = -vr;
222     csf = log(CSF_C2*(1./3.)*vr);
223     if (csf < 0) csf = -csf;
224     csf = CSF_S1 + CSF_S2*csf*csf*csf;
225     csf *= CSF_C0*CSF_C2*4.*PI*PI*CSF_C1*CSF_C1*cpd*cpd;
226     csf *= exp(-CSF_C1*4.*PI/CSF_P1*(CSF_C2*vr + 2.)*cpd);
227     /* compute visible error */
228     eest = eest*csf/ndthresh - 1.;
229     if (eest <= FTINY)
230     return;
231     /* scale by saleincy */
232     eest *= cs->hls;
233     /* worth the bother? */
234     if (eest <= .01)
235     return;
236 greg 3.7 /* put into map */
237 greg 3.1 for ( ; y0 < y1; y0++) {
238     float *em0 = cerrmap + fndx(x0, y0);
239     register float *emp = em0 + (x1-x0);
240     while (emp-- > em0)
241 greg 3.7 if (eest > *emp)
242     *emp = eest;
243 greg 3.1 }
244     cerrzero = 0;
245     }
246    
247     static void
248 schorsch 3.6 subconspicuity( /* compute subportion of conspicuity */
249     int x0,
250     int y0,
251     int x1,
252     int y1,
253     struct ConspSum *cs
254     )
255 greg 3.1 {
256     struct ConspSum mysum;
257     int i;
258    
259 schorsch 3.5 if ((x0 >= x1) | (y0 >= y1))
260 greg 3.1 error(CONSISTENCY, "bad call to subconspicuity");
261    
262     clr_consp(&mysum); /* prepare sum */
263    
264 schorsch 3.5 if ((x0 == x1-1) & (y0 == y1-1)) { /* single pixel */
265 greg 3.1 double hls;
266     register int n = fndx(x0, y0);
267     if (sbuffer[n]) {
268     copycolor(mysum.vsum, cbuffer[n]);
269     copycolor(mysum.v2sum, val2map[n]);
270     mysum.nsamp = sbuffer[n];
271     }
272     if ((mysum.xmsum = xmbuffer[n]) == MO_UNK)
273     mysum.xmsum = 0;
274     else
275     mysum.ymsum = ymbuffer[n];
276     mysum.npix = 1;
277     /* max. hls in fovea */
278     mysum.hls = obj_prio(obuffer[n]);
279     if (x0 >= fhsiz) {
280     hls = obj_prio(obuffer[fndx(x0-fhsiz,y0)]);
281     if (hls > mysum.hls) mysum.hls = hls;
282     }
283     if (x0 < hres-fhsiz) {
284     hls = obj_prio(obuffer[fndx(x0+fhsiz,y0)]);
285     if (hls > mysum.hls) mysum.hls = hls;
286     }
287     if (y0 >= fvsiz) {
288     hls = obj_prio(obuffer[fndx(x0,y0-fvsiz)]);
289     if (hls > mysum.hls) mysum.hls = hls;
290     }
291     if (y0 < vres-fvsiz) {
292     hls = obj_prio(obuffer[fndx(x0,y0+fvsiz)]);
293     if (hls > mysum.hls) mysum.hls = hls;
294     }
295     } else if (x0 == x1-1) { /* vertical pair */
296     for (i = y0 ; i < y1; i++)
297     subconspicuity(x0, i, x1, i+1, &mysum);
298     } else if (y0 == y1-1) { /* horizontal pair */
299     for (i = x0 ; i < x1; i++)
300     subconspicuity(i, y0, i+1, y1, &mysum);
301     } else { /* rectangle */
302     subconspicuity(x0, y0, (x0+x1)>>1, (y0+y1)>>1, &mysum);
303     subconspicuity((x0+x1)>>1, y0, x1, (y0+y1)>>1, &mysum);
304     subconspicuity(x0, (y0+y1)>>1, (x0+x1)>>1, y1, &mysum);
305     subconspicuity((x0+x1)>>1, (y0+y1)>>1, x1, y1, &mysum);
306     }
307     /* update conspicuity */
308     est_consp(x0, y0, x1, y1, &mysum);
309     /* sum into return value */
310     sum_consp(cs, &mysum);
311     }
312    
313 schorsch 3.6 extern void
314     conspicuity(void) /* compute conspicuous error map */
315 greg 3.1 {
316     int fhres, fvres;
317     int fx, fy;
318     /* reuse previous z-buffer */
319     cerrmap = (float *)zprev;
320 schorsch 3.3 memset((void *)cerrmap, '\0', sizeof(float)*hres*vres);
321 greg 3.1 cerrzero = 1;
322     /* compute base pixel frequency */
323     pixel_deg = .5*(hres/vw.horiz + vres/vw.vert);
324     /* compute foveal resolution */
325     fhres = vw.horiz/FOV_DEG + 0.5;
326     if (fhres <= 0) fhres = 1;
327     else if (fhres > hres) fhres = hres;
328     fvres = vw.vert/FOV_DEG + 0.5;
329     if (fvres <= 0) fvres = 1;
330     else if (fvres > vres) fvres = vres;
331     fhsiz = hres/fhres;
332     fvsiz = vres/fvres;
333     /* call our foveal subroutine */
334     for (fy = fvres; fy--; )
335     for (fx = fhres; fx--; )
336     subconspicuity(hres*fx/fhres, vres*fy/fvres,
337     hres*(fx+1)/fhres, vres*(fy+1)/fvres,
338     NULL);
339     }
340    
341    
342     /*
343     * The following structure is used to collect data on the
344     * initial error in the ambient value estimate, in order
345     * to correct for it in the subsequent frames.
346     */
347     static struct AmbSum {
348     double diffsum[3]; /* sum of (correct - ambval) */
349     long nsamps; /* number of values in sum */
350     } *asump = NULL;
351    
352    
353     static int
354 schorsch 3.6 ppri_cmp( /* pixel priority comparison */
355     const void *pp1,
356     const void *pp2
357     )
358 greg 3.1 {
359     double se1 = cerrmap[*(const int *)pp1];
360     double se2 = cerrmap[*(const int *)pp2];
361     int adiff;
362     /* higher conspicuity to front */
363     if (se1 < se2) return(1);
364     if (se1 > se2) return(-1);
365     /* else higher error to front */
366     adiff = (int)abuffer[*(const int *)pp1] -
367     (int)abuffer[*(const int *)pp2];
368     if (adiff)
369     return(adiff);
370     /* else fewer samples to front */
371     return((int)sbuffer[*(const int *)pp1] -
372     (int)sbuffer[*(const int *)pp2]);
373     }
374    
375    
376     static int
377 schorsch 3.6 ray_refine( /* refine the given pixel by tracing a ray */
378     register int n
379     )
380 greg 3.1 {
381     RAY ir;
382     COLOR ctmp;
383     int i;
384    
385     if (n < 0) { /* fetching stragglers */
386     if (nprocs <= 1 || !ray_presult(&ir, 0))
387     return(-1);
388     n = ir.rno;
389     } else { /* else tracing a new ray */
390     double hv[2];
391     if (sbuffer[n] >= 255) /* reached limit? */
392     return(-1);
393     sample_pos(hv, n%hres, n/hres, sbuffer[n]);
394     ir.rmax = viewray(ir.rorg, ir.rdir, &vw, hv[0], hv[1]);
395     if (ir.rmax < -FTINY)
396     return(-1);
397     if (nprocs > 1) {
398     int rval;
399 greg 3.8 rayorigin(&ir, PRIMARY, NULL, NULL);
400 greg 3.1 ir.rno = n;
401     rval = ray_pqueue(&ir);
402     if (!rval)
403     return(-1);
404     if (rval < 0)
405     quit(1);
406     n = ir.rno;
407     } else
408     ray_trace(&ir);
409     }
410     if (abuffer[n] == ALOWQ && asump != NULL) {
411     if (sbuffer[n] != 1)
412     error(CONSISTENCY, "bad code in ray_refine");
413     if (getambcolor(ctmp, obuffer[n]) &&
414 schorsch 3.5 (colval(ctmp,RED) > 0.01) &
415     (colval(ctmp,GRN) > 0.01) &
416     (colval(ctmp,BLU) > 0.01)) {
417 greg 3.1 for (i = 0; i < 3; i++)
418     asump->diffsum[i] +=
419     (colval(ir.rcol,i) - colval(cbuffer[n],i))
420     / colval(ctmp,i);
421     asump->nsamps++;
422     }
423     sbuffer[n] = 0;
424     }
425     setcolor(ctmp,
426     colval(ir.rcol,RED)*colval(ir.rcol,RED),
427     colval(ir.rcol,GRN)*colval(ir.rcol,GRN),
428     colval(ir.rcol,BLU)*colval(ir.rcol,BLU));
429     if (!sbuffer[n]) { /* first sample */
430     copycolor(cbuffer[n], ir.rcol);
431     copycolor(val2map[n], ctmp);
432     abuffer[n] = AHIGHQ;
433     sbuffer[n] = 1;
434     } else { /* else sum in sample */
435     addcolor(cbuffer[n], ir.rcol);
436     addcolor(val2map[n], ctmp);
437     sbuffer[n]++;
438     }
439     return(n);
440     }
441    
442    
443     static long
444 schorsch 3.6 refine_rays( /* compute refinement rays */
445     long nrays
446     )
447 greg 3.1 {
448     int *pord;
449     int ntodo;
450     long rdone;
451     int i;
452     /* skip if nothing significant */
453     if (ndtset && cerrzero)
454 greg 3.4 return(0);
455 greg 3.1 /* initialize priority list */
456     pord = (int *)malloc(sizeof(int)*hres*vres);
457     for (i = hres*vres; i--; )
458     pord[i] = i;
459     /* sort our priorities */
460     ntodo = hres*vres;
461     if (nrays < ntodo)
462     qsort((void *)pord, hres*vres, sizeof(int), ppri_cmp);
463     i = 0;
464     /* trace rays in list */
465     for (rdone = 0; rdone < nrays; rdone++) {
466     if (ndtset && i >= 1000 && cerrmap[pord[i]] <= FTINY)
467     ntodo = i;
468     if (i >= ntodo) { /* redo conspicuity & priority */
469     while (ray_refine(-1) >= 0)
470     ;
471     conspicuity();
472     if (ndtset && cerrzero)
473     break;
474     qsort((void *)pord, hres*vres, sizeof(int), ppri_cmp);
475     ntodo = hres*vres/8;
476     i = 0;
477     }
478     /* sample next pixel */
479     ray_refine(pord[i++]);
480     }
481     /* clean up and return */
482     while (ray_refine(-1) >= 0)
483     ;
484     free((void *)pord);
485     return(rdone);
486     }
487    
488    
489 schorsch 3.6 extern int
490     refine_frame( /* refine current frame */
491     int pass
492     )
493 greg 3.1 {
494     static double rtime_used = 0;
495     static long ray_cnt = 0;
496     static double ctime_used = 0;
497     static int csp_cnt = 0;
498 schorsch 3.5 int timed = (fcur > fbeg) | (pass > 0) | (quickstart);
499 greg 3.1 double time_start, rtime_start, time_done;
500     struct AmbSum myAmbSum;
501     long rays_todo, nr;
502     register int n;
503     /* IBR refinement? */
504 schorsch 3.5 if ((pass == 0) & (fcur > fbeg))
505 greg 3.1 return(refine_first());
506     /* any time left? */
507     time_start = getTime();
508     if (timed) {
509     if (time_start >= frm_stop)
510     goto nomore;
511     if (csp_cnt > 0 && time_start + ctime_used/csp_cnt >= frm_stop)
512     goto nomore;
513     }
514     asump = NULL; /* use resampling to update ambval? */
515     if (!curparams->ambounce && hirendparams.ambounce) {
516     myAmbSum.diffsum[RED] =
517     myAmbSum.diffsum[GRN] =
518     myAmbSum.diffsum[BLU] = 0;
519     myAmbSum.nsamps = 0;
520     asump = &myAmbSum;
521     }
522     /* initialize value-squared map */
523     if (val2map == NULL) {
524     val2map = cprev; /* OK to reuse at this point */
525     n = (asump == NULL) ? hres*vres : 0;
526     while (n--)
527     if (sbuffer[n])
528     setcolor(val2map[n],
529     colval(cbuffer[n],RED)*colval(cbuffer[n],RED),
530     colval(cbuffer[n],GRN)*colval(cbuffer[n],GRN),
531     colval(cbuffer[n],BLU)*colval(cbuffer[n],BLU));
532     else
533     setcolor(val2map[n], 0., 0., 0.);
534     }
535     /* compute conspicuity */
536     if (!silent) {
537     printf("\tComputing conspicuity map\n");
538     fflush(stdout);
539     }
540     conspicuity();
541     csp_cnt++;
542     #if 0
543     if (pass == 1) {
544     char fnm[256];
545     sprintf(fnm, vval(BASENAME), fcur);
546     strcat(fnm, "_incmap.pic");
547     write_map(cerrmap, fnm);
548     }
549     #endif
550     /* get ray start time */
551     rtime_start = getTime();
552     ctime_used += rtime_start - time_start;
553     if (timed && rtime_start >= frm_stop)
554     return(0); /* error done but out of time */
555     if (rtime_used <= FTINY) {
556     if (quickstart)
557     rays_todo = 1000;
558     else
559     rays_todo = hres*vres;
560     } else {
561     rays_todo = (long)((frm_stop - rtime_start) *
562     ray_cnt / rtime_used);
563     if (rays_todo < 1000)
564     return(0); /* let's call it a frame */
565     }
566     /* set higher rendering quality */
567     if (twolevels && curparams != &hirendparams) {
568     ray_restore(curparams = &hirendparams);
569     if (nprocs > 1) { /* need to update children */
570     if (!silent) {
571     printf("\tRestarting %d processes\n", nprocs);
572     fflush(stdout);
573     }
574     ray_pclose(0);
575     ray_popen(nprocs);
576     }
577     }
578     /* compute refinement rays */
579     if (!silent) {
580     printf("\tRefinement pass %d...",
581 schorsch 3.6 pass+1); /*, rays_todo); */
582 greg 3.1 fflush(stdout);
583     }
584     if (asump != NULL) /* flag low-quality samples */
585     for (n = hres*vres; n--; )
586     if (sbuffer[n])
587     abuffer[n] = ALOWQ;
588     /* trace those rays */
589     nr = refine_rays(rays_todo);
590     if (!silent)
591 schorsch 3.6 printf("traced %ld HQ rays\n", nr);
592 greg 3.1 if (nr <= 0)
593     return(0);
594     /* update timing stats */
595     while (ray_cnt >= 1L<<20) {
596     ray_cnt >>= 1;
597     rtime_used *= .5;
598     }
599     ray_cnt += nr;
600     time_done = getTime();
601     rtime_used += time_done - rtime_start;
602     if (!timed && time_done > frm_stop)
603     frm_stop = time_done;
604     /* update ambient value */
605     if (asump != NULL && asump->nsamps >= 1000) {
606     double sf = 1./(double)asump->nsamps;
607     for (n = 3; n--; ) {
608     asump->diffsum[n] *= sf;
609     asump->diffsum[n] += colval(lorendparams.ambval,n);
610     if (asump->diffsum[n] < 0) asump->diffsum[n] = 0;
611     }
612     setcolor(lorendparams.ambval,
613     asump->diffsum[RED],
614     asump->diffsum[GRN],
615     asump->diffsum[BLU]);
616     if (!silent)
617     printf("\tUpdated parameter: -av %f %f %f\n",
618     asump->diffsum[RED],
619     asump->diffsum[GRN],
620     asump->diffsum[BLU]);
621     asump = NULL;
622     }
623     return(1);
624     nomore:
625     /* make sure error map is updated */
626 schorsch 3.5 if ((fcur == fbeg) | (pass > 1))
627 greg 3.1 comp_frame_error();
628     return(0);
629     }