35 |
|
#define ambsamp(h,i,j) (h)->sa[(i)*(h)->ns + (j)] |
36 |
|
|
37 |
|
typedef struct { |
38 |
< |
FVECT r_i, r_i1, e_i, rI2_eJ2; |
39 |
< |
double nf, I1, I2; |
38 |
> |
FVECT r_i, r_i1, e_i, rcp, rI2_eJ2; |
39 |
> |
double I1, I2; |
40 |
|
} FFTRI; /* vectors and coefficients for Hessian calculation */ |
41 |
|
|
42 |
|
|
70 |
|
d = 1.0/(n*n); |
71 |
|
scalecolor(hp->acoef, d); |
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) |
73 |
> |
hp->uy[0] = 0.5 - frandom(); |
74 |
> |
hp->uy[1] = 0.5 - frandom(); |
75 |
> |
hp->uy[2] = 0.5 - frandom(); |
76 |
> |
for (i = 3; i--; ) |
77 |
> |
if ((-0.6 < r->ron[i]) & (r->ron[i] < 0.6)) |
78 |
|
break; |
79 |
< |
if (i >= 3) |
80 |
< |
error(CONSISTENCY, "bad ray direction in inithemi()"); |
79 |
> |
if (i < 0) |
80 |
> |
error(CONSISTENCY, "bad ray direction in inithemi"); |
81 |
|
hp->uy[i] = 1.0; |
82 |
|
VCROSS(hp->ux, hp->uy, r->ron); |
83 |
|
normalize(hp->ux); |
141 |
|
static void |
142 |
|
comp_fftri(FFTRI *ftp, FVECT ap0, FVECT ap1, FVECT rop) |
143 |
|
{ |
144 |
< |
FVECT vcp; |
145 |
< |
double dot_e, dot_er, rdot_r, rdot_r1, J2; |
144 |
> |
double rdot_cp, dot_e, dot_er, rdot_r, rdot_r1, J2; |
145 |
|
int i; |
146 |
|
|
147 |
|
VSUB(ftp->r_i, ap0, rop); |
148 |
|
VSUB(ftp->r_i1, ap1, rop); |
149 |
|
VSUB(ftp->e_i, ap1, ap0); |
150 |
< |
VCROSS(vcp, ftp->e_i, ftp->r_i); |
151 |
< |
ftp->nf = 1.0/DOT(vcp,vcp); |
150 |
> |
VCROSS(ftp->rcp, ftp->r_i, ftp->r_i1); |
151 |
> |
rdot_cp = 1.0/DOT(ftp->rcp,ftp->rcp); |
152 |
|
dot_e = DOT(ftp->e_i,ftp->e_i); |
153 |
|
dot_er = DOT(ftp->e_i, ftp->r_i); |
154 |
|
rdot_r = 1.0/DOT(ftp->r_i,ftp->r_i); |
155 |
|
rdot_r1 = 1.0/DOT(ftp->r_i1,ftp->r_i1); |
156 |
|
ftp->I1 = acos( DOT(ftp->r_i, ftp->r_i1) * sqrt(rdot_r*rdot_r1) ) * |
157 |
< |
sqrt( ftp->nf ); |
157 |
> |
sqrt( rdot_cp ); |
158 |
|
ftp->I2 = ( DOT(ftp->e_i, ftp->r_i1)*rdot_r1 - dot_er*rdot_r + |
159 |
< |
dot_e*ftp->I1 )*0.5*ftp->nf; |
159 |
> |
dot_e*ftp->I1 )*0.5*rdot_cp; |
160 |
|
J2 = ( 0.5*(rdot_r - rdot_r1) - dot_er*ftp->I2 ) / dot_e; |
161 |
|
for (i = 3; i--; ) |
162 |
|
ftp->rI2_eJ2[i] = ftp->I2*ftp->r_i[i] + J2*ftp->e_i[i]; |
180 |
|
static void |
181 |
|
comp_hessian(FVECT hess[3], FFTRI *ftp, FVECT nrm) |
182 |
|
{ |
183 |
< |
FVECT vcp; |
183 |
> |
FVECT ncp; |
184 |
|
FVECT m1[3], m2[3], m3[3], m4[3]; |
185 |
|
double d1, d2, d3, d4; |
186 |
|
double I3, J3, K3; |
190 |
|
d2 = 1.0/DOT(ftp->r_i1,ftp->r_i1); |
191 |
|
d3 = 1.0/DOT(ftp->e_i,ftp->e_i); |
192 |
|
d4 = DOT(ftp->e_i, ftp->r_i); |
193 |
< |
I3 = 0.25*ftp->nf*( DOT(ftp->e_i, ftp->r_i1)*d2*d2 - d4*d1*d1 + |
194 |
< |
3.0/d3*ftp->I2 ); |
193 |
> |
I3 = ( DOT(ftp->e_i, ftp->r_i1)*d2*d2 - d4*d1*d1 + 3.0/d3*ftp->I2 ) |
194 |
> |
/ ( 4.0*DOT(ftp->rcp,ftp->rcp) ); |
195 |
|
J3 = 0.25*d3*(d1*d1 - d2*d2) - d4*d3*I3; |
196 |
|
K3 = d3*(ftp->I2 - I3/d1 - 2.0*d4*J3); |
197 |
|
/* intermediate matrices */ |
198 |
< |
VCROSS(vcp, nrm, ftp->e_i); |
199 |
< |
compose_matrix(m1, vcp, ftp->rI2_eJ2); |
198 |
> |
VCROSS(ncp, nrm, ftp->e_i); |
199 |
> |
compose_matrix(m1, ncp, ftp->rI2_eJ2); |
200 |
|
compose_matrix(m2, ftp->r_i, ftp->r_i); |
201 |
|
compose_matrix(m3, ftp->e_i, ftp->e_i); |
202 |
|
compose_matrix(m4, ftp->r_i, ftp->e_i); |
203 |
< |
VCROSS(vcp, ftp->r_i, ftp->e_i); |
205 |
< |
d1 = DOT(nrm, vcp); |
203 |
> |
d1 = DOT(nrm, ftp->rcp); |
204 |
|
d2 = -d1*ftp->I2; |
205 |
|
d1 *= 2.0; |
206 |
|
for (i = 3; i--; ) /* final matrix sum */ |
244 |
|
static void |
245 |
|
comp_gradient(FVECT grad, FFTRI *ftp, FVECT nrm) |
246 |
|
{ |
247 |
< |
FVECT vcp; |
247 |
> |
FVECT ncp; |
248 |
|
double f1; |
249 |
|
int i; |
250 |
|
|
251 |
< |
VCROSS(vcp, ftp->r_i, ftp->r_i1); |
252 |
< |
f1 = 2.0*DOT(nrm, vcp); |
255 |
< |
VCROSS(vcp, nrm, ftp->e_i); |
251 |
> |
f1 = 2.0*DOT(nrm, ftp->rcp); |
252 |
> |
VCROSS(ncp, nrm, ftp->e_i); |
253 |
|
for (i = 3; i--; ) |
254 |
< |
grad[i] = (-0.5/PI)*( ftp->I1*vcp[i] + f1*ftp->rI2_eJ2[i] ); |
254 |
> |
grad[i] = (-0.5/PI)*( ftp->I1*ncp[i] + f1*ftp->rI2_eJ2[i] ); |
255 |
|
} |
256 |
|
|
257 |
|
|
319 |
|
hess2[0][1] = DOT(uv[0], b); |
320 |
|
hess2[1][0] = DOT(uv[1], a); |
321 |
|
hess2[1][1] = DOT(uv[1], b); |
322 |
< |
/* compute eigenvalues */ |
323 |
< |
if ( quadratic(evalue, 1.0, -hess2[0][0]-hess2[1][1], |
324 |
< |
hess2[0][0]*hess2[1][1]-hess2[0][1]*hess2[1][0]) != 2 || |
325 |
< |
(evalue[0] = fabs(evalue[0])) <= FTINY*FTINY || |
326 |
< |
(evalue[1] = fabs(evalue[1])) <= FTINY*FTINY ) |
322 |
> |
/* compute eigenvalue(s) */ |
323 |
> |
i = quadratic(evalue, 1.0, -hess2[0][0]-hess2[1][1], |
324 |
> |
hess2[0][0]*hess2[1][1]-hess2[0][1]*hess2[1][0]); |
325 |
> |
if (i == 1) /* double-root (circle) */ |
326 |
> |
evalue[1] = evalue[0]; |
327 |
> |
if (!i || ((evalue[0] = fabs(evalue[0])) <= FTINY*FTINY) | |
328 |
> |
((evalue[1] = fabs(evalue[1])) <= FTINY*FTINY) ) |
329 |
|
error(INTERNAL, "bad eigenvalue calculation"); |
330 |
|
|
331 |
|
if (evalue[0] > evalue[1]) { |
416 |
|
rev_hessian(hesscol); |
417 |
|
add2hessian(hessian, hessrow[j], hessdia, hesscol, backg); |
418 |
|
} |
419 |
< |
if (gradient != NULL) { |
419 |
> |
if (gradrow != NULL) { |
420 |
|
comp_gradient(graddia, &fftr, hp->rp->ron); |
421 |
|
rev_gradient(gradcol); |
422 |
|
add2gradient(gradient, gradrow[j], graddia, gradcol, backg); |
478 |
|
VSUB(vd, ap->p, hp->rp->rop); |
479 |
|
/* brightness over cosine factor */ |
480 |
|
gfact = colval(ap->v,CIEY) / DOT(hp->rp->ron, vd); |
481 |
< |
/* -sine = -proj_radius/vd_length */ |
482 |
< |
dgsum[0] += DOT(uv[1], vd) * gfact; |
483 |
< |
dgsum[1] -= DOT(uv[0], vd) * gfact; |
481 |
> |
/* sine = proj_radius/vd_length */ |
482 |
> |
dgsum[0] -= DOT(uv[1], vd) * gfact; |
483 |
> |
dgsum[1] += DOT(uv[0], vd) * gfact; |
484 |
|
} |
485 |
|
dg[0] = dgsum[0] / (hp->ns*hp->ns); |
486 |
|
dg[1] = dgsum[1] / (hp->ns*hp->ns); |
501 |
|
AMBHEMI *hp = inithemi(rcol, r, wt); |
502 |
|
int cnt = 0; |
503 |
|
FVECT my_uv[2]; |
504 |
< |
double d, acol[3]; |
504 |
> |
double d, K, acol[3]; |
505 |
|
struct s_ambsamp *ap; |
506 |
|
int i, j; |
507 |
|
/* check/initialize */ |
534 |
|
free(hp); |
535 |
|
return(-1); /* no radius or gradient calc. */ |
536 |
|
} |
537 |
< |
if (bright(acol) > FTINY) /* normalize Y values */ |
538 |
< |
d = cnt/bright(acol); |
539 |
< |
else |
537 |
> |
if (bright(acol) > FTINY) { /* normalize Y values */ |
538 |
> |
d = 0.99*cnt/bright(acol); |
539 |
> |
K = 0.01; |
540 |
> |
} else { /* geometric Hessian fall-back */ |
541 |
|
d = 0.0; |
542 |
+ |
K = 1.0; |
543 |
+ |
pg = NULL; |
544 |
+ |
dg = NULL; |
545 |
+ |
} |
546 |
|
ap = hp->sa; /* relative Y channel from here on... */ |
547 |
|
for (i = hp->ns*hp->ns; i--; ap++) |
548 |
< |
colval(ap->v,CIEY) = bright(ap->v)*d + 0.01; |
548 |
> |
colval(ap->v,CIEY) = bright(ap->v)*d + K; |
549 |
|
|
550 |
|
if (uv == NULL) /* make sure we have axis pointers */ |
551 |
|
uv = my_uv; |
556 |
|
ambdirgrad(hp, uv, dg); |
557 |
|
|
558 |
|
if (ra != NULL) { /* scale/clamp radii */ |
559 |
+ |
if (pg != NULL) { |
560 |
+ |
if (ra[0]*(d = fabs(pg[0])) > 1.0) |
561 |
+ |
ra[0] = 1.0/d; |
562 |
+ |
if (ra[1]*(d = fabs(pg[1])) > 1.0) |
563 |
+ |
ra[1] = 1.0/d; |
564 |
+ |
if (ra[0] > ra[1]) |
565 |
+ |
ra[0] = ra[1]; |
566 |
+ |
} |
567 |
|
if (ra[0] < minarad) { |
568 |
|
ra[0] = minarad; |
569 |
|
if (ra[1] < minarad) |
570 |
|
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 |
– |
} |
571 |
|
} |
572 |
|
ra[0] *= d = 1.0/sqrt(sqrt(wt)); |
573 |
|
if ((ra[1] *= d) > 2.0*ra[0]) |
576 |
|
ra[1] = maxarad; |
577 |
|
if (ra[0] > maxarad) |
578 |
|
ra[0] = maxarad; |
579 |
+ |
} |
580 |
+ |
if (pg != NULL) { /* cap gradient if necessary */ |
581 |
+ |
d = pg[0]*pg[0]*ra[0]*ra[0] + pg[1]*pg[1]*ra[1]*ra[1]; |
582 |
+ |
if (d > 1.0) { |
583 |
+ |
d = 1.0/sqrt(d); |
584 |
+ |
pg[0] *= d; |
585 |
+ |
pg[1] *= d; |
586 |
+ |
} |
587 |
|
} |
588 |
|
} |
589 |
|
free(hp); /* clean up and return */ |