diff --git a/pmapdata.c b/pmapdata.c index 3b18f51..68d1c18 100644 --- a/pmapdata.c +++ b/pmapdata.c @@ -1,993 +1,993 @@ #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", SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment") (c) Tokyo University of Science, supported by the JSPS Grants-in-Aid for Scientific Research (KAKENHI JP19KK0115, "Three-Dimensional Light Flow") (c) Tokyo University of Science, supported by the KAJIMA Foundation under the project title: "Reflections from Building Façades and their Glare Potential on the Built Environment -- Application of the Photon Flow Method using Annual Simulation" ========================================================================= $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 "pmaproi.h" #include "pmapcontrib.h" #include "otypes.h" #include "random.h" PhotonMap *photonMaps [NUM_PMAP_TYPES] = { NULL, NULL, NULL, NULL, NULL, NULL, NULL, #ifdef PMAP_PHOTONFLOW NULL, NULL, #endif NULL }; /* Include routines to handle underlying point cloud data structure */ #ifdef PMAP_OOC #include "pmapooc.c" #else #include "pmapkdt.c" #include "pmaptkdt.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_PATHFILT /* Init path lookup table and its key buffer */ pmap -> pathLUT = NULL; pmap -> pathLUTKeys = NULL; pmap -> numPathLUTKeys = 0; #endif /* Init stats */ setcolor(pmap -> photonFlux, 0, 0, 0); pmap -> CoG [0] = pmap -> CoG [1] = pmap -> CoG [2] = 0; pmap -> CoGdist = 0; pmap -> maxDist2 = 0; /* Init transient photon mapping stuff */ pmap -> velocity = photonVelocity; pmap -> minPathLen = pmap -> maxPathLen = pmap -> avgPathLen = 0; /* 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]; /* Init storage */ pmap -> heap = NULL; pmap -> heapBuf = NULL; pmap -> heapBufLen = 0; /* Init precomputed contrib photon stuff */ pmap -> contribMode = 0; pmap -> preCompContribTab = NULL; pmap -> contribHeap = NULL; pmap -> contribHeapBuf = NULL; pmap -> preCompContrib = NULL; pmap -> rcOpts = NULL; pmap -> numContribSrc = 0; pmap -> contribSrc = NULL; /* Mark photon contribution source as unused */ pmap -> lastContribSrc.srcIdx = pmap -> lastContribSrc.srcBin = -1; #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 W33nd0z? */ fdFlags = fcntl(fileno(pmap -> heap), F_GETFL); fcntl(fileno(pmap -> heap), F_SETFL, fdFlags | O_APPEND); #endif /* ftruncate(fileno(pmap -> heap), 0); */ } #ifdef PMAP_CONTRIB if (PMAP_CONTRIB_BINNING(pmap)) { /* Allocate precomputed contribution heap buffer */ if (!pmap -> contribHeap) { /* Open heap file for precomputed contributions */ - mktemp(strcpy(pmap -> contribHeapFname, PMAP_TMPFNAME)); + mktemp(strcpy(pmap -> contribHeapFname, PMAP_TMPFNAME)); if (!(pmap -> contribHeap = fopen(pmap -> contribHeapFname, "w+b"))) error(SYSTEM, "failed opening precomputed contribution " "heap file in initPhotonHeap()" ); #ifdef F_SETFL /* XXX is there an alternative needed for W33nd0z? */ fdFlags = fcntl(fileno(pmap -> contribHeap), F_GETFL); fcntl(fileno(pmap -> contribHeap), F_SETFL, fdFlags | O_APPEND); #endif /* ftruncate(fileno(pmap -> contribHeap), 0); */ } } #endif } void flushPhotonHeap (PhotonMap *pmap) { int fd, fdContrib = -1; unsigned long len, lenContrib = 0; if (!pmap) error(INTERNAL, "undefined photon map in flushPhotonHeap()"); if (!pmap -> heap || !pmap -> heapBuf) { /* Silently ignore undefined heap error(INTERNAL, "undefined photon heap in flushPhotonHeap()"); */ return; } fd = fileno(pmap -> heap); len = pmap -> heapBufLen * sizeof(Photon); #ifdef PMAP_CONTRIB if (PMAP_CONTRIB_BINNING(pmap)) { if (!pmap -> contribHeap || !pmap -> contribHeapBuf) { /* Silently ignore undefined heap error(INTERNAL, "undefined contribution heap in flushPhotonHeap()"); */ return; } /* Do same for the contribution photon buffa if binning is enabled. * Note the photon and contribution heap contain the same number of * entries to ensure photons are flushed in sync along with their * corresponding contributions to their respective heap files */ fdContrib = fileno(pmap -> contribHeap); lenContrib = pmap -> heapBufLen * ((PreComputedContrib*)pmap -> preCompContrib) -> contribSize; } #endif #if NIX /* Lock heap file(s) prior to writing */ shmLock(fd, F_WRLCK); if (lenContrib) shmLock(fdContrib, F_WRLCK); #endif #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); #if NIX fflush(stderr); #endif #endif /* Atomically write write block to photon heap file */ /* !!! Unbuffered I/O via pwrite() avoids potential race conditions * !!! and buffer corruption which can occur with lseek()/fseek() * !!! followed by write()/fwrite(). */ /*if (pwrite(fd, pmap -> heapBuf, len, lseek(fd, 0, SEEK_END)) != len) */ if (write(fd, pmap -> heapBuf, len) != len) { /* Clean up temp file */ fclose(pmap -> heap); unlink(pmap -> heapFname); error(SYSTEM, "can't append to photon heap file in flushPhotonHeap()"); } #ifdef PMAP_CONTRIB /* Atomically write block to contribution heap file */ if (fdContrib > -1 && write(fdContrib, pmap -> contribHeapBuf, lenContrib) != lenContrib ) { /* Clean up temp file */ fclose(pmap -> contribHeap); unlink(pmap -> contribHeapFname); error(SYSTEM, "can't append to contrib heap file in flushPhotonHeap()"); } #endif #if NIX if (fsync(fd) || fdContrib > -1 && fsync(fdContrib)) error(SYSTEM, "failed to fsync heap in flushPhotonHeap()"); /* Release lock on heap file(s) */ shmLock(fd, F_UNLCK); if (fdContrib > -1) shmLock(fdContrib, F_UNLCK); #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, const char *contribs) { 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; /* Store photon if within a region of interest (for ze Ecksperts!) */ if (!photonInROI(ray)) 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) ); #endif setPhotonFlux(&photon, photonFlux); /* Set photon position and flags */ VCOPY(photon.pos, ray -> rop); photon.flags = 0; photon.caustic = PMAP_CAUSTICRAY(ray); /* Set photon's subprocess index (to consolidate contribution sources after photon distribution completes) */ photon.proc = PMAP_GETRAYPROC(ray); switch (pmap -> type) { #ifdef PMAP_CONTRIB case PMAP_TYPE_CONTRIB: /* Contrib photon before precomputation (i.e. in forward pass); set contribution source from last index in contrib source array. Note the contribution source bin has already been set by newPhotonContribSrc(). */ photon.aux.contribSrc = pmap -> numContribSrc; /* Photon will be stored; set contribution source index, * thereby marking it as having spawned photon(s) */ if (ray -> rsrc > PMAP_MAXSRCIDX) error(INTERNAL, "contribution source index overflow"); else pmap -> lastContribSrc.srcIdx = ray -> rsrc; break; case PMAP_TYPE_CONTRIB_CHILD: /* Contrib photon being precomputed; set contribution source from index passed in ray. NOTE: This is currently only for info and consistency with the above, but not used by the lookup routines */ if (ray -> rsrc > PMAP_MAXSRCIDX) error(INTERNAL, "contribution source index overflow"); else photon.aux.contribSrc = ray -> rsrc; break; #endif case PMAP_TYPE_TRANSIENT: #ifdef PMAP_PHOTONFLOW case PMAP_TYPE_TRANSLIGHTFLOW: #endif /* Set cumulative path length for transient photon, corresponding to its time of flight. The cumulative path length is obtained relative to the maximum path length photonMaxPathDist. The last intersection distance is subtracted as this isn't factored in until photonRay() generates a new photon ray. */ photon.aux.pathLen = photonMaxPathDist - ray -> rmax + ray -> rot; break; default: /* Set photon's path ID from ray serial num (supplemented by photon.proc below) */ photon.aux.pathID = ray -> rno; } /* Set normal */ for (i = 0; i <= 2; i++) photon.norm [i] = 127.0 * (isVolumePmap(pmap) #ifdef PMAP_PHOTONFLOW || isLightFlowPmap(pmap) #endif ? ray -> rdir [i] : ray -> ron [i] ); if (!pmap -> heapBuf) { /* Lazily allocate photon heap buffa */ #if NIX /* Randomise buffa size to temporally decorellate flushes in * multiprocessing mode */ srandom(randSeed + getpid()); pmap -> heapBufSize = PMAP_HEAPBUFSIZE * 0.5 * (1 + frandom()); #else /* Randomisation disabled for single processes on WIN; also useful * for reproducability during debugging */ pmap -> heapBufSize = PMAP_HEAPBUFSIZE; #endif #ifdef PMAP_CONTRIB if (contribs && pmap -> preCompContrib) /* Adapt buffa size to accommodate precomputed contributions */ pmap -> heapBufSize *= (float)sizeof(Photon) / ((PreComputedContrib*)pmap -> preCompContrib) -> contribSize; #endif if (!(pmap -> heapBuf = calloc(pmap -> heapBufSize, sizeof(Photon)))) error(SYSTEM, "out of memory allocating heap buffer in newPhoton()" ); pmap -> heapBufLen = 0; } /* Photon initialised; now append to heap buffa */ memcpy(pmap -> heapBuf + pmap -> heapBufLen, &photon, sizeof(Photon)); #ifdef PMAP_CONTRIB if (contribs && pmap -> preCompContrib) { const unsigned long contribSize = ((PreComputedContrib*)pmap -> preCompContrib) -> contribSize; if (!pmap -> contribHeapBuf) { /* Lazily allocate contribution heap buffa of same size as photon * heap buffa */ pmap -> contribHeapBuf = calloc(pmap -> heapBufSize, contribSize); if (!pmap -> contribHeapBuf) error(SYSTEM, "out of memory allocating precomputed " "contribution heap buffer in newPhoton()" ); } /* Append precomputed contribs to heap buffa */ memcpy(pmap -> contribHeapBuf + pmap -> heapBufLen * contribSize, contribs, contribSize ); } #endif if (++pmap -> heapBufLen >= pmap -> heapBufSize) /* Photon (and contrib, if applicable) heap buffa(s) full, flush to heap file(s) */ 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, PhotonContribSourceIdx *contribSrcOfs, 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, avgPathLen = 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 sprintf(errmsg, "Reading %lu at %lu... ", pmap -> heapBufSize, ftell(pmap->heap) ); #endif pmap -> heapBufLen = fread(pmap -> heapBuf, sizeof(Photon), pmap -> heapBufSize, pmap -> heap ); #ifdef DEBUG_PMAP sprintf(errmsg, "%s Got %lu\n", errmsg, pmap -> heapBufLen); eputs(errmsg); #if NIX fflush(stderr); #endif #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 (isContribPmap(pmap) && contribSrcOfs) /* Linearise contrib photon origin index from subprocess index * using the per-subprocess offsets in contribSrcOfs */ p -> aux.contribSrc += contribSrcOfs [p -> proc]; else if (isTransientPmap(pmap) #ifdef PMAP_PHOTONFLOW || isTransLightFlowPmap(pmap) #endif ) { /* Update min and max path length */ if (p -> aux.pathLen < pmap -> minPathLen) pmap -> minPathLen = p -> aux.pathLen; else if (p -> aux.pathLen > pmap -> maxPathLen) pmap -> maxPathLen = p -> aux.pathLen; /* Update avg photon path length */ avgPathLen += p -> aux.pathLen; } /* 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, #if 1 photonFlux [isContribPmap(pmap) ? photonSrcIdx(pmap, p) : 0] #else photonFlux [0] #endif ); 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()" ); /* Clean up temp files */ fclose(pmap -> heap); fclose(nuHeap); unlink(pmap -> heapFname); unlink(nuHeapFname); } nCheck += pmap -> heapBufLen; } #if NIX /* !!! Flush pending writes to new heap !!! */ fflush(nuHeap); #endif #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; /* Average photon path lengths */ pmap -> avgPathLen = avgPathLen /= 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); if (isTransientPmap(pmap) #ifdef PMAP_PHOTONFLOW || isTransLightFlowPmap(pmap) #endif ) { /* Add temporal distance (in units of photon path length) for transient photons */ d [0] = p -> aux.pathLen - avgPathLen; CoGdist += d [0] * d [0]; } } } 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 if (isTransientPmap(pmap) #ifdef PMAP_PHOTONFLOW || isTransLightFlowPmap(pmap) #endif ) kdT_TransBuildPhotonMap(pmap); 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; const RREAL *norm = ray -> ron, *pos = ray -> rop; 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; if (isVolumePmap(pmap) #ifdef PMAP_PHOTONFLOW || isLightFlowPmap(pmap) && sphericalIrrad #endif ) { /* Volume photons are not associated with a normal or intersection point; null the normal and set origin as lookup point. */ /* If a lightflow photon map is used and sphericalIrrad = 1, the mean spherical irradiance is evaluated, so normals are irrelevant and not filtered. */ norm = NULL; pos = ray -> rorg; } /* 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_OOC OOC_FindPhotons(pmap, pos, norm); #else if (isTransientPmap(pmap) #ifdef PMAP_PHOTONFLOW || isTransLightFlowPmap(pmap) #endif ) kdT_TransFindPhotons(pmap, pos, norm); else kdT_FindPhotons(pmap, pos, norm); #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, void *photonIdx ) { /* 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, (OOC_DataIdx*)photonIdx ); #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_PATHFILT 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 /* Contrib stuff */ if (isContribPmap(pmap) && pmap -> preCompContribTab) { /* Parent contrib photon map; clean up child photon maps containing precomputed contribs via lu_done() */ lu_done(pmap -> preCompContribTab); free(pmap -> preCompContribTab); pmap -> preCompContribTab = NULL; if (pmap -> contribHeapBuf) { free(pmap-> contribHeapBuf); pmap -> contribHeapBuf = NULL; } if (pmap -> contribHeap) { fclose(pmap -> contribHeap); pmap -> contribHeap = NULL; } } } diff --git a/pmcontrib3.c b/pmcontrib3.c index bb44554..461dfc4 100644 --- a/pmcontrib3.c +++ b/pmcontrib3.c @@ -1,252 +1,265 @@ #ifndef lint static const char RCSid[] = "$Id$"; #endif /* ========================================================================= Photon map routines for precomputed light source contributions. These routines build and save precomputed contribution photon maps, and are used by mkpmap. Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) (c) Lucerne University of Applied Sciences and Arts, supported by the Swiss National Science Foundation (SNSF #147053, "Daylight Redirecting Components", SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment") (c) Tokyo University of Science, supported by the KAJIMA Foundation under the project title: "Reflections from Building Façades and their Glare Potential on the Built Environment -- Application of the Photon Flow Method using Annual Simulation" ========================================================================= $Id$ */ #include "pmapcontrib.h" #ifdef PMAP_CONTRIB #include "pmapdiag.h" #include "pmapio.h" #include #if NIX #define __USE_XOPEN_EXTENDED /* Enables nftw() under Linux */ #include #endif /* ------------------ CONTRIB PMAP BUILDING STUFF --------------------- */ #if NIX static int ftwFunc (const char *path, const struct stat *statData, int typeFlags, struct FTW *ftwData ) /* Callback for nftw() to delete a directory entry */ { return remove(path); /* = unlink()/rmdir() depending on type */ } #endif static int buildPreCompContribPmap (const LUENT *preCompContribNode, void *p) -/* Build per-modifier precomputed photon map. Returns 0 on success. */ +/* Build per-modifier precomputed photon map. Empty photo maps are + (nonfatally) skipped. Returns 0 on success. */ { PhotonMap *preCompContribPmap = (PhotonMap*)( preCompContribNode -> data ); PreComputedContrib *preCompContrib = (PreComputedContrib*)( preCompContribPmap -> preCompContrib ); const char *dirName = (char*)p; unsigned fnLen, wfnLen; /* Flush photon/contrib heaps */ flushPhotonHeap(preCompContribPmap); fflush(preCompContribPmap -> contribHeap); - if (verbose) { - sprintf(errmsg, "Building precomputed contribution photon map " - "for modifier %s\n", preCompContribNode -> key + /* Get heap file size, skip if zero (= no photons for this modifier) */ + if (fseek(preCompContribPmap -> heap, 0, SEEK_END) < 0) + error(SYSTEM, "failed seek to end of photon heap in buildPhotonMap()"); + if (!ftell(preCompContribPmap -> heap)) { + sprintf(errmsg, "No contribution photons precomputed for modifier %s", + preCompContribNode -> key ); - eputs(errmsg); + error(WARNING, errmsg); + } + else { + if (verbose) { + sprintf(errmsg, "Building precomputed contribution photon map " + "for modifier %s\n", preCompContribNode -> key + ); + eputs(errmsg); #if NIX - fflush(stderr); + fflush(stderr); #endif - } - - /* Allocate & set output filenames to subdirectory and modifier */ - fnLen = strlen(dirName) + strlen(preCompContribNode -> key) + - strlen(PMAP_CONTRIB_FILE) + 1; - wfnLen = strlen(dirName) + strlen(preCompContribNode -> key) + - strlen(PMAP_CONTRIB_WAVELETFILE) + 1; - if (!(preCompContribPmap -> fileName = malloc(fnLen)) || - !(preCompContrib -> waveletFname = malloc(wfnLen)) - ) - error(SYSTEM, "out of memory allocating filename in " - "buildPreCompContribPmap()" - ); - - snprintf(preCompContribPmap -> fileName, fnLen, - PMAP_CONTRIB_FILE, dirName, preCompContribNode -> key - ); - snprintf(preCompContrib -> waveletFname, wfnLen, - PMAP_CONTRIB_WAVELETFILE, dirName, preCompContribNode -> key - ); - - if (preCompContrib -> nBins > 1) { - /* Binning enabled; open wavelet coefficient file; this is where - the sorted contributions are written to by contribPhotonSort() - (via buildPhotonMap()) */ - if (!(preCompContrib -> waveletFile = fopen( - preCompContrib -> waveletFname, "wb" - )) - ) { - sprintf(errmsg, "can't open wavelet coefficient file %s", - preCompContrib -> waveletFname + } + + /* Allocate & set output filenames to subdirectory and modifier */ + fnLen = strlen(dirName) + strlen(preCompContribNode -> key) + + strlen(PMAP_CONTRIB_FILE) + 1; + wfnLen = strlen(dirName) + strlen(preCompContribNode -> key) + + strlen(PMAP_CONTRIB_WAVELETFILE) + 1; + if (!(preCompContribPmap -> fileName = malloc(fnLen)) || + !(preCompContrib -> waveletFname = malloc(wfnLen)) + ) + error(SYSTEM, "out of memory allocating filename in " + "buildPreCompContribPmap()" ); - error(SYSTEM, errmsg); + + snprintf(preCompContribPmap -> fileName, fnLen, + PMAP_CONTRIB_FILE, dirName, preCompContribNode -> key + ); + snprintf(preCompContrib -> waveletFname, wfnLen, + PMAP_CONTRIB_WAVELETFILE, dirName, preCompContribNode -> key + ); + + if (preCompContrib -> nBins > 1) { + /* Binning enabled; open wavelet coefficient file; this is where + the sorted contributions are written to by contribPhotonSort() + (via buildPhotonMap()) */ + if (!(preCompContrib -> waveletFile = fopen( + preCompContrib -> waveletFname, "wb" + )) + ) { + sprintf(errmsg, "can't open wavelet coefficient file %s", + preCompContrib -> waveletFname + ); + error(SYSTEM, errmsg); + } } + + /* Rebuild underlying data structure, destroying heap; photon flux is now + * finalised, so no flux normalisation (note NULL flux param) */ + buildPhotonMap(preCompContribPmap, NULL, NULL, 1); + + /* Free primary and transposed wavelet matrices; no longer needed */ + freeWaveletMatrix2(preCompContrib -> waveletMatrix, + preCompContrib -> coeffDim + ); + freeWaveletMatrix2(preCompContrib -> tWaveletMatrix, + preCompContrib -> coeffDim + ); + + /* Free thresholded coefficients; no longer needed */ + free(preCompContrib -> compressedCoeffs); + + preCompContrib -> waveletMatrix = + preCompContrib -> tWaveletMatrix = NULL; + preCompContrib -> compressedCoeffs = NULL; } - /* Rebuild underlying data structure, destroying heap; photon flux is now - * finalised, so no flux normalisation (note NULL flux param) */ - buildPhotonMap(preCompContribPmap, NULL, NULL, 1); - - /* Free primary and transposed wavelet matrices; no longer needed */ - freeWaveletMatrix2(preCompContrib -> waveletMatrix, - preCompContrib -> coeffDim - ); - freeWaveletMatrix2(preCompContrib -> tWaveletMatrix, - preCompContrib -> coeffDim - ); - - /* Free thresholded coefficients; no longer needed */ - free(preCompContrib -> compressedCoeffs); - - preCompContrib -> waveletMatrix = - preCompContrib -> tWaveletMatrix = NULL; - preCompContrib -> compressedCoeffs = NULL; - return 0; } void buildContribPhotonMap (const PhotonMap *pmap) /* Build contribution pmap and its per-modifier precomputed chilren */ { char dirName [1024]; struct stat dirStat; /* Set subdirectory containing per-modified child pmaps from filename */ snprintf(dirName, sizeof(dirName), PMAP_CONTRIB_DIR, pmap -> fileName); /* Cleanup previous directory contents if necessary */ if (!stat(dirName, &dirStat)) { /* File/dir exists */ if (S_ISREG(dirStat.st_mode)) /* Found regular file; delete and skip rest */ unlink(dirName); else if (S_ISDIR(dirStat.st_mode)) { /* Found subdirectory; delete its contents, skipping symlinks */ #if NIX if (nftw(dirName, ftwFunc, 1, FTW_DEPTH | FTW_MOUNT | FTW_PHYS) < 0 ) { sprintf(errmsg, "failed cleanup of output directory %s", dirName ); error(SYSTEM, errmsg); } #else /* Apparently there's no nftw() under Wind0z, so just skip * cleanup. Tuff luck, Wind0z Weeniez! There's probably some * replacement for this, but we couldn't be buggered... */ #endif } else { /* Found neither regular file nor directory; whatever it is, * just stuff it and quit! */ sprintf(errmsg, "cannot remove existing output file %s", dirName ); error(SYSTEM, errmsg); } } /* Create new directory for per-modifier contribution photon maps */ if (mkdir(dirName, S_IFDIR | S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH) < 0) { sprintf(errmsg, "error creating output directory %s", dirName); error(SYSTEM, errmsg); } - /* Build per-modifier precomputed pmaps from their contribution heaps. - !!! NOTE: buildPreCompContribPmap() bails out if any child contrib - !!! pmaps are empty; should this trigger a warning instead? */ + /* Build per-modifier precomputed pmaps from their contribution heaps. */ lu_doall(pmap -> preCompContribTab, buildPreCompContribPmap, dirName); } /* ------------------ CONTRIB PMAP SAVING STUFF --------------------- */ static int savePreCompContrib (const LUENT *preCompContribNode, void *p) +/* Save child contrib photon map for current modifier */ { PhotonMap *preCompContribPmap = (PhotonMap*)( preCompContribNode -> data ); PhotonMapArgz *argz = (PhotonMapArgz*)p; - /* Save child contrib photon map for current modifier; savePhotonMap() - * will not recurse and call saveContribPhotonMap() again because - * preCompContribPmap -> preCompContribTab == NULL */ - savePhotonMap(preCompContribPmap, preCompContribPmap -> fileName, - argz -> argc, argz -> argv - ); + /* Skip empty photon maps */ + if (preCompContribPmap -> numPhotons) + /* NOTE: savePhotonMap() will not recurse and call + saveContribPhotonMap() again because preCompContribPmap -> + preCompContribTab == NULL */ + savePhotonMap(preCompContribPmap, preCompContribPmap -> fileName, + argz -> argc, argz -> argv + ); return 0; } void saveContribPhotonMap (const PhotonMap *pmap, const char *fname, int argc, char **argv ) /* Save contribution pmap and its per-modifier precomputed children. NOTE: fname accepted to conform with savePhotonMap(), but currently ignored; this was already set in preComputedContrib() */ { PhotonMapArgz argz = {argc, argv}; char rcOptFname [1024]; FILE *rcOptFile; /* Save rcontrib options file */ if (!pmap -> rcOpts) error(INTERNAL, "undefined rcontrib options in saveContribPhotonMap()" ); if (!pmap -> fileName) error(INTERNAL, "undefined photon map filename in saveContribPhotonMap()" ); snprintf(rcOptFname, sizeof(rcOptFname), PMAP_CONTRIB_OPTFILE, pmap -> fileName ); if (!(rcOptFile = fopen(rcOptFname, "w")) || !fprintf(rcOptFile, "%s\n", pmap -> rcOpts) ) { sprintf(errmsg, "cannot write rcontrib options to %s", rcOptFname); error(SYSTEM, errmsg); } fclose(rcOptFile); /* Now iterate over per-modifier child pmaps */ lu_doall(pmap -> preCompContribTab, savePreCompContrib, &argz); } #endif /* PMAP_CONTRIB */ diff --git a/pmcontrib4.c b/pmcontrib4.c index 1cf6019..6b42031 100644 --- a/pmcontrib4.c +++ b/pmcontrib4.c @@ -1,645 +1,663 @@ #ifndef lint static const char RCSid[] = "$Id: pmcontrib2.c,v 2.5 2018/11/08 00:54:07 greg Exp $"; #endif /* ========================================================================= Photon map routines for precomputed light source contributions. These routines interface to rcontrib. Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) (c) Lucerne University of Applied Sciences and Arts, supported by the Swiss National Science Foundation (SNSF #147053, "Daylight Redirecting Components", SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment") + (c) Tokyo University of Science, + supported by the KAJIMA Foundation under the project title: + "Reflections from Building Façades and their Glare Potential on the + Built Environment -- Application of the Photon Flow Method using + Annual Simulation" ========================================================================= $Id$ */ #include "pmapcontrib.h" #ifdef PMAP_CONTRIB #include "pmapmat.h" #include "pmapsrc.h" #include "pmaprand.h" #include "pmapio.h" #include "pmapdiag.h" #include "otypes.h" #include "otspecial.h" #include "random.h" #include "pmcontribcache.h" #include "oocnn.h" void initPmapContribTab (LUTAB *contribTab, int *contribMode) /* Set contribution table (interface to rcmain.c) */ { pmapContribTab = contribTab; pmapContribMode = contribMode; } /* ------------------ CONTRIB PMAP LOADING STUFF --------------------- */ static int checkContribModifier (const LUENT *modContNode, void *p) /* Check for valid modifier in LUT entry */ { MODCONT *modCont = (MODCONT*)modContNode -> data; char *modName = (char*)modCont -> modname; OBJECT modObj = modifier(modName); const PhotonMap *pmap = (PhotonMap*)p; if (modObj == OVOID) { sprintf(errmsg, "invalid modifier %s", modName); error(USER, errmsg); } if (!islight(objptr(modObj) -> otype)) { sprintf(errmsg, "%s is not a light source modifier", modName); error(USER, errmsg); } return 0; } void initPmapContrib (PhotonMap *pmap) /* Set up photon map contributions and get binning parameters */ { unsigned t; /* Avoid incompatible photon maps */ for (t = 0; t < NUM_PMAP_TYPES; t++) if (photonMaps [t] && t != PMAP_TYPE_CONTRIB) { sprintf(errmsg, "%s photon map does not support contributions", pmapName [t] ); error(USER, errmsg); } if (!pmapContribTab) error(INTERNAL, "contribution LUT not initialised in initPmapContrib()" ); /* Check for valid contribution modifiers */ lu_doall(pmapContribTab, checkContribModifier, pmap); /* Set callback to collect contributions */ if (pmap) { pmap -> lookup = getPreCompPhotonContrib; /* XXX: Need tighter tolerance for contribution photon lookups? */ pmap -> gatherTolerance = 1.0; } } static int loadPreCompContrib (const LUENT *modContNode, void *p) /* Load child contribution photon map for current contrib modifier entry */ { MODCONT *modCont = (MODCONT*)modContNode -> data; char *modName = (char*)modCont -> modname; PhotonMap *preCompContribPmap; const PhotonMap *parentPmap = (const PhotonMap*)p; LUENT *preCompContribNode; PreComputedContrib *preCompContrib; - char dirName [1024]; + char dirName [1024], fileName [2048]; + struct stat fileStat; /* Allocate new LUT entry for child precomputed contribution pmap */ preCompContribNode = lu_find(parentPmap -> preCompContribTab, modName); if (!preCompContribNode) error(SYSTEM, "out of memory allocating LUT entry in loadPreCompContrib()" ); preCompContribNode -> key = modName; + /* Set subdirectory from parent photon map's filename, set child pmap's + filename from subdir and modifier */ + snprintf(dirName, sizeof(dirName), PMAP_CONTRIB_DIR, + parentPmap -> fileName + ); + snprintf(fileName, sizeof(dirName) + 4, PMAP_CONTRIB_FILE, + dirName, modName + ); + + /* Check if child photon map file exists, else skip */ + if (stat(fileName, &fileStat)) { + sprintf(errmsg, "no child photon map found for modifier %s, " + "contributions omitted", modName + ); + error(WARNING, errmsg); + preCompContribNode -> data = NULL; + return 0; + } + /* Allocate child precomputed contribution photon map */ preCompContribNode -> data = (char*)( preCompContribPmap = malloc(sizeof(PhotonMap)) ); if (!preCompContribPmap) error(SYSTEM, "out of memory allocating precomputed contribution " "photon map in loadPreCompContrib()" ); - /* Set subdirectory from parent photon map's filename */ - snprintf(dirName, sizeof(dirName), PMAP_CONTRIB_DIR, - parentPmap -> fileName - ); - /* Allocate & set child pmap's filename from subdir and modifier */ preCompContribPmap -> fileName = malloc(strlen(parentPmap -> fileName) + strlen(modName) + strlen(PMAP_CONTRIB_FILE) ); if (!preCompContribPmap -> fileName) error(SYSTEM, "out of memory allocating filename in " "loadPreCompContrib()" - ); - snprintf(preCompContribPmap -> fileName, sizeof(dirName) + 4, - PMAP_CONTRIB_FILE, dirName, modName - ); + ); + strcpy(preCompContribPmap -> fileName, fileName); /* Load child precomp contrib photon map for current modifier; loadPhotonMap() will not recurse and call loadContribPhotonMap() again because preCompContribPmap -> preCompContribTab == NULL */ - loadPhotonMap(preCompContribPmap, preCompContribPmap -> fileName); + loadPhotonMap(preCompContribPmap, preCompContribPmap -> fileName); /* Pass parent pmap's lookup parameters on to child */ preCompContribPmap -> minGather = parentPmap -> minGather; preCompContribPmap -> maxGather = parentPmap -> maxGather; - preCompContribPmap -> gatherTolerance = parentPmap -> gatherTolerance; + preCompContribPmap -> gatherTolerance = parentPmap -> gatherTolerance; /* Override contrib/coefficient mode if it doesn't match precomputed */ if (*pmapContribMode != preCompContribPmap -> contribMode) { sprintf(errmsg, "photon map contains precomputed %s, overriding rcontrib setting", preCompContribPmap -> contribMode ? "contributions" : "coefficients" ); error(WARNING, errmsg); *pmapContribMode = preCompContribPmap -> contribMode; } /* NOTE: preCompContribPmap -> preCompContrib is initialised by * loadPhotonMap() */ preCompContrib = (PreComputedContrib*)( preCompContribPmap -> preCompContrib ); /* Set number of bins and wavelet coefficients */ preCompContrib -> nBins = preCompContrib -> scDim * preCompContrib -> scDim; preCompContrib -> nCoeffs = preCompContrib -> coeffDim * preCompContrib -> coeffDim; /* Set wavelet coefficient filename */ preCompContrib -> waveletFname = malloc(strlen(dirName) + strlen(modName) + strlen(PMAP_CONTRIB_WAVELETFILE) ); if (!preCompContribPmap -> fileName) error(SYSTEM, "out of memory allocating filename in " "loadPreCompContrib()" ); snprintf(preCompContrib -> waveletFname, sizeof(dirName) + 5, PMAP_CONTRIB_WAVELETFILE, dirName, modName ); return 0; } void loadContribPhotonMap (PhotonMap *pmap, const char *fname) /* Load contribution pmap and its per-modifier precomputed children */ { LUTAB lutInit = LU_SINIT(NULL, freePreCompContribNode); /* Allocate & init LUT containing per-modifier child contrib pmaps */ if (!(pmap -> preCompContribTab = malloc(sizeof(lutInit)))) error(SYSTEM, "out of memory allocating precomputed contribution " "LUT in loadContribPhotonMap()" ); memcpy(pmap -> preCompContribTab, &lutInit, sizeof(lutInit)); /* NOTE: filename already set in loadPhotonMap() pmap -> fileName = (char*)fname; */ /* Init contrib photon map for light source contributions */ initPmapContrib(pmap); /* Load child contribution photon maps for each modifier in contribTab */ lu_doall(pmapContribTab, loadPreCompContrib, pmap); } /* ------------------ CONTRIB PMAP DECODING STUFF --------------------- */ static int decodeContribs (PreComputedContrib *preCompContrib) /* Decode and decompress mRGBE-encoded wavelet coefficients in preCompContrib -> mrgbeCoeffs, apply inverse wavelet transform and return decoded contributions in the wavelet coefficient matrix preCompContrib -> waveletMatrix. NOTE: THE WAVELET COEFFICIENT MATRIX IS ASSUMED TO BE ZEROED, with the wavelet approximation coefficients in the upper left of the matrix (i.e. preCompContrib -> waveletMatrix [0..approxDim-1][0..approxDim-1], where approxDim = WAVELET_PADD2_APPROXDIM). Returns 0 on success, else -1. */ { unsigned c, i, j, coeffIdx, scDim, coeffDim, nCoeffs, nCompressedCoeffs; WaveletMatrix2 waveletMatrix; mRGBERange *mrgbeRange; mRGBE *mrgbeCoeffs; DCOLOR fCoeff; if (!preCompContrib || !preCompContrib -> mrgbeCoeffs || !(waveletMatrix = preCompContrib -> waveletMatrix) || !preCompContrib -> tWaveletMatrix ) /* Should be initialised by caller! */ return -1; scDim = preCompContrib -> scDim; coeffDim = preCompContrib -> coeffDim; nCoeffs = preCompContrib -> nCoeffs; nCompressedCoeffs = preCompContrib -> nCompressedCoeffs; mrgbeCoeffs = preCompContrib -> mrgbeCoeffs; mrgbeRange = &preCompContrib -> mrgbeRange; /* Init mRGBE coefficient normalisation from range */ if (!mRGBEinit(mrgbeRange, mrgbeRange -> min, mrgbeRange -> max)) return -1; /* Decode mRGBE detail coeffs and populate wavelet coefficient matrix, based on the encoded incremental coefficient index. This omits the approximation coefficients in the upper left of the coefficient matrix (these should have already been set by the caller). NOTE: The detail coefficients in the matrix are assumed to be initialised to zero to account for those that were dropped by thresholding. */ for (c = coeffIdx = 0; c < nCompressedCoeffs; c++) { coeffIdx += mRGBEdecode(mrgbeCoeffs [c], mrgbeRange, fCoeff); i = PMAP_CONTRIB_LIN2X(coeffDim, coeffIdx); j = PMAP_CONTRIB_LIN2Y(coeffDim, coeffIdx); #ifdef PMAP_CONTRIB_DBG /* Check for invalid decoded coefficients */ if (i < WAVELET_PADD2_APPROXDIM && j < WAVELET_PADD2_APPROXDIM || coeffIdx >= nCoeffs ) error(CONSISTENCY, "wavelet coefficient index out of range in decodeContribs()" ); #endif copycolor(waveletMatrix [i] [j], fCoeff); } /* Do inverse 2D wavelet transform on preCompContrib -> waveletMatrix */ if (!padWaveletInvXform2(waveletMatrix, preCompContrib -> tWaveletMatrix, coeffDim, scDim ) ) error(INTERNAL, "failed inverse wavelet transform in decodeContribs()" ); #if 0 #ifdef PMAP_CONTRIB_DBG for (i = 0; i < scDim; i++) { for (j = 0; j < scDim; j++) { #if 0 /* Dump decoded contribs for debugging */ printf("% 7.3g\t", colorAvg(waveletMatrix [i][j])); #endif /* Warn if coefficient is negative; this _can_ happen due to loss of precision by the mRGBE encoding */ if (min(min(waveletMatrix [i][j][0], waveletMatrix [i][j][1]), waveletMatrix [i][j][2] ) < 0 ) error(WARNING, "negative contributions in getPreCompContrib()"); } #if 0 putchar('\n'); } putchar('\n'); #else } #endif #endif #endif return 0; } static void getPreCompContribByPhoton (const Photon *photon, OOC_DataIdx photonIdx, PreComputedContrib *preCompContrib, DCOLOR *binnedContribs ) /* Fetch and decode mRGBE-encoded wavelet coefficients for given photon from wavelet file at offset photonIdx, perform inverse wavelet xform, and return the reconstructed binned contributions in binnedContribs (which is assumed to be variable). Returns 0 on success, else -1. */ /* NOTE: binnedContribs IS ASSUMED TO POINT TO CACHE PAGE DATA! */ { unsigned i, j, k, coeffIdx, scDim, coeffDim, nCompressedCoeffs; COLOR fCoeff; COLR rgbe32 [WAVELET_PADD2_NUMAPPROX + 2]; WaveletMatrix2 waveletMatrix = NULL; WaveletCoeff3 *coeffPtr; if (!binnedContribs || !preCompContrib) /* Shouldn't happen */ error(INTERNAL, "contributions not initialised in getPreCompContrib()" ); if (preCompContrib -> nBins <= 1) /* No binning, so nothing to decode */ return; scDim = preCompContrib -> scDim; coeffDim = preCompContrib -> coeffDim; nCompressedCoeffs = preCompContrib -> nCompressedCoeffs; /* Lazily initialise preCompContrib */ if (!preCompContrib -> waveletFile) { /* Lazily open wavelet coefficient file */ preCompContrib -> waveletFile = fopen( preCompContrib -> waveletFname, "rb" ); if (!preCompContrib -> waveletFile) { sprintf(errmsg, "can't open wavelet coefficient file %s " "in getPreCompContrib()", preCompContrib -> waveletFname ); error(SYSTEM, errmsg); } /* Set record size of encoded coeffs in wavelet file */ preCompContrib -> contribSize = PMAP_CONTRIB_ENCSIZE( nCompressedCoeffs ); /* Lazily allocate buffer for mRGBE wavelet coeffs */ preCompContrib -> mrgbeCoeffs = calloc(nCompressedCoeffs, sizeof(mRGBE) ); if (!preCompContrib -> mrgbeCoeffs) error(SYSTEM, "out of memory allocating decoded wavelet " "coefficients in getPreCompContrib()" ); /* Lazily allocate primary and transposed wavelet matrices */ preCompContrib -> waveletMatrix = allocWaveletMatrix2(coeffDim); preCompContrib -> tWaveletMatrix = allocWaveletMatrix2(coeffDim); if (!preCompContrib -> waveletMatrix || !preCompContrib -> tWaveletMatrix ) error(SYSTEM, "out of memory allocating wavelet coefficient " "matrix in getPreCompContrib()" ); } /* Seek to photon index in wavelet file and read its associated coefficients */ if (fseek(preCompContrib -> waveletFile, photonIdx * preCompContrib -> contribSize, SEEK_SET ) < 0 ) error(SYSTEM, "can't seek in wavelet coefficient file " "in getPreCompContrib()" ); /* Read 32-bit encoded wavelet approximation coefficients and mRGBE * range; omit mRGBE minimum if only 1 compressed bin (=maximum * compression), since it's implicitly zero in this case. */ if (getbinary(rgbe32, sizeof(COLR), WAVELET_PADD2_NUMAPPROX + 1 + (nCompressedCoeffs > 1), preCompContrib -> waveletFile ) != WAVELET_PADD2_NUMAPPROX + 1 + (nCompressedCoeffs > 1) ) error(SYSTEM, "can't read approximation coefficients from " "wavelet coefficient file in getPreCompContrib()" ); if (getbinary(preCompContrib -> mrgbeCoeffs, sizeof(mRGBE), nCompressedCoeffs, preCompContrib -> waveletFile ) != nCompressedCoeffs ) error(SYSTEM, "can't read detail coefficients from " "wavelet coefficient file in getPreCompContrib()" ); /* Get mRGBE range (NOTE min/max are reversed in the coeff file) Note that direct assignment isn't possible as colr_color() returns float, not double. */ colr_color(fCoeff, rgbe32 [WAVELET_PADD2_NUMAPPROX]); copycolor(preCompContrib -> mrgbeRange.max, fCoeff); if (nCompressedCoeffs > 1) { colr_color(fCoeff, rgbe32 [WAVELET_PADD2_NUMAPPROX + 1]); copycolor(preCompContrib -> mrgbeRange.min, fCoeff); } else { /* Single compressed bin; mRGBE minimum is implicitly zero */ setcolor(preCompContrib -> mrgbeRange.min, 0, 0, 0); } /* Zero wavelet coefficient matrix and set approximation coefficients in * the upper left of the matrix, as expected by the decodeContribs(). */ waveletMatrix = preCompContrib -> waveletMatrix; zeroCoeffs2(waveletMatrix, coeffDim); for (i = 0; i < WAVELET_PADD2_APPROXDIM; i++) { /* Calc row index (=mult) in outer loop, increment in inner */ coeffIdx = PMAP_CONTRIB_XY2LIN(WAVELET_PADD2_APPROXDIM, i, 0); for (j = 0; j < WAVELET_PADD2_APPROXDIM; j++, coeffIdx++) { /* Direct assignment to wavelet matrix via colr_color() isn't * possible as the latter returns float, not double */ colr_color(fCoeff, rgbe32 [coeffIdx]); /* HACK: depending on the boundary extension mode, some approx * coeffs may be NEGATIVE (!), so set sign from extra bit in * 32-bit RGBE mantissa. */ for (k = 0; k < 3; k++) waveletMatrix [i] [j] [k] = PMAP_CONTRIB_GET_RGBE32_SGN( rgbe32 [coeffIdx] [k], fCoeff [k] ); } } /* All set, now decode mRGBE coeffs and invert wavelet transform */ if (decodeContribs(preCompContrib)) error(INTERNAL, "failed contribution decoding/decompression " "in getPreCompContrib()" ); /* Copy decoded contributions from wavelet coefficient matrix */ for (i = 0; i < scDim; i++) { /* Calc row pointer (=mult) in outer loop, increment in inner */ coeffPtr = binnedContribs + PMAP_CONTRIB_XY2LIN(scDim, i, 0); #ifdef PMAP_CONTRIB_LOG /* Logarithmic contribs; copy & apply log per element */ for (j = 0; j < scDim; j++, coeffPtr++) for (k = 0; k < 3; k++) (*coeffPtr) [k] = PMAP_CONTRIB_IVAL(waveletMatrix [i] [j] [k]); #else /* Linear contribs; copy matrix rows */ memcpy(*coeffPtr, waveletMatrix [i], scDim * WAVELET_COEFFSIZE); #endif } } typedef struct { const RAY *ray; RREAL rayCoeff [3]; COLORV *totalContrib; const PhotonMap *parentPmap; } PreCompContribRCData; static int getPreCompContribByMod (const LUENT *preCompContribTabNode, void *p ) /* Fetch and decode precomputed contributions from closest photon in pmap * referenced by preCompContribTabNode, and accumulate in pmapContribTab * for the corresponding modifier. * THIS VERSION FETCHES THE SINGLE CLOSEST PRECOMP. PHOTON */ { PhotonMap *preCompContribPmap = (PhotonMap*)( preCompContribTabNode -> data ); PreComputedContrib *preCompContrib; const char *modName = preCompContribTabNode -> key; MODCONT *modCont; PreCompContribRCData *rcData = (PreCompContribRCData*)p; LUENT *contribTabNode; Photon photon; OOC_DataIdx photonIdx; COLOR photonContrib; unsigned b, k; if (!preCompContribPmap || !rcData || !(preCompContrib = (PreComputedContrib*)( preCompContribPmap -> preCompContrib )) ) { sprintf(errmsg, "missing precomputed contributions for modifier %s", modName ); error(INTERNAL, errmsg); } /* Find rcontrib's LUT entry for this modifier */ contribTabNode = lu_find(pmapContribTab, modName); if (!contribTabNode || !contribTabNode -> key || !(modCont = (MODCONT*)contribTabNode -> data) ) { sprintf(errmsg, "missing contribution LUT entry for modifier %s", modName ); error(INTERNAL, errmsg); } /* Check if rcontrib bins match precomputed; bail out if not, as this should be set by the photon map's rcontrib options file */ if (modCont -> nbins != preCompContrib -> nBins) { sprintf(errmsg, "bin count mismatch for modifier %s (expected %d); " "missing option @" PMAP_CONTRIB_OPTFILE " ?", modName, preCompContrib -> nBins, rcData -> parentPmap -> fileName ); error(USER, errmsg); } /* Fetch closest photon and its contributions */ if (find1Photon(preCompContribPmap, rcData -> ray, &photon, &photonIdx ) ) { if (preCompContrib -> nBins > 1) { /* BINNING ENABLED; lazily init contrib cache if necessary */ if (initContribCache(preCompContrib)) { /* CACHE ENABLED; fetch cache page for found photon */ DCOLOR *cacheBins; if (!getContribCache(preCompContrib, photonIdx, &cacheBins)) { /* PAGE NOT CACHED; decode contribs into new cache page */ getPreCompContribByPhoton(&photon, photonIdx, preCompContrib, cacheBins ); } else; /* PAGE CACHED; (re)use decoded contribs from cache! */ for (b = 0; b < preCompContrib -> nBins; b++) { /* Accumulate binned contributions into rcontrib's LUT entry for this modifier, weighted by ray coefficient. Also sum total contributions. */ /* NOTE: Using multcolor() on the decoded contribs would modify them in the cache, so use a temp variable here. */ for (k = 0; k < 3; k++) photonContrib [k] = cacheBins [b] [k] * rcData -> rayCoeff [k]; addcolor(modCont -> cbin [b], photonContrib); addcolor(rcData -> totalContrib, photonContrib); } } else { /* CACHE DISABLED; decode contribs into temp array on stack */ DCOLOR tempBins [preCompContrib -> nBins]; getPreCompContribByPhoton(&photon, photonIdx, preCompContrib, tempBins ); for (b = 0; b < preCompContrib -> nBins; b++) { /* Accumulate binned contributions into rcontrib's LUT entry for this modifier, weighted by ray coefficient. Also sum total contributions. */ for (k = 0; k < 3; k++) photonContrib [k] = tempBins [b] [k] * rcData -> rayCoeff [k]; addcolor(modCont -> cbin [b], photonContrib); addcolor(rcData -> totalContrib, photonContrib); } } } else { /* NO BINNING; get single contribution directly from photon, scale by ray coefficient and sum to total contribs */ getPhotonFlux(&photon, photonContrib); multcolor(photonContrib, rcData -> rayCoeff); addcolor(modCont -> cbin [0], photonContrib); addcolor(rcData -> totalContrib, photonContrib); } } return 0; } void getPreCompPhotonContrib (PhotonMap *pmap, RAY *ray, COLOR totalContrib ) /* Get precomputed light source contributions at ray -> rop for all specified modifiers and accumulate them in pmapContribTab (which maps to rcontrib's binning LUT). Also return total precomputed contribs. */ { PreCompContribRCData rcData; /* Ignore sources */ if (ray -> ro && islight(objptr(ray -> ro -> omod) -> otype)) return; /* Get cumulative path coefficient up to photon lookup point */ raycontrib(rcData.rayCoeff, ray, PRIMARY); setcolor(totalContrib, 0, 0, 0); /* Rcontrib results passed to per-modifier child contrib pmaps */ rcData.ray = ray; rcData.totalContrib = totalContrib; rcData.parentPmap = pmap; /* Iterate over child contribution photon maps for each modifier and * collect their contributions */ lu_doall(pmap -> preCompContribTab, getPreCompContribByMod, &rcData); } #endif /* PMAP_CONTRIB */