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

# Content
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 #include "copyright.h"
13
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 }