Back to index

radiance  4R0+20100331
aniso.c
Go to the documentation of this file.
00001 #ifndef lint
00002 static const char RCSid[] = "$Id: aniso.c,v 2.43 2005/04/19 01:15:06 greg Exp $";
00003 #endif
00004 /*
00005  *  Shading functions for anisotropic materials.
00006  */
00007 
00008 #include "copyright.h"
00009 
00010 #include  "ray.h"
00011 #include  "ambient.h"
00012 #include  "otypes.h"
00013 #include  "rtotypes.h"
00014 #include  "source.h"
00015 #include  "func.h"
00016 #include  "random.h"
00017 
00018 #ifndef  MAXITER
00019 #define  MAXITER     10            /* maximum # specular ray attempts */
00020 #endif
00021 
00022 /*
00023  *     This routine implements the anisotropic Gaussian
00024  *  model described by Ward in Siggraph `92 article.
00025  *     We orient the surface towards the incoming ray, so a single
00026  *  surface can be used to represent an infinitely thin object.
00027  *
00028  *  Arguments for MAT_PLASTIC2 and MAT_METAL2 are:
00029  *  4+ ux     uy     uz     funcfile      [transform...]
00030  *  0
00031  *  6  red    grn    blu    specular-frac.       u-facet-slope v-facet-slope
00032  *
00033  *  Real arguments for MAT_TRANS2 are:
00034  *  8  red    grn    blu    rspec  u-rough       v-rough       trans  tspec
00035  */
00036 
00037                             /* specularity flags */
00038 #define  SP_REFL     01            /* has reflected specular component */
00039 #define  SP_TRAN     02            /* has transmitted specular */
00040 #define  SP_FLAT     04            /* reflecting surface is flat */
00041 #define  SP_RBLT     010           /* reflection below sample threshold */
00042 #define  SP_TBLT     020           /* transmission below threshold */
00043 #define  SP_BADU     040           /* bad u direction calculation */
00044 
00045 typedef struct {
00046        OBJREC  *mp;         /* material pointer */
00047        RAY  *rp;            /* ray pointer */
00048        short  specfl;              /* specularity flags, defined above */
00049        COLOR  mcolor;              /* color of this material */
00050        COLOR  scolor;              /* color of specular component */
00051        FVECT  vrefl;        /* vector in reflected direction */
00052        FVECT  prdir;        /* vector in transmitted direction */
00053        FVECT  u, v;         /* u and v vectors orienting anisotropy */
00054        double  u_alpha;     /* u roughness */
00055        double  v_alpha;     /* v roughness */
00056        double  rdiff, rspec;       /* reflected specular, diffuse */
00057        double  trans;              /* transmissivity */
00058        double  tdiff, tspec;       /* transmitted specular, diffuse */
00059        FVECT  pnorm;        /* perturbed surface normal */
00060        double  pdot;        /* perturbed dot product */
00061 }  ANISODAT;         /* anisotropic material data */
00062 
00063 static srcdirf_t diraniso;
00064 static void getacoords(RAY  *r, ANISODAT  *np);
00065 static void agaussamp(RAY  *r, ANISODAT  *np);
00066 
00067 
00068 static void
00069 diraniso(            /* compute source contribution */
00070        COLOR  cval,                /* returned coefficient */
00071        void  *nnp,          /* material data */
00072        FVECT  ldir,                /* light source direction */
00073        double  omega               /* light source size */
00074 )
00075 {
00076        register ANISODAT *np = nnp;
00077        double  ldot;
00078        double  dtmp, dtmp1, dtmp2;
00079        FVECT  h;
00080        double  au2, av2;
00081        COLOR  ctmp;
00082 
00083        setcolor(cval, 0.0, 0.0, 0.0);
00084 
00085        ldot = DOT(np->pnorm, ldir);
00086 
00087        if (ldot < 0.0 ? np->trans <= FTINY : np->trans >= 1.0-FTINY)
00088               return;              /* wrong side */
00089 
00090        if (ldot > FTINY && np->rdiff > FTINY) {
00091               /*
00092                *  Compute and add diffuse reflected component to returned
00093                *  color.  The diffuse reflected component will always be
00094                *  modified by the color of the material.
00095                */
00096               copycolor(ctmp, np->mcolor);
00097               dtmp = ldot * omega * np->rdiff * (1.0/PI);
00098               scalecolor(ctmp, dtmp);
00099               addcolor(cval, ctmp);
00100        }
00101        if (ldot > FTINY && (np->specfl&(SP_REFL|SP_BADU)) == SP_REFL) {
00102               /*
00103                *  Compute specular reflection coefficient using
00104                *  anisotropic gaussian distribution model.
00105                */
00106                                           /* add source width if flat */
00107               if (np->specfl & SP_FLAT)
00108                      au2 = av2 = omega * (0.25/PI);
00109               else
00110                      au2 = av2 = 0.0;
00111               au2 += np->u_alpha*np->u_alpha;
00112               av2 += np->v_alpha*np->v_alpha;
00113                                           /* half vector */
00114               h[0] = ldir[0] - np->rp->rdir[0];
00115               h[1] = ldir[1] - np->rp->rdir[1];
00116               h[2] = ldir[2] - np->rp->rdir[2];
00117                                           /* ellipse */
00118               dtmp1 = DOT(np->u, h);
00119               dtmp1 *= dtmp1 / au2;
00120               dtmp2 = DOT(np->v, h);
00121               dtmp2 *= dtmp2 / av2;
00122                                           /* gaussian */
00123               dtmp = DOT(np->pnorm, h);
00124               dtmp = (dtmp1 + dtmp2) / (dtmp*dtmp);
00125               dtmp = exp(-dtmp) / (4.0*PI * np->pdot * sqrt(au2*av2));
00126                                           /* worth using? */
00127               if (dtmp > FTINY) {
00128                      copycolor(ctmp, np->scolor);
00129                      dtmp *= omega;
00130                      scalecolor(ctmp, dtmp);
00131                      addcolor(cval, ctmp);
00132               }
00133        }
00134        if (ldot < -FTINY && np->tdiff > FTINY) {
00135               /*
00136                *  Compute diffuse transmission.
00137                */
00138               copycolor(ctmp, np->mcolor);
00139               dtmp = -ldot * omega * np->tdiff * (1.0/PI);
00140               scalecolor(ctmp, dtmp);
00141               addcolor(cval, ctmp);
00142        }
00143        if (ldot < -FTINY && (np->specfl&(SP_TRAN|SP_BADU)) == SP_TRAN) {
00144               /*
00145                *  Compute specular transmission.  Specular transmission
00146                *  is always modified by material color.
00147                */
00148                                           /* roughness + source */
00149               au2 = av2 = omega * (1.0/PI);
00150               au2 += np->u_alpha*np->u_alpha;
00151               av2 += np->v_alpha*np->v_alpha;
00152                                           /* "half vector" */
00153               h[0] = ldir[0] - np->prdir[0];
00154               h[1] = ldir[1] - np->prdir[1];
00155               h[2] = ldir[2] - np->prdir[2];
00156               dtmp = DOT(h,h);
00157               if (dtmp > FTINY*FTINY) {
00158                      dtmp1 = DOT(h,np->pnorm);
00159                      dtmp = 1.0 - dtmp1*dtmp1/dtmp;
00160                      if (dtmp > FTINY*FTINY) {
00161                             dtmp1 = DOT(h,np->u);
00162                             dtmp1 *= dtmp1 / au2;
00163                             dtmp2 = DOT(h,np->v);
00164                             dtmp2 *= dtmp2 / av2;
00165                             dtmp = (dtmp1 + dtmp2) / dtmp;
00166                      }
00167               } else
00168                      dtmp = 0.0;
00169                                           /* gaussian */
00170               dtmp = exp(-dtmp) / (PI * np->pdot * sqrt(au2*av2));
00171                                           /* worth using? */
00172               if (dtmp > FTINY) {
00173                      copycolor(ctmp, np->mcolor);
00174                      dtmp *= np->tspec * omega;
00175                      scalecolor(ctmp, dtmp);
00176                      addcolor(cval, ctmp);
00177               }
00178        }
00179 }
00180 
00181 
00182 extern int
00183 m_aniso(                    /* shade ray that hit something anisotropic */
00184        register OBJREC  *m,
00185        register RAY  *r
00186 )
00187 {
00188        ANISODAT  nd;
00189        COLOR  ctmp;
00190        register int  i;
00191                                           /* easy shadow test */
00192        if (r->crtype & SHADOW)
00193               return(1);
00194 
00195        if (m->oargs.nfargs != (m->otype == MAT_TRANS2 ? 8 : 6))
00196               objerror(m, USER, "bad number of real arguments");
00197                                           /* check for back side */
00198        if (r->rod < 0.0) {
00199               if (!backvis && m->otype != MAT_TRANS2) {
00200                      raytrans(r);
00201                      return(1);
00202               }
00203               raytexture(r, m->omod);
00204               flipsurface(r);                    /* reorient if backvis */
00205        } else
00206               raytexture(r, m->omod);
00207                                           /* get material color */
00208        nd.mp = m;
00209        nd.rp = r;
00210        setcolor(nd.mcolor, m->oargs.farg[0],
00211                         m->oargs.farg[1],
00212                         m->oargs.farg[2]);
00213                                           /* get roughness */
00214        nd.specfl = 0;
00215        nd.u_alpha = m->oargs.farg[4];
00216        nd.v_alpha = m->oargs.farg[5];
00217        if (nd.u_alpha < FTINY || nd.v_alpha <= FTINY)
00218               objerror(m, USER, "roughness too small");
00219 
00220        nd.pdot = raynormal(nd.pnorm, r);  /* perturb normal */
00221        if (nd.pdot < .001)
00222               nd.pdot = .001;                    /* non-zero for diraniso() */
00223        multcolor(nd.mcolor, r->pcol);            /* modify material color */
00224                                           /* get specular component */
00225        if ((nd.rspec = m->oargs.farg[3]) > FTINY) {
00226               nd.specfl |= SP_REFL;
00227                                           /* compute specular color */
00228               if (m->otype == MAT_METAL2)
00229                      copycolor(nd.scolor, nd.mcolor);
00230               else
00231                      setcolor(nd.scolor, 1.0, 1.0, 1.0);
00232               scalecolor(nd.scolor, nd.rspec);
00233                                           /* check threshold */
00234               if (specthresh >= nd.rspec-FTINY)
00235                      nd.specfl |= SP_RBLT;
00236                                           /* compute refl. direction */
00237               for (i = 0; i < 3; i++)
00238                      nd.vrefl[i] = r->rdir[i] + 2.0*nd.pdot*nd.pnorm[i];
00239               if (DOT(nd.vrefl, r->ron) <= FTINY)       /* penetration? */
00240                      for (i = 0; i < 3; i++)            /* safety measure */
00241                             nd.vrefl[i] = r->rdir[i] + 2.*r->rod*r->ron[i];
00242        }
00243                                           /* compute transmission */
00244        if (m->otype == MAT_TRANS2) {
00245               nd.trans = m->oargs.farg[6]*(1.0 - nd.rspec);
00246               nd.tspec = nd.trans * m->oargs.farg[7];
00247               nd.tdiff = nd.trans - nd.tspec;
00248               if (nd.tspec > FTINY) {
00249                      nd.specfl |= SP_TRAN;
00250                                                  /* check threshold */
00251                      if (specthresh >= nd.tspec-FTINY)
00252                             nd.specfl |= SP_TBLT;
00253                      if (DOT(r->pert,r->pert) <= FTINY*FTINY) {
00254                             VCOPY(nd.prdir, r->rdir);
00255                      } else {
00256                             for (i = 0; i < 3; i++)            /* perturb */
00257                                    nd.prdir[i] = r->rdir[i] - r->pert[i];
00258                             if (DOT(nd.prdir, r->ron) < -FTINY)
00259                                    normalize(nd.prdir); /* OK */
00260                             else
00261                                    VCOPY(nd.prdir, r->rdir);
00262                      }
00263               }
00264        } else
00265               nd.tdiff = nd.tspec = nd.trans = 0.0;
00266 
00267                                           /* diffuse reflection */
00268        nd.rdiff = 1.0 - nd.trans - nd.rspec;
00269 
00270        if (r->ro != NULL && isflat(r->ro->otype))
00271               nd.specfl |= SP_FLAT;
00272 
00273        getacoords(r, &nd);                /* set up coordinates */
00274 
00275        if (nd.specfl & (SP_REFL|SP_TRAN) && !(nd.specfl & SP_BADU))
00276               agaussamp(r, &nd);
00277 
00278        if (nd.rdiff > FTINY) {            /* ambient from this side */
00279               copycolor(ctmp, nd.mcolor); /* modified by material color */
00280               if (nd.specfl & SP_RBLT)
00281                      scalecolor(ctmp, 1.0-nd.trans);
00282               else
00283                      scalecolor(ctmp, nd.rdiff);
00284               multambient(ctmp, r, nd.pnorm);
00285               addcolor(r->rcol, ctmp);    /* add to returned color */
00286        }
00287        if (nd.tdiff > FTINY) {            /* ambient from other side */
00288               FVECT  bnorm;
00289 
00290               flipsurface(r);
00291               bnorm[0] = -nd.pnorm[0];
00292               bnorm[1] = -nd.pnorm[1];
00293               bnorm[2] = -nd.pnorm[2];
00294               copycolor(ctmp, nd.mcolor); /* modified by color */
00295               if (nd.specfl & SP_TBLT)
00296                      scalecolor(ctmp, nd.trans);
00297               else
00298                      scalecolor(ctmp, nd.tdiff);
00299               multambient(ctmp, r, bnorm);
00300               addcolor(r->rcol, ctmp);
00301               flipsurface(r);
00302        }
00303                                    /* add direct component */
00304        direct(r, diraniso, &nd);
00305 
00306        return(1);
00307 }
00308 
00309 
00310 static void
00311 getacoords(          /* set up coordinate system */
00312        RAY  *r,
00313        register ANISODAT  *np
00314 )
00315 {
00316        register MFUNC  *mf;
00317        register int  i;
00318 
00319        mf = getfunc(np->mp, 3, 0x7, 1);
00320        setfunc(np->mp, r);
00321        errno = 0;
00322        for (i = 0; i < 3; i++)
00323               np->u[i] = evalue(mf->ep[i]);
00324        if (errno == EDOM || errno == ERANGE) {
00325               objerror(np->mp, WARNING, "compute error");
00326               np->specfl |= SP_BADU;
00327               return;
00328        }
00329        if (mf->f != &unitxf)
00330               multv3(np->u, np->u, mf->f->xfm);
00331        fcross(np->v, np->pnorm, np->u);
00332        if (normalize(np->v) == 0.0) {
00333               objerror(np->mp, WARNING, "illegal orientation vector");
00334               np->specfl |= SP_BADU;
00335               return;
00336        }
00337        fcross(np->u, np->v, np->pnorm);
00338 }
00339 
00340 
00341 static void
00342 agaussamp(           /* sample anisotropic gaussian specular */
00343        RAY  *r,
00344        register ANISODAT  *np
00345 )
00346 {
00347        RAY  sr;
00348        FVECT  h;
00349        double  rv[2];
00350        double  d, sinp, cosp;
00351        int  niter;
00352        register int  i;
00353                                    /* compute reflection */
00354        if ((np->specfl & (SP_REFL|SP_RBLT)) == SP_REFL &&
00355                      rayorigin(&sr, SPECULAR, r, np->scolor) == 0) {
00356               dimlist[ndims++] = (int)np->mp;
00357               for (niter = 0; niter < MAXITER; niter++) {
00358                      if (niter)
00359                             d = frandom();
00360                      else
00361                             d = urand(ilhash(dimlist,ndims)+samplendx);
00362                      multisamp(rv, 2, d);
00363                      d = 2.0*PI * rv[0];
00364                      cosp = tcos(d) * np->u_alpha;
00365                      sinp = tsin(d) * np->v_alpha;
00366                      d = sqrt(cosp*cosp + sinp*sinp);
00367                      cosp /= d;
00368                      sinp /= d;
00369                      rv[1] = 1.0 - specjitter*rv[1];
00370                      if (rv[1] <= FTINY)
00371                             d = 1.0;
00372                      else
00373                             d = sqrt(-log(rv[1]) /
00374                                    (cosp*cosp/(np->u_alpha*np->u_alpha) +
00375                                     sinp*sinp/(np->v_alpha*np->v_alpha)));
00376                      for (i = 0; i < 3; i++)
00377                             h[i] = np->pnorm[i] +
00378                                    d*(cosp*np->u[i] + sinp*np->v[i]);
00379                      d = -2.0 * DOT(h, r->rdir) / (1.0 + d*d);
00380                      for (i = 0; i < 3; i++)
00381                             sr.rdir[i] = r->rdir[i] + d*h[i];
00382                      if (DOT(sr.rdir, r->ron) > FTINY) {
00383                             rayvalue(&sr);
00384                             multcolor(sr.rcol, sr.rcoef);
00385                             addcolor(r->rcol, sr.rcol);
00386                             break;
00387                      }
00388               }
00389               ndims--;
00390        }
00391                                    /* compute transmission */
00392        copycolor(sr.rcoef, np->mcolor);          /* modify by material color */
00393        scalecolor(sr.rcoef, np->tspec);
00394        if ((np->specfl & (SP_TRAN|SP_TBLT)) == SP_TRAN &&
00395                      rayorigin(&sr, SPECULAR, r, sr.rcoef) == 0) {
00396               dimlist[ndims++] = (int)np->mp;
00397               for (niter = 0; niter < MAXITER; niter++) {
00398                      if (niter)
00399                             d = frandom();
00400                      else
00401                             d = urand(ilhash(dimlist,ndims)+1823+samplendx);
00402                      multisamp(rv, 2, d);
00403                      d = 2.0*PI * rv[0];
00404                      cosp = tcos(d) * np->u_alpha;
00405                      sinp = tsin(d) * np->v_alpha;
00406                      d = sqrt(cosp*cosp + sinp*sinp);
00407                      cosp /= d;
00408                      sinp /= d;
00409                      rv[1] = 1.0 - specjitter*rv[1];
00410                      if (rv[1] <= FTINY)
00411                             d = 1.0;
00412                      else
00413                             d = sqrt(-log(rv[1]) /
00414                                    (cosp*cosp/(np->u_alpha*np->u_alpha) +
00415                                     sinp*sinp/(np->v_alpha*np->v_alpha)));
00416                      for (i = 0; i < 3; i++)
00417                             sr.rdir[i] = np->prdir[i] +
00418                                           d*(cosp*np->u[i] + sinp*np->v[i]);
00419                      if (DOT(sr.rdir, r->ron) < -FTINY) {
00420                             normalize(sr.rdir);  /* OK, normalize */
00421                             rayvalue(&sr);
00422                             multcolor(sr.rcol, sr.rcoef);
00423                             addcolor(r->rcol, sr.rcol);
00424                             break;
00425                      }
00426               }
00427               ndims--;
00428        }
00429 }