src/util/ranimove2.c

#ifndef lint
static const char       RCSid[] = "$Id: ranimove2.c,v 3.8 2005/04/19 01:15:07 greg Exp $";
#endif
/*
 *  ranimove2.c
 *
 *  Frame refinement routines for ranimove(1).
 *
 *  Created by Gregory Ward on Wed Jan 08 2003.
 */

#include "copyright.h"

#include <string.h>

#include "ranimove.h"
#include "random.h"


#define HL_ERR  0.32            /* highlight error threshold */

int     cerrzero;               /* is cerrmap all zeroes? */

static int ppri_cmp(const void *pp1, const void *pp2);
static int ray_refine(int       n);
static long refine_rays(long    nrays);


extern int
refine_first(void)                      /* initial refinement pass */
{
        int     *esamp = (int *)zprev;  /* OK to reuse */
        int     hl_erri = errori(HL_ERR);
        int     nextra = 0;
        int     x, y, xp, yp;
        int     neigh;
        register int    n, np;

        if (sizeof(int) < sizeof(*zprev))
                error(CONSISTENCY, "code error in refine_first");
        if (!silent) {
                printf("\tFirst refinement pass...");
                fflush(stdout);
        }
        memset((void *)esamp, '\0', sizeof(int)*hres*vres);
        /*
         * In our initial pass, we look for lower error pixels from
         * the same objects in the previous frame, and copy them here.
         */
        for (y = vres; y--; )
            for (x = hres; x--; ) {
                n = fndx(x, y);
                if (obuffer[n] == OVOID)
                        continue;
                if (xmbuffer[n] == MO_UNK)
                        continue;
                xp = x + xmbuffer[n];
                if ((xp < 0) | (xp >= hres))
                        continue;
                yp = y + ymbuffer[n];
                if ((yp < 0) | (yp >= vres))
                        continue;
                np = fndx(xp, yp);
                                        /* make sure we hit same object */
                if (oprev[np] != obuffer[n])
                        continue;
                                        /* is previous frame error lower? */
                if (aprev[np] < AMIN + ATIDIFF)
                        continue;
                if (aprev[np] <= abuffer[n] + ATIDIFF)
                        continue;
                                        /* shadow & highlight detection */
                if (abuffer[n] > hl_erri &&
                                getclosest(&neigh, 1, x, y) &&
                                bigdiff(cbuffer[neigh], cprev[np],
                                        HL_ERR*(.9+.2*frandom())))
                        continue;
                abuffer[n] = aprev[np] - ATIDIFF;
                copycolor(cbuffer[n], cprev[np]);
                esamp[n] = 1;           /* record extrapolated sample */
                nextra++;
            }
        for (n = hres*vres; n--; )      /* update sample counts */
                if (esamp[n])
                        sbuffer[n] = 1;
        if (!silent)
                printf("extrapolated %d pixels\n", nextra);
        return(1);
}


/*
 * We use a recursive computation of the conspicuity
 * map to avoid associated memory costs and simplify
 * coding.  We create a virtual image pyramid, pooling
 * variance calculations, etc.  The top of the pyramid
 * corresponds to the foveal resolution, as there should
 * not be any interesting mechanisms above this level.
 */

#define CSF_C0  1.14
#define CSF_C1  0.67
#define CSF_C2  1.7
#define CSF_S1  6.1
#define CSF_S2  7.3
#define CSF_P1  45.9
#define CSF_PC  (30./45.9*CSF_P1)
#define CSF_VR0 0.15
#define CSF_VRC 80.

struct ConspSum {
        COLOR           vsum;           /* value sum */
        COLOR           v2sum;          /* value^2 sum */
        long            nsamp;          /* number of samples */
        long            xmsum;          /* x-motion sum */
        long            ymsum;          /* y-motion sum */
        int             npix;           /* number of pixels */
        double          hls;            /* high-level saliency */
};

static double   pixel_deg;      /* base pixel frequency */
static int      fhsiz, fvsiz;   /* foveal subimage size */

static void clr_consp(struct ConspSum *cs);
static void sum_consp(struct ConspSum *cdest, struct ConspSum *cs);
static void est_consp(int x0, int y0, int x1, int y1, struct ConspSum *cs);
static void subconspicuity(int x0, int y0, int x1, int y1, struct ConspSum *cs);

static void
clr_consp(                      /* initialize a conspicuity sum */
        register struct ConspSum        *cs
)
{
        if (cs == NULL)
                return;
        setcolor(cs->vsum, 0., 0., 0.);
        setcolor(cs->v2sum, 0., 0., 0.);
        cs->nsamp = 0;
        cs->xmsum = cs->ymsum = 0;
        cs->npix = 0;
        cs->hls = 0;
}

static void
sum_consp(              /* sum in conspicuity result */
        register struct ConspSum        *cdest,
        register struct ConspSum        *cs
)
{
        if ((cdest == NULL) | (cs == NULL))
                return;
        addcolor(cdest->vsum, cs->vsum);
        addcolor(cdest->v2sum, cs->v2sum);
        cdest->nsamp += cs->nsamp;
        cdest->xmsum += cs->xmsum;
        cdest->ymsum += cs->ymsum;
        cdest->npix += cs->npix;
        if (cs->hls > cdest->hls)
                cdest->hls = cs->hls;
}

static void
est_consp(      /* estimate error conspicuity & update */
        int     x0,
        int     y0,
        int     x1,
        int     y1,
        register struct ConspSum        *cs
)
{
        double  rad2, mtn2, cpd, vm, vr, csf, eest;
                                                /* do we care? */
        if (cs->hls <= FTINY)
                return;
                                                /* get relative error */
        if (cs->nsamp < NSAMPOK) {
                int     neigh[NSAMPOK];         /* gather neighbors */
                eest = comperr(neigh,
                        getclosest(neigh, NSAMPOK, (x0+x1)>>1, (y0+y1)>>1),
                                cs->nsamp);
        } else
                eest = estimaterr(cs->vsum, cs->v2sum, cs->nsamp, cs->nsamp);
        
        if ((x0 == x1-1) & (y0 == y1-1)) {      /* update pixel error */
                int     n = fndx(x0, y0);
                int     ai;
                int     ne;
                if (sbuffer[n] >= 255) {
                        abuffer[n] = ADISTANT;
                } else {
                        ai = errori(eest);
                        if (ai < AMIN) ai = AMIN;
                        else if (ai >= ADISTANT) ai = ADISTANT-1;
                        abuffer[n] = ai;
                                                /* can't improve on closest */
                        if (!cs->nsamp && getclosest(&ne, 1, x0, y0) &&
                                        abuffer[ne] < ai &&
                                        abuffer[ne] >= AMIN)
                                abuffer[n] = abuffer[ne];
                }
        }
                                                /* compute radius^2 */
        rad2 = 0.125*((x1-x0)*(x1-x0) + (y1-y0)*(y1-y0));

                                                /* average motion^2 */
        mtn2 = (double)cs->xmsum*cs->xmsum + (double)cs->ymsum*cs->ymsum;
        mtn2 /= (double)(cs->npix*cs->npix);
                                                /* motion blur hides us? */
        if (mblur*mblur*mtn2 >= 4.*rad2)
                return;
                                                /* too small to see? */
        cpd = pixel_deg * pixel_deg / rad2;
        if (cpd > CSF_PC*CSF_PC)
                return;
        cpd = sqrt(cpd);
                                                /* compute CSF [Daley98] */
        vm = rate * sqrt(mtn2) / pixel_deg;
        vr = cs->hls/hlsmax*vm + CSF_VR0;       /* use hls tracking eff. */
        if (vr > CSF_VRC) vr = CSF_VRC;
        vr = vm - vr;
        if (vr < 0) vr = -vr;
        csf = log(CSF_C2*(1./3.)*vr);
        if (csf < 0) csf = -csf;
        csf = CSF_S1 + CSF_S2*csf*csf*csf;
        csf *= CSF_C0*CSF_C2*4.*PI*PI*CSF_C1*CSF_C1*cpd*cpd;
        csf *= exp(-CSF_C1*4.*PI/CSF_P1*(CSF_C2*vr + 2.)*cpd);
                                                /* compute visible error */
        eest = eest*csf/ndthresh - 1.;
        if (eest <= FTINY)
                return;
                                                /* scale by saleincy */
        eest *= cs->hls;
                                                /* worth the bother? */
        if (eest <= .01)
                return;
                                                /* put into map */
        for ( ; y0 < y1; y0++) {
                float   *em0 = cerrmap + fndx(x0, y0);
                register float  *emp = em0 + (x1-x0);
                while (emp-- > em0)
                        if (eest > *emp)
                                *emp = eest;
        }
        cerrzero = 0;
}

static void
subconspicuity( /* compute subportion of conspicuity */
        int     x0,
        int     y0,
        int     x1,
        int     y1,
        struct ConspSum *cs
)
{
        struct ConspSum mysum;
        int     i;

        if ((x0 >= x1) | (y0 >= y1))
                error(CONSISTENCY, "bad call to subconspicuity");

        clr_consp(&mysum);                      /* prepare sum */

        if ((x0 == x1-1) & (y0 == y1-1)) {      /* single pixel */
                double  hls;
                register int    n = fndx(x0, y0);
                if (sbuffer[n]) {
                        copycolor(mysum.vsum, cbuffer[n]);
                        copycolor(mysum.v2sum, val2map[n]);
                        mysum.nsamp = sbuffer[n];
                }
                if ((mysum.xmsum = xmbuffer[n]) == MO_UNK)
                        mysum.xmsum = 0;
                else
                        mysum.ymsum = ymbuffer[n];
                mysum.npix = 1;
                                                /* max. hls in fovea */
                mysum.hls = obj_prio(obuffer[n]);
                if (x0 >= fhsiz) {
                        hls = obj_prio(obuffer[fndx(x0-fhsiz,y0)]);
                        if (hls > mysum.hls) mysum.hls = hls;
                }
                if (x0 < hres-fhsiz) {
                        hls = obj_prio(obuffer[fndx(x0+fhsiz,y0)]);
                        if (hls > mysum.hls) mysum.hls = hls;
                }
                if (y0 >= fvsiz) {
                        hls = obj_prio(obuffer[fndx(x0,y0-fvsiz)]);
                        if (hls > mysum.hls) mysum.hls = hls;
                }
                if (y0 < vres-fvsiz) {
                        hls = obj_prio(obuffer[fndx(x0,y0+fvsiz)]);
                        if (hls > mysum.hls) mysum.hls = hls;
                }
        } else if (x0 == x1-1) {                /* vertical pair */
                for (i = y0 ; i < y1; i++)
                        subconspicuity(x0, i, x1, i+1, &mysum);
        } else if (y0 == y1-1) {                /* horizontal pair */
                for (i = x0 ; i < x1; i++)
                        subconspicuity(i, y0, i+1, y1, &mysum);
        } else {                                /* rectangle */
                subconspicuity(x0, y0, (x0+x1)>>1, (y0+y1)>>1, &mysum);
                subconspicuity((x0+x1)>>1, y0, x1, (y0+y1)>>1, &mysum);
                subconspicuity(x0, (y0+y1)>>1, (x0+x1)>>1, y1, &mysum);
                subconspicuity((x0+x1)>>1, (y0+y1)>>1, x1, y1, &mysum);
        }
                                                /* update conspicuity */
        est_consp(x0, y0, x1, y1, &mysum);
                                                /* sum into return value */
        sum_consp(cs, &mysum);
}

extern void
conspicuity(void)                       /* compute conspicuous error map */
{
        int     fhres, fvres;
        int     fx, fy;
                                        /* reuse previous z-buffer */
        cerrmap = (float *)zprev;
        memset((void *)cerrmap, '\0', sizeof(float)*hres*vres);
        cerrzero = 1;
                                        /* compute base pixel frequency */
        pixel_deg = .5*(hres/vw.horiz + vres/vw.vert);
                                        /* compute foveal resolution */
        fhres = vw.horiz/FOV_DEG + 0.5;
        if (fhres <= 0) fhres = 1;
        else if (fhres > hres) fhres = hres;
        fvres = vw.vert/FOV_DEG + 0.5;
        if (fvres <= 0) fvres = 1;
        else if (fvres > vres) fvres = vres;
        fhsiz = hres/fhres;
        fvsiz = vres/fvres;
                                        /* call our foveal subroutine */
        for (fy = fvres; fy--; )
                for (fx = fhres; fx--; )
                        subconspicuity(hres*fx/fhres, vres*fy/fvres,
                                        hres*(fx+1)/fhres, vres*(fy+1)/fvres,
                                        NULL);
}


/*
 * The following structure is used to collect data on the
 * initial error in the ambient value estimate, in order
 * to correct for it in the subsequent frames.
 */
static struct AmbSum {
        double          diffsum[3];     /* sum of (correct - ambval) */
        long            nsamps;         /* number of values in sum */
}               *asump = NULL;


static int
ppri_cmp(               /* pixel priority comparison */
        const void *pp1,
        const void *pp2
)
{
        double  se1 = cerrmap[*(const int *)pp1];
        double  se2 = cerrmap[*(const int *)pp2];
        int     adiff;
                                        /* higher conspicuity to front */
        if (se1 < se2) return(1);
        if (se1 > se2) return(-1);
                                        /* else higher error to front */
        adiff = (int)abuffer[*(const int *)pp1] -
                        (int)abuffer[*(const int *)pp2];
        if (adiff)
                return(adiff);
                                        /* else fewer samples to front */
        return((int)sbuffer[*(const int *)pp1] -
                        (int)sbuffer[*(const int *)pp2]);
}


static int
ray_refine(                     /* refine the given pixel by tracing a ray */
        register int    n
)
{
        RAY     ir;
        COLOR   ctmp;
        int     i;

        if (n < 0) {                    /* fetching stragglers */
                if (nprocs <= 1 || !ray_presult(&ir, 0))
                        return(-1);
                n = ir.rno;
        } else {                        /* else tracing a new ray */
                double  hv[2];
                if (sbuffer[n] >= 255)  /* reached limit? */
                        return(-1);
                sample_pos(hv, n%hres, n/hres, sbuffer[n]);
                ir.rmax = viewray(ir.rorg, ir.rdir, &vw, hv[0], hv[1]);
                if (ir.rmax < -FTINY)
                        return(-1);
                if (nprocs > 1) {
                        int     rval;
                        rayorigin(&ir, PRIMARY, NULL, NULL);
                        ir.rno = n;
                        rval = ray_pqueue(&ir);
                        if (!rval)
                                return(-1);
                        if (rval < 0)
                                quit(1);
                        n = ir.rno;
                } else
                        ray_trace(&ir);
        }
        if (abuffer[n] == ALOWQ && asump != NULL) {
                if (sbuffer[n] != 1)
                        error(CONSISTENCY, "bad code in ray_refine");
                if (getambcolor(ctmp, obuffer[n]) &&
                                (colval(ctmp,RED) > 0.01) &
                                (colval(ctmp,GRN) > 0.01) &
                                (colval(ctmp,BLU) > 0.01)) {
                        for (i = 0; i < 3; i++)
                                asump->diffsum[i] +=
                                    (colval(ir.rcol,i) - colval(cbuffer[n],i))
                                        / colval(ctmp,i);
                        asump->nsamps++;
                }
                sbuffer[n] = 0;
        }
        setcolor(ctmp,
                colval(ir.rcol,RED)*colval(ir.rcol,RED),
                colval(ir.rcol,GRN)*colval(ir.rcol,GRN),
                colval(ir.rcol,BLU)*colval(ir.rcol,BLU));
        if (!sbuffer[n]) {                      /* first sample */
                copycolor(cbuffer[n], ir.rcol);
                copycolor(val2map[n], ctmp);
                abuffer[n] = AHIGHQ;
                sbuffer[n] = 1;
        } else {                                /* else sum in sample */
                addcolor(cbuffer[n], ir.rcol);
                addcolor(val2map[n], ctmp);
                sbuffer[n]++;
        }
        return(n);
}


static long
refine_rays(            /* compute refinement rays */
        long    nrays
)
{
        int     *pord;
        int     ntodo;
        long    rdone;
        int     i;
                                        /* skip if nothing significant */
        if (ndtset && cerrzero)
                return(0);
                                        /* initialize priority list */
        pord = (int *)malloc(sizeof(int)*hres*vres);
        for (i = hres*vres; i--; )
                pord[i] = i;
                                        /* sort our priorities */
        ntodo = hres*vres;
        if (nrays < ntodo)
                qsort((void *)pord, hres*vres, sizeof(int), ppri_cmp);
        i = 0;
                                        /* trace rays in list */
        for (rdone = 0; rdone < nrays; rdone++) {
                if (ndtset && i >= 1000 && cerrmap[pord[i]] <= FTINY)
                        ntodo = i;
                if (i >= ntodo) {       /* redo conspicuity & priority */
                        while (ray_refine(-1) >= 0)
                                ;
                        conspicuity();
                        if (ndtset && cerrzero)
                                break;
                        qsort((void *)pord, hres*vres, sizeof(int), ppri_cmp);
                        ntodo = hres*vres/8;
                        i = 0;
                }
                                        /* sample next pixel */
                ray_refine(pord[i++]);
        }
                                        /* clean up and return */
        while (ray_refine(-1) >= 0)
                ;
        free((void *)pord);
        return(rdone);
}


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