--- ray/src/rt/ambcomp.c 2014/04/23 06:04:17 2.29 +++ ray/src/rt/ambcomp.c 2014/04/30 18:27:14 2.40 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: ambcomp.c,v 2.29 2014/04/23 06:04:17 greg Exp $"; +static const char RCSid[] = "$Id: ambcomp.c,v 2.40 2014/04/30 18:27:14 greg Exp $"; #endif /* * Routines to compute "ambient" values using Monte Carlo @@ -28,15 +28,15 @@ 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 */ #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 */ @@ -70,14 +70,14 @@ inithemi( /* initialize sampling hemisphere */ d = 1.0/(n*n); scalecolor(hp->acoef, d); /* make tangent plane axes */ - hp->uy[0] = 0.1 - 0.2*frandom(); - hp->uy[1] = 0.1 - 0.2*frandom(); - hp->uy[2] = 0.1 - 0.2*frandom(); - for (i = 0; i < 3; i++) - if (r->ron[i] < 0.6 && r->ron[i] > -0.6) + 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); @@ -103,11 +103,8 @@ ambsample( /* sample an ambient direction */ setcolor(ar.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(NULL); /* no sample taken */ - } + if (rayorigin(&ar, AMBIENT, hp->rp, ar.rcoef) < 0) + goto badsample; if (ambacc > FTINY) { multcolor(ar.rcoef, hp->acoef); scalecolor(ar.rcoef, 1./AVGREFL); @@ -124,38 +121,45 @@ ambsample( /* sample an ambient direction */ 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 */ - VSUM(ap->p, ar.rorg, ar.rdir, 20.0*maxarad); - else - VCOPY(ap->p, ar.rop); return(ap); +badsample: + setcolor(ap->v, 0., 0., 0.); + VCOPY(ap->p, hp->rp->rop); + return(NULL); } /* 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(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->I1 = acos( DOT(ftp->r_i, ftp->r_i1) / sqrt(dot_r*dot_r1) ) * - sqrt( ftp->nf ); - 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.5/dot_e*( 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]; } @@ -176,7 +180,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; @@ -186,20 +190,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/d3*ftp->I2 ); + 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[j] = 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 */ @@ -207,7 +208,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; } } @@ -243,16 +244,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] ); } @@ -320,11 +319,13 @@ 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 || - (evalue[1] = fabs(evalue[1])) <= 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]) { @@ -415,7 +416,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); @@ -454,9 +455,9 @@ ambHessian( /* anisotropic radii & pos. gradient */ if (ra != NULL) /* extract eigenvectors & radii */ eigenvectors(uv, ra, hessian); - if (pg != NULL) { /* tangential position gradient/PI */ - pg[0] = DOT(gradient, uv[0]) / PI; - pg[1] = DOT(gradient, uv[1]) / PI; + if (pg != NULL) { /* tangential position gradient */ + pg[0] = DOT(gradient, uv[0]); + pg[1] = DOT(gradient, uv[1]); } } @@ -477,9 +478,9 @@ ambdirgrad(AMBHEMI *hp, FVECT uv[2], float dg[2]) VSUB(vd, ap->p, hp->rp->rop); /* brightness over cosine factor */ gfact = colval(ap->v,CIEY) / DOT(hp->rp->ron, vd); - /* -sine = -proj_radius/vd_length */ - dgsum[0] += DOT(uv[1], vd) * gfact; - dgsum[1] -= DOT(uv[0], vd) * gfact; + /* sine = proj_radius/vd_length */ + 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); @@ -500,7 +501,7 @@ doambient( /* compute ambient component */ AMBHEMI *hp = inithemi(rcol, r, wt); int cnt = 0; FVECT my_uv[2]; - double d, acol[3]; + double d, K, acol[3]; struct s_ambsamp *ap; int i, j; /* check/initialize */ @@ -533,14 +534,18 @@ doambient( /* compute ambient component */ free(hp); return(-1); /* no radius or gradient calc. */ } - multcolor(acol, hp->acoef); /* normalize Y values */ - if ((d = bright(acol)) > FTINY) - d = 1.0/d; - else + 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 + 0.0314; + colval(ap->v,CIEY) = bright(ap->v)*d + K; if (uv == NULL) /* make sure we have axis pointers */ uv = my_uv; @@ -551,6 +556,14 @@ doambient( /* compute ambient component */ ambdirgrad(hp, uv, dg); 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) @@ -563,6 +576,14 @@ doambient( /* compute ambient component */ 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 */