--- ray/src/rt/pmap.c 2016/09/26 20:19:30 2.12 +++ ray/src/rt/pmap.c 2021/02/19 02:10:35 2.18 @@ -1,22 +1,25 @@ #ifndef lint -static const char RCSid[] = "$Id: pmap.c,v 2.12 2016/09/26 20:19:30 greg Exp $"; +static const char RCSid[] = "$Id: pmap.c,v 2.18 2021/02/19 02:10:35 rschregle Exp $"; #endif + /* ====================================================================== Photon map main module 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.12 2016/09/26 20:19:30 greg Exp $ + $Id: pmap.c,v 2.18 2021/02/19 02:10:35 rschregle Exp $ */ - #include "pmap.h" #include "pmapmat.h" #include "pmapsrc.h" @@ -25,13 +28,15 @@ static const char RCSid[] = "$Id: pmap.c,v 2.12 2016/0 #include "pmapbias.h" #include "pmapdiag.h" #include "otypes.h" +#include "otspecial.h" #include -#include -#include -#include +#if NIX + #include + #include + #include +#endif - void savePmaps (const PhotonMap **pmaps, int argc, char **argv) { unsigned t; @@ -49,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; @@ -59,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])); @@ -68,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); @@ -81,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; @@ -105,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; } @@ -127,7 +150,15 @@ void tracePhoton (RAY *ray) /* Follow photon as it bounces around... */ { long mod; - OBJREC* mat; + OBJREC *mat, *port = NULL; + + if (!ray -> parent) { + /* !!! PHOTON PORT REJECTION SAMPLING HACK: get photon port for + * !!! primary ray from ray -> ro, then reset the latter to NULL so + * !!! as not to interfere with localhit() */ + port = ray -> ro; + ray -> ro = NULL; + } if (ray -> rlvl > photonMaxBounce) { #ifdef PMAP_RUNAWAY_WARN @@ -138,10 +169,31 @@ void tracePhoton (RAY *ray) if (colorAvg(ray -> cext) > FTINY && !photonParticipate(ray)) return; - + if (localhit(ray, &thescene)) { mod = ray -> ro -> omod; - + + /* 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 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 + * constitutes rejection sampling and thereby compensates any bias + * incurred by the overestimated flux. */ +#ifdef PMAP_PORTREJECT_WARN + sprintf(errmsg, "photon outside port %s", ray -> ro -> oname); + error(WARNING, errmsg); +#endif + return; + } + if ((ray -> clipset && inset(ray -> clipset, mod)) || mod == OVOID) { /* Transfer ray if modifier is VOID or clipped within antimatta */ RAY tray; @@ -173,11 +225,14 @@ static void preComputeGlobal (PhotonMap *pmap) repComplete = numPreComp = finalGather * pmap -> numPhotons; - if (photonRepTime) { - sprintf(errmsg, "Precomputing irradiance for %ld global photons...\n", + if (verbose) { + sprintf(errmsg, + "\nPrecomputing irradiance for %ld global photons\n", numPreComp); eputs(errmsg); +#if NIX fflush(stderr); +#endif } /* Copy photon map for precomputed photons */ @@ -204,7 +259,7 @@ static void preComputeGlobal (PhotonMap *pmap) for (i = 0; i < numPreComp; i++) { /* Get random photon from stratified distribution in source heap to - * avoid duplicates and clutering */ + * avoid duplicates and clustering */ pIdx = firstPhoton(pmap) + (unsigned long)((i + pmapRandom(pmap -> randState)) / finalGather); @@ -242,9 +297,11 @@ static void preComputeGlobal (PhotonMap *pmap) deletePhotons(pmap); memcpy(pmap, &nuPmap, sizeof(PhotonMap)); - if (photonRepTime) { - eputs("Rebuilding precomputed photon map...\n"); + if (verbose) { + eputs("\nRebuilding precomputed photon map\n"); +#if NIX fflush(stderr); +#endif } /* Rebuild underlying data structure, destroying heap */ @@ -263,7 +320,7 @@ typedef struct { void distribPhotons (PhotonMap **pmaps, unsigned numProc) { EmissionMap emap; - char errmsg2 [128], shmFname [255]; + char errmsg2 [128], shmFname [PMAP_TMPFNLEN]; unsigned t, srcIdx, proc; double totalFlux = 0; int shmFile, stat, pid; @@ -295,41 +352,51 @@ void distribPhotons (PhotonMap **pmaps, unsigned numPr initPhotonHeap(pmaps [t]); /* Per-subprocess target count */ pmaps [t] -> distribTarget /= numProc; + + if (!pmaps [t] -> distribTarget) + error(INTERNAL, "no photons to distribute in distribPhotons"); } initPhotonEmissionFuncs(); initPhotonScatterFuncs(); - /* Get photon ports if specified */ - if (ambincl == 1) - getPhotonPorts(); + /* Get photon ports from modifier list */ + getPhotonPorts(photonPortList); /* Get photon sensor modifiers */ getPhotonSensors(photonSensorList); +#if NIX /* Set up shared mem for photon counters (zeroed by ftruncate) */ -#if 0 - snprintf(shmFname, 255, PMAP_SHMFNAME, getpid()); - shmFile = shm_open(shmFname, O_CREAT | O_RDWR, S_IRUSR | S_IWUSR); -#else - strcpy(shmFname, PMAP_SHMFNAME); + strcpy(shmFname, PMAP_TMPFNAME); shmFile = mkstemp(shmFname); -#endif - if (shmFile < 0) - error(SYSTEM, "failed opening shared memory file in distribPhotons"); + if (shmFile < 0 || ftruncate(shmFile, sizeof(*photonCnt)) < 0) + error(SYSTEM, "failed shared mem init in distribPhotons"); - if (ftruncate(shmFile, sizeof(*photonCnt)) < 0) - error(SYSTEM, "failed setting shared memory size in distribPhotons"); - photonCnt = mmap(NULL, sizeof(*photonCnt), PROT_READ | PROT_WRITE, MAP_SHARED, shmFile, 0); if (photonCnt == MAP_FAILED) - error(SYSTEM, "failed mapping shared memory in distribPhotons"); + error(SYSTEM, "failed mapping shared memory in distribPhotons"); +#else + /* Allocate photon counters statically on Windoze */ + if (!(photonCnt = malloc(sizeof(PhotonCnt)))) + error(SYSTEM, "failed trivial malloc in distribPhotons"); + photonCnt -> numEmitted = photonCnt -> numComplete = 0; +#endif /* NIX */ - if (photonRepTime) - eputs("\n"); + if (verbose) { + sprintf(errmsg, "\nIntegrating flux from %d sources", nsources); + + if (photonPorts) { + sprintf(errmsg2, " via %d ports", numPhotonPorts); + strcat(errmsg, errmsg2); + } + + strcat(errmsg, "\n"); + eputs(errmsg); + } /* =================================================================== * FLUX INTEGRATION - Get total photon flux from light sources @@ -343,20 +410,22 @@ void distribPhotons (PhotonMap **pmaps, unsigned numPr : NULL; photonPartition [emap.src -> so -> otype] (&emap); - if (photonRepTime) { - sprintf(errmsg, "Integrating flux from source %s ", + if (verbose) { + sprintf(errmsg, "\tIntegrating flux from source %s ", source [srcIdx].so -> oname); - + if (emap.port) { sprintf(errmsg2, "via port %s ", photonPorts [portCnt].so -> oname); strcat(errmsg, errmsg2); } - - sprintf(errmsg2, "(%lu partitions)...\n", emap.numPartitions); + + sprintf(errmsg2, "(%lu partitions)\n", emap.numPartitions); strcat(errmsg, errmsg2); eputs(errmsg); +#if NIX fflush(stderr); +#endif } for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions; @@ -371,12 +440,25 @@ void distribPhotons (PhotonMap **pmaps, unsigned numPr if (totalFlux < FTINY) error(USER, "zero flux from light sources"); + + /* Record start time for progress reports */ + repStartTime = time(NULL); + if (verbose) { + sprintf(errmsg, "\nPhoton distribution @ %d procs\n", numProc); + eputs(errmsg); + } + /* MAIN LOOP */ for (proc = 0; proc < numProc; proc++) { +#if NIX if (!(pid = fork())) { - /* SUBPROCESS ENTERS HERE. - All opened and memory mapped files are inherited */ + /* SUBPROCESS ENTERS HERE; open and mmapped files inherited */ +#else + if (1) { + /* No subprocess under Windoze */ +#endif + /* Local photon counters for this subprocess */ unsigned passCnt = 0, prePassCnt = 0; unsigned long lastNumPhotons [NUM_PMAP_TYPES]; unsigned long localNumEmitted = 0; /* Num photons emitted by this @@ -384,12 +466,23 @@ void distribPhotons (PhotonMap **pmaps, unsigned numPr /* Seed RNGs from PID for decorellated photon distribution */ pmapSeed(randSeed + proc, partState); - pmapSeed(randSeed + proc, emitState); - pmapSeed(randSeed + proc, cntState); - pmapSeed(randSeed + proc, mediumState); - pmapSeed(randSeed + proc, scatterState); - pmapSeed(randSeed + proc, rouletteState); - + pmapSeed(randSeed + (proc + 1) % numProc, emitState); + pmapSeed(randSeed + (proc + 2) % numProc, cntState); + 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; @@ -409,9 +502,7 @@ void distribPhotons (PhotonMap **pmaps, unsigned numPr * iterations; make it clear to user which photon maps are * missing so (s)he can check geometry and materials */ if (++prePassCnt > maxPreDistrib) { - sprintf(errmsg, - "proc %d, source %s: too many prepasses", - proc, source [srcIdx].so -> oname); + sprintf(errmsg, "proc %d: too many prepasses", proc); for (t = 0; t < NUM_PMAP_TYPES; t++) if (pmaps [t] && !pmaps [t] -> numPhotons) { @@ -483,7 +574,7 @@ void distribPhotons (PhotonMap **pmaps, unsigned numPr break; } - /* Update shared completion counter for prog.report by parent */ + /* Update shared completion counter for progreport by parent */ photonCnt -> numComplete += numEmit; /* PHOTON DISTRIBUTION LOOP */ @@ -496,25 +587,29 @@ void distribPhotons (PhotonMap **pmaps, unsigned numPr ? photonPorts + portCnt : NULL; photonPartition [emap.src -> so -> otype] (&emap); - if (photonRepTime && !proc) { + if (verbose && !proc) { + /* Output from subproc 0 only to avoid race condition + * on console I/O */ if (!passCnt) - sprintf(errmsg, "PREPASS %d on source %s ", + sprintf(errmsg, "\tPREPASS %d on source %s ", prePassCnt, source [srcIdx].so -> oname); else - sprintf(errmsg, "MAIN PASS on source %s ", + sprintf(errmsg, "\tMAIN PASS on source %s ", source [srcIdx].so -> oname); - + if (emap.port) { sprintf(errmsg2, "via port %s ", photonPorts [portCnt].so -> oname); strcat(errmsg, errmsg2); } - + sprintf(errmsg2, "(%lu partitions)\n", emap.numPartitions); strcat(errmsg, errmsg2); eputs(errmsg); +#if NIX fflush(stderr); +#endif } for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions; @@ -529,8 +624,7 @@ void distribPhotons (PhotonMap **pmaps, unsigned numPr /* Number of photons to emit from ziss partishunn -- * proportional to flux; photon ray weight and scalar - * flux are uniform (the latter only varying in RGB). - * */ + * flux are uniform (latter only varying in RGB). */ partNumEmit = numEmit * colorAvg(emap.partFlux) / totalFlux; partEmitCnt = (unsigned long)partNumEmit; @@ -551,9 +645,17 @@ void distribPhotons (PhotonMap **pmaps, unsigned numPr /* Emit photon based on PDF and trace through scene * until absorbed/leaked */ emitPhoton(&emap, &photonRay); +#if 1 + if (emap.port) + /* !!! PHOTON PORT REJECTION SAMPLING HACK: set + * !!! photon port as fake hit object for + * !!! primary ray to check for intersection in + * !!! tracePhoton() */ + photonRay.ro = emap.port -> so; +#endif tracePhoton(&photonRay); } - + /* Update shared global photon count for each pmap */ for (t = 0; t < NUM_PMAP_TYPES; t++) if (pmaps [t]) { @@ -561,6 +663,15 @@ void distribPhotons (PhotonMap **pmaps, unsigned numPr pmaps [t] -> numPhotons - lastNumPhotons [t]; lastNumPhotons [t] = pmaps [t] -> numPhotons; } +#if !NIX + /* Synchronous progress report on Windoze */ + if (!proc && photonRepTime > 0 && + time(NULL) >= repLastTime + photonRepTime) { + repEmitted = repProgress = photonCnt -> numEmitted; + repComplete = photonCnt -> numComplete; + pmapDistribReport(); + } +#endif } portCnt++; @@ -576,55 +687,40 @@ void distribPhotons (PhotonMap **pmaps, unsigned numPr break; } - if (t >= NUM_PMAP_TYPES) { + if (t >= NUM_PMAP_TYPES) /* No empty photon maps found; now do pass 2 */ passCnt++; -#if 0 - if (photonRepTime) - eputs("\n"); -#endif - } } while (passCnt < 2); - /* Unmap shared photon counters */ -#if 0 - munmap(photonCnt, sizeof(*photonCnt)); - close(shmFile); -#endif - - /* Flush heap buffa for every pmap one final time; this is required - * to prevent data corruption! */ + /* Flush heap buffa for every pmap one final time; + * avoids potential data corruption! */ for (t = 0; t < NUM_PMAP_TYPES; t++) if (pmaps [t]) { -#if 0 - eputs("Final flush\n"); -#endif flushPhotonHeap(pmaps [t]); - fclose(pmaps [t] -> heap); + /* Heap file closed automatically on exit + fclose(pmaps [t] -> heap); */ #ifdef DEBUG_PMAP - sprintf(errmsg, "Proc %d: total %ld photons\n", getpid(), + sprintf(errmsg, "Proc %d: total %ld photons\n", proc, pmaps [t] -> numPhotons); eputs(errmsg); #endif } - +#if NIX + /* Terminate subprocess */ exit(0); +#endif } else if (pid < 0) error(SYSTEM, "failed to fork subprocess in distribPhotons"); } +#if NIX /* PARENT PROCESS CONTINUES HERE */ - /* Record start time and enable progress report signal handler */ - repStartTime = time(NULL); #ifdef SIGCONT + /* Enable progress report signal handler */ signal(SIGCONT, pmapDistribReport); -#endif - - if (photonRepTime) - eputs("\n"); - - /* Wait for subprocesses to complete while reporting progress */ +#endif + /* Wait for subprocesses complete while reporting progress */ proc = numProc; while (proc) { while (waitpid(-1, &stat, WNOHANG) > 0) { @@ -637,21 +733,19 @@ void distribPhotons (PhotonMap **pmaps, unsigned numPr /* Nod off for a bit and update progress */ sleep(1); - /* Update progress report from shared subprocess counters */ + + /* Asynchronous progress report from shared subprocess counters */ repEmitted = repProgress = photonCnt -> numEmitted; - repComplete = photonCnt -> numComplete; + repComplete = photonCnt -> numComplete; + repProgress = repComplete = 0; for (t = 0; t < NUM_PMAP_TYPES; t++) if ((pm = pmaps [t])) { -#if 0 - /* Get photon count from heapfile size for progress update */ - fseek(pm -> heap, 0, SEEK_END); - pm -> numPhotons = ftell(pm -> heap) / sizeof(Photon); */ -#else /* Get global photon count from shmem updated by subprocs */ - pm -> numPhotons = photonCnt -> numPhotons [t]; -#endif + repProgress += pm -> numPhotons = photonCnt -> numPhotons [t]; + repComplete += pm -> distribTarget; } + repComplete *= numProc; if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime) pmapDistribReport(); @@ -659,41 +753,52 @@ void distribPhotons (PhotonMap **pmaps, unsigned numPr else signal(SIGCONT, pmapDistribReport); #endif } +#endif /* NIX */ /* =================================================================== * POST-DISTRIBUTION - Set photon flux and build data struct for photon * storage, etc. * =================================================================== */ #ifdef SIGCONT + /* Reset signal handler */ signal(SIGCONT, SIG_DFL); #endif free(emap.samples); - /* Set photon flux (repProgress is total num emitted) */ + /* Set photon flux */ totalFlux /= photonCnt -> numEmitted; - +#if NIX /* Photon counters no longer needed, unmap shared memory */ munmap(photonCnt, sizeof(*photonCnt)); close(shmFile); -#if 0 - shm_unlink(shmFname); -#else unlink(shmFname); +#else + free(photonCnt); #endif - + if (verbose) + eputs("\n"); + for (t = 0; t < NUM_PMAP_TYPES; t++) if (pmaps [t]) { - if (photonRepTime) { - sprintf(errmsg, "\nBuilding %s photon map...\n", pmapName [t]); + if (verbose) { + sprintf(errmsg, "Building %s photon map\n", pmapName [t]); eputs(errmsg); +#if NIX fflush(stderr); +#endif } /* Build underlying data structure; heap is destroyed */ buildPhotonMap(pmaps [t], &totalFlux, NULL, numProc); } - + /* Precompute photon irradiance if necessary */ - if (preCompPmap) + if (preCompPmap) { + if (verbose) + eputs("\n"); preComputeGlobal(preCompPmap); + } + + if (verbose) + eputs("\n"); }