diff --git a/pmapcontrib.c b/pmapcontrib.c index 1923138..4b253d5 100644 --- a/pmapcontrib.c +++ b/pmapcontrib.c @@ -1,1685 +1,1688 @@ #ifndef lin 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 handle contribution binning, compression and encoding, and are used by mkpmap. Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) (c) Lucerne University of Applied Sciences and Arts, supported by the Swiss National Science Foundation (SNSF #147053, "Daylight Redirecting Components", SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment") (c) Tokyo University of Science, supported by the KAJIMA Foundation under the project title: "Reflections from Building Façades and their Glare Potential on the Built Environment -- Application of the Photon Flow Method using Annual Simulation" ========================================================================= $Id: pmapcontrib.c,v 2.20 2021/02/16 20:06:06 greg Exp $ */ #include "pmapcontrib.h" #ifdef PMAP_CONTRIB #include "pmapdiag.h" #include "pmaprand.h" #include "pmapmat.h" #include "pmapsrc.h" #include "pmutil.h" #include "otspecial.h" #include "otypes.h" #include "lookup.h" #if NIX #include <sys/mman.h> #include <sys/wait.h> #endif /* The following are convenient placeholders interfacing to mkpmap and rcontrib. They can be externally set via initPmapContribTab() and referenced within the contrib pmap modules. These variables can then be safely ignored by rtrace/rpict/rvu, without annoying linking errors. */ /* Global pointer to rcontrib's contribution binning LUT */ LUTAB *pmapContribTab = NULL; /* Contribution/coefficient mode flag */ int *pmapContribMode; #ifdef PMAP_CONTRIB_DBG /* Sum/num of mRGBE encoding errors for wavelet coeffs */ static WAVELET_COEFF mrgbeDiffs = 0; static unsigned long mrgbeDiffCnt = 0; #endif static int ray2bin (const RAY *ray, unsigned scDim) /* Map ray dir (pointing away from origin) to its 1D bin in an (scDim x scDim) Shirley-Chiu square. Returns -1 if mapped bin is invalid (e.g. behind plane) */ { 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 func variables */ /* XXX: Note inversion of normal and up vectors to match rcontrib's convention, since the input ray points away from its origin */ 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 */ /* normalize(diskPlane); NEED TO NORMALISE? */ diskPlane [0] *= scRHS; diskx = DOT(ray -> rdir, diskPlane); disky = DOT(ray -> rdir, scUp); /* Apply Shirley-Chiu mapping of (diskx, disky) to square */ disk2square(scCoord, diskx, disky); /* Map Shirley-Chiu square coords to 1D bin */ scBin [0] = scCoord [0] * scDim; scBin [1] = scCoord [1] * scDim; return PMAP_CONTRIB_XY2LIN(scDim, scBin [0], scBin [1]); } else return -1; } /* ------------------ CONTRIBSRC STUFF --------------------- */ #ifndef PMAP_CONTRIB_TEST MODCONT *addContribModifier (LUTAB *contribTab, unsigned *numContribs, char *mod, char *binParm, 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; unsigned numCoeffs; if (binCnt > 1) { /* Binning enabled; warn if potential coefficient index overflow in mRGBE encoding. This requires getting the number of wavelet coefficients generated by the transform a priori. */ if (!(numCoeffs = padWaveletXform2(NULL, NULL, sqrt(binCnt), NULL))) { sprintf(errmsg, "can't determine number of wavelet coefficients " "for modifier %s", mod ); error(INTERNAL, errmsg); } if (numCoeffs * numCoeffs > PMAP_CONTRIB_MAXCOEFFS) { sprintf(errmsg, "mRGBE data field may overflow for modifier %s; reduce -bn " "and/or compression if contribution precomputation fails", mod ); error(WARNING, errmsg); } } 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 (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, 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, 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) /* Not interested in this source (modifier not in contrib LUT) */ return -1; if (srcCont -> nbins <= 1) /* Binning disabled; return bin 0 */ return 0; /* 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; #if 0 worldfunc(RCCONTEXT, &srcRay); #endif set_eparams((char*)srcCont -> params); if ((bin = ray2bin(&srcRay, sqrt(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 contrib * 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 )) error(SYSTEM, "failed writing photon contrib source in " "newPhotonContribSource()" ); pmap -> numContribSrc++; if (!pmap -> numContribSrc || pmap -> numContribSrc > PMAP_MAXCONTRIBSRC ) error(INTERNAL, "contribution source overflow"); } if (contribSrcRay) { /* 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 contrib source will be used, since the ray will be lost once a photon is spawned */ pmap -> lastContribSrc.srcBin = contribSourceBin( pmapContribTab, 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 source 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 contrib source seek " "in buildContribSources()" ); 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; } /* ----------------- CONTRIB PMAP ENCODING STUFF -------------------- */ static void coeffSwap (PreComputedContribCoeff *c1, PreComputedContribCoeff *c2 ) { PreComputedContribCoeff tCoeff; tCoeff.coeff = c1 -> coeff; tCoeff.idx = c1 -> idx; c1 -> coeff = c2 -> coeff; c1 -> idx = c2 -> idx; c2 -> coeff = tCoeff.coeff; c2 -> idx = tCoeff.idx; } static int coeffPartition (PreComputedContribCoeff *coeffs, unsigned medianPos, unsigned left, unsigned right ) /* REVERSE partition coefficients by magnitude, such that coeffs [left..medianPos] >= coeffs [medianPos+1..right]. Returns median's position. */ { unsigned long l, r, m; WAVELET_COEFF lVal, rVal, mVal, tVal; #define COEFFVAL(c) (DOT((c) -> coeff, (c) -> coeff)) if (left < right) { /* Select pivot from median-of-three and move to photons [right] (a.k.a. Lomuto partitioning) */ l = left; r = right; m = l + ((r - l) >> 1); /* Avoids overflow vs. (l+r) >> 1 */ lVal = COEFFVAL(coeffs + l); rVal = COEFFVAL(coeffs + r); mVal = COEFFVAL(coeffs + m); if (mVal > lVal) { coeffSwap(coeffs + m, coeffs + l); tVal = mVal; mVal = lVal; lVal = tVal; } if (rVal > lVal) { coeffSwap(coeffs + r, coeffs + l); tVal = rVal; rVal = lVal; lVal = tVal; } if (mVal > rVal) { coeffSwap(coeffs + m, coeffs + r); tVal = mVal; mVal = rVal; rVal = tVal; } /* Pivot with key rVal is now in coeffs [right] */ /* l & r converge, swapping coefficients out of (reversed) order with respect to pivot. The convergence point l = r is the pivot's final position */ while (l < r) { while (l < r && COEFFVAL(coeffs + l) > rVal) l++; while (r > l && COEFFVAL(coeffs + r) <= rVal) r--; /* Coeffs out of order, must swap */ if (l < r) coeffSwap(coeffs + l, coeffs + r); }; /* Now l == r is pivot's final position; swap these coeffs */ coeffSwap(coeffs + l, coeffs + right); /* Recurse in partition containing median */ if (l > medianPos) return coeffPartition(coeffs, medianPos, left, l - 1); else if (l < medianPos) return coeffPartition(coeffs, medianPos, l + 1, right); else /* l == medianPos, partitioning done */ return l; } else return left; } static int coeffIdxCompare (const void *c1, const void *c2) /* Comparison function to sort thresholded coefficients by index */ /* TODO: scale coeffs by 2^-l to adapt to resolution/frequency l? (this drops higher-frequency coeffs first). */ { const PreComputedContribCoeff *tcoeff1 = c1, *tcoeff2 = c2; const unsigned v1 = tcoeff1 -> idx, v2 = tcoeff2 -> idx; if (v1 < v2) return -1; else if (v1 > v2) return 1; else return 0; } static int thresholdContribs (PreComputedContrib *preCompContrib) /* Threshold wavelet detail coefficients in preCompContrib -> waveletMatrix [approxDim..coeffDim-1] [approxDim..coeffDim-1] (where approxDim = WAVELET_PADD2_APPROXDIM) by keeping the (preCompContrib -> nCompressedCoeffs) largest of these and returning them in preCompContrib -> compressedCoeffs along with their original matrix indices. NOTE: The wavelet approximation coefficients preCompContrib -> waveletMatrix [0..approxDim-1] [0..approxDim-1] are excluded from thresholding to minimise compression artefacts. Returns 0 on success, else -1. */ { unsigned i, j, coeffDim, coeffIdx, nNzDetailCoeffs, nCompressedCoeffs, numThresh; WaveletMatrix2 waveletMatrix; PreComputedContribCoeff *threshCoeffs, *threshCoeffPtr; WAVELET_COEFF *coeffPtr; if (!preCompContrib || !(coeffDim = preCompContrib -> coeffDim) || !(threshCoeffs = preCompContrib -> compressedCoeffs) || !(waveletMatrix = preCompContrib -> waveletMatrix) ) /* Should be initialised by caller! */ return -1; /* Set up coefficient buffer for compression (thresholding), skipping * the approximation coefficients in the upper left of waveletMatrix, * which are not thresholded. Also skip zero coefficients (resulting * from padding), since these are already implicitly thresholded. The * original 2D matrix indices are linearised to 1D and saved with each * coefficient to restore the original sparse coefficient matrix. */ for (i = 0, threshCoeffPtr = threshCoeffs; i < coeffDim; i++) { /* Calc row pointer (=mult) in outer loop, inc in inner */ coeffIdx = PMAP_CONTRIB_XY2LIN(coeffDim, i, 0); for (j = 0; j < coeffDim; j++, coeffIdx++) { if ((i >= WAVELET_PADD2_APPROXDIM || j >= WAVELET_PADD2_APPROXDIM ) && !coeffIsZero(waveletMatrix [i] [j]) ) { /* Nonzero detail coefficient; set up pointer to coeff (sorts faster than 3 doubles) and save original (linearised) matrix index */ threshCoeffPtr -> idx = coeffIdx; threshCoeffPtr++ -> coeff = (WAVELET_COEFF*)&( waveletMatrix [i] [j] ); } } } /* Num of nonzero detail coeffs in buffer, number actually expected */ numThresh = threshCoeffPtr - threshCoeffs; nNzDetailCoeffs = WAVELET_PADD2_NUMDETAIL( preCompContrib -> nNonZeroCoeffs ); nCompressedCoeffs = preCompContrib -> nCompressedCoeffs; /* If fewer nonzero detail coeff are in the threshold buffer than * anticipated, the loop below fills the remainder of the threshold * buffer with duplicates of a coefficient in the lower right of the * matrix, which is padding and guaranteed to be zero. This condition * can occur if the wavelet transform actually generated genuine zero * detail coefficients. Infact it's quite common if the wavelet * transformed contributions are quite sparse. */ i = j = coeffDim - 1; coeffIdx = PMAP_CONTRIB_XY2LIN(coeffDim, i, j); for (coeffPtr = (WAVELET_COEFF*)(waveletMatrix [i][j]); numThresh < nNzDetailCoeffs; numThresh++ ) { threshCoeffPtr -> idx = coeffIdx; threshCoeffPtr++ -> coeff = coeffPtr; } /* Partition coeffs in reverse order, such that all coeffs in threshCoeffs [0..nCompressedCoeffs-1] are larger than those in threshCoeffs [nCompressedCoeffs..nNzDetailCoeffs-1] */ coeffPartition(threshCoeffs, nCompressedCoeffs - 1, 0, nNzDetailCoeffs - 1 ); #ifdef PMAP_CONTRIB_DBG /* Check coefficient partitioning */ threshCoeffPtr = threshCoeffs + nCompressedCoeffs - 1; for (i = 0; i < nCompressedCoeffs; i++) if (COEFFVAL(threshCoeffs + i) < COEFFVAL(threshCoeffPtr)) error(CONSISTENCY, "invalid wavelet coefficient partitioning " "in thresholdContribs()" ); for (; i < nNzDetailCoeffs; i++) if (COEFFVAL(threshCoeffs + i) > COEFFVAL(threshCoeffPtr)) error(CONSISTENCY, "invalid wavelet coefficient partitioning " "in thresholdContribs()" ); #endif /* Sort partition containing nCompressedCoeffs coefficients by index * (ignoring the remaining coefficients since these are now dropped * due to tresholding) */ qsort(threshCoeffs, nCompressedCoeffs, sizeof(PreComputedContribCoeff), coeffIdxCompare ); return 0; } static int encodeContribs (PreComputedContrib *preCompContrib, float compressRatio ) /* Apply wavelet transform to input matrix preCompContrib -> waveletMatrix and compress according to compressRatio, storing thresholded and mRGBE-encoded coefficients in preCompContrib -> mrgbeCoeffs. Note that the approximation coefficients in the upper left of the matrix are not encoded, and returned in preCompContrib -> waveletMatrix [0..WAVELET_PADD2_APPROXDIM-1] [0..WAVELET_PADD2_APPROXDIM-1]. Returns 0 on success. */ { unsigned i, j, k, scDim, lastCoeffIdx; WaveletMatrix2 waveletMatrix, tWaveletMatrix; PreComputedContribCoeff *threshCoeffs; mRGBERange *mrgbeRange; mRGBE *mrgbeCoeffs; WAVELET_COEFF absCoeff; #ifdef PMAP_CONTRIB_DBG WaveletCoeff3 decCoeff; unsigned decIdx; #endif if (!preCompContrib || !preCompContrib -> mrgbeCoeffs || !preCompContrib -> compressedCoeffs || !(waveletMatrix = preCompContrib -> waveletMatrix) || !(tWaveletMatrix = preCompContrib -> tWaveletMatrix) || !(scDim = preCompContrib -> scDim) ) /* Should be initialised by the caller! */ return -1; #ifdef PMAP_CONTRIB_DBG for (i = 0; i < scDim; i++) { for (j = 0; j < scDim; j++) { for (k = 0; k < 3; k++) { #if 0 /* Set contributions to bins for debugging */ waveletMatrix [i] [j] [k] = PMAP_CONTRIB_VAL( PMAP_CONTRIB_XY2LIN(scDim, i, j) ); #elif 0 /* Replace contribs with "bump" function */ waveletMatrix [i] [j] [k] = PMAP_CONTRIB_VAL( (1. - fabs(i - scDim/2. + 0.5) / (scDim/2. - 0.5)) * (1. - fabs(j - scDim/2. + 0.5) / (scDim/2. - 0.5)) ); #elif 0 /* Inject reference contribs from rcontrib classic */ #include "rc-ref.c" if (preCompContrib -> nBins != PMAP_CONTRIB_REFSIZE) { sprintf(errmsg, "reference contribs require %d bins", PMAP_CONTRIB_REFSIZE ); error(USER, errmsg); } waveletMatrix [i] [j] [k] = PMAP_CONTRIB_VAL( refContrib [PMAP_CONTRIB_XY2LIN(scDim, i, j)] [k] ); #endif } #if 0 /* Dump contribs prior to encoding for debugging */ printf("% 7.3g\t", colorAvg(waveletMatrix [i] [j])); } putchar('\n'); } putchar('\n'); #else } } #endif #endif /* Do 2D wavelet transform on preCompContrib -> waveletMatrix */ if (padWaveletXform2(waveletMatrix, tWaveletMatrix, scDim, NULL) != preCompContrib -> coeffDim ) error(INTERNAL, "failed wavelet transform in encodeContribs()"); /* Compress wavelet detail coeffs by thresholding */ if (thresholdContribs(preCompContrib) < 0) error(INTERNAL, "failed wavelet compression in encodeContribs()"); threshCoeffs = preCompContrib -> compressedCoeffs; /* Init per-channel coefficient range for mRGBE encoding */ mrgbeRange = &preCompContrib -> mrgbeRange; setcolor(mrgbeRange -> min, FHUGE, FHUGE, FHUGE); setcolor(mrgbeRange -> max, 0, 0, 0); /* Update per-channel coefficient range */ for (i = 0; i < preCompContrib -> nCompressedCoeffs; i++) for (k = 0; k < 3; k++) { absCoeff = fabs(threshCoeffs [i].coeff [k]); if (absCoeff < mrgbeRange -> min [k]) mrgbeRange -> min [k] = absCoeff; if (absCoeff > mrgbeRange -> max [k]) mrgbeRange -> max [k] = absCoeff; } if (preCompContrib -> nCompressedCoeffs == 1) /* Maximum compression with just 1 (!) compressed detail coeff (is * this even useful?), so mRGBE range is undefined since min & max * coincide; set minimum to 0, maximum to the single remaining * coefficient */ setcolor(mrgbeRange -> min, 0, 0, 0); /* Init mRGBE coefficient normalisation from range */ if (!mRGBEinit(mrgbeRange, mrgbeRange -> min, mrgbeRange -> max)) return -1; mrgbeCoeffs = preCompContrib -> mrgbeCoeffs; /* Encode wavelet detail coefficients to mRGBE */ for (i = lastCoeffIdx = 0; i < preCompContrib -> nCompressedCoeffs; lastCoeffIdx = threshCoeffs [i++].idx ) { /* HACK: To reduce the likelihood of overflowing the mRGBE data * field with the coefficient index, it is expressed incrementally * w.r.t. the previously encoded coefficient's index, instead of * absolutely. This implies threshCoeffs must be sorted by * coefficient index to ensure increments are positive, and to * minimise their magnitude. * Note that an overflow cannot be predicted beforehand, e.g. by * mkpmap when parsing the number of bins, as this depends on the * sparseness of the wavelet coefficients (which in turn depends on * the frequency distribution of the binned contributions), and the * fraction of those that are dropped (i.e. the compression ratio). * */ mrgbeCoeffs [i] = mRGBEencode(threshCoeffs [i].coeff, mrgbeRange, threshCoeffs [i].idx - lastCoeffIdx ); if (!mrgbeCoeffs [i].all) error(INTERNAL, "failed mRGBE encoding in encodeContribs()"); #ifdef PMAP_CONTRIB_DBG /* mRGBE encoding sanity check */ decIdx = mRGBEdecode(mrgbeCoeffs [i], mrgbeRange, decCoeff); mrgbeDiffs += sqrt(dist2(decCoeff, threshCoeffs [i].coeff)); mrgbeDiffCnt++; if (/*decIdx > preCompContrib -> nCompressedCoeffs || */ decIdx != threshCoeffs [i].idx - lastCoeffIdx || sqrt(dist2(decCoeff, threshCoeffs [i].coeff)) > 0.1 * colorAvg(mrgbeRange -> max) ) { sprintf(errmsg, "failed sanity check in encodeContribs()\n" "Encoded: [%.3g %.3g %.3g %d]\nDecoded: [%.3g %.3g %.3g %d]", threshCoeffs [i].coeff [0], threshCoeffs [i].coeff [1], threshCoeffs [i].coeff [2], threshCoeffs [i].idx - lastCoeffIdx, decCoeff [0], decCoeff [1], decCoeff [2], decIdx ); error(CONSISTENCY, errmsg); } for (k = 0; k < 3; k++) { absCoeff = fabs(threshCoeffs [i].coeff [k]) > 1e-6 ? decCoeff [k] / threshCoeffs [i].coeff [k] : 0; if (absCoeff < 0) error(CONSISTENCY, "coefficient sign reversal in encodeContribs()" ); } #endif } return 0; } static MODCONT *getPhotonContrib (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 *modCont = (MODCONT*)lu_find( pmapContribTab, srcMod -> oname ) -> data; Photon *photon; COLOR flux; DCOLOR *cbin; PhotonSearchQueueNode *sqn; float r2, norm; unsigned i, bin, numZeroBins; setcolor(irrad, 0, 0, 0); if (!modCont) /* Modifier not in LUT */ return NULL; /* Zero bins for this source modifier */ memset(modCont -> cbin, 0, sizeof(DCOLOR) * modCont -> nbins); if (!pmap -> maxGather) /* Zero bandwidth */ return NULL; /* Lookup photons; 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 from source %s found at (%.3f, %.3f, %.3f)", source [ray->rsrc].so -> oname, ray -> rop [0], ray -> rop [1], ray -> rop [2] ); error(WARNING, errmsg); #endif return NULL; } /* Avg radius^2 between furthest two photons to improve accuracy */ sqn = pmap -> squeue.node + 1; r2 = max(sqn -> dist2, (sqn + 1) -> dist2); r2 = 0.25 * (pmap -> maxDist2 + r2 + 2 * sqrt(pmap -> maxDist2 * r2)); #ifdef PMAP_EPANECHNIKOV /* Normalise accumulated flux by Epanechnikov kernel integral in 2D (see Eq. 4.4, p.76 in Silverman, "Density Estimation for Statistics and Data Analysis", 1st Ed., 1986, and Wann Jensen, "Realistic Image Synthesis using Photon Mapping"), include RADIANCE-specific lambertian factor PI */ norm = 2 / (PI * PI * r2); #else /* Normalise accumulated flux by search area PI * r^2, including RADIANCE-specific lambertian factor PI */ norm = 1 / (PI * PI * r2); #endif /* Skip the extra photon */ for (i = 1, numZeroBins = modCont -> nbins; i < pmap -> squeue.tail; i++, sqn++ ) { /* Get photon's contribution to density estimate */ photon = getNearestPhoton(&pmap -> squeue, sqn -> idx); getPhotonFlux(photon, flux); scalecolor(flux, norm); #ifdef PMAP_EPANECHNIKOV /* Apply Epanechnikov kernel to photon flux based on distance */ scalecolor(flux, 1 - sqn -> dist2 / r2); #endif addcolor(irrad, flux); cbin = &modCont -> cbin [photonSrcBin(pmap, photon)]; /* Bin was empty (zero) before accumulating flux; decrement corresponding counter */ numZeroBins -= coeffIsZero(*cbin); addcolor(*cbin, flux); } #ifdef PMAP_CONTRIB_WARNZEROBINS /* Warn if too many empty bins */ if (numZeroBins > PMAP_CONTRIB_MAXZEROBINS * modCont -> nbins) { sprintf(errmsg, "empty bin ratio > %.2f at (%.3f, %.3f, %.3f) " "for modifier %s; contributions may be biased", PMAP_CONTRIB_MAXZEROBINS, ray -> rop [0], ray -> rop [1], ray -> rop [2], modCont -> modname ); error(WARNING, errmsg); } #endif return modCont; } void freePreCompContribNode (void *p) /* Clean up precomputed contribution LUT entry */ { PhotonMap *preCompContribPmap = (PhotonMap*)p; PreComputedContrib *preCompContrib; if (preCompContribPmap) { preCompContrib = (PreComputedContrib*)( preCompContribPmap -> preCompContrib ); if (preCompContrib) { /* Free primary and transposed wavelet matrices */ freeWaveletMatrix2(preCompContrib -> waveletMatrix, preCompContrib -> coeffDim ); freeWaveletMatrix2(preCompContrib -> tWaveletMatrix, preCompContrib -> coeffDim ); /* Free thresholded coefficients */ free(preCompContrib -> compressedCoeffs); /* Free mRGBE encoded coefficients */ free(preCompContrib -> mrgbeCoeffs); free(preCompContrib -> waveletFname); if (preCompContrib -> cache) { /* Free contribution cache */ OOC_DeleteCache(preCompContrib -> cache); free(preCompContrib -> cache); } #if NIX if (preCompContrib -> numPhotonsShm) { /* Release shared memory for photon counter */ munmap(preCompContrib -> numPhotonsShm, sizeof(unsigned long)); close(preCompContrib -> shmFile); unlink(preCompContrib -> shmFname); } #endif #ifdef PMAP_CONTRIB_BINHISTO free(preCompContrib -> binHisto); #endif } /* Clean up precomputed contrib photon map */ deletePhotons(preCompContribPmap); free(preCompContribPmap); } } void initPreComputedContrib (PreComputedContrib *preCompContrib) /* Initialise precomputed contribution container in photon map */ { preCompContrib -> waveletFname = NULL; preCompContrib -> waveletFile = NULL; preCompContrib -> waveletMatrix = preCompContrib -> tWaveletMatrix = NULL; preCompContrib -> compressedCoeffs = NULL; setcolor(preCompContrib -> mrgbeRange.min, 0, 0, 0); setcolor(preCompContrib -> mrgbeRange.max, 0, 0, 0); setcolor(preCompContrib -> mrgbeRange.norm, 0, 0, 0); preCompContrib -> mrgbeCoeffs = NULL; preCompContrib -> scDim = preCompContrib -> nBins = preCompContrib -> coeffDim = preCompContrib -> nCoeffs = preCompContrib -> nNonZeroCoeffs = preCompContrib -> nCompressedCoeffs = preCompContrib -> contribSize = 0; preCompContrib -> cache = NULL; preCompContrib -> numPhotonsShm = NULL; preCompContrib -> shmFile = -1; #ifdef PMAP_CONTRIB_BINHISTO preCompContrib -> binHisto = NULL; #endif } static int allocPreCompContribPmap (const LUENT *modContNode, void *p) /* Allocate & init child contrib photon map for modContNode's modifier * and its contributions for parent photon map p. */ { unsigned scDim, nBins, coeffDim, nCoeffs, nCompressedCoeffs, nNzDetailCoeffs; PhotonMap *preCompContribPmap, *parentPmap = (PhotonMap*)p; LUENT *preCompContribNode; PreComputedContrib *preCompContrib; const MODCONT *modCont = (MODCONT*)modContNode -> data; if (!modCont || !parentPmap) { error(INTERNAL, "undefined parent photon map or modifier LUT node" " in allocPreCompContribPmap()" ); return -1; } preCompContribNode = lu_find(parentPmap -> preCompContribTab, modCont -> modname ); if (!preCompContribNode) error(SYSTEM, "out of memory allocating LUT entry in " "allocPreCompContribPmap()" ); preCompContribNode -> key = (char*)modCont -> modname; preCompContribNode -> data = (char*)(preCompContribPmap = malloc(sizeof(PhotonMap)) ); if (!preCompContribPmap) error(SYSTEM, "out of memory allocating photon map " "in allocPreCompContribPmap()" ); initPhotonMap(preCompContribPmap, PMAP_TYPE_CONTRIB_CHILD); /* Set output filename from parent photon map */ preCompContribPmap -> fileName = parentPmap -> fileName; /* Set contrib/coeff mode from parent photon map */ preCompContribPmap -> contribMode = parentPmap -> contribMode; preCompContrib = (PreComputedContrib*)( preCompContribPmap -> preCompContrib = malloc(sizeof(PreComputedContrib)) ); if (!preCompContrib) error(SYSTEM, "out of memory allocating contributions " "in allocPreCompContribPmap()" ); initPreComputedContrib(preCompContrib); #if NIX /* Allocate and map shared memory for photon counter */ strcpy(preCompContrib -> shmFname, PMAP_TMPFNAME); preCompContrib -> shmFile = mkstemp(preCompContrib -> shmFname); if (preCompContrib -> shmFile < 0 || ftruncate(preCompContrib -> shmFile, sizeof(unsigned long)) < 0 ) error(SYSTEM, "failed opening shared memory file in allocPreCompContribPmap()" ); preCompContrib -> numPhotonsShm = mmap(NULL, sizeof(unsigned long), PROT_READ | PROT_WRITE, MAP_SHARED, preCompContrib -> shmFile, 0 ); if (preCompContrib -> numPhotonsShm == MAP_FAILED) error(SYSTEM, "failed mapping shared memory in allocPreCompContribPmap()" ); #endif /* Set Shirley-Chiu square & wavelet matrix dimensions (number of bins resp. coefficients). */ preCompContrib -> scDim = scDim = sqrt( preCompContrib -> nBins = nBins = modCont -> nbins ); if (scDim * scDim != nBins) /* Shouldn't happen; mkpmap took care of dis */ error(INTERNAL, "contribution bin count not integer square " "in allocPreCompContribPmap()" ); /* Allocate and initialise photon and contribution heaps */ initPhotonHeap(preCompContribPmap); if (nBins > 1) { /* BINNING ENABLED; set up wavelet xform & compression. Get dimensions of wavelet coefficient matrix after boundary padding, and number of nonzero coeffs. The number of compressed (detail) coeffs is fixed and based on the number of NONZERO coeffs (minus approx, see NOTE below) since zero coeffs are considered thresholded by default. !!! NOTE: THE APPROX COEFFS IN THE UPPER LEFT OF !!! THE MATRIX ARE _NEVER_ THRESHOLDED TO MINIMISE !!! COMPRESSION ARTEFACTS! THEY ARE ESSENTIAL FOR !!! RECONSTRUCTING THE ORIGINAL CONTRIBS. */ preCompContrib -> coeffDim = coeffDim = padWaveletXform2( NULL, NULL, scDim, &preCompContrib -> nNonZeroCoeffs ); preCompContrib -> nCoeffs = nCoeffs = coeffDim * coeffDim; nNzDetailCoeffs = WAVELET_PADD2_NUMDETAIL( preCompContrib -> nNonZeroCoeffs ); nCompressedCoeffs = (1 - compressRatio) * nNzDetailCoeffs; if (nCompressedCoeffs < 1) { error(WARNING, "maximum contribution compression, " "clamping number of wavelet coefficients in " "allocPreCompContribPmap()" ); nCompressedCoeffs = 1; } if (nCompressedCoeffs >= preCompContrib -> nCoeffs) { error(WARNING, "minimum contribution compression, " "clamping number of wavelet coefficients in " "allocPreCompContribPmap()" ); nCompressedCoeffs = nNzDetailCoeffs; } preCompContrib -> nCompressedCoeffs = nCompressedCoeffs; /* Lazily allocate primary & transposed wavelet coeff matrix */ preCompContrib -> waveletMatrix = allocWaveletMatrix2(coeffDim); preCompContrib -> tWaveletMatrix = allocWaveletMatrix2(coeffDim); if (!preCompContrib -> waveletMatrix || !preCompContrib -> tWaveletMatrix ) error(SYSTEM, "out of memory allocating wavelet coefficient " "matrix in allocPreCompContribPmap()" ); /* Lazily allocate thresholded detail coefficients */ preCompContrib -> compressedCoeffs = calloc(nNzDetailCoeffs, sizeof(PreComputedContribCoeff) ); /* Lazily alloc mRGBE-encoded compressed wavelet coeffs */ preCompContrib -> mrgbeCoeffs = calloc(nCompressedCoeffs, sizeof(mRGBE) ); if (!preCompContrib -> compressedCoeffs || !preCompContrib -> mrgbeCoeffs ) error(SYSTEM, "out of memory allocating compressed " "contribution buffer in allocPreCompContribPmap()" ); /* Set size of compressed contributions, in bytes */ preCompContrib -> contribSize = PMAP_CONTRIB_ENCSIZE( nCompressedCoeffs ); #ifdef PMAP_CONTRIB_BINHISTO preCompContrib -> binHisto = calloc(nBins, sizeof(unsigned long)); if (!preCompContrib -> binHisto) error(SYSTEM, "out of memory allocating contribution bin " "histogram in allocPreCompContribPmap()" ); else memset(preCompContrib -> binHisto, 0, (size_t)nBins * sizeof(unsigned long) ); preCompContrib -> minZero = nBins; preCompContrib -> maxZero = preCompContrib -> sumZero = 0; #endif } return 0; } static int sumNumPhotons (const LUENT *preCompContribNode, void *p) /* Set number of photons for this modifier's child contrib photon map * from the corresponding shared photon counter, and sum over all * modifiers in (unsigned long*)p */ { unsigned long *numPhotons = (unsigned long*)p; PhotonMap *preCompContribPmap = (PhotonMap*)( preCompContribNode -> data ); const PreComputedContrib *preCompContrib; if (!preCompContribPmap) return 0; else preCompContrib = (PreComputedContrib*)( preCompContribPmap -> preCompContrib ); #if NIX /* Lock shared photon counter file for reading */ shmLock(preCompContrib -> shmFile, F_RDLCK); #endif *numPhotons += (preCompContribPmap -> numPhotons = *preCompContrib -> numPhotonsShm ); #if NIX /* Release lock on shard photon counter file */ shmLock(preCompContrib -> shmFile, F_UNLCK); #endif return 0; } static int flushContribPhotonHeaps (const LUENT *preCompContribNode, void *p ) /* Flush photon heaps for this modifier's precomputed contrib pmap */ { PhotonMap *preCompContribPmap = (PhotonMap*)( preCompContribNode -> data ); flushPhotonHeap(preCompContribPmap); return 0; } void preComputeContrib (PhotonMap *pmap, unsigned numProc) /* Precompute contributions for a random subset of (finalGather * pmap -> numPhotons) photons, and return the per-modifier precomputed contribution photon maps in LUT preCompContribTab, discarding the original photons. */ { unsigned long p, numPreComp; unsigned proc; int stat; pid_t pid, procPids [PMAP_MAXPROC]; PhotonMap nuPmap; LUTAB lutInit = LU_SINIT(NULL, freePreCompContribNode); /* Init new parent photon map and set output filename and rcontrib options */ initPhotonMap(&nuPmap, pmap -> type); nuPmap.fileName = pmap -> fileName; nuPmap.rcOpts = pmap -> rcOpts; /* Set contrib/coeff mode in new parent photon map */ nuPmap.contribMode = pmap -> contribMode; /* Allocate and init LUT containing per-modifier child photon maps */ if (!(nuPmap.preCompContribTab = malloc(sizeof(LUTAB)))) error(SYSTEM, "out of memory allocating LUT in preComputeContrib()" ); memcpy(nuPmap.preCompContribTab, &lutInit, sizeof(LUTAB)); /* Allocate and init per-modifier child photon maps */ lu_doall(pmapContribTab, allocPreCompContribPmap, &nuPmap); /* Record start time, baby */ repStartTime = repLastTime = time(NULL); repComplete = numPreComp = finalGather * pmap -> numPhotons; repProgress = 0; if (verbose) { sprintf(errmsg, "\nPrecomputing contributions for %ld photons\n", numPreComp ); eputs(errmsg); #if NIX fflush(stderr); #endif } 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 int photonOk; unsigned i, j, k, coeffIdx, scDim, nBins, coeffDim, nCoeffs, nCompressedCoeffs; const double pInc = finalGather > FTINY ? 1 / finalGather : 0, numPhotonsPerProc = (double)pmap -> numPhotons / numProc, numPreCompPerProc = (double)numPreComp / numProc; PhotonIdx pIdx; Photon photon; MODCONT *modCont = NULL; LUENT *preCompContribNode; PhotonMap *preCompContribPmap; PreComputedContrib *preCompContrib; WaveletMatrix2 waveletMatrix; WaveletCoeff3 *coeffPtr; RAY ray; unsigned long lastpIdx = 0; photonRay(NULL, &ray, PRIMARY, NULL); ray.ro = NULL; /* Decorrelate subprocess RNGs */ pmapSeed(randSeed + proc, pmap -> randState); #ifdef PMAP_CONTRIB_DBG /* Output child process PID after random delay to prevent * scrambling console output with sibling subprocs */ usleep(1e6 * pmapRandom(rouletteState)); sprintf(errmsg, "*DBG* Proc %d: PID = %d " "(waiting 10 sec to attach debugger...)\n", proc, getpid() ); eputs(errmsg); #if NIX fflush(stderr); #endif /* Give debugga time to attach to subproc */ sleep(10); #endif for (p = 0; p < (unsigned long)numPreCompPerProc; p++) { do { /* Get random photon from stratified distribution in * source heap to avoid duplicates and clustering. */ pIdx = firstPhoton(pmap) + (unsigned long)(p * pInc + pmapRandom(pmap -> randState) * max(pInc - 1, 0) + proc * numPhotonsPerProc ); getPhoton(pmap, pIdx, &photon); #ifdef PMAP_CONTRIB_DBG /* Shouldn't happen unless there's a bug in the pIdx calculation above. Skipped if modCont == NULL, i.e. if the last photon lookup failed */ if (modCont && pIdx == lastpIdx) { sprintf(errmsg, "duplicate contrib photon " "(proc %d, pIdx %d)", proc, pIdx ); error(WARNING, errmsg); } lastpIdx = pIdx; #endif /* 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 (i = 0; i < 3; i++) ray.ron [i] = photon.norm [i] / 127.0; /* Get contributions at photon position; retry with another photon if no contribs found */ modCont = getPhotonContrib(pmap, &ray, ray.rcol); } while (!modCont); /* Get preallocated LUT node for this modifier */ preCompContribNode = lu_find(nuPmap.preCompContribTab, modCont -> modname ); if (!preCompContribNode || !preCompContribNode -> key || !preCompContribNode -> data ) { sprintf(errmsg, "missing LUT entry for modifier %s" " in preComputeContrib()", modCont -> modname ); error(CONSISTENCY, errmsg); } preCompContribPmap = (PhotonMap*)preCompContribNode -> data; preCompContrib = (PreComputedContrib*)( preCompContribPmap -> preCompContrib ); scDim = preCompContrib -> scDim; nBins = preCompContrib -> nBins; nCoeffs = preCompContrib -> nCoeffs; nCompressedCoeffs = preCompContrib -> nCompressedCoeffs; waveletMatrix = preCompContrib -> waveletMatrix; if (nBins > 1) { /* Binning enabled */ char contribBuf0 [preCompContrib -> contribSize], *contribBuf; COLR *rgbe32 = (COLR*)contribBuf0; /* copy binned contribs to wavelet * input matrix before compressing and encoding */ for (i = 0; i < scDim; i++) { /* Calculate row pointer (=multiplication) in outer * loop, increment in inner loop */ coeffPtr = modCont -> cbin + PMAP_CONTRIB_XY2LIN(scDim, i, 0); #ifdef PMAP_CONTRIB_LOG /* Logarithmic contribs; copy & apply per element */ for (j = 0; j < scDim; j++, coeffPtr++) for (k = 0; k < 3; k++) waveletMatrix [i] [j] [k] = PMAP_CONTRIB_VAL((*coeffPtr) [k]); #else /* Linear contribs; copy matrix rows */ memcpy(waveletMatrix [i], *coeffPtr, scDim * WAVELET_COEFFSIZE ); #endif } /* Compress and encode contribs */ if (encodeContribs(preCompContrib, compressRatio)) { sprintf(errmsg, "failed contribution " "compression/encoding for modifier %s " "in preComputeContrib()", modCont -> modname ); error(INTERNAL, errmsg); } /* Encode wavelet approx coeffs in the upper left of * the wavelet matrix to 32-bit RGBE, and buffer them. * These coeffs are not thresholded to minimise * compression artefacts. */ for (i = 0; i < WAVELET_PADD2_APPROXDIM; i++) { for (j = 0; j < WAVELET_PADD2_APPROXDIM; j++, rgbe32++) { setcolr(*rgbe32, fabs(waveletMatrix [i] [j] [0]), fabs(waveletMatrix [i] [j] [1]), fabs(waveletMatrix [i] [j] [2]) ); /* HACK: depending on the boundary extension mode, * some approx coeffs can be NEGATIVE (!), which * 32-bit RGBE doesn't accommodate. To get around * this, we sacrifice a bit in each mantissa for the * sign. */ for (k = 0; k < 3; k++) (*rgbe32) [k] = PMAP_CONTRIB_SET_RGBE32_SGN( (*rgbe32) [k], waveletMatrix [i] [j] [k] ); } } /* Encode wavelet detail coeff range to 32-bit RGBE, and * buffer it after approx coeffs. * NOTE: mRGBE range maximum is buffered first to * facilitate handling the special (but pointless) case of * a single compressed detail coeff; if so, the mRGBE * minimum is omitted, since implicitly zero. */ setcolr(*(rgbe32++), preCompContrib -> mrgbeRange.max [0], preCompContrib -> mrgbeRange.max [1], preCompContrib -> mrgbeRange.max [2] ); if (nCompressedCoeffs > 1) setcolr(*(rgbe32++), preCompContrib -> mrgbeRange.min [0], preCompContrib -> mrgbeRange.min [1], preCompContrib -> mrgbeRange.min [2] ); contribBuf = (char*)rgbe32; /* Now buffa mRGBE encoded, compressed detail coeffs. * NOTE: we can't copy to preCompContribPmap -> * contribHeapBuf directly because it isn't allocated * until the 1st call to newPhoton() below. IMHO this is * cleaner and more consistent with the existing code than * allocating it here. */ memcpy(contribBuf, preCompContrib -> mrgbeCoeffs, nCompressedCoeffs * sizeof(mRGBE) ); /* Append precomputed photon to heap from ray, check if * photon was accepted (i.e. newPhoton() returned 0). * Note the locally buffered encoded contribs are passed * as an extra parameter. */ ray.rsrc = proc; photonOk = !newPhoton(preCompContribPmap, &ray, contribBuf0 ); #ifdef PMAP_CONTRIB_BINHISTO if (preCompContrib -> binHisto) { unsigned long numZero = 0; /* Update bin histogram for last photon lookup */ for (i = 1; i < pmap -> squeue.tail; i++) { Photon *photon = getNearestPhoton(&pmap -> squeue, pmap -> squeue.node [i].idx ); ++preCompContrib -> binHisto [ photonSrcBin(pmap, photon) ]; } for (i = 0; i < preCompContrib -> nBins; i++) numZero += (intens(modCont -> cbin [i]) < FTINY); preCompContrib -> sumZero += numZero; if (numZero < preCompContrib -> minZero) preCompContrib -> minZero = numZero; if (numZero > preCompContrib -> maxZero) preCompContrib -> maxZero = numZero; } #endif } else { /* Binning disabled; append precomputed photon to heap * from ray _without_ encoded contribs, check if photon * was accepted (i.e. newPhoton() returned 0). */ photonOk = !newPhoton(preCompContribPmap, &ray, NULL); } /* Increment parent photon map's photon counter if photon was * accepted by newPhoton(). */ nuPmap.numPhotons += photonOk; #if NIX /* Asynchronously increment shared photon counter for this * modifier if photon was accepted by newPhoton(). * NOTE: write lock on memmapped file! */ if (photonOk) { shmLock(preCompContrib -> shmFile, F_WRLCK); *preCompContrib -> numPhotonsShm += photonOk; shmLock(preCompContrib -> shmFile, F_UNLCK); } #else /* Synchronously report progress on W33nD0z */ if (!proc && photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime ) { repProgress = nuPmap.numPhotons; pmapPreCompReport(); } #endif } /* End precomp. photon loop */ /* Flush heap buffers to finalise them */ lu_doall(nuPmap.preCompContribTab, flushContribPhotonHeaps, NULL); #ifdef PMAP_CONTRIB_BINHISTO /* Print binning histogram/population count for 1st subproc */ if (!proc && preCompContrib -> binHisto) { printf("Binning histogram for modifier %s (avg/bin):\n", preCompContribNode -> key ); for (i = 0; i < preCompContrib -> nBins; i++) printf("%f\n", (float)preCompContrib -> binHisto [i] / preCompContrib -> nBins / nuPmap.numPhotons); printf("Empty bin ratio (avg, min, max):\t%f\t%f\t%f\n", (float)preCompContrib -> sumZero / preCompContrib -> nBins / nuPmap.numPhotons, (float)preCompContrib -> minZero / preCompContrib -> nBins, (float)preCompContrib -> maxZero / preCompContrib -> nBins ); } #endif #if NIX #ifdef PMAP_CONTRIB_DBG if (mrgbeDiffCnt) { sprintf(errmsg, "*DBG* preComputeContrib: Avg mRGBE encoding diff = %g\n", mrgbeDiffs / mrgbeDiffCnt ); eputs(errmsg); } sprintf(errmsg, "*DBG* Subproc %d exiting normally\n", proc); eputs(errmsg); + #if NIX + fflush(stderr); + #endif #endif /* Terminate subprocess normally */ exit(0); #endif } /* if (!fork()) */ else { /* PARENT PROCESS CONTINUES FORKING LOOP HERE */ if (pid < 0) error(SYSTEM, "failed to fork subprocess in preComputeContrib()" ); else /* Save child process ID */ procPids [proc] = pid; } } /* End forking loop */ #if NIX /* PARENT PROCESS CONTINUES HERE AFTER FORKING ALL SUBPROCS */ #ifdef SIGCONT /* Enable progress report signal handler */ signal(SIGCONT, pmapPreCompReport); #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 & quit */ for (proc = 0; proc < numProc; proc++) kill(procPids [proc], SIGKILL); cleanUpPmaps(photonMaps); error(USER, "failed precomputing contribution photons"); } --proc; } /* Nod off for a bit and update progress */ sleep(1); /* Update num photons from shared photon counters for all mods */ nuPmap.numPhotons = 0; lu_doall(nuPmap.preCompContribTab, sumNumPhotons, &nuPmap.numPhotons ); repProgress = nuPmap.numPhotons; if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime) pmapPreCompReport(); #ifdef SIGCONT else signal(SIGCONT, pmapPreCompReport); #endif } #endif /* NIX */ #ifdef SIGCONT /* Reset signal handler */ signal(SIGCONT, SIG_DFL); #endif /* Trash original pmap and its leaf file, and replace with new parent, which now acts as container for the per-modifier child pmaps */ unlink(pmap -> store.leafFname); deletePhotons(pmap); memcpy(pmap, &nuPmap, sizeof(PhotonMap)); /* Bail out if no photons could be precomputed */ if (!pmap -> numPhotons) error(USER, "no contribution photons precomputed; try increasing -am" ); } #else /* ------------------------------------------------------------------- U N I T T E S T S ------------------------------------------------------------------- */ #include <stdio.h> #include <random.h> #include "func.h" int main (int argc, char *argv []) { unsigned i, scdim, nsamp, numTheta = 0, numPhi = 0; double t, p; RAY ray; int bin; if (argc < 3) { fprintf(stderr, "%s <scdim> <nsamp> [<var>=<value>; ..]\n", argv [0] ); fprintf(stderr, "Default variable defs: %s\n", PMAP_CONTRIB_SCDEFAULTS ); fputs("\nMissing Shirley-Chiu dimension scDim>1, " "number of samples nsamp\n", stderr ); return -1; } if (!(scdim = atoi(argv [1]))) { fputs("Invalid Shirley-Chiu dimension\n", stderr); return -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(); calcontext(RCCONTEXT); #endif /* Compile default orientation variables for contribution binning */ scompile(PMAP_CONTRIB_SCDEFAULTS, NULL, 0); /* Compile custom orientation variables from command line */ for (i = 3; i < argc; i++) scompile(argv [i], NULL, 0); puts("x\ty\tz\t\ttheta\tphi\t\tbin"); 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, scdim); printf("%.3f\t%.3f\t%.3f\t-->\t%.3f\t%.3f\t-->\t%d\n", ray.rdir [0], ray.rdir [1], ray.rdir [2], t, p, bin ); } return 0; } #endif #endif /* PMAP_CONTRIB */ diff --git a/pmapdata.c b/pmapdata.c index 68d1c18..7db313d 100644 --- a/pmapdata.c +++ b/pmapdata.c @@ -1,993 +1,999 @@ #ifndef lint static const char RCSid[] = "$Id: pmapdata.c,v 2.23 2021/04/14 11:26:25 rschregle Exp $"; #endif /* ========================================================================= Photon map types and high level interface to nearest neighbour lookups in underlying point cloud data structure. The default data structure is an in-core kd-tree (see pmapkdt.{h,c}). This can be overriden with the PMAP_OOC compiletime switch, which enables an out-of-core octree (see oococt.{h,c}). Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) (c) Fraunhofer Institute for Solar Energy Systems, supported by the German Research Foundation (DFG LU-204/10-2, "Fassadenintegrierte Regelsysteme" (FARESYS)) (c) Lucerne University of Applied Sciences and Arts, supported by the Swiss National Science Foundation (SNSF #147053, "Daylight Redirecting Components", SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment") (c) Tokyo University of Science, supported by the JSPS Grants-in-Aid for Scientific Research (KAKENHI JP19KK0115, "Three-Dimensional Light Flow") (c) Tokyo University of Science, supported by the KAJIMA Foundation under the project title: "Reflections from Building Façades and their Glare Potential on the Built Environment -- Application of the Photon Flow Method using Annual Simulation" ========================================================================= $Id: pmapdata.c,v 2.23 2021/04/14 11:26:25 rschregle Exp $ */ #include "pmapdata.h" #include "pmapdens.h" #include "pmaprand.h" #include "pmapmat.h" #include "pmaproi.h" #include "pmapcontrib.h" #include "otypes.h" #include "random.h" PhotonMap *photonMaps [NUM_PMAP_TYPES] = { NULL, NULL, NULL, NULL, NULL, NULL, NULL, #ifdef PMAP_PHOTONFLOW NULL, NULL, #endif NULL }; /* Include routines to handle underlying point cloud data structure */ #ifdef PMAP_OOC #include "pmapooc.c" #else #include "pmapkdt.c" #include "pmaptkdt.c" #endif /* Ambient include/exclude set (from ambient.c) */ #ifndef MAXASET #define MAXASET 4095 #endif extern OBJECT ambset [MAXASET+1]; /* Callback to print photon attributes acc. to user defined format */ int (*printPhoton)(RAY *r, Photon *p, PhotonMap *pm); /* PHOTON MAP BUILD ROUTINES ------------------------------------------ */ void initPhotonMap (PhotonMap *pmap, PhotonMapType t) /* Init photon map 'n' stuff... */ { if (!pmap) return; pmap -> numPhotons = 0; pmap -> biasCompHist = NULL; pmap -> maxPos [0] = pmap -> maxPos [1] = pmap -> maxPos [2] = -FHUGE; pmap -> minPos [0] = pmap -> minPos [1] = pmap -> minPos [2] = FHUGE; pmap -> minGathered = pmap -> maxGathered = pmap -> totalGathered = 0; pmap -> gatherTolerance = gatherTolerance; pmap -> minError = pmap -> maxError = pmap -> rmsError = 0; pmap -> numDensity = 0; pmap -> distribRatio = 1; pmap -> type = t; pmap -> squeue.node = NULL; pmap -> squeue.len = 0; #ifdef PMAP_PATHFILT /* Init path lookup table and its key buffer */ pmap -> pathLUT = NULL; pmap -> pathLUTKeys = NULL; pmap -> numPathLUTKeys = 0; #endif /* Init stats */ setcolor(pmap -> photonFlux, 0, 0, 0); pmap -> CoG [0] = pmap -> CoG [1] = pmap -> CoG [2] = 0; pmap -> CoGdist = 0; pmap -> maxDist2 = 0; /* Init transient photon mapping stuff */ pmap -> velocity = photonVelocity; pmap -> minPathLen = pmap -> maxPathLen = pmap -> avgPathLen = 0; /* Init local RNG state */ pmap -> randState [0] = 10243; pmap -> randState [1] = 39829; pmap -> randState [2] = 9433; pmapSeed(randSeed, pmap -> randState); /* Set up type-specific photon lookup callback */ pmap -> lookup = pmapLookup [t]; /* Init storage */ pmap -> heap = NULL; pmap -> heapBuf = NULL; pmap -> heapBufLen = 0; /* Init precomputed contrib photon stuff */ pmap -> contribMode = 0; pmap -> preCompContribTab = NULL; pmap -> contribHeap = NULL; pmap -> contribHeapBuf = NULL; pmap -> preCompContrib = NULL; pmap -> rcOpts = NULL; pmap -> numContribSrc = 0; pmap -> contribSrc = NULL; /* Mark photon contribution source as unused */ pmap -> lastContribSrc.srcIdx = pmap -> lastContribSrc.srcBin = -1; #ifdef PMAP_OOC OOC_Null(&pmap -> store); #else kdT_Null(&pmap -> store); #endif } void initPhotonHeap (PhotonMap *pmap) { int fdFlags; if (!pmap) error(INTERNAL, "undefined photon map in initPhotonHeap()"); if (!pmap -> heap) { /* Open heap file */ mktemp(strcpy(pmap -> heapFname, PMAP_TMPFNAME)); if (!(pmap -> heap = fopen(pmap -> heapFname, "w+b"))) error(SYSTEM, "failed opening heap file in initPhotonHeap()"); #ifdef F_SETFL /* XXX is there an alternative needed for W33nd0z? */ fdFlags = fcntl(fileno(pmap -> heap), F_GETFL); fcntl(fileno(pmap -> heap), F_SETFL, fdFlags | O_APPEND); #endif /* ftruncate(fileno(pmap -> heap), 0); */ } #ifdef PMAP_CONTRIB if (PMAP_CONTRIB_BINNING(pmap)) { /* Allocate precomputed contribution heap buffer */ if (!pmap -> contribHeap) { /* Open heap file for precomputed contributions */ mktemp(strcpy(pmap -> contribHeapFname, PMAP_TMPFNAME)); if (!(pmap -> contribHeap = fopen(pmap -> contribHeapFname, "w+b"))) error(SYSTEM, "failed opening precomputed contribution " "heap file in initPhotonHeap()" ); #ifdef F_SETFL /* XXX is there an alternative needed for W33nd0z? */ fdFlags = fcntl(fileno(pmap -> contribHeap), F_GETFL); fcntl(fileno(pmap -> contribHeap), F_SETFL, fdFlags | O_APPEND); #endif /* ftruncate(fileno(pmap -> contribHeap), 0); */ } } #endif } void flushPhotonHeap (PhotonMap *pmap) { int fd, fdContrib = -1; unsigned long len, lenContrib = 0; if (!pmap) error(INTERNAL, "undefined photon map in flushPhotonHeap()"); if (!pmap -> heap || !pmap -> heapBuf) { /* Silently ignore undefined heap error(INTERNAL, "undefined photon heap in flushPhotonHeap()"); */ return; } fd = fileno(pmap -> heap); len = pmap -> heapBufLen * sizeof(Photon); #ifdef PMAP_CONTRIB if (PMAP_CONTRIB_BINNING(pmap)) { if (!pmap -> contribHeap || !pmap -> contribHeapBuf) { /* Silently ignore undefined heap error(INTERNAL, "undefined contribution heap in flushPhotonHeap()"); */ return; } /* Do same for the contribution photon buffa if binning is enabled. * Note the photon and contribution heap contain the same number of * entries to ensure photons are flushed in sync along with their * corresponding contributions to their respective heap files */ fdContrib = fileno(pmap -> contribHeap); lenContrib = pmap -> heapBufLen * ((PreComputedContrib*)pmap -> preCompContrib) -> contribSize; } #endif #if NIX /* Lock heap file(s) prior to writing */ shmLock(fd, F_WRLCK); if (lenContrib) shmLock(fdContrib, F_WRLCK); #endif #ifdef DEBUG_PMAP sprintf(errmsg, "Proc %d: flushing %ld photons from pos %ld\n", getpid(), pmap -> heapBufLen, lseek(fd, 0, SEEK_END) / sizeof(Photon) ); eputs(errmsg); #if NIX fflush(stderr); #endif #endif /* Atomically write write block to photon heap file */ /* !!! Unbuffered I/O via pwrite() avoids potential race conditions * !!! and buffer corruption which can occur with lseek()/fseek() * !!! followed by write()/fwrite(). */ /*if (pwrite(fd, pmap -> heapBuf, len, lseek(fd, 0, SEEK_END)) != len) */ if (write(fd, pmap -> heapBuf, len) != len) { /* Clean up temp file */ fclose(pmap -> heap); unlink(pmap -> heapFname); error(SYSTEM, "can't append to photon heap file in flushPhotonHeap()"); } #ifdef PMAP_CONTRIB /* Atomically write block to contribution heap file */ if (fdContrib > -1 && write(fdContrib, pmap -> contribHeapBuf, lenContrib) != lenContrib ) { /* Clean up temp file */ fclose(pmap -> contribHeap); unlink(pmap -> contribHeapFname); error(SYSTEM, "can't append to contrib heap file in flushPhotonHeap()"); } #endif #if NIX if (fsync(fd) || fdContrib > -1 && fsync(fdContrib)) error(SYSTEM, "failed to fsync heap in flushPhotonHeap()"); /* Release lock on heap file(s) */ shmLock(fd, F_UNLCK); if (fdContrib > -1) shmLock(fdContrib, F_UNLCK); #endif pmap -> heapBufLen = 0; } #ifdef DEBUG_PMAP static int checkPhotonHeap (FILE *file) /* Check heap for nonsensical or duplicate photons */ { Photon p, lastp; int i, dup; rewind(file); memset(&lastp, 0, sizeof(lastp)); while (fread(&p, sizeof(p), 1, file)) { dup = 1; for (i = 0; i <= 2; i++) { if (p.pos [i] < thescene.cuorg [i] || p.pos [i] > thescene.cuorg [i] + thescene.cusize ) { sprintf(errmsg, "corrupt photon in heap at [%f, %f, %f]\n", p.pos [0], p.pos [1], p.pos [2] ); error(WARNING, errmsg); } dup &= p.pos [i] == lastp.pos [i]; } if (dup) { sprintf(errmsg, "consecutive duplicate photon in heap " "at [%f, %f, %f]\n", p.pos [0], p.pos [1], p.pos [2] ); error(WARNING, errmsg); } } return 0; } #endif int newPhoton (PhotonMap* pmap, const RAY* ray, const char *contribs) { unsigned i; Photon photon; COLOR photonFlux; /* Account for distribution ratio */ if (!pmap || pmapRandom(pmap -> randState) > pmap -> distribRatio) return -1; /* Don't store on sources */ if (ray -> robj > -1 && islight(objptr(ray -> ro -> omod) -> otype)) return -1; /* Ignore photon if modifier in/outside exclude/include set */ if (ambincl != -1 && ray -> ro && ambincl != inset(ambset, ray->ro->omod) ) return -1; /* Store photon if within a region of interest (for ze Ecksperts!) */ if (!photonInROI(ray)) return -1; /* Adjust flux according to distribution ratio and ray weight */ copycolor(photonFlux, ray -> rcol); #if 0 /* Factored out ray -> rweight as deprecated (?) for pmap, and infact erroneously attenuates volume photon flux based on extinction, which is already factored in by photonParticipate() */ scalecolor(photonFlux, ray -> rweight / (pmap -> distribRatio ? pmap -> distribRatio : 1) ); #else scalecolor(photonFlux, 1.0 / (pmap -> distribRatio ? pmap -> distribRatio : 1) ); #endif setPhotonFlux(&photon, photonFlux); /* Set photon position and flags */ VCOPY(photon.pos, ray -> rop); photon.flags = 0; photon.caustic = PMAP_CAUSTICRAY(ray); /* Set photon's subprocess index (to consolidate contribution sources after photon distribution completes) */ photon.proc = PMAP_GETRAYPROC(ray); switch (pmap -> type) { #ifdef PMAP_CONTRIB case PMAP_TYPE_CONTRIB: /* Contrib photon before precomputation (i.e. in forward pass); set contribution source from last index in contrib source array. Note the contribution source bin has already been set by newPhotonContribSrc(). */ photon.aux.contribSrc = pmap -> numContribSrc; /* Photon will be stored; set contribution source index, * thereby marking it as having spawned photon(s) */ if (ray -> rsrc > PMAP_MAXSRCIDX) error(INTERNAL, "contribution source index overflow"); else pmap -> lastContribSrc.srcIdx = ray -> rsrc; break; case PMAP_TYPE_CONTRIB_CHILD: /* Contrib photon being precomputed; set contribution source from index passed in ray. NOTE: This is currently only for info and consistency with the above, but not used by the lookup routines */ if (ray -> rsrc > PMAP_MAXSRCIDX) error(INTERNAL, "contribution source index overflow"); else photon.aux.contribSrc = ray -> rsrc; break; #endif case PMAP_TYPE_TRANSIENT: #ifdef PMAP_PHOTONFLOW case PMAP_TYPE_TRANSLIGHTFLOW: #endif /* Set cumulative path length for transient photon, corresponding to its time of flight. The cumulative path length is obtained relative to the maximum path length photonMaxPathDist. The last intersection distance is subtracted as this isn't factored in until photonRay() generates a new photon ray. */ photon.aux.pathLen = photonMaxPathDist - ray -> rmax + ray -> rot; break; default: /* Set photon's path ID from ray serial num (supplemented by photon.proc below) */ photon.aux.pathID = ray -> rno; } /* Set normal */ for (i = 0; i <= 2; i++) photon.norm [i] = 127.0 * (isVolumePmap(pmap) #ifdef PMAP_PHOTONFLOW || isLightFlowPmap(pmap) #endif ? ray -> rdir [i] : ray -> ron [i] ); if (!pmap -> heapBuf) { /* Lazily allocate photon heap buffa */ #if NIX /* Randomise buffa size to temporally decorellate flushes in * multiprocessing mode */ srandom(randSeed + getpid()); pmap -> heapBufSize = PMAP_HEAPBUFSIZE * 0.5 * (1 + frandom()); #else /* Randomisation disabled for single processes on WIN; also useful * for reproducability during debugging */ pmap -> heapBufSize = PMAP_HEAPBUFSIZE; #endif #ifdef PMAP_CONTRIB if (contribs && pmap -> preCompContrib) /* Adapt buffa size to accommodate precomputed contributions */ pmap -> heapBufSize *= (float)sizeof(Photon) / ((PreComputedContrib*)pmap -> preCompContrib) -> contribSize; #endif if (!(pmap -> heapBuf = calloc(pmap -> heapBufSize, sizeof(Photon)))) error(SYSTEM, "out of memory allocating heap buffer in newPhoton()" ); pmap -> heapBufLen = 0; } /* Photon initialised; now append to heap buffa */ memcpy(pmap -> heapBuf + pmap -> heapBufLen, &photon, sizeof(Photon)); #ifdef PMAP_CONTRIB if (contribs && pmap -> preCompContrib) { const unsigned long contribSize = ((PreComputedContrib*)pmap -> preCompContrib) -> contribSize; if (!pmap -> contribHeapBuf) { /* Lazily allocate contribution heap buffa of same size as photon * heap buffa */ pmap -> contribHeapBuf = calloc(pmap -> heapBufSize, contribSize); if (!pmap -> contribHeapBuf) error(SYSTEM, "out of memory allocating precomputed " "contribution heap buffer in newPhoton()" ); } /* Append precomputed contribs to heap buffa */ memcpy(pmap -> contribHeapBuf + pmap -> heapBufLen * contribSize, contribs, contribSize ); } #endif if (++pmap -> heapBufLen >= pmap -> heapBufSize) /* Photon (and contrib, if applicable) heap buffa(s) full, flush to heap file(s) */ flushPhotonHeap(pmap); pmap -> numPhotons++; /* Print photon attributes */ if (printPhoton) /* Non-const kludge */ printPhoton((RAY*)ray, &photon, pmap); return 0; } void buildPhotonMap (PhotonMap *pmap, double *photonFlux, PhotonContribSourceIdx *contribSrcOfs, unsigned nproc ) { unsigned long n, nCheck = 0; unsigned i; Photon *p; COLOR flux; char nuHeapFname [sizeof(PMAP_TMPFNAME)]; FILE *nuHeap; /* Need double here to reduce summation errors */ double avgFlux [3] = {0, 0, 0}, CoG [3] = {0, 0, 0}, CoGdist = 0, avgPathLen = 0; FVECT d; if (!pmap) error(INTERNAL, "undefined photon map in buildPhotonMap()"); /* Get number of photons from heapfile size */ if (fseek(pmap -> heap, 0, SEEK_END) < 0) error(SYSTEM, "failed seek to end of photon heap in buildPhotonMap()"); pmap -> numPhotons = ftell(pmap -> heap) / sizeof(Photon); - if (!pmap -> numPhotons) + if (pmap -> numPhotons < 2) + /* Bail out if no photons (obviously). There is also little point in + having only one photon, and this will infact also cause numerical + issues when calculating the photon distribution bounds later on, + since its extent then collapses to zero. */ error(INTERNAL, "empty photon map in buildPhotonMap()"); if (!pmap -> heap) error(INTERNAL, "no heap in buildPhotonMap()"); #ifdef DEBUG_PMAP eputs("Checking photon heap consistency...\n"); checkPhotonHeap(pmap -> heap); sprintf(errmsg, "Heap contains %ld photons\n", pmap -> numPhotons); eputs(errmsg); #endif /* Allocate heap buffa */ if (!pmap -> heapBuf) { pmap -> heapBufSize = PMAP_HEAPBUFSIZE; pmap -> heapBuf = calloc(pmap -> heapBufSize, sizeof(Photon)); if (!pmap -> heapBuf) error(SYSTEM, "failed to allocate postprocessed photon heap in buildPhotonMap()" ); } /* We REALLY don't need yet another @%&*! heap just to hold the scaled * photons, but can't think of a quicker fix... */ mktemp(strcpy(nuHeapFname, PMAP_TMPFNAME)); if (!(nuHeap = fopen(nuHeapFname, "w+b"))) error(SYSTEM, "failed to open postprocessed photon heap in buildPhotonMap()" ); rewind(pmap -> heap); #ifdef DEBUG_PMAP eputs("Postprocessing photons...\n"); #endif while (!feof(pmap -> heap)) { #ifdef DEBUG_PMAP sprintf(errmsg, "Reading %lu at %lu... ", pmap -> heapBufSize, ftell(pmap->heap) ); #endif pmap -> heapBufLen = fread(pmap -> heapBuf, sizeof(Photon), pmap -> heapBufSize, pmap -> heap ); #ifdef DEBUG_PMAP sprintf(errmsg, "%s Got %lu\n", errmsg, pmap -> heapBufLen); eputs(errmsg); #if NIX fflush(stderr); #endif #endif if (ferror(pmap -> heap)) error(SYSTEM, "failed to read photon heap in buildPhotonMap()"); for (n = pmap -> heapBufLen, p = pmap -> heapBuf; n; n--, p++) { /* Update min and max pos and set photon flux */ for (i = 0; i <= 2; i++) { if (p -> pos [i] < pmap -> minPos [i]) pmap -> minPos [i] = p -> pos [i]; - else if (p -> pos [i] > pmap -> maxPos [i]) + if (p -> pos [i] > pmap -> maxPos [i]) pmap -> maxPos [i] = p -> pos [i]; - /* Update centre of gravity with photon position */ CoG [i] += p -> pos [i]; } if (isContribPmap(pmap) && contribSrcOfs) /* Linearise contrib photon origin index from subprocess index * using the per-subprocess offsets in contribSrcOfs */ p -> aux.contribSrc += contribSrcOfs [p -> proc]; else if (isTransientPmap(pmap) #ifdef PMAP_PHOTONFLOW || isTransLightFlowPmap(pmap) #endif ) { /* Update min and max path length */ if (p -> aux.pathLen < pmap -> minPathLen) pmap -> minPathLen = p -> aux.pathLen; else if (p -> aux.pathLen > pmap -> maxPathLen) pmap -> maxPathLen = p -> aux.pathLen; /* Update avg photon path length */ avgPathLen += p -> aux.pathLen; } /* Scale photon's flux (hitherto normalised to 1 over RGB); in * case of a contrib photon map, this is done per light source, * and photonFlux is assumed to be an array */ getPhotonFlux(p, flux); if (photonFlux) { scalecolor(flux, #if 1 photonFlux [isContribPmap(pmap) ? photonSrcIdx(pmap, p) : 0] #else photonFlux [0] #endif ); setPhotonFlux(p, flux); } /* Update average photon flux; need a double here */ addcolor(avgFlux, flux); } /* Write modified photons to new heap */ fwrite(pmap -> heapBuf, sizeof(Photon), pmap -> heapBufLen, nuHeap); if (ferror(nuHeap)) { error(SYSTEM, "failed postprocessing photon flux in buildPhotonMap()" ); /* Clean up temp files */ fclose(pmap -> heap); fclose(nuHeap); unlink(pmap -> heapFname); unlink(nuHeapFname); } nCheck += pmap -> heapBufLen; } #if NIX /* !!! Flush pending writes to new heap !!! */ fflush(nuHeap); #endif #ifdef DEBUG_PMAP if (nCheck < pmap -> numPhotons) error(INTERNAL, "truncated photon heap in buildPhotonMap"); #endif /* Finalise average photon flux */ scalecolor(avgFlux, 1.0 / pmap -> numPhotons); copycolor(pmap -> photonFlux, avgFlux); /* Average photon positions to get centre of gravity */ for (i = 0; i < 3; i++) pmap -> CoG [i] = CoG [i] /= pmap -> numPhotons; /* Average photon path lengths */ pmap -> avgPathLen = avgPathLen /= pmap -> numPhotons; rewind(pmap -> heap); /* Compute average photon distance to centre of gravity */ while (!feof(pmap -> heap)) { pmap -> heapBufLen = fread(pmap -> heapBuf, sizeof(Photon), pmap -> heapBufSize, pmap -> heap ); for (n = pmap -> heapBufLen, p = pmap -> heapBuf; n; n--, p++) { VSUB(d, p -> pos, CoG); CoGdist += DOT(d, d); if (isTransientPmap(pmap) #ifdef PMAP_PHOTONFLOW || isTransLightFlowPmap(pmap) #endif ) { /* Add temporal distance (in units of photon path length) for transient photons */ d [0] = p -> aux.pathLen - avgPathLen; CoGdist += d [0] * d [0]; } } } pmap -> CoGdist = CoGdist /= pmap -> numPhotons; /* Swap heaps, discarding unscaled photons */ fclose(pmap -> heap); unlink(pmap -> heapFname); pmap -> heap = nuHeap; strcpy(pmap -> heapFname, nuHeapFname); #ifdef PMAP_OOC OOC_BuildPhotonMap(pmap, nproc); #else if (isTransientPmap(pmap) #ifdef PMAP_PHOTONFLOW || isTransLightFlowPmap(pmap) #endif ) kdT_TransBuildPhotonMap(pmap); else kdT_BuildPhotonMap(pmap); #endif /* Trash heap and its buffa */ free(pmap -> heapBuf); fclose(pmap -> heap); unlink(pmap -> heapFname); pmap -> heap = NULL; pmap -> heapBuf = NULL; } /* PHOTON MAP SEARCH ROUTINES ------------------------------------------ */ /* Dynamic max photon search radius increase and reduction factors */ #define PMAP_MAXDIST_INC 4 #define PMAP_MAXDIST_DEC 0.9 /* Num successful lookups before reducing in max search radius */ #define PMAP_MAXDIST_CNT 1000 /* Threshold below which we assume increasing max radius won't help */ #define PMAP_SHORT_LOOKUP_THRESH 1 /* Coefficient for adaptive maximum search radius */ #define PMAP_MAXDIST_COEFF 100 void findPhotons (PhotonMap* pmap, const RAY* ray) { int redo = 0; const RREAL *norm = ray -> ron, *pos = ray -> rop; if (!pmap -> squeue.len) { /* Lazy init priority queue */ #ifdef PMAP_OOC OOC_InitFindPhotons(pmap); #else kdT_InitFindPhotons(pmap); #endif pmap -> minGathered = pmap -> maxGather; pmap -> maxGathered = pmap -> minGather; pmap -> totalGathered = 0; pmap -> numLookups = pmap -> numShortLookups = 0; pmap -> shortLookupPct = 0; pmap -> minError = FHUGE; pmap -> maxError = -FHUGE; pmap -> rmsError = 0; /* SQUARED max search radius limit is based on avg photon distance to * centre of gravity, unless fixed by user (maxDistFix > 0) */ pmap -> maxDist0 = pmap -> maxDist2Limit = maxDistFix > 0 ? maxDistFix * maxDistFix : PMAP_MAXDIST_COEFF * pmap -> squeue.len * pmap -> CoGdist / pmap -> numPhotons; } do { pmap -> squeue.tail = 0; pmap -> maxDist2 = pmap -> maxDist0; if (isVolumePmap(pmap) #ifdef PMAP_PHOTONFLOW || isLightFlowPmap(pmap) && sphericalIrrad #endif ) { /* Volume photons are not associated with a normal or intersection point; null the normal and set origin as lookup point. */ /* If a lightflow photon map is used and sphericalIrrad = 1, the mean spherical irradiance is evaluated, so normals are irrelevant and not filtered. */ norm = NULL; pos = ray -> rorg; } /* XXX/NOTE: status returned by XXX_FindPhotons() is currently ignored; if no photons are found, an empty queue is returned under the assumption all photons are too distant to contribute significant flux. */ #ifdef PMAP_OOC OOC_FindPhotons(pmap, pos, norm); #else if (isTransientPmap(pmap) #ifdef PMAP_PHOTONFLOW || isTransLightFlowPmap(pmap) #endif ) kdT_TransFindPhotons(pmap, pos, norm); else kdT_FindPhotons(pmap, pos, norm); #endif #ifdef PMAP_LOOKUP_INFO fprintf(stderr, "%d/%d %s photons found within radius %.3f " "at (%.2f,%.2f,%.2f) on %s\n", pmap -> squeue.tail, pmap -> squeue.len, pmapName [pmap -> type], sqrt(pmap -> maxDist2), ray -> rop [0], ray -> rop [1], ray -> rop [2], ray -> ro ? ray -> ro -> oname : "<null>" ); #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 : "<null>" ); error(WARNING, errmsg); #endif /* Bail out after warning if maxDist is fixed */ if (maxDistFix > 0) return; if (pmap -> maxDist0 < pmap -> maxDist2Limit) { /* Increase max search radius if below limit & redo search */ pmap -> maxDist0 *= PMAP_MAXDIST_INC; #ifdef PMAP_LOOKUP_REDO redo = 1; #endif #ifdef PMAP_LOOKUP_WARN sprintf(errmsg, redo ? "restarting photon lookup with max radius %.1e" : "max photon lookup radius adjusted to %.1e", pmap -> maxDist0 ); error(WARNING, errmsg); #endif } #ifdef PMAP_LOOKUP_REDO else { sprintf(errmsg, "max photon lookup radius clamped to %.1e", pmap -> maxDist0 ); error(WARNING, errmsg); } #endif } /* Reset successful lookup counter */ pmap -> numLookups = 0; } else { /* Skip search radius reduction if maxDist is fixed */ if (maxDistFix > 0) return; /* Increment successful lookup counter and reduce max search radius if * wraparound */ pmap -> numLookups = (pmap -> numLookups + 1) % PMAP_MAXDIST_CNT; if (!pmap -> numLookups) pmap -> maxDist0 *= PMAP_MAXDIST_DEC; redo = 0; } } while (redo); } Photon *find1Photon (PhotonMap *pmap, const RAY* ray, Photon *photon, void *photonIdx ) { - /* Init (squared) search radius to avg photon dist to centre of gravity */ - float maxDist2_0 = pmap -> CoGdist; - int found = 0; + /* Init (squared) search radius to avg photon dist to centre of gravity, + or fixed user-specified radius */ + float maxDist2_0 = maxDistFix > 0 + ? maxDistFix * maxDistFix + : PMAP_MAXDIST_COEFF * pmap -> CoGdist; + int found = 0; #ifdef PMAP_LOOKUP_REDO #define REDO 1 #else #define REDO 0 #endif do { pmap -> maxDist2 = maxDist2_0; #ifdef PMAP_OOC found = OOC_Find1Photon(pmap, ray -> rop, ray -> ron, photon, (OOC_DataIdx*)photonIdx ); #else found = kdT_Find1Photon(pmap, ray -> rop, ray -> ron, photon); #endif if (found < 0) { /* Expand search radius to retry */ maxDist2_0 *= 2; #ifdef PMAP_LOOKUP_WARN sprintf(errmsg, "failed 1-NN photon lookup" #ifdef PMAP_LOOKUP_REDO ", retrying with search radius %.2f", maxDist2_0 #endif ); error(WARNING, errmsg); #endif } } while (REDO && found < 0); /* Return photon buffer containing valid photon, else NULL */ return found < 0 ? NULL : photon; } void getPhoton (PhotonMap *pmap, PhotonIdx idx, Photon *photon) { #ifdef PMAP_OOC if (OOC_GetPhoton(pmap, idx, photon)) #else if (kdT_GetPhoton(pmap, idx, photon)) #endif error(INTERNAL, "failed photon lookup"); } Photon *getNearestPhoton (const PhotonSearchQueue *squeue, PhotonIdx idx) { #ifdef PMAP_OOC return OOC_GetNearestPhoton(squeue, idx); #else return kdT_GetNearestPhoton(squeue, idx); #endif } PhotonIdx firstPhoton (const PhotonMap *pmap) { #ifdef PMAP_OOC return OOC_FirstPhoton(pmap); #else return kdT_FirstPhoton(pmap); #endif } /* PHOTON MAP CLEANUP ROUTINES ------------------------------------------ */ void deletePhotons (PhotonMap* pmap) { #ifdef PMAP_OOC OOC_Delete(&pmap -> store); #else kdT_Delete(&pmap -> store); #endif free(pmap -> squeue.node); free(pmap -> biasCompHist); pmap -> numPhotons = pmap -> minGather = pmap -> maxGather = pmap -> squeue.len = pmap -> squeue.tail = 0; #ifdef PMAP_PATHFILT if (pmap -> pathLUT) { lu_done(pmap -> pathLUT); free(pmap -> pathLUT); } if (pmap -> pathLUTKeys) { unsigned k; for (k = 0; k < pmap->squeue.len + 1; free(pmap->pathLUTKeys [k++])); free(pmap -> pathLUTKeys); } #endif /* Contrib stuff */ if (isContribPmap(pmap) && pmap -> preCompContribTab) { /* Parent contrib photon map; clean up child photon maps containing precomputed contribs via lu_done() */ lu_done(pmap -> preCompContribTab); free(pmap -> preCompContribTab); pmap -> preCompContribTab = NULL; if (pmap -> contribHeapBuf) { free(pmap-> contribHeapBuf); pmap -> contribHeapBuf = NULL; } if (pmap -> contribHeap) { fclose(pmap -> contribHeap); pmap -> contribHeap = NULL; } } } diff --git a/pmapooc.c b/pmapooc.c index d07b4b9..4ece9d8 100644 --- a/pmapooc.c +++ b/pmapooc.c @@ -1,364 +1,364 @@ /* ====================================================================== Photon map interface to out-of-core octree Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) (c) Lucerne University of Applied Sciences and Arts, supported by the Swiss National Science Foundation (SNSF #147053, "Daylight Redirecting Components", SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment") (c) Tokyo University of Science, supported by the JSPS Grants-in-Aid for Scientific Research (KAKENHI JP19KK0115, "Three-Dimensional Light Flow") ====================================================================== $Id: pmapooc.c,v 1.7 2021/03/22 23:00:00 rschregle Exp $ */ #include "pmapdata.h" /* Includes pmapooc.h */ #include "source.h" #include "otspecial.h" #include "oocsort.h" #include "oocbuild.h" #include "pmcontribsort.h" /* PHOTON MAP BUILD ROUTINES ------------------------------------------ */ /* Returns photon position as sorting key for OOC_Octree & friends * (notably for Morton code generation). * !!! XXX: Uses type conversion from float via TEMPORARY storage; * !!! THIS IS NOT THREAD SAFE! * !!! RETURNED KEY PERSISTS ONLY IF COPIED BEFORE NEXT CALL! */ RREAL *OOC_PhotonKey (const void *p) { static FVECT photonPos; /* Temp storage for type conversion */ VCOPY(photonPos, ((Photon*)p) -> pos); return photonPos; } #ifdef DEBUG_OOC static void OOC_CheckKeys (FILE *file, const OOC_Octree *oct) { Photon p, lastp; RREAL *key; MortonIdx zIdx, lastzIdx = 0; rewind(file); memset(&lastp, 0, sizeof(lastp)); while (fread(&p, sizeof(p), 1, file) > 0) { key = OOC_PhotonKey(&p); zIdx = OOC_KEY2MORTON(key, oct); if (zIdx < lastzIdx) error(INTERNAL, "photons not sorted"); if (zIdx == lastzIdx) { sprintf(errmsg, "identical key %021ld " "for [%.9f, %.9f, %.9f]\tand [%.9f, %.9f, %.9f]", zIdx, lastp.pos [0], lastp.pos [1], lastp.pos [2], p.pos [0], p.pos [1], p.pos [2] ); error(WARNING, errmsg); } lastzIdx = zIdx; memcpy(&lastp, &p, sizeof(p)); } } #endif void OOC_BuildPhotonMap (struct PhotonMap *pmap, unsigned numProc) { FILE *leafFile; char *leafFname; FVECT d, octOrg; int i; RREAL octSize = 0; /* Determine octree size and origin from pmap extent and init octree */ VCOPY(octOrg, pmap -> minPos); VSUB(d, pmap -> maxPos, pmap -> minPos); for (i = 0; i < 3; i++) if (octSize < d [i]) octSize = d [i]; - + if (octSize < FTINY) error(INTERNAL, "zero octree size in OOC_BuildPhotonMap"); - + /* Derive leaf filename from photon map and open file */ leafFname = malloc( strlen(pmap -> fileName) + strlen(PMAP_OOC_LEAFSUFFIX) + 1 ); if (!leafFname) error(SYSTEM, "failed leaf filename alloc in OOC_BuildPhotonMap"); strcpy(leafFname, pmap -> fileName); strcat(leafFname, PMAP_OOC_LEAFSUFFIX); if (!(leafFile = fopen(leafFname, "w+b"))) error(SYSTEM, "failed to open leaf file in OOC_BuildPhotonMap"); #ifdef DEBUG_OOC eputs("Sorting photons by Morton code...\n"); #endif /* Sort photons in heapfile by Morton code and write to leaf file; redirect to contrib sorting routine if contrib child pmap and we have precomputed contributions, else we need a regular sort before we can precompute */ if (isContribChildPmap(pmap) && pmap -> preCompContrib ? contribPhotonSort(pmap, leafFile, octOrg, octSize, PMAP_OOC_NUMBLK, PMAP_OOC_BLKSIZE, numProc ) : OOC_Sort(pmap -> heap, leafFile, PMAP_OOC_NUMBLK, PMAP_OOC_BLKSIZE, numProc, sizeof(Photon), octOrg, octSize, OOC_PhotonKey ) ) error(INTERNAL, "failed out-of-core photon sort in OOC_BuildPhotonMap"); /* Init and build octree */ OOC_Init(&pmap -> store, sizeof(Photon), octOrg, octSize, OOC_PhotonKey, leafFile, leafFname ); #ifdef DEBUG_OOC eputs("Checking leaf file consistency...\n"); OOC_CheckKeys(leafFile, &pmap -> store); eputs("Building out-of-core octree...\n"); #endif if (!OOC_Build(&pmap -> store, PMAP_OOC_LEAFMAX, PMAP_OOC_MAXDEPTH)) error(INTERNAL, "failed out-of-core octree build in OOC_BuildPhotonMap"); #ifdef DEBUG_OOC eputs("Checking out-of-core octree consistency...\n"); if (OOC_Check(&pmap -> store, OOC_ROOT(&pmap -> store), octOrg, octSize, 0 )) error(INTERNAL, "inconsistent out-of-core octree; Time4Harakiri"); #endif } /* PHOTON MAP I/O ROUTINES ------------------------------------------ */ int OOC_SavePhotons (const struct PhotonMap *pmap, FILE *out) { return OOC_SaveOctree(&pmap -> store, out); } int OOC_LoadPhotons (struct PhotonMap *pmap, FILE *nodeFile) { FILE *leafFile; char leafFname [1024]; /* Derive leaf filename from photon map and open file */ strncpy(leafFname, pmap -> fileName, sizeof(leafFname)); strncat(leafFname, PMAP_OOC_LEAFSUFFIX, sizeof(leafFname) - strlen(leafFname) - 1 ); if (!(leafFile = fopen(leafFname, "r"))) error(SYSTEM, "failed to open leaf file in OOC_LoadPhotons"); if (OOC_LoadOctree(&pmap -> store, nodeFile, OOC_PhotonKey, leafFile)) return -1; #ifdef DEBUG_OOC /* Check octree for consistency */ if (OOC_Check(&pmap -> store, OOC_ROOT(&pmap -> store), pmap -> store.org, pmap -> store.size, 0 )) return -1; #endif return 0; } /* PHOTON MAP SEARCH ROUTINES ------------------------------------------ */ void OOC_InitFindPhotons (struct PhotonMap *pmap) { if (OOC_InitNearest(&pmap -> squeue, pmap -> maxGather + 1, sizeof(Photon) )) error(SYSTEM, "can't allocate photon search queue"); #ifdef PMAP_PATHFILT initPhotonPaths(pmap); #endif } static void OOC_InitPhotonCache (struct PhotonMap *pmap) /* Initialise OOC photon cache */ { static char warn = 1; if (!pmap -> store.cache && !pmap -> numDensity) { if (pmapCacheSize > 0) { const unsigned pageSize = pmapCachePageSize * pmap -> maxGather, numPages = pmapCacheSize / pageSize; /* Allocate & init I/O cache in octree */ pmap -> store.cache = malloc(sizeof(OOC_Cache)); if (!pmap -> store.cache || OOC_CacheInit(pmap -> store.cache, numPages, pageSize, sizeof(Photon) )) { error(SYSTEM, "failed OOC photon map cache init"); } } else if (warn) { error(WARNING, "OOC photon map cache DISABLED"); warn = 0; } } } int OOC_FindPhotons (struct PhotonMap *pmap, const FVECT pos, const FVECT norm ) { OOC_SearchFilterCallback filt; OOC_SearchFilterState filtState; #ifdef PMAP_PATHFILT OOC_SearchAttribCallback paths, *pathsPtr = NULL; #endif float res, n [3]; /* Lazily init OOC cache */ if (!pmap -> store.cache) OOC_InitPhotonCache(pmap); /* Set up filter callback */ if (norm) VCOPY(n, norm); filtState.pmap = pmap; filtState.pos = pos; filtState.norm = norm ? n : NULL; filt.state = &filtState; filt.filtFunc = filterPhoton; /* Get sought light source modifier from pmap -> lastContribSrc if precomputing contribution photon */ filtState.srcMod = isContribPmap(pmap) ? findmaterial(source [pmap -> lastContribSrc.srcIdx].so) : NULL; #ifdef PMAP_PATHFILT /* Set up path ID callback to filter for photons from unique paths */ paths.state = pmap; paths.findFunc = findPhotonPath; paths.addFunc = addPhotonPath; paths.delFunc = deletePhotonPath; paths.checkFunc = checkPhotonPaths; pathsPtr = &paths; resetPhotonPaths(pmap); res = OOC_FindNearest( &pmap -> store, OOC_ROOT(&pmap -> store), 0, pmap -> store.org, pmap -> store.size, pos, &filt, pathsPtr, &pmap -> squeue, pmap -> maxDist2 ); #else res = OOC_FindNearest( &pmap -> store, OOC_ROOT(&pmap -> store), 0, pmap -> store.org, pmap -> store.size, pos, &filt, NULL, &pmap -> squeue, pmap -> maxDist2 ); #endif /* PMAP_PATHFILT */ if (res < 0) error(INTERNAL, "failed k-NN photon lookup in OOC_FindPhotons"); /* Get (maximum distance)^2 from search queue head */ pmap -> maxDist2 = pmap -> squeue.node [0].dist2; /* Return success or failure (empty queue => none found) */ return pmap -> squeue.tail ? 0 : -1; } int OOC_Find1Photon (struct PhotonMap* pmap, const FVECT pos, const FVECT norm, Photon *photon, OOC_DataIdx *photonIdx ) { OOC_SearchFilterCallback filt; OOC_SearchFilterState filtState; float n [3], maxDist2; /* Lazily init OOC cache */ if (!pmap -> store.cache) OOC_InitPhotonCache(pmap); /* Set up filter callback */ if (norm) VCOPY(n, norm); filtState.pmap = pmap; filtState.norm = norm ? n : NULL; filtState.srcMod = NULL; filt.state = &filtState; filt.filtFunc = filterPhoton; *photonIdx = (OOC_DataIdx)(-1); maxDist2 = OOC_Find1Nearest(&pmap -> store, OOC_ROOT(&pmap -> store), 0, pmap -> store.org, pmap -> store.size,pos, &filt, photon, photonIdx, pmap -> maxDist2 ); if (maxDist2 < 0) error(INTERNAL, "failed 1-NN photon lookup in OOC_Find1Photon"); if (maxDist2 >= pmap -> maxDist2) /* No photon found => failed */ return -1; else { /* Set photon distance => success */ pmap -> maxDist2 = maxDist2; return 0; } } /* PHOTON ACCESS ROUTINES ------------------------------------------ */ int OOC_GetPhoton (struct PhotonMap *pmap, PhotonIdx idx, Photon *photon) { return OOC_GetData(&pmap -> store, idx, photon); } Photon *OOC_GetNearestPhoton (const PhotonSearchQueue *squeue, PhotonIdx idx) { return OOC_GetNearest(squeue, idx); } PhotonIdx OOC_FirstPhoton (const struct PhotonMap* pmap) { return 0; } diff --git a/pmcontrib3.c b/pmcontrib3.c index 1536f2c..d60fd39 100644 --- a/pmcontrib3.c +++ b/pmcontrib3.c @@ -1,266 +1,281 @@ #ifndef lint static const char RCSid[] = "$Id$"; #endif /* ========================================================================= Photon map routines for precomputed light source contributions. These routines build and save precomputed contribution photon maps, and are used by mkpmap. Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) (c) Lucerne University of Applied Sciences and Arts, supported by the Swiss National Science Foundation (SNSF #147053, "Daylight Redirecting Components", SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment") (c) Tokyo University of Science, supported by the KAJIMA Foundation under the project title: "Reflections from Building Façades and their Glare Potential on the Built Environment -- Application of the Photon Flow Method using Annual Simulation" ========================================================================= $Id$ */ #include "pmapcontrib.h" #ifdef PMAP_CONTRIB #include "pmapdiag.h" #include "pmapio.h" #include <sys/stat.h> #if NIX #define __USE_XOPEN_EXTENDED /* Enables nftw() under Linux */ #include <ftw.h> #endif /* ------------------ CONTRIB PMAP BUILDING STUFF --------------------- */ #if NIX static int ftwFunc (const char *path, const struct stat *statData, int typeFlags, struct FTW *ftwData ) /* Callback for nftw() to delete a directory entry */ { return remove(path); /* = unlink()/rmdir() depending on type */ } #endif static int buildPreCompContribPmap (const LUENT *preCompContribNode, void *p) /* Build per-modifier precomputed photon map. Empty photon maps are (nonfatally) skipped. Returns 0 on success. */ { PhotonMap *preCompContribPmap = (PhotonMap*)( preCompContribNode -> data ); PreComputedContrib *preCompContrib = (PreComputedContrib*)( preCompContribPmap -> preCompContrib ); const char *dirName = (char*)p; unsigned fnLen, wfnLen; /* Flush photon/contrib heaps */ flushPhotonHeap(preCompContribPmap); fflush(preCompContribPmap -> contribHeap); - /* Get heap file size, skip if zero (= no photons for this modifier) */ + /* Get heap file size */ if (fseek(preCompContribPmap -> heap, 0, SEEK_END) < 0) - error(SYSTEM, "failed seek to end of photon heap in buildPhotonMap()"); - if (!ftell(preCompContribPmap -> heap)) { - sprintf(errmsg, "No contribution photons precomputed for modifier %s", + error(SYSTEM, + "failed seek to end of photon heap in buildPreCompContribPmap()" + ); + if (ftell(preCompContribPmap -> heap) / sizeof(Photon) < 2) { + /* Skip if no photons for this modifier, as this causes + buildPhotonMap() to fail. Also skip if only one, since pointless + and will cause problems building the underlying data structure. */ + sprintf(errmsg, "no contribution photons precomputed for modifier %s", preCompContribNode -> key ); error(WARNING, errmsg); + + /* Trash heap and its buffa; this would have been taken care of by + buildPhotonMap(). */ + fclose(preCompContribPmap -> heap); + preCompContribPmap -> heap = NULL; + unlink(preCompContribPmap -> heapFname); + preCompContribPmap -> numPhotons = 0; } else { if (verbose) { sprintf(errmsg, "Building precomputed contribution photon map " "for modifier %s\n", preCompContribNode -> key ); eputs(errmsg); #if NIX fflush(stderr); #endif } /* Allocate & set output filenames to subdirectory and modifier */ fnLen = strlen(dirName) + strlen(preCompContribNode -> key) + strlen(PMAP_CONTRIB_FILE) + 1; wfnLen = strlen(dirName) + strlen(preCompContribNode -> key) + strlen(PMAP_CONTRIB_WAVELETFILE) + 1; if (!(preCompContribPmap -> fileName = malloc(fnLen)) || !(preCompContrib -> waveletFname = malloc(wfnLen)) ) error(SYSTEM, "out of memory allocating filename in " "buildPreCompContribPmap()" ); snprintf(preCompContribPmap -> fileName, fnLen, PMAP_CONTRIB_FILE, dirName, preCompContribNode -> key ); snprintf(preCompContrib -> waveletFname, wfnLen, PMAP_CONTRIB_WAVELETFILE, dirName, preCompContribNode -> key ); if (preCompContrib -> nBins > 1) { /* Binning enabled; open wavelet coefficient file; this is where the sorted contributions are written to by contribPhotonSort() (via buildPhotonMap()) */ if (!(preCompContrib -> waveletFile = fopen( preCompContrib -> waveletFname, "wb" )) ) { sprintf(errmsg, "can't open wavelet coefficient file %s", preCompContrib -> waveletFname ); error(SYSTEM, errmsg); } } /* Rebuild underlying data structure, destroying heap; photon flux is now * finalised, so no flux normalisation (note NULL flux param) */ buildPhotonMap(preCompContribPmap, NULL, NULL, 1); - - /* Free primary and transposed wavelet matrices; no longer needed */ + } + + if (preCompContrib -> nBins > 1) { + /* Binning enbled; free primary and transposed wavelet matrices, + as no longer needed */ freeWaveletMatrix2(preCompContrib -> waveletMatrix, preCompContrib -> coeffDim ); freeWaveletMatrix2(preCompContrib -> tWaveletMatrix, preCompContrib -> coeffDim ); - /* Free thresholded coefficients; no longer needed */ + /* Free thresholded coefficients, as no longer needed */ free(preCompContrib -> compressedCoeffs); preCompContrib -> waveletMatrix = preCompContrib -> tWaveletMatrix = NULL; preCompContrib -> compressedCoeffs = NULL; } return 0; } void buildContribPhotonMap (const PhotonMap *pmap) /* Build contribution pmap and its per-modifier precomputed chilren */ { char dirName [1024]; struct stat dirStat; /* Set subdirectory containing per-modified child pmaps from filename */ snprintf(dirName, sizeof(dirName), PMAP_CONTRIB_DIR, pmap -> fileName); /* Cleanup previous directory contents if necessary */ if (!stat(dirName, &dirStat)) { /* File/dir exists */ if (S_ISREG(dirStat.st_mode)) /* Found regular file; delete and skip rest */ unlink(dirName); else if (S_ISDIR(dirStat.st_mode)) { /* Found subdirectory; delete its contents, skipping symlinks */ #if NIX if (nftw(dirName, ftwFunc, 1, FTW_DEPTH | FTW_MOUNT | FTW_PHYS) < 0 ) { sprintf(errmsg, "failed cleanup of output directory %s", dirName ); error(SYSTEM, errmsg); } #else /* Apparently there's no nftw() under Wind0z, so just skip * cleanup. Tuff luck, Wind0z Weeniez! There's probably some * replacement for this, but we couldn't be buggered... */ #endif } else { /* Found neither regular file nor directory; whatever it is, * just stuff it and quit! */ sprintf(errmsg, "cannot remove existing output file %s", dirName ); error(SYSTEM, errmsg); } } /* Create new directory for per-modifier contribution photon maps */ if (mkdir(dirName, S_IFDIR | S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH) < 0) { sprintf(errmsg, "error creating output directory %s", dirName); error(SYSTEM, errmsg); } /* Build per-modifier precomputed pmaps from their contribution heaps. */ lu_doall(pmap -> preCompContribTab, buildPreCompContribPmap, dirName); } /* ------------------ CONTRIB PMAP SAVING STUFF --------------------- */ static int savePreCompContrib (const LUENT *preCompContribNode, void *p) /* Save child contrib photon map for current modifier; empty photon maps are (nonfatally) skipped. */ { PhotonMap *preCompContribPmap = (PhotonMap*)( preCompContribNode -> data ); PhotonMapArgz *argz = (PhotonMapArgz*)p; /* Skip empty photon maps */ if (preCompContribPmap -> numPhotons) /* NOTE: savePhotonMap() will not recurse and call saveContribPhotonMap() again because preCompContribPmap -> preCompContribTab == NULL */ savePhotonMap(preCompContribPmap, preCompContribPmap -> fileName, argz -> argc, argz -> argv ); return 0; } void saveContribPhotonMap (const PhotonMap *pmap, const char *fname, int argc, char **argv ) /* Save contribution pmap and its per-modifier precomputed children. NOTE: fname accepted to conform with savePhotonMap(), but currently ignored; this was already set in preComputedContrib() */ { PhotonMapArgz argz = {argc, argv}; char rcOptFname [1024]; FILE *rcOptFile; /* Save rcontrib options file */ if (!pmap -> rcOpts) error(INTERNAL, "undefined rcontrib options in saveContribPhotonMap()" ); if (!pmap -> fileName) error(INTERNAL, "undefined photon map filename in saveContribPhotonMap()" ); snprintf(rcOptFname, sizeof(rcOptFname), PMAP_CONTRIB_OPTFILE, pmap -> fileName ); if (!(rcOptFile = fopen(rcOptFname, "w")) || !fprintf(rcOptFile, "%s\n", pmap -> rcOpts) ) { sprintf(errmsg, "cannot write rcontrib options to %s", rcOptFname); error(SYSTEM, errmsg); } fclose(rcOptFile); /* Now iterate over per-modifier child pmaps */ lu_doall(pmap -> preCompContribTab, savePreCompContrib, &argz); } #endif /* PMAP_CONTRIB */