ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/ranimove2.c
Revision: 3.11
Committed: Sun Mar 29 02:55:01 2020 UTC (2 years, 6 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad5R3, HEAD
Changes since 3.10: +3 -2 lines
Log Message:
Fixed seg fault caused by reading past end of array

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: ranimove2.c,v 3.10 2016/04/18 22:39:13 greg Exp $";
3 #endif
4 /*
5 * ranimove2.c
6 *
7 * Frame refinement routines for ranimove(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 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
28
29 int
30 refine_first(void) /* initial refinement pass */
31 {
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 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 memset((void *)esamp, '\0', sizeof(int)*hres*vres);
46 /*
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 if ((xp < 0) | (xp >= hres))
59 continue;
60 yp = y + ymbuffer[n];
61 if ((yp < 0) | (yp >= vres))
62 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 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 static void
130 clr_consp( /* initialize a conspicuity sum */
131 struct ConspSum *cs
132 )
133 {
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 sum_consp( /* sum in conspicuity result */
146 struct ConspSum *cdest,
147 struct ConspSum *cs
148 )
149 {
150 if ((cdest == NULL) | (cs == NULL))
151 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 est_consp( /* estimate error conspicuity & update */
164 int x0,
165 int y0,
166 int x1,
167 int y1,
168 struct ConspSum *cs
169 )
170 {
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 if ((x0 == x1-1) & (y0 == y1-1)) { /* update pixel error */
185 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 /* put into map */
237 for ( ; y0 < y1; y0++) {
238 float *em0 = cerrmap + fndx(x0, y0);
239 float *emp = em0 + (x1-x0);
240 while (emp-- > em0)
241 if (eest > *emp)
242 *emp = eest;
243 }
244 cerrzero = 0;
245 }
246
247 static void
248 subconspicuity( /* compute subportion of conspicuity */
249 int x0,
250 int y0,
251 int x1,
252 int y1,
253 struct ConspSum *cs
254 )
255 {
256 struct ConspSum mysum;
257 int i;
258
259 if ((x0 >= x1) | (y0 >= y1))
260 error(CONSISTENCY, "bad call to subconspicuity");
261
262 clr_consp(&mysum); /* prepare sum */
263
264 if ((x0 == x1-1) & (y0 == y1-1)) { /* single pixel */
265 double hls;
266 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 void
314 conspicuity(void) /* compute conspicuous error map */
315 {
316 int fhres, fvres;
317 int fx, fy;
318 /* reuse previous z-buffer */
319 cerrmap = (float *)zprev;
320 memset((void *)cerrmap, '\0', sizeof(float)*hres*vres);
321 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 ppri_cmp( /* pixel priority comparison */
355 const void *pp1,
356 const void *pp2
357 )
358 {
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 ray_refine( /* refine the given pixel by tracing a ray */
378 int n
379 )
380 {
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 rayorigin(&ir, PRIMARY, NULL, NULL);
400 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 (colval(ctmp,RED) > 0.01) &
415 (colval(ctmp,GRN) > 0.01) &
416 (colval(ctmp,BLU) > 0.01)) {
417 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 refine_rays( /* compute refinement rays */
445 long nrays
446 )
447 {
448 int *pord;
449 int ntodo;
450 long rdone;
451 int i;
452 /* skip if nothing significant */
453 if (ndtset && cerrzero)
454 return(0);
455 /* 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) & (i < ntodo) &&
467 cerrmap[pord[i]] <= FTINY)
468 ntodo = i;
469 if (i >= ntodo) { /* redo conspicuity & priority */
470 while (ray_refine(-1) >= 0)
471 ;
472 conspicuity();
473 if (ndtset && cerrzero)
474 break;
475 qsort((void *)pord, hres*vres, sizeof(int), ppri_cmp);
476 ntodo = hres*vres/8;
477 i = 0;
478 }
479 /* sample next pixel */
480 ray_refine(pord[i++]);
481 }
482 /* clean up and return */
483 while (ray_refine(-1) >= 0)
484 ;
485 free((void *)pord);
486 return(rdone);
487 }
488
489
490 int
491 refine_frame( /* refine current frame */
492 int pass
493 )
494 {
495 static double rtime_used = 0;
496 static long ray_cnt = 0;
497 static double ctime_used = 0;
498 static int csp_cnt = 0;
499 int timed = (fcur > fbeg) | (pass > 0) | (quickstart);
500 double time_start, rtime_start, time_done;
501 struct AmbSum myAmbSum;
502 long rays_todo, nr;
503 int n;
504 /* IBR refinement? */
505 if ((pass == 0) & (fcur > fbeg))
506 return(refine_first());
507 /* any time left? */
508 time_start = getTime();
509 if (timed) {
510 if (time_start >= frm_stop)
511 goto nomore;
512 if (csp_cnt > 0 && time_start + ctime_used/csp_cnt >= frm_stop)
513 goto nomore;
514 }
515 asump = NULL; /* use resampling to update ambval? */
516 if (!curparams->ambounce && hirendparams.ambounce) {
517 myAmbSum.diffsum[RED] =
518 myAmbSum.diffsum[GRN] =
519 myAmbSum.diffsum[BLU] = 0;
520 myAmbSum.nsamps = 0;
521 asump = &myAmbSum;
522 }
523 /* initialize value-squared map */
524 if (val2map == NULL) {
525 val2map = cprev; /* OK to reuse at this point */
526 n = (asump == NULL) ? hres*vres : 0;
527 while (n--)
528 if (sbuffer[n])
529 setcolor(val2map[n],
530 colval(cbuffer[n],RED)*colval(cbuffer[n],RED),
531 colval(cbuffer[n],GRN)*colval(cbuffer[n],GRN),
532 colval(cbuffer[n],BLU)*colval(cbuffer[n],BLU));
533 else
534 setcolor(val2map[n], 0., 0., 0.);
535 }
536 /* compute conspicuity */
537 if (!silent) {
538 printf("\tComputing conspicuity map\n");
539 fflush(stdout);
540 }
541 conspicuity();
542 csp_cnt++;
543 #if 0
544 if (pass == 1) {
545 char fnm[256];
546 sprintf(fnm, vval(BASENAME), fcur);
547 strcat(fnm, "_incmap.pic");
548 write_map(cerrmap, fnm);
549 }
550 #endif
551 /* get ray start time */
552 rtime_start = getTime();
553 ctime_used += rtime_start - time_start;
554 if (timed && rtime_start >= frm_stop)
555 return(0); /* error done but out of time */
556 if (rtime_used <= FTINY) {
557 if (quickstart)
558 rays_todo = 1000;
559 else
560 rays_todo = hres*vres;
561 } else {
562 rays_todo = (long)((frm_stop - rtime_start) *
563 ray_cnt / rtime_used);
564 if (rays_todo < 1000)
565 return(0); /* let's call it a frame */
566 }
567 /* set higher rendering quality */
568 if (twolevels && curparams != &hirendparams) {
569 ray_restore(curparams = &hirendparams);
570 if (nprocs > 1) { /* need to update children */
571 if (!silent) {
572 printf("\tRestarting %d processes\n", nprocs);
573 fflush(stdout);
574 }
575 ray_pclose(0);
576 ray_popen(nprocs);
577 }
578 }
579 /* compute refinement rays */
580 if (!silent) {
581 printf("\tRefinement pass %d...",
582 pass+1); /*, rays_todo); */
583 fflush(stdout);
584 }
585 if (asump != NULL) /* flag low-quality samples */
586 for (n = hres*vres; n--; )
587 if (sbuffer[n])
588 abuffer[n] = ALOWQ;
589 /* trace those rays */
590 nr = refine_rays(rays_todo);
591 if (!silent)
592 printf("traced %ld HQ rays\n", nr);
593 if (nr <= 0)
594 return(0);
595 /* update timing stats */
596 while (ray_cnt >= 1L<<20) {
597 ray_cnt >>= 1;
598 rtime_used *= .5;
599 }
600 ray_cnt += nr;
601 time_done = getTime();
602 rtime_used += time_done - rtime_start;
603 if (!timed && time_done > frm_stop)
604 frm_stop = time_done;
605 /* update ambient value */
606 if (asump != NULL && asump->nsamps >= 1000) {
607 double sf = 1./(double)asump->nsamps;
608 for (n = 3; n--; ) {
609 asump->diffsum[n] *= sf;
610 asump->diffsum[n] += colval(lorendparams.ambval,n);
611 if (asump->diffsum[n] < 0) asump->diffsum[n] = 0;
612 }
613 setcolor(lorendparams.ambval,
614 asump->diffsum[RED],
615 asump->diffsum[GRN],
616 asump->diffsum[BLU]);
617 if (!silent)
618 printf("\tUpdated parameter: -av %f %f %f\n",
619 asump->diffsum[RED],
620 asump->diffsum[GRN],
621 asump->diffsum[BLU]);
622 asump = NULL;
623 }
624 return(1);
625 nomore:
626 /* make sure error map is updated */
627 if ((fcur == fbeg) | (pass > 1))
628 comp_frame_error();
629 return(0);
630 }