--- ray/src/rt/pmap.c 2018/03/20 19:55:33 2.14 +++ ray/src/rt/pmap.c 2021/02/19 02:10:35 2.18 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: pmap.c,v 2.14 2018/03/20 19:55:33 rschregle Exp $"; +static const char RCSid[] = "$Id: pmap.c,v 2.18 2021/02/19 02:10:35 rschregle Exp $"; #endif @@ -9,11 +9,14 @@ static const char RCSid[] = "$Id: pmap.c,v 2.14 2018/0 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") (c) Lucerne University of Applied Sciences and Arts, - supported by the Swiss National Science Foundation (SNSF, #147053) + supported by the Swiss National Science Foundation + (SNSF #147053, "Daylight Redirecting Components") ====================================================================== - $Id: pmap.c,v 2.14 2018/03/20 19:55:33 rschregle Exp $ + $Id: pmap.c,v 2.18 2021/02/19 02:10:35 rschregle Exp $ */ @@ -25,6 +28,7 @@ static const char RCSid[] = "$Id: pmap.c,v 2.14 2018/0 #include "pmapbias.h" #include "pmapdiag.h" #include "otypes.h" +#include "otspecial.h" #include #if NIX #include @@ -50,9 +54,10 @@ static int photonParticipate (RAY *ray) or 0 if absorbed and $*%&ed. Analogon to rayparticipate(). */ { int i; - RREAL cosTheta, cosPhi, du, dv; + RREAL xi1, cosTheta, phi, du, dv; const float cext = colorAvg(ray -> cext), - albedo = colorAvg(ray -> albedo); + albedo = colorAvg(ray -> albedo), + gecc = ray -> gecc, gecc2 = sqr(gecc); FVECT u, v; COLOR cvext; @@ -60,6 +65,16 @@ static int photonParticipate (RAY *ray) ray -> rmax = -log(pmapRandom(mediumState)) / cext; while (!localhit(ray, &thescene)) { + if (!incube(&thescene, ray -> rop)) { + /* Terminate photon if it has leaked from the scene */ +#ifdef DEBUG_PMAP + fprintf(stderr, + "Volume photon leaked from scene at [%.3f %.3f %.3f]\n", + ray -> rop [0], ray -> rop [1], ray -> rop [2]); +#endif + return 0; + } + setcolor(cvext, exp(-ray -> rmax * ray -> cext [0]), exp(-ray -> rmax * ray -> cext [1]), exp(-ray -> rmax * ray -> cext [2])); @@ -69,7 +84,13 @@ static int photonParticipate (RAY *ray) colorNorm(ray -> rcol); VCOPY(ray -> rorg, ray -> rop); +#if 0 if (albedo > FTINY && ray -> rlvl > 0) +#else + /* Store volume photons unconditionally in mist to also account for + direct inscattering from sources */ + if (albedo > FTINY) +#endif /* Add to volume photon map */ newPhoton(volumePmap, ray); @@ -82,17 +103,17 @@ static int photonParticipate (RAY *ray) scalecolor(ray -> rcol, 1 / albedo); /* Scatter photon */ - cosTheta = ray -> gecc <= FTINY ? 2 * pmapRandom(scatterState) - 1 - : 1 / (2 * ray -> gecc) * - (1 + ray -> gecc * ray -> gecc - - (1 - ray -> gecc * ray -> gecc) / - (1 - ray -> gecc + 2 * ray -> gecc * - pmapRandom(scatterState))); - - cosPhi = cos(2 * PI * pmapRandom(scatterState)); - du = dv = sqrt(1 - cosTheta * cosTheta); /* sin(theta) */ - du *= cosPhi; - dv *= sqrt(1 - cosPhi * cosPhi); /* sin(phi) */ + xi1 = pmapRandom(scatterState); + cosTheta = ray -> gecc <= FTINY + ? 2 * xi1 - 1 + : 0.5 / gecc * + (1 + gecc2 - sqr((1 - gecc2) / + (1 + gecc * (2 * xi1 - 1)))); + + phi = 2 * PI * pmapRandom(scatterState); + du = dv = sqrt(1 - sqr(cosTheta)); /* sin(theta) */ + du *= cos(phi); + dv *= sin(phi); /* Get axes u & v perpendicular to photon direction */ i = 0; @@ -106,19 +127,20 @@ static int photonParticipate (RAY *ray) for (i = 0; i < 3; i++) ray -> rdir [i] = du * u [i] + dv * v [i] + cosTheta * ray -> rdir [i]; + ray -> rlvl++; ray -> rmax = -log(pmapRandom(mediumState)) / cext; } - + + /* Passed through medium until intersecting local object */ setcolor(cvext, exp(-ray -> rot * ray -> cext [0]), exp(-ray -> rot * ray -> cext [1]), exp(-ray -> rot * ray -> cext [2])); /* Modify ray color and normalise */ multcolor(ray -> rcol, cvext); - colorNorm(ray -> rcol); - - /* Passed through medium */ + colorNorm(ray -> rcol); + return 1; } @@ -151,10 +173,15 @@ void tracePhoton (RAY *ray) if (localhit(ray, &thescene)) { mod = ray -> ro -> omod; - if (port && ray -> ro != port) { + /* XXX: Is port -> omod != mod sufficient here? Probably not... */ + if ( + port && ray -> ro != port && + findmaterial(port) != findmaterial(ray -> ro) + ) { /* !!! PHOTON PORT REJECTION SAMPLING HACK !!! - * Terminate photon if emitted from port without intersecting it; - * this can happen when the port's partitions extend beyond its + * Terminate photon if emitted from port without intersecting it or + * its other associated surfaces or same material. + * This can happen when the port's partitions extend beyond its * actual geometry, e.g. with polygons. Since the total flux * relayed by the port is based on the (in this case) larger * partition area, it is overestimated; terminating these photons @@ -444,7 +471,18 @@ void distribPhotons (PhotonMap **pmaps, unsigned numPr pmapSeed(randSeed + (proc + 3) % numProc, mediumState); pmapSeed(randSeed + (proc + 4) % numProc, scatterState); pmapSeed(randSeed + (proc + 5) % numProc, rouletteState); - + +#ifdef DEBUG_PMAP + /* Output child process PID after random delay to prevent corrupted + * console output due to race condition */ + usleep(1e6 * pmapRandom(rouletteState)); + fprintf(stderr, "Proc %d: PID = %d " + "(waiting 10 sec to attach debugger...)\n", + proc, getpid()); + /* Allow time for debugger to attach to child process */ + sleep(10); +#endif + for (t = 0; t < NUM_PMAP_TYPES; t++) lastNumPhotons [t] = 0;