ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/ranimove2.c
Revision: 3.2
Committed: Tue Feb 25 02:47:24 2003 UTC (21 years ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad3R5
Changes since 3.1: +1 -56 lines
Log Message:
Replaced inline copyright notice with #include "copyright.h"

File Contents

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