src/rt/ambient.c

/* Copyright (c) 1986 Regents of the University of California */

#ifndef lint
static char SCCSid[] = "$SunId$ LBL";
#endif

/*
 *  ambient.c - routines dealing with ambient (inter-reflected) component.
 *
 *  The macro AMBFLUSH (if defined) is the number of ambient values
 *      to wait before flushing to the ambient file.
 *
 *     5/9/86
 */

#include  "ray.h"

#include  "octree.h"

#include  "otypes.h"

#include  "random.h"

#define  OCTSCALE       0.5     /* ceil((valid rad.)/(cube size)) */

extern CUBE  thescene;          /* contains space boundaries */

extern COLOR  ambval;           /* global ambient component */
extern double  ambacc;          /* ambient accuracy */
extern int  ambres;             /* ambient resolution */
extern int  ambdiv;             /* number of divisions for calculation */
extern int  ambssamp;           /* number of super-samples */
extern int  ambounce;           /* number of ambient bounces */
extern char  *amblist[];        /* ambient include/exclude list */
extern int  ambincl;            /* include == 1, exclude == 0 */

#define  MAXASET        511     /* maximum number of elements in ambient set */
OBJECT  ambset[MAXASET+1]={0};  /* ambient include/exclude set */

double  maxarad;                /* maximum ambient radius */
double  minarad;                /* minimum ambient radius */

/*
 * Since we've defined our vectors as float below to save space,
 * watch out for changes in the definitions of VCOPY() and DOT().
 */
typedef struct ambval {
        float  pos[3];          /* position in space */
        float  dir[3];          /* normal direction */
        int  lvl;               /* recursion level of parent ray */
        float  weight;          /* weight of parent ray */
        COLOR  val;             /* computed ambient value */
        float  rad;             /* validity radius */
        struct ambval  *next;   /* next in list */
}  AMBVAL;                      /* ambient value */

typedef struct ambtree {
        AMBVAL  *alist;         /* ambient value list */
        struct ambtree  *kid;   /* 8 child nodes */
}  AMBTREE;                     /* ambient octree */

typedef struct {
        float  k;               /* error contribution per sample */
        COLOR  v;               /* ray sum */
        int  n;                 /* number of samples */
        short  t, p;            /* theta, phi indices */
}  AMBSAMP;                     /* ambient sample */

static AMBTREE  atrunk;         /* our ambient trunk node */

static FILE  *ambfp = NULL;     /* ambient file pointer */

#define  newambval()    (AMBVAL *)bmalloc(sizeof(AMBVAL))

#define  newambtree()   (AMBTREE *)calloc(8, sizeof(AMBTREE))

double  sumambient(), doambient(), makeambient();


setambient(afile)                       /* initialize calculation */
char  *afile;
{
        long  ftell();
        AMBVAL  amb;

        maxarad = thescene.cusize / 2.0;                /* maximum radius */
                                                        /* minimum radius */
        minarad = ambres > 0 ? thescene.cusize/ambres : 0.0;

                                        /* open ambient file */
        if (afile != NULL)
                if ((ambfp = fopen(afile, "r+")) != NULL) {
                        while (fread((char *)&amb,sizeof(AMBVAL),1,ambfp) == 1)
                                avinsert(&amb, &atrunk, thescene.cuorg,
                                                thescene.cusize);
                                                        /* align */
                        fseek(ambfp, -(ftell(ambfp)%sizeof(AMBVAL)), 1);
                } else if ((ambfp = fopen(afile, "w")) == NULL) {
                        sprintf(errmsg, "cannot open ambient file \"%s\"",
                                        afile);
                        error(SYSTEM, errmsg);
                }
}


ambnotify(obj)                  /* record new modifier */
OBJECT  obj;
{
        static int  hitlimit = 0;
        register OBJREC  *o = objptr(obj);
        register char  **amblp;

        if (hitlimit || !ismodifier(o->otype))
                return;
        for (amblp = amblist; *amblp != NULL; amblp++)
                if (!strcmp(o->oname, *amblp)) {
                        if (ambset[0] >= MAXASET) {
                                error(WARNING, "too many modifiers in ambient list");
                                hitlimit++;
                                return;         /* should this be fatal? */
                        }
                        insertelem(ambset, obj);
                        return;
                }
}


ambient(acol, r)                /* compute ambient component for ray */
COLOR  acol;
register RAY  *r;
{
        static int  rdepth = 0;                 /* ambient recursion */
        double  wsum;

        rdepth++;                               /* increment level */

        if (ambdiv <= 0)                        /* no ambient calculation */
                goto dumbamb;
                                                /* check number of bounces */
        if (rdepth > ambounce)
                goto dumbamb;
                                                /* check ambient list */
        if (ambincl != -1 && r->ro != NULL &&
                        ambincl != inset(ambset, r->ro->omod))
                goto dumbamb;

        if (ambacc <= FTINY) {                  /* no ambient storage */
                if (doambient(acol, r) == 0.0)
                        goto dumbamb;
                goto done;
        }
                                                /* get ambient value */
        setcolor(acol, 0.0, 0.0, 0.0);
        wsum = sumambient(acol, r, &atrunk, thescene.cuorg, thescene.cusize);
        if (wsum > FTINY)
                scalecolor(acol, 1.0/wsum);
        else if (makeambient(acol, r) == 0.0)
                goto dumbamb;
        goto done;

dumbamb:                                        /* return global value */
        copycolor(acol, ambval);
done:                                           /* must finish here! */
        rdepth--;
}


double
sumambient(acol, r, at, c0, s)          /* get interpolated ambient value */
COLOR  acol;
register RAY  *r;
AMBTREE  *at;
FVECT  c0;
double  s;
{
        extern double  sqrt();
        double  d, e1, e2, wt, wsum;
        COLOR  ct;
        FVECT  ck0;
        int  i;
        register int  j;
        register AMBVAL  *av;
                                        /* do this node */
        wsum = 0.0;
        for (av = at->alist; av != NULL; av = av->next) {
                /*
                 *  Ray strength test.
                 */
                if (av->lvl > r->rlvl || av->weight < r->rweight-FTINY)
                        continue;
                /*
                 *  Ambient radius test.
                 */
                e1 = 0.0;
                for (j = 0; j < 3; j++) {
                        d = av->pos[j] - r->rop[j];
                        e1 += d * d;
                }
                e1 /= av->rad * av->rad;
                if (e1 > ambacc*ambacc*1.21)
                        continue;
                /*
                 *  Normal direction test.
                 */
                e2 = (1.0 - DOT(av->dir, r->ron)) * r->rweight;
                if (e2 < 0.0) e2 = 0.0;
                if (e1 + e2 > ambacc*ambacc*1.21)
                        continue;
                /*
                 *  Ray behind test.
                 */
                d = 0.0;
                for (j = 0; j < 3; j++)
                        d += (r->rop[j] - av->pos[j]) *
                                        (av->dir[j] + r->ron[j]);
                if (d*0.5 < -minarad*ambacc)
                        continue;
                /*
                 *  Jittering final test reduces image artifacts.
                 */
                wt = sqrt(e1) + sqrt(e2);
                wt *= .9 + .2*frandom();
                if (wt > ambacc)
                        continue;
                if (wt <= 1e-3)
                        wt = 1e3;
                else
                        wt = 1.0 / wt;
                wsum += wt;
                copycolor(ct, av->val);
                scalecolor(ct, wt);
                addcolor(acol, ct);
        }
        if (at->kid == NULL)
                return(wsum);
                                        /* do children */
        s *= 0.5;
        for (i = 0; i < 8; i++) {
                for (j = 0; j < 3; j++) {
                        ck0[j] = c0[j];
                        if (1<<j & i)
                                ck0[j] += s;
                        if (r->rop[j] < ck0[j] - OCTSCALE*s)
                                break;
                        if (r->rop[j] > ck0[j] + (1.0+OCTSCALE)*s)
                                break;
                }
                if (j == 3)
                        wsum += sumambient(acol, r, at->kid+i, ck0, s);
        }
        return(wsum);
}


double
makeambient(acol, r)            /* make a new ambient value */
COLOR  acol;
register RAY  *r;
{
        AMBVAL  amb;

        amb.rad = doambient(acol, r);           /* compute ambient */
        if (amb.rad == 0.0)
                return(0.0);
                                                /* store it */
        VCOPY(amb.pos, r->rop);
        VCOPY(amb.dir, r->ron);
        amb.lvl = r->rlvl;
        amb.weight = r->rweight;
        copycolor(amb.val, acol);
                                                /* insert into tree */
        avinsert(&amb, &atrunk, thescene.cuorg, thescene.cusize);
        avsave(&amb);                           /* write to file */
        return(amb.rad);
}


double
doambient(acol, r)                      /* compute ambient component */
COLOR  acol;
register RAY  *r;
{
        extern int  ambcmp();
        extern double  sin(), cos(), sqrt();
        int  hlist[4];
        double  phi, xd, yd, zd;
        double  b, b2;
        register AMBSAMP  *div;
        AMBSAMP  dnew;
        RAY  ar;
        FVECT  ux, uy;
        double  arad;
        int  ndivs, nt, np, ns, ne, i, j;
        register int  k;

        setcolor(acol, 0.0, 0.0, 0.0);
                                        /* set number of divisions */
        nt = sqrt(ambdiv * r->rweight * 0.5) + 0.5;
        np = 2 * nt;
        ndivs = nt * np;
                                        /* check first */
        if (ndivs == 0 || rayorigin(&ar, r, AMBIENT, 0.5) < 0)
                return(0.0);
                                        /* set number of super-samples */
        ns = ambssamp * r->rweight + 0.5;
        if (ns > 0) {
                div = (AMBSAMP *)malloc(ndivs*sizeof(AMBSAMP));
                if (div == NULL)
                        error(SYSTEM, "out of memory in doambient");
        } else
                div = NULL;
                                        /* make axes */
        uy[0] = uy[1] = uy[2] = 0.0;
        for (k = 0; k < 3; k++)
                if (r->ron[k] < 0.6 && r->ron[k] > -0.6)
                        break;
        uy[k] = 1.0;
        fcross(ux, r->ron, uy);
        normalize(ux);
        fcross(uy, ux, r->ron);
                                        /* set up urand */
        hlist[0] = r->rno;
                                                /* sample divisions */
        arad = 0.0;
        ne = 0;
        for (i = 0; i < nt; i++) {
                hlist[1] = i;
                for (j = 0; j < np; j++) {
                        rayorigin(&ar, r, AMBIENT, 0.5);        /* pretested */
                        hlist[2] = j;
                        hlist[3] = 0;
                        zd = sqrt((i+urand(ilhash(hlist,4)))/nt);
                        hlist[3] = 1;
                        phi = 2.0*PI * (j+urand(ilhash(hlist,4)))/np;
                        xd = cos(phi) * zd;
                        yd = sin(phi) * zd;
                        zd = sqrt(1.0 - zd*zd);
                        for (k = 0; k < 3; k++)
                                ar.rdir[k] = xd*ux[k]+yd*uy[k]+zd*r->ron[k];
                        dimlist[ndims++] = i*np + j + 38813;
                        rayvalue(&ar);
                        ndims--;
                        if (ar.rot < FHUGE)
                                arad += 1.0 / ar.rot;
                        if (div != NULL) {              /* save division */
                                div[ne].k = 0.0;
                                copycolor(div[ne].v, ar.rcol);
                                div[ne].n = 0;
                                div[ne].t = i; div[ne].p = j;
                                                        /* sum errors */
                                b = bright(ar.rcol);
                                if (i > 0) {            /* from above */
                                        b2 = bright(div[ne-np].v) - b;
                                        b2 *= b2 * 0.25;
                                        div[ne].k += b2;
                                        div[ne].n++;
                                        div[ne-np].k += b2;
                                        div[ne-np].n++;
                                }
                                if (j > 0) {            /* from behind */
                                        b2 = bright(div[ne-1].v) - b;
                                        b2 *= b2 * 0.25;
                                        div[ne].k += b2;
                                        div[ne].n++;
                                        div[ne-1].k += b2;
                                        div[ne-1].n++;
                                }
                                if (j == np-1) {        /* around */
                                        b2 = bright(div[ne-(np-1)].v) - b;
                                        b2 *= b2 * 0.25;
                                        div[ne].k += b2;
                                        div[ne].n++;
                                        div[ne-(np-1)].k += b2;
                                        div[ne-(np-1)].n++;
                                }
                                ne++;
                        } else
                                addcolor(acol, ar.rcol);
                }
        }
        for (k = 0; k < ne; k++) {              /* compute errors */
                if (div[k].n > 1)
                        div[k].k /= div[k].n;
                div[k].n = 1;
        }
                                                /* sort the divisions */
        qsort(div, ne, sizeof(AMBSAMP), ambcmp);
                                                /* skim excess */
        while (ne > ns) {
                ne--;
                addcolor(acol, div[ne].v);
        }
                                                /* super-sample */
        for (i = ns; i > 0; i--) {
                rayorigin(&ar, r, AMBIENT, 0.5);        /* pretested */
                hlist[1] = div[0].t;
                hlist[2] = div[0].p;
                hlist[3] = 0;
                zd = sqrt((div[0].t+urand(ilhash(hlist,4)+div[0].n))/nt);
                hlist[3] = 1;
                phi = 2.0*PI * (div[0].p+urand(ilhash(hlist,4)+div[0].n))/np;
                xd = cos(phi) * zd;
                yd = sin(phi) * zd;
                zd = sqrt(1.0 - zd*zd);
                for (k = 0; k < 3; k++)
                        ar.rdir[k] = xd*ux[k]+yd*uy[k]+zd*r->ron[k];
                dimlist[ndims++] = div[0].t*np + div[0].p + 38813;
                rayvalue(&ar);
                ndims--;
                rayvalue(&ar);
                if (ar.rot < FHUGE)
                        arad += 1.0 / ar.rot;
                                                /* recompute error */
                copycolor(dnew.v, div[0].v);
                addcolor(dnew.v, ar.rcol);
                dnew.n = div[0].n + 1;
                dnew.t = div[0].t; dnew.p = div[0].p;
                b2 = bright(dnew.v)/dnew.n - bright(ar.rcol);
                b2 = b2*b2 + div[0].k*(div[0].n*div[0].n);
                dnew.k = b2/(dnew.n*dnew.n);
                                                /* reinsert */
                for (k = 0; k < ne-1 && dnew.k < div[k+1].k; k++)
                        copystruct(&div[k], &div[k+1]);
                copystruct(&div[k], &dnew);

                if (ne >= i) {          /* extract darkest division */
                        ne--;
                        if (div[ne].n > 1) {
                                b = 1.0/div[ne].n;
                                scalecolor(div[ne].v, b);
                                div[ne].n = 1;
                        }
                        addcolor(acol, div[ne].v);
                }
        }
        scalecolor(acol, 1.0/ndivs);
        if (arad <= FTINY)
                arad = FHUGE;
        else
                arad = (ndivs+ns) / arad / sqrt(r->rweight);
        if (arad > maxarad)
                arad = maxarad;
        else if (arad < minarad)
                arad = minarad;
        if (div != NULL)
                free((char *)div);
        return(arad);
}


static int
ambcmp(d1, d2)                          /* decreasing order */
AMBSAMP  *d1, *d2;
{
        if (d1->k < d2->k)
                return(1);
        if (d1->k > d2->k)
                return(-1);
        return(0);
}


static
avsave(av)                              /* save an ambient value */
AMBVAL  *av;
{
#ifdef  AMBFLUSH
        static int  nunflshed = 0;
#endif
        if (ambfp == NULL)
                return;
        if (fwrite((char *)av, sizeof(AMBVAL), 1, ambfp) != 1)
                goto writerr;
#ifdef  AMBFLUSH
        if (++nunflshed >= AMBFLUSH) {
                if (fflush(ambfp) == EOF)
                        goto writerr;
                nunflshed = 0;
        }
#endif
        return;
writerr:
        error(SYSTEM, "error writing ambient file");
}


static
avinsert(aval, at, c0, s)               /* insert ambient value in a tree */
AMBVAL  *aval;
register AMBTREE  *at;
FVECT  c0;
double  s;
{
        FVECT  ck0;
        int  branch;
        register AMBVAL  *av;
        register int  i;

        if ((av = newambval()) == NULL)
                goto memerr;
        copystruct(av, aval);
        VCOPY(ck0, c0);
        while (s*(OCTSCALE/2) > av->rad*ambacc) {
                if (at->kid == NULL)
                        if ((at->kid = newambtree()) == NULL)
                                goto memerr;
                s *= 0.5;
                branch = 0;
                for (i = 0; i < 3; i++)
                        if (av->pos[i] > ck0[i] + s) {
                                ck0[i] += s;
                                branch |= 1 << i;
                        }
                at = at->kid + branch;
        }
        av->next = at->alist;
        at->alist = av;
        return;
memerr:
        error(SYSTEM, "out of memory in avinsert");
}
Revision:	1.13
Committed:	Thu Jun 6 10:48:50 1991 UTC (32 years, 10 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 1.12:	+16 -8 lines
Log Message:	some minor improvements
#	User	Rev	Content
1	greg	1.1	/* Copyright (c) 1986 Regents of the University of California */
2
3			#ifndef lint
4			static char SCCSid[] = "$SunId$ LBL";
5			#endif
6
7			/*
8			* ambient.c - routines dealing with ambient (inter-reflected) component.
9			*
10			* The macro AMBFLUSH (if defined) is the number of ambient values
11			* to wait before flushing to the ambient file.
12			*
13			* 5/9/86
14			*/
15
16			#include "ray.h"
17
18			#include "octree.h"
19
20	greg	1.8	#include "otypes.h"
21
22	greg	1.1	#include "random.h"
23
24			#define OCTSCALE 0.5 /* ceil((valid rad.)/(cube size)) */
25
26			extern CUBE thescene; /* contains space boundaries */
27
28			extern COLOR ambval; /* global ambient component */
29			extern double ambacc; /* ambient accuracy */
30			extern int ambres; /* ambient resolution */
31			extern int ambdiv; /* number of divisions for calculation */
32			extern int ambssamp; /* number of super-samples */
33			extern int ambounce; /* number of ambient bounces */
34			extern char amblist[]; / ambient include/exclude list */
35			extern int ambincl; /* include == 1, exclude == 0 */
36
37	greg	1.11	#define MAXASET 511 /* maximum number of elements in ambient set */
38			OBJECT ambset[MAXASET+1]={0}; /* ambient include/exclude set */
39	greg	1.1
40			double maxarad; /* maximum ambient radius */
41			double minarad; /* minimum ambient radius */
42
43	greg	1.13	/*
44			* Since we've defined our vectors as float below to save space,
45			* watch out for changes in the definitions of VCOPY() and DOT().
46			*/
47	greg	1.1	typedef struct ambval {
48	greg	1.13	float pos[3]; /* position in space */
49			float dir[3]; /* normal direction */
50	greg	1.1	int lvl; /* recursion level of parent ray */
51			float weight; /* weight of parent ray */
52			COLOR val; /* computed ambient value */
53			float rad; /* validity radius */
54			struct ambval next; / next in list */
55			} AMBVAL; /* ambient value */
56
57			typedef struct ambtree {
58			AMBVAL alist; / ambient value list */
59			struct ambtree kid; / 8 child nodes */
60			} AMBTREE; /* ambient octree */
61
62			typedef struct {
63			float k; /* error contribution per sample */
64			COLOR v; /* ray sum */
65			int n; /* number of samples */
66			short t, p; /* theta, phi indices */
67			} AMBSAMP; /* ambient sample */
68
69			static AMBTREE atrunk; /* our ambient trunk node */
70
71			static FILE ambfp = NULL; / ambient file pointer */
72
73			#define newambval() (AMBVAL *)bmalloc(sizeof(AMBVAL))
74
75			#define newambtree() (AMBTREE *)calloc(8, sizeof(AMBTREE))
76
77			double sumambient(), doambient(), makeambient();
78
79
80			setambient(afile) /* initialize calculation */
81			char *afile;
82			{
83			long ftell();
84			AMBVAL amb;
85	greg	1.8
86	greg	1.1	maxarad = thescene.cusize / 2.0; /* maximum radius */
87			/* minimum radius */
88			minarad = ambres > 0 ? thescene.cusize/ambres : 0.0;
89
90			/* open ambient file */
91			if (afile != NULL)
92			if ((ambfp = fopen(afile, "r+")) != NULL) {
93	greg	1.9	while (fread((char *)&amb,sizeof(AMBVAL),1,ambfp) == 1)
94	greg	1.1	avinsert(&amb, &atrunk, thescene.cuorg,
95			thescene.cusize);
96			/* align */
97			fseek(ambfp, -(ftell(ambfp)%sizeof(AMBVAL)), 1);
98			} else if ((ambfp = fopen(afile, "w")) == NULL) {
99			sprintf(errmsg, "cannot open ambient file \"%s\"",
100			afile);
101			error(SYSTEM, errmsg);
102	greg	1.8	}
103			}
104
105
106			ambnotify(obj) /* record new modifier */
107			OBJECT obj;
108			{
109	greg	1.11	static int hitlimit = 0;
110	greg	1.8	register OBJREC *o = objptr(obj);
111			register char **amblp;
112
113	greg	1.11	if (hitlimit \|\| !ismodifier(o->otype))
114	greg	1.8	return;
115			for (amblp = amblist; *amblp != NULL; amblp++)
116			if (!strcmp(o->oname, *amblp)) {
117	greg	1.11	if (ambset[0] >= MAXASET) {
118			error(WARNING, "too many modifiers in ambient list");
119			hitlimit++;
120			return; /* should this be fatal? */
121			}
122	greg	1.8	insertelem(ambset, obj);
123			return;
124	greg	1.1	}
125			}
126
127
128			ambient(acol, r) /* compute ambient component for ray */
129			COLOR acol;
130			register RAY *r;
131			{
132			static int rdepth = 0; /* ambient recursion */
133			double wsum;
134
135			rdepth++; /* increment level */
136
137			if (ambdiv <= 0) /* no ambient calculation */
138			goto dumbamb;
139			/* check number of bounces */
140			if (rdepth > ambounce)
141			goto dumbamb;
142			/* check ambient list */
143			if (ambincl != -1 && r->ro != NULL &&
144			ambincl != inset(ambset, r->ro->omod))
145			goto dumbamb;
146
147			if (ambacc <= FTINY) { /* no ambient storage */
148			if (doambient(acol, r) == 0.0)
149			goto dumbamb;
150			goto done;
151			}
152			/* get ambient value */
153			setcolor(acol, 0.0, 0.0, 0.0);
154			wsum = sumambient(acol, r, &atrunk, thescene.cuorg, thescene.cusize);
155			if (wsum > FTINY)
156			scalecolor(acol, 1.0/wsum);
157			else if (makeambient(acol, r) == 0.0)
158			goto dumbamb;
159			goto done;
160
161			dumbamb: /* return global value */
162			copycolor(acol, ambval);
163			done: /* must finish here! */
164			rdepth--;
165			}
166
167
168			double
169			sumambient(acol, r, at, c0, s) /* get interpolated ambient value */
170			COLOR acol;
171			register RAY *r;
172			AMBTREE *at;
173			FVECT c0;
174			double s;
175			{
176			extern double sqrt();
177			double d, e1, e2, wt, wsum;
178			COLOR ct;
179			FVECT ck0;
180			int i;
181			register int j;
182			register AMBVAL *av;
183	greg	1.7	/* do this node */
184	greg	1.1	wsum = 0.0;
185			for (av = at->alist; av != NULL; av = av->next) {
186			/*
187			* Ray strength test.
188			*/
189			if (av->lvl > r->rlvl \|\| av->weight < r->rweight-FTINY)
190			continue;
191			/*
192			* Ambient radius test.
193			*/
194			e1 = 0.0;
195			for (j = 0; j < 3; j++) {
196			d = av->pos[j] - r->rop[j];
197			e1 += d * d;
198			}
199			e1 /= av->rad * av->rad;
200			if (e1 > ambaccambacc1.21)
201			continue;
202			/*
203			* Normal direction test.
204			*/
205			e2 = (1.0 - DOT(av->dir, r->ron)) * r->rweight;
206			if (e2 < 0.0) e2 = 0.0;
207			if (e1 + e2 > ambaccambacc1.21)
208			continue;
209			/*
210			* Ray behind test.
211			*/
212			d = 0.0;
213			for (j = 0; j < 3; j++)
214			d += (r->rop[j] - av->pos[j]) *
215			(av->dir[j] + r->ron[j]);
216	greg	1.13	if (d0.5 < -minaradambacc)
217	greg	1.1	continue;
218			/*
219			* Jittering final test reduces image artifacts.
220			*/
221			wt = sqrt(e1) + sqrt(e2);
222	greg	1.7	wt = .9 + .2frandom();
223	greg	1.6	if (wt > ambacc)
224	greg	1.1	continue;
225			if (wt <= 1e-3)
226			wt = 1e3;
227			else
228			wt = 1.0 / wt;
229			wsum += wt;
230			copycolor(ct, av->val);
231			scalecolor(ct, wt);
232			addcolor(acol, ct);
233	greg	1.7	}
234			if (at->kid == NULL)
235			return(wsum);
236			/* do children */
237			s *= 0.5;
238			for (i = 0; i < 8; i++) {
239			for (j = 0; j < 3; j++) {
240			ck0[j] = c0[j];
241			if (1<<j & i)
242			ck0[j] += s;
243			if (r->rop[j] < ck0[j] - OCTSCALE*s)
244			break;
245			if (r->rop[j] > ck0[j] + (1.0+OCTSCALE)*s)
246			break;
247			}
248			if (j == 3)
249			wsum += sumambient(acol, r, at->kid+i, ck0, s);
250	greg	1.1	}
251			return(wsum);
252			}
253
254
255			double
256			makeambient(acol, r) /* make a new ambient value */
257			COLOR acol;
258			register RAY *r;
259			{
260			AMBVAL amb;
261
262			amb.rad = doambient(acol, r); /* compute ambient */
263			if (amb.rad == 0.0)
264			return(0.0);
265			/* store it */
266			VCOPY(amb.pos, r->rop);
267			VCOPY(amb.dir, r->ron);
268			amb.lvl = r->rlvl;
269			amb.weight = r->rweight;
270			copycolor(amb.val, acol);
271			/* insert into tree */
272			avinsert(&amb, &atrunk, thescene.cuorg, thescene.cusize);
273			avsave(&amb); /* write to file */
274			return(amb.rad);
275			}
276
277
278			double
279			doambient(acol, r) /* compute ambient component */
280			COLOR acol;
281			register RAY *r;
282			{
283			extern int ambcmp();
284			extern double sin(), cos(), sqrt();
285	greg	1.12	int hlist[4];
286	greg	1.1	double phi, xd, yd, zd;
287	greg	1.3	double b, b2;
288	greg	1.1	register AMBSAMP *div;
289			AMBSAMP dnew;
290			RAY ar;
291			FVECT ux, uy;
292			double arad;
293			int ndivs, nt, np, ns, ne, i, j;
294			register int k;
295
296			setcolor(acol, 0.0, 0.0, 0.0);
297			/* set number of divisions */
298			nt = sqrt(ambdiv * r->rweight * 0.5) + 0.5;
299			np = 2 * nt;
300			ndivs = nt * np;
301			/* check first */
302			if (ndivs == 0 \|\| rayorigin(&ar, r, AMBIENT, 0.5) < 0)
303			return(0.0);
304			/* set number of super-samples */
305			ns = ambssamp * r->rweight + 0.5;
306			if (ns > 0) {
307			div = (AMBSAMP )malloc(ndivssizeof(AMBSAMP));
308			if (div == NULL)
309			error(SYSTEM, "out of memory in doambient");
310	greg	1.13	} else
311			div = NULL;
312	greg	1.1	/* make axes */
313			uy[0] = uy[1] = uy[2] = 0.0;
314			for (k = 0; k < 3; k++)
315			if (r->ron[k] < 0.6 && r->ron[k] > -0.6)
316			break;
317			uy[k] = 1.0;
318			fcross(ux, r->ron, uy);
319			normalize(ux);
320			fcross(uy, ux, r->ron);
321	greg	1.12	/* set up urand */
322			hlist[0] = r->rno;
323	greg	1.1	/* sample divisions */
324			arad = 0.0;
325			ne = 0;
326	greg	1.12	for (i = 0; i < nt; i++) {
327			hlist[1] = i;
328	greg	1.1	for (j = 0; j < np; j++) {
329			rayorigin(&ar, r, AMBIENT, 0.5); /* pretested */
330	greg	1.12	hlist[2] = j;
331			hlist[3] = 0;
332			zd = sqrt((i+urand(ilhash(hlist,4)))/nt);
333			hlist[3] = 1;
334			phi = 2.0PI (j+urand(ilhash(hlist,4)))/np;
335	greg	1.1	xd = cos(phi) * zd;
336			yd = sin(phi) * zd;
337			zd = sqrt(1.0 - zd*zd);
338			for (k = 0; k < 3; k++)
339			ar.rdir[k] = xdux[k]+yduy[k]+zd*r->ron[k];
340	greg	1.12	dimlist[ndims++] = i*np + j + 38813;
341	greg	1.1	rayvalue(&ar);
342	greg	1.12	ndims--;
343	greg	1.1	if (ar.rot < FHUGE)
344			arad += 1.0 / ar.rot;
345	greg	1.13	if (div != NULL) { /* save division */
346	greg	1.1	div[ne].k = 0.0;
347			copycolor(div[ne].v, ar.rcol);
348			div[ne].n = 0;
349			div[ne].t = i; div[ne].p = j;
350			/* sum errors */
351	greg	1.3	b = bright(ar.rcol);
352	greg	1.1	if (i > 0) { /* from above */
353	greg	1.3	b2 = bright(div[ne-np].v) - b;
354			b2 = b2 0.25;
355			div[ne].k += b2;
356	greg	1.1	div[ne].n++;
357	greg	1.3	div[ne-np].k += b2;
358	greg	1.1	div[ne-np].n++;
359			}
360			if (j > 0) { /* from behind */
361	greg	1.3	b2 = bright(div[ne-1].v) - b;
362			b2 = b2 0.25;
363			div[ne].k += b2;
364	greg	1.1	div[ne].n++;
365	greg	1.3	div[ne-1].k += b2;
366	greg	1.1	div[ne-1].n++;
367			}
368			if (j == np-1) { /* around */
369	greg	1.3	b2 = bright(div[ne-(np-1)].v) - b;
370			b2 = b2 0.25;
371			div[ne].k += b2;
372	greg	1.1	div[ne].n++;
373	greg	1.3	div[ne-(np-1)].k += b2;
374	greg	1.1	div[ne-(np-1)].n++;
375			}
376			ne++;
377			} else
378			addcolor(acol, ar.rcol);
379			}
380	greg	1.12	}
381	greg	1.1	for (k = 0; k < ne; k++) { /* compute errors */
382			if (div[k].n > 1)
383			div[k].k /= div[k].n;
384			div[k].n = 1;
385			}
386			/* sort the divisions */
387			qsort(div, ne, sizeof(AMBSAMP), ambcmp);
388			/* skim excess */
389			while (ne > ns) {
390			ne--;
391			addcolor(acol, div[ne].v);
392			}
393			/* super-sample */
394			for (i = ns; i > 0; i--) {
395			rayorigin(&ar, r, AMBIENT, 0.5); /* pretested */
396	greg	1.12	hlist[1] = div[0].t;
397			hlist[2] = div[0].p;
398			hlist[3] = 0;
399			zd = sqrt((div[0].t+urand(ilhash(hlist,4)+div[0].n))/nt);
400			hlist[3] = 1;
401			phi = 2.0PI (div[0].p+urand(ilhash(hlist,4)+div[0].n))/np;
402	greg	1.1	xd = cos(phi) * zd;
403			yd = sin(phi) * zd;
404			zd = sqrt(1.0 - zd*zd);
405			for (k = 0; k < 3; k++)
406			ar.rdir[k] = xdux[k]+yduy[k]+zd*r->ron[k];
407	greg	1.12	dimlist[ndims++] = div[0].t*np + div[0].p + 38813;
408			rayvalue(&ar);
409			ndims--;
410	greg	1.1	rayvalue(&ar);
411			if (ar.rot < FHUGE)
412			arad += 1.0 / ar.rot;
413			/* recompute error */
414			copycolor(dnew.v, div[0].v);
415			addcolor(dnew.v, ar.rcol);
416			dnew.n = div[0].n + 1;
417			dnew.t = div[0].t; dnew.p = div[0].p;
418	greg	1.3	b2 = bright(dnew.v)/dnew.n - bright(ar.rcol);
419	greg	1.5	b2 = b2b2 + div[0].k(div[0].n*div[0].n);
420			dnew.k = b2/(dnew.n*dnew.n);
421	greg	1.1	/* reinsert */
422			for (k = 0; k < ne-1 && dnew.k < div[k+1].k; k++)
423	greg	1.9	copystruct(&div[k], &div[k+1]);
424			copystruct(&div[k], &dnew);
425	greg	1.1
426			if (ne >= i) { /* extract darkest division */
427			ne--;
428	greg	1.13	if (div[ne].n > 1) {
429			b = 1.0/div[ne].n;
430			scalecolor(div[ne].v, b);
431			div[ne].n = 1;
432			}
433	greg	1.1	addcolor(acol, div[ne].v);
434			}
435			}
436			scalecolor(acol, 1.0/ndivs);
437			if (arad <= FTINY)
438			arad = FHUGE;
439			else
440			arad = (ndivs+ns) / arad / sqrt(r->rweight);
441			if (arad > maxarad)
442			arad = maxarad;
443			else if (arad < minarad)
444			arad = minarad;
445	greg	1.13	if (div != NULL)
446	greg	1.1	free((char *)div);
447			return(arad);
448			}
449
450
451			static int
452			ambcmp(d1, d2) /* decreasing order */
453			AMBSAMP d1, d2;
454			{
455			if (d1->k < d2->k)
456			return(1);
457			if (d1->k > d2->k)
458			return(-1);
459			return(0);
460			}
461
462
463			static
464			avsave(av) /* save an ambient value */
465			AMBVAL *av;
466			{
467			#ifdef AMBFLUSH
468			static int nunflshed = 0;
469			#endif
470			if (ambfp == NULL)
471			return;
472	greg	1.9	if (fwrite((char *)av, sizeof(AMBVAL), 1, ambfp) != 1)
473	greg	1.1	goto writerr;
474			#ifdef AMBFLUSH
475			if (++nunflshed >= AMBFLUSH) {
476			if (fflush(ambfp) == EOF)
477			goto writerr;
478			nunflshed = 0;
479			}
480			#endif
481			return;
482			writerr:
483			error(SYSTEM, "error writing ambient file");
484			}
485
486
487			static
488			avinsert(aval, at, c0, s) /* insert ambient value in a tree */
489			AMBVAL *aval;
490			register AMBTREE *at;
491			FVECT c0;
492			double s;
493			{
494			FVECT ck0;
495			int branch;
496			register AMBVAL *av;
497			register int i;
498
499			if ((av = newambval()) == NULL)
500			goto memerr;
501	greg	1.9	copystruct(av, aval);
502	greg	1.1	VCOPY(ck0, c0);
503			while (s(OCTSCALE/2) > av->radambacc) {
504			if (at->kid == NULL)
505			if ((at->kid = newambtree()) == NULL)
506			goto memerr;
507			s *= 0.5;
508			branch = 0;
509			for (i = 0; i < 3; i++)
510			if (av->pos[i] > ck0[i] + s) {
511			ck0[i] += s;
512			branch \|= 1 << i;
513			}
514			at = at->kid + branch;
515			}
516			av->next = at->alist;
517			at->alist = av;
518			return;
519			memerr:
520			error(SYSTEM, "out of memory in avinsert");
521			}