ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/bsdf.c
(Generate patch)

Comparing ray/src/common/bsdf.c (file contents):
Revision 2.3 by greg, Sat Jul 3 02:50:22 2010 UTC vs.
Revision 2.58 by greg, Thu May 14 19:20:13 2020 UTC

# Line 2 | Line 2
2   static const char RCSid[] = "$Id$";
3   #endif
4   /*
5 < * Routines for handling BSDF data
5 > *  bsdf.c
6 > *  
7 > *  Definitions for bidirectional scattering distribution functions.
8 > *
9 > *  Created by Greg Ward on 1/10/11.
10 > *
11   */
12  
13 < #include "standard.h"
14 < #include "bsdf.h"
15 < #include "paths.h"
16 < #include "ezxml.h"
13 > #define _USE_MATH_DEFINES
14 > #include <stdio.h>
15 > #include <stdlib.h>
16 > #include <string.h>
17 > #include <math.h>
18   #include <ctype.h>
19 + #include "ezxml.h"
20 + #include "hilbert.h"
21 + #include "bsdf.h"
22 + #include "bsdf_m.h"
23 + #include "bsdf_t.h"
24  
25 < #define MAXLATS         46              /* maximum number of latitudes */
25 > /* English ASCII strings corresponding to ennumerated errors */
26 > const char              *SDerrorEnglish[] = {
27 >                                "No error",
28 >                                "Memory error",
29 >                                "File input/output error",
30 >                                "File format error",
31 >                                "Illegal argument",
32 >                                "Invalid data",
33 >                                "Unsupported feature",
34 >                                "Internal program error",
35 >                                "Unknown error"
36 >                        };
37  
38 < /* BSDF angle specification */
39 < typedef struct {
18 <        char    name[64];               /* basis name */
19 <        int     nangles;                /* total number of directions */
20 <        struct {
21 <                float   tmin;                   /* starting theta */
22 <                short   nphis;                  /* number of phis (0 term) */
23 <        }       lat[MAXLATS+1];         /* latitudes */
24 < } ANGLE_BASIS;
38 > /* Pointer to error list in preferred language */
39 > const char              **SDerrorList = SDerrorEnglish;
40  
41 < #define MAXABASES       5               /* limit on defined bases */
41 > /* Additional information on last error (ASCII English) */
42 > char                    SDerrorDetail[256];
43  
44 < static ANGLE_BASIS      abase_list[MAXABASES] = {
45 <        {
30 <                "LBNL/Klems Full", 145,
31 <                { {-5., 1},
32 <                {5., 8},
33 <                {15., 16},
34 <                {25., 20},
35 <                {35., 24},
36 <                {45., 24},
37 <                {55., 24},
38 <                {65., 16},
39 <                {75., 12},
40 <                {90., 0} }
41 <        }, {
42 <                "LBNL/Klems Half", 73,
43 <                { {-6.5, 1},
44 <                {6.5, 8},
45 <                {19.5, 12},
46 <                {32.5, 16},
47 <                {46.5, 20},
48 <                {61.5, 12},
49 <                {76.5, 4},
50 <                {90., 0} }
51 <        }, {
52 <                "LBNL/Klems Quarter", 41,
53 <                { {-9., 1},
54 <                {9., 8},
55 <                {27., 12},
56 <                {46., 12},
57 <                {66., 8},
58 <                {90., 0} }
59 <        }
60 < };
44 > /* Empty distribution for getCDist() calls that fail for some reason */
45 > const SDCDst            SDemptyCD;
46  
47 < static int      nabases = 3;    /* current number of defined bases */
47 > /* Cache of loaded BSDFs */
48 > struct SDCache_s        *SDcacheList = NULL;
49  
50 < #define  FEQ(a,b)       ((a)-(b) <= 1e-7 && (b)-(a) <= 1e-7)
50 > /* Retain BSDFs in cache list? */
51 > int                     SDretainSet = SDretainNone;
52  
53 < // returns the name of the given tag
54 < #ifdef ezxml_name
68 < #undef ezxml_name
69 < static char *
70 < ezxml_name(ezxml_t xml)
71 < {
72 <        if (xml == NULL)
73 <                return(NULL);
74 <        return(xml->name);
75 < }
76 < #endif
53 > /* Maximum cache size for any given BSDF? */
54 > unsigned long           SDmaxCache = 0;         /* 0 == unlimited */
55  
56 < // returns the given tag's character content or empty string if none
57 < #ifdef ezxml_txt
58 < #undef ezxml_txt
81 < static char *
82 < ezxml_txt(ezxml_t xml)
56 > /* Report any error to the indicated stream */
57 > SDError
58 > SDreportError(SDError ec, FILE *fp)
59   {
60 <        if (xml == NULL)
61 <                return("");
62 <        return(xml->txt);
60 >        if (!ec)
61 >                return SDEnone;
62 >        if ((ec < SDEnone) | (ec > SDEunknown)) {
63 >                SDerrorDetail[0] = '\0';
64 >                ec = SDEunknown;
65 >        }
66 >        if (fp == NULL)
67 >                return ec;
68 >        fputs(SDerrorList[ec], fp);
69 >        if (SDerrorDetail[0]) {
70 >                fputs(": ", fp);
71 >                fputs(SDerrorDetail, fp);
72 >        }
73 >        fputc('\n', fp);
74 >        if (fp != stderr)
75 >                fflush(fp);
76 >        return ec;
77   }
88 #endif
78  
79 <
80 < static int
81 < ab_getvec(              /* get vector for this angle basis index */
93 <        FVECT v,
94 <        int ndx,
95 <        void *p
79 > static double
80 > to_meters(              /* return factor to convert given unit to meters */
81 >        const char *unit
82   )
83   {
84 <        ANGLE_BASIS  *ab = (ANGLE_BASIS *)p;
85 <        int     li;
86 <        double  pol, azi, d;
87 <        
88 <        if ((ndx < 0) | (ndx >= ab->nangles))
89 <                return(0);
90 <        for (li = 0; ndx >= ab->lat[li].nphis; li++)
91 <                ndx -= ab->lat[li].nphis;
106 <        pol = PI/180.*0.5*(ab->lat[li].tmin + ab->lat[li+1].tmin);
107 <        azi = 2.*PI*ndx/ab->lat[li].nphis;
108 <        v[2] = d = cos(pol);
109 <        d = sqrt(1. - d*d);     /* sin(pol) */
110 <        v[0] = cos(azi)*d;
111 <        v[1] = sin(azi)*d;
112 <        return(1);
84 >        if (unit == NULL) return(1.);           /* safe assumption? */
85 >        if (!strcasecmp(unit, "Meter")) return(1.);
86 >        if (!strcasecmp(unit, "Foot")) return(.3048);
87 >        if (!strcasecmp(unit, "Inch")) return(.0254);
88 >        if (!strcasecmp(unit, "Centimeter")) return(.01);
89 >        if (!strcasecmp(unit, "Millimeter")) return(.001);
90 >        sprintf(SDerrorDetail, "Unknown dimensional unit '%s'", unit);
91 >        return(-1.);
92   }
93  
94 <
95 < static int
96 < ab_getndx(              /* get index corresponding to the given vector */
118 <        FVECT v,
119 <        void *p
120 < )
94 > /* Load geometric dimensions and description (if any) */
95 > static SDError
96 > SDloadGeometry(SDData *sd, ezxml_t wtl)
97   {
98 <        ANGLE_BASIS  *ab = (ANGLE_BASIS *)p;
99 <        int     li, ndx;
100 <        double  pol, azi, d;
98 >        ezxml_t         node, matl, geom;
99 >        double          cfact;
100 >        const char      *fmt = NULL, *mgfstr;
101  
102 <        if ((v[2] < -1.0) | (v[2] > 1.0))
103 <                return(-1);
104 <        pol = 180.0/PI*acos(v[2]);
105 <        azi = 180.0/PI*atan2(v[1], v[0]);
106 <        if (azi < 0.0) azi += 360.0;
107 <        for (li = 1; ab->lat[li].tmin <= pol; li++)
108 <                if (!ab->lat[li].nphis)
109 <                        return(-1);
110 <        --li;
111 <        ndx = (int)((1./360.)*azi*ab->lat[li].nphis + 0.5);
112 <        if (ndx >= ab->lat[li].nphis) ndx = 0;
113 <        while (li--)
114 <                ndx += ab->lat[li].nphis;
115 <        return(ndx);
102 >        SDerrorDetail[0] = '\0';
103 >        sd->matn[0] = '\0'; sd->makr[0] = '\0';
104 >        sd->dim[0] = sd->dim[1] = sd->dim[2] = 0;
105 >        matl = ezxml_child(wtl, "Material");
106 >        if (matl != NULL) {                     /* get material info. */
107 >                if ((node = ezxml_child(matl, "Name")) != NULL) {
108 >                        strncpy(sd->matn, ezxml_txt(node), SDnameLn);
109 >                        if (sd->matn[SDnameLn-1])
110 >                                strcpy(sd->matn+(SDnameLn-4), "...");
111 >                }
112 >                if ((node = ezxml_child(matl, "Manufacturer")) != NULL) {
113 >                        strncpy(sd->makr, ezxml_txt(node), SDnameLn);
114 >                        if (sd->makr[SDnameLn-1])
115 >                                strcpy(sd->makr+(SDnameLn-4), "...");
116 >                }
117 >                if ((node = ezxml_child(matl, "Width")) != NULL)
118 >                        sd->dim[0] = atof(ezxml_txt(node)) *
119 >                                        to_meters(ezxml_attr(node, "unit"));
120 >                if ((node = ezxml_child(matl, "Height")) != NULL)
121 >                        sd->dim[1] = atof(ezxml_txt(node)) *
122 >                                        to_meters(ezxml_attr(node, "unit"));
123 >                if ((node = ezxml_child(matl, "Thickness")) != NULL)
124 >                        sd->dim[2] = atof(ezxml_txt(node)) *
125 >                                        to_meters(ezxml_attr(node, "unit"));
126 >                if ((sd->dim[0] < 0) | (sd->dim[1] < 0) | (sd->dim[2] < 0)) {
127 >                        if (!SDerrorDetail[0])
128 >                                sprintf(SDerrorDetail, "Negative dimension in \"%s\"",
129 >                                                        sd->name);
130 >                        return SDEdata;
131 >                }
132 >        }
133 >        sd->mgf = NULL;
134 >        geom = ezxml_child(wtl, "Geometry");
135 >        if (geom == NULL)                       /* no actual geometry? */
136 >                return SDEnone;
137 >        fmt = ezxml_attr(geom, "format");
138 >        if (fmt != NULL && strcasecmp(fmt, "MGF")) {
139 >                sprintf(SDerrorDetail,
140 >                        "Unrecognized geometry format '%s' in \"%s\"",
141 >                                        fmt, sd->name);
142 >                return SDEsupport;
143 >        }
144 >        if ((node = ezxml_child(geom, "MGFblock")) == NULL ||
145 >                        (mgfstr = ezxml_txt(node)) == NULL)
146 >                return SDEnone;
147 >        while (isspace(*mgfstr))
148 >                ++mgfstr;
149 >        if (!*mgfstr)
150 >                return SDEnone;
151 >        cfact = to_meters(ezxml_attr(node, "unit"));
152 >        if (cfact <= 0)
153 >                return SDEformat;
154 >        sd->mgf = (char *)malloc(strlen(mgfstr)+32);
155 >        if (sd->mgf == NULL) {
156 >                strcpy(SDerrorDetail, "Out of memory in SDloadGeometry");
157 >                return SDEmemory;
158 >        }
159 >        if (cfact < 0.99 || cfact > 1.01)
160 >                sprintf(sd->mgf, "xf -s %.5f\n%s\nxf\n", cfact, mgfstr);
161 >        else
162 >                strcpy(sd->mgf, mgfstr);
163 >        return SDEnone;
164   }
165  
166 + /* Load a BSDF struct from the given file (free first and keep name) */
167 + SDError
168 + SDloadFile(SDData *sd, const char *fname)
169 + {
170 +        SDError         lastErr;
171 +        ezxml_t         fl, wtl;
172  
173 < static double
174 < ab_getohm(              /* get solid angle for this angle basis index */
175 <        int ndx,
176 <        void *p
177 < )
173 >        if ((sd == NULL) | (fname == NULL || !*fname))
174 >                return SDEargument;
175 >                                /* free old data, keeping name */
176 >        SDfreeBSDF(sd);
177 >                                /* parse XML file */
178 >        fl = ezxml_parse_file(fname);
179 >        if (fl == NULL) {
180 >                sprintf(SDerrorDetail, "Cannot open BSDF \"%s\"", fname);
181 >                return SDEfile;
182 >        }
183 >        if (ezxml_error(fl)[0]) {
184 >                sprintf(SDerrorDetail, "BSDF \"%s\" %s", fname, ezxml_error(fl));
185 >                ezxml_free(fl);
186 >                return SDEformat;
187 >        }
188 >        if (strcmp(ezxml_name(fl), "WindowElement")) {
189 >                sprintf(SDerrorDetail,
190 >                        "BSDF \"%s\": top level node not 'WindowElement'",
191 >                                sd->name);
192 >                ezxml_free(fl);
193 >                return SDEformat;
194 >        }
195 >        wtl = ezxml_child(fl, "FileType");
196 >        if (wtl != NULL && strcmp(ezxml_txt(wtl), "BSDF")) {
197 >                sprintf(SDerrorDetail,
198 >                        "XML \"%s\": wrong FileType (must be 'BSDF')",
199 >                                sd->name);
200 >                ezxml_free(fl);
201 >                return SDEformat;
202 >        }
203 >        wtl = ezxml_child(ezxml_child(fl, "Optical"), "Layer");
204 >        if (wtl == NULL) {
205 >                sprintf(SDerrorDetail, "BSDF \"%s\": no optical layers",
206 >                                sd->name);
207 >                ezxml_free(fl);
208 >                return SDEformat;
209 >        }
210 >                                /* load geometry if present */
211 >        lastErr = SDloadGeometry(sd, wtl);
212 >        if (lastErr) {
213 >                ezxml_free(fl);
214 >                return lastErr;
215 >        }
216 >                                /* try loading variable resolution data */
217 >        lastErr = SDloadTre(sd, wtl);
218 >                                /* check our result */
219 >        if (lastErr == SDEsupport)      /* try matrix BSDF if not tree data */
220 >                lastErr = SDloadMtx(sd, wtl);
221 >                
222 >                                /* done with XML file */
223 >        ezxml_free(fl);
224 >        
225 >        if (lastErr) {          /* was there a load error? */
226 >                SDfreeBSDF(sd);
227 >                return lastErr;
228 >        }
229 >                                /* remove any insignificant components */
230 >        if (sd->rf != NULL && sd->rf->maxHemi <= .001) {
231 >                SDfreeSpectralDF(sd->rf); sd->rf = NULL;
232 >        }
233 >        if (sd->rb != NULL && sd->rb->maxHemi <= .001) {
234 >                SDfreeSpectralDF(sd->rb); sd->rb = NULL;
235 >        }
236 >        if (sd->tf != NULL && sd->tf->maxHemi <= .001) {
237 >                SDfreeSpectralDF(sd->tf); sd->tf = NULL;
238 >        }
239 >        if (sd->tb != NULL && sd->tb->maxHemi <= .001) {
240 >                SDfreeSpectralDF(sd->tb); sd->tb = NULL;
241 >        }
242 >                                /* return success */
243 >        return SDEnone;
244 > }
245 >
246 > /* Allocate new spectral distribution function */
247 > SDSpectralDF *
248 > SDnewSpectralDF(int nc)
249   {
250 <        ANGLE_BASIS  *ab = (ANGLE_BASIS *)p;
150 <        int     li;
151 <        double  theta, theta1;
250 >        SDSpectralDF    *df;
251          
252 <        if ((ndx < 0) | (ndx >= ab->nangles))
253 <                return(0);
254 <        for (li = 0; ndx >= ab->lat[li].nphis; li++)
156 <                ndx -= ab->lat[li].nphis;
157 <        theta1 = PI/180. * ab->lat[li+1].tmin;
158 <        if (ab->lat[li].nphis == 1) {           /* special case */
159 <                if (ab->lat[li].tmin > FTINY)
160 <                        error(USER, "unsupported BSDF coordinate system");
161 <                return(2.*PI*(1. - cos(theta1)));
252 >        if (nc <= 0) {
253 >                strcpy(SDerrorDetail, "Zero component spectral DF request");
254 >                return NULL;
255          }
256 <        theta = PI/180. * ab->lat[li].tmin;
257 <        return(2.*PI*(cos(theta) - cos(theta1))/(double)ab->lat[li].nphis);
256 >        df = (SDSpectralDF *)malloc(sizeof(SDSpectralDF) +
257 >                                        (nc-1)*sizeof(SDComponent));
258 >        if (df == NULL) {
259 >                sprintf(SDerrorDetail,
260 >                                "Cannot allocate %d component spectral DF", nc);
261 >                return NULL;
262 >        }
263 >        df->minProjSA = .0;
264 >        df->maxHemi = .0;
265 >        df->ncomp = nc;
266 >        memset(df->comp, 0, nc*sizeof(SDComponent));
267 >        return df;
268   }
269  
270 <
271 < static int
272 < ab_getvecR(             /* get reverse vector for this angle basis index */
170 <        FVECT v,
171 <        int ndx,
172 <        void *p
173 < )
270 > /* Add component(s) to spectral distribution function */
271 > SDSpectralDF *
272 > SDaddComponent(SDSpectralDF *odf, int nadd)
273   {
274 <        if (!ab_getvec(v, ndx, p))
176 <                return(0);
274 >        SDSpectralDF    *df;
275  
276 <        v[0] = -v[0];
277 <        v[1] = -v[1];
278 <        v[2] = -v[2];
279 <
280 <        return(1);
276 >        if (odf == NULL)
277 >                return SDnewSpectralDF(nadd);
278 >        if (nadd <= 0)
279 >                return odf;
280 >        df = (SDSpectralDF *)realloc(odf, sizeof(SDSpectralDF) +
281 >                                (odf->ncomp+nadd-1)*sizeof(SDComponent));
282 >        if (df == NULL) {
283 >                sprintf(SDerrorDetail,
284 >                        "Cannot add %d component(s) to spectral DF", nadd);
285 >                SDfreeSpectralDF(odf);
286 >                return NULL;
287 >        }
288 >        memset(df->comp+df->ncomp, 0, nadd*sizeof(SDComponent));
289 >        df->ncomp += nadd;
290 >        return df;
291   }
292  
293 <
294 < static int
295 < ab_getndxR(             /* get index corresponding to the reverse vector */
188 <        FVECT v,
189 <        void *p
190 < )
293 > /* Free cached cumulative distributions for BSDF component */
294 > void
295 > SDfreeCumulativeCache(SDSpectralDF *df)
296   {
297 <        FVECT  v2;
298 <        
194 <        v2[0] = -v[0];
195 <        v2[1] = -v[1];
196 <        v2[2] = -v[2];
297 >        int     n;
298 >        SDCDst  *cdp;
299  
300 <        return ab_getndx(v2, p);
300 >        if (df == NULL)
301 >                return;
302 >        for (n = df->ncomp; n-- > 0; )
303 >                while ((cdp = df->comp[n].cdList) != NULL) {
304 >                        df->comp[n].cdList = cdp->next;
305 >                        free(cdp);
306 >                }
307   }
308  
309 <
310 < static void
311 < load_angle_basis(       /* load BSDF angle basis */
204 <        ezxml_t wab
205 < )
309 > /* Free a spectral distribution function */
310 > void
311 > SDfreeSpectralDF(SDSpectralDF *df)
312   {
313 <        char    *abname = ezxml_txt(ezxml_child(wab, "AngleBasisName"));
314 <        ezxml_t wbb;
315 <        int     i;
210 <        
211 <        if (abname == NULL || !*abname)
313 >        int     n;
314 >
315 >        if (df == NULL)
316                  return;
317 <        for (i = nabases; i--; )
318 <                if (!strcmp(abname, abase_list[i].name))
319 <                        return;         /* XXX assume it's the same */
320 <        if (nabases >= MAXABASES)
321 <                error(INTERNAL, "too many angle bases");
218 <        strcpy(abase_list[nabases].name, abname);
219 <        abase_list[nabases].nangles = 0;
220 <        for (i = 0, wbb = ezxml_child(wab, "AngleBasisBlock");
221 <                        wbb != NULL; i++, wbb = wbb->next) {
222 <                if (i >= MAXLATS)
223 <                        error(INTERNAL, "too many latitudes in custom basis");
224 <                abase_list[nabases].lat[i+1].tmin = atof(ezxml_txt(
225 <                                ezxml_child(ezxml_child(wbb,
226 <                                        "ThetaBounds"), "UpperTheta")));
227 <                if (!i)
228 <                        abase_list[nabases].lat[i].tmin =
229 <                                        -abase_list[nabases].lat[i+1].tmin;
230 <                else if (!FEQ(atof(ezxml_txt(ezxml_child(ezxml_child(wbb,
231 <                                        "ThetaBounds"), "LowerTheta"))),
232 <                                abase_list[nabases].lat[i].tmin))
233 <                        error(WARNING, "theta values disagree in custom basis");
234 <                abase_list[nabases].nangles +=
235 <                        abase_list[nabases].lat[i].nphis =
236 <                                atoi(ezxml_txt(ezxml_child(wbb, "nPhis")));
237 <        }
238 <        abase_list[nabases++].lat[i].nphis = 0;
317 >        SDfreeCumulativeCache(df);
318 >        for (n = df->ncomp; n-- > 0; )
319 >                if (df->comp[n].dist != NULL)
320 >                        (*df->comp[n].func->freeSC)(df->comp[n].dist);
321 >        free(df);
322   }
323  
324 <
325 < static void
326 < load_bsdf_data(         /* load BSDF distribution for this wavelength */
244 <        struct BSDF_data *dp,
245 <        ezxml_t wdb
246 < )
324 > /* Shorten file path to useable BSDF name, removing suffix */
325 > void
326 > SDclipName(char *res, const char *fname)
327   {
328 <        char  *cbasis = ezxml_txt(ezxml_child(wdb,"ColumnAngleBasis"));
249 <        char  *rbasis = ezxml_txt(ezxml_child(wdb,"RowAngleBasis"));
250 <        char  *sdata;
251 <        int  i;
328 >        const char      *cp, *dot = NULL;
329          
330 <        if ((cbasis == NULL || !*cbasis) | (rbasis == NULL || !*rbasis)) {
331 <                error(WARNING, "missing column/row basis for BSDF");
330 >        for (cp = fname; *cp; cp++)
331 >                if (*cp == '.')
332 >                        dot = cp;
333 >        if ((dot == NULL) | (dot < fname+2))
334 >                dot = cp;
335 >        if (dot - fname >= SDnameLn)
336 >                fname = dot - SDnameLn + 1;
337 >        while (fname < dot)
338 >                *res++ = *fname++;
339 >        *res = '\0';
340 > }
341 >
342 > /* Initialize an unused BSDF struct (simply clears to zeroes) */
343 > void
344 > SDclearBSDF(SDData *sd, const char *fname)
345 > {
346 >        if (sd == NULL)
347                  return;
348 <        }
349 <        /* XXX need to add routines for loading in foreign bases */
258 <        for (i = nabases; i--; )
259 <                if (!strcmp(cbasis, abase_list[i].name)) {
260 <                        dp->ninc = abase_list[i].nangles;
261 <                        dp->ib_priv = (void *)&abase_list[i];
262 <                        dp->ib_vec = ab_getvecR;
263 <                        dp->ib_ndx = ab_getndxR;
264 <                        dp->ib_ohm = ab_getohm;
265 <                        break;
266 <                }
267 <        if (i < 0) {
268 <                sprintf(errmsg, "unsupported ColumnAngleBasis '%s'", cbasis);
269 <                error(WARNING, errmsg);
348 >        memset(sd, 0, sizeof(SDData));
349 >        if (fname == NULL)
350                  return;
351 <        }
352 <        for (i = nabases; i--; )
353 <                if (!strcmp(rbasis, abase_list[i].name)) {
354 <                        dp->nout = abase_list[i].nangles;
355 <                        dp->ob_priv = (void *)&abase_list[i];
356 <                        dp->ob_vec = ab_getvec;
357 <                        dp->ob_ndx = ab_getndx;
358 <                        dp->ob_ohm = ab_getohm;
279 <                        break;
280 <                }
281 <        if (i < 0) {
282 <                sprintf(errmsg, "unsupported RowAngleBasis '%s'", cbasis);
283 <                error(WARNING, errmsg);
351 >        SDclipName(sd->name, fname);
352 > }
353 >
354 > /* Free data associated with BSDF struct */
355 > void
356 > SDfreeBSDF(SDData *sd)
357 > {
358 >        if (sd == NULL)
359                  return;
360 +        if (sd->mgf != NULL) {
361 +                free(sd->mgf);
362 +                sd->mgf = NULL;
363          }
364 <                                /* read BSDF data */
365 <        sdata  = ezxml_txt(ezxml_child(wdb,"ScatteringData"));
366 <        if (sdata == NULL || !*sdata) {
289 <                error(WARNING, "missing BSDF ScatteringData");
290 <                return;
364 >        if (sd->rf != NULL) {
365 >                SDfreeSpectralDF(sd->rf);
366 >                sd->rf = NULL;
367          }
368 <        dp->bsdf = (float *)malloc(sizeof(float)*dp->ninc*dp->nout);
369 <        if (dp->bsdf == NULL)
370 <                error(SYSTEM, "out of memory in load_bsdf_data");
295 <        for (i = 0; i < dp->ninc*dp->nout; i++) {
296 <                char  *sdnext = fskip(sdata);
297 <                if (sdnext == NULL) {
298 <                        error(WARNING, "bad/missing BSDF ScatteringData");
299 <                        free(dp->bsdf); dp->bsdf = NULL;
300 <                        return;
301 <                }
302 <                while (*sdnext && isspace(*sdnext))
303 <                        sdnext++;
304 <                if (*sdnext == ',') sdnext++;
305 <                dp->bsdf[i] = atof(sdata);
306 <                sdata = sdnext;
368 >        if (sd->rb != NULL) {
369 >                SDfreeSpectralDF(sd->rb);
370 >                sd->rb = NULL;
371          }
372 <        while (isspace(*sdata))
373 <                sdata++;
374 <        if (*sdata) {
311 <                sprintf(errmsg, "%d extra characters after BSDF ScatteringData",
312 <                                strlen(sdata));
313 <                error(WARNING, errmsg);
372 >        if (sd->tf != NULL) {
373 >                SDfreeSpectralDF(sd->tf);
374 >                sd->tf = NULL;
375          }
376 +        if (sd->tb != NULL) {
377 +                SDfreeSpectralDF(sd->tb);
378 +                sd->tb = NULL;
379 +        }
380 +        sd->rLambFront.cieY = .0;
381 +        sd->rLambFront.spec.flags = 0;
382 +        sd->rLambBack.cieY = .0;
383 +        sd->rLambBack.spec.flags = 0;
384 +        sd->tLamb.cieY = .0;
385 +        sd->tLamb.spec.flags = 0;
386   }
387  
388 <
389 < static int
390 < check_bsdf_data(        /* check that BSDF data is sane */
320 <        struct BSDF_data *dp
321 < )
388 > /* Find writeable BSDF by name, or allocate new cache entry if absent */
389 > SDData *
390 > SDgetCache(const char *bname)
391   {
392 <        double          *omega_iarr, *omega_oarr;
393 <        double          dom, contrib, hemi_total;
325 <        int             nneg;
326 <        FVECT           v;
327 <        int             i, o;
392 >        struct SDCache_s        *sdl;
393 >        char                    sdnam[SDnameLn];
394  
395 <        if (dp == NULL || dp->bsdf == NULL)
396 <                return(0);
331 <        omega_iarr = (double *)calloc(dp->ninc, sizeof(double));
332 <        omega_oarr = (double *)calloc(dp->nout, sizeof(double));
333 <        if ((omega_iarr == NULL) | (omega_oarr == NULL))
334 <                error(SYSTEM, "out of memory in check_bsdf_data");
335 <                                        /* incoming projected solid angles */
336 <        hemi_total = .0;
337 <        for (i = dp->ninc; i--; ) {
338 <                dom = getBSDF_incohm(dp,i);
339 <                if (dom <= .0) {
340 <                        error(WARNING, "zero/negative incoming solid angle");
341 <                        continue;
342 <                }
343 <                if (!getBSDF_incvec(v,dp,i) || v[2] > FTINY) {
344 <                        error(WARNING, "illegal incoming BSDF direction");
345 <                        free(omega_iarr); free(omega_oarr);
346 <                        return(0);
347 <                }
348 <                hemi_total += omega_iarr[i] = dom * -v[2];
349 <        }
350 <        if ((hemi_total > 1.02*PI) | (hemi_total < 0.98*PI)) {
351 <                sprintf(errmsg, "incoming BSDF hemisphere off by %.1f%%",
352 <                                100.*(hemi_total/PI - 1.));
353 <                error(WARNING, errmsg);
354 <        }
355 <        dom = PI / hemi_total;          /* fix normalization */
356 <        for (i = dp->ninc; i--; )
357 <                omega_iarr[i] *= dom;
358 <                                        /* outgoing projected solid angles */
359 <        hemi_total = .0;
360 <        for (o = dp->nout; o--; ) {
361 <                dom = getBSDF_outohm(dp,o);
362 <                if (dom <= .0) {
363 <                        error(WARNING, "zero/negative outgoing solid angle");
364 <                        continue;
365 <                }
366 <                if (!getBSDF_outvec(v,dp,o) || v[2] < -FTINY) {
367 <                        error(WARNING, "illegal outgoing BSDF direction");
368 <                        free(omega_iarr); free(omega_oarr);
369 <                        return(0);
370 <                }
371 <                hemi_total += omega_oarr[o] = dom * v[2];
372 <        }
373 <        if ((hemi_total > 1.02*PI) | (hemi_total < 0.98*PI)) {
374 <                sprintf(errmsg, "outgoing BSDF hemisphere off by %.1f%%",
375 <                                100.*(hemi_total/PI - 1.));
376 <                error(WARNING, errmsg);
377 <        }
378 <        dom = PI / hemi_total;          /* fix normalization */
379 <        for (o = dp->nout; o--; )
380 <                omega_oarr[o] *= dom;
381 <        nneg = 0;                       /* check outgoing totals */
382 <        for (i = 0; i < dp->ninc; i++) {
383 <                hemi_total = .0;
384 <                for (o = dp->nout; o--; ) {
385 <                        double  f = BSDF_value(dp,i,o);
386 <                        if (f >= .0)
387 <                                hemi_total += f*omega_oarr[o];
388 <                        else {
389 <                                nneg += (f < -FTINY);
390 <                                BSDF_value(dp,i,o) = .0f;
391 <                        }
392 <                }
393 <                if (hemi_total > 1.02) {
394 <                        sprintf(errmsg,
395 <                        "incoming BSDF direction %d passes %.1f%% of light",
396 <                                        i, 100.*hemi_total);
397 <                        error(WARNING, errmsg);
398 <                }
399 <        }
400 <        if (nneg) {
401 <                sprintf(errmsg, "%d negative BSDF values (ignored)", nneg);
402 <                error(WARNING, errmsg);
403 <        }
404 <                                        /* reverse roles and check again */
405 <        for (o = 0; o < dp->nout; o++) {
406 <                hemi_total = .0;
407 <                for (i = dp->ninc; i--; )
408 <                        hemi_total += BSDF_value(dp,i,o) * omega_iarr[i];
395 >        if (bname == NULL)
396 >                return NULL;
397  
398 <                if (hemi_total > 1.02) {
399 <                        sprintf(errmsg,
400 <                        "outgoing BSDF direction %d collects %.1f%% of light",
401 <                                        o, 100.*hemi_total);
402 <                        error(WARNING, errmsg);
398 >        SDclipName(sdnam, bname);
399 >        for (sdl = SDcacheList; sdl != NULL; sdl = sdl->next)
400 >                if (!strcmp(sdl->bsdf.name, sdnam)) {
401 >                        sdl->refcnt++;
402 >                        return &sdl->bsdf;
403                  }
404 <        }
405 <        free(omega_iarr); free(omega_oarr);
406 <        return(1);
404 >
405 >        sdl = (struct SDCache_s *)calloc(1, sizeof(struct SDCache_s));
406 >        if (sdl == NULL)
407 >                return NULL;
408 >
409 >        strcpy(sdl->bsdf.name, sdnam);
410 >        sdl->next = SDcacheList;
411 >        SDcacheList = sdl;
412 >
413 >        sdl->refcnt = 1;
414 >        return &sdl->bsdf;
415   }
416  
417 < struct BSDF_data *
418 < load_BSDF(              /* load BSDF data from file */
419 <        char *fname
420 < )
417 > /* Get loaded BSDF from cache (or load and cache it on first call) */
418 > /* Report any problem to stderr and return NULL on failure */
419 > const SDData *
420 > SDcacheFile(const char *fname)
421   {
422 <        char                    *path;
423 <        ezxml_t                 fl, wtl, wld, wdb;
428 <        struct BSDF_data        *dp;
422 >        SDData          *sd;
423 >        SDError         ec;
424          
425 <        path = getpath(fname, getrlibpath(), R_OK);
426 <        if (path == NULL) {
427 <                sprintf(errmsg, "cannot find BSDF file \"%s\"", fname);
428 <                error(WARNING, errmsg);
429 <                return(NULL);
425 >        if (fname == NULL || !*fname)
426 >                return NULL;
427 >        SDerrorDetail[0] = '\0';
428 >        /* PLACE MUTEX LOCK HERE FOR THREAD-SAFE */
429 >        if ((sd = SDgetCache(fname)) == NULL) {
430 >                SDreportError(SDEmemory, stderr);
431 >                return NULL;
432          }
433 <        fl = ezxml_parse_file(path);
434 <        if (fl == NULL) {
435 <                sprintf(errmsg, "cannot open BSDF \"%s\"", path);
436 <                error(WARNING, errmsg);
440 <                return(NULL);
433 >        if (!SDisLoaded(sd) && (ec = SDloadFile(sd, fname))) {
434 >                SDreportError(ec, stderr);
435 >                SDfreeCache(sd);
436 >                sd = NULL;
437          }
438 <        if (ezxml_error(fl)[0]) {
439 <                sprintf(errmsg, "BSDF \"%s\" %s", path, ezxml_error(fl));
440 <                error(WARNING, errmsg);
441 <                ezxml_free(fl);
442 <                return(NULL);
438 >        /* END MUTEX LOCK */
439 >        return sd;
440 > }
441 >
442 > /* Free a BSDF from our cache (clear all if NULL) */
443 > void
444 > SDfreeCache(const SDData *sd)
445 > {
446 >        struct SDCache_s        *sdl, *sdLast = NULL;
447 >
448 >        if (sd == NULL) {               /* free entire list */
449 >                while ((sdl = SDcacheList) != NULL) {
450 >                        SDcacheList = sdl->next;
451 >                        SDfreeBSDF(&sdl->bsdf);
452 >                        free(sdl);
453 >                }
454 >                return;
455          }
456 <        if (strcmp(ezxml_name(fl), "WindowElement")) {
457 <                sprintf(errmsg,
458 <                        "BSDF \"%s\": top level node not 'WindowElement'",
459 <                                path);
460 <                error(WARNING, errmsg);
461 <                ezxml_free(fl);
462 <                return(NULL);
456 >        for (sdl = SDcacheList; sdl != NULL; sdl = (sdLast=sdl)->next)
457 >                if (&sdl->bsdf == sd)
458 >                        break;
459 >        if (sdl == NULL || (sdl->refcnt -= (sdl->refcnt > 0)))
460 >                return;                 /* missing or still in use */
461 >                                        /* keep unreferenced data? */
462 >        if (SDisLoaded(sd) && SDretainSet) {
463 >                if (SDretainSet == SDretainAll)
464 >                        return;         /* keep everything */
465 >                                        /* else free cumulative data */
466 >                SDfreeCumulativeCache(sd->rf);
467 >                SDfreeCumulativeCache(sd->rb);
468 >                SDfreeCumulativeCache(sd->tf);
469 >                SDfreeCumulativeCache(sd->tb);
470 >                return;
471          }
472 <        wtl = ezxml_child(ezxml_child(fl, "Optical"), "Layer");
473 <        load_angle_basis(ezxml_child(ezxml_child(wtl,
474 <                                "DataDefinition"), "AngleBasis"));
475 <        dp = (struct BSDF_data *)calloc(1, sizeof(struct BSDF_data));
476 <        for (wld = ezxml_child(wtl, "WavelengthData");
477 <                                wld != NULL; wld = wld->next) {
478 <                if (strcmp(ezxml_txt(ezxml_child(wld,"Wavelength")), "Visible"))
479 <                        continue;
480 <                wdb = ezxml_child(wld, "WavelengthDataBlock");
481 <                if (wdb == NULL) continue;
482 <                if (strcmp(ezxml_txt(ezxml_child(wdb,"WavelengthDataDirection")),
483 <                                        "Transmission Front"))
484 <                        continue;
485 <                load_bsdf_data(dp, wdb);        /* load front BTDF */
486 <                break;                          /* ignore the rest */
472 >                                        /* remove from list and free */
473 >        if (sdLast == NULL)
474 >                SDcacheList = sdl->next;
475 >        else
476 >                sdLast->next = sdl->next;
477 >        SDfreeBSDF(&sdl->bsdf);
478 >        free(sdl);
479 > }
480 >
481 > /* Sample an individual BSDF component */
482 > SDError
483 > SDsampComponent(SDValue *sv, FVECT ioVec, double randX, SDComponent *sdc)
484 > {
485 >        float           coef[SDmaxCh];
486 >        SDError         ec;
487 >        FVECT           inVec;
488 >        const SDCDst    *cd;
489 >        double          d;
490 >        int             n;
491 >                                        /* check arguments */
492 >        if ((sv == NULL) | (ioVec == NULL) | (sdc == NULL))
493 >                return SDEargument;
494 >                                        /* get cumulative distribution */
495 >        VCOPY(inVec, ioVec);
496 >        sv->cieY = 0;
497 >        cd = (*sdc->func->getCDist)(inVec, sdc);
498 >        if (cd != NULL)
499 >                sv->cieY = cd->cTotal;
500 >        if (sv->cieY <= 1e-6) {         /* nothing to sample? */
501 >                sv->spec = c_dfcolor;
502 >                memset(ioVec, 0, sizeof(FVECT));
503 >                return SDEnone;
504          }
505 <        ezxml_free(fl);                         /* done with XML file */
506 <        if (!check_bsdf_data(dp)) {
507 <                sprintf(errmsg, "bad/missing BTDF data in \"%s\"", path);
508 <                error(WARNING, errmsg);
509 <                free_BSDF(dp);
510 <                dp = NULL;
505 >                                        /* compute sample direction */
506 >        ec = (*sdc->func->sampCDist)(ioVec, randX, cd);
507 >        if (ec)
508 >                return ec;
509 >                                        /* get BSDF color */
510 >        n = (*sdc->func->getBSDFs)(coef, ioVec, inVec, sdc);
511 >        if (n <= 0) {
512 >                strcpy(SDerrorDetail, "BSDF sample value error");
513 >                return SDEinternal;
514          }
515 <        return(dp);
515 >        sv->spec = sdc->cspec[0];
516 >        d = coef[0];
517 >        while (--n) {
518 >                c_cmix(&sv->spec, d, &sv->spec, coef[n], &sdc->cspec[n]);
519 >                d += coef[n];
520 >        }
521 >        c_ccvt(&sv->spec, C_CSXY);      /* make sure (x,y) is set */
522 >        return SDEnone;
523   }
524  
525 + #define MS_MAXDIM       15
526  
527 + /* Convert 1-dimensional random variable to N-dimensional */
528   void
529 < free_BSDF(              /* free BSDF data structure */
485 <        struct BSDF_data *b
486 < )
529 > SDmultiSamp(double t[], int n, double randX)
530   {
531 <        if (b == NULL)
531 >        unsigned        nBits;
532 >        double          scale;
533 >        bitmask_t       ndx, coord[MS_MAXDIM];
534 >
535 >        if (n <= 0)                     /* check corner cases */
536                  return;
537 <        if (b->bsdf != NULL)
538 <                free(b->bsdf);
539 <        free(b);
537 >        if (randX < 0) randX = 0;
538 >        else if (randX >= 1.) randX = 0.999999999999999;
539 >        if (n == 1) {
540 >                t[0] = randX;
541 >                return;
542 >        }
543 >        while (n > MS_MAXDIM)           /* punt for higher dimensions */
544 >                t[--n] = rand()*(1./(RAND_MAX+.5));
545 >        nBits = (8*sizeof(bitmask_t) - 1) / n;
546 >        ndx = randX * (double)((bitmask_t)1 << (nBits*n));
547 >                                        /* get coordinate on Hilbert curve */
548 >        hilbert_i2c(n, nBits, ndx, coord);
549 >                                        /* convert back to [0,1) range */
550 >        scale = 1. / (double)((bitmask_t)1 << nBits);
551 >        while (n--)
552 >                t[n] = scale * ((double)coord[n] + rand()*(1./(RAND_MAX+.5)));
553   }
554  
555 + #undef MS_MAXDIM
556  
557 < int
558 < r_BSDF_incvec(          /* compute random input vector at given location */
559 <        FVECT v,
499 <        struct BSDF_data *b,
500 <        int i,
501 <        double rv,
502 <        MAT4 xm
503 < )
557 > /* Generate diffuse hemispherical sample */
558 > static void
559 > SDdiffuseSamp(FVECT outVec, int outFront, double randX)
560   {
561 <        FVECT   pert;
562 <        double  rad;
563 <        int     j;
564 <        
565 <        if (!getBSDF_incvec(v, b, i))
566 <                return(0);
567 <        rad = sqrt(getBSDF_incohm(b, i) / PI);
512 <        multisamp(pert, 3, rv);
513 <        for (j = 0; j < 3; j++)
514 <                v[j] += rad*(2.*pert[j] - 1.);
515 <        if (xm != NULL)
516 <                multv3(v, v, xm);
517 <        return(normalize(v) != 0.0);
561 >                                        /* convert to position on hemisphere */
562 >        SDmultiSamp(outVec, 2, randX);
563 >        SDsquare2disk(outVec, outVec[0], outVec[1]);
564 >        outVec[2] = 1. - outVec[0]*outVec[0] - outVec[1]*outVec[1];
565 >        outVec[2] = sqrt(outVec[2]*(outVec[2]>0));
566 >        if (!outFront)                  /* going out back? */
567 >                outVec[2] = -outVec[2];
568   }
569  
570 + /* Query projected solid angle coverage for non-diffuse BSDF direction */
571 + SDError
572 + SDsizeBSDF(double *projSA, const FVECT v1, const RREAL *v2,
573 +                                int qflags, const SDData *sd)
574 + {
575 +        SDSpectralDF    *rdf, *tdf;
576 +        SDError         ec;
577 +        int             i;
578 +                                        /* check arguments */
579 +        if ((projSA == NULL) | (v1 == NULL) | (sd == NULL))
580 +                return SDEargument;
581 +                                        /* initialize extrema */
582 +        switch (qflags) {
583 +        case SDqueryMax:
584 +                projSA[0] = .0;
585 +                break;
586 +        case SDqueryMin+SDqueryMax:
587 +                projSA[1] = .0;
588 +                /* fall through */
589 +        case SDqueryMin:
590 +                projSA[0] = 10.;
591 +                break;
592 +        case 0:
593 +                return SDEargument;
594 +        }
595 +        if (v1[2] > 0) {                /* front surface query? */
596 +                rdf = sd->rf;
597 +                tdf = (sd->tf != NULL) ? sd->tf : sd->tb;
598 +        } else {
599 +                rdf = sd->rb;
600 +                tdf = (sd->tb != NULL) ? sd->tb : sd->tf;
601 +        }
602 +        if (v2 != NULL) {               /* bidirectional? */
603 +                if (v1[2] > 0 ^ v2[2] > 0)
604 +                        rdf = NULL;
605 +                else
606 +                        tdf = NULL;
607 +        }
608 +        ec = SDEdata;                   /* run through components */
609 +        for (i = (rdf==NULL) ? 0 : rdf->ncomp; i--; ) {
610 +                ec = (*rdf->comp[i].func->queryProjSA)(projSA, v1, v2,
611 +                                                qflags, &rdf->comp[i]);
612 +                if (ec)
613 +                        return ec;
614 +        }
615 +        for (i = (tdf==NULL) ? 0 : tdf->ncomp; i--; ) {
616 +                ec = (*tdf->comp[i].func->queryProjSA)(projSA, v1, v2,
617 +                                                qflags, &tdf->comp[i]);
618 +                if (ec)
619 +                        return ec;
620 +        }
621 +        if (ec) {                       /* all diffuse? */
622 +                projSA[0] = M_PI;
623 +                if (qflags == SDqueryMin+SDqueryMax)
624 +                        projSA[1] = M_PI;
625 +        } else if (qflags == SDqueryMin+SDqueryMax && projSA[0] > projSA[1])
626 +                projSA[0] = projSA[1];
627 +        return SDEnone;
628 + }
629  
630 < int
631 < r_BSDF_outvec(          /* compute random output vector at given location */
632 <        FVECT v,
524 <        struct BSDF_data *b,
525 <        int o,
526 <        double rv,
527 <        MAT4 xm
528 < )
630 > /* Return BSDF for the given incident and scattered ray vectors */
631 > SDError
632 > SDevalBSDF(SDValue *sv, const FVECT outVec, const FVECT inVec, const SDData *sd)
633   {
634 <        FVECT   pert;
635 <        double  rad;
636 <        int     j;
637 <        
638 <        if (!getBSDF_outvec(v, b, o))
639 <                return(0);
640 <        rad = sqrt(getBSDF_outohm(b, o) / PI);
641 <        multisamp(pert, 3, rv);
642 <        for (j = 0; j < 3; j++)
643 <                v[j] += rad*(2.*pert[j] - 1.);
644 <        if (xm != NULL)
645 <                multv3(v, v, xm);
646 <        return(normalize(v) != 0.0);
634 >        int             inFront, outFront;
635 >        SDSpectralDF    *sdf;
636 >        float           coef[SDmaxCh];
637 >        int             nch, i;
638 >                                        /* check arguments */
639 >        if ((sv == NULL) | (outVec == NULL) | (inVec == NULL) | (sd == NULL))
640 >                return SDEargument;
641 >                                        /* whose side are we on? */
642 >        inFront = (inVec[2] > 0);
643 >        outFront = (outVec[2] > 0);
644 >                                        /* start with diffuse portion */
645 >        if (inFront & outFront) {
646 >                *sv = sd->rLambFront;
647 >                sdf = sd->rf;
648 >        } else if (!(inFront | outFront)) {
649 >                *sv = sd->rLambBack;
650 >                sdf = sd->rb;
651 >        } else if (inFront) {
652 >                *sv = sd->tLamb;
653 >                sdf = (sd->tf != NULL) ? sd->tf : sd->tb;
654 >        } else /* outFront & !inFront */ {
655 >                *sv = sd->tLamb;
656 >                sdf = (sd->tb != NULL) ? sd->tb : sd->tf;
657 >        }
658 >        sv->cieY *= 1./M_PI;
659 >                                        /* add non-diffuse components */
660 >        i = (sdf != NULL) ? sdf->ncomp : 0;
661 >        while (i-- > 0) {
662 >                nch = (*sdf->comp[i].func->getBSDFs)(coef, outVec, inVec,
663 >                                                        &sdf->comp[i]);
664 >                while (nch-- > 0) {
665 >                        c_cmix(&sv->spec, sv->cieY, &sv->spec,
666 >                                        coef[nch], &sdf->comp[i].cspec[nch]);
667 >                        sv->cieY += coef[nch];
668 >                }
669 >        }
670 >        c_ccvt(&sv->spec, C_CSXY);      /* make sure (x,y) is set */
671 >        return SDEnone;
672   }
673  
674 <
675 < static int
676 < addrot(                 /* compute rotation (x,y,z) => (xp,yp,zp) */
548 <        char *xfarg[],
549 <        FVECT xp,
550 <        FVECT yp,
551 <        FVECT zp
552 < )
674 > /* Compute directional hemispherical scattering at this incident angle */
675 > double
676 > SDdirectHemi(const FVECT inVec, int sflags, const SDData *sd)
677   {
678 <        static char     bufs[3][16];
679 <        int     bn = 0;
680 <        char    **xfp = xfarg;
681 <        double  theta;
678 >        double          hsum;
679 >        SDSpectralDF    *rdf, *tdf;
680 >        const SDCDst    *cd;
681 >        int             i;
682 >                                        /* check arguments */
683 >        if ((inVec == NULL) | (sd == NULL))
684 >                return .0;
685 >                                        /* gather diffuse components */
686 >        if (inVec[2] > 0) {
687 >                hsum = sd->rLambFront.cieY;
688 >                rdf = sd->rf;
689 >                tdf = (sd->tf != NULL) ? sd->tf : sd->tb;
690 >        } else /* !inFront */ {
691 >                hsum = sd->rLambBack.cieY;
692 >                rdf = sd->rb;
693 >                tdf = (sd->tb != NULL) ? sd->tb : sd->tf;
694 >        }
695 >        if ((sflags & SDsampDf+SDsampR) != SDsampDf+SDsampR)
696 >                hsum = .0;
697 >        if ((sflags & SDsampDf+SDsampT) == SDsampDf+SDsampT)
698 >                hsum += sd->tLamb.cieY;
699 >                                        /* gather non-diffuse components */
700 >        i = (((sflags & SDsampSp+SDsampR) == SDsampSp+SDsampR) &
701 >                        (rdf != NULL)) ? rdf->ncomp : 0;
702 >        while (i-- > 0) {               /* non-diffuse reflection */
703 >                cd = (*rdf->comp[i].func->getCDist)(inVec, &rdf->comp[i]);
704 >                if (cd != NULL)
705 >                        hsum += cd->cTotal;
706 >        }
707 >        i = (((sflags & SDsampSp+SDsampT) == SDsampSp+SDsampT) &
708 >                        (tdf != NULL)) ? tdf->ncomp : 0;
709 >        while (i-- > 0) {               /* non-diffuse transmission */
710 >                cd = (*tdf->comp[i].func->getCDist)(inVec, &tdf->comp[i]);
711 >                if (cd != NULL)
712 >                        hsum += cd->cTotal;
713 >        }
714 >        return hsum;
715 > }
716  
717 <        if (yp[2]*yp[2] + zp[2]*zp[2] < 2.*FTINY*FTINY) {
718 <                /* Special case for X' along Z-axis */
719 <                theta = -atan2(yp[0], yp[1]);
720 <                *xfp++ = "-ry";
721 <                *xfp++ = xp[2] < 0.0 ? "90" : "-90";
722 <                *xfp++ = "-rz";
723 <                sprintf(bufs[bn], "%f", theta*(180./PI));
724 <                *xfp++ = bufs[bn++];
725 <                return(xfp - xfarg);
717 > /* Sample BSDF direction based on the given random variable */
718 > SDError
719 > SDsampBSDF(SDValue *sv, FVECT ioVec, double randX, int sflags, const SDData *sd)
720 > {
721 >        SDError         ec;
722 >        FVECT           inVec;
723 >        int             inFront;
724 >        SDSpectralDF    *rdf, *tdf;
725 >        double          rdiff;
726 >        float           coef[SDmaxCh];
727 >        int             i, j, n, nr;
728 >        SDComponent     *sdc;
729 >        const SDCDst    **cdarr = NULL;
730 >                                        /* check arguments */
731 >        if ((sv == NULL) | (ioVec == NULL) | (sd == NULL) |
732 >                        (randX < 0) | (randX >= 1.))
733 >                return SDEargument;
734 >                                        /* whose side are we on? */
735 >        VCOPY(inVec, ioVec);
736 >        inFront = (inVec[2] > 0);
737 >                                        /* remember diffuse portions */
738 >        if (inFront) {
739 >                *sv = sd->rLambFront;
740 >                rdf = sd->rf;
741 >                tdf = (sd->tf != NULL) ? sd->tf : sd->tb;
742 >        } else /* !inFront */ {
743 >                *sv = sd->rLambBack;
744 >                rdf = sd->rb;
745 >                tdf = (sd->tb != NULL) ? sd->tb : sd->tf;
746          }
747 <        theta = atan2(yp[2], zp[2]);
748 <        if (!FEQ(theta,0.0)) {
749 <                *xfp++ = "-rx";
750 <                sprintf(bufs[bn], "%f", theta*(180./PI));
751 <                *xfp++ = bufs[bn++];
747 >        if ((sflags & SDsampDf+SDsampR) != SDsampDf+SDsampR)
748 >                sv->cieY = .0;
749 >        rdiff = sv->cieY;
750 >        if ((sflags & SDsampDf+SDsampT) == SDsampDf+SDsampT)
751 >                sv->cieY += sd->tLamb.cieY;
752 >                                        /* gather non-diffuse components */
753 >        i = nr = (((sflags & SDsampSp+SDsampR) == SDsampSp+SDsampR) &
754 >                        (rdf != NULL)) ? rdf->ncomp : 0;
755 >        j = (((sflags & SDsampSp+SDsampT) == SDsampSp+SDsampT) &
756 >                        (tdf != NULL)) ? tdf->ncomp : 0;
757 >        n = i + j;
758 >        if (n > 0 && (cdarr = (const SDCDst **)malloc(n*sizeof(SDCDst *))) == NULL)
759 >                return SDEmemory;
760 >        while (j-- > 0) {               /* non-diffuse transmission */
761 >                cdarr[i+j] = (*tdf->comp[j].func->getCDist)(inVec, &tdf->comp[j]);
762 >                if (cdarr[i+j] == NULL)
763 >                        cdarr[i+j] = &SDemptyCD;
764 >                sv->cieY += cdarr[i+j]->cTotal;
765          }
766 <        theta = asin(-xp[2]);
767 <        if (!FEQ(theta,0.0)) {
768 <                *xfp++ = "-ry";
769 <                sprintf(bufs[bn], " %f", theta*(180./PI));
770 <                *xfp++ = bufs[bn++];
766 >        while (i-- > 0) {               /* non-diffuse reflection */
767 >                cdarr[i] = (*rdf->comp[i].func->getCDist)(inVec, &rdf->comp[i]);
768 >                if (cdarr[i] == NULL)
769 >                        cdarr[i] = &SDemptyCD;
770 >                sv->cieY += cdarr[i]->cTotal;
771          }
772 <        theta = atan2(xp[1], xp[0]);
773 <        if (!FEQ(theta,0.0)) {
774 <                *xfp++ = "-rz";
775 <                sprintf(bufs[bn], "%f", theta*(180./PI));
585 <                *xfp++ = bufs[bn++];
772 >        if (sv->cieY <= 1e-6) {         /* anything to sample? */
773 >                sv->cieY = .0;
774 >                memset(ioVec, 0, sizeof(FVECT));
775 >                return SDEnone;
776          }
777 <        *xfp = NULL;
778 <        return(xfp - xfarg);
777 >                                        /* scale random variable */
778 >        randX *= sv->cieY;
779 >                                        /* diffuse reflection? */
780 >        if (randX < rdiff) {
781 >                SDdiffuseSamp(ioVec, inFront, randX/rdiff);
782 >                goto done;
783 >        }
784 >        randX -= rdiff;
785 >                                        /* diffuse transmission? */
786 >        if ((sflags & SDsampDf+SDsampT) == SDsampDf+SDsampT) {
787 >                if (randX < sd->tLamb.cieY) {
788 >                        sv->spec = sd->tLamb.spec;
789 >                        SDdiffuseSamp(ioVec, !inFront, randX/sd->tLamb.cieY);
790 >                        goto done;
791 >                }
792 >                randX -= sd->tLamb.cieY;
793 >        }
794 >                                        /* else one of cumulative dist. */
795 >        for (i = 0; i < n && randX >= cdarr[i]->cTotal; i++)
796 >                randX -= cdarr[i]->cTotal;
797 >        if (i >= n)
798 >                return SDEinternal;
799 >                                        /* compute sample direction */
800 >        sdc = (i < nr) ? &rdf->comp[i] : &tdf->comp[i-nr];
801 >        ec = (*sdc->func->sampCDist)(ioVec, randX/cdarr[i]->cTotal, cdarr[i]);
802 >        if (ec)
803 >                return ec;
804 >                                        /* compute color */
805 >        j = (*sdc->func->getBSDFs)(coef, ioVec, inVec, sdc);
806 >        if (j <= 0) {
807 >                sprintf(SDerrorDetail, "BSDF \"%s\" sampling value error",
808 >                                sd->name);
809 >                return SDEinternal;
810 >        }
811 >        sv->spec = sdc->cspec[0];
812 >        rdiff = coef[0];
813 >        while (--j) {
814 >                c_cmix(&sv->spec, rdiff, &sv->spec, coef[j], &sdc->cspec[j]);
815 >                rdiff += coef[j];
816 >        }
817 > done:
818 >        if (cdarr != NULL)
819 >                free(cdarr);
820 >        c_ccvt(&sv->spec, C_CSXY);      /* make sure (x,y) is set */
821 >        return SDEnone;
822   }
823  
824 + /* Compute World->BSDF transform from surface normal and up (Y) vector */
825 + SDError
826 + SDcompXform(RREAL vMtx[3][3], const FVECT sNrm, const FVECT uVec)
827 + {
828 +        if ((vMtx == NULL) | (sNrm == NULL) | (uVec == NULL))
829 +                return SDEargument;
830 +        VCOPY(vMtx[2], sNrm);
831 +        if (normalize(vMtx[2]) == 0)
832 +                return SDEargument;
833 +        fcross(vMtx[0], uVec, vMtx[2]);
834 +        if (normalize(vMtx[0]) == 0)
835 +                return SDEargument;
836 +        fcross(vMtx[1], vMtx[2], vMtx[0]);
837 +        return SDEnone;
838 + }
839  
840 < int
841 < getBSDF_xfm(            /* compute BSDF orient. -> world orient. transform */
842 <        MAT4 xm,
595 <        FVECT nrm,
596 <        UpDir ud
597 < )
840 > /* Compute inverse transform */
841 > SDError
842 > SDinvXform(RREAL iMtx[3][3], RREAL vMtx[3][3])
843   {
844 <        char    *xfargs[7];
845 <        XF      myxf;
601 <        FVECT   updir, xdest, ydest;
844 >        RREAL   mTmp[3][3];
845 >        double  d;
846  
847 <        updir[0] = updir[1] = updir[2] = 0.;
848 <        switch (ud) {
849 <        case UDzneg:
850 <                updir[2] = -1.;
851 <                break;
852 <        case UDyneg:
853 <                updir[1] = -1.;
854 <                break;
855 <        case UDxneg:
856 <                updir[0] = -1.;
613 <                break;
614 <        case UDxpos:
615 <                updir[0] = 1.;
616 <                break;
617 <        case UDypos:
618 <                updir[1] = 1.;
619 <                break;
620 <        case UDzpos:
621 <                updir[2] = 1.;
622 <                break;
623 <        case UDunknown:
624 <                return(0);
847 >        if ((iMtx == NULL) | (vMtx == NULL))
848 >                return SDEargument;
849 >                                        /* compute determinant */
850 >        mTmp[0][0] = vMtx[2][2]*vMtx[1][1] - vMtx[2][1]*vMtx[1][2];
851 >        mTmp[0][1] = vMtx[2][1]*vMtx[0][2] - vMtx[2][2]*vMtx[0][1];
852 >        mTmp[0][2] = vMtx[1][2]*vMtx[0][1] - vMtx[1][1]*vMtx[0][2];
853 >        d = vMtx[0][0]*mTmp[0][0] + vMtx[1][0]*mTmp[0][1] + vMtx[2][0]*mTmp[0][2];
854 >        if (d == 0) {
855 >                strcpy(SDerrorDetail, "Zero determinant in matrix inversion");
856 >                return SDEargument;
857          }
858 <        fcross(xdest, updir, nrm);
859 <        if (normalize(xdest) == 0.0)
860 <                return(0);
861 <        fcross(ydest, nrm, xdest);
862 <        xf(&myxf, addrot(xfargs, xdest, ydest, nrm), xfargs);
863 <        copymat4(xm, myxf.xfm);
864 <        return(1);
858 >        d = 1./d;                       /* invert matrix */
859 >        mTmp[0][0] *= d; mTmp[0][1] *= d; mTmp[0][2] *= d;
860 >        mTmp[1][0] = d*(vMtx[2][0]*vMtx[1][2] - vMtx[2][2]*vMtx[1][0]);
861 >        mTmp[1][1] = d*(vMtx[2][2]*vMtx[0][0] - vMtx[2][0]*vMtx[0][2]);
862 >        mTmp[1][2] = d*(vMtx[1][0]*vMtx[0][2] - vMtx[1][2]*vMtx[0][0]);
863 >        mTmp[2][0] = d*(vMtx[2][1]*vMtx[1][0] - vMtx[2][0]*vMtx[1][1]);
864 >        mTmp[2][1] = d*(vMtx[2][0]*vMtx[0][1] - vMtx[2][1]*vMtx[0][0]);
865 >        mTmp[2][2] = d*(vMtx[1][1]*vMtx[0][0] - vMtx[1][0]*vMtx[0][1]);
866 >        memcpy(iMtx, mTmp, sizeof(mTmp));
867 >        return SDEnone;
868 > }
869 >
870 > /* Transform and normalize direction (column) vector */
871 > SDError
872 > SDmapDir(FVECT resVec, RREAL vMtx[3][3], const FVECT inpVec)
873 > {
874 >        FVECT   vTmp;
875 >
876 >        if ((resVec == NULL) | (inpVec == NULL))
877 >                return SDEargument;
878 >        if (vMtx == NULL) {             /* assume they just want to normalize */
879 >                if (resVec != inpVec)
880 >                        VCOPY(resVec, inpVec);
881 >                return (normalize(resVec) > 0) ? SDEnone : SDEargument;
882 >        }
883 >        vTmp[0] = DOT(vMtx[0], inpVec);
884 >        vTmp[1] = DOT(vMtx[1], inpVec);
885 >        vTmp[2] = DOT(vMtx[2], inpVec);
886 >        if (normalize(vTmp) == 0)
887 >                return SDEargument;
888 >        VCOPY(resVec, vTmp);
889 >        return SDEnone;
890   }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines