diff --git a/pmapcontrib.c b/pmapcontrib.c index 3104a0b..53da282 100644 --- a/pmapcontrib.c +++ b/pmapcontrib.c @@ -1,1115 +1,1120 @@ #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 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") ========================================================================= $Id: pmapcontrib.c,v 2.20 2021/02/16 20:06:06 greg Exp $ */ -/* TODO: Isolate sanity checks with PMAP_CONTRIB_DBG once stable! */ +/* TODO: + -- Rubbish decoded coeffs and invalid mRGBE range in decodeContribs() + -- Sporadic attempt to save empty/invalid pmap in savePhotonMap() +*/ #include "pmapcontrib.h" #include "pmapdiag.h" #include "pmaprand.h" #include "pmapmat.h" #include "pmapsrc.h" #include "pmutil.h" #include "otspecial.h" #include "otypes.h" #include "lookup.h" #ifdef PMAP_CONTRIB /* 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; 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 & PMAP_CONTRIB_SCDIM(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 = PMAP_CONTRIB_SCDIM(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; } /* ------------------ 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; 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; /* 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, 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 ); } 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 -------------------- */ - typedef struct { - WAVELET_COEFF *coeff; - unsigned bin; - } EncodedPreComputedContrib; + static int coeffCompare (const void *c1, const void *c2) /* Comparison function to REVERSE sort thresholded coefficients */ { const EncodedPreComputedContrib *tcoeff1 = c1, *tcoeff2 = c2; /* Use dot product as magnitude to compare _absolute_ values */ const WAVELET_COEFF v1 = DOT(tcoeff1 -> coeff, tcoeff1 -> coeff), v2 = DOT(tcoeff2 -> coeff, tcoeff2 -> coeff); if (v1 < v2) return 1; else if (v1 > v2) return -1; else return 0; } static int thresholdContribs (MODCONT *binnedContribs, PreComputedContrib *preCompContrib ) /* Threshold binned detail coefficients in binnedContribs -> cbin [1..] by keeping the (preCompContrib -> nCompressedBins) largest and returning these in preContrib -> compressedBins along with their original bin order. NOTE: binnedContribs -> cbin [0] is the average wavelet coefficient and excluded from thresholding. Returns 0 on success. */ { unsigned b; EncodedPreComputedContrib *threshCoeffs; if (!binnedContribs || !preCompContrib || - !preCompContrib -> compressedBins + !(threshCoeffs = preCompContrib -> compressedBins) ) /* Should be initialised by caller! */ return -1; /* Set up coefficient buffer for compression (thresholding), skipping * first coeff since it's the average and therefore not thresholded */ - threshCoeffs = (EncodedPreComputedContrib*)( - preCompContrib -> compressedBins - ); for (b = 0; b < binnedContribs -> nbins - 1; b++) { /* Set up pointers to coeffs (sorted faster than 3 floats/doubles) and remember original bin position prior to sorting. NOTE: The encoded bin position preserves the offset for the 1st coefficient i.e. the average! */ threshCoeffs [b].coeff = (WAVELET_COEFF*)( &(binnedContribs -> cbin [b + 1]) ); threshCoeffs [b].bin = b + 1; } /* REVERSE sort coeffs by magnitude; the non-thresholded coeffs will then start at threshCoeffs [0] */ qsort(threshCoeffs, binnedContribs -> nbins - 1, sizeof(EncodedPreComputedContrib), coeffCompare ); return 0; } static int encodeContribs (MODCONT *binnedContribs, PreComputedContrib *preCompContrib, float compressRatio ) /* Apply wavelet transform to binnedContribs -> cbin and compress according to compressRatio, storing thresholded and mRGBE-encoded coefficients in preCompContrib -> mrgbeCoeffs. Note that the average coefficient is not encoded, and returned in binnedContribs -> cbin [0]. Returns 0 on success. */ { unsigned b, i; EncodedPreComputedContrib *threshCoeffs; mRGBERange *mrgbeRange; mRGBE *mrgbeCoeffs; WAVELET_COEFF absCoeff; #ifdef PMAP_CONTRIB_DBG WaveletCoeff3 decCoeff; unsigned decBin; #endif if (!binnedContribs || !preCompContrib || !preCompContrib -> compressedBins || !preCompContrib -> mrgbeCoeffs ) /* Should be initialised by the caller! */ return -1; + #ifdef PMAP_CONTRIB_DBG + for (b = 0; b < binnedContribs -> nbins; b++) { + #if 1 + /* Set contributions to bins for debugging */ + setcolor(binnedContribs -> cbin [b], b, b, b); + #endif + + /* Dump contribs prior to encoding for debugging */ + printf("%.3g\t", colorAvg(binnedContribs -> cbin [b])); + } + + putchar('\n'); + #endif + /* Do 2D wavelet transform on preCompContrib -> waveletMatrix (which actually maps to binnedContribs -> cbin) */ if (waveletXform2(preCompContrib -> waveletMatrix, preCompContrib -> tWaveletMatrix, preCompContrib -> l ) < 0 ) error(INTERNAL, "failed wavelet transform in encodeContribs()"); /* Compress wavelet detail coeffs by thresholding */ if (thresholdContribs(binnedContribs, preCompContrib) < 0) error(INTERNAL, "failed wavelet compression in encodeContribs()"); - threshCoeffs = (EncodedPreComputedContrib*)( - preCompContrib -> compressedBins - ); + threshCoeffs = preCompContrib -> compressedBins; /* 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 (b = 0; b < preCompContrib -> nCompressedBins; b++) { for (i = 0; i < 3; i++) { #ifdef PMAP_CONTRIB_DBG #if 0 - /* Encode bin as coeff instead of actual coeff itself */ + /* Replace wavelet coefficient with bin for debugging, ordering + bins linearly */ threshCoeffs [b].coeff [i] = b + 1; + threshCoeffs [b].bin = b + 1; #endif #endif absCoeff = fabs(threshCoeffs [b].coeff [i]); if (absCoeff < mrgbeRange -> min [i]) mrgbeRange -> min [i] = absCoeff; if (absCoeff > mrgbeRange -> max [i]) mrgbeRange -> max [i] = absCoeff; } } if (preCompContrib -> nCompressedBins == 1) /* Maximum compression with just 1 (!) compressed bin (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 */ - mRGBEinit(mrgbeRange, mrgbeRange -> min, mrgbeRange -> max); + if (!mRGBEinit(mrgbeRange, mrgbeRange -> min, mrgbeRange -> max)) + return -1; + mrgbeCoeffs = preCompContrib -> mrgbeCoeffs; /* Encode wavelet detail coefficients to mRGBE */ for (b = 0; b < preCompContrib -> nCompressedBins; b++) { mrgbeCoeffs [b] = mRGBEencode(threshCoeffs [b].coeff, mrgbeRange, threshCoeffs [b].bin ); #ifdef PMAP_CONTRIB_DBG /* Encoding sanity check */ decBin = mRGBEdecode(mrgbeCoeffs [b], mrgbeRange, decCoeff); #if 0 if (decBin != threshCoeffs [b].bin || sqrt(dist2(decCoeff, threshCoeffs [b].coeff)) > 1.7 ) error(CONSISTENCY, "failed sanity check in encodeContribs()"); #endif for (i = 0; i < 3; i++) if (decCoeff [i] / threshCoeffs [b].coeff [i] < 0) error(CONSISTENCY, "coefficient sign reversal in encodeContribs()" ); #endif } return 0; } static void initContribHeap (PhotonMap *pmap) /* Initialise precomputed contribution heap */ { int fdFlags; 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 initContribHeap()" ); #ifdef F_SETFL /* XXX is there an alternative needed for Wind0z? */ fdFlags = fcntl(fileno(pmap -> contribHeap), F_GETFL); fcntl(fileno(pmap -> contribHeap), F_SETFL, fdFlags | O_APPEND); #endif /* ftruncate(fileno(pmap -> heap), 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 *contrib = (MODCONT*)lu_find( pmapContribTab, srcMod -> oname ) -> data; Photon *photon; COLOR flux; PhotonSearchQueueNode *sqn; float r2, norm; 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; } /* 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)); /* XXX: OMIT extra normalisation factor 1/PI for ambient calc? */ #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 ; 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 photon distance */ scalecolor(flux, 1 - sqn -> dist2 / r2); #endif addcolor(irrad, flux); - #ifdef PMAP_CONTRIB_DBG - /* Set photon flux to bin for debugging */ - setcolor(contrib -> cbin [photonSrcBin(pmap, photon)], - photonSrcBin(pmap, photon), - photonSrcBin(pmap, photon), - photonSrcBin(pmap, photon) - ); - #endif addcolor(contrib -> cbin [photonSrcBin(pmap, photon)], flux); } return contrib; } 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 */ free(preCompContrib -> waveletMatrix); freeWaveletMatrix(preCompContrib -> tWaveletMatrix, preCompContrib -> l ); /* Free thresholded coefficients */ free(preCompContrib -> compressedBins); /* Free mRGBE encoded coefficients */ free(preCompContrib -> mrgbeCoeffs); free(preCompContrib -> waveletFname); } /* 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 -> compressedBins = 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 -> nBins = preCompContrib -> nCompressedBins = preCompContrib -> l = preCompContrib -> scDim = preCompContrib -> contribSize = 0; } void preComputeContrib (PhotonMap *pmap) /* 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 i, numPreComp; unsigned b, j; long nCompressedBins; PhotonIdx pIdx; Photon photon; RAY ray; MODCONT *binnedContribs; COLR mrgbeRange32 [2]; LUENT *preCompContribNode; PhotonMap *preCompContribPmap, nuPmap; PreComputedContrib *preCompContrib; LUTAB lutInit = LU_SINIT(NULL, freePreCompContribNode); #ifdef PMAP_CONTRIB_DBG RAY dbgRay; #endif if (verbose) { sprintf(errmsg, "\nPrecomputing contributions for %ld photons\n", numPreComp ); eputs(errmsg); #if NIX fflush(stderr); #endif } /* Init new parent photon map and set output filename */ initPhotonMap(&nuPmap, pmap -> type); nuPmap.fileName = pmap -> fileName; /* 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)); /* Record start time, baby */ repStartTime = time(NULL); #ifdef SIGCONT signal(SIGCONT, pmapPreCompReport); #endif repComplete = numPreComp = finalGather * pmap -> numPhotons; repProgress = 0; photonRay(NULL, &ray, PRIMARY, NULL); ray.ro = NULL; for (i = 0; i < numPreComp; i++) { /* Get random photon from stratified distribution in source heap * to avoid duplicates and clustering */ pIdx = firstPhoton(pmap) + (unsigned long)( (i + pmapRandom(pmap -> randState)) / finalGather ); getPhoton(pmap, pIdx, &photon); /* Init dummy photon ray with intersection 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; #ifndef PMAP_CONTRIB_DBG /* Get contributions at photon position; retry with another photon if no contribs found */ binnedContribs = getPhotonContrib(pmap, &ray, ray.rcol); #else /* Get all contribs from same photon for debugging */ getPhoton(pmap, 0, &photon); memcpy(&dbgRay, &ray, sizeof(RAY)); VCOPY(dbgRay.rop, photon.pos); dbgRay.rsrc = photonSrcIdx(pmap, &photon); for (j = 0; j < 3; j++) dbgRay.ron [j] = photon.norm [j] / 127.0; binnedContribs = getPhotonContrib(pmap, &dbgRay, dbgRay.rcol); #endif if (binnedContribs) { /* Found contributions at photon position, so generate * precomputed photon */ if (!(preCompContribNode = lu_find(nuPmap.preCompContribTab, binnedContribs -> modname) ) ) error(SYSTEM, "out of memory allocating LUT entry in " "preComputeContrib()" ); if (!preCompContribNode -> key) { /* New LUT node for precomputed contribs for this modifier */ preCompContribNode -> key = (char*)binnedContribs -> modname; preCompContribNode -> data = (char*)( preCompContribPmap = malloc(sizeof(PhotonMap)) ); if (preCompContribPmap) { /* Init new child photon map and its contributions */ initPhotonMap(preCompContribPmap, PMAP_TYPE_CONTRIB_CHILD ); initPhotonHeap(preCompContribPmap); initContribHeap(preCompContribPmap); preCompContrib = (PreComputedContrib*)( preCompContribPmap -> preCompContrib = malloc(sizeof(PreComputedContrib)) ); initPreComputedContrib(preCompContrib); } if (!preCompContribPmap || !preCompContrib) error(SYSTEM, "out of memory allocating new photon map " "in preComputeContrib()" ); /* Set output filename from parent photon map */ preCompContribPmap -> fileName = nuPmap.fileName; /* Set contrib/coeff mode from parent photon map */ preCompContribPmap -> contribMode = nuPmap.contribMode; /* Get Shirley-Chiu square & wavelet matrix resolution, bins per dimension, and number of bins */ preCompContrib -> l = logDim(binnedContribs -> nbins); preCompContrib -> scDim = PMAP_CONTRIB_SCDIM( preCompContrib -> l ); preCompContrib -> nBins = preCompContrib -> scDim * preCompContrib -> scDim; if (preCompContrib -> nBins != binnedContribs -> nbins) /* Shouldn't happen */ error(INTERNAL, "bin count mismatch in preComputeContrib()" ); if (preCompContrib -> nBins > 1) { /* BINNING ENABLED; set up wavelet xform & compression. Number of compressed coeffs/bins is fixed (subtract one as the average coefficient is not thresholded) */ nCompressedBins = (preCompContrib -> nBins - 1) * (1 - compressRatio); if (nCompressedBins < 1) { error(WARNING, "maximum contribution compression, clamping number " "of bins in preComputeContrib()" ); nCompressedBins = 1; } if (nCompressedBins >= preCompContrib -> nBins) { error(WARNING, "minimum contribution compression, clamping number " "of bins in preComputeContrib()" ); nCompressedBins = preCompContrib -> nBins - 1; } preCompContrib -> nCompressedBins = nCompressedBins; /* Lazily allocate 2D primary wavelet matrix and map its rows to linear array binnedContribs -> cbin */ if (!(preCompContrib -> waveletMatrix = calloc( preCompContrib -> scDim, sizeof(WaveletCoeff3*)) ) ) error(SYSTEM, "out of memory allocating primary " "wavelet matrix in preComputeContrib()" ); for (b = 0; b < preCompContrib -> scDim; b++) /* Point to each row in existing 1D contrib array */ preCompContrib -> waveletMatrix [b] = &binnedContribs -> cbin [b * preCompContrib -> scDim]; /* Lazily allocate transposed wavelet matrix */ preCompContrib -> tWaveletMatrix = allocWaveletMatrix(preCompContrib -> l); if (!preCompContrib -> tWaveletMatrix) error(SYSTEM, "out of memory allocating transposed wavelet " "matrix in preComputeContrib()" ); /* Lazily allocate thresholded detail coefficients (minus the average coeff) */ preCompContrib -> compressedBins = calloc( preCompContrib -> nBins - 1, sizeof(EncodedPreComputedContrib) ); /* Lazily alloc mRGBE-encoded compressed wavelet coeffs */ preCompContrib -> mrgbeCoeffs = calloc( preCompContrib -> nCompressedBins, sizeof(mRGBE) ); if (!preCompContrib -> compressedBins || !preCompContrib -> mrgbeCoeffs ) error(SYSTEM, "out of memory allocating compressed " "coefficients in preComputeContrib()" ); /* Set size of compressed contributions, in bytes */ preCompContrib -> contribSize = PMAP_CONTRIB_ENCSIZE( preCompContrib -> nCompressedBins ); } } else { /* LUT node already exists for this modifier; use its data */ preCompContribPmap = (PhotonMap*)preCompContribNode -> data; preCompContrib = (PreComputedContrib*)( preCompContribPmap -> preCompContrib ); } if (binnedContribs -> nbins > 1) { /* Binning enabled; compress & encode binned contribs */ if (encodeContribs(binnedContribs, preCompContrib, compressRatio ) ) error(INTERNAL, "failed contribution " "compression/encoding in preComputeContrib()" ); /* Dump encoded bins to contrib heap file, prepended by mRGBE range in 32-bit RGBE. NOTE: minimum and maximum are reversed in the file to facilitate handling of single compressed bins (see below) */ setcolr(mrgbeRange32 [0], preCompContrib -> mrgbeRange.max [0], preCompContrib -> mrgbeRange.max [1], preCompContrib -> mrgbeRange.max [2] ); setcolr(mrgbeRange32 [1], preCompContrib -> mrgbeRange.min [0], preCompContrib -> mrgbeRange.min [1], preCompContrib -> mrgbeRange.min [2] ); /* Dump only mRGBE maximum if 1 compressed bin (maximum * compression), since minimum is implicitly zero and can * be omitted to save space */ putbinary(mrgbeRange32, sizeof(COLR), 1 + (preCompContrib -> nCompressedBins > 1), preCompContribPmap -> contribHeap ); if (putbinary(preCompContrib -> mrgbeCoeffs, sizeof(mRGBE), preCompContrib -> nCompressedBins, preCompContribPmap -> contribHeap ) != preCompContrib -> nCompressedBins ) error(SYSTEM, "failed writing to coefficients to " "contribution heap in preComputeContrib()" ); } /* Set photon flux to coarsest average wavelet coefficient NOTE: binnedContrib -> cbin [0] == preCompContrib -> waveletMatrix [0][0] == preCompContrib -> tWaveletMatrix [0][0]. IF BINNING IS DISABLED, this is the total contribution from all directions */ copycolor(ray.rcol, binnedContribs -> cbin [0]); /* HACK: signal newPhoton() to set precomputed photon's contribution source from ray -> rsrc */ /* XXX: DO WE STILL NEED THIS SHIT AFTER PRECOMP, SINCE CHILD PMAPS ARE PER-MODIFIER ANYWAY? */ preCompContribPmap -> lastContribSrc.srcIdx = -2; /* Append photon to new heap from ray and increment total count for all modifiers in parent photon map */ newPhoton(preCompContribPmap, &ray); nuPmap.numPhotons++; } /* Update progress */ repProgress++; if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime) pmapPreCompReport(); #ifdef SIGCONT else signal(SIGCONT, pmapPreCompReport); #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)); #ifdef SIGCONT signal(SIGCONT, SIG_DFL); #endif /* Build per-modifier precomputed photon maps from their contribution heaps */ lu_doall(pmap -> preCompContribTab, buildPreCompContribPmap, NULL); } #else /* ------------------------------------------------------------------- U N I T T E S T S ------------------------------------------------------------------- */ #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 defs: %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 = PMAP_CONTRIB_SCBINS(l); 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 #endif /* PMAP_CONTRIB */ diff --git a/pmapcontrib.h b/pmapcontrib.h index ff43136..7db1424 100644 --- a/pmapcontrib.h +++ b/pmapcontrib.h @@ -1,198 +1,204 @@ /* RCSid $Id: pmapcontrib.h,v 2.5 2016/05/17 17:39:47 rschregle Exp $ */ /* ========================================================================= Photon map routines for precomputed light source contributions; these routines interface to mkpmap and rcontrib. Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) (c) Lucerne University of Applied Sciences and Arts, supported by the Swiss National Science Foundation (SNSF #147053, "Daylight Redirecting Components", SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment") ========================================================================= $Id: pmapcontrib.h,v 2.5 2016/05/17 17:39:47 rschregle Exp $ */ #ifndef _PMAPCONTRIB_H #define _PMAPCONTRIB_H #if defined(PMAP_OOC) && !defined(PMAP_CONTRIB) #define PMAP_CONTRIB #endif #include "pmapdata.h" #include "wavelet2.h" #include "mrgbe.h" #ifndef MAXPROCESS #include "rcontrib.h" #endif #ifndef MAXMODLIST /* Max number of contributing sources */ #define MAXMODLIST 1024 #endif /* Filename templates for per-modifier contrib photon maps and wavelet coeffs; these are held in a separate subdirectory PMAP_CONTRIB_DIR */ #define PMAP_CONTRIB_DIR "%s.rc" #define PMAP_CONTRIB_FILE "%s/%s.pm" #define PMAP_CONTRIB_WAVELETFILE "%s/%s.wvt" /* The following variables can be specified to override the orientation of the Shirley-Chiu mapping (see also disk2square.cal). We use the built-in functions in disk2square.c for efficiency and in order to not depend on an external function file. These variables merely mimic the function file interace. RHS : right-hand-coordinate system (-1 for left-hand) rNx, rNy, rNz : surface normal Ux, Uy, Uz : up vector (defines phi = 0) */ #define PMAP_CONTRIB_SCRHS "RHS" #define PMAP_CONTRIB_SCNORMX "rNx" #define PMAP_CONTRIB_SCNORMY "rNy" #define PMAP_CONTRIB_SCNORMZ "rNz" #define PMAP_CONTRIB_SCUPX "Ux" #define PMAP_CONTRIB_SCUPY "Uy" #define PMAP_CONTRIB_SCUPZ "Uz" #define PMAP_CONTRIB_SCDEFAULTS ( \ "RHS=1; rNx=0; rNy=0; rNz=1; Ux=0; Uy=1; Uz=0;" \ ) /* Maximum Shirley-Chiu binning resolution l*/ #define PMAP_CONTRIB_MAXBINRES (mRGBE_DATABITS >> 1) /* Shirley-Chiu square dimensions and number of bins for resolution l */ #define PMAP_CONTRIB_SCDIM(l) (1 << (l)) #define PMAP_CONTRIB_SCBINS(l) (1 << ((l) << 1)) /* Size of encoded & compressed coefficients in wavelet file (in bytes), * as a function of the number of compressed bins nComp. NOTE: the * minimum is omitted from the mRGBE range if nComp == 1, since it is * implicitly zero! */ #define PMAP_CONTRIB_ENCSIZE(nComp) ( \ sizeof(COLR) * (1 + ((nComp) > 1)) + sizeof(mRGBE) * (nComp) \ ) - + + + /* Struct for wavelet coeff thresholding; saves original bin order */ + typedef struct { + WAVELET_COEFF *coeff; + unsigned bin; + } EncodedPreComputedContrib; typedef struct { - char *waveletFname; - FILE *waveletFile; - WaveletMatrix waveletMatrix, tWaveletMatrix; - /* Coeff compression/encoding-decoding/decompression buffer */ - void *compressedBins; - mRGBERange mrgbeRange; - mRGBE *mrgbeCoeffs; - unsigned l, scDim, nBins, nCompressedBins; - unsigned long contribSize; + char *waveletFname; + FILE *waveletFile; + WaveletMatrix waveletMatrix, tWaveletMatrix; + /* Wavelet coeff compression/encoding buffer */ + EncodedPreComputedContrib *compressedBins; + mRGBERange mrgbeRange; + mRGBE *mrgbeCoeffs; + unsigned l, scDim, nBins, nCompressedBins; + unsigned long contribSize; } PreComputedContrib; /* 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 */ extern LUTAB *pmapContribTab; /* Contribution/coefficient mode flag */ extern int *pmapContribMode; MODCONT *addContribModifier(LUTAB *contribTab, unsigned *numContribs, char *mod, char *binParm, int binCnt ); /* Add light source modifier mod to contribution lookup table contribsTab, and update numContribs. Return initialised contribution data for this modifier. */ 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. */ void freePreCompContribNode (void *p); /* Clean up precomputed contribution LUT entry */ void initPmapContribTab (LUTAB *contribTab, int *contribMode); /* Set contribution table and contrib/coeff mode flag (interface to rcmain.c, see also pmapContribTab above) */ void initPmapContrib (PhotonMap *pmap); /* Initialise contribution photon map and check modifiers. NOTE: pmapContribTab must be set before via initPmapContribTab() */ 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. */ 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. */ void initPreComputedContrib (PreComputedContrib *preCompContrib); /* Initialise precomputed contribution container in photon map */ void preComputeContrib (PhotonMap *pmap); /* Precompute contributions for a random subset of (finalGather * pmap -> numPhotons) photons, and init the per-modifier precomputed contribution photon maps in LUT pmap -> preCompContribTab, discarding the original photons. */ void distribPhotonContrib (PhotonMap *pmap, LUTAB *contribTab, unsigned numContribs, int *contribMode, unsigned numProc ); /* Emit photons with binned light source contributions, precompute their contributions and build photon map */ int buildPreCompContribPmap (const LUENT *preCompContribNode, void *p); /* Build per-modifier precomputed photon map. Returns 0 on success. */ void saveContribPhotonMap (const PhotonMap *pmap, const char *fname, int argc, char **argv ); /* Save contribution pmap and its per-modifier precomputed children */ void loadContribPhotonMap (PhotonMap *pmap, const char *fname); /* Load contribution pmap and its per-modifier precomputed children */ void getPreCompPhotonContrib (PhotonMap *pmap, RAY *ray, COLOR totalContrib ); /* Get precomputed light source contributions at ray -> rop for all specified modifiers and accumulate them in pmapContribTab (which maps to rcontrib's binning LUT). Also return total precomputed contribs. */ #endif diff --git a/pmcontrib4.c b/pmcontrib4.c index debbc9b..f85e6a7 100644 --- a/pmcontrib4.c +++ b/pmcontrib4.c @@ -1,571 +1,572 @@ #ifndef lint static const char RCSid[] = "$Id: pmcontrib2.c,v 2.5 2018/11/08 00:54:07 greg Exp $"; #endif /* ========================================================================= Photon map routines for precomputed light source contributions. These routines interface to rcontrib. Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) (c) Lucerne University of Applied Sciences and Arts, supported by the Swiss National Science Foundation (SNSF #147053, "Daylight Redirecting Components", SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment") ========================================================================= $Id$ */ #include "pmapcontrib.h" #ifdef PMAP_CONTRIB #include "pmapmat.h" #include "pmapsrc.h" #include "pmaprand.h" #include "pmapio.h" #include "pmapdiag.h" #include "otypes.h" #include "otspecial.h" #include "oocnn.h" void initPmapContribTab (LUTAB *contribTab, int *contribMode) /* Set contribution table (interface to rcmain.c) */ { pmapContribTab = contribTab; pmapContribMode = contribMode; } /* ------------------ CONTRIB PMAP LOADING STUFF --------------------- */ static int checkContribModifier (const LUENT *modContNode, void *p) /* Check for valid modifier in LUT entry */ { MODCONT *modCont = (MODCONT*)modContNode -> data; char *modName = (char*)modCont -> modname; OBJECT modObj = modifier(modName); const PhotonMap *pmap = (PhotonMap*)p; if (modObj == OVOID) { sprintf(errmsg, "invalid modifier %s", modName); error(USER, errmsg); } if (!islight(objptr(modObj) -> otype)) { sprintf(errmsg, "%s is not a light source modifier", modName); error(USER, errmsg); } /* Warn if number of bins exceeds lookup bandwidth (as this is * guaranteed to result in empty bins!) */ if (pmap -> maxGather < modCont -> nbins) { sprintf(errmsg, "contribution photon lookup bandwidth too low " "for modifier %s, some bins will be zero", modName ); error(WARNING, errmsg); } return 0; } void initPmapContrib (PhotonMap *pmap) /* Set up photon map contributions and get binning parameters */ { unsigned t; /* Avoid incompatible photon maps */ for (t = 0; t < NUM_PMAP_TYPES; t++) if (photonMaps [t] && t != PMAP_TYPE_CONTRIB) { sprintf(errmsg, "%s photon map does not support contributions", pmapName [t] ); error(USER, errmsg); } if (!pmapContribTab) error(INTERNAL, "contribution LUT not initialised in initPmapContrib()" ); /* Check for valid contribution modifiers */ lu_doall(pmapContribTab, checkContribModifier, pmap); /* Set callback to collect contributions */ if (pmap) { pmap -> lookup = getPreCompPhotonContrib; /* XXX: Need tighter tolerance for contribution photon lookups? */ pmap -> gatherTolerance = 1.0; } } static int loadPreCompContrib (const LUENT *modContNode, void *p) /* Load child contribution photon map for current contrib modifier entry */ { MODCONT *modCont = (MODCONT*)modContNode -> data; char *modName = (char*)modCont -> modname; PhotonMap *preCompContribPmap; const PhotonMap *parentPmap = (const PhotonMap*)p; LUENT *preCompContribNode; PreComputedContrib *preCompContrib; char dirName [1024]; /* Allocate new LUT entry for child precomputed contribution pmap */ preCompContribNode = lu_find(parentPmap -> preCompContribTab, modName); if (!preCompContribNode) error(SYSTEM, "out of memory allocating LUT entry in loadPreCompContrib()" ); preCompContribNode -> key = modName; /* Allocate child precomputed contribution photon map */ preCompContribNode -> data = (char*)( preCompContribPmap = malloc(sizeof(PhotonMap)) ); if (!preCompContribPmap) error(SYSTEM, "out of memory allocating precomputed contribution " "photon map in loadPreCompContrib()" ); /* Set subdirectory from parent photon map's filename */ snprintf(dirName, sizeof(dirName), PMAP_CONTRIB_DIR, parentPmap -> fileName ); /* Allocate & set child pmap's filename from subdir and modifier */ preCompContribPmap -> fileName = malloc(strlen(parentPmap -> fileName) + strlen(modName) + strlen(PMAP_CONTRIB_FILE) ); if (!preCompContribPmap -> fileName) error(SYSTEM, "out of memory allocating filename in " "loadPreCompContrib()" ); snprintf(preCompContribPmap -> fileName, sizeof(dirName) + 4, PMAP_CONTRIB_FILE, dirName, modName ); /* Load child precomp contrib photon map for current modifier; loadPhotonMap() will not recurse and call loadContribPhotonMap() again because preCompContribPmap -> preCompContribTab == NULL */ loadPhotonMap(preCompContribPmap, preCompContribPmap -> fileName); /* Pass parent pmap's lookup parameters on to child */ preCompContribPmap -> minGather = parentPmap -> minGather; preCompContribPmap -> maxGather = parentPmap -> maxGather; preCompContribPmap -> gatherTolerance = parentPmap -> gatherTolerance; /* Override contrib/coeff mode if it doesn't match precomputed */ if (*pmapContribMode != preCompContribPmap -> contribMode) { sprintf(errmsg, "photon map contains precomputed %s, overriding rcontrib setting", preCompContribPmap -> contribMode ? "contributions" : "coefficients" ); error(WARNING, errmsg); *pmapContribMode = preCompContribPmap -> contribMode; } /* NOTE: preCompContribPmap -> preCompContrib is initialised by * loadPhotonMap() */ preCompContrib = (PreComputedContrib*)( preCompContribPmap -> preCompContrib ); /* Set bins per dimension and number of bins */ preCompContrib -> scDim = PMAP_CONTRIB_SCDIM(preCompContrib -> l); preCompContrib -> nBins = preCompContrib -> scDim * preCompContrib -> scDim; /* Set wavelet coeff filename */ preCompContrib -> waveletFname = malloc(strlen(dirName) + strlen(modName) + strlen(PMAP_CONTRIB_WAVELETFILE) ); if (!preCompContribPmap -> fileName) error(SYSTEM, "out of memory allocating filename in " "loadPreCompContrib()" ); snprintf(preCompContrib -> waveletFname, sizeof(dirName) + 5, PMAP_CONTRIB_WAVELETFILE, dirName, modName ); return 0; } void loadContribPhotonMap (PhotonMap *pmap, const char *fname) /* Load contribution pmap and its per-modifier precomputed children */ { LUTAB lutInit = LU_SINIT(NULL, freePreCompContribNode); /* Allocate & init LUT containing per-modifier child contrib pmaps */ if (!(pmap -> preCompContribTab = malloc(sizeof(lutInit)))) error(SYSTEM, "out of memory allocating precomputed contribution " "LUT in loadContribPhotonMap()" ); memcpy(pmap -> preCompContribTab, &lutInit, sizeof(lutInit)); /* NOTE: filename already set in loadPhotonMap() pmap -> fileName = (char*)fname; */ /* Init contrib photon map for light source contributions */ initPmapContrib(pmap); /* Load child contribution photon maps for each modifier in contribTab */ lu_doall(pmapContribTab, loadPreCompContrib, pmap); } static int decodeContribs (PreComputedContrib *preCompContrib, DCOLOR *binnedContribs ) /* Decode and decompress mRGBE-encoded wavelet coefficients in preCompContrib -> mrgbeCoeffs, apply inverse wavelet transform and return decoded contributions in binnedContribs. NOTE: The average wavelet coefficient is expected as input in the upper left of the wavelet matrix (i.e. preCompContrib->waveletMatrix [0][0]). Returns 0 on success, else -1. */ { unsigned b, i, coeffBin; mRGBERange *mrgbeRange; mRGBE *mrgbeCoeffs; DCOLOR fCoeff; if (!binnedContribs || !preCompContrib || !binnedContribs || - !preCompContrib -> mrgbeCoeffs || !preCompContrib -> compressedBins + !preCompContrib -> mrgbeCoeffs ) /* Should be initialised by caller! */ return -1; mrgbeCoeffs = preCompContrib -> mrgbeCoeffs; mrgbeRange = &preCompContrib -> mrgbeRange; /* Init mRGBE coefficient normalisation from range */ - mRGBEinit(mrgbeRange, mrgbeRange -> min, mrgbeRange -> max); + if (!mRGBEinit(mrgbeRange, mrgbeRange -> min, mrgbeRange -> max)) + return -1; /* Zero all detail coeffs (omitting average at binnedContribs [0]) */ for (b = 1; b < preCompContrib -> nBins; b++) setcolor(binnedContribs [b], 0, 0, 0); /* Decode mRGBE to wavelet detail coefficients, omitting average at * binnedContrib [0]. * NOTE: The encoded bin preserves the offset for the omitted 1st * coefficient (i.e. the average), so no shifting done on decoding. */ for (b = 0; b < preCompContrib -> nCompressedBins; b++) { coeffBin = mRGBEdecode(mrgbeCoeffs [b], mrgbeRange, fCoeff); #ifdef PMAP_CONTRIB_DBG /* Check for invalid decoded bins */ if (coeffBin < 1 || coeffBin >= preCompContrib -> nBins) error(CONSISTENCY, "bin out of range in decodeContribs()"); #endif copycolor(binnedContribs [coeffBin], fCoeff); } /* Do inverse 2D wavelet transform on preCompContrib -> waveletMatrix (which actually maps to binnedContribs) */ if (waveletInvXform2(preCompContrib -> waveletMatrix, preCompContrib -> tWaveletMatrix, preCompContrib -> l ) < 0 ) error(INTERNAL, "failed inverse wavelet transform in decodeContribs()" ); #ifdef PMAP_CONTRIB_DBG - for (b = 0; b < preCompContrib -> nBins; b++) - /* Sanity check; decoded contribs should be positive! */ - for (i = 0; i < 3; i++) - if (binnedContribs [b][i] < 0) - error(CONSISTENCY, - "negative contributions in getPreCompContrib()" - ); + for (b = 0; b < preCompContrib -> nBins; b++) { + /* Warn if coefficient is negative; this _can_ happen due to loss of + * precision by the mRGBE encoding */ + if (min(min(binnedContribs [b][0], binnedContribs [b][1]), + binnedContribs [b][2] + ) < 0 + ) + error(WARNING, "negative contributions in getPreCompContrib()"); + + /* Dump decoded contribs for debugging */ + printf("%.3g\t", colorAvg(binnedContribs [b])); + } + + putchar('\n'); #endif return 0; } static void getPreCompContribByPhoton (const Photon *photon, OOC_DataIdx photonIdx, PreComputedContrib *preCompContrib, DCOLOR *binnedContribs ) /* Fetch and decode mRGBE-encoded wavelet coefficients for given photon from wavelet file at offset photonIdx, perform inverse wavelet xform, and return the reconstructed binned contributions in binnedContribs (which is assumed to be variable). Returns 0 on success, else -1. */ /* XXX: HERE binnedContribs POINTS TO CACHE PAGE DATA! */ { unsigned b; COLR mrgbeRange32 [2]; COLOR fCoeff; if (!binnedContribs || !preCompContrib) /* Shouldn't happen */ error(INTERNAL, "contributions not initialised in getPreCompContrib()" ); if (preCompContrib -> nBins <= 1) /* No binning, so nothing to decode */ return; /* Lazily initialise preCompContrib */ if (!preCompContrib -> waveletFile) { /* Lazily open wavelet coefficient file */ preCompContrib -> waveletFile = fopen( preCompContrib -> waveletFname, "rb" ); if (!preCompContrib -> waveletFile) { sprintf(errmsg, "failed to open wavelet coefficient file %s " "in getPreCompContrib()", preCompContrib -> waveletFname ); error(SYSTEM, errmsg); } /* Set record size of encoded coeffs in wavelet file */ preCompContrib -> contribSize = PMAP_CONTRIB_ENCSIZE( preCompContrib -> nCompressedBins ); /* Lazily allocate buffer for mRGBE wavelet coeffs */ preCompContrib -> mrgbeCoeffs = calloc( preCompContrib -> nCompressedBins, sizeof(mRGBE) ); - /* Lazily allocate buffer for decoded detail coefficients - (minus avg coeff, which is stored in corresponding photon) */ - preCompContrib -> compressedBins = calloc( - preCompContrib -> nBins - 1, sizeof(DCOLOR) - ); - if (!preCompContrib -> mrgbeCoeffs || - !preCompContrib -> compressedBins - ) + if (!preCompContrib -> mrgbeCoeffs) error(SYSTEM, "out of memory allocating decoded coefficients " "in getPreCompContrib()" ); /* Lazily allocate primary wavelet matrix rows */ preCompContrib -> waveletMatrix = calloc( preCompContrib -> scDim, sizeof(WaveletCoeff3*) ); if (!preCompContrib -> waveletMatrix) error(SYSTEM, "out of memory allocating primary wavelet " "matrix in getPreCompContrib()" ); /* Lazily allocate transposed wavelet matrix */ preCompContrib -> tWaveletMatrix = allocWaveletMatrix(preCompContrib -> l); if (!preCompContrib -> tWaveletMatrix) error(SYSTEM, "out of memory allocating transposed wavelet " " matrix in getPreCompContrib()" ); } /* Seek to photon index in wavelet file and read mRGBE normalisation and encoded bins */ if (fseek(preCompContrib -> waveletFile, photonIdx * preCompContrib -> contribSize, SEEK_SET ) < 0 ) error(SYSTEM, "failed seek in wavelet coefficient file " "in getPreCompContrib()" ); /* Read only mRGBE maximum if 1 compressed bin (maximum compression), * since minimum is implicitly zero and therefore omitted */ getbinary(mrgbeRange32, sizeof(COLR), 1 + (preCompContrib -> nCompressedBins > 1), preCompContrib -> waveletFile ); if (getbinary(preCompContrib -> mrgbeCoeffs, sizeof(mRGBE), preCompContrib -> nCompressedBins, preCompContrib -> waveletFile ) != preCompContrib -> nCompressedBins ) error(SYSTEM, "failed reading from wavelet coefficient file " "in getPreCompContrib()" ); /* Get mRGBE range (NOTE min/max are reversed in the coeff file) */ colr_color(fCoeff, mrgbeRange32 [0]); copycolor(preCompContrib -> mrgbeRange.max, fCoeff); if (preCompContrib -> nCompressedBins > 1) { colr_color(fCoeff, mrgbeRange32 [1]); copycolor(preCompContrib -> mrgbeRange.min, fCoeff); } else { /* Single compressed bin; mRGBE minimum is implicitly zero */ setcolor(preCompContrib -> mrgbeRange.min, 0, 0, 0); } /* Map primary wavelet matrix rows to linear output array binnedContribs; do this every time under the assumption that binnedContribs is not static! */ for (b = 0; b < preCompContrib -> scDim; b++) /* Point to each row in cache linear output array */ preCompContrib -> waveletMatrix [b] = &binnedContribs [b * preCompContrib -> scDim]; /* Set average wavelet coefficient from photon flux; this coefficient is expected by the inverse transform routine in upper left of the wavelet matrix (i.e. preCompContrib -> waveletMatrix [0][0] = binnedContribs [0]) */ getPhotonFlux((Photon*)photon, fCoeff); copycolor(preCompContrib -> waveletMatrix [0][0], fCoeff); /* All set, now decode mRGBE coeffs and invert wavelet transform */ if (decodeContribs(preCompContrib, binnedContribs)) error(INTERNAL, "failed contribution decoding/decompression " "in getPreCompContrib()" ); } typedef struct { const RAY *ray; RREAL rayCoeff [3]; COLORV *totalContrib; } PreCompContribRCData; static int getPreCompContribByMod (const LUENT *preCompContribTabNode, void *p ) /* Fetch and decode precomputed contributions from closest photon in pmap * referenced by preCompContribTabNode, and accumulate in pmapContribTab * for the corresponding modifier */ { PhotonMap *preCompContribPmap = (PhotonMap*)( preCompContribTabNode -> data ); PreComputedContrib *preCompContrib; const char *modName = preCompContribTabNode -> key; MODCONT *modCont; PreCompContribRCData *rcData = (PreCompContribRCData*)p; LUENT *contribTabNode; Photon photon; OOC_DataIdx photonIdx; COLOR photonContrib; unsigned b; if (!preCompContribPmap || !rcData || !(preCompContrib = (PreComputedContrib*)( preCompContribPmap -> preCompContrib )) ) { sprintf(errmsg, "missing precomputed contributions for modifier %s", modName ); error(INTERNAL, errmsg); } /* Find rcontrib's LUT entry for this modifier */ contribTabNode = lu_find(pmapContribTab, modName); if (!contribTabNode || !contribTabNode -> key || !(modCont = (MODCONT*)contribTabNode -> data) ) { sprintf(errmsg, "missing contribution LUT entry for modifier %s", modName ); error(INTERNAL, errmsg); } /* Reallocate rcontrib bins if they don't match precomputed */ if (modCont -> nbins != preCompContrib -> nBins) { contribTabNode -> data = realloc(modCont, sizeof(MODCONT) + sizeof(DCOLOR) * (preCompContrib -> nBins - 1) ); if (!contribTabNode -> data) { sprintf(errmsg, "out of memory reallocating contribution bins " "for modifier %s", modCont -> modname ); error(SYSTEM, errmsg); } else { sprintf(errmsg, "bin count mismatch for modifier %s, resizing", modCont -> modname ); error(WARNING, errmsg); } modCont = (MODCONT*)contribTabNode -> data; modCont -> nbins = preCompContrib -> nBins; memset(modCont -> cbin, 0, sizeof(DCOLOR) * modCont -> nbins); } /* Fetch closest photon and its contributions */ if (find1Photon(preCompContribPmap, rcData -> ray, &photon, &photonIdx )) { /* XXX: TEMPORARY BUFFER FOR DECODED BINS, PENDING CONTRIB CACHE */ DCOLOR binnedContribs [preCompContrib -> nBins]; if (preCompContrib -> nBins > 1) { /* BINNING ENABLED; fetch and decode bins for found photon */ getPreCompContribByPhoton(&photon, photonIdx, preCompContrib, binnedContribs ); } else { /* NO BINNING; get single contribution directly from photon */ getPhotonFlux(&photon, photonContrib); copycolor(binnedContribs [0], photonContrib); } - for (b = 0; b <= preCompContrib -> nBins; b++) { + for (b = 0; b < preCompContrib -> nBins; b++) { /* Accumulate binned contributions into rcontrib's LUT entry for this modifier, weighted by ray coefficient. Also sum total contributions. */ multcolor(binnedContribs [b], rcData -> rayCoeff); addcolor(modCont -> cbin [b], binnedContribs [b]); addcolor(rcData -> totalContrib, binnedContribs [b]); } } return 0; } void getPreCompPhotonContrib (PhotonMap *pmap, RAY *ray, COLOR totalContrib ) /* Get precomputed light source contributions at ray -> rop for all specified modifiers and accumulate them in pmapContribTab (which maps to rcontrib's binning LUT). Also return total precomputed contribs. */ { PreCompContribRCData rcData; /* Ignore sources */ if (ray -> ro && islight(objptr(ray -> ro -> omod) -> otype)) return; /* Get cumulative path coefficient up to photon lookup point */ raycontrib(rcData.rayCoeff, ray, PRIMARY); setcolor(totalContrib, 0, 0, 0); /* Rcontrib results passed to per-modifier child contrib pmaps */ rcData.ray = ray; rcData.totalContrib = totalContrib; /* Iterate over child contribution photon maps for each modifier and * collect their contributions */ lu_doall(pmap -> preCompContribTab, getPreCompContribByMod, &rcData); } #endif /* PMAP_CONTRIB */