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