Back to index

radiance  4R0+20100331
bsdf.c
Go to the documentation of this file.
00001 #ifndef lint
00002 static const char RCSid[] = "$Id: bsdf.c,v 2.2 2009/06/19 06:49:25 greg Exp $";
00003 #endif
00004 /*
00005  * Routines for handling BSDF data
00006  */
00007 
00008 #include "standard.h"
00009 #include "bsdf.h"
00010 #include "paths.h"
00011 #include "ezxml.h"
00012 #include <ctype.h>
00013 
00014 #define MAXLATS             46            /* maximum number of latitudes */
00015 
00016 /* BSDF angle specification */
00017 typedef struct {
00018        char   name[64];            /* basis name */
00019        int    nangles;             /* total number of directions */
00020        struct {
00021               float  tmin;                /* starting theta */
00022               short  nphis;               /* number of phis (0 term) */
00023        }      lat[MAXLATS+1];             /* latitudes */
00024 } ANGLE_BASIS;
00025 
00026 #define       MAXABASES     3             /* limit on defined bases */
00027 
00028 static ANGLE_BASIS   abase_list[MAXABASES] = {
00029        {
00030               "LBNL/Klems Full", 145,
00031               { {-5., 1},
00032               {5., 8},
00033               {15., 16},
00034               {25., 20},
00035               {35., 24},
00036               {45., 24},
00037               {55., 24},
00038               {65., 16},
00039               {75., 12},
00040               {90., 0} }
00041        }, {
00042               "LBNL/Klems Half", 73,
00043               { {-6.5, 1},
00044               {6.5, 8},
00045               {19.5, 12},
00046               {32.5, 16},
00047               {46.5, 20},
00048               {61.5, 12},
00049               {76.5, 4},
00050               {90., 0} }
00051        }, {
00052               "LBNL/Klems Quarter", 41,
00053               { {-9., 1},
00054               {9., 8},
00055               {27., 12},
00056               {46., 12},
00057               {66., 8},
00058               {90., 0} }
00059        }
00060 };
00061 
00062 static int    nabases = 3;  /* current number of defined bases */
00063 
00064 
00065 static int
00066 ab_getvec(           /* get vector for this angle basis index */
00067        FVECT v,
00068        int ndx,
00069        void *p
00070 )
00071 {
00072        ANGLE_BASIS  *ab = (ANGLE_BASIS *)p;
00073        int    li;
00074        double pol, azi, d;
00075        
00076        if ((ndx < 0) | (ndx >= ab->nangles))
00077               return(0);
00078        for (li = 0; ndx >= ab->lat[li].nphis; li++)
00079               ndx -= ab->lat[li].nphis;
00080        pol = PI/180.*0.5*(ab->lat[li].tmin + ab->lat[li+1].tmin);
00081        azi = 2.*PI*ndx/ab->lat[li].nphis;
00082        v[2] = d = cos(pol);
00083        d = sqrt(1. - d*d);  /* sin(pol) */
00084        v[0] = cos(azi)*d;
00085        v[1] = sin(azi)*d;
00086        return(1);
00087 }
00088 
00089 
00090 static int
00091 ab_getndx(           /* get index corresponding to the given vector */
00092        FVECT v,
00093        void *p
00094 )
00095 {
00096        ANGLE_BASIS  *ab = (ANGLE_BASIS *)p;
00097        int    li, ndx;
00098        double pol, azi, d;
00099 
00100        if ((v[2] < -1.0) | (v[2] > 1.0))
00101               return(-1);
00102        pol = 180.0/PI*acos(v[2]);
00103        azi = 180.0/PI*atan2(v[1], v[0]);
00104        if (azi < 0.0) azi += 360.0;
00105        for (li = 1; ab->lat[li].tmin <= pol; li++)
00106               if (!ab->lat[li].nphis)
00107                      return(-1);
00108        --li;
00109        ndx = (int)((1./360.)*azi*ab->lat[li].nphis + 0.5);
00110        if (ndx >= ab->lat[li].nphis) ndx = 0;
00111        while (li--)
00112               ndx += ab->lat[li].nphis;
00113        return(ndx);
00114 }
00115 
00116 
00117 static double
00118 ab_getohm(           /* get solid angle for this angle basis index */
00119        int ndx,
00120        void *p
00121 )
00122 {
00123        ANGLE_BASIS  *ab = (ANGLE_BASIS *)p;
00124        int    li;
00125        double theta, theta1;
00126        
00127        if ((ndx < 0) | (ndx >= ab->nangles))
00128               return(0);
00129        for (li = 0; ndx >= ab->lat[li].nphis; li++)
00130               ndx -= ab->lat[li].nphis;
00131        theta1 = PI/180. * ab->lat[li+1].tmin;
00132        if (ab->lat[li].nphis == 1) {             /* special case */
00133               if (ab->lat[li].tmin > FTINY)
00134                      error(USER, "unsupported BSDF coordinate system");
00135               return(2.*PI*(1. - cos(theta1)));
00136        }
00137        theta = PI/180. * ab->lat[li].tmin;
00138        return(2.*PI*(cos(theta) - cos(theta1))/(double)ab->lat[li].nphis);
00139 }
00140 
00141 
00142 static int
00143 ab_getvecR(          /* get reverse vector for this angle basis index */
00144        FVECT v,
00145        int ndx,
00146        void *p
00147 )
00148 {
00149        if (!ab_getvec(v, ndx, p))
00150               return(0);
00151 
00152        v[0] = -v[0];
00153        v[1] = -v[1];
00154        v[2] = -v[2];
00155 
00156        return(1);
00157 }
00158 
00159 
00160 static int
00161 ab_getndxR(          /* get index corresponding to the reverse vector */
00162        FVECT v,
00163        void *p
00164 )
00165 {
00166        FVECT  v2;
00167        
00168        v2[0] = -v[0];
00169        v2[1] = -v[1];
00170        v2[2] = -v[2];
00171 
00172        return ab_getndx(v2, p);
00173 }
00174 
00175 
00176 static void
00177 load_bsdf_data(             /* load BSDF distribution for this wavelength */
00178        struct BSDF_data *dp,
00179        ezxml_t wdb
00180 )
00181 {
00182        char  *cbasis = ezxml_txt(ezxml_child(wdb,"ColumnAngleBasis"));
00183        char  *rbasis = ezxml_txt(ezxml_child(wdb,"RowAngleBasis"));
00184        char  *sdata;
00185        int  i;
00186        
00187        if ((cbasis == NULL) | (rbasis == NULL)) {
00188               error(WARNING, "missing column/row basis for BSDF");
00189               return;
00190        }
00191        /* XXX need to add routines for loading in foreign bases */
00192        for (i = nabases; i--; )
00193               if (!strcmp(cbasis, abase_list[i].name)) {
00194                      dp->ninc = abase_list[i].nangles;
00195                      dp->ib_priv = (void *)&abase_list[i];
00196                      dp->ib_vec = ab_getvecR;
00197                      dp->ib_ndx = ab_getndxR;
00198                      dp->ib_ohm = ab_getohm;
00199                      break;
00200               }
00201        if (i < 0) {
00202               sprintf(errmsg, "unsupported ColumnAngleBasis '%s'", cbasis);
00203               error(WARNING, errmsg);
00204               return;
00205        }
00206        for (i = nabases; i--; )
00207               if (!strcmp(rbasis, abase_list[i].name)) {
00208                      dp->nout = abase_list[i].nangles;
00209                      dp->ob_priv = (void *)&abase_list[i];
00210                      dp->ob_vec = ab_getvec;
00211                      dp->ob_ndx = ab_getndx;
00212                      dp->ob_ohm = ab_getohm;
00213                      break;
00214               }
00215        if (i < 0) {
00216               sprintf(errmsg, "unsupported RowAngleBasis '%s'", cbasis);
00217               error(WARNING, errmsg);
00218               return;
00219        }
00220                             /* read BSDF data */
00221        sdata  = ezxml_txt(ezxml_child(wdb,"ScatteringData"));
00222        if (sdata == NULL) {
00223               error(WARNING, "missing BSDF ScatteringData");
00224               return;
00225        }
00226        dp->bsdf = (float *)malloc(sizeof(float)*dp->ninc*dp->nout);
00227        if (dp->bsdf == NULL)
00228               error(SYSTEM, "out of memory in load_bsdf_data");
00229        for (i = 0; i < dp->ninc*dp->nout; i++) {
00230               char  *sdnext = fskip(sdata);
00231               if (sdnext == NULL) {
00232                      error(WARNING, "bad/missing BSDF ScatteringData");
00233                      free(dp->bsdf); dp->bsdf = NULL;
00234                      return;
00235               }
00236               while (*sdnext && isspace(*sdnext))
00237                      sdnext++;
00238               if (*sdnext == ',') sdnext++;
00239               dp->bsdf[i] = atof(sdata);
00240               sdata = sdnext;
00241        }
00242        while (isspace(*sdata))
00243               sdata++;
00244        if (*sdata) {
00245               sprintf(errmsg, "%d extra characters after BSDF ScatteringData",
00246                             strlen(sdata));
00247               error(WARNING, errmsg);
00248        }
00249 }
00250 
00251 
00252 static int
00253 check_bsdf_data(     /* check that BSDF data is sane */
00254        struct BSDF_data *dp
00255 )
00256 {
00257        double        *omega_iarr, *omega_oarr;
00258        double        dom, contrib, hemi_total;
00259        int           nneg;
00260        FVECT         v;
00261        int           i, o;
00262 
00263        if (dp == NULL || dp->bsdf == NULL)
00264               return(0);
00265        omega_iarr = (double *)calloc(dp->ninc, sizeof(double));
00266        omega_oarr = (double *)calloc(dp->nout, sizeof(double));
00267        if ((omega_iarr == NULL) | (omega_oarr == NULL))
00268               error(SYSTEM, "out of memory in check_bsdf_data");
00269                                    /* incoming projected solid angles */
00270        hemi_total = .0;
00271        for (i = dp->ninc; i--; ) {
00272               dom = getBSDF_incohm(dp,i);
00273               if (dom <= .0) {
00274                      error(WARNING, "zero/negative incoming solid angle");
00275                      continue;
00276               }
00277               if (!getBSDF_incvec(v,dp,i) || v[2] > FTINY) {
00278                      error(WARNING, "illegal incoming BSDF direction");
00279                      free(omega_iarr); free(omega_oarr);
00280                      return(0);
00281               }
00282               hemi_total += omega_iarr[i] = dom * -v[2];
00283        }
00284        if ((hemi_total > 1.02*PI) | (hemi_total < 0.98*PI)) {
00285               sprintf(errmsg, "incoming BSDF hemisphere off by %.1f%%",
00286                             100.*(hemi_total/PI - 1.));
00287               error(WARNING, errmsg);
00288        }
00289        dom = PI / hemi_total;             /* fix normalization */
00290        for (i = dp->ninc; i--; )
00291               omega_iarr[i] *= dom;
00292                                    /* outgoing projected solid angles */
00293        hemi_total = .0;
00294        for (o = dp->nout; o--; ) {
00295               dom = getBSDF_outohm(dp,o);
00296               if (dom <= .0) {
00297                      error(WARNING, "zero/negative outgoing solid angle");
00298                      continue;
00299               }
00300               if (!getBSDF_outvec(v,dp,o) || v[2] < -FTINY) {
00301                      error(WARNING, "illegal outgoing BSDF direction");
00302                      free(omega_iarr); free(omega_oarr);
00303                      return(0);
00304               }
00305               hemi_total += omega_oarr[o] = dom * v[2];
00306        }
00307        if ((hemi_total > 1.02*PI) | (hemi_total < 0.98*PI)) {
00308               sprintf(errmsg, "outgoing BSDF hemisphere off by %.1f%%",
00309                             100.*(hemi_total/PI - 1.));
00310               error(WARNING, errmsg);
00311        }
00312        dom = PI / hemi_total;             /* fix normalization */
00313        for (o = dp->nout; o--; )
00314               omega_oarr[o] *= dom;
00315        nneg = 0;                   /* check outgoing totals */
00316        for (i = 0; i < dp->ninc; i++) {
00317               hemi_total = .0;
00318               for (o = dp->nout; o--; ) {
00319                      double f = BSDF_value(dp,i,o);
00320                      if (f >= .0)
00321                             hemi_total += f*omega_oarr[o];
00322                      else {
00323                             nneg += (f < -FTINY);
00324                             BSDF_value(dp,i,o) = .0f;
00325                      }
00326               }
00327               if (hemi_total > 1.02) {
00328                      sprintf(errmsg,
00329                      "incoming BSDF direction %d passes %.1f%% of light",
00330                                    i, 100.*hemi_total);
00331                      error(WARNING, errmsg);
00332               }
00333        }
00334        if (nneg) {
00335               sprintf(errmsg, "%d negative BSDF values (ignored)", nneg);
00336               error(WARNING, errmsg);
00337        }
00338                                    /* reverse roles and check again */
00339        for (o = 0; o < dp->nout; o++) {
00340               hemi_total = .0;
00341               for (i = dp->ninc; i--; )
00342                      hemi_total += BSDF_value(dp,i,o) * omega_iarr[i];
00343 
00344               if (hemi_total > 1.02) {
00345                      sprintf(errmsg,
00346                      "outgoing BSDF direction %d collects %.1f%% of light",
00347                                    o, 100.*hemi_total);
00348                      error(WARNING, errmsg);
00349               }
00350        }
00351        free(omega_iarr); free(omega_oarr);
00352        return(1);
00353 }
00354 
00355 struct BSDF_data *
00356 load_BSDF(           /* load BSDF data from file */
00357        char *fname
00358 )
00359 {
00360        char                 *path;
00361        ezxml_t                     fl, wtl, wld, wdb;
00362        struct BSDF_data     *dp;
00363        
00364        path = getpath(fname, getrlibpath(), R_OK);
00365        if (path == NULL) {
00366               sprintf(errmsg, "cannot find BSDF file \"%s\"", fname);
00367               error(WARNING, errmsg);
00368               return(NULL);
00369        }
00370        fl = ezxml_parse_file(path);
00371        if (fl == NULL) {
00372               sprintf(errmsg, "cannot open BSDF \"%s\"", path);
00373               error(WARNING, errmsg);
00374               return(NULL);
00375        }
00376        if (ezxml_error(fl)[0]) {
00377               sprintf(errmsg, "BSDF \"%s\" %s", path, ezxml_error(fl));
00378               error(WARNING, errmsg);
00379               ezxml_free(fl);
00380               return(NULL);
00381        }
00382        if (strcmp(ezxml_name(fl), "WindowElement")) {
00383               sprintf(errmsg,
00384                      "BSDF \"%s\": top level node not 'WindowElement'",
00385                             path);
00386               error(WARNING, errmsg);
00387               ezxml_free(fl);
00388               return(NULL);
00389        }
00390        wtl = ezxml_child(ezxml_child(fl, "Optical"), "Layer");
00391        dp = (struct BSDF_data *)calloc(1, sizeof(struct BSDF_data));
00392        for (wld = ezxml_child(wtl, "WavelengthData");
00393                             wld != NULL; wld = wld->next) {
00394               if (strcmp(ezxml_txt(ezxml_child(wld,"Wavelength")), "Visible"))
00395                      continue;
00396               wdb = ezxml_child(wld, "WavelengthDataBlock");
00397               if (wdb == NULL) continue;
00398               if (strcmp(ezxml_txt(ezxml_child(wdb,"WavelengthDataDirection")),
00399                                    "Transmission Front"))
00400                      continue;
00401               load_bsdf_data(dp, wdb);    /* load front BTDF */
00402               break;                      /* ignore the rest */
00403        }
00404        ezxml_free(fl);                           /* done with XML file */
00405        if (!check_bsdf_data(dp)) {
00406               sprintf(errmsg, "bad/missing BTDF data in \"%s\"", path);
00407               error(WARNING, errmsg);
00408               free_BSDF(dp);
00409               dp = NULL;
00410        }
00411        return(dp);
00412 }
00413 
00414 
00415 void
00416 free_BSDF(           /* free BSDF data structure */
00417        struct BSDF_data *b
00418 )
00419 {
00420        if (b == NULL)
00421               return;
00422        if (b->bsdf != NULL)
00423               free(b->bsdf);
00424        free(b);
00425 }
00426 
00427 
00428 int
00429 r_BSDF_incvec(              /* compute random input vector at given location */
00430        FVECT v,
00431        struct BSDF_data *b,
00432        int i,
00433        double rv,
00434        MAT4 xm
00435 )
00436 {
00437        FVECT  pert;
00438        double rad;
00439        int    j;
00440        
00441        if (!getBSDF_incvec(v, b, i))
00442               return(0);
00443        rad = sqrt(getBSDF_incohm(b, i) / PI);
00444        multisamp(pert, 3, rv);
00445        for (j = 0; j < 3; j++)
00446               v[j] += rad*(2.*pert[j] - 1.);
00447        if (xm != NULL)
00448               multv3(v, v, xm);
00449        return(normalize(v) != 0.0);
00450 }
00451 
00452 
00453 int
00454 r_BSDF_outvec(              /* compute random output vector at given location */
00455        FVECT v,
00456        struct BSDF_data *b,
00457        int o,
00458        double rv,
00459        MAT4 xm
00460 )
00461 {
00462        FVECT  pert;
00463        double rad;
00464        int    j;
00465        
00466        if (!getBSDF_outvec(v, b, o))
00467               return(0);
00468        rad = sqrt(getBSDF_outohm(b, o) / PI);
00469        multisamp(pert, 3, rv);
00470        for (j = 0; j < 3; j++)
00471               v[j] += rad*(2.*pert[j] - 1.);
00472        if (xm != NULL)
00473               multv3(v, v, xm);
00474        return(normalize(v) != 0.0);
00475 }
00476 
00477 
00478 #define  FEQ(a,b)    ((a)-(b) <= 1e-7 && (b)-(a) <= 1e-7)
00479 
00480 static int
00481 addrot(                     /* compute rotation (x,y,z) => (xp,yp,zp) */
00482        char *xfarg[],
00483        FVECT xp,
00484        FVECT yp,
00485        FVECT zp
00486 )
00487 {
00488        static char   bufs[3][16];
00489        int    bn = 0;
00490        char   **xfp = xfarg;
00491        double theta;
00492 
00493        if (yp[2]*yp[2] + zp[2]*zp[2] < 2.*FTINY*FTINY) {
00494               /* Special case for X' along Z-axis */
00495               theta = -atan2(yp[0], yp[1]);
00496               *xfp++ = "-ry";
00497               *xfp++ = xp[2] < 0.0 ? "90" : "-90";
00498               *xfp++ = "-rz";
00499               sprintf(bufs[bn], "%f", theta*(180./PI));
00500               *xfp++ = bufs[bn++];
00501               return(xfp - xfarg);
00502        }
00503        theta = atan2(yp[2], zp[2]);
00504        if (!FEQ(theta,0.0)) {
00505               *xfp++ = "-rx";
00506               sprintf(bufs[bn], "%f", theta*(180./PI));
00507               *xfp++ = bufs[bn++];
00508        }
00509        theta = asin(-xp[2]);
00510        if (!FEQ(theta,0.0)) {
00511               *xfp++ = "-ry";
00512               sprintf(bufs[bn], " %f", theta*(180./PI));
00513               *xfp++ = bufs[bn++];
00514        }
00515        theta = atan2(xp[1], xp[0]);
00516        if (!FEQ(theta,0.0)) {
00517               *xfp++ = "-rz";
00518               sprintf(bufs[bn], "%f", theta*(180./PI));
00519               *xfp++ = bufs[bn++];
00520        }
00521        *xfp = NULL;
00522        return(xfp - xfarg);
00523 }
00524 
00525 
00526 int
00527 getBSDF_xfm(         /* compute BSDF orient. -> world orient. transform */
00528        MAT4 xm,
00529        FVECT nrm,
00530        UpDir ud
00531 )
00532 {
00533        char   *xfargs[7];
00534        XF     myxf;
00535        FVECT  updir, xdest, ydest;
00536 
00537        updir[0] = updir[1] = updir[2] = 0.;
00538        switch (ud) {
00539        case UDzneg:
00540               updir[2] = -1.;
00541               break;
00542        case UDyneg:
00543               updir[1] = -1.;
00544               break;
00545        case UDxneg:
00546               updir[0] = -1.;
00547               break;
00548        case UDxpos:
00549               updir[0] = 1.;
00550               break;
00551        case UDypos:
00552               updir[1] = 1.;
00553               break;
00554        case UDzpos:
00555               updir[2] = 1.;
00556               break;
00557        case UDunknown:
00558               return(0);
00559        }
00560        fcross(xdest, updir, nrm);
00561        if (normalize(xdest) == 0.0)
00562               return(0);
00563        fcross(ydest, nrm, xdest);
00564        xf(&myxf, addrot(xfargs, xdest, ydest, nrm), xfargs);
00565        copymat4(xm, myxf.xfm);
00566        return(1);
00567 }