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.2 by greg, Fri Jun 19 06:49:25 2009 UTC vs.
Revision 2.61 by greg, Tue Dec 7 23:49:50 2021 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       3               /* 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 + /* Retain BSDFs in cache list? */
51 + int                     SDretainSet = SDretainNone;
52  
53 < static int
54 < ab_getvec(              /* get vector for this angle basis index */
55 <        FVECT v,
56 <        int ndx,
57 <        void *p
58 < )
53 > /* Maximum cache size for any given BSDF? */
54 > unsigned long           SDmaxCache = 0;         /* 0 == unlimited */
55 >
56 > /* Report any error to the indicated stream */
57 > SDError
58 > SDreportError(SDError ec, FILE *fp)
59   {
60 <        ANGLE_BASIS  *ab = (ANGLE_BASIS *)p;
61 <        int     li;
62 <        double  pol, azi, d;
63 <        
64 <        if ((ndx < 0) | (ndx >= ab->nangles))
65 <                return(0);
66 <        for (li = 0; ndx >= ab->lat[li].nphis; li++)
67 <                ndx -= ab->lat[li].nphis;
68 <        pol = PI/180.*0.5*(ab->lat[li].tmin + ab->lat[li+1].tmin);
69 <        azi = 2.*PI*ndx/ab->lat[li].nphis;
70 <        v[2] = d = cos(pol);
71 <        d = sqrt(1. - d*d);     /* sin(pol) */
72 <        v[0] = cos(azi)*d;
73 <        v[1] = sin(azi)*d;
74 <        return(1);
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   }
78  
79 <
80 < static int
81 < ab_getndx(              /* get index corresponding to the given vector */
92 <        FVECT v,
93 <        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, ndx;
86 <        double  pol, azi, d;
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 <        if ((v[2] < -1.0) | (v[2] > 1.0))
95 <                return(-1);
96 <        pol = 180.0/PI*acos(v[2]);
97 <        azi = 180.0/PI*atan2(v[1], v[0]);
98 <        if (azi < 0.0) azi += 360.0;
99 <        for (li = 1; ab->lat[li].tmin <= pol; li++)
100 <                if (!ab->lat[li].nphis)
101 <                        return(-1);
102 <        --li;
103 <        ndx = (int)((1./360.)*azi*ab->lat[li].nphis + 0.5);
104 <        if (ndx >= ab->lat[li].nphis) ndx = 0;
105 <        while (li--)
106 <                ndx += ab->lat[li].nphis;
107 <        return(ndx);
94 > /* Load geometric dimensions and description (if any) */
95 > static SDError
96 > SDloadGeometry(SDData *sd, ezxml_t wtl)
97 > {
98 >        ezxml_t         node, matl, geom;
99 >        double          cfact;
100 >        const char      *fmt = NULL, *mgfstr;
101 >
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;
124 <        int     li;
125 <        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++)
130 <                ndx -= ab->lat[li].nphis;
131 <        theta1 = PI/180. * ab->lat[li+1].tmin;
132 <        if (ab->lat[li].nphis == 1) {           /* special case */
133 <                if (ab->lat[li].tmin > FTINY)
134 <                        error(USER, "unsupported BSDF coordinate system");
135 <                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 */
144 <        FVECT v,
145 <        int ndx,
146 <        void *p
147 < )
270 > /* Add component(s) to spectral distribution function */
271 > SDSpectralDF *
272 > SDaddComponent(SDSpectralDF *odf, int nadd)
273   {
274 <        if (!ab_getvec(v, ndx, p))
150 <                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 */
162 <        FVECT v,
163 <        void *p
164 < )
293 > /* Free cached cumulative distributions for BSDF component */
294 > void
295 > SDfreeCumulativeCache(SDSpectralDF *df)
296   {
297 <        FVECT  v2;
298 <        
168 <        v2[0] = -v[0];
169 <        v2[1] = -v[1];
170 <        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 + /* Free a spectral distribution function */
310 + void
311 + SDfreeSpectralDF(SDSpectralDF *df)
312 + {
313 +        int     n;
314  
315 < static void
316 < load_bsdf_data(         /* load BSDF distribution for this wavelength */
317 <        struct BSDF_data *dp,
318 <        ezxml_t wdb
319 < )
315 >        if (df == NULL)
316 >                return;
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 > /* 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"));
183 <        char  *rbasis = ezxml_txt(ezxml_child(wdb,"RowAngleBasis"));
184 <        char  *sdata;
185 <        int  i;
328 >        const char      *cp, *dot = NULL;
329          
330 <        if ((cbasis == NULL) | (rbasis == NULL)) {
331 <                error(WARNING, "missing column/row basis for BSDF");
330 >        for (cp = fname; *cp; cp++)
331 >                if (*cp == '.')
332 >                        dot = cp;
333 >                else if (*cp == '/')
334 >                        dot = NULL;
335 >        if ((dot == NULL) | (dot < fname+2))
336 >                dot = cp;
337 >        if (dot - fname >= SDnameLn)
338 >                fname = dot - SDnameLn + 1;
339 >        while (fname < dot)
340 >                *res++ = *fname++;
341 >        *res = '\0';
342 > }
343 >
344 > /* Initialize an unused BSDF struct (simply clears to zeroes) */
345 > void
346 > SDclearBSDF(SDData *sd, const char *fname)
347 > {
348 >        if (sd == NULL)
349                  return;
350 <        }
351 <        /* XXX need to add routines for loading in foreign bases */
192 <        for (i = nabases; i--; )
193 <                if (!strcmp(cbasis, abase_list[i].name)) {
194 <                        dp->ninc = abase_list[i].nangles;
195 <                        dp->ib_priv = (void *)&abase_list[i];
196 <                        dp->ib_vec = ab_getvecR;
197 <                        dp->ib_ndx = ab_getndxR;
198 <                        dp->ib_ohm = ab_getohm;
199 <                        break;
200 <                }
201 <        if (i < 0) {
202 <                sprintf(errmsg, "unsupported ColumnAngleBasis '%s'", cbasis);
203 <                error(WARNING, errmsg);
350 >        memset(sd, 0, sizeof(SDData));
351 >        if (fname == NULL)
352                  return;
353 <        }
354 <        for (i = nabases; i--; )
355 <                if (!strcmp(rbasis, abase_list[i].name)) {
356 <                        dp->nout = abase_list[i].nangles;
357 <                        dp->ob_priv = (void *)&abase_list[i];
358 <                        dp->ob_vec = ab_getvec;
359 <                        dp->ob_ndx = ab_getndx;
360 <                        dp->ob_ohm = ab_getohm;
213 <                        break;
214 <                }
215 <        if (i < 0) {
216 <                sprintf(errmsg, "unsupported RowAngleBasis '%s'", cbasis);
217 <                error(WARNING, errmsg);
353 >        SDclipName(sd->name, fname);
354 > }
355 >
356 > /* Free data associated with BSDF struct */
357 > void
358 > SDfreeBSDF(SDData *sd)
359 > {
360 >        if (sd == NULL)
361                  return;
362 +        if (sd->mgf != NULL) {
363 +                free(sd->mgf);
364 +                sd->mgf = NULL;
365          }
366 <                                /* read BSDF data */
367 <        sdata  = ezxml_txt(ezxml_child(wdb,"ScatteringData"));
368 <        if (sdata == NULL) {
223 <                error(WARNING, "missing BSDF ScatteringData");
224 <                return;
366 >        if (sd->rf != NULL) {
367 >                SDfreeSpectralDF(sd->rf);
368 >                sd->rf = NULL;
369          }
370 <        dp->bsdf = (float *)malloc(sizeof(float)*dp->ninc*dp->nout);
371 <        if (dp->bsdf == NULL)
372 <                error(SYSTEM, "out of memory in load_bsdf_data");
229 <        for (i = 0; i < dp->ninc*dp->nout; i++) {
230 <                char  *sdnext = fskip(sdata);
231 <                if (sdnext == NULL) {
232 <                        error(WARNING, "bad/missing BSDF ScatteringData");
233 <                        free(dp->bsdf); dp->bsdf = NULL;
234 <                        return;
235 <                }
236 <                while (*sdnext && isspace(*sdnext))
237 <                        sdnext++;
238 <                if (*sdnext == ',') sdnext++;
239 <                dp->bsdf[i] = atof(sdata);
240 <                sdata = sdnext;
370 >        if (sd->rb != NULL) {
371 >                SDfreeSpectralDF(sd->rb);
372 >                sd->rb = NULL;
373          }
374 <        while (isspace(*sdata))
375 <                sdata++;
376 <        if (*sdata) {
245 <                sprintf(errmsg, "%d extra characters after BSDF ScatteringData",
246 <                                strlen(sdata));
247 <                error(WARNING, errmsg);
374 >        if (sd->tf != NULL) {
375 >                SDfreeSpectralDF(sd->tf);
376 >                sd->tf = NULL;
377          }
378 +        if (sd->tb != NULL) {
379 +                SDfreeSpectralDF(sd->tb);
380 +                sd->tb = NULL;
381 +        }
382 +        sd->rLambFront.cieY = .0;
383 +        sd->rLambFront.spec.flags = 0;
384 +        sd->rLambBack.cieY = .0;
385 +        sd->rLambBack.spec.flags = 0;
386 +        sd->tLambFront.cieY = .0;
387 +        sd->tLambFront.spec.flags = 0;
388 +        sd->tLambBack.cieY = .0;
389 +        sd->tLambBack.spec.flags = 0;
390   }
391  
392 <
393 < static int
394 < check_bsdf_data(        /* check that BSDF data is sane */
254 <        struct BSDF_data *dp
255 < )
392 > /* Find writeable BSDF by name, or allocate new cache entry if absent */
393 > SDData *
394 > SDgetCache(const char *bname)
395   {
396 <        double          *omega_iarr, *omega_oarr;
397 <        double          dom, contrib, hemi_total;
259 <        int             nneg;
260 <        FVECT           v;
261 <        int             i, o;
396 >        struct SDCache_s        *sdl;
397 >        char                    sdnam[SDnameLn];
398  
399 <        if (dp == NULL || dp->bsdf == NULL)
400 <                return(0);
265 <        omega_iarr = (double *)calloc(dp->ninc, sizeof(double));
266 <        omega_oarr = (double *)calloc(dp->nout, sizeof(double));
267 <        if ((omega_iarr == NULL) | (omega_oarr == NULL))
268 <                error(SYSTEM, "out of memory in check_bsdf_data");
269 <                                        /* incoming projected solid angles */
270 <        hemi_total = .0;
271 <        for (i = dp->ninc; i--; ) {
272 <                dom = getBSDF_incohm(dp,i);
273 <                if (dom <= .0) {
274 <                        error(WARNING, "zero/negative incoming solid angle");
275 <                        continue;
276 <                }
277 <                if (!getBSDF_incvec(v,dp,i) || v[2] > FTINY) {
278 <                        error(WARNING, "illegal incoming BSDF direction");
279 <                        free(omega_iarr); free(omega_oarr);
280 <                        return(0);
281 <                }
282 <                hemi_total += omega_iarr[i] = dom * -v[2];
283 <        }
284 <        if ((hemi_total > 1.02*PI) | (hemi_total < 0.98*PI)) {
285 <                sprintf(errmsg, "incoming BSDF hemisphere off by %.1f%%",
286 <                                100.*(hemi_total/PI - 1.));
287 <                error(WARNING, errmsg);
288 <        }
289 <        dom = PI / hemi_total;          /* fix normalization */
290 <        for (i = dp->ninc; i--; )
291 <                omega_iarr[i] *= dom;
292 <                                        /* outgoing projected solid angles */
293 <        hemi_total = .0;
294 <        for (o = dp->nout; o--; ) {
295 <                dom = getBSDF_outohm(dp,o);
296 <                if (dom <= .0) {
297 <                        error(WARNING, "zero/negative outgoing solid angle");
298 <                        continue;
299 <                }
300 <                if (!getBSDF_outvec(v,dp,o) || v[2] < -FTINY) {
301 <                        error(WARNING, "illegal outgoing BSDF direction");
302 <                        free(omega_iarr); free(omega_oarr);
303 <                        return(0);
304 <                }
305 <                hemi_total += omega_oarr[o] = dom * v[2];
306 <        }
307 <        if ((hemi_total > 1.02*PI) | (hemi_total < 0.98*PI)) {
308 <                sprintf(errmsg, "outgoing BSDF hemisphere off by %.1f%%",
309 <                                100.*(hemi_total/PI - 1.));
310 <                error(WARNING, errmsg);
311 <        }
312 <        dom = PI / hemi_total;          /* fix normalization */
313 <        for (o = dp->nout; o--; )
314 <                omega_oarr[o] *= dom;
315 <        nneg = 0;                       /* check outgoing totals */
316 <        for (i = 0; i < dp->ninc; i++) {
317 <                hemi_total = .0;
318 <                for (o = dp->nout; o--; ) {
319 <                        double  f = BSDF_value(dp,i,o);
320 <                        if (f >= .0)
321 <                                hemi_total += f*omega_oarr[o];
322 <                        else {
323 <                                nneg += (f < -FTINY);
324 <                                BSDF_value(dp,i,o) = .0f;
325 <                        }
326 <                }
327 <                if (hemi_total > 1.02) {
328 <                        sprintf(errmsg,
329 <                        "incoming BSDF direction %d passes %.1f%% of light",
330 <                                        i, 100.*hemi_total);
331 <                        error(WARNING, errmsg);
332 <                }
333 <        }
334 <        if (nneg) {
335 <                sprintf(errmsg, "%d negative BSDF values (ignored)", nneg);
336 <                error(WARNING, errmsg);
337 <        }
338 <                                        /* reverse roles and check again */
339 <        for (o = 0; o < dp->nout; o++) {
340 <                hemi_total = .0;
341 <                for (i = dp->ninc; i--; )
342 <                        hemi_total += BSDF_value(dp,i,o) * omega_iarr[i];
399 >        if (bname == NULL)
400 >                return NULL;
401  
402 <                if (hemi_total > 1.02) {
403 <                        sprintf(errmsg,
404 <                        "outgoing BSDF direction %d collects %.1f%% of light",
405 <                                        o, 100.*hemi_total);
406 <                        error(WARNING, errmsg);
402 >        SDclipName(sdnam, bname);
403 >        for (sdl = SDcacheList; sdl != NULL; sdl = sdl->next)
404 >                if (!strcmp(sdl->bsdf.name, sdnam)) {
405 >                        sdl->refcnt++;
406 >                        return &sdl->bsdf;
407                  }
408 <        }
409 <        free(omega_iarr); free(omega_oarr);
410 <        return(1);
408 >
409 >        sdl = (struct SDCache_s *)calloc(1, sizeof(struct SDCache_s));
410 >        if (sdl == NULL)
411 >                return NULL;
412 >
413 >        strcpy(sdl->bsdf.name, sdnam);
414 >        sdl->next = SDcacheList;
415 >        SDcacheList = sdl;
416 >
417 >        sdl->refcnt = 1;
418 >        return &sdl->bsdf;
419   }
420  
421 < struct BSDF_data *
422 < load_BSDF(              /* load BSDF data from file */
423 <        char *fname
424 < )
421 > /* Get loaded BSDF from cache (or load and cache it on first call) */
422 > /* Report any problem to stderr and return NULL on failure */
423 > const SDData *
424 > SDcacheFile(const char *fname)
425   {
426 <        char                    *path;
427 <        ezxml_t                 fl, wtl, wld, wdb;
362 <        struct BSDF_data        *dp;
426 >        SDData          *sd;
427 >        SDError         ec;
428          
429 <        path = getpath(fname, getrlibpath(), R_OK);
430 <        if (path == NULL) {
431 <                sprintf(errmsg, "cannot find BSDF file \"%s\"", fname);
432 <                error(WARNING, errmsg);
433 <                return(NULL);
429 >        if (fname == NULL || !*fname)
430 >                return NULL;
431 >        SDerrorDetail[0] = '\0';
432 >        /* PLACE MUTEX LOCK HERE FOR THREAD-SAFE */
433 >        if ((sd = SDgetCache(fname)) == NULL) {
434 >                SDreportError(SDEmemory, stderr);
435 >                return NULL;
436          }
437 <        fl = ezxml_parse_file(path);
438 <        if (fl == NULL) {
439 <                sprintf(errmsg, "cannot open BSDF \"%s\"", path);
440 <                error(WARNING, errmsg);
374 <                return(NULL);
437 >        if (!SDisLoaded(sd) && (ec = SDloadFile(sd, fname))) {
438 >                SDreportError(ec, stderr);
439 >                SDfreeCache(sd);
440 >                sd = NULL;
441          }
442 <        if (ezxml_error(fl)[0]) {
443 <                sprintf(errmsg, "BSDF \"%s\" %s", path, ezxml_error(fl));
444 <                error(WARNING, errmsg);
445 <                ezxml_free(fl);
446 <                return(NULL);
442 >        /* END MUTEX LOCK */
443 >        return sd;
444 > }
445 >
446 > /* Free a BSDF from our cache (clear all if NULL) */
447 > void
448 > SDfreeCache(const SDData *sd)
449 > {
450 >        struct SDCache_s        *sdl, *sdLast = NULL;
451 >
452 >        if (sd == NULL) {               /* free entire list */
453 >                while ((sdl = SDcacheList) != NULL) {
454 >                        SDcacheList = sdl->next;
455 >                        SDfreeBSDF(&sdl->bsdf);
456 >                        free(sdl);
457 >                }
458 >                return;
459          }
460 <        if (strcmp(ezxml_name(fl), "WindowElement")) {
461 <                sprintf(errmsg,
462 <                        "BSDF \"%s\": top level node not 'WindowElement'",
463 <                                path);
464 <                error(WARNING, errmsg);
465 <                ezxml_free(fl);
466 <                return(NULL);
460 >        for (sdl = SDcacheList; sdl != NULL; sdl = (sdLast=sdl)->next)
461 >                if (&sdl->bsdf == sd)
462 >                        break;
463 >        if (sdl == NULL || (sdl->refcnt -= (sdl->refcnt > 0)))
464 >                return;                 /* missing or still in use */
465 >                                        /* keep unreferenced data? */
466 >        if (SDisLoaded(sd) && SDretainSet) {
467 >                if (SDretainSet == SDretainAll)
468 >                        return;         /* keep everything */
469 >                                        /* else free cumulative data */
470 >                SDfreeCumulativeCache(sd->rf);
471 >                SDfreeCumulativeCache(sd->rb);
472 >                SDfreeCumulativeCache(sd->tf);
473 >                SDfreeCumulativeCache(sd->tb);
474 >                return;
475          }
476 <        wtl = ezxml_child(ezxml_child(fl, "Optical"), "Layer");
477 <        dp = (struct BSDF_data *)calloc(1, sizeof(struct BSDF_data));
478 <        for (wld = ezxml_child(wtl, "WavelengthData");
479 <                                wld != NULL; wld = wld->next) {
480 <                if (strcmp(ezxml_txt(ezxml_child(wld,"Wavelength")), "Visible"))
481 <                        continue;
482 <                wdb = ezxml_child(wld, "WavelengthDataBlock");
483 <                if (wdb == NULL) continue;
484 <                if (strcmp(ezxml_txt(ezxml_child(wdb,"WavelengthDataDirection")),
485 <                                        "Transmission Front"))
486 <                        continue;
487 <                load_bsdf_data(dp, wdb);        /* load front BTDF */
488 <                break;                          /* ignore the rest */
476 >                                        /* remove from list and free */
477 >        if (sdLast == NULL)
478 >                SDcacheList = sdl->next;
479 >        else
480 >                sdLast->next = sdl->next;
481 >        SDfreeBSDF(&sdl->bsdf);
482 >        free(sdl);
483 > }
484 >
485 > /* Sample an individual BSDF component */
486 > SDError
487 > SDsampComponent(SDValue *sv, FVECT ioVec, double randX, SDComponent *sdc)
488 > {
489 >        float           coef[SDmaxCh];
490 >        SDError         ec;
491 >        FVECT           inVec;
492 >        const SDCDst    *cd;
493 >        double          d;
494 >        int             n;
495 >                                        /* check arguments */
496 >        if ((sv == NULL) | (ioVec == NULL) | (sdc == NULL))
497 >                return SDEargument;
498 >                                        /* get cumulative distribution */
499 >        VCOPY(inVec, ioVec);
500 >        sv->cieY = 0;
501 >        cd = (*sdc->func->getCDist)(inVec, sdc);
502 >        if (cd != NULL)
503 >                sv->cieY = cd->cTotal;
504 >        if (sv->cieY <= 1e-6) {         /* nothing to sample? */
505 >                sv->spec = c_dfcolor;
506 >                memset(ioVec, 0, sizeof(FVECT));
507 >                return SDEnone;
508          }
509 <        ezxml_free(fl);                         /* done with XML file */
510 <        if (!check_bsdf_data(dp)) {
511 <                sprintf(errmsg, "bad/missing BTDF data in \"%s\"", path);
512 <                error(WARNING, errmsg);
513 <                free_BSDF(dp);
514 <                dp = NULL;
509 >                                        /* compute sample direction */
510 >        ec = (*sdc->func->sampCDist)(ioVec, randX, cd);
511 >        if (ec)
512 >                return ec;
513 >                                        /* get BSDF color */
514 >        n = (*sdc->func->getBSDFs)(coef, inVec, ioVec, sdc);
515 >        if (n <= 0) {
516 >                strcpy(SDerrorDetail, "BSDF sample value error");
517 >                return SDEinternal;
518          }
519 <        return(dp);
519 >        sv->spec = sdc->cspec[0];
520 >        d = coef[0];
521 >        while (--n) {
522 >                c_cmix(&sv->spec, d, &sv->spec, coef[n], &sdc->cspec[n]);
523 >                d += coef[n];
524 >        }
525 >        c_ccvt(&sv->spec, C_CSXY);      /* make sure (x,y) is set */
526 >        return SDEnone;
527   }
528  
529 + #define MS_MAXDIM       15
530  
531 + /* Convert 1-dimensional random variable to N-dimensional */
532   void
533 < free_BSDF(              /* free BSDF data structure */
417 <        struct BSDF_data *b
418 < )
533 > SDmultiSamp(double t[], int n, double randX)
534   {
535 <        if (b == NULL)
535 >        unsigned        nBits;
536 >        double          scale;
537 >        bitmask_t       ndx, coord[MS_MAXDIM];
538 >
539 >        if (n <= 0)                     /* check corner cases */
540                  return;
541 <        if (b->bsdf != NULL)
542 <                free(b->bsdf);
543 <        free(b);
541 >        if (randX < 0) randX = 0;
542 >        else if (randX >= 1.) randX = 0.999999999999999;
543 >        if (n == 1) {
544 >                t[0] = randX;
545 >                return;
546 >        }
547 >        while (n > MS_MAXDIM)           /* punt for higher dimensions */
548 >                t[--n] = rand()*(1./(RAND_MAX+.5));
549 >        nBits = (8*sizeof(bitmask_t) - 1) / n;
550 >        ndx = randX * (double)((bitmask_t)1 << (nBits*n));
551 >                                        /* get coordinate on Hilbert curve */
552 >        hilbert_i2c(n, nBits, ndx, coord);
553 >                                        /* convert back to [0,1) range */
554 >        scale = 1. / (double)((bitmask_t)1 << nBits);
555 >        while (n--)
556 >                t[n] = scale * ((double)coord[n] + rand()*(1./(RAND_MAX+.5)));
557   }
558  
559 + #undef MS_MAXDIM
560  
561 < int
562 < r_BSDF_incvec(          /* compute random input vector at given location */
563 <        FVECT v,
431 <        struct BSDF_data *b,
432 <        int i,
433 <        double rv,
434 <        MAT4 xm
435 < )
561 > /* Generate diffuse hemispherical sample */
562 > static void
563 > SDdiffuseSamp(FVECT ioVec, int outFront, double randX)
564   {
565 <        FVECT   pert;
566 <        double  rad;
567 <        int     j;
568 <        
569 <        if (!getBSDF_incvec(v, b, i))
570 <                return(0);
571 <        rad = sqrt(getBSDF_incohm(b, i) / PI);
444 <        multisamp(pert, 3, rv);
445 <        for (j = 0; j < 3; j++)
446 <                v[j] += rad*(2.*pert[j] - 1.);
447 <        if (xm != NULL)
448 <                multv3(v, v, xm);
449 <        return(normalize(v) != 0.0);
565 >                                        /* convert to position on hemisphere */
566 >        SDmultiSamp(ioVec, 2, randX);
567 >        SDsquare2disk(ioVec, ioVec[0], ioVec[1]);
568 >        ioVec[2] = 1. - ioVec[0]*ioVec[0] - ioVec[1]*ioVec[1];
569 >        ioVec[2] = sqrt(ioVec[2]*(ioVec[2]>0));
570 >        if (!outFront)                  /* going out back? */
571 >                ioVec[2] = -ioVec[2];
572   }
573  
574 + /* Query projected solid angle coverage for non-diffuse BSDF direction */
575 + SDError
576 + SDsizeBSDF(double *projSA, const FVECT v1, const RREAL *v2,
577 +                                int qflags, const SDData *sd)
578 + {
579 +        SDSpectralDF    *rdf, *tdf;
580 +        SDError         ec;
581 +        int             i;
582 +                                        /* check arguments */
583 +        if ((projSA == NULL) | (v1 == NULL) | (sd == NULL))
584 +                return SDEargument;
585 +                                        /* initialize extrema */
586 +        switch (qflags) {
587 +        case SDqueryMax:
588 +                projSA[0] = .0;
589 +                break;
590 +        case SDqueryMin+SDqueryMax:
591 +                projSA[1] = .0;
592 +                /* fall through */
593 +        case SDqueryMin:
594 +                projSA[0] = 10.;
595 +                break;
596 +        case 0:
597 +                return SDEargument;
598 +        }
599 +        if (v1[2] > 0) {                /* front surface query? */
600 +                rdf = sd->rf;
601 +                tdf = (sd->tf != NULL) ? sd->tf : sd->tb;
602 +        } else {
603 +                rdf = sd->rb;
604 +                tdf = (sd->tb != NULL) ? sd->tb : sd->tf;
605 +        }
606 +        if (v2 != NULL) {               /* bidirectional? */
607 +                if (v1[2] > 0 ^ v2[2] > 0)
608 +                        rdf = NULL;
609 +                else
610 +                        tdf = NULL;
611 +        }
612 +        ec = SDEdata;                   /* run through components */
613 +        for (i = (rdf==NULL) ? 0 : rdf->ncomp; i--; ) {
614 +                ec = (*rdf->comp[i].func->queryProjSA)(projSA, v1, v2,
615 +                                                qflags, &rdf->comp[i]);
616 +                if (ec)
617 +                        return ec;
618 +        }
619 +        for (i = (tdf==NULL) ? 0 : tdf->ncomp; i--; ) {
620 +                ec = (*tdf->comp[i].func->queryProjSA)(projSA, v1, v2,
621 +                                                qflags, &tdf->comp[i]);
622 +                if (ec)
623 +                        return ec;
624 +        }
625 +        if (ec) {                       /* all diffuse? */
626 +                projSA[0] = M_PI;
627 +                if (qflags == SDqueryMin+SDqueryMax)
628 +                        projSA[1] = M_PI;
629 +        } else if (qflags == SDqueryMin+SDqueryMax && projSA[0] > projSA[1])
630 +                projSA[0] = projSA[1];
631 +        return SDEnone;
632 + }
633  
634 < int
635 < r_BSDF_outvec(          /* compute random output vector at given location */
636 <        FVECT v,
456 <        struct BSDF_data *b,
457 <        int o,
458 <        double rv,
459 <        MAT4 xm
460 < )
634 > /* Return BSDF for the given incident and scattered ray vectors */
635 > SDError
636 > SDevalBSDF(SDValue *sv, const FVECT inVec, const FVECT outVec, const SDData *sd)
637   {
638 <        FVECT   pert;
639 <        double  rad;
640 <        int     j;
641 <        
642 <        if (!getBSDF_outvec(v, b, o))
643 <                return(0);
644 <        rad = sqrt(getBSDF_outohm(b, o) / PI);
645 <        multisamp(pert, 3, rv);
646 <        for (j = 0; j < 3; j++)
647 <                v[j] += rad*(2.*pert[j] - 1.);
648 <        if (xm != NULL)
649 <                multv3(v, v, xm);
650 <        return(normalize(v) != 0.0);
638 >        int             inFront, outFront;
639 >        SDSpectralDF    *sdf;
640 >        float           coef[SDmaxCh];
641 >        int             nch, i;
642 >                                        /* check arguments */
643 >        if ((sv == NULL) | (outVec == NULL) | (inVec == NULL) | (sd == NULL))
644 >                return SDEargument;
645 >                                        /* whose side are we on? */
646 >        inFront = (inVec[2] > 0);
647 >        outFront = (outVec[2] > 0);
648 >                                        /* start with diffuse portion */
649 >        if (inFront & outFront) {
650 >                *sv = sd->rLambFront;
651 >                sdf = sd->rf;
652 >        } else if (!(inFront | outFront)) {
653 >                *sv = sd->rLambBack;
654 >                sdf = sd->rb;
655 >        } else if (inFront) {
656 >                *sv = sd->tLambFront;
657 >                sdf = (sd->tf != NULL) ? sd->tf : sd->tb;
658 >        } else /* outFront & !inFront */ {
659 >                *sv = sd->tLambBack;
660 >                sdf = (sd->tb != NULL) ? sd->tb : sd->tf;
661 >        }
662 >        sv->cieY *= 1./M_PI;
663 >                                        /* add non-diffuse components */
664 >        i = (sdf != NULL) ? sdf->ncomp : 0;
665 >        while (i-- > 0) {
666 >                nch = (*sdf->comp[i].func->getBSDFs)(coef, inVec, outVec,
667 >                                                        &sdf->comp[i]);
668 >                while (nch-- > 0) {
669 >                        c_cmix(&sv->spec, sv->cieY, &sv->spec,
670 >                                        coef[nch], &sdf->comp[i].cspec[nch]);
671 >                        sv->cieY += coef[nch];
672 >                }
673 >        }
674 >        c_ccvt(&sv->spec, C_CSXY);      /* make sure (x,y) is set */
675 >        return SDEnone;
676   }
677  
678 <
679 < #define  FEQ(a,b)       ((a)-(b) <= 1e-7 && (b)-(a) <= 1e-7)
680 <
480 < static int
481 < addrot(                 /* compute rotation (x,y,z) => (xp,yp,zp) */
482 <        char *xfarg[],
483 <        FVECT xp,
484 <        FVECT yp,
485 <        FVECT zp
486 < )
678 > /* Compute directional hemispherical scattering at this incident angle */
679 > double
680 > SDdirectHemi(const FVECT inVec, int sflags, const SDData *sd)
681   {
682 <        static char     bufs[3][16];
683 <        int     bn = 0;
684 <        char    **xfp = xfarg;
685 <        double  theta;
682 >        double          hsum;
683 >        SDSpectralDF    *rdf, *tdf;
684 >        const SDCDst    *cd;
685 >        int             i;
686 >                                        /* check arguments */
687 >        if ((inVec == NULL) | (sd == NULL))
688 >                return .0;
689 >                                        /* gather diffuse components */
690 >        if (inVec[2] > 0) {
691 >                hsum = sd->rLambFront.cieY;
692 >                rdf = sd->rf;
693 >                tdf = (sd->tf != NULL) ? sd->tf : sd->tb;
694 >        } else /* !inFront */ {
695 >                hsum = sd->rLambBack.cieY;
696 >                rdf = sd->rb;
697 >                tdf = (sd->tb != NULL) ? sd->tb : sd->tf;
698 >        }
699 >        if ((sflags & SDsampDf+SDsampR) != SDsampDf+SDsampR)
700 >                hsum = .0;
701 >        if ((sflags & SDsampDf+SDsampT) == SDsampDf+SDsampT)
702 >                hsum += (inVec[2] > 0) ?
703 >                                sd->tLambFront.cieY : sd->tLambBack.cieY;
704 >                                        /* gather non-diffuse components */
705 >        i = (((sflags & SDsampSp+SDsampR) == SDsampSp+SDsampR) &
706 >                        (rdf != NULL)) ? rdf->ncomp : 0;
707 >        while (i-- > 0) {               /* non-diffuse reflection */
708 >                cd = (*rdf->comp[i].func->getCDist)(inVec, &rdf->comp[i]);
709 >                if (cd != NULL)
710 >                        hsum += cd->cTotal;
711 >        }
712 >        i = (((sflags & SDsampSp+SDsampT) == SDsampSp+SDsampT) &
713 >                        (tdf != NULL)) ? tdf->ncomp : 0;
714 >        while (i-- > 0) {               /* non-diffuse transmission */
715 >                cd = (*tdf->comp[i].func->getCDist)(inVec, &tdf->comp[i]);
716 >                if (cd != NULL)
717 >                        hsum += cd->cTotal;
718 >        }
719 >        return hsum;
720 > }
721  
722 <        if (yp[2]*yp[2] + zp[2]*zp[2] < 2.*FTINY*FTINY) {
723 <                /* Special case for X' along Z-axis */
724 <                theta = -atan2(yp[0], yp[1]);
725 <                *xfp++ = "-ry";
726 <                *xfp++ = xp[2] < 0.0 ? "90" : "-90";
727 <                *xfp++ = "-rz";
728 <                sprintf(bufs[bn], "%f", theta*(180./PI));
729 <                *xfp++ = bufs[bn++];
730 <                return(xfp - xfarg);
722 > /* Sample BSDF direction based on the given random variable */
723 > SDError
724 > SDsampBSDF(SDValue *sv, FVECT ioVec, double randX, int sflags, const SDData *sd)
725 > {
726 >        SDError         ec;
727 >        FVECT           inVec;
728 >        int             inFront;
729 >        SDSpectralDF    *rdf, *tdf;
730 >        double          rdiff;
731 >        float           coef[SDmaxCh];
732 >        int             i, j, n, nr;
733 >        SDComponent     *sdc;
734 >        const SDCDst    **cdarr = NULL;
735 >                                        /* check arguments */
736 >        if ((sv == NULL) | (ioVec == NULL) | (sd == NULL) |
737 >                        (randX < 0) | (randX >= 1.))
738 >                return SDEargument;
739 >                                        /* whose side are we on? */
740 >        VCOPY(inVec, ioVec);
741 >        inFront = (inVec[2] > 0);
742 >                                        /* remember diffuse portions */
743 >        if (inFront) {
744 >                *sv = sd->rLambFront;
745 >                rdf = sd->rf;
746 >                tdf = (sd->tf != NULL) ? sd->tf : sd->tb;
747 >        } else /* !inFront */ {
748 >                *sv = sd->rLambBack;
749 >                rdf = sd->rb;
750 >                tdf = (sd->tb != NULL) ? sd->tb : sd->tf;
751          }
752 <        theta = atan2(yp[2], zp[2]);
753 <        if (!FEQ(theta,0.0)) {
754 <                *xfp++ = "-rx";
755 <                sprintf(bufs[bn], "%f", theta*(180./PI));
756 <                *xfp++ = bufs[bn++];
752 >        if ((sflags & SDsampDf+SDsampR) != SDsampDf+SDsampR)
753 >                sv->cieY = .0;
754 >        rdiff = sv->cieY;
755 >        if ((sflags & SDsampDf+SDsampT) == SDsampDf+SDsampT)
756 >                sv->cieY += inFront ? sd->tLambFront.cieY : sd->tLambBack.cieY;
757 >                                        /* gather non-diffuse components */
758 >        i = nr = (((sflags & SDsampSp+SDsampR) == SDsampSp+SDsampR) &
759 >                        (rdf != NULL)) ? rdf->ncomp : 0;
760 >        j = (((sflags & SDsampSp+SDsampT) == SDsampSp+SDsampT) &
761 >                        (tdf != NULL)) ? tdf->ncomp : 0;
762 >        n = i + j;
763 >        if (n > 0 && (cdarr = (const SDCDst **)malloc(n*sizeof(SDCDst *))) == NULL)
764 >                return SDEmemory;
765 >        while (j-- > 0) {               /* non-diffuse transmission */
766 >                cdarr[i+j] = (*tdf->comp[j].func->getCDist)(inVec, &tdf->comp[j]);
767 >                if (cdarr[i+j] == NULL)
768 >                        cdarr[i+j] = &SDemptyCD;
769 >                sv->cieY += cdarr[i+j]->cTotal;
770          }
771 <        theta = asin(-xp[2]);
772 <        if (!FEQ(theta,0.0)) {
773 <                *xfp++ = "-ry";
774 <                sprintf(bufs[bn], " %f", theta*(180./PI));
775 <                *xfp++ = bufs[bn++];
771 >        while (i-- > 0) {               /* non-diffuse reflection */
772 >                cdarr[i] = (*rdf->comp[i].func->getCDist)(inVec, &rdf->comp[i]);
773 >                if (cdarr[i] == NULL)
774 >                        cdarr[i] = &SDemptyCD;
775 >                sv->cieY += cdarr[i]->cTotal;
776          }
777 <        theta = atan2(xp[1], xp[0]);
778 <        if (!FEQ(theta,0.0)) {
779 <                *xfp++ = "-rz";
780 <                sprintf(bufs[bn], "%f", theta*(180./PI));
519 <                *xfp++ = bufs[bn++];
777 >        if (sv->cieY <= 1e-6) {         /* anything to sample? */
778 >                sv->cieY = .0;
779 >                memset(ioVec, 0, sizeof(FVECT));
780 >                return SDEnone;
781          }
782 <        *xfp = NULL;
783 <        return(xfp - xfarg);
782 >                                        /* scale random variable */
783 >        randX *= sv->cieY;
784 >                                        /* diffuse reflection? */
785 >        if (randX < rdiff) {
786 >                SDdiffuseSamp(ioVec, inFront, randX/rdiff);
787 >                goto done;
788 >        }
789 >        randX -= rdiff;
790 >                                        /* diffuse transmission? */
791 >        if ((sflags & SDsampDf+SDsampT) == SDsampDf+SDsampT) {
792 >                const SDValue   *sdt = inFront ? &sd->tLambFront : &sd->tLambBack;
793 >                if (randX < sdt->cieY) {
794 >                        sv->spec = sdt->spec;
795 >                        SDdiffuseSamp(ioVec, !inFront, randX/sdt->cieY);
796 >                        goto done;
797 >                }
798 >                randX -= sdt->cieY;
799 >        }
800 >                                        /* else one of cumulative dist. */
801 >        for (i = 0; i < n && randX >= cdarr[i]->cTotal; i++)
802 >                randX -= cdarr[i]->cTotal;
803 >        if (i >= n)
804 >                return SDEinternal;
805 >                                        /* compute sample direction */
806 >        sdc = (i < nr) ? &rdf->comp[i] : &tdf->comp[i-nr];
807 >        ec = (*sdc->func->sampCDist)(ioVec, randX/cdarr[i]->cTotal, cdarr[i]);
808 >        if (ec)
809 >                return ec;
810 >                                        /* compute color */
811 >        j = (*sdc->func->getBSDFs)(coef, inVec, ioVec, sdc);
812 >        if (j <= 0) {
813 >                sprintf(SDerrorDetail, "BSDF \"%s\" sampling value error",
814 >                                sd->name);
815 >                return SDEinternal;
816 >        }
817 >        sv->spec = sdc->cspec[0];
818 >        rdiff = coef[0];
819 >        while (--j) {
820 >                c_cmix(&sv->spec, rdiff, &sv->spec, coef[j], &sdc->cspec[j]);
821 >                rdiff += coef[j];
822 >        }
823 > done:
824 >        if (cdarr != NULL)
825 >                free(cdarr);
826 >        c_ccvt(&sv->spec, C_CSXY);      /* make sure (x,y) is set */
827 >        return SDEnone;
828   }
829  
830 + /* Compute World->BSDF transform from surface normal and up (Y) vector */
831 + SDError
832 + SDcompXform(RREAL vMtx[3][3], const FVECT sNrm, const FVECT uVec)
833 + {
834 +        if ((vMtx == NULL) | (sNrm == NULL) | (uVec == NULL))
835 +                return SDEargument;
836 +        VCOPY(vMtx[2], sNrm);
837 +        if (normalize(vMtx[2]) == 0)
838 +                return SDEargument;
839 +        fcross(vMtx[0], uVec, vMtx[2]);
840 +        if (normalize(vMtx[0]) == 0)
841 +                return SDEargument;
842 +        fcross(vMtx[1], vMtx[2], vMtx[0]);
843 +        return SDEnone;
844 + }
845  
846 < int
847 < getBSDF_xfm(            /* compute BSDF orient. -> world orient. transform */
848 <        MAT4 xm,
529 <        FVECT nrm,
530 <        UpDir ud
531 < )
846 > /* Compute inverse transform */
847 > SDError
848 > SDinvXform(RREAL iMtx[3][3], RREAL vMtx[3][3])
849   {
850 <        char    *xfargs[7];
851 <        XF      myxf;
535 <        FVECT   updir, xdest, ydest;
850 >        RREAL   mTmp[3][3];
851 >        double  d;
852  
853 <        updir[0] = updir[1] = updir[2] = 0.;
854 <        switch (ud) {
855 <        case UDzneg:
856 <                updir[2] = -1.;
857 <                break;
858 <        case UDyneg:
859 <                updir[1] = -1.;
860 <                break;
861 <        case UDxneg:
862 <                updir[0] = -1.;
547 <                break;
548 <        case UDxpos:
549 <                updir[0] = 1.;
550 <                break;
551 <        case UDypos:
552 <                updir[1] = 1.;
553 <                break;
554 <        case UDzpos:
555 <                updir[2] = 1.;
556 <                break;
557 <        case UDunknown:
558 <                return(0);
853 >        if ((iMtx == NULL) | (vMtx == NULL))
854 >                return SDEargument;
855 >                                        /* compute determinant */
856 >        mTmp[0][0] = vMtx[2][2]*vMtx[1][1] - vMtx[2][1]*vMtx[1][2];
857 >        mTmp[0][1] = vMtx[2][1]*vMtx[0][2] - vMtx[2][2]*vMtx[0][1];
858 >        mTmp[0][2] = vMtx[1][2]*vMtx[0][1] - vMtx[1][1]*vMtx[0][2];
859 >        d = vMtx[0][0]*mTmp[0][0] + vMtx[1][0]*mTmp[0][1] + vMtx[2][0]*mTmp[0][2];
860 >        if (d == 0) {
861 >                strcpy(SDerrorDetail, "Zero determinant in matrix inversion");
862 >                return SDEargument;
863          }
864 <        fcross(xdest, updir, nrm);
865 <        if (normalize(xdest) == 0.0)
866 <                return(0);
867 <        fcross(ydest, nrm, xdest);
868 <        xf(&myxf, addrot(xfargs, xdest, ydest, nrm), xfargs);
869 <        copymat4(xm, myxf.xfm);
870 <        return(1);
864 >        d = 1./d;                       /* invert matrix */
865 >        mTmp[0][0] *= d; mTmp[0][1] *= d; mTmp[0][2] *= d;
866 >        mTmp[1][0] = d*(vMtx[2][0]*vMtx[1][2] - vMtx[2][2]*vMtx[1][0]);
867 >        mTmp[1][1] = d*(vMtx[2][2]*vMtx[0][0] - vMtx[2][0]*vMtx[0][2]);
868 >        mTmp[1][2] = d*(vMtx[1][0]*vMtx[0][2] - vMtx[1][2]*vMtx[0][0]);
869 >        mTmp[2][0] = d*(vMtx[2][1]*vMtx[1][0] - vMtx[2][0]*vMtx[1][1]);
870 >        mTmp[2][1] = d*(vMtx[2][0]*vMtx[0][1] - vMtx[2][1]*vMtx[0][0]);
871 >        mTmp[2][2] = d*(vMtx[1][1]*vMtx[0][0] - vMtx[1][0]*vMtx[0][1]);
872 >        memcpy(iMtx, mTmp, sizeof(mTmp));
873 >        return SDEnone;
874 > }
875 >
876 > /* Transform and normalize direction (column) vector */
877 > SDError
878 > SDmapDir(FVECT resVec, RREAL vMtx[3][3], const FVECT inpVec)
879 > {
880 >        FVECT   vTmp;
881 >
882 >        if ((resVec == NULL) | (inpVec == NULL))
883 >                return SDEargument;
884 >        if (vMtx == NULL) {             /* assume they just want to normalize */
885 >                if (resVec != inpVec)
886 >                        VCOPY(resVec, inpVec);
887 >                return (normalize(resVec) > 0) ? SDEnone : SDEargument;
888 >        }
889 >        vTmp[0] = DOT(vMtx[0], inpVec);
890 >        vTmp[1] = DOT(vMtx[1], inpVec);
891 >        vTmp[2] = DOT(vMtx[2], inpVec);
892 >        if (normalize(vTmp) == 0)
893 >                return SDEargument;
894 >        VCOPY(resVec, vTmp);
895 >        return SDEnone;
896   }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines