ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/bsdf_t.c
Revision: 3.46
Committed: Tue Apr 17 18:11:16 2018 UTC (6 years, 6 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 3.45: +7 -6 lines
Log Message:
Fixed issues with zero transmission values in chromatic tensor trees

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: bsdf_t.c,v 3.45 2018/01/05 21:00:24 greg Exp $";
3 #endif
4 /*
5 * bsdf_t.c
6 *
7 * Definitions for variable-resolution BSDF trees
8 *
9 * Created by Greg Ward on 2/2/11.
10 *
11 */
12
13 #define _USE_MATH_DEFINES
14 #include "rtio.h"
15 #include <stdlib.h>
16 #include <math.h>
17 #include <ctype.h>
18 #include "ezxml.h"
19 #include "bsdf.h"
20 #include "bsdf_t.h"
21 #include "hilbert.h"
22
23 /* Callback function type for SDtraverseTre() */
24 typedef int SDtreCallback(float val, const double *cmin, double csiz,
25 void *cptr);
26 /* reference width maximum (1.0) */
27 static const unsigned iwbits = sizeof(unsigned)*4;
28 static const unsigned iwmax = 1<<(sizeof(unsigned)*4);
29 /* maximum cumulative value */
30 static const unsigned cumlmax = ~0;
31 /* constant z-vector */
32 static const FVECT zvec = {.0, .0, 1.};
33 /* quantization value */
34 static double quantum = 1./256.;
35 /* our RGB primaries */
36 static C_COLOR tt_RGB_prim[3];
37 static float tt_RGB_coef[3];
38
39 static const double czero[SD_MAXDIM];
40
41 enum {tt_Y, tt_u, tt_v}; /* tree components (tt_Y==0) */
42
43 /* Struct used for our distribution-building callback */
44 typedef struct {
45 short nic; /* number of input coordinates */
46 short rev; /* reversing query */
47 unsigned alen; /* current array length */
48 unsigned nall; /* number of allocated entries */
49 unsigned wmin; /* minimum square size so far */
50 unsigned wmax; /* maximum square size */
51 struct outdir_s {
52 unsigned hent; /* entering Hilbert index */
53 int wid; /* this square size */
54 float bsdf; /* BSDF for this square */
55 } *darr; /* output direction array */
56 } SDdistScaffold;
57
58 /* Allocate a new scattering distribution node */
59 static SDNode *
60 SDnewNode(int nd, int lg)
61 {
62 SDNode *st;
63
64 if (nd <= 0) {
65 strcpy(SDerrorDetail, "Zero dimension BSDF node request");
66 return NULL;
67 }
68 if (nd > SD_MAXDIM) {
69 sprintf(SDerrorDetail, "Illegal BSDF dimension (%d > %d)",
70 nd, SD_MAXDIM);
71 return NULL;
72 }
73 if (lg < 0) {
74 st = (SDNode *)malloc(sizeof(SDNode) +
75 sizeof(st->u.t[0])*((1<<nd) - 1));
76 if (st == NULL) {
77 sprintf(SDerrorDetail,
78 "Cannot allocate %d branch BSDF tree", 1<<nd);
79 return NULL;
80 }
81 memset(st->u.t, 0, sizeof(st->u.t[0])<<nd);
82 } else {
83 st = (SDNode *)malloc(sizeof(SDNode) +
84 sizeof(st->u.v[0])*((1 << nd*lg) - 1));
85 if (st == NULL) {
86 sprintf(SDerrorDetail,
87 "Cannot allocate %d BSDF leaves", 1 << nd*lg);
88 return NULL;
89 }
90 }
91 st->ndim = nd;
92 st->log2GR = lg;
93 return st;
94 }
95
96 /* Free an SD tree */
97 static void
98 SDfreeTre(SDNode *st)
99 {
100 int n;
101
102 if (st == NULL)
103 return;
104 for (n = (st->log2GR < 0) << st->ndim; n--; )
105 SDfreeTre(st->u.t[n]);
106 free(st);
107 }
108
109 /* Free a variable-resolution BSDF */
110 static void
111 SDFreeBTre(void *p)
112 {
113 SDTre *sdt = (SDTre *)p;
114
115 if (sdt == NULL)
116 return;
117 SDfreeTre(sdt->stc[tt_Y]);
118 SDfreeTre(sdt->stc[tt_u]);
119 SDfreeTre(sdt->stc[tt_v]);
120 free(sdt);
121 }
122
123 /* Fill branch's worth of grid values from subtree */
124 static void
125 fill_grid_branch(float *dptr, const float *sptr, int nd, int shft)
126 {
127 unsigned n = 1 << (shft-1);
128
129 if (!--nd) { /* end on the line */
130 memcpy(dptr, sptr, sizeof(*dptr)*n);
131 return;
132 }
133 while (n--) /* recurse on each slice */
134 fill_grid_branch(dptr + (n << shft*nd),
135 sptr + (n << (shft-1)*nd), nd, shft);
136 }
137
138 /* Get pointer at appropriate offset for the given branch */
139 static float *
140 grid_branch_start(SDNode *st, int n)
141 {
142 unsigned skipsiz = 1 << (st->log2GR - 1);
143 float *vptr = st->u.v;
144 int i;
145
146 for (i = st->ndim; i--; skipsiz <<= st->log2GR)
147 if (1<<i & n)
148 vptr += skipsiz;
149 return vptr;
150 }
151
152 /* Simplify (consolidate) a tree by flattening uniform depth regions */
153 static SDNode *
154 SDsimplifyTre(SDNode *st)
155 {
156 int match, n;
157
158 if (st == NULL) /* check for invalid tree */
159 return NULL;
160 if (st->log2GR >= 0) /* grid just returns unaltered */
161 return st;
162 match = 1; /* check if grids below match */
163 for (n = 0; n < 1<<st->ndim; n++) {
164 if ((st->u.t[n] = SDsimplifyTre(st->u.t[n])) == NULL)
165 return NULL; /* propogate error up call stack */
166 match &= (st->u.t[n]->log2GR == st->u.t[0]->log2GR);
167 }
168 if (match && (match = st->u.t[0]->log2GR) >= 0) {
169 SDNode *stn = SDnewNode(st->ndim, match + 1);
170 if (stn == NULL) /* out of memory? */
171 return st;
172 /* transfer values to new grid */
173 for (n = 1 << st->ndim; n--; )
174 fill_grid_branch(grid_branch_start(stn, n),
175 st->u.t[n]->u.v, stn->ndim, stn->log2GR);
176 SDfreeTre(st); /* free old tree */
177 st = stn; /* return new one */
178 }
179 return st;
180 }
181
182 /* Assign the given voxel in tree (produces no grid nodes) */
183 static SDNode *
184 SDsetVoxel(SDNode *sroot, int nd, const double *tmin, const double tsiz, float val)
185 {
186 double ctrk[SD_MAXDIM];
187 double csiz = 1.;
188 SDNode *st;
189 int i, n;
190 /* check arguments */
191 for (i = nd; i-- > 0; )
192 if ((tmin[i] < .0) | (tmin[i] >= 1.-FTINY))
193 break;
194 if ((i >= 0) | (nd <= 0) | (tsiz <= FTINY) | (tsiz > 1.+FTINY) |
195 (sroot != NULL && sroot->ndim != nd)) {
196 SDfreeTre(sroot);
197 return NULL;
198 }
199 if (tsiz >= 1.-FTINY) { /* special case when tree is a leaf */
200 SDfreeTre(sroot);
201 if ((sroot = SDnewNode(nd, 0)) != NULL)
202 sroot->u.v[0] = val;
203 return sroot;
204 }
205 /* make sure we have branching root */
206 if (sroot != NULL && sroot->log2GR >= 0) {
207 SDfreeTre(sroot); sroot = NULL;
208 }
209 if (sroot == NULL && (sroot = SDnewNode(nd, -1)) == NULL)
210 return NULL;
211 st = sroot; /* climb/grow tree */
212 memset(ctrk, 0, sizeof(ctrk));
213 for ( ; ; ) {
214 csiz *= .5; /* find appropriate branch */
215 n = 0;
216 for (i = nd; i--; )
217 if (ctrk[i]+csiz <= tmin[i]+FTINY) {
218 ctrk[i] += csiz;
219 n |= 1 << i;
220 }
221 /* reached desired voxel? */
222 if (csiz <= tsiz+FTINY) {
223 SDfreeTre(st->u.t[n]);
224 st = st->u.t[n] = SDnewNode(nd, 0);
225 break;
226 }
227 /* else grow tree as needed */
228 if (st->u.t[n] != NULL && st->u.t[n]->log2GR >= 0) {
229 SDfreeTre(st->u.t[n]); st->u.t[n] = NULL;
230 }
231 if (st->u.t[n] == NULL)
232 st->u.t[n] = SDnewNode(nd, -1);
233 if ((st = st->u.t[n]) == NULL)
234 break;
235 }
236 if (st == NULL) {
237 SDfreeTre(sroot);
238 return NULL;
239 }
240 st->u.v[0] = val; /* assign leaf and return root */
241 return sroot;
242 }
243
244 /* Find smallest leaf in tree */
245 static double
246 SDsmallestLeaf(const SDNode *st)
247 {
248 if (st->log2GR < 0) { /* tree branches */
249 double lmin = 1.;
250 int n;
251 for (n = 1<<st->ndim; n--; ) {
252 double lsiz = SDsmallestLeaf(st->u.t[n]);
253 if (lsiz < lmin)
254 lmin = lsiz;
255 }
256 return .5*lmin;
257 }
258 /* leaf grid width */
259 return 1. / (double)(1 << st->log2GR);
260 }
261
262 /* Add up N-dimensional hypercube array values over the given box */
263 static double
264 SDiterSum(const float *va, int nd, int shft, const int *imin, const int *imax)
265 {
266 const unsigned skipsiz = 1 << --nd*shft;
267 double sum = .0;
268 int i;
269
270 va += *imin * skipsiz;
271
272 if (skipsiz == 1)
273 for (i = *imin; i < *imax; i++)
274 sum += *va++;
275 else
276 for (i = *imin; i < *imax; i++, va += skipsiz)
277 sum += SDiterSum(va, nd, shft, imin+1, imax+1);
278 return sum;
279 }
280
281 /* Average BSDF leaves over an orthotope defined by the unit hypercube */
282 static double
283 SDavgTreBox(const SDNode *st, const double *bmin, const double *bmax)
284 {
285 unsigned n;
286 int i;
287
288 if (!st)
289 return .0;
290 /* check box limits */
291 for (i = st->ndim; i--; ) {
292 if (bmin[i] >= 1.)
293 return .0;
294 if (bmax[i] <= 0)
295 return .0;
296 if (bmin[i] >= bmax[i])
297 return .0;
298 }
299 if (st->log2GR < 0) { /* iterate on subtree */
300 double sum = .0, wsum = 1e-20;
301 double sbmin[SD_MAXDIM], sbmax[SD_MAXDIM], w;
302 for (n = 1 << st->ndim; n--; ) {
303 w = 1.;
304 for (i = st->ndim; i--; ) {
305 sbmin[i] = 2.*bmin[i];
306 sbmax[i] = 2.*bmax[i];
307 if (n & 1<<i) {
308 sbmin[i] -= 1.;
309 sbmax[i] -= 1.;
310 }
311 if (sbmin[i] < .0) sbmin[i] = .0;
312 if (sbmax[i] > 1.) sbmax[i] = 1.;
313 if (sbmin[i] >= sbmax[i]) {
314 w = .0;
315 break;
316 }
317 w *= sbmax[i] - sbmin[i];
318 }
319 if (w > 1e-10) {
320 sum += w * SDavgTreBox(st->u.t[n], sbmin, sbmax);
321 wsum += w;
322 }
323 }
324 return sum / wsum;
325 } else { /* iterate over leaves */
326 int imin[SD_MAXDIM], imax[SD_MAXDIM];
327
328 n = 1;
329 for (i = st->ndim; i--; ) {
330 imin[i] = (bmin[i] <= 0) ? 0 :
331 (int)((1 << st->log2GR)*bmin[i]);
332 imax[i] = (bmax[i] >= 1.) ? (1 << st->log2GR) :
333 (int)((1 << st->log2GR)*bmax[i] + .999999);
334 n *= imax[i] - imin[i];
335 }
336 if (n)
337 return SDiterSum(st->u.v, st->ndim,
338 st->log2GR, imin, imax) / (double)n;
339 }
340 return .0;
341 }
342
343 /* Recursive call for SDtraverseTre() */
344 static int
345 SDdotravTre(const SDNode *st, const double *pos, int cmask,
346 SDtreCallback *cf, void *cptr,
347 const double *cmin, double csiz)
348 {
349 int rv, rval = 0;
350 double bmin[SD_MAXDIM];
351 int i, n;
352 /* paranoia */
353 if (st == NULL)
354 return 0;
355 /* in branches? */
356 if (st->log2GR < 0) {
357 unsigned skipmask = 0;
358 csiz *= .5;
359 for (i = st->ndim; i--; )
360 if (1<<i & cmask) {
361 if (pos[i] < cmin[i] + csiz)
362 for (n = 1 << st->ndim; n--; ) {
363 if (n & 1<<i)
364 skipmask |= 1<<n;
365 }
366 else
367 for (n = 1 << st->ndim; n--; ) {
368 if (!(n & 1<<i))
369 skipmask |= 1<<n;
370 }
371 }
372 for (n = 1 << st->ndim; n--; ) {
373 if (1<<n & skipmask)
374 continue;
375 for (i = st->ndim; i--; )
376 if (1<<i & n)
377 bmin[i] = cmin[i] + csiz;
378 else
379 bmin[i] = cmin[i];
380
381 rval += rv = SDdotravTre(st->u.t[n], pos, cmask,
382 cf, cptr, bmin, csiz);
383 if (rv < 0)
384 return rv;
385 }
386 } else { /* else traverse leaves */
387 int clim[SD_MAXDIM][2];
388 int cpos[SD_MAXDIM];
389
390 if (st->log2GR == 0) /* short cut */
391 return (*cf)(st->u.v[0], cmin, csiz, cptr);
392
393 csiz /= (double)(1 << st->log2GR);
394 /* assign coord. ranges */
395 for (i = st->ndim; i--; )
396 if (1<<i & cmask) {
397 clim[i][0] = (pos[i] - cmin[i])/csiz;
398 /* check overflow from f.p. error */
399 clim[i][0] -= clim[i][0] >> st->log2GR;
400 clim[i][1] = clim[i][0] + 1;
401 } else {
402 clim[i][0] = 0;
403 clim[i][1] = 1 << st->log2GR;
404 }
405 #if (SD_MAXDIM == 4)
406 bmin[0] = cmin[0] + csiz*clim[0][0];
407 for (cpos[0] = clim[0][0]; cpos[0] < clim[0][1]; cpos[0]++) {
408 bmin[1] = cmin[1] + csiz*clim[1][0];
409 for (cpos[1] = clim[1][0]; cpos[1] < clim[1][1]; cpos[1]++) {
410 bmin[2] = cmin[2] + csiz*clim[2][0];
411 if (st->ndim == 3) {
412 cpos[2] = clim[2][0];
413 n = cpos[0];
414 for (i = 1; i < 3; i++)
415 n = (n << st->log2GR) + cpos[i];
416 for ( ; cpos[2] < clim[2][1]; cpos[2]++) {
417 rval += rv = (*cf)(st->u.v[n++], bmin, csiz, cptr);
418 if (rv < 0)
419 return rv;
420 bmin[2] += csiz;
421 }
422 } else {
423 for (cpos[2] = clim[2][0]; cpos[2] < clim[2][1]; cpos[2]++) {
424 bmin[3] = cmin[3] + csiz*(cpos[3] = clim[3][0]);
425 n = cpos[0];
426 for (i = 1; i < 4; i++)
427 n = (n << st->log2GR) + cpos[i];
428 for ( ; cpos[3] < clim[3][1]; cpos[3]++) {
429 rval += rv = (*cf)(st->u.v[n++], bmin, csiz, cptr);
430 if (rv < 0)
431 return rv;
432 bmin[3] += csiz;
433 }
434 bmin[2] += csiz;
435 }
436 }
437 bmin[1] += csiz;
438 }
439 bmin[0] += csiz;
440 }
441 #else
442 _!_ "broken code segment!"
443 #endif
444 }
445 return rval;
446 }
447
448 /* Traverse a tree, visiting nodes in a slice that fits partial position */
449 static int
450 SDtraverseTre(const SDNode *st, const double *pos, int cmask,
451 SDtreCallback *cf, void *cptr)
452 {
453 int i;
454 /* check arguments */
455 if ((st == NULL) | (cf == NULL))
456 return -1;
457 for (i = st->ndim; i--; )
458 if (1<<i & cmask && (pos[i] < 0) | (pos[i] >= 1.))
459 return -1;
460
461 return SDdotravTre(st, pos, cmask, cf, cptr, czero, 1.);
462 }
463
464 /* Look up tree value at the given grid position */
465 static float
466 SDlookupTre(const SDNode *st, const double *pos, double *hcube)
467 {
468 double spos[SD_MAXDIM];
469 int i, n, t;
470 /* initialize voxel return */
471 if (hcube) {
472 hcube[i = st->ndim] = 1.;
473 while (i--)
474 hcube[i] = .0;
475 }
476 /* climb the tree */
477 while (st != NULL && st->log2GR < 0) {
478 n = 0; /* move to appropriate branch */
479 if (hcube) hcube[st->ndim] *= .5;
480 for (i = st->ndim; i--; ) {
481 spos[i] = 2.*pos[i];
482 t = (spos[i] >= 1.);
483 n |= t<<i;
484 spos[i] -= (double)t;
485 if (hcube) hcube[i] += (double)t * hcube[st->ndim];
486 }
487 st = st->u.t[n]; /* avoids tail recursion */
488 pos = spos;
489 }
490 if (st == NULL) /* should never happen? */
491 return .0;
492 if (st->log2GR == 0) /* short cut */
493 return st->u.v[0];
494 n = t = 0; /* find grid array index */
495 for (i = st->ndim; i--; ) {
496 n += (int)((1<<st->log2GR)*pos[i]) << t;
497 t += st->log2GR;
498 }
499 if (hcube) { /* compute final hypercube */
500 hcube[st->ndim] /= (double)(1<<st->log2GR);
501 for (i = st->ndim; i--; )
502 hcube[i] += floor((1<<st->log2GR)*pos[i])*hcube[st->ndim];
503 }
504 return st->u.v[n]; /* no interpolation */
505 }
506
507 /* Convert CIE (Y,u',v') color to our RGB */
508 static void
509 SDyuv2rgb(double yval, double uprime, double vprime, float rgb[3])
510 {
511 const double dfact = 1./(6.*uprime - 16.*vprime + 12.);
512 C_COLOR cxy;
513
514 c_cset(&cxy, 9.*uprime*dfact, 4.*vprime*dfact);
515 c_toSharpRGB(&cxy, yval, rgb);
516 }
517
518 /* Query BSDF value and sample hypercube for the given vectors */
519 static int
520 SDqueryTre(const SDTre *sdt, float *coef,
521 const FVECT outVec, const FVECT inVec, double *hc)
522 {
523 const RREAL *vtmp;
524 float yval;
525 FVECT rOutVec;
526 double gridPos[4];
527
528 if (sdt->stc[tt_Y] == NULL) /* paranoia, I hope */
529 return 0;
530
531 switch (sdt->sidef) { /* whose side are you on? */
532 case SD_FREFL:
533 if ((outVec[2] < 0) | (inVec[2] < 0))
534 return 0;
535 break;
536 case SD_BREFL:
537 if ((outVec[2] > 0) | (inVec[2] > 0))
538 return 0;
539 break;
540 case SD_FXMIT:
541 if (outVec[2] > 0) {
542 if (inVec[2] > 0)
543 return 0;
544 vtmp = outVec; outVec = inVec; inVec = vtmp;
545 } else if (inVec[2] < 0)
546 return 0;
547 break;
548 case SD_BXMIT:
549 if (inVec[2] > 0) {
550 if (outVec[2] > 0)
551 return 0;
552 vtmp = outVec; outVec = inVec; inVec = vtmp;
553 } else if (outVec[2] < 0)
554 return 0;
555 break;
556 default:
557 return 0;
558 }
559 /* convert vector coordinates */
560 if (sdt->stc[tt_Y]->ndim == 3) {
561 spinvector(rOutVec, outVec, zvec, -atan2(-inVec[1],-inVec[0]));
562 gridPos[0] = (.5-FTINY) -
563 .5*sqrt(inVec[0]*inVec[0] + inVec[1]*inVec[1]);
564 SDdisk2square(gridPos+1, rOutVec[0], rOutVec[1]);
565 } else if (sdt->stc[tt_Y]->ndim == 4) {
566 SDdisk2square(gridPos, -inVec[0], -inVec[1]);
567 SDdisk2square(gridPos+2, outVec[0], outVec[1]);
568 } else
569 return 0; /* should be internal error */
570 /* get BSDF value */
571 yval = SDlookupTre(sdt->stc[tt_Y], gridPos, hc);
572 if (coef == NULL) /* just getting hypercube? */
573 return 1;
574 if (sdt->stc[tt_u] == NULL || sdt->stc[tt_v] == NULL) {
575 *coef = yval;
576 return 1; /* no color */
577 }
578 /* else decode color */
579 SDyuv2rgb(yval, SDlookupTre(sdt->stc[tt_u], gridPos, NULL),
580 SDlookupTre(sdt->stc[tt_v], gridPos, NULL), coef);
581 coef[0] *= tt_RGB_coef[0];
582 coef[1] *= tt_RGB_coef[1];
583 coef[2] *= tt_RGB_coef[2];
584 return 3;
585 }
586
587 /* Compute non-diffuse component for variable-resolution BSDF */
588 static int
589 SDgetTreBSDF(float coef[SDmaxCh], const FVECT outVec,
590 const FVECT inVec, SDComponent *sdc)
591 {
592 /* check arguments */
593 if ((coef == NULL) | (outVec == NULL) | (inVec == NULL) | (sdc == NULL)
594 || sdc->dist == NULL)
595 return 0;
596 /* get nearest BSDF value */
597 return SDqueryTre((SDTre *)sdc->dist, coef, outVec, inVec, NULL);
598 }
599
600 /* Callback to build cumulative distribution using SDtraverseTre() */
601 static int
602 build_scaffold(float val, const double *cmin, double csiz, void *cptr)
603 {
604 SDdistScaffold *sp = (SDdistScaffold *)cptr;
605 int wid = csiz*(double)iwmax + .5;
606 double revcmin[2];
607 bitmask_t bmin[2], bmax[2];
608
609 if (sp->rev) { /* need to reverse sense? */
610 revcmin[0] = 1. - cmin[0] - csiz;
611 revcmin[1] = 1. - cmin[1] - csiz;
612 cmin = revcmin;
613 } else {
614 cmin += sp->nic; /* else skip to output coords */
615 }
616 if (wid < sp->wmin) /* new minimum width? */
617 sp->wmin = wid;
618 if (wid > sp->wmax) /* new maximum? */
619 sp->wmax = wid;
620 if (sp->alen >= sp->nall) { /* need more space? */
621 struct outdir_s *ndarr;
622 sp->nall = (int)(1.5*sp->nall) + 256;
623 ndarr = (struct outdir_s *)realloc(sp->darr,
624 sizeof(struct outdir_s)*sp->nall);
625 if (ndarr == NULL) {
626 sprintf(SDerrorDetail,
627 "Cannot grow scaffold to %u entries", sp->nall);
628 return -1; /* abort build */
629 }
630 sp->darr = ndarr;
631 }
632 /* find Hilbert entry index */
633 bmin[0] = cmin[0]*(double)iwmax + .5;
634 bmin[1] = cmin[1]*(double)iwmax + .5;
635 bmax[0] = bmin[0] + wid-1;
636 bmax[1] = bmin[1] + wid-1;
637 hilbert_box_vtx(2, sizeof(bitmask_t), iwbits, 1, bmin, bmax);
638 sp->darr[sp->alen].hent = hilbert_c2i(2, iwbits, bmin);
639 sp->darr[sp->alen].wid = wid;
640 sp->darr[sp->alen].bsdf = val;
641 sp->alen++; /* on to the next entry */
642 return 0;
643 }
644
645 /* Scaffold comparison function for qsort -- ascending Hilbert index */
646 static int
647 sscmp(const void *p1, const void *p2)
648 {
649 unsigned h1 = (*(const struct outdir_s *)p1).hent;
650 unsigned h2 = (*(const struct outdir_s *)p2).hent;
651
652 if (h1 > h2)
653 return 1;
654 if (h1 < h2)
655 return -1;
656 return 0;
657 }
658
659 /* Create a new cumulative distribution for the given input direction */
660 static SDTreCDst *
661 make_cdist(const SDTre *sdt, const double *invec, int rev)
662 {
663 SDdistScaffold myScaffold;
664 double pos[4];
665 int cmask;
666 SDTreCDst *cd;
667 struct outdir_s *sp;
668 double scale, cursum;
669 int i;
670 /* initialize scaffold */
671 myScaffold.wmin = iwmax;
672 myScaffold.wmax = 0;
673 myScaffold.nic = sdt->stc[tt_Y]->ndim - 2;
674 myScaffold.rev = rev;
675 myScaffold.alen = 0;
676 myScaffold.nall = 512;
677 myScaffold.darr = (struct outdir_s *)malloc(sizeof(struct outdir_s) *
678 myScaffold.nall);
679 if (myScaffold.darr == NULL)
680 return NULL;
681 /* set up traversal */
682 cmask = (1<<myScaffold.nic) - 1;
683 for (i = myScaffold.nic; i--; )
684 pos[i+2*rev] = invec[i];
685 cmask <<= 2*rev;
686 /* grow the distribution */
687 if (SDtraverseTre(sdt->stc[tt_Y], pos, cmask,
688 build_scaffold, &myScaffold) < 0) {
689 free(myScaffold.darr);
690 return NULL;
691 }
692 /* allocate result holder */
693 cd = (SDTreCDst *)malloc(sizeof(SDTreCDst) +
694 sizeof(cd->carr[0])*myScaffold.alen);
695 if (cd == NULL) {
696 sprintf(SDerrorDetail,
697 "Cannot allocate %u entry cumulative distribution",
698 myScaffold.alen);
699 free(myScaffold.darr);
700 return NULL;
701 }
702 cd->isodist = (myScaffold.nic == 1);
703 /* sort the distribution */
704 qsort(myScaffold.darr, cd->calen = myScaffold.alen,
705 sizeof(struct outdir_s), sscmp);
706
707 /* record input range */
708 scale = myScaffold.wmin / (double)iwmax;
709 for (i = myScaffold.nic; i--; ) {
710 cd->clim[i][0] = floor(pos[i+2*rev]/scale) * scale;
711 cd->clim[i][1] = cd->clim[i][0] + scale;
712 }
713 if (cd->isodist) { /* avoid issue in SDqueryTreProjSA() */
714 cd->clim[1][0] = cd->clim[0][0];
715 cd->clim[1][1] = cd->clim[0][1];
716 }
717 cd->max_psa = myScaffold.wmax / (double)iwmax;
718 cd->max_psa *= cd->max_psa * M_PI;
719 if (rev)
720 cd->sidef = (sdt->sidef==SD_BXMIT) ? SD_FXMIT : SD_BXMIT;
721 else
722 cd->sidef = sdt->sidef;
723 cd->cTotal = 1e-20; /* compute directional total */
724 sp = myScaffold.darr;
725 for (i = myScaffold.alen; i--; sp++)
726 cd->cTotal += sp->bsdf * (double)sp->wid * sp->wid;
727 cursum = .0; /* go back and get cumulative values */
728 scale = (double)cumlmax / cd->cTotal;
729 sp = myScaffold.darr;
730 for (i = 0; i < cd->calen; i++, sp++) {
731 cd->carr[i].hndx = sp->hent;
732 cd->carr[i].cuml = scale*cursum + .5;
733 cursum += sp->bsdf * (double)sp->wid * sp->wid;
734 }
735 cd->carr[i].hndx = ~0; /* make final entry */
736 cd->carr[i].cuml = cumlmax;
737 cd->cTotal *= M_PI/(double)iwmax/iwmax;
738 /* all done, clean up and return */
739 free(myScaffold.darr);
740 return cd;
741 }
742
743 /* Find or allocate a cumulative distribution for the given incoming vector */
744 const SDCDst *
745 SDgetTreCDist(const FVECT inVec, SDComponent *sdc)
746 {
747 const SDTre *sdt;
748 double inCoord[2];
749 int i;
750 int mode;
751 SDTreCDst *cd, *cdlast;
752 /* check arguments */
753 if ((inVec == NULL) | (sdc == NULL) ||
754 (sdt = (SDTre *)sdc->dist) == NULL)
755 return NULL;
756 switch (mode = sdt->sidef) { /* check direction */
757 case SD_FREFL:
758 if (inVec[2] < 0)
759 return NULL;
760 break;
761 case SD_BREFL:
762 if (inVec[2] > 0)
763 return NULL;
764 break;
765 case SD_FXMIT:
766 if (inVec[2] < 0)
767 mode = SD_BXMIT;
768 break;
769 case SD_BXMIT:
770 if (inVec[2] > 0)
771 mode = SD_FXMIT;
772 break;
773 default:
774 return NULL;
775 }
776 if (sdt->stc[tt_Y]->ndim == 3) { /* isotropic BSDF? */
777 if (mode != sdt->sidef) /* XXX unhandled reciprocity */
778 return &SDemptyCD;
779 inCoord[0] = (.5-FTINY) -
780 .5*sqrt(inVec[0]*inVec[0] + inVec[1]*inVec[1]);
781 } else if (sdt->stc[tt_Y]->ndim == 4) {
782 if (mode != sdt->sidef) /* use reciprocity? */
783 SDdisk2square(inCoord, inVec[0], inVec[1]);
784 else
785 SDdisk2square(inCoord, -inVec[0], -inVec[1]);
786 } else
787 return NULL; /* should be internal error */
788 /* quantize to avoid f.p. errors */
789 for (i = sdt->stc[tt_Y]->ndim - 2; i--; )
790 inCoord[i] = floor(inCoord[i]/quantum)*quantum + .5*quantum;
791 cdlast = NULL; /* check for direction in cache list */
792 /* PLACE MUTEX LOCK HERE FOR THREAD-SAFE */
793 for (cd = (SDTreCDst *)sdc->cdList; cd != NULL;
794 cdlast = cd, cd = cd->next) {
795 if (cd->sidef != mode)
796 continue;
797 for (i = sdt->stc[tt_Y]->ndim - 2; i--; )
798 if ((cd->clim[i][0] > inCoord[i]) |
799 (inCoord[i] >= cd->clim[i][1]))
800 break;
801 if (i < 0)
802 break; /* means we have a match */
803 }
804 if (cd == NULL) /* need to create new entry? */
805 cdlast = cd = make_cdist(sdt, inCoord, mode != sdt->sidef);
806 if (cdlast != NULL) { /* move entry to head of cache list */
807 cdlast->next = cd->next;
808 cd->next = (SDTreCDst *)sdc->cdList;
809 sdc->cdList = (SDCDst *)cd;
810 }
811 /* END MUTEX LOCK */
812 return (SDCDst *)cd; /* ready to go */
813 }
814
815 /* Query solid angle for vector(s) */
816 static SDError
817 SDqueryTreProjSA(double *psa, const FVECT v1, const RREAL *v2,
818 int qflags, SDComponent *sdc)
819 {
820 double myPSA[2];
821 /* check arguments */
822 if ((psa == NULL) | (v1 == NULL) | (sdc == NULL) ||
823 sdc->dist == NULL)
824 return SDEargument;
825 /* get projected solid angle(s) */
826 if (v2 != NULL) {
827 const SDTre *sdt = (SDTre *)sdc->dist;
828 double hcube[SD_MAXDIM+1];
829 if (!SDqueryTre(sdt, NULL, v1, v2, hcube)) {
830 strcpy(SDerrorDetail, "Bad call to SDqueryTreProjSA");
831 return SDEinternal;
832 }
833 myPSA[0] = hcube[sdt->stc[tt_Y]->ndim];
834 myPSA[1] = myPSA[0] *= myPSA[0] * M_PI;
835 } else {
836 const SDTreCDst *cd = (const SDTreCDst *)SDgetTreCDist(v1, sdc);
837 if (cd == NULL)
838 myPSA[0] = myPSA[1] = 0;
839 else {
840 myPSA[0] = M_PI * (cd->clim[0][1] - cd->clim[0][0]) *
841 (cd->clim[1][1] - cd->clim[1][0]);
842 myPSA[1] = cd->max_psa;
843 }
844 }
845 switch (qflags) { /* record based on flag settings */
846 case SDqueryVal:
847 *psa = myPSA[0];
848 break;
849 case SDqueryMax:
850 if (myPSA[1] > *psa)
851 *psa = myPSA[1];
852 break;
853 case SDqueryMin+SDqueryMax:
854 if (myPSA[1] > psa[1])
855 psa[1] = myPSA[1];
856 /* fall through */
857 case SDqueryMin:
858 if ((myPSA[0] > 0) & (myPSA[0] < psa[0]))
859 psa[0] = myPSA[0];
860 break;
861 }
862 return SDEnone;
863 }
864
865 /* Sample cumulative distribution */
866 static SDError
867 SDsampTreCDist(FVECT ioVec, double randX, const SDCDst *cdp)
868 {
869 const unsigned nBitsC = 4*sizeof(bitmask_t);
870 const unsigned nExtraBits = 8*(sizeof(bitmask_t)-sizeof(unsigned));
871 const SDTreCDst *cd = (const SDTreCDst *)cdp;
872 const unsigned target = randX*cumlmax;
873 bitmask_t hndx, hcoord[2];
874 double gpos[3], rotangle;
875 int i, iupper, ilower;
876 /* check arguments */
877 if ((ioVec == NULL) | (cd == NULL))
878 return SDEargument;
879 if (!cd->sidef)
880 return SDEnone; /* XXX should never happen */
881 if (ioVec[2] > 0) {
882 if ((cd->sidef != SD_FREFL) & (cd->sidef != SD_FXMIT))
883 return SDEargument;
884 } else if ((cd->sidef != SD_BREFL) & (cd->sidef != SD_BXMIT))
885 return SDEargument;
886 /* binary search to find position */
887 ilower = 0; iupper = cd->calen;
888 while ((i = (iupper + ilower) >> 1) != ilower)
889 if (target >= cd->carr[i].cuml)
890 ilower = i;
891 else
892 iupper = i;
893 /* localize random position */
894 randX = (randX*cumlmax - cd->carr[ilower].cuml) /
895 (double)(cd->carr[iupper].cuml - cd->carr[ilower].cuml);
896 /* index in longer Hilbert curve */
897 hndx = (randX*cd->carr[iupper].hndx + (1.-randX)*cd->carr[ilower].hndx)
898 * (double)((bitmask_t)1 << nExtraBits);
899 /* convert Hilbert index to vector */
900 hilbert_i2c(2, nBitsC, hndx, hcoord);
901 for (i = 2; i--; )
902 gpos[i] = ((double)hcoord[i] + rand()*(1./(RAND_MAX+.5))) /
903 (double)((bitmask_t)1 << nBitsC);
904 SDsquare2disk(gpos, gpos[0], gpos[1]);
905 /* compute Z-coordinate */
906 gpos[2] = 1. - gpos[0]*gpos[0] - gpos[1]*gpos[1];
907 gpos[2] = sqrt(gpos[2]*(gpos[2]>0));
908 /* emit from back? */
909 if ((cd->sidef == SD_BREFL) | (cd->sidef == SD_FXMIT))
910 gpos[2] = -gpos[2];
911 if (cd->isodist) { /* rotate isotropic sample */
912 rotangle = atan2(-ioVec[1],-ioVec[0]);
913 spinvector(ioVec, gpos, zvec, rotangle);
914 } else
915 VCOPY(ioVec, gpos);
916 return SDEnone;
917 }
918
919 /* Advance pointer to the next non-white character in the string (or nul) */
920 static int
921 next_token(char **spp)
922 {
923 while (isspace(**spp))
924 ++*spp;
925 return **spp;
926 }
927
928 /* Advance pointer past matching token (or any token if c==0) */
929 #define eat_token(spp,c) ((next_token(spp)==(c)) ^ !(c) ? *(*(spp))++ : 0)
930
931 /* Count words from this point in string to '}' */
932 static int
933 count_values(char *cp)
934 {
935 int n = 0;
936
937 while (next_token(&cp) != '}' && *cp) {
938 while (!isspace(*cp) & (*cp != ',') & (*cp != '}'))
939 if (!*++cp)
940 break;
941 ++n;
942 eat_token(&cp, ',');
943 }
944 return n;
945 }
946
947 /* Load an array of real numbers, returning total */
948 static int
949 load_values(char **spp, float *va, int n)
950 {
951 float *v = va;
952 char *svnext;
953
954 while (n-- > 0 && (svnext = fskip(*spp)) != NULL) {
955 if ((*v++ = atof(*spp)) < 0)
956 v[-1] = 0;
957 *spp = svnext;
958 eat_token(spp, ',');
959 }
960 return v - va;
961 }
962
963 /* Load BSDF tree data */
964 static SDNode *
965 load_tree_data(char **spp, int nd)
966 {
967 SDNode *st;
968 int n;
969
970 if (!eat_token(spp, '{')) {
971 strcpy(SDerrorDetail, "Missing '{' in tensor tree");
972 return NULL;
973 }
974 if (next_token(spp) == '{') { /* tree branches */
975 st = SDnewNode(nd, -1);
976 if (st == NULL)
977 return NULL;
978 for (n = 0; n < 1<<nd; n++)
979 if ((st->u.t[n] = load_tree_data(spp, nd)) == NULL) {
980 SDfreeTre(st);
981 return NULL;
982 }
983 } else { /* else load value grid */
984 int bsiz;
985 n = count_values(*spp); /* see how big the grid is */
986 for (bsiz = 0; bsiz < 8*sizeof(size_t); bsiz += nd)
987 if (1<<bsiz == n)
988 break;
989 if (bsiz >= 8*sizeof(size_t)) {
990 strcpy(SDerrorDetail, "Illegal value count in tensor tree");
991 return NULL;
992 }
993 st = SDnewNode(nd, bsiz/nd);
994 if (st == NULL)
995 return NULL;
996 if (load_values(spp, st->u.v, n) != n) {
997 strcpy(SDerrorDetail, "Real format error in tensor tree");
998 SDfreeTre(st);
999 return NULL;
1000 }
1001 }
1002 if (!eat_token(spp, '}')) {
1003 strcpy(SDerrorDetail, "Missing '}' in tensor tree");
1004 SDfreeTre(st);
1005 return NULL;
1006 }
1007 eat_token(spp, ',');
1008 return st;
1009 }
1010
1011 /* Compute min. proj. solid angle and max. direct hemispherical scattering */
1012 static SDError
1013 get_extrema(SDSpectralDF *df)
1014 {
1015 SDNode *st = (*(SDTre *)df->comp[0].dist).stc[tt_Y];
1016 double stepWidth, dhemi, bmin[4], bmax[4];
1017
1018 stepWidth = SDsmallestLeaf(st);
1019 if (quantum > stepWidth) /* adjust quantization factor */
1020 quantum = stepWidth;
1021 df->minProjSA = M_PI*stepWidth*stepWidth;
1022 if (stepWidth < .03125)
1023 stepWidth = .03125; /* 1/32 resolution good enough */
1024 df->maxHemi = .0;
1025 if (st->ndim == 3) { /* isotropic BSDF */
1026 bmin[1] = bmin[2] = .0;
1027 bmax[1] = bmax[2] = 1.;
1028 for (bmin[0] = .0; bmin[0] < .5-FTINY; bmin[0] += stepWidth) {
1029 bmax[0] = bmin[0] + stepWidth;
1030 dhemi = SDavgTreBox(st, bmin, bmax);
1031 if (dhemi > df->maxHemi)
1032 df->maxHemi = dhemi;
1033 }
1034 } else if (st->ndim == 4) { /* anisotropic BSDF */
1035 bmin[2] = bmin[3] = .0;
1036 bmax[2] = bmax[3] = 1.;
1037 for (bmin[0] = .0; bmin[0] < 1.-FTINY; bmin[0] += stepWidth) {
1038 bmax[0] = bmin[0] + stepWidth;
1039 for (bmin[1] = .0; bmin[1] < 1.-FTINY; bmin[1] += stepWidth) {
1040 bmax[1] = bmin[1] + stepWidth;
1041 dhemi = SDavgTreBox(st, bmin, bmax);
1042 if (dhemi > df->maxHemi)
1043 df->maxHemi = dhemi;
1044 }
1045 }
1046 } else
1047 return SDEinternal;
1048 /* correct hemispherical value */
1049 df->maxHemi *= M_PI;
1050 return SDEnone;
1051 }
1052
1053 /* Load BSDF distribution for this wavelength */
1054 static SDError
1055 load_bsdf_data(SDData *sd, ezxml_t wdb, int ct, int ndim)
1056 {
1057 SDSpectralDF *df;
1058 SDTre *sdt;
1059 char *sdata;
1060 /* allocate BSDF component */
1061 sdata = ezxml_txt(ezxml_child(wdb, "WavelengthDataDirection"));
1062 if (!sdata)
1063 return SDEnone;
1064 /*
1065 * Remember that front and back are reversed from WINDOW 6 orientations
1066 */
1067 if (!strcasecmp(sdata, "Transmission Front")) {
1068 if (sd->tb == NULL && (sd->tb = SDnewSpectralDF(1)) == NULL)
1069 return SDEmemory;
1070 df = sd->tb;
1071 } else if (!strcasecmp(sdata, "Transmission Back")) {
1072 if (sd->tf == NULL && (sd->tf = SDnewSpectralDF(1)) == NULL)
1073 return SDEmemory;
1074 df = sd->tf;
1075 } else if (!strcasecmp(sdata, "Reflection Front")) {
1076 if (sd->rb == NULL && (sd->rb = SDnewSpectralDF(1)) == NULL)
1077 return SDEmemory;
1078 df = sd->rb;
1079 } else if (!strcasecmp(sdata, "Reflection Back")) {
1080 if (sd->rf == NULL && (sd->rf = SDnewSpectralDF(1)) == NULL)
1081 return SDEmemory;
1082 df = sd->rf;
1083 } else
1084 return SDEnone;
1085 /* get angle bases */
1086 sdata = ezxml_txt(ezxml_child(wdb,"AngleBasis"));
1087 if (!sdata || strcasecmp(sdata, "LBNL/Shirley-Chiu")) {
1088 sprintf(SDerrorDetail, "%s angle basis for BSDF '%s'",
1089 !sdata ? "Missing" : "Unsupported", sd->name);
1090 return !sdata ? SDEformat : SDEsupport;
1091 }
1092 if (df->comp[0].dist == NULL) { /* need to allocate BSDF tree? */
1093 sdt = (SDTre *)malloc(sizeof(SDTre));
1094 if (sdt == NULL)
1095 return SDEmemory;
1096 if (df == sd->rf)
1097 sdt->sidef = SD_FREFL;
1098 else if (df == sd->rb)
1099 sdt->sidef = SD_BREFL;
1100 else if (df == sd->tf)
1101 sdt->sidef = SD_FXMIT;
1102 else /* df == sd->tb */
1103 sdt->sidef = SD_BXMIT;
1104 sdt->stc[tt_Y] = sdt->stc[tt_u] = sdt->stc[tt_v] = NULL;
1105 df->comp[0].dist = sdt;
1106 df->comp[0].func = &SDhandleTre;
1107 } else {
1108 sdt = (SDTre *)df->comp[0].dist;
1109 if (sdt->stc[ct] != NULL) {
1110 SDfreeTre(sdt->stc[ct]);
1111 sdt->stc[ct] = NULL;
1112 }
1113 }
1114 /* read BSDF data */
1115 sdata = ezxml_txt(ezxml_child(wdb, "ScatteringData"));
1116 if (!sdata || !next_token(&sdata)) {
1117 sprintf(SDerrorDetail, "Missing BSDF ScatteringData in '%s'",
1118 sd->name);
1119 return SDEformat;
1120 }
1121 sdt->stc[ct] = load_tree_data(&sdata, ndim);
1122 if (sdt->stc[ct] == NULL)
1123 return SDEformat;
1124 if (next_token(&sdata)) { /* check for unconsumed characters */
1125 sprintf(SDerrorDetail,
1126 "Extra characters at end of ScatteringData in '%s'",
1127 sd->name);
1128 return SDEformat;
1129 }
1130 /* flatten branches where possible */
1131 sdt->stc[ct] = SDsimplifyTre(sdt->stc[ct]);
1132 if (sdt->stc[ct] == NULL)
1133 return SDEinternal;
1134 /* compute global quantities for Y */
1135 return (ct == tt_Y) ? get_extrema(df) : SDEnone;
1136 }
1137
1138 /* Find minimum value in tree */
1139 static float
1140 SDgetTreMin(const SDNode *st)
1141 {
1142 float vmin = FHUGE;
1143 int n;
1144
1145 if (st->log2GR < 0) {
1146 for (n = 1<<st->ndim; n--; ) {
1147 float v = SDgetTreMin(st->u.t[n]);
1148 if (v < vmin)
1149 vmin = v;
1150 }
1151 } else {
1152 for (n = 1<<(st->ndim*st->log2GR); n--; )
1153 if (st->u.v[n] < vmin)
1154 vmin = st->u.v[n];
1155 }
1156 return vmin;
1157 }
1158
1159 /* Subtract the given value from all tree nodes */
1160 static void
1161 SDsubtractTreVal(SDNode *st, float val)
1162 {
1163 int n;
1164
1165 if (st->log2GR < 0) {
1166 for (n = 1<<st->ndim; n--; )
1167 SDsubtractTreVal(st->u.t[n], val);
1168 } else {
1169 for (n = 1<<(st->ndim*st->log2GR); n--; )
1170 if ((st->u.v[n] -= val) < 0)
1171 st->u.v[n] = .0f;
1172 }
1173 }
1174
1175 /* Subtract minimum Y value from BSDF */
1176 static double
1177 subtract_min_Y(SDNode *st)
1178 {
1179 const float vmaxmin = 1.5/M_PI;
1180 float vmin;
1181 /* be sure to skip unused portion */
1182 if (st->ndim == 3) {
1183 int n;
1184 vmin = vmaxmin;
1185 if (st->log2GR < 0) {
1186 for (n = 0; n < 8; n += 2) {
1187 float v = SDgetTreMin(st->u.t[n]);
1188 if (v < vmin)
1189 vmin = v;
1190 }
1191 } else if (st->log2GR) {
1192 for (n = 1 << (3*st->log2GR - 1); n--; )
1193 if (st->u.v[n] < vmin)
1194 vmin = st->u.v[n];
1195 } else
1196 vmin = st->u.v[0];
1197 } else /* anisotropic covers entire tree */
1198 vmin = SDgetTreMin(st);
1199
1200 if ((vmin >= vmaxmin) | (vmin <= .01/M_PI))
1201 return .0; /* not worth bothering about */
1202
1203 SDsubtractTreVal(st, vmin);
1204
1205 return M_PI * vmin; /* return hemispherical value */
1206 }
1207
1208 /* Struct used in callback to find RGB extrema */
1209 typedef struct {
1210 SDNode **stc; /* original Y, u' & v' trees */
1211 float rgb[3]; /* RGB value */
1212 SDNode *new_stu, *new_stv; /* replacement u' & v' trees */
1213 } SDextRGBs;
1214
1215 /* Callback to find minimum RGB from Y value plus CIE (u',v') trees */
1216 static int
1217 get_min_RGB(float yval, const double *cmin, double csiz, void *cptr)
1218 {
1219 SDextRGBs *mp = (SDextRGBs *)cptr;
1220 double cmax[SD_MAXDIM];
1221 float rgb[3];
1222
1223 if (mp->stc[tt_Y]->ndim == 3) {
1224 if (cmin[0] + .5*csiz >= .5)
1225 return 0; /* ignore dead half of isotropic */
1226 } else
1227 cmax[3] = cmin[3] + csiz;
1228 cmax[0] = cmin[0] + csiz;
1229 cmax[1] = cmin[1] + csiz;
1230 cmax[2] = cmin[2] + csiz;
1231 /* average RGB color over voxel */
1232 SDyuv2rgb(yval, SDavgTreBox(mp->stc[tt_u], cmin, cmax),
1233 SDavgTreBox(mp->stc[tt_v], cmin, cmax), rgb);
1234 /* track smallest components */
1235 if (rgb[0] < mp->rgb[0]) mp->rgb[0] = rgb[0];
1236 if (rgb[1] < mp->rgb[1]) mp->rgb[1] = rgb[1];
1237 if (rgb[2] < mp->rgb[2]) mp->rgb[2] = rgb[2];
1238 return 0;
1239 }
1240
1241 /* Callback to build adjusted u' tree */
1242 static int
1243 adjust_utree(float uprime, const double *cmin, double csiz, void *cptr)
1244 {
1245 SDextRGBs *mp = (SDextRGBs *)cptr;
1246 double cmax[SD_MAXDIM];
1247 double yval;
1248 float rgb[3];
1249 C_COLOR clr;
1250
1251 if (mp->stc[tt_Y]->ndim == 3) {
1252 if (cmin[0] + .5*csiz >= .5)
1253 return 0; /* ignore dead half of isotropic */
1254 } else
1255 cmax[3] = cmin[3] + csiz;
1256 cmax[0] = cmin[0] + csiz;
1257 cmax[1] = cmin[1] + csiz;
1258 cmax[2] = cmin[2] + csiz;
1259 /* average RGB color over voxel */
1260 SDyuv2rgb(yval=SDavgTreBox(mp->stc[tt_Y], cmin, cmax), uprime,
1261 SDavgTreBox(mp->stc[tt_v], cmin, cmax), rgb);
1262 /* subtract minimum (& clamp) */
1263 if ((rgb[0] -= mp->rgb[0]) < 1e-5*yval) rgb[0] = 1e-5*yval;
1264 if ((rgb[1] -= mp->rgb[1]) < 1e-5*yval) rgb[1] = 1e-5*yval;
1265 if ((rgb[2] -= mp->rgb[2]) < 1e-5*yval) rgb[2] = 1e-5*yval;
1266 c_fromSharpRGB(rgb, &clr); /* compute new u' for adj. RGB */
1267 uprime = 4.*clr.cx/(-2.*clr.cx + 12.*clr.cy + 3.);
1268 /* assign in new u' tree */
1269 mp->new_stu = SDsetVoxel(mp->new_stu, mp->stc[tt_Y]->ndim,
1270 cmin, csiz, uprime);
1271 return -(mp->new_stu == NULL);
1272 }
1273
1274 /* Callback to build adjusted v' tree */
1275 static int
1276 adjust_vtree(float vprime, const double *cmin, double csiz, void *cptr)
1277 {
1278 SDextRGBs *mp = (SDextRGBs *)cptr;
1279 double cmax[SD_MAXDIM];
1280 double yval;
1281 float rgb[3];
1282 C_COLOR clr;
1283
1284 if (mp->stc[tt_Y]->ndim == 3) {
1285 if (cmin[0] + .5*csiz >= .5)
1286 return 0; /* ignore dead half of isotropic */
1287 } else
1288 cmax[3] = cmin[3] + csiz;
1289 cmax[0] = cmin[0] + csiz;
1290 cmax[1] = cmin[1] + csiz;
1291 cmax[2] = cmin[2] + csiz;
1292 /* average RGB color over voxel */
1293 SDyuv2rgb(yval=SDavgTreBox(mp->stc[tt_Y], cmin, cmax),
1294 SDavgTreBox(mp->stc[tt_u], cmin, cmax),
1295 vprime, rgb);
1296 /* subtract minimum (& clamp) */
1297 if ((rgb[0] -= mp->rgb[0]) < 1e-5*yval) rgb[0] = 1e-5*yval;
1298 if ((rgb[1] -= mp->rgb[1]) < 1e-5*yval) rgb[1] = 1e-5*yval;
1299 if ((rgb[2] -= mp->rgb[2]) < 1e-5*yval) rgb[2] = 1e-5*yval;
1300 c_fromSharpRGB(rgb, &clr); /* compute new v' for adj. RGB */
1301 vprime = 9.*clr.cy/(-2.*clr.cx + 12.*clr.cy + 3.);
1302 /* assign in new v' tree */
1303 mp->new_stv = SDsetVoxel(mp->new_stv, mp->stc[tt_Y]->ndim,
1304 cmin, csiz, vprime);
1305 return -(mp->new_stv == NULL);
1306 }
1307
1308 /* Subtract minimum (diffuse) color and return luminance & CIE (x,y) */
1309 static double
1310 subtract_min_RGB(C_COLOR *cs, SDNode *stc[])
1311 {
1312 SDextRGBs my_min;
1313 double ymin;
1314
1315 my_min.stc = stc;
1316 my_min.rgb[0] = my_min.rgb[1] = my_min.rgb[2] = FHUGE;
1317 my_min.new_stu = my_min.new_stv = NULL;
1318 /* get minimum RGB value */
1319 SDtraverseTre(stc[tt_Y], NULL, 0, get_min_RGB, &my_min);
1320 /* convert to C_COLOR */
1321 ymin = c_fromSharpRGB(my_min.rgb, cs);
1322 if ((ymin >= .5*FHUGE) | (ymin <= .01/M_PI))
1323 return .0; /* close to zero or no tree */
1324 /* adjust u' & v' trees */
1325 SDtraverseTre(stc[tt_u], NULL, 0, adjust_utree, &my_min);
1326 SDtraverseTre(stc[tt_v], NULL, 0, adjust_vtree, &my_min);
1327 SDfreeTre(stc[tt_u]); SDfreeTre(stc[tt_v]);
1328 stc[tt_u] = SDsimplifyTre(my_min.new_stu);
1329 stc[tt_v] = SDsimplifyTre(my_min.new_stv);
1330 /* subtract Y & return hemispherical */
1331 SDsubtractTreVal(stc[tt_Y], ymin);
1332
1333 return M_PI * ymin;
1334 }
1335
1336 /* Extract and separate diffuse portion of BSDF */
1337 static void
1338 extract_diffuse(SDValue *dv, SDSpectralDF *df)
1339 {
1340 int n;
1341 SDTre *sdt;
1342
1343 if (df == NULL || df->ncomp <= 0) {
1344 dv->spec = c_dfcolor;
1345 dv->cieY = .0;
1346 return;
1347 }
1348 sdt = (SDTre *)df->comp[0].dist;
1349 /* subtract minimum color/grayscale */
1350 if (sdt->stc[tt_u] != NULL && sdt->stc[tt_v] != NULL) {
1351 int i = 3*(tt_RGB_coef[1] < .001);
1352 while (i--) { /* initialize on first call */
1353 float rgb[3];
1354 rgb[0] = rgb[1] = rgb[2] = .0f; rgb[i] = 1.f;
1355 tt_RGB_coef[i] = c_fromSharpRGB(rgb, &tt_RGB_prim[i]);
1356 }
1357 memcpy(df->comp[0].cspec, tt_RGB_prim, sizeof(tt_RGB_prim));
1358 dv->cieY = subtract_min_RGB(&dv->spec, sdt->stc);
1359 } else {
1360 df->comp[0].cspec[0] = dv->spec = c_dfcolor;
1361 dv->cieY = subtract_min_Y(sdt->stc[tt_Y]);
1362 }
1363 df->maxHemi -= dv->cieY; /* adjust maximum hemispherical */
1364 /* make sure everything is set */
1365 c_ccvt(&dv->spec, C_CSXY+C_CSSPEC);
1366 }
1367
1368 /* Load a variable-resolution BSDF tree from an open XML file */
1369 SDError
1370 SDloadTre(SDData *sd, ezxml_t wtl)
1371 {
1372 SDError ec;
1373 ezxml_t wld, wdb;
1374 int rank;
1375 char *txt;
1376 /* basic checks and tensor rank */
1377 txt = ezxml_txt(ezxml_child(ezxml_child(wtl,
1378 "DataDefinition"), "IncidentDataStructure"));
1379 if (txt == NULL || !*txt) {
1380 sprintf(SDerrorDetail,
1381 "BSDF \"%s\": missing IncidentDataStructure",
1382 sd->name);
1383 return SDEformat;
1384 }
1385 if (!strcasecmp(txt, "TensorTree3"))
1386 rank = 3;
1387 else if (!strcasecmp(txt, "TensorTree4"))
1388 rank = 4;
1389 else {
1390 sprintf(SDerrorDetail,
1391 "BSDF \"%s\": unsupported IncidentDataStructure",
1392 sd->name);
1393 return SDEsupport;
1394 }
1395 /* load BSDF components */
1396 for (wld = ezxml_child(wtl, "WavelengthData");
1397 wld != NULL; wld = wld->next) {
1398 const char *cnm = ezxml_txt(ezxml_child(wld,"Wavelength"));
1399 int ct = -1;
1400 if (!strcasecmp(cnm, "Visible"))
1401 ct = tt_Y;
1402 else if (!strcasecmp(cnm, "CIE-u"))
1403 ct = tt_u;
1404 else if (!strcasecmp(cnm, "CIE-v"))
1405 ct = tt_v;
1406 else
1407 continue;
1408 for (wdb = ezxml_child(wld, "WavelengthDataBlock");
1409 wdb != NULL; wdb = wdb->next)
1410 if ((ec = load_bsdf_data(sd, wdb, ct, rank)) != SDEnone)
1411 return ec;
1412 }
1413 /* separate diffuse components */
1414 extract_diffuse(&sd->rLambFront, sd->rf);
1415 extract_diffuse(&sd->rLambBack, sd->rb);
1416 if (sd->tf != NULL)
1417 extract_diffuse(&sd->tLamb, sd->tf);
1418 if (sd->tb != NULL)
1419 extract_diffuse(&sd->tLamb, sd->tb);
1420 /* return success */
1421 return SDEnone;
1422 }
1423
1424 /* Variable resolution BSDF methods */
1425 const SDFunc SDhandleTre = {
1426 &SDgetTreBSDF,
1427 &SDqueryTreProjSA,
1428 &SDgetTreCDist,
1429 &SDsampTreCDist,
1430 &SDFreeBTre,
1431 };