diff --git a/pmapcontrib.c b/pmapcontrib.c index ac2496d..025019a 100644 --- a/pmapcontrib.c +++ b/pmapcontrib.c @@ -1,615 +1,626 @@ #ifndef lint static const char RCSid[] = "$Id: pmapcontrib.c,v 2.20 2021/02/16 20:06:06 greg Exp $"; #endif /* ========================================================================= Photon map routines for precomputed light source contributions; these routines interface to 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") ========================================================================= $Id: pmapcontrib.c,v 2.20 2021/02/16 20:06:06 greg Exp $ */ #include "pmapcontrib.h" #include "pmapdiag.h" #include "pmaprand.h" #include "pmapmat.h" #include "pmapsrc.h" #include "pmapio.h" #include "pmutil.h" #include "otspecial.h" #include "otypes.h" extern void SDdisk2square(double sq[2], double diskx, double disky); static int logDim (unsigned size) /* Return log2(sqrt(size)) = l, where size = (2^l)(2^l). Is there a faster way (e.g. binary search)? */ { unsigned i, sz, dim; if (!size) return 0; for (i = 0, dim = 0, sz = size; sz >>= 2; dim++); return dim; } static int xy2bin (unsigned l, int x, int y) /* Serialise 2D coordinates in range (2^l) x (2^l) to 1D bin. Returns -1 if coordinates invalid */ { return x < 0 || y < 0 ? -1 : (x << l) + y; } static void bin2xy (unsigned l, int bin, int *x, int *y) /* Deserialise 1D bin to 2D coordinates in range (2^l) x (2^l). Returns -1 if bin invalid */ { /* Obviously this is faster than integer division/modulo */ *x = bin < 0 ? -1 : bin >> l; *y = bin < 0 ? -1 : bin & (1 << l) - 1; } static int ray2bin (const RAY *ray, unsigned nbins) /* Map ray dir (pointing away from origin) to its 1D Shirley-Chiu bin, where nbins = (2^l) x (2^l) for l > 1. Returns -1 if mapped bin is invalid (e.g. behind plane) */ { const unsigned l = logDim(nbins), scDim = 1 << l; static int scRHS, varInit = 0; static FVECT scNorm, scUp; unsigned scBin [2]; FVECT diskPlane; RREAL dz, diskx, disky, rad, diskd2, scCoord [2]; if (!varInit) { /* Lazy init shirley-Chiu mapping orientation from function variables */ scRHS = varvalue(PMAP_CONTRIB_SCRHS); scNorm [0] = varvalue(PMAP_CONTRIB_SCNORMX); scNorm [1] = varvalue(PMAP_CONTRIB_SCNORMY); scNorm [2] = varvalue(PMAP_CONTRIB_SCNORMZ); scUp [0] = varvalue(PMAP_CONTRIB_SCUPX); scUp [1] = varvalue(PMAP_CONTRIB_SCUPY); scUp [2] = varvalue(PMAP_CONTRIB_SCUPZ); varInit ^= 1; } /* Map incident direction to disk */ dz = DOT(ray -> rdir, scNorm); /* normal proj */ if (dz > 0) { fcross(diskPlane, scUp, scNorm); /* y-axis in plane, perp to up */ diskPlane [0] *= scRHS; diskx = DOT(ray -> rdir, diskPlane); disky = DOT(ray -> rdir, scUp) - dz * DOT(ray -> rdir, scUp); diskd2 = diskx * diskx + disky * disky; /* in-plane length^2 */ rad = dz>FTINY ? sqrt((1 - dz*dz) / diskd2) : 0; /* radial factor */ disky *= rad; /* diskx, disky now in range [-1, 1] */ /* Apply Shirley-Chiu mapping of (diskx, disky) to square */ SDdisk2square(scCoord, diskx, disky); /* Map Shirley-Chiu square coords to 1D bin */ scBin [0] = scCoord [0] * scDim; scBin [1] = scCoord [1] * scDim; return xy2bin(l, scBin [0], scBin [1]); } else return -1; } #ifndef PMAP_CONTRIB_TEST MODCONT *addContribModifier (LUTAB *contribTab, unsigned *numContribs, char *mod, char *binParm, char *binVal, int binCnt ) /* Add light source modifier mod to contribution lookup table contribTab, and update numContribs. Return initialised contribution data for this modifier. This code adapted from rcontrib.c:addmodifier(). */ { LUENT *lutEntry = lu_find(contribTab, mod); MODCONT *contrib; EPNODE *eBinVal; if (lutEntry -> data) { /* Reject duplicate modifiers */ sprintf(errmsg, "duplicate light source modifier %s", mod); error(USER, errmsg); } if (*numContribs >= MAXMODLIST) { sprintf(errmsg, "too many modifiers (max. %d)", MAXMODLIST); error(INTERNAL, errmsg); } lutEntry -> key = mod; if (!binVal) /* Fall back to single bin if no value specified */ binVal = "0"; eBinVal = eparse(binVal); if (eBinVal -> type == NUM) { /* Bin value is constant */ binCnt = (int)(evalue(eBinVal) + 1.5); if (binCnt != 1) { sprintf(errmsg, "invalid non-zero constant for bin %s", binVal); error(USER, errmsg); } } else if (binCnt <= 0) { sprintf(errmsg, "undefined/invalid bin count for modifier %s", mod); error(USER, errmsg); } /* Allocate and init contributions */ contrib = (MODCONT *)malloc( sizeof(MODCONT) + sizeof(DCOLOR) * (binCnt - 1) ); if (!contrib) error(SYSTEM, "out of memory in addContribModifier()"); contrib -> outspec = NULL; contrib -> modname = mod; contrib -> params = binParm; contrib -> nbins = binCnt; contrib -> binv = eBinVal; contrib -> bin0 = 0; memset(contrib -> cbin, 0, sizeof(DCOLOR) * binCnt); lutEntry -> data = lutEntry -> data = (char *)contrib; ++(*numContribs); return(contrib); } void addContribModfile (LUTAB *contribTab, unsigned *numContribs, char *modFile, char *binParm, char *binVal, int binCnt ) /* Add light source modifiers from file modFile to contribution lookup table * contribTab, and update numContribs. * NOTE: This code is adapted from rcontrib.c */ { char *mods [MAXMODLIST]; int i; /* Find file and store strings */ i = wordfile(mods, MAXMODLIST, getpath(modFile, getrlibpath(), R_OK)); if (i < 0) { sprintf(errmsg, "can't open modifier file %s", modFile); error(SYSTEM, errmsg); } if (*numContribs + i >= MAXMODLIST - 1) { sprintf(errmsg, "too many modifiers (max. %d) in file %s", MAXMODLIST - 1, modFile ); error(INTERNAL, errmsg); } for (i = 0; mods [i]; i++) /* Add each modifier */ addContribModifier(contribTab, numContribs, mods [i], binParm, binVal, binCnt ); } static int contribSourceBin (LUTAB *contribs, const RAY *ray) /* Map contribution source ray to its bin for light source ray -> rsrc, using Shirley-Chiu disk-to-square mapping. Return invalid bin -1 if mapping failed. */ { const SRCREC *src; const OBJREC *srcMod; const MODCONT *srcCont; RAY srcRay; int bin, i; /* Check we have a valid ray and contribution LUT */ if (!ray || !contribs) return -1; src = &source [ray -> rsrc]; srcMod = findmaterial(src -> so); srcCont = (MODCONT*)lu_find(contribs, srcMod -> oname) -> data; if (!srcCont) /* We're not interested in this source (modifier not in contrib LUT) */ return -1; /* Set up shadow ray pointing to source for disk2square mapping */ rayorigin(&srcRay, SHADOW, NULL, NULL); srcRay.rsrc = ray -> rsrc; VCOPY(srcRay.rorg, ray -> rop); for (i = 0; i < 3; i++) srcRay.rdir [i] = -ray -> rdir [i]; if (!(src -> sflags & SDISTANT ? sourcehit(&srcRay) : (*ofun[srcMod -> otype].funp)(srcMod, &srcRay) )) /* (Redundant?) sanity check for valid source ray? */ return -1; worldfunc(RCCONTEXT, &srcRay); set_eparams((char*)srcCont -> params); if ((bin = ray2bin(&srcRay, srcCont -> nbins)) < 0) error(WARNING, "Ignoring invalid bin in contribSourceBin()"); return bin; } PhotonContribSourceIdx newPhotonContribSource (PhotonMap *pmap, const RAY *contribSrcRay, FILE *contribSrcHeap ) /* Add contribution source ray for emitted contribution photon and save * light source index and binned direction. The current contribution source * is stored in pmap -> lastContribSrc. If the previous contribution source * spawned photons (i.e. has srcIdx >= 0), it's appended to contribSrcHeap. * If contribSrcRay == NULL, the current contribution source is still * flushed, but no new source is set. Returns updated contribution source * counter pmap -> numContribSrc. */ { if (!pmap || !contribSrcHeap) return 0; /* Check if last contribution source has spawned photons (srcIdx >= 0, * see newPhoton()), in which case we save it to the heap file before * clobbering it. (Note this is short-term storage, so we doan' need * da portable I/O stuff here). */ if (pmap -> lastContribSrc.srcIdx >= 0) { - if (!fwrite(&pmap -> lastContribSrc, - sizeof(PhotonContribSource), 1, contribSrcHeap + if (!fwrite(&pmap -> lastContribSrc, sizeof(PhotonContribSource), + 1, contribSrcHeap )) - error(SYSTEM, - "failed writing photon contrib source in newPhotonContribSource()" + error(SYSTEM, "failed writing photon contrib source in " + "newPhotonContribSource()" ); pmap -> numContribSrc++; if (!pmap -> numContribSrc || pmap -> numContribSrc > PMAP_MAXCONTRIBSRC ); } /* Mark this contribution source unused with a negative source index until its path spawns a photon (see newPhoton()) */ pmap -> lastContribSrc.srcIdx = -1; /* Map ray to bin in anticipation that this contribution source will be used, since the ray will be lost once a photon is spawned */ pmap -> lastContribSrc.srcBin = contribSourceBin( pmap -> contribTab, contribSrcRay ); + + if (pmap -> lastContribSrc.srcBin < 0) { + /* Warn if invalid bin, but trace photon nonetheless. It will + count as emitted to prevent bias, but will not be stored in + newPhoton(), as it contributes zero flux */ + sprintf(errmsg, "invalid bin for light source %s, " + "contribution photons will be discarded", + source [contribSrcRay -> rsrc].so -> oname + ); + error(WARNING, errmsg); + } return pmap -> numContribSrc; } PhotonContribSourceIdx buildContribSources (PhotonMap *pmap, FILE **contribSrcHeap, char **contribSrcHeapFname, PhotonContribSourceIdx *contribSrcOfs, unsigned numHeaps ) /* Consolidate per-subprocess contribution sources heaps into array * pmap -> contribSrc. Returns offset for contribution source index * linearisation in pmap -> numContribSrc. The heap files in contribSrcHeap * are closed on return. */ { PhotonContribSourceIdx heapLen; unsigned heap; if (!pmap || !contribSrcHeap || !contribSrcOfs || !numHeaps) return 0; pmap -> numContribSrc = 0; for (heap = 0; heap < numHeaps; heap++) { contribSrcOfs [heap] = pmap -> numContribSrc; if (fseek(contribSrcHeap [heap], 0, SEEK_END) < 0) error(SYSTEM, - "failed photon contribution source seek in buildContribSources()" + "failed photon contrib source seek in buildContribSources()" ); - pmap -> numContribSrc += heapLen = ftell(contribSrcHeap [heap]) / - sizeof(PhotonContribSource); + pmap -> numContribSrc += heapLen = + ftell(contribSrcHeap [heap]) / sizeof(PhotonContribSource); if (!(pmap -> contribSrc = realloc(pmap -> contribSrc, pmap -> numContribSrc * sizeof(PhotonContribSource) ))) error(SYSTEM, "failed photon contrib source alloc in buildContribSources()" ); rewind(contribSrcHeap [heap]); if (fread(pmap -> contribSrc + contribSrcOfs [heap], sizeof(PhotonContribSource), heapLen, contribSrcHeap [heap] ) != heapLen) error(SYSTEM, "failed reading photon contrib source in buildContribSources()" ); fclose(contribSrcHeap [heap]); unlink(contribSrcHeapFname [heap]); } return pmap -> numContribSrc; - } + } static MODCONT *photonContrib (PhotonMap *pmap, RAY *ray, COLOR irrad) /* Locate photons near ray -> rop which originated from the light source modifier indexed by ray -> rsrc, and accumulate their contributions in pmap -> srcContrib. Return total contributions in irrad and pointer to binned contributions (or NULL if none were found). */ { const OBJREC *srcMod = findmaterial(source [ray -> rsrc].so); MODCONT *contrib = (MODCONT*)lu_find( pmap -> contribTab, srcMod -> oname ) -> data; Photon *photon; COLOR flux; PhotonSearchQueueNode *sqn; float r2, invArea; unsigned i; /* Zero bins for this source modifier */ memset(contrib -> cbin, 0, sizeof(DCOLOR) * contrib -> nbins); setcolor(irrad, 0, 0, 0); if (!contrib || !pmap -> maxGather) /* Modifier not in LUT or zero bandwidth */ return NULL; /* Lookup photons */ pmap -> squeue.tail = 0; /* Pass light source index to filter in findPhotons() via pmap -> lastContribSrc */ pmap -> lastContribSrc.srcIdx = ray -> rsrc; findPhotons(pmap, ray); /* Need at least 2 photons */ if (pmap -> squeue.tail < 2) { #ifdef PMAP_NONEFOUND sprintf(errmsg, "no photons found on %s at (%.3f, %.3f, %.3f)", ray -> ro ? ray -> ro -> oname : "", ray -> rop [0], ray -> rop [1], ray -> rop [2] ); error(WARNING, errmsg); #endif return NULL; } /* Average radius^2 between furthest two photons to improve accuracy and get inverse search area 1 / (PI * r^2). */ /* XXX: OMIT extra normalisation factor 1/PI for ambient calculation? */ sqn = pmap -> squeue.node + 1; r2 = max(sqn -> dist2, (sqn + 1) -> dist2); r2 = 0.25 * (pmap -> maxDist2 + r2 + 2 * sqrt(pmap -> maxDist2 * r2)); invArea = 1 / (PI * PI * r2); /* Skip the extra photon */ for (i = 1 ; i < pmap -> squeue.tail; i++, sqn++) { /* Get photon's contribution to density estimate */ photon = getNearestPhoton(&pmap -> squeue, sqn -> idx); getPhotonFlux(photon, flux); scalecolor(flux, invArea); #ifdef PMAP_EPANECHNIKOV /* Apply Epanechnikov kernel to photon flux based on photon distance */ scalecolor(flux, 2 * (1 - sqn -> dist2 / r2)); #endif addcolor(irrad, flux); addcolor(contrib -> cbin [photonSrcBin(pmap, photon)], flux); } return contrib; } void preComputeContrib (PhotonMap *pmap) /* Precompute contributions for a random subset of finalGather * pmap -> numPhotons photons, and build the precomputed contribution photon map, discarding the original photons. */ /* !!! NOTE: PRECOMPUTATION WITH OOC CURRENTLY WITHOUT CACHE !!! */ { unsigned long i, numPreComp; unsigned j; PhotonIdx pIdx; Photon photon; RAY ray; MODCONT *binnedContribs; PhotonMap nuPmap; repComplete = numPreComp = finalGather * pmap -> numPhotons; if (verbose) { sprintf(errmsg, "\nPrecomputing contributions for %ld photons\n", numPreComp ); eputs(errmsg); #if NIX fflush(stderr); #endif } /* Copy photon map for precomputed photons */ memcpy(&nuPmap, pmap, sizeof(PhotonMap)); /* Zero counters, init new heap and extents */ nuPmap.numPhotons = 0; initPhotonHeap(&nuPmap); for (j = 0; j < 3; j++) { nuPmap.minPos [j] = FHUGE; nuPmap.maxPos [j] = -FHUGE; } /* Record start time, baby */ repStartTime = time(NULL); #ifdef SIGCONT signal(SIGCONT, pmapPreCompReport); #endif repProgress = 0; photonRay(NULL, &ray, PRIMARY, NULL); ray.ro = NULL; for (i = 0; i < numPreComp; i++) { do { /* Get random photon from stratified distribution in source heap * to avoid duplicates and clustering */ pIdx = firstPhoton(pmap) + (unsigned long)( (i + pmapRandom(pmap -> randState)) / finalGather ); getPhoton(pmap, pIdx, &photon); /* Init dummy photon ray with intersection and normal at photon position. Set emitting light source index from origin. */ VCOPY(ray.rop, photon.pos); ray.rsrc = photonSrcIdx(pmap, &photon); for (j = 0; j < 3; j++) ray.ron [j] = photon.norm [j] / 127.0; /* Get contributions at photon position; retry with another photon * if no contribs found */ binnedContribs = photonContrib(pmap, &ray, ray.rcol); } while (!binnedContribs); /* Append photon to new heap from ray */ newPhoton(&nuPmap, &ray); /* TODO: Apply wavelet xform to binned contribs, apply lossy compression, encode to RGBE */ #if 0 encodeContrib(binnedContribs, pmap -> compRatio); #endif /* Update progress */ repProgress++; if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime) pmapPreCompReport(); #ifdef SIGCONT else signal(SIGCONT, pmapPreCompReport); #endif } /* Flush heap */ flushPhotonHeap(&nuPmap); #ifdef SIGCONT signal(SIGCONT, SIG_DFL); #endif /* Trash original pmap, replace with precomputed one */ deletePhotons(pmap); memcpy(pmap, &nuPmap, sizeof(PhotonMap)); if (verbose) { eputs("\nRebuilding precomputed contribution photon map\n"); #if NIX fflush(stderr); #endif } /* Rebuild underlying data structure, destroying heap */ buildPhotonMap(pmap, NULL, NULL, 1); /* TODO: Save wavelet coeffs here? */ } #else #include #include #include "func.h" int main (int argc, char *argv []) { unsigned i, l, nbins, nsamp, numTheta = 0, numPhi = 0; double t, p; RAY ray; int bin; if (argc < 3) { fprintf(stderr, "%s [= ..]\n", argv [0]); fprintf(stderr, "Default variable values: %s\n", PMAP_CONTRIB_SCDEFAULTS ); fputs("\nMissing resolution l>1, number of samples nsamp\n", stderr ); return -1; } if (!(l = atoi(argv [1]))) { fputs("Invalid resolution\n", stderr); return -1; } else nbins = 1 << (l << 1); if (!(nsamp = atoi(argv [2]))) { fputs("Invalid num samples\n", stderr); return -1; } else { numTheta = (int)(sqrt(nsamp) / 2); numPhi = 4* numTheta; } /* (Doan' need to) Init cal func routines for binning */ #if 0 initfunc(); setcontext(RCCONTEXT); #endif /* Compile default orientation variables for contribution binning */ scompile(PMAP_CONTRIB_SCDEFAULTS, NULL, 0); /* Compile custom orientation variabls from command line */ for (i = 3; i < argc; i++) scompile(argv [i], NULL, 0); for (i = 0; i < nsamp; i++) { #if 0 /* Random */ t = 0.5 * PI * drand48(); p = 2 * PI * drand48(); #else /* Stratified */ t = 0.5 * PI * ((i % numTheta) + drand48()) / numTheta; p = 2.0 * PI * ((i / numTheta) + drand48()) / numPhi; #endif ray.rdir [0] = sin(t) * cos(p); ray.rdir [1] = sin(t) * sin(p); ray.rdir [2] = cos(t); bin = ray2bin(&ray, nbins); printf("%.3f\t%.3f\t%.3f\t-->\t%d\n", ray.rdir [0], ray.rdir [1], ray.rdir [2], bin ); } return 0; } #endif diff --git a/pmapdata.c b/pmapdata.c index 40e14f8..30f8488 100644 --- a/pmapdata.c +++ b/pmapdata.c @@ -1,826 +1,834 @@ #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 "pmaproi.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; /* 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 */ 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; +#if 0 /* NOTE: Invalid bins are now handled in distribPhotonContrib() + immediately after returning from this routine */ + if (photonSrcBin(pmap, &photon) < 0) + /* Photon has an invalid bin, and therefore implicitly zero + contributions; don't bother storing! */ + return -1; +#endif + /* Photon will be stored; set contribution source index, + thereby marking it as having spawned photon(s) */ 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. */ 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 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]; 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, 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/pmcontrib2.c b/pmcontrib2.c index a393d23..4c7f0e6 100644 --- a/pmcontrib2.c +++ b/pmcontrib2.c @@ -1,722 +1,732 @@ #ifndef lint static const char RCSid[] = "$Id$"; #endif /* ========================================================================= Photon map routines for precomputed light source contributions; this module contains the main photon distribution routine as part of 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") ========================================================================= $Id$ */ #include "pmapcontrib.h" #include "pmapdiag.h" #include "pmaprand.h" #include "pmapmat.h" #include "pmaproi.h" #include "pmapsrc.h" #include "otspecial.h" #if NIX #include #include #endif /* Defs for photon emission counter array passed by sub-processes to parent * via shared memory */ typedef unsigned long PhotonContribCnt; /* Indices for photon emission counter array: num photons stored and num * emitted per source */ #define PHOTONCNT_NUMPHOT 0 #define PHOTONCNT_NUMEMIT(n) (1 + n) /* Photon distribution counter update interval expressed as bitmask; counters shared among subling subprocesses will only be updated in multiples of PMAP_CNTUPDATE in order to reduce contention */ #define PMAP_CNTUPDATE 0xffL void distribPhotonContrib (PhotonMap* pmap, LUTAB *contribTab, unsigned numContribs, unsigned numProc ) { EmissionMap emap; char errmsg2 [128], shmFname [PMAP_TMPFNLEN]; unsigned srcIdx, proc; int shmFile, stat, pid; double *srcFlux, /* Emitted flux per light source */ srcDistribTarget; /* Target photon count per source */ PhotonContribCnt *photonCnt, /* Photon emission counter array */ lastPhotonCnt [PHOTONCNT_NUMEMIT(nsources)]; unsigned photonCntSize = (sizeof(PhotonContribCnt) * PHOTONCNT_NUMEMIT(nsources) ); FILE **contribSrcHeap = NULL; char **contribSrcHeapFname = NULL; PhotonContribSourceIdx *contribSrcOfs = NULL; pid_t procPids [PMAP_MAXPROC]; if (!pmap) error(USER, "no contribution photon map specified"); if (!nsources) error(USER, "no light sources"); if (nsources > MAXMODLIST) error(USER, "too many light sources"); if (!contribTab || !numContribs) error(USER, "no modifiers specified for contribution photon map"); /* Allocate photon flux per light source; this differs for every * source as all sources contribute the same number of distributed * photons (srcDistribTarget), hence the number of photons emitted per * source does not correlate with its emitted flux. The resulting flux * per photon is therefore adjusted individually for each source. */ if (!(srcFlux = calloc(nsources, sizeof(double)))) error(SYSTEM, "can't allocate source flux in distribPhotonContrib"); /* =================================================================== * INITIALISATION - Set up emission and scattering funcs * =================================================================== */ emap.samples = NULL; emap.src = NULL; emap.maxPartitions = MAXSPART; emap.partitions = (unsigned char*)malloc(emap.maxPartitions >> 1); if (!emap.partitions) error(USER, "can't allocate source partitions in distribPhotonContrib"); /* Initialise contrib photon map */ initPhotonMap(pmap, PMAP_TYPE_CONTRIB); initPmapContrib(contribTab); initPhotonHeap(pmap); initPhotonEmissionFuncs(); initPhotonScatterFuncs(); /* Per-subprocess / per-source target counts */ pmap -> distribTarget /= numProc; srcDistribTarget = nsources ? (double)pmap -> distribTarget / nsources : 0; if (!pmap -> distribTarget) error(INTERNAL, "no photons to distribute in distribPhotonContrib"); /* Get photon ports from modifier list */ getPhotonPorts(photonPortList); /* Get photon sensor modifiers */ getPhotonSensors(photonSensorList); /* Get polyhedral regions of interest */ getPolyROIs(pmapROImodList); #if NIX /* Set up shared mem for photon counters (zeroed by ftruncate) */ strcpy(shmFname, PMAP_TMPFNAME); shmFile = mkstemp(shmFname); if (shmFile < 0 || ftruncate(shmFile, photonCntSize) < 0) error(SYSTEM, "failed shared mem init in distribPhotonContrib"); photonCnt = mmap(NULL, photonCntSize, PROT_READ | PROT_WRITE, MAP_SHARED, shmFile, 0 ); if (photonCnt == MAP_FAILED) error(SYSTEM, "failed shsared mem mapping in distribPhotonContrib"); #else /* Allocate photon counters statically on Windoze */ if (!(photonCnt = malloc(photonCntSize))) error(SYSTEM, "failed trivial malloc in distribPhotonContrib"); for (srcIdx = 0; srcIdx < PHOTONCNT_NUMEMIT(nsources); srcIdx++) photonCnt [srcIdx] = 0; #endif /* NIX */ /* Zero overflow tracking counters */ memset(&lastPhotonCnt, 0, sizeof(lastPhotonCnt)); if (verbose) { sprintf(errmsg, "\nIntegrating flux from %d sources", nsources); if (photonPorts) { sprintf(errmsg2, " via %d ports", numPhotonPorts); strcat(errmsg, errmsg2); } strcat(errmsg, "\n"); eputs(errmsg); } /* ============================================================= * FLUX INTEGRATION - Get total flux emitted from sources/ports * ============================================================= */ for (srcIdx = 0; srcIdx < nsources; srcIdx++) { unsigned portCnt = 0; const OBJREC *srcMod = findmaterial(source [srcIdx].so); srcFlux [srcIdx] = 0; /* Skip this source if its contributions are not sought */ if (!lu_find(pmap -> contribTab, srcMod -> oname) -> data) { sprintf(errmsg, "ignoring contributions from source %s", source [srcIdx].so -> oname ); error(WARNING, errmsg); continue; } emap.src = source + srcIdx; do { /* Need at least one iteration if no ports! */ emap.port = emap.src -> sflags & SDISTANT ? photonPorts + portCnt : NULL; photonPartition [emap.src -> so -> otype] (&emap); if (verbose) { sprintf(errmsg, "\tIntegrating flux from source %s ", source [srcIdx].so -> oname ); if (emap.port) { sprintf(errmsg2, "via port %s ", photonPorts [portCnt].so -> oname ); strcat(errmsg, errmsg2); } sprintf(errmsg2, "(%lu partitions)\n", emap.numPartitions); strcat(errmsg, errmsg2); eputs(errmsg); #if NIX fflush(stderr); #endif } for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions; emap.partitionCnt++ ) { initPhotonEmission(&emap, pdfSamples); srcFlux [srcIdx] += colorAvg(emap.partFlux); } portCnt++; } while (portCnt < numPhotonPorts); if (srcFlux [srcIdx] < FTINY) { sprintf(errmsg, "source %s has zero emission", source [srcIdx].so -> oname ); error(WARNING, errmsg); } } /* Allocate & init per-subprocess contribution source heap files */ contribSrcHeap = calloc(numProc, sizeof(FILE*)); contribSrcHeapFname = calloc(numProc, sizeof(char*)); contribSrcOfs = calloc(numProc, sizeof(PhotonContribSourceIdx)); if (!contribSrcHeap || !contribSrcHeapFname || !contribSrcOfs) error(SYSTEM, "failed contribution source heap allocation " "in distribPhotonContrib()" ); for (proc = 0; proc < numProc; proc++) { contribSrcHeapFname [proc] = malloc(PMAP_TMPFNLEN); if (!contribSrcHeapFname [proc]) error(SYSTEM, "failed contribution source heap file allocation " "in distribPhotonContrib()" ); mktemp(strcpy(contribSrcHeapFname [proc], PMAP_TMPFNAME)); if (!(contribSrcHeap [proc] = fopen(contribSrcHeapFname [proc], "w+b"))) error(SYSTEM, "failed opening contribution source heap file " "in distribPhotonContrib()" ); } /* Record start time for progress reports */ repStartTime = time(NULL); if (verbose) { sprintf(errmsg, "\nPhoton distribution @ %d procs\n", numProc); eputs(errmsg); } /* MAIN LOOP */ for (proc = 0; proc < numProc; proc++) { #if NIX if (!(pid = fork())) { /* SUBPROCESS ENTERS HERE; opened and mmapped files inherited */ #else if (1) { /* No subprocess under Windoze */ #endif /* Local photon counters for this subprocess */ unsigned long lastNumPhotons = 0, localNumEmitted = 0; double photonFluxSum = 0; /* Accum. photon flux */ /* Seed RNGs from PID for decorellated photon distribution */ pmapSeed(randSeed + proc, partState); pmapSeed(randSeed + (proc + 1) % numProc, emitState); pmapSeed(randSeed + (proc + 2) % numProc, cntState); pmapSeed(randSeed + (proc + 3) % numProc, mediumState); pmapSeed(randSeed + (proc + 4) % numProc, scatterState); pmapSeed(randSeed + (proc + 5) % numProc, rouletteState); #ifdef PMAP_SIGUSR { double partNumEmit; unsigned long partEmitCnt; double srcPhotonFlux, avgPhotonFlux; unsigned portCnt, passCnt, prePassCnt; float srcPreDistrib; double srcNumEmit; /* # to emit from source */ unsigned long srcNumDistrib; /* # stored */ void sigUsrDiags() /* Loop diags via SIGUSR1 */ { sprintf(errmsg, "********************* Proc %d Diags *********************\n" "srcIdx = %d (%s)\nportCnt = %d (%s)\npassCnt = %d\n" "srcFlux = %f\nsrcPhotonFlux = %f\navgPhotonFlux = %f\n" "partNumEmit = %f\npartEmitCnt = %lu\n\n", proc, srcIdx, findmaterial(source [srcIdx].so) -> oname, portCnt, photonPorts [portCnt].so -> oname, passCnt, srcFlux [srcIdx], srcPhotonFlux, avgPhotonFlux, partNumEmit, partEmitCnt ); eputs(errmsg); fflush(stderr); } } #endif #ifdef PMAP_SIGUSR signal(SIGUSR1, sigUsrDiags); #endif #ifdef DEBUG_PMAP /* Output child process PID after random delay to prevent corrupted * console output due to race condition */ usleep(1e6 * pmapRandom(rouletteState)); fprintf(stderr, "Proc %d: PID = %d (waiting 10 sec to attach debugger...)\n", proc, getpid() ); /* Allow time for debugger to attach to child process */ sleep(10); #endif /* ============================================================= * 2-PASS PHOTON DISTRIBUTION * Pass 1 (pre): emit fraction of target photon count * Pass 2 (main): based on outcome of pass 1, estimate remaining * number of photons to emit to approximate target * count * ============================================================= */ for (srcIdx = 0; srcIdx < nsources; srcIdx++) { const unsigned numEmitIdx = PHOTONCNT_NUMEMIT(srcIdx); #ifndef PMAP_SIGUSR unsigned portCnt, passCnt = 0, prePassCnt = 0; float srcPreDistrib = preDistrib; double srcNumEmit = 0; /* # to emit from source */ unsigned long srcNumDistrib = pmap -> numPhotons; /* #stored */ #else passCnt = prePassCnt = 0; srcPreDistrib = preDistrib; srcNumEmit = 0; /* # to emit from source */ srcNumDistrib = pmap -> numPhotons; /* # stored */ #endif if (srcFlux [srcIdx] < FTINY) /* Source has zero emission or was skipped in prepass because its contributions are not sought */ continue; while (passCnt < 2) { if (!passCnt) { /* INIT PASS 1 */ if (++prePassCnt > maxPreDistrib) { /* Warn if no photons contributed after sufficient * iterations; only output from subprocess 0 to reduce * console clutter */ if (!proc) { sprintf(errmsg, "source %s: too many prepasses, skipped", source [srcIdx].so -> oname ); error(WARNING, errmsg); } break; } /* Num to emit is fraction of target count */ srcNumEmit = srcPreDistrib * srcDistribTarget; } else { /* INIT PASS 2 */ #ifndef PMAP_SIGUSR double srcPhotonFlux, avgPhotonFlux; #endif /* Based on the outcome of the predistribution we can now * figure out how many more photons we have to emit from * the current source to meet the target count, * srcDistribTarget. This value is clamped to 0 in case * the target has already been exceeded in pass 1. * srcNumEmit and srcNumDistrib is the number of photons * emitted and distributed (stored) from the current * source in pass 1, respectively. */ srcNumDistrib = pmap -> numPhotons - srcNumDistrib; srcNumEmit *= srcNumDistrib ? max(srcDistribTarget/srcNumDistrib, 1) - 1 : 0; if (!srcNumEmit) /* No photons left to distribute in main pass */ break; srcPhotonFlux = srcFlux [srcIdx] / srcNumEmit; avgPhotonFlux = photonFluxSum / (srcIdx + 1); if (avgPhotonFlux > FTINY && srcPhotonFlux / avgPhotonFlux < FTINY ) { /* Skip source if its photon flux is grossly below the * running average, indicating negligible contributions * at the expense of excessive distribution time; only * output from subproc 0 to reduce console clutter */ if (!proc) { sprintf(errmsg, "source %s: itsy bitsy photon flux, skipped", source [srcIdx].so -> oname ); error(WARNING, errmsg); } srcNumEmit = 0; /* Or just break??? */ } /* Update sum of photon flux per light source */ photonFluxSum += srcPhotonFlux; } portCnt = 0; do { /* Need at least one iteration if no ports! */ emap.src = source + srcIdx; emap.port = emap.src -> sflags & SDISTANT ? photonPorts + portCnt : NULL; photonPartition [emap.src -> so -> otype] (&emap); if (verbose && !proc) { /* Output from subproc 0 only to avoid race condition * on console I/O */ if (!passCnt) sprintf(errmsg, "\tPREPASS %d on source %s ", prePassCnt, source [srcIdx].so -> oname ); else sprintf(errmsg, "\tMAIN PASS on source %s ", source [srcIdx].so -> oname ); if (emap.port) { sprintf(errmsg2, "via port %s ", photonPorts [portCnt].so -> oname ); strcat(errmsg, errmsg2); } sprintf(errmsg2, "(%lu partitions)\n", emap.numPartitions ); strcat(errmsg, errmsg2); eputs(errmsg); #if NIX fflush(stderr); #endif } for (emap.partitionCnt = 0; emap.partitionCnt < emap.numPartitions; emap.partitionCnt++ ) { #ifndef PMAP_SIGUSR double partNumEmit; unsigned long partEmitCnt; #endif /* Get photon origin within current source partishunn * and build emission map */ photonOrigin [emap.src -> so -> otype] (&emap); initPhotonEmission(&emap, pdfSamples); /* Number of photons to emit from ziss partishunn; * scale according to its normalised contribushunn to * the emitted source flux */ partNumEmit = (srcNumEmit * colorAvg(emap.partFlux) / srcFlux [srcIdx] ); partEmitCnt = (unsigned long)partNumEmit; /* Probabilistically account for fractional photons */ if (pmapRandom(cntState) < partNumEmit - partEmitCnt) partEmitCnt++; /* Integer counter avoids FP rounding errors during * iteration */ while (partEmitCnt--) { RAY photonRay; /* Emit photon according to PDF (if any). If * accepted, allocate associated contribution * origin, and trace through scene until * absorbed/leaked; emitPhoton() sets the emitting * light source index in photonRay */ /* NOTE: rejection sampling skips incrementing the * emission counter (see below), thus compensating * for the rejected photons by increasing the photon * flux in proportion to the lower effective * emission count. * BUG: THIS INTERFERES WITH THE PROGRESS COUNTER * REPORTED TO THE PARENT, AND WITH THE * PREDISTRIBUTION PASS --> PHOTON DISTRIBUTION WILL * FINISH EARLY, WITH FEWER PHOTONS THAN TARGETTED! */ - if (!emitPhoton(&emap, &photonRay)) continue; + if (!emitPhoton(&emap, &photonRay)) + continue; newPhotonContribSource(pmap, &photonRay, contribSrcHeap [proc] ); + + /* Skip photon if it has an invalid bin. It will + implicitly contribute zero flux, so don't bother + tracing and storing it. However, it counts as + emitted (see enclosing while() loop) to avoid + bias. */ + if (pmap -> lastContribSrc.srcBin < 0) + continue; + /* Set subprocess index in photonRay for post- * distrib contribution source index linearisation; * this is propagated with the contrib source index * in photonRay and set for photon hits by * newPhoton() */ PMAP_SETRAYPROC(&photonRay, proc); tracePhoton(&photonRay); /* Update local emission counter */ localNumEmitted++; if (!(localNumEmitted & PMAP_CNTUPDATE)) { /* Update global counters shared with siblings in intervals to reduce overhead/contention */ lastPhotonCnt [numEmitIdx] = photonCnt [numEmitIdx]; photonCnt [numEmitIdx] += PMAP_CNTUPDATE; /* Check for overflow using previous values */ if (photonCnt [numEmitIdx] < lastPhotonCnt [numEmitIdx] ) { sprintf(errmsg, "photon emission counter " "overflow (was: %ld, is: %ld)", lastPhotonCnt [numEmitIdx], photonCnt [numEmitIdx] ); error(INTERNAL, errmsg); } /* Differentially increment photon counter */ lastNumPhotons = pmap -> numPhotons; photonCnt [PHOTONCNT_NUMPHOT] += pmap -> numPhotons - lastNumPhotons; /* Check for photon counter overflow (this could only happen before an emission counter overflow if the scene has an absurdly high albedo and/or very dense geometry) */ if (photonCnt [PHOTONCNT_NUMPHOT] < lastNumPhotons ) { sprintf(errmsg, "photon counter overlow " "(was: %ld, is: %ld)", lastNumPhotons, photonCnt [PHOTONCNT_NUMPHOT] ); error(INTERNAL, errmsg); } } } #if !NIX /* Synchronous progress report on Windoze */ if (!proc && photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime ) { unsigned s; repComplete = pmap -> distribTarget * numProc; repProgress = photonCnt [PHOTONCNT_NUMPHOT]; for (repEmitted = 0, s = 0; s < nsources; s++) repEmitted += photonCnt [numEmitIdx]; pmapDistribReport(); } #endif } portCnt++; } while (portCnt < numPhotonPorts); if (pmap -> numPhotons == srcNumDistrib) { /* Double predistrib factor in case no photons were stored * for this source and redo pass 1 */ srcPreDistrib *= 2; } else { /* Now do pass 2 */ passCnt++; } } } /* Flush heap buffa one final time to prevent data corruption */ flushPhotonHeap(pmap); /* Flush last contribution origin to origin heap file */ newPhotonContribSource(pmap, NULL, contribSrcHeap [proc]); /* Heap files closed automatically on exit fclose(pmap -> heap); fclose(orgHeap [proc]); */ #ifdef DEBUG_PMAP sprintf( errmsg, "Proc %d total %ld photons\n", proc, pmap -> numPhotons ); eputs(errmsg); fflush(stderr); #endif #ifdef PMAP_SIGUSR signal(SIGUSR1, SIG_DFL); #endif #if NIX /* Terminate subprocess */ exit(0); #endif } else { /* PARENT PROCESS CONTINUES LOOP ITERATION HERE */ if (pid < 0) error(SYSTEM, "failed to fork subprocess in distribPhotonContrib()" ); else /* Saves child process IDs */ procPids [proc] = pid; } } #if NIX /* PARENT PROCESS CONTINUES HERE */ #ifdef SIGCONT /* Enable progress report signal handler */ signal(SIGCONT, pmapDistribReport); #endif /* Wait for subprocesses to complete while reporting progress */ proc = numProc; while (proc) { while (waitpid(-1, &stat, WNOHANG) > 0) { /* Subprocess exited; check status */ if (!WIFEXITED(stat) || WEXITSTATUS(stat)) { /* Exited with error; terminate siblings, clean up & bail out */ for (proc = 0; proc < numProc; proc++) kill(procPids [proc], SIGKILL); /* Unmap shared memory, squish mapped file */ munmap(photonCnt, sizeof(*photonCnt)); close(shmFile); unlink(shmFname); error(USER, "failed photon distribution"); } --proc; } /* Nod off for a bit and update progress */ sleep(1); /* Asynchronous progress report from shared subprocess counters */ repComplete = pmap -> distribTarget * numProc; repProgress = photonCnt [PHOTONCNT_NUMPHOT]; for (repEmitted = 0, srcIdx = 0; srcIdx < nsources; srcIdx++) repEmitted += photonCnt [PHOTONCNT_NUMEMIT(srcIdx)]; /* Get global photon count from shmem updated by subprocs */ pmap -> numPhotons = photonCnt [PHOTONCNT_NUMPHOT]; if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime) pmapDistribReport(); #ifdef SIGCONT else signal(SIGCONT, pmapDistribReport); #endif } #endif /* NIX */ /* ================================================================ * POST-DISTRIBUTION - Set photon flux and build kd-tree, etc. * ================================================================ */ #ifdef SIGCONT /* Reset signal handler */ signal(SIGCONT, SIG_DFL); #endif free(emap.samples); if (!pmap -> numPhotons) error(USER, "empty contribution photon map"); /* Load per-subprocess contribution sources into pmap -> contribSrc */ /* Dumb compilers apparently need the char** cast */ pmap -> numContribSrc = buildContribSources(pmap, contribSrcHeap, (char**)contribSrcHeapFname, contribSrcOfs, numProc ); if (!pmap -> numContribSrc) error(INTERNAL, "no contribution sources in contribution photon map"); /* Set photon flux per source */ /* TODO: HOW/WHERE DO WE HANDLE COEFFICIENT MODE??? */ for (srcIdx = 0; srcIdx < nsources; srcIdx++) srcFlux [srcIdx] /= photonCnt [PHOTONCNT_NUMEMIT(srcIdx)]; #if NIX /* Photon counters no longer needed, unmap shared memory */ munmap(photonCnt, sizeof(*photonCnt)); close(shmFile); unlink(shmFname); #else free(photonCnt); #endif if (verbose) { eputs("\nBuilding contribution photon map...\n"); #if NIX fflush(stderr); #endif } /* Build underlying data structure; heap is destroyed */ buildPhotonMap(pmap, srcFlux, contribSrcOfs, numProc); /* Precompute binned photon contributions */ if (verbose) eputs("\n"); preComputeContrib(contribPmap); /* Free per-subprocess origin heap files */ for (proc = 0; proc < numProc; proc++) free(contribSrcHeapFname [proc]); free(contribSrcHeapFname); free(contribSrcHeap); free(contribSrcOfs); if (verbose) eputs("\n"); }