home *** CD-ROM | disk | FTP | other *** search
/ DP Tool Club 8 / CDASC08.ISO / NEWS / RADIANCE / SRC / RT / ANISO.C < prev    next >
C/C++ Source or Header  |  1993-10-07  |  12KB  |  408 lines
  1. /* Copyright (c) 1992 Regents of the University of California */
  2.  
  3. #ifndef lint
  4. static char SCCSid[] = "@(#)aniso.c 2.22 10/16/92 LBL";
  5. #endif
  6.  
  7. /*
  8.  *  Shading functions for anisotropic materials.
  9.  */
  10.  
  11. #include  "ray.h"
  12.  
  13. #include  "otypes.h"
  14.  
  15. #include  "func.h"
  16.  
  17. #include  "random.h"
  18.  
  19. extern double  specthresh;        /* specular sampling threshold */
  20. extern double  specjitter;        /* specular sampling jitter */
  21.  
  22. /*
  23.  *    This routine implements the anisotropic Gaussian
  24.  *  model described by Ward in Siggraph `92 article.
  25.  *    We orient the surface towards the incoming ray, so a single
  26.  *  surface can be used to represent an infinitely thin object.
  27.  *
  28.  *  Arguments for MAT_PLASTIC2 and MAT_METAL2 are:
  29.  *  4+ ux    uy    uz    funcfile    [transform...]
  30.  *  0
  31.  *  6  red    grn    blu    specular-frac.    u-facet-slope v-facet-slope
  32.  *
  33.  *  Real arguments for MAT_TRANS2 are:
  34.  *  8  red     grn    blu    rspec    u-rough    v-rough    trans    tspec
  35.  */
  36.  
  37. #define  BSPEC(m)    (6.0)        /* specularity parameter b */
  38.  
  39.                 /* specularity flags */
  40. #define  SP_REFL    01        /* has reflected specular component */
  41. #define  SP_TRAN    02        /* has transmitted specular */
  42. #define  SP_FLAT    04        /* reflecting surface is flat */
  43. #define  SP_RBLT    010        /* reflection below sample threshold */
  44. #define  SP_TBLT    020        /* transmission below threshold */
  45. #define  SP_BADU    040        /* bad u direction calculation */
  46.  
  47. typedef struct {
  48.     OBJREC  *mp;        /* material pointer */
  49.     RAY  *rp;        /* ray pointer */
  50.     short  specfl;        /* specularity flags, defined above */
  51.     COLOR  mcolor;        /* color of this material */
  52.     COLOR  scolor;        /* color of specular component */
  53.     FVECT  vrefl;        /* vector in reflected direction */
  54.     FVECT  prdir;        /* vector in transmitted direction */
  55.     FVECT  u, v;        /* u and v vectors orienting anisotropy */
  56.     double  u_alpha;    /* u roughness */
  57.     double  v_alpha;    /* v roughness */
  58.     double  rdiff, rspec;    /* reflected specular, diffuse */
  59.     double  trans;        /* transmissivity */
  60.     double  tdiff, tspec;    /* transmitted specular, diffuse */
  61.     FVECT  pnorm;        /* perturbed surface normal */
  62.     double  pdot;        /* perturbed dot product */
  63. }  ANISODAT;        /* anisotropic material data */
  64.  
  65.  
  66. diraniso(cval, np, ldir, omega)        /* compute source contribution */
  67. COLOR  cval;            /* returned coefficient */
  68. register ANISODAT  *np;        /* material data */
  69. FVECT  ldir;            /* light source direction */
  70. double  omega;            /* light source size */
  71. {
  72.     double  ldot;
  73.     double  dtmp, dtmp1, dtmp2;
  74.     FVECT  h;
  75.     double  au2, av2;
  76.     COLOR  ctmp;
  77.  
  78.     setcolor(cval, 0.0, 0.0, 0.0);
  79.  
  80.     ldot = DOT(np->pnorm, ldir);
  81.  
  82.     if (ldot < 0.0 ? np->trans <= FTINY : np->trans >= 1.0-FTINY)
  83.         return;        /* wrong side */
  84.  
  85.     if (ldot > FTINY && np->rdiff > FTINY) {
  86.         /*
  87.          *  Compute and add diffuse reflected component to returned
  88.          *  color.  The diffuse reflected component will always be
  89.          *  modified by the color of the material.
  90.          */
  91.         copycolor(ctmp, np->mcolor);
  92.         dtmp = ldot * omega * np->rdiff / PI;
  93.         scalecolor(ctmp, dtmp);
  94.         addcolor(cval, ctmp);
  95.     }
  96.     if (ldot > FTINY && (np->specfl&(SP_REFL|SP_BADU)) == SP_REFL) {
  97.         /*
  98.          *  Compute specular reflection coefficient using
  99.          *  anisotropic gaussian distribution model.
  100.          */
  101.                         /* add source width if flat */
  102.         if (np->specfl & SP_FLAT)
  103.             au2 = av2 = omega/(4.0*PI);
  104.         else
  105.             au2 = av2 = 0.0;
  106.         au2 += np->u_alpha*np->u_alpha;
  107.         av2 += np->v_alpha*np->v_alpha;
  108.                         /* half vector */
  109.         h[0] = ldir[0] - np->rp->rdir[0];
  110.         h[1] = ldir[1] - np->rp->rdir[1];
  111.         h[2] = ldir[2] - np->rp->rdir[2];
  112.         normalize(h);
  113.                         /* ellipse */
  114.         dtmp1 = DOT(np->u, h);
  115.         dtmp1 *= dtmp1 / au2;
  116.         dtmp2 = DOT(np->v, h);
  117.         dtmp2 *= dtmp2 / av2;
  118.                         /* gaussian */
  119.         dtmp = (dtmp1 + dtmp2) / (1.0 + DOT(np->pnorm, h));
  120.         dtmp = exp(-2.0*dtmp) * (1.0/4.0/PI)
  121.                 * sqrt(ldot/(np->pdot*au2*av2));
  122.                         /* worth using? */
  123.         if (dtmp > FTINY) {
  124.             copycolor(ctmp, np->scolor);
  125.             dtmp *= omega;
  126.             scalecolor(ctmp, dtmp);
  127.             addcolor(cval, ctmp);
  128.         }
  129.     }
  130.     if (ldot < -FTINY && np->tdiff > FTINY) {
  131.         /*
  132.          *  Compute diffuse transmission.
  133.          */
  134.         copycolor(ctmp, np->mcolor);
  135.         dtmp = -ldot * omega * np->tdiff / PI;
  136.         scalecolor(ctmp, dtmp);
  137.         addcolor(cval, ctmp);
  138.     }
  139.     if (ldot < -FTINY && (np->specfl&(SP_TRAN|SP_BADU)) == SP_TRAN) {
  140.         /*
  141.          *  Compute specular transmission.  Specular transmission
  142.          *  is always modified by material color.
  143.          */
  144.                         /* roughness + source */
  145.         au2 = av2 = omega / PI;
  146.         au2 += np->u_alpha*np->u_alpha;
  147.         av2 += np->v_alpha*np->v_alpha;
  148.                         /* "half vector" */
  149.         h[0] = ldir[0] - np->prdir[0];
  150.         h[1] = ldir[1] - np->prdir[1];
  151.         h[2] = ldir[2] - np->prdir[2];
  152.         dtmp = DOT(h,h);
  153.         if (dtmp > FTINY*FTINY) {
  154.             dtmp1 = DOT(h,np->pnorm);
  155.             dtmp = 1.0 - dtmp1*dtmp1/dtmp;
  156.             if (dtmp > FTINY*FTINY) {
  157.                 dtmp1 = DOT(h,np->u);
  158.                 dtmp1 = dtmp1*dtmp1 / au2;
  159.                 dtmp2 = DOT(h,np->v);
  160.                 dtmp2 = dtmp2*dtmp2 / av2;
  161.                 dtmp = (dtmp1 + dtmp2) / dtmp;
  162.             }
  163.         } else
  164.             dtmp = 0.0;
  165.                         /* gaussian */
  166.         dtmp = exp(-dtmp) * (1.0/PI)
  167.                 * sqrt(-ldot/(np->pdot*au2*av2));
  168.                         /* worth using? */
  169.         if (dtmp > FTINY) {
  170.             copycolor(ctmp, np->mcolor);
  171.             dtmp *= np->tspec * omega;
  172.             scalecolor(ctmp, dtmp);
  173.             addcolor(cval, ctmp);
  174.         }
  175.     }
  176. }
  177.  
  178.  
  179. m_aniso(m, r)            /* shade ray that hit something anisotropic */
  180. register OBJREC  *m;
  181. register RAY  *r;
  182. {
  183.     ANISODAT  nd;
  184.     double  dtmp;
  185.     COLOR  ctmp;
  186.     register int  i;
  187.                         /* easy shadow test */
  188.     if (r->crtype & SHADOW)
  189.         return;
  190.  
  191.     if (m->oargs.nfargs != (m->otype == MAT_TRANS2 ? 8 : 6))
  192.         objerror(m, USER, "bad number of real arguments");
  193.     nd.mp = m;
  194.     nd.rp = r;
  195.                         /* get material color */
  196.     setcolor(nd.mcolor, m->oargs.farg[0],
  197.                m->oargs.farg[1],
  198.                m->oargs.farg[2]);
  199.                         /* get roughness */
  200.     nd.specfl = 0;
  201.     nd.u_alpha = m->oargs.farg[4];
  202.     nd.v_alpha = m->oargs.farg[5];
  203.     if (nd.u_alpha < FTINY || nd.v_alpha <= FTINY)
  204.         objerror(m, USER, "roughness too small");
  205.                         /* reorient if necessary */
  206.     if (r->rod < 0.0)
  207.         flipsurface(r);
  208.                         /* get modifiers */
  209.     raytexture(r, m->omod);
  210.     nd.pdot = raynormal(nd.pnorm, r);    /* perturb normal */
  211.     if (nd.pdot < .001)
  212.         nd.pdot = .001;            /* non-zero for diraniso() */
  213.     multcolor(nd.mcolor, r->pcol);        /* modify material color */
  214.                         /* get specular component */
  215.     if ((nd.rspec = m->oargs.farg[3]) > FTINY) {
  216.         nd.specfl |= SP_REFL;
  217.                         /* compute specular color */
  218.         if (m->otype == MAT_METAL2)
  219.             copycolor(nd.scolor, nd.mcolor);
  220.         else
  221.             setcolor(nd.scolor, 1.0, 1.0, 1.0);
  222.         scalecolor(nd.scolor, nd.rspec);
  223.                         /* improved model */
  224.         dtmp = exp(-BSPEC(m)*nd.pdot);
  225.         for (i = 0; i < 3; i++)
  226.             colval(nd.scolor,i) += (1.0-colval(nd.scolor,i))*dtmp;
  227.         nd.rspec += (1.0-nd.rspec)*dtmp;
  228.                         /* check threshold */
  229.         if (specthresh > FTINY &&
  230.                 (specthresh >= 1.-FTINY ||
  231.                 specthresh + .05 - .1*frandom() > nd.rspec))
  232.             nd.specfl |= SP_RBLT;
  233.                         /* compute refl. direction */
  234.         for (i = 0; i < 3; i++)
  235.             nd.vrefl[i] = r->rdir[i] + 2.0*nd.pdot*nd.pnorm[i];
  236.         if (DOT(nd.vrefl, r->ron) <= FTINY)    /* penetration? */
  237.             for (i = 0; i < 3; i++)        /* safety measure */
  238.                 nd.vrefl[i] = r->rdir[i] + 2.*r->rod*r->ron[i];
  239.     }
  240.                         /* compute transmission */
  241.     if (m->otype == MAT_TRANS2) {
  242.         nd.trans = m->oargs.farg[6]*(1.0 - nd.rspec);
  243.         nd.tspec = nd.trans * m->oargs.farg[7];
  244.         nd.tdiff = nd.trans - nd.tspec;
  245.         if (nd.tspec > FTINY) {
  246.             nd.specfl |= SP_TRAN;
  247.                             /* check threshold */
  248.             if (specthresh > FTINY &&
  249.                     (specthresh >= 1.-FTINY ||
  250.                 specthresh + .05 - .1*frandom() > nd.tspec))
  251.                 nd.specfl |= SP_TBLT;
  252.             if (DOT(r->pert,r->pert) <= FTINY*FTINY) {
  253.                 VCOPY(nd.prdir, r->rdir);
  254.             } else {
  255.                 for (i = 0; i < 3; i++)        /* perturb */
  256.                     nd.prdir[i] = r->rdir[i] - r->pert[i];
  257.                 if (DOT(nd.prdir, r->ron) < -FTINY)
  258.                     normalize(nd.prdir);    /* OK */
  259.                 else
  260.                     VCOPY(nd.prdir, r->rdir);
  261.             }
  262.         }
  263.     } else
  264.         nd.tdiff = nd.tspec = nd.trans = 0.0;
  265.  
  266.                         /* diffuse reflection */
  267.     nd.rdiff = 1.0 - nd.trans - nd.rspec;
  268.  
  269.     if (r->ro != NULL && (r->ro->otype == OBJ_FACE ||
  270.             r->ro->otype == OBJ_RING))
  271.         nd.specfl |= SP_FLAT;
  272.  
  273.     getacoords(r, &nd);            /* set up coordinates */
  274.  
  275.     if (nd.specfl & (SP_REFL|SP_TRAN) && !(nd.specfl & SP_BADU))
  276.         agaussamp(r, &nd);
  277.  
  278.     if (nd.rdiff > FTINY) {        /* ambient from this side */
  279.         ambient(ctmp, r);
  280.         if (nd.specfl & SP_RBLT)
  281.             scalecolor(ctmp, 1.0-nd.trans);
  282.         else
  283.             scalecolor(ctmp, nd.rdiff);
  284.         multcolor(ctmp, nd.mcolor);    /* modified by material color */
  285.         addcolor(r->rcol, ctmp);    /* add to returned color */
  286.     }
  287.     if (nd.tdiff > FTINY) {        /* ambient from other side */
  288.         flipsurface(r);
  289.         ambient(ctmp, r);
  290.         if (nd.specfl & SP_TBLT)
  291.             scalecolor(ctmp, nd.trans);
  292.         else
  293.             scalecolor(ctmp, nd.tdiff);
  294.         multcolor(ctmp, nd.mcolor);    /* modified by color */
  295.         addcolor(r->rcol, ctmp);
  296.         flipsurface(r);
  297.     }
  298.                     /* add direct component */
  299.     direct(r, diraniso, &nd);
  300. }
  301.  
  302.  
  303. static
  304. getacoords(r, np)        /* set up coordinate system */
  305. RAY  *r;
  306. register ANISODAT  *np;
  307. {
  308.     register MFUNC  *mf;
  309.     register int  i;
  310.  
  311.     mf = getfunc(np->mp, 3, 0x7, 1);
  312.     setfunc(np->mp, r);
  313.     errno = 0;
  314.     for (i = 0; i < 3; i++)
  315.         np->u[i] = evalue(mf->ep[i]);
  316.     if (errno) {
  317.         objerror(np->mp, WARNING, "compute error");
  318.         np->specfl |= SP_BADU;
  319.         return;
  320.     }
  321.     if (mf->f != &unitxf)
  322.         multv3(np->u, np->u, mf->f->xfm);
  323.     fcross(np->v, np->pnorm, np->u);
  324.     if (normalize(np->v) == 0.0) {
  325.         objerror(np->mp, WARNING, "illegal orientation vector");
  326.         np->specfl |= SP_BADU;
  327.         return;
  328.     }
  329.     fcross(np->u, np->v, np->pnorm);
  330. }
  331.  
  332.  
  333. static
  334. agaussamp(r, np)        /* sample anisotropic gaussian specular */
  335. RAY  *r;
  336. register ANISODAT  *np;
  337. {
  338.     RAY  sr;
  339.     FVECT  h;
  340.     double  rv[2];
  341.     double  d, sinp, cosp;
  342.     register int  i;
  343.                     /* compute reflection */
  344.     if ((np->specfl & (SP_REFL|SP_RBLT)) == SP_REFL &&
  345.             rayorigin(&sr, r, SPECULAR, np->rspec) == 0) {
  346.         dimlist[ndims++] = (int)np->mp;
  347.         d = urand(ilhash(dimlist,ndims)+samplendx);
  348.         multisamp(rv, 2, d);
  349.         d = 2.0*PI * rv[0];
  350.         cosp = cos(d) * np->u_alpha;
  351.         sinp = sin(d) * np->v_alpha;
  352.         d = sqrt(cosp*cosp + sinp*sinp);
  353.         cosp /= d;
  354.         sinp /= d;
  355.         rv[1] = 1.0 - specjitter*rv[1];
  356.         if (rv[1] <= FTINY)
  357.             d = 1.0;
  358.         else
  359.             d = sqrt(-log(rv[1]) /
  360.                 (cosp*cosp/(np->u_alpha*np->u_alpha) +
  361.                  sinp*sinp/(np->v_alpha*np->v_alpha)));
  362.         for (i = 0; i < 3; i++)
  363.             h[i] = np->pnorm[i] +
  364.                 d*(cosp*np->u[i] + sinp*np->v[i]);
  365.         d = -2.0 * DOT(h, r->rdir) / (1.0 + d*d);
  366.         for (i = 0; i < 3; i++)
  367.             sr.rdir[i] = r->rdir[i] + d*h[i];
  368.         if (DOT(sr.rdir, r->ron) <= FTINY)    /* penetration? */
  369.             VCOPY(sr.rdir, np->vrefl);    /* jitter no good */
  370.         rayvalue(&sr);
  371.         multcolor(sr.rcol, np->scolor);
  372.         addcolor(r->rcol, sr.rcol);
  373.         ndims--;
  374.     }
  375.                     /* compute transmission */
  376.     if ((np->specfl & (SP_TRAN|SP_TBLT)) == SP_TRAN &&
  377.             rayorigin(&sr, r, SPECULAR, np->tspec) == 0) {
  378.         dimlist[ndims++] = (int)np->mp;
  379.         d = urand(ilhash(dimlist,ndims)+1823+samplendx);
  380.         multisamp(rv, 2, d);
  381.         d = 2.0*PI * rv[0];
  382.         cosp = cos(d) * np->u_alpha;
  383.         sinp = sin(d) * np->v_alpha;
  384.         d = sqrt(cosp*cosp + sinp*sinp);
  385.         cosp /= d;
  386.         sinp /= d;
  387.         rv[1] = 1.0 - specjitter*rv[1];
  388.         if (rv[1] <= FTINY)
  389.             d = 1.0;
  390.         else
  391.             d = sqrt(-log(rv[1]) /
  392.                 (cosp*cosp/(np->u_alpha*np->u_alpha) +
  393.                  sinp*sinp/(np->v_alpha*np->u_alpha)));
  394.         for (i = 0; i < 3; i++)
  395.             sr.rdir[i] = np->prdir[i] +
  396.                     d*(cosp*np->u[i] + sinp*np->v[i]);
  397.         if (DOT(sr.rdir, r->ron) < -FTINY)
  398.             normalize(sr.rdir);        /* OK, normalize */
  399.         else
  400.             VCOPY(sr.rdir, np->prdir);    /* else no jitter */
  401.         rayvalue(&sr);
  402.         scalecolor(sr.rcol, np->tspec);
  403.         multcolor(sr.rcol, np->mcolor);        /* modify by color */
  404.         addcolor(r->rcol, sr.rcol);
  405.         ndims--;
  406.     }
  407. }
  408.