32 |
|
} sa[1]; /* sample array (extends struct) */ |
33 |
|
} AMBHEMI; /* ambient sample hemisphere */ |
34 |
|
|
35 |
< |
#define ambsamp(h,i,j) (h)->sa[(i)*(h)->ns + (j)] |
35 |
> |
typedef struct s_ambsamp AMBSAMP; |
36 |
|
|
37 |
+ |
#define ambsam(h,i,j) (h)->sa[(i)*(h)->ns + (j)] |
38 |
+ |
|
39 |
|
typedef struct { |
40 |
< |
FVECT r_i, r_i1, e_i, rI2_eJ2; |
41 |
< |
double nf, I1, I2; |
40 |
> |
FVECT r_i, r_i1, e_i, rcp, rI2_eJ2; |
41 |
> |
double I1, I2; |
42 |
|
} FFTRI; /* vectors and coefficients for Hessian calculation */ |
43 |
|
|
44 |
|
|
61 |
|
if (n < i) |
62 |
|
n = i; |
63 |
|
/* allocate sampling array */ |
64 |
< |
hp = (AMBHEMI *)malloc(sizeof(AMBHEMI) + |
63 |
< |
sizeof(struct s_ambsamp)*(n*n - 1)); |
64 |
> |
hp = (AMBHEMI *)malloc(sizeof(AMBHEMI) + sizeof(AMBSAMP)*(n*n - 1)); |
65 |
|
if (hp == NULL) |
66 |
|
return(NULL); |
67 |
|
hp->rp = r; |
71 |
|
d = 1.0/(n*n); |
72 |
|
scalecolor(hp->acoef, d); |
73 |
|
/* make tangent plane axes */ |
74 |
< |
hp->uy[0] = 0.1 - 0.2*frandom(); |
75 |
< |
hp->uy[1] = 0.1 - 0.2*frandom(); |
76 |
< |
hp->uy[2] = 0.1 - 0.2*frandom(); |
77 |
< |
for (i = 0; i < 3; i++) |
78 |
< |
if (r->ron[i] < 0.6 && r->ron[i] > -0.6) |
74 |
> |
hp->uy[0] = 0.5 - frandom(); |
75 |
> |
hp->uy[1] = 0.5 - frandom(); |
76 |
> |
hp->uy[2] = 0.5 - frandom(); |
77 |
> |
for (i = 3; i--; ) |
78 |
> |
if ((-0.6 < r->ron[i]) & (r->ron[i] < 0.6)) |
79 |
|
break; |
80 |
< |
if (i >= 3) |
81 |
< |
error(CONSISTENCY, "bad ray direction in inithemi()"); |
80 |
> |
if (i < 0) |
81 |
> |
error(CONSISTENCY, "bad ray direction in inithemi"); |
82 |
|
hp->uy[i] = 1.0; |
83 |
|
VCROSS(hp->ux, hp->uy, r->ron); |
84 |
|
normalize(hp->ux); |
88 |
|
} |
89 |
|
|
90 |
|
|
91 |
< |
static struct s_ambsamp * |
92 |
< |
ambsample( /* sample an ambient direction */ |
93 |
< |
AMBHEMI *hp, |
93 |
< |
int i, |
94 |
< |
int j |
95 |
< |
) |
91 |
> |
/* Sample ambient division and apply weighting coefficient */ |
92 |
> |
static int |
93 |
> |
getambsamp(RAY *arp, AMBHEMI *hp, int i, int j, int n) |
94 |
|
{ |
95 |
< |
struct s_ambsamp *ap = &ambsamp(hp,i,j); |
96 |
< |
RAY ar; |
99 |
< |
double spt[2], zd; |
100 |
< |
int ii; |
95 |
> |
int hlist[3], ii; |
96 |
> |
double spt[2], zd; |
97 |
|
/* ambient coefficient for weight */ |
98 |
|
if (ambacc > FTINY) |
99 |
< |
setcolor(ar.rcoef, AVGREFL, AVGREFL, AVGREFL); |
99 |
> |
setcolor(arp->rcoef, AVGREFL, AVGREFL, AVGREFL); |
100 |
|
else |
101 |
< |
copycolor(ar.rcoef, hp->acoef); |
102 |
< |
if (rayorigin(&ar, AMBIENT, hp->rp, ar.rcoef) < 0) |
103 |
< |
goto badsample; |
101 |
> |
copycolor(arp->rcoef, hp->acoef); |
102 |
> |
if (rayorigin(arp, AMBIENT, hp->rp, arp->rcoef) < 0) |
103 |
> |
return(0); |
104 |
|
if (ambacc > FTINY) { |
105 |
< |
multcolor(ar.rcoef, hp->acoef); |
106 |
< |
scalecolor(ar.rcoef, 1./AVGREFL); |
105 |
> |
multcolor(arp->rcoef, hp->acoef); |
106 |
> |
scalecolor(arp->rcoef, 1./AVGREFL); |
107 |
|
} |
108 |
< |
/* generate hemispherical sample */ |
109 |
< |
SDsquare2disk(spt, (i+.1+.8*frandom())/hp->ns, |
110 |
< |
(j+.1+.8*frandom())/hp->ns ); |
108 |
> |
hlist[0] = hp->rp->rno; |
109 |
> |
hlist[1] = i; |
110 |
> |
hlist[2] = j; |
111 |
> |
multisamp(spt, 2, urand(ilhash(hlist,3)+n)); |
112 |
> |
if (!n) { /* avoid border samples for n==0 */ |
113 |
> |
if ((spt[0] < 0.1) | (spt[0] > 0.9)) |
114 |
> |
spt[0] = 0.1 + 0.8*frandom(); |
115 |
> |
if ((spt[1] < 0.1) | (spt[1] > 0.9)) |
116 |
> |
spt[1] = 0.1 + 0.8*frandom(); |
117 |
> |
} |
118 |
> |
SDsquare2disk(spt, (i+spt[0])/hp->ns, (j+spt[1])/hp->ns); |
119 |
|
zd = sqrt(1. - spt[0]*spt[0] - spt[1]*spt[1]); |
120 |
|
for (ii = 3; ii--; ) |
121 |
< |
ar.rdir[ii] = spt[0]*hp->ux[ii] + |
121 |
> |
arp->rdir[ii] = spt[0]*hp->ux[ii] + |
122 |
|
spt[1]*hp->uy[ii] + |
123 |
|
zd*hp->rp->ron[ii]; |
124 |
< |
checknorm(ar.rdir); |
124 |
> |
checknorm(arp->rdir); |
125 |
|
dimlist[ndims++] = i*hp->ns + j + 90171; |
126 |
< |
rayvalue(&ar); /* evaluate ray */ |
127 |
< |
ndims--; |
128 |
< |
if (ar.rt > 20.0*maxarad) /* limit vertex distance */ |
129 |
< |
ar.rt = 20.0*maxarad; |
126 |
> |
rayvalue(arp); /* evaluate ray */ |
127 |
> |
ndims--; /* apply coefficient */ |
128 |
> |
multcolor(arp->rcol, arp->rcoef); |
129 |
> |
return(1); |
130 |
> |
} |
131 |
> |
|
132 |
> |
|
133 |
> |
static AMBSAMP * |
134 |
> |
ambsample( /* initial ambient division sample */ |
135 |
> |
AMBHEMI *hp, |
136 |
> |
int i, |
137 |
> |
int j |
138 |
> |
) |
139 |
> |
{ |
140 |
> |
AMBSAMP *ap = &ambsam(hp,i,j); |
141 |
> |
RAY ar; |
142 |
> |
/* generate hemispherical sample */ |
143 |
> |
if (!getambsamp(&ar, hp, i, j, 0)) |
144 |
> |
goto badsample; |
145 |
> |
/* limit vertex distance */ |
146 |
> |
if (ar.rt > 10.0*thescene.cusize) |
147 |
> |
ar.rt = 10.0*thescene.cusize; |
148 |
|
else if (ar.rt <= FTINY) /* should never happen! */ |
149 |
|
goto badsample; |
150 |
|
VSUM(ap->p, ar.rorg, ar.rdir, ar.rt); |
129 |
– |
multcolor(ar.rcol, ar.rcoef); /* apply coefficient */ |
151 |
|
copycolor(ap->v, ar.rcol); |
152 |
|
return(ap); |
153 |
|
badsample: |
157 |
|
} |
158 |
|
|
159 |
|
|
160 |
+ |
/* Estimate errors based on ambient division differences */ |
161 |
+ |
static float * |
162 |
+ |
getambdiffs(AMBHEMI *hp) |
163 |
+ |
{ |
164 |
+ |
float *earr = calloc(hp->ns*hp->ns, sizeof(float)); |
165 |
+ |
float *ep; |
166 |
+ |
AMBSAMP *ap; |
167 |
+ |
double b, d2; |
168 |
+ |
int i, j; |
169 |
+ |
|
170 |
+ |
if (earr == NULL) /* out of memory? */ |
171 |
+ |
return(NULL); |
172 |
+ |
/* compute squared neighbor diffs */ |
173 |
+ |
for (ap = hp->sa, ep = earr, i = 0; i < hp->ns; i++) |
174 |
+ |
for (j = 0; j < hp->ns; j++, ap++, ep++) { |
175 |
+ |
b = bright(ap[0].v); |
176 |
+ |
if (i) { /* from above */ |
177 |
+ |
d2 = b - bright(ap[-hp->ns].v); |
178 |
+ |
d2 *= d2; |
179 |
+ |
ep[0] += d2; |
180 |
+ |
ep[-hp->ns] += d2; |
181 |
+ |
} |
182 |
+ |
if (j) { /* from behind */ |
183 |
+ |
d2 = b - bright(ap[-1].v); |
184 |
+ |
d2 *= d2; |
185 |
+ |
ep[0] += d2; |
186 |
+ |
ep[-1] += d2; |
187 |
+ |
} |
188 |
+ |
} |
189 |
+ |
/* correct for number of neighbors */ |
190 |
+ |
earr[0] *= 2.f; |
191 |
+ |
earr[hp->ns-1] *= 2.f; |
192 |
+ |
earr[(hp->ns-1)*hp->ns] *= 2.f; |
193 |
+ |
earr[(hp->ns-1)*hp->ns + hp->ns-1] *= 2.f; |
194 |
+ |
for (i = 1; i < hp->ns-1; i++) { |
195 |
+ |
earr[i*hp->ns] *= 4./3.; |
196 |
+ |
earr[i*hp->ns + hp->ns-1] *= 4./3.; |
197 |
+ |
} |
198 |
+ |
for (j = 1; j < hp->ns-1; j++) { |
199 |
+ |
earr[j] *= 4./3.; |
200 |
+ |
earr[(hp->ns-1)*hp->ns + j] *= 4./3.; |
201 |
+ |
} |
202 |
+ |
return(earr); |
203 |
+ |
} |
204 |
+ |
|
205 |
+ |
|
206 |
+ |
/* Perform super-sampling on hemisphere (introduces bias) */ |
207 |
+ |
static void |
208 |
+ |
ambsupersamp(double acol[3], AMBHEMI *hp, int cnt) |
209 |
+ |
{ |
210 |
+ |
float *earr = getambdiffs(hp); |
211 |
+ |
double e2sum = 0; |
212 |
+ |
AMBSAMP *ap; |
213 |
+ |
RAY ar; |
214 |
+ |
COLOR asum; |
215 |
+ |
float *ep; |
216 |
+ |
int i, j, n; |
217 |
+ |
|
218 |
+ |
if (earr == NULL) /* just skip calc. if no memory */ |
219 |
+ |
return; |
220 |
+ |
/* add up estimated variances */ |
221 |
+ |
for (ep = earr + hp->ns*hp->ns; ep-- > earr; ) |
222 |
+ |
e2sum += *ep; |
223 |
+ |
ep = earr; /* perform super-sampling */ |
224 |
+ |
for (ap = hp->sa, i = 0; i < hp->ns; i++) |
225 |
+ |
for (j = 0; j < hp->ns; j++, ap++) { |
226 |
+ |
int nss = *ep/e2sum*cnt + frandom(); |
227 |
+ |
setcolor(asum, 0., 0., 0.); |
228 |
+ |
for (n = 1; n <= nss; n++) { |
229 |
+ |
if (!getambsamp(&ar, hp, i, j, n)) { |
230 |
+ |
nss = n-1; |
231 |
+ |
break; |
232 |
+ |
} |
233 |
+ |
addcolor(asum, ar.rcol); |
234 |
+ |
} |
235 |
+ |
if (nss) { /* update returned ambient value */ |
236 |
+ |
const double ssf = 1./(nss + 1); |
237 |
+ |
for (n = 3; n--; ) |
238 |
+ |
acol[n] += ssf*colval(asum,n) + |
239 |
+ |
(ssf - 1.)*colval(ap->v,n); |
240 |
+ |
} |
241 |
+ |
e2sum -= *ep++; /* update remainders */ |
242 |
+ |
cnt -= nss; |
243 |
+ |
} |
244 |
+ |
free(earr); |
245 |
+ |
} |
246 |
+ |
|
247 |
+ |
|
248 |
|
/* Compute vectors and coefficients for Hessian/gradient calcs */ |
249 |
|
static void |
250 |
|
comp_fftri(FFTRI *ftp, FVECT ap0, FVECT ap1, FVECT rop) |
251 |
|
{ |
252 |
< |
FVECT vcp; |
144 |
< |
double dot_e, dot_er, rdot_r, rdot_r1, J2; |
252 |
> |
double rdot_cp, dot_e, dot_er, rdot_r, rdot_r1, J2; |
253 |
|
int i; |
254 |
|
|
255 |
|
VSUB(ftp->r_i, ap0, rop); |
256 |
|
VSUB(ftp->r_i1, ap1, rop); |
257 |
|
VSUB(ftp->e_i, ap1, ap0); |
258 |
< |
VCROSS(vcp, ftp->e_i, ftp->r_i); |
259 |
< |
ftp->nf = 1.0/DOT(vcp,vcp); |
258 |
> |
VCROSS(ftp->rcp, ftp->r_i, ftp->r_i1); |
259 |
> |
rdot_cp = 1.0/DOT(ftp->rcp,ftp->rcp); |
260 |
|
dot_e = DOT(ftp->e_i,ftp->e_i); |
261 |
|
dot_er = DOT(ftp->e_i, ftp->r_i); |
262 |
|
rdot_r = 1.0/DOT(ftp->r_i,ftp->r_i); |
263 |
|
rdot_r1 = 1.0/DOT(ftp->r_i1,ftp->r_i1); |
264 |
|
ftp->I1 = acos( DOT(ftp->r_i, ftp->r_i1) * sqrt(rdot_r*rdot_r1) ) * |
265 |
< |
sqrt( ftp->nf ); |
265 |
> |
sqrt( rdot_cp ); |
266 |
|
ftp->I2 = ( DOT(ftp->e_i, ftp->r_i1)*rdot_r1 - dot_er*rdot_r + |
267 |
< |
dot_e*ftp->I1 )*0.5*ftp->nf; |
267 |
> |
dot_e*ftp->I1 )*0.5*rdot_cp; |
268 |
|
J2 = ( 0.5*(rdot_r - rdot_r1) - dot_er*ftp->I2 ) / dot_e; |
269 |
|
for (i = 3; i--; ) |
270 |
|
ftp->rI2_eJ2[i] = ftp->I2*ftp->r_i[i] + J2*ftp->e_i[i]; |
288 |
|
static void |
289 |
|
comp_hessian(FVECT hess[3], FFTRI *ftp, FVECT nrm) |
290 |
|
{ |
291 |
< |
FVECT vcp; |
291 |
> |
FVECT ncp; |
292 |
|
FVECT m1[3], m2[3], m3[3], m4[3]; |
293 |
|
double d1, d2, d3, d4; |
294 |
|
double I3, J3, K3; |
298 |
|
d2 = 1.0/DOT(ftp->r_i1,ftp->r_i1); |
299 |
|
d3 = 1.0/DOT(ftp->e_i,ftp->e_i); |
300 |
|
d4 = DOT(ftp->e_i, ftp->r_i); |
301 |
< |
I3 = 0.25*ftp->nf*( DOT(ftp->e_i, ftp->r_i1)*d2*d2 - d4*d1*d1 + |
302 |
< |
3.0/d3*ftp->I2 ); |
301 |
> |
I3 = ( DOT(ftp->e_i, ftp->r_i1)*d2*d2 - d4*d1*d1 + 3.0/d3*ftp->I2 ) |
302 |
> |
/ ( 4.0*DOT(ftp->rcp,ftp->rcp) ); |
303 |
|
J3 = 0.25*d3*(d1*d1 - d2*d2) - d4*d3*I3; |
304 |
|
K3 = d3*(ftp->I2 - I3/d1 - 2.0*d4*J3); |
305 |
|
/* intermediate matrices */ |
306 |
< |
VCROSS(vcp, nrm, ftp->e_i); |
307 |
< |
compose_matrix(m1, vcp, ftp->rI2_eJ2); |
306 |
> |
VCROSS(ncp, nrm, ftp->e_i); |
307 |
> |
compose_matrix(m1, ncp, ftp->rI2_eJ2); |
308 |
|
compose_matrix(m2, ftp->r_i, ftp->r_i); |
309 |
|
compose_matrix(m3, ftp->e_i, ftp->e_i); |
310 |
|
compose_matrix(m4, ftp->r_i, ftp->e_i); |
311 |
< |
VCROSS(vcp, ftp->r_i, ftp->e_i); |
204 |
< |
d1 = DOT(nrm, vcp); |
311 |
> |
d1 = DOT(nrm, ftp->rcp); |
312 |
|
d2 = -d1*ftp->I2; |
313 |
|
d1 *= 2.0; |
314 |
|
for (i = 3; i--; ) /* final matrix sum */ |
352 |
|
static void |
353 |
|
comp_gradient(FVECT grad, FFTRI *ftp, FVECT nrm) |
354 |
|
{ |
355 |
< |
FVECT vcp; |
355 |
> |
FVECT ncp; |
356 |
|
double f1; |
357 |
|
int i; |
358 |
|
|
359 |
< |
VCROSS(vcp, ftp->r_i, ftp->r_i1); |
360 |
< |
f1 = 2.0*DOT(nrm, vcp); |
254 |
< |
VCROSS(vcp, nrm, ftp->e_i); |
359 |
> |
f1 = 2.0*DOT(nrm, ftp->rcp); |
360 |
> |
VCROSS(ncp, nrm, ftp->e_i); |
361 |
|
for (i = 3; i--; ) |
362 |
< |
grad[i] = (-0.5/PI)*( ftp->I1*vcp[i] + f1*ftp->rI2_eJ2[i] ); |
362 |
> |
grad[i] = (-0.5/PI)*( ftp->I1*ncp[i] + f1*ftp->rI2_eJ2[i] ); |
363 |
|
} |
364 |
|
|
365 |
|
|
386 |
|
|
387 |
|
/* Return brightness of furthest ambient sample */ |
388 |
|
static COLORV |
389 |
< |
back_ambval(struct s_ambsamp *ap1, struct s_ambsamp *ap2, |
284 |
< |
struct s_ambsamp *ap3, FVECT orig) |
389 |
> |
back_ambval(AMBSAMP *ap1, AMBSAMP *ap2, AMBSAMP *ap3, FVECT orig) |
390 |
|
{ |
391 |
|
COLORV vback; |
392 |
|
FVECT vec; |
426 |
|
hess2[0][1] = DOT(uv[0], b); |
427 |
|
hess2[1][0] = DOT(uv[1], a); |
428 |
|
hess2[1][1] = DOT(uv[1], b); |
429 |
< |
/* compute eigenvalues */ |
430 |
< |
if ( quadratic(evalue, 1.0, -hess2[0][0]-hess2[1][1], |
431 |
< |
hess2[0][0]*hess2[1][1]-hess2[0][1]*hess2[1][0]) != 2 || |
432 |
< |
(evalue[0] = fabs(evalue[0])) <= FTINY*FTINY || |
433 |
< |
(evalue[1] = fabs(evalue[1])) <= FTINY*FTINY ) |
429 |
> |
/* compute eigenvalue(s) */ |
430 |
> |
i = quadratic(evalue, 1.0, -hess2[0][0]-hess2[1][1], |
431 |
> |
hess2[0][0]*hess2[1][1]-hess2[0][1]*hess2[1][0]); |
432 |
> |
if (i == 1) /* double-root (circle) */ |
433 |
> |
evalue[1] = evalue[0]; |
434 |
> |
if (!i || ((evalue[0] = fabs(evalue[0])) <= FTINY*FTINY) | |
435 |
> |
((evalue[1] = fabs(evalue[1])) <= FTINY*FTINY) ) |
436 |
|
error(INTERNAL, "bad eigenvalue calculation"); |
437 |
|
|
438 |
|
if (evalue[0] > evalue[1]) { |
492 |
|
} |
493 |
|
/* compute first row of edges */ |
494 |
|
for (j = 0; j < hp->ns-1; j++) { |
495 |
< |
comp_fftri(&fftr, ambsamp(hp,0,j).p, |
496 |
< |
ambsamp(hp,0,j+1).p, hp->rp->rop); |
495 |
> |
comp_fftri(&fftr, ambsam(hp,0,j).p, |
496 |
> |
ambsam(hp,0,j+1).p, hp->rp->rop); |
497 |
|
if (hessrow != NULL) |
498 |
|
comp_hessian(hessrow[j], &fftr, hp->rp->ron); |
499 |
|
if (gradrow != NULL) |
503 |
|
for (i = 0; i < hp->ns-1; i++) { |
504 |
|
FVECT hesscol[3]; /* compute first vertical edge */ |
505 |
|
FVECT gradcol; |
506 |
< |
comp_fftri(&fftr, ambsamp(hp,i,0).p, |
507 |
< |
ambsamp(hp,i+1,0).p, hp->rp->rop); |
506 |
> |
comp_fftri(&fftr, ambsam(hp,i,0).p, |
507 |
> |
ambsam(hp,i+1,0).p, hp->rp->rop); |
508 |
|
if (hessrow != NULL) |
509 |
|
comp_hessian(hesscol, &fftr, hp->rp->ron); |
510 |
|
if (gradrow != NULL) |
513 |
|
FVECT hessdia[3]; /* compute triangle contributions */ |
514 |
|
FVECT graddia; |
515 |
|
COLORV backg; |
516 |
< |
backg = back_ambval(&ambsamp(hp,i,j), &ambsamp(hp,i,j+1), |
517 |
< |
&ambsamp(hp,i+1,j), hp->rp->rop); |
516 |
> |
backg = back_ambval(&ambsam(hp,i,j), &ambsam(hp,i,j+1), |
517 |
> |
&ambsam(hp,i+1,j), hp->rp->rop); |
518 |
|
/* diagonal (inner) edge */ |
519 |
< |
comp_fftri(&fftr, ambsamp(hp,i,j+1).p, |
520 |
< |
ambsamp(hp,i+1,j).p, hp->rp->rop); |
519 |
> |
comp_fftri(&fftr, ambsam(hp,i,j+1).p, |
520 |
> |
ambsam(hp,i+1,j).p, hp->rp->rop); |
521 |
|
if (hessrow != NULL) { |
522 |
|
comp_hessian(hessdia, &fftr, hp->rp->ron); |
523 |
|
rev_hessian(hesscol); |
524 |
|
add2hessian(hessian, hessrow[j], hessdia, hesscol, backg); |
525 |
|
} |
526 |
< |
if (gradient != NULL) { |
526 |
> |
if (gradrow != NULL) { |
527 |
|
comp_gradient(graddia, &fftr, hp->rp->ron); |
528 |
|
rev_gradient(gradcol); |
529 |
|
add2gradient(gradient, gradrow[j], graddia, gradcol, backg); |
530 |
|
} |
531 |
|
/* initialize edge in next row */ |
532 |
< |
comp_fftri(&fftr, ambsamp(hp,i+1,j+1).p, |
533 |
< |
ambsamp(hp,i+1,j).p, hp->rp->rop); |
532 |
> |
comp_fftri(&fftr, ambsam(hp,i+1,j+1).p, |
533 |
> |
ambsam(hp,i+1,j).p, hp->rp->rop); |
534 |
|
if (hessrow != NULL) |
535 |
|
comp_hessian(hessrow[j], &fftr, hp->rp->ron); |
536 |
|
if (gradrow != NULL) |
537 |
|
comp_gradient(gradrow[j], &fftr, hp->rp->ron); |
538 |
|
/* new column edge & paired triangle */ |
539 |
< |
backg = back_ambval(&ambsamp(hp,i,j+1), &ambsamp(hp,i+1,j+1), |
540 |
< |
&ambsamp(hp,i+1,j), hp->rp->rop); |
541 |
< |
comp_fftri(&fftr, ambsamp(hp,i,j+1).p, ambsamp(hp,i+1,j+1).p, |
539 |
> |
backg = back_ambval(&ambsam(hp,i,j+1), &ambsam(hp,i+1,j+1), |
540 |
> |
&ambsam(hp,i+1,j), hp->rp->rop); |
541 |
> |
comp_fftri(&fftr, ambsam(hp,i,j+1).p, ambsam(hp,i+1,j+1).p, |
542 |
|
hp->rp->rop); |
543 |
|
if (hessrow != NULL) { |
544 |
|
comp_hessian(hesscol, &fftr, hp->rp->ron); |
573 |
|
static void |
574 |
|
ambdirgrad(AMBHEMI *hp, FVECT uv[2], float dg[2]) |
575 |
|
{ |
576 |
< |
struct s_ambsamp *ap; |
577 |
< |
double dgsum[2]; |
578 |
< |
int n; |
579 |
< |
FVECT vd; |
580 |
< |
double gfact; |
576 |
> |
AMBSAMP *ap; |
577 |
> |
double dgsum[2]; |
578 |
> |
int n; |
579 |
> |
FVECT vd; |
580 |
> |
double gfact; |
581 |
|
|
582 |
|
dgsum[0] = dgsum[1] = 0.0; /* sum values times -tan(theta) */ |
583 |
|
for (ap = hp->sa, n = hp->ns*hp->ns; n--; ap++) { |
585 |
|
VSUB(vd, ap->p, hp->rp->rop); |
586 |
|
/* brightness over cosine factor */ |
587 |
|
gfact = colval(ap->v,CIEY) / DOT(hp->rp->ron, vd); |
588 |
< |
/* -sine = -proj_radius/vd_length */ |
589 |
< |
dgsum[0] += DOT(uv[1], vd) * gfact; |
590 |
< |
dgsum[1] -= DOT(uv[0], vd) * gfact; |
588 |
> |
/* sine = proj_radius/vd_length */ |
589 |
> |
dgsum[0] -= DOT(uv[1], vd) * gfact; |
590 |
> |
dgsum[1] += DOT(uv[0], vd) * gfact; |
591 |
|
} |
592 |
|
dg[0] = dgsum[0] / (hp->ns*hp->ns); |
593 |
|
dg[1] = dgsum[1] / (hp->ns*hp->ns); |
605 |
|
float dg[2] /* returned (optional) */ |
606 |
|
) |
607 |
|
{ |
608 |
< |
AMBHEMI *hp = inithemi(rcol, r, wt); |
609 |
< |
int cnt = 0; |
610 |
< |
FVECT my_uv[2]; |
611 |
< |
double d, acol[3]; |
612 |
< |
struct s_ambsamp *ap; |
613 |
< |
int i, j; |
608 |
> |
AMBHEMI *hp = inithemi(rcol, r, wt); |
609 |
> |
int cnt = 0; |
610 |
> |
FVECT my_uv[2]; |
611 |
> |
double d, K, acol[3]; |
612 |
> |
AMBSAMP *ap; |
613 |
> |
int i, j; |
614 |
|
/* check/initialize */ |
615 |
|
if (hp == NULL) |
616 |
|
return(0); |
635 |
|
free(hp); |
636 |
|
return(0); /* no valid samples */ |
637 |
|
} |
638 |
+ |
if (cnt < hp->ns*hp->ns) { /* incomplete sampling? */ |
639 |
+ |
copycolor(rcol, acol); |
640 |
+ |
free(hp); |
641 |
+ |
return(-1); /* return value w/o Hessian */ |
642 |
+ |
} |
643 |
+ |
cnt = ambssamp*wt + 0.5; /* perform super-sampling? */ |
644 |
+ |
if (cnt > 0) |
645 |
+ |
ambsupersamp(acol, hp, cnt); |
646 |
|
copycolor(rcol, acol); /* final indirect irradiance/PI */ |
647 |
< |
if (cnt < hp->ns*hp->ns || /* incomplete sampling? */ |
533 |
< |
(ra == NULL) & (pg == NULL) & (dg == NULL)) { |
647 |
> |
if ((ra == NULL) & (pg == NULL) & (dg == NULL)) { |
648 |
|
free(hp); |
649 |
|
return(-1); /* no radius or gradient calc. */ |
650 |
|
} |
651 |
< |
if (bright(acol) > FTINY) /* normalize Y values */ |
652 |
< |
d = cnt/bright(acol); |
653 |
< |
else |
651 |
> |
if (bright(acol) > FTINY) { /* normalize Y values */ |
652 |
> |
d = 0.99*cnt/bright(acol); |
653 |
> |
K = 0.01; |
654 |
> |
} else { /* geometric Hessian fall-back */ |
655 |
|
d = 0.0; |
656 |
+ |
K = 1.0; |
657 |
+ |
pg = NULL; |
658 |
+ |
dg = NULL; |
659 |
+ |
} |
660 |
|
ap = hp->sa; /* relative Y channel from here on... */ |
661 |
|
for (i = hp->ns*hp->ns; i--; ap++) |
662 |
< |
colval(ap->v,CIEY) = bright(ap->v)*d + 0.01; |
662 |
> |
colval(ap->v,CIEY) = bright(ap->v)*d + K; |
663 |
|
|
664 |
|
if (uv == NULL) /* make sure we have axis pointers */ |
665 |
|
uv = my_uv; |
670 |
|
ambdirgrad(hp, uv, dg); |
671 |
|
|
672 |
|
if (ra != NULL) { /* scale/clamp radii */ |
673 |
+ |
if (pg != NULL) { |
674 |
+ |
if (ra[0]*(d = fabs(pg[0])) > 1.0) |
675 |
+ |
ra[0] = 1.0/d; |
676 |
+ |
if (ra[1]*(d = fabs(pg[1])) > 1.0) |
677 |
+ |
ra[1] = 1.0/d; |
678 |
+ |
if (ra[0] > ra[1]) |
679 |
+ |
ra[0] = ra[1]; |
680 |
+ |
} |
681 |
|
if (ra[0] < minarad) { |
682 |
|
ra[0] = minarad; |
683 |
|
if (ra[1] < minarad) |
690 |
|
ra[1] = maxarad; |
691 |
|
if (ra[0] > maxarad) |
692 |
|
ra[0] = maxarad; |
693 |
+ |
} |
694 |
+ |
if (pg != NULL) { /* cap gradient if necessary */ |
695 |
+ |
d = pg[0]*pg[0]*ra[0]*ra[0] + pg[1]*pg[1]*ra[1]*ra[1]; |
696 |
+ |
if (d > 1.0) { |
697 |
+ |
d = 1.0/sqrt(d); |
698 |
+ |
pg[0] *= d; |
699 |
+ |
pg[1] *= d; |
700 |
+ |
} |
701 |
|
} |
702 |
|
} |
703 |
|
free(hp); /* clean up and return */ |