--- ray/src/rt/aniso.c 1992/05/16 08:37:03 2.17 +++ ray/src/rt/aniso.c 1994/12/21 09:51:44 2.28 @@ -19,9 +19,13 @@ static char SCCSid[] = "$SunId$ LBL"; extern double specthresh; /* specular sampling threshold */ extern double specjitter; /* specular sampling jitter */ +extern int backvis; /* back faces visible? */ + +static agaussamp(), getacoords(); + /* - * This anisotropic reflection model uses a variant on the - * exponential Gaussian used in normal.c. + * This routine implements the anisotropic Gaussian + * model described by Ward in Siggraph `92 article. * We orient the surface towards the incoming ray, so a single * surface can be used to represent an infinitely thin object. * @@ -34,8 +38,6 @@ extern double specjitter; /* specular sampling jitte * 8 red grn blu rspec u-rough v-rough trans tspec */ -#define BSPEC(m) (6.0) /* specularity parameter b */ - /* specularity flags */ #define SP_REFL 01 /* has reflected specular component */ #define SP_TRAN 02 /* has transmitted specular */ @@ -53,8 +55,8 @@ typedef struct { FVECT vrefl; /* vector in reflected direction */ FVECT prdir; /* vector in transmitted direction */ FVECT u, v; /* u and v vectors orienting anisotropy */ - double u_alpha2; /* u roughness squared */ - double v_alpha2; /* v roughness squared */ + double u_alpha; /* u roughness */ + double v_alpha; /* v roughness */ double rdiff, rspec; /* reflected specular, diffuse */ double trans; /* transmissivity */ double tdiff, tspec; /* transmitted specular, diffuse */ @@ -103,21 +105,21 @@ double omega; /* light source size */ au2 = av2 = omega/(4.0*PI); else au2 = av2 = 0.0; - au2 += np->u_alpha2; - av2 += np->v_alpha2; + au2 += np->u_alpha*np->u_alpha; + av2 += np->v_alpha*np->v_alpha; /* half vector */ h[0] = ldir[0] - np->rp->rdir[0]; h[1] = ldir[1] - np->rp->rdir[1]; h[2] = ldir[2] - np->rp->rdir[2]; - normalize(h); /* ellipse */ dtmp1 = DOT(np->u, h); dtmp1 *= dtmp1 / au2; dtmp2 = DOT(np->v, h); dtmp2 *= dtmp2 / av2; /* gaussian */ - dtmp = (dtmp1 + dtmp2) / (1.0 + DOT(np->pnorm, h)); - dtmp = exp(-2.0*dtmp) * 1.0/(4.0*PI) + dtmp = DOT(np->pnorm, h); + dtmp = (dtmp1 + dtmp2) / (dtmp*dtmp); + dtmp = exp(-dtmp) * (0.25/PI) * sqrt(ldot/(np->pdot*au2*av2)); /* worth using? */ if (dtmp > FTINY) { @@ -143,25 +145,27 @@ double omega; /* light source size */ */ /* roughness + source */ au2 = av2 = omega / PI; - au2 += np->u_alpha2; - av2 += np->v_alpha2; + au2 += np->u_alpha*np->u_alpha; + av2 += np->v_alpha*np->v_alpha; /* "half vector" */ h[0] = ldir[0] - np->prdir[0]; h[1] = ldir[1] - np->prdir[1]; h[2] = ldir[2] - np->prdir[2]; - dtmp = DOT(h,np->pnorm); - dtmp = DOT(h,h) - dtmp*dtmp; + dtmp = DOT(h,h); if (dtmp > FTINY*FTINY) { - dtmp1 = DOT(h,np->u); - dtmp1 = dtmp1*dtmp1 / (au2*dtmp); - dtmp2 = DOT(h,np->v); - dtmp2 = dtmp2*dtmp2 / (av2*dtmp); - dtmp = 2. - 2.*DOT(ldir,np->prdir); - dtmp *= dtmp1 + dtmp2; + dtmp1 = DOT(h,np->pnorm); + dtmp = 1.0 - dtmp1*dtmp1/dtmp; + if (dtmp > FTINY*FTINY) { + dtmp1 = DOT(h,np->u); + dtmp1 *= dtmp1 / au2; + dtmp2 = DOT(h,np->v); + dtmp2 *= dtmp2 / av2; + dtmp = (dtmp1 + dtmp2) / dtmp; + } } else dtmp = 0.0; /* gaussian */ - dtmp = exp(-dtmp) * 1.0/(4.0*PI) + dtmp = exp(-dtmp) * (1.0/PI) * sqrt(-ldot/(np->pdot*au2*av2)); /* worth using? */ if (dtmp > FTINY) { @@ -179,12 +183,11 @@ register OBJREC *m; register RAY *r; { ANISODAT nd; - double dtmp; COLOR ctmp; register int i; /* easy shadow test */ if (r->crtype & SHADOW) - return; + return(1); if (m->oargs.nfargs != (m->otype == MAT_TRANS2 ? 8 : 6)) objerror(m, USER, "bad number of real arguments"); @@ -196,15 +199,18 @@ register RAY *r; m->oargs.farg[2]); /* get roughness */ nd.specfl = 0; - nd.u_alpha2 = m->oargs.farg[4]; - nd.u_alpha2 *= nd.u_alpha2; - nd.v_alpha2 = m->oargs.farg[5]; - nd.v_alpha2 *= nd.v_alpha2; - if (nd.u_alpha2 < FTINY*FTINY || nd.v_alpha2 <= FTINY*FTINY) + nd.u_alpha = m->oargs.farg[4]; + nd.v_alpha = m->oargs.farg[5]; + if (nd.u_alpha < FTINY || nd.v_alpha <= FTINY) objerror(m, USER, "roughness too small"); - /* reorient if necessary */ - if (r->rod < 0.0) - flipsurface(r); + /* check for back side */ + if (r->rod < 0.0) { + if (!backvis && m->otype != MAT_TRANS2) { + raytrans(r); + return(1); + } + flipsurface(r); /* reorient if backvis */ + } /* get modifiers */ raytexture(r, m->omod); nd.pdot = raynormal(nd.pnorm, r); /* perturb normal */ @@ -220,15 +226,8 @@ register RAY *r; else setcolor(nd.scolor, 1.0, 1.0, 1.0); scalecolor(nd.scolor, nd.rspec); - /* improved model */ - dtmp = exp(-BSPEC(m)*nd.pdot); - for (i = 0; i < 3; i++) - colval(nd.scolor,i) += (1.0-colval(nd.scolor,i))*dtmp; - nd.rspec += (1.0-nd.rspec)*dtmp; /* check threshold */ - if (specthresh > FTINY && - (specthresh >= 1.-FTINY || - specthresh + .05 - .1*frandom() > nd.rspec)) + if (specthresh >= nd.rspec-FTINY) nd.specfl |= SP_RBLT; /* compute refl. direction */ for (i = 0; i < 3; i++) @@ -245,9 +244,7 @@ register RAY *r; if (nd.tspec > FTINY) { nd.specfl |= SP_TRAN; /* check threshold */ - if (specthresh > FTINY && - (specthresh >= 1.-FTINY || - specthresh + .05 - .1*frandom() > nd.tspec)) + if (specthresh >= nd.tspec-FTINY) nd.specfl |= SP_TBLT; if (DOT(r->pert,r->pert) <= FTINY*FTINY) { VCOPY(nd.prdir, r->rdir); @@ -297,6 +294,8 @@ register RAY *r; } /* add direct component */ direct(r, diraniso, &nd); + + return(1); } @@ -347,9 +346,9 @@ register ANISODAT *np; d = urand(ilhash(dimlist,ndims)+samplendx); multisamp(rv, 2, d); d = 2.0*PI * rv[0]; - cosp = cos(d); - sinp = sin(d); - d = sqrt(np->u_alpha2*cosp*cosp + np->v_alpha2*sinp*sinp); + cosp = cos(d) * np->u_alpha; + sinp = sin(d) * np->v_alpha; + d = sqrt(cosp*cosp + sinp*sinp); cosp /= d; sinp /= d; rv[1] = 1.0 - specjitter*rv[1]; @@ -357,8 +356,8 @@ register ANISODAT *np; d = 1.0; else d = sqrt(-log(rv[1]) / - (cosp*cosp/np->u_alpha2 + - sinp*sinp/np->v_alpha2)); + (cosp*cosp/(np->u_alpha*np->u_alpha) + + sinp*sinp/(np->v_alpha*np->v_alpha))); for (i = 0; i < 3; i++) h[i] = np->pnorm[i] + d*(cosp*np->u[i] + sinp*np->v[i]); @@ -379,15 +378,18 @@ register ANISODAT *np; d = urand(ilhash(dimlist,ndims)+1823+samplendx); multisamp(rv, 2, d); d = 2.0*PI * rv[0]; - cosp = cos(d); - sinp = sin(d); + cosp = cos(d) * np->u_alpha; + sinp = sin(d) * np->v_alpha; + d = sqrt(cosp*cosp + sinp*sinp); + cosp /= d; + sinp /= d; rv[1] = 1.0 - specjitter*rv[1]; if (rv[1] <= FTINY) d = 1.0; else d = sqrt(-log(rv[1]) / - (cosp*cosp*4./np->u_alpha2 + - sinp*sinp*4./np->v_alpha2)); + (cosp*cosp/(np->u_alpha*np->u_alpha) + + sinp*sinp/(np->v_alpha*np->u_alpha))); for (i = 0; i < 3; i++) sr.rdir[i] = np->prdir[i] + d*(cosp*np->u[i] + sinp*np->v[i]);