ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/src/common/bsdf_m.c
Revision: 3.29
Committed: Thu Apr 2 16:40:32 2015 UTC (10 years, 7 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 3.28: +172 -70 lines
Log Message: 
Support for color BSDF rendering using Klems representation

File Contents

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