diff --git a/pmap.c b/pmap.c index c2c24f1..7936d15 100644 --- a/pmap.c +++ b/pmap.c @@ -1,935 +1,944 @@ #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 "pmapio.h" #include "pmapdens.h" #include "pmapbias.h" #include "pmapdiag.h" #include "otypes.h" #include "otspecial.h" #include #if NIX #include #include #include #endif /* Photon distribution counter update interval expressed as bitmask; counters shared among subling subprocesses will only be updated in multiples of PMAP_CNTUPDATE in order to reduce contention */ #define PMAP_CNTUPDATE 0xffL 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, argc, argv); } } static int photonParticipate (RAY *ray) /* Trace photon through participating medium. Returns 1 if passed through, or 0 if absorbed and $*%&ed. Analogon to rayparticipate(). */ { int i; RREAL xi1, cosTheta, phi, du, dv; - const float cext = colorAvg(ray -> cext), - albedo = colorAvg(ray -> albedo), + const float 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; +#ifdef DEBUG_PHOTONSPERPATH unsigned long photonsPerPath = 0, lastNumPhotons; +#endif /* Mean free distance until interaction with medium */ ray -> rmax = -log(pmapRandom(mediumState)) / cext; while (!localhit(ray, &thescene)) { if (!incube(&thescene, ray -> rop)) { /* Terminate photon if it has leaked from the scene */ #ifdef DEBUG_PMAP_LEAKED fprintf( stderr, "Volume photon leaked from scene at [%.3f %.3f %.3f]\n", ray -> rop [0], ray -> rop [1], ray -> rop [2] ); #endif #ifdef DEBUG_PHOTONSPERPATH printf("Photons/path =\t%ld\n", photonsPerPath); #endif return 0; } - setcolor(cvext, + /* 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); - + +#ifdef PMAP_PHOTONFLOW + if (lightFlowPhotonMapping) + /* Correct normalised flux for lightflow photons based on + extinction (= mean num photons deposited per unit path length) */ + scalecolor(ray -> rcol, invCext); +#endif #if 0 if (albedo > FTINY && ray -> rlvl > 0) { #else - /* Store volume photons unconditionally in mist to also account for - direct inscattering from sources */ + /* Store volume photons unconditionally in mist to also account for + direct inscattering from sources */ if (albedo > FTINY) { #endif /* Add to volume photon map (implicitly overidden by lightFlowPmap in photon flow mode via volumePmap macro) */ +#ifdef DEBUG_PHOTONSPERPATH lastNumPhotons = volumePmap -> numPhotons; +#endif newPhoton(volumePmap, ray); +#ifdef DEBUG_PHOTONSPERPATH photonsPerPath += volumePmap -> numPhotons - lastNumPhotons; +#endif } /* Absorbed? */ if (pmapRandom(rouletteState) > albedo) return 0; - + /* Colour bleeding without attenuation (?) */ multcolor(ray -> rcol, ray -> albedo); - scalecolor(ray -> rcol, 1 / albedo); + scalecolor(ray -> rcol, invAlbedo); /* Scatter photon */ 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]; + ray -> rdir [i] = + du * u [i] + dv * v [i] + cosTheta * ray -> rdir [i]; ray -> rlvl++; ray -> rmax = -log(pmapRandom(mediumState)) / cext; /* fprintf(stderr, "%.3f\n", ray -> rmax); */ } - /* Passed through medium until intersecting local object */ + /* 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); 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; if (localhit(ray, &thescene)) { mod = ray -> ro -> omod; /* XXX: Is port -> omod != mod sufficient here? Probably not... */ if ( port && ray -> ro != port && findmaterial(port) != findmaterial(ray -> ro) ) { /* !!! PHOTON PORT REJECTION SAMPLING HACK !!! * Terminate photon if emitted from port without intersecting it or * its other associated surfaces or same material. * This can happen when the port's partitions extend beyond its * actual geometry, e.g. with polygons. Since the total flux * relayed by the port is based on the (in this case) larger * partition area, it is overestimated; terminating these photons * constitutes rejection sampling and thereby compensates any bias * incurred by the overestimated flux. */ #ifdef PMAP_PORTREJECT_WARN sprintf(errmsg, "photon outside port %s", ray -> ro -> oname); error(WARNING, errmsg); #endif return; } if ((ray -> clipset && inset(ray -> clipset, mod)) || mod == OVOID) { /* Transfer ray if modifier is VOID or clipped within antimatta */ RAY tray; 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 i, numPreComp; unsigned j; 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 = time(NULL); #ifdef SIGCONT signal(SIGCONT, pmapPreCompReport); #endif repProgress = 0; 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); /* 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); /* 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, lastPhotonCnt; 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); #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"); /* Zero overflow tracking counters */ memset(&lastPhotonCnt, 0, sizeof(lastPhotonCnt)); #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 lastNumPhotons [NUM_PMAP_TYPES]; unsigned long localNumEmitted = 0; /* Num photons emitted by this subprocess alone */ /* 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); } 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 = localNumEmitted * 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++; #if 0 /* Update local and shared (global) emission counter */ photonCnt -> numEmitted += partEmitCnt; localNumEmitted += partEmitCnt; /* Update shared completion counter for progress report by parent */ photonCnt -> numComplete += partEmitCnt; /* Check for overflowed emission/completion counters */ if ( photonCnt -> numEmitted < lastPhotonCnt.numEmitted || photonCnt -> numComplete < lastPhotonCnt.numComplete ) error(INTERNAL, "photon emission counter overflow"); else { /* Checked ok, update overflow tracking counters */ lastPhotonCnt.numEmitted = photonCnt -> numEmitted; lastPhotonCnt.numComplete = photonCnt -> numComplete; } #endif /* 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() */ /* TODO: Move this into emitPhoton()? */ photonRay.ro = emap.port -> so; #endif /* Set subprocess ID in photonRay; extends path ID. Used by newPhoton() as photon attributes. */ PMAP_SETRAYPROC(&photonRay, proc); tracePhoton(&photonRay); #if 1 /* Update local emission counter*/ localNumEmitted++; if (!(localNumEmitted & PMAP_CNTUPDATE)) { /* Update shared (global) counters with lower frequency to reduce overhead and contention with siblings */ photonCnt -> numEmitted += PMAP_CNTUPDATE; photonCnt -> numComplete += PMAP_CNTUPDATE; /* Check for overflow using previous values */ /* XXX: NEED LOCKING HERE??? */ if ( photonCnt -> numEmitted < lastPhotonCnt.numEmitted || photonCnt -> numComplete < lastPhotonCnt.numComplete ) error(INTERNAL, "photon emission counter overflow"); else { /* Checked ok, update overflow tracking counters */ lastPhotonCnt.numEmitted = photonCnt -> numEmitted; lastPhotonCnt.numComplete = photonCnt -> numComplete; } /* 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; /* 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) */ /* XXX: NEED LOCKING HERE??? */ if ( photonCnt -> numPhotons [t] < lastPhotonCnt.numPhotons [t] ) error( INTERNAL, "photon counter overflow" ); else lastPhotonCnt.numPhotons [t] = photonCnt -> numPhotons [t]; } } #endif } #if 0 /* 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; /* 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] < lastPhotonCnt.numPhotons [t] ) error( INTERNAL, "photon counter overflow" ); else lastPhotonCnt.numPhotons [t] = photonCnt -> numPhotons [t]; } #endif #if !NIX /* Synchronous progress report on Windoze */ if (!proc && photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime) { repEmitted = repProgress = photonCnt -> numEmitted; repComplete = photonCnt -> numComplete; pmapDistribReport(); } #endif } portCnt++; } 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 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 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 */ repEmitted = repProgress = photonCnt -> numEmitted; repComplete = photonCnt -> numComplete; repProgress = repComplete = 0; for (t = 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; 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; - - #ifdef PMAP_PHOTONFLOW - /* Adjust flux for mean volume photons deposited per unit path length */ - /* totalFlux /= */ - #endif -#if NIX +#if NIX /* Photon counters no longer needed, unmap shared memory */ munmap(photonCnt, sizeof(*photonCnt)); close(shmFile); unlink(shmFname); #else free(photonCnt); -#endif +#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"); } diff --git a/pmapdata.c b/pmapdata.c index e7c9a0f..1671726 100644 --- a/pmapdata.c +++ b/pmapdata.c @@ -1,751 +1,753 @@ #ifndef lint static const char RCSid[] = "$Id: pmapdata.c,v 2.23 2021/04/14 11:26:25 rschregle Exp $"; #endif /* ========================================================================= Photon map types and high level interface to nearest neighbour lookups in underlying point cloud data structure. The default data structure is an in-core kd-tree (see pmapkdt.{h,c}). This can be overriden with the PMAP_OOC compiletime switch, which enables an out-of-core octree (see oococt.{h,c}). 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: pmapdata.c,v 2.23 2021/04/14 11:26:25 rschregle Exp $ */ #include "pmapdata.h" #include "pmapdens.h" #include "pmaprand.h" #include "pmapmat.h" #include "otypes.h" #include "random.h" PhotonMap *photonMaps [NUM_PMAP_TYPES] = { NULL, NULL, NULL, NULL, NULL, NULL, #ifdef PMAP_PHOTONFLOW NULL #endif }; /* Include routines to handle underlying point cloud data structure */ #ifdef PMAP_OOC #include "pmapooc.c" #else #include "pmapkdt.c" #endif /* Ambient include/exclude set (from ambient.c) */ #ifndef MAXASET #define MAXASET 4095 #endif extern OBJECT ambset [MAXASET+1]; /* Callback to print photon attributes acc. to user defined format */ int (*printPhoton)(RAY *r, Photon *p, PhotonMap *pm); /* PHOTON MAP BUILD ROUTINES ------------------------------------------ */ void initPhotonMap (PhotonMap *pmap, PhotonMapType t) /* Init photon map 'n' stuff... */ { if (!pmap) return; pmap -> numPhotons = 0; pmap -> biasCompHist = NULL; pmap -> maxPos [0] = pmap -> maxPos [1] = pmap -> maxPos [2] = -FHUGE; pmap -> minPos [0] = pmap -> minPos [1] = pmap -> minPos [2] = FHUGE; pmap -> minGathered = pmap -> maxGathered = pmap -> totalGathered = 0; pmap -> gatherTolerance = gatherTolerance; pmap -> minError = pmap -> maxError = pmap -> rmsError = 0; pmap -> numDensity = 0; pmap -> distribRatio = 1; pmap -> type = t; pmap -> squeue.node = NULL; pmap -> squeue.len = 0; #ifdef PMAP_PHOTONFLOW /* Init path lookup table and its key buffer */ pmap -> pathLUT = NULL; pmap -> pathLUTKeys = NULL; pmap -> numPathLUTKeys = 0; #endif /* Init local RNG state */ pmap -> randState [0] = 10243; pmap -> randState [1] = 39829; pmap -> randState [2] = 9433; pmapSeed(randSeed, pmap -> randState); /* Set up type-specific photon lookup callback */ pmap -> lookup = pmapLookup [t]; /* Mark primary photon ray as unused */ pmap -> lastPrimary.srcIdx = -1; pmap -> numPrimary = 0; pmap -> primaries = NULL; /* Init storage */ pmap -> heap = NULL; pmap -> heapBuf = NULL; pmap -> heapBufLen = 0; #ifdef PMAP_OOC OOC_Null(&pmap -> store); #else kdT_Null(&pmap -> store); #endif } void initPhotonHeap (PhotonMap *pmap) { int fdFlags; if (!pmap) error(INTERNAL, "undefined photon map in initPhotonHeap"); if (!pmap -> heap) { /* Open heap file */ mktemp(strcpy(pmap -> heapFname, PMAP_TMPFNAME)); if (!(pmap -> heap = fopen(pmap -> heapFname, "w+b"))) error(SYSTEM, "failed opening heap file in initPhotonHeap"); #ifdef F_SETFL /* XXX is there an alternative needed for Windows? */ fdFlags = fcntl(fileno(pmap -> heap), F_GETFL); fcntl(fileno(pmap -> heap), F_SETFL, fdFlags | O_APPEND); #endif/* ftruncate(fileno(pmap -> heap), 0); */ } } void flushPhotonHeap (PhotonMap *pmap) { int fd; const unsigned long len = pmap -> heapBufLen * sizeof(Photon); if (!pmap) error(INTERNAL, "undefined photon map in flushPhotonHeap"); if (!pmap -> heap || !pmap -> heapBuf) { /* Silently ignore undefined heap error(INTERNAL, "undefined heap in flushPhotonHeap"); */ return; } /* Atomically seek and write block to heap */ /* !!! Unbuffered I/O via pwrite() avoids potential race conditions * !!! and buffer corruption which can occur with lseek()/fseek() * !!! followed by write()/fwrite(). */ fd = fileno(pmap -> heap); #ifdef DEBUG_PMAP sprintf(errmsg, "Proc %d: flushing %ld photons from pos %ld\n", getpid(), pmap -> heapBufLen, lseek(fd, 0, SEEK_END) / sizeof(Photon)); eputs(errmsg); #endif /*if (pwrite(fd, pmap -> heapBuf, len, lseek(fd, 0, SEEK_END)) != len) */ if (write(fd, pmap -> heapBuf, len) != len) error(SYSTEM, "failed append to heap file in flushPhotonHeap"); #if NIX if (fsync(fd)) error(SYSTEM, "failed fsync in flushPhotonHeap"); #endif pmap -> heapBufLen = 0; } #ifdef DEBUG_PMAP static int checkPhotonHeap (FILE *file) /* Check heap for nonsensical or duplicate photons */ { Photon p, lastp; int i, dup; rewind(file); memset(&lastp, 0, sizeof(lastp)); 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]; } 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 int newPhoton (PhotonMap* pmap, const RAY* ray) { unsigned i; Photon photon; COLOR photonFlux; /* Account for distribution ratio */ if (!pmap || pmapRandom(pmap -> randState) > pmap -> distribRatio) return -1; /* Don't store on sources */ if (ray -> robj > -1 && islight(objptr(ray -> ro -> omod) -> otype)) return -1; /* Ignore photon if modifier in/outside exclude/include set */ if (ambincl != -1 && ray -> ro && ambincl != inset(ambset, ray -> ro -> omod)) return -1; if (pmapNumROI && pmapROI) { unsigned inROI = 0; FVECT photonDist; /* Store photon if within a region of interest (for ze Ecksperts!) Note size of spherical ROI is squared. */ for (i = 0; !inROI && i < pmapNumROI; i++) { VSUB(photonDist, ray -> rop, pmapROI [i].pos); inROI = ( PMAP_ROI_ISSPHERE(pmapROI + i) ? DOT(photonDist, photonDist) <= pmapROI [i].siz [0] : fabs(photonDist [0]) <= pmapROI [i].siz [0] && fabs(photonDist [1]) <= pmapROI [i].siz [1] && fabs(photonDist [2]) <= pmapROI [i].siz [2] ); } if (!inROI) return -1; } /* Adjust flux according to distribution ratio and ray weight */ copycolor(photonFlux, ray -> rcol); #if 0 /* Factored out ray -> rweight as deprecated (?) for pmap, and infact erroneously attenuates volume photon flux based on extinction, which is already factored in by photonParticipate() */ scalecolor(photonFlux, ray -> rweight / (pmap -> distribRatio ? pmap -> distribRatio : 1)); #else - scalecolor(photonFlux, - 1.0 / (pmap -> distribRatio ? pmap -> distribRatio : 1)); + scalecolor( + photonFlux, + 1.0 / (pmap -> distribRatio ? pmap -> distribRatio : 1) + ); #endif setPhotonFlux(&photon, photonFlux); /* Set photon position and flags */ VCOPY(photon.pos, ray -> rop); photon.flags = 0; photon.caustic = PMAP_CAUSTICRAY(ray); /* Set contrib photon's primary ray and subprocess index (the latter * to linearise the primary ray indices after photon distribution is * complete). Also set primary ray's source index, thereby marking it * as used. */ if (isContribPmap(pmap)) { photon.primary = pmap -> numPrimary; photon.proc = PMAP_GETRAYPROC(ray); pmap -> lastPrimary.srcIdx = ray -> rsrc; } /* Set photon's path ID from ray and subprocess */ - else { + else { photon.primary = ray -> rno; photon.proc = PMAP_GETRAYPROC(ray); } - + /* Set normal */ for (i = 0; i <= 2; i++) photon.norm [i] = 127.0 * (isVolumePmap(pmap) ? ray -> rdir [i] : ray -> ron [i] ); if (!pmap -> heapBuf) { /* Lazily allocate heap buffa */ #if NIX /* Randomise buffa size to temporally decorellate flushes in * multiprocessing mode */ srandom(randSeed + getpid()); pmap -> heapBufSize = PMAP_HEAPBUFSIZE * (0.5 + frandom()); #else /* Randomisation disabled for single processes on WIN; also useful - * for reproducability during debugging */ + * for reproducability during debugging */ pmap -> heapBufSize = PMAP_HEAPBUFSIZE; -#endif +#endif if (!(pmap -> heapBuf = calloc(pmap -> heapBufSize, sizeof(Photon)))) error(SYSTEM, "failed heap buffer allocation in newPhoton"); - pmap -> heapBufLen = 0; + pmap -> heapBufLen = 0; } /* Photon initialised; now append to heap buffa */ memcpy(pmap -> heapBuf + pmap -> heapBufLen, &photon, sizeof(Photon)); - + if (++pmap -> heapBufLen >= pmap -> heapBufSize) /* Heap buffa full, flush to heap file */ flushPhotonHeap(pmap); pmap -> numPhotons++; - + /* Print photon attributes */ if (printPhoton) /* Non-const kludge */ printPhoton((RAY*)ray, &photon, pmap); - + return 0; } void buildPhotonMap (PhotonMap *pmap, double *photonFlux, PhotonPrimaryIdx *primaryOfs, unsigned nproc) { unsigned long n, nCheck = 0; unsigned i; Photon *p; COLOR flux; char nuHeapFname [sizeof(PMAP_TMPFNAME)]; FILE *nuHeap; /* Need double here to reduce summation errors */ double avgFlux [3] = {0, 0, 0}, CoG [3] = {0, 0, 0}, CoGdist = 0; FVECT d; if (!pmap) error(INTERNAL, "undefined photon map in buildPhotonMap"); /* Get number of photons from heapfile size */ if (fseek(pmap -> heap, 0, SEEK_END) < 0) error(SYSTEM, "failed seek to end of photon heap in buildPhotonMap"); pmap -> numPhotons = ftell(pmap -> heap) / sizeof(Photon); if (!pmap -> numPhotons) error(INTERNAL, "empty photon map in buildPhotonMap"); if (!pmap -> heap) error(INTERNAL, "no heap in buildPhotonMap"); #ifdef DEBUG_PMAP eputs("Checking photon heap consistency...\n"); checkPhotonHeap(pmap -> heap); sprintf(errmsg, "Heap contains %ld photons\n", pmap -> numPhotons); eputs(errmsg); #endif /* Allocate heap buffa */ if (!pmap -> heapBuf) { pmap -> heapBufSize = PMAP_HEAPBUFSIZE; pmap -> heapBuf = calloc(pmap -> heapBufSize, sizeof(Photon)); if (!pmap -> heapBuf) error(SYSTEM, "failed to allocate postprocessed photon heap in" "buildPhotonMap"); } /* We REALLY don't need yet another @%&*! heap just to hold the scaled * photons, but can't think of a quicker fix... */ mktemp(strcpy(nuHeapFname, PMAP_TMPFNAME)); if (!(nuHeap = fopen(nuHeapFname, "w+b"))) error(SYSTEM, "failed to open postprocessed photon heap in " "buildPhotonMap"); rewind(pmap -> heap); #ifdef DEBUG_PMAP eputs("Postprocessing photons...\n"); #endif while (!feof(pmap -> heap)) { #ifdef DEBUG_PMAP printf("Reading %lu at %lu... ", pmap -> heapBufSize, ftell(pmap->heap)); #endif pmap -> heapBufLen = fread(pmap -> heapBuf, sizeof(Photon), pmap -> heapBufSize, pmap -> heap); #ifdef DEBUG_PMAP printf("Got %lu\n", pmap -> heapBufLen); #endif if (ferror(pmap -> heap)) error(SYSTEM, "failed to read photon heap in buildPhotonMap"); for (n = pmap -> heapBufLen, p = pmap -> heapBuf; n; n--, p++) { /* Update min and max pos and set photon flux */ for (i = 0; i <= 2; i++) { if (p -> pos [i] < pmap -> minPos [i]) pmap -> minPos [i] = p -> pos [i]; else if (p -> pos [i] > pmap -> maxPos [i]) pmap -> maxPos [i] = p -> pos [i]; /* Update centre of gravity with photon position */ CoG [i] += p -> pos [i]; } if (primaryOfs) /* Linearise photon primary index from subprocess index using the * per-subprocess offsets in primaryOfs */ p -> primary += primaryOfs [p -> proc]; /* Scale photon's flux (hitherto normalised to 1 over RGB); in * case of a contrib photon map, this is done per light source, * and photonFlux is assumed to be an array */ getPhotonFlux(p, flux); if (photonFlux) { scalecolor(flux, photonFlux [isContribPmap(pmap) ? photonSrcIdx(pmap, p) : 0]); setPhotonFlux(p, flux); } /* Update average photon flux; need a double here */ addcolor(avgFlux, flux); } /* Write modified photons to new heap */ fwrite(pmap -> heapBuf, sizeof(Photon), pmap -> heapBufLen, nuHeap); if (ferror(nuHeap)) error(SYSTEM, "failed postprocessing photon flux in " "buildPhotonMap"); nCheck += pmap -> heapBufLen; } #ifdef DEBUG_PMAP if (nCheck < pmap -> numPhotons) error(INTERNAL, "truncated photon heap in buildPhotonMap"); #endif /* Finalise average photon flux */ scalecolor(avgFlux, 1.0 / pmap -> numPhotons); copycolor(pmap -> photonFlux, avgFlux); /* Average photon positions to get centre of gravity */ for (i = 0; i < 3; i++) pmap -> CoG [i] = CoG [i] /= pmap -> numPhotons; rewind(pmap -> heap); /* Compute average photon distance to centre of gravity */ while (!feof(pmap -> heap)) { pmap -> heapBufLen = fread(pmap -> heapBuf, sizeof(Photon), pmap -> heapBufSize, pmap -> heap); for (n = pmap -> heapBufLen, p = pmap -> heapBuf; n; n--, p++) { VSUB(d, p -> pos, CoG); CoGdist += DOT(d, d); } } pmap -> CoGdist = CoGdist /= pmap -> numPhotons; /* Swap heaps, discarding unscaled photons */ fclose(pmap -> heap); unlink(pmap -> heapFname); pmap -> heap = nuHeap; strcpy(pmap -> heapFname, nuHeapFname); #ifdef PMAP_OOC OOC_BuildPhotonMap(pmap, nproc); #else kdT_BuildPhotonMap(pmap); #endif /* Trash heap and its buffa */ free(pmap -> heapBuf); fclose(pmap -> heap); unlink(pmap -> heapFname); pmap -> heap = NULL; pmap -> heapBuf = NULL; } /* PHOTON MAP SEARCH ROUTINES ------------------------------------------ */ /* Dynamic max photon search radius increase and reduction factors */ #define PMAP_MAXDIST_INC 4 #define PMAP_MAXDIST_DEC 0.9 /* Num successful lookups before reducing in max search radius */ #define PMAP_MAXDIST_CNT 1000 /* Threshold below which we assume increasing max radius won't help */ #define PMAP_SHORT_LOOKUP_THRESH 1 /* Coefficient for adaptive maximum search radius */ #define PMAP_MAXDIST_COEFF 100 void findPhotons (PhotonMap* pmap, const RAY* ray) { int redo = 0; if (!pmap -> squeue.len) { /* Lazy init priority queue */ #ifdef PMAP_OOC OOC_InitFindPhotons(pmap); #else kdT_InitFindPhotons(pmap); #endif pmap -> minGathered = pmap -> maxGather; pmap -> maxGathered = pmap -> minGather; pmap -> totalGathered = 0; pmap -> numLookups = pmap -> numShortLookups = 0; pmap -> shortLookupPct = 0; pmap -> minError = FHUGE; pmap -> maxError = -FHUGE; pmap -> rmsError = 0; /* SQUARED max search radius limit is based on avg photon distance to * centre of gravity, unless fixed by user (maxDistFix > 0) */ pmap -> maxDist0 = pmap -> maxDist2Limit = maxDistFix > 0 ? maxDistFix * maxDistFix : PMAP_MAXDIST_COEFF * pmap -> squeue.len * pmap -> CoGdist / pmap -> numPhotons; } do { pmap -> squeue.tail = 0; pmap -> maxDist2 = pmap -> maxDist0; /* XXX/NOTE: status returned by XXX_FindPhotons() is currently ignored; if no photons are found, an empty queue is returned under the assumption all photons are too distant to contribute significant flux. */ #ifdef PMAP_PHOTONFLOW if (isVolumePmap(pmap) && !isLightFlowPmap(pmap)) { #else if (isVolumePmap(pmap)) { #endif /* Search position is ray -> rorg for volume photons, since we have no intersection point. Photon normals are ignored; these are incident directions for the volume scattering func). */ #ifdef PMAP_OOC OOC_FindPhotons(pmap, ray -> rorg, NULL); #else kdT_FindPhotons(pmap, ray -> rorg, NULL); #endif } else { #ifdef PMAP_OOC OOC_FindPhotons(pmap, ray -> rop, ray -> ron); #else kdT_FindPhotons(pmap, ray -> rop, ray -> ron); #endif } #ifdef PMAP_LOOKUP_INFO fprintf(stderr, "%d/%d %s photons found within radius %.3f " "at (%.2f,%.2f,%.2f) on %s\n", pmap -> squeue.tail, pmap -> squeue.len, pmapName [pmap -> type], sqrt(pmap -> maxDist2), ray -> rop [0], ray -> rop [1], ray -> rop [2], ray -> ro ? ray -> ro -> oname : ""); #endif if (pmap -> squeue.tail < pmap -> squeue.len * pmap -> gatherTolerance) { /* Short lookup; too few photons found */ if (pmap -> squeue.tail > PMAP_SHORT_LOOKUP_THRESH) { /* Ignore short lookups which return fewer than * PMAP_SHORT_LOOKUP_THRESH photons under the assumption there * really are no photons in the vicinity, and increasing the max * search radius therefore won't help */ #ifdef PMAP_LOOKUP_WARN sprintf(errmsg, "%d/%d %s photons found at (%.2f,%.2f,%.2f) on %s", pmap -> squeue.tail, pmap -> squeue.len, pmapName [pmap -> type], ray -> rop [0], ray -> rop [1], ray -> rop [2], ray -> ro ? ray -> ro -> oname : ""); error(WARNING, errmsg); #endif /* Bail out after warning if maxDist is fixed */ if (maxDistFix > 0) return; if (pmap -> maxDist0 < pmap -> maxDist2Limit) { /* Increase max search radius if below limit & redo search */ pmap -> maxDist0 *= PMAP_MAXDIST_INC; #ifdef PMAP_LOOKUP_REDO redo = 1; #endif #ifdef PMAP_LOOKUP_WARN sprintf(errmsg, redo ? "restarting photon lookup with max radius %.1e" : "max photon lookup radius adjusted to %.1e", pmap -> maxDist0); error(WARNING, errmsg); #endif } #ifdef PMAP_LOOKUP_REDO else { sprintf(errmsg, "max photon lookup radius clamped to %.1e", pmap -> maxDist0); error(WARNING, errmsg); } #endif } /* Reset successful lookup counter */ pmap -> numLookups = 0; } else { /* Skip search radius reduction if maxDist is fixed */ if (maxDistFix > 0) return; /* Increment successful lookup counter and reduce max search radius if * wraparound */ pmap -> numLookups = (pmap -> numLookups + 1) % PMAP_MAXDIST_CNT; if (!pmap -> numLookups) pmap -> maxDist0 *= PMAP_MAXDIST_DEC; redo = 0; } } while (redo); } Photon *find1Photon (PhotonMap *pmap, const RAY* ray, Photon *photon) { /* Init (squared) search radius to avg photon dist to centre of gravity */ float maxDist2_0 = pmap -> CoGdist; int found = 0; #ifdef PMAP_LOOKUP_REDO #define REDO 1 #else #define REDO 0 #endif do { pmap -> maxDist2 = maxDist2_0; #ifdef PMAP_OOC found = OOC_Find1Photon(pmap, ray -> rop, ray -> ron, photon); #else found = kdT_Find1Photon(pmap, ray -> rop, ray -> ron, photon); #endif if (found < 0) { /* Expand search radius to retry */ maxDist2_0 *= 2; #ifdef PMAP_LOOKUP_WARN sprintf(errmsg, "failed 1-NN photon lookup" #ifdef PMAP_LOOKUP_REDO ", retrying with search radius %.2f", maxDist2_0 #endif ); error(WARNING, errmsg); #endif } } while (REDO && found < 0); /* Return photon buffer containing valid photon, else NULL */ return found < 0 ? NULL : photon; } void getPhoton (PhotonMap *pmap, PhotonIdx idx, Photon *photon) { #ifdef PMAP_OOC if (OOC_GetPhoton(pmap, idx, photon)) #else if (kdT_GetPhoton(pmap, idx, photon)) #endif error(INTERNAL, "failed photon lookup"); } Photon *getNearestPhoton (const PhotonSearchQueue *squeue, PhotonIdx idx) { #ifdef PMAP_OOC return OOC_GetNearestPhoton(squeue, idx); #else return kdT_GetNearestPhoton(squeue, idx); #endif } PhotonIdx firstPhoton (const PhotonMap *pmap) { #ifdef PMAP_OOC return OOC_FirstPhoton(pmap); #else return kdT_FirstPhoton(pmap); #endif } /* PHOTON MAP CLEANUP ROUTINES ------------------------------------------ */ void deletePhotons (PhotonMap* pmap) { #ifdef PMAP_OOC OOC_Delete(&pmap -> store); #else kdT_Delete(&pmap -> store); #endif free(pmap -> squeue.node); free(pmap -> biasCompHist); pmap -> numPhotons = pmap -> minGather = pmap -> maxGather = pmap -> squeue.len = pmap -> squeue.tail = 0; #ifdef PMAP_PHOTONFLOW if (pmap -> pathLUT) { lu_done(pmap -> pathLUT); free(pmap -> pathLUT); } if (pmap -> pathLUTKeys) { unsigned k; for (k = 0; k < pmap->squeue.len + 1; free(pmap->pathLUTKeys [k++])); free(pmap -> pathLUTKeys); } #endif }