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 (19 years, 5 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

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: ranimove2.c,v 3.2 2003/02/25 02:47:24 greg Exp $";
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 <string.h>
15
16 #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 memset((void *)esamp, '\0', sizeof(int)*hres*vres);
42 /*
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 memset((void *)cerrmap, '\0', sizeof(float)*hres*vres);
300 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 }