Back to index

radiance  4R0+20100331
dielectric.c
Go to the documentation of this file.
00001 #ifndef lint
00002 static const char    RCSid[] = "$Id: dielectric.c,v 2.20 2005/04/19 01:15:06 greg Exp $";
00003 #endif
00004 /*
00005  *  dielectric.c - shading function for transparent materials.
00006  */
00007 
00008 #include "copyright.h"
00009 
00010 #include  "ray.h"
00011 #include  "otypes.h"
00012 #include  "rtotypes.h"
00013 
00014 #ifdef  DISPERSE
00015 #include  "source.h"
00016 static int disperse(OBJREC *m,RAY *r,FVECT vt,double tr,COLOR cet,COLOR abt);
00017 static int lambda(OBJREC  *m, FVECT  v2, FVECT  dv, FVECT  lr);
00018 #endif
00019 
00020 static double mylog(double  x);
00021 
00022 
00023 /*
00024  *     Explicit calculations for Fresnel's equation are performed,
00025  *  but only one square root computation is necessary.
00026  *     The index of refraction is given as a Hartmann equation
00027  *  with lambda0 equal to zero.  If the slope of Hartmann's
00028  *  equation is non-zero, the material disperses light upon
00029  *  refraction.  This condition is examined on rays traced to
00030  *  light sources.  If a ray is exiting a dielectric material, we
00031  *  check the sources to see if any would cause bright color to be
00032  *  directed to the viewer due to dispersion.  This gives colorful
00033  *  sparkle to crystals, etc.  (Only if DISPERSE is defined!)
00034  *
00035  *  Arguments for MAT_DIELECTRIC are:
00036  *     red    grn    blu    rndx   Hartmann
00037  *
00038  *  Arguments for MAT_INTERFACE are:
00039  *     red1   grn1   blu1   rndx1  red2   grn2   blu2   rndx2
00040  *
00041  *  The primaries are material transmission per unit length.
00042  *  MAT_INTERFACE uses dielectric1 for inside and dielectric2 for
00043  *  outside.
00044  */
00045 
00046 
00047 #define  MLAMBDA     500           /* mean lambda */
00048 #define  MAXLAMBDA   779           /* maximum lambda */
00049 #define  MINLAMBDA   380           /* minimum lambda */
00050 
00051 #define  MINCOS             0.997         /* minimum dot product for dispersion */
00052 
00053 
00054 static double
00055 mylog(        /* special log for extinction coefficients */
00056        double  x
00057 )
00058 {
00059        if (x < 1e-40)
00060               return(-100.);
00061        if (x >= 1.)
00062               return(0.);
00063        return(log(x));
00064 }
00065 
00066 
00067 extern int
00068 m_dielectric( /* color a ray which hit a dielectric interface */
00069        OBJREC  *m,
00070        register RAY  *r
00071 )
00072 {
00073        double  cos1, cos2, nratio;
00074        COLOR  ctrans;
00075        COLOR  talb;
00076        int  hastexture;
00077        double  transdist, transtest=0;
00078        double  mirdist, mirtest=0;
00079        int    flatsurface;
00080        double  refl, trans;
00081        FVECT  dnorm;
00082        double  d1, d2;
00083        RAY  p;
00084        register int  i;
00085 
00086        if (m->oargs.nfargs != (m->otype==MAT_DIELECTRIC ? 5 : 8))
00087               objerror(m, USER, "bad arguments");
00088 
00089        raytexture(r, m->omod);                   /* get modifiers */
00090 
00091        if ( (hastexture = DOT(r->pert,r->pert) > FTINY*FTINY) )
00092               cos1 = raynormal(dnorm, r); /* perturb normal */
00093        else {
00094               VCOPY(dnorm, r->ron);
00095               cos1 = r->rod;
00096        }
00097        flatsurface = !hastexture && r->ro != NULL && isflat(r->ro->otype);
00098 
00099                                           /* index of refraction */
00100        if (m->otype == MAT_DIELECTRIC)
00101               nratio = m->oargs.farg[3] + m->oargs.farg[4]/MLAMBDA;
00102        else
00103               nratio = m->oargs.farg[3] / m->oargs.farg[7];
00104        
00105        if (cos1 < 0.0) {                  /* inside */
00106               hastexture = -hastexture;
00107               cos1 = -cos1;
00108               dnorm[0] = -dnorm[0];
00109               dnorm[1] = -dnorm[1];
00110               dnorm[2] = -dnorm[2];
00111               setcolor(r->cext, -mylog(m->oargs.farg[0]*colval(r->pcol,RED)),
00112                              -mylog(m->oargs.farg[1]*colval(r->pcol,GRN)),
00113                              -mylog(m->oargs.farg[2]*colval(r->pcol,BLU)));
00114               setcolor(r->albedo, 0., 0., 0.);
00115               r->gecc = 0.;
00116               if (m->otype == MAT_INTERFACE) {
00117                      setcolor(ctrans,
00118                             -mylog(m->oargs.farg[4]*colval(r->pcol,RED)),
00119                             -mylog(m->oargs.farg[5]*colval(r->pcol,GRN)),
00120                             -mylog(m->oargs.farg[6]*colval(r->pcol,BLU)));
00121                      setcolor(talb, 0., 0., 0.);
00122               } else {
00123                      copycolor(ctrans, cextinction);
00124                      copycolor(talb, salbedo);
00125               }
00126        } else {                           /* outside */
00127               nratio = 1.0 / nratio;
00128 
00129               setcolor(ctrans, -mylog(m->oargs.farg[0]*colval(r->pcol,RED)),
00130                              -mylog(m->oargs.farg[1]*colval(r->pcol,GRN)),
00131                              -mylog(m->oargs.farg[2]*colval(r->pcol,BLU)));
00132               setcolor(talb, 0., 0., 0.);
00133               if (m->otype == MAT_INTERFACE) {
00134                      setcolor(r->cext,
00135                             -mylog(m->oargs.farg[4]*colval(r->pcol,RED)),
00136                             -mylog(m->oargs.farg[5]*colval(r->pcol,GRN)),
00137                             -mylog(m->oargs.farg[6]*colval(r->pcol,BLU)));
00138                      setcolor(r->albedo, 0., 0., 0.);
00139                      r->gecc = 0.;
00140               }
00141        }
00142 
00143        d2 = 1.0 - nratio*nratio*(1.0 - cos1*cos1);      /* compute cos theta2 */
00144 
00145        if (d2 < FTINY)                    /* total reflection */
00146 
00147               refl = 1.0;
00148 
00149        else {                      /* refraction occurs */
00150                                    /* compute Fresnel's equations */
00151               cos2 = sqrt(d2);
00152               d1 = cos1;
00153               d2 = nratio*cos2;
00154               d1 = (d1 - d2) / (d1 + d2);
00155               refl = d1 * d1;
00156 
00157               d1 = 1.0 / cos1;
00158               d2 = nratio / cos2;
00159               d1 = (d1 - d2) / (d1 + d2);
00160               refl += d1 * d1;
00161 
00162               refl *= 0.5;
00163               trans = 1.0 - refl;
00164 
00165               trans *= nratio*nratio;            /* solid angle ratio */
00166 
00167               setcolor(p.rcoef, trans, trans, trans);
00168 
00169               if (rayorigin(&p, REFRACTED, r, p.rcoef) == 0) {
00170 
00171                                           /* compute refracted ray */
00172                      d1 = nratio*cos1 - cos2;
00173                      for (i = 0; i < 3; i++)
00174                             p.rdir[i] = nratio*r->rdir[i] + d1*dnorm[i];
00175                                           /* accidental reflection? */
00176                      if (hastexture &&
00177                             DOT(p.rdir,r->ron)*hastexture >= -FTINY) {
00178                             d1 *= (double)hastexture;
00179                             for (i = 0; i < 3; i++)     /* ignore texture */
00180                                    p.rdir[i] = nratio*r->rdir[i] +
00181                                                  d1*r->ron[i];
00182                             normalize(p.rdir);   /* not exact */
00183                      }
00184 #ifdef  DISPERSE
00185                      if (m->otype != MAT_DIELECTRIC
00186                                    || r->rod > 0.0
00187                                    || r->crtype & SHADOW
00188                                    || !directvis
00189                                    || m->oargs.farg[4] == 0.0
00190                                    || !disperse(m, r, p.rdir,
00191                                                  trans, ctrans, talb))
00192 #endif
00193                      {
00194                             copycolor(p.cext, ctrans);
00195                             copycolor(p.albedo, talb);
00196                             rayvalue(&p);
00197                             multcolor(p.rcol, p.rcoef);
00198                             addcolor(r->rcol, p.rcol);
00199                                           /* virtual distance */
00200                             if (flatsurface ||
00201                                    (1.-FTINY <= nratio &&
00202                                           nratio <= 1.+FTINY)) {
00203                                    transtest = 2*bright(p.rcol);
00204                                    transdist = r->rot + p.rt;
00205                             }
00206                      }
00207               }
00208        }
00209        setcolor(p.rcoef, refl, refl, refl);
00210 
00211        if (!(r->crtype & SHADOW) &&
00212                      rayorigin(&p, REFLECTED, r, p.rcoef) == 0) {
00213 
00214                                    /* compute reflected ray */
00215               for (i = 0; i < 3; i++)
00216                      p.rdir[i] = r->rdir[i] + 2.0*cos1*dnorm[i];
00217                                    /* accidental penetration? */
00218               if (hastexture && DOT(p.rdir,r->ron)*hastexture <= FTINY)
00219                      for (i = 0; i < 3; i++)            /* ignore texture */
00220                             p.rdir[i] = r->rdir[i] + 2.0*r->rod*r->ron[i];
00221 
00222               rayvalue(&p);               /* reflected ray value */
00223 
00224               multcolor(p.rcol, p.rcoef); /* color contribution */
00225               addcolor(r->rcol, p.rcol);
00226                                           /* virtual distance */
00227               if (flatsurface) {
00228                      mirtest = 2*bright(p.rcol);
00229                      mirdist = r->rot + p.rt;
00230               }
00231        }
00232                             /* check distance to return */
00233        d1 = bright(r->rcol);
00234        if (transtest > d1)
00235               r->rt = transdist;
00236        else if (mirtest > d1)
00237               r->rt = mirdist;
00238                             /* rayvalue() computes absorption */
00239        return(1);
00240 }
00241 
00242 
00243 #ifdef  DISPERSE
00244 
00245 static int
00246 disperse(  /* check light sources for dispersion */
00247        OBJREC  *m,
00248        RAY  *r,
00249        FVECT  vt,
00250        double  tr,
00251        COLOR  cet,
00252        COLOR  abt
00253 )
00254 {
00255        RAY  sray;
00256        const RAY  *entray;
00257        FVECT  v1, v2, n1, n2;
00258        FVECT  dv, v2Xdv;
00259        double  v2Xdvv2Xdv;
00260        int  success = 0;
00261        SRCINDEX  si;
00262        FVECT  vtmp1, vtmp2;
00263        double  dtmp1, dtmp2;
00264        int  l1, l2;
00265        COLOR  ctmp;
00266        int  i;
00267        
00268        /*
00269         *     This routine computes dispersion to the first order using
00270         *  the following assumptions:
00271         *
00272         *     1) The dependency of the index of refraction on wavelength
00273         *            is approximated by Hartmann's equation with lambda0
00274         *            equal to zero.
00275         *     2) The entry and exit locations are constant with respect
00276         *            to dispersion.
00277         *
00278         *     The second assumption permits us to model dispersion without
00279         *  having to sample refracted directions.  We assume that the
00280         *  geometry inside the material is constant, and concern ourselves
00281         *  only with the relationship between the entering and exiting ray.
00282         *  We compute the first derivatives of the entering and exiting
00283         *  refraction with respect to the index of refraction.  This 
00284         *  is then used in a first order Taylor series to determine the
00285         *  index of refraction necessary to send the exiting ray to each
00286         *  light source.
00287         *     If an exiting ray hits a light source within the refraction
00288         *  boundaries, we sum all the frequencies over the disc of the
00289         *  light source to determine the resulting color.  A smaller light
00290         *  source will therefore exhibit a sharper spectrum.
00291         */
00292 
00293        if (!(r->crtype & REFRACTED)) {           /* ray started in material */
00294               VCOPY(v1, r->rdir);
00295               n1[0] = -r->rdir[0]; n1[1] = -r->rdir[1]; n1[2] = -r->rdir[2];
00296        } else {
00297                                           /* find entry point */
00298               for (entray = r; entray->rtype != REFRACTED;
00299                             entray = entray->parent)
00300                      ;
00301               entray = entray->parent;
00302               if (entray->crtype & REFRACTED)           /* too difficult */
00303                      return(0);
00304               VCOPY(v1, entray->rdir);
00305               VCOPY(n1, entray->ron);
00306        }
00307        VCOPY(v2, vt);                     /* exiting ray */
00308        VCOPY(n2, r->ron);
00309 
00310                                    /* first order dispersion approx. */
00311        dtmp1 = DOT(n1, v1);
00312        dtmp2 = DOT(n2, v2);
00313        for (i = 0; i < 3; i++)
00314               dv[i] = v1[i] + v2[i] - n1[i]/dtmp1 - n2[i]/dtmp2;
00315               
00316        if (DOT(dv, dv) <= FTINY)   /* null effect */
00317               return(0);
00318                                    /* compute plane normal */
00319        fcross(v2Xdv, v2, dv);
00320        v2Xdvv2Xdv = DOT(v2Xdv, v2Xdv);
00321 
00322                                    /* check sources */
00323        initsrcindex(&si);
00324        while (srcray(&sray, r, &si)) {
00325 
00326               if (DOT(sray.rdir, v2) < MINCOS)
00327                      continue;                   /* bad source */
00328                                           /* adjust source ray */
00329 
00330               dtmp1 = DOT(v2Xdv, sray.rdir) / v2Xdvv2Xdv;
00331               sray.rdir[0] -= dtmp1 * v2Xdv[0];
00332               sray.rdir[1] -= dtmp1 * v2Xdv[1];
00333               sray.rdir[2] -= dtmp1 * v2Xdv[2];
00334 
00335               l1 = lambda(m, v2, dv, sray.rdir); /* mean lambda */
00336 
00337               if (l1 > MAXLAMBDA || l1 < MINLAMBDA)     /* not visible */
00338                      continue;
00339                                           /* trace source ray */
00340               copycolor(sray.cext, cet);
00341               copycolor(sray.albedo, abt);
00342               normalize(sray.rdir);
00343               rayvalue(&sray);
00344               if (bright(sray.rcol) <= FTINY)    /* missed it */
00345                      continue;
00346               
00347               /*
00348                *     Compute spectral sum over diameter of source.
00349                *  First find directions for rays going to opposite
00350                *  sides of source, then compute wavelengths for each.
00351                */
00352                
00353               fcross(vtmp1, v2Xdv, sray.rdir);
00354               dtmp1 = sqrt(si.dom  / v2Xdvv2Xdv / PI);
00355 
00356                                                  /* compute first ray */
00357               for (i = 0; i < 3; i++)
00358                      vtmp2[i] = sray.rdir[i] + dtmp1*vtmp1[i];
00359 
00360               l1 = lambda(m, v2, dv, vtmp2);            /* first lambda */
00361               if (l1 < 0)
00362                      continue;
00363                                                  /* compute second ray */
00364               for (i = 0; i < 3; i++)
00365                      vtmp2[i] = sray.rdir[i] - dtmp1*vtmp1[i];
00366 
00367               l2 = lambda(m, v2, dv, vtmp2);            /* second lambda */
00368               if (l2 < 0)
00369                      continue;
00370                                    /* compute color from spectrum */
00371               if (l1 < l2)
00372                      spec_rgb(ctmp, l1, l2);
00373               else
00374                      spec_rgb(ctmp, l2, l1);
00375               multcolor(ctmp, sray.rcol);
00376               scalecolor(ctmp, tr);
00377               addcolor(r->rcol, ctmp);
00378               success++;
00379        }
00380        return(success);
00381 }
00382 
00383 
00384 static int
00385 lambda(                     /* compute lambda for material */
00386        register OBJREC  *m,
00387        FVECT  v2,
00388        FVECT  dv,
00389        FVECT  lr
00390 )
00391 {
00392        FVECT  lrXdv, v2Xlr;
00393        double  dtmp, denom;
00394        int  i;
00395 
00396        fcross(lrXdv, lr, dv);
00397        for (i = 0; i < 3; i++)     
00398               if (lrXdv[i] > FTINY || lrXdv[i] < -FTINY)
00399                      break;
00400        if (i >= 3)
00401               return(-1);
00402 
00403        fcross(v2Xlr, v2, lr);
00404 
00405        dtmp = m->oargs.farg[4] / MLAMBDA;
00406        denom = dtmp + v2Xlr[i]/lrXdv[i] * (m->oargs.farg[3] + dtmp);
00407 
00408        if (denom < FTINY)
00409               return(-1);
00410 
00411        return(m->oargs.farg[4] / denom);
00412 }
00413 
00414 #endif  /* DISPERSE */