ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/bsdf_m.c
Revision: 3.24
Committed: Mon Sep 10 18:06:08 2012 UTC (11 years, 7 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 3.23: +9 -9 lines
Log Message: 
Hopeful fixes to tensor tree reciprocity and orientation

File Contents

# User Rev Content
1 greg 3.2 #ifndef lint
2 greg 3.24 static const char RCSid[] = "$Id: bsdf_m.c,v 3.23 2012/09/02 15:33:15 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 <stdlib.h>
17     #include <math.h>
18     #include <ctype.h>
19     #include "ezxml.h"
20     #include "bsdf.h"
21     #include "bsdf_m.h"
22    
23     /* Function return codes */
24     #define RC_GOOD 1
25     #define RC_FAIL 0
26     #define RC_FORMERR (-1)
27     #define RC_DATERR (-2)
28     #define RC_UNSUPP (-3)
29     #define RC_INTERR (-4)
30     #define RC_MEMERR (-5)
31    
32     #define MAXLATS 46 /* maximum number of latitudes */
33    
34     /* BSDF angle specification */
35     typedef struct {
36     char name[64]; /* basis name */
37     int nangles; /* total number of directions */
38     struct {
39     float tmin; /* starting theta */
40     int nphis; /* number of phis (0 term) */
41     } lat[MAXLATS+1]; /* latitudes */
42     } ANGLE_BASIS;
43    
44     #define MAXABASES 7 /* limit on defined bases */
45    
46     static ANGLE_BASIS abase_list[MAXABASES] = {
47     {
48     "LBNL/Klems Full", 145,
49 greg 3.22 { {0., 1},
50 greg 3.1 {5., 8},
51     {15., 16},
52     {25., 20},
53     {35., 24},
54     {45., 24},
55     {55., 24},
56     {65., 16},
57     {75., 12},
58     {90., 0} }
59     }, {
60     "LBNL/Klems Half", 73,
61 greg 3.22 { {0., 1},
62 greg 3.1 {6.5, 8},
63     {19.5, 12},
64     {32.5, 16},
65     {46.5, 20},
66     {61.5, 12},
67     {76.5, 4},
68     {90., 0} }
69     }, {
70     "LBNL/Klems Quarter", 41,
71 greg 3.22 { {0., 1},
72 greg 3.1 {9., 8},
73     {27., 12},
74     {46., 12},
75     {66., 8},
76     {90., 0} }
77     }
78     };
79    
80     static int nabases = 3; /* current number of defined bases */
81    
82     static int
83     fequal(double a, double b)
84     {
85 greg 3.10 if (b != 0)
86 greg 3.1 a = a/b - 1.;
87     return (a <= 1e-6) & (a >= -1e-6);
88     }
89    
90 greg 3.9 /* Returns the given tag's character content or empty string if none */
91 greg 3.1 #ifdef ezxml_txt
92     #undef ezxml_txt
93     static char *
94     ezxml_txt(ezxml_t xml)
95     {
96     if (xml == NULL)
97     return "";
98     return xml->txt;
99     }
100     #endif
101    
102     /* Convert error to standard BSDF code */
103     static SDError
104     convert_errcode(int ec)
105     {
106     switch (ec) {
107     case RC_GOOD:
108     return SDEnone;
109     case RC_FORMERR:
110     return SDEformat;
111     case RC_DATERR:
112     return SDEdata;
113     case RC_UNSUPP:
114     return SDEsupport;
115     case RC_INTERR:
116     return SDEinternal;
117     case RC_MEMERR:
118     return SDEmemory;
119     }
120     return SDEunknown;
121     }
122    
123     /* Allocate a BSDF matrix of the given size */
124     static SDMat *
125     SDnewMatrix(int ni, int no)
126     {
127     SDMat *sm;
128    
129     if ((ni <= 0) | (no <= 0)) {
130     strcpy(SDerrorDetail, "Empty BSDF matrix request");
131     return NULL;
132     }
133     sm = (SDMat *)malloc(sizeof(SDMat) + (ni*no - 1)*sizeof(float));
134     if (sm == NULL) {
135     sprintf(SDerrorDetail, "Cannot allocate %dx%d BSDF matrix",
136     ni, no);
137     return NULL;
138     }
139     memset(sm, 0, sizeof(SDMat)-sizeof(float));
140     sm->ninc = ni;
141     sm->nout = no;
142    
143     return sm;
144     }
145    
146     /* Free a BSDF matrix */
147     #define SDfreeMatrix free
148    
149 greg 3.8 /* get vector for this angle basis index (front exiting) */
150 greg 3.1 static int
151 greg 3.9 fo_getvec(FVECT v, double ndxr, void *p)
152 greg 3.1 {
153 greg 3.9 ANGLE_BASIS *ab = (ANGLE_BASIS *)p;
154     int ndx = (int)ndxr;
155     double randX = ndxr - ndx;
156 greg 3.1 double rx[2];
157     int li;
158     double pol, azi, d;
159    
160 greg 3.9 if ((ndxr < 0) | (ndx >= ab->nangles))
161 greg 3.1 return RC_FAIL;
162     for (li = 0; ndx >= ab->lat[li].nphis; li++)
163     ndx -= ab->lat[li].nphis;
164     SDmultiSamp(rx, 2, randX);
165     pol = M_PI/180.*( (1.-rx[0])*ab->lat[li].tmin +
166     rx[0]*ab->lat[li+1].tmin );
167     azi = 2.*M_PI*(ndx + rx[1] - .5)/ab->lat[li].nphis;
168     v[2] = d = cos(pol);
169     d = sqrt(1. - d*d); /* sin(pol) */
170     v[0] = cos(azi)*d;
171     v[1] = sin(azi)*d;
172     return RC_GOOD;
173     }
174    
175 greg 3.8 /* get index corresponding to the given vector (front exiting) */
176 greg 3.1 static int
177 greg 3.8 fo_getndx(const FVECT v, void *p)
178 greg 3.1 {
179 greg 3.9 ANGLE_BASIS *ab = (ANGLE_BASIS *)p;
180 greg 3.1 int li, ndx;
181 greg 3.17 double pol, azi;
182 greg 3.1
183     if (v == NULL)
184     return -1;
185 greg 3.10 if ((v[2] < 0) | (v[2] > 1.))
186 greg 3.1 return -1;
187     pol = 180.0/M_PI*acos(v[2]);
188     azi = 180.0/M_PI*atan2(v[1], v[0]);
189     if (azi < 0.0) azi += 360.0;
190     for (li = 1; ab->lat[li].tmin <= pol; li++)
191     if (!ab->lat[li].nphis)
192     return -1;
193     --li;
194     ndx = (int)((1./360.)*azi*ab->lat[li].nphis + 0.5);
195     if (ndx >= ab->lat[li].nphis) ndx = 0;
196     while (li--)
197     ndx += ab->lat[li].nphis;
198     return ndx;
199     }
200    
201     /* compute square of real value */
202     static double sq(double x) { return x*x; }
203    
204 greg 3.8 /* get projected solid angle for this angle basis index (universal) */
205 greg 3.1 static double
206 greg 3.8 io_getohm(int ndx, void *p)
207 greg 3.1 {
208     static int last_li = -1;
209     static double last_ohm;
210     ANGLE_BASIS *ab = (ANGLE_BASIS *)p;
211     int li;
212     double theta, theta1;
213    
214     if ((ndx < 0) | (ndx >= ab->nangles))
215     return -1.;
216     for (li = 0; ndx >= ab->lat[li].nphis; li++)
217     ndx -= ab->lat[li].nphis;
218     if (li == last_li) /* cached latitude? */
219     return last_ohm;
220     last_li = li;
221 greg 3.21 theta = M_PI/180. * ab->lat[li].tmin;
222 greg 3.1 theta1 = M_PI/180. * ab->lat[li+1].tmin;
223     return last_ohm = M_PI*(sq(cos(theta)) - sq(cos(theta1))) /
224     (double)ab->lat[li].nphis;
225     }
226    
227 greg 3.8 /* get vector for this angle basis index (back incident) */
228 greg 3.1 static int
229 greg 3.9 bi_getvec(FVECT v, double ndxr, void *p)
230 greg 3.1 {
231 greg 3.9 if (!fo_getvec(v, ndxr, p))
232 greg 3.1 return RC_FAIL;
233    
234     v[0] = -v[0];
235     v[1] = -v[1];
236     v[2] = -v[2];
237    
238     return RC_GOOD;
239     }
240    
241 greg 3.8 /* get index corresponding to the vector (back incident) */
242 greg 3.1 static int
243 greg 3.8 bi_getndx(const FVECT v, void *p)
244 greg 3.1 {
245     FVECT v2;
246    
247     v2[0] = -v[0];
248     v2[1] = -v[1];
249     v2[2] = -v[2];
250    
251 greg 3.8 return fo_getndx(v2, p);
252     }
253    
254     /* get vector for this angle basis index (back exiting) */
255     static int
256 greg 3.9 bo_getvec(FVECT v, double ndxr, void *p)
257 greg 3.8 {
258 greg 3.9 if (!fo_getvec(v, ndxr, p))
259 greg 3.8 return RC_FAIL;
260    
261     v[2] = -v[2];
262    
263     return RC_GOOD;
264     }
265    
266     /* get index corresponding to the vector (back exiting) */
267     static int
268     bo_getndx(const FVECT v, void *p)
269     {
270     FVECT v2;
271    
272     v2[0] = v[0];
273     v2[1] = v[1];
274     v2[2] = -v[2];
275    
276     return fo_getndx(v2, p);
277     }
278    
279     /* get vector for this angle basis index (front incident) */
280     static int
281 greg 3.9 fi_getvec(FVECT v, double ndxr, void *p)
282 greg 3.8 {
283 greg 3.9 if (!fo_getvec(v, ndxr, p))
284 greg 3.8 return RC_FAIL;
285    
286     v[0] = -v[0];
287     v[1] = -v[1];
288    
289     return RC_GOOD;
290     }
291    
292     /* get index corresponding to the vector (front incident) */
293     static int
294     fi_getndx(const FVECT v, void *p)
295     {
296     FVECT v2;
297    
298     v2[0] = -v[0];
299     v2[1] = -v[1];
300     v2[2] = v[2];
301    
302     return fo_getndx(v2, p);
303 greg 3.1 }
304    
305     /* load custom BSDF angle basis */
306     static int
307     load_angle_basis(ezxml_t wab)
308     {
309     char *abname = ezxml_txt(ezxml_child(wab, "AngleBasisName"));
310     ezxml_t wbb;
311     int i;
312    
313     if (!abname || !*abname)
314     return RC_FAIL;
315     for (i = nabases; i--; )
316     if (!strcasecmp(abname, abase_list[i].name))
317     return RC_GOOD; /* assume it's the same */
318     if (nabases >= MAXABASES) {
319     sprintf(SDerrorDetail, "Out of angle bases reading '%s'",
320     abname);
321     return RC_INTERR;
322     }
323     strcpy(abase_list[nabases].name, abname);
324     abase_list[nabases].nangles = 0;
325     for (i = 0, wbb = ezxml_child(wab, "AngleBasisBlock");
326     wbb != NULL; i++, wbb = wbb->next) {
327     if (i >= MAXLATS) {
328     sprintf(SDerrorDetail, "Too many latitudes for '%s'",
329     abname);
330     return RC_INTERR;
331     }
332     abase_list[nabases].lat[i+1].tmin = atof(ezxml_txt(
333     ezxml_child(ezxml_child(wbb,
334     "ThetaBounds"), "UpperTheta")));
335     if (!i)
336 greg 3.21 abase_list[nabases].lat[0].tmin = 0;
337 greg 3.1 else if (!fequal(atof(ezxml_txt(ezxml_child(ezxml_child(wbb,
338     "ThetaBounds"), "LowerTheta"))),
339     abase_list[nabases].lat[i].tmin)) {
340     sprintf(SDerrorDetail, "Theta values disagree in '%s'",
341 greg 3.12 abname);
342 greg 3.1 return RC_DATERR;
343     }
344     abase_list[nabases].nangles +=
345     abase_list[nabases].lat[i].nphis =
346     atoi(ezxml_txt(ezxml_child(wbb, "nPhis")));
347     if (abase_list[nabases].lat[i].nphis <= 0 ||
348     (abase_list[nabases].lat[i].nphis == 1 &&
349     abase_list[nabases].lat[i].tmin > FTINY)) {
350     sprintf(SDerrorDetail, "Illegal phi count in '%s'",
351 greg 3.12 abname);
352 greg 3.1 return RC_DATERR;
353     }
354     }
355     abase_list[nabases++].lat[i].nphis = 0;
356     return RC_GOOD;
357     }
358    
359     /* compute min. proj. solid angle and max. direct hemispherical scattering */
360     static int
361     get_extrema(SDSpectralDF *df)
362     {
363     SDMat *dp = (SDMat *)df->comp[0].dist;
364     double *ohma;
365     int i, o;
366     /* initialize extrema */
367     df->minProjSA = M_PI;
368     df->maxHemi = .0;
369     ohma = (double *)malloc(dp->nout*sizeof(double));
370     if (ohma == NULL)
371     return RC_MEMERR;
372     /* get outgoing solid angles */
373     for (o = dp->nout; o--; )
374     if ((ohma[o] = mBSDF_outohm(dp,o)) < df->minProjSA)
375     df->minProjSA = ohma[o];
376     /* compute hemispherical sums */
377     for (i = dp->ninc; i--; ) {
378     double hemi = .0;
379     for (o = dp->nout; o--; )
380     hemi += ohma[o] * mBSDF_value(dp, i, o);
381     if (hemi > df->maxHemi)
382     df->maxHemi = hemi;
383     }
384     free(ohma);
385     /* need incoming solid angles, too? */
386 greg 3.5 if ((dp->ib_ohm != dp->ob_ohm) | (dp->ib_priv != dp->ob_priv)) {
387 greg 3.1 double ohm;
388     for (i = dp->ninc; i--; )
389     if ((ohm = mBSDF_incohm(dp,i)) < df->minProjSA)
390     df->minProjSA = ohm;
391     }
392     return (df->maxHemi <= 1.01);
393     }
394    
395     /* load BSDF distribution for this wavelength */
396     static int
397     load_bsdf_data(SDData *sd, ezxml_t wdb, int rowinc)
398     {
399     SDSpectralDF *df;
400     SDMat *dp;
401     char *sdata;
402     int inbi, outbi;
403     int i;
404     /* allocate BSDF component */
405     sdata = ezxml_txt(ezxml_child(wdb, "WavelengthDataDirection"));
406 greg 3.15 if (!sdata)
407     return RC_FAIL;
408 greg 3.8 /*
409     * Remember that front and back are reversed from WINDOW 6 orientations
410     */
411 greg 3.23 if (!strcasecmp(sdata, "Transmission Front")) {
412 greg 3.24 if (sd->tb != NULL)
413     SDfreeSpectralDF(sd->tb);
414     if ((sd->tb = SDnewSpectralDF(1)) == NULL)
415     return RC_MEMERR;
416     df = sd->tb;
417     } else if (!strcasecmp(sdata, "Transmission Back")) {
418 greg 3.1 if (sd->tf != NULL)
419     SDfreeSpectralDF(sd->tf);
420     if ((sd->tf = SDnewSpectralDF(1)) == NULL)
421     return RC_MEMERR;
422     df = sd->tf;
423     } else if (!strcasecmp(sdata, "Reflection Front")) {
424 greg 3.24 if (sd->rb != NULL)
425 greg 3.6 SDfreeSpectralDF(sd->rb);
426     if ((sd->rb = SDnewSpectralDF(1)) == NULL)
427     return RC_MEMERR;
428     df = sd->rb;
429     } else if (!strcasecmp(sdata, "Reflection Back")) {
430 greg 3.24 if (sd->rf != NULL)
431 greg 3.1 SDfreeSpectralDF(sd->rf);
432     if ((sd->rf = SDnewSpectralDF(1)) == NULL)
433     return RC_MEMERR;
434     df = sd->rf;
435     } else
436     return RC_FAIL;
437 greg 3.4 /* XXX should also check "ScatteringDataType" for consistency? */
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.1 df->comp[0].cspec[0] = c_dfcolor; /* XXX monochrome for now */
495     df->comp[0].dist = dp;
496     df->comp[0].func = &SDhandleMtx;
497     /* read BSDF data */
498 greg 3.15 sdata = ezxml_txt(ezxml_child(wdb, "ScatteringData"));
499 greg 3.1 if (!sdata || !*sdata) {
500     sprintf(SDerrorDetail, "Missing BSDF ScatteringData in '%s'",
501     sd->name);
502     return RC_FORMERR;
503     }
504     for (i = 0; i < dp->ninc*dp->nout; i++) {
505 greg 3.3 char *sdnext = fskip(sdata);
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     if (rowinc) {
516     int r = i/dp->nout;
517 greg 3.16 int c = i - r*dp->nout;
518 greg 3.1 mBSDF_value(dp,r,c) = atof(sdata);
519     } else
520     dp->bsdf[i] = atof(sdata);
521     sdata = sdnext;
522     }
523     return get_extrema(df);
524     }
525    
526     /* Subtract minimum (diffuse) scattering amount from BSDF */
527     static double
528     subtract_min(SDMat *sm)
529     {
530     float minv = sm->bsdf[0];
531     int n = sm->ninc*sm->nout;
532     int i;
533    
534     for (i = n; --i; )
535     if (sm->bsdf[i] < minv)
536     minv = sm->bsdf[i];
537 greg 3.15
538     if (minv <= FTINY)
539     return .0;
540    
541 greg 3.1 for (i = n; i--; )
542     sm->bsdf[i] -= minv;
543    
544     return minv*M_PI; /* be sure to include multiplier */
545     }
546    
547     /* Extract and separate diffuse portion of BSDF */
548     static void
549     extract_diffuse(SDValue *dv, SDSpectralDF *df)
550     {
551     int n;
552    
553     if (df == NULL || df->ncomp <= 0) {
554     dv->spec = c_dfcolor;
555     dv->cieY = .0;
556     return;
557     }
558     dv->spec = df->comp[0].cspec[0];
559     dv->cieY = subtract_min((SDMat *)df->comp[0].dist);
560     /* in case of multiple components */
561     for (n = df->ncomp; --n; ) {
562     double ymin = subtract_min((SDMat *)df->comp[n].dist);
563     c_cmix(&dv->spec, dv->cieY, &dv->spec, ymin, &df->comp[n].cspec[0]);
564     dv->cieY += ymin;
565     }
566 greg 3.15 df->maxHemi -= dv->cieY; /* adjust maximum hemispherical */
567 greg 3.4 /* make sure everything is set */
568 greg 3.1 c_ccvt(&dv->spec, C_CSXY+C_CSSPEC);
569     }
570    
571     /* Load a BSDF matrix from an open XML file */
572     SDError
573 greg 3.4 SDloadMtx(SDData *sd, ezxml_t wtl)
574 greg 3.1 {
575 greg 3.15 ezxml_t wld, wdb;
576     int rowIn;
577     char *txt;
578     int rval;
579     /* basic checks and data ordering */
580 greg 3.4 txt = ezxml_txt(ezxml_child(ezxml_child(wtl,
581     "DataDefinition"), "IncidentDataStructure"));
582     if (txt == NULL || !*txt) {
583 greg 3.1 sprintf(SDerrorDetail,
584 greg 3.4 "BSDF \"%s\": missing IncidentDataStructure",
585 greg 3.1 sd->name);
586     return SDEformat;
587     }
588     if (!strcasecmp(txt, "Rows"))
589     rowIn = 1;
590     else if (!strcasecmp(txt, "Columns"))
591     rowIn = 0;
592     else {
593     sprintf(SDerrorDetail,
594     "BSDF \"%s\": unsupported IncidentDataStructure",
595     sd->name);
596     return SDEsupport;
597     }
598 greg 3.18 /* get angle bases */
599     for (wld = ezxml_child(ezxml_child(wtl, "DataDefinition"), "AngleBasis");
600     wld != NULL; wld = wld->next) {
601     rval = load_angle_basis(wld);
602     if (rval < 0)
603     return convert_errcode(rval);
604     }
605 greg 3.15 /* load BSDF components */
606 greg 3.1 for (wld = ezxml_child(wtl, "WavelengthData");
607     wld != NULL; wld = wld->next) {
608     if (strcasecmp(ezxml_txt(ezxml_child(wld,"Wavelength")),
609     "Visible"))
610     continue; /* just visible for now */
611     for (wdb = ezxml_child(wld, "WavelengthDataBlock");
612     wdb != NULL; wdb = wdb->next)
613     if ((rval = load_bsdf_data(sd, wdb, rowIn)) < 0)
614 greg 3.4 return convert_errcode(rval);
615 greg 3.1 }
616 greg 3.15 /* separate diffuse components */
617 greg 3.1 extract_diffuse(&sd->rLambFront, sd->rf);
618     extract_diffuse(&sd->rLambBack, sd->rb);
619 greg 3.23 extract_diffuse(&sd->tLamb, (sd->tf != NULL) ? sd->tf : sd->tb);
620 greg 3.15 /* return success */
621 greg 3.1 return SDEnone;
622     }
623    
624     /* Get Matrix BSDF value */
625     static int
626     SDgetMtxBSDF(float coef[SDmaxCh], const FVECT outVec,
627 greg 3.12 const FVECT inVec, SDComponent *sdc)
628 greg 3.1 {
629 greg 3.12 const SDMat *dp;
630 greg 3.1 int i_ndx, o_ndx;
631 greg 3.12 /* check arguments */
632     if ((coef == NULL) | (outVec == NULL) | (inVec == NULL) | (sdc == NULL)
633     || (dp = (SDMat *)sdc->dist) == NULL)
634     return 0;
635 greg 3.1 /* get angle indices */
636     i_ndx = mBSDF_incndx(dp, inVec);
637     o_ndx = mBSDF_outndx(dp, outVec);
638     /* try reciprocity if necessary */
639 greg 3.19 if ((i_ndx < 0) & (o_ndx < 0)) {
640 greg 3.1 i_ndx = mBSDF_incndx(dp, outVec);
641     o_ndx = mBSDF_outndx(dp, inVec);
642     }
643     if ((i_ndx < 0) | (o_ndx < 0))
644     return 0; /* nothing from this component */
645     coef[0] = mBSDF_value(dp, i_ndx, o_ndx);
646     return 1; /* XXX monochrome for now */
647     }
648    
649 greg 3.12 /* Query solid angle for vector(s) */
650 greg 3.1 static SDError
651 greg 3.10 SDqueryMtxProjSA(double *psa, const FVECT v1, const RREAL *v2,
652 greg 3.12 int qflags, SDComponent *sdc)
653 greg 3.1 {
654 greg 3.12 const SDMat *dp;
655 greg 3.5 double inc_psa, out_psa;
656     /* check arguments */
657 greg 3.12 if ((psa == NULL) | (v1 == NULL) | (sdc == NULL) ||
658     (dp = (SDMat *)sdc->dist) == NULL)
659 greg 3.1 return SDEargument;
660 greg 3.10 if (v2 == NULL)
661     v2 = v1;
662 greg 3.5 /* get projected solid angles */
663 greg 3.10 out_psa = mBSDF_outohm(dp, mBSDF_outndx(dp, v1));
664     inc_psa = mBSDF_incohm(dp, mBSDF_incndx(dp, v2));
665 greg 3.12 if ((v1 != v2) & (out_psa <= 0) & (inc_psa <= 0)) {
666 greg 3.11 inc_psa = mBSDF_outohm(dp, mBSDF_outndx(dp, v2));
667     out_psa = mBSDF_incohm(dp, mBSDF_incndx(dp, v1));
668     }
669 greg 3.5
670     switch (qflags) { /* record based on flag settings */
671     case SDqueryMax:
672     if (inc_psa > psa[0])
673     psa[0] = inc_psa;
674     if (out_psa > psa[0])
675     psa[0] = out_psa;
676     break;
677     case SDqueryMin+SDqueryMax:
678 greg 3.13 if (inc_psa > psa[1])
679 greg 3.5 psa[1] = inc_psa;
680 greg 3.13 if (out_psa > psa[1])
681 greg 3.5 psa[1] = out_psa;
682     /* fall through */
683 greg 3.12 case SDqueryVal:
684     if (qflags == SDqueryVal)
685     psa[0] = M_PI;
686 greg 3.14 /* fall through */
687     case SDqueryMin:
688 greg 3.10 if ((inc_psa > 0) & (inc_psa < psa[0]))
689 greg 3.1 psa[0] = inc_psa;
690 greg 3.10 if ((out_psa > 0) & (out_psa < psa[0]))
691 greg 3.5 psa[0] = out_psa;
692     break;
693 greg 3.1 }
694 greg 3.5 /* make sure it's legal */
695 greg 3.10 return (psa[0] <= 0) ? SDEinternal : SDEnone;
696 greg 3.1 }
697    
698     /* Compute new cumulative distribution from BSDF */
699     static int
700     make_cdist(SDMatCDst *cd, const FVECT inVec, SDMat *dp, int rev)
701     {
702     const unsigned maxval = ~0;
703     double *cmtab, scale;
704     int o;
705    
706     cmtab = (double *)malloc((cd->calen+1)*sizeof(double));
707     if (cmtab == NULL)
708     return 0;
709     cmtab[0] = .0;
710     for (o = 0; o < cd->calen; o++) {
711     if (rev)
712     cmtab[o+1] = mBSDF_value(dp, o, cd->indx) *
713     (*dp->ib_ohm)(o, dp->ib_priv);
714     else
715     cmtab[o+1] = mBSDF_value(dp, cd->indx, o) *
716     (*dp->ob_ohm)(o, dp->ob_priv);
717     cmtab[o+1] += cmtab[o];
718     }
719     cd->cTotal = cmtab[cd->calen];
720     scale = (double)maxval / cd->cTotal;
721     cd->carr[0] = 0;
722     for (o = 1; o < cd->calen; o++)
723     cd->carr[o] = scale*cmtab[o] + .5;
724     cd->carr[cd->calen] = maxval;
725     free(cmtab);
726     return 1;
727     }
728    
729     /* Get cumulative distribution for matrix BSDF */
730     static const SDCDst *
731     SDgetMtxCDist(const FVECT inVec, SDComponent *sdc)
732     {
733 greg 3.12 SDMat *dp;
734 greg 3.1 int reverse;
735     SDMatCDst myCD;
736     SDMatCDst *cd, *cdlast;
737 greg 3.5 /* check arguments */
738 greg 3.12 if ((inVec == NULL) | (sdc == NULL) ||
739     (dp = (SDMat *)sdc->dist) == NULL)
740 greg 3.1 return NULL;
741     memset(&myCD, 0, sizeof(myCD));
742     myCD.indx = mBSDF_incndx(dp, inVec);
743     if (myCD.indx >= 0) {
744     myCD.ob_priv = dp->ob_priv;
745     myCD.ob_vec = dp->ob_vec;
746     myCD.calen = dp->nout;
747     reverse = 0;
748 greg 3.19 } else { /* try reciprocity */
749 greg 3.1 myCD.indx = mBSDF_outndx(dp, inVec);
750     if (myCD.indx < 0)
751     return NULL;
752     myCD.ob_priv = dp->ib_priv;
753     myCD.ob_vec = dp->ib_vec;
754     myCD.calen = dp->ninc;
755     reverse = 1;
756     }
757     cdlast = NULL; /* check for it in cache list */
758 greg 3.14 for (cd = (SDMatCDst *)sdc->cdList; cd != NULL;
759 greg 3.20 cdlast = cd, cd = cd->next)
760 greg 3.1 if (cd->indx == myCD.indx && (cd->calen == myCD.calen) &
761     (cd->ob_priv == myCD.ob_priv) &
762     (cd->ob_vec == myCD.ob_vec))
763     break;
764     if (cd == NULL) { /* need to allocate new entry */
765     cd = (SDMatCDst *)malloc(sizeof(SDMatCDst) +
766 greg 3.14 sizeof(myCD.carr[0])*myCD.calen);
767 greg 3.1 if (cd == NULL)
768     return NULL;
769     *cd = myCD; /* compute cumulative distribution */
770     if (!make_cdist(cd, inVec, dp, reverse)) {
771     free(cd);
772     return NULL;
773     }
774     cdlast = cd;
775     }
776     if (cdlast != NULL) { /* move entry to head of cache list */
777     cdlast->next = cd->next;
778 greg 3.20 cd->next = (SDMatCDst *)sdc->cdList;
779 greg 3.1 sdc->cdList = (SDCDst *)cd;
780     }
781     return (SDCDst *)cd; /* ready to go */
782     }
783    
784     /* Sample cumulative distribution */
785     static SDError
786 greg 3.12 SDsampMtxCDist(FVECT ioVec, double randX, const SDCDst *cdp)
787 greg 3.1 {
788     const unsigned maxval = ~0;
789     const SDMatCDst *mcd = (const SDMatCDst *)cdp;
790     const unsigned target = randX*maxval;
791     int i, iupper, ilower;
792 greg 3.5 /* check arguments */
793 greg 3.12 if ((ioVec == NULL) | (mcd == NULL))
794 greg 3.5 return SDEargument;
795 greg 3.1 /* binary search to find index */
796     ilower = 0; iupper = mcd->calen;
797     while ((i = (iupper + ilower) >> 1) != ilower)
798 greg 3.18 if (target >= mcd->carr[i])
799 greg 3.1 ilower = i;
800     else
801     iupper = i;
802     /* localize random position */
803     randX = (randX*maxval - mcd->carr[ilower]) /
804     (double)(mcd->carr[iupper] - mcd->carr[ilower]);
805     /* convert index to vector */
806 greg 3.12 if ((*mcd->ob_vec)(ioVec, i+randX, mcd->ob_priv))
807 greg 3.1 return SDEnone;
808 greg 3.12 strcpy(SDerrorDetail, "Matrix BSDF sampling fault");
809 greg 3.1 return SDEinternal;
810     }
811    
812     /* Fixed resolution BSDF methods */
813     SDFunc SDhandleMtx = {
814     &SDgetMtxBSDF,
815     &SDqueryMtxProjSA,
816     &SDgetMtxCDist,
817     &SDsampMtxCDist,
818     &SDfreeMatrix,
819     };