| 21 |
|
#include "ambient.h" |
| 22 |
|
#include "random.h" |
| 23 |
|
|
| 24 |
< |
#ifdef NEWAMB |
| 24 |
> |
#ifndef OLDAMB |
| 25 |
|
|
| 26 |
|
extern void SDsquare2disk(double ds[2], double seedx, double seedy); |
| 27 |
|
|
| 33 |
|
|
| 34 |
|
typedef struct { |
| 35 |
|
RAY *rp; /* originating ray sample */ |
| 36 |
– |
FVECT ux, uy; /* tangent axis unit vectors */ |
| 36 |
|
int ns; /* number of samples per axis */ |
| 37 |
+ |
int sampOK; /* acquired full sample set? */ |
| 38 |
|
COLOR acoef; /* division contribution coefficient */ |
| 39 |
+ |
double acol[3]; /* accumulated color */ |
| 40 |
+ |
FVECT ux, uy; /* tangent axis unit vectors */ |
| 41 |
|
AMBSAMP sa[1]; /* sample array (extends struct) */ |
| 42 |
|
} AMBHEMI; /* ambient sample hemisphere */ |
| 43 |
|
|
| 50 |
|
} FFTRI; /* vectors and coefficients for Hessian calculation */ |
| 51 |
|
|
| 52 |
|
|
| 53 |
< |
static AMBHEMI * |
| 54 |
< |
inithemi( /* initialize sampling hemisphere */ |
| 55 |
< |
COLOR ac, |
| 56 |
< |
RAY *r, |
| 57 |
< |
double wt |
| 53 |
> |
static int |
| 54 |
> |
ambsample( /* initial ambient division sample */ |
| 55 |
> |
AMBHEMI *hp, |
| 56 |
> |
int i, |
| 57 |
> |
int j, |
| 58 |
> |
int n |
| 59 |
|
) |
| 60 |
|
{ |
| 61 |
< |
AMBHEMI *hp; |
| 62 |
< |
double d; |
| 60 |
< |
int n, i; |
| 61 |
< |
/* set number of divisions */ |
| 62 |
< |
if (ambacc <= FTINY && |
| 63 |
< |
wt > (d = 0.8*intens(ac)*r->rweight/(ambdiv*minweight))) |
| 64 |
< |
wt = d; /* avoid ray termination */ |
| 65 |
< |
n = sqrt(ambdiv * wt) + 0.5; |
| 66 |
< |
i = 1 + 5*(ambacc > FTINY); /* minimum number of samples */ |
| 67 |
< |
if (n < i) |
| 68 |
< |
n = i; |
| 69 |
< |
/* allocate sampling array */ |
| 70 |
< |
hp = (AMBHEMI *)malloc(sizeof(AMBHEMI) + sizeof(AMBSAMP)*(n*n - 1)); |
| 71 |
< |
if (hp == NULL) |
| 72 |
< |
return(NULL); |
| 73 |
< |
hp->rp = r; |
| 74 |
< |
hp->ns = n; |
| 75 |
< |
/* assign coefficient */ |
| 76 |
< |
copycolor(hp->acoef, ac); |
| 77 |
< |
d = 1.0/(n*n); |
| 78 |
< |
scalecolor(hp->acoef, d); |
| 79 |
< |
/* make tangent plane axes */ |
| 80 |
< |
hp->uy[0] = 0.5 - frandom(); |
| 81 |
< |
hp->uy[1] = 0.5 - frandom(); |
| 82 |
< |
hp->uy[2] = 0.5 - frandom(); |
| 83 |
< |
for (i = 3; i--; ) |
| 84 |
< |
if ((-0.6 < r->ron[i]) & (r->ron[i] < 0.6)) |
| 85 |
< |
break; |
| 86 |
< |
if (i < 0) |
| 87 |
< |
error(CONSISTENCY, "bad ray direction in inithemi"); |
| 88 |
< |
hp->uy[i] = 1.0; |
| 89 |
< |
VCROSS(hp->ux, hp->uy, r->ron); |
| 90 |
< |
normalize(hp->ux); |
| 91 |
< |
VCROSS(hp->uy, r->ron, hp->ux); |
| 92 |
< |
/* we're ready to sample */ |
| 93 |
< |
return(hp); |
| 94 |
< |
} |
| 95 |
< |
|
| 96 |
< |
|
| 97 |
< |
/* Sample ambient division and apply weighting coefficient */ |
| 98 |
< |
static int |
| 99 |
< |
getambsamp(RAY *arp, AMBHEMI *hp, int i, int j, int n) |
| 100 |
< |
{ |
| 61 |
> |
AMBSAMP *ap = &ambsam(hp,i,j); |
| 62 |
> |
RAY ar; |
| 63 |
|
int hlist[3], ii; |
| 64 |
|
double spt[2], zd; |
| 65 |
+ |
/* generate hemispherical sample */ |
| 66 |
|
/* ambient coefficient for weight */ |
| 67 |
|
if (ambacc > FTINY) |
| 68 |
< |
setcolor(arp->rcoef, AVGREFL, AVGREFL, AVGREFL); |
| 68 |
> |
setcolor(ar.rcoef, AVGREFL, AVGREFL, AVGREFL); |
| 69 |
|
else |
| 70 |
< |
copycolor(arp->rcoef, hp->acoef); |
| 71 |
< |
if (rayorigin(arp, AMBIENT, hp->rp, arp->rcoef) < 0) |
| 70 |
> |
copycolor(ar.rcoef, hp->acoef); |
| 71 |
> |
if (rayorigin(&ar, AMBIENT, hp->rp, ar.rcoef) < 0) |
| 72 |
|
return(0); |
| 73 |
|
if (ambacc > FTINY) { |
| 74 |
< |
multcolor(arp->rcoef, hp->acoef); |
| 75 |
< |
scalecolor(arp->rcoef, 1./AVGREFL); |
| 74 |
> |
multcolor(ar.rcoef, hp->acoef); |
| 75 |
> |
scalecolor(ar.rcoef, 1./AVGREFL); |
| 76 |
|
} |
| 77 |
|
hlist[0] = hp->rp->rno; |
| 78 |
|
hlist[1] = j; |
| 79 |
|
hlist[2] = i; |
| 80 |
|
multisamp(spt, 2, urand(ilhash(hlist,3)+n)); |
| 81 |
< |
if (!n) { /* avoid border samples for n==0 */ |
| 81 |
> |
/* avoid coincident samples */ |
| 82 |
> |
if (!n && (0 < i) & (i < hp->ns-1) && |
| 83 |
> |
(0 < j) & (j < hp->ns-1)) { |
| 84 |
|
if ((spt[0] < 0.1) | (spt[0] >= 0.9)) |
| 85 |
|
spt[0] = 0.1 + 0.8*frandom(); |
| 86 |
|
if ((spt[1] < 0.1) | (spt[1] >= 0.9)) |
| 89 |
|
SDsquare2disk(spt, (j+spt[1])/hp->ns, (i+spt[0])/hp->ns); |
| 90 |
|
zd = sqrt(1. - spt[0]*spt[0] - spt[1]*spt[1]); |
| 91 |
|
for (ii = 3; ii--; ) |
| 92 |
< |
arp->rdir[ii] = spt[0]*hp->ux[ii] + |
| 92 |
> |
ar.rdir[ii] = spt[0]*hp->ux[ii] + |
| 93 |
|
spt[1]*hp->uy[ii] + |
| 94 |
|
zd*hp->rp->ron[ii]; |
| 95 |
< |
checknorm(arp->rdir); |
| 95 |
> |
checknorm(ar.rdir); |
| 96 |
|
dimlist[ndims++] = AI(hp,i,j) + 90171; |
| 97 |
< |
rayvalue(arp); /* evaluate ray */ |
| 98 |
< |
ndims--; /* apply coefficient */ |
| 99 |
< |
multcolor(arp->rcol, arp->rcoef); |
| 97 |
> |
rayvalue(&ar); /* evaluate ray */ |
| 98 |
> |
ndims--; |
| 99 |
> |
if (ar.rt <= FTINY) |
| 100 |
> |
return(0); /* should never happen */ |
| 101 |
> |
multcolor(ar.rcol, ar.rcoef); /* apply coefficient */ |
| 102 |
> |
if (ar.rt*ap->d < 1.0) /* new/closer distance? */ |
| 103 |
> |
ap->d = 1.0/ar.rt; |
| 104 |
> |
if (!n) { /* record first vertex & value */ |
| 105 |
> |
if (ar.rt > 10.0*thescene.cusize) |
| 106 |
> |
ar.rt = 10.0*thescene.cusize; |
| 107 |
> |
VSUM(ap->p, ar.rorg, ar.rdir, ar.rt); |
| 108 |
> |
copycolor(ap->v, ar.rcol); |
| 109 |
> |
} else { /* else update recorded value */ |
| 110 |
> |
hp->acol[RED] -= colval(ap->v,RED); |
| 111 |
> |
hp->acol[GRN] -= colval(ap->v,GRN); |
| 112 |
> |
hp->acol[BLU] -= colval(ap->v,BLU); |
| 113 |
> |
zd = 1.0/(double)(n+1); |
| 114 |
> |
scalecolor(ar.rcol, zd); |
| 115 |
> |
zd *= (double)n; |
| 116 |
> |
scalecolor(ap->v, zd); |
| 117 |
> |
addcolor(ap->v, ar.rcol); |
| 118 |
> |
} |
| 119 |
> |
addcolor(hp->acol, ap->v); /* add to our sum */ |
| 120 |
|
return(1); |
| 121 |
|
} |
| 122 |
|
|
| 123 |
|
|
| 139 |
– |
static AMBSAMP * |
| 140 |
– |
ambsample( /* initial ambient division sample */ |
| 141 |
– |
AMBHEMI *hp, |
| 142 |
– |
int i, |
| 143 |
– |
int j |
| 144 |
– |
) |
| 145 |
– |
{ |
| 146 |
– |
AMBSAMP *ap = &ambsam(hp,i,j); |
| 147 |
– |
RAY ar; |
| 148 |
– |
/* generate hemispherical sample */ |
| 149 |
– |
if (!getambsamp(&ar, hp, i, j, 0) || ar.rt <= FTINY) { |
| 150 |
– |
memset(ap, 0, sizeof(AMBSAMP)); |
| 151 |
– |
return(NULL); |
| 152 |
– |
} |
| 153 |
– |
ap->d = 1.0/ar.rt; /* limit vertex distance */ |
| 154 |
– |
if (ar.rt > 10.0*thescene.cusize) |
| 155 |
– |
ar.rt = 10.0*thescene.cusize; |
| 156 |
– |
VSUM(ap->p, ar.rorg, ar.rdir, ar.rt); |
| 157 |
– |
copycolor(ap->v, ar.rcol); |
| 158 |
– |
return(ap); |
| 159 |
– |
} |
| 160 |
– |
|
| 161 |
– |
|
| 124 |
|
/* Estimate errors based on ambient division differences */ |
| 125 |
|
static float * |
| 126 |
|
getambdiffs(AMBHEMI *hp) |
| 175 |
|
|
| 176 |
|
/* Perform super-sampling on hemisphere (introduces bias) */ |
| 177 |
|
static void |
| 178 |
< |
ambsupersamp(double acol[3], AMBHEMI *hp, int cnt) |
| 178 |
> |
ambsupersamp(AMBHEMI *hp, int cnt) |
| 179 |
|
{ |
| 180 |
|
float *earr = getambdiffs(hp); |
| 181 |
|
double e2rem = 0; |
| 182 |
|
AMBSAMP *ap; |
| 221 |
– |
RAY ar; |
| 222 |
– |
double asum[3]; |
| 183 |
|
float *ep; |
| 184 |
|
int i, j, n, nss; |
| 185 |
|
|
| 194 |
|
if (e2rem <= FTINY) |
| 195 |
|
goto done; /* nothing left to do */ |
| 196 |
|
nss = *ep/e2rem*cnt + frandom(); |
| 197 |
< |
asum[0] = asum[1] = asum[2] = 0.0; |
| 198 |
< |
for (n = 1; n <= nss; n++) { |
| 199 |
< |
if (!getambsamp(&ar, hp, i, j, n)) { |
| 240 |
< |
nss = n-1; |
| 241 |
< |
break; |
| 242 |
< |
} |
| 243 |
< |
addcolor(asum, ar.rcol); |
| 244 |
< |
} |
| 245 |
< |
if (nss) { /* update returned ambient value */ |
| 246 |
< |
const double ssf = 1./(nss + 1.); |
| 247 |
< |
for (n = 3; n--; ) |
| 248 |
< |
acol[n] += ssf*asum[n] + |
| 249 |
< |
(ssf - 1.)*colval(ap->v,n); |
| 250 |
< |
} |
| 251 |
< |
e2rem -= *ep++; /* update remainders */ |
| 252 |
< |
cnt -= nss; |
| 197 |
> |
for (n = 1; n <= nss && ambsample(hp,i,j,n); n++) |
| 198 |
> |
--cnt; |
| 199 |
> |
e2rem -= *ep++; /* update remainder */ |
| 200 |
|
} |
| 201 |
|
done: |
| 202 |
|
free(earr); |
| 203 |
|
} |
| 204 |
|
|
| 205 |
|
|
| 206 |
+ |
static AMBHEMI * |
| 207 |
+ |
samp_hemi( /* sample indirect hemisphere */ |
| 208 |
+ |
COLOR rcol, |
| 209 |
+ |
RAY *r, |
| 210 |
+ |
double wt |
| 211 |
+ |
) |
| 212 |
+ |
{ |
| 213 |
+ |
AMBHEMI *hp; |
| 214 |
+ |
double d; |
| 215 |
+ |
int n, i, j; |
| 216 |
+ |
/* set number of divisions */ |
| 217 |
+ |
if (ambacc <= FTINY && |
| 218 |
+ |
wt > (d = 0.8*intens(rcol)*r->rweight/(ambdiv*minweight))) |
| 219 |
+ |
wt = d; /* avoid ray termination */ |
| 220 |
+ |
n = sqrt(ambdiv * wt) + 0.5; |
| 221 |
+ |
i = 1 + 5*(ambacc > FTINY); /* minimum number of samples */ |
| 222 |
+ |
if (n < i) |
| 223 |
+ |
n = i; |
| 224 |
+ |
/* allocate sampling array */ |
| 225 |
+ |
hp = (AMBHEMI *)malloc(sizeof(AMBHEMI) + sizeof(AMBSAMP)*(n*n - 1)); |
| 226 |
+ |
if (hp == NULL) |
| 227 |
+ |
error(SYSTEM, "out of memory in samp_hemi"); |
| 228 |
+ |
hp->rp = r; |
| 229 |
+ |
hp->ns = n; |
| 230 |
+ |
hp->acol[RED] = hp->acol[GRN] = hp->acol[BLU] = 0.0; |
| 231 |
+ |
memset(hp->sa, 0, sizeof(AMBSAMP)*n*n); |
| 232 |
+ |
hp->sampOK = 0; |
| 233 |
+ |
/* assign coefficient */ |
| 234 |
+ |
copycolor(hp->acoef, rcol); |
| 235 |
+ |
d = 1.0/(n*n); |
| 236 |
+ |
scalecolor(hp->acoef, d); |
| 237 |
+ |
/* make tangent plane axes */ |
| 238 |
+ |
if (!getperpendicular(hp->ux, r->ron, 1)) |
| 239 |
+ |
error(CONSISTENCY, "bad ray direction in samp_hemi"); |
| 240 |
+ |
VCROSS(hp->uy, r->ron, hp->ux); |
| 241 |
+ |
/* sample divisions */ |
| 242 |
+ |
for (i = hp->ns; i--; ) |
| 243 |
+ |
for (j = hp->ns; j--; ) |
| 244 |
+ |
hp->sampOK += ambsample(hp, i, j, 0); |
| 245 |
+ |
copycolor(rcol, hp->acol); |
| 246 |
+ |
if (!hp->sampOK) { /* utter failure? */ |
| 247 |
+ |
free(hp); |
| 248 |
+ |
return(NULL); |
| 249 |
+ |
} |
| 250 |
+ |
if (hp->sampOK < hp->ns*hp->ns) { |
| 251 |
+ |
hp->sampOK *= -1; /* soft failure */ |
| 252 |
+ |
return(hp); |
| 253 |
+ |
} |
| 254 |
+ |
n = ambssamp*wt + 0.5; |
| 255 |
+ |
if (n > 8) { /* perform super-sampling? */ |
| 256 |
+ |
ambsupersamp(hp, n); |
| 257 |
+ |
copycolor(rcol, hp->acol); |
| 258 |
+ |
} |
| 259 |
+ |
return(hp); /* all is well */ |
| 260 |
+ |
} |
| 261 |
+ |
|
| 262 |
+ |
|
| 263 |
|
/* Return brightness of farthest ambient sample */ |
| 264 |
|
static double |
| 265 |
|
back_ambval(AMBHEMI *hp, const int n1, const int n2, const int n3) |
| 600 |
|
ambcorral(AMBHEMI *hp, FVECT uv[2], const double r0, const double r1) |
| 601 |
|
{ |
| 602 |
|
const double max_d = 1.0/(minarad*ambacc + 0.001); |
| 603 |
< |
const double ang_res = 0.5*PI/(hp->ns-1); |
| 604 |
< |
const double ang_step = ang_res/((int)(16/PI*ang_res) + (1+FTINY)); |
| 603 |
> |
const double ang_res = 0.5*PI/hp->ns; |
| 604 |
> |
const double ang_step = ang_res/((int)(16/PI*ang_res) + 1.01); |
| 605 |
|
double avg_d = 0; |
| 606 |
|
uint32 flgs = 0; |
| 607 |
+ |
FVECT vec; |
| 608 |
+ |
double u, v; |
| 609 |
+ |
double ang, a1; |
| 610 |
|
int i, j; |
| 611 |
|
/* don't bother for a few samples */ |
| 612 |
|
if (hp->ns < 12) |
| 624 |
|
for (i = 0; i < hp->ns; i++) |
| 625 |
|
for (j = 0; j < hp->ns; j += !i|(i==hp->ns-1) ? 1 : hp->ns-1) { |
| 626 |
|
AMBSAMP *ap = &ambsam(hp,i,j); |
| 620 |
– |
FVECT vec; |
| 621 |
– |
double u, v; |
| 622 |
– |
double ang, a1; |
| 623 |
– |
int abp; |
| 627 |
|
if ((ap->d <= FTINY) | (ap->d >= max_d)) |
| 628 |
|
continue; /* too far or too near */ |
| 629 |
|
VSUB(vec, ap->p, hp->rp->rop); |
| 630 |
< |
u = DOT(vec, uv[0]) * ap->d; |
| 631 |
< |
v = DOT(vec, uv[1]) * ap->d; |
| 632 |
< |
if ((r0*r0*u*u + r1*r1*v*v) * ap->d*ap->d <= 1.0) |
| 630 |
> |
u = DOT(vec, uv[0]); |
| 631 |
> |
v = DOT(vec, uv[1]); |
| 632 |
> |
if ((r0*r0*u*u + r1*r1*v*v) * ap->d*ap->d <= u*u + v*v) |
| 633 |
|
continue; /* occluder outside ellipse */ |
| 634 |
|
ang = atan2a(v, u); /* else set direction flags */ |
| 635 |
< |
for (a1 = ang-.5*ang_res; a1 <= ang+.5*ang_res; a1 += ang_step) |
| 635 |
> |
for (a1 = ang-ang_res; a1 <= ang+ang_res; a1 += ang_step) |
| 636 |
|
flgs |= 1L<<(int)(16/PI*(a1 + 2.*PI*(a1 < 0))); |
| 637 |
|
} |
| 638 |
+ |
/* add low-angle incident (< 20deg) */ |
| 639 |
+ |
if (fabs(hp->rp->rod) <= 0.342) { |
| 640 |
+ |
u = -DOT(hp->rp->rdir, uv[0]); |
| 641 |
+ |
v = -DOT(hp->rp->rdir, uv[1]); |
| 642 |
+ |
if ((r0*r0*u*u + r1*r1*v*v) > hp->rp->rot*hp->rp->rot) { |
| 643 |
+ |
ang = atan2a(v, u); |
| 644 |
+ |
ang += 2.*PI*(ang < 0); |
| 645 |
+ |
ang *= 16/PI; |
| 646 |
+ |
if ((ang < .5) | (ang >= 31.5)) |
| 647 |
+ |
flgs |= 0x80000001; |
| 648 |
+ |
else |
| 649 |
+ |
flgs |= 3L<<(int)(ang-.5); |
| 650 |
+ |
} |
| 651 |
+ |
} |
| 652 |
|
return(flgs); |
| 653 |
|
} |
| 654 |
|
|
| 665 |
|
uint32 *crlp /* returned (optional) */ |
| 666 |
|
) |
| 667 |
|
{ |
| 668 |
< |
AMBHEMI *hp = inithemi(rcol, r, wt); |
| 652 |
< |
int cnt; |
| 668 |
> |
AMBHEMI *hp = samp_hemi(rcol, r, wt); |
| 669 |
|
FVECT my_uv[2]; |
| 670 |
< |
double d, K, acol[3]; |
| 670 |
> |
double d, K; |
| 671 |
|
AMBSAMP *ap; |
| 672 |
< |
int i, j; |
| 673 |
< |
/* check/initialize */ |
| 658 |
< |
if (hp == NULL) |
| 659 |
< |
return(0); |
| 672 |
> |
int i; |
| 673 |
> |
/* clear return values */ |
| 674 |
|
if (uv != NULL) |
| 675 |
|
memset(uv, 0, sizeof(FVECT)*2); |
| 676 |
|
if (ra != NULL) |
| 681 |
|
dg[0] = dg[1] = 0.0; |
| 682 |
|
if (crlp != NULL) |
| 683 |
|
*crlp = 0; |
| 684 |
< |
/* sample the hemisphere */ |
| 685 |
< |
acol[0] = acol[1] = acol[2] = 0.0; |
| 686 |
< |
cnt = 0; |
| 687 |
< |
for (i = hp->ns; i--; ) |
| 688 |
< |
for (j = hp->ns; j--; ) |
| 689 |
< |
if ((ap = ambsample(hp, i, j)) != NULL) { |
| 690 |
< |
addcolor(acol, ap->v); |
| 677 |
< |
++cnt; |
| 678 |
< |
} |
| 679 |
< |
if (!cnt) { |
| 680 |
< |
setcolor(rcol, 0.0, 0.0, 0.0); |
| 681 |
< |
free(hp); |
| 682 |
< |
return(0); /* no valid samples */ |
| 684 |
> |
if (hp == NULL) /* sampling falure? */ |
| 685 |
> |
return(0); |
| 686 |
> |
|
| 687 |
> |
if ((ra == NULL) & (pg == NULL) & (dg == NULL) || |
| 688 |
> |
(hp->sampOK < 0) | (hp->ns < 6)) { |
| 689 |
> |
free(hp); /* Hessian not requested/possible */ |
| 690 |
> |
return(-1); /* value-only return value */ |
| 691 |
|
} |
| 692 |
< |
if (cnt < hp->ns*hp->ns) { /* incomplete sampling? */ |
| 685 |
< |
copycolor(rcol, acol); |
| 686 |
< |
free(hp); |
| 687 |
< |
return(-1); /* return value w/o Hessian */ |
| 688 |
< |
} |
| 689 |
< |
cnt = ambssamp*wt + 0.5; /* perform super-sampling? */ |
| 690 |
< |
if (cnt > 8) |
| 691 |
< |
ambsupersamp(acol, hp, cnt); |
| 692 |
< |
copycolor(rcol, acol); /* final indirect irradiance/PI */ |
| 693 |
< |
if ((ra == NULL) & (pg == NULL) & (dg == NULL)) { |
| 694 |
< |
free(hp); |
| 695 |
< |
return(-1); /* no radius or gradient calc. */ |
| 696 |
< |
} |
| 697 |
< |
if ((d = bright(acol)) > FTINY) { /* normalize Y values */ |
| 692 |
> |
if ((d = bright(rcol)) > FTINY) { /* normalize Y values */ |
| 693 |
|
d = 0.99*(hp->ns*hp->ns)/d; |
| 694 |
|
K = 0.01; |
| 695 |
|
} else { /* or fall back on geometric Hessian */ |
| 724 |
|
if (ra[1] < minarad) |
| 725 |
|
ra[1] = minarad; |
| 726 |
|
} |
| 727 |
< |
ra[0] *= d = 1.0/sqrt(sqrt(wt)); |
| 727 |
> |
ra[0] *= d = 1.0/sqrt(wt); |
| 728 |
|
if ((ra[1] *= d) > 2.0*ra[0]) |
| 729 |
|
ra[1] = 2.0*ra[0]; |
| 730 |
|
if (ra[1] > maxarad) { |
