--- ray/src/rt/pmapmat.c 2015/02/24 19:39:27 2.1 +++ ray/src/rt/pmapmat.c 2015/05/21 05:54:54 2.5 @@ -4,10 +4,11 @@ Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) (c) Fraunhofer Institute for Solar Energy Systems, - Lucerne University of Applied Sciences & Arts + (c) Lucerne University of Applied Sciences and Arts, + supported by the Swiss National Science Foundation (SNSF, #147053) ================================================================== - $Id: pmapmat.c,v 2.1 2015/02/24 19:39:27 greg Exp $ + $Id: pmapmat.c,v 2.5 2015/05/21 05:54:54 greg Exp $ */ @@ -113,7 +114,7 @@ void photonRay (const RAY *rayIn, RAY *rayOut, VCOPY(rayOut -> rdir, rayIn -> rdir); } else if (fluxAtten) { - /* Attenuate and normalised flux for scattered rays */ + /* Attenuate and normalise flux for scattered rays */ multcolor(rayOut -> rcol, fluxAtten); colorNorm(rayOut -> rcol); } @@ -198,12 +199,7 @@ static int isoSpecPhotonScatter (NORMDAT *nd, RAY *ray int niter, i = 0; /* Set up sample coordinates */ - do { - v [0] = v [1] = v [2] = 0; - v [i++] = 1; - fcross(u, v, nd -> pnorm); - } while (normalize(u) < FTINY); - + getperpendicular(u, nd -> pnorm, 1); fcross(v, nd -> pnorm, u); if (nd -> specfl & SP_REFL) { @@ -267,12 +263,7 @@ static void diffPhotonScatter (FVECT normal, RAY* rayO int i = 0; /* Set up sample coordinates */ - do { - v [0] = v [1] = v [2] = 0; - v [i++] = 1; - fcross(u, v, normal); - } while (normalize(u) < FTINY); - + getperpendicular(u, normal, 1); fcross(v, normal, u); /* Convert theta & phi to cartesian */ @@ -323,7 +314,7 @@ static int normalPhotonScatter (OBJREC *mat, RAY *rayI nd.specfl |= SP_FLAT; /* Perturb normal */ - if ((hastexture = DOT(rayIn -> pert, rayIn -> pert)) > sqr(FTINY)) + if ((hastexture = (DOT(rayIn -> pert, rayIn -> pert) > sqr(FTINY)) )) nd.pdot = raynormal(nd.pnorm, rayIn); else { VCOPY(nd.pnorm, rayIn -> ron); @@ -752,7 +743,7 @@ static int dielectricPhotonScatter (OBJREC *mat, RAY * /* get modifiers */ raytexture(rayIn, mat -> omod); - if ((hastexture = DOT(rayIn -> pert, rayIn -> pert)) > FTINY * FTINY) + if ((hastexture = (DOT(rayIn -> pert, rayIn -> pert) > FTINY * FTINY))) /* Perturb normal */ cos1 = raynormal(dnorm, rayIn); else { @@ -899,7 +890,7 @@ static int glassPhotonScatter (OBJREC *mat, RAY *rayIn /* reorient if necessary */ if (rayIn -> rod < 0) flipsurface(rayIn); - if ((hastexture = DOT(rayIn -> pert, rayIn -> pert)) > FTINY * FTINY) + if ((hastexture = (DOT(rayIn -> pert, rayIn -> pert) > FTINY * FTINY) )) pdot = raynormal(pnorm, rayIn); else { VCOPY(pnorm, rayIn -> ron); @@ -1383,16 +1374,303 @@ static int pattexPhotonScatter (OBJREC *mat, RAY *rayI +#if 0 + static int bsdfPhotonScatter (OBJREC *mat, RAY *rayIn) + /* Generate new photon ray for BSDF modifier and recurse. */ + { + int hitFront; + SDError err; + FVECT upvec; + MFUNC *mf; + BSDFDAT nd; + RAY rayOut; + + /* Following code adapted from m_bsdf() */ + /* Check arguments */ + if (mat -> oargs.nsargs < 6 || mat -> oargs.nfargs > 9 || + mat -> oargs.nfargs % 3) + objerror(mat, USER, "bad # arguments"); + + hitFront = (rayIn -> rod > 0); + + /* Load cal file */ + mf = getfunc(mat, 5, 0x1d, 1); + + /* Get thickness */ + nd.thick = evalue(mf -> ep [0]); + if ((-FTINY <= nd.thick) & (nd.thick <= FTINY)) + nd.thick = .0; + + if (nd.thick != .0 || (!hitFront && !backvis)) { + /* Proxy geometry present, so use it instead and transfer ray */ + photonRay(rayIn, &rayOut, PMAP_XFER, NULL); + tracePhoton(&rayOut); + + return 0; + } + + /* Get BSDF data */ + nd.sd = loadBSDF(mat -> oargs.sarg [1]); + + /* Diffuse reflectance */ + if (hitFront) { + if (mat -> oargs.nfargs < 3) + setcolor(nd.rdiff, .0, .0, .0); + else setcolor(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, mat -> oargs.farg [3], mat -> oargs.farg [4], + mat -> oargs.farg [5]); + + /* Diffuse transmittance */ + if (mat -> oargs.nfargs < 9) + setcolor(nd.tdiff, .0, .0, .0); + else setcolor(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 [1]); + upvec [1] = evalue(mf -> ep [2]); + upvec [2] = evalue(mf -> ep [3]); + + 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]; + } + + /* Following code adapted from SDsampBSDF() */ + { + SDSpectralDF *rdf, *tdf; + SDValue bsdfVal; + double xi, rhoDiff = 0; + float coef [SDmaxCh]; + int i, j, n, nr; + SDComponent *sdc; + const SDCDst **cdarr = NULL; + + /* Get diffuse albedo (?) */ + if (hitFront) { + bsdfVal = nd.sd -> rLambFront; + rdf = nd.sd -> rf; + tdf = nd.sd -> tf ? nd.sd -> tf : nd.sd -> tb; + } + else { + bsdfVal = nd.sd -> rLambBack; + rdf = nd.sd -> rb; + tdf = nd.sd -> tb ? nd.sd -> tb : nd.sd -> tf; + } + + rhoDiff = bsdfVal.cieY; + bsdfVal.cieY += nd.sd -> tLamb.cieY; + + /* Allocate non-diffuse sampling */ + i = nr = rdf ? rdf -> ncomp : 0; + j = tdf ? tdf -> ncomp : 0; + n = i + j; + + if (n > 0 && !(cdarr = (const SDCDst**)malloc(n * sizeof(SDCDst*)))) + objerror(mat, USER, transSDError(SDEmemory)); + + while (j-- > 0) { + /* Sum up non-diffuse transmittance */ + cdarr [i + j] = (*tdf -> comp [j].func -> getCDist)(nd.vray, &tdf -> comp [j]); + + if (!cdarr [i + j]) + cdarr [i + j] = &SDemptyCD; + else bsdfVal.cieY += cdarr [i + j] -> cTotal; + } + + while (i-- > 0) { + /* Sum up non-diffuse reflectance */ + cdarr [i] = (*rdf -> comp [i].func -> getCDist)(nd.vray, &rdf -> comp [i]); + + if (!cdarr [i]) + cdarr [i] = &SDemptyCD; + else bsdfVal.cieY += cdarr [i] -> cTotal; + } + + if (bsdfVal.cieY <= FTINY) { + /* Don't bother sampling, just absorb photon */ + if (cdarr) + free(cdarr); + return 0; + } + + /* Insert direct and indirect photon hits if diffuse component */ + if (rhoDiff > FTINY || nd.sd -> tLamb.cieY > FTINY) + addPhotons(rayIn); + + xi = pmapRandom(rouletteState); + + if ((xi -= rhoDiff) <= 0) { + /* Diffuse reflection */ + photonRay(rayIn, &rayOut, PMAP_DIFFREFL, nd.rdiff); + diffPhotonScatter(nd.pnorm, &rayOut); + } + else if ((xi -= nd.sd -> tLamb.cieY) <= 0) { + /* Diffuse transmission */ + flipsurface(rayIn); + photonRay(rayIn, &rayOut, PMAP_DIFFTRANS, nd.tdiff); + bsdfVal.spec = nd.sd -> tLamb.spec; + diffPhotonScatter(nd.pnorm, &rayOut); + } + else { + int rayOutType; + COLOR bsdfRGB; + + /* Non-diffuse CDF inversion (?) */ + for (i = 0; i < n && (xi -= cdarr [i] -> cTotal) > 0; i++); + + if (i >= n) { + /* Absorbed -- photon went Deer Hunter */ + if (cdarr) + free(cdarr); + return 0; + } + + if (i < nr) { + /* Non-diffuse reflection */ + sdc = &rdf -> comp [i]; + rayOutType = PMAP_SPECREFL; + } + else { + /* Non-diffuse transmission */ + sdc = &tdf -> comp [i - nr]; + rayOutType = PMAP_SPECTRANS; + } + + /* Generate non-diff sample dir */ + VCOPY(rayOut.rdir, nd.vray); + err = (*sdc -> func -> sampCDist) + (rayOut.rdir, pmapRandom(scatterState), cdarr [i]); + if (err) + objerror(mat, USER, transSDError(SDEinternal)); + + /* Get colour */ + j = (*sdc -> func -> getBSDFs)(coef, rayOut.rdir, nd.vray, sdc); + + if (j <= 0) { + sprintf(SDerrorDetail, "BSDF \"%s\" sampling value error", + nd.sd -> name); + objerror(mat, USER, transSDError(SDEinternal)); + } + + bsdfVal.spec = sdc -> cspec [0]; + rhoDiff = coef [0]; + + while (--j) { + c_cmix(&bsdfVal.spec, rhoDiff, &bsdfVal.spec, coef [j], + &sdc -> cspec [j]); + rhoDiff += coef [j]; + } + + /* ? */ + c_ccvt(&bsdfVal.spec, C_CSXY + C_CSSPEC); + ccy2rgb(&bsdfVal.spec, bsdfVal.cieY, bsdfRGB); + + /* Xform outgoing dir to world coords */ + if ((err = SDmapDir(rayOut.rdir, nd.fromloc, rayOut.rdir))) { + objerror(mat, USER, transSDError(err)); + return 0; + } + + photonRay(rayIn, &rayOut, rayOutType, bsdfRGB); + } + + if (cdarr) + free(cdarr); + } + + /* Clean up BSDF */ + SDfreeCache(nd.sd); + + tracePhoton(&rayOut); + return 0; + } +#else + +/* + The following code is + (c) Lucerne University of Applied Sciences and Arts, + supported by the Swiss National Science Foundation (SNSF, #147053) +*/ + static int bsdfPhotonScatter (OBJREC *mat, RAY *rayIn) /* Generate new photon ray for BSDF modifier and recurse. */ { int hitFront; SDError err; + SDValue bsdfVal; FVECT upvec; MFUNC *mf; BSDFDAT nd; RAY rayOut; - + COLOR bsdfRGB; + double prDiff, ptDiff, prDiffSD, ptDiffSD, prSpecSD, ptSpecSD, + albedo, xi, xi2; + const double patAlb = colorAvg(rayIn -> pcol); + /* Following code adapted from m_bsdf() */ /* Check arguments */ if (mat -> oargs.nsargs < 6 || mat -> oargs.nfargs > 9 || @@ -1420,7 +1698,7 @@ static int bsdfPhotonScatter (OBJREC *mat, RAY *rayIn) /* Get BSDF data */ nd.sd = loadBSDF(mat -> oargs.sarg [1]); - /* Diffuse reflectance */ + /* Extra diffuse reflectance from material def */ if (hitFront) { if (mat -> oargs.nfargs < 3) setcolor(nd.rdiff, .0, .0, .0); @@ -1439,7 +1717,7 @@ static int bsdfPhotonScatter (OBJREC *mat, RAY *rayIn) else setcolor(nd.rdiff, mat -> oargs.farg [3], mat -> oargs.farg [4], mat -> oargs.farg [5]); - /* Diffuse transmittance */ + /* Extra diffuse transmittance from material def */ if (mat -> oargs.nfargs < 9) setcolor(nd.tdiff, .0, .0, .0); else setcolor(nd.tdiff, mat -> oargs.farg [6], mat -> oargs.farg [7], @@ -1506,156 +1784,112 @@ static int bsdfPhotonScatter (OBJREC *mat, RAY *rayIn) nd.pnorm [1] = -nd.pnorm [1]; nd.pnorm [2] = -nd.pnorm [2]; } - - /* Following code adapted from SDsampBSDF() */ - { - SDSpectralDF *rdf, *tdf; - SDValue bsdfVal; - double xi, rhoDiff = 0; - float coef [SDmaxCh]; - int i, j, n, nr; - SDComponent *sdc; - const SDCDst **cdarr = NULL; - - /* Get diffuse albedo (?) */ - if (hitFront) { - bsdfVal = nd.sd -> rLambFront; - rdf = nd.sd -> rf; - tdf = nd.sd -> tf ? nd.sd -> tf : nd.sd -> tb; - } - else { - bsdfVal = nd.sd -> rLambBack; - rdf = nd.sd -> rb; - tdf = nd.sd -> tb ? nd.sd -> tb : nd.sd -> tf; - } - - rhoDiff = bsdfVal.cieY; - bsdfVal.cieY += nd.sd -> tLamb.cieY; - - /* Allocate non-diffuse sampling */ - i = nr = rdf ? rdf -> ncomp : 0; - j = tdf ? tdf -> ncomp : 0; - n = i + j; - - if (n > 0 && !(cdarr = (const SDCDst**)malloc(n * sizeof(SDCDst*)))) - objerror(mat, USER, transSDError(SDEmemory)); + + /* 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 + * 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"); + */ - while (j-- > 0) { - /* Sum up non-diffuse transmittance */ - cdarr [i + j] = (*tdf -> comp [j].func -> getCDist)(nd.vray, &tdf -> comp [j]); - - if (!cdarr [i + j]) - cdarr [i + j] = &SDemptyCD; - else bsdfVal.cieY += cdarr [i + j] -> cTotal; - } + /* Insert direct and indirect photon hits if diffuse component */ + if (prDiff + ptDiff + prDiffSD + ptDiffSD > FTINY) + addPhotons(rayIn); + + xi = xi2 = pmapRandom(rouletteState); - while (i-- > 0) { - /* Sum up non-diffuse reflectance */ - cdarr [i] = (*rdf -> comp [i].func -> getCDist)(nd.vray, &rdf -> comp [i]); + if (xi > albedo) + /* Absorbtion */ + return 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) */ + flipsurface(rayIn); + photonRay(rayIn, &rayOut, PMAP_DIFFTRANS, nd.tdiff); + diffPhotonScatter(nd.pnorm, &rayOut); + } + + else { /* Sample SDF */ + if ((xi -= prDiffSD) <= 0) { + /* Diffuse SDF reflection (constant component) */ + if ((err = SDsampBSDF(&bsdfVal, nd.vray, xi2, + SDsampDf | SDsampR, nd.sd))) + objerror(mat, USER, transSDError(err)); - if (!cdarr [i]) - cdarr [i] = &SDemptyCD; - else bsdfVal.cieY += cdarr [i] -> cTotal; + /* Apply pattern to spectral component */ + ccy2rgb(&bsdfVal.spec, bsdfVal.cieY, bsdfRGB); + multcolor(bsdfRGB, rayIn -> pcol); + photonRay(rayIn, &rayOut, PMAP_DIFFREFL, bsdfRGB); } - - if (bsdfVal.cieY <= FTINY) { - /* Don't bother sampling, just absorb photon */ - if (cdarr) - free(cdarr); - return 0; - } - - /* Insert direct and indirect photon hits if diffuse component */ - if (rhoDiff > FTINY || nd.sd -> tLamb.cieY > FTINY) - addPhotons(rayIn); - - xi = pmapRandom(rouletteState); - - if ((xi -= rhoDiff) <= 0) { - /* Diffuse reflection */ - photonRay(rayIn, &rayOut, PMAP_DIFFREFL, nd.rdiff); - diffPhotonScatter(nd.pnorm, &rayOut); - } - else if ((xi -= nd.sd -> tLamb.cieY) <= 0) { - /* Diffuse transmission */ - flipsurface(rayIn); - photonRay(rayIn, &rayOut, PMAP_DIFFTRANS, nd.tdiff); - bsdfVal.spec = nd.sd -> tLamb.spec; - diffPhotonScatter(nd.pnorm, &rayOut); - } - else { - int rayOutType; - COLOR bsdfRGB; - - /* Non-diffuse CDF inversion (?) */ - for (i = 0; i < n && (xi -= cdarr [i] -> cTotal) > 0; i++); - - if (i >= n) { - /* Absorbed -- photon went Deer Hunter */ - if (cdarr) - free(cdarr); - return 0; - } - if (i < nr) { - /* Non-diffuse reflection */ - sdc = &rdf -> comp [i]; - rayOutType = PMAP_SPECREFL; - } - else { - /* Non-diffuse transmission */ - sdc = &tdf -> comp [i - nr]; - rayOutType = PMAP_SPECTRANS; - } + else if ((xi -= ptDiffSD) <= 0) { + /* Diffuse SDF transmission (constant component) */ + if ((err = SDsampBSDF(&bsdfVal, nd.vray, xi2, + SDsampDf | SDsampT, nd.sd))) + objerror(mat, USER, transSDError(err)); - /* Generate non-diff sample dir */ - VCOPY(rayOut.rdir, nd.vray); - err = (*sdc -> func -> sampCDist) - (rayOut.rdir, pmapRandom(scatterState), cdarr [i]); - if (err) - objerror(mat, USER, transSDError(SDEinternal)); - - /* Get colour */ - j = (*sdc -> func -> getBSDFs)(coef, rayOut.rdir, nd.vray, sdc); - - if (j <= 0) { - sprintf(SDerrorDetail, "BSDF \"%s\" sampling value error", - nd.sd -> name); - objerror(mat, USER, transSDError(SDEinternal)); - } - - bsdfVal.spec = sdc -> cspec [0]; - rhoDiff = coef [0]; - - while (--j) { - c_cmix(&bsdfVal.spec, rhoDiff, &bsdfVal.spec, coef [j], - &sdc -> cspec [j]); - rhoDiff += coef [j]; - } - - /* ? */ - c_ccvt(&bsdfVal.spec, C_CSXY + C_CSSPEC); + /* Apply pattern to spectral component */ ccy2rgb(&bsdfVal.spec, bsdfVal.cieY, bsdfRGB); - - /* Xform outgoing dir to world coords */ - if ((err = SDmapDir(rayOut.rdir, nd.fromloc, rayOut.rdir))) { + multcolor(bsdfRGB, rayIn -> pcol); + addcolor(bsdfRGB, nd.tdiff); + flipsurface(rayIn); /* Necessary? */ + photonRay(rayIn, &rayOut, PMAP_DIFFTRANS, bsdfRGB); + } + + else if ((xi -= prSpecSD) <= 0) { + /* Non-diffuse ("specular") SDF reflection */ + if ((err = SDsampBSDF(&bsdfVal, nd.vray, xi2, + SDsampSp | SDsampR, nd.sd))) objerror(mat, USER, transSDError(err)); - return 0; - } - photonRay(rayIn, &rayOut, rayOutType, bsdfRGB); + ccy2rgb(&bsdfVal.spec, bsdfVal.cieY, bsdfRGB); + photonRay(rayIn, &rayOut, PMAP_SPECREFL, bsdfRGB); } - if (cdarr) - free(cdarr); + else { + /* Non-diffuse ("specular") SDF transmission */ + if ((err = SDsampBSDF(&bsdfVal, nd.vray, xi2, + SDsampSp | SDsampT, nd.sd))) + objerror(mat, USER, transSDError(err)); + + /* Apply pattern to spectral component */ + ccy2rgb(&bsdfVal.spec, bsdfVal.cieY, bsdfRGB); + multcolor(bsdfRGB, rayIn -> pcol); + flipsurface(rayIn); /* Necessary? */ + photonRay(rayIn, &rayOut, PMAP_SPECTRANS, bsdfRGB); + } + + /* Xform outgoing dir to world coords */ + if ((err = SDmapDir(rayOut.rdir, nd.fromloc, nd.vray))) { + objerror(mat, USER, transSDError(err)); + return 0; + } } - - /* Clean up BSDF */ + + /* Clean up */ SDfreeCache(nd.sd); tracePhoton(&rayOut); return 0; } +#endif