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 |
|
|
28 |
|
typedef struct { |
29 |
|
COLOR v; /* hemisphere sample value */ |
30 |
< |
float d; /* reciprocal distance (1/rt) */ |
30 |
> |
float d; /* reciprocal distance */ |
31 |
|
FVECT p; /* intersection point */ |
32 |
|
} AMBSAMP; /* sample value */ |
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 |
> |
ambcollision( /* proposed direciton collides? */ |
55 |
> |
AMBHEMI *hp, |
56 |
> |
int i, |
57 |
> |
int j, |
58 |
> |
FVECT dv |
59 |
|
) |
60 |
|
{ |
61 |
< |
AMBHEMI *hp; |
62 |
< |
double d; |
63 |
< |
int n, i; |
64 |
< |
/* set number of divisions */ |
65 |
< |
if (ambacc <= FTINY && |
66 |
< |
wt > (d = 0.8*intens(ac)*r->rweight/(ambdiv*minweight))) |
67 |
< |
wt = d; /* avoid ray termination */ |
68 |
< |
n = sqrt(ambdiv * wt) + 0.5; |
69 |
< |
i = 1 + 5*(ambacc > FTINY); /* minimum number of samples */ |
70 |
< |
if (n < i) |
71 |
< |
n = i; |
72 |
< |
/* allocate sampling array */ |
73 |
< |
hp = (AMBHEMI *)malloc(sizeof(AMBHEMI) + sizeof(AMBSAMP)*(n*n - 1)); |
74 |
< |
if (hp == NULL) |
75 |
< |
return(NULL); |
76 |
< |
hp->rp = r; |
77 |
< |
hp->ns = n; |
78 |
< |
/* assign coefficient */ |
79 |
< |
copycolor(hp->acoef, ac); |
80 |
< |
d = 1.0/(n*n); |
81 |
< |
scalecolor(hp->acoef, d); |
82 |
< |
/* make tangent plane axes */ |
83 |
< |
hp->uy[0] = 0.5 - frandom(); |
84 |
< |
hp->uy[1] = 0.5 - frandom(); |
85 |
< |
hp->uy[2] = 0.5 - frandom(); |
86 |
< |
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); |
61 |
> |
double cos_thresh; |
62 |
> |
int ii, jj; |
63 |
> |
/* min. spacing = 1/4th division */ |
64 |
> |
cos_thresh = (PI/4.)/(double)hp->ns; |
65 |
> |
cos_thresh = 1. - .5*cos_thresh*cos_thresh; |
66 |
> |
/* check existing neighbors */ |
67 |
> |
for (ii = i-1; ii <= i+1; ii++) { |
68 |
> |
if (ii < 0) continue; |
69 |
> |
if (ii >= hp->ns) break; |
70 |
> |
for (jj = j-1; jj <= j+1; jj++) { |
71 |
> |
AMBSAMP *ap; |
72 |
> |
FVECT avec; |
73 |
> |
double dprod; |
74 |
> |
if (jj < 0) continue; |
75 |
> |
if (jj >= hp->ns) break; |
76 |
> |
if ((ii==i) & (jj==j)) continue; |
77 |
> |
ap = &ambsam(hp,ii,jj); |
78 |
> |
if (ap->d <= .5/FHUGE) |
79 |
> |
continue; /* no one home */ |
80 |
> |
VSUB(avec, ap->p, hp->rp->rop); |
81 |
> |
dprod = DOT(avec, dv); |
82 |
> |
if (dprod >= cos_thresh*VLEN(avec)) |
83 |
> |
return(1); /* collision */ |
84 |
> |
} |
85 |
> |
} |
86 |
> |
return(0); /* nothing to worry about */ |
87 |
|
} |
88 |
|
|
89 |
|
|
97 |
– |
/* Sample ambient division and apply weighting coefficient */ |
90 |
|
static int |
91 |
< |
getambsamp(RAY *arp, AMBHEMI *hp, int i, int j, int n) |
91 |
> |
ambsample( /* initial ambient division sample */ |
92 |
> |
AMBHEMI *hp, |
93 |
> |
int i, |
94 |
> |
int j, |
95 |
> |
int n |
96 |
> |
) |
97 |
|
{ |
98 |
+ |
AMBSAMP *ap = &ambsam(hp,i,j); |
99 |
+ |
RAY ar; |
100 |
|
int hlist[3], ii; |
101 |
|
double spt[2], zd; |
102 |
+ |
/* generate hemispherical sample */ |
103 |
|
/* ambient coefficient for weight */ |
104 |
|
if (ambacc > FTINY) |
105 |
< |
setcolor(arp->rcoef, AVGREFL, AVGREFL, AVGREFL); |
105 |
> |
setcolor(ar.rcoef, AVGREFL, AVGREFL, AVGREFL); |
106 |
|
else |
107 |
< |
copycolor(arp->rcoef, hp->acoef); |
108 |
< |
if (rayorigin(arp, AMBIENT, hp->rp, arp->rcoef) < 0) |
107 |
> |
copycolor(ar.rcoef, hp->acoef); |
108 |
> |
if (rayorigin(&ar, AMBIENT, hp->rp, ar.rcoef) < 0) |
109 |
|
return(0); |
110 |
|
if (ambacc > FTINY) { |
111 |
< |
multcolor(arp->rcoef, hp->acoef); |
112 |
< |
scalecolor(arp->rcoef, 1./AVGREFL); |
111 |
> |
multcolor(ar.rcoef, hp->acoef); |
112 |
> |
scalecolor(ar.rcoef, 1./AVGREFL); |
113 |
|
} |
114 |
|
hlist[0] = hp->rp->rno; |
115 |
|
hlist[1] = j; |
116 |
|
hlist[2] = i; |
117 |
|
multisamp(spt, 2, urand(ilhash(hlist,3)+n)); |
118 |
< |
if (!n) { /* avoid border samples for n==0 */ |
119 |
< |
if ((spt[0] < 0.1) | (spt[0] >= 0.9)) |
120 |
< |
spt[0] = 0.1 + 0.8*frandom(); |
121 |
< |
if ((spt[1] < 0.1) | (spt[1] >= 0.9)) |
122 |
< |
spt[1] = 0.1 + 0.8*frandom(); |
123 |
< |
} |
118 |
> |
resample: |
119 |
|
SDsquare2disk(spt, (j+spt[1])/hp->ns, (i+spt[0])/hp->ns); |
120 |
|
zd = sqrt(1. - spt[0]*spt[0] - spt[1]*spt[1]); |
121 |
|
for (ii = 3; ii--; ) |
122 |
< |
arp->rdir[ii] = spt[0]*hp->ux[ii] + |
122 |
> |
ar.rdir[ii] = spt[0]*hp->ux[ii] + |
123 |
|
spt[1]*hp->uy[ii] + |
124 |
|
zd*hp->rp->ron[ii]; |
125 |
< |
checknorm(arp->rdir); |
125 |
> |
checknorm(ar.rdir); |
126 |
> |
/* avoid coincident samples */ |
127 |
> |
if (!n && ambcollision(hp, i, j, ar.rdir)) { |
128 |
> |
spt[0] = frandom(); spt[1] = frandom(); |
129 |
> |
goto resample; /* reject this sample */ |
130 |
> |
} |
131 |
|
dimlist[ndims++] = AI(hp,i,j) + 90171; |
132 |
< |
rayvalue(arp); /* evaluate ray */ |
133 |
< |
ndims--; /* apply coefficient */ |
134 |
< |
multcolor(arp->rcol, arp->rcoef); |
132 |
> |
rayvalue(&ar); /* evaluate ray */ |
133 |
> |
ndims--; |
134 |
> |
zd = raydistance(&ar); |
135 |
> |
if (zd <= FTINY) |
136 |
> |
return(0); /* should never happen */ |
137 |
> |
multcolor(ar.rcol, ar.rcoef); /* apply coefficient */ |
138 |
> |
if (zd*ap->d < 1.0) /* new/closer distance? */ |
139 |
> |
ap->d = 1.0/zd; |
140 |
> |
if (!n) { /* record first vertex & value */ |
141 |
> |
if (zd > 10.0*thescene.cusize + 1000.) |
142 |
> |
zd = 10.0*thescene.cusize + 1000.; |
143 |
> |
VSUM(ap->p, ar.rorg, ar.rdir, zd); |
144 |
> |
copycolor(ap->v, ar.rcol); |
145 |
> |
} else { /* else update recorded value */ |
146 |
> |
hp->acol[RED] -= colval(ap->v,RED); |
147 |
> |
hp->acol[GRN] -= colval(ap->v,GRN); |
148 |
> |
hp->acol[BLU] -= colval(ap->v,BLU); |
149 |
> |
zd = 1.0/(double)(n+1); |
150 |
> |
scalecolor(ar.rcol, zd); |
151 |
> |
zd *= (double)n; |
152 |
> |
scalecolor(ap->v, zd); |
153 |
> |
addcolor(ap->v, ar.rcol); |
154 |
> |
} |
155 |
> |
addcolor(hp->acol, ap->v); /* add to our sum */ |
156 |
|
return(1); |
157 |
|
} |
158 |
|
|
159 |
|
|
160 |
< |
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 |
< |
|
162 |
< |
/* Estimate errors based on ambient division differences */ |
160 |
> |
/* Estimate variance based on ambient division differences */ |
161 |
|
static float * |
162 |
|
getambdiffs(AMBHEMI *hp) |
163 |
|
{ |
164 |
+ |
const double normf = 1./bright(hp->acoef); |
165 |
|
float *earr = (float *)calloc(hp->ns*hp->ns, sizeof(float)); |
166 |
|
float *ep; |
167 |
|
AMBSAMP *ap; |
168 |
< |
double b, d2; |
168 |
> |
double b, b1, d2; |
169 |
|
int i, j; |
170 |
|
|
171 |
|
if (earr == NULL) /* out of memory? */ |
172 |
|
return(NULL); |
173 |
< |
/* compute squared neighbor diffs */ |
173 |
> |
/* sum squared neighbor diffs */ |
174 |
|
for (ap = hp->sa, ep = earr, i = 0; i < hp->ns; i++) |
175 |
|
for (j = 0; j < hp->ns; j++, ap++, ep++) { |
176 |
|
b = bright(ap[0].v); |
177 |
|
if (i) { /* from above */ |
178 |
< |
d2 = b - bright(ap[-hp->ns].v); |
179 |
< |
d2 *= d2; |
178 |
> |
b1 = bright(ap[-hp->ns].v); |
179 |
> |
d2 = b - b1; |
180 |
> |
d2 *= d2*normf/(b + b1); |
181 |
|
ep[0] += d2; |
182 |
|
ep[-hp->ns] += d2; |
183 |
|
} |
184 |
|
if (!j) continue; |
185 |
|
/* from behind */ |
186 |
< |
d2 = b - bright(ap[-1].v); |
187 |
< |
d2 *= d2; |
186 |
> |
b1 = bright(ap[-1].v); |
187 |
> |
d2 = b - b1; |
188 |
> |
d2 *= d2*normf/(b + b1); |
189 |
|
ep[0] += d2; |
190 |
|
ep[-1] += d2; |
191 |
|
if (!i) continue; |
192 |
|
/* diagonal */ |
193 |
< |
d2 = b - bright(ap[-hp->ns-1].v); |
194 |
< |
d2 *= d2; |
193 |
> |
b1 = bright(ap[-hp->ns-1].v); |
194 |
> |
d2 = b - b1; |
195 |
> |
d2 *= d2*normf/(b + b1); |
196 |
|
ep[0] += d2; |
197 |
|
ep[-hp->ns-1] += d2; |
198 |
|
} |
215 |
|
|
216 |
|
/* Perform super-sampling on hemisphere (introduces bias) */ |
217 |
|
static void |
218 |
< |
ambsupersamp(double acol[3], AMBHEMI *hp, int cnt) |
218 |
> |
ambsupersamp(AMBHEMI *hp, int cnt) |
219 |
|
{ |
220 |
|
float *earr = getambdiffs(hp); |
221 |
|
double e2rem = 0; |
220 |
– |
AMBSAMP *ap; |
221 |
– |
RAY ar; |
222 |
– |
double asum[3]; |
222 |
|
float *ep; |
223 |
|
int i, j, n, nss; |
224 |
|
|
228 |
|
for (ep = earr + hp->ns*hp->ns; ep > earr; ) |
229 |
|
e2rem += *--ep; |
230 |
|
ep = earr; /* perform super-sampling */ |
231 |
< |
for (ap = hp->sa, i = 0; i < hp->ns; i++) |
232 |
< |
for (j = 0; j < hp->ns; j++, ap++) { |
231 |
> |
for (i = 0; i < hp->ns; i++) |
232 |
> |
for (j = 0; j < hp->ns; j++) { |
233 |
|
if (e2rem <= FTINY) |
234 |
|
goto done; /* nothing left to do */ |
235 |
|
nss = *ep/e2rem*cnt + frandom(); |
236 |
< |
asum[0] = asum[1] = asum[2] = 0.0; |
237 |
< |
for (n = 1; n <= nss; n++) { |
238 |
< |
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; |
236 |
> |
for (n = 1; n <= nss && ambsample(hp,i,j,n); n++) |
237 |
> |
if (!--cnt) goto done; |
238 |
> |
e2rem -= *ep++; /* update remainder */ |
239 |
|
} |
240 |
|
done: |
241 |
|
free(earr); |
242 |
|
} |
243 |
|
|
244 |
|
|
245 |
+ |
static AMBHEMI * |
246 |
+ |
samp_hemi( /* sample indirect hemisphere */ |
247 |
+ |
COLOR rcol, |
248 |
+ |
RAY *r, |
249 |
+ |
double wt |
250 |
+ |
) |
251 |
+ |
{ |
252 |
+ |
AMBHEMI *hp; |
253 |
+ |
double d; |
254 |
+ |
int n, i, j; |
255 |
+ |
/* insignificance check */ |
256 |
+ |
if (bright(rcol) <= FTINY) |
257 |
+ |
return(NULL); |
258 |
+ |
/* set number of divisions */ |
259 |
+ |
if (ambacc <= FTINY && |
260 |
+ |
wt > (d = 0.8*intens(rcol)*r->rweight/(ambdiv*minweight))) |
261 |
+ |
wt = d; /* avoid ray termination */ |
262 |
+ |
n = sqrt(ambdiv * wt) + 0.5; |
263 |
+ |
i = 1 + 5*(ambacc > FTINY); /* minimum number of samples */ |
264 |
+ |
if (n < i) |
265 |
+ |
n = i; |
266 |
+ |
/* allocate sampling array */ |
267 |
+ |
hp = (AMBHEMI *)malloc(sizeof(AMBHEMI) + sizeof(AMBSAMP)*(n*n - 1)); |
268 |
+ |
if (hp == NULL) |
269 |
+ |
error(SYSTEM, "out of memory in samp_hemi"); |
270 |
+ |
hp->rp = r; |
271 |
+ |
hp->ns = n; |
272 |
+ |
hp->acol[RED] = hp->acol[GRN] = hp->acol[BLU] = 0.0; |
273 |
+ |
memset(hp->sa, 0, sizeof(AMBSAMP)*n*n); |
274 |
+ |
hp->sampOK = 0; |
275 |
+ |
/* assign coefficient */ |
276 |
+ |
copycolor(hp->acoef, rcol); |
277 |
+ |
d = 1.0/(n*n); |
278 |
+ |
scalecolor(hp->acoef, d); |
279 |
+ |
/* make tangent plane axes */ |
280 |
+ |
if (!getperpendicular(hp->ux, r->ron, 1)) |
281 |
+ |
error(CONSISTENCY, "bad ray direction in samp_hemi"); |
282 |
+ |
VCROSS(hp->uy, r->ron, hp->ux); |
283 |
+ |
/* sample divisions */ |
284 |
+ |
for (i = hp->ns; i--; ) |
285 |
+ |
for (j = hp->ns; j--; ) |
286 |
+ |
hp->sampOK += ambsample(hp, i, j, 0); |
287 |
+ |
copycolor(rcol, hp->acol); |
288 |
+ |
if (!hp->sampOK) { /* utter failure? */ |
289 |
+ |
free(hp); |
290 |
+ |
return(NULL); |
291 |
+ |
} |
292 |
+ |
if (hp->sampOK < hp->ns*hp->ns) { |
293 |
+ |
hp->sampOK *= -1; /* soft failure */ |
294 |
+ |
return(hp); |
295 |
+ |
} |
296 |
+ |
if (hp->sampOK < 64) |
297 |
+ |
return(hp); /* insufficient for super-sampling */ |
298 |
+ |
n = ambssamp*wt + 0.5; |
299 |
+ |
if (n > 8) { /* perform super-sampling? */ |
300 |
+ |
ambsupersamp(hp, n); |
301 |
+ |
copycolor(rcol, hp->acol); |
302 |
+ |
} |
303 |
+ |
return(hp); /* all is well */ |
304 |
+ |
} |
305 |
+ |
|
306 |
+ |
|
307 |
|
/* Return brightness of farthest ambient sample */ |
308 |
|
static double |
309 |
|
back_ambval(AMBHEMI *hp, const int n1, const int n2, const int n3) |
644 |
|
ambcorral(AMBHEMI *hp, FVECT uv[2], const double r0, const double r1) |
645 |
|
{ |
646 |
|
const double max_d = 1.0/(minarad*ambacc + 0.001); |
647 |
< |
const double ang_res = 0.5*PI/(hp->ns-1); |
648 |
< |
const double ang_step = ang_res/((int)(16/PI*ang_res) + (1+FTINY)); |
647 |
> |
const double ang_res = 0.5*PI/hp->ns; |
648 |
> |
const double ang_step = ang_res/((int)(16/PI*ang_res) + 1.01); |
649 |
|
double avg_d = 0; |
650 |
|
uint32 flgs = 0; |
651 |
|
FVECT vec; |
653 |
|
double ang, a1; |
654 |
|
int i, j; |
655 |
|
/* don't bother for a few samples */ |
656 |
< |
if (hp->ns < 12) |
656 |
> |
if (hp->ns < 8) |
657 |
|
return(0); |
658 |
|
/* check distances overhead */ |
659 |
|
for (i = hp->ns*3/4; i-- > hp->ns>>2; ) |
671 |
|
if ((ap->d <= FTINY) | (ap->d >= max_d)) |
672 |
|
continue; /* too far or too near */ |
673 |
|
VSUB(vec, ap->p, hp->rp->rop); |
674 |
< |
u = DOT(vec, uv[0]) * ap->d; |
675 |
< |
v = DOT(vec, uv[1]) * ap->d; |
676 |
< |
if ((r0*r0*u*u + r1*r1*v*v) * ap->d*ap->d <= 1.0) |
674 |
> |
u = DOT(vec, uv[0]); |
675 |
> |
v = DOT(vec, uv[1]); |
676 |
> |
if ((r0*r0*u*u + r1*r1*v*v) * ap->d*ap->d <= u*u + v*v) |
677 |
|
continue; /* occluder outside ellipse */ |
678 |
|
ang = atan2a(v, u); /* else set direction flags */ |
679 |
< |
for (a1 = ang-.5*ang_res; a1 <= ang+.5*ang_res; a1 += ang_step) |
679 |
> |
for (a1 = ang-ang_res; a1 <= ang+ang_res; a1 += ang_step) |
680 |
|
flgs |= 1L<<(int)(16/PI*(a1 + 2.*PI*(a1 < 0))); |
681 |
|
} |
634 |
– |
/* add low-angle incident (< 20deg) */ |
635 |
– |
if (fabs(hp->rp->rod) <= 0.342) { |
636 |
– |
u = -DOT(hp->rp->rdir, uv[0]); |
637 |
– |
v = -DOT(hp->rp->rdir, uv[1]); |
638 |
– |
if ((r0*r0*u*u + r1*r1*v*v) > hp->rp->rot*hp->rp->rot) { |
639 |
– |
ang = atan2a(v, u); |
640 |
– |
ang += 2.*PI*(ang < 0); |
641 |
– |
ang *= 16/PI; |
642 |
– |
if ((ang < .5) | (ang >= 31.5)) |
643 |
– |
flgs |= 0x80000001; |
644 |
– |
else |
645 |
– |
flgs |= 3L<<(int)(ang-.5); |
646 |
– |
} |
647 |
– |
} |
682 |
|
return(flgs); |
683 |
|
} |
684 |
|
|
695 |
|
uint32 *crlp /* returned (optional) */ |
696 |
|
) |
697 |
|
{ |
698 |
< |
AMBHEMI *hp = inithemi(rcol, r, wt); |
665 |
< |
int cnt; |
698 |
> |
AMBHEMI *hp = samp_hemi(rcol, r, wt); |
699 |
|
FVECT my_uv[2]; |
700 |
< |
double d, K, acol[3]; |
700 |
> |
double d, K; |
701 |
|
AMBSAMP *ap; |
702 |
< |
int i, j; |
703 |
< |
/* check/initialize */ |
671 |
< |
if (hp == NULL) |
672 |
< |
return(0); |
702 |
> |
int i; |
703 |
> |
/* clear return values */ |
704 |
|
if (uv != NULL) |
705 |
|
memset(uv, 0, sizeof(FVECT)*2); |
706 |
|
if (ra != NULL) |
711 |
|
dg[0] = dg[1] = 0.0; |
712 |
|
if (crlp != NULL) |
713 |
|
*crlp = 0; |
714 |
< |
/* sample the hemisphere */ |
715 |
< |
acol[0] = acol[1] = acol[2] = 0.0; |
716 |
< |
cnt = 0; |
717 |
< |
for (i = hp->ns; i--; ) |
718 |
< |
for (j = hp->ns; j--; ) |
719 |
< |
if ((ap = ambsample(hp, i, j)) != NULL) { |
720 |
< |
addcolor(acol, ap->v); |
690 |
< |
++cnt; |
691 |
< |
} |
692 |
< |
if ((hp->ns < 4) | (cnt < hp->ns*hp->ns)) { |
693 |
< |
free(hp); /* inadequate sampling */ |
694 |
< |
copycolor(rcol, acol); |
695 |
< |
return(-cnt); /* value-only result */ |
714 |
> |
if (hp == NULL) /* sampling falure? */ |
715 |
> |
return(0); |
716 |
> |
|
717 |
> |
if ((ra == NULL) & (pg == NULL) & (dg == NULL) || |
718 |
> |
(hp->sampOK < 0) | (hp->ns < 6)) { |
719 |
> |
free(hp); /* Hessian not requested/possible */ |
720 |
> |
return(-1); /* value-only return value */ |
721 |
|
} |
722 |
< |
cnt = ambssamp*wt + 0.5; /* perform super-sampling? */ |
698 |
< |
if (cnt > 8) |
699 |
< |
ambsupersamp(acol, hp, cnt); |
700 |
< |
copycolor(rcol, acol); /* final indirect irradiance/PI */ |
701 |
< |
if ((ra == NULL) & (pg == NULL) & (dg == NULL)) { |
702 |
< |
free(hp); |
703 |
< |
return(-1); /* no Hessian or gradients requested */ |
704 |
< |
} |
705 |
< |
if ((d = bright(acol)) > FTINY) { /* normalize Y values */ |
722 |
> |
if ((d = bright(rcol)) > FTINY) { /* normalize Y values */ |
723 |
|
d = 0.99*(hp->ns*hp->ns)/d; |
724 |
|
K = 0.01; |
725 |
|
} else { /* or fall back on geometric Hessian */ |
754 |
|
if (ra[1] < minarad) |
755 |
|
ra[1] = minarad; |
756 |
|
} |
757 |
< |
ra[0] *= d = 1.0/sqrt(sqrt(wt)); |
757 |
> |
ra[0] *= d = 1.0/sqrt(wt); |
758 |
|
if ((ra[1] *= d) > 2.0*ra[0]) |
759 |
|
ra[1] = 2.0*ra[0]; |
760 |
|
if (ra[1] > maxarad) { |
763 |
|
ra[0] = maxarad; |
764 |
|
} |
765 |
|
/* flag encroached directions */ |
766 |
< |
if ((wt >= 0.89*AVGREFL) & (crlp != NULL)) |
766 |
> |
if (crlp != NULL) |
767 |
|
*crlp = ambcorral(hp, uv, ra[0]*ambacc, ra[1]*ambacc); |
768 |
|
if (pg != NULL) { /* cap gradient if necessary */ |
769 |
|
d = pg[0]*pg[0]*ra[0]*ra[0] + pg[1]*pg[1]*ra[1]*ra[1]; |
866 |
|
ndims--; |
867 |
|
multcolor(ar.rcol, ar.rcoef); /* apply coefficient */ |
868 |
|
addcolor(dp->v, ar.rcol); |
869 |
< |
/* use rt to improve gradient calc */ |
870 |
< |
if (ar.rt > FTINY && ar.rt < FHUGE) |
871 |
< |
dp->r += 1.0/ar.rt; |
869 |
> |
/* use rxt to improve gradient calc */ |
870 |
> |
if (ar.rxt > FTINY && ar.rxt < FHUGE) |
871 |
> |
dp->r += 1.0/ar.rxt; |
872 |
|
/* (re)initialize error */ |
873 |
|
if (dp->n++) { |
874 |
|
b2 = bright(dp->v)/dp->n - bright(ar.rcol); |