25 |
|
|
26 |
|
extern void SDsquare2disk(double ds[2], double seedx, double seedy); |
27 |
|
|
28 |
– |
/* vertex direction bit positions */ |
29 |
– |
#define VDB_xy 0 |
30 |
– |
#define VDB_y 01 |
31 |
– |
#define VDB_x 02 |
32 |
– |
#define VDB_Xy 03 |
33 |
– |
#define VDB_xY 04 |
34 |
– |
#define VDB_X 05 |
35 |
– |
#define VDB_Y 06 |
36 |
– |
#define VDB_XY 07 |
37 |
– |
/* get opposite vertex direction bit */ |
38 |
– |
#define VDB_OPP(f) (~(f) & 07) |
39 |
– |
/* adjacent triangle vertex flags */ |
40 |
– |
static const int adjacent_trifl[8] = { |
41 |
– |
0, /* forbidden diagonal */ |
42 |
– |
1<<VDB_x|1<<VDB_y|1<<VDB_Xy, |
43 |
– |
1<<VDB_y|1<<VDB_x|1<<VDB_xY, |
44 |
– |
1<<VDB_y|1<<VDB_Xy|1<<VDB_X, |
45 |
– |
1<<VDB_x|1<<VDB_xY|1<<VDB_Y, |
46 |
– |
1<<VDB_Xy|1<<VDB_X|1<<VDB_Y, |
47 |
– |
1<<VDB_xY|1<<VDB_Y|1<<VDB_X, |
48 |
– |
0, /* forbidden diagonal */ |
49 |
– |
}; |
50 |
– |
|
28 |
|
typedef struct { |
29 |
|
COLOR v; /* hemisphere sample value */ |
30 |
|
float d; /* reciprocal distance (1/rt) */ |
39 |
|
AMBSAMP sa[1]; /* sample array (extends struct) */ |
40 |
|
} AMBHEMI; /* ambient sample hemisphere */ |
41 |
|
|
42 |
< |
#define ambndx(h,i,j) ((i)*(h)->ns + (j)) |
43 |
< |
#define ambsam(h,i,j) (h)->sa[ambndx(h,i,j)] |
42 |
> |
#define AI(h,i,j) ((i)*(h)->ns + (j)) |
43 |
> |
#define ambsam(h,i,j) (h)->sa[AI(h,i,j)] |
44 |
|
|
45 |
|
typedef struct { |
46 |
|
FVECT r_i, r_i1, e_i, rcp, rI2_eJ2; |
47 |
|
double I1, I2; |
71 |
– |
int valid; |
48 |
|
} FFTRI; /* vectors and coefficients for Hessian calculation */ |
49 |
|
|
50 |
|
|
75 |
– |
/* Get index for adjacent vertex */ |
76 |
– |
static int |
77 |
– |
adjacent_verti(AMBHEMI *hp, int i, int j, int dbit) |
78 |
– |
{ |
79 |
– |
int i0 = i*hp->ns + j; |
80 |
– |
|
81 |
– |
switch (dbit) { |
82 |
– |
case VDB_y: return(i0 - hp->ns); |
83 |
– |
case VDB_x: return(i0 - 1); |
84 |
– |
case VDB_Xy: return(i0 - hp->ns + 1); |
85 |
– |
case VDB_xY: return(i0 + hp->ns - 1); |
86 |
– |
case VDB_X: return(i0 + 1); |
87 |
– |
case VDB_Y: return(i0 + hp->ns); |
88 |
– |
/* the following should never occur */ |
89 |
– |
case VDB_xy: return(i0 - hp->ns - 1); |
90 |
– |
case VDB_XY: return(i0 + hp->ns + 1); |
91 |
– |
} |
92 |
– |
return(-1); |
93 |
– |
} |
94 |
– |
|
95 |
– |
|
96 |
– |
/* Get vertex direction bit for the opposite edge to complete triangle */ |
97 |
– |
static int |
98 |
– |
vdb_edge(int db1, int db2) |
99 |
– |
{ |
100 |
– |
switch (db1) { |
101 |
– |
case VDB_x: return(db2==VDB_y ? VDB_Xy : VDB_Y); |
102 |
– |
case VDB_y: return(db2==VDB_x ? VDB_xY : VDB_X); |
103 |
– |
case VDB_X: return(db2==VDB_Xy ? VDB_y : VDB_xY); |
104 |
– |
case VDB_Y: return(db2==VDB_xY ? VDB_x : VDB_Xy); |
105 |
– |
case VDB_xY: return(db2==VDB_x ? VDB_y : VDB_X); |
106 |
– |
case VDB_Xy: return(db2==VDB_y ? VDB_x : VDB_Y); |
107 |
– |
} |
108 |
– |
error(INTERNAL, "forbidden diagonal in vdb_edge()"); |
109 |
– |
return(-1); |
110 |
– |
} |
111 |
– |
|
112 |
– |
|
51 |
|
static AMBHEMI * |
52 |
|
inithemi( /* initialize sampling hemisphere */ |
53 |
|
COLOR ac, |
128 |
|
spt[1]*hp->uy[ii] + |
129 |
|
zd*hp->rp->ron[ii]; |
130 |
|
checknorm(arp->rdir); |
131 |
< |
dimlist[ndims++] = ambndx(hp,i,j) + 90171; |
131 |
> |
dimlist[ndims++] = AI(hp,i,j) + 90171; |
132 |
|
rayvalue(arp); /* evaluate ray */ |
133 |
|
ndims--; /* apply coefficient */ |
134 |
|
multcolor(arp->rcol, arp->rcoef); |
181 |
|
ep[0] += d2; |
182 |
|
ep[-hp->ns] += d2; |
183 |
|
} |
184 |
< |
if (j) { /* from behind */ |
185 |
< |
d2 = b - bright(ap[-1].v); |
186 |
< |
d2 *= d2; |
187 |
< |
ep[0] += d2; |
188 |
< |
ep[-1] += d2; |
189 |
< |
} |
184 |
> |
if (!j) continue; |
185 |
> |
/* from behind */ |
186 |
> |
d2 = b - bright(ap[-1].v); |
187 |
> |
d2 *= d2; |
188 |
> |
ep[0] += d2; |
189 |
> |
ep[-1] += d2; |
190 |
> |
if (!i) continue; |
191 |
> |
/* diagonal */ |
192 |
> |
d2 = b - bright(ap[-hp->ns-1].v); |
193 |
> |
d2 *= d2; |
194 |
> |
ep[0] += d2; |
195 |
> |
ep[-hp->ns-1] += d2; |
196 |
|
} |
197 |
|
/* correct for number of neighbors */ |
198 |
< |
earr[0] *= 2.f; |
199 |
< |
earr[hp->ns-1] *= 2.f; |
200 |
< |
earr[(hp->ns-1)*hp->ns] *= 2.f; |
201 |
< |
earr[(hp->ns-1)*hp->ns + hp->ns-1] *= 2.f; |
198 |
> |
earr[0] *= 8./3.; |
199 |
> |
earr[hp->ns-1] *= 8./3.; |
200 |
> |
earr[(hp->ns-1)*hp->ns] *= 8./3.; |
201 |
> |
earr[(hp->ns-1)*hp->ns + hp->ns-1] *= 8./3.; |
202 |
|
for (i = 1; i < hp->ns-1; i++) { |
203 |
< |
earr[i*hp->ns] *= 4./3.; |
204 |
< |
earr[i*hp->ns + hp->ns-1] *= 4./3.; |
203 |
> |
earr[i*hp->ns] *= 8./5.; |
204 |
> |
earr[i*hp->ns + hp->ns-1] *= 8./5.; |
205 |
|
} |
206 |
|
for (j = 1; j < hp->ns-1; j++) { |
207 |
< |
earr[j] *= 4./3.; |
208 |
< |
earr[(hp->ns-1)*hp->ns + j] *= 4./3.; |
207 |
> |
earr[j] *= 8./5.; |
208 |
> |
earr[(hp->ns-1)*hp->ns + j] *= 8./5.; |
209 |
|
} |
210 |
|
return(earr); |
211 |
|
} |
216 |
|
ambsupersamp(double acol[3], AMBHEMI *hp, int cnt) |
217 |
|
{ |
218 |
|
float *earr = getambdiffs(hp); |
219 |
< |
double e2sum = 0; |
219 |
> |
double e2rem = 0; |
220 |
|
AMBSAMP *ap; |
221 |
|
RAY ar; |
222 |
< |
COLOR asum; |
222 |
> |
double asum[3]; |
223 |
|
float *ep; |
224 |
< |
int i, j, n; |
224 |
> |
int i, j, n, nss; |
225 |
|
|
226 |
|
if (earr == NULL) /* just skip calc. if no memory */ |
227 |
|
return; |
228 |
< |
/* add up estimated variances */ |
229 |
< |
for (ep = earr + hp->ns*hp->ns; ep-- > earr; ) |
230 |
< |
e2sum += *ep; |
228 |
> |
/* accumulate estimated variances */ |
229 |
> |
for (ep = earr + hp->ns*hp->ns; ep > earr; ) |
230 |
> |
e2rem += *--ep; |
231 |
|
ep = earr; /* perform super-sampling */ |
232 |
|
for (ap = hp->sa, i = 0; i < hp->ns; i++) |
233 |
|
for (j = 0; j < hp->ns; j++, ap++) { |
234 |
< |
int nss = *ep/e2sum*cnt + frandom(); |
235 |
< |
setcolor(asum, 0., 0., 0.); |
234 |
> |
if (e2rem <= FTINY) |
235 |
> |
goto done; /* nothing left to do */ |
236 |
> |
nss = *ep/e2rem*cnt + frandom(); |
237 |
> |
asum[0] = asum[1] = asum[2] = 0.0; |
238 |
|
for (n = 1; n <= nss; n++) { |
239 |
|
if (!getambsamp(&ar, hp, i, j, n)) { |
240 |
|
nss = n-1; |
243 |
|
addcolor(asum, ar.rcol); |
244 |
|
} |
245 |
|
if (nss) { /* update returned ambient value */ |
246 |
< |
const double ssf = 1./(nss + 1); |
246 |
> |
const double ssf = 1./(nss + 1.); |
247 |
|
for (n = 3; n--; ) |
248 |
< |
acol[n] += ssf*colval(asum,n) + |
248 |
> |
acol[n] += ssf*asum[n] + |
249 |
|
(ssf - 1.)*colval(ap->v,n); |
250 |
|
} |
251 |
< |
e2sum -= *ep++; /* update remainders */ |
251 |
> |
e2rem -= *ep++; /* update remainders */ |
252 |
|
cnt -= nss; |
253 |
|
} |
254 |
+ |
done: |
255 |
|
free(earr); |
256 |
|
} |
257 |
|
|
258 |
|
|
312 |
– |
/* Compute vertex flags, indicating farthest in each direction */ |
313 |
– |
static uby8 * |
314 |
– |
vertex_flags(AMBHEMI *hp) |
315 |
– |
{ |
316 |
– |
uby8 *vflags = (uby8 *)calloc(hp->ns*hp->ns, sizeof(uby8)); |
317 |
– |
uby8 *vf; |
318 |
– |
AMBSAMP *ap; |
319 |
– |
int i, j; |
320 |
– |
|
321 |
– |
if (vflags == NULL) |
322 |
– |
error(SYSTEM, "out of memory in vertex_flags()"); |
323 |
– |
vf = vflags; |
324 |
– |
ap = hp->sa; /* compute farthest along first row */ |
325 |
– |
for (j = 0; j < hp->ns-1; j++, vf++, ap++) |
326 |
– |
if (ap[0].d <= ap[1].d) |
327 |
– |
vf[0] |= 1<<VDB_X; |
328 |
– |
else |
329 |
– |
vf[1] |= 1<<VDB_x; |
330 |
– |
++vf; ++ap; |
331 |
– |
/* flag subsequent rows */ |
332 |
– |
for (i = 1; i < hp->ns; i++) { |
333 |
– |
for (j = 0; j < hp->ns-1; j++, vf++, ap++) { |
334 |
– |
if (ap[0].d <= ap[-hp->ns].d) /* row before */ |
335 |
– |
vf[0] |= 1<<VDB_y; |
336 |
– |
else |
337 |
– |
vf[-hp->ns] |= 1<<VDB_Y; |
338 |
– |
if (ap[0].d <= ap[1-hp->ns].d) /* diagonal we care about */ |
339 |
– |
vf[0] |= 1<<VDB_Xy; |
340 |
– |
else |
341 |
– |
vf[1-hp->ns] |= 1<<VDB_xY; |
342 |
– |
if (ap[0].d <= ap[1].d) /* column after */ |
343 |
– |
vf[0] |= 1<<VDB_X; |
344 |
– |
else |
345 |
– |
vf[1] |= 1<<VDB_x; |
346 |
– |
} |
347 |
– |
if (ap[0].d <= ap[-hp->ns].d) /* final column edge */ |
348 |
– |
vf[0] |= 1<<VDB_y; |
349 |
– |
else |
350 |
– |
vf[-hp->ns] |= 1<<VDB_Y; |
351 |
– |
++vf; ++ap; |
352 |
– |
} |
353 |
– |
return(vflags); |
354 |
– |
} |
355 |
– |
|
356 |
– |
|
259 |
|
/* Return brightness of farthest ambient sample */ |
260 |
|
static double |
261 |
< |
back_ambval(AMBHEMI *hp, int i, int j, int dbit1, int dbit2, const uby8 *vflags) |
261 |
> |
back_ambval(AMBHEMI *hp, const int n1, const int n2, const int n3) |
262 |
|
{ |
263 |
< |
const int v0 = ambndx(hp,i,j); |
264 |
< |
const int tflags = (1<<dbit1 | 1<<dbit2); |
265 |
< |
int v1, v2; |
266 |
< |
|
267 |
< |
if ((vflags[v0] & tflags) == tflags) /* is v0 the farthest? */ |
268 |
< |
return(colval(hp->sa[v0].v,CIEY)); |
269 |
< |
v1 = adjacent_verti(hp, i, j, dbit1); |
270 |
< |
if (vflags[v0] & 1<<dbit2) /* v1 farthest if v0>v2 */ |
369 |
< |
return(colval(hp->sa[v1].v,CIEY)); |
370 |
< |
v2 = adjacent_verti(hp, i, j, dbit2); |
371 |
< |
if (vflags[v0] & 1<<dbit1) /* v2 farthest if v0>v1 */ |
372 |
< |
return(colval(hp->sa[v2].v,CIEY)); |
373 |
< |
/* else check if v1>v2 */ |
374 |
< |
if (vflags[v1] & 1<<vdb_edge(dbit1,dbit2)) |
375 |
< |
return(colval(hp->sa[v1].v,CIEY)); |
376 |
< |
return(colval(hp->sa[v2].v,CIEY)); |
263 |
> |
if (hp->sa[n1].d <= hp->sa[n2].d) { |
264 |
> |
if (hp->sa[n1].d <= hp->sa[n3].d) |
265 |
> |
return(colval(hp->sa[n1].v,CIEY)); |
266 |
> |
return(colval(hp->sa[n3].v,CIEY)); |
267 |
> |
} |
268 |
> |
if (hp->sa[n2].d <= hp->sa[n3].d) |
269 |
> |
return(colval(hp->sa[n2].v,CIEY)); |
270 |
> |
return(colval(hp->sa[n3].v,CIEY)); |
271 |
|
} |
272 |
|
|
273 |
|
|
274 |
|
/* Compute vectors and coefficients for Hessian/gradient calcs */ |
275 |
|
static void |
276 |
< |
comp_fftri(FFTRI *ftp, AMBHEMI *hp, int i, int j, int dbit, const uby8 *vflags) |
276 |
> |
comp_fftri(FFTRI *ftp, AMBHEMI *hp, const int n0, const int n1) |
277 |
|
{ |
278 |
< |
const int i0 = ambndx(hp,i,j); |
279 |
< |
double rdot_cp, dot_e, dot_er, rdot_r, rdot_r1, J2; |
386 |
< |
int i1, ii; |
278 |
> |
double rdot_cp, dot_e, dot_er, rdot_r, rdot_r1, J2; |
279 |
> |
int ii; |
280 |
|
|
281 |
< |
ftp->valid = 0; /* check if we can skip this edge */ |
282 |
< |
ii = adjacent_trifl[dbit]; |
283 |
< |
if ((vflags[i0] & ii) == ii) /* cancels if vertex used as value */ |
391 |
< |
return; |
392 |
< |
i1 = adjacent_verti(hp, i, j, dbit); |
393 |
< |
ii = adjacent_trifl[VDB_OPP(dbit)]; |
394 |
< |
if ((vflags[i1] & ii) == ii) /* on either end (for both triangles) */ |
395 |
< |
return; |
396 |
< |
/* else go ahead with calculation */ |
397 |
< |
VSUB(ftp->r_i, hp->sa[i0].p, hp->rp->rop); |
398 |
< |
VSUB(ftp->r_i1, hp->sa[i1].p, hp->rp->rop); |
399 |
< |
VSUB(ftp->e_i, hp->sa[i1].p, hp->sa[i0].p); |
281 |
> |
VSUB(ftp->r_i, hp->sa[n0].p, hp->rp->rop); |
282 |
> |
VSUB(ftp->r_i1, hp->sa[n1].p, hp->rp->rop); |
283 |
> |
VSUB(ftp->e_i, hp->sa[n1].p, hp->sa[n0].p); |
284 |
|
VCROSS(ftp->rcp, ftp->r_i, ftp->r_i1); |
285 |
|
rdot_cp = 1.0/DOT(ftp->rcp,ftp->rcp); |
286 |
|
dot_e = DOT(ftp->e_i,ftp->e_i); |
294 |
|
J2 = ( 0.5*(rdot_r - rdot_r1) - dot_er*ftp->I2 ) / dot_e; |
295 |
|
for (ii = 3; ii--; ) |
296 |
|
ftp->rI2_eJ2[ii] = ftp->I2*ftp->r_i[ii] + J2*ftp->e_i[ii]; |
413 |
– |
ftp->valid++; |
297 |
|
} |
298 |
|
|
299 |
|
|
319 |
|
double d1, d2, d3, d4; |
320 |
|
double I3, J3, K3; |
321 |
|
int i, j; |
439 |
– |
|
440 |
– |
if (!ftp->valid) { /* preemptive test */ |
441 |
– |
memset(hess, 0, sizeof(FVECT)*3); |
442 |
– |
return; |
443 |
– |
} |
322 |
|
/* compute intermediate coefficients */ |
323 |
|
d1 = 1.0/DOT(ftp->r_i,ftp->r_i); |
324 |
|
d2 = 1.0/DOT(ftp->r_i1,ftp->r_i1); |
382 |
|
double f1; |
383 |
|
int i; |
384 |
|
|
507 |
– |
if (!ftp->valid) { /* preemptive test */ |
508 |
– |
memset(grad, 0, sizeof(FVECT)); |
509 |
– |
return; |
510 |
– |
} |
385 |
|
f1 = 2.0*DOT(nrm, ftp->rcp); |
386 |
|
VCROSS(ncp, nrm, ftp->e_i); |
387 |
|
for (i = 3; i--; ) |
411 |
|
|
412 |
|
|
413 |
|
/* Compute anisotropic radii and eigenvector directions */ |
414 |
< |
static int |
414 |
> |
static void |
415 |
|
eigenvectors(FVECT uv[2], float ra[2], FVECT hessian[3]) |
416 |
|
{ |
417 |
|
double hess2[2][2]; |
433 |
|
if (i == 1) /* double-root (circle) */ |
434 |
|
evalue[1] = evalue[0]; |
435 |
|
if (!i || ((evalue[0] = fabs(evalue[0])) <= FTINY*FTINY) | |
436 |
< |
((evalue[1] = fabs(evalue[1])) <= FTINY*FTINY) ) |
437 |
< |
error(INTERNAL, "bad eigenvalue calculation"); |
438 |
< |
|
436 |
> |
((evalue[1] = fabs(evalue[1])) <= FTINY*FTINY) ) { |
437 |
> |
ra[0] = ra[1] = maxarad; |
438 |
> |
return; |
439 |
> |
} |
440 |
|
if (evalue[0] > evalue[1]) { |
441 |
|
ra[0] = sqrt(sqrt(4.0/evalue[0])); |
442 |
|
ra[1] = sqrt(sqrt(4.0/evalue[1])); |
471 |
|
static char memerrmsg[] = "out of memory in ambHessian()"; |
472 |
|
FVECT (*hessrow)[3] = NULL; |
473 |
|
FVECT *gradrow = NULL; |
599 |
– |
uby8 *vflags; |
474 |
|
FVECT hessian[3]; |
475 |
|
FVECT gradient; |
476 |
|
FFTRI fftr; |
492 |
|
error(SYSTEM, memerrmsg); |
493 |
|
memset(gradient, 0, sizeof(gradient)); |
494 |
|
} |
621 |
– |
/* get vertex position flags */ |
622 |
– |
vflags = vertex_flags(hp); |
495 |
|
/* compute first row of edges */ |
496 |
|
for (j = 0; j < hp->ns-1; j++) { |
497 |
< |
comp_fftri(&fftr, hp, 0, j, VDB_X, vflags); |
497 |
> |
comp_fftri(&fftr, hp, AI(hp,0,j), AI(hp,0,j+1)); |
498 |
|
if (hessrow != NULL) |
499 |
|
comp_hessian(hessrow[j], &fftr, hp->rp->ron); |
500 |
|
if (gradrow != NULL) |
504 |
|
for (i = 0; i < hp->ns-1; i++) { |
505 |
|
FVECT hesscol[3]; /* compute first vertical edge */ |
506 |
|
FVECT gradcol; |
507 |
< |
comp_fftri(&fftr, hp, i, 0, VDB_Y, vflags); |
507 |
> |
comp_fftri(&fftr, hp, AI(hp,i,0), AI(hp,i+1,0)); |
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 |
|
double backg; |
516 |
< |
backg = back_ambval(hp, i, j, VDB_X, VDB_Y, vflags); |
516 |
> |
backg = back_ambval(hp, AI(hp,i,j), |
517 |
> |
AI(hp,i,j+1), AI(hp,i+1,j)); |
518 |
|
/* diagonal (inner) edge */ |
519 |
< |
comp_fftri(&fftr, hp, i, j+1, VDB_xY, vflags); |
519 |
> |
comp_fftri(&fftr, hp, AI(hp,i,j+1), AI(hp,i+1,j)); |
520 |
|
if (hessrow != NULL) { |
521 |
|
comp_hessian(hessdia, &fftr, hp->rp->ron); |
522 |
|
rev_hessian(hesscol); |
528 |
|
add2gradient(gradient, gradrow[j], graddia, gradcol, backg); |
529 |
|
} |
530 |
|
/* initialize edge in next row */ |
531 |
< |
comp_fftri(&fftr, hp, i+1, j+1, VDB_x, vflags); |
531 |
> |
comp_fftri(&fftr, hp, AI(hp,i+1,j+1), AI(hp,i+1,j)); |
532 |
|
if (hessrow != NULL) |
533 |
|
comp_hessian(hessrow[j], &fftr, hp->rp->ron); |
534 |
|
if (gradrow != NULL) |
535 |
|
comp_gradient(gradrow[j], &fftr, hp->rp->ron); |
536 |
|
/* new column edge & paired triangle */ |
537 |
< |
backg = back_ambval(hp, i+1, j+1, VDB_x, VDB_y, vflags); |
538 |
< |
comp_fftri(&fftr, hp, i, j+1, VDB_Y, vflags); |
537 |
> |
backg = back_ambval(hp, AI(hp,i+1,j+1), |
538 |
> |
AI(hp,i+1,j), AI(hp,i,j+1)); |
539 |
> |
comp_fftri(&fftr, hp, AI(hp,i,j+1), AI(hp,i+1,j+1)); |
540 |
|
if (hessrow != NULL) { |
541 |
|
comp_hessian(hesscol, &fftr, hp->rp->ron); |
542 |
|
rev_hessian(hessdia); |
556 |
|
/* release row buffers */ |
557 |
|
if (hessrow != NULL) free(hessrow); |
558 |
|
if (gradrow != NULL) free(gradrow); |
685 |
– |
free(vflags); |
559 |
|
|
560 |
|
if (ra != NULL) /* extract eigenvectors & radii */ |
561 |
|
eigenvectors(uv, ra, hessian); |
595 |
|
static uint32 |
596 |
|
ambcorral(AMBHEMI *hp, FVECT uv[2], const double r0, const double r1) |
597 |
|
{ |
598 |
< |
uint32 flgs = 0; |
599 |
< |
int i, j; |
600 |
< |
/* circle around perimeter */ |
598 |
> |
const double max_d = 1.0/(minarad*ambacc + 0.001); |
599 |
> |
const double ang_res = 0.5*PI/(hp->ns-1); |
600 |
> |
const double ang_step = ang_res/((int)(16/PI*ang_res) + (1+FTINY)); |
601 |
> |
double avg_d = 0; |
602 |
> |
uint32 flgs = 0; |
603 |
> |
int i, j; |
604 |
> |
/* don't bother for a few samples */ |
605 |
> |
if (hp->ns < 12) |
606 |
> |
return(0); |
607 |
> |
/* check distances overhead */ |
608 |
> |
for (i = hp->ns*3/4; i-- > hp->ns>>2; ) |
609 |
> |
for (j = hp->ns*3/4; j-- > hp->ns>>2; ) |
610 |
> |
avg_d += ambsam(hp,i,j).d; |
611 |
> |
avg_d *= 4.0/(hp->ns*hp->ns); |
612 |
> |
if (avg_d*r0 >= 1.0) /* ceiling too low for corral? */ |
613 |
> |
return(0); |
614 |
> |
if (avg_d >= max_d) /* insurance */ |
615 |
> |
return(0); |
616 |
> |
/* else circle around perimeter */ |
617 |
|
for (i = 0; i < hp->ns; i++) |
618 |
|
for (j = 0; j < hp->ns; j += !i|(i==hp->ns-1) ? 1 : hp->ns-1) { |
619 |
|
AMBSAMP *ap = &ambsam(hp,i,j); |
620 |
|
FVECT vec; |
621 |
|
double u, v; |
622 |
< |
double ang; |
622 |
> |
double ang, a1; |
623 |
|
int abp; |
624 |
< |
if (ap->d <= FTINY) |
625 |
< |
continue; |
624 |
> |
if ((ap->d <= FTINY) | (ap->d >= max_d)) |
625 |
> |
continue; /* too far or too near */ |
626 |
|
VSUB(vec, ap->p, hp->rp->rop); |
627 |
|
u = DOT(vec, uv[0]) * ap->d; |
628 |
|
v = DOT(vec, uv[1]) * ap->d; |
629 |
|
if ((r0*r0*u*u + r1*r1*v*v) * ap->d*ap->d <= 1.0) |
630 |
|
continue; /* occluder outside ellipse */ |
631 |
|
ang = atan2a(v, u); /* else set direction flags */ |
632 |
< |
ang += 2.0*PI*(ang < 0); |
633 |
< |
ang *= 16./PI; |
745 |
< |
if ((ang < .5) | (ang >= 31.5)) |
746 |
< |
flgs |= 0x80000001; |
747 |
< |
else |
748 |
< |
flgs |= 3L<<(int)(ang-.5); |
632 |
> |
for (a1 = ang-.5*ang_res; a1 <= ang+.5*ang_res; a1 += ang_step) |
633 |
> |
flgs |= 1L<<(int)(16/PI*(a1 + 2.*PI*(a1 < 0))); |
634 |
|
} |
635 |
|
return(flgs); |
636 |
|
} |
687 |
|
return(-1); /* return value w/o Hessian */ |
688 |
|
} |
689 |
|
cnt = ambssamp*wt + 0.5; /* perform super-sampling? */ |
690 |
< |
if (cnt > 0) |
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)) { |
701 |
|
K = 1.0; |
702 |
|
pg = NULL; |
703 |
|
dg = NULL; |
704 |
+ |
crlp = NULL; |
705 |
|
} |
706 |
|
ap = hp->sa; /* relative Y channel from here on... */ |
707 |
|
for (i = hp->ns*hp->ns; i--; ap++) |
737 |
|
if (ra[0] > maxarad) |
738 |
|
ra[0] = maxarad; |
739 |
|
} |
740 |
< |
if (crlp != NULL) /* flag encroached directions */ |
740 |
> |
/* flag encroached directions */ |
741 |
> |
if ((wt >= 0.5-FTINY) & (crlp != NULL)) |
742 |
|
*crlp = ambcorral(hp, uv, ra[0]*ambacc, ra[1]*ambacc); |
743 |
|
if (pg != NULL) { /* cap gradient if necessary */ |
744 |
|
d = pg[0]*pg[0]*ra[0]*ra[0] + pg[1]*pg[1]*ra[1]*ra[1]; |