diff --git a/m_bsdf.c b/m_bsdf.c index fd24770..78d5e9d 100644 --- a/m_bsdf.c +++ b/m_bsdf.c @@ -1,813 +1,818 @@ #ifndef lint static const char RCSid[] = "$Id: m_bsdf.c,v 2.70 2022/03/09 00:27:25 greg Exp $"; #endif /* * Shading for materials with BSDFs taken from XML data files */ #include "copyright.h" #include "ray.h" #include "otypes.h" #include "ambient.h" #include "source.h" #include "func.h" #include "bsdf.h" #include "random.h" #include "pmapmat.h" /* * Arguments to this material include optional diffuse colors. * String arguments include the BSDF and function files. * For the MAT_BSDF type, a non-zero thickness causes the useful behavior * of translating transmitted rays this distance beneath the surface * (opposite the surface normal) to bypass any intervening geometry. * Translation only affects scattered, non-source-directed samples. * A non-zero thickness has the further side-effect that an unscattered * (view) ray will pass right through our material, making the BSDF * surface invisible and showing the proxied geometry instead. Thickness * has the further effect of turning off reflection on the reverse side so * rays heading in the opposite direction pass unimpeded through the BSDF * surface. A paired surface may be placed on the opposide side of * the detail geometry, less than this thickness away, if a two-way * proxy is desired. Note that the sign of the thickness is important. * A positive thickness hides geometry behind the BSDF surface and uses * front reflectance and transmission properties. A negative thickness * hides geometry in front of the surface when rays hit from behind, * and applies only the transmission and backside reflectance properties. * Reflection is ignored on the hidden side, as those rays pass through. * For the MAT_ABSDF type, we check for a strong "through" component. * Such a component will cause direct rays to pass through unscattered. * A separate test prevents over-counting by dropping samples that are * too close to this "through" direction. BSDFs with such a through direction * will also have a view component, meaning they are somewhat see-through. * A MAT_BSDF type with zero thickness behaves the same as a MAT_ABSDF * type with no strong through component. * The "up" vector for the BSDF is given by three variables, defined * (along with the thickness) by the named function file, or '.' if none. * Together with the surface normal, this defines the local coordinate * system for the BSDF. * We do not reorient the surface, so if the BSDF has no back-side * reflectance and none is given in the real arguments, a BSDF surface * with zero thickness will appear black when viewed from behind * unless backface visibility is on, when it becomes invisible. * The diffuse arguments are added to components in the BSDF file, * not multiplied. However, patterns affect this material as a multiplier * on everything except non-diffuse reflection. * * Arguments for MAT_ABSDF are: * 5+ BSDFfile ux uy uz funcfile transform * 0 * 0|3|6|9 rdf gdf bdf * rdb gdb bdb * rdt gdt bdt * * Arguments for MAT_BSDF are: * 6+ thick BSDFfile ux uy uz funcfile transform * 0 * 0|3|6|9 rdf gdf bdf * rdb gdb bdb * rdt gdt bdt */ /* * Note that our reverse ray-tracing process means that the positions * of incoming and outgoing vectors may be reversed in our calls * to the BSDF library. This is usually fine, since the bidirectional nature * of the BSDF (that's what the 'B' stands for) means it all works out. */ typedef struct { OBJREC *mp; /* material pointer */ RAY *pr; /* intersected ray */ FVECT pnorm; /* perturbed surface normal */ FVECT vray; /* local outgoing (return) vector */ double sr_vpsa[2]; /* sqrt of BSDF projected solid angle extrema */ RREAL toloc[3][3]; /* world to local BSDF coords */ RREAL fromloc[3][3]; /* local BSDF coords to world */ double thick; /* surface thickness */ COLOR cthru; /* "through" component for MAT_ABSDF */ COLOR cthru_surr; /* surround for "through" component */ SDData *sd; /* loaded BSDF data */ COLOR rdiff; /* diffuse reflection */ COLOR runsamp; /* BSDF hemispherical reflection */ COLOR tdiff; /* diffuse transmission */ COLOR tunsamp; /* BSDF hemispherical transmission */ } BSDFDAT; /* BSDF material data */ #define cvt_sdcolor(cv, svp) ccy2rgb(&(svp)->spec, (svp)->cieY, cv) typedef struct { double vy; /* brightness (for sorting) */ FVECT tdir; /* through sample direction (normalized) */ COLOR vcol; /* BTDF color */ } PEAKSAMP; /* BTDF peak sample */ /* Comparison function to put near-peak values in descending order */ static int cmp_psamp(const void *p1, const void *p2) { double diff = (*(const PEAKSAMP *)p1).vy - (*(const PEAKSAMP *)p2).vy; if (diff > 0) return(-1); if (diff < 0) return(1); return(0); } /* Compute "through" component color for MAT_ABSDF */ static void compute_through(BSDFDAT *ndp) { #define NDIR2CHECK 29 static const float dir2check[NDIR2CHECK][2] = { {0, 0}, {-0.6, 0}, {0, 0.6}, {0, -0.6}, {0.6, 0}, {-0.6, 0.6}, {-0.6, -0.6}, {0.6, 0.6}, {0.6, -0.6}, {-1.2, 0}, {0, 1.2}, {0, -1.2}, {1.2, 0}, {-1.2, 1.2}, {-1.2, -1.2}, {1.2, 1.2}, {1.2, -1.2}, {-1.8, 0}, {0, 1.8}, {0, -1.8}, {1.8, 0}, {-1.8, 1.8}, {-1.8, -1.8}, {1.8, 1.8}, {1.8, -1.8}, {-2.4, 0}, {0, 2.4}, {0, -2.4}, {2.4, 0}, }; PEAKSAMP psamp[NDIR2CHECK]; SDSpectralDF *dfp; FVECT pdir; double tomega, srchrad; double tomsum, tomsurr; COLOR vpeak, vsurr, btdiff; double vypeak; int i, ns; SDError ec; if (ndp->pr->rod > 0) dfp = (ndp->sd->tf != NULL) ? ndp->sd->tf : ndp->sd->tb; else dfp = (ndp->sd->tb != NULL) ? ndp->sd->tb : ndp->sd->tf; if (dfp == NULL) return; /* no specular transmission */ if (bright(ndp->pr->pcol) <= FTINY) return; /* pattern is black, here */ srchrad = sqrt(dfp->minProjSA); /* else evaluate peak */ for (i = 0; i < NDIR2CHECK; i++) { SDValue sv; psamp[i].tdir[0] = -ndp->vray[0] + dir2check[i][0]*srchrad; psamp[i].tdir[1] = -ndp->vray[1] + dir2check[i][1]*srchrad; psamp[i].tdir[2] = -ndp->vray[2]; normalize(psamp[i].tdir); ec = SDevalBSDF(&sv, ndp->vray, psamp[i].tdir, ndp->sd); if (ec) goto baderror; cvt_sdcolor(psamp[i].vcol, &sv); psamp[i].vy = sv.cieY; } qsort(psamp, NDIR2CHECK, sizeof(PEAKSAMP), cmp_psamp); if (psamp[0].vy <= FTINY) return; /* zero BTDF here */ setcolor(vpeak, 0, 0, 0); setcolor(vsurr, 0, 0, 0); vypeak = tomsum = tomsurr = 0; /* combine top unique values */ ns = 0; for (i = 0; i < NDIR2CHECK; i++) { if (i && psamp[i].vy == psamp[i-1].vy) continue; /* assume duplicate sample */ ec = SDsizeBSDF(&tomega, ndp->vray, psamp[i].tdir, SDqueryMin, ndp->sd); if (ec) goto baderror; scalecolor(psamp[i].vcol, tomega); /* not part of peak? */ if (tomega > 1.5*dfp->minProjSA || vypeak > 8.*psamp[i].vy*ns) { if (!i) return; /* abort */ addcolor(vsurr, psamp[i].vcol); tomsurr += tomega; continue; } addcolor(vpeak, psamp[i].vcol); tomsum += tomega; vypeak += psamp[i].vy; ++ns; } if (tomsurr < 0.2*tomsum) /* insufficient surround? */ return; scalecolor(vsurr, 1./tomsurr); /* surround is avg. BTDF */ if (ndp->vray[2] > 0) /* get diffuse BTDF */ cvt_sdcolor(btdiff, &ndp->sd->tLambFront); else cvt_sdcolor(btdiff, &ndp->sd->tLambBack); scalecolor(btdiff, (1./PI)); for (i = 3; i--; ) { /* remove diffuse contrib. */ if ((colval(vpeak,i) -= tomsum*colval(btdiff,i)) < 0) colval(vpeak,i) = 0; if ((colval(vsurr,i) -= colval(btdiff,i)) < 0) colval(vsurr,i) = 0; } if (bright(vpeak) < .0005) /* < 0.05% specular? */ return; multcolor(vsurr, ndp->pr->pcol); /* modify by color */ multcolor(vpeak, ndp->pr->pcol); copycolor(ndp->cthru, vpeak); copycolor(ndp->cthru_surr, vsurr); return; baderror: objerror(ndp->mp, USER, transSDError(ec)); #undef NDIR2CHECK } /* Jitter ray sample according to projected solid angle and specjitter */ static void bsdf_jitter(FVECT vres, BSDFDAT *ndp, double sr_psa) { VCOPY(vres, ndp->vray); if (specjitter < 1.) sr_psa *= specjitter; if (sr_psa <= FTINY) return; vres[0] += sr_psa*(.5 - frandom()); vres[1] += sr_psa*(.5 - frandom()); normalize(vres); } /* Get BSDF specular for direct component, returning true if OK to proceed */ static int direct_specular_OK(COLOR cval, FVECT ldir, double omega, BSDFDAT *ndp) { int nsamp = 1; int scnt = 0; FVECT vsrc, vjit; double tomega, tomega2; double sf, tsr, sd[2]; COLOR csmp, cdiff; double diffY; SDValue sv; SDError ec; int i; /* in case we fail */ setcolor(cval, 0, 0, 0); /* transform source direction */ if (SDmapDir(vsrc, ndp->toloc, ldir) != SDEnone) return(0); /* check indirect over-counting */ if ((vsrc[2] > 0) ^ (ndp->vray[2] > 0) && bright(ndp->cthru) > FTINY) { double dx = vsrc[0] + ndp->vray[0]; double dy = vsrc[1] + ndp->vray[1]; SDSpectralDF *dfp = (ndp->pr->rod > 0) ? ((ndp->sd->tf != NULL) ? ndp->sd->tf : ndp->sd->tb) : ((ndp->sd->tb != NULL) ? ndp->sd->tb : ndp->sd->tf) ; tomega = omega*fabs(vsrc[2]); if (dx*dx + dy*dy <= (2.5*4./PI)*(tomega + dfp->minProjSA + 2.*sqrt(tomega*dfp->minProjSA))) { if (bright(ndp->cthru_surr) <= FTINY) return(0); copycolor(cval, ndp->cthru_surr); return(1); /* return non-zero surround BTDF */ } } /* will discount diffuse portion */ switch ((vsrc[2] > 0)<<1 | (ndp->vray[2] > 0)) { case 3: if (ndp->sd->rf == NULL) return(0); /* all diffuse */ sv = ndp->sd->rLambFront; break; case 0: if (ndp->sd->rb == NULL) return(0); /* all diffuse */ sv = ndp->sd->rLambBack; break; case 1: if ((ndp->sd->tf == NULL) & (ndp->sd->tb == NULL)) return(0); /* all diffuse */ sv = ndp->sd->tLambFront; break; case 2: if ((ndp->sd->tf == NULL) & (ndp->sd->tb == NULL)) return(0); /* all diffuse */ sv = ndp->sd->tLambBack; break; } if (sv.cieY > FTINY) { diffY = sv.cieY *= 1./PI; cvt_sdcolor(cdiff, &sv); } else { diffY = 0; setcolor(cdiff, 0, 0, 0); } ec = SDsizeBSDF(&tomega, ndp->vray, vsrc, SDqueryMin, ndp->sd); if (ec) goto baderror; /* check if sampling BSDF */ if ((tsr = sqrt(tomega)) > 0) { nsamp = 4.*specjitter*ndp->pr->rweight + .5; nsamp += !nsamp; } /* jitter to fuzz BSDF cells */ for (i = nsamp; i--; ) { bsdf_jitter(vjit, ndp, tsr); /* compute BSDF */ ec = SDevalBSDF(&sv, vjit, vsrc, ndp->sd); if (ec) goto baderror; if (sv.cieY - diffY <= FTINY) continue; /* no specular part */ /* check for variable resolution */ ec = SDsizeBSDF(&tomega2, vjit, vsrc, SDqueryMin, ndp->sd); if (ec) goto baderror; if (tomega2 < .12*tomega) continue; /* not safe to include */ cvt_sdcolor(csmp, &sv); addcolor(cval, csmp); ++scnt; } if (!scnt) /* no valid specular samples? */ return(0); sf = 1./scnt; /* weighted average BSDF */ scalecolor(cval, sf); /* subtract diffuse contribution */ for (i = 3*(diffY > FTINY); i--; ) if ((colval(cval,i) -= colval(cdiff,i)) < 0) colval(cval,i) = 0; return(1); baderror: objerror(ndp->mp, USER, transSDError(ec)); return(0); /* gratis return */ } /* Compute source contribution for BSDF (reflected & transmitted) */ static void dir_bsdf( COLOR cval, /* returned coefficient */ void *nnp, /* material data */ FVECT ldir, /* light source direction */ double omega /* light source size */ ) { BSDFDAT *np = (BSDFDAT *)nnp; double ldot; double dtmp; COLOR ctmp; setcolor(cval, 0, 0, 0); ldot = DOT(np->pnorm, ldir); if ((-FTINY <= ldot) & (ldot <= FTINY)) return; if (ldot > 0 && bright(np->rdiff) > FTINY) { /* * Compute diffuse reflected component */ copycolor(ctmp, np->rdiff); dtmp = ldot * omega * (1./PI); scalecolor(ctmp, dtmp); addcolor(cval, ctmp); } if (ldot < 0 && bright(np->tdiff) > FTINY) { /* * Compute diffuse transmission */ copycolor(ctmp, np->tdiff); dtmp = -ldot * omega * (1./PI); scalecolor(ctmp, dtmp); addcolor(cval, ctmp); } if (ambRayInPmap(np->pr)) return; /* specular already in photon map */ /* * Compute specular scattering coefficient using BSDF */ if (!direct_specular_OK(ctmp, ldir, omega, np)) return; if (ldot < 0) { /* pattern for specular transmission */ multcolor(ctmp, np->pr->pcol); dtmp = -ldot * omega; } else dtmp = ldot * omega; scalecolor(ctmp, dtmp); addcolor(cval, ctmp); } /* Compute source contribution for BSDF (reflected only) */ static void dir_brdf( COLOR cval, /* returned coefficient */ void *nnp, /* material data */ FVECT ldir, /* light source direction */ double omega /* light source size */ ) { BSDFDAT *np = (BSDFDAT *)nnp; double ldot; double dtmp; COLOR ctmp, ctmp1, ctmp2; setcolor(cval, 0, 0, 0); ldot = DOT(np->pnorm, ldir); - + if (ldot <= FTINY) return; if (bright(np->rdiff) > FTINY) { /* * Compute diffuse reflected component */ copycolor(ctmp, np->rdiff); dtmp = ldot * omega * (1./PI); scalecolor(ctmp, dtmp); addcolor(cval, ctmp); } if (ambRayInPmap(np->pr)) return; /* specular already in photon map */ /* * Compute specular reflection coefficient using BSDF */ if (!direct_specular_OK(ctmp, ldir, omega, np)) return; dtmp = ldot * omega; scalecolor(ctmp, dtmp); addcolor(cval, ctmp); } /* Compute source contribution for BSDF (transmitted only) */ static void dir_btdf( COLOR cval, /* returned coefficient */ void *nnp, /* material data */ FVECT ldir, /* light source direction */ double omega /* light source size */ ) { BSDFDAT *np = (BSDFDAT *)nnp; double ldot; double dtmp; COLOR ctmp; setcolor(cval, 0, 0, 0); ldot = DOT(np->pnorm, ldir); if (ldot >= -FTINY) return; if (bright(np->tdiff) > FTINY) { /* * Compute diffuse transmission */ copycolor(ctmp, np->tdiff); dtmp = -ldot * omega * (1./PI); scalecolor(ctmp, dtmp); addcolor(cval, ctmp); } if (ambRayInPmap(np->pr)) return; /* specular already in photon map */ /* * Compute specular scattering coefficient using BSDF */ if (!direct_specular_OK(ctmp, ldir, omega, np)) return; /* full pattern on transmission */ multcolor(ctmp, np->pr->pcol); dtmp = -ldot * omega; scalecolor(ctmp, dtmp); addcolor(cval, ctmp); } /* Sample separate BSDF component */ static int sample_sdcomp(BSDFDAT *ndp, SDComponent *dcp, int xmit) { const int hasthru = (xmit && !(ndp->pr->crtype & (SPECULAR|AMBIENT)) && bright(ndp->cthru) > FTINY); int nstarget = 1; int nsent = 0; int n; SDError ec; SDValue bsv; double xrand; FVECT vsmp, vinc; RAY sr; /* multiple samples? */ if (specjitter > 1.5) { nstarget = specjitter*ndp->pr->rweight + .5; nstarget += !nstarget; } /* run through our samples */ for (n = 0; n < nstarget; n++) { if (nstarget == 1) { /* stratify random variable */ xrand = urand(ilhash(dimlist,ndims)+samplendx); if (specjitter < 1.) xrand = .5 + specjitter*(xrand-.5); } else { xrand = (n + frandom())/(double)nstarget; } SDerrorDetail[0] = '\0'; /* sample direction & coef. */ bsdf_jitter(vsmp, ndp, ndp->sr_vpsa[0]); VCOPY(vinc, vsmp); /* to compare after */ ec = SDsampComponent(&bsv, vsmp, xrand, dcp); if (ec) objerror(ndp->mp, USER, transSDError(ec)); if (bsv.cieY <= FTINY) /* zero component? */ break; if (hasthru) { /* check for view ray */ double dx = vinc[0] + vsmp[0]; double dy = vinc[1] + vsmp[1]; if (dx*dx + dy*dy <= ndp->sr_vpsa[0]*ndp->sr_vpsa[0]) continue; /* exclude view sample */ } /* map non-view sample->world */ if (SDmapDir(sr.rdir, ndp->fromloc, vsmp) != SDEnone) break; /* spawn a specular ray */ if (nstarget > 1) bsv.cieY /= (double)nstarget; cvt_sdcolor(sr.rcoef, &bsv); /* use sample color */ if (xmit) /* apply pattern on transmit */ multcolor(sr.rcoef, ndp->pr->pcol); if (rayorigin(&sr, SPECULAR, ndp->pr, sr.rcoef) < 0) { if (!n & (nstarget > 1)) { n = nstarget; /* avoid infinitue loop */ nstarget = nstarget*sr.rweight/minweight; if (n == nstarget) break; n = -1; /* moved target */ } continue; /* try again */ } if (xmit && ndp->thick != 0) /* need to offset origin? */ VSUM(sr.rorg, sr.rorg, ndp->pr->ron, -ndp->thick); rayvalue(&sr); /* send & evaluate sample */ multcolor(sr.rcol, sr.rcoef); addcolor(ndp->pr->rcol, sr.rcol); ++nsent; } return(nsent); } /* Sample non-diffuse components of BSDF */ static int sample_sdf(BSDFDAT *ndp, int sflags) { int hasthru = (sflags == SDsampSpT && !(ndp->pr->crtype & (SPECULAR|AMBIENT)) && bright(ndp->cthru) > FTINY); int n, ntotal = 0; double b = 0; SDSpectralDF *dfp; COLORV *unsc; if (sflags == SDsampSpT) { unsc = ndp->tunsamp; if (ndp->pr->rod > 0) dfp = (ndp->sd->tf != NULL) ? ndp->sd->tf : ndp->sd->tb; else dfp = (ndp->sd->tb != NULL) ? ndp->sd->tb : ndp->sd->tf; } else /* sflags == SDsampSpR */ { unsc = ndp->runsamp; if (ndp->pr->rod > 0) dfp = ndp->sd->rf; else dfp = ndp->sd->rb; } setcolor(unsc, 0, 0, 0); if (dfp == NULL) /* no specular component? */ return(0); if (hasthru) { /* separate view sample? */ RAY tr; if (rayorigin(&tr, TRANS, ndp->pr, ndp->cthru) == 0) { VCOPY(tr.rdir, ndp->pr->rdir); rayvalue(&tr); multcolor(tr.rcol, tr.rcoef); addcolor(ndp->pr->rcol, tr.rcol); ndp->pr->rxt = ndp->pr->rot + raydistance(&tr); ++ntotal; b = bright(ndp->cthru); } else hasthru = 0; } if (dfp->maxHemi - b <= FTINY) { /* have specular to sample? */ b = 0; } else { FVECT vjit; bsdf_jitter(vjit, ndp, ndp->sr_vpsa[1]); b = SDdirectHemi(vjit, sflags, ndp->sd) - b; if (b < 0) b = 0; } if (b <= specthresh+FTINY) { /* below sampling threshold? */ if (b > FTINY) { /* XXX no color from BSDF */ if (sflags == SDsampSpT) { copycolor(unsc, ndp->pr->pcol); scalecolor(unsc, b); } else /* no pattern on reflection */ setcolor(unsc, b, b, b); } return(ntotal); } dimlist[ndims] = (int)(size_t)ndp->mp; /* else sample specular */ ndims += 2; for (n = dfp->ncomp; n--; ) { /* loop over components */ dimlist[ndims-1] = n + 9438; ntotal += sample_sdcomp(ndp, &dfp->comp[n], sflags==SDsampSpT); } ndims -= 2; return(ntotal); } /* Color a ray that hit a BSDF material */ int m_bsdf(OBJREC *m, RAY *r) { int hasthick = (m->otype == MAT_BSDF); int hitfront; COLOR ctmp; SDError ec; FVECT upvec, vtmp; MFUNC *mf; BSDFDAT nd; /* check arguments */ if ((m->oargs.nsargs < hasthick+5) | (m->oargs.nfargs > 9) | (m->oargs.nfargs % 3)) objerror(m, USER, "bad # arguments"); /* record surface struck */ hitfront = (r->rod > 0); /* load cal file */ mf = hasthick ? getfunc(m, 5, 0x1d, 1) : getfunc(m, 4, 0xe, 1) ; setfunc(m, r); nd.thick = 0; /* set thickness */ if (hasthick) { nd.thick = evalue(mf->ep[0]); if ((-FTINY <= nd.thick) & (nd.thick <= FTINY)) nd.thick = 0; } /* check backface visibility */ if (!hitfront & !backvis) { raytrans(r); return(1); } + /* block shadow rays when photons are directly visualized + (global pmap and -ab -1, or caustic or contrib pmaps) */ + if ( ( (globalPmap && ambounce<0) || causticPmap || contribPmap ) + && (r->crtype & SHADOW) ) + return(1); /* check other rays to pass */ if (nd.thick != 0 && (r->crtype & SHADOW || !(r->crtype & (SPECULAR|AMBIENT)) || (nd.thick > 0) ^ hitfront)) { raytrans(r); /* hide our proxy */ return(1); } if (hasthick && r->crtype & SHADOW) /* early shadow check #1 */ return(1); nd.mp = m; nd.pr = r; /* get BSDF data */ nd.sd = loadBSDF(m->oargs.sarg[hasthick]); /* early shadow check #2 */ if (r->crtype & SHADOW && (nd.sd->tf == NULL) & (nd.sd->tb == NULL)) { SDfreeCache(nd.sd); return(1); } /* diffuse components */ if (hitfront) { cvt_sdcolor(nd.rdiff, &nd.sd->rLambFront); if (m->oargs.nfargs >= 3) { setcolor(ctmp, m->oargs.farg[0], m->oargs.farg[1], m->oargs.farg[2]); addcolor(nd.rdiff, ctmp); } cvt_sdcolor(nd.tdiff, &nd.sd->tLambFront); } else { cvt_sdcolor(nd.rdiff, &nd.sd->rLambBack); if (m->oargs.nfargs >= 6) { setcolor(ctmp, m->oargs.farg[3], m->oargs.farg[4], m->oargs.farg[5]); addcolor(nd.rdiff, ctmp); } cvt_sdcolor(nd.tdiff, &nd.sd->tLambBack); } if (m->oargs.nfargs >= 9) { /* add diffuse transmittance? */ setcolor(ctmp, m->oargs.farg[6], m->oargs.farg[7], m->oargs.farg[8]); addcolor(nd.tdiff, ctmp); } /* get modifiers */ raytexture(r, m->omod); /* modify diffuse values */ multcolor(nd.rdiff, r->pcol); multcolor(nd.tdiff, r->pcol); /* get up vector */ upvec[0] = evalue(mf->ep[hasthick+0]); upvec[1] = evalue(mf->ep[hasthick+1]); upvec[2] = evalue(mf->ep[hasthick+2]); /* return to world coords */ if (mf->fxp != &unitxf) { multv3(upvec, upvec, mf->fxp->xfm); nd.thick *= mf->fxp->sca; } if (r->rox != NULL) { multv3(upvec, upvec, r->rox->f.xfm); nd.thick *= r->rox->f.sca; } raynormal(nd.pnorm, r); /* compute local BSDF xform */ ec = SDcompXform(nd.toloc, nd.pnorm, upvec); if (!ec) { nd.vray[0] = -r->rdir[0]; nd.vray[1] = -r->rdir[1]; nd.vray[2] = -r->rdir[2]; ec = SDmapDir(nd.vray, nd.toloc, nd.vray); } if (ec) { objerror(m, WARNING, "Illegal orientation vector"); SDfreeCache(nd.sd); return(1); } setcolor(nd.cthru, 0, 0, 0); /* consider through component */ setcolor(nd.cthru_surr, 0, 0, 0); if (m->otype == MAT_ABSDF) { compute_through(&nd); if (r->crtype & SHADOW) { RAY tr; /* attempt to pass shadow ray */ SDfreeCache(nd.sd); if (rayorigin(&tr, TRANS, r, nd.cthru) < 0) return(1); /* no through component */ VCOPY(tr.rdir, r->rdir); rayvalue(&tr); /* transmit with scaling */ multcolor(tr.rcol, tr.rcoef); copycolor(r->rcol, tr.rcol); return(1); /* we're done */ } } ec = SDinvXform(nd.fromloc, nd.toloc); if (!ec) /* determine BSDF resolution */ ec = SDsizeBSDF(nd.sr_vpsa, nd.vray, NULL, SDqueryMin+SDqueryMax, nd.sd); if (ec) objerror(m, USER, transSDError(ec)); nd.sr_vpsa[0] = sqrt(nd.sr_vpsa[0]); nd.sr_vpsa[1] = sqrt(nd.sr_vpsa[1]); if (!hitfront) { /* perturb normal towards hit */ nd.pnorm[0] = -nd.pnorm[0]; nd.pnorm[1] = -nd.pnorm[1]; nd.pnorm[2] = -nd.pnorm[2]; } /* sample reflection */ sample_sdf(&nd, SDsampSpR); /* sample transmission */ sample_sdf(&nd, SDsampSpT); /* compute indirect diffuse */ copycolor(ctmp, nd.rdiff); addcolor(ctmp, nd.runsamp); if (bright(ctmp) > FTINY) { /* ambient from reflection */ if (!hitfront) flipsurface(r); multambient(ctmp, r, nd.pnorm); addcolor(r->rcol, ctmp); if (!hitfront) flipsurface(r); } copycolor(ctmp, nd.tdiff); addcolor(ctmp, nd.tunsamp); if (bright(ctmp) > FTINY) { /* ambient from other side */ FVECT bnorm; if (hitfront) flipsurface(r); bnorm[0] = -nd.pnorm[0]; bnorm[1] = -nd.pnorm[1]; bnorm[2] = -nd.pnorm[2]; if (nd.thick != 0) { /* proxy with offset? */ VCOPY(vtmp, r->rop); VSUM(r->rop, vtmp, r->ron, nd.thick); multambient(ctmp, r, bnorm); VCOPY(r->rop, vtmp); } else multambient(ctmp, r, bnorm); addcolor(r->rcol, ctmp); if (hitfront) flipsurface(r); } /* add direct component */ if ((bright(nd.tdiff) <= FTINY) & (nd.sd->tf == NULL) & (nd.sd->tb == NULL)) { direct(r, dir_brdf, &nd); /* reflection only */ } else if (nd.thick == 0) { direct(r, dir_bsdf, &nd); /* thin surface scattering */ } else { direct(r, dir_brdf, &nd); /* reflection first */ VCOPY(vtmp, r->rop); /* offset for transmitted */ VSUM(r->rop, vtmp, r->ron, -nd.thick); direct(r, dir_btdf, &nd); /* separate transmission */ VCOPY(r->rop, vtmp); } /* clean up */ SDfreeCache(nd.sd); return(1); } diff --git a/pmapdata.c b/pmapdata.c index bd25b4b..4893339 100644 --- a/pmapdata.c +++ b/pmapdata.c @@ -1,979 +1,979 @@ #ifndef lint static const char RCSid[] = "$Id: pmapdata.c,v 2.23 2021/04/14 11:26:25 rschregle Exp $"; #endif -/* +/* ========================================================================= Photon map types and high level interface to nearest neighbour lookups in underlying point cloud data structure. - + The default data structure is an in-core kd-tree (see pmapkdt.{h,c}). This can be overriden with the PMAP_OOC compiletime switch, which enables an out-of-core octree (see oococt.{h,c}). Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) (c) Fraunhofer Institute for Solar Energy Systems, - supported by the German Research Foundation - (DFG LU-204/10-2, "Fassadenintegrierte Regelsysteme" (FARESYS)) + supported by the German Research Foundation + (DFG LU-204/10-2, "Fassadenintegrierte Regelsysteme" (FARESYS)) (c) Lucerne University of Applied Sciences and Arts, - supported by the Swiss National Science Foundation + supported by the Swiss National Science Foundation (SNSF #147053, "Daylight Redirecting Components", SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment") (c) Tokyo University of Science, - supported by the JSPS Grants-in-Aid for Scientific Research + supported by the JSPS Grants-in-Aid for Scientific Research (KAKENHI JP19KK0115, "Three-Dimensional Light Flow") ========================================================================= - + $Id: pmapdata.c,v 2.23 2021/04/14 11:26:25 rschregle Exp $ */ #include "pmapdata.h" #include "pmapdens.h" #include "pmaprand.h" #include "pmapmat.h" +#include "pmapsrc.h" #include "pmaproi.h" #include "pmapcontrib.h" #include "otypes.h" #include "random.h" PhotonMap *photonMaps [NUM_PMAP_TYPES] = { NULL, NULL, NULL, NULL, NULL, NULL, NULL, #ifdef PMAP_PHOTONFLOW NULL, NULL, #endif NULL }; /* Include routines to handle underlying point cloud data structure */ #ifdef PMAP_OOC #include "pmapooc.c" #else #include "pmapkdt.c" #include "pmaptkdt.c" #endif /* Ambient include/exclude set (from ambient.c) */ #ifndef MAXASET #define MAXASET 4095 #endif extern OBJECT ambset [MAXASET+1]; /* Callback to print photon attributes acc. to user defined format */ int (*printPhoton)(RAY *r, Photon *p, PhotonMap *pm); /* PHOTON MAP BUILD ROUTINES ------------------------------------------ */ void initPhotonMap (PhotonMap *pmap, PhotonMapType t) /* Init photon map 'n' stuff... */ { - if (!pmap) + if (!pmap) return; pmap -> numPhotons = 0; pmap -> biasCompHist = NULL; pmap -> maxPos [0] = pmap -> maxPos [1] = pmap -> maxPos [2] = -FHUGE; pmap -> minPos [0] = pmap -> minPos [1] = pmap -> minPos [2] = FHUGE; pmap -> minGathered = pmap -> maxGathered = pmap -> totalGathered = 0; pmap -> gatherTolerance = gatherTolerance; pmap -> minError = pmap -> maxError = pmap -> rmsError = 0; pmap -> numDensity = 0; pmap -> distribRatio = 1; pmap -> type = t; pmap -> squeue.node = NULL; pmap -> squeue.len = 0; #ifdef PMAP_PATHFILT /* Init path lookup table and its key buffer */ pmap -> pathLUT = NULL; pmap -> pathLUTKeys = NULL; pmap -> numPathLUTKeys = 0; #endif /* Init stats */ setcolor(pmap -> photonFlux, 0, 0, 0); pmap -> CoG [0] = pmap -> CoG [1] = pmap -> CoG [2] = 0; pmap -> CoGdist = 0; pmap -> maxDist2 = 0; - + /* Init transient photon mapping stuff */ pmap -> velocity = photonVelocity; pmap -> minPathLen = pmap -> maxPathLen = pmap -> avgPathLen = 0; /* Init local RNG state */ pmap -> randState [0] = 10243; pmap -> randState [1] = 39829; pmap -> randState [2] = 9433; pmapSeed(randSeed, pmap -> randState); /* Set up type-specific photon lookup callback */ pmap -> lookup = pmapLookup [t]; /* Init storage */ pmap -> heap = NULL; pmap -> heapBuf = NULL; pmap -> heapBufLen = 0; /* Init precomputed contrib photon stuff */ pmap -> contribMode = 0; pmap -> preCompContribTab = NULL; pmap -> contribHeap = NULL; pmap -> contribHeapBuf = NULL; pmap -> preCompContrib = NULL; pmap -> rcOpts = NULL; - pmap -> numContribSrc = 0; + pmap -> numContribSrc = 0; pmap -> contribSrc = NULL; /* Mark photon contribution source as unused */ pmap -> lastContribSrc.srcIdx = pmap -> lastContribSrc.srcBin = -1; #ifdef PMAP_OOC OOC_Null(&pmap -> store); #else kdT_Null(&pmap -> store); #endif } void initPhotonHeap (PhotonMap *pmap) { int fdFlags; - + if (!pmap) error(INTERNAL, "undefined photon map in initPhotonHeap()"); - + if (!pmap -> heap) { /* Open heap file */ mktemp(strcpy(pmap -> heapFname, PMAP_TMPFNAME)); if (!(pmap -> heap = fopen(pmap -> heapFname, "w+b"))) error(SYSTEM, "failed opening heap file in initPhotonHeap()"); #ifdef F_SETFL /* XXX is there an alternative needed for W33nd0z? */ fdFlags = fcntl(fileno(pmap -> heap), F_GETFL); fcntl(fileno(pmap -> heap), F_SETFL, fdFlags | O_APPEND); #endif /* ftruncate(fileno(pmap -> heap), 0); */ } #ifdef PMAP_CONTRIB if (PMAP_CONTRIB_BINNING(pmap)) { /* Allocate precomputed contribution heap buffer */ if (!pmap -> contribHeap) { /* Open heap file for precomputed contributions */ - mktemp(strcpy(pmap -> contribHeapFname, PMAP_TMPFNAME)); + mktemp(strcpy(pmap -> contribHeapFname, PMAP_TMPFNAME)); if (!(pmap -> contribHeap = fopen(pmap -> contribHeapFname, "w+b"))) error(SYSTEM, "failed opening precomputed contribution " "heap file in initPhotonHeap()" ); #ifdef F_SETFL /* XXX is there an alternative needed for W33nd0z? */ fdFlags = fcntl(fileno(pmap -> contribHeap), F_GETFL); fcntl(fileno(pmap -> contribHeap), F_SETFL, fdFlags | O_APPEND); #endif /* ftruncate(fileno(pmap -> contribHeap), 0); */ } } #endif } void flushPhotonHeap (PhotonMap *pmap) { int fd, fdContrib = -1; unsigned long len, lenContrib = 0; if (!pmap) error(INTERNAL, "undefined photon map in flushPhotonHeap()"); if (!pmap -> heap || !pmap -> heapBuf) { - /* Silently ignore undefined heap + /* Silently ignore undefined heap error(INTERNAL, "undefined photon heap in flushPhotonHeap()"); */ return; } fd = fileno(pmap -> heap); len = pmap -> heapBufLen * sizeof(Photon); -#ifdef PMAP_CONTRIB +#ifdef PMAP_CONTRIB if (PMAP_CONTRIB_BINNING(pmap)) { if (!pmap -> contribHeap || !pmap -> contribHeapBuf) { - /* Silently ignore undefined heap + /* Silently ignore undefined heap error(INTERNAL, "undefined contribution heap in flushPhotonHeap()"); */ return; } - /* Do same for the contribution photon buffa if binning is enabled. + /* Do same for the contribution photon buffa if binning is enabled. * Note the photon and contribution heap contain the same number of * entries to ensure photons are flushed in sync along with their * corresponding contributions to their respective heap files */ fdContrib = fileno(pmap -> contribHeap); - lenContrib = pmap -> heapBufLen * + lenContrib = pmap -> heapBufLen * ((PreComputedContrib*)pmap -> preCompContrib) -> contribSize; } #endif - + #if NIX /* Lock heap file(s) prior to writing */ shmLock(fd, F_WRLCK); if (lenContrib) shmLock(fdContrib, F_WRLCK); #endif - + #ifdef DEBUG_PMAP - sprintf(errmsg, "Proc %d: flushing %ld photons from pos %ld\n", + sprintf(errmsg, "Proc %d: flushing %ld photons from pos %ld\n", getpid(), pmap -> heapBufLen, lseek(fd, 0, SEEK_END) / sizeof(Photon) - ); + ); eputs(errmsg); #endif /* Atomically write write block to photon heap file */ /* !!! Unbuffered I/O via pwrite() avoids potential race conditions * !!! and buffer corruption which can occur with lseek()/fseek() - * !!! followed by write()/fwrite(). */ + * !!! followed by write()/fwrite(). */ /*if (pwrite(fd, pmap -> heapBuf, len, lseek(fd, 0, SEEK_END)) != len) */ if (write(fd, pmap -> heapBuf, len) != len) { /* Clean up temp file */ fclose(pmap -> heap); unlink(pmap -> heapFname); error(SYSTEM, "can't append to photon heap file in flushPhotonHeap()"); } - + #ifdef PMAP_CONTRIB /* Atomically write block to contribution heap file */ if (fdContrib > -1 && write(fdContrib, pmap -> contribHeapBuf, lenContrib) != lenContrib ) { /* Clean up temp file */ fclose(pmap -> contribHeap); unlink(pmap -> contribHeapFname); error(SYSTEM, "can't append to contrib heap file in flushPhotonHeap()"); } #endif - + #if NIX if (fsync(fd) || fdContrib > -1 && fsync(fdContrib)) error(SYSTEM, "failed to fsync heap in flushPhotonHeap()"); /* Release lock on heap file(s) */ shmLock(fd, F_UNLCK); if (fdContrib > -1) shmLock(fdContrib, F_UNLCK); #endif - + pmap -> heapBufLen = 0; } #ifdef DEBUG_PMAP static int checkPhotonHeap (FILE *file) /* Check heap for nonsensical or duplicate photons */ { Photon p, lastp; int i, dup; - + rewind(file); memset(&lastp, 0, sizeof(lastp)); - + while (fread(&p, sizeof(p), 1, file)) { dup = 1; - + for (i = 0; i <= 2; i++) { - if (p.pos [i] < thescene.cuorg [i] || + if (p.pos [i] < thescene.cuorg [i] || p.pos [i] > thescene.cuorg [i] + thescene.cusize - ) { - sprintf(errmsg, "corrupt photon in heap at [%f, %f, %f]\n", + ) { + sprintf(errmsg, "corrupt photon in heap at [%f, %f, %f]\n", p.pos [0], p.pos [1], p.pos [2] ); error(WARNING, errmsg); } dup &= p.pos [i] == lastp.pos [i]; } if (dup) { sprintf(errmsg, "consecutive duplicate photon in heap " "at [%f, %f, %f]\n", p.pos [0], p.pos [1], p.pos [2] ); error(WARNING, errmsg); } } return 0; } #endif int newPhoton (PhotonMap* pmap, const RAY* ray, const char *contribs) { unsigned i; Photon photon; COLOR photonFlux; /* Account for distribution ratio */ - if (!pmap || pmapRandom(pmap -> randState) > pmap -> distribRatio) + if (!pmap || pmapRandom(pmap -> randState) > pmap -> distribRatio) return -1; /* Don't store on sources */ - if (ray -> robj > -1 && islight(objptr(ray -> ro -> omod) -> otype)) + if (ray -> robj > -1 && islight(objptr(ray -> ro -> omod) -> otype)) return -1; /* Ignore photon if modifier in/outside exclude/include set */ - if (ambincl != -1 && ray -> ro && + if (ambincl != -1 && ray -> ro && ambincl != inset(ambset, ray->ro->omod) ) return -1; /* Store photon if within a region of interest (for ze Ecksperts!) */ if (!photonInROI(ray)) return -1; /* Adjust flux according to distribution ratio and ray weight */ copycolor(photonFlux, ray -> rcol); #if 0 /* Factored out ray -> rweight as deprecated (?) for pmap, and infact erroneously attenuates volume photon flux based on extinction, which is already factored in by photonParticipate() */ - scalecolor(photonFlux, + scalecolor(photonFlux, ray -> rweight / (pmap -> distribRatio ? pmap -> distribRatio : 1) ); #else - scalecolor(photonFlux, + scalecolor(photonFlux, 1.0 / (pmap -> distribRatio ? pmap -> distribRatio : 1) ); #endif setPhotonFlux(&photon, photonFlux); /* Set photon position and flags */ VCOPY(photon.pos, ray -> rop); photon.flags = 0; photon.caustic = PMAP_CAUSTICRAY(ray); /* Set photon's subprocess index (to consolidate contribution sources after photon distribution completes) */ photon.proc = PMAP_GETRAYPROC(ray); switch (pmap -> type) { #ifdef PMAP_CONTRIB case PMAP_TYPE_CONTRIB: - /* Contrib photon before precomputation (i.e. in forward + /* Contrib photon before precomputation (i.e. in forward pass); set contribution source from last index in contrib - source array. Note the contribution source bin has already + source array. Note the contribution source bin has already been set by newPhotonContribSrc(). */ photon.aux.contribSrc = pmap -> numContribSrc; - /* Photon will be stored; set contribution source index, - * thereby marking it as having spawned photon(s) */ + /* Photon will be stored; set contribution source index, + * thereby marking it as having spawned photon(s) */ if (ray -> rsrc > PMAP_MAXSRCIDX) error(INTERNAL, "contribution source index overflow"); - else pmap -> lastContribSrc.srcIdx = ray -> rsrc; + else pmap -> lastContribSrc.srcIdx = ray -> rsrc; break; case PMAP_TYPE_CONTRIB_CHILD: - /* Contrib photon being precomputed; set contribution + /* Contrib photon being precomputed; set contribution source from index passed in ray. NOTE: This is currently only for info and consistency with the above, but not used by the lookup routines */ if (ray -> rsrc > PMAP_MAXSRCIDX) error(INTERNAL, "contribution source index overflow"); else photon.aux.contribSrc = ray -> rsrc; break; #endif case PMAP_TYPE_TRANSIENT: #ifdef PMAP_PHOTONFLOW case PMAP_TYPE_TRANSLIGHTFLOW: #endif - /* Set cumulative path length for transient photon, corresponding + /* Set cumulative path length for transient photon, corresponding to its time of flight. The cumulative path length is obtained - relative to the maximum path length photonMaxPathDist. + relative to the maximum path length photonMaxPathDist. The last intersection distance is subtracted as this isn't factored in until photonRay() generates a new photon ray. */ photon.aux.pathLen = photonMaxPathDist - ray -> rmax + ray -> rot; break; default: - /* Set photon's path ID from ray serial num + /* Set photon's path ID from ray serial num (supplemented by photon.proc below) */ photon.aux.pathID = ray -> rno; } - + /* Set normal */ for (i = 0; i <= 2; i++) - photon.norm [i] = 127.0 * (isVolumePmap(pmap) + photon.norm [i] = 127.0 * (isVolumePmap(pmap) #ifdef PMAP_PHOTONFLOW || isLightFlowPmap(pmap) #endif ? ray -> rdir [i] : ray -> ron [i] ); if (!pmap -> heapBuf) { /* Lazily allocate photon heap buffa */ #if NIX /* Randomise buffa size to temporally decorellate flushes in * multiprocessing mode */ srandom(randSeed + getpid()); pmap -> heapBufSize = PMAP_HEAPBUFSIZE * 0.5 * (1 + frandom()); #else /* Randomisation disabled for single processes on WIN; also useful * for reproducability during debugging */ pmap -> heapBufSize = PMAP_HEAPBUFSIZE; #endif -#ifdef PMAP_CONTRIB +#ifdef PMAP_CONTRIB if (contribs && pmap -> preCompContrib) /* Adapt buffa size to accommodate precomputed contributions */ - pmap -> heapBufSize *= (float)sizeof(Photon) / + pmap -> heapBufSize *= (float)sizeof(Photon) / ((PreComputedContrib*)pmap -> preCompContrib) -> contribSize; #endif if (!(pmap -> heapBuf = calloc(pmap -> heapBufSize, sizeof(Photon)))) - error(SYSTEM, + error(SYSTEM, "out of memory allocating heap buffer in newPhoton()" ); pmap -> heapBufLen = 0; } /* Photon initialised; now append to heap buffa */ memcpy(pmap -> heapBuf + pmap -> heapBufLen, &photon, sizeof(Photon)); #ifdef PMAP_CONTRIB if (contribs && pmap -> preCompContrib) { - const unsigned long contribSize = + const unsigned long contribSize = ((PreComputedContrib*)pmap -> preCompContrib) -> contribSize; - + if (!pmap -> contribHeapBuf) { /* Lazily allocate contribution heap buffa of same size as photon * heap buffa */ pmap -> contribHeapBuf = calloc(pmap -> heapBufSize, contribSize); if (!pmap -> contribHeapBuf) error(SYSTEM, "out of memory allocating precomputed " "contribution heap buffer in newPhoton()" ); } - + /* Append precomputed contribs to heap buffa */ memcpy(pmap -> contribHeapBuf + pmap -> heapBufLen * contribSize, contribs, contribSize ); } #endif if (++pmap -> heapBufLen >= pmap -> heapBufSize) - /* Photon (and contrib, if applicable) heap buffa(s) full, - flush to heap file(s) */ + /* Photon (and contrib, if applicable) heap buffa(s) full, + flush to heap file(s) */ flushPhotonHeap(pmap); pmap -> numPhotons++; /* Print photon attributes */ if (printPhoton) /* Non-const kludge */ printPhoton((RAY*)ray, &photon, pmap); return 0; } -void buildPhotonMap (PhotonMap *pmap, double *photonFlux, +void buildPhotonMap (PhotonMap *pmap, double *photonFlux, PhotonContribSourceIdx *contribSrcOfs, unsigned nproc ) { unsigned long n, nCheck = 0; unsigned i; Photon *p; COLOR flux; char nuHeapFname [sizeof(PMAP_TMPFNAME)]; FILE *nuHeap; /* Need double here to reduce summation errors */ double avgFlux [3] = {0, 0, 0}, CoG [3] = {0, 0, 0}, CoGdist = 0, avgPathLen = 0; FVECT d; if (!pmap) error(INTERNAL, "undefined photon map in buildPhotonMap()"); /* Get number of photons from heapfile size */ if (fseek(pmap -> heap, 0, SEEK_END) < 0) error(SYSTEM, "failed seek to end of photon heap in buildPhotonMap()"); pmap -> numPhotons = ftell(pmap -> heap) / sizeof(Photon); if (!pmap -> numPhotons) error(INTERNAL, "empty photon map in buildPhotonMap()"); if (!pmap -> heap) error(INTERNAL, "no heap in buildPhotonMap()"); #ifdef DEBUG_PMAP eputs("Checking photon heap consistency...\n"); checkPhotonHeap(pmap -> heap); sprintf(errmsg, "Heap contains %ld photons\n", pmap -> numPhotons); eputs(errmsg); #endif /* Allocate heap buffa */ if (!pmap -> heapBuf) { pmap -> heapBufSize = PMAP_HEAPBUFSIZE; pmap -> heapBuf = calloc(pmap -> heapBufSize, sizeof(Photon)); if (!pmap -> heapBuf) - error(SYSTEM, + error(SYSTEM, "failed to allocate postprocessed photon heap in buildPhotonMap()" ); } /* We REALLY don't need yet another @%&*! heap just to hold the scaled * photons, but can't think of a quicker fix... */ mktemp(strcpy(nuHeapFname, PMAP_TMPFNAME)); if (!(nuHeap = fopen(nuHeapFname, "w+b"))) - error(SYSTEM, + error(SYSTEM, "failed to open postprocessed photon heap in buildPhotonMap()" ); rewind(pmap -> heap); -#ifdef DEBUG_PMAP +#ifdef DEBUG_PMAP eputs("Postprocessing photons...\n"); #endif while (!feof(pmap -> heap)) { -#ifdef DEBUG_PMAP - printf("Reading %lu at %lu... ", +#ifdef DEBUG_PMAP + printf("Reading %lu at %lu... ", pmap -> heapBufSize, ftell(pmap->heap) ); #endif - pmap -> heapBufLen = fread(pmap -> heapBuf, sizeof(Photon), + pmap -> heapBufLen = fread(pmap -> heapBuf, sizeof(Photon), pmap -> heapBufSize, pmap -> heap ); #ifdef DEBUG_PMAP printf("Got %lu\n", pmap -> heapBufLen); #endif if (ferror(pmap -> heap)) error(SYSTEM, "failed to read photon heap in buildPhotonMap()"); for (n = pmap -> heapBufLen, p = pmap -> heapBuf; n; n--, p++) { /* Update min and max pos and set photon flux */ for (i = 0; i <= 2; i++) { - if (p -> pos [i] < pmap -> minPos [i]) + if (p -> pos [i] < pmap -> minPos [i]) pmap -> minPos [i] = p -> pos [i]; - else if (p -> pos [i] > pmap -> maxPos [i]) + else if (p -> pos [i] > pmap -> maxPos [i]) pmap -> maxPos [i] = p -> pos [i]; /* Update centre of gravity with photon position */ CoG [i] += p -> pos [i]; } if (isContribPmap(pmap) && contribSrcOfs) - /* Linearise contrib photon origin index from subprocess index + /* Linearise contrib photon origin index from subprocess index * using the per-subprocess offsets in contribSrcOfs */ p -> aux.contribSrc += contribSrcOfs [p -> proc]; else if (isTransientPmap(pmap) #ifdef PMAP_PHOTONFLOW || isTransLightFlowPmap(pmap) #endif ) { /* Update min and max path length */ if (p -> aux.pathLen < pmap -> minPathLen) pmap -> minPathLen = p -> aux.pathLen; else if (p -> aux.pathLen > pmap -> maxPathLen) pmap -> maxPathLen = p -> aux.pathLen; /* Update avg photon path length */ avgPathLen += p -> aux.pathLen; } /* Scale photon's flux (hitherto normalised to 1 over RGB); in * case of a contrib photon map, this is done per light source, * and photonFlux is assumed to be an array */ getPhotonFlux(p, flux); if (photonFlux) { - scalecolor(flux, + scalecolor(flux, #if 1 photonFlux [isContribPmap(pmap) ? photonSrcIdx(pmap, p) : 0] #else photonFlux [0] #endif ); setPhotonFlux(p, flux); } /* Update average photon flux; need a double here */ addcolor(avgFlux, flux); } /* Write modified photons to new heap */ fwrite(pmap -> heapBuf, sizeof(Photon), pmap -> heapBufLen, nuHeap); if (ferror(nuHeap)) { - error(SYSTEM, + error(SYSTEM, "failed postprocessing photon flux in buildPhotonMap()" ); /* Clean up temp files */ fclose(pmap -> heap); fclose(nuHeap); unlink(pmap -> heapFname); unlink(nuHeapFname); } nCheck += pmap -> heapBufLen; } - + /* !!! Flush pending writes to new heap !!! */ fflush(nuHeap); #ifdef DEBUG_PMAP if (nCheck < pmap -> numPhotons) error(INTERNAL, "truncated photon heap in buildPhotonMap"); #endif - + /* Finalise average photon flux */ scalecolor(avgFlux, 1.0 / pmap -> numPhotons); copycolor(pmap -> photonFlux, avgFlux); /* Average photon positions to get centre of gravity */ for (i = 0; i < 3; i++) pmap -> CoG [i] = CoG [i] /= pmap -> numPhotons; /* Average photon path lengths */ pmap -> avgPathLen = avgPathLen /= pmap -> numPhotons; rewind(pmap -> heap); /* Compute average photon distance to centre of gravity */ while (!feof(pmap -> heap)) { - pmap -> heapBufLen = fread(pmap -> heapBuf, sizeof(Photon), + pmap -> heapBufLen = fread(pmap -> heapBuf, sizeof(Photon), pmap -> heapBufSize, pmap -> heap ); for (n = pmap -> heapBufLen, p = pmap -> heapBuf; n; n--, p++) { VSUB(d, p -> pos, CoG); CoGdist += DOT(d, d); - + if (isTransientPmap(pmap) #ifdef PMAP_PHOTONFLOW || isTransLightFlowPmap(pmap) #endif ) { - /* Add temporal distance (in units of photon path length) + /* Add temporal distance (in units of photon path length) for transient photons */ d [0] = p -> aux.pathLen - avgPathLen; CoGdist += d [0] * d [0]; } } } pmap -> CoGdist = CoGdist /= pmap -> numPhotons; /* Swap heaps, discarding unscaled photons */ fclose(pmap -> heap); unlink(pmap -> heapFname); pmap -> heap = nuHeap; strcpy(pmap -> heapFname, nuHeapFname); - + #ifdef PMAP_OOC OOC_BuildPhotonMap(pmap, nproc); #else if (isTransientPmap(pmap) #ifdef PMAP_PHOTONFLOW || isTransLightFlowPmap(pmap) #endif ) kdT_TransBuildPhotonMap(pmap); else kdT_BuildPhotonMap(pmap); #endif /* Trash heap and its buffa */ free(pmap -> heapBuf); fclose(pmap -> heap); unlink(pmap -> heapFname); pmap -> heap = NULL; pmap -> heapBuf = NULL; } /* PHOTON MAP SEARCH ROUTINES ------------------------------------------ */ /* Dynamic max photon search radius increase and reduction factors */ #define PMAP_MAXDIST_INC 4 #define PMAP_MAXDIST_DEC 0.9 /* Num successful lookups before reducing in max search radius */ #define PMAP_MAXDIST_CNT 1000 /* Threshold below which we assume increasing max radius won't help */ #define PMAP_SHORT_LOOKUP_THRESH 1 /* Coefficient for adaptive maximum search radius */ #define PMAP_MAXDIST_COEFF 100 void findPhotons (PhotonMap* pmap, const RAY* ray) { int redo = 0; const RREAL *norm = ray -> ron, *pos = ray -> rop; if (!pmap -> squeue.len) { /* Lazy init priority queue */ #ifdef PMAP_OOC OOC_InitFindPhotons(pmap); #else kdT_InitFindPhotons(pmap); #endif pmap -> minGathered = pmap -> maxGather; pmap -> maxGathered = pmap -> minGather; pmap -> totalGathered = 0; pmap -> numLookups = pmap -> numShortLookups = 0; pmap -> shortLookupPct = 0; pmap -> minError = FHUGE; pmap -> maxError = -FHUGE; pmap -> rmsError = 0; /* SQUARED max search radius limit is based on avg photon distance to * centre of gravity, unless fixed by user (maxDistFix > 0) */ - pmap -> maxDist0 = pmap -> maxDist2Limit = maxDistFix > 0 + pmap -> maxDist0 = pmap -> maxDist2Limit = maxDistFix > 0 ? maxDistFix * maxDistFix - : PMAP_MAXDIST_COEFF * pmap -> squeue.len * + : PMAP_MAXDIST_COEFF * pmap -> squeue.len * pmap -> CoGdist / pmap -> numPhotons; } do { pmap -> squeue.tail = 0; pmap -> maxDist2 = pmap -> maxDist0; - if (isVolumePmap(pmap) -#ifdef PMAP_PHOTONFLOW + if (isVolumePmap(pmap) +#ifdef PMAP_PHOTONFLOW || isLightFlowPmap(pmap) && sphericalIrrad #endif ) { /* Volume photons are not associated with a normal or intersection point; null the normal and set origin as lookup point. */ - /* If a lightflow photon map is used and sphericalIrrad = 1, + /* If a lightflow photon map is used and sphericalIrrad = 1, the mean spherical irradiance is evaluated, so normals are irrelevant and not filtered. */ norm = NULL; pos = ray -> rorg; } - /* XXX/NOTE: status returned by XXX_FindPhotons() is currently - ignored; if no photons are found, an empty queue is returned + /* XXX/NOTE: status returned by XXX_FindPhotons() is currently + ignored; if no photons are found, an empty queue is returned under the assumption all photons are too distant to contribute significant flux. */ #ifdef PMAP_OOC OOC_FindPhotons(pmap, pos, norm); #else if (isTransientPmap(pmap) #ifdef PMAP_PHOTONFLOW || isTransLightFlowPmap(pmap) #endif ) kdT_TransFindPhotons(pmap, pos, norm); else kdT_FindPhotons(pmap, pos, norm); #endif #ifdef PMAP_LOOKUP_INFO - fprintf(stderr, "%d/%d %s photons found within radius %.3f " - "at (%.2f,%.2f,%.2f) on %s\n", pmap -> squeue.tail, + fprintf(stderr, "%d/%d %s photons found within radius %.3f " + "at (%.2f,%.2f,%.2f) on %s\n", pmap -> squeue.tail, pmap -> squeue.len, pmapName [pmap -> type], sqrt(pmap -> maxDist2), - ray -> rop [0], ray -> rop [1], ray -> rop [2], + ray -> rop [0], ray -> rop [1], ray -> rop [2], ray -> ro ? ray -> ro -> oname : "" ); #endif if (pmap -> squeue.tail < pmap -> squeue.len * pmap->gatherTolerance) { /* Short lookup; too few photons found */ if (pmap -> squeue.tail > PMAP_SHORT_LOOKUP_THRESH) { /* Ignore short lookups which return fewer than * PMAP_SHORT_LOOKUP_THRESH photons under the assumption there * really are no photons in the vicinity, and increasing the max * search radius therefore won't help */ #ifdef PMAP_LOOKUP_WARN - sprintf(errmsg, + sprintf(errmsg, "%d/%d %s photons found at (%.2f,%.2f,%.2f) on %s", pmap -> squeue.tail, pmap->squeue.len, pmapName [pmap->type], - ray -> rop [0], ray -> rop [1], ray -> rop [2], + ray -> rop [0], ray -> rop [1], ray -> rop [2], ray -> ro ? ray -> ro -> oname : "" ); error(WARNING, errmsg); #endif /* Bail out after warning if maxDist is fixed */ if (maxDistFix > 0) return; if (pmap -> maxDist0 < pmap -> maxDist2Limit) { /* Increase max search radius if below limit & redo search */ pmap -> maxDist0 *= PMAP_MAXDIST_INC; #ifdef PMAP_LOOKUP_REDO redo = 1; #endif #ifdef PMAP_LOOKUP_WARN - sprintf(errmsg, redo + sprintf(errmsg, redo ? "restarting photon lookup with max radius %.1e" : "max photon lookup radius adjusted to %.1e", pmap -> maxDist0 ); error(WARNING, errmsg); #endif } #ifdef PMAP_LOOKUP_REDO else { sprintf(errmsg, "max photon lookup radius clamped to %.1e", pmap -> maxDist0 ); error(WARNING, errmsg); } #endif } /* Reset successful lookup counter */ pmap -> numLookups = 0; } else { /* Skip search radius reduction if maxDist is fixed */ if (maxDistFix > 0) return; /* Increment successful lookup counter and reduce max search radius if * wraparound */ pmap -> numLookups = (pmap -> numLookups + 1) % PMAP_MAXDIST_CNT; if (!pmap -> numLookups) pmap -> maxDist0 *= PMAP_MAXDIST_DEC; - + redo = 0; } - + } while (redo); } -Photon *find1Photon (PhotonMap *pmap, const RAY* ray, Photon *photon, +Photon *find1Photon (PhotonMap *pmap, const RAY* ray, Photon *photon, void *photonIdx ) { /* Init (squared) search radius to avg photon dist to centre of gravity */ float maxDist2_0 = pmap -> CoGdist; - int found = 0; + int found = 0; #ifdef PMAP_LOOKUP_REDO #define REDO 1 #else #define REDO 0 #endif do { pmap -> maxDist2 = maxDist2_0; #ifdef PMAP_OOC found = OOC_Find1Photon(pmap, ray -> rop, ray -> ron, photon, (OOC_DataIdx*)photonIdx ); #else found = kdT_Find1Photon(pmap, ray -> rop, ray -> ron, photon); #endif if (found < 0) { /* Expand search radius to retry */ maxDist2_0 *= 2; #ifdef PMAP_LOOKUP_WARN sprintf(errmsg, "failed 1-NN photon lookup" #ifdef PMAP_LOOKUP_REDO ", retrying with search radius %.2f", maxDist2_0 #endif ); error(WARNING, errmsg); #endif } } while (REDO && found < 0); /* Return photon buffer containing valid photon, else NULL */ return found < 0 ? NULL : photon; } void getPhoton (PhotonMap *pmap, PhotonIdx idx, Photon *photon) { #ifdef PMAP_OOC if (OOC_GetPhoton(pmap, idx, photon)) #else if (kdT_GetPhoton(pmap, idx, photon)) #endif error(INTERNAL, "failed photon lookup"); } Photon *getNearestPhoton (const PhotonSearchQueue *squeue, PhotonIdx idx) { #ifdef PMAP_OOC return OOC_GetNearestPhoton(squeue, idx); #else return kdT_GetNearestPhoton(squeue, idx); #endif } PhotonIdx firstPhoton (const PhotonMap *pmap) { #ifdef PMAP_OOC return OOC_FirstPhoton(pmap); #else return kdT_FirstPhoton(pmap); #endif } /* PHOTON MAP CLEANUP ROUTINES ------------------------------------------ */ void deletePhotons (PhotonMap* pmap) { #ifdef PMAP_OOC OOC_Delete(&pmap -> store); #else kdT_Delete(&pmap -> store); #endif free(pmap -> squeue.node); free(pmap -> biasCompHist); pmap -> numPhotons = pmap -> minGather = pmap -> maxGather = pmap -> squeue.len = pmap -> squeue.tail = 0; #ifdef PMAP_PATHFILT if (pmap -> pathLUT) { lu_done(pmap -> pathLUT); free(pmap -> pathLUT); } if (pmap -> pathLUTKeys) { unsigned k; for (k = 0; k < pmap->squeue.len + 1; free(pmap->pathLUTKeys [k++])); free(pmap -> pathLUTKeys); } #endif /* Contrib stuff */ if (isContribPmap(pmap) && pmap -> preCompContribTab) { /* Parent contrib photon map; clean up child photon maps containing precomputed contribs via lu_done() */ lu_done(pmap -> preCompContribTab); free(pmap -> preCompContribTab); pmap -> preCompContribTab = NULL; - + if (pmap -> contribHeapBuf) { free(pmap-> contribHeapBuf); pmap -> contribHeapBuf = NULL; } if (pmap -> contribHeap) { fclose(pmap -> contribHeap); pmap -> contribHeap = NULL; } } } - diff --git a/pmapmat.c b/pmapmat.c index c0c04eb..f7e47e5 100644 --- a/pmapmat.c +++ b/pmapmat.c @@ -1,2075 +1,2087 @@ #ifndef lint static const char RCSid[] = "$Id: pmapmat.c,v 2.24 2021/02/22 13:27:49 rschregle Exp $"; #endif -/* +/* ====================================================================== - Photon map support routines for scattering by materials. + Photon map support routines for scattering by materials. Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) (c) Fraunhofer Institute for Solar Energy Systems, - supported by the German Research Foundation + supported by the German Research Foundation (DFG LU-204/10-2, "Fassadenintegrierte Regelsysteme FARESYS") (c) Lucerne University of Applied Sciences and Arts, - supported by the Swiss National Science Foundation + supported by the Swiss National Science Foundation (SNSF #147053, "Daylight Redirecting Components") ====================================================================== - + */ #include "pmapmat.h" #include "pmapdata.h" #include "pmaprand.h" +#include "pmapsrc.h" +#include "otspecial.h" #include "otypes.h" #include "data.h" #include "func.h" #include "bsdf.h" #include /* Stuff ripped off from material modules */ #define MAXITER 10 #define SP_REFL 01 -#define SP_TRAN 02 +#define SP_TRAN 02 #define SP_PURE 04 #define SP_FLAT 010 #define SP_BADU 040 #define MLAMBDA 500 #define RINDEX 1.52 #define FRESNE(ci) (exp(-5.85*(ci)) - 0.00287989916) /* Get antimatter sensor orientation flags. HACK: these replace the last character in the object name! */ #define PMAP_GETSENSFLAGS(obj) ((obj) -> oname [strlen((obj) -> oname) - 1]) typedef struct { OBJREC *mp; RAY *rp; short specfl; COLOR mcolor, scolor; FVECT vrefl, prdir, pnorm; double alpha2, rdiff, rspec, trans, tdiff, tspec, pdot; } NORMDAT; typedef struct { OBJREC *mp; RAY *rp; short specfl; COLOR mcolor, scolor; FVECT vrefl, prdir, u, v, pnorm; double u_alpha, v_alpha, rdiff, rspec, trans, tdiff, tspec, pdot; } ANISODAT; typedef struct { OBJREC *mp; RAY *pr; DATARRAY *dp; COLOR mcolor; COLOR rdiff; COLOR tdiff; double rspec; double trans; double tspec; FVECT pnorm; double pdot; } BRDFDAT; typedef struct { OBJREC *mp; RAY *pr; FVECT pnorm; FVECT vray; double sr_vpsa [2]; RREAL toloc [3][3]; RREAL fromloc [3][3]; double thick; SDData *sd; COLOR runsamp; COLOR rdiff; COLOR tunsamp; COLOR tdiff; } BSDFDAT; extern const SDCDst SDemptyCD; /* Per-material scattering function dispatch table; return value is usually * zero, indicating photon termination */ int (*photonScatter [NUMOTYPE]) (OBJREC*, RAY*); /* List of antimatter sensor modifier names and associated object set */ char *photonSensorList [MAXSET + 1] = {NULL}; static OBJECT photonSensorSet [MAXSET + 1] = {0}; /* ================ General support routines ================ */ -void photonRay (const RAY *rayIn, RAY *rayOut, int rayOutType, +void photonRay (const RAY *rayIn, RAY *rayOut, int rayOutType, COLOR fluxAtten ) /* Spawn a new photon ray from a previous one; this is effectively a - * customised rayorigin(). + * customised rayorigin(). * A SPECULAR rayOutType flags this photon as _caustic_ for subsequent hits. * It is preserved for transferred rays (of type PMAP_XFER). * fluxAtten specifies the RGB attenuation of the photon flux effected by * the scattering material. The outgoing flux is then normalised to maintain * a uniform average of 1 over RGB. If fluxAtten == NULL, the flux remains * unchanged for the outgoing photon. fluxAtten is ignored for transferred * rays. * The ray direction is preserved for transferred rays, and undefined for * scattered rays and must be subsequently set by the caller. */ { rayorigin(rayOut, rayOutType, rayIn, NULL); - + if (rayIn) { /* Transfer flux */ copycolor(rayOut -> rcol, rayIn -> rcol); - + /* Copy caustic flag & direction for transferred rays */ if (rayOutType == PMAP_XFER) { /* rayOut -> rtype |= rayIn -> rtype & SPECULAR; */ rayOut -> rtype |= rayIn -> rtype; VCOPY(rayOut -> rdir, rayIn -> rdir); } else if (fluxAtten) { /* Attenuate and normalise flux for scattered rays */ multcolor(rayOut -> rcol, fluxAtten); colorNorm(rayOut -> rcol); } /* Propagate index of emitting light source */ rayOut -> rsrc = rayIn -> rsrc; - + /* Update maximum photon path distance; this is passed as aft plane distance to localhit() in tracePhoton(); it has no effect if negative, but limits the ray distance if positive. */ rayOut -> rmax = rayIn -> rmax - rayIn -> rot; } } static void addPhotons (const RAY *r) /* Insert photon hits, where applicable */ { if (!r -> rlvl) /* Add direct photon at primary hitpoint */ newPhoton(directPmap, r, NULL); - else { - /* Add global or precomputed photon at indirect hitpoint */ + + if (r->rlvl || weaksrcmat(source[(r)->rsrc].so->omod) ) + /* Add global or precomputed photon at indirect hitpoint, + primary global photon if from glow or illum */ newPhoton(preCompPmap ? preCompPmap : globalPmap, r, NULL); + if (PMAP_CAUSTICRAY(r)) /* Store caustic photon if specular flag set */ - if (PMAP_CAUSTICRAY(r)) - newPhoton(causticPmap, r, NULL); - + newPhoton(causticPmap, r, NULL); + + if (r->rlvl || weaksrcmat(source[(r)->rsrc].so->omod) ) /* Store in contribution photon map */ newPhoton(contribPmap, r, NULL); - } - + /* Add transient photon at all hitpoints */ newPhoton(transientPmap, r, NULL); } void getPhotonSensors (char **sensorList) /* Find antimatter geometry declared as photon sensors */ { OBJECT i; OBJREC *obj; char **lp; - + /* Init sensor set */ photonSensorSet [0] = 0; - + if (!sensorList [0]) return; - + for (i = 0; i < nobjects; i++) { obj = objptr(i); - + /* Insert object in sensor set if it's in the specified sensor list * and of type antimatter. Note the string comparison omits the last * character, since this contains the orientation flags. */ for (lp = sensorList; *lp; lp++) { if (!strncmp(obj -> oname, *lp, strlen(*lp) - 1)) { if (obj -> otype != MAT_CLIP) { sprintf(errmsg, "photon sensor modifier %s is not antimatter", obj -> oname ); error(USER, errmsg); } - + if (photonSensorSet [0] >= AMBLLEN) error(USER, "too many photon sensor modifiers"); - + /* Set orientation flags in object */ PMAP_GETSENSFLAGS(obj) = (*lp) [strlen(*lp) - 1]; insertelem(photonSensorSet, i); } } } - + if (!photonSensorSet [0]) error(USER, "no photon sensors found"); } -/* ================ Material specific scattering routines ================ */ +/* ================ Material specific scattering routines ================ */ static int isoSpecPhotonScatter (NORMDAT *nd, RAY *rayOut) -/* Generate direction for isotropically specularly reflected +/* Generate direction for isotropically specularly reflected or transmitted ray. Returns 1 if successful. */ { FVECT u, v, h; RAY *rayIn = nd -> rp; double d, d2, sinp, cosp; int niter, i = 0; - - /* Set up sample coordinates */ + + /* Set up sample coordinates */ getperpendicular(u, nd -> pnorm, 1); fcross(v, nd -> pnorm, u); - + if (nd -> specfl & SP_REFL) { /* Specular reflection; make MAXITER attempts at getting a ray */ - + for (niter = 0; niter < MAXITER; niter++) { d = 2 * PI * pmapRandom(scatterState); cosp = cos(d); sinp = sin(d); d2 = pmapRandom(scatterState); d = d2 <= FTINY ? 1 : sqrt(nd -> alpha2 * -log(d2)); - + for (i = 0; i < 3; i++) h [i] = nd -> pnorm [i] + d * (cosp * u [i] + sinp * v [i]); - + d = -2 * DOT(h, rayIn -> rdir) / (1 + d * d); VSUM(rayOut -> rdir, rayIn -> rdir, h, d); - - if (DOT(rayOut -> rdir, rayIn -> ron) > FTINY) + + if (DOT(rayOut -> rdir, rayIn -> ron) > FTINY) return 1; } - + return 0; } - + else { /* Specular transmission; make MAXITER attempts at getting a ray */ - + for (niter = 0; niter < MAXITER; niter++) { d = 2 * PI * pmapRandom(scatterState); cosp = cos(d); sinp = sin(d); d2 = pmapRandom(scatterState); d = d2 <= FTINY ? 1 : sqrt(-log(d2) * nd -> alpha2); for (i = 0; i < 3; i++) - rayOut -> rdir [i] = nd -> prdir [i] + + rayOut -> rdir [i] = nd -> prdir [i] + d * (cosp * u [i] + sinp * v [i]); if (DOT(rayOut -> rdir, rayIn -> ron) < -FTINY) { normalize(rayOut -> rdir); return 1; } } - + return 0; } } static void diffPhotonScatter (FVECT normal, RAY* rayOut) /* Generate cosine-weighted direction for diffuse ray */ { const RREAL cosThetaSqr = pmapRandom(scatterState), cosTheta = sqrt(cosThetaSqr), sinTheta = sqrt(1 - cosThetaSqr), phi = 2 * PI * pmapRandom(scatterState), du = cos(phi) * sinTheta, dv = sin(phi) * sinTheta; FVECT u, v; int i = 0; /* Set up sample coordinates */ getperpendicular(u, normal, 1); fcross(v, normal, u); - + /* Convert theta & phi to cartesian */ for (i = 0; i < 3; i++) rayOut -> rdir [i] = du * u [i] + dv * v [i] + cosTheta * normal [i]; normalize(rayOut -> rdir); } static int normalPhotonScatter (OBJREC *mat, RAY *rayIn) /* Generate new photon ray for isotropic material and recurse */ { NORMDAT nd; int i, hastexture; float xi, albedo, prdiff, ptdiff, prspec, ptspec; double d, fresnel; RAY rayOut; if (mat -> oargs.nfargs != (mat -> otype == MAT_TRANS ? 7 : 5)) objerror(mat, USER, "bad number of arguments"); - + /* Check for back side; reorient if back is visible */ if (rayIn -> rod < 0) - if (!backvis && mat -> otype != MAT_TRANS) + if (!backvis && mat -> otype != MAT_TRANS) return 0; else { /* Get modifiers */ raytexture(rayIn, mat -> omod); flipsurface(rayIn); } else raytexture(rayIn, mat -> omod); - + nd.mp = mat; nd.rp = rayIn; - + /* Get material color */ copycolor(nd.mcolor, mat -> oargs.farg); - + /* Get roughness */ nd.specfl = 0; nd.alpha2 = mat -> oargs.farg [4]; - - if ((nd.alpha2 *= nd.alpha2) <= FTINY) + + if ((nd.alpha2 *= nd.alpha2) <= FTINY) nd.specfl |= SP_PURE; - - if (rayIn -> ro != NULL && isflat(rayIn -> ro -> otype)) - nd.specfl |= SP_FLAT; - + + if (rayIn -> ro != NULL && isflat(rayIn -> ro -> otype)) + nd.specfl |= SP_FLAT; + /* Perturb normal */ if ((hastexture = (DOT(rayIn -> pert, rayIn -> pert) > sqr(FTINY)) )) nd.pdot = raynormal(nd.pnorm, rayIn); else { VCOPY(nd.pnorm, rayIn -> ron); nd.pdot = rayIn -> rod; } - + nd.pdot = max(nd.pdot, .001); - + /* Modify material color */ multcolor(nd.mcolor, rayIn -> pcol); nd.rspec = mat -> oargs.farg [3]; - + /* Approximate Fresnel term */ if (nd.specfl & SP_PURE && nd.rspec > FTINY) { fresnel = FRESNE(rayIn -> rod); nd.rspec += fresnel * (1 - nd.rspec); - } + } else fresnel = 0; - + /* Transmission params */ if (mat -> otype == MAT_TRANS) { nd.trans = mat -> oargs.farg [5] * (1 - nd.rspec); nd.tspec = nd.trans * mat -> oargs.farg [6]; - nd.tdiff = nd.trans - nd.tspec; - } + nd.tdiff = nd.trans - nd.tspec; + } else nd.tdiff = nd.tspec = nd.trans = 0; - + /* Specular reflection params */ if (nd.rspec > FTINY) { /* Specular color */ - if (mat -> otype != MAT_METAL) + if (mat -> otype != MAT_METAL) setcolor(nd.scolor, nd.rspec, nd.rspec, nd.rspec); else if (fresnel > FTINY) { d = nd.rspec * (1 - fresnel); - for (i = 0; i < 3; i++) + for (i = 0; i < 3; i++) nd.scolor [i] = fresnel + nd.mcolor [i] * d; } else { copycolor(nd.scolor, nd.mcolor); scalecolor(nd.scolor, nd.rspec); } } else setcolor(nd.scolor, 0, 0, 0); - + /* Diffuse reflection params */ nd.rdiff = 1 - nd.trans - nd.rspec; - + /* Set up probabilities */ prdiff = ptdiff = ptspec = colorAvg(nd.mcolor); prdiff *= nd.rdiff; ptdiff *= nd.tdiff; prspec = colorAvg(nd.scolor); ptspec *= nd.tspec; albedo = prdiff + ptdiff + prspec + ptspec; - + /* Insert direct and indirect photon hits if diffuse component */ if (prdiff > FTINY || ptdiff > FTINY) addPhotons(rayIn); - + xi = pmapRandom(rouletteState); - if (xi > albedo) + if (xi > albedo) /* Absorbed */ return 0; if (xi > (albedo -= prspec)) { /* Specular reflection */ + int rotype = 0; + rotype|=REFLECTED; nd.specfl |= SP_REFL; - + if (nd.specfl & SP_PURE) { /* Perfect specular reflection */ for (i = 0; i < 3; i++) { /* Reflected ray */ nd.vrefl [i] = rayIn -> rdir [i] + 2 * nd.pdot * nd.pnorm [i]; } /* Penetration? */ if (hastexture && DOT(nd.vrefl, rayIn -> ron) <= FTINY) for (i = 0; i < 3; i++) { /* Safety measure */ - nd.vrefl [i] = rayIn -> rdir [i] + + nd.vrefl [i] = rayIn -> rdir [i] + 2 * rayIn -> rod * rayIn -> ron [i]; } VCOPY(rayOut.rdir, nd.vrefl); } - - else if (!isoSpecPhotonScatter(&nd, &rayOut)) - return 0; - - photonRay(rayIn, &rayOut, PMAP_SPECREFL, nd.scolor); + + else { + if (!isoSpecPhotonScatter(&nd, &rayOut)) + return 0; + rotype|=SPECULAR; + } + photonRay(rayIn, &rayOut, rotype, nd.scolor); } - + else if (xi > (albedo -= ptspec)) { /* Specular transmission */ + int rotype = 0; + rotype |= TRANS; nd.specfl |= SP_TRAN; - + if (hastexture) { /* Perturb */ for (i = 0; i < 3; i++) nd.prdir [i] = rayIn -> rdir [i] - rayIn -> pert [i]; - + if (DOT(nd.prdir, rayIn -> ron) < -FTINY) normalize(nd.prdir); else VCOPY(nd.prdir, rayIn -> rdir); } else VCOPY(nd.prdir, rayIn -> rdir); - + if ((nd.specfl & (SP_TRAN | SP_PURE)) == (SP_TRAN | SP_PURE)) /* Perfect specular transmission */ VCOPY(rayOut.rdir, nd.prdir); - else if (!isoSpecPhotonScatter(&nd, &rayOut)) - return 0; - - photonRay(rayIn, &rayOut, PMAP_SPECTRANS, nd.mcolor); + else { + if (!isoSpecPhotonScatter(&nd, &rayOut)) + return 0; + rotype |= SPECULAR; + } + photonRay(rayIn, &rayOut, rotype, nd.mcolor); } - + else if (xi > (albedo -= prdiff)) { /* Diffuse reflection */ photonRay(rayIn, &rayOut, PMAP_DIFFREFL, nd.mcolor); diffPhotonScatter(hastexture ? nd.pnorm : rayIn -> ron, &rayOut); } - + else { /* Diffuse transmission */ flipsurface(rayIn); photonRay(rayIn, &rayOut, PMAP_DIFFTRANS, nd.mcolor); if (hastexture) { FVECT bnorm; bnorm [0] = -nd.pnorm [0]; bnorm [1] = -nd.pnorm [1]; bnorm [2] = -nd.pnorm [2]; diffPhotonScatter(bnorm, &rayOut); - } + } else diffPhotonScatter(rayIn -> ron, &rayOut); } - + tracePhoton(&rayOut); return 0; } static void getacoords (ANISODAT *nd) /* Set up coordinate system for anisotropic sampling; cloned from aniso.c */ { MFUNC *mf; int i; mf = getfunc(nd -> mp, 3, 0x7, 1); setfunc(nd -> mp, nd -> rp); errno = 0; for (i = 0; i < 3; i++) nd -> u [i] = evalue(mf -> ep [i]); - + if (errno == EDOM || errno == ERANGE) nd -> u [0] = nd -> u [1] = nd -> u [2] = 0.0; - + if (mf -> fxp != &unitxf) multv3(nd -> u, nd -> u, mf -> fxp -> xfm); fcross(nd -> v, nd -> pnorm, nd -> u); - + if (normalize(nd -> v) == 0.0) { if (fabs(nd -> u_alpha - nd -> v_alpha) > 0.001) objerror(nd -> mp, WARNING, "illegal orientation vector"); getperpendicular(nd -> u, nd -> pnorm, 1); fcross(nd -> v, nd -> pnorm, nd -> u); - nd -> u_alpha = nd -> v_alpha = + nd -> u_alpha = nd -> v_alpha = sqrt(0.5 * (sqr(nd -> u_alpha) + sqr(nd -> v_alpha))); - } + } else fcross(nd -> u, nd -> v, nd -> pnorm); } static int anisoSpecPhotonScatter (ANISODAT *nd, RAY *rayOut) -/* Generate direction for anisotropically specularly reflected +/* Generate direction for anisotropically specularly reflected or transmitted ray. Returns 1 if successful. */ { FVECT h; double d, d2, sinp, cosp; int niter, i; RAY *rayIn = nd -> rp; - if (rayIn -> ro != NULL && isflat(rayIn -> ro -> otype)) + if (rayIn -> ro != NULL && isflat(rayIn -> ro -> otype)) nd -> specfl |= SP_FLAT; - + /* set up coordinates */ getacoords(nd); - + if (rayOut -> rtype & TRANS) { /* Specular transmission */ - if (DOT(rayIn -> pert, rayIn -> pert) <= sqr(FTINY)) + if (DOT(rayIn -> pert, rayIn -> pert) <= sqr(FTINY)) VCOPY(nd -> prdir, rayIn -> rdir); else { /* perturb */ - for (i = 0; i < 3; i++) + for (i = 0; i < 3; i++) nd -> prdir [i] = rayIn -> rdir [i] - rayIn -> pert [i]; - - if (DOT(nd -> prdir, rayIn -> ron) < -FTINY) + + if (DOT(nd -> prdir, rayIn -> ron) < -FTINY) normalize(nd -> prdir); else VCOPY(nd -> prdir, rayIn -> rdir); } - + /* Make MAXITER attempts at getting a ray */ for (niter = 0; niter < MAXITER; niter++) { d = 2 * PI * pmapRandom(scatterState); cosp = cos(d) * nd -> u_alpha; sinp = sin(d) * nd -> v_alpha; d = sqrt(sqr(cosp) + sqr(sinp)); cosp /= d; sinp /= d; d2 = pmapRandom(scatterState); - d = d2 <= FTINY + d = d2 <= FTINY ? 1 - : sqrt(-log(d2) / - (sqr(cosp) / sqr(nd -> u_alpha) + + : sqrt(-log(d2) / + (sqr(cosp) / sqr(nd -> u_alpha) + sqr(sinp) / (nd -> v_alpha * nd -> u_alpha) ) ); for (i = 0; i < 3; i++) - rayOut -> rdir [i] = nd -> prdir [i] + + rayOut -> rdir [i] = nd -> prdir [i] + d * (cosp * nd -> u [i] + sinp * nd -> v [i]); if (DOT(rayOut -> rdir, rayIn -> ron) < -FTINY) { normalize(rayOut -> rdir); return 1; } } - + return 0; } - + else { /* Specular reflection */ - + /* Make MAXITER attempts at getting a ray */ for (niter = 0; niter < MAXITER; niter++) { d = 2 * PI * pmapRandom(scatterState); cosp = cos(d) * nd -> u_alpha; sinp = sin(d) * nd -> v_alpha; d = sqrt(sqr(cosp) + sqr(sinp)); cosp /= d; sinp /= d; d2 = pmapRandom(scatterState); - d = d2 <= FTINY - ? 1 - : sqrt(-log(d2) / + d = d2 <= FTINY + ? 1 + : sqrt(-log(d2) / (sqr(cosp) / sqr(nd -> u_alpha) + sqr(sinp) / (nd->v_alpha * nd->v_alpha) ) ); for (i = 0; i < 3; i++) - h [i] = nd -> pnorm [i] + + h [i] = nd -> pnorm [i] + d * (cosp * nd -> u [i] + sinp * nd -> v [i]); - + d = -2 * DOT(h, rayIn -> rdir) / (1 + d * d); VSUM(rayOut -> rdir, rayIn -> rdir, h, d); - - if (DOT(rayOut -> rdir, rayIn -> ron) > FTINY) + + if (DOT(rayOut -> rdir, rayIn -> ron) > FTINY) return 1; } - + return 0; } } static int anisoPhotonScatter (OBJREC *mat, RAY *rayIn) /* Generate new photon ray for anisotropic material and recurse */ { ANISODAT nd; float xi, albedo, prdiff, ptdiff, prspec, ptspec; RAY rayOut; - + if (mat -> oargs.nfargs != (mat -> otype == MAT_TRANS2 ? 8 : 6)) objerror(mat, USER, "bad number of real arguments"); - + nd.mp = mat; nd.rp = rayIn; - + /* get material color */ copycolor(nd.mcolor, mat -> oargs.farg); - + /* get roughness */ nd.specfl = 0; nd.u_alpha = mat -> oargs.farg [4]; nd.v_alpha = mat -> oargs.farg [5]; if (nd.u_alpha < FTINY || nd.v_alpha <= FTINY) objerror(mat, USER, "roughness too small"); - + /* check for back side; reorient if back is visible */ if (rayIn -> rod < 0) - if (!backvis && mat -> otype != MAT_TRANS2) + if (!backvis && mat -> otype != MAT_TRANS2) return 0; else { /* get modifiers */ raytexture(rayIn, mat -> omod); flipsurface(rayIn); } else raytexture(rayIn, mat -> omod); - + /* perturb normal */ nd.pdot = max(raynormal(nd.pnorm, rayIn), .001); - + /* modify material color */ multcolor(nd.mcolor, rayIn -> pcol); nd.rspec = mat -> oargs.farg [3]; - + /* transmission params */ if (mat -> otype == MAT_TRANS2) { nd.trans = mat -> oargs.farg [6] * (1 - nd.rspec); nd.tspec = nd.trans * mat -> oargs.farg [7]; nd.tdiff = nd.trans - nd.tspec; - if (nd.tspec > FTINY) + if (nd.tspec > FTINY) nd.specfl |= SP_TRAN; - } + } else nd.tdiff = nd.tspec = nd.trans = 0; - + /* specular reflection params */ if (nd.rspec > FTINY) { nd.specfl |= SP_REFL; - + /* compute specular color */ - if (mat -> otype == MAT_METAL2) + if (mat -> otype == MAT_METAL2) copycolor(nd.scolor, nd.mcolor); else setcolor(nd.scolor, 1, 1, 1); - + scalecolor(nd.scolor, nd.rspec); } else setcolor(nd.scolor, 0, 0, 0); - + /* diffuse reflection params */ nd.rdiff = 1 - nd.trans - nd.rspec; - + /* Set up probabilities */ prdiff = ptdiff = ptspec = colorAvg(nd.mcolor); prdiff *= nd.rdiff; ptdiff *= nd.tdiff; prspec = colorAvg(nd.scolor); ptspec *= nd.tspec; albedo = prdiff + ptdiff + prspec + ptspec; - + /* Insert direct and indirect photon hits if diffuse component */ if (prdiff > FTINY || ptdiff > FTINY) addPhotons(rayIn); - + xi = pmapRandom(rouletteState); - - if (xi > albedo) + + if (xi > albedo) /* Absorbed */ return 0; - + if (xi > (albedo -= prspec)) /* Specular reflection */ if (!(nd.specfl & SP_BADU)) { photonRay(rayIn, &rayOut, PMAP_SPECREFL, nd.scolor); - if (!anisoSpecPhotonScatter(&nd, &rayOut)) + if (!anisoSpecPhotonScatter(&nd, &rayOut)) return 0; } else return 0; - + else if (xi > (albedo -= ptspec)) /* Specular transmission */ - + if (!(nd.specfl & SP_BADU)) { /* Specular transmission */ photonRay(rayIn, &rayOut, PMAP_SPECTRANS, nd.mcolor); - if (!anisoSpecPhotonScatter(&nd, &rayOut)) + if (!anisoSpecPhotonScatter(&nd, &rayOut)) return 0; } else return 0; - + else if (xi > (albedo -= prdiff)) { - /* Diffuse reflection */ + /* Diffuse reflection */ photonRay(rayIn, &rayOut, PMAP_DIFFREFL, nd.mcolor); diffPhotonScatter(nd.pnorm, &rayOut); } - + else { /* Diffuse transmission */ FVECT bnorm; flipsurface(rayIn); bnorm [0] = -nd.pnorm [0]; bnorm [1] = -nd.pnorm [1]; bnorm [2] = -nd.pnorm [2]; - + photonRay(rayIn, &rayOut, PMAP_DIFFTRANS, nd.mcolor); diffPhotonScatter(bnorm, &rayOut); } - + tracePhoton(&rayOut); return 0; } static double mylog (double x) /* special log for extinction coefficients; cloned from dielectric.c */ { if (x < 1e-40) return(-100.); - + if (x >= 1.) return(0.); - + return(log(x)); } static int dielectricPhotonScatter (OBJREC *mat, RAY *rayIn) /* Generate new photon ray for dielectric material and recurse */ { double cos1, cos2, nratio, d1, d2, refl; COLOR ctrans, talb; FVECT dnorm; int hastexture, i; RAY rayOut; if (mat -> oargs.nfargs != (mat -> otype == MAT_DIELECTRIC ? 5 : 8)) objerror(mat, USER, "bad arguments"); - + /* get modifiers */ raytexture(rayIn, mat -> omod); - + if ((hastexture = (DOT(rayIn -> pert, rayIn -> pert) > sqr(FTINY)))) /* Perturb normal */ cos1 = raynormal(dnorm, rayIn); else { VCOPY(dnorm, rayIn -> ron); cos1 = rayIn -> rod; } - + /* index of refraction */ - nratio = mat -> otype == MAT_DIELECTRIC + nratio = mat -> otype == MAT_DIELECTRIC ? mat->oargs.farg[3] + mat->oargs.farg[4] / MLAMBDA : mat->oargs.farg[3] / mat->oargs.farg[7]; if (cos1 < 0) { /* inside */ hastexture = -hastexture; cos1 = -cos1; dnorm [0] = -dnorm [0]; dnorm [1] = -dnorm [1]; dnorm [2] = -dnorm [2]; - setcolor(rayIn -> cext, + setcolor(rayIn -> cext, -mylog(mat -> oargs.farg [0] * rayIn -> pcol [0]), -mylog(mat -> oargs.farg [1] * rayIn -> pcol [1]), -mylog(mat -> oargs.farg [2] * rayIn -> pcol [2]) ); setcolor(rayIn -> albedo, 0, 0, 0); rayIn -> gecc = 0; - + if (mat -> otype == MAT_INTERFACE) { - setcolor(ctrans, + setcolor(ctrans, -mylog(mat -> oargs.farg [4] * rayIn -> pcol [0]), -mylog(mat -> oargs.farg [5] * rayIn -> pcol [1]), -mylog(mat -> oargs.farg [6] * rayIn -> pcol [2]) ); setcolor(talb, 0, 0, 0); - } + } else { copycolor(ctrans, cextinction); copycolor(talb, salbedo); } - } - + } + else { /* outside */ nratio = 1.0 / nratio; - setcolor(ctrans, + setcolor(ctrans, -mylog(mat -> oargs.farg [0] * rayIn -> pcol [0]), -mylog(mat -> oargs.farg [1] * rayIn -> pcol [1]), -mylog(mat -> oargs.farg [2] * rayIn -> pcol [2]) ); setcolor(talb, 0, 0, 0); - + if (mat -> otype == MAT_INTERFACE) { setcolor(rayIn -> cext, -mylog(mat -> oargs.farg [4] * rayIn -> pcol [0]), -mylog(mat -> oargs.farg [5] * rayIn -> pcol [1]), -mylog(mat -> oargs.farg [6] * rayIn -> pcol [2]) ); setcolor(rayIn -> albedo, 0, 0, 0); rayIn -> gecc = 0; } - } + } /* compute cos theta2 */ d2 = 1 - sqr(nratio) * (1 - sqr(cos1)); - + if (d2 < FTINY) { /* Total reflection */ refl = cos2 = 1.0; } else { /* Refraction, compute Fresnel's equations */ cos2 = sqrt(d2); d1 = cos1; d2 = nratio * cos2; d1 = (d1 - d2) / (d1 + d2); refl = sqr(d1); d1 = 1 / cos1; d2 = nratio / cos2; d1 = (d1 - d2) / (d1 + d2); refl += sqr(d1); refl *= 0.5; } - + if (pmapRandom(rouletteState) > refl) { /* Refraction */ photonRay(rayIn, &rayOut, PMAP_REFRACT, NULL); d1 = nratio * cos1 - cos2; - + for (i = 0; i < 3; i++) rayOut.rdir [i] = nratio * rayIn -> rdir [i] + d1 * dnorm [i]; - + if (hastexture && DOT(rayOut.rdir, rayIn->ron)*hastexture >= -FTINY) { d1 *= hastexture; for (i = 0; i < 3; i++) - rayOut.rdir [i] = nratio * rayIn -> rdir [i] + + rayOut.rdir [i] = nratio * rayIn -> rdir [i] + d1 * rayIn -> ron [i]; - + normalize(rayOut.rdir); } copycolor(rayOut.cext, ctrans); copycolor(rayOut.albedo, talb); } - + else { /* Reflection */ - photonRay(rayIn, &rayOut, PMAP_SPECREFL, NULL); + photonRay(rayIn, &rayOut, REFLECTED, NULL); VSUM(rayOut.rdir, rayIn -> rdir, dnorm, 2 * cos1); - + if (hastexture && DOT(rayOut.rdir, rayIn->ron) * hastexture <= FTINY) for (i = 0; i < 3; i++) - rayOut.rdir [i] = rayIn -> rdir [i] + + rayOut.rdir [i] = rayIn -> rdir [i] + 2 * rayIn -> rod * rayIn -> ron [i]; } - + /* Ray is modified by medium defined by cext and albedo in * photonParticipate() */ tracePhoton(&rayOut); - + return 0; } static int glassPhotonScatter (OBJREC *mat, RAY *rayIn) /* Generate new photon ray for glass material and recurse */ { float albedo, xi, ptrans; COLOR mcolor, refl, trans; double pdot, cos2, d, r1e, r1m, rindex = 0.0; FVECT pnorm, pdir; int hastexture, i; RAY rayOut; /* check arguments */ if (mat -> oargs.nfargs == 3) rindex = RINDEX; else if (mat -> oargs.nfargs == 4) rindex = mat -> oargs.farg [3]; else objerror(mat, USER, "bad arguments"); - copycolor(mcolor, mat -> oargs.farg); - + copycolor(mcolor, mat -> oargs.farg); + /* get modifiers */ raytexture(rayIn, mat -> omod); - + /* reorient if necessary */ - if (rayIn -> rod < 0) + if (rayIn -> rod < 0) flipsurface(rayIn); if ((hastexture = (DOT(rayIn -> pert, rayIn -> pert) > sqr(FTINY)))) pdot = raynormal(pnorm, rayIn); else { VCOPY(pnorm, rayIn -> ron); pdot = rayIn -> rod; } - + /* Modify material color */ multcolor(mcolor, rayIn -> pcol); - + /* angular transmission */ cos2 = sqrt((1 - 1 / sqr(rindex)) + sqr(pdot / rindex)); - setcolor(mcolor, - pow(mcolor [0], 1 / cos2), - pow(mcolor [1], 1 / cos2), + setcolor(mcolor, + pow(mcolor [0], 1 / cos2), + pow(mcolor [1], 1 / cos2), pow(mcolor [2], 1 / cos2) ); - + /* compute reflection */ r1e = (pdot - rindex * cos2) / (pdot + rindex * cos2); r1e *= r1e; r1m = (1 / pdot - rindex / cos2) / (1 / pdot + rindex / cos2); r1m *= r1m; - + for (i = 0; i < 3; i++) { double r1ed2, r1md2, d2; - + d = mcolor [i]; d2 = sqr(d); r1ed2 = sqr(r1e) * d2; r1md2 = sqr(r1m) * d2; - + /* compute transmittance */ - trans [i] = 0.5 * d * (sqr(1 - r1e) / + trans [i] = 0.5 * d * (sqr(1 - r1e) / (1 - r1ed2) + sqr(1 - r1m) / (1 - r1md2) ); /* compute reflectance */ refl [i] = 0.5 * (r1e * (1 + (1 - 2 * r1e) * d2) / (1 - r1ed2) + r1m * (1 + (1 - 2 * r1m) * d2) / (1 - r1md2) ); } - + /* Set up probabilities */ ptrans = colorAvg(trans); albedo = colorAvg(refl) + ptrans; xi = pmapRandom(rouletteState); - if (xi > albedo) + if (xi > albedo) /* Absorbed */ return 0; if (xi > (albedo -= ptrans)) { /* Transmitted */ if (hastexture) { /* perturb direction */ VSUM(pdir, rayIn -> rdir, rayIn -> pert, 2 * (1 - rindex)); if (normalize(pdir) == 0) { objerror(mat, WARNING, "bad perturbation"); VCOPY(pdir, rayIn -> rdir); } - } + } else VCOPY(pdir, rayIn -> rdir); - + VCOPY(rayOut.rdir, pdir); - photonRay(rayIn, &rayOut, PMAP_SPECTRANS, mcolor); + photonRay(rayIn, &rayOut, TRANS, mcolor); } - + else { /* reflected ray */ VSUM(rayOut.rdir, rayIn -> rdir, pnorm, 2 * pdot); - photonRay(rayIn, &rayOut, PMAP_SPECREFL, mcolor); + photonRay(rayIn, &rayOut, REFLECTED, mcolor); } tracePhoton(&rayOut); return 0; } static int aliasPhotonScatter (OBJREC *mat, RAY *rayIn) /* Transfer photon scattering to alias target */ { OBJECT aliasObj; OBJREC aliasRec, *aliasPtr; - + /* Straight replacement? */ if (!mat -> oargs.nsargs) { /* Skip void modifier! */ - if (mat -> omod != OVOID) { + if (mat -> omod != OVOID) { mat = objptr(mat -> omod); photonScatter [mat -> otype] (mat, rayIn); } - + return 0; } - + /* Else replace alias */ if (mat -> oargs.nsargs != 1) objerror(mat, INTERNAL, "bad # string arguments"); - + aliasPtr = mat; aliasObj = objndx(aliasPtr); - + /* Follow alias trail */ do { - aliasObj = aliasPtr -> oargs.nsargs == 1 + aliasObj = aliasPtr -> oargs.nsargs == 1 ? lastmod(aliasObj, aliasPtr -> oargs.sarg [0]) : aliasPtr -> omod; if (aliasObj < 0) objerror(aliasPtr, USER, "bad reference"); - + aliasPtr = objptr(aliasObj); } while (aliasPtr -> otype == MOD_ALIAS); /* Copy alias object */ aliasRec = *aliasPtr; - + /* Substitute modifier */ aliasRec.omod = mat -> omod; - + /* Replacement scattering routine */ photonScatter [aliasRec.otype] (&aliasRec, rayIn); /* Avoid potential memory leak? */ if (aliasRec.os != aliasPtr -> os) { if (aliasPtr -> os) free_os(aliasPtr); aliasPtr -> os = aliasRec.os; } return 0; } static int clipPhotonScatter (OBJREC *mat, RAY *rayIn) /* Generate new photon ray for antimatter material and recurse */ { OBJECT obj = objndx(mat), mod, cset [MAXSET + 1], *modset; int entering, inside = 0, i; const RAY *rp; RAY rayOut; - if ((modset = (OBJECT*)mat -> os) == NULL) { + if ((modset = (OBJECT*)mat -> os) == NULL) { if (mat -> oargs.nsargs < 1 || mat -> oargs.nsargs > MAXSET) objerror(mat, USER, "bad # arguments"); - + modset = (OBJECT*)malloc((mat -> oargs.nsargs + 1) * sizeof(OBJECT)); - - if (modset == NULL) + + if (modset == NULL) error(SYSTEM, "out of memory in clipPhotonScatter"); modset [0] = 0; - + for (i = 0; i < mat -> oargs.nsargs; i++) { if (!strcmp(mat -> oargs.sarg [i], VOIDID)) continue; - + if ((mod = lastmod(obj, mat -> oargs.sarg [i])) == OVOID) { sprintf(errmsg, "unknown modifier \"%s\"", mat->oargs.sarg[i]); objerror(mat, WARNING, errmsg); continue; } - + if (inset(modset, mod)) { objerror(mat, WARNING, "duplicate modifier"); continue; } - + insertelem(modset, mod); } - + mat -> os = (char*)modset; } - - if (rayIn -> clipset != NULL) + + if (rayIn -> clipset != NULL) setcopy(cset, rayIn -> clipset); else cset [0] = 0; - + entering = rayIn -> rod > 0; - + if (inset(photonSensorSet, obj)) { - /* Material defined as virtual sensor; store photon if incident + /* Material defined as virtual sensor; store photon if incident * from side(s) corresponding to orientation flags. Note PMAPSENSBI * subsumes both PMAP_SENSFWD and PMAP_SENSBWD. */ const char sensFlags = PMAP_GETSENSFLAGS(mat); if (entering && sensFlags & PMAP_SENSFWD || !entering && sensFlags & PMAP_SENSBWD ) if (entering) addPhotons(rayIn); else { - /* Flip normal if incident from back, then flip back after + /* Flip normal if incident from back, then flip back after saving photon */ flipsurface(rayIn); addPhotons(rayIn); flipsurface(rayIn); } } for (i = modset [0]; i > 0; i--) { if (entering) { if (!inset(cset, modset [i])) { - if (cset [0] >= MAXSET) + if (cset [0] >= MAXSET) error(INTERNAL, "set overflow in clipPhotonScatter"); insertelem(cset, modset [i]); } - } - else if (inset(cset, modset [i])) + } + else if (inset(cset, modset [i])) deletelem(cset, modset [i]); } - + rayIn -> newcset = cset; - + if (strcmp(mat -> oargs.sarg [0], VOIDID)) { for (rp = rayIn; rp -> parent != NULL; rp = rp -> parent) { - if ( !(rp -> rtype & RAYREFL) && rp->parent->ro != NULL && + if ( !(rp -> rtype & RAYREFL) && rp->parent->ro != NULL && inset(modset, rp -> parent -> ro -> omod)) { if (rp -> parent -> rod > 0) inside++; else inside--; } } - + if (inside > 0) { flipsurface(rayIn); mat = objptr(lastmod(obj, mat -> oargs.sarg [0])); photonScatter [mat -> otype] (mat, rayIn); return 0; } } - + /* Else transfer ray */ photonRay(rayIn, &rayOut, PMAP_XFER, NULL); tracePhoton(&rayOut); - + return 0; } static int mirrorPhotonScatter (OBJREC *mat, RAY *rayIn) /* Generate new photon ray for mirror material and recurse */ { RAY rayOut; int rpure = 1, i; FVECT pnorm; double pdot; float albedo; COLOR mcolor; /* check arguments */ if (mat -> oargs.nfargs != 3 || mat -> oargs.nsargs > 1) objerror(mat, USER, "bad number of arguments"); - + /* back is black */ - if (rayIn -> rod < 0) + if (rayIn -> rod < 0) return 0; - + /* get modifiers */ raytexture(rayIn, mat -> omod); - + /* assign material color */ - copycolor(mcolor, mat -> oargs.farg); + copycolor(mcolor, mat -> oargs.farg); multcolor(mcolor, rayIn -> pcol); - + /* Set up probabilities */ albedo = colorAvg(mcolor); - - if (pmapRandom(rouletteState) > albedo) + + if (pmapRandom(rouletteState) > albedo) /* Absorbed */ return 0; - + /* compute reflected ray */ - photonRay(rayIn, &rayOut, PMAP_SPECREFL, mcolor); + photonRay(rayIn, &rayOut, REFLECTED, mcolor); if (DOT(rayIn -> pert, rayIn -> pert) > sqr(FTINY)) { /* use textures */ - pdot = raynormal(pnorm, rayIn); - + pdot = raynormal(pnorm, rayIn); + for (i = 0; i < 3; i++) rayOut.rdir [i] = rayIn -> rdir [i] + 2 * pdot * pnorm [i]; - + rpure = 0; } - + /* Check for penetration */ if (rpure || DOT(rayOut.rdir, rayIn -> ron) <= FTINY) for (i = 0; i < 3; i++) - rayOut.rdir [i] = rayIn -> rdir [i] + + rayOut.rdir [i] = rayIn -> rdir [i] + 2 * rayIn -> rod * rayIn -> ron [i]; tracePhoton(&rayOut); return 0; } static int mistPhotonScatter (OBJREC *mat, RAY *rayIn) /* Generate new photon ray within mist and recurse */ { COLOR mext; RREAL re, ge, be; RAY rayOut; /* check arguments */ - if (mat -> oargs.nfargs > 7) + if (mat -> oargs.nfargs > 7) objerror(mat, USER, "bad arguments"); - + if (mat -> oargs.nfargs > 2) { /* compute extinction */ copycolor(mext, mat -> oargs.farg); /* get modifiers */ raytexture(rayIn, mat -> omod); multcolor(mext, rayIn -> pcol); - } + } else setcolor(mext, 0, 0, 0); - + photonRay(rayIn, &rayOut, PMAP_XFER, NULL); - + if (rayIn -> rod > 0) { /* entering ray */ addcolor(rayOut.cext, mext); - - if (mat -> oargs.nfargs > 5) + + if (mat -> oargs.nfargs > 5) copycolor(rayOut.albedo, mat -> oargs.farg + 3); - if (mat -> oargs.nfargs > 6) + if (mat -> oargs.nfargs > 6) rayOut.gecc = mat -> oargs.farg [6]; - } - + } + else { /* leaving ray */ re = max(rayIn -> cext [0] - mext [0], cextinction [0]); ge = max(rayIn -> cext [1] - mext [1], cextinction [1]); be = max(rayIn -> cext [2] - mext [2], cextinction [2]); setcolor(rayOut.cext, re, ge, be); - - if (mat -> oargs.nfargs > 5) + + if (mat -> oargs.nfargs > 5) copycolor(rayOut.albedo, salbedo); if (mat -> oargs.nfargs > 6) rayOut.gecc = seccg; } - + tracePhoton(&rayOut); return 0; } static int mx_dataPhotonScatter (OBJREC *mat, RAY *rayIn) /* Pass photon on to materials selected by mixture data */ { OBJECT obj; double coef, pt [MAXDIM]; DATARRAY *dp; OBJECT mod [2]; MFUNC *mf; int i; - if (mat -> oargs.nsargs < 6) + if (mat -> oargs.nsargs < 6) objerror(mat, USER, "bad # arguments"); - + obj = objndx(mat); - + for (i = 0; i < 2; i++) - if (!strcmp(mat -> oargs.sarg [i], VOIDID)) + if (!strcmp(mat -> oargs.sarg [i], VOIDID)) mod [i] = OVOID; else if ((mod [i] = lastmod(obj, mat -> oargs.sarg [i])) == OVOID) { sprintf(errmsg, "undefined modifier \"%s\"", mat->oargs.sarg[i]); objerror(mat, USER, errmsg); } - + dp = getdata(mat -> oargs.sarg [3]); i = (1 << dp -> nd) - 1; mf = getfunc(mat, 4, i << 5, 0); setfunc(mat, rayIn); errno = 0; - + for (i = 0; i < dp -> nd; i++) { pt [i] = evalue(mf -> ep [i]); - + if (errno) { objerror(mat, WARNING, "compute error"); return 0; } } - + coef = datavalue(dp, pt); errno = 0; coef = funvalue(mat -> oargs.sarg [2], 1, &coef); - - if (errno) + + if (errno) objerror(mat, WARNING, "compute error"); else { OBJECT mxMod = mod [pmapRandom(rouletteState) < coef ? 0 : 1]; - + if (mxMod != OVOID) { mat = objptr(mxMod); photonScatter [mat -> otype] (mat, rayIn); } else { /* Transfer ray if no modifier */ RAY rayOut; - + photonRay(rayIn, &rayOut, PMAP_XFER, NULL); tracePhoton(&rayOut); } } - + return 0; } static int mx_pdataPhotonScatter (OBJREC *mat, RAY *rayIn) /* Pass photon on to materials selected by mixture picture */ { OBJECT obj; double col [3], coef, pt [MAXDIM]; DATARRAY *dp; OBJECT mod [2]; MFUNC *mf; int i; - if (mat -> oargs.nsargs < 7) + if (mat -> oargs.nsargs < 7) objerror(mat, USER, "bad # arguments"); - + obj = objndx(mat); - + for (i = 0; i < 2; i++) - if (!strcmp(mat -> oargs.sarg [i], VOIDID)) + if (!strcmp(mat -> oargs.sarg [i], VOIDID)) mod [i] = OVOID; else if ((mod [i] = lastmod(obj, mat -> oargs.sarg [i])) == OVOID) { sprintf(errmsg, "undefined modifier \"%s\"", mat->oargs.sarg[i]); objerror(mat, USER, errmsg); } - + dp = getpict(mat -> oargs.sarg [3]); mf = getfunc(mat, 4, 0x3 << 5, 0); setfunc(mat, rayIn); errno = 0; pt [1] = evalue(mf -> ep [0]); pt [0] = evalue(mf -> ep [1]); - + if (errno) { objerror(mat, WARNING, "compute error"); return 0; } - - for (i = 0; i < 3; i++) + + for (i = 0; i < 3; i++) col [i] = datavalue(dp + i, pt); - + errno = 0; coef = funvalue(mat -> oargs.sarg [2], 3, col); - - if (errno) + + if (errno) objerror(mat, WARNING, "compute error"); else { OBJECT mxMod = mod [pmapRandom(rouletteState) < coef ? 0 : 1]; - + if (mxMod != OVOID) { mat = objptr(mxMod); photonScatter [mat -> otype] (mat, rayIn); } else { /* Transfer ray if no modifier */ RAY rayOut; - + photonRay(rayIn, &rayOut, PMAP_XFER, NULL); tracePhoton(&rayOut); } } return 0; } static int mx_funcPhotonScatter (OBJREC *mat, RAY *rayIn) /* Pass photon on to materials selected by mixture function */ { OBJECT obj, mod [2]; int i; double coef; MFUNC *mf; - if (mat -> oargs.nsargs < 4) + if (mat -> oargs.nsargs < 4) objerror(mat, USER, "bad # arguments"); - + obj = objndx(mat); - + for (i = 0; i < 2; i++) if (!strcmp(mat -> oargs.sarg [i], VOIDID)) mod [i] = OVOID; else if ((mod [i] = lastmod(obj, mat -> oargs.sarg [i])) == OVOID) { sprintf(errmsg, "undefined modifier \"%s\"", mat->oargs.sarg[i]); objerror(mat, USER, errmsg); } - + mf = getfunc(mat, 3, 0x4, 0); setfunc(mat, rayIn); errno = 0; - + /* bound coefficient */ coef = min(1, max(0, evalue(mf -> ep [0]))); - - if (errno) + + if (errno) objerror(mat, WARNING, "compute error"); - else { + else { OBJECT mxMod = mod [pmapRandom(rouletteState) < coef ? 0 : 1]; - + if (mxMod != OVOID) { mat = objptr(mxMod); photonScatter [mat -> otype] (mat, rayIn); } else { /* Transfer ray if no modifier */ RAY rayOut; - + photonRay(rayIn, &rayOut, PMAP_XFER, NULL); tracePhoton(&rayOut); } } - + return 0; } static int pattexPhotonScatter (OBJREC *mat, RAY *rayIn) /* Generate new photon ray for pattern or texture modifier and recurse. This code is brought to you by Henkel! :^) */ { RAY rayOut; - + /* Get pattern */ ofun [mat -> otype].funp(mat, rayIn); if (mat -> omod != OVOID) { /* Scatter using modifier (if any) */ mat = objptr(mat -> omod); photonScatter [mat -> otype] (mat, rayIn); } else { /* Transfer ray if no modifier */ photonRay(rayIn, &rayOut, PMAP_XFER, NULL); tracePhoton(&rayOut); } - + return 0; } static int setbrdfunc(BRDFDAT *bd) /* Set up brdf function and variables; ripped off from m_brdf.c */ { FVECT v; - + if (setfunc(bd -> mp, bd -> pr) == 0) return 0; /* (Re)Assign func variables */ multv3(v, bd -> pnorm, funcxf.xfm); varset("NxP", '=', v [0] / funcxf.sca); varset("NyP", '=', v [1] / funcxf.sca); varset("NzP", '=', v [2] / funcxf.sca); - varset("RdotP", '=', + varset("RdotP", '=', bd -> pdot <= -1. ? -1. : bd -> pdot >= 1. ? 1. : bd -> pdot); varset("CrP", '=', colval(bd -> mcolor, RED)); varset("CgP", '=', colval(bd -> mcolor, GRN)); varset("CbP", '=', colval(bd -> mcolor, BLU)); - + return 1; } static int brdfPhotonScatter (OBJREC *mat, RAY *rayIn) /* Generate new photon ray for BRTDfunc material and recurse. Only ideal reflection and transmission are sampled for the specular componentent. */ { int hitfront = 1, hastexture, i; BRDFDAT nd; RAY rayOut; COLOR rspecCol, tspecCol; double prDiff, ptDiff, prSpec, ptSpec, albedo, xi; MFUNC *mf; FVECT bnorm; /* Check argz */ if (mat -> oargs.nsargs < 10 || mat -> oargs.nfargs < 9) objerror(mat, USER, "bad # arguments"); - + nd.mp = mat; nd.pr = rayIn; /* Dummiez */ nd.rspec = nd.tspec = 1.0; nd.trans = 0.5; /* Diffuz reflektanz */ if (rayIn -> rod > 0.0) setcolor(nd.rdiff, mat -> oargs.farg[0], mat -> oargs.farg [1], mat -> oargs.farg [2]); else setcolor(nd.rdiff, mat-> oargs.farg [3], mat -> oargs.farg [4], mat -> oargs.farg [5]); /* Diffuz tranzmittanz */ setcolor(nd.tdiff, mat -> oargs.farg [6], mat -> oargs.farg [7], mat -> oargs.farg [8]); /* Get modz */ raytexture(rayIn, mat -> omod); hastexture = (DOT(rayIn -> pert, rayIn -> pert) > sqr(FTINY)); if (hastexture) { /* Perturb normal */ nd.pdot = raynormal(nd.pnorm, rayIn); - } + } else { VCOPY(nd.pnorm, rayIn -> ron); nd.pdot = rayIn -> rod; } if (rayIn -> rod < 0.0) { /* Orient perturbed valuz */ nd.pdot = -nd.pdot; for (i = 0; i < 3; i++) { nd.pnorm [i] = -nd.pnorm [i]; rayIn -> pert [i] = -rayIn -> pert [i]; } - + hitfront = 0; } - + /* Get pattern kolour, modify diffuz valuz */ copycolor(nd.mcolor, rayIn -> pcol); multcolor(nd.rdiff, nd.mcolor); multcolor(nd.tdiff, nd.mcolor); /* Load cal file, evaluate spekula refl/tranz varz */ nd.dp = NULL; mf = getfunc(mat, 9, 0x3f, 0); setbrdfunc(&nd); errno = 0; - setcolor(rspecCol, + setcolor(rspecCol, evalue(mf->ep[0]), evalue(mf->ep[1]), evalue(mf->ep[2])); - setcolor(tspecCol, + setcolor(tspecCol, evalue(mf->ep[3]), evalue(mf->ep[4]), evalue(mf->ep[5])); if (errno == EDOM || errno == ERANGE) objerror(mat, WARNING, "compute error"); else { /* Set up probz */ prDiff = colorAvg(nd.rdiff); ptDiff = colorAvg(nd.tdiff); prSpec = colorAvg(rspecCol); ptSpec = colorAvg(tspecCol); albedo = prDiff + ptDiff + prSpec + ptSpec; } /* Insert direct and indirect photon hitz if diffuz komponent */ if (prDiff > FTINY || ptDiff > FTINY) addPhotons(rayIn); /* Stochastically sample absorption or scattering evenz */ if ((xi = pmapRandom(rouletteState)) > albedo) /* Absorbed */ return 0; if (xi > (albedo -= prSpec)) { /* Ideal spekula reflekzion */ - photonRay(rayIn, &rayOut, PMAP_SPECREFL, rspecCol); + photonRay(rayIn, &rayOut, REFLECTED, rspecCol); VSUM(rayOut.rdir, rayIn -> rdir, nd.pnorm, 2 * nd.pdot); checknorm(rayOut.rdir); } else if (xi > (albedo -= ptSpec)) { /* Ideal spekula tranzmission */ - photonRay(rayIn, &rayOut, PMAP_SPECTRANS, tspecCol); + photonRay(rayIn, &rayOut, TRANS, tspecCol); if (hastexture) { /* Perturb direkzion */ VSUB(rayOut.rdir, rayIn -> rdir, rayIn -> pert); if (normalize(rayOut.rdir) == 0.0) { objerror(mat, WARNING, "illegal perturbation"); VCOPY(rayOut.rdir, rayIn -> rdir); } } else VCOPY(rayOut.rdir, rayIn -> rdir); } else if (xi > (albedo -= prDiff)) { /* Diffuz reflekzion */ if (!hitfront) flipsurface(rayIn); photonRay(rayIn, &rayOut, PMAP_DIFFREFL, nd.mcolor); diffPhotonScatter(nd.pnorm, &rayOut); } else { /* Diffuz tranzmission */ if (hitfront) flipsurface(rayIn); photonRay(rayIn, &rayOut, PMAP_DIFFTRANS, nd.mcolor); bnorm [0] = -nd.pnorm [0]; bnorm [1] = -nd.pnorm [1]; bnorm [2] = -nd.pnorm [2]; - diffPhotonScatter(bnorm, &rayOut); + diffPhotonScatter(bnorm, &rayOut); } tracePhoton(&rayOut); return 0; } int brdf2PhotonScatter (OBJREC *mat, RAY *rayIn) /* Generate new photon ray for procedural or data driven BRDF material and recurse. Only diffuse reflection and transmission are sampled. */ { BRDFDAT nd; RAY rayOut; double dtmp, prDiff, ptDiff, albedo, xi; MFUNC *mf; FVECT bnorm; /* Check argz */ - if (mat -> oargs.nsargs < (hasdata(mat -> otype) ? 4 : 2) || - mat -> oargs.nfargs < (mat -> otype == MAT_TFUNC || + if (mat -> oargs.nsargs < (hasdata(mat -> otype) ? 4 : 2) || + mat -> oargs.nfargs < (mat -> otype == MAT_TFUNC || mat -> otype == MAT_TDATA ? 6 : 4)) objerror(mat, USER, "bad # arguments"); - + if (rayIn -> rod < 0.0) { /* Hit backside; reorient if visible, else transfer photon */ if (!backvis) { photonRay(rayIn, &rayOut, PMAP_XFER, NULL); tracePhoton(&rayOut); return 0; } - + raytexture(rayIn, mat -> omod); flipsurface(rayIn); - } + } else raytexture(rayIn, mat -> omod); nd.mp = mat; nd.pr = rayIn; - + /* Material kolour */ - setcolor(nd.mcolor, mat -> oargs.farg [0], mat -> oargs.farg [1], + setcolor(nd.mcolor, mat -> oargs.farg [0], mat -> oargs.farg [1], mat -> oargs.farg [2]); /* Spekula komponent */ nd.rspec = mat -> oargs.farg [3]; - + /* Tranzmittanz */ if (mat -> otype == MAT_TFUNC || mat -> otype == MAT_TDATA) { nd.trans = mat -> oargs.farg [4] * (1.0 - nd.rspec); nd.tspec = nd.trans * mat -> oargs.farg [5]; dtmp = nd.trans - nd.tspec; setcolor(nd.tdiff, dtmp, dtmp, dtmp); - } + } else { nd.tspec = nd.trans = 0.0; setcolor(nd.tdiff, 0.0, 0.0, 0.0); } - + /* Reflektanz */ dtmp = 1.0 - nd.trans - nd.rspec; setcolor(nd.rdiff, dtmp, dtmp, dtmp); /* Perturb normal */ nd.pdot = raynormal(nd.pnorm, rayIn); /* Modify material kolour */ multcolor(nd.mcolor, rayIn -> pcol); multcolor(nd.rdiff, nd.mcolor); multcolor(nd.tdiff, nd.mcolor); - + /* Load auxiliary filez */ if (hasdata(mat -> otype)) { nd.dp = getdata(mat -> oargs.sarg [1]); getfunc(mat, 2, 0, 0); - } + } else { nd.dp = NULL; getfunc(mat, 1, 0, 0); } /* Set up probz */ prDiff = colorAvg(nd.rdiff); ptDiff = colorAvg(nd.tdiff); albedo = prDiff + ptDiff; /* Insert direct and indirect photon hitz if diffuz komponent */ if (prDiff > FTINY || ptDiff > FTINY) addPhotons(rayIn); /* Stochastically sample absorption or scattering evenz */ if ((xi = pmapRandom(rouletteState)) > albedo) /* Absorbed */ return 0; if (xi > (albedo -= prDiff)) { /* Diffuz reflekzion */ photonRay(rayIn, &rayOut, PMAP_DIFFREFL, nd.rdiff); diffPhotonScatter(nd.pnorm, &rayOut); } else { /* Diffuz tranzmission */ flipsurface(rayIn); photonRay(rayIn, &rayOut, PMAP_DIFFTRANS, nd.tdiff); bnorm [0] = -nd.pnorm [0]; bnorm [1] = -nd.pnorm [1]; bnorm [2] = -nd.pnorm [2]; - diffPhotonScatter(bnorm, &rayOut); + diffPhotonScatter(bnorm, &rayOut); } tracePhoton(&rayOut); return 0; } -/* +/* ====================================================================== The following code is (c) Lucerne University of Applied Sciences and Arts, - supported by the Swiss National Science Foundation + supported by the Swiss National Science Foundation (SNSF #147053, "Daylight Redirecting Components") ====================================================================== */ static int bsdfPhotonScatter (OBJREC *mat, RAY *rayIn) /* Generate new photon ray for BSDF modifier and recurse. */ { int hasthick = (mat->otype == MAT_BSDF); int hitFront; SDError err; SDValue bsdfVal; FVECT upvec; MFUNC *mf; BSDFDAT nd; RAY rayOut; COLOR bsdfRGB; int transmitted; - double prDiff, ptDiff, prDiffSD, ptDiffSD, prSpecSD, ptSpecSD, + double prDiff, ptDiff, prDiffSD, ptDiffSD, prSpecSD, ptSpecSD, albedo, xi; const double patAlb = bright(rayIn -> pcol); - + /* Following code adapted from m_bsdf() */ /* Check arguments */ if ( - mat -> oargs.nsargs < hasthick+5 || + mat -> oargs.nsargs < hasthick+5 || mat -> oargs.nfargs > 9 || mat -> oargs.nfargs % 3 ) objerror(mat, USER, "bad # arguments"); - + hitFront = (rayIn -> rod > 0); /* Load cal file */ mf = hasthick ? getfunc(mat, 5, 0x1d, 1) : getfunc(mat, 4, 0xe, 1); /* Get thickness */ nd.thick = 0; if (hasthick) { nd.thick = evalue(mf -> ep [0]); if ((-FTINY <= nd.thick) & (nd.thick <= FTINY)) nd.thick = .0; } /* Get BSDF data */ nd.sd = loadBSDF(mat -> oargs.sarg [hasthick]); - + /* Extra diffuse reflectance from material def */ if (hitFront) { if (mat -> oargs.nfargs < 3) setcolor(nd.rdiff, .0, .0, .0); else setcolor( - nd.rdiff, + nd.rdiff, mat -> oargs.farg [0], mat -> oargs.farg [1], mat -> oargs.farg [2] ); - } + } else if (mat -> oargs.nfargs < 6) { /* Check for absorbing backside */ if (!backvis && !nd.sd -> rb && !nd.sd -> tf) { SDfreeCache(nd.sd); return 0; } - + setcolor(nd.rdiff, .0, .0, .0); - } + } else setcolor( - nd.rdiff, + nd.rdiff, mat -> oargs.farg [3], mat -> oargs.farg [4], mat -> oargs.farg [5] ); /* Extra diffuse transmittance from material def */ if (mat -> oargs.nfargs < 9) setcolor(nd.tdiff, .0, .0, .0); else setcolor( - nd.tdiff, + nd.tdiff, mat -> oargs.farg [6], mat -> oargs.farg [7], mat -> oargs.farg [8] ); - + nd.mp = mat; nd.pr = rayIn; /* Get modifiers */ raytexture(rayIn, mat -> omod); - + /* Modify diffuse values */ multcolor(nd.rdiff, rayIn -> pcol); multcolor(nd.tdiff, rayIn -> pcol); /* Get up vector & xform to world coords */ upvec [0] = evalue(mf -> ep [hasthick+0]); upvec [1] = evalue(mf -> ep [hasthick+1]); upvec [2] = evalue(mf -> ep [hasthick+2]); - + if (mf -> fxp != &unitxf) { multv3(upvec, upvec, mf -> fxp -> xfm); nd.thick *= mf -> fxp -> sca; } - + if (rayIn -> rox) { multv3(upvec, upvec, rayIn -> rox -> f.xfm); nd.thick *= rayIn -> rox -> f.sca; } - + /* Perturb normal */ raynormal(nd.pnorm, rayIn); - + /* Xform incident dir to local BSDF coords */ err = SDcompXform(nd.toloc, nd.pnorm, upvec); - + if (!err) { nd.vray [0] = -rayIn -> rdir [0]; nd.vray [1] = -rayIn -> rdir [1]; nd.vray [2] = -rayIn -> rdir [2]; err = SDmapDir(nd.vray, nd.toloc, nd.vray); } - + if (!err) err = SDinvXform(nd.fromloc, nd.toloc); - + if (err) { objerror(mat, WARNING, "Illegal orientation vector"); return 0; } - + /* Determine BSDF resolution */ err = SDsizeBSDF( nd.sr_vpsa, nd.vray, NULL, SDqueryMin + SDqueryMax, nd.sd ); - + if (err) objerror(mat, USER, transSDError(err)); - + nd.sr_vpsa [0] = sqrt(nd.sr_vpsa [0]); nd.sr_vpsa [1] = sqrt(nd.sr_vpsa [1]); /* Orient perturbed normal towards incident side */ if (!hitFront) { nd.pnorm [0] = -nd.pnorm [0]; nd.pnorm [1] = -nd.pnorm [1]; nd.pnorm [2] = -nd.pnorm [2]; } /* Get scatter probabilities (weighted by pattern except for spec refl) * prDiff, ptDiff: extra diffuse component in material def * prDiffSD, ptDiffSD: diffuse (constant) component in SDF - * prSpecSD, ptSpecSD: non-diffuse ("specular") component in SDF + * prSpecSD, ptSpecSD: non-diffuse ("specular") component in SDF * albedo: sum of above, inverse absorption probability */ prDiff = colorAvg(nd.rdiff); ptDiff = colorAvg(nd.tdiff); prDiffSD = patAlb * SDdirectHemi(nd.vray, SDsampDf | SDsampR, nd.sd); ptDiffSD = patAlb * SDdirectHemi(nd.vray, SDsampDf | SDsampT, nd.sd); prSpecSD = SDdirectHemi(nd.vray, SDsampSp | SDsampR, nd.sd); ptSpecSD = patAlb * SDdirectHemi(nd.vray, SDsampSp | SDsampT, nd.sd); albedo = prDiff + ptDiff + prDiffSD + ptDiffSD + prSpecSD + ptSpecSD; - /* + /* if (albedo > 1) objerror(mat, WARNING, "Invalid albedo"); */ - + /* Insert direct and indirect photon hits if diffuse component */ if (prDiff + ptDiff + prDiffSD + ptDiffSD > FTINY) - addPhotons(rayIn); + addPhotons(rayIn); xi = pmapRandom(rouletteState); - + if (xi > albedo) /* Absorbtion */ return 0; - + transmitted = 0; if ((xi -= prDiff) <= 0) { /* Diffuse reflection (extra component in material def) */ photonRay(rayIn, &rayOut, PMAP_DIFFREFL, nd.rdiff); diffPhotonScatter(nd.pnorm, &rayOut); } - + else if ((xi -= ptDiff) <= 0) { /* Diffuse transmission (extra component in material def) */ photonRay(rayIn, &rayOut, PMAP_DIFFTRANS, nd.tdiff); diffPhotonScatter(nd.pnorm, &rayOut); transmitted = 1; } - + else { /* Sample SDF */ if ((xi -= prDiffSD) <= 0) { /* Diffuse SDF reflection (constant component) */ - if ((err = SDsampBSDF(&bsdfVal, nd.vray, + if ((err = SDsampBSDF(&bsdfVal, nd.vray, pmapRandom(scatterState), SDsampDf | SDsampR, nd.sd ))) objerror(mat, USER, transSDError(err)); - + /* Apply pattern to spectral component */ ccy2rgb(&bsdfVal.spec, bsdfVal.cieY, bsdfRGB); multcolor(bsdfRGB, rayIn -> pcol); photonRay(rayIn, &rayOut, PMAP_DIFFREFL, bsdfRGB); } else if ((xi -= ptDiffSD) <= 0) { /* Diffuse SDF transmission (constant component) */ - if ((err = SDsampBSDF(&bsdfVal, nd.vray, + if ((err = SDsampBSDF(&bsdfVal, nd.vray, pmapRandom(scatterState), SDsampDf | SDsampT, nd.sd ))) objerror(mat, USER, transSDError(err)); - + /* Apply pattern to spectral component */ ccy2rgb(&bsdfVal.spec, bsdfVal.cieY, bsdfRGB); multcolor(bsdfRGB, rayIn -> pcol); addcolor(bsdfRGB, nd.tdiff); photonRay(rayIn, &rayOut, PMAP_DIFFTRANS, bsdfRGB); transmitted = 1; } else if ((xi -= prSpecSD) <= 0) { /* Non-diffuse ("specular") SDF reflection */ - if ((err = SDsampBSDF(&bsdfVal, nd.vray, + if ((err = SDsampBSDF(&bsdfVal, nd.vray, pmapRandom(scatterState), SDsampSp | SDsampR, nd.sd ))) objerror(mat, USER, transSDError(err)); - + ccy2rgb(&bsdfVal.spec, bsdfVal.cieY, bsdfRGB); photonRay(rayIn, &rayOut, PMAP_SPECREFL, bsdfRGB); } - + else { /* Non-diffuse ("specular") SDF transmission */ - if ((err = SDsampBSDF(&bsdfVal, nd.vray, + if ((err = SDsampBSDF(&bsdfVal, nd.vray, pmapRandom(scatterState), SDsampSp | SDsampT, nd.sd ))) objerror(mat, USER, transSDError(err)); /* Apply pattern to spectral component */ ccy2rgb(&bsdfVal.spec, bsdfVal.cieY, bsdfRGB); multcolor(bsdfRGB, rayIn -> pcol); photonRay(rayIn, &rayOut, PMAP_SPECTRANS, bsdfRGB); transmitted = 1; - } - + } + /* Xform outgoing dir to world coords */ if ((err = SDmapDir(rayOut.rdir, nd.fromloc, nd.vray))) { objerror(mat, USER, transSDError(err)); return 0; } } - + /* Clean up */ SDfreeCache(nd.sd); /* Offset outgoing photon origin by thickness to bypass proxy geometry */ if (transmitted && nd.thick != 0) VSUM(rayOut.rorg, rayOut.rorg, rayIn -> ron, -nd.thick); tracePhoton(&rayOut); return 0; } static int lightPhotonScatter (OBJREC* mat, RAY* ray) /* Light sources doan' reflect, mang */ { return 0; } void initPhotonScatterFuncs () /* Init photonScatter[] dispatch table */ { int i; /* Catch-all for inconsistencies */ - for (i = 0; i < NUMOTYPE; i++) + for (i = 0; i < NUMOTYPE; i++) photonScatter [i] = o_default; photonScatter [MAT_LIGHT] = photonScatter [MAT_ILLUM] = - photonScatter [MAT_GLOW] = photonScatter [MAT_SPOT] = + photonScatter [MAT_GLOW] = photonScatter [MAT_SPOT] = lightPhotonScatter; photonScatter [MAT_PLASTIC] = photonScatter [MAT_METAL] = photonScatter [MAT_TRANS] = normalPhotonScatter; - + photonScatter [MAT_PLASTIC2] = photonScatter [MAT_METAL2] = photonScatter [MAT_TRANS2] = anisoPhotonScatter; - - photonScatter [MAT_DIELECTRIC] = photonScatter [MAT_INTERFACE] = + + photonScatter [MAT_DIELECTRIC] = photonScatter [MAT_INTERFACE] = dielectricPhotonScatter; photonScatter [MAT_MIST] = mistPhotonScatter; photonScatter [MAT_GLASS] = glassPhotonScatter; photonScatter [MAT_CLIP] = clipPhotonScatter; photonScatter [MAT_MIRROR] = mirrorPhotonScatter; photonScatter [MIX_FUNC] = mx_funcPhotonScatter; photonScatter [MIX_DATA] = mx_dataPhotonScatter; photonScatter [MIX_PICT]= mx_pdataPhotonScatter; photonScatter [PAT_BDATA] = photonScatter [PAT_CDATA] = photonScatter [PAT_BFUNC] = photonScatter [PAT_CFUNC] = - photonScatter [PAT_CPICT] = photonScatter [TEX_FUNC] = + photonScatter [PAT_CPICT] = photonScatter [TEX_FUNC] = photonScatter [TEX_DATA] = pattexPhotonScatter; photonScatter [MOD_ALIAS] = aliasPhotonScatter; photonScatter [MAT_BRTDF] = brdfPhotonScatter; - + photonScatter [MAT_PFUNC] = photonScatter [MAT_MFUNC] = - photonScatter [MAT_PDATA] = photonScatter [MAT_MDATA] = - photonScatter [MAT_TFUNC] = photonScatter [MAT_TDATA] = + photonScatter [MAT_PDATA] = photonScatter [MAT_MDATA] = + photonScatter [MAT_TFUNC] = photonScatter [MAT_TDATA] = brdf2PhotonScatter; - photonScatter [MAT_BSDF] = photonScatter [MAT_ABSDF] = + photonScatter [MAT_BSDF] = photonScatter [MAT_ABSDF] = bsdfPhotonScatter; } diff --git a/pmapmat.h b/pmapmat.h index 54e7452..e4d1101 100644 --- a/pmapmat.h +++ b/pmapmat.h @@ -1,100 +1,133 @@ /* RCSid $Id: pmapmat.h,v 2.14 2021/02/22 13:27:49 rschregle Exp $ */ -/* +/* ====================================================================== - Photon map support routines for scattering by materials. - + Photon map support routines for scattering by materials. + Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) (c) Fraunhofer Institute for Solar Energy Systems, - supported by the German Research Foundation + supported by the German Research Foundation (DFG LU-204/10-2, "Fassadenintegrierte Regelsysteme FARESYS") (c) Lucerne University of Applied Sciences and Arts, - supported by the Swiss National Science Foundation + supported by the Swiss National Science Foundation (SNSF #147053, "Daylight Redirecting Components") ====================================================================== - + */ #ifndef PMAPMAT_H #define PMAPMAT_H #include "pmap.h" - /* - Check for paths already accounted for in photon map to avoid + /* Assumption: Caustic is anything that cannot be reached by backward + ray-tracing or stochastic sampling -> CAUSTICFLAGS=TRUE */ + #define ISWEAK(r) ( findmaterial(objptr(source[(r)->rsrc].so->omod))->otype==MAT_GLOW||findmaterial(objptr(source[(r)->rsrc].so->omod))->otype==MAT_ILLUM) + #define CAUSTICFLAGS(r) ( ( ( (r)->rtype&SPECULAR || \ + (r)->rtype&REFRACTED || \ + (r)->rtype&REFLECTED || \ + (r)->rtype&TRANS && \ + ( ( (r)->rlvl?(r)->parent->rtype&(SPECULAR|REFRACTED|REFLECTED|AMBIENT) : 0) || ISWEAK(r)) \ + ) && !((r)->rtype&AMBIENT) \ + ) != 0 \ + ) + /* Check for paths already accounted for in photon map to avoid double-counting during backward raytracing. - - ambRayInPmap(): Check for DIFFUSE -> (DIFFUSE|SPECULAR) -> * - subpaths. These occur during the backward pass - when an ambient ray spawns an indirect diffuse or - specular ray. These subpaths are already + + ambRayInPmap(): Check for DIFFUSE -> (DIFFUSE|SPECULAR) -> * + subpaths. These occur during the backward pass + when an ambient ray spawns an indirect diffuse or + specular ray. These subpaths are already accounted for in the global photon map. */ - #define ambRayInPmap(r) ( \ - (r) -> crtype & AMBIENT && photonMapping && \ - (ambounce < 0 || (r) -> rlvl > 1) \ + #define ambRayInPmap(r) ( photonMapping && \ + ( ambounce < 0 || (r)->rlvl > 1 || causticPmap || contribPmap ) && \ + (r)->rtype==AMBIENT ) + + /*Theseareinthecausticphotonmap*/ + #define specularRayInPmap(r) ( (causticPmap || contribPmap ) && \ + ( (r)->crtype == (PRIMARY|AMBIENT|SPECULAR) || \ + (r)->crtype == (PRIMARY|SHADOW|SPECULAR)) || \ + (r)->crtype == (PRIMARY|AMBIENT|TRANS) || \ + (r)->crtype == (PRIMARY|SHADOW|TRANS)) \ + ) + #define ambientRayInPmap(m,r) ( photonMapping && ambounce && \ + (r)->rtype == SHADOW && \ + (r)->crtype == (PRIMARY|AMBIENT|SHADOW) && \ + (r)->parent->crtype != (PRIMARY) && \ + (r)->rlvl == 1 && 0 \ + ) + #define blockAmb(r) ( photonMapping && (ambounce<0) ) + #define viewRayInPmap (m,r) (photonMapping && (ambounce<0) &&0 ) + + /* Block sampling of glow through data-driven BSDF */ + #define glowignore(m,r) (r->rlvl ? \ + ( ( ( globalPmap || contribPmap ) && ambounce<0 ) && \ + (m->otype == MAT_GLOW || m->otype == MAT_ILLUM) && \ + !(r->parent->crtype&(SPECULAR|AMBIENT)) && \ + ( (objptr(r->parent->ro->omod)->otype == MAT_BSDF) || \ + (objptr(r->parent->ro->omod)->otype == MAT_ABSDF) \ + ) \ + ) : 0 \ ) /* !!! shadowRayInPmap() NOW DISABLED !!! */ #define shadowRayInPmap(r) 0 /* srcRayInPmap(): Check whether a source ray transferred through * medium (e.g. glass/dielectric) is already * accounted for in the photon map. This is used by * badcomponent() in source.c when checking source * hits (glows and lights, hence ambient and shadow * rays). */ - #define srcRayInPmap(r) ( \ - ((globalPmap && ambounce < 0) || causticPmap || contribPmap) && \ - (r) -> crtype & (AMBIENT | SHADOW) && \ - (r) -> rtype & (TRANS | REFRACTED) \ - ) + #define srcRayInPmap(r) (0) - /* Check if scattered ray spawns a caustic photon; + /* Check if scattered ray spawns a caustic photon; * !!! NOTE this returns a single bit as boolean value (0|1), rather * !!! than the original short int, hence the explicit test against zero. * !!! This allows the macro the be used in a conditional statement * !!! and when setting a photon's caustic flag in newPhoton(). */ - #define PMAP_CAUSTICRAY(r) (((r) -> rtype & SPECULAR) != 0) - + #define PMAP_CAUSTICRAY(r) ( ( CAUSTICFLAGS(r) ) != 0 ) /* Scattered photon ray types for photonRay() */ #define PMAP_DIFFREFL (REFLECTED | AMBIENT) - #define PMAP_DIFFTRANS (REFLECTED | AMBIENT | TRANS) + #define PMAP_DIFFTRANS (AMBIENT | TRANS) + #define PMAP_REFL (REFLECTED) #define PMAP_SPECREFL (REFLECTED | SPECULAR) - #define PMAP_SPECTRANS (REFLECTED | SPECULAR | TRANS) - #define PMAP_REFRACT (REFRACTED | SPECULAR) + #define PMAP_TRANS (TRANS) + #define PMAP_SPECTRANS (SPECULAR | TRANS) + #define PMAP_REFRACT (REFRACTED) #define PMAP_XFER (TRANS) /* Antimatter sensor flags (orientation relative to surface normal): * Forward, backward, both (bidirectional). */ #define PMAP_SENSFWD 1 #define PMAP_SENSBWD 2 #define PMAP_SENSBI (PMAP_SENSFWD | PMAP_SENSBWD) /* Dispatch table for photon scattering functions */ extern int (*photonScatter []) (OBJREC*, RAY*); /* List of antimatter sensor modifier names */ extern char *photonSensorList [MAXSET + 1]; /* Spawn a new photon ray from a previous one; this is effectively a * customised rayorigin(). */ - void photonRay (const RAY *rayIn, RAY *rayOut, int rayOutType, + void photonRay (const RAY *rayIn, RAY *rayOut, int rayOutType, COLOR fluxAtten ); /* Init photonScatter[] dispatch table with material specific scattering routines */ void initPhotonScatterFuncs (); - /* Find antimatter geometry declared as photon sensors */ + /* Find antimatter geometry declared as photon sensors */ void getPhotonSensors (char **sensorList); - + #endif diff --git a/source.c b/source.c index 9fd5c5b..345d136 100644 --- a/source.c +++ b/source.c @@ -1,790 +1,796 @@ #ifndef lint static const char RCSid[] = "$Id: source.c,v 2.78 2022/03/30 16:40:03 greg Exp $"; #endif /* * source.c - routines dealing with illumination sources. * * External symbols declared in source.h */ #include "ray.h" #include "otypes.h" #include "otspecial.h" #include "rtotypes.h" #include "source.h" #include "random.h" #include "pmapsrc.h" #include "pmapmat.h" #ifndef MAXSSAMP #define MAXSSAMP 16 /* maximum samples per ray */ #endif /* * Structures used by direct() */ typedef struct { int sno; /* source number */ FVECT dir; /* source direction */ COLOR coef; /* material coefficient */ COLOR val; /* contribution */ } CONTRIB; /* direct contribution */ typedef struct { int sndx; /* source index (to CONTRIB array) */ float brt; /* brightness (for comparison) */ } CNTPTR; /* contribution pointer */ static CONTRIB *srccnt = NULL; /* source contributions in direct() */ static CNTPTR *cntord = NULL; /* source ordering in direct() */ static int maxcntr = 0; /* size of contribution arrays */ static int cntcmp(const void *p1, const void *p2); void marksources(void) /* find and mark source objects */ { int indirect = 0; int i; OBJREC *o, *m; int ns; /* call us only once! */ if (nsources) error(CONSISTENCY, "Multiple calls to marksources!"); /* initialize dispatch table */ initstypes(); /* find direct sources */ for (i = 0; i < nsceneobjs; i++) { - + o = objptr(i); if (!issurface(o->otype) || o->omod == OVOID) continue; /* find material */ m = findmaterial(o); if (m == NULL) continue; if (m->otype == MAT_CLIP) { markclip(m); /* special case for antimatter */ continue; } if (!islight(m->otype)) continue; /* not source modifier */ - + if (m->oargs.nfargs != (m->otype == MAT_GLOW ? 4 : m->otype == MAT_SPOT ? 7 : 3)) objerror(m, USER, "bad # arguments"); if (m->oargs.farg[0] <= FTINY && (m->oargs.farg[1] <= FTINY) & (m->oargs.farg[2] <= FTINY)) continue; /* don't bother */ if (m->otype == MAT_GLOW && o->otype != OBJ_SOURCE && m->oargs.farg[3] <= FTINY) { indirect += (ambounce > 0); continue; /* don't track these */ } if (sfun[o->otype].of == NULL || sfun[o->otype].of->setsrc == NULL) objerror(o, USER, "illegal material"); if ((ns = newsource()) < 0) goto memerr; setsource(&source[ns], o); if (m->otype == MAT_GLOW) { source[ns].sflags |= SPROX; source[ns].sl.prox = m->oargs.farg[3]; if (source[ns].sflags & SDISTANT) { source[ns].sflags |= SSKIP; indirect += (ambounce > 0); } } else if (m->otype == MAT_SPOT) { if (source[ns].sflags & SDISTANT) objerror(o, WARNING, "distant source is a spotlight"); source[ns].sflags |= SSPOT; if ((source[ns].sl.s = makespot(m)) == NULL) goto memerr; if (source[ns].sflags & SFLAT && !checkspot(source[ns].sl.s,source[ns].snorm)) { objerror(o, WARNING, "invalid spotlight direction"); source[ns].sflags |= SSKIP; } } maxcntr += !(source[ns].sflags & SSKIP); } if (!maxcntr) { if (!indirect) error(WARNING, "no light sources found"); return; /* no direct calculation, it seems */ } #if SHADCACHE for (ns = 0; ns < nsources; ns++) /* initialize obstructor cache */ initobscache(ns); #endif /* PMAP: disable virtual sources */ if (!photonMapping) markvirtuals(); /* find and add virtual sources */ - + /* allocate our contribution arrays */ maxcntr += MAXSPART; /* start with this many */ srccnt = (CONTRIB *)malloc(maxcntr*sizeof(CONTRIB)); cntord = (CNTPTR *)malloc(maxcntr*sizeof(CNTPTR)); if ((srccnt != NULL) & (cntord != NULL)) return; memerr: error(SYSTEM, "out of memory in marksources"); } void distantsources(void) /* only mark distant sources */ { int i; OBJREC *o, *m; int ns; /* call us only once! */ if (nsources) error(CONSISTENCY, "Multiple calls to distantsources!"); /* initialize dispatch table */ initstypes(); /* sources needed for sourcehit() */ for (i = 0; i < nsceneobjs; i++) { - + o = objptr(i); if ((o->otype != OBJ_SOURCE) | (o->omod == OVOID)) continue; /* find material */ m = findmaterial(o); if (m == NULL) continue; if (!islight(m->otype)) continue; /* not source modifier */ - + if (m->oargs.nfargs != (m->otype == MAT_GLOW ? 4 : m->otype == MAT_SPOT ? 7 : 3)) objerror(m, USER, "bad # arguments"); if (m->oargs.farg[0] <= FTINY && (m->oargs.farg[1] <= FTINY) & (m->oargs.farg[2] <= FTINY)) continue; /* don't bother */ if (sfun[o->otype].of == NULL || sfun[o->otype].of->setsrc == NULL) objerror(o, USER, "illegal material"); if ((ns = newsource()) < 0) error(SYSTEM, "out of memory in distantsources"); setsource(&source[ns], o); if (m->otype == MAT_GLOW) { source[ns].sflags |= SPROX|SSKIP; source[ns].sl.prox = m->oargs.farg[3]; } else if (m->otype == MAT_SPOT) objerror(o, WARNING, "distant source is a spotlight"); } } void freesources(void) /* free all source structures */ { if (nsources > 0) { #if SHADCACHE while (nsources--) freeobscache(&source[nsources]); #endif free((void *)source); source = NULL; nsources = 0; } markclip(NULL); if (maxcntr <= 0) return; free((void *)srccnt); srccnt = NULL; free((void *)cntord); cntord = NULL; maxcntr = 0; } int srcray( /* send a ray to a source, return domega */ RAY *sr, /* returned source ray */ RAY *r, /* ray which hit object */ SRCINDEX *si /* source sample index */ ) { double d; /* distance to source */ SRCREC *srcp; rayorigin(sr, SHADOW, r, NULL); /* ignore limits */ if (r == NULL) sr->rmax = 0.0; while ((d = nextssamp(sr, si)) != 0.0) { sr->rsrc = si->sn; /* remember source */ srcp = source + si->sn; if (srcp->sflags & SDISTANT) { if (srcp->sflags & SSPOT && spotout(sr, srcp->sl.s)) continue; return(1); /* sample OK */ } /* local source */ /* check proximity */ if (srcp->sflags & SPROX && d > srcp->sl.prox) continue; /* check angle */ if (srcp->sflags & SSPOT) { if (spotout(sr, srcp->sl.s)) continue; /* adjust solid angle */ si->dom *= d*d; d += srcp->sl.s->flen; si->dom /= d*d; } return(1); /* sample OK */ } return(0); /* no more samples */ } void srcvalue( /* punch ray to source and compute value */ RAY *r ) { SRCREC *sp; sp = &source[r->rsrc]; if (sp->sflags & SVIRTUAL) { /* virtual source */ /* check intersection */ if (!(*ofun[sp->so->otype].funp)(sp->so, r)) return; if (!rayshade(r, r->ro->omod)) /* compute contribution */ goto nomat; rayparticipate(r); return; } /* compute intersection */ if (sp->sflags & SDISTANT ? sourcehit(r) : (*ofun[sp->so->otype].funp)(sp->so, r)) { if (sp->sa.success >= 0) sp->sa.success++; if (!rayshade(r, r->ro->omod)) /* compute contribution */ goto nomat; rayparticipate(r); return; } /* we missed our mark! */ if (sp->sa.success < 0) return; /* bitched already */ sp->sa.success -= AIMREQT; if (sp->sa.success >= 0) return; /* leniency */ sprintf(errmsg, "aiming failure for light source \"%s\"", sp->so->oname); error(WARNING, errmsg); /* issue warning */ return; nomat: objerror(r->ro, USER, "material not found"); } static int transillum( /* check if material is transparent illum */ OBJREC *m ) { m = findmaterial(m); if (m == NULL) return(1); if (m->otype != MAT_ILLUM) return(0); return(!m->oargs.nsargs || !strcmp(m->oargs.sarg[0], VOIDID)); } int sourcehit( /* check to see if ray hit distant source */ RAY *r ) { int glowsrc = -1; int transrc = -1; int first, last; int i; if (r->rsrc >= 0) { /* check only one if aimed */ first = last = r->rsrc; } else { /* otherwise check all */ first = 0; last = nsources-1; } for (i = first; i <= last; i++) { if ((source[i].sflags & (SDISTANT|SVIRTUAL)) != SDISTANT) continue; /* * Check to see if ray is within * solid angle of source. */ if (2.*PI*(1. - DOT(source[i].sloc,r->rdir)) > source[i].ss2) continue; /* is it what we aimed for? */ if (i == r->rsrc) { r->ro = source[i].so; break; } /* * If it's a glow or transparent illum, just remember it. */ if (source[i].sflags & SSKIP) { if (glowsrc < 0) glowsrc = i; continue; } if (transillum(source[i].so)) { if (transrc < 0) transrc = i; continue; } r->ro = source[i].so; /* otherwise, use first hit */ break; } /* * Do we need fallback? */ if (r->ro == NULL) { if (transrc >= 0 && r->crtype & (AMBIENT|SPECULAR)) return(0); /* avoid overcounting */ if (glowsrc >= 0) r->ro = source[glowsrc].so; else return(0); /* nothing usable */ } /* * Assign object index */ r->robj = objndx(r->ro); return(1); } static int cntcmp( /* contribution compare (descending) */ const void *p1, const void *p2 ) { const CNTPTR *sc1 = (const CNTPTR *)p1; const CNTPTR *sc2 = (const CNTPTR *)p2; if (sc1->brt > sc2->brt) return(-1); if (sc1->brt < sc2->brt) return(1); return(0); } void direct( /* add direct component */ RAY *r, /* ray that hit surface */ srcdirf_t *f, /* direct component coefficient function */ void *p /* data for f */ ) { int sn; CONTRIB *scp; SRCINDEX si; int nshadcheck, ncnts; int nhits; double prob, ourthresh, hwt; RAY sr; - + /* PMAP: Factor in direct photons (primarily for debugging/validation) */ /* PMAP: Also primary hook in lightflow mode, since this accounts for direct _and_ indirect component */ if (directPhotonMapping #ifdef PMAP_PHOTONFLOW || lightFlowPhotonMapping #endif ) { (*f)(r -> rcol, p, r -> ron, PI); multDirectPmap(r); return; } /* NOTE: srccnt and cntord global so no recursion */ if (maxcntr <= 0) return; /* no direct?! */ initsrcindex(&si); /* potential contributions */ for (sn = 0; srcray(&sr, r, &si); sn++) { if (sn >= maxcntr) { maxcntr = sn + MAXSPART; srccnt = (CONTRIB *)realloc((void *)srccnt, maxcntr*sizeof(CONTRIB)); cntord = (CNTPTR *)realloc((void *)cntord, maxcntr*sizeof(CNTPTR)); if ((srccnt == NULL) | (cntord == NULL)) error(SYSTEM, "out of memory in direct"); } cntord[sn].sndx = sn; scp = srccnt + sn; scp->sno = sr.rsrc; #if SHADCACHE /* check shadow cache */ if (si.np == 1 && srcblocked(&sr)) { cntord[sn].brt = 0.0; continue; } #endif /* compute coefficient */ (*f)(scp->coef, p, sr.rdir, si.dom); cntord[sn].brt = intens(scp->coef); if (cntord[sn].brt <= 0.0) continue; VCOPY(scp->dir, sr.rdir); copycolor(sr.rcoef, scp->coef); /* compute potential */ sr.revf = srcvalue; rayvalue(&sr); multcolor(sr.rcol, sr.rcoef); copycolor(scp->val, sr.rcol); cntord[sn].brt = bright(sr.rcol); } /* sort contributions */ qsort(cntord, sn, sizeof(CNTPTR), cntcmp); { /* find last */ int l, m; ncnts = l = sn; sn = 0; while ((m = (sn + ncnts) >> 1) != l) { if (cntord[m].brt > 0.0) sn = m; else ncnts = m; l = m; } } if (ncnts == 0) return; /* no contributions! */ /* accumulate tail */ for (sn = ncnts-1; sn > 0; sn--) cntord[sn-1].brt += cntord[sn].brt; /* compute number to check */ nshadcheck = pow((double)ncnts, shadcert) + .5; /* modify threshold */ ourthresh = shadthresh / r->rweight; /* test for shadows */ for (nhits = 0, hwt = 0.0, sn = 0; sn < ncnts; hwt += (double)source[scp->sno].nhits / (double)source[scp->sno].ntests, sn++) { /* check threshold */ if (sn >= MINSHADCNT && (sn+nshadcheck>=ncnts ? cntord[sn].brt : cntord[sn].brt-cntord[sn+nshadcheck].brt) < ourthresh*bright(r->rcol)) break; scp = srccnt + cntord[sn].sndx; /* test for hit */ rayorigin(&sr, SHADOW, r, NULL); copycolor(sr.rcoef, scp->coef); VCOPY(sr.rdir, scp->dir); sr.rsrc = scp->sno; /* keep statistics */ if (source[scp->sno].ntests++ > 0xfffffff0) { source[scp->sno].ntests >>= 1; source[scp->sno].nhits >>= 1; } if (localhit(&sr, &thescene) && ( sr.ro != source[scp->sno].so || source[scp->sno].sflags & SFOLLOW )) { /* follow entire path */ raycont(&sr); if (trace != NULL) (*trace)(&sr); /* trace execution */ if (bright(sr.rcol) <= FTINY) { #if SHADCACHE if ((scp <= srccnt || scp[-1].sno != scp->sno) && (scp >= srccnt+ncnts-1 || scp[1].sno != scp->sno)) srcblocker(&sr); #endif continue; /* missed! */ } rayparticipate(&sr); multcolor(sr.rcol, sr.rcoef); copycolor(scp->val, sr.rcol); } else if (trace != NULL && (source[scp->sno].sflags & (SDISTANT|SVIRTUAL|SFOLLOW)) == (SDISTANT|SFOLLOW) && sourcehit(&sr) && rayshade(&sr, sr.ro->omod)) { (*trace)(&sr); /* trace execution */ /* skip call to rayparticipate() & scp->val update */ } /* add contribution if hit */ addcolor(r->rcol, scp->val); nhits++; source[scp->sno].nhits++; } /* source hit rate */ if (hwt > FTINY) hwt = (double)nhits / hwt; else hwt = 0.5; #ifdef DEBUG sprintf(errmsg, "%d tested, %d untested, %f conditional hit rate\n", sn, ncnts-sn, hwt); eputs(errmsg); #endif /* add in untested sources */ for ( ; sn < ncnts; sn++) { scp = srccnt + cntord[sn].sndx; prob = hwt * (double)source[scp->sno].nhits / (double)source[scp->sno].ntests; if (prob < 1.0) scalecolor(scp->val, prob); addcolor(r->rcol, scp->val); } } void srcscatter( /* compute source scattering into ray */ RAY *r ) { int oldsampndx; int nsamps; RAY sr; SRCINDEX si; double t, d; double re, ge, be; COLOR cvext; int i, j; if (r->rot >= FHUGE*.99 || r->gecc >= 1.-FTINY) return; /* this can never work */ /* PMAP: do unconditional inscattering for volume photons */ if (!volumePhotonMapping && (r->slights == NULL || r->slights[0] == 0)) return; - + if (ssampdist <= FTINY || (nsamps = r->rot/ssampdist + .5) < 1) nsamps = 1; #if MAXSSAMP else if (nsamps > MAXSSAMP) nsamps = MAXSSAMP; #endif oldsampndx = samplendx; samplendx = random()&0x7fff; /* randomize */ for (i = volumePhotonMapping ? 1 : r->slights[0]; i > 0; i--) { /* for each source OR once if volume photon map enabled */ for (j = 0; j < nsamps; j++) { /* for each sample position */ samplendx++; t = r->rot * (j+frandom())/nsamps; /* extinction */ re = t*colval(r->cext,RED); ge = t*colval(r->cext,GRN); be = t*colval(r->cext,BLU); setcolor(cvext, re > 92. ? 0. : exp(-re), ge > 92. ? 0. : exp(-ge), be > 92. ? 0. : exp(-be)); if (intens(cvext) <= FTINY) break; /* too far away */ sr.rorg[0] = r->rorg[0] + r->rdir[0]*t; sr.rorg[1] = r->rorg[1] + r->rdir[1]*t; sr.rorg[2] = r->rorg[2] + r->rdir[2]*t; - + if (!volumePhotonMapping) { initsrcindex(&si); /* sample ray to this source */ si.sn = r->slights[i]; nopart(&si, &sr); if (!srcray(&sr, NULL, &si) || sr.rsrc != r->slights[i]) continue; /* no path */ #if SHADCACHE if (srcblocked(&sr)) /* check shadow cache */ continue; #endif copycolor(sr.cext, r->cext); copycolor(sr.albedo, r->albedo); sr.gecc = r->gecc; sr.slights = r->slights; rayvalue(&sr); /* eval. source ray */ if (bright(sr.rcol) <= FTINY) { #if SHADCACHE srcblocker(&sr); /* add blocker to cache */ #endif continue; } if (r->gecc <= FTINY) /* compute P(theta) */ d = 1.; else { d = DOT(r->rdir, sr.rdir); d = 1. + r->gecc*r->gecc - 2.*r->gecc*d; d = (1. - r->gecc*r->gecc) / (d*sqrt(d)); } /* other factors */ d *= si.dom * r->rot / (4.*PI*nsamps); scalecolor(sr.rcol, d); } else { /* PMAP: Add ambient inscattering from - * volume photons; note we reverse the + * volume photons; note we reverse the * incident ray direction since we're * now in *backward* raytracing mode! */ sr.rdir [0] = -r -> rdir [0]; sr.rdir [1] = -r -> rdir [1]; sr.rdir [2] = -r -> rdir [2]; sr.gecc = r -> gecc; inscatterVolumePmap(&sr, sr.rcol); scalecolor(sr.rcol, r -> rot / nsamps); } multcolor(sr.rcol, r->cext); multcolor(sr.rcol, r->albedo); multcolor(sr.rcol, cvext); addcolor(r->rcol, sr.rcol); /* add it in */ } } samplendx = oldsampndx; } /**************************************************************** * The following macros were separated from the m_light() routine * because they are very nasty and difficult to understand. */ /* illumblock * * * We cannot allow an illum to pass to another illum, because that * would almost certainly constitute overcounting. * However, we do allow an illum to pass to another illum * that is actually going to relay to a virtual light source. * We also prevent an illum from passing to a glow; this provides a * convenient mechanism for defining detailed light source * geometry behind (or inside) an effective radiator. */ -static int -weaksrcmat(OBJREC *m) /* identify material */ +int +weaksrcmat(OBJECT obj) /* identify material */ { - m = findmaterial(m); + OBJREC *m = findmaterial(objptr(obj)); + if (m == NULL) return(0); return((m->otype==MAT_ILLUM) | (m->otype==MAT_GLOW)); } #define illumblock(m, r) (!(source[r->rsrc].sflags&SVIRTUAL) && \ r->rod > 0.0 && \ - weaksrcmat(source[r->rsrc].so)) + weaksrcmat(source[r->rsrc].so->omod)) /* wrongsource * * * This source is the wrong source (ie. overcounted) if we are * aimed to a different source than the one we hit and the one * we hit is not an illum that should be passed. */ #define wrongsource(m, r) (r->rsrc>=0 && source[r->rsrc].so!=r->ro && \ (m->otype!=MAT_ILLUM || illumblock(m,r))) /* distglow * * * A distant glow is an object that sometimes acts as a light source, * but is too far away from the test point to be one in this case. * (Glows with negative radii should NEVER participate in illumination.) */ #define distglow(m, r, d) (m->otype==MAT_GLOW && \ m->oargs.farg[3] >= -FTINY && \ d > m->oargs.farg[3]) /* badcomponent * * * We must avoid counting light sources in the ambient calculation, * since the direct component is handled separately. Therefore, any * ambient ray which hits an active light source must be discarded. * The same is true for stray specular samples, since the specular * contribution from light sources is calculated separately. */ /* PMAP: Also avoid counting sources via transferred ambient rays (e.g. * through glass) when photon mapping is enabled, as these indirect - * components are already accounted for. + * components are already accounted for. */ #define badcomponent(m, r) (srcRayInPmap(r) || \ (r->crtype&(AMBIENT|SPECULAR) && \ !(r->crtype&SHADOW || r->rod < 0.0 || \ /* not 100% correct */ distglow(m, r, r->rot)))) /* passillum * * * An illum passes to another material type when we didn't hit it * on purpose (as part of a direct calculation), or it is relaying * a virtual light source. */ #define passillum(m, r) (m->otype==MAT_ILLUM && \ (r->rsrc<0 || source[r->rsrc].so!=r->ro || \ source[r->rsrc].sflags&SVIRTUAL)) /* srcignore * * * The -dv flag is normally on for sources to be visible. */ #define srcignore(m, r) !(directvis || r->crtype&SHADOW || \ distglow(m, r, raydist(r,PRIMARY))) int m_light( /* ray hit a light source */ OBJREC *m, RAY *r ) { /* check for over-counting */ if (badcomponent(m, r)) { setcolor(r->rcoef, 0.0, 0.0, 0.0); return(1); } if (wrongsource(m, r)) { setcolor(r->rcoef, 0.0, 0.0, 0.0); return(1); + } + /* check for glow through BSDF with pmap */ + if (glowignore(m, r)) { + setcolor(r->rcoef, 0.0, 0.0, 0.0); + return(1); } /* check for passed illum */ if (passillum(m, r)) { if (m->oargs.nsargs && strcmp(m->oargs.sarg[0], VOIDID)) return(rayshade(r,lastmod(objndx(m),m->oargs.sarg[0]))); raytrans(r); return(1); } /* check for invisibility */ if (srcignore(m, r)) { setcolor(r->rcoef, 0.0, 0.0, 0.0); return(1); } /* otherwise treat as source */ /* check for behind */ if (r->rod < 0.0) return(1); /* check for outside spot */ if (m->otype==MAT_SPOT && spotout(r, makespot(m))) return(1); /* get distribution pattern */ raytexture(r, m->omod); /* get source color */ setcolor(r->rcol, m->oargs.farg[0], m->oargs.farg[1], m->oargs.farg[2]); /* modify value */ multcolor(r->rcol, r->pcol); return(1); } diff --git a/source.h b/source.h index 3e7b025..f4a6a1b 100644 --- a/source.h +++ b/source.h @@ -1,209 +1,209 @@ /* RCSid $Id: source.h,v 2.22 2020/12/17 03:30:37 greg Exp $ */ /* * source.h - header file for ray tracing sources. * * Include after ray.h */ #ifndef _RAD_SOURCE_H_ #define _RAD_SOURCE_H_ #include #ifdef __cplusplus extern "C" { #endif #ifndef AIMREQT #define AIMREQT 100 /* required aim success/failure */ #endif #ifndef SHADCACHE #define SHADCACHE 20 /* shadow cache resolution */ #endif #ifndef MINSHADCNT #define MINSHADCNT 2 /* test at least this many shadows */ #endif #define SDISTANT 01 /* source distant flag */ #define SSKIP 02 /* source skip flag */ #define SPROX 04 /* source proximity flag */ #define SSPOT 010 /* source spotlight flag */ #define SVIRTUAL 020 /* source virtual flag */ #define SFLAT 040 /* source flat flag */ #define SCIR 0100 /* source circular flag */ #define SCYL 0200 /* source cylindrical flag */ #define SFOLLOW 0400 /* source follow path flag */ typedef struct { FVECT aim; /* aim direction or center */ float siz; /* output solid angle or area */ float flen; /* focal length (negative if distant source) */ } SPOT; /* spotlight */ typedef struct { union { struct { FVECT u, v; /* unit vectors */ } f; /* flat source indexing */ struct { FVECT o; /* origin position */ RREAL e1, e2; /* 1/extent */ int ax; /* major direction */ } d; /* distant source indexing */ } p; /* indexing parameters */ OBJECT obs[1]; /* cache obstructors (extends struct) */ } OBSCACHE; /* obstructor cache */ typedef struct { FVECT sloc; /* direction or position of source */ FVECT ss[3]; /* source dimension vectors, U, V, and W */ float srad; /* maximum source radius */ float ss2; /* solid angle or projected area */ OBJREC *so; /* source destination object */ struct { float prox; /* proximity */ SPOT *s; /* spot */ } sl; /* localized source information */ union { long success; /* successes - AIMREQT*failures */ struct { short pn; /* projection number */ int sn; /* next source to aim for */ } sv; /* virtual source */ } sa; /* source aiming information */ unsigned long ntests, nhits; /* shadow tests and hits */ #ifdef SHADCACHE OBSCACHE *obscache; /* obstructor cache */ #endif int sflags; /* source flags */ } SRCREC; /* light source */ #define MAXSPART 64 /* maximum partitions per source */ #define SU 0 /* U vector or partition */ #define SV 1 /* V vector or partition */ #define SW 2 /* W vector or partition */ #define S0 3 /* leaf partition */ #define snorm ss[SW] /* normal vector for flat source */ typedef struct { double dom; /* solid angle of partition */ int sn; /* source number */ short np; /* number of partitions */ short sp; /* this partition number */ unsigned char spt[MAXSPART/2]; /* source partitioning */ } SRCINDEX; /* source index structure */ #define initsrcindex(s) ((s)->sn = (s)->sp = -1, (s)->np = 0) #define clrpart(pt) memset((char *)(pt), '\0', MAXSPART/2) #define setpart(pt,i,v) ((pt)[(i)>>2] |= (v)<<(((i)&3)<<1)) #define spart(pt,pi) ((pt)[(pi)>>2] >> (((pi)&3)<<1) & 3) /* * Special support functions for sources */ /* * Virtual source materials must define the following. * * vproj(pm, op, sp, i) Compute i'th virtual projection * of source sp in object op and assign * the 4x4 transformation matrix pm. * Return 1 on success, 0 if no i'th projection. * * nproj The number of projections. The value of * i passed to vproj runs from 0 to nproj-1. */ typedef struct { int (*vproj)(); /* project virtual sources */ int nproj; /* number of possible projections */ } VSMATERIAL; /* virtual source material functions */ typedef struct { void (*setsrc)(); /* set light source for object */ void (*partit)(); /* partition light source object */ double (*getpleq)(); /* plane equation for surface */ double (*getdisk)(); /* maximum disk for surface */ } SOBJECT; /* source object functions */ typedef union { VSMATERIAL *mf; /* material functions */ SOBJECT *of; /* object functions */ } SRCFUNC; /* source functions */ extern SRCFUNC sfun[]; /* source dispatch table */ extern SRCREC *source; /* our source list */ extern int nsources; /* the number of sources */ #define sflatform(sn,dir) -DOT(source[sn].snorm, dir) #define getplaneq(c,o) (*sfun[(o)->otype].of->getpleq)(c,o) #define getmaxdisk(c,o) (*sfun[(o)->otype].of->getdisk)(c,o) #define setsource(s,o) (*sfun[(o)->otype].of->setsrc)(s,o) /* defined in source.c */ extern void marksources(void); extern void distantsources(void); extern void freesources(void); extern int srcray(RAY *sr, RAY *r, SRCINDEX *si); extern void srcvalue(RAY *r); extern int sourcehit(RAY *r); typedef void srcdirf_t(COLOR cv, void *np, FVECT ldir, double omega); extern void direct(RAY *r, srcdirf_t *f, void *p); extern void srcscatter(RAY *r); extern int m_light(OBJREC *m, RAY *r); +extern int weaksrcmat(OBJECT obj); /* defined in srcobstr.c */ extern void initobscache(int sn); extern int srcblocker(RAY *r); extern int srcblocked(RAY *r); extern void freeobscache(SRCREC *s); extern void markclip(OBJREC *m); /* defined in srcsamp.c */ extern double nextssamp(RAY *r, SRCINDEX *si); extern int skipparts(int ct[3], int sz[3], int pp[2], unsigned char *pt); extern void nopart(SRCINDEX *si, RAY *r); extern void cylpart(SRCINDEX *si, RAY *r); extern void flatpart(SRCINDEX *si, RAY *r); extern double scylform(int sn, FVECT dir); /* defined in srcsupp.c */ extern void initstypes(void); extern int newsource(void); extern void setflatss(SRCREC *src); extern void fsetsrc(SRCREC *src, OBJREC *so); extern void ssetsrc(SRCREC *src, OBJREC *so); extern void sphsetsrc(SRCREC *src, OBJREC *so); extern void rsetsrc(SRCREC *src, OBJREC *so); extern void cylsetsrc(SRCREC *src, OBJREC *so); extern SPOT *makespot(OBJREC *m); extern int spotout(RAY *r, SPOT *s); extern double fgetmaxdisk(FVECT ocent, OBJREC *op); extern double rgetmaxdisk(FVECT ocent, OBJREC *op); extern double fgetplaneq(FVECT nvec, OBJREC *op); extern double rgetplaneq(FVECT nvec, OBJREC *op); extern int commonspot(SPOT *sp1, SPOT *sp2, FVECT org); extern int commonbeam(SPOT *sp1, SPOT *sp2, FVECT org); extern int checkspot(SPOT *sp, FVECT nrm); extern double spotdisk(FVECT oc, OBJREC *op, SPOT *sp, FVECT pos); extern double beamdisk(FVECT oc, OBJREC *op, SPOT *sp, FVECT dir); extern double intercircle(FVECT cc, FVECT c1, FVECT c2, double r1s, double r2s); /* defined in virtuals.c */ extern void markvirtuals(void); extern void addvirtuals(int sn, int nr); extern void vproject(OBJREC *o, int sn, int n); extern OBJREC *vsmaterial(OBJREC *o); extern int makevsrc(OBJREC *op, int sn, MAT4 pm); extern double getdisk(FVECT oc, OBJREC *op, int sn); extern int vstestvis(int f, OBJREC *o, FVECT oc, double or2, int sn); extern void virtverb(int sn, FILE *fp); #ifdef __cplusplus } #endif #endif /* _RAD_SOURCE_H_ */ -