--- ray/src/rt/pmap.c 2015/09/01 16:27:52 2.10 +++ ray/src/rt/pmap.c 2016/09/26 20:19:30 2.12 @@ -1,16 +1,18 @@ #ifndef lint -static const char RCSid[] = "$Id: pmap.c,v 2.10 2015/09/01 16:27:52 greg Exp $"; +static const char RCSid[] = "$Id: pmap.c,v 2.12 2016/09/26 20:19:30 greg Exp $"; #endif + /* - ================================================================== + ====================================================================== Photon map main module Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) (c) Fraunhofer Institute for Solar Energy Systems, (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) + ====================================================================== + $Id: pmap.c,v 2.12 2016/09/26 20:19:30 greg Exp $ */ @@ -25,78 +27,11 @@ static const char RCSid[] = "$Id: pmap.c,v 2.10 2015/0 #include "otypes.h" #include #include +#include +#include -extern char *octname; - -static char PmapRevision [] = "$Revision: 2.10 $"; - - - -/* 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 (pm -> fileName && octname && - !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; @@ -108,20 +43,6 @@ void savePmaps (const PhotonMap **pmaps, int argc, cha } - -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, @@ -147,12 +68,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); @@ -237,101 +159,124 @@ 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) { - sprintf(errmsg, - "Precomputing irradiance for %ld global photons...\n", - nuHeapSize); + sprintf(errmsg, "Precomputing irradiance for %ld global photons...\n", + numPreComp); eputs(errmsg); fflush(stderr); } - 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 clutering */ + 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 + /* Trash original pmap, replace with precomputed one */ + deletePhotons(pmap); + memcpy(pmap, &nuPmap, sizeof(PhotonMap)); + if (photonRepTime) { - eputs("Rebuilding global photon heap...\n"); + eputs("Rebuilding precomputed photon map...\n"); fflush(stderr); } - - 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 [255]; + unsigned t, srcIdx, proc; + double totalFlux = 0; + int shmFile, stat, pid; + PhotonMap *pm; + PhotonCnt *photonCnt; 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 @@ -340,11 +285,17 @@ 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(pmaps [t], t); + if (pmaps [t]) { + initPhotonMap(pmaps [t], t); + /* Open photon heapfile */ + initPhotonHeap(pmaps [t]); + /* Per-subprocess target count */ + pmaps [t] -> distribTarget /= numProc; + } initPhotonEmissionFuncs(); initPhotonScatterFuncs(); @@ -356,25 +307,38 @@ void distribPhotons (PhotonMap **pmaps) /* 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); - + /* 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); + shmFile = mkstemp(shmFname); +#endif + + if (shmFile < 0) + error(SYSTEM, "failed opening shared memory file 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"); + 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); @@ -408,208 +372,315 @@ 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"); - - for (t = 0; t < NUM_PMAP_TYPES; t++) - if (pmaps [t] && !pmaps [t] -> heapEnd) { - 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; + /* MAIN LOOP */ + for (proc = 0; proc < numProc; proc++) { + if (!(pid = fork())) { + /* SUBPROCESS ENTERS HERE. + All opened and memory mapped files are inherited */ + 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 (pmaps [t]) - numEmit = min(pmaps [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 = pmaps [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, emitState); + pmapSeed(randSeed + proc, cntState); + pmapSeed(randSeed + proc, mediumState); + pmapSeed(randSeed + proc, scatterState); + pmapSeed(randSeed + proc, 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 = pmaps [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, source %s: too many prepasses", + proc, source [srcIdx].so -> oname); + + 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 prog.report 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 (photonRepTime && !proc) { + 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); + } + + sprintf(errmsg2, "(%lu partitions)\n", + emap.numPartitions); + strcat(errmsg, errmsg2); + eputs(errmsg); + fflush(stderr); + } - /* 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 (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++; + + /* 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); + 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; + } + } - /* 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++; +#if 0 + if (photonRepTime) + eputs("\n"); +#endif } - - portCnt++; - } while (portCnt < numPhotonPorts); + } 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! */ + 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); +#ifdef DEBUG_PMAP + sprintf(errmsg, "Proc %d: total %ld photons\n", getpid(), + pmaps [t] -> numPhotons); + eputs(errmsg); +#endif + } + + exit(0); } + else if (pid < 0) + error(SYSTEM, "failed to fork subprocess in distribPhotons"); + } + + /* PARENT PROCESS CONTINUES HERE */ + /* Record start time and enable progress report signal handler */ + repStartTime = time(NULL); +#ifdef SIGCONT + signal(SIGCONT, pmapDistribReport); +#endif + + if (photonRepTime) + eputs("\n"); + + /* 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); + /* Update progress report from shared subprocess counters */ + repEmitted = repProgress = photonCnt -> numEmitted; + repComplete = photonCnt -> numComplete; + for (t = 0; t < NUM_PMAP_TYPES; t++) - if (pmaps [t] && !pmaps [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])) { +#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 } - - if (t >= NUM_PMAP_TYPES) { - /* No empty photon maps found; now do pass 2 */ - passCnt++; - if (photonRepTime) - eputs("\n"); - } - } while (passCnt < 2); + if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime) + pmapDistribReport(); +#ifdef SIGCONT + else signal(SIGCONT, pmapDistribReport); +#endif + } + /* =================================================================== - * 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 + signal(SIGCONT, SIG_DFL); +#endif free(emap.samples); /* Set photon flux (repProgress is total num emitted) */ - totalFlux /= repProgress; + totalFlux /= photonCnt -> numEmitted; + /* Photon counters no longer needed, unmap shared memory */ + munmap(photonCnt, sizeof(*photonCnt)); + close(shmFile); +#if 0 + shm_unlink(shmFname); +#else + unlink(shmFname); +#endif + for (t = 0; t < NUM_PMAP_TYPES; t++) if (pmaps [t]) { if (photonRepTime) { @@ -617,151 +688,12 @@ void distribPhotons (PhotonMap **pmaps) eputs(errmsg); fflush(stderr); } - - balancePhotons(pmaps [t], &totalFlux); + + /* Build underlying data structure; heap is destroyed */ + buildPhotonMap(pmaps [t], &totalFlux, NULL, numProc); } - + /* Precompute photon irradiance if necessary */ if (preCompPmap) 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); }