28 |
|
COLOR acoef; /* division contribution coefficient */ |
29 |
|
struct s_ambsamp { |
30 |
|
COLOR v; /* hemisphere sample value */ |
31 |
< |
float p[3]; /* intersection point */ |
31 |
> |
FVECT p; /* intersection point */ |
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)] |
36 |
|
|
37 |
|
typedef struct { |
38 |
< |
FVECT r_i, r_i1, e_i; |
39 |
< |
double nf, I1, I2, J2; |
38 |
> |
FVECT r_i, r_i1, e_i, rI2_eJ2; |
39 |
> |
double nf, I1, I2; |
40 |
|
} FFTRI; /* vectors and coefficients for Hessian calculation */ |
41 |
|
|
42 |
|
|
69 |
|
copycolor(hp->acoef, ac); |
70 |
|
d = 1.0/(n*n); |
71 |
|
scalecolor(hp->acoef, d); |
72 |
< |
/* make tangent axes */ |
73 |
< |
hp->uy[0] = hp->uy[1] = hp->uy[2] = 0.0; |
72 |
> |
/* make tangent plane axes */ |
73 |
> |
hp->uy[0] = 0.1 - 0.2*frandom(); |
74 |
> |
hp->uy[1] = 0.1 - 0.2*frandom(); |
75 |
> |
hp->uy[2] = 0.1 - 0.2*frandom(); |
76 |
|
for (i = 0; i < 3; i++) |
77 |
|
if (r->ron[i] < 0.6 && r->ron[i] > -0.6) |
78 |
|
break; |
87 |
|
} |
88 |
|
|
89 |
|
|
90 |
< |
static int |
90 |
> |
static struct s_ambsamp * |
91 |
|
ambsample( /* sample an ambient direction */ |
92 |
|
AMBHEMI *hp, |
93 |
|
int i, |
96 |
|
{ |
97 |
|
struct s_ambsamp *ap = &ambsamp(hp,i,j); |
98 |
|
RAY ar; |
97 |
– |
int hlist[3]; |
99 |
|
double spt[2], zd; |
100 |
|
int ii; |
101 |
|
/* ambient coefficient for weight */ |
103 |
|
setcolor(ar.rcoef, AVGREFL, AVGREFL, AVGREFL); |
104 |
|
else |
105 |
|
copycolor(ar.rcoef, hp->acoef); |
106 |
< |
if (rayorigin(&ar, AMBIENT, hp->rp, ar.rcoef) < 0) { |
107 |
< |
setcolor(ap->v, 0., 0., 0.); |
107 |
< |
VCOPY(ap->p, hp->rp->rop); |
108 |
< |
return(0); /* no sample taken */ |
109 |
< |
} |
106 |
> |
if (rayorigin(&ar, AMBIENT, hp->rp, ar.rcoef) < 0) |
107 |
> |
goto badsample; |
108 |
|
if (ambacc > FTINY) { |
109 |
|
multcolor(ar.rcoef, hp->acoef); |
110 |
|
scalecolor(ar.rcoef, 1./AVGREFL); |
111 |
|
} |
112 |
|
/* generate hemispherical sample */ |
113 |
|
SDsquare2disk(spt, (i+.1+.8*frandom())/hp->ns, |
114 |
< |
(j+.1+.8*frandom())/hp->ns); |
114 |
> |
(j+.1+.8*frandom())/hp->ns ); |
115 |
|
zd = sqrt(1. - spt[0]*spt[0] - spt[1]*spt[1]); |
116 |
|
for (ii = 3; ii--; ) |
117 |
|
ar.rdir[ii] = spt[0]*hp->ux[ii] + |
121 |
|
dimlist[ndims++] = i*hp->ns + j + 90171; |
122 |
|
rayvalue(&ar); /* evaluate ray */ |
123 |
|
ndims--; |
124 |
+ |
/* limit vertex distance */ |
125 |
+ |
if (ar.rt > 10.0*thescene.cusize) |
126 |
+ |
ar.rt = 10.0*thescene.cusize; |
127 |
+ |
else if (ar.rt <= FTINY) /* should never happen! */ |
128 |
+ |
goto badsample; |
129 |
+ |
VSUM(ap->p, ar.rorg, ar.rdir, ar.rt); |
130 |
|
multcolor(ar.rcol, ar.rcoef); /* apply coefficient */ |
131 |
|
copycolor(ap->v, ar.rcol); |
132 |
< |
if (ar.rt > 20.0*maxarad) /* limit vertex distance */ |
133 |
< |
ar.rt = 20.0*maxarad; |
134 |
< |
VSUM(ap->p, ar.rorg, ar.rdir, ar.rt); |
135 |
< |
return(1); |
132 |
> |
return(ap); |
133 |
> |
badsample: |
134 |
> |
setcolor(ap->v, 0., 0., 0.); |
135 |
> |
VCOPY(ap->p, hp->rp->rop); |
136 |
> |
return(NULL); |
137 |
|
} |
138 |
|
|
139 |
|
|
140 |
|
/* Compute vectors and coefficients for Hessian/gradient calcs */ |
141 |
|
static void |
142 |
< |
comp_fftri(FFTRI *ftp, float ap0[3], float ap1[3], FVECT rop) |
142 |
> |
comp_fftri(FFTRI *ftp, FVECT ap0, FVECT ap1, FVECT rop) |
143 |
|
{ |
144 |
< |
FVECT v1; |
145 |
< |
double dot_e, dot_er, dot_r, dot_r1; |
144 |
> |
FVECT vcp; |
145 |
> |
double dot_e, dot_er, rdot_r, rdot_r1, J2; |
146 |
> |
int i; |
147 |
|
|
148 |
|
VSUB(ftp->r_i, ap0, rop); |
149 |
|
VSUB(ftp->r_i1, ap1, rop); |
150 |
|
VSUB(ftp->e_i, ap1, ap0); |
151 |
< |
VCROSS(v1, ftp->e_i, ftp->r_i); |
152 |
< |
ftp->nf = 1.0/DOT(v1,v1); |
147 |
< |
VCROSS(v1, ftp->r_i, ftp->r_i1); |
148 |
< |
ftp->I1 = sqrt(DOT(v1,v1)*ftp->nf); |
151 |
> |
VCROSS(vcp, ftp->e_i, ftp->r_i); |
152 |
> |
ftp->nf = 1.0/DOT(vcp,vcp); |
153 |
|
dot_e = DOT(ftp->e_i,ftp->e_i); |
154 |
|
dot_er = DOT(ftp->e_i, ftp->r_i); |
155 |
< |
dot_r = DOT(ftp->r_i,ftp->r_i); |
156 |
< |
dot_r1 = DOT(ftp->r_i1,ftp->r_i1); |
157 |
< |
ftp->I2 = ( DOT(ftp->e_i, ftp->r_i1)/dot_r1 - dot_er/dot_r + |
155 |
> |
rdot_r = 1.0/DOT(ftp->r_i,ftp->r_i); |
156 |
> |
rdot_r1 = 1.0/DOT(ftp->r_i1,ftp->r_i1); |
157 |
> |
ftp->I1 = acos( DOT(ftp->r_i, ftp->r_i1) * sqrt(rdot_r*rdot_r1) ) * |
158 |
> |
sqrt( ftp->nf ); |
159 |
> |
ftp->I2 = ( DOT(ftp->e_i, ftp->r_i1)*rdot_r1 - dot_er*rdot_r + |
160 |
|
dot_e*ftp->I1 )*0.5*ftp->nf; |
161 |
< |
ftp->J2 = 0.25*ftp->nf*( 1.0/dot_r - 1.0/dot_r1 ) - |
162 |
< |
dot_er/dot_e*ftp->I2; |
161 |
> |
J2 = ( 0.5*(rdot_r - rdot_r1) - dot_er*ftp->I2 ) / dot_e; |
162 |
> |
for (i = 3; i--; ) |
163 |
> |
ftp->rI2_eJ2[i] = ftp->I2*ftp->r_i[i] + J2*ftp->e_i[i]; |
164 |
|
} |
165 |
|
|
166 |
|
|
167 |
< |
/* Compose matrix from two vectors */ |
167 |
> |
/* Compose 3x3 matrix from two vectors */ |
168 |
|
static void |
169 |
|
compose_matrix(FVECT mat[3], FVECT va, FVECT vb) |
170 |
|
{ |
181 |
|
static void |
182 |
|
comp_hessian(FVECT hess[3], FFTRI *ftp, FVECT nrm) |
183 |
|
{ |
184 |
< |
FVECT v1, v2; |
184 |
> |
FVECT vcp; |
185 |
|
FVECT m1[3], m2[3], m3[3], m4[3]; |
186 |
|
double d1, d2, d3, d4; |
187 |
|
double I3, J3, K3; |
192 |
|
d3 = 1.0/DOT(ftp->e_i,ftp->e_i); |
193 |
|
d4 = DOT(ftp->e_i, ftp->r_i); |
194 |
|
I3 = 0.25*ftp->nf*( DOT(ftp->e_i, ftp->r_i1)*d2*d2 - d4*d1*d1 + |
195 |
< |
3.0*ftp->I2*d3 ); |
195 |
> |
3.0/d3*ftp->I2 ); |
196 |
|
J3 = 0.25*d3*(d1*d1 - d2*d2) - d4*d3*I3; |
197 |
|
K3 = d3*(ftp->I2 - I3/d1 - 2.0*d4*J3); |
198 |
|
/* intermediate matrices */ |
199 |
< |
VCROSS(v1, nrm, ftp->e_i); |
200 |
< |
for (j = 3; j--; ) |
194 |
< |
v2[i] = ftp->I2*ftp->r_i[j] + ftp->J2*ftp->e_i[j]; |
195 |
< |
compose_matrix(m1, v1, v2); |
199 |
> |
VCROSS(vcp, nrm, ftp->e_i); |
200 |
> |
compose_matrix(m1, vcp, ftp->rI2_eJ2); |
201 |
|
compose_matrix(m2, ftp->r_i, ftp->r_i); |
202 |
|
compose_matrix(m3, ftp->e_i, ftp->e_i); |
203 |
|
compose_matrix(m4, ftp->r_i, ftp->e_i); |
204 |
< |
VCROSS(v1, ftp->r_i, ftp->e_i); |
205 |
< |
d1 = DOT(nrm, v1); |
204 |
> |
VCROSS(vcp, ftp->r_i, ftp->e_i); |
205 |
> |
d1 = DOT(nrm, vcp); |
206 |
|
d2 = -d1*ftp->I2; |
207 |
|
d1 *= 2.0; |
208 |
|
for (i = 3; i--; ) /* final matrix sum */ |
210 |
|
hess[i][j] = m1[i][j] + d1*( I3*m2[i][j] + K3*m3[i][j] + |
211 |
|
2.0*J3*m4[i][j] ); |
212 |
|
hess[i][j] += d2*(i==j); |
213 |
< |
hess[i][j] *= -1.0/PI; |
213 |
> |
hess[i][j] *= 1.0/PI; |
214 |
|
} |
215 |
|
} |
216 |
|
|
254 |
|
f1 = 2.0*DOT(nrm, vcp); |
255 |
|
VCROSS(vcp, nrm, ftp->e_i); |
256 |
|
for (i = 3; i--; ) |
257 |
< |
grad[i] = (0.5/PI)*( ftp->I1*vcp[i] + |
253 |
< |
f1*(ftp->I2*ftp->r_i[i] + ftp->J2*ftp->e_i[i]) ); |
257 |
> |
grad[i] = (-0.5/PI)*( ftp->I1*vcp[i] + f1*ftp->rI2_eJ2[i] ); |
258 |
|
} |
259 |
|
|
260 |
|
|
290 |
|
|
291 |
|
VSUB(vec, ap1->p, orig); |
292 |
|
d2best = DOT(vec,vec); |
293 |
< |
vback = ap1->v[CIEY]; |
293 |
> |
vback = colval(ap1->v,CIEY); |
294 |
|
VSUB(vec, ap2->p, orig); |
295 |
|
d2 = DOT(vec,vec); |
296 |
|
if (d2 > d2best) { |
297 |
|
d2best = d2; |
298 |
< |
vback = ap2->v[CIEY]; |
298 |
> |
vback = colval(ap2->v,CIEY); |
299 |
|
} |
300 |
|
VSUB(vec, ap3->p, orig); |
301 |
|
d2 = DOT(vec,vec); |
302 |
|
if (d2 > d2best) |
303 |
< |
return(ap3->v[CIEY]); |
303 |
> |
return(colval(ap3->v,CIEY)); |
304 |
|
return(vback); |
305 |
|
} |
306 |
|
|
323 |
|
hess2[1][0] = DOT(uv[1], a); |
324 |
|
hess2[1][1] = DOT(uv[1], b); |
325 |
|
/* compute eigenvalues */ |
326 |
< |
if (quadratic(evalue, 1.0, -hess2[0][0]-hess2[1][1], |
326 |
> |
if ( quadratic(evalue, 1.0, -hess2[0][0]-hess2[1][1], |
327 |
|
hess2[0][0]*hess2[1][1]-hess2[0][1]*hess2[1][0]) != 2 || |
328 |
< |
(evalue[0] = fabs(evalue[0])) <= FTINY*FTINY*FTINY || |
329 |
< |
(evalue[1] = fabs(evalue[1])) <= FTINY*FTINY*FTINY) |
328 |
> |
(evalue[0] = fabs(evalue[0])) <= FTINY*FTINY || |
329 |
> |
(evalue[1] = fabs(evalue[1])) <= FTINY*FTINY ) |
330 |
|
error(INTERNAL, "bad eigenvalue calculation"); |
331 |
|
|
332 |
|
if (evalue[0] > evalue[1]) { |
333 |
< |
ra[0] = 1.0/sqrt(sqrt(evalue[0])); |
334 |
< |
ra[1] = 1.0/sqrt(sqrt(evalue[1])); |
333 |
> |
ra[0] = sqrt(sqrt(4.0/evalue[0])); |
334 |
> |
ra[1] = sqrt(sqrt(4.0/evalue[1])); |
335 |
|
slope1 = evalue[1]; |
336 |
|
} else { |
337 |
< |
ra[0] = 1.0/sqrt(sqrt(evalue[1])); |
338 |
< |
ra[1] = 1.0/sqrt(sqrt(evalue[0])); |
337 |
> |
ra[0] = sqrt(sqrt(4.0/evalue[1])); |
338 |
> |
ra[1] = sqrt(sqrt(4.0/evalue[0])); |
339 |
|
slope1 = evalue[0]; |
340 |
|
} |
341 |
|
/* compute unit eigenvectors */ |
356 |
|
ambHessian( /* anisotropic radii & pos. gradient */ |
357 |
|
AMBHEMI *hp, |
358 |
|
FVECT uv[2], /* returned */ |
359 |
< |
float ra[2], /* returned */ |
360 |
< |
float pg[2] /* returned */ |
359 |
> |
float ra[2], /* returned (optional) */ |
360 |
> |
float pg[2] /* returned (optional) */ |
361 |
|
) |
362 |
|
{ |
363 |
|
static char memerrmsg[] = "out of memory in ambHessian()"; |
372 |
|
VCOPY(uv[1], hp->uy); |
373 |
|
/* clock-wise vertex traversal from sample POV */ |
374 |
|
if (ra != NULL) { /* initialize Hessian row buffer */ |
375 |
< |
hessrow = (FVECT (*)[3])malloc(sizeof(FVECT)*3*hp->ns); |
375 |
> |
hessrow = (FVECT (*)[3])malloc(sizeof(FVECT)*3*(hp->ns-1)); |
376 |
|
if (hessrow == NULL) |
377 |
|
error(SYSTEM, memerrmsg); |
378 |
|
memset(hessian, 0, sizeof(hessian)); |
379 |
|
} else if (pg == NULL) /* bogus call? */ |
380 |
|
return; |
381 |
|
if (pg != NULL) { /* initialize form factor row buffer */ |
382 |
< |
gradrow = (FVECT *)malloc(sizeof(FVECT)*hp->ns); |
382 |
> |
gradrow = (FVECT *)malloc(sizeof(FVECT)*(hp->ns-1)); |
383 |
|
if (gradrow == NULL) |
384 |
|
error(SYSTEM, memerrmsg); |
385 |
|
memset(gradient, 0, sizeof(gradient)); |
456 |
|
|
457 |
|
if (ra != NULL) /* extract eigenvectors & radii */ |
458 |
|
eigenvectors(uv, ra, hessian); |
459 |
< |
if (pg != NULL) { /* project position gradient */ |
459 |
> |
if (pg != NULL) { /* tangential position gradient */ |
460 |
|
pg[0] = DOT(gradient, uv[0]); |
461 |
|
pg[1] = DOT(gradient, uv[1]); |
462 |
|
} |
468 |
|
ambdirgrad(AMBHEMI *hp, FVECT uv[2], float dg[2]) |
469 |
|
{ |
470 |
|
struct s_ambsamp *ap; |
471 |
+ |
double dgsum[2]; |
472 |
|
int n; |
473 |
+ |
FVECT vd; |
474 |
+ |
double gfact; |
475 |
|
|
476 |
< |
dg[0] = dg[1] = 0; |
476 |
> |
dgsum[0] = dgsum[1] = 0.0; /* sum values times -tan(theta) */ |
477 |
|
for (ap = hp->sa, n = hp->ns*hp->ns; n--; ap++) { |
471 |
– |
FVECT vd; |
472 |
– |
double gfact; |
478 |
|
/* use vector for azimuth + 90deg */ |
479 |
|
VSUB(vd, ap->p, hp->rp->rop); |
480 |
< |
/* brightness with tangent factor */ |
481 |
< |
gfact = ap->v[CIEY] / DOT(hp->rp->ron, vd); |
482 |
< |
/* sine = proj_radius/vd_length */ |
483 |
< |
dg[0] -= DOT(uv[1], vd) * gfact ; |
484 |
< |
dg[1] += DOT(uv[0], vd) * gfact; |
480 |
> |
/* brightness over cosine factor */ |
481 |
> |
gfact = colval(ap->v,CIEY) / DOT(hp->rp->ron, vd); |
482 |
> |
/* -sine = -proj_radius/vd_length */ |
483 |
> |
dgsum[0] += DOT(uv[1], vd) * gfact; |
484 |
> |
dgsum[1] -= DOT(uv[0], vd) * gfact; |
485 |
|
} |
486 |
+ |
dg[0] = dgsum[0] / (hp->ns*hp->ns); |
487 |
+ |
dg[1] = dgsum[1] / (hp->ns*hp->ns); |
488 |
|
} |
489 |
|
|
490 |
|
|
499 |
|
float dg[2] /* returned (optional) */ |
500 |
|
) |
501 |
|
{ |
502 |
+ |
AMBHEMI *hp = inithemi(rcol, r, wt); |
503 |
|
int cnt = 0; |
504 |
|
FVECT my_uv[2]; |
497 |
– |
AMBHEMI *hp; |
505 |
|
double d, acol[3]; |
506 |
|
struct s_ambsamp *ap; |
507 |
|
int i, j; |
508 |
< |
/* initialize */ |
509 |
< |
if ((hp = inithemi(rcol, r, wt)) == NULL) |
508 |
> |
/* check/initialize */ |
509 |
> |
if (hp == NULL) |
510 |
|
return(0); |
511 |
|
if (uv != NULL) |
512 |
|
memset(uv, 0, sizeof(FVECT)*2); |
520 |
|
acol[0] = acol[1] = acol[2] = 0.0; |
521 |
|
for (i = hp->ns; i--; ) |
522 |
|
for (j = hp->ns; j--; ) |
523 |
< |
if (ambsample(hp, i, j)) { |
517 |
< |
ap = &ambsamp(hp,i,j); |
523 |
> |
if ((ap = ambsample(hp, i, j)) != NULL) { |
524 |
|
addcolor(acol, ap->v); |
525 |
|
++cnt; |
526 |
|
} |
529 |
|
free(hp); |
530 |
|
return(0); /* no valid samples */ |
531 |
|
} |
532 |
< |
d = 1.0 / cnt; /* final indirect irradiance/PI */ |
527 |
< |
acol[0] *= d; acol[1] *= d; acol[2] *= d; |
528 |
< |
copycolor(rcol, acol); |
532 |
> |
copycolor(rcol, acol); /* final indirect irradiance/PI */ |
533 |
|
if (cnt < hp->ns*hp->ns || /* incomplete sampling? */ |
534 |
|
(ra == NULL) & (pg == NULL) & (dg == NULL)) { |
535 |
|
free(hp); |
536 |
|
return(-1); /* no radius or gradient calc. */ |
537 |
|
} |
538 |
< |
d = 0.01 * bright(rcol); /* add in 1% before Hessian comp. */ |
539 |
< |
if (d < FTINY) d = FTINY; |
540 |
< |
ap = hp->sa; /* using Y channel from here on... */ |
538 |
> |
if (bright(acol) > FTINY) /* normalize Y values */ |
539 |
> |
d = cnt/bright(acol); |
540 |
> |
else |
541 |
> |
d = 0.0; |
542 |
> |
ap = hp->sa; /* relative Y channel from here on... */ |
543 |
|
for (i = hp->ns*hp->ns; i--; ap++) |
544 |
< |
colval(ap->v,CIEY) = bright(ap->v) + d; |
544 |
> |
colval(ap->v,CIEY) = bright(ap->v)*d + 0.01; |
545 |
|
|
546 |
|
if (uv == NULL) /* make sure we have axis pointers */ |
547 |
|
uv = my_uv; |
548 |
|
/* compute radii & pos. gradient */ |
549 |
|
ambHessian(hp, uv, ra, pg); |
550 |
+ |
|
551 |
|
if (dg != NULL) /* compute direction gradient */ |
552 |
|
ambdirgrad(hp, uv, dg); |
553 |
< |
if (ra != NULL) { /* adjust/clamp radii */ |
554 |
< |
d = sqrt(sqrt((4.0/PI)*bright(rcol)/wt)); |
555 |
< |
if ((ra[0] *= d) > maxarad) |
556 |
< |
ra[0] = maxarad; |
553 |
> |
|
554 |
> |
if (ra != NULL) { /* scale/clamp radii */ |
555 |
> |
if (ra[0] < minarad) { |
556 |
> |
ra[0] = minarad; |
557 |
> |
if (ra[1] < minarad) |
558 |
> |
ra[1] = minarad; |
559 |
> |
/* cap gradient if necessary */ |
560 |
> |
if (pg != NULL) { |
561 |
> |
d = pg[0]*pg[0]*ra[0]*ra[0] + |
562 |
> |
pg[1]*pg[1]*ra[1]*ra[1]; |
563 |
> |
if (d > 1.0) { |
564 |
> |
d = 1.0/sqrt(d); |
565 |
> |
pg[0] *= d; |
566 |
> |
pg[1] *= d; |
567 |
> |
} |
568 |
> |
} |
569 |
> |
} |
570 |
> |
ra[0] *= d = 1.0/sqrt(sqrt(wt)); |
571 |
|
if ((ra[1] *= d) > 2.0*ra[0]) |
572 |
|
ra[1] = 2.0*ra[0]; |
573 |
+ |
if (ra[1] > maxarad) { |
574 |
+ |
ra[1] = maxarad; |
575 |
+ |
if (ra[0] > maxarad) |
576 |
+ |
ra[0] = maxarad; |
577 |
+ |
} |
578 |
|
} |
579 |
|
free(hp); /* clean up and return */ |
580 |
|
return(1); |