--- ray/src/rt/pmapcontrib.c 2015/04/21 19:16:51 2.2 +++ ray/src/rt/pmapcontrib.c 2018/03/16 21:00:09 2.16 @@ -1,18 +1,21 @@ +#ifndef lint +static const char RCSid[] = "$Id: pmapcontrib.c,v 2.16 2018/03/16 21:00:09 rschregle Exp $"; +#endif + /* - ================================================================== - Photon map support for light source contributions + ====================================================================== + Photon map for light source contributions 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: pmapcontrib.c,v 2.2 2015/04/21 19:16:51 greg Exp $ + $Id: pmapcontrib.c,v 2.16 2018/03/16 21:00:09 rschregle Exp $ */ #include "pmapcontrib.h" -#include "pmap.h" #include "pmapmat.h" #include "pmapsrc.h" #include "pmaprand.h" @@ -20,472 +23,713 @@ #include "pmapdiag.h" #include "rcontrib.h" #include "otypes.h" +#if NIX + #include + #include +#endif - -extern int contrib; /* coeff/contrib flag */ - - - -static void setPmapContribParams (PhotonMap *pmap, LUTAB *srcContrib) -/* Set parameters for light source contributions */ +static PhotonPrimaryIdx newPhotonPrimary (PhotonMap *pmap, + const RAY *primRay, + FILE *primHeap) +/* Add primary ray for emitted photon and save light source index, origin on + * source, and emitted direction; used by contrib photons. The current + * primary is stored in pmap -> lastPrimary. If the previous primary + * contributed photons (has srcIdx >= 0), it's appended to primHeap. If + * primRay == NULL, the current primary is still flushed, but no new primary + * is set. Returns updated primary counter pmap -> numPrimary. */ { - /* Set light source modifier list and appropriate callback to extract - * their contributions from the photon map */ - if (pmap) { - pmap -> srcContrib = srcContrib; - pmap -> lookup = photonContrib; - /* Ensure we get all requested photon contribs during lookups */ - pmap -> gatherTolerance = 1.0; + if (!pmap || !primHeap) + return 0; + + /* Check if last primary ray has spawned photons (srcIdx >= 0, see + * newPhoton()), in which case we save it to the primary heap file + * before clobbering it */ + if (pmap -> lastPrimary.srcIdx >= 0) { + if (!fwrite(&pmap -> lastPrimary, sizeof(PhotonPrimary), 1, primHeap)) + error(SYSTEM, "failed writing photon primary in newPhotonPrimary"); + + pmap -> numPrimary++; + if (pmap -> numPrimary > PMAP_MAXPRIMARY) + error(INTERNAL, "photon primary overflow in newPhotonPrimary"); } -} + /* Mark unused with negative source index until path spawns a photon (see + * newPhoton()) */ + pmap -> lastPrimary.srcIdx = -1; + + if (primRay) { + FVECT dvec; - -static void checkPmapContribs (const PhotonMap *pmap, LUTAB *srcContrib) -/* Check modifiers for light source contributions */ -{ - const PhotonPrimary *primary = pmap -> primary; - OBJREC *srcMod; - unsigned long i, found = 0; - - /* Make sure at least one of the modifiers is actually in the pmap, - * otherwise findPhotons() winds up in an infinite loop! */ - for (i = pmap -> primarySize; i; --i, ++primary) { - if (primary -> srcIdx < 0 || primary -> srcIdx >= nsources) - error(INTERNAL, "invalid light source index in photon map"); - - srcMod = objptr(source [primary -> srcIdx].so -> omod); - if ((MODCONT*)lu_find(srcContrib, srcMod -> oname) -> data) - ++found; +#ifdef PMAP_PRIMARYDIR + /* Reverse incident direction to point to light source */ + dvec [0] = -primRay -> rdir [0]; + dvec [1] = -primRay -> rdir [1]; + dvec [2] = -primRay -> rdir [2]; + pmap -> lastPrimary.dir = encodedir(dvec); +#endif +#ifdef PMAP_PRIMARYPOS + VCOPY(pmap -> lastPrimary.pos, primRay -> rop); +#endif } - if (!found) - error(USER, "modifiers not in photon map"); + return pmap -> numPrimary; } - - -void initPmapContrib (LUTAB *srcContrib, unsigned numSrcContrib) + + +#ifdef DEBUG_PMAP +static int checkPrimaryHeap (FILE *file) +/* Check heap for ordered primaries */ { - unsigned t; + Photon p, lastp; + int i, dup; - for (t = 0; t < NUM_PMAP_TYPES; t++) - if (photonMaps [t] && t != PMAP_TYPE_CONTRIB) { - sprintf(errmsg, "%s photon map does not support contributions", - pmapName [t]); - error(USER, errmsg); - } + rewind(file); + memset(&lastp, 0, sizeof(lastp)); - /* Get params */ - setPmapContribParams(contribPmap, srcContrib); - - if (contribPhotonMapping) { - if (contribPmap -> maxGather < numSrcContrib) { - /* Adjust density estimate bandwidth if lower than modifier - * count, otherwise contributions are missing */ - error(WARNING, "contrib density estimate bandwidth too low, " - "adjusting to modifier count"); - contribPmap -> maxGather = numSrcContrib; + while (fread(&p, sizeof(p), 1, file)) { + dup = 1; + + for (i = 0; i <= 2; 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", + p.pos [0], p.pos [1], p.pos [2]); + error(WARNING, errmsg); + } + + dup &= p.pos [i] == lastp.pos [i]; } - /* Sanity check */ - checkPmapContribs(contribPmap, srcContrib); + 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 -void photonContrib (PhotonMap *pmap, RAY *ray, COLOR irrad) -/* Sum up light source contributions from photons in pmap->srcContrib */ +static PhotonPrimaryIdx buildPrimaries (PhotonMap *pmap, FILE **primaryHeap, + char **primaryHeapFname, + PhotonPrimaryIdx *primaryOfs, + unsigned numHeaps) +/* Consolidate per-subprocess photon primary heaps into the primary array + * pmap -> primaries. Returns offset for primary index linearisation in + * numPrimary. The heap files in primaryHeap are closed on return. */ { - unsigned i; - PhotonSQNode *sq; - float r, invArea; - RREAL rayCoeff [3]; - FVECT rdir, rop; - - setcolor(irrad, 0, 0, 0); - - if (!pmap -> maxGather) - return; - - /* Ignore sources */ - if (ray -> ro) - if (islight(objptr(ray -> ro -> omod) -> otype)) - return; - - /* Set context for binning function evaluation and get cumulative path - * coefficient up to photon lookup point */ - worldfunc(RCCONTEXT, ray); - raycontrib(rayCoeff, ray, PRIMARY); - - /* Save incident ray's direction and hitpoint */ - VCOPY(rdir, ray -> rdir); - VCOPY(rop, ray -> rop); - - /* Lookup photons */ - pmap -> squeueEnd = 0; - findPhotons(pmap, ray); + PhotonPrimaryIdx heapLen; + unsigned heap; - /* Need at least 2 photons */ - if (pmap -> squeueEnd < 2) { - #ifdef PMAP_NONEFOUND - sprintf(errmsg, "no photons found on %s at (%.3f, %.3f, %.3f)", - ray -> ro ? ray -> ro -> oname : "", - ray -> rop [0], ray -> rop [1], ray -> rop [2]); - error(WARNING, errmsg); - #endif + if (!pmap || !primaryHeap || !primaryOfs || !numHeaps) + return 0; - return; - } - - /* Average (squared) radius between furthest two photons to improve - * accuracy and get inverse search area 1 / (PI * r^2), with extra - * normalisation factor 1 / PI for ambient calculation */ - sq = pmap -> squeue + 1; - r = max(sq -> dist, (sq + 1) -> dist); - r = 0.25 * (pmap -> maxDist + r + 2 * sqrt(pmap -> maxDist * r)); - invArea = 1 / (PI * PI * r); + pmap -> numPrimary = 0; - /* Skip the extra photon */ - for (i = 1 ; i < pmap -> squeueEnd; i++, sq++) { - COLOR flux; + for (heap = 0; heap < numHeaps; heap++) { + primaryOfs [heap] = pmap -> numPrimary; - /* Get photon's contribution to density estimate */ - getPhotonFlux(sq -> photon, flux); - scalecolor(flux, invArea); -#ifdef PMAP_EPANECHNIKOV - /* Apply Epanechnikov kernel to photon flux (dists are squared) */ - scalecolor(flux, 2 * (1 - sq -> dist / r)); -#endif - addcolor(irrad, flux); - - if (pmap -> srcContrib) { - const PhotonPrimary *primary = pmap -> primary + - sq -> photon -> primary; - OBJREC *srcMod = objptr(source [primary -> srcIdx].so -> omod); - MODCONT *srcContrib = (MODCONT*)lu_find(pmap -> srcContrib, - srcMod -> oname) -> data; - - if (srcContrib) { - /* Photon's emitting light source has modifier whose - * contributions are sought */ - int srcBin; + if (fseek(primaryHeap [heap], 0, SEEK_END) < 0) + error(SYSTEM, "failed photon primary seek in buildPrimaries"); + pmap -> numPrimary += heapLen = ftell(primaryHeap [heap]) / + sizeof(PhotonPrimary); - /* Set incident dir and origin of photon's primary ray on - * light source for dummy shadow ray, and evaluate binning - * function */ - VCOPY(ray -> rdir, primary -> dir); - VCOPY(ray -> rop, primary -> org); - srcBin = evalue(srcContrib -> binv) + .5; + pmap -> primaries = realloc(pmap -> primaries, + pmap -> numPrimary * + sizeof(PhotonPrimary)); + if (!pmap -> primaries) + error(SYSTEM, "failed photon primary alloc in buildPrimaries"); - if (srcBin < 0 || srcBin >= srcContrib -> nbins) { - error(WARNING, "bad bin number (ignored)"); - continue; - } - - if (!contrib) { - /* Ray coefficient mode; normalise by light source radiance - * after applying distrib pattern */ - int j; - raytexture(ray, srcMod -> omod); - setcolor(ray -> rcol, srcMod -> oargs.farg [0], - srcMod -> oargs.farg [1], srcMod -> oargs.farg [2]); - multcolor(ray -> rcol, ray -> pcol); - for (j = 0; j < 3; j++) - flux [j] = ray -> rcol [j] ? flux [j] / ray -> rcol [j] - : 0; - } - - multcolor(flux, rayCoeff); - addcolor(srcContrib -> cbin [srcBin], flux); - } - else fprintf(stderr, "Skipped contrib from %s\n", srcMod -> oname); - } + rewind(primaryHeap [heap]); + if (fread(pmap -> primaries + primaryOfs [heap], sizeof(PhotonPrimary), + heapLen, primaryHeap [heap]) != heapLen) + error(SYSTEM, "failed reading photon primaries in buildPrimaries"); + + fclose(primaryHeap [heap]); + unlink(primaryHeapFname [heap]); } - /* Restore incident ray's direction and hitpoint */ - VCOPY(ray -> rdir, rdir); - VCOPY(ray -> rop, rop); - - return; -} + return pmap -> numPrimary; +} -void distribPhotonContrib (PhotonMap* pm) -{ - EmissionMap emap; - char errmsg2 [128]; - unsigned srcIdx; - double *srcFlux; /* Emitted flux per light source */ - const double srcDistribTarget = /* Target photon count per source */ - nsources ? (double)pm -> distribTarget / nsources : 0; +/* Defs for photon emission counter array passed by sub-processes to parent + * via shared memory */ +typedef unsigned long PhotonContribCnt; +/* Indices for photon emission counter array: num photons stored and num + * emitted per source */ +#define PHOTONCNT_NUMPHOT 0 +#define PHOTONCNT_NUMEMIT(n) (1 + n) + + + + + + +void distribPhotonContrib (PhotonMap* pm, unsigned numProc) +{ + EmissionMap emap; + char errmsg2 [128], shmFname [PMAP_TMPFNLEN]; + unsigned srcIdx, proc; + int shmFile, stat, pid; + double *srcFlux, /* Emitted flux per light source */ + srcDistribTarget; /* Target photon count per source */ + PhotonContribCnt *photonCnt; /* Photon emission counter array */ + unsigned photonCntSize = sizeof(PhotonContribCnt) * + PHOTONCNT_NUMEMIT(nsources); + FILE **primaryHeap = NULL; + char **primaryHeapFname = NULL; + PhotonPrimaryIdx *primaryOfs = NULL; + if (!pm) - error(USER, "no photon map defined"); + error(USER, "no photon map defined in distribPhotonContrib"); if (!nsources) - error(USER, "no light sources"); - + error(USER, "no light sources in distribPhotonContrib"); + + if (nsources > MAXMODLIST) + error(USER, "too many light sources in distribPhotonContrib"); + /* Allocate photon flux per light source; this differs for every * source as all sources contribute the same number of distributed * photons (srcDistribTarget), hence the number of photons emitted per * source does not correlate with its emitted flux. The resulting flux * per photon is therefore adjusted individually for each source. */ if (!(srcFlux = calloc(nsources, sizeof(double)))) - error(SYSTEM, "cannot allocate source flux"); + error(SYSTEM, "can't allocate source flux in distribPhotonContrib"); - /* ================================================================ - * INITIALISASHUNN - Set up emisshunn and scattering funcs - * ================================================================ */ + /* =================================================================== + * INITIALISATION - Set up emission and scattering funcs + * =================================================================== */ emap.samples = NULL; emap.src = NULL; emap.maxPartitions = MAXSPART; emap.partitions = (unsigned char*)malloc(emap.maxPartitions >> 1); if (!emap.partitions) - error(USER, "can't allocate source partitions"); + error(USER, "can't allocate source partitions in distribPhotonContrib"); + /* Initialise contrib photon map */ initPhotonMap(pm, PMAP_TYPE_CONTRIB); + initPhotonHeap(pm); initPhotonEmissionFuncs(); initPhotonScatterFuncs(); + /* Per-subprocess / per-source target counts */ + pm -> distribTarget /= numProc; + srcDistribTarget = nsources ? (double)pm -> distribTarget / nsources : 0; + + if (!pm -> distribTarget) + error(INTERNAL, "no photons to distribute in distribPhotonContrib"); + /* Get photon ports if specified */ if (ambincl == 1) getPhotonPorts(); /* Get photon sensor modifiers */ getPhotonSensors(photonSensorList); + +#if NIX + /* Set up shared mem for photon counters (zeroed by ftruncate) */ + strcpy(shmFname, PMAP_TMPFNAME); + shmFile = mkstemp(shmFname); - /* Seed RNGs for photon distribution */ - pmapSeed(randSeed, partState); - pmapSeed(randSeed, emitState); - pmapSeed(randSeed, cntState); - pmapSeed(randSeed, mediumState); - pmapSeed(randSeed, scatterState); - pmapSeed(randSeed, rouletteState); + if (shmFile < 0 || ftruncate(shmFile, photonCntSize) < 0) + error(SYSTEM, "failed shared mem init in distribPhotonContrib"); - /* Record start time and enable progress report signal handler */ - repStartTime = time(NULL); - signal(SIGCONT, pmapDistribReport); + photonCnt = mmap(NULL, photonCntSize, PROT_READ | PROT_WRITE, + MAP_SHARED, shmFile, 0); + + if (photonCnt == MAP_FAILED) + error(SYSTEM, "failed shared mem mapping in distribPhotonContrib"); +#else + /* Allocate photon counters statically on Windoze */ + if (!(photonCnt = malloc(photonCntSize))) + error(SYSTEM, "failed trivial malloc in distribPhotonContrib"); - for (srcIdx = 0; srcIdx < nsources; srcIdx++) { - unsigned portCnt = 0, passCnt = 0, prePassCnt = 0; - double srcNumEmit = 0; /* # photons to emit from source */ - unsigned long srcNumDistrib = pm -> heapEnd; /* # photons stored */ + for (srcIdx = 0; srcIdx < PHOTONCNT_NUMEMIT(nsources); srcIdx++) + photonCnt [srcIdx] = 0; +#endif /* NIX */ - srcFlux [srcIdx] = repProgress = 0; + 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 flux emitted from sources/ports + * ============================================================= */ + for (srcIdx = 0; srcIdx < nsources; srcIdx++) { + unsigned portCnt = 0; + srcFlux [srcIdx] = 0; emap.src = source + srcIdx; - if (photonRepTime) - eputs("\n"); - - /* ============================================================= - * FLUX INTEGRATION - Get total flux emitted from light source - * ============================================================= */ - do { - emap.port = emap.src -> sflags & SDISTANT - ? photonPorts + portCnt : NULL; + do { /* Need at least one iteration if no ports! */ + emap.port = emap.src -> sflags & SDISTANT ? photonPorts + portCnt + : NULL; photonPartition [emap.src -> so -> otype] (&emap); - - if (photonRepTime) { - sprintf(errmsg, "Integrating flux from source %s (mod %s) ", - source [srcIdx].so -> oname, - objptr(source [srcIdx].so -> omod) -> oname); - + + 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; + for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions; emap.partitionCnt++) { initPhotonEmission(&emap, pdfSamples); srcFlux [srcIdx] += colorAvg(emap.partFlux); } portCnt++; - } while (portCnt < numPhotonPorts); - + } while (portCnt < numPhotonPorts); + if (srcFlux [srcIdx] < FTINY) { sprintf(errmsg, "source %s has zero emission", source [srcIdx].so -> oname); error(WARNING, errmsg); } - else { - /* ========================================================== - * 2-PASS PHOTON DISTRIBUTION - * Pass 1 (pre): emit fraction of target photon count - * Pass 2 (main): based on outcome of pass 1, estimate - * remaining number of photons to emit to - * approximate target count - * ========================================================== */ - do { - if (!passCnt) { - /* INIT PASS 1 */ - if (++prePassCnt > maxPreDistrib) { - /* Warn if no photons contributed after sufficient - * iterations */ - sprintf(errmsg, "too many prepasses, no photons " - "from source %s", source [srcIdx].so -> oname); - error(WARNING, errmsg); - break; - } + } + + /* Allocate & init per-subprocess primary heap files */ + primaryHeap = calloc(numProc, sizeof(FILE*)); + primaryHeapFname = calloc(numProc, sizeof(char*)); + primaryOfs = calloc(numProc, sizeof(PhotonPrimaryIdx)); + if (!primaryHeap || !primaryHeapFname || !primaryOfs) + error(SYSTEM, "failed primary heap allocation in " + "distribPhotonContrib"); + + for (proc = 0; proc < numProc; proc++) { + primaryHeapFname [proc] = malloc(PMAP_TMPFNLEN); + if (!primaryHeapFname [proc]) + error(SYSTEM, "failed primary heap file allocation in " + "distribPhotonContrib"); - /* Num to emit is fraction of target count */ - srcNumEmit = preDistrib * srcDistribTarget; - } + mktemp(strcpy(primaryHeapFname [proc], PMAP_TMPFNAME)); + if (!(primaryHeap [proc] = fopen(primaryHeapFname [proc], "w+b"))) + error(SYSTEM, "failed opening primary heap file in " + "distribPhotonContrib"); + } - else { - /* INIT PASS 2 */ - /* Based on the outcome of the predistribution we can now - * figure out how many more photons we have to emit from - * the current source to meet the target count, - * srcDistribTarget. This value is clamped to 0 in case - * the target has already been exceeded in pass 1. - * srcNumEmit and srcNumDistrib is the number of photons - * emitted and distributed (stored) from the current - * source in pass 1, respectively. */ - srcNumDistrib = pm -> heapEnd - srcNumDistrib; - srcNumEmit *= srcNumDistrib - ? max(srcDistribTarget/srcNumDistrib, 1) - 1 - : 0; + /* Record start time for progress reports */ + repStartTime = time(NULL); - if (!srcNumEmit) - /* No photons left to distribute in main pass */ - break; - } - - /* Set completion count for progress report */ - repComplete = srcNumEmit + repProgress; - portCnt = 0; - - do { - emap.port = emap.src -> sflags & SDISTANT - ? photonPorts + portCnt : NULL; - photonPartition [emap.src -> so -> otype] (&emap); + if (verbose) { + sprintf(errmsg, "\nPhoton distribution @ %d procs\n", numProc); + eputs(errmsg); + } - if (photonRepTime) { - if (!passCnt) - sprintf(errmsg, "PREPASS %d on source %s (mod %s) ", - prePassCnt, source [srcIdx].so -> oname, - objptr(source[srcIdx].so->omod) -> oname); - else - sprintf(errmsg, "MAIN PASS on source %s (mod %s) ", - source [srcIdx].so -> oname, - objptr(source[srcIdx].so->omod) -> oname); - - if (emap.port) { - sprintf(errmsg2, "via port %s ", - photonPorts [portCnt].so -> oname); - strcat(errmsg, errmsg2); + /* MAIN LOOP */ + for (proc = 0; proc < numProc; proc++) { +#if NIX + if (!(pid = fork())) { + /* SUBPROCESS ENTERS HERE; opened and mmapped files inherited */ +#else + if (1) { + /* No subprocess under Windoze */ +#endif + /* Local photon counters for this subprocess */ + unsigned long lastNumPhotons = 0, localNumEmitted = 0; + double photonFluxSum = 0; /* Accum. photon flux */ + + /* Seed RNGs from PID for decorellated photon distribution */ + pmapSeed(randSeed + proc, partState); + 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 PMAP_SIGUSR + double partNumEmit; + unsigned long partEmitCnt; + double srcPhotonFlux, avgPhotonFlux; + unsigned portCnt, passCnt, prePassCnt; + float srcPreDistrib; + double srcNumEmit; /* # to emit from source */ + unsigned long srcNumDistrib; /* # stored */ + + void sigUsrDiags() + /* Loop diags via SIGUSR1 */ + { + sprintf(errmsg, + "********************* Proc %d Diags *********************\n" + "srcIdx = %d (%s)\nportCnt = %d (%s)\npassCnt = %d\n" + "srcFlux = %f\nsrcPhotonFlux = %f\navgPhotonFlux = %f\n" + "partNumEmit = %f\npartEmitCnt = %lu\n\n", + proc, srcIdx, findmaterial(source [srcIdx].so) -> oname, + portCnt, photonPorts [portCnt].so -> oname, + passCnt, srcFlux [srcIdx], srcPhotonFlux, avgPhotonFlux, + partNumEmit, partEmitCnt); + eputs(errmsg); + fflush(stderr); + } +#endif + +#if PMAP_SIGUSR + signal(SIGUSR1, sigUsrDiags); +#endif + + /* 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\n", proc, getpid()); + /* Allow time for debugger to attach to child process */ + sleep(10); + + /* ============================================================= + * 2-PASS PHOTON DISTRIBUTION + * Pass 1 (pre): emit fraction of target photon count + * Pass 2 (main): based on outcome of pass 1, estimate remaining + * number of photons to emit to approximate target + * count + * ============================================================= */ + for (srcIdx = 0; srcIdx < nsources; srcIdx++) { +#ifndef PMAP_SIGUSR + unsigned portCnt, passCnt = 0, prePassCnt = 0; + float srcPreDistrib = preDistrib; + double srcNumEmit = 0; /* # to emit from source */ + unsigned long srcNumDistrib = pm -> numPhotons; /* # stored */ +#else + passCnt = prePassCnt = 0; + srcPreDistrib = preDistrib; + srcNumEmit = 0; /* # to emit from source */ + srcNumDistrib = pm -> numPhotons; /* # stored */ +#endif + + if (srcFlux [srcIdx] < FTINY) + continue; + + while (passCnt < 2) { + if (!passCnt) { + /* INIT PASS 1 */ + if (++prePassCnt > maxPreDistrib) { + /* Warn if no photons contributed after sufficient + * iterations; only output from subprocess 0 to reduce + * console clutter */ + if (!proc) { + sprintf(errmsg, + "source %s: too many prepasses, skipped", + source [srcIdx].so -> oname); + error(WARNING, errmsg); + } + + break; } - sprintf(errmsg2, "(%lu partitions)...\n", - emap.numPartitions); - strcat(errmsg, errmsg2); - eputs(errmsg); - fflush(stderr); + /* Num to emit is fraction of target count */ + srcNumEmit = srcPreDistrib * srcDistribTarget; } - - for (emap.partitionCnt = 0; - emap.partitionCnt < emap.numPartitions; - emap.partitionCnt++) { - double partNumEmit; - unsigned long partEmitCnt; + else { + /* INIT PASS 2 */ +#ifndef PMAP_SIGUSR + double srcPhotonFlux, avgPhotonFlux; +#endif - /* Get photon origin within current source partishunn - * and build emission map */ - photonOrigin [emap.src -> so -> otype] (&emap); - initPhotonEmission(&emap, pdfSamples); - - /* Number of photons to emit from ziss partishunn; - * scale according to its normalised contribushunn to - * the emitted source flux */ - partNumEmit = srcNumEmit * colorAvg(emap.partFlux) / - srcFlux [srcIdx]; - partEmitCnt = (unsigned long)partNumEmit; - - /* Probabilistically account for fractional photons */ - if (pmapRandom(cntState) < partNumEmit - partEmitCnt) - partEmitCnt++; + /* Based on the outcome of the predistribution we can now + * figure out how many more photons we have to emit from + * the current source to meet the target count, + * srcDistribTarget. This value is clamped to 0 in case + * the target has already been exceeded in pass 1. + * srcNumEmit and srcNumDistrib is the number of photons + * emitted and distributed (stored) from the current + * source in pass 1, respectively. */ + srcNumDistrib = pm -> numPhotons - srcNumDistrib; + srcNumEmit *= srcNumDistrib + ? max(srcDistribTarget/srcNumDistrib, 1) - 1 + : 0; + + if (!srcNumEmit) + /* No photons left to distribute in main pass */ + break; - /* Integer counter avoids FP rounding errors */ - while (partEmitCnt--) { - RAY photonRay; + srcPhotonFlux = srcFlux [srcIdx] / srcNumEmit; + avgPhotonFlux = photonFluxSum / (srcIdx + 1); - /* Emit photon according to PDF (if any), allocate - * associated primary ray, and trace through scene - * until absorbed/leaked */ - emitPhoton(&emap, &photonRay); - addPhotonPrimary(pm, &photonRay); - tracePhoton(&photonRay); + if (avgPhotonFlux > FTINY && + srcPhotonFlux / avgPhotonFlux < FTINY) { + /* Skip source if its photon flux is grossly below the + * running average, indicating negligible contributions + * at the expense of excessive distribution time; only + * output from subproc 0 to reduce console clutter */ + if (!proc) { + sprintf(errmsg, + "source %s: itsy bitsy photon flux, skipped", + source [srcIdx].so -> oname); + error(WARNING, errmsg); + } + + srcNumEmit = 0; /* Or just break??? */ + } + + /* Update sum of photon flux per light source */ + photonFluxSum += srcPhotonFlux; + } + + portCnt = 0; + do { /* Need at least one iteration if no ports! */ + emap.src = source + srcIdx; + emap.port = emap.src -> sflags & SDISTANT + ? photonPorts + portCnt : NULL; + photonPartition [emap.src -> so -> otype] (&emap); + + if (verbose && !proc) { + /* Output from subproc 0 only to avoid race condition + * on console I/O */ + if (!passCnt) + sprintf(errmsg, "\tPREPASS %d on source %s ", + prePassCnt, source [srcIdx].so -> oname); + else + 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; + emap.partitionCnt++) { +#ifndef PMAP_SIGUSR + double partNumEmit; + unsigned long partEmitCnt; +#endif - /* Record progress */ - repProgress++; + /* Get photon origin within current source partishunn + * and build emission map */ + photonOrigin [emap.src -> so -> otype] (&emap); + initPhotonEmission(&emap, pdfSamples); - if (photonRepTime > 0 && - time(NULL) >= repLastTime + photonRepTime) + /* Number of photons to emit from ziss partishunn; + * scale according to its normalised contribushunn to + * the emitted source flux */ + partNumEmit = srcNumEmit * colorAvg(emap.partFlux) / + srcFlux [srcIdx]; + partEmitCnt = (unsigned long)partNumEmit; + + /* Probabilistically account for fractional photons */ + if (pmapRandom(cntState) < partNumEmit - partEmitCnt) + partEmitCnt++; + + /* Update local and shared global emission counter */ + photonCnt [PHOTONCNT_NUMEMIT(srcIdx)] += partEmitCnt; + localNumEmitted += partEmitCnt; + + /* Integer counter avoids FP rounding errors during + * iteration */ + while (partEmitCnt--) { + RAY photonRay; + + /* Emit photon according to PDF (if any), allocate + * associated primary ray, and trace through scene + * until absorbed/leaked; emitPhoton() sets the + * emitting light source index in photonRay */ + 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 + newPhotonPrimary(pm, &photonRay, primaryHeap[proc]); + /* Set subprocess index in photonRay for post- + * distrib primary index linearisation; this is + * propagated with the primary index in photonRay + * and set for photon hits by newPhoton() */ + PMAP_SETRAYPROC(&photonRay, proc); + tracePhoton(&photonRay); + } + + /* Update shared global photon count */ + photonCnt [PHOTONCNT_NUMPHOT] += pm -> numPhotons - + lastNumPhotons; + lastNumPhotons = pm -> numPhotons; +#if !NIX + /* Synchronous progress report on Windoze */ + if (!proc && photonRepTime > 0 && + time(NULL) >= repLastTime + photonRepTime) { + unsigned s; + repComplete = pm -> distribTarget * numProc; + repProgress = photonCnt [PHOTONCNT_NUMPHOT]; + + for (repEmitted = 0, s = 0; s < nsources; s++) + repEmitted += photonCnt [PHOTONCNT_NUMEMIT(s)]; + pmapDistribReport(); - #ifndef BSD - else signal(SIGCONT, pmapDistribReport); - #endif + } +#endif } - } - - portCnt++; - } while (portCnt < numPhotonPorts); - if (pm -> heapEnd == srcNumDistrib) - /* Double preDistrib in case no photons were stored - * for this source and redo pass 1 */ - preDistrib *= 2; - else { - /* Now do pass 2 */ - passCnt++; - if (photonRepTime) - eputs("\n"); + portCnt++; + } while (portCnt < numPhotonPorts); + + if (pm -> numPhotons == srcNumDistrib) { + /* Double predistrib factor in case no photons were stored + * for this source and redo pass 1 */ + srcPreDistrib *= 2; + } + else { + /* Now do pass 2 */ + passCnt++; + } } - } while (passCnt < 2); - - /* Flux per photon emitted from this source; repProgress is the - * number of emitted photons after both passes */ - srcFlux [srcIdx] = repProgress ? srcFlux [srcIdx] / repProgress - : 0; + } + + /* Flush heap buffa one final time to prevent data corruption */ + flushPhotonHeap(pm); + /* Flush final photon primary to primary heap file */ + newPhotonPrimary(pm, NULL, primaryHeap [proc]); + /* Heap files closed automatically on exit + fclose(pm -> heap); + fclose(primaryHeap [proc]); */ + +#ifdef DEBUG_PMAP + sprintf(errmsg, "Proc %d total %ld photons\n", proc, + pm -> numPhotons); + eputs(errmsg); + fflush(stderr); +#endif + +#ifdef PMAP_SIGUSR + signal(SIGUSR1, SIG_DFL); +#endif + +#if NIX + /* Terminate subprocess */ + exit(0); +#endif } + else if (pid < 0) + error(SYSTEM, "failed to fork subprocess in distribPhotonContrib"); } +#if NIX + /* PARENT PROCESS CONTINUES HERE */ +#ifdef SIGCONT + /* Enable progress report signal handler */ + signal(SIGCONT, pmapDistribReport); +#endif + /* Wait for subprocesses to complete while reporting progress */ + proc = numProc; + while (proc) { + while (waitpid(-1, &stat, WNOHANG) > 0) { + /* Subprocess exited; check status */ + if (!WIFEXITED(stat) || WEXITSTATUS(stat)) + error(USER, "failed photon distribution"); + + --proc; + } + + /* Nod off for a bit and update progress */ + sleep(1); + + /* Asynchronous progress report from shared subprocess counters */ + repComplete = pm -> distribTarget * numProc; + repProgress = photonCnt [PHOTONCNT_NUMPHOT]; + + for (repEmitted = 0, srcIdx = 0; srcIdx < nsources; srcIdx++) + repEmitted += photonCnt [PHOTONCNT_NUMEMIT(srcIdx)]; + + /* Get global photon count from shmem updated by subprocs */ + pm -> numPhotons = photonCnt [PHOTONCNT_NUMPHOT]; + + if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime) + pmapDistribReport(); +#ifdef SIGCONT + else signal(SIGCONT, pmapDistribReport); +#endif + } +#endif /* NIX */ + /* ================================================================ * POST-DISTRIBUTION - Set photon flux and build kd-tree, etc. * ================================================================ */ +#ifdef SIGCONT + /* Reset signal handler */ signal(SIGCONT, SIG_DFL); +#endif free(emap.samples); - if (!pm -> heapEnd) - error(USER, "empty photon map"); + if (!pm -> numPhotons) + error(USER, "empty contribution photon map"); - /* Check for valid primary photon rays */ - if (!pm -> primary) + /* Load per-subprocess primary rays into pm -> primary array */ + /* Dumb compilers apparently need the char** cast */ + pm -> numPrimary = buildPrimaries(pm, primaryHeap, + (char**)primaryHeapFname, + primaryOfs, numProc); + if (!pm -> numPrimary) error(INTERNAL, "no primary rays in contribution photon map"); - - if (pm -> primary [pm -> primaryEnd].srcIdx < 0) - /* Last primary ray is unused, so decrement counter */ - pm -> primaryEnd--; - if (photonRepTime) { - eputs("\nBuilding contrib photon heap...\n"); + /* Set photon flux per source */ + for (srcIdx = 0; srcIdx < nsources; srcIdx++) + srcFlux [srcIdx] /= photonCnt [PHOTONCNT_NUMEMIT(srcIdx)]; +#if NIX + /* Photon counters no longer needed, unmap shared memory */ + munmap(photonCnt, sizeof(*photonCnt)); + close(shmFile); + unlink(shmFname); +#else + free(photonCnt); +#endif + + if (verbose) { + eputs("\nBuilding contribution photon map...\n"); +#if NIX fflush(stderr); +#endif } + + /* Build underlying data structure; heap is destroyed */ + buildPhotonMap(pm, srcFlux, primaryOfs, numProc); + + /* Free per-subprocess primary heap files */ + for (proc = 0; proc < numProc; proc++) + free(primaryHeapFname [proc]); - balancePhotons(pm, srcFlux); + free(primaryHeapFname); + free(primaryHeap); + free(primaryOfs); + + if (verbose) + eputs("\n"); }