--- ray/src/rt/ambcomp.c 2014/04/19 02:39:44 2.27 +++ ray/src/rt/ambcomp.c 2014/04/30 23:38:58 2.41 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: ambcomp.c,v 2.27 2014/04/19 02:39:44 greg Exp $"; +static const char RCSid[] = "$Id: ambcomp.c,v 2.41 2014/04/30 23:38:58 greg Exp $"; #endif /* * Routines to compute "ambient" values using Monte Carlo @@ -28,15 +28,17 @@ typedef struct { COLOR acoef; /* division contribution coefficient */ struct s_ambsamp { COLOR v; /* hemisphere sample value */ - float p[3]; /* intersection point */ + FVECT p; /* intersection point */ } sa[1]; /* sample array (extends struct) */ } AMBHEMI; /* ambient sample hemisphere */ +typedef struct s_ambsamp AMBSAMP; + #define ambsamp(h,i,j) (h)->sa[(i)*(h)->ns + (j)] typedef struct { - FVECT r_i, r_i1, e_i; - double nf, I1, I2, J2; + FVECT r_i, r_i1, e_i, rcp, rI2_eJ2; + double I1, I2; } FFTRI; /* vectors and coefficients for Hessian calculation */ @@ -59,8 +61,7 @@ inithemi( /* initialize sampling hemisphere */ if (n < i) n = i; /* allocate sampling array */ - hp = (AMBHEMI *)malloc(sizeof(AMBHEMI) + - sizeof(struct s_ambsamp)*(n*n - 1)); + hp = (AMBHEMI *)malloc(sizeof(AMBHEMI) + sizeof(AMBSAMP)*(n*n - 1)); if (hp == NULL) return(NULL); hp->rp = r; @@ -69,13 +70,15 @@ inithemi( /* initialize sampling hemisphere */ copycolor(hp->acoef, ac); d = 1.0/(n*n); scalecolor(hp->acoef, d); - /* make tangent axes */ - hp->uy[0] = hp->uy[1] = hp->uy[2] = 0.0; - for (i = 0; i < 3; i++) - if (r->ron[i] < 0.6 && r->ron[i] > -0.6) + /* make tangent plane axes */ + hp->uy[0] = 0.5 - frandom(); + hp->uy[1] = 0.5 - frandom(); + hp->uy[2] = 0.5 - frandom(); + for (i = 3; i--; ) + if ((-0.6 < r->ron[i]) & (r->ron[i] < 0.6)) break; - if (i >= 3) - error(CONSISTENCY, "bad ray direction in inithemi()"); + if (i < 0) + error(CONSISTENCY, "bad ray direction in inithemi"); hp->uy[i] = 1.0; VCROSS(hp->ux, hp->uy, r->ron); normalize(hp->ux); @@ -85,79 +88,193 @@ inithemi( /* initialize sampling hemisphere */ } +/* Prepare ambient division sample */ static int -ambsample( /* sample an ambient direction */ - AMBHEMI *hp, - int i, - int j -) +prepambsamp(RAY *arp, AMBHEMI *hp, int i, int j, int n) { - struct s_ambsamp *ap = &ambsamp(hp,i,j); - RAY ar; - int hlist[3]; - double spt[2], zd; - int ii; + int hlist[3], ii; + double spt[2], zd; /* ambient coefficient for weight */ if (ambacc > FTINY) - setcolor(ar.rcoef, AVGREFL, AVGREFL, AVGREFL); + setcolor(arp->rcoef, AVGREFL, AVGREFL, AVGREFL); else - copycolor(ar.rcoef, hp->acoef); - if (rayorigin(&ar, AMBIENT, hp->rp, ar.rcoef) < 0) { - setcolor(ap->v, 0., 0., 0.); - VCOPY(ap->p, hp->rp->rop); - return(0); /* no sample taken */ - } + copycolor(arp->rcoef, hp->acoef); + if (rayorigin(arp, AMBIENT, hp->rp, arp->rcoef) < 0) + return(0); if (ambacc > FTINY) { - multcolor(ar.rcoef, hp->acoef); - scalecolor(ar.rcoef, 1./AVGREFL); + multcolor(arp->rcoef, hp->acoef); + scalecolor(arp->rcoef, 1./AVGREFL); } - /* generate hemispherical sample */ - SDsquare2disk(spt, (i+.1+.8*frandom())/hp->ns, - (j+.1+.8*frandom())/hp->ns); + hlist[0] = hp->rp->rno; + hlist[1] = i; + hlist[2] = j; + multisamp(spt, 2, urand(ilhash(hlist,3)+n)); + if (!n) { /* avoid border samples for n==0 */ + if ((spt[0] < 0.1) | (spt[0] > 0.9)) + spt[0] = 0.1 + 0.8*frandom(); + if ((spt[1] < 0.1) | (spt[1] > 0.9)) + spt[1] = 0.1 + 0.8*frandom(); + } + SDsquare2disk(spt, (i+spt[0])/hp->ns, (j+spt[1])/hp->ns); zd = sqrt(1. - spt[0]*spt[0] - spt[1]*spt[1]); for (ii = 3; ii--; ) - ar.rdir[ii] = spt[0]*hp->ux[ii] + + arp->rdir[ii] = spt[0]*hp->ux[ii] + spt[1]*hp->uy[ii] + zd*hp->rp->ron[ii]; - checknorm(ar.rdir); + checknorm(arp->rdir); + return(1); +} + + +static AMBSAMP * +ambsample( /* sample an ambient direction */ + AMBHEMI *hp, + int i, + int j +) +{ + AMBSAMP *ap = &ambsamp(hp,i,j); + RAY ar; + /* generate hemispherical sample */ + if (!prepambsamp(&ar, hp, i, j, 0)) + goto badsample; dimlist[ndims++] = i*hp->ns + j + 90171; rayvalue(&ar); /* evaluate ray */ ndims--; + /* limit vertex distance */ + if (ar.rt > 10.0*thescene.cusize) + ar.rt = 10.0*thescene.cusize; + else if (ar.rt <= FTINY) /* should never happen! */ + goto badsample; + VSUM(ap->p, ar.rorg, ar.rdir, ar.rt); multcolor(ar.rcol, ar.rcoef); /* apply coefficient */ copycolor(ap->v, ar.rcol); - if (ar.rt > 20.0*maxarad) /* limit vertex distance */ - ar.rt = 20.0*maxarad; - VSUM(ap->p, ar.rorg, ar.rdir, ar.rt); - return(1); + return(ap); +badsample: + setcolor(ap->v, 0., 0., 0.); + VCOPY(ap->p, hp->rp->rop); + return(NULL); } +/* Estimate errors based on ambient division differences */ +static float * +getambdiffs(AMBHEMI *hp) +{ + float *earr = calloc(hp->ns*hp->ns, sizeof(float)); + float *ep; + double b, d2; + int i, j; + + if (earr == NULL) /* out of memory? */ + return(NULL); + /* compute squared neighbor diffs */ + for (ep = earr, i = 0; i < hp->ns; i++) + for (j = 0; j < hp->ns; j++, ep++) { + b = bright(ambsamp(hp,i,j).v); + if (i) { /* from above */ + d2 = b - bright(ambsamp(hp,i-1,j).v); + d2 *= d2; + ep[0] += d2; + ep[-hp->ns] += d2; + } + if (j) { /* from behind */ + d2 = b - bright(ambsamp(hp,i,j-1).v); + d2 *= d2; + ep[0] += d2; + ep[-1] += d2; + } + } + /* correct for number of neighbors */ + earr[0] *= 2.f; + earr[hp->ns-1] *= 2.f; + earr[(hp->ns-1)*hp->ns] *= 2.f; + earr[(hp->ns-1)*hp->ns + hp->ns-1] *= 2.f; + for (i = 1; i < hp->ns-1; i++) { + earr[i*hp->ns] *= 4./3.; + earr[i*hp->ns + hp->ns-1] *= 4./3.; + } + for (j = 1; j < hp->ns-1; j++) { + earr[j] *= 4./3.; + earr[(hp->ns-1)*hp->ns + j] *= 4./3.; + } + return(earr); +} + + +/* Perform super-sampling on hemisphere */ +static void +ambsupersamp(double acol[3], AMBHEMI *hp, int cnt) +{ + float *earr = getambdiffs(hp); + double e2sum = 0; + AMBSAMP *ap; + RAY ar; + COLOR asum; + float *ep; + int i, j, n; + + if (earr == NULL) /* just skip calc. if no memory */ + return; + /* add up estimated variances */ + for (ep = earr + hp->ns*hp->ns; ep-- > earr; ) + e2sum += *ep; + ep = earr; /* perform super-sampling */ + for (ap = hp->sa, i = 0; i < hp->ns; i++) + for (j = 0; j < hp->ns; j++, ap++) { + int nss = *ep/e2sum*cnt + frandom(); + setcolor(asum, 0., 0., 0.); + for (n = 1; n <= nss; n++) { + if (!prepambsamp(&ar, hp, i, j, n)) { + nss = n-1; + break; + } + dimlist[ndims++] = i*hp->ns + j + 90171; + rayvalue(&ar); /* evaluate super-sample */ + ndims--; + multcolor(ar.rcol, ar.rcoef); + addcolor(asum, ar.rcol); + } + if (nss) { /* update returned ambient value */ + const double ssf = 1./(nss + 1); + for (n = 3; n--; ) + acol[n] += ssf*colval(asum,n) + + (ssf - 1.)*colval(ap->v,n); + } + e2sum -= *ep++; /* update remainders */ + cnt -= nss; + } + free(earr); +} + + /* Compute vectors and coefficients for Hessian/gradient calcs */ static void -comp_fftri(FFTRI *ftp, float ap0[3], float ap1[3], FVECT rop) +comp_fftri(FFTRI *ftp, FVECT ap0, FVECT ap1, FVECT rop) { - FVECT v1; - double dot_e, dot_er, dot_r, dot_r1; + double rdot_cp, dot_e, dot_er, rdot_r, rdot_r1, J2; + int i; VSUB(ftp->r_i, ap0, rop); VSUB(ftp->r_i1, ap1, rop); VSUB(ftp->e_i, ap1, ap0); - VCROSS(v1, ftp->e_i, ftp->r_i); - ftp->nf = 1.0/DOT(v1,v1); - VCROSS(v1, ftp->r_i, ftp->r_i1); - ftp->I1 = sqrt(DOT(v1,v1)*ftp->nf); + VCROSS(ftp->rcp, ftp->r_i, ftp->r_i1); + rdot_cp = 1.0/DOT(ftp->rcp,ftp->rcp); dot_e = DOT(ftp->e_i,ftp->e_i); dot_er = DOT(ftp->e_i, ftp->r_i); - dot_r = DOT(ftp->r_i,ftp->r_i); - dot_r1 = DOT(ftp->r_i1,ftp->r_i1); - ftp->I2 = ( DOT(ftp->e_i, ftp->r_i1)/dot_r1 - dot_er/dot_r + - dot_e*ftp->I1 )*0.5*ftp->nf; - ftp->J2 = 0.25*ftp->nf*( 1.0/dot_r - 1.0/dot_r1 ) - - dot_er/dot_e*ftp->I2; + rdot_r = 1.0/DOT(ftp->r_i,ftp->r_i); + rdot_r1 = 1.0/DOT(ftp->r_i1,ftp->r_i1); + ftp->I1 = acos( DOT(ftp->r_i, ftp->r_i1) * sqrt(rdot_r*rdot_r1) ) * + sqrt( rdot_cp ); + ftp->I2 = ( DOT(ftp->e_i, ftp->r_i1)*rdot_r1 - dot_er*rdot_r + + dot_e*ftp->I1 )*0.5*rdot_cp; + J2 = ( 0.5*(rdot_r - rdot_r1) - dot_er*ftp->I2 ) / dot_e; + for (i = 3; i--; ) + ftp->rI2_eJ2[i] = ftp->I2*ftp->r_i[i] + J2*ftp->e_i[i]; } -/* Compose matrix from two vectors */ +/* Compose 3x3 matrix from two vectors */ static void compose_matrix(FVECT mat[3], FVECT va, FVECT vb) { @@ -174,7 +291,7 @@ compose_matrix(FVECT mat[3], FVECT va, FVECT vb) static void comp_hessian(FVECT hess[3], FFTRI *ftp, FVECT nrm) { - FVECT v1, v2; + FVECT ncp; FVECT m1[3], m2[3], m3[3], m4[3]; double d1, d2, d3, d4; double I3, J3, K3; @@ -184,20 +301,17 @@ comp_hessian(FVECT hess[3], FFTRI *ftp, FVECT nrm) d2 = 1.0/DOT(ftp->r_i1,ftp->r_i1); d3 = 1.0/DOT(ftp->e_i,ftp->e_i); d4 = DOT(ftp->e_i, ftp->r_i); - I3 = 0.25*ftp->nf*( DOT(ftp->e_i, ftp->r_i1)*d2*d2 - d4*d1*d1 + - 3.0*ftp->I2*d3 ); + I3 = ( DOT(ftp->e_i, ftp->r_i1)*d2*d2 - d4*d1*d1 + 3.0/d3*ftp->I2 ) + / ( 4.0*DOT(ftp->rcp,ftp->rcp) ); J3 = 0.25*d3*(d1*d1 - d2*d2) - d4*d3*I3; K3 = d3*(ftp->I2 - I3/d1 - 2.0*d4*J3); /* intermediate matrices */ - VCROSS(v1, nrm, ftp->e_i); - for (j = 3; j--; ) - v2[i] = ftp->I2*ftp->r_i[j] + ftp->J2*ftp->e_i[j]; - compose_matrix(m1, v1, v2); + VCROSS(ncp, nrm, ftp->e_i); + compose_matrix(m1, ncp, ftp->rI2_eJ2); compose_matrix(m2, ftp->r_i, ftp->r_i); compose_matrix(m3, ftp->e_i, ftp->e_i); compose_matrix(m4, ftp->r_i, ftp->e_i); - VCROSS(v1, ftp->r_i, ftp->e_i); - d1 = DOT(nrm, v1); + d1 = DOT(nrm, ftp->rcp); d2 = -d1*ftp->I2; d1 *= 2.0; for (i = 3; i--; ) /* final matrix sum */ @@ -205,7 +319,7 @@ comp_hessian(FVECT hess[3], FFTRI *ftp, FVECT nrm) hess[i][j] = m1[i][j] + d1*( I3*m2[i][j] + K3*m3[i][j] + 2.0*J3*m4[i][j] ); hess[i][j] += d2*(i==j); - hess[i][j] *= -1.0/PI; + hess[i][j] *= 1.0/PI; } } @@ -241,16 +355,14 @@ add2hessian(FVECT hess[3], FVECT ehess1[3], static void comp_gradient(FVECT grad, FFTRI *ftp, FVECT nrm) { - FVECT vcp; + FVECT ncp; double f1; int i; - VCROSS(vcp, ftp->r_i, ftp->r_i1); - f1 = 2.0*DOT(nrm, vcp); - VCROSS(vcp, nrm, ftp->e_i); + f1 = 2.0*DOT(nrm, ftp->rcp); + VCROSS(ncp, nrm, ftp->e_i); for (i = 3; i--; ) - grad[i] = (0.5/PI)*( ftp->I1*vcp[i] + - f1*(ftp->I2*ftp->r_i[i] + ftp->J2*ftp->e_i[i]) ); + grad[i] = (-0.5/PI)*( ftp->I1*ncp[i] + f1*ftp->rI2_eJ2[i] ); } @@ -277,8 +389,7 @@ add2gradient(FVECT grad, FVECT egrad1, FVECT egrad2, F /* Return brightness of furthest ambient sample */ static COLORV -back_ambval(struct s_ambsamp *ap1, struct s_ambsamp *ap2, - struct s_ambsamp *ap3, FVECT orig) +back_ambval(AMBSAMP *ap1, AMBSAMP *ap2, AMBSAMP *ap3, FVECT orig) { COLORV vback; FVECT vec; @@ -286,17 +397,17 @@ back_ambval(struct s_ambsamp *ap1, struct s_ambsamp *a VSUB(vec, ap1->p, orig); d2best = DOT(vec,vec); - vback = ap1->v[CIEY]; + vback = colval(ap1->v,CIEY); VSUB(vec, ap2->p, orig); d2 = DOT(vec,vec); if (d2 > d2best) { d2best = d2; - vback = ap2->v[CIEY]; + vback = colval(ap2->v,CIEY); } VSUB(vec, ap3->p, orig); d2 = DOT(vec,vec); if (d2 > d2best) - return(ap3->v[CIEY]); + return(colval(ap3->v,CIEY)); return(vback); } @@ -318,20 +429,22 @@ eigenvectors(FVECT uv[2], float ra[2], FVECT hessian[3 hess2[0][1] = DOT(uv[0], b); hess2[1][0] = DOT(uv[1], a); hess2[1][1] = DOT(uv[1], b); - /* compute eigenvalues */ - if (quadratic(evalue, 1.0, -hess2[0][0]-hess2[1][1], - hess2[0][0]*hess2[1][1]-hess2[0][1]*hess2[1][0]) != 2 || - (evalue[0] = fabs(evalue[0])) <= FTINY*FTINY*FTINY || - (evalue[1] = fabs(evalue[1])) <= FTINY*FTINY*FTINY) + /* compute eigenvalue(s) */ + i = quadratic(evalue, 1.0, -hess2[0][0]-hess2[1][1], + hess2[0][0]*hess2[1][1]-hess2[0][1]*hess2[1][0]); + if (i == 1) /* double-root (circle) */ + evalue[1] = evalue[0]; + if (!i || ((evalue[0] = fabs(evalue[0])) <= FTINY*FTINY) | + ((evalue[1] = fabs(evalue[1])) <= FTINY*FTINY) ) error(INTERNAL, "bad eigenvalue calculation"); if (evalue[0] > evalue[1]) { - ra[0] = 1.0/sqrt(sqrt(evalue[0])); - ra[1] = 1.0/sqrt(sqrt(evalue[1])); + ra[0] = sqrt(sqrt(4.0/evalue[0])); + ra[1] = sqrt(sqrt(4.0/evalue[1])); slope1 = evalue[1]; } else { - ra[0] = 1.0/sqrt(sqrt(evalue[1])); - ra[1] = 1.0/sqrt(sqrt(evalue[0])); + ra[0] = sqrt(sqrt(4.0/evalue[1])); + ra[1] = sqrt(sqrt(4.0/evalue[0])); slope1 = evalue[0]; } /* compute unit eigenvectors */ @@ -352,8 +465,8 @@ static void ambHessian( /* anisotropic radii & pos. gradient */ AMBHEMI *hp, FVECT uv[2], /* returned */ - float ra[2], /* returned */ - float pg[2] /* returned */ + float ra[2], /* returned (optional) */ + float pg[2] /* returned (optional) */ ) { static char memerrmsg[] = "out of memory in ambHessian()"; @@ -368,14 +481,14 @@ ambHessian( /* anisotropic radii & pos. gradient */ VCOPY(uv[1], hp->uy); /* clock-wise vertex traversal from sample POV */ if (ra != NULL) { /* initialize Hessian row buffer */ - hessrow = (FVECT (*)[3])malloc(sizeof(FVECT)*3*hp->ns); + hessrow = (FVECT (*)[3])malloc(sizeof(FVECT)*3*(hp->ns-1)); if (hessrow == NULL) error(SYSTEM, memerrmsg); memset(hessian, 0, sizeof(hessian)); } else if (pg == NULL) /* bogus call? */ return; if (pg != NULL) { /* initialize form factor row buffer */ - gradrow = (FVECT *)malloc(sizeof(FVECT)*hp->ns); + gradrow = (FVECT *)malloc(sizeof(FVECT)*(hp->ns-1)); if (gradrow == NULL) error(SYSTEM, memerrmsg); memset(gradient, 0, sizeof(gradient)); @@ -413,7 +526,7 @@ ambHessian( /* anisotropic radii & pos. gradient */ rev_hessian(hesscol); add2hessian(hessian, hessrow[j], hessdia, hesscol, backg); } - if (gradient != NULL) { + if (gradrow != NULL) { comp_gradient(graddia, &fftr, hp->rp->ron); rev_gradient(gradcol); add2gradient(gradient, gradrow[j], graddia, gradcol, backg); @@ -452,7 +565,7 @@ ambHessian( /* anisotropic radii & pos. gradient */ if (ra != NULL) /* extract eigenvectors & radii */ eigenvectors(uv, ra, hessian); - if (pg != NULL) { /* project position gradient */ + if (pg != NULL) { /* tangential position gradient */ pg[0] = DOT(gradient, uv[0]); pg[1] = DOT(gradient, uv[1]); } @@ -463,21 +576,24 @@ ambHessian( /* anisotropic radii & pos. gradient */ static void ambdirgrad(AMBHEMI *hp, FVECT uv[2], float dg[2]) { - struct s_ambsamp *ap; - int n; + AMBSAMP *ap; + double dgsum[2]; + int n; + FVECT vd; + double gfact; - dg[0] = dg[1] = 0; + dgsum[0] = dgsum[1] = 0.0; /* sum values times -tan(theta) */ for (ap = hp->sa, n = hp->ns*hp->ns; n--; ap++) { - FVECT vd; - double gfact; /* use vector for azimuth + 90deg */ VSUB(vd, ap->p, hp->rp->rop); - /* brightness with tangent factor */ - gfact = ap->v[CIEY] / DOT(hp->rp->ron, vd); + /* brightness over cosine factor */ + gfact = colval(ap->v,CIEY) / DOT(hp->rp->ron, vd); /* sine = proj_radius/vd_length */ - dg[0] -= DOT(uv[1], vd) * gfact ; - dg[1] += DOT(uv[0], vd) * gfact; + dgsum[0] -= DOT(uv[1], vd) * gfact; + dgsum[1] += DOT(uv[0], vd) * gfact; } + dg[0] = dgsum[0] / (hp->ns*hp->ns); + dg[1] = dgsum[1] / (hp->ns*hp->ns); } @@ -492,14 +608,14 @@ doambient( /* compute ambient component */ float dg[2] /* returned (optional) */ ) { - int cnt = 0; - FVECT my_uv[2]; - AMBHEMI *hp; - double d, acol[3]; - struct s_ambsamp *ap; - int i, j; - /* initialize */ - if ((hp = inithemi(rcol, r, wt)) == NULL) + AMBHEMI *hp = inithemi(rcol, r, wt); + int cnt = 0; + FVECT my_uv[2]; + double d, K, acol[3]; + AMBSAMP *ap; + int i, j; + /* check/initialize */ + if (hp == NULL) return(0); if (uv != NULL) memset(uv, 0, sizeof(FVECT)*2); @@ -513,8 +629,7 @@ doambient( /* compute ambient component */ acol[0] = acol[1] = acol[2] = 0.0; for (i = hp->ns; i--; ) for (j = hp->ns; j--; ) - if (ambsample(hp, i, j)) { - ap = &ambsamp(hp,i,j); + if ((ap = ambsample(hp, i, j)) != NULL) { addcolor(acol, ap->v); ++cnt; } @@ -523,32 +638,70 @@ doambient( /* compute ambient component */ free(hp); return(0); /* no valid samples */ } - d = 1.0 / cnt; /* final indirect irradiance/PI */ - acol[0] *= d; acol[1] *= d; acol[2] *= d; - copycolor(rcol, acol); - if (cnt < hp->ns*hp->ns || /* incomplete sampling? */ - (ra == NULL) & (pg == NULL) & (dg == NULL)) { + if (cnt < hp->ns*hp->ns) { /* incomplete sampling? */ + copycolor(rcol, acol); free(hp); + return(-1); /* return value w/o Hessian */ + } + cnt = ambssamp*wt + 0.5; /* perform super-sampling? */ + if (cnt > 0) + ambsupersamp(acol, hp, cnt); + copycolor(rcol, acol); /* final indirect irradiance/PI */ + if ((ra == NULL) & (pg == NULL) & (dg == NULL)) { + free(hp); return(-1); /* no radius or gradient calc. */ } - d = 0.01 * bright(rcol); /* add in 1% before Hessian comp. */ - if (d < FTINY) d = FTINY; - ap = hp->sa; /* using Y channel from here on... */ + if (bright(acol) > FTINY) { /* normalize Y values */ + d = 0.99*cnt/bright(acol); + K = 0.01; + } else { /* geometric Hessian fall-back */ + d = 0.0; + K = 1.0; + pg = NULL; + dg = NULL; + } + ap = hp->sa; /* relative Y channel from here on... */ for (i = hp->ns*hp->ns; i--; ap++) - colval(ap->v,CIEY) = bright(ap->v) + d; + colval(ap->v,CIEY) = bright(ap->v)*d + K; if (uv == NULL) /* make sure we have axis pointers */ uv = my_uv; /* compute radii & pos. gradient */ ambHessian(hp, uv, ra, pg); + if (dg != NULL) /* compute direction gradient */ ambdirgrad(hp, uv, dg); - if (ra != NULL) { /* adjust/clamp radii */ - d = sqrt(sqrt((4.0/PI)*bright(rcol)/wt)); - if ((ra[0] *= d) > maxarad) - ra[0] = maxarad; + + if (ra != NULL) { /* scale/clamp radii */ + if (pg != NULL) { + if (ra[0]*(d = fabs(pg[0])) > 1.0) + ra[0] = 1.0/d; + if (ra[1]*(d = fabs(pg[1])) > 1.0) + ra[1] = 1.0/d; + if (ra[0] > ra[1]) + ra[0] = ra[1]; + } + if (ra[0] < minarad) { + ra[0] = minarad; + if (ra[1] < minarad) + ra[1] = minarad; + } + ra[0] *= d = 1.0/sqrt(sqrt(wt)); if ((ra[1] *= d) > 2.0*ra[0]) ra[1] = 2.0*ra[0]; + if (ra[1] > maxarad) { + ra[1] = maxarad; + if (ra[0] > maxarad) + ra[0] = maxarad; + } + if (pg != NULL) { /* cap gradient if necessary */ + d = pg[0]*pg[0]*ra[0]*ra[0] + pg[1]*pg[1]*ra[1]*ra[1]; + if (d > 1.0) { + d = 1.0/sqrt(d); + pg[0] *= d; + pg[1] *= d; + } + } } free(hp); /* clean up and return */ return(1);