diff --git a/pmapdata.c b/pmapdata.c index 8c3e64f..2540df8 100644 --- a/pmapdata.c +++ b/pmapdata.c @@ -1,836 +1,840 @@ #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") ========================================================================= $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, NULL #ifdef PMAP_PHOTONFLOW , NULL, NULL #endif }; /* 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 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]; /* Mark photon contribution source as unused */ pmap -> lastContribSrc.srcIdx = pmap -> lastContribSrc.srcBin = -1; pmap -> numContribSrc = 0; pmap -> contribSrc = 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()"); /* Clean up temp file */ fclose(pmap -> heap); unlink(pmap -> heapFname); } #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) ); #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 */ photon.proc = PMAP_GETRAYPROC(ray); switch (pmap -> type) { case PMAP_TYPE_CONTRIB: /* Set contrib photon's origin and subprocess index (the latter to * linearise the origin array after photon distribution completes). * Also set contribution source index, thereby marking it as used. * Note the contribution source bin has already been set by * newPhotonContribSrc(). */ photon.aux.contribSrc = pmap -> numContribSrc; pmap -> lastContribSrc.srcIdx = ray -> rsrc; break; 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. */ + #if 1 photon.aux.pathLen = photonMaxPathDist - ray -> rmax + ray -> rot; + #else + photon.aux.pathLen = photonMaxPathDist - ray -> rmax; + #endif 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 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 */ pmap -> heapBufSize = PMAP_HEAPBUFSIZE; #endif if (!(pmap -> heapBuf = calloc(pmap -> heapBufSize, sizeof(Photon)))) error(SYSTEM, "failed heap buffer allocation in newPhoton"); 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, 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 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 (isContribPmap(pmap) && contribSrcOfs) /* Linearise contrib photon origin index from subprocess index * using the per-subprocess offsets in contribSrcOfs */ p -> aux.contribSrc += contribSrcOfs [p -> proc]; /* 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, 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()" ); /* Clean up temp files */ fclose(pmap -> heap); fclose(nuHeap); unlink(pmap -> heapFname); unlink(nuHeapFname); } 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; /* 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) { /* 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_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 } diff --git a/pmapsrc.c b/pmapsrc.c index b96d926..4894888 100644 --- a/pmapsrc.c +++ b/pmapsrc.c @@ -1,1070 +1,1061 @@ #ifndef lint static const char RCSid[] = "$Id: pmapsrc.c,v 2.20 2021/04/09 17:42:34 rschregle Exp $"; #endif /* ====================================================================== Photon map support routines for emission from light sources 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: pmapsrc.c,v 2.20 2021/04/09 17:42:34 rschregle Exp $" */ #include "pmapsrc.h" #include "pmap.h" #include "pmaprand.h" #include "face.h" #include "otypes.h" #include "otspecial.h" #ifdef PMAP_PHOTONFLOW #include "pmapdens.h" #endif /* List of photon port modifier names */ char *photonPortList [MAXSET + 1] = {NULL}; /* Photon port objects (with modifiers in photonPortMods) */ SRCREC *photonPorts = NULL; unsigned numPhotonPorts = 0; int (*photonPartition [NUMOTYPE]) (EmissionMap*); int (*photonOrigin [NUMOTYPE]) (EmissionMap*); /* PHOTON PORT SUPPORT ROUTINES ------------------------------------------ */ /* Get/set photon port orientation flags from/in source flags. * HACK: the port orientation flags are embedded in the source flags and * shifted so they won't clobber the latter, since these are interpreted * by the *PhotonPartition() and *PhotonOrigin() routines! */ #define PMAP_SETPORTFLAGS(portdir) ((int)(portdir) << 14) #define PMAP_GETPORTFLAGS(sflags) ((sflags) >> 14) /* Set number of source partitions. * HACK: this is doubled if the source acts as a bidirectionally * emitting photon port, resulting in alternating front/backside partitions, * although essentially each partition is just used twice with opposing * normals. */ #define PMAP_SETNUMPARTITIONS(emap) ( \ (emap) -> numPartitions <<= ( \ (emap) -> port && \ PMAP_GETPORTFLAGS((emap) -> port -> sflags) == PMAP_PORTBI \ ) \ ) /* Get current source partition and numer of partitions * HACK: halve the number partitions if the source acts as a bidrectionally * emitting photon port, since each partition is used twice with opposing * normals. */ #define PMAP_GETNUMPARTITIONS(emap) (\ (emap) -> numPartitions >> ( \ (emap) -> port && \ PMAP_GETPORTFLAGS((emap) -> port -> sflags) == PMAP_PORTBI \ ) \ ) #define PMAP_GETPARTITION(emap) ( \ (emap) -> partitionCnt >> ( \ (emap) -> port && \ PMAP_GETPORTFLAGS((emap) -> port -> sflags) == PMAP_PORTBI \ ) \ ) void getPhotonPorts (char **portList) /* Find geometry declared as photon ports from modifiers in portList */ { OBJECT i; OBJREC *obj, *mat; int mLen; char **lp; /* Init photon port objects */ photonPorts = NULL; if (!portList [0]) return; for (i = numPhotonPorts = 0; i < nobjects; i++) { obj = objptr(i); mat = findmaterial(obj); /* Check if object is a surface and NOT a light source (duh) and * resolve its material (if any) via any aliases, then check for * inclusion in modifier list; note use of strncmp() to ignore port * flags */ if (issurface(obj -> otype) && mat && !islight(mat -> otype)) { mLen = strlen(mat -> oname); for (lp = portList; *lp && strncmp(mat -> oname, *lp, mLen); lp++); if (*lp) { /* Add photon port */ photonPorts = (SRCREC*)realloc( photonPorts, (numPhotonPorts + 1) * sizeof(SRCREC) ); if (!photonPorts) error(USER, "can't allocate photon ports"); photonPorts [numPhotonPorts].so = obj; /* Extract port orientation flags and embed in source flags. * Note setsource() combines (i.e. ORs) these with the actual * source flags below. */ photonPorts [numPhotonPorts].sflags = PMAP_SETPORTFLAGS((*lp) [mLen]); if (!sfun [obj -> otype].of || !sfun[obj -> otype].of -> setsrc) objerror(obj, USER, "illegal photon port"); setsource(photonPorts + numPhotonPorts, obj); numPhotonPorts++; } } } if (!numPhotonPorts) error(USER, "no valid photon ports found"); } static void setPhotonPortNormal (EmissionMap* emap) /* Set normal for current photon port partition (if defined) based on its * orientation */ { int i, portFlags; if (emap -> port) { /* Extract photon port orientation flags, set surface normal as follows: -- Port oriented forwards --> flip surface normal to point outwards, since normal points inwards per mkillum convention) -- Port oriented backwards --> surface normal is NOT flipped, since it already points inwards. -- Port is bidirectionally/bilaterally oriented --> flip normal based on the parity of the current partition emap -> partitionCnt. In this case, photon emission alternates between port front/back faces for consecutive partitions. */ portFlags = PMAP_GETPORTFLAGS(emap -> port -> sflags); if ( portFlags == PMAP_PORTFWD || portFlags == PMAP_PORTBI && !(emap -> partitionCnt & 1) ) for (i = 0; i < 3; i++) emap -> ws [i] = -emap -> ws [i]; } } /* SOURCE / PHOTON PORT PARTITIONING ROUTINES----------------------------- */ static int flatPhotonPartition2 ( EmissionMap* emap, unsigned long mp, FVECT cent, FVECT u, FVECT v, double du2, double dv2 ) /* Recursive part of flatPhotonPartition(..) */ { FVECT newct, newax; unsigned long npl, npu; if (mp > emap -> maxPartitions) { /* Enlarge partition array */ emap -> maxPartitions <<= 1; emap -> partitions = (unsigned char*)realloc( emap -> partitions, emap -> maxPartitions >> 1 ); if (!emap -> partitions) error(USER, "can't allocate source partitions"); memset( emap -> partitions + (emap -> maxPartitions >> 2), 0, emap -> maxPartitions >> 2 ); } if (du2 * dv2 <= 1) { /* hit limit? */ setpart(emap -> partitions, emap -> partitionCnt, S0); emap -> partitionCnt++; return 1; } if (du2 > dv2) { /* subdivide in U */ setpart(emap -> partitions, emap -> partitionCnt, SU); emap -> partitionCnt++; newax [0] = 0.5 * u [0]; newax [1] = 0.5 * u [1]; newax [2] = 0.5 * u [2]; u = newax; du2 *= 0.25; } else { /* subdivide in V */ setpart(emap -> partitions, emap -> partitionCnt, SV); emap -> partitionCnt++; newax [0] = 0.5 * v [0]; newax [1] = 0.5 * v [1]; newax [2] = 0.5 * v [2]; v = newax; dv2 *= 0.25; } /* lower half */ newct [0] = cent [0] - newax [0]; newct [1] = cent [1] - newax [1]; newct [2] = cent [2] - newax [2]; npl = flatPhotonPartition2(emap, mp << 1, newct, u, v, du2, dv2); /* upper half */ newct [0] = cent [0] + newax [0]; newct [1] = cent [1] + newax [1]; newct [2] = cent [2] + newax [2]; npu = flatPhotonPartition2(emap, mp << 1, newct, u, v, du2, dv2); /* return total */ return npl + npu; } static int flatPhotonPartition (EmissionMap* emap) /* Partition flat source for photon emission */ { RREAL *vp; double du2, dv2; memset(emap -> partitions, 0, emap -> maxPartitions >> 1); emap -> partArea = srcsizerat * thescene.cusize; emap -> partArea *= emap -> partArea; vp = emap -> src -> ss [SU]; du2 = DOT(vp, vp) / emap -> partArea; vp = emap -> src -> ss [SV]; dv2 = DOT(vp, vp) / emap -> partArea; emap -> partitionCnt = 0; emap -> numPartitions = flatPhotonPartition2( emap, 1, emap -> src -> sloc, emap -> src -> ss [SU], emap -> src -> ss [SV], du2, dv2 ); emap -> partitionCnt = 0; emap -> partArea = emap -> src -> ss2 / emap -> numPartitions; return 1; } static int sourcePhotonPartition (EmissionMap* emap) /* Partition scene cube faces or photon port for photon emission from distant source */ { if (emap -> port) { /* Relay partitioning to photon port */ SRCREC *src = emap -> src; emap -> src = emap -> port; photonPartition [emap -> src -> so -> otype] (emap); PMAP_SETNUMPARTITIONS(emap); emap -> src = src; } else { /* No photon ports defined; partition scene cube faces (SUBOPTIMAL) */ memset(emap -> partitions, 0, emap -> maxPartitions >> 1); setpart(emap -> partitions, 0, S0); emap -> partitionCnt = 0; emap -> numPartitions = 1 / srcsizerat; emap -> numPartitions *= emap -> numPartitions; emap -> partArea = sqr(thescene.cusize) / emap -> numPartitions; emap -> numPartitions *= 6; } return 1; } static int spherePhotonPartition (EmissionMap* emap) /* Partition spherical source into equal solid angles using uniform mapping for photon emission */ { unsigned numTheta, numPhi; memset(emap -> partitions, 0, emap -> maxPartitions >> 1); setpart(emap -> partitions, 0, S0); emap -> partArea = 4 * PI * sqr(emap -> src -> srad); emap -> numPartitions = emap -> partArea / sqr(srcsizerat * thescene.cusize); numTheta = max(sqrt(2 * emap -> numPartitions / PI) + 0.5, 1); numPhi = 0.5 * PI * numTheta + 0.5; emap -> numPartitions = (unsigned long)numTheta * numPhi; emap -> partitionCnt = 0; emap -> partArea /= emap -> numPartitions; return 1; } static int cylPhotonPartition2 ( EmissionMap* emap, unsigned long mp, FVECT cent, FVECT axis, double d2 ) /* Recursive part of cyPhotonPartition(..) */ { FVECT newct, newax; unsigned long npl, npu; if (mp > emap -> maxPartitions) { /* Enlarge partition array */ emap -> maxPartitions <<= 1; emap -> partitions = (unsigned char*)realloc( emap -> partitions, emap -> maxPartitions >> 1 ); if (!emap -> partitions) error(USER, "can't allocate source partitions"); memset( emap -> partitions + (emap -> maxPartitions >> 2), 0, emap -> maxPartitions >> 2 ); } if (d2 <= 1) { /* hit limit? */ setpart(emap -> partitions, emap -> partitionCnt, S0); emap -> partitionCnt++; return 1; } /* subdivide */ setpart(emap -> partitions, emap -> partitionCnt, SU); emap -> partitionCnt++; newax [0] = 0.5 * axis [0]; newax [1] = 0.5 * axis [1]; newax [2] = 0.5 * axis [2]; d2 *= 0.25; /* lower half */ newct [0] = cent [0] - newax [0]; newct [1] = cent [1] - newax [1]; newct [2] = cent [2] - newax [2]; npl = cylPhotonPartition2(emap, mp << 1, newct, newax, d2); /* upper half */ newct [0] = cent [0] + newax [0]; newct [1] = cent [1] + newax [1]; newct [2] = cent [2] + newax [2]; npu = cylPhotonPartition2(emap, mp << 1, newct, newax, d2); /* return total */ return npl + npu; } static int cylPhotonPartition (EmissionMap* emap) /* Partition cylindrical source for photon emission */ { double d2; memset(emap -> partitions, 0, emap -> maxPartitions >> 1); d2 = srcsizerat * thescene.cusize; d2 = PI * emap -> src -> ss2 / (2 * emap -> src -> srad * sqr(d2)); d2 *= d2 * DOT(emap -> src -> ss [SU], emap -> src -> ss [SU]); emap -> partitionCnt = 0; emap -> numPartitions = cylPhotonPartition2( emap, 1, emap -> src -> sloc, emap -> src -> ss [SU], d2 ); emap -> partitionCnt = 0; emap -> partArea = PI * emap -> src -> ss2 / emap -> numPartitions; return 1; } /* PHOTON ORIGIN ROUTINES ------------------------------------------------ */ static int flatPhotonOrigin (EmissionMap* emap) /* Init emission map with photon origin and associated surface axes on flat light source surface. Also sets source aperture and sampling hemisphere axes at origin. Returns 1 if photon lies within source boundary, else 0 (signalling rejection to trim non-quadrilateral geometry). Some code fragments are adapted from srcsamp.c:nextssamp(). */ { int i, cent[3], size[3], parr[2]; FVECT vpos; const SRCREC *src = emap -> src; /* Get surface axes */ VCOPY(emap -> us, src -> ss [SU]); normalize(emap -> us); VCOPY(emap -> ws, src -> ss [SW]); /* Flip normal emap -> ws if port and required by its orientation */ setPhotonPortNormal(emap); fcross(emap -> vs, emap -> ws, emap -> us); /* Get hemisphere axes & aperture */ if (src -> sflags & SSPOT) { VCOPY(emap -> wh, src -> sl.s -> aim); i = 0; do { emap -> vh [0] = emap -> vh [1] = emap -> vh [2] = 0; emap -> vh [i++] = 1; fcross(emap -> uh, emap -> vh, emap -> wh); } while (normalize(emap -> uh) < FTINY); fcross(emap -> vh, emap -> wh, emap -> uh); emap -> cosThetaMax = 1 - src -> sl.s -> siz / (2 * PI); } else { VCOPY(emap -> uh, emap -> us); VCOPY(emap -> vh, emap -> vs); VCOPY(emap -> wh, emap -> ws); emap -> cosThetaMax = 0; } cent [SU] = cent [SV] = cent [SW] = 0; size [SU] = size [SV] = size [SW] = emap -> maxPartitions; parr [0] = 0; parr [1] = PMAP_GETPARTITION(emap); if (!skipparts(cent, size, parr, emap -> partitions)) error(CONSISTENCY, "bad source partition in flatPhotonOrigin"); vpos [SU] = (1 - 2 * pmapRandom(partState)) * size [SU] / emap -> maxPartitions; vpos [SV] = (1 - 2 * pmapRandom(partState)) * size [SV] / emap -> maxPartitions; vpos [SW] = 0; for (i = 0; i < 3; i++) vpos [i] += (double)cent [i] / emap -> maxPartitions; if (src -> sflags & SCIR) { /* Disk source; trim origin by correcting partition extent (note implicit trim [SW] = 0) */ FVECT trim; trim [SU] = 1.12837917 * sqrt(1 - 0.5 * vpos [SV] * vpos [SV]); trim [SV] = 1.12837917 * sqrt(1 - 0.5 * vpos [SU] * vpos [SU]); for (i = 0; i < 3; i++) emap -> photonOrg [i] = src -> sloc [i] + trim [SU] * vpos [SU] * src -> ss [SU][i] + trim [SV] * vpos [SV] * src -> ss [SV][i]; /* Accept origin unconditionally */ return 1; } else { /* Polygonal source */ FACE *inYaFace = getface(src -> so); if (inYaFace -> nv == 3) { /* Triangle; sample separately; code adapted from rfluxmtx.c */ vpos [SU] = 0.5 * (vpos [SU] + 1); vpos [SV] = 0.5 * (vpos [SV] + 1); vpos [SU] *= vpos [SV] = sqrt(vpos [SV]); vpos [SV] = 1 - vpos [SV]; for (i = 0; i < 3; i++) { const RREAL v0 = VERTEX(inYaFace, 0) [i]; emap -> photonOrg [i] = v0 + vpos [SU] * (VERTEX(inYaFace, 1) [i] - v0) + vpos [SV] * (VERTEX(inYaFace, 2) [i] - v0); } /* Accept photon origin unconditionally */ return 1; } else if (inYaFace -> nv == 4) { /* Quadrilateral; use standard sampling from srcsamp.c */ for (i = 0; i < 3; i++) emap -> photonOrg [i] = emap -> src -> sloc [i] + vpos [SU] * emap -> src -> ss [SU][i] + vpos [SV] * emap -> src -> ss [SV][i] + vpos [SW] * emap -> src -> ss [SW][i]; /* Accept photon origin unconditionally */ return 1; } else { /* Arbitrary polygon */ for (i = 0; i < 3; i++) emap -> photonOrg [i] = emap -> src -> sloc [i] + vpos [SU] * src -> srad * emap -> us [i] + vpos [SV] * src -> srad * emap -> vs [i]; /* Accept origin only if within polygon bounds (= rejection sampling) */ return inface(emap -> photonOrg, inYaFace); } } } static int spherePhotonOrigin (EmissionMap* emap) /* Init emission map with photon origin and associated surface axes on spherical light source. Also sets source aperture and sampling hemisphere axes at origin */ { int i = 0; unsigned numTheta, numPhi, t, p; RREAL cosTheta, sinTheta, phi; /* Get current partition */ numTheta = max(sqrt(2 * PMAP_GETNUMPARTITIONS(emap) / PI) + 0.5, 1); numPhi = 0.5 * PI * numTheta + 0.5; t = PMAP_GETPARTITION(emap) / numPhi; p = PMAP_GETPARTITION(emap) - t * numPhi; emap -> ws [2] = cosTheta = 1 - 2 * (t + pmapRandom(partState)) / numTheta; sinTheta = sqrt(1 - sqr(cosTheta)); phi = 2 * PI * (p + pmapRandom(partState)) / numPhi; emap -> ws [0] = cos(phi) * sinTheta; emap -> ws [1] = sin(phi) * sinTheta; /* Flip normal emap -> ws if port and required by its orientation */ setPhotonPortNormal(emap); /* Get surface axes us & vs perpendicular to ws */ do { emap -> vs [0] = emap -> vs [1] = emap -> vs [2] = 0; emap -> vs [i++] = 1; fcross(emap -> us, emap -> vs, emap -> ws); } while (normalize(emap -> us) < FTINY); fcross(emap -> vs, emap -> ws, emap -> us); /* Get origin */ for (i = 0; i < 3; i++) emap -> photonOrg [i] = emap -> src -> sloc [i] + emap -> src -> srad * emap -> ws [i]; /* Get hemisphere axes & aperture */ if (emap -> src -> sflags & SSPOT) { VCOPY(emap -> wh, emap -> src -> sl.s -> aim); i = 0; do { emap -> vh [0] = emap -> vh [1] = emap -> vh [2] = 0; emap -> vh [i++] = 1; fcross(emap -> uh, emap -> vh, emap -> wh); } while (normalize(emap -> uh) < FTINY); fcross(emap -> vh, emap -> wh, emap -> uh); emap -> cosThetaMax = 1 - emap -> src -> sl.s -> siz / (2 * PI); } else { VCOPY(emap -> uh, emap -> us); VCOPY(emap -> vh, emap -> vs); VCOPY(emap -> wh, emap -> ws); emap -> cosThetaMax = 0; } return 1; } static int sourcePhotonOrigin (EmissionMap* emap) /* Init emission map with photon origin and associated surface axes on scene cube face for distant light source. Also sets source aperture (solid angle) and sampling hemisphere axes at origin */ { unsigned long i, partsPerDim, partsPerFace; unsigned face; RREAL du, dv; int acceptOrigin = 1; if (emap -> port) { /* Relay to photon port; get origin on its surface */ SRCREC *src = emap -> src; emap -> src = emap -> port; acceptOrigin = photonOrigin [emap -> src -> so -> otype] (emap); emap -> src = src; } else { /* No ports defined, so get origin on scene cube face (SUBOPTIMAL) */ /* Get current face from partition number */ partsPerDim = 1 / srcsizerat; partsPerFace = sqr(partsPerDim); face = emap -> partitionCnt / partsPerFace; if (!(emap -> partitionCnt % partsPerFace)) { /* Skipped to a new face; get its normal */ emap -> ws [0] = emap -> ws [1] = emap -> ws [2] = 0; emap -> ws [face >> 1] = face & 1 ? 1 : -1; /* Get surface axes us & vs perpendicular to ws */ face >>= 1; emap -> vs [0] = emap -> vs [1] = emap -> vs [2] = 0; emap -> vs [(face + (emap -> ws [face] > 0 ? 2 : 1)) % 3] = 1; fcross(emap -> us, emap -> vs, emap -> ws); } /* Get jittered offsets within face from partition number (in range [-0.5, 0.5]) */ i = emap -> partitionCnt % partsPerFace; du = (i / partsPerDim + pmapRandom(partState)) / partsPerDim - 0.5; dv = (i % partsPerDim + pmapRandom(partState)) / partsPerDim - 0.5; /* Jittered destination point within partition */ for (i = 0; i < 3; i++) emap -> photonOrg [i] = thescene.cuorg [i] + thescene.cusize * ( 0.5 + du * emap -> us [i] + dv * emap -> vs [i] + 0.5 * emap -> ws [i] ); } /* Get hemisphere axes & aperture */ VCOPY(emap -> wh, emap -> src -> sloc); i = 0; do { emap -> vh [0] = emap -> vh [1] = emap -> vh [2] = 0; emap -> vh [i++] = 1; fcross(emap -> uh, emap -> vh, emap -> wh); } while (normalize(emap -> uh) < FTINY); fcross(emap -> vh, emap -> wh, emap -> uh); /* Get aperture */ emap -> cosThetaMax = 1 - emap -> src -> ss2 / (2 * PI); emap -> cosThetaMax = min(1, max(-1, emap -> cosThetaMax)); return acceptOrigin; } static int cylPhotonOrigin (EmissionMap* emap) /* Init emission map with photon origin and associated surface axes on cylindrical light source surface. Also sets source aperture and sampling hemisphere axes at origin */ { int i, cent[3], size[3], parr[2]; FVECT v; cent [0] = cent [1] = cent [2] = 0; size [0] = size [1] = size [2] = emap -> maxPartitions; parr [0] = 0; parr [1] = PMAP_GETPARTITION(emap); if (!skipparts(cent, size, parr, emap -> partitions)) error(CONSISTENCY, "bad source partition in cylPhotonOrigin"); v [SU] = 0; v [SV] = (1 - 2 * pmapRandom(partState)) * (double)size [1] / emap -> maxPartitions; v [SW] = (1 - 2 * pmapRandom(partState)) * (double)size [2] / emap -> maxPartitions; normalize(v); v [SU] = (1 - 2 * pmapRandom(partState)) * (double)size [1] / emap -> maxPartitions; for (i = 0; i < 3; i++) v [i] += (double)cent [i] / emap -> maxPartitions; /* Get surface axes */ for (i = 0; i < 3; i++) emap -> photonOrg [i] = emap -> ws [i] = ( v [SV] * emap -> src -> ss [SV][i] + v [SW] * emap -> src -> ss [SW][i] ) / 0.8559; /* Flip normal emap -> ws if port and required by its orientation */ setPhotonPortNormal(emap); normalize(emap -> ws); VCOPY(emap -> us, emap -> src -> ss [SU]); normalize(emap -> us); fcross(emap -> vs, emap -> ws, emap -> us); /* Get origin */ for (i = 0; i < 3; i++) emap -> photonOrg [i] += v [SU] * emap -> src -> ss [SU][i] + emap -> src -> sloc [i]; /* Get hemisphere axes & aperture */ if (emap -> src -> sflags & SSPOT) { VCOPY(emap -> wh, emap -> src -> sl.s -> aim); i = 0; do { emap -> vh [0] = emap -> vh [1] = emap -> vh [2] = 0; emap -> vh [i++] = 1; fcross(emap -> uh, emap -> vh, emap -> wh); } while (normalize(emap -> uh) < FTINY); fcross(emap -> vh, emap -> wh, emap -> uh); emap -> cosThetaMax = 1 - emap -> src -> sl.s -> siz / (2 * PI); } else { VCOPY(emap -> uh, emap -> us); VCOPY(emap -> vh, emap -> vs); VCOPY(emap -> wh, emap -> ws); emap -> cosThetaMax = 0; } return 1; } /* PHOTON EMISSION ROUTINES ---------------------------------------------- */ static int defaultEmissionFunc (EmissionMap* emap) /* Default behaviour when no emission funcs defined for this source type */ { objerror(emap -> src -> so, INTERNAL, "undefined photon emission function" ); return 0; } void initPhotonEmissionFuncs () /* Init photonPartition[] and photonOrigin[] dispatch tables */ { int i; for (i = 0; i < NUMOTYPE; i++) photonPartition [i] = photonOrigin [i] = defaultEmissionFunc; photonPartition [OBJ_FACE] = photonPartition [OBJ_RING] = flatPhotonPartition; photonPartition [OBJ_SOURCE] = sourcePhotonPartition; photonPartition [OBJ_SPHERE] = spherePhotonPartition; photonPartition [OBJ_CYLINDER] = cylPhotonPartition; photonOrigin [OBJ_FACE] = photonOrigin [OBJ_RING] = flatPhotonOrigin; photonOrigin [OBJ_SOURCE] = sourcePhotonOrigin; photonOrigin [OBJ_SPHERE] = spherePhotonOrigin; photonOrigin [OBJ_CYLINDER] = cylPhotonOrigin; } void initPhotonEmission (EmissionMap *emap, float numPdfSamples) /* Initialise photon emission from partitioned light source emap -> src; * this involves integrating the flux emitted from the current photon * origin emap -> photonOrg and setting up a PDF to sample the emission * distribution with numPdfSamples samples */ { unsigned i, t, p; double phi, cosTheta, sinTheta, du, dv, dOmega; EmissionSample* sample; const OBJREC* mod = findmaterial(emap -> src -> so); static RAY r; photonOrigin [emap -> src -> so -> otype] (emap); setcolor(emap -> partFlux, 0, 0, 0); emap -> cdf = 0; emap -> numSamples = 0; cosTheta = DOT(emap -> ws, emap -> wh); if (cosTheta <= 0 && sqrt(1-sqr(cosTheta)) <= emap -> cosThetaMax + FTINY) /* Aperture completely below surface; no emission from current origin */ return; if (mod -> omod == OVOID && !emap -> port && (cosTheta >= 1 - FTINY || ( emap -> src -> sflags & SDISTANT && acos(cosTheta) + acos(emap -> cosThetaMax) <= 0.5 * PI ))) { /* Source is unmodified and has no port (which requires testing for occlusion), and is either local with normal aligned aperture or distant with aperture above surface --> get flux & PDF via analytical solution */ setcolor(emap -> partFlux, mod -> oargs.farg [0], mod -> oargs.farg [1], mod -> oargs.farg [2] ); /* Multiply radiance by projected Omega * dA to get flux */ scalecolor(emap -> partFlux, PI * cosTheta * (1 - sqr(max(emap -> cosThetaMax, 0))) * emap -> partArea ); } else { /* Source is either modified, has a port, is local with off-normal aperture, or distant with aperture partly below surface --> get flux via numerical integration using aperture samples with constant differential solid angle */ /* Figga out numba of aperture samples for integration; numTheta / numPhi ratio is 1 / PI */ emap -> dCosTheta = (1 - emap -> cosThetaMax); du = sqrt(pdfSamples * 2 * emap -> dCosTheta); emap -> numTheta = max(du + 0.5, 1); emap -> numPhi = max(PI * du + 0.5, 1); emap -> dCosTheta /= emap -> numTheta; emap -> dPhi = 2 * PI / emap -> numPhi; /* Differential solid angle is constant for this mapping */ dOmega = emap -> dCosTheta * emap -> dPhi; /* Allocate PDF, baby */ sample = emap -> samples = (EmissionSample*)realloc(emap -> samples, sizeof(EmissionSample) * emap -> numTheta * emap -> numPhi ); if (!emap -> samples) error(USER, "can't allocate emission PDF"); VCOPY(r.rorg, emap -> photonOrg); VCOPY(r.rop, emap -> photonOrg); r.rmax = 0; for (t = 0; t < emap -> numTheta; t++) { for (p = 0; p < emap -> numPhi; p++) { /* This uniform mapping handles 0 <= thetaMax <= PI */ cosTheta = 1 - (t + pmapRandom(emitState)) * emap -> dCosTheta; sinTheta = sqrt(1 - sqr(cosTheta)); phi = (p + pmapRandom(emitState)) * emap -> dPhi; du = cos(phi) * sinTheta; dv = sin(phi) * sinTheta; rayorigin(&r, PRIMARY, NULL, NULL); for (i = 0; i < 3; i++) r.rdir [i] = (du * emap -> uh [i] + dv * emap -> vh [i] + cosTheta * emap -> wh [i] ); /* Sample behind surface? */ VCOPY(r.ron, emap -> ws); if ((r.rod = DOT(r.rdir, r.ron)) <= 0) continue; /* Get radiance emitted in this direction; to get flux we multiply by cos(theta_surface), dOmega, and dA. Ray is directed towards light source for raytexture(). */ if (!(emap -> src -> sflags & SDISTANT)) for (i = 0; i < 3; i++) r.rdir [i] = -r.rdir [i]; /* Port occluded in this direction? */ if (emap -> port && localhit(&r, &thescene)) continue; raytexture(&r, mod -> omod); setcolor(r.rcol, mod -> oargs.farg [0], mod -> oargs.farg [1], mod -> oargs.farg [2] ); multcolor(r.rcol, r.pcol); /* Multiply by projection cos(theta_surface) */ scalecolor(r.rcol, r.rod); /* Add PDF sample if nonzero. TODO: importance info for photon emission could go here.*/ if (colorAvg(r.rcol)) { copycolor(sample -> pdf, r.rcol); sample -> cdf = emap -> cdf += colorAvg(sample -> pdf); sample -> theta = t; sample++ -> phi = p; emap -> numSamples++; addcolor(emap -> partFlux, r.rcol); } } } /* Multiply by differential angle & area, dOmega * dA */ scalecolor(emap -> partFlux, dOmega * emap -> partArea); } } #define vomitPhoton emitPhoton #define bluarrrghPhoton vomitPhoton int emitPhoton (const EmissionMap* emap, RAY* ray) /* Emit photon from current partition emap -> partitionCnt based on emission distribution. Returns new photon ray and an acceptance/rejection flag: 1 if origin lies within the emitting source/port, else 0 if it lies outside, in which case the emitted photon should be rejected. */ { unsigned long i, lo, hi; const EmissionSample* sample = emap -> samples; RREAL du, dv, cosTheta, cosThetaSqr, sinTheta, phi; const OBJREC* mod = findmaterial(emap -> src -> so); /* Choose a new origin within current partition for every emitted photon to break up clustering artifacts */ if (!photonOrigin [emap -> src -> so -> otype] ((EmissionMap*)emap)) /* Origin lies outside source/port geometry; signal its rejection */ return 0; /* If we have a local glow source with a maximum radius, then restrict our photon to the specified distance, otherwise we set the limit imposed by photonMaxPathDist (or no limit if 0) */ if (mod -> otype == MAT_GLOW && !(emap -> src -> sflags & SDISTANT) && mod -> oargs.farg[3] > FTINY ) ray -> rmax = mod -> oargs.farg[3]; else ray -> rmax = photonMaxPathDist; rayorigin(ray, PRIMARY, NULL, NULL); /* Init directional samples */ if (!emap -> numSamples) { /* Source is unmodified and has no port, and either local with normal aligned aperture, or distant with aperture above surface --> use cosine weighted distribution */ cosThetaSqr = (1 - pmapRandom(emitState) * (1 - sqr(max(emap -> cosThetaMax, 0))) ); cosTheta = sqrt(cosThetaSqr); sinTheta = sqrt(1 - cosThetaSqr); phi = 2 * PI * pmapRandom(emitState); setcolor( ray -> rcol, mod -> oargs.farg [0], mod -> oargs.farg [1], mod -> oargs.farg [2] ); } else { /* Source is either modified, has a port, is local with off-normal aperture, or distant with aperture partly below surface --> choose direction from constructed cumulative distribution function with Monte Carlo inversion using binary search. */ du = pmapRandom(emitState) * emap -> cdf; lo = 1; hi = emap -> numSamples; while (hi > lo) { i = (lo + hi) >> 1; sample = emap -> samples + i - 1; if (sample -> cdf >= du) hi = i; if (sample -> cdf < du) lo = i + 1; } /* Finalise found sample */ i = (lo + hi) >> 1; sample = emap -> samples + i - 1; /* This is a uniform mapping, mon */ cosTheta = ( 1 - (sample -> theta + pmapRandom(emitState)) * emap -> dCosTheta ); sinTheta = sqrt(1 - sqr(cosTheta)); phi = (sample -> phi + pmapRandom(emitState)) * emap -> dPhi; copycolor(ray -> rcol, sample -> pdf); } /* Normalize photon flux so that average over RGB is 1 */ colorNorm(ray -> rcol); VCOPY(ray -> rorg, emap -> photonOrg); du = cos(phi) * sinTheta; dv = sin(phi) * sinTheta; for (i = 0; i < 3; i++) ray -> rdir [i] = ( du * emap -> uh [i] + dv * emap -> vh [i] + cosTheta * emap -> wh [i] ); if (emap -> src -> sflags & SDISTANT) /* Distant source; reverse ray direction to point into the scene. */ for (i = 0; i < 3; i++) ray -> rdir [i] = -ray -> rdir [i]; if (emap -> port) { /* Photon emitted from port; move origin just behind port so it will be scattered */ for (i = 0; i < 3; i++) ray -> rorg [i] -= 2 * FTINY * ray -> rdir [i]; /* !!! PHOTON PORT REJECTION SAMPLING HACK: set photon port as fake !!! hit object for primary ray to check for its intersection in !!! tracePhoton() */ ray -> ro = emap -> port -> so; } /* Assign emitting light source index */ ray -> rsrc = emap -> src - source; /* Signal photon's acceptance */ return 1; } /* SOURCE CONTRIBS FROM DIRECT / VOLUME PHOTONS -------------------------- */ void multDirectPmap (RAY *r) /* Factor irradiance from direct photons into r -> rcol; interface to * direct() */ { COLOR photonIrrad; #ifdef PMAP_PHOTONFLOW if (lightFlowPhotonMapping) /* Photon flow mode: evaluate volume photon map as lightfield * (including direct component), ignoring participating medium. */ (lightFlowPmap -> lookup)(lightFlowPmap, r, photonIrrad); else #endif if (transientPhotonMapping) /* Lookup transient photons (direct & indirect) */ (transientPmap -> lookup)(transientPmap, r, photonIrrad); else /* Lookup direct photon irradiance */ (directPmap -> lookup)(directPmap, r, photonIrrad); /* Multiply with coefficient in ray */ multcolor(r -> rcol, photonIrrad); return; } void inscatterVolumePmap (RAY *r, COLOR inscatter) /* Add inscattering from volume photon map; interface to srcscatter() */ { /* Add ambient in-scattering via lookup callback */ -#ifdef PMAP_PHOTONFLOW - if (lightFlowPhotonMapping) { - /* Get inscatter from (transient) lightflow photons; force evaluation - as mean spherical irradiance */ - sphericalIrrad = 1; - (lightFlowPmap -> lookup)(lightFlowPmap, r, inscatter); - } - else -#endif - (volumePmap -> lookup)(volumePmap, r, inscatter); + (volumePmap -> lookup)(volumePmap, r, inscatter); }