96 |
|
{ |
97 |
|
struct s_ambsamp *ap = &ambsamp(hp,i,j); |
98 |
|
RAY ar; |
99 |
– |
int hlist[3]; |
99 |
|
double spt[2], zd; |
100 |
|
int ii; |
101 |
|
/* ambient coefficient for weight */ |
127 |
|
multcolor(ar.rcol, ar.rcoef); /* apply coefficient */ |
128 |
|
copycolor(ap->v, ar.rcol); |
129 |
|
if (ar.rt > 20.0*maxarad) /* limit vertex distance */ |
130 |
< |
ar.rt = 20.0*maxarad; |
131 |
< |
VSUM(ap->p, ar.rorg, ar.rdir, ar.rt); |
130 |
> |
VSUM(ap->p, ar.rorg, ar.rdir, 20.0*maxarad); |
131 |
> |
else |
132 |
> |
VCOPY(ap->p, ar.rop); |
133 |
|
return(ap); |
134 |
|
} |
135 |
|
|
146 |
|
VSUB(ftp->e_i, ap1, ap0); |
147 |
|
VCROSS(v1, ftp->e_i, ftp->r_i); |
148 |
|
ftp->nf = 1.0/DOT(v1,v1); |
149 |
– |
VCROSS(v1, ftp->r_i, ftp->r_i1); |
150 |
– |
ftp->I1 = sqrt(DOT(v1,v1)*ftp->nf); |
149 |
|
dot_e = DOT(ftp->e_i,ftp->e_i); |
150 |
|
dot_er = DOT(ftp->e_i, ftp->r_i); |
151 |
|
dot_r = DOT(ftp->r_i,ftp->r_i); |
152 |
|
dot_r1 = DOT(ftp->r_i1,ftp->r_i1); |
153 |
+ |
ftp->I1 = acos( DOT(ftp->r_i, ftp->r_i1) / sqrt(dot_r*dot_r1) ) * |
154 |
+ |
sqrt( ftp->nf ); |
155 |
|
ftp->I2 = ( DOT(ftp->e_i, ftp->r_i1)/dot_r1 - dot_er/dot_r + |
156 |
|
dot_e*ftp->I1 )*0.5*ftp->nf; |
157 |
< |
ftp->J2 = 0.25*ftp->nf*( 1.0/dot_r - 1.0/dot_r1 ) - |
157 |
> |
ftp->J2 = 0.5/dot_e*( 1.0/dot_r - 1.0/dot_r1 ) - |
158 |
|
dot_er/dot_e*ftp->I2; |
159 |
|
} |
160 |
|
|
187 |
|
d3 = 1.0/DOT(ftp->e_i,ftp->e_i); |
188 |
|
d4 = DOT(ftp->e_i, ftp->r_i); |
189 |
|
I3 = 0.25*ftp->nf*( DOT(ftp->e_i, ftp->r_i1)*d2*d2 - d4*d1*d1 + |
190 |
< |
3.0*d3*ftp->I2 ); |
190 |
> |
3.0/d3*ftp->I2 ); |
191 |
|
J3 = 0.25*d3*(d1*d1 - d2*d2) - d4*d3*I3; |
192 |
|
K3 = d3*(ftp->I2 - I3/d1 - 2.0*d4*J3); |
193 |
|
/* intermediate matrices */ |
288 |
|
|
289 |
|
VSUB(vec, ap1->p, orig); |
290 |
|
d2best = DOT(vec,vec); |
291 |
< |
vback = ap1->v[CIEY]; |
291 |
> |
vback = colval(ap1->v,CIEY); |
292 |
|
VSUB(vec, ap2->p, orig); |
293 |
|
d2 = DOT(vec,vec); |
294 |
|
if (d2 > d2best) { |
295 |
|
d2best = d2; |
296 |
< |
vback = ap2->v[CIEY]; |
296 |
> |
vback = colval(ap2->v,CIEY); |
297 |
|
} |
298 |
|
VSUB(vec, ap3->p, orig); |
299 |
|
d2 = DOT(vec,vec); |
300 |
|
if (d2 > d2best) |
301 |
< |
return(ap3->v[CIEY]); |
301 |
> |
return(colval(ap3->v,CIEY)); |
302 |
|
return(vback); |
303 |
|
} |
304 |
|
|
328 |
|
error(INTERNAL, "bad eigenvalue calculation"); |
329 |
|
|
330 |
|
if (evalue[0] > evalue[1]) { |
331 |
< |
ra[0] = 1.0/sqrt(sqrt(evalue[0])); |
332 |
< |
ra[1] = 1.0/sqrt(sqrt(evalue[1])); |
331 |
> |
ra[0] = sqrt(sqrt(4.0/evalue[0])); |
332 |
> |
ra[1] = sqrt(sqrt(4.0/evalue[1])); |
333 |
|
slope1 = evalue[1]; |
334 |
|
} else { |
335 |
< |
ra[0] = 1.0/sqrt(sqrt(evalue[1])); |
336 |
< |
ra[1] = 1.0/sqrt(sqrt(evalue[0])); |
335 |
> |
ra[0] = sqrt(sqrt(4.0/evalue[1])); |
336 |
> |
ra[1] = sqrt(sqrt(4.0/evalue[0])); |
337 |
|
slope1 = evalue[0]; |
338 |
|
} |
339 |
|
/* compute unit eigenvectors */ |
454 |
|
|
455 |
|
if (ra != NULL) /* extract eigenvectors & radii */ |
456 |
|
eigenvectors(uv, ra, hessian); |
457 |
< |
if (pg != NULL) { /* project position gradient */ |
458 |
< |
pg[0] = DOT(gradient, uv[0]); |
459 |
< |
pg[1] = DOT(gradient, uv[1]); |
457 |
> |
if (pg != NULL) { /* tangential position gradient/PI */ |
458 |
> |
pg[0] = DOT(gradient, uv[0]) / PI; |
459 |
> |
pg[1] = DOT(gradient, uv[1]) / PI; |
460 |
|
} |
461 |
|
} |
462 |
|
|
466 |
|
ambdirgrad(AMBHEMI *hp, FVECT uv[2], float dg[2]) |
467 |
|
{ |
468 |
|
struct s_ambsamp *ap; |
469 |
+ |
double dgsum[2]; |
470 |
|
int n; |
471 |
|
FVECT vd; |
472 |
|
double gfact; |
473 |
|
|
474 |
< |
dg[0] = dg[1] = 0; |
474 |
> |
dgsum[0] = dgsum[1] = 0.0; /* sum values times -tan(theta) */ |
475 |
|
for (ap = hp->sa, n = hp->ns*hp->ns; n--; ap++) { |
476 |
|
/* use vector for azimuth + 90deg */ |
477 |
|
VSUB(vd, ap->p, hp->rp->rop); |
478 |
< |
/* brightness with tangent factor */ |
479 |
< |
gfact = ap->v[CIEY] / DOT(hp->rp->ron, vd); |
480 |
< |
/* sine = proj_radius/vd_length */ |
481 |
< |
dg[0] -= DOT(uv[1], vd) * gfact; |
482 |
< |
dg[1] += DOT(uv[0], vd) * gfact; |
478 |
> |
/* brightness over cosine factor */ |
479 |
> |
gfact = colval(ap->v,CIEY) / DOT(hp->rp->ron, vd); |
480 |
> |
/* -sine = -proj_radius/vd_length */ |
481 |
> |
dgsum[0] += DOT(uv[1], vd) * gfact; |
482 |
> |
dgsum[1] -= DOT(uv[0], vd) * gfact; |
483 |
|
} |
484 |
+ |
dg[0] = dgsum[0] / (hp->ns*hp->ns); |
485 |
+ |
dg[1] = dgsum[1] / (hp->ns*hp->ns); |
486 |
|
} |
487 |
|
|
488 |
|
|
527 |
|
free(hp); |
528 |
|
return(0); /* no valid samples */ |
529 |
|
} |
530 |
< |
d = 1.0 / cnt; /* final indirect irradiance/PI */ |
528 |
< |
acol[0] *= d; acol[1] *= d; acol[2] *= d; |
529 |
< |
copycolor(rcol, acol); |
530 |
> |
copycolor(rcol, acol); /* final indirect irradiance/PI */ |
531 |
|
if (cnt < hp->ns*hp->ns || /* incomplete sampling? */ |
532 |
|
(ra == NULL) & (pg == NULL) & (dg == NULL)) { |
533 |
|
free(hp); |
534 |
|
return(-1); /* no radius or gradient calc. */ |
535 |
|
} |
536 |
< |
d = 0.01 * bright(rcol); /* add in 1% before Hessian comp. */ |
537 |
< |
if (d < FTINY) d = FTINY; |
538 |
< |
ap = hp->sa; /* using Y channel from here on... */ |
536 |
> |
multcolor(acol, hp->acoef); /* normalize Y values */ |
537 |
> |
if ((d = bright(acol)) > FTINY) |
538 |
> |
d = 1.0/d; |
539 |
> |
else |
540 |
> |
d = 0.0; |
541 |
> |
ap = hp->sa; /* relative Y channel from here on... */ |
542 |
|
for (i = hp->ns*hp->ns; i--; ap++) |
543 |
< |
colval(ap->v,CIEY) = bright(ap->v) + d; |
543 |
> |
colval(ap->v,CIEY) = bright(ap->v)*d + 0.0314; |
544 |
|
|
545 |
|
if (uv == NULL) /* make sure we have axis pointers */ |
546 |
|
uv = my_uv; |
547 |
|
/* compute radii & pos. gradient */ |
548 |
|
ambHessian(hp, uv, ra, pg); |
549 |
+ |
|
550 |
|
if (dg != NULL) /* compute direction gradient */ |
551 |
|
ambdirgrad(hp, uv, dg); |
552 |
+ |
|
553 |
|
if (ra != NULL) { /* scale/clamp radii */ |
554 |
< |
d = sqrt(sqrt((4.0/PI)*bright(rcol)/wt)); |
555 |
< |
ra[0] *= d; |
554 |
> |
if (ra[0] < minarad) { |
555 |
> |
ra[0] = minarad; |
556 |
> |
if (ra[1] < minarad) |
557 |
> |
ra[1] = minarad; |
558 |
> |
} |
559 |
> |
ra[0] *= d = 1.0/sqrt(sqrt(wt)); |
560 |
|
if ((ra[1] *= d) > 2.0*ra[0]) |
561 |
|
ra[1] = 2.0*ra[0]; |
562 |
|
if (ra[1] > maxarad) { |