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; |
87 |
|
SDsquare2disk(spt, (j+spt[1])/hp->ns, (i+spt[0])/hp->ns); |
88 |
|
zd = sqrt(1. - spt[0]*spt[0] - spt[1]*spt[1]); |
89 |
|
for (ii = 3; ii--; ) |
90 |
< |
arp->rdir[ii] = spt[0]*hp->ux[ii] + |
90 |
> |
ar.rdir[ii] = spt[0]*hp->ux[ii] + |
91 |
|
spt[1]*hp->uy[ii] + |
92 |
|
zd*hp->rp->ron[ii]; |
93 |
< |
checknorm(arp->rdir); |
93 |
> |
checknorm(ar.rdir); |
94 |
|
dimlist[ndims++] = AI(hp,i,j) + 90171; |
95 |
< |
rayvalue(arp); /* evaluate ray */ |
96 |
< |
ndims--; /* apply coefficient */ |
97 |
< |
multcolor(arp->rcol, arp->rcoef); |
95 |
> |
rayvalue(&ar); /* evaluate ray */ |
96 |
> |
ndims--; |
97 |
> |
if (ar.rt <= FTINY) |
98 |
> |
return(0); /* should never happen */ |
99 |
> |
multcolor(ar.rcol, ar.rcoef); /* apply coefficient */ |
100 |
> |
if (ar.rt*ap->d < 1.0) /* new/closer distance? */ |
101 |
> |
ap->d = 1.0/ar.rt; |
102 |
> |
if (!n) { /* record first vertex & value */ |
103 |
> |
if (ar.rt > 10.0*thescene.cusize) |
104 |
> |
ar.rt = 10.0*thescene.cusize; |
105 |
> |
VSUM(ap->p, ar.rorg, ar.rdir, ar.rt); |
106 |
> |
copycolor(ap->v, ar.rcol); |
107 |
> |
} else { /* else update recorded value */ |
108 |
> |
hp->acol[RED] -= colval(ap->v,RED); |
109 |
> |
hp->acol[GRN] -= colval(ap->v,GRN); |
110 |
> |
hp->acol[BLU] -= colval(ap->v,BLU); |
111 |
> |
zd = 1.0/(double)(n+1); |
112 |
> |
scalecolor(ar.rcol, zd); |
113 |
> |
zd *= (double)n; |
114 |
> |
scalecolor(ap->v, zd); |
115 |
> |
addcolor(ap->v, ar.rcol); |
116 |
> |
} |
117 |
> |
addcolor(hp->acol, ap->v); /* add to our sum */ |
118 |
|
return(1); |
119 |
|
} |
120 |
|
|
121 |
|
|
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 |
– |
|
122 |
|
/* Estimate errors based on ambient division differences */ |
123 |
|
static float * |
124 |
|
getambdiffs(AMBHEMI *hp) |
173 |
|
|
174 |
|
/* Perform super-sampling on hemisphere (introduces bias) */ |
175 |
|
static void |
176 |
< |
ambsupersamp(double acol[3], AMBHEMI *hp, int cnt) |
176 |
> |
ambsupersamp(AMBHEMI *hp, int cnt) |
177 |
|
{ |
178 |
|
float *earr = getambdiffs(hp); |
179 |
|
double e2rem = 0; |
180 |
|
AMBSAMP *ap; |
181 |
|
RAY ar; |
222 |
– |
double asum[3]; |
182 |
|
float *ep; |
183 |
|
int i, j, n, nss; |
184 |
|
|
193 |
|
if (e2rem <= FTINY) |
194 |
|
goto done; /* nothing left to do */ |
195 |
|
nss = *ep/e2rem*cnt + frandom(); |
196 |
< |
asum[0] = asum[1] = asum[2] = 0.0; |
197 |
< |
for (n = 1; n <= nss; n++) { |
198 |
< |
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; |
196 |
> |
for (n = 1; n <= nss && ambsample(hp,i,j,n); n++) |
197 |
> |
--cnt; |
198 |
> |
e2rem -= *ep++; /* update remainder */ |
199 |
|
} |
200 |
|
done: |
201 |
|
free(earr); |
202 |
|
} |
203 |
|
|
204 |
|
|
205 |
+ |
static AMBHEMI * |
206 |
+ |
samp_hemi( /* sample indirect hemisphere */ |
207 |
+ |
COLOR rcol, |
208 |
+ |
RAY *r, |
209 |
+ |
double wt |
210 |
+ |
) |
211 |
+ |
{ |
212 |
+ |
AMBHEMI *hp; |
213 |
+ |
double d; |
214 |
+ |
int n, i, j; |
215 |
+ |
/* set number of divisions */ |
216 |
+ |
if (ambacc <= FTINY && |
217 |
+ |
wt > (d = 0.8*intens(rcol)*r->rweight/(ambdiv*minweight))) |
218 |
+ |
wt = d; /* avoid ray termination */ |
219 |
+ |
n = sqrt(ambdiv * wt) + 0.5; |
220 |
+ |
i = 1 + 5*(ambacc > FTINY); /* minimum number of samples */ |
221 |
+ |
if (n < i) |
222 |
+ |
n = i; |
223 |
+ |
/* allocate sampling array */ |
224 |
+ |
hp = (AMBHEMI *)malloc(sizeof(AMBHEMI) + sizeof(AMBSAMP)*(n*n - 1)); |
225 |
+ |
if (hp == NULL) |
226 |
+ |
error(SYSTEM, "out of memory in samp_hemi"); |
227 |
+ |
hp->rp = r; |
228 |
+ |
hp->ns = n; |
229 |
+ |
hp->acol[RED] = hp->acol[GRN] = hp->acol[BLU] = 0.0; |
230 |
+ |
memset(hp->sa, 0, sizeof(AMBSAMP)*n*n); |
231 |
+ |
hp->sampOK = 0; |
232 |
+ |
/* assign coefficient */ |
233 |
+ |
copycolor(hp->acoef, rcol); |
234 |
+ |
d = 1.0/(n*n); |
235 |
+ |
scalecolor(hp->acoef, d); |
236 |
+ |
/* make tangent plane axes */ |
237 |
+ |
hp->uy[0] = 0.5 - frandom(); |
238 |
+ |
hp->uy[1] = 0.5 - frandom(); |
239 |
+ |
hp->uy[2] = 0.5 - frandom(); |
240 |
+ |
for (i = 3; i--; ) |
241 |
+ |
if ((-0.6 < r->ron[i]) & (r->ron[i] < 0.6)) |
242 |
+ |
break; |
243 |
+ |
if (i < 0) |
244 |
+ |
error(CONSISTENCY, "bad ray direction in samp_hemi"); |
245 |
+ |
hp->uy[i] = 1.0; |
246 |
+ |
VCROSS(hp->ux, hp->uy, r->ron); |
247 |
+ |
normalize(hp->ux); |
248 |
+ |
VCROSS(hp->uy, r->ron, hp->ux); |
249 |
+ |
/* sample divisions */ |
250 |
+ |
for (i = hp->ns; i--; ) |
251 |
+ |
for (j = hp->ns; j--; ) |
252 |
+ |
hp->sampOK += ambsample(hp, i, j, 0); |
253 |
+ |
copycolor(rcol, hp->acol); |
254 |
+ |
if (!hp->sampOK) { /* utter failure? */ |
255 |
+ |
free(hp); |
256 |
+ |
return(NULL); |
257 |
+ |
} |
258 |
+ |
if (hp->sampOK < hp->ns*hp->ns) { |
259 |
+ |
hp->sampOK *= -1; /* soft failure */ |
260 |
+ |
return(hp); |
261 |
+ |
} |
262 |
+ |
n = ambssamp*wt + 0.5; |
263 |
+ |
if (n > 8) { /* perform super-sampling? */ |
264 |
+ |
ambsupersamp(hp, n); |
265 |
+ |
copycolor(rcol, hp->acol); |
266 |
+ |
} |
267 |
+ |
return(hp); /* all is well */ |
268 |
+ |
} |
269 |
+ |
|
270 |
+ |
|
271 |
|
/* Return brightness of farthest ambient sample */ |
272 |
|
static double |
273 |
|
back_ambval(AMBHEMI *hp, const int n1, const int n2, const int n3) |
612 |
|
const double ang_step = ang_res/((int)(16/PI*ang_res) + (1+FTINY)); |
613 |
|
double avg_d = 0; |
614 |
|
uint32 flgs = 0; |
615 |
+ |
FVECT vec; |
616 |
+ |
double u, v; |
617 |
+ |
double ang, a1; |
618 |
|
int i, j; |
619 |
|
/* don't bother for a few samples */ |
620 |
|
if (hp->ns < 12) |
632 |
|
for (i = 0; i < hp->ns; i++) |
633 |
|
for (j = 0; j < hp->ns; j += !i|(i==hp->ns-1) ? 1 : hp->ns-1) { |
634 |
|
AMBSAMP *ap = &ambsam(hp,i,j); |
620 |
– |
FVECT vec; |
621 |
– |
double u, v; |
622 |
– |
double ang, a1; |
623 |
– |
int abp; |
635 |
|
if ((ap->d <= FTINY) | (ap->d >= max_d)) |
636 |
|
continue; /* too far or too near */ |
637 |
|
VSUB(vec, ap->p, hp->rp->rop); |
638 |
< |
u = DOT(vec, uv[0]) * ap->d; |
639 |
< |
v = DOT(vec, uv[1]) * ap->d; |
640 |
< |
if ((r0*r0*u*u + r1*r1*v*v) * ap->d*ap->d <= 1.0) |
638 |
> |
u = DOT(vec, uv[0]); |
639 |
> |
v = DOT(vec, uv[1]); |
640 |
> |
if ((r0*r0*u*u + r1*r1*v*v) * ap->d*ap->d <= u*u + v*v) |
641 |
|
continue; /* occluder outside ellipse */ |
642 |
|
ang = atan2a(v, u); /* else set direction flags */ |
643 |
|
for (a1 = ang-.5*ang_res; a1 <= ang+.5*ang_res; a1 += ang_step) |
644 |
|
flgs |= 1L<<(int)(16/PI*(a1 + 2.*PI*(a1 < 0))); |
645 |
|
} |
646 |
+ |
/* add low-angle incident (< 20deg) */ |
647 |
+ |
if (fabs(hp->rp->rod) <= 0.342) { |
648 |
+ |
u = -DOT(hp->rp->rdir, uv[0]); |
649 |
+ |
v = -DOT(hp->rp->rdir, uv[1]); |
650 |
+ |
if ((r0*r0*u*u + r1*r1*v*v) > hp->rp->rot*hp->rp->rot) { |
651 |
+ |
ang = atan2a(v, u); |
652 |
+ |
ang += 2.*PI*(ang < 0); |
653 |
+ |
ang *= 16/PI; |
654 |
+ |
if ((ang < .5) | (ang >= 31.5)) |
655 |
+ |
flgs |= 0x80000001; |
656 |
+ |
else |
657 |
+ |
flgs |= 3L<<(int)(ang-.5); |
658 |
+ |
} |
659 |
+ |
} |
660 |
|
return(flgs); |
661 |
|
} |
662 |
|
|
673 |
|
uint32 *crlp /* returned (optional) */ |
674 |
|
) |
675 |
|
{ |
676 |
< |
AMBHEMI *hp = inithemi(rcol, r, wt); |
652 |
< |
int cnt; |
676 |
> |
AMBHEMI *hp = samp_hemi(rcol, r, wt); |
677 |
|
FVECT my_uv[2]; |
678 |
< |
double d, K, acol[3]; |
678 |
> |
double d, K; |
679 |
|
AMBSAMP *ap; |
680 |
< |
int i, j; |
681 |
< |
/* check/initialize */ |
658 |
< |
if (hp == NULL) |
659 |
< |
return(0); |
680 |
> |
int i; |
681 |
> |
/* clear return values */ |
682 |
|
if (uv != NULL) |
683 |
|
memset(uv, 0, sizeof(FVECT)*2); |
684 |
|
if (ra != NULL) |
689 |
|
dg[0] = dg[1] = 0.0; |
690 |
|
if (crlp != NULL) |
691 |
|
*crlp = 0; |
692 |
< |
/* sample the hemisphere */ |
693 |
< |
acol[0] = acol[1] = acol[2] = 0.0; |
694 |
< |
cnt = 0; |
695 |
< |
for (i = hp->ns; i--; ) |
696 |
< |
for (j = hp->ns; j--; ) |
697 |
< |
if ((ap = ambsample(hp, i, j)) != NULL) { |
698 |
< |
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 */ |
692 |
> |
if (hp == NULL) /* sampling falure? */ |
693 |
> |
return(0); |
694 |
> |
|
695 |
> |
if ((ra == NULL) & (pg == NULL) & (dg == NULL) || |
696 |
> |
(hp->sampOK < 0) | (hp->ns < 4)) { |
697 |
> |
free(hp); /* Hessian not requested/possible */ |
698 |
> |
return(-1); /* value-only return value */ |
699 |
|
} |
700 |
< |
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 */ |
700 |
> |
if ((d = bright(rcol)) > FTINY) { /* normalize Y values */ |
701 |
|
d = 0.99*(hp->ns*hp->ns)/d; |
702 |
|
K = 0.01; |
703 |
|
} else { /* or fall back on geometric Hessian */ |
732 |
|
if (ra[1] < minarad) |
733 |
|
ra[1] = minarad; |
734 |
|
} |
735 |
< |
ra[0] *= d = 1.0/sqrt(sqrt(wt)); |
735 |
> |
ra[0] *= d = 1.0/sqrt(wt); |
736 |
|
if ((ra[1] *= d) > 2.0*ra[0]) |
737 |
|
ra[1] = 2.0*ra[0]; |
738 |
|
if (ra[1] > maxarad) { |