ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/m_bsdf.c
Revision: 2.64
Committed: Wed Aug 25 16:12:21 2021 UTC (2 years, 8 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.63: +11 -4 lines
Log Message:
fix: Improved duplicate sample testing, which improves peak extraction

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.64 static const char RCSid[] = "$Id: m_bsdf.c,v 2.63 2021/03/27 20:08:58 greg Exp $";
3 greg 2.1 #endif
4     /*
5     * Shading for materials with BSDFs taken from XML data files
6     */
7    
8     #include "copyright.h"
9    
10     #include "ray.h"
11 greg 2.50 #include "otypes.h"
12 greg 2.1 #include "ambient.h"
13     #include "source.h"
14     #include "func.h"
15     #include "bsdf.h"
16     #include "random.h"
17 greg 2.30 #include "pmapmat.h"
18 greg 2.1
19     /*
20 greg 2.50 * Arguments to this material include optional diffuse colors.
21 greg 2.1 * String arguments include the BSDF and function files.
22 greg 2.50 * For the MAT_BSDF type, a non-zero thickness causes the useful behavior
23 greg 2.5 * of translating transmitted rays this distance beneath the surface
24     * (opposite the surface normal) to bypass any intervening geometry.
25     * Translation only affects scattered, non-source-directed samples.
26     * A non-zero thickness has the further side-effect that an unscattered
27 greg 2.35 * (view) ray will pass right through our material, making the BSDF
28     * surface invisible and showing the proxied geometry instead. Thickness
29     * has the further effect of turning off reflection on the reverse side so
30     * rays heading in the opposite direction pass unimpeded through the BSDF
31 greg 2.5 * surface. A paired surface may be placed on the opposide side of
32     * the detail geometry, less than this thickness away, if a two-way
33     * proxy is desired. Note that the sign of the thickness is important.
34     * A positive thickness hides geometry behind the BSDF surface and uses
35     * front reflectance and transmission properties. A negative thickness
36     * hides geometry in front of the surface when rays hit from behind,
37     * and applies only the transmission and backside reflectance properties.
38     * Reflection is ignored on the hidden side, as those rays pass through.
39 greg 2.52 * For the MAT_ABSDF type, we check for a strong "through" component.
40 greg 2.50 * Such a component will cause direct rays to pass through unscattered.
41 greg 2.40 * A separate test prevents over-counting by dropping samples that are
42     * too close to this "through" direction. BSDFs with such a through direction
43     * will also have a view component, meaning they are somewhat see-through.
44 greg 2.52 * A MAT_BSDF type with zero thickness behaves the same as a MAT_ABSDF
45 greg 2.50 * type with no strong through component.
46 greg 2.1 * The "up" vector for the BSDF is given by three variables, defined
47     * (along with the thickness) by the named function file, or '.' if none.
48     * Together with the surface normal, this defines the local coordinate
49     * system for the BSDF.
50     * We do not reorient the surface, so if the BSDF has no back-side
51 greg 2.5 * reflectance and none is given in the real arguments, a BSDF surface
52     * with zero thickness will appear black when viewed from behind
53 greg 2.35 * unless backface visibility is on, when it becomes invisible.
54 greg 2.5 * The diffuse arguments are added to components in the BSDF file,
55 greg 2.1 * not multiplied. However, patterns affect this material as a multiplier
56     * on everything except non-diffuse reflection.
57     *
58 greg 2.52 * Arguments for MAT_ABSDF are:
59 greg 2.50 * 5+ BSDFfile ux uy uz funcfile transform
60     * 0
61     * 0|3|6|9 rdf gdf bdf
62     * rdb gdb bdb
63     * rdt gdt bdt
64     *
65 greg 2.1 * Arguments for MAT_BSDF are:
66     * 6+ thick BSDFfile ux uy uz funcfile transform
67     * 0
68 greg 2.8 * 0|3|6|9 rdf gdf bdf
69 greg 2.1 * rdb gdb bdb
70     * rdt gdt bdt
71     */
72    
73 greg 2.4 /*
74     * Note that our reverse ray-tracing process means that the positions
75     * of incoming and outgoing vectors may be reversed in our calls
76 greg 2.35 * to the BSDF library. This is usually fine, since the bidirectional nature
77 greg 2.4 * of the BSDF (that's what the 'B' stands for) means it all works out.
78     */
79    
80 greg 2.1 typedef struct {
81     OBJREC *mp; /* material pointer */
82     RAY *pr; /* intersected ray */
83     FVECT pnorm; /* perturbed surface normal */
84 greg 2.4 FVECT vray; /* local outgoing (return) vector */
85 greg 2.9 double sr_vpsa[2]; /* sqrt of BSDF projected solid angle extrema */
86 greg 2.1 RREAL toloc[3][3]; /* world to local BSDF coords */
87     RREAL fromloc[3][3]; /* local BSDF coords to world */
88     double thick; /* surface thickness */
89 greg 2.52 COLOR cthru; /* "through" component for MAT_ABSDF */
90 greg 2.61 COLOR cthru_surr; /* surround for "through" component */
91 greg 2.1 SDData *sd; /* loaded BSDF data */
92 greg 2.31 COLOR rdiff; /* diffuse reflection */
93 greg 2.39 COLOR runsamp; /* BSDF hemispherical reflection */
94 greg 2.31 COLOR tdiff; /* diffuse transmission */
95 greg 2.39 COLOR tunsamp; /* BSDF hemispherical transmission */
96 greg 2.1 } BSDFDAT; /* BSDF material data */
97    
98     #define cvt_sdcolor(cv, svp) ccy2rgb(&(svp)->spec, (svp)->cieY, cv)
99    
100 greg 2.58 typedef struct {
101     double vy; /* brightness (for sorting) */
102     FVECT tdir; /* through sample direction (normalized) */
103     COLOR vcol; /* BTDF color */
104     } PEAKSAMP; /* BTDF peak sample */
105    
106     /* Comparison function to put near-peak values in descending order */
107     static int
108     cmp_psamp(const void *p1, const void *p2)
109     {
110     double diff = (*(const PEAKSAMP *)p1).vy - (*(const PEAKSAMP *)p2).vy;
111     if (diff > 0) return(-1);
112     if (diff < 0) return(1);
113     return(0);
114     }
115    
116 greg 2.52 /* Compute "through" component color for MAT_ABSDF */
117 greg 2.34 static void
118     compute_through(BSDFDAT *ndp)
119     {
120 greg 2.60 #define NDIR2CHECK 29
121 greg 2.34 static const float dir2check[NDIR2CHECK][2] = {
122 greg 2.60 {0, 0}, {-0.6, 0}, {0, 0.6},
123     {0, -0.6}, {0.6, 0}, {-0.6, 0.6},
124     {-0.6, -0.6}, {0.6, 0.6}, {0.6, -0.6},
125     {-1.2, 0}, {0, 1.2}, {0, -1.2},
126     {1.2, 0}, {-1.2, 1.2}, {-1.2, -1.2},
127     {1.2, 1.2}, {1.2, -1.2}, {-1.8, 0},
128     {0, 1.8}, {0, -1.8}, {1.8, 0},
129     {-1.8, 1.8}, {-1.8, -1.8}, {1.8, 1.8},
130     {1.8, -1.8}, {-2.4, 0}, {0, 2.4},
131     {0, -2.4}, {2.4, 0},
132 greg 2.34 };
133 greg 2.64 #define neighbors(i,j) \
134     ((dir2check[i][0]-dir2check[j][0])*(dir2check[i][0]-dir2check[j][0]) + \
135     (dir2check[i][1]-dir2check[j][1])*(dir2check[i][1]-dir2check[j][1]) <= 0.73)
136 greg 2.58 const double peak_over = 1.5;
137     PEAKSAMP psamp[NDIR2CHECK];
138 greg 2.34 SDSpectralDF *dfp;
139     FVECT pdir;
140     double tomega, srchrad;
141 greg 2.61 double tomsum, tomsurr;
142     COLOR vpeak, vsurr;
143     double vypeak;
144 greg 2.64 int i, j, ns;
145 greg 2.34 SDError ec;
146    
147     if (ndp->pr->rod > 0)
148     dfp = (ndp->sd->tf != NULL) ? ndp->sd->tf : ndp->sd->tb;
149     else
150     dfp = (ndp->sd->tb != NULL) ? ndp->sd->tb : ndp->sd->tf;
151    
152     if (dfp == NULL)
153     return; /* no specular transmission */
154     if (bright(ndp->pr->pcol) <= FTINY)
155     return; /* pattern is black, here */
156 greg 2.58 srchrad = sqrt(dfp->minProjSA); /* else evaluate peak */
157     for (i = 0; i < NDIR2CHECK; i++) {
158     SDValue sv;
159     psamp[i].tdir[0] = -ndp->vray[0] + dir2check[i][0]*srchrad;
160     psamp[i].tdir[1] = -ndp->vray[1] + dir2check[i][1]*srchrad;
161     psamp[i].tdir[2] = -ndp->vray[2];
162     normalize(psamp[i].tdir);
163     ec = SDevalBSDF(&sv, psamp[i].tdir, ndp->vray, ndp->sd);
164     if (ec)
165     goto baderror;
166     cvt_sdcolor(psamp[i].vcol, &sv);
167 greg 2.61 psamp[i].vy = sv.cieY;
168 greg 2.58 }
169     qsort(psamp, NDIR2CHECK, sizeof(PEAKSAMP), cmp_psamp);
170 greg 2.61 if (psamp[0].vy <= FTINY)
171     return; /* zero area */
172 greg 2.40 setcolor(vpeak, 0, 0, 0);
173 greg 2.61 setcolor(vsurr, 0, 0, 0);
174     vypeak = tomsum = tomsurr = 0; /* combine top unique values */
175     ns = 0;
176 greg 2.34 for (i = 0; i < NDIR2CHECK; i++) {
177 greg 2.64 for (j = i; j--; ) /* check for duplicate sample */
178     if (psamp[j].vy == psamp[i].vy && neighbors(i,j))
179     break;
180     if (j >= 0)
181     continue; /* skip duplicate */
182 greg 2.61
183 greg 2.58 ec = SDsizeBSDF(&tomega, psamp[i].tdir, ndp->vray,
184     SDqueryMin, ndp->sd);
185 greg 2.34 if (ec)
186     goto baderror;
187 greg 2.61 /* not really a peak? */
188     if (tomega > 1.5*dfp->minProjSA ||
189     vypeak > 8.*psamp[i].vy*ns) {
190     if (!i) return; /* abort */
191     scalecolor(psamp[i].vcol, tomega);
192     addcolor(vsurr, psamp[i].vcol);
193     tomsurr += tomega;
194 greg 2.58 continue;
195 greg 2.34 }
196 greg 2.58 scalecolor(psamp[i].vcol, tomega);
197     addcolor(vpeak, psamp[i].vcol);
198     tomsum += tomega;
199     vypeak += psamp[i].vy;
200     ++ns;
201     }
202 greg 2.61 if (vypeak*tomsurr < peak_over*bright(vsurr)*ns)
203 greg 2.58 return; /* peak not peaky enough */
204 greg 2.62 if ((vypeak/ns - (ndp->vray[2] > 0 ? ndp->sd->tLambFront.cieY
205     : ndp->sd->tLambBack.cieY)*(1./PI))*tomsum <= .001)
206 greg 2.40 return; /* < 0.1% transmission */
207 greg 2.58 copycolor(ndp->cthru, vpeak); /* already scaled by omega */
208 greg 2.34 multcolor(ndp->cthru, ndp->pr->pcol); /* modify by pattern */
209 greg 2.61 if (tomsurr > FTINY) { /* surround contribution? */
210     scalecolor(vsurr, 1./tomsurr); /* this one is avg. BTDF */
211     copycolor(ndp->cthru_surr, vsurr);
212     multcolor(ndp->cthru_surr, ndp->pr->pcol);
213     }
214 greg 2.34 return;
215     baderror:
216     objerror(ndp->mp, USER, transSDError(ec));
217 greg 2.64 #undef neighbors
218 greg 2.34 #undef NDIR2CHECK
219     }
220    
221 greg 2.4 /* Jitter ray sample according to projected solid angle and specjitter */
222     static void
223 greg 2.15 bsdf_jitter(FVECT vres, BSDFDAT *ndp, double sr_psa)
224 greg 2.4 {
225     VCOPY(vres, ndp->vray);
226     if (specjitter < 1.)
227     sr_psa *= specjitter;
228     if (sr_psa <= FTINY)
229     return;
230     vres[0] += sr_psa*(.5 - frandom());
231     vres[1] += sr_psa*(.5 - frandom());
232     normalize(vres);
233     }
234    
235 greg 2.33 /* Get BSDF specular for direct component, returning true if OK to proceed */
236 greg 2.7 static int
237 greg 2.33 direct_specular_OK(COLOR cval, FVECT ldir, double omega, BSDFDAT *ndp)
238 greg 2.7 {
239 greg 2.43 int nsamp;
240     double wtot = 0;
241 greg 2.13 FVECT vsrc, vsmp, vjit;
242 greg 2.36 double tomega, tomega2;
243 greg 2.15 double sf, tsr, sd[2];
244 greg 2.32 COLOR csmp, cdiff;
245     double diffY;
246 greg 2.7 SDValue sv;
247     SDError ec;
248 greg 2.13 int i;
249 greg 2.37 /* in case we fail */
250 greg 2.40 setcolor(cval, 0, 0, 0);
251 greg 2.7 /* transform source direction */
252     if (SDmapDir(vsrc, ndp->toloc, ldir) != SDEnone)
253     return(0);
254 greg 2.32 /* will discount diffuse portion */
255     switch ((vsrc[2] > 0)<<1 | (ndp->vray[2] > 0)) {
256     case 3:
257     if (ndp->sd->rf == NULL)
258     return(0); /* all diffuse */
259     sv = ndp->sd->rLambFront;
260     break;
261     case 0:
262     if (ndp->sd->rb == NULL)
263     return(0); /* all diffuse */
264     sv = ndp->sd->rLambBack;
265     break;
266 greg 2.62 case 1:
267 greg 2.32 if ((ndp->sd->tf == NULL) & (ndp->sd->tb == NULL))
268     return(0); /* all diffuse */
269 greg 2.62 sv = ndp->sd->tLambFront;
270     break;
271     case 2:
272     if ((ndp->sd->tf == NULL) & (ndp->sd->tb == NULL))
273     return(0); /* all diffuse */
274     sv = ndp->sd->tLambBack;
275 greg 2.32 break;
276     }
277 greg 2.33 if (sv.cieY > FTINY) {
278     diffY = sv.cieY *= 1./PI;
279 greg 2.32 cvt_sdcolor(cdiff, &sv);
280     } else {
281 greg 2.40 diffY = 0;
282     setcolor(cdiff, 0, 0, 0);
283 greg 2.32 }
284 greg 2.53 /* need projected solid angle */
285 greg 2.37 omega *= fabs(vsrc[2]);
286 greg 2.13 /* check indirect over-counting */
287 greg 2.40 if ((vsrc[2] > 0) ^ (ndp->vray[2] > 0) && bright(ndp->cthru) > FTINY) {
288 greg 2.53 double dx = vsrc[0] + ndp->vray[0];
289     double dy = vsrc[1] + ndp->vray[1];
290     SDSpectralDF *dfp = (ndp->pr->rod > 0) ?
291     ((ndp->sd->tf != NULL) ? ndp->sd->tf : ndp->sd->tb) :
292     ((ndp->sd->tb != NULL) ? ndp->sd->tb : ndp->sd->tf) ;
293    
294     if (dx*dx + dy*dy <= (2.5*4./PI)*(omega + dfp->minProjSA +
295 greg 2.61 2.*sqrt(omega*dfp->minProjSA))) {
296     if (bright(ndp->cthru_surr) <= FTINY)
297     return(0);
298     copycolor(cval, ndp->cthru_surr);
299     return(1); /* return non-zero surround BTDF */
300     }
301 greg 2.7 }
302 greg 2.53 ec = SDsizeBSDF(&tomega, ndp->vray, vsrc, SDqueryMin, ndp->sd);
303     if (ec)
304     goto baderror;
305 greg 2.37 /* assign number of samples */
306 greg 2.15 sf = specjitter * ndp->pr->rweight;
307 greg 2.40 if (tomega <= 0)
308 greg 2.24 nsamp = 1;
309     else if (25.*tomega <= omega)
310 greg 2.15 nsamp = 100.*sf + .5;
311     else
312     nsamp = 4.*sf*omega/tomega + .5;
313     nsamp += !nsamp;
314 greg 2.37 sf = sqrt(omega); /* sample our source area */
315 greg 2.15 tsr = sqrt(tomega);
316 greg 2.13 for (i = nsamp; i--; ) {
317     VCOPY(vsmp, vsrc); /* jitter query directions */
318     if (nsamp > 1) {
319     multisamp(sd, 2, (i + frandom())/(double)nsamp);
320     vsmp[0] += (sd[0] - .5)*sf;
321     vsmp[1] += (sd[1] - .5)*sf;
322 greg 2.36 normalize(vsmp);
323 greg 2.13 }
324 greg 2.15 bsdf_jitter(vjit, ndp, tsr);
325 greg 2.37 /* compute BSDF */
326     ec = SDevalBSDF(&sv, vjit, vsmp, ndp->sd);
327     if (ec)
328     goto baderror;
329     if (sv.cieY - diffY <= FTINY)
330     continue; /* no specular part */
331 greg 2.36 /* check for variable resolution */
332     ec = SDsizeBSDF(&tomega2, vjit, vsmp, SDqueryMin, ndp->sd);
333     if (ec)
334     goto baderror;
335     if (tomega2 < .12*tomega)
336     continue; /* not safe to include */
337 greg 2.13 cvt_sdcolor(csmp, &sv);
338 greg 2.59 #if 0
339     if (sf < 2.5*tsr) { /* weight by BSDF for small sources */
340 greg 2.44 scalecolor(csmp, sv.cieY);
341     wtot += sv.cieY;
342     } else
343 greg 2.59 #endif
344     wtot += 1.;
345 greg 2.43 addcolor(cval, csmp);
346 greg 2.13 }
347 greg 2.43 if (wtot <= FTINY) /* no valid specular samples? */
348 greg 2.37 return(0);
349    
350 greg 2.43 sf = 1./wtot; /* weighted average BSDF */
351 greg 2.13 scalecolor(cval, sf);
352 greg 2.32 /* subtract diffuse contribution */
353     for (i = 3*(diffY > FTINY); i--; )
354 greg 2.40 if ((colval(cval,i) -= colval(cdiff,i)) < 0)
355     colval(cval,i) = 0;
356 greg 2.32 return(1);
357 greg 2.13 baderror:
358     objerror(ndp->mp, USER, transSDError(ec));
359 greg 2.17 return(0); /* gratis return */
360 greg 2.7 }
361    
362 greg 2.5 /* Compute source contribution for BSDF (reflected & transmitted) */
363 greg 2.1 static void
364 greg 2.5 dir_bsdf(
365 greg 2.1 COLOR cval, /* returned coefficient */
366     void *nnp, /* material data */
367     FVECT ldir, /* light source direction */
368     double omega /* light source size */
369     )
370     {
371 greg 2.3 BSDFDAT *np = (BSDFDAT *)nnp;
372 greg 2.1 double ldot;
373     double dtmp;
374     COLOR ctmp;
375    
376 greg 2.40 setcolor(cval, 0, 0, 0);
377 greg 2.1
378     ldot = DOT(np->pnorm, ldir);
379     if ((-FTINY <= ldot) & (ldot <= FTINY))
380     return;
381    
382 greg 2.9 if (ldot > 0 && bright(np->rdiff) > FTINY) {
383 greg 2.1 /*
384 greg 2.39 * Compute diffuse reflected component
385 greg 2.1 */
386     copycolor(ctmp, np->rdiff);
387     dtmp = ldot * omega * (1./PI);
388     scalecolor(ctmp, dtmp);
389     addcolor(cval, ctmp);
390     }
391 greg 2.9 if (ldot < 0 && bright(np->tdiff) > FTINY) {
392 greg 2.1 /*
393 greg 2.39 * Compute diffuse transmission
394 greg 2.1 */
395     copycolor(ctmp, np->tdiff);
396     dtmp = -ldot * omega * (1.0/PI);
397     scalecolor(ctmp, dtmp);
398     addcolor(cval, ctmp);
399     }
400 greg 2.30 if (ambRayInPmap(np->pr))
401     return; /* specular already in photon map */
402 greg 2.1 /*
403 greg 2.39 * Compute specular scattering coefficient using BSDF
404 greg 2.1 */
405 greg 2.33 if (!direct_specular_OK(ctmp, ldir, omega, np))
406 greg 2.1 return;
407 greg 2.31 if (ldot < 0) { /* pattern for specular transmission */
408 greg 2.1 multcolor(ctmp, np->pr->pcol);
409     dtmp = -ldot * omega;
410 greg 2.31 } else
411     dtmp = ldot * omega;
412 greg 2.1 scalecolor(ctmp, dtmp);
413     addcolor(cval, ctmp);
414     }
415    
416 greg 2.5 /* Compute source contribution for BSDF (reflected only) */
417     static void
418     dir_brdf(
419     COLOR cval, /* returned coefficient */
420     void *nnp, /* material data */
421     FVECT ldir, /* light source direction */
422     double omega /* light source size */
423     )
424     {
425     BSDFDAT *np = (BSDFDAT *)nnp;
426     double ldot;
427     double dtmp;
428     COLOR ctmp, ctmp1, ctmp2;
429    
430 greg 2.40 setcolor(cval, 0, 0, 0);
431 greg 2.5
432     ldot = DOT(np->pnorm, ldir);
433    
434     if (ldot <= FTINY)
435     return;
436    
437     if (bright(np->rdiff) > FTINY) {
438     /*
439 greg 2.39 * Compute diffuse reflected component
440 greg 2.5 */
441     copycolor(ctmp, np->rdiff);
442     dtmp = ldot * omega * (1./PI);
443     scalecolor(ctmp, dtmp);
444     addcolor(cval, ctmp);
445     }
446 greg 2.30 if (ambRayInPmap(np->pr))
447     return; /* specular already in photon map */
448 greg 2.5 /*
449 greg 2.39 * Compute specular reflection coefficient using BSDF
450 greg 2.5 */
451 greg 2.33 if (!direct_specular_OK(ctmp, ldir, omega, np))
452 greg 2.5 return;
453     dtmp = ldot * omega;
454     scalecolor(ctmp, dtmp);
455     addcolor(cval, ctmp);
456     }
457    
458     /* Compute source contribution for BSDF (transmitted only) */
459     static void
460     dir_btdf(
461     COLOR cval, /* returned coefficient */
462     void *nnp, /* material data */
463     FVECT ldir, /* light source direction */
464     double omega /* light source size */
465     )
466     {
467     BSDFDAT *np = (BSDFDAT *)nnp;
468     double ldot;
469     double dtmp;
470     COLOR ctmp;
471    
472 greg 2.40 setcolor(cval, 0, 0, 0);
473 greg 2.5
474     ldot = DOT(np->pnorm, ldir);
475    
476     if (ldot >= -FTINY)
477     return;
478    
479     if (bright(np->tdiff) > FTINY) {
480     /*
481 greg 2.39 * Compute diffuse transmission
482 greg 2.5 */
483     copycolor(ctmp, np->tdiff);
484     dtmp = -ldot * omega * (1.0/PI);
485     scalecolor(ctmp, dtmp);
486     addcolor(cval, ctmp);
487     }
488 greg 2.30 if (ambRayInPmap(np->pr))
489     return; /* specular already in photon map */
490 greg 2.5 /*
491 greg 2.39 * Compute specular scattering coefficient using BSDF
492 greg 2.5 */
493 greg 2.33 if (!direct_specular_OK(ctmp, ldir, omega, np))
494 greg 2.5 return;
495     /* full pattern on transmission */
496     multcolor(ctmp, np->pr->pcol);
497     dtmp = -ldot * omega;
498     scalecolor(ctmp, dtmp);
499     addcolor(cval, ctmp);
500     }
501    
502 greg 2.1 /* Sample separate BSDF component */
503     static int
504 greg 2.40 sample_sdcomp(BSDFDAT *ndp, SDComponent *dcp, int xmit)
505 greg 2.1 {
506 greg 2.47 const int hasthru = (xmit &&
507     !(ndp->pr->crtype & (SPECULAR|AMBIENT))
508     && bright(ndp->cthru) > FTINY);
509 greg 2.41 int nstarget = 1;
510     int nsent = 0;
511     int n;
512     SDError ec;
513     SDValue bsv;
514     double xrand;
515     FVECT vsmp, vinc;
516     RAY sr;
517 greg 2.1 /* multiple samples? */
518     if (specjitter > 1.5) {
519     nstarget = specjitter*ndp->pr->rweight + .5;
520 greg 2.14 nstarget += !nstarget;
521 greg 2.1 }
522 greg 2.11 /* run through our samples */
523 greg 2.40 for (n = 0; n < nstarget; n++) {
524 greg 2.15 if (nstarget == 1) { /* stratify random variable */
525 greg 2.11 xrand = urand(ilhash(dimlist,ndims)+samplendx);
526 greg 2.15 if (specjitter < 1.)
527     xrand = .5 + specjitter*(xrand-.5);
528     } else {
529 greg 2.40 xrand = (n + frandom())/(double)nstarget;
530 greg 2.15 }
531 greg 2.11 SDerrorDetail[0] = '\0'; /* sample direction & coef. */
532 greg 2.15 bsdf_jitter(vsmp, ndp, ndp->sr_vpsa[0]);
533 greg 2.40 VCOPY(vinc, vsmp); /* to compare after */
534 greg 2.11 ec = SDsampComponent(&bsv, vsmp, xrand, dcp);
535 greg 2.1 if (ec)
536 greg 2.2 objerror(ndp->mp, USER, transSDError(ec));
537 greg 2.11 if (bsv.cieY <= FTINY) /* zero component? */
538 greg 2.1 break;
539 greg 2.40 if (hasthru) { /* check for view ray */
540     double dx = vinc[0] + vsmp[0];
541     double dy = vinc[1] + vsmp[1];
542     if (dx*dx + dy*dy <= ndp->sr_vpsa[0]*ndp->sr_vpsa[0])
543     continue; /* exclude view sample */
544     }
545     /* map non-view sample->world */
546 greg 2.4 if (SDmapDir(sr.rdir, ndp->fromloc, vsmp) != SDEnone)
547 greg 2.1 break;
548     /* spawn a specular ray */
549     if (nstarget > 1)
550     bsv.cieY /= (double)nstarget;
551 greg 2.11 cvt_sdcolor(sr.rcoef, &bsv); /* use sample color */
552 greg 2.40 if (xmit) /* apply pattern on transmit */
553 greg 2.1 multcolor(sr.rcoef, ndp->pr->pcol);
554     if (rayorigin(&sr, SPECULAR, ndp->pr, sr.rcoef) < 0) {
555 greg 2.48 if (!n & (nstarget > 1)) {
556 greg 2.49 n = nstarget; /* avoid infinitue loop */
557 greg 2.48 nstarget = nstarget*sr.rweight/minweight;
558 greg 2.49 if (n == nstarget) break;
559 greg 2.48 n = -1; /* moved target */
560     }
561     continue; /* try again */
562 greg 2.1 }
563 greg 2.40 if (xmit && ndp->thick != 0) /* need to offset origin? */
564 greg 2.5 VSUM(sr.rorg, sr.rorg, ndp->pr->ron, -ndp->thick);
565 greg 2.1 rayvalue(&sr); /* send & evaluate sample */
566     multcolor(sr.rcol, sr.rcoef);
567     addcolor(ndp->pr->rcol, sr.rcol);
568 greg 2.40 ++nsent;
569 greg 2.1 }
570     return(nsent);
571     }
572    
573     /* Sample non-diffuse components of BSDF */
574     static int
575     sample_sdf(BSDFDAT *ndp, int sflags)
576     {
577 greg 2.46 int hasthru = (sflags == SDsampSpT &&
578 greg 2.47 !(ndp->pr->crtype & (SPECULAR|AMBIENT))
579     && bright(ndp->cthru) > FTINY);
580 greg 2.1 int n, ntotal = 0;
581 greg 2.40 double b = 0;
582 greg 2.1 SDSpectralDF *dfp;
583     COLORV *unsc;
584    
585     if (sflags == SDsampSpT) {
586 greg 2.39 unsc = ndp->tunsamp;
587 greg 2.22 if (ndp->pr->rod > 0)
588     dfp = (ndp->sd->tf != NULL) ? ndp->sd->tf : ndp->sd->tb;
589     else
590     dfp = (ndp->sd->tb != NULL) ? ndp->sd->tb : ndp->sd->tf;
591 greg 2.1 } else /* sflags == SDsampSpR */ {
592 greg 2.39 unsc = ndp->runsamp;
593 greg 2.31 if (ndp->pr->rod > 0)
594 greg 2.1 dfp = ndp->sd->rf;
595 greg 2.31 else
596 greg 2.1 dfp = ndp->sd->rb;
597     }
598 greg 2.40 setcolor(unsc, 0, 0, 0);
599 greg 2.1 if (dfp == NULL) /* no specular component? */
600     return(0);
601 greg 2.40
602     if (hasthru) { /* separate view sample? */
603     RAY tr;
604     if (rayorigin(&tr, TRANS, ndp->pr, ndp->cthru) == 0) {
605     VCOPY(tr.rdir, ndp->pr->rdir);
606     rayvalue(&tr);
607     multcolor(tr.rcol, tr.rcoef);
608     addcolor(ndp->pr->rcol, tr.rcol);
609 greg 2.56 ndp->pr->rxt = ndp->pr->rot + raydistance(&tr);
610 greg 2.40 ++ntotal;
611     b = bright(ndp->cthru);
612     } else
613     hasthru = 0;
614     }
615 greg 2.43 if (dfp->maxHemi - b <= FTINY) { /* have specular to sample? */
616 greg 2.40 b = 0;
617     } else {
618     FVECT vjit;
619     bsdf_jitter(vjit, ndp, ndp->sr_vpsa[1]);
620     b = SDdirectHemi(vjit, sflags, ndp->sd) - b;
621     if (b < 0) b = 0;
622     }
623     if (b <= specthresh+FTINY) { /* below sampling threshold? */
624     if (b > FTINY) { /* XXX no color from BSDF */
625 greg 2.1 if (sflags == SDsampSpT) {
626 greg 2.39 copycolor(unsc, ndp->pr->pcol);
627 greg 2.40 scalecolor(unsc, b);
628 greg 2.1 } else /* no pattern on reflection */
629 greg 2.40 setcolor(unsc, b, b, b);
630 greg 2.1 }
631 greg 2.40 return(ntotal);
632 greg 2.1 }
633 greg 2.41 dimlist[ndims] = (int)(size_t)ndp->mp; /* else sample specular */
634     ndims += 2;
635 greg 2.1 for (n = dfp->ncomp; n--; ) { /* loop over components */
636     dimlist[ndims-1] = n + 9438;
637     ntotal += sample_sdcomp(ndp, &dfp->comp[n], sflags==SDsampSpT);
638     }
639     ndims -= 2;
640     return(ntotal);
641     }
642    
643     /* Color a ray that hit a BSDF material */
644     int
645     m_bsdf(OBJREC *m, RAY *r)
646     {
647 greg 2.50 int hasthick = (m->otype == MAT_BSDF);
648 greg 2.6 int hitfront;
649 greg 2.1 COLOR ctmp;
650     SDError ec;
651 greg 2.5 FVECT upvec, vtmp;
652 greg 2.1 MFUNC *mf;
653     BSDFDAT nd;
654     /* check arguments */
655 greg 2.50 if ((m->oargs.nsargs < hasthick+5) | (m->oargs.nfargs > 9) |
656 greg 2.1 (m->oargs.nfargs % 3))
657     objerror(m, USER, "bad # arguments");
658 greg 2.6 /* record surface struck */
659 greg 2.9 hitfront = (r->rod > 0);
660 greg 2.1 /* load cal file */
661 greg 2.50 mf = hasthick ? getfunc(m, 5, 0x1d, 1)
662     : getfunc(m, 4, 0xe, 1) ;
663 greg 2.25 setfunc(m, r);
664 greg 2.50 nd.thick = 0; /* set thickness */
665     if (hasthick) {
666     nd.thick = evalue(mf->ep[0]);
667     if ((-FTINY <= nd.thick) & (nd.thick <= FTINY))
668     nd.thick = 0;
669     }
670 greg 2.26 /* check backface visibility */
671     if (!hitfront & !backvis) {
672     raytrans(r);
673     return(1);
674     }
675 greg 2.5 /* check other rays to pass */
676 greg 2.34 if (nd.thick != 0 && (r->crtype & SHADOW ||
677     !(r->crtype & (SPECULAR|AMBIENT)) ||
678 greg 2.29 (nd.thick > 0) ^ hitfront)) {
679 greg 2.5 raytrans(r); /* hide our proxy */
680 greg 2.1 return(1);
681     }
682 greg 2.51 if (hasthick && r->crtype & SHADOW) /* early shadow check #1 */
683     return(1);
684 greg 2.31 nd.mp = m;
685     nd.pr = r;
686 greg 2.5 /* get BSDF data */
687 greg 2.50 nd.sd = loadBSDF(m->oargs.sarg[hasthick]);
688 greg 2.51 /* early shadow check #2 */
689 greg 2.55 if (r->crtype & SHADOW && (nd.sd->tf == NULL) & (nd.sd->tb == NULL)) {
690     SDfreeCache(nd.sd);
691 greg 2.34 return(1);
692 greg 2.55 }
693 greg 2.63 /* diffuse components */
694 greg 2.6 if (hitfront) {
695 greg 2.31 cvt_sdcolor(nd.rdiff, &nd.sd->rLambFront);
696     if (m->oargs.nfargs >= 3) {
697     setcolor(ctmp, m->oargs.farg[0],
698 greg 2.1 m->oargs.farg[1],
699     m->oargs.farg[2]);
700 greg 2.31 addcolor(nd.rdiff, ctmp);
701     }
702 greg 2.63 cvt_sdcolor(nd.tdiff, &nd.sd->tLambFront);
703 greg 2.1 } else {
704 greg 2.31 cvt_sdcolor(nd.rdiff, &nd.sd->rLambBack);
705     if (m->oargs.nfargs >= 6) {
706     setcolor(ctmp, m->oargs.farg[3],
707 greg 2.1 m->oargs.farg[4],
708     m->oargs.farg[5]);
709 greg 2.31 addcolor(nd.rdiff, ctmp);
710     }
711 greg 2.63 cvt_sdcolor(nd.tdiff, &nd.sd->tLambBack);
712 greg 2.1 }
713 greg 2.63 if (m->oargs.nfargs >= 9) { /* add diffuse transmittance? */
714 greg 2.31 setcolor(ctmp, m->oargs.farg[6],
715 greg 2.1 m->oargs.farg[7],
716     m->oargs.farg[8]);
717 greg 2.31 addcolor(nd.tdiff, ctmp);
718     }
719 greg 2.1 /* get modifiers */
720     raytexture(r, m->omod);
721     /* modify diffuse values */
722     multcolor(nd.rdiff, r->pcol);
723     multcolor(nd.tdiff, r->pcol);
724     /* get up vector */
725 greg 2.50 upvec[0] = evalue(mf->ep[hasthick+0]);
726     upvec[1] = evalue(mf->ep[hasthick+1]);
727     upvec[2] = evalue(mf->ep[hasthick+2]);
728 greg 2.1 /* return to world coords */
729 greg 2.21 if (mf->fxp != &unitxf) {
730     multv3(upvec, upvec, mf->fxp->xfm);
731     nd.thick *= mf->fxp->sca;
732 greg 2.1 }
733 greg 2.23 if (r->rox != NULL) {
734     multv3(upvec, upvec, r->rox->f.xfm);
735     nd.thick *= r->rox->f.sca;
736     }
737 greg 2.1 raynormal(nd.pnorm, r);
738     /* compute local BSDF xform */
739     ec = SDcompXform(nd.toloc, nd.pnorm, upvec);
740     if (!ec) {
741 greg 2.4 nd.vray[0] = -r->rdir[0];
742     nd.vray[1] = -r->rdir[1];
743     nd.vray[2] = -r->rdir[2];
744     ec = SDmapDir(nd.vray, nd.toloc, nd.vray);
745 greg 2.20 }
746 greg 2.19 if (ec) {
747     objerror(m, WARNING, "Illegal orientation vector");
748 greg 2.55 SDfreeCache(nd.sd);
749 greg 2.19 return(1);
750 greg 2.1 }
751 greg 2.50 setcolor(nd.cthru, 0, 0, 0); /* consider through component */
752 greg 2.61 setcolor(nd.cthru_surr, 0, 0, 0);
753 greg 2.52 if (m->otype == MAT_ABSDF) {
754 greg 2.50 compute_through(&nd);
755     if (r->crtype & SHADOW) {
756     RAY tr; /* attempt to pass shadow ray */
757 greg 2.55 SDfreeCache(nd.sd);
758 greg 2.50 if (rayorigin(&tr, TRANS, r, nd.cthru) < 0)
759     return(1); /* no through component */
760     VCOPY(tr.rdir, r->rdir);
761     rayvalue(&tr); /* transmit with scaling */
762     multcolor(tr.rcol, tr.rcoef);
763     copycolor(r->rcol, tr.rcol);
764     return(1); /* we're done */
765     }
766 greg 2.34 }
767     ec = SDinvXform(nd.fromloc, nd.toloc);
768     if (!ec) /* determine BSDF resolution */
769     ec = SDsizeBSDF(nd.sr_vpsa, nd.vray, NULL,
770     SDqueryMin+SDqueryMax, nd.sd);
771 greg 2.20 if (ec)
772     objerror(m, USER, transSDError(ec));
773    
774 greg 2.9 nd.sr_vpsa[0] = sqrt(nd.sr_vpsa[0]);
775     nd.sr_vpsa[1] = sqrt(nd.sr_vpsa[1]);
776 greg 2.6 if (!hitfront) { /* perturb normal towards hit */
777 greg 2.1 nd.pnorm[0] = -nd.pnorm[0];
778     nd.pnorm[1] = -nd.pnorm[1];
779     nd.pnorm[2] = -nd.pnorm[2];
780     }
781     /* sample reflection */
782     sample_sdf(&nd, SDsampSpR);
783     /* sample transmission */
784     sample_sdf(&nd, SDsampSpT);
785     /* compute indirect diffuse */
786 greg 2.39 copycolor(ctmp, nd.rdiff);
787     addcolor(ctmp, nd.runsamp);
788     if (bright(ctmp) > FTINY) { /* ambient from reflection */
789 greg 2.6 if (!hitfront)
790 greg 2.1 flipsurface(r);
791     multambient(ctmp, r, nd.pnorm);
792     addcolor(r->rcol, ctmp);
793 greg 2.6 if (!hitfront)
794 greg 2.1 flipsurface(r);
795     }
796 greg 2.39 copycolor(ctmp, nd.tdiff);
797     addcolor(ctmp, nd.tunsamp);
798     if (bright(ctmp) > FTINY) { /* ambient from other side */
799 greg 2.1 FVECT bnorm;
800 greg 2.6 if (hitfront)
801 greg 2.1 flipsurface(r);
802     bnorm[0] = -nd.pnorm[0];
803     bnorm[1] = -nd.pnorm[1];
804     bnorm[2] = -nd.pnorm[2];
805 greg 2.9 if (nd.thick != 0) { /* proxy with offset? */
806 greg 2.5 VCOPY(vtmp, r->rop);
807 greg 2.18 VSUM(r->rop, vtmp, r->ron, nd.thick);
808 greg 2.5 multambient(ctmp, r, bnorm);
809     VCOPY(r->rop, vtmp);
810     } else
811     multambient(ctmp, r, bnorm);
812 greg 2.1 addcolor(r->rcol, ctmp);
813 greg 2.6 if (hitfront)
814 greg 2.1 flipsurface(r);
815     }
816     /* add direct component */
817 greg 2.22 if ((bright(nd.tdiff) <= FTINY) & (nd.sd->tf == NULL) &
818     (nd.sd->tb == NULL)) {
819 greg 2.5 direct(r, dir_brdf, &nd); /* reflection only */
820 greg 2.9 } else if (nd.thick == 0) {
821 greg 2.5 direct(r, dir_bsdf, &nd); /* thin surface scattering */
822     } else {
823     direct(r, dir_brdf, &nd); /* reflection first */
824     VCOPY(vtmp, r->rop); /* offset for transmitted */
825     VSUM(r->rop, vtmp, r->ron, -nd.thick);
826 greg 2.6 direct(r, dir_btdf, &nd); /* separate transmission */
827 greg 2.5 VCOPY(r->rop, vtmp);
828     }
829 greg 2.1 /* clean up */
830     SDfreeCache(nd.sd);
831     return(1);
832     }