ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/bsdf_m.c
Revision: 3.39
Committed: Wed Sep 11 00:24:03 2019 UTC (5 years, 1 month ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad5R3
Changes since 3.38: +9 -9 lines
Log Message: 
Reversed order of matrix look-up macros to be more consistent

File Contents

# User Rev Content
1 greg 3.2 #ifndef lint
2 greg 3.39 static const char RCSid[] = "$Id: bsdf_m.c,v 3.38 2018/05/10 22:55:35 greg Exp $";
3 greg 3.2 #endif
4 greg 3.1 /*
5     * bsdf_m.c
6     *
7     * Definitions supporting BSDF matrices
8     *
9     * Created by Greg Ward on 2/2/11.
10     * Copyright 2011 Anyhere Software. All rights reserved.
11     *
12     */
13    
14 greg 3.16 #define _USE_MATH_DEFINES
15 greg 3.3 #include "rtio.h"
16 greg 3.1 #include <math.h>
17     #include <ctype.h>
18     #include "ezxml.h"
19     #include "bsdf.h"
20     #include "bsdf_m.h"
21    
22     /* Function return codes */
23     #define RC_GOOD 1
24     #define RC_FAIL 0
25     #define RC_FORMERR (-1)
26     #define RC_DATERR (-2)
27     #define RC_UNSUPP (-3)
28     #define RC_INTERR (-4)
29     #define RC_MEMERR (-5)
30    
31 greg 3.25 ANGLE_BASIS abase_list[MAXABASES] = {
32 greg 3.1 {
33     "LBNL/Klems Full", 145,
34 greg 3.22 { {0., 1},
35 greg 3.1 {5., 8},
36     {15., 16},
37     {25., 20},
38     {35., 24},
39     {45., 24},
40     {55., 24},
41     {65., 16},
42     {75., 12},
43     {90., 0} }
44     }, {
45     "LBNL/Klems Half", 73,
46 greg 3.22 { {0., 1},
47 greg 3.1 {6.5, 8},
48     {19.5, 12},
49     {32.5, 16},
50     {46.5, 20},
51     {61.5, 12},
52     {76.5, 4},
53     {90., 0} }
54     }, {
55     "LBNL/Klems Quarter", 41,
56 greg 3.22 { {0., 1},
57 greg 3.1 {9., 8},
58     {27., 12},
59     {46., 12},
60     {66., 8},
61     {90., 0} }
62     }
63     };
64    
65 greg 3.25 int nabases = 3; /* current number of defined bases */
66 greg 3.1
67 greg 3.30 C_COLOR mtx_RGB_prim[3]; /* our RGB primaries */
68     float mtx_RGB_coef[3]; /* corresponding Y coefficients */
69 greg 3.29
70     enum {mtx_Y, mtx_X, mtx_Z}; /* matrix components (mtx_Y==0) */
71    
72     /* check if two real values are near enough to equal */
73 greg 3.1 static int
74     fequal(double a, double b)
75     {
76 greg 3.10 if (b != 0)
77 greg 3.1 a = a/b - 1.;
78     return (a <= 1e-6) & (a >= -1e-6);
79     }
80    
81 greg 3.29 /* convert error to standard BSDF code */
82 greg 3.1 static SDError
83     convert_errcode(int ec)
84     {
85     switch (ec) {
86     case RC_GOOD:
87     return SDEnone;
88     case RC_FORMERR:
89     return SDEformat;
90     case RC_DATERR:
91     return SDEdata;
92     case RC_UNSUPP:
93     return SDEsupport;
94     case RC_INTERR:
95     return SDEinternal;
96     case RC_MEMERR:
97     return SDEmemory;
98     }
99     return SDEunknown;
100     }
101    
102 greg 3.29 /* allocate a BSDF matrix of the given size */
103 greg 3.1 static SDMat *
104     SDnewMatrix(int ni, int no)
105     {
106     SDMat *sm;
107    
108     if ((ni <= 0) | (no <= 0)) {
109     strcpy(SDerrorDetail, "Empty BSDF matrix request");
110     return NULL;
111     }
112     sm = (SDMat *)malloc(sizeof(SDMat) + (ni*no - 1)*sizeof(float));
113     if (sm == NULL) {
114     sprintf(SDerrorDetail, "Cannot allocate %dx%d BSDF matrix",
115     ni, no);
116     return NULL;
117     }
118     memset(sm, 0, sizeof(SDMat)-sizeof(float));
119     sm->ninc = ni;
120     sm->nout = no;
121    
122     return sm;
123     }
124    
125     /* Free a BSDF matrix */
126 greg 3.29 void
127     SDfreeMatrix(void *ptr)
128     {
129     SDMat *mp = (SDMat *)ptr;
130    
131     if (mp->chroma != NULL) free(mp->chroma);
132     free(ptr);
133     }
134 greg 3.1
135 greg 3.32 /* compute square of real value */
136     static double sq(double x) { return x*x; }
137    
138 greg 3.25 /* Get vector for this angle basis index (front exiting) */
139     int
140 greg 3.9 fo_getvec(FVECT v, double ndxr, void *p)
141 greg 3.1 {
142 greg 3.9 ANGLE_BASIS *ab = (ANGLE_BASIS *)p;
143     int ndx = (int)ndxr;
144     double randX = ndxr - ndx;
145 greg 3.1 double rx[2];
146     int li;
147 greg 3.32 double azi, d;
148 greg 3.1
149 greg 3.9 if ((ndxr < 0) | (ndx >= ab->nangles))
150 greg 3.1 return RC_FAIL;
151     for (li = 0; ndx >= ab->lat[li].nphis; li++)
152     ndx -= ab->lat[li].nphis;
153     SDmultiSamp(rx, 2, randX);
154 greg 3.33 d = (1. - rx[0])*sq(cos(M_PI/180.*ab->lat[li].tmin)) +
155     rx[0]*sq(cos(M_PI/180.*ab->lat[li+1].tmin));
156 greg 3.32 v[2] = d = sqrt(d); /* cos(pol) */
157 greg 3.1 azi = 2.*M_PI*(ndx + rx[1] - .5)/ab->lat[li].nphis;
158     d = sqrt(1. - d*d); /* sin(pol) */
159     v[0] = cos(azi)*d;
160     v[1] = sin(azi)*d;
161     return RC_GOOD;
162     }
163    
164 greg 3.25 /* Get index corresponding to the given vector (front exiting) */
165     int
166 greg 3.8 fo_getndx(const FVECT v, void *p)
167 greg 3.1 {
168 greg 3.9 ANGLE_BASIS *ab = (ANGLE_BASIS *)p;
169 greg 3.1 int li, ndx;
170 greg 3.17 double pol, azi;
171 greg 3.1
172     if (v == NULL)
173     return -1;
174 greg 3.10 if ((v[2] < 0) | (v[2] > 1.))
175 greg 3.1 return -1;
176 greg 3.27 pol = 180.0/M_PI*Acos(v[2]);
177 greg 3.1 azi = 180.0/M_PI*atan2(v[1], v[0]);
178     if (azi < 0.0) azi += 360.0;
179     for (li = 1; ab->lat[li].tmin <= pol; li++)
180     if (!ab->lat[li].nphis)
181     return -1;
182     --li;
183     ndx = (int)((1./360.)*azi*ab->lat[li].nphis + 0.5);
184     if (ndx >= ab->lat[li].nphis) ndx = 0;
185     while (li--)
186     ndx += ab->lat[li].nphis;
187     return ndx;
188     }
189    
190 greg 3.25 /* Get projected solid angle for this angle basis index (universal) */
191     double
192 greg 3.8 io_getohm(int ndx, void *p)
193 greg 3.1 {
194 greg 3.36 static void *last_p = NULL;
195 greg 3.1 static int last_li = -1;
196     static double last_ohm;
197     ANGLE_BASIS *ab = (ANGLE_BASIS *)p;
198     int li;
199     double theta, theta1;
200    
201     if ((ndx < 0) | (ndx >= ab->nangles))
202     return -1.;
203     for (li = 0; ndx >= ab->lat[li].nphis; li++)
204     ndx -= ab->lat[li].nphis;
205 greg 3.36 if ((p == last_p) & (li == last_li)) /* cached latitude? */
206 greg 3.1 return last_ohm;
207 greg 3.36 last_p = p;
208 greg 3.1 last_li = li;
209 greg 3.21 theta = M_PI/180. * ab->lat[li].tmin;
210 greg 3.1 theta1 = M_PI/180. * ab->lat[li+1].tmin;
211     return last_ohm = M_PI*(sq(cos(theta)) - sq(cos(theta1))) /
212     (double)ab->lat[li].nphis;
213     }
214    
215 greg 3.25 /* Get vector for this angle basis index (back incident) */
216     int
217 greg 3.9 bi_getvec(FVECT v, double ndxr, void *p)
218 greg 3.1 {
219 greg 3.9 if (!fo_getvec(v, ndxr, p))
220 greg 3.1 return RC_FAIL;
221    
222     v[0] = -v[0];
223     v[1] = -v[1];
224     v[2] = -v[2];
225    
226     return RC_GOOD;
227     }
228    
229 greg 3.25 /* Get index corresponding to the vector (back incident) */
230     int
231 greg 3.8 bi_getndx(const FVECT v, void *p)
232 greg 3.1 {
233     FVECT v2;
234    
235     v2[0] = -v[0];
236     v2[1] = -v[1];
237     v2[2] = -v[2];
238    
239 greg 3.8 return fo_getndx(v2, p);
240     }
241    
242 greg 3.25 /* Get vector for this angle basis index (back exiting) */
243     int
244 greg 3.9 bo_getvec(FVECT v, double ndxr, void *p)
245 greg 3.8 {
246 greg 3.9 if (!fo_getvec(v, ndxr, p))
247 greg 3.8 return RC_FAIL;
248    
249     v[2] = -v[2];
250    
251     return RC_GOOD;
252     }
253    
254 greg 3.25 /* Get index corresponding to the vector (back exiting) */
255     int
256 greg 3.8 bo_getndx(const FVECT v, void *p)
257     {
258     FVECT v2;
259    
260     v2[0] = v[0];
261     v2[1] = v[1];
262     v2[2] = -v[2];
263    
264     return fo_getndx(v2, p);
265     }
266    
267 greg 3.25 /* Get vector for this angle basis index (front incident) */
268     int
269 greg 3.9 fi_getvec(FVECT v, double ndxr, void *p)
270 greg 3.8 {
271 greg 3.9 if (!fo_getvec(v, ndxr, p))
272 greg 3.8 return RC_FAIL;
273    
274     v[0] = -v[0];
275     v[1] = -v[1];
276    
277     return RC_GOOD;
278     }
279    
280 greg 3.25 /* Get index corresponding to the vector (front incident) */
281     int
282 greg 3.8 fi_getndx(const FVECT v, void *p)
283     {
284     FVECT v2;
285    
286     v2[0] = -v[0];
287     v2[1] = -v[1];
288     v2[2] = v[2];
289    
290     return fo_getndx(v2, p);
291 greg 3.1 }
292    
293 greg 3.29 /* Get color or grayscale value for BSDF for the given direction pair */
294     int
295     mBSDF_color(float coef[], const SDMat *dp, int i, int o)
296     {
297     C_COLOR cxy;
298    
299 greg 3.39 coef[0] = mBSDF_value(dp, o, i);
300 greg 3.29 if (dp->chroma == NULL)
301     return 1; /* grayscale */
302    
303 greg 3.39 c_decodeChroma(&cxy, mBSDF_chroma(dp,o,i));
304 greg 3.29 c_toSharpRGB(&cxy, coef[0], coef);
305     coef[0] *= mtx_RGB_coef[0];
306     coef[1] *= mtx_RGB_coef[1];
307     coef[2] *= mtx_RGB_coef[2];
308     return 3; /* RGB color */
309     }
310    
311 greg 3.1 /* load custom BSDF angle basis */
312     static int
313     load_angle_basis(ezxml_t wab)
314     {
315     char *abname = ezxml_txt(ezxml_child(wab, "AngleBasisName"));
316     ezxml_t wbb;
317     int i;
318    
319     if (!abname || !*abname)
320     return RC_FAIL;
321     for (i = nabases; i--; )
322     if (!strcasecmp(abname, abase_list[i].name))
323     return RC_GOOD; /* assume it's the same */
324     if (nabases >= MAXABASES) {
325     sprintf(SDerrorDetail, "Out of angle bases reading '%s'",
326     abname);
327     return RC_INTERR;
328     }
329     strcpy(abase_list[nabases].name, abname);
330     abase_list[nabases].nangles = 0;
331     for (i = 0, wbb = ezxml_child(wab, "AngleBasisBlock");
332     wbb != NULL; i++, wbb = wbb->next) {
333     if (i >= MAXLATS) {
334     sprintf(SDerrorDetail, "Too many latitudes for '%s'",
335     abname);
336     return RC_INTERR;
337     }
338     abase_list[nabases].lat[i+1].tmin = atof(ezxml_txt(
339     ezxml_child(ezxml_child(wbb,
340     "ThetaBounds"), "UpperTheta")));
341     if (!i)
342 greg 3.21 abase_list[nabases].lat[0].tmin = 0;
343 greg 3.1 else if (!fequal(atof(ezxml_txt(ezxml_child(ezxml_child(wbb,
344     "ThetaBounds"), "LowerTheta"))),
345     abase_list[nabases].lat[i].tmin)) {
346     sprintf(SDerrorDetail, "Theta values disagree in '%s'",
347 greg 3.12 abname);
348 greg 3.1 return RC_DATERR;
349     }
350     abase_list[nabases].nangles +=
351     abase_list[nabases].lat[i].nphis =
352     atoi(ezxml_txt(ezxml_child(wbb, "nPhis")));
353     if (abase_list[nabases].lat[i].nphis <= 0 ||
354     (abase_list[nabases].lat[i].nphis == 1 &&
355     abase_list[nabases].lat[i].tmin > FTINY)) {
356     sprintf(SDerrorDetail, "Illegal phi count in '%s'",
357 greg 3.12 abname);
358 greg 3.1 return RC_DATERR;
359     }
360     }
361     abase_list[nabases++].lat[i].nphis = 0;
362     return RC_GOOD;
363     }
364    
365     /* compute min. proj. solid angle and max. direct hemispherical scattering */
366     static int
367     get_extrema(SDSpectralDF *df)
368     {
369     SDMat *dp = (SDMat *)df->comp[0].dist;
370     double *ohma;
371     int i, o;
372     /* initialize extrema */
373     df->minProjSA = M_PI;
374     df->maxHemi = .0;
375     ohma = (double *)malloc(dp->nout*sizeof(double));
376     if (ohma == NULL)
377     return RC_MEMERR;
378     /* get outgoing solid angles */
379     for (o = dp->nout; o--; )
380     if ((ohma[o] = mBSDF_outohm(dp,o)) < df->minProjSA)
381     df->minProjSA = ohma[o];
382     /* compute hemispherical sums */
383     for (i = dp->ninc; i--; ) {
384     double hemi = .0;
385     for (o = dp->nout; o--; )
386 greg 3.39 hemi += ohma[o] * mBSDF_value(dp, o, i);
387 greg 3.1 if (hemi > df->maxHemi)
388     df->maxHemi = hemi;
389     }
390     free(ohma);
391     /* need incoming solid angles, too? */
392 greg 3.5 if ((dp->ib_ohm != dp->ob_ohm) | (dp->ib_priv != dp->ob_priv)) {
393 greg 3.1 double ohm;
394     for (i = dp->ninc; i--; )
395     if ((ohm = mBSDF_incohm(dp,i)) < df->minProjSA)
396     df->minProjSA = ohm;
397     }
398     return (df->maxHemi <= 1.01);
399     }
400    
401     /* load BSDF distribution for this wavelength */
402     static int
403 greg 3.29 load_bsdf_data(SDData *sd, ezxml_t wdb, int ct, int rowinc)
404 greg 3.1 {
405     SDSpectralDF *df;
406     SDMat *dp;
407     char *sdata;
408     int inbi, outbi;
409     int i;
410     /* allocate BSDF component */
411     sdata = ezxml_txt(ezxml_child(wdb, "WavelengthDataDirection"));
412 greg 3.15 if (!sdata)
413     return RC_FAIL;
414 greg 3.8 /*
415     * Remember that front and back are reversed from WINDOW 6 orientations
416     */
417 greg 3.23 if (!strcasecmp(sdata, "Transmission Front")) {
418 greg 3.29 if (sd->tb == NULL && (sd->tb = SDnewSpectralDF(3)) == NULL)
419 greg 3.24 return RC_MEMERR;
420     df = sd->tb;
421     } else if (!strcasecmp(sdata, "Transmission Back")) {
422 greg 3.29 if (sd->tf == NULL && (sd->tf = SDnewSpectralDF(3)) == NULL)
423 greg 3.1 return RC_MEMERR;
424     df = sd->tf;
425     } else if (!strcasecmp(sdata, "Reflection Front")) {
426 greg 3.29 if (sd->rb == NULL && (sd->rb = SDnewSpectralDF(3)) == NULL)
427 greg 3.6 return RC_MEMERR;
428     df = sd->rb;
429     } else if (!strcasecmp(sdata, "Reflection Back")) {
430 greg 3.29 if (sd->rf == NULL && (sd->rf = SDnewSpectralDF(3)) == NULL)
431 greg 3.1 return RC_MEMERR;
432     df = sd->rf;
433     } else
434     return RC_FAIL;
435 greg 3.29 /* free previous matrix if any */
436     if (df->comp[ct].dist != NULL) {
437     SDfreeMatrix(df->comp[ct].dist);
438     df->comp[ct].dist = NULL;
439     }
440 greg 3.1 /* get angle bases */
441     sdata = ezxml_txt(ezxml_child(wdb,"ColumnAngleBasis"));
442     if (!sdata || !*sdata) {
443     sprintf(SDerrorDetail, "Missing column basis for BSDF '%s'",
444     sd->name);
445     return RC_FORMERR;
446     }
447     for (inbi = nabases; inbi--; )
448 greg 3.4 if (!strcasecmp(sdata, abase_list[inbi].name))
449 greg 3.1 break;
450     if (inbi < 0) {
451 greg 3.4 sprintf(SDerrorDetail, "Undefined ColumnAngleBasis '%s'", sdata);
452 greg 3.1 return RC_FORMERR;
453     }
454     sdata = ezxml_txt(ezxml_child(wdb,"RowAngleBasis"));
455     if (!sdata || !*sdata) {
456     sprintf(SDerrorDetail, "Missing row basis for BSDF '%s'",
457     sd->name);
458     return RC_FORMERR;
459     }
460     for (outbi = nabases; outbi--; )
461 greg 3.4 if (!strcasecmp(sdata, abase_list[outbi].name))
462 greg 3.1 break;
463     if (outbi < 0) {
464 greg 3.4 sprintf(SDerrorDetail, "Undefined RowAngleBasis '%s'", sdata);
465 greg 3.1 return RC_FORMERR;
466     }
467     /* allocate BSDF matrix */
468     dp = SDnewMatrix(abase_list[inbi].nangles, abase_list[outbi].nangles);
469     if (dp == NULL)
470     return RC_MEMERR;
471 greg 3.5 dp->ib_priv = &abase_list[inbi];
472     dp->ob_priv = &abase_list[outbi];
473 greg 3.1 if (df == sd->tf) {
474 greg 3.23 dp->ib_vec = &fi_getvec;
475     dp->ib_ndx = &fi_getndx;
476     dp->ob_vec = &bo_getvec;
477     dp->ob_ndx = &bo_getndx;
478     } else if (df == sd->tb) {
479     dp->ib_vec = &bi_getvec;
480     dp->ib_ndx = &bi_getndx;
481     dp->ob_vec = &fo_getvec;
482     dp->ob_ndx = &fo_getndx;
483 greg 3.1 } else if (df == sd->rf) {
484 greg 3.8 dp->ib_vec = &fi_getvec;
485     dp->ib_ndx = &fi_getndx;
486     dp->ob_vec = &fo_getvec;
487     dp->ob_ndx = &fo_getndx;
488 greg 3.1 } else /* df == sd->rb */ {
489 greg 3.8 dp->ib_vec = &bi_getvec;
490     dp->ib_ndx = &bi_getndx;
491     dp->ob_vec = &bo_getvec;
492     dp->ob_ndx = &bo_getndx;
493 greg 3.1 }
494 greg 3.8 dp->ib_ohm = &io_getohm;
495     dp->ob_ohm = &io_getohm;
496 greg 3.29 df->comp[ct].dist = dp;
497     df->comp[ct].func = &SDhandleMtx;
498 greg 3.1 /* read BSDF data */
499 greg 3.15 sdata = ezxml_txt(ezxml_child(wdb, "ScatteringData"));
500 greg 3.1 if (!sdata || !*sdata) {
501     sprintf(SDerrorDetail, "Missing BSDF ScatteringData in '%s'",
502     sd->name);
503     return RC_FORMERR;
504     }
505     for (i = 0; i < dp->ninc*dp->nout; i++) {
506 greg 3.28 char *sdnext = fskip(sdata);
507     double val;
508 greg 3.1 if (sdnext == NULL) {
509     sprintf(SDerrorDetail,
510     "Bad/missing BSDF ScatteringData in '%s'",
511     sd->name);
512     return RC_FORMERR;
513     }
514 greg 3.15 while (isspace(*sdnext))
515 greg 3.1 sdnext++;
516     if (*sdnext == ',') sdnext++;
517 greg 3.28 if ((val = atof(sdata)) < 0)
518     val = 0; /* don't allow negative values */
519 greg 3.1 if (rowinc) {
520     int r = i/dp->nout;
521 greg 3.16 int c = i - r*dp->nout;
522 greg 3.39 mBSDF_value(dp,c,r) = val;
523 greg 3.1 } else
524 greg 3.28 dp->bsdf[i] = val;
525 greg 3.1 sdata = sdnext;
526     }
527 greg 3.29 return (ct == mtx_Y) ? get_extrema(df) : RC_GOOD;
528     }
529    
530     /* copy our RGB (x,y) primary chromaticities */
531     static void
532     copy_RGB_prims(C_COLOR cspec[])
533     {
534     if (mtx_RGB_coef[1] < .001) { /* need to initialize */
535     int i = 3;
536     while (i--) {
537     float rgb[3];
538     rgb[0] = rgb[1] = rgb[2] = .0f;
539     rgb[i] = 1.f;
540     mtx_RGB_coef[i] = c_fromSharpRGB(rgb, &mtx_RGB_prim[i]);
541     }
542     }
543     memcpy(cspec, mtx_RGB_prim, sizeof(mtx_RGB_prim));
544 greg 3.1 }
545    
546 greg 3.29 /* encode chromaticity if XYZ -- reduce to one channel in any case */
547     static SDSpectralDF *
548     encode_chroma(SDSpectralDF *df)
549     {
550     SDMat *mpx, *mpy, *mpz;
551     int n;
552    
553     if (df == NULL || df->ncomp != 3)
554     return df;
555    
556     mpy = (SDMat *)df->comp[mtx_Y].dist;
557     if (mpy == NULL) {
558     free(df);
559     return NULL;
560     }
561     mpx = (SDMat *)df->comp[mtx_X].dist;
562     mpz = (SDMat *)df->comp[mtx_Z].dist;
563     if (mpx == NULL || (mpx->ninc != mpy->ninc) | (mpx->nout != mpy->nout))
564     goto done;
565     if (mpz == NULL || (mpz->ninc != mpy->ninc) | (mpz->nout != mpy->nout))
566     goto done;
567     mpy->chroma = (C_CHROMA *)malloc(sizeof(C_CHROMA)*mpy->ninc*mpy->nout);
568     if (mpy->chroma == NULL)
569     goto done; /* XXX punt */
570     /* encode chroma values */
571     for (n = mpy->ninc*mpy->nout; n--; ) {
572     const double sum = mpx->bsdf[n] + mpy->bsdf[n] + mpz->bsdf[n];
573     C_COLOR cxy;
574     if (sum > .0)
575     c_cset(&cxy, mpx->bsdf[n]/sum, mpy->bsdf[n]/sum);
576     else
577     c_cset(&cxy, 1./3., 1./3.);
578     mpy->chroma[n] = c_encodeChroma(&cxy);
579     }
580     done: /* free X & Z channels */
581     if (mpx != NULL) SDfreeMatrix(mpx);
582     if (mpz != NULL) SDfreeMatrix(mpz);
583     if (mpy->chroma == NULL) /* grayscale after all? */
584     df->comp[0].cspec[0] = c_dfcolor;
585     else /* else copy RGB primaries */
586     copy_RGB_prims(df->comp[0].cspec);
587     df->ncomp = 1; /* return resized struct */
588     return (SDSpectralDF *)realloc(df, sizeof(SDSpectralDF));
589     }
590    
591     /* subtract minimum (diffuse) scattering amount from BSDF */
592 greg 3.1 static double
593 greg 3.29 subtract_min(C_COLOR *cs, SDMat *sm)
594 greg 3.1 {
595 greg 3.29 const int ncomp = 1 + 2*(sm->chroma != NULL);
596 greg 3.31 float min_coef[3], ymin, coef[3];
597 greg 3.29 int i, o, c;
598 greg 3.1
599 greg 3.29 min_coef[0] = min_coef[1] = min_coef[2] = FHUGE;
600     for (i = 0; i < sm->ninc; i++)
601     for (o = 0; o < sm->nout; o++) {
602     c = mBSDF_color(coef, sm, i, o);
603     while (c--)
604     if (coef[c] < min_coef[c])
605     min_coef[c] = coef[c];
606     }
607 greg 3.31 ymin = 0;
608 greg 3.29 for (c = ncomp; c--; )
609 greg 3.31 ymin += min_coef[c];
610     if (ymin <= .01/M_PI) /* not worth bothering about? */
611 greg 3.15 return .0;
612 greg 3.29 if (ncomp == 1) { /* subtract grayscale minimum */
613     for (i = sm->ninc*sm->nout; i--; )
614 greg 3.31 sm->bsdf[i] -= ymin;
615 greg 3.29 *cs = c_dfcolor;
616 greg 3.31 return M_PI*ymin;
617 greg 3.29 }
618     /* else subtract colored minimum */
619     for (i = 0; i < sm->ninc; i++)
620     for (o = 0; o < sm->nout; o++) {
621     C_COLOR cxy;
622     c = mBSDF_color(coef, sm, i, o);
623     while (c--)
624     coef[c] = (coef[c] - min_coef[c]) /
625     mtx_RGB_coef[c];
626 greg 3.30 if (c_fromSharpRGB(coef, &cxy) > 1e-5)
627 greg 3.39 mBSDF_chroma(sm,o,i) = c_encodeChroma(&cxy);
628     mBSDF_value(sm,o,i) -= ymin;
629 greg 3.29 }
630     /* return colored minimum */
631 greg 3.30 for (i = 3; i--; )
632     coef[i] = min_coef[i]/mtx_RGB_coef[i];
633     c_fromSharpRGB(coef, cs);
634 greg 3.1
635 greg 3.31 return M_PI*ymin;
636 greg 3.1 }
637    
638 greg 3.29 /* Extract and separate diffuse portion of BSDF & convert color */
639     static SDSpectralDF *
640 greg 3.1 extract_diffuse(SDValue *dv, SDSpectralDF *df)
641     {
642    
643 greg 3.29 df = encode_chroma(df); /* reduce XYZ to Y + chroma */
644 greg 3.1 if (df == NULL || df->ncomp <= 0) {
645     dv->spec = c_dfcolor;
646     dv->cieY = .0;
647 greg 3.29 return df;
648 greg 3.1 }
649 greg 3.29 /* subtract minimum value */
650     dv->cieY = subtract_min(&dv->spec, (SDMat *)df->comp[0].dist);
651 greg 3.15 df->maxHemi -= dv->cieY; /* adjust maximum hemispherical */
652 greg 3.38
653     c_ccvt(&dv->spec, C_CSXY); /* make sure (x,y) is set */
654 greg 3.29 return df;
655 greg 3.1 }
656    
657     /* Load a BSDF matrix from an open XML file */
658     SDError
659 greg 3.4 SDloadMtx(SDData *sd, ezxml_t wtl)
660 greg 3.1 {
661 greg 3.15 ezxml_t wld, wdb;
662     int rowIn;
663     char *txt;
664     int rval;
665     /* basic checks and data ordering */
666 greg 3.4 txt = ezxml_txt(ezxml_child(ezxml_child(wtl,
667     "DataDefinition"), "IncidentDataStructure"));
668     if (txt == NULL || !*txt) {
669 greg 3.1 sprintf(SDerrorDetail,
670 greg 3.4 "BSDF \"%s\": missing IncidentDataStructure",
671 greg 3.1 sd->name);
672     return SDEformat;
673     }
674     if (!strcasecmp(txt, "Rows"))
675     rowIn = 1;
676     else if (!strcasecmp(txt, "Columns"))
677     rowIn = 0;
678     else {
679     sprintf(SDerrorDetail,
680     "BSDF \"%s\": unsupported IncidentDataStructure",
681     sd->name);
682     return SDEsupport;
683     }
684 greg 3.18 /* get angle bases */
685     for (wld = ezxml_child(ezxml_child(wtl, "DataDefinition"), "AngleBasis");
686     wld != NULL; wld = wld->next) {
687     rval = load_angle_basis(wld);
688     if (rval < 0)
689     return convert_errcode(rval);
690     }
691 greg 3.15 /* load BSDF components */
692 greg 3.1 for (wld = ezxml_child(wtl, "WavelengthData");
693     wld != NULL; wld = wld->next) {
694 greg 3.29 const char *cnm = ezxml_txt(ezxml_child(wld,"Wavelength"));
695     int ct = -1;
696     if (!strcasecmp(cnm, "Visible"))
697     ct = mtx_Y;
698     else if (!strcasecmp(cnm, "CIE-X"))
699     ct = mtx_X;
700     else if (!strcasecmp(cnm, "CIE-Z"))
701     ct = mtx_Z;
702     else
703     continue;
704 greg 3.1 for (wdb = ezxml_child(wld, "WavelengthDataBlock");
705     wdb != NULL; wdb = wdb->next)
706 greg 3.29 if ((rval = load_bsdf_data(sd, wdb, ct, rowIn)) < 0)
707 greg 3.4 return convert_errcode(rval);
708 greg 3.1 }
709 greg 3.15 /* separate diffuse components */
710 greg 3.29 sd->rf = extract_diffuse(&sd->rLambFront, sd->rf);
711     sd->rb = extract_diffuse(&sd->rLambBack, sd->rb);
712 greg 3.26 if (sd->tf != NULL)
713 greg 3.29 sd->tf = extract_diffuse(&sd->tLamb, sd->tf);
714 greg 3.26 if (sd->tb != NULL)
715 greg 3.29 sd->tb = extract_diffuse(&sd->tLamb, sd->tb);
716 greg 3.15 /* return success */
717 greg 3.1 return SDEnone;
718     }
719    
720     /* Get Matrix BSDF value */
721     static int
722     SDgetMtxBSDF(float coef[SDmaxCh], const FVECT outVec,
723 greg 3.12 const FVECT inVec, SDComponent *sdc)
724 greg 3.1 {
725 greg 3.12 const SDMat *dp;
726 greg 3.1 int i_ndx, o_ndx;
727 greg 3.12 /* check arguments */
728     if ((coef == NULL) | (outVec == NULL) | (inVec == NULL) | (sdc == NULL)
729     || (dp = (SDMat *)sdc->dist) == NULL)
730     return 0;
731 greg 3.1 /* get angle indices */
732     i_ndx = mBSDF_incndx(dp, inVec);
733     o_ndx = mBSDF_outndx(dp, outVec);
734     /* try reciprocity if necessary */
735 greg 3.19 if ((i_ndx < 0) & (o_ndx < 0)) {
736 greg 3.1 i_ndx = mBSDF_incndx(dp, outVec);
737     o_ndx = mBSDF_outndx(dp, inVec);
738     }
739     if ((i_ndx < 0) | (o_ndx < 0))
740     return 0; /* nothing from this component */
741 greg 3.29
742     return mBSDF_color(coef, dp, i_ndx, o_ndx);
743 greg 3.1 }
744    
745 greg 3.12 /* Query solid angle for vector(s) */
746 greg 3.1 static SDError
747 greg 3.10 SDqueryMtxProjSA(double *psa, const FVECT v1, const RREAL *v2,
748 greg 3.12 int qflags, SDComponent *sdc)
749 greg 3.1 {
750 greg 3.12 const SDMat *dp;
751 greg 3.5 double inc_psa, out_psa;
752     /* check arguments */
753 greg 3.12 if ((psa == NULL) | (v1 == NULL) | (sdc == NULL) ||
754     (dp = (SDMat *)sdc->dist) == NULL)
755 greg 3.1 return SDEargument;
756 greg 3.10 if (v2 == NULL)
757     v2 = v1;
758 greg 3.5 /* get projected solid angles */
759 greg 3.10 out_psa = mBSDF_outohm(dp, mBSDF_outndx(dp, v1));
760     inc_psa = mBSDF_incohm(dp, mBSDF_incndx(dp, v2));
761 greg 3.12 if ((v1 != v2) & (out_psa <= 0) & (inc_psa <= 0)) {
762 greg 3.11 inc_psa = mBSDF_outohm(dp, mBSDF_outndx(dp, v2));
763     out_psa = mBSDF_incohm(dp, mBSDF_incndx(dp, v1));
764     }
765 greg 3.5
766     switch (qflags) { /* record based on flag settings */
767     case SDqueryMax:
768     if (inc_psa > psa[0])
769     psa[0] = inc_psa;
770     if (out_psa > psa[0])
771     psa[0] = out_psa;
772     break;
773     case SDqueryMin+SDqueryMax:
774 greg 3.13 if (inc_psa > psa[1])
775 greg 3.5 psa[1] = inc_psa;
776 greg 3.13 if (out_psa > psa[1])
777 greg 3.5 psa[1] = out_psa;
778     /* fall through */
779 greg 3.12 case SDqueryVal:
780     if (qflags == SDqueryVal)
781     psa[0] = M_PI;
782 greg 3.14 /* fall through */
783     case SDqueryMin:
784 greg 3.10 if ((inc_psa > 0) & (inc_psa < psa[0]))
785 greg 3.1 psa[0] = inc_psa;
786 greg 3.10 if ((out_psa > 0) & (out_psa < psa[0]))
787 greg 3.5 psa[0] = out_psa;
788     break;
789 greg 3.1 }
790 greg 3.5 /* make sure it's legal */
791 greg 3.10 return (psa[0] <= 0) ? SDEinternal : SDEnone;
792 greg 3.1 }
793    
794     /* Compute new cumulative distribution from BSDF */
795     static int
796     make_cdist(SDMatCDst *cd, const FVECT inVec, SDMat *dp, int rev)
797     {
798     const unsigned maxval = ~0;
799     double *cmtab, scale;
800     int o;
801    
802     cmtab = (double *)malloc((cd->calen+1)*sizeof(double));
803     if (cmtab == NULL)
804     return 0;
805     cmtab[0] = .0;
806     for (o = 0; o < cd->calen; o++) {
807     if (rev)
808 greg 3.39 cmtab[o+1] = mBSDF_value(dp, cd->indx, o) *
809 greg 3.1 (*dp->ib_ohm)(o, dp->ib_priv);
810     else
811 greg 3.39 cmtab[o+1] = mBSDF_value(dp, o, cd->indx) *
812 greg 3.1 (*dp->ob_ohm)(o, dp->ob_priv);
813     cmtab[o+1] += cmtab[o];
814     }
815     cd->cTotal = cmtab[cd->calen];
816     scale = (double)maxval / cd->cTotal;
817     cd->carr[0] = 0;
818     for (o = 1; o < cd->calen; o++)
819     cd->carr[o] = scale*cmtab[o] + .5;
820     cd->carr[cd->calen] = maxval;
821     free(cmtab);
822     return 1;
823     }
824    
825     /* Get cumulative distribution for matrix BSDF */
826     static const SDCDst *
827     SDgetMtxCDist(const FVECT inVec, SDComponent *sdc)
828     {
829 greg 3.12 SDMat *dp;
830 greg 3.1 int reverse;
831     SDMatCDst myCD;
832     SDMatCDst *cd, *cdlast;
833 greg 3.5 /* check arguments */
834 greg 3.12 if ((inVec == NULL) | (sdc == NULL) ||
835     (dp = (SDMat *)sdc->dist) == NULL)
836 greg 3.1 return NULL;
837     memset(&myCD, 0, sizeof(myCD));
838     myCD.indx = mBSDF_incndx(dp, inVec);
839     if (myCD.indx >= 0) {
840     myCD.ob_priv = dp->ob_priv;
841     myCD.ob_vec = dp->ob_vec;
842     myCD.calen = dp->nout;
843     reverse = 0;
844 greg 3.19 } else { /* try reciprocity */
845 greg 3.1 myCD.indx = mBSDF_outndx(dp, inVec);
846     if (myCD.indx < 0)
847     return NULL;
848     myCD.ob_priv = dp->ib_priv;
849     myCD.ob_vec = dp->ib_vec;
850     myCD.calen = dp->ninc;
851     reverse = 1;
852     }
853     cdlast = NULL; /* check for it in cache list */
854 greg 3.34 /* PLACE MUTEX LOCK HERE FOR THREAD-SAFE */
855 greg 3.14 for (cd = (SDMatCDst *)sdc->cdList; cd != NULL;
856 greg 3.20 cdlast = cd, cd = cd->next)
857 greg 3.1 if (cd->indx == myCD.indx && (cd->calen == myCD.calen) &
858     (cd->ob_priv == myCD.ob_priv) &
859     (cd->ob_vec == myCD.ob_vec))
860     break;
861     if (cd == NULL) { /* need to allocate new entry */
862     cd = (SDMatCDst *)malloc(sizeof(SDMatCDst) +
863 greg 3.14 sizeof(myCD.carr[0])*myCD.calen);
864 greg 3.1 if (cd == NULL)
865     return NULL;
866     *cd = myCD; /* compute cumulative distribution */
867     if (!make_cdist(cd, inVec, dp, reverse)) {
868     free(cd);
869     return NULL;
870     }
871     cdlast = cd;
872     }
873     if (cdlast != NULL) { /* move entry to head of cache list */
874     cdlast->next = cd->next;
875 greg 3.20 cd->next = (SDMatCDst *)sdc->cdList;
876 greg 3.1 sdc->cdList = (SDCDst *)cd;
877     }
878 greg 3.34 /* END MUTEX LOCK */
879 greg 3.1 return (SDCDst *)cd; /* ready to go */
880     }
881    
882     /* Sample cumulative distribution */
883     static SDError
884 greg 3.12 SDsampMtxCDist(FVECT ioVec, double randX, const SDCDst *cdp)
885 greg 3.1 {
886     const unsigned maxval = ~0;
887     const SDMatCDst *mcd = (const SDMatCDst *)cdp;
888     const unsigned target = randX*maxval;
889     int i, iupper, ilower;
890 greg 3.5 /* check arguments */
891 greg 3.12 if ((ioVec == NULL) | (mcd == NULL))
892 greg 3.5 return SDEargument;
893 greg 3.1 /* binary search to find index */
894     ilower = 0; iupper = mcd->calen;
895     while ((i = (iupper + ilower) >> 1) != ilower)
896 greg 3.18 if (target >= mcd->carr[i])
897 greg 3.1 ilower = i;
898     else
899     iupper = i;
900     /* localize random position */
901     randX = (randX*maxval - mcd->carr[ilower]) /
902     (double)(mcd->carr[iupper] - mcd->carr[ilower]);
903     /* convert index to vector */
904 greg 3.12 if ((*mcd->ob_vec)(ioVec, i+randX, mcd->ob_priv))
905 greg 3.1 return SDEnone;
906 greg 3.12 strcpy(SDerrorDetail, "Matrix BSDF sampling fault");
907 greg 3.1 return SDEinternal;
908     }
909    
910     /* Fixed resolution BSDF methods */
911 greg 3.37 const SDFunc SDhandleMtx = {
912 greg 3.1 &SDgetMtxBSDF,
913     &SDqueryMtxProjSA,
914     &SDgetMtxCDist,
915     &SDsampMtxCDist,
916     &SDfreeMatrix,
917     };