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, |
107 |
|
if (rayorigin(&ar, AMBIENT, hp->rp, ar.rcoef) < 0) { |
108 |
|
setcolor(ap->v, 0., 0., 0.); |
109 |
|
VCOPY(ap->p, hp->rp->rop); |
110 |
< |
return(0); /* no sample taken */ |
110 |
> |
return(NULL); /* no sample taken */ |
111 |
|
} |
112 |
|
if (ambacc > FTINY) { |
113 |
|
multcolor(ar.rcoef, hp->acoef); |
115 |
|
} |
116 |
|
/* generate hemispherical sample */ |
117 |
|
SDsquare2disk(spt, (i+.1+.8*frandom())/hp->ns, |
118 |
< |
(j+.1+.8*frandom())/hp->ns); |
118 |
> |
(j+.1+.8*frandom())/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] + |
130 |
|
if (ar.rt > 20.0*maxarad) /* limit vertex distance */ |
131 |
|
ar.rt = 20.0*maxarad; |
132 |
|
VSUM(ap->p, ar.rorg, ar.rdir, ar.rt); |
133 |
< |
return(1); |
133 |
> |
return(ap); |
134 |
|
} |
135 |
|
|
136 |
|
|
159 |
|
} |
160 |
|
|
161 |
|
|
162 |
< |
/* Compose matrix from two vectors */ |
162 |
> |
/* Compose 3x3 matrix from two vectors */ |
163 |
|
static void |
164 |
|
compose_matrix(FVECT mat[3], FVECT va, FVECT vb) |
165 |
|
{ |
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*ftp->I2*d3 ); |
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 */ |
194 |
|
VCROSS(v1, nrm, ftp->e_i); |
195 |
|
for (j = 3; j--; ) |
196 |
< |
v2[i] = ftp->I2*ftp->r_i[j] + ftp->J2*ftp->e_i[j]; |
196 |
> |
v2[j] = ftp->I2*ftp->r_i[j] + ftp->J2*ftp->e_i[j]; |
197 |
|
compose_matrix(m1, v1, v2); |
198 |
|
compose_matrix(m2, ftp->r_i, ftp->r_i); |
199 |
|
compose_matrix(m3, ftp->e_i, ftp->e_i); |
321 |
|
hess2[1][0] = DOT(uv[1], a); |
322 |
|
hess2[1][1] = DOT(uv[1], b); |
323 |
|
/* compute eigenvalues */ |
324 |
< |
if (quadratic(evalue, 1.0, -hess2[0][0]-hess2[1][1], |
324 |
> |
if ( quadratic(evalue, 1.0, -hess2[0][0]-hess2[1][1], |
325 |
|
hess2[0][0]*hess2[1][1]-hess2[0][1]*hess2[1][0]) != 2 || |
326 |
< |
(evalue[0] = fabs(evalue[0])) <= FTINY*FTINY*FTINY || |
327 |
< |
(evalue[1] = fabs(evalue[1])) <= FTINY*FTINY*FTINY) |
326 |
> |
(evalue[0] = fabs(evalue[0])) <= FTINY*FTINY || |
327 |
> |
(evalue[1] = fabs(evalue[1])) <= FTINY*FTINY ) |
328 |
|
error(INTERNAL, "bad eigenvalue calculation"); |
329 |
|
|
330 |
|
if (evalue[0] > evalue[1]) { |
354 |
|
ambHessian( /* anisotropic radii & pos. gradient */ |
355 |
|
AMBHEMI *hp, |
356 |
|
FVECT uv[2], /* returned */ |
357 |
< |
float ra[2], /* returned */ |
358 |
< |
float pg[2] /* returned */ |
357 |
> |
float ra[2], /* returned (optional) */ |
358 |
> |
float pg[2] /* returned (optional) */ |
359 |
|
) |
360 |
|
{ |
361 |
|
static char memerrmsg[] = "out of memory in ambHessian()"; |
370 |
|
VCOPY(uv[1], hp->uy); |
371 |
|
/* clock-wise vertex traversal from sample POV */ |
372 |
|
if (ra != NULL) { /* initialize Hessian row buffer */ |
373 |
< |
hessrow = (FVECT (*)[3])malloc(sizeof(FVECT)*3*hp->ns); |
373 |
> |
hessrow = (FVECT (*)[3])malloc(sizeof(FVECT)*3*(hp->ns-1)); |
374 |
|
if (hessrow == NULL) |
375 |
|
error(SYSTEM, memerrmsg); |
376 |
|
memset(hessian, 0, sizeof(hessian)); |
377 |
|
} else if (pg == NULL) /* bogus call? */ |
378 |
|
return; |
379 |
|
if (pg != NULL) { /* initialize form factor row buffer */ |
380 |
< |
gradrow = (FVECT *)malloc(sizeof(FVECT)*hp->ns); |
380 |
> |
gradrow = (FVECT *)malloc(sizeof(FVECT)*(hp->ns-1)); |
381 |
|
if (gradrow == NULL) |
382 |
|
error(SYSTEM, memerrmsg); |
383 |
|
memset(gradient, 0, sizeof(gradient)); |
467 |
|
{ |
468 |
|
struct s_ambsamp *ap; |
469 |
|
int n; |
470 |
+ |
FVECT vd; |
471 |
+ |
double gfact; |
472 |
|
|
473 |
|
dg[0] = dg[1] = 0; |
474 |
|
for (ap = hp->sa, n = hp->ns*hp->ns; n--; ap++) { |
471 |
– |
FVECT vd; |
472 |
– |
double gfact; |
475 |
|
/* use vector for azimuth + 90deg */ |
476 |
|
VSUB(vd, ap->p, hp->rp->rop); |
477 |
|
/* brightness with tangent factor */ |
478 |
|
gfact = ap->v[CIEY] / DOT(hp->rp->ron, vd); |
479 |
|
/* sine = proj_radius/vd_length */ |
480 |
< |
dg[0] -= DOT(uv[1], vd) * gfact ; |
480 |
> |
dg[0] -= DOT(uv[1], vd) * gfact; |
481 |
|
dg[1] += DOT(uv[0], vd) * gfact; |
482 |
|
} |
483 |
|
} |
494 |
|
float dg[2] /* returned (optional) */ |
495 |
|
) |
496 |
|
{ |
497 |
+ |
AMBHEMI *hp = inithemi(rcol, r, wt); |
498 |
|
int cnt = 0; |
499 |
|
FVECT my_uv[2]; |
497 |
– |
AMBHEMI *hp; |
500 |
|
double d, acol[3]; |
501 |
|
struct s_ambsamp *ap; |
502 |
|
int i, j; |
503 |
< |
/* initialize */ |
504 |
< |
if ((hp = inithemi(rcol, r, wt)) == NULL) |
503 |
> |
/* check/initialize */ |
504 |
> |
if (hp == NULL) |
505 |
|
return(0); |
506 |
|
if (uv != NULL) |
507 |
|
memset(uv, 0, sizeof(FVECT)*2); |
515 |
|
acol[0] = acol[1] = acol[2] = 0.0; |
516 |
|
for (i = hp->ns; i--; ) |
517 |
|
for (j = hp->ns; j--; ) |
518 |
< |
if (ambsample(hp, i, j)) { |
517 |
< |
ap = &ambsamp(hp,i,j); |
518 |
> |
if ((ap = ambsample(hp, i, j)) != NULL) { |
519 |
|
addcolor(acol, ap->v); |
520 |
|
++cnt; |
521 |
|
} |
544 |
|
ambHessian(hp, uv, ra, pg); |
545 |
|
if (dg != NULL) /* compute direction gradient */ |
546 |
|
ambdirgrad(hp, uv, dg); |
547 |
< |
if (ra != NULL) { /* adjust/clamp radii */ |
547 |
> |
if (ra != NULL) { /* scale/clamp radii */ |
548 |
|
d = sqrt(sqrt((4.0/PI)*bright(rcol)/wt)); |
549 |
< |
if ((ra[0] *= d) > maxarad) |
549 |
< |
ra[0] = maxarad; |
549 |
> |
ra[0] *= d; |
550 |
|
if ((ra[1] *= d) > 2.0*ra[0]) |
551 |
|
ra[1] = 2.0*ra[0]; |
552 |
+ |
if (ra[1] > maxarad) { |
553 |
+ |
ra[1] = maxarad; |
554 |
+ |
if (ra[0] > maxarad) |
555 |
+ |
ra[0] = maxarad; |
556 |
+ |
} |
557 |
|
} |
558 |
|
free(hp); /* clean up and return */ |
559 |
|
return(1); |