diff --git a/pmap.c b/pmap.c index 025b19c..9b55883 100644 --- a/pmap.c +++ b/pmap.c @@ -1,933 +1,947 @@ #ifndef lint static const char RCSid[] = "$Id: pmap.c,v 2.18 2021/02/19 02:10:35 rschregle Exp $"; #endif /* ====================================================================== Photon map main module Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) (c) Fraunhofer Institute for Solar Energy Systems, supported by the German Research Foundation (DFG LU-204/10-2, "Fassadenintegrierte Regelsysteme FARESYS") (c) Lucerne University of Applied Sciences and Arts, supported by the Swiss National Science Foundation (SNSF #147053, "Daylight Redirecting Components") (c) Tokyo University of Science, supported by the JSPS Grants-in-Aid for Scientific Research (KAKENHI JP19KK0115, "Three-Dimensional Light Flow") ====================================================================== $Id: pmap.c,v 2.18 2021/02/19 02:10:35 rschregle Exp $ */ #include "pmap.h" #include "pmapmat.h" #include "pmapsrc.h" #include "pmaprand.h" #include "pmaproi.h" #include "pmapio.h" #include "pmapdens.h" #include "pmapbias.h" #include "pmapdiag.h" #include "pmutil.h" #include "otypes.h" #include "otspecial.h" #include #if NIX #include #include #include #include #endif static int photonParticipate (RAY *ray) -/* Trace photon through participating medium. Returns 1 if passed through, - or 0 if absorbed and $*%&ed. Analogon to rayparticipate(). */ +/* Trace photon through participating medium (including dielectric + and interface). Returns 1 if passed through, or 0 if absorbed. + Analogon to rayparticipate(). */ { int i; RREAL xi1, cosTheta, phi, du, dv, rmax0; const RREAL cext = colorAvg(ray -> cext), invCext = 1 / cext, albedo = colorAvg(ray -> albedo), invAlbedo = 1 / albedo, gecc = ray -> gecc, gecc2 = sqr(gecc); FVECT u, v; COLOR cvext; /* Save incoming ray's max distance; this is required since rmax is set to the mean free distance for each path segment (see below) */ rmax0 = ray -> rmax; /* Mean free distance until interaction with medium */ ray -> rmax = -log(pmapRandom(mediumState)) / cext; while (!localhit(ray, &thescene)) { /* No geometry hit within mean free dist ray->rmax; photon interacts with medium */ if (!incube(&thescene, ray -> rop)) { /* Terminate photon if it has leaked from the scene */ #ifdef DEBUG_PMAP_LEAKED fprintf( stderr, "Volume photon leaked from scene at [%.3f %.3f %.3f]\n", ray -> rop [0], ray -> rop [1], ray -> rop [2] ); #endif return 0; } /* RGB attenuation over mean free distance */ setcolor(cvext, exp(-ray -> rmax * ray -> cext [0]), exp(-ray -> rmax * ray -> cext [1]), exp(-ray -> rmax * ray -> cext [2]) ); /* Modify ray color and normalise */ multcolor(ray -> rcol, cvext); colorNorm(ray -> rcol); /* Set ray origin for next interaction */ VCOPY(ray -> rorg, ray -> rop); /* Update ray's path length; terminate photon if maximum reached. See also photonRay(). */ rmax0 -= ray -> rot; if (photonMaxPathDist > 0 && rmax0 < 0) return 0; - /* Store volume photon if nonzero albedo; this also accounts for - direct inscattering from light sources */ - if (albedo > FTINY) { + /* Store volume photon if nonzero albedo (or unconditionally if + lightflow photon); this also accounts for direct inscattering + from light sources */ + if (albedo > FTINY || lightFlowPhotonMapping) { /* Set ray's path length to store with photon */ ray -> rmax = rmax0; -#ifdef PMAP_PHOTONFLOW if (lightFlowPhotonMapping) { /* Temporarily correct normalised flux for lightflow photons based on extinction (= mean num photons deposited per unit path length). */ scalecolor(ray -> rcol, invCext); newPhoton(lightFlowPmap, ray, NULL); /* Undo correction for next iteration after storing photon */ scalecolor(ray -> rcol, cext); } else -#endif newPhoton(volumePmap, ray, NULL); } /* Colour bleeding without attenuation (?) */ multcolor(ray -> rcol, ray -> albedo); scalecolor(ray -> rcol, invAlbedo); - if (pmapRandom(rouletteState) < albedo) { - /* Scatter photon */ + if (pmapRandom(rouletteState) < albedo && !lightFlowPhotonMapping) { + /* Scatter volume photon; lightflow photons are NOT scattered! + + NOTE: "albedo" here specifically refers to the SCATTERING ALBEDO + sigma_s / sigma_t (see Greg's email "Albedo within dielectric / + interface material?" dated 24 Jun 2021), where sigma_s is the + is the scattering probability, and sigma_t is the extinction. + The extinction sigma_t = sigma_a + sigma_s is itself the + sum of the absorption and scattering probabilities, sigma_a + resp. sigma_s. */ xi1 = pmapRandom(scatterState); cosTheta = ray -> gecc <= FTINY ? 2 * xi1 - 1 : 0.5 / gecc * ( 1 + gecc2 - sqr((1 - gecc2) / (1 + gecc * (2 * xi1 - 1))) ); phi = 2 * PI * pmapRandom(scatterState); du = dv = sqrt(1 - sqr(cosTheta)); /* sin(theta) */ du *= cos(phi); dv *= sin(phi); /* Get axes u & v perpendicular to photon direction */ i = 0; do { v [0] = v [1] = v [2] = 0; v [i++] = 1; fcross(u, v, ray -> rdir); } while (normalize(u) < FTINY); fcross(v, ray -> rdir, u); for (i = 0; i < 3; i++) ray -> rdir [i] = du * u [i] + dv * v [i] + cosTheta * ray -> rdir [i]; } +#if 0 + /* Photons should be terminated if not scattered; however, this + unconditionally terminates all photon in dielectrics, since their + scattering albedo is implicitly zero */ + else return 0; +#endif ray -> rlvl++; ray -> rmax = -log(pmapRandom(mediumState)) / cext; /* fprintf(stderr, "%.3f\n", ray -> rmax); */ } /* Passed through medium until intersecting local object */ setcolor(cvext, exp(-ray -> rot * ray -> cext [0]), exp(-ray -> rot * ray -> cext [1]), exp(-ray -> rot * ray -> cext [2]) ); /* Modify ray color and normalise */ multcolor(ray -> rcol, cvext); colorNorm(ray -> rcol); /* Restore outgoing ray's max distance */ ray -> rmax = rmax0; return 1; } void tracePhoton (RAY *ray) /* Follow photon as it bounces around... */ { long mod; 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; /* NOTE: localhit() limits intersections to the aft plane's distance ray->rmax. This mechanism is (ab)used here to limit the photon path length by initialising ray->rmax to photonMaxPathDist; see emitPhoton(). With unlimited path length (photonMaxPathDist=0), ray->rmax becomes negative (see photonRay()), and the aft plane has no effect. */ if (localhit(ray, &thescene)) { mod = ray -> ro -> omod; if (port && ray -> ro != port) { /* Terminate photon if emitted from port without intersecting it */ #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; photonRay(ray, &tray, PMAP_XFER, NULL); tracePhoton(&tray); } else { /* Scatter for modifier material */ mat = objptr(mod); photonScatter [mat -> otype] (mat, ray); } } } static void preComputeGlobal (PhotonMap *pmap) /* 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 p, numPreComp; unsigned j; const double pInc = finalGather > FTINY ? 1 / finalGather : 0; PhotonIdx pIdx; Photon photon; RAY ray; PhotonMap nuPmap; repComplete = numPreComp = finalGather * pmap -> numPhotons; if (verbose) { sprintf(errmsg, "\nPrecomputing irradiance for %ld global photons\n", numPreComp ); eputs(errmsg); #if NIX fflush(stderr); #endif } /* Copy photon map for precomputed photons */ 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 = repLastTime = time(NULL); #ifdef SIGCONT signal(SIGCONT, pmapPreCompReport); #endif repProgress = 0; photonRay(NULL, &ray, PRIMARY, NULL); ray.ro = NULL; for (p = 0; p < numPreComp; p++) { /* Get random photon from stratified distribution in source heap to * avoid duplicates and clustering. * NOTE: multiplication with pInc expanded to prevent rounding error * which would occasionally cause the last photon to be selected * for precomputation again! */ pIdx = firstPhoton(pmap) + (unsigned long)( p * pInc + pmapRandom(pmap -> randState) * pInc ); getPhoton(pmap, pIdx, &photon); /* 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; /* Get density estimate at photon position */ photonDensity(pmap, &ray, ray.rcol); /* Append photon to new heap from ray */ newPhoton(&nuPmap, &ray, NULL); /* Update progress */ repProgress++; if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime) pmapPreCompReport(); #ifdef SIGCONT else signal(SIGCONT, pmapPreCompReport); #endif } /* Flush heap */ flushPhotonHeap(&nuPmap); #ifdef SIGCONT signal(SIGCONT, SIG_DFL); #endif /* 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 } /* Rebuild underlying data structure, destroying heap */ buildPhotonMap(pmap, NULL, NULL, 1); } typedef struct { unsigned long numPhotons [NUM_PMAP_TYPES], numEmitted, numComplete; } PhotonCnt; void distribPhotons (PhotonMap **pmaps, unsigned numProc) { EmissionMap emap; char errmsg2 [128], shmFname [PMAP_TMPFNLEN]; unsigned t, srcIdx, proc; double totalFlux = 0; int shmFile, stat, pid; PhotonMap *pm; PhotonCnt *photonCnt; pid_t procPids [PMAP_MAXPROC]; for (t = 0; t < NUM_PMAP_TYPES && !pmaps [t]; t++); if (t >= NUM_PMAP_TYPES) error(USER, "no photon maps defined in distribPhotons"); if (!nsources) error(USER, "no light sources in distribPhotons"); /* =================================================================== * INITIALISATION - Set up emission and scattering funcs * =================================================================== */ emap.samples = NULL; emap.maxPartitions = MAXSPART; emap.partitions = (unsigned char*)malloc(emap.maxPartitions >> 1); if (!emap.partitions) error(INTERNAL, "can't allocate source partitions in distribPhotons"); /* Initialise all defined photon maps */ for (t = 0; t < NUM_PMAP_TYPES; 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 from modifier list */ getPhotonPorts(photonPortList); /* Get photon sensor modifiers */ getPhotonSensors(photonSensorList); /* Get polyhedral regions of interest */ getPolyROIs(pmapROImodList); #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); } /* =================================================================== * FLUX INTEGRATION - Get total photon flux from light sources * =================================================================== */ for (srcIdx = 0; srcIdx < nsources; srcIdx++) { unsigned portCnt = 0; emap.src = source + srcIdx; do { /* Set next photon port if defined; note this loop iterates once if * no ports are defined. */ emap.port = emap.src -> sflags & SDISTANT ? photonPorts + portCnt : NULL; photonPartition [emap.src -> so -> otype] (&emap); 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); strcat(errmsg, errmsg2); eputs(errmsg); #if NIX fflush(stderr); #endif } for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions; emap.partitionCnt++ ) { initPhotonEmission(&emap, pdfSamples); totalFlux += colorAvg(emap.partFlux); } portCnt++; } while (portCnt < numPhotonPorts); } if (totalFlux < FTINY) error(USER, "zero flux from light sources"); /* Record start time for progress reports */ repStartTime = time(NULL); if (verbose) { sprintf(errmsg, "\nPhoton distribution @ %d procs\n", numProc); eputs(errmsg); } /* MAIN LOOP */ for (proc = 0; proc < numProc; proc++) { #if NIX if (!(pid = fork())) { /* SUBPROCESS ENTERS HERE; 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 preIncPhotonCnt, lastNumPhotons [NUM_PMAP_TYPES], numEmitted = 0, lastNumEmitted = 0; /* Decorrelate shared photon counter updates to reduce * contention by setting unique update intervals per subproc */ const unsigned photonCntUpdate = PMAP_CNTUPDATE + proc; /* 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 DEBUG_PMAP /* Output child process PID after random delay to prevent corrupted * console output due to race condition */ usleep(1e6 * pmapRandom(rouletteState)); fprintf(stderr, "Proc %d: PID = %d (waiting 10 sec to attach debugger...)\n", proc, getpid() ); /* Allow time for debugger to attach to child process */ sleep(10); #endif for (t = 0; t < NUM_PMAP_TYPES; t++) lastNumPhotons [t] = 0; /* ============================================================= * 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 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); } sprintf(errmsg2, "; increase -apD"); strcat(errmsg, errmsg2); error(USER, errmsg); break; } /* Num to emit is fraction of minimum target count */ numEmit = FHUGE; 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. */ 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 ); } /* 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 = numEmitted * maxDistribRatio) < FTINY) /* No photons left to distribute in main pass */ break; } /* PHOTON DISTRIBUTION LOOP */ for (srcIdx = 0; srcIdx < nsources; srcIdx++) { unsigned portCnt = 0; emap.src = source + srcIdx; do { /* Set next photon port if defined; note this loop iterates * once if no ports are defined. */ 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++ ) { 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++; /* Integer counter avoids FP rounding errors during * iteration */ while (partEmitCnt--) { RAY photonRay; /* Emit photon based on PDF, skip if rejected. */ /* NOTE: rejection sampling skips incrementing the * emission counter (see below), thus compensating * for the rejected photons by increasing the photon * flux in proportion to the lower effective * emission count. * BUG: THIS INTERFERES WITH THE PROGRESS COUNTER * REPORTED TO THE PARENT, AND WITH THE * PREDISTRIBUTION PASS --> PHOTON DISTRIBUTION WILL * FINISH EARLY, WITH FEWER PHOTONS THAN TARGETTED! */ if (!emitPhoton(&emap, &photonRay)) continue; /* Set subprocess ID in photonRay; extends path ID. Used by newPhoton() as photon attributes. */ PMAP_SETRAYPROC(&photonRay, proc); tracePhoton(&photonRay); /* Update local emission counter*/ numEmitted++; if (!(partEmitCnt && (numEmitted & photonCntUpdate))) { /* Update global counters shared with siblings in intervals to reduce overhead & contention; ALSO WHEN EMITTING FINAL PHOTON FROM THIS PARTITION, OTHERWISE COUNTERS WILL SKIP LAST UPDATE, BIASING THE PHOTON FLUX! */ #if NIX /* Shared photon counter file must be locked! */ shmLock(shmFile, F_WRLCK); #endif /* Differentially increment number emitted using * local counter numEmitted */ preIncPhotonCnt = photonCnt -> numEmitted; photonCnt -> numEmitted += numEmitted - lastNumEmitted; photonCnt -> numComplete += numEmitted - lastNumEmitted; lastNumEmitted = numEmitted; /* Check overflow using pre-increment values */ if (photonCnt -> numEmitted < preIncPhotonCnt) { sprintf(errmsg, "photon emission counter " "overflow (was: %ld, is: %ld)", preIncPhotonCnt, photonCnt -> numEmitted ); error(INTERNAL, errmsg); } /* Update global photon count for each pmap */ for (t = 0; t < NUM_PMAP_TYPES; t++) { if (pmaps [t]) { preIncPhotonCnt = photonCnt->numPhotons [t]; /* Differentially increment using local counter lastNumPhotons */ photonCnt -> numPhotons [t] += pmaps [t] -> numPhotons - lastNumPhotons [t]; lastNumPhotons [t] = pmaps [t] -> numPhotons; /* Check for photon counter overflow (this * could only happen before an emission * counter overflow if the scene has an * absurdly high albedo and/or very dense * geometry) */ if (photonCnt -> numPhotons [t] < preIncPhotonCnt) { sprintf(errmsg, "photon counter " "overflow (was: %ld, is: %ld)", preIncPhotonCnt, photonCnt -> numPhotons [t] ); error(INTERNAL, errmsg); } } } #if NIX /* Release lock on shared photon counter file! */ shmLock(shmFile, F_UNLCK); #endif } } #if !NIX /* Synchronous progress report on Windoze */ if (!proc && photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime) { repEmitted = repProgress = photonCnt -> numEmitted; repComplete = photonCnt -> numComplete; for (t = 0; t < NUM_PMAP_TYPES; t++) if (pmaps [t]) { repProgress += pmaps [t] -> numPhotons; repComplete += pmaps [t] -> distribTarget; } 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 } #if NIX /* Terminate subprocess */ exit(0); #endif } else { /* PARENT PROCESS CONTINUES LOOP ITERATION HERE */ if (pid < 0) error(SYSTEM, "failed to fork subprocess in distribPhotons"); else /* Save child process IDs */ procPids [proc] = pid; } } #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)) { /* Exited with error; terminate siblings, clean up & bail out */ for (proc = 0; proc < numProc; proc++) kill(procPids [proc], SIGKILL); /* Unmap shared memory, squish mapped file */ munmap(photonCnt, sizeof(*photonCnt)); close(shmFile); unlink(shmFname); error(USER, "failed photon distribution"); } --proc; } /* Nod off for a bit and update progress */ sleep(1); /* Asynchronous progress report from shared subprocess counters */ shmLock(shmFile, F_RDLCK); repEmitted = repProgress = photonCnt -> numEmitted; repComplete = photonCnt -> numComplete; for (t = repProgress = repComplete = 0; t < NUM_PMAP_TYPES; t++) if ((pm = pmaps [t])) { /* Get global photon count from shmem updated by subprocs */ repProgress += pm -> numPhotons = photonCnt -> numPhotons [t]; repComplete += pm -> distribTarget; } repComplete *= numProc; shmLock(shmFile, F_UNLCK); if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime) pmapDistribReport(); #ifdef SIGCONT else signal(SIGCONT, pmapDistribReport); #endif } #endif /* NIX */ /* =================================================================== * POST-DISTRIBUTION - Set photon flux and build data struct for photon * storage, etc. * =================================================================== */ #ifdef SIGCONT /* Reset signal handler */ signal(SIGCONT, SIG_DFL); #endif free(emap.samples); /* Set photon flux */ 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 (pmaps [t]) { if (verbose) { sprintf(errmsg, "Building %s photon map\n", pmapName [t]); eputs(errmsg); #if NIX fflush(stderr); #endif } /* Build underlying data structure; heap is destroyed */ buildPhotonMap(pmaps [t], &totalFlux, NULL, numProc); } /* Precompute photon irradiance if necessary */ if (preCompPmap) { if (verbose) eputs("\n"); preComputeGlobal(preCompPmap); } if (verbose) eputs("\n"); }