--- ray/src/rt/pmap.c 2015/04/22 20:28:16 2.3 +++ ray/src/rt/pmap.c 2018/03/20 19:55:33 2.14 @@ -1,17 +1,22 @@ +#ifndef lint +static const char RCSid[] = "$Id: pmap.c,v 2.14 2018/03/20 19:55:33 rschregle Exp $"; +#endif + + /* - ================================================================== + ====================================================================== Photon map main module 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: pmap.c,v 2.3 2015/04/22 20:28:16 rschregle Exp $ + $Id: pmap.c,v 2.14 2018/03/20 19:55:33 rschregle Exp $ */ - #include "pmap.h" #include "pmapmat.h" #include "pmapsrc.h" @@ -21,103 +26,24 @@ #include "pmapdiag.h" #include "otypes.h" #include -#include +#if NIX + #include + #include + #include +#endif - -extern char *octname; - -static char PmapRevision [] = "$Revision: 2.3 $"; - - - -/* Photon map lookup functions per type */ -void (*pmapLookup [NUM_PMAP_TYPES])(PhotonMap*, RAY*, COLOR) = { - photonDensity, photonPreCompDensity, photonDensity, volumePhotonDensity, - photonDensity, NULL -}; - - - -void colorNorm (COLOR c) -/* Normalise colour channels to average of 1 */ -{ - const float avg = colorAvg(c); - - if (!avg) - return; - - c [0] /= avg; - c [1] /= avg; - c [2] /= avg; -} - - - -void loadPmaps (PhotonMap **pmaps, const PhotonMapParams *parm) -{ - unsigned t; - struct stat octstat, pmstat; - PhotonMap *pm; - PhotonMapType type; - - for (t = 0; t < NUM_PMAP_TYPES; t++) - if (setPmapParam(&pm, parm + t)) { - /* Check if photon map newer than octree */ - if (!stat(pm -> fileName, &pmstat) && !stat(octname, &octstat) && - octstat.st_mtime > pmstat.st_mtime) { - sprintf(errmsg, "photon map in file %s may be stale", - pm -> fileName); - error(USER, errmsg); - } - - /* Load photon map from file and get its type */ - if ((type = loadPhotonMap(pm, pm -> fileName)) == PMAP_TYPE_NONE) - error(USER, "failed loading photon map"); - - /* Assign to appropriate photon map type (deleting previously - * loaded photon map of same type if necessary) */ - if (pmaps [type]) { - deletePhotons(pmaps [type]); - free(pmaps [type]); - } - pmaps [type] = pm; - - /* Check for invalid density estimate bandwidth */ - if (pm -> maxGather > pm -> heapSize) { - error(WARNING, "adjusting density estimate bandwidth"); - pm -> minGather = pm -> maxGather = pm -> heapSize; - } - } -} - - - void savePmaps (const PhotonMap **pmaps, int argc, char **argv) { unsigned t; for (t = 0; t < NUM_PMAP_TYPES; t++) { if (pmaps [t]) - savePhotonMap(pmaps [t], pmaps [t] -> fileName, t, argc, argv); + savePhotonMap(pmaps [t], pmaps [t] -> fileName, argc, argv); } } - -void cleanUpPmaps (PhotonMap **pmaps) -{ - unsigned t; - - for (t = 0; t < NUM_PMAP_TYPES; t++) { - if (pmaps [t]) { - deletePhotons(pmaps [t]); - free(pmaps [t]); - } - } -} - - static int photonParticipate (RAY *ray) /* Trace photon through participating medium. Returns 1 if passed through, @@ -143,12 +69,13 @@ static int photonParticipate (RAY *ray) colorNorm(ray -> rcol); VCOPY(ray -> rorg, ray -> rop); - if (albedo > FTINY) + if (albedo > FTINY && ray -> rlvl > 0) /* Add to volume photon map */ - if (ray -> rlvl > 0) addPhoton(volumePmap, ray); + newPhoton(volumePmap, ray); /* Absorbed? */ - if (pmapRandom(rouletteState) > albedo) return 0; + if (pmapRandom(rouletteState) > albedo) + return 0; /* Colour bleeding without attenuation (?) */ multcolor(ray -> rcol, ray -> albedo); @@ -201,19 +128,45 @@ 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 error(WARNING, "runaway photon!"); +#endif return; } - + if (colorAvg(ray -> cext) > FTINY && !photonParticipate(ray)) return; - + if (localhit(ray, &thescene)) { mod = ray -> ro -> omod; - + + if (port && ray -> ro != port) { + /* !!! 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 + * 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; @@ -231,101 +184,129 @@ void tracePhoton (RAY *ray) static void preComputeGlobal (PhotonMap *pmap) -/* Precompute irradiance from global photons for final gathering using - the first finalGather * pmap -> heapSize photons in the heap. Returns - new heap with precomputed photons. */ +/* Precompute irradiance from global photons for final gathering for + a random subset of finalGather * pmap -> numPhotons photons, and builds + the photon map, discarding the original photons. */ +/* !!! NOTE: PRECOMPUTATION WITH OOC CURRENTLY WITHOUT CACHE !!! */ { - unsigned long i, nuHeapSize; - unsigned j; - Photon *nuHeap, *p; - COLOR irrad; - RAY ray; - float nuMinPos [3], nuMaxPos [3]; + unsigned long i, numPreComp; + unsigned j; + PhotonIdx pIdx; + Photon photon; + RAY ray; + PhotonMap nuPmap; - repComplete = nuHeapSize = finalGather * pmap -> heapSize; + repComplete = numPreComp = finalGather * pmap -> numPhotons; - if (photonRepTime) { + if (verbose) { sprintf(errmsg, - "Precomputing irradiance for %ld global photons...\n", - nuHeapSize); + "\nPrecomputing irradiance for %ld global photons\n", + numPreComp); eputs(errmsg); +#if NIX fflush(stderr); +#endif } - p = nuHeap = (Photon*)malloc(nuHeapSize * sizeof(Photon)); - if (!nuHeap) - error(USER, "can't allocate photon heap"); - - for (j = 0; j <= 2; j++) { - nuMinPos [j] = FHUGE; - nuMaxPos [j] = -FHUGE; - } + /* Copy photon map for precomputed photons */ + memcpy(&nuPmap, pmap, sizeof(PhotonMap)); + + /* Zero counters, init new heap and extents */ + nuPmap.numPhotons = 0; + initPhotonHeap(&nuPmap); + for (j = 0; j < 3; j++) { + nuPmap.minPos [j] = FHUGE; + nuPmap.maxPos [j] = -FHUGE; + } + /* Record start time, baby */ repStartTime = time(NULL); - #ifdef SIGCONT - signal(SIGCONT, pmapPreCompReport); - #endif +#ifdef SIGCONT + signal(SIGCONT, pmapPreCompReport); +#endif repProgress = 0; - memcpy(nuHeap, pmap -> heap, nuHeapSize * sizeof(Photon)); - for (i = 0, p = nuHeap; i < nuHeapSize; i++, p++) { - ray.ro = NULL; - VCOPY(ray.rop, p -> pos); + photonRay(NULL, &ray, PRIMARY, NULL); + ray.ro = NULL; + + for (i = 0; i < numPreComp; i++) { + /* Get random photon from stratified distribution in source heap to + * avoid duplicates and clustering */ + pIdx = firstPhoton(pmap) + + (unsigned long)((i + pmapRandom(pmap -> randState)) / + finalGather); + getPhoton(pmap, pIdx, &photon); - /* Update min and max positions & set ray normal */ - for (j = 0; j < 3; j++) { - if (p -> pos [j] < nuMinPos [j]) nuMinPos [j] = p -> pos [j]; - if (p -> pos [j] > nuMaxPos [j]) nuMaxPos [j] = p -> pos [j]; - ray.ron [j] = p -> norm [j] / 127.0; - } + /* Init dummy photon ray with intersection at photon position */ + VCOPY(ray.rop, photon.pos); + for (j = 0; j < 3; j++) + ray.ron [j] = photon.norm [j] / 127.0; - photonDensity(pmap, &ray, irrad); - setPhotonFlux(p, irrad); + /* Get density estimate at photon position */ + photonDensity(pmap, &ray, ray.rcol); + + /* Append photon to new heap from ray */ + newPhoton(&nuPmap, &ray); + + /* Update progress */ repProgress++; if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime) pmapPreCompReport(); - #ifdef SIGCONT - else signal(SIGCONT, pmapPreCompReport); - #endif +#ifdef SIGCONT + else signal(SIGCONT, pmapPreCompReport); +#endif } - #ifdef SIGCONT - signal(SIGCONT, SIG_DFL); - #endif + /* Flush heap */ + flushPhotonHeap(&nuPmap); - /* Replace & rebuild heap */ - free(pmap -> heap); - pmap -> heap = nuHeap; - pmap -> heapSize = pmap -> heapEnd = nuHeapSize; - VCOPY(pmap -> minPos, nuMinPos); - VCOPY(pmap -> maxPos, nuMaxPos); +#ifdef SIGCONT + signal(SIGCONT, SIG_DFL); +#endif - if (photonRepTime) { - eputs("Rebuilding global photon heap...\n"); + /* Trash original pmap, replace with precomputed one */ + deletePhotons(pmap); + memcpy(pmap, &nuPmap, sizeof(PhotonMap)); + + if (verbose) { + eputs("\nRebuilding precomputed photon map\n"); +#if NIX fflush(stderr); +#endif } - - balancePhotons(pmap, NULL); + + /* Rebuild underlying data structure, destroying heap */ + buildPhotonMap(pmap, NULL, NULL, 1); } -void distribPhotons (PhotonMap **pmaps) +typedef struct { + unsigned long numPhotons [NUM_PMAP_TYPES], + numEmitted, numComplete; +} PhotonCnt; + + + +void distribPhotons (PhotonMap **pmaps, unsigned numProc) { - EmissionMap emap; - char errmsg2 [128]; - unsigned t, srcIdx, passCnt = 0, prePassCnt = 0; - double totalFlux = 0; - PhotonMap *pm; + EmissionMap emap; + char errmsg2 [128], shmFname [PMAP_TMPFNLEN]; + unsigned t, srcIdx, proc; + double totalFlux = 0; + int shmFile, stat, pid; + PhotonMap *pm; + PhotonCnt *photonCnt; - for (t = 0; t < NUM_PMAP_TYPES && !photonMaps [t]; t++); + for (t = 0; t < NUM_PMAP_TYPES && !pmaps [t]; t++); + if (t >= NUM_PMAP_TYPES) - error(USER, "no photon maps defined"); + error(USER, "no photon maps defined in distribPhotons"); if (!nsources) - error(USER, "no light sources"); + error(USER, "no light sources in distribPhotons"); /* =================================================================== * INITIALISATION - Set up emission and scattering funcs @@ -334,59 +315,90 @@ void distribPhotons (PhotonMap **pmaps) emap.maxPartitions = MAXSPART; emap.partitions = (unsigned char*)malloc(emap.maxPartitions >> 1); if (!emap.partitions) - error(INTERNAL, "can't allocate source partitions"); + error(INTERNAL, "can't allocate source partitions in distribPhotons"); /* Initialise all defined photon maps */ for (t = 0; t < NUM_PMAP_TYPES; t++) - initPhotonMap(photonMaps [t], t); + if (pmaps [t]) { + initPhotonMap(pmaps [t], t); + /* Open photon heapfile */ + 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); - /* Seed RNGs for photon distribution */ - pmapSeed(randSeed, partState); - pmapSeed(randSeed, emitState); - pmapSeed(randSeed, cntState); - pmapSeed(randSeed, mediumState); - pmapSeed(randSeed, scatterState); - pmapSeed(randSeed, rouletteState); +#if NIX + /* Set up shared mem for photon counters (zeroed by ftruncate) */ + strcpy(shmFname, PMAP_TMPFNAME); + shmFile = mkstemp(shmFname); + + if (shmFile < 0 || ftruncate(shmFile, sizeof(*photonCnt)) < 0) + error(SYSTEM, "failed shared mem init 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"); +#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 (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); + } - if (photonRepTime) - eputs("\n"); - /* =================================================================== * FLUX INTEGRATION - Get total photon flux from light sources * =================================================================== */ - for (srcIdx = 0; srcIdx < nsources; srcIdx++) { + for (srcIdx = 0; srcIdx < nsources; srcIdx++) { unsigned portCnt = 0; emap.src = source + srcIdx; - do { + 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 ", + 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; @@ -401,361 +413,354 @@ void distribPhotons (PhotonMap **pmaps) if (totalFlux < FTINY) error(USER, "zero flux from light sources"); - - /* Record start time and enable progress report signal handler */ - repStartTime = time(NULL); - #ifdef SIGCONT - signal(SIGCONT, pmapDistribReport); - #endif - repProgress = prePassCnt = 0; - - if (photonRepTime) - eputs("\n"); - - /* =================================================================== - * 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 { - double numEmit; - if (!passCnt) { - /* INIT PASS 1 */ - /* Skip if no photons contributed after sufficient iterations; make - * it clear to user which photon maps are missing so (s)he can - * check the scene geometry and materials */ - if (++prePassCnt > maxPreDistrib) { - sprintf(errmsg, "too many prepasses"); + /* Record start time for progress reports */ + repStartTime = time(NULL); - for (t = 0; t < NUM_PMAP_TYPES; t++) - if (photonMaps [t] && !photonMaps [t] -> heapEnd) { - sprintf(errmsg2, ", no %s photons stored", pmapName [t]); - strcat(errmsg, errmsg2); - } - - error(USER, errmsg); - break; - } + if (verbose) { + sprintf(errmsg, "\nPhoton distribution @ %d procs\n", numProc); + eputs(errmsg); + } - /* Num to emit is fraction of minimum target count */ - numEmit = FHUGE; + /* MAIN LOOP */ + for (proc = 0; proc < numProc; proc++) { +#if NIX + if (!(pid = fork())) { + /* 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 + subprocess alone */ - for (t = 0; t < NUM_PMAP_TYPES; t++) - if (photonMaps [t]) - numEmit = min(photonMaps [t] -> distribTarget, numEmit); - - numEmit *= preDistrib; - } - - else { - /* INIT PASS 2 */ - /* Based on the outcome of the predistribution we can now estimate - * how many more photons we have to emit for each photon map to - * meet its respective target count. This value is clamped to 0 in - * case the target has already been exceeded in the pass 1. Note - * repProgress is the number of photons emitted thus far, while - * heapEnd is the number of photons stored in each photon map. */ - double maxDistribRatio = 0; - - /* Set the distribution ratio for each map; this indicates how many - * photons of each respective type are stored per emitted photon, - * and is used as probability for storing a photon by addPhoton(). - * Since this biases the photon density, addPhoton() promotes the - * flux of stored photons to compensate. */ - for (t = 0; t < NUM_PMAP_TYPES; t++) - if ((pm = photonMaps [t])) { - pm -> distribRatio = (double)pm -> distribTarget / - pm -> heapEnd - 1; - - /* Check if photon map "overflowed", i.e. exceeded its target - * count in the prepass; correcting the photon flux via the - * distribution ratio is no longer possible, as no more - * photons of this type will be stored, so notify the user - * rather than deliver incorrect results. - * In future we should handle this more intelligently by - * using the photonFlux in each photon map to individually - * correct the flux after distribution. */ - if (pm -> distribRatio <= FTINY) { - sprintf(errmsg, - "%s photon map overflow in prepass, reduce -apD", - pmapName [t]); - error(INTERNAL, errmsg); - } + /* 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); - maxDistribRatio = max(pm -> distribRatio, maxDistribRatio); - } - - /* Normalise distribution ratios and calculate number of photons to - * emit in main pass */ for (t = 0; t < NUM_PMAP_TYPES; t++) - if ((pm = photonMaps [t])) - pm -> distribRatio /= maxDistribRatio; - - if ((numEmit = repProgress * maxDistribRatio) < FTINY) - /* No photons left to distribute in main pass */ - break; - } - - /* Set completion count for progress report */ - repComplete = numEmit + repProgress; - - /* PHOTON DISTRIBUTION LOOP */ - for (srcIdx = 0; srcIdx < nsources; srcIdx++) { - unsigned portCnt = 0; - emap.src = source + srcIdx; - + lastNumPhotons [t] = 0; + + /* ============================================================= + * 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 { - emap.port = emap.src -> sflags & SDISTANT ? photonPorts + portCnt - : NULL; - photonPartition [emap.src -> so -> otype] (&emap); + double numEmit; - if (photonRepTime) { - if (!passCnt) - sprintf(errmsg, "PREPASS %d on source %s ", - prePassCnt, source [srcIdx].so -> oname); - else - sprintf(errmsg, "MAIN PASS on source %s ", - source [srcIdx].so -> oname); - - if (emap.port) { - sprintf(errmsg2, "via port %s ", - photonPorts [portCnt].so -> oname); - strcat(errmsg, errmsg2); + if (!passCnt) { + /* INIT PASS 1 */ + /* Skip if no photons contributed after sufficient + * 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: too many prepasses", proc); + + for (t = 0; t < NUM_PMAP_TYPES; t++) + if (pmaps [t] && !pmaps [t] -> numPhotons) { + sprintf(errmsg2, ", no %s photons stored", + pmapName [t]); + strcat(errmsg, errmsg2); + } + + error(USER, errmsg); + break; } + + /* Num to emit is fraction of minimum target count */ + numEmit = FHUGE; - sprintf(errmsg2, "(%lu partitions)...\n", emap.numPartitions); - strcat(errmsg, errmsg2); - eputs(errmsg); - fflush(stderr); + for (t = 0; t < NUM_PMAP_TYPES; t++) + if (pmaps [t]) + numEmit = min(pmaps [t] -> distribTarget, numEmit); + + numEmit *= preDistrib; } - - for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions; - emap.partitionCnt++) { - double partNumEmit; - unsigned long partEmitCnt; + else { + /* INIT PASS 2 */ + /* Based on the outcome of the predistribution we can now + * estimate how many more photons we have to emit for each + * photon map to meet its respective target count. This + * value is clamped to 0 in case the target has already been + * exceeded in the pass 1. */ + double maxDistribRatio = 0; + + /* Set the distribution ratio for each map; this indicates + * how many photons of each respective type are stored per + * emitted photon, and is used as probability for storing a + * photon by newPhoton(). Since this biases the photon + * density, newPhoton() promotes the flux of stored photons + * to compensate. */ + for (t = 0; t < NUM_PMAP_TYPES; t++) + if ((pm = pmaps [t])) { + pm -> distribRatio = (double)pm -> distribTarget / + pm -> numPhotons - 1; + + /* Check if photon map "overflowed", i.e. exceeded its + * target count in the prepass; correcting the photon + * flux via the distribution ratio is no longer + * possible, as no more photons of this type will be + * stored, so notify the user rather than deliver + * incorrect results. In future we should handle this + * more intelligently by using the photonFlux in each + * photon map to individually correct the flux after + * distribution. */ + if (pm -> distribRatio <= FTINY) { + sprintf(errmsg, "%s photon map overflow in " + "prepass, reduce -apD", pmapName [t]); + error(INTERNAL, errmsg); + } + + maxDistribRatio = max(pm -> distribRatio, + maxDistribRatio); + } - /* 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 -- - * proportional to flux; photon ray weight and scalar flux - * are uniform (the latter only varying in RGB). */ - partNumEmit = numEmit * colorAvg(emap.partFlux) / totalFlux; - partEmitCnt = (unsigned long)partNumEmit; - - /* Probabilistically account for fractional photons */ - if (pmapRandom(cntState) < partNumEmit - partEmitCnt) - partEmitCnt++; + /* Normalise distribution ratios and calculate number of + * photons to emit in main pass */ + for (t = 0; t < NUM_PMAP_TYPES; t++) + if ((pm = pmaps [t])) + pm -> distribRatio /= maxDistribRatio; + + if ((numEmit = localNumEmitted * maxDistribRatio) < FTINY) + /* No photons left to distribute in main pass */ + break; + } - /* Integer counter avoids FP rounding errors */ - while (partEmitCnt--) { - RAY photonRay; + /* Update shared completion counter for progreport by parent */ + photonCnt -> numComplete += numEmit; + + /* PHOTON DISTRIBUTION LOOP */ + for (srcIdx = 0; srcIdx < nsources; srcIdx++) { + unsigned portCnt = 0; + emap.src = source + srcIdx; + + 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 (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 + } - /* Emit photon based on PDF and trace through scene until - * absorbed/leaked */ - emitPhoton(&emap, &photonRay); - tracePhoton(&photonRay); + for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions; + emap.partitionCnt++) { + double partNumEmit; + unsigned long partEmitCnt; + + /* 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 -- + * proportional to flux; photon ray weight and scalar + * flux are uniform (latter only varying in RGB). */ + partNumEmit = numEmit * colorAvg(emap.partFlux) / + totalFlux; + partEmitCnt = (unsigned long)partNumEmit; + + /* Probabilistically account for fractional photons */ + if (pmapRandom(cntState) < partNumEmit - partEmitCnt) + partEmitCnt++; + + /* Update local and shared (global) emission counter */ + photonCnt -> numEmitted += partEmitCnt; + localNumEmitted += partEmitCnt; + + /* Integer counter avoids FP rounding errors during + * iteration */ + while (partEmitCnt--) { + RAY photonRay; + + /* 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]) { + photonCnt -> numPhotons [t] += + 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 + } - /* Record progress */ - repProgress++; - - if (photonRepTime > 0 && - time(NULL) >= repLastTime + photonRepTime) - pmapDistribReport(); - #ifdef SIGCONT - else signal(SIGCONT, pmapDistribReport); - #endif + portCnt++; + } while (portCnt < numPhotonPorts); + } + + for (t = 0; t < NUM_PMAP_TYPES; t++) + if (pmaps [t] && !pmaps [t] -> numPhotons) { + /* Double preDistrib in case a photon map is empty and + * redo pass 1 --> possibility of infinite loop for + * pathological scenes (e.g. absorbing materials) */ + preDistrib *= 2; + break; } + + if (t >= NUM_PMAP_TYPES) + /* No empty photon maps found; now do pass 2 */ + passCnt++; + } while (passCnt < 2); + + /* Flush heap buffa for every pmap one final time; + * avoids potential data corruption! */ + for (t = 0; t < NUM_PMAP_TYPES; t++) + if (pmaps [t]) { + flushPhotonHeap(pmaps [t]); + /* Heap file closed automatically on exit + fclose(pmaps [t] -> heap); */ +#ifdef DEBUG_PMAP + sprintf(errmsg, "Proc %d: total %ld photons\n", proc, + pmaps [t] -> numPhotons); + eputs(errmsg); +#endif } - - portCnt++; - } while (portCnt < numPhotonPorts); +#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 */ +#ifdef SIGCONT + /* Enable progress report signal handler */ + signal(SIGCONT, pmapDistribReport); +#endif + /* Wait for subprocesses 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 */ + repEmitted = repProgress = photonCnt -> numEmitted; + repComplete = photonCnt -> numComplete; + + repProgress = repComplete = 0; for (t = 0; t < NUM_PMAP_TYPES; t++) - if (photonMaps [t] && !photonMaps [t] -> heapEnd) { - /* Double preDistrib in case a photon map is empty and redo - * pass 1 --> possibility of infinite loop for pathological - * scenes (e.g. absorbing materials) */ - preDistrib *= 2; - break; + if ((pm = pmaps [t])) { + /* Get global photon count from shmem updated by subprocs */ + repProgress += pm -> numPhotons = photonCnt -> numPhotons [t]; + repComplete += pm -> distribTarget; } - - if (t >= NUM_PMAP_TYPES) { - /* No empty photon maps found; now do pass 2 */ - passCnt++; - if (photonRepTime) - eputs("\n"); - } - } while (passCnt < 2); + repComplete *= numProc; + 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. + * POST-DISTRIBUTION - Set photon flux and build data struct for photon + * storage, etc. * =================================================================== */ - #ifdef SIGCONT - signal(SIGCONT, SIG_DFL); - #endif +#ifdef SIGCONT + /* Reset signal handler */ + signal(SIGCONT, SIG_DFL); +#endif free(emap.samples); - /* Set photon flux (repProgress is total num emitted) */ - totalFlux /= repProgress; - + /* Set photon flux */ + totalFlux /= photonCnt -> numEmitted; +#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("\n"); + for (t = 0; t < NUM_PMAP_TYPES; t++) - if (photonMaps [t]) { - if (photonRepTime) { - sprintf(errmsg, "\nBuilding %s photon map...\n", pmapName [t]); + if (pmaps [t]) { + if (verbose) { + sprintf(errmsg, "Building %s photon map\n", pmapName [t]); eputs(errmsg); +#if NIX fflush(stderr); +#endif } - - balancePhotons(photonMaps [t], &totalFlux); + + /* 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); -} - - - -void photonDensity (PhotonMap *pmap, RAY *ray, COLOR irrad) -/* Photon density estimate. Returns irradiance at ray -> rop. */ -{ - unsigned i; - PhotonSQNode *sq; - float r; - COLOR flux; - - setcolor(irrad, 0, 0, 0); - - if (!pmap -> maxGather) - return; - - /* Ignore sources */ - if (ray -> ro) - if (islight(objptr(ray -> ro -> omod) -> otype)) - return; - - pmap -> squeueEnd = 0; - findPhotons(pmap, ray); + } - /* 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 - - return; - } - - if (pmap -> minGather == pmap -> maxGather) { - /* No bias compensation. Just do a plain vanilla estimate */ - sq = pmap -> squeue + 1; - - /* Average radius between furthest two photons to improve accuracy */ - r = max(sq -> dist, (sq + 1) -> dist); - r = 0.25 * (pmap -> maxDist + r + 2 * sqrt(pmap -> maxDist * r)); - - /* Skip the extra photon */ - for (i = 1 ; i < pmap -> squeueEnd; i++, sq++) { - getPhotonFlux(sq -> photon, flux); -#ifdef PMAP_EPANECHNIKOV - /* Apply Epanechnikov kernel to photon flux (dists are squared) */ - scalecolor(flux, 2 * (1 - sq -> dist / r)); -#endif - addcolor(irrad, flux); - } - - /* Divide by search area PI * r^2, 1 / PI required as ambient - normalisation factor */ - scalecolor(irrad, 1 / (PI * PI * r)); - - return; - } - else - /* Apply bias compensation to density estimate */ - biasComp(pmap, irrad); -} - - - -void photonPreCompDensity (PhotonMap *pmap, RAY *r, COLOR irrad) -/* Returns precomputed photon density estimate at ray -> rop. */ -{ - Photon *p; - - setcolor(irrad, 0, 0, 0); - - /* Ignore sources */ - if (r -> ro && islight(objptr(r -> ro -> omod) -> otype)) - return; - - if ((p = find1Photon(preCompPmap, r))) - getPhotonFlux(p, irrad); -} - - - -void volumePhotonDensity (PhotonMap *pmap, RAY *ray, COLOR irrad) -/* Photon volume density estimate. Returns irradiance at ray -> rop. */ -{ - unsigned i; - PhotonSQNode *sq; - float gecc2, r, ph; - COLOR flux; - - setcolor(irrad, 0, 0, 0); - - if (!pmap -> maxGather) - return; - - pmap -> squeueEnd = 0; - findPhotons(pmap, ray); - - /* Need at least 2 photons */ - if (pmap -> squeueEnd < 2) - return; - - if (pmap -> minGather == pmap -> maxGather) { - /* No bias compensation. Just do a plain vanilla estimate */ - gecc2 = ray -> gecc * ray -> gecc; - sq = pmap -> squeue + 1; - - /* Average radius between furthest two photons to improve accuracy */ - r = max(sq -> dist, (sq + 1) -> dist); - r = 0.25 * (pmap -> maxDist + r + 2 * sqrt(pmap -> maxDist * r)); - - /* Skip the extra photon */ - for (i = 1 ; i < pmap -> squeueEnd; i++, sq++) { - /* Compute phase function for inscattering from photon */ - if (gecc2 <= FTINY) - ph = 1; - else { - ph = DOT(ray -> rdir, sq -> photon -> norm) / 127; - ph = 1 + gecc2 - 2 * ray -> gecc * ph; - ph = (1 - gecc2) / (ph * sqrt(ph)); - } - - getPhotonFlux(sq -> photon, flux); - scalecolor(flux, ph); - addcolor(irrad, flux); - } - - /* Divide by search volume 4 / 3 * PI * r^3 and phase function - normalization factor 1 / (4 * PI) */ - scalecolor(irrad, 3 / (16 * PI * PI * r * sqrt(r))); - - return; - } - - else - /* Apply bias compensation to density estimate */ - volumeBiasComp(pmap, ray, irrad); + if (verbose) + eputs("\n"); }