ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/ranimove2.c
Revision: 3.3
Committed: Mon Jun 30 14:59:13 2003 UTC (20 years, 8 months ago) by schorsch
Content type: text/plain
Branch: MAIN
Changes since 3.2: +5 -3 lines
Log Message:
Replaced most outdated BSD function calls with their posix equivalents, and cleaned up a few other platform dependencies.

File Contents

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