diff --git a/oocnn.c b/oocnn.c index 542ffdb..77c17f2 100644 --- a/oocnn.c +++ b/oocnn.c @@ -1,485 +1,492 @@ #ifndef lint static const char RCSid[] = "$Id: oocnn.c,v 1.2 2020/11/10 01:09:47 u-no-hoo Exp u-no-hoo $"; #endif /* ========================================================================= k-nearest neighbour lookup routines for out-of-core octree data structure 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") (c) Tokyo University of Science, supported by the JSPS Grants-in-Aid for Scientific Research (KAKENHI JP19KK0115, "Three-Dimensional Light Flow") ========================================================================= $Id: oocnn.c,v 1.2 2020/11/10 01:09:47 u-no-hoo Exp u-no-hoo $ */ #if !defined(_WIN32) && !defined(_WIN64) && defined(PMAP_OOC) /* No Windoze support for now */ #include "oocnn.h" #include "oocsort.h" #include #include #include /* ========================= SEARCH QUEUE STUFF ========================= */ #ifdef DEBUG_OOC_NN static int OOC_CheckSearchQueue ( OOC_SearchQueue *squeue, const FVECT key, RREAL *(*keyFunc)(const void*), unsigned root ) /* Priority queue sanity check */ { unsigned kid; const OOC_SearchQueueNode *sqn = squeue -> node; void *rec; float d2; if (root < squeue -> tail) { rec = OOC_GetNearest(squeue, sqn [root].idx); d2 = dist2(keyFunc(rec), key); if (sqn [root].dist2 != d2) return -1; if ((kid = (root << 1) + 1) < squeue -> tail) { if (sqn [kid].dist2 > sqn [root].dist2) return -1; else return OOC_CheckSearchQueue(squeue, key, keyFunc, kid); } if (++kid < squeue -> tail) { if (sqn [kid].dist2 > sqn [root].dist2) return -1; else return OOC_CheckSearchQueue(squeue, key, keyFunc, kid); } } return 0; } #endif static float OOC_PutNearest ( - OOC_SearchQueue *squeue, float d2, void *rec, + OOC_SearchQueue *squeue, float d2, void *rec, OOC_DataIdx recIdx, const OOC_SearchAttribCallback *attrib ) /* Insert data record with SQUARED distance to query point into search * priority queue, maintaining the most distant record at the queue head. * If the queue is full, the new record is only inserted if it is closer * than the queue head. * * The search priority queue is represented as a linearised binary tree with * the root corresponding to the queue head, and the tail corresponding to * the last leaf. * * The optional attrib callback maps record attributes to their queue nodes * by calling attrib->findFunc() if attrib != NULL. It enables finding a * previously queued data record and replacing it by rec if the latter's * SQUARED distance d2 is lower. * * The record is copied into the queue's local record buffa for post-search * retrieval to minimise redundant disk access. Note that it suffices to * only rearrange the corresponding indices in the queue nodes when * restoring the priority queue after every insertion, rather than moving * the actual records. * * This function returns the new maximum SQUARED distance at the head if the * search queue is full. Otherwise it returns -1, indicating a maximum for * the entire queue is as yet undefined. */ { void *tRec = NULL; OOC_SearchQueueNode *sqn = squeue -> node, **attribSqn = NULL, **tAttribSqn = NULL; - unsigned root, kid, kid2, rootIdx; + unsigned root, kid, kid2, rootIdx, rootRecIdx; /* Start with undefined max distance */ float d2max = -1; #ifdef PMAP_PATHFILT /* Find previously enqueued record with same attribute if callback defined */ #ifdef PMAP_DEBUGPATHS printf("-------- Enqueue --------\n"); #endif attribSqn = (OOC_SearchQueueNode**)( attrib ? attrib -> findFunc(rec, attrib -> state) : NULL ); #ifdef PMAP_DEBUGPATHS printf(attribSqn ? "found\n" : "not found\n"); #endif #endif if (!attribSqn && squeue -> tail < squeue -> len) { /* No duplicate attribute in queue, and queue not full; insert record at tail and resort towards head */ kid = squeue -> tail++; while (kid) { root = (kid - 1) >> 1; /* Compare with parent and swap if necessary (bubble up), else terminate */ if (d2 > sqn [root].dist2) { - sqn [kid].dist2 = sqn [root].dist2; - sqn [kid].idx = sqn [root].idx; + sqn [kid].dist2 = sqn [root].dist2; + sqn [kid].idx = sqn [root].idx; + sqn [kid].recIdx = sqn [root].recIdx; #ifdef PMAP_PATHFILT if (attrib) { /* Update root's attribute map entry to refer to kid */ tRec = OOC_GetNearest(squeue, sqn [root].idx); tAttribSqn = (OOC_SearchQueueNode**)( attrib -> findFunc(tRec, attrib -> state) ); #ifdef PMAP_DEBUGPATHS printf("update\n"); #endif if (!tAttribSqn) { /* Consistency error: a previously enqueued record must * have an entry in the attribute map! */ fputs("OOC_PutNearest: corrupt attribute map\n", stderr); exit(-1); } *tAttribSqn = &sqn [kid]; } #endif kid = root; } else break; } /* Priority queue restored; assign tail position as linear index into * record buffa squeue -> nnRec and append record */ - sqn [kid].dist2 = d2; - sqn [kid].idx = squeue -> tail - 1; + sqn [kid].dist2 = d2; + sqn [kid].idx = squeue -> tail - 1; + sqn [kid].recIdx = recIdx; memcpy(OOC_GetNearest(squeue, sqn [kid].idx), rec, squeue -> recSize); #ifdef PMAP_PATHFILT if (attrib) /* Add new record to attribute map */ attrib -> addFunc(rec, &sqn [kid], attrib -> state); #endif } else { /* Search queue full or duplicate attribute in queue and new node may be closer; set root to replaced node, otherwise to queue head */ root = attribSqn ? *attribSqn - sqn : 0; if (d2 < sqn [root].dist2) { /* New record closer than maximum at root; replace and resort towards leaves (=search queue tail) */ - rootIdx = sqn [root].idx; + rootIdx = sqn [root].idx; + rootRecIdx = sqn [root].recIdx; #ifdef PMAP_PATHFILT if (attrib && !attribSqn) { /* New record will replace root and has unique attribute; delete * root from the attribute map */ tRec = OOC_GetNearest(squeue, sqn [root].idx); attrib -> delFunc(tRec, attrib -> state); } #endif while ((kid = (root << 1) + 1) < squeue -> tail) { /* Compare with larger kid & swap if necessary (bubble down), else terminate */ if ( (kid2 = (kid + 1)) < squeue -> tail && sqn [kid2].dist2 > sqn [kid].dist2 ) kid = kid2; if (d2 < sqn [kid].dist2) { - sqn [root].dist2 = sqn [kid].dist2; - sqn [root].idx = sqn [kid].idx; + sqn [root].dist2 = sqn [kid].dist2; + sqn [root].idx = sqn [kid].idx; + sqn [root].recIdx = sqn [kid].recIdx; #ifdef PMAP_PATHFILT if (attrib) { /* Update kid's attribute map entry to refer to root */ tRec = OOC_GetNearest(squeue, sqn [kid].idx); tAttribSqn = (OOC_SearchQueueNode**)( attrib -> findFunc(tRec, attrib -> state) ); #ifdef PMAP_DEBUGPATHS printf("update\n"); #endif if (!tAttribSqn) { /* Consistency error: a previously enqueued record must * have an entry in the attribute map! */ fputs("OOC_PutNearest: corrupt attribute map\n", stderr); exit(-1); } *tAttribSqn = &sqn [root]; } #endif } else break; root = kid; } /* Reassign head's / replaced node's previous buffa index and overwrite corresponding record */ sqn [root].dist2 = d2; sqn [root].idx = rootIdx; + sqn [root].recIdx = rootRecIdx; memcpy(OOC_GetNearest(squeue, sqn[root].idx), rec, squeue->recSize); #ifdef PMAP_PATHFILT if (attrib) { if (!attribSqn) /* Add new record to attribute map */ attrib -> addFunc(rec, &sqn [root], attrib -> state); else { #ifdef PMAP_DEBUGPATHS attrib -> findFunc(rec, attrib -> state); printf("update\n"); #endif /* Update new record's existing attribute map entry */ *attribSqn = &sqn [root]; } } #endif } } if (squeue -> tail >= squeue -> len) /* Search queue full; update max distance^2 from head node */ d2max = sqn [0].dist2; #if defined(PMAP_PATHFILT) && defined(PMAP_DEBUGPATHS) if (attrib) /* Run sanity check on attribute map */ attrib -> checkFunc(attrib -> state); #endif return d2max; } int OOC_InitNearest (OOC_SearchQueue *squeue, unsigned len, unsigned recSize) { squeue -> len = len; squeue -> recSize = recSize; squeue -> node = calloc(len, sizeof(OOC_SearchQueueNode)); squeue -> nnRec = calloc(len, recSize); if (!squeue -> node || !squeue -> nnRec) { perror("OOC_InitNearest: failed search queue allocation"); return -1; } return 0; } void *OOC_GetNearest (const OOC_SearchQueue *squeue, unsigned idx) { return squeue -> nnRec + idx * squeue -> recSize; } /* ============================ OCTREE STUFF ============================ */ static float OOC_BBoxDist2 (const FVECT bbOrg, RREAL bbSiz, const FVECT key) /* Return minimum *SQUARED* distance between key and bounding box defined by * bbOrg and bbSiz; a distance of 0 implies the key lies INSIDE the bbox */ { /* Explicit comparison with bbox corners */ int i; FVECT d; for (i = 0; i < 3; i++) { d [i] = key [i] - bbOrg [i]; d [i] = d [i] < 0 ? -d [i] : d [i] - bbSiz; /* Set to 0 if inside bbox */ if (d [i] < 0) d [i] = 0; } return DOT(d, d); } float OOC_FindNearest ( OOC_Octree *oct, OOC_Node *node, OOC_DataIdx dataIdx, const FVECT org, float size, const FVECT key, const OOC_SearchFilterCallback *filter, const OOC_SearchAttribCallback *attrib, OOC_SearchQueue *squeue, float maxDist2 ) { const float kidSize = size * 0.5; unsigned i, kid, kid0; float d2; char rec [oct -> recSize]; FVECT kidOrg; OOC_DataIdx kidDataIdx, recIdx; OOC_Node *kidNode; if (!node) /* No current node -- shouldn't happen */ return -1; /* Start with suboctant closest to key */ for (kid0 = 0, i = 0; i < 3; i++) kid0 |= (key [i] > org [i] + kidSize) << i; for (i = 0; i < 8; i++) { kid = kid0 ^ i; kidNode = node; kidDataIdx = dataIdx + OOC_GetKid(oct, &kidNode, kid); /* Prune empty suboctant */ if ((!kidNode && !OOC_ISLEAF(node)) || (OOC_ISLEAF(node) && !node -> leaf.num [kid]) ) continue; /* Set up suboctant */ VCOPY(kidOrg, org); OOC_OCTORIGIN(kidOrg, kid, kidSize); /* Prune suboctant if not overlapped by maxDist2 */ if (OOC_BBoxDist2(kidOrg, kidSize, key) > maxDist2) continue; if (kidNode) { /* Internal node; recurse into non-empty suboctant */ maxDist2 = OOC_FindNearest(oct, kidNode, kidDataIdx, kidOrg, kidSize, key, filter, attrib, squeue, maxDist2 ); if (maxDist2 < 0) /* Bail out on error */ break; } else if (OOC_ISLEAF(node)) { /* Leaf node; do linear check of all records in suboctant */ for (recIdx = kidDataIdx; recIdx < kidDataIdx + node -> leaf.num [kid]; recIdx++ ) { if (OOC_GetData(oct, recIdx, rec)) return -1; if (!filter || filter -> filtFunc(rec, filter -> state)) { /* Insert record in search queue if it passes filter * and its SQUARED dist to key is below maxDist2 */ if ((d2 = dist2(key, oct -> key(rec))) < maxDist2) { - if ((d2 = OOC_PutNearest(squeue, d2, rec, attrib)) >= 0) + if ((d2 = OOC_PutNearest(squeue, d2, rec, recIdx, attrib) + ) >= 0 + ) /* Update maxDist2 if search queue is full */ maxDist2 = d2; #ifdef DEBUG_OOC_NN if (OOC_CheckSearchQueue(squeue, key, oct->key, 0)) { fprintf(stderr, "OOC_PutNearest: corrupt search queue\n" ); return -1; } #endif } } } } } /* XXX: Note this will not correspond to the maximum distance in a partially filled search queue, since it hasn't been finalised; this distance should therefore be obtained from the search queue head by the caller when this func returns! */ return maxDist2; } float OOC_Find1Nearest ( OOC_Octree *oct, OOC_Node *node, OOC_DataIdx dataIdx, const FVECT org, float size, const FVECT key, const OOC_SearchFilterCallback *filter, void *nnRec, OOC_DataIdx *nnIdx, float maxDist2 ) { const float kidSize = size * 0.5; unsigned i, kid, kid0; float d2; char rec [oct -> recSize]; FVECT kidOrg; OOC_DataIdx kidDataIdx, recIdx; OOC_Node *kidNode; if (!node) /* No current node -- shouldn't happen */ return -1; /* Start with suboctant closest to key */ for (kid0 = 0, i = 0; i < 3; i++) kid0 |= (key [i] > org [i] + kidSize) << i; for (i = 0; i < 8; i++) { kid = kid0 ^ i; kidNode = node; kidDataIdx = dataIdx + OOC_GetKid(oct, &kidNode, kid); /* Prune empty suboctant */ if ((!kidNode && !OOC_ISLEAF(node)) || (OOC_ISLEAF(node) && !node -> leaf.num [kid]) ) continue; /* Set up suboctant */ VCOPY(kidOrg, org); OOC_OCTORIGIN(kidOrg, kid, kidSize); /* Prune suboctant if not overlapped by maxDist2 */ if (OOC_BBoxDist2(kidOrg, kidSize, key) > maxDist2) continue; if (kidNode) { /* Internal node; recurse into non-empty suboctant */ maxDist2 = OOC_Find1Nearest(oct, kidNode, kidDataIdx, kidOrg, kidSize, key, filter, nnRec, nnIdx, maxDist2 ); if (maxDist2 < 0) /* Bail out on error */ break; } else if (OOC_ISLEAF(node)) /* Leaf node; do linear check of all records in suboctant */ for (recIdx = kidDataIdx; recIdx < kidDataIdx + node -> leaf.num [kid]; recIdx++ ) { if (OOC_GetData(oct, recIdx, rec)) return -1; if (!filter || filter -> filtFunc(rec, filter -> state)) /* If record passes filter, update closest record, its index, and max SQUARED distance to its key */ if ((d2 = dist2(key, oct -> key(rec))) < maxDist2) { memcpy(nnRec, rec, oct -> recSize); maxDist2 = d2; if (nnIdx) *nnIdx = recIdx; } } } return maxDist2; } #endif /* NIX / PMAP_OOC */ diff --git a/oocnn.h b/oocnn.h index f747d3b..37bb3af 100644 --- a/oocnn.h +++ b/oocnn.h @@ -1,125 +1,126 @@ /* ========================================================================= k-nearest neighbour lookup routines for out-of-core octree data structure 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") (c) Tokyo University of Science, supported by the JSPS Grants-in-Aid for Scientific Research (KAKENHI JP19KK0115, "Three-Dimensional Light Flow") ========================================================================= $Id: oocnn.h,v 1.2 2020/11/10 01:15:09 u-no-hoo Exp u-no-hoo $ */ #ifndef OOC_SEARCH_H #define OOC_SEARCH_H #include "oococt.h" /* Priority queue node for octree lookups */ typedef struct { float dist2; /* Record's distance^2 to query point */ unsigned idx; /* Record's index in queue buffer */ + OOC_DataIdx recIdx; /* Record's absolute index */ } OOC_SearchQueueNode; /* Priority queue for octree lookups. Note that we explicitly store the * NN candidates in a local in-core buffer nnRec for later retrieval * without incurring additional disk I/O */ typedef struct { OOC_SearchQueueNode *node; unsigned len, tail, recSize; void *nnRec; } OOC_SearchQueue; /* Requires above typedefs */ #include "pmapfilt.h" float OOC_FindNearest ( OOC_Octree *oct, OOC_Node *node, OOC_DataIdx dataIdx, const FVECT org, float size, const FVECT key, const OOC_SearchFilterCallback *filter, const OOC_SearchAttribCallback *attrib, OOC_SearchQueue *queue, float maxDist2 ); /* Find k nearest neighbours (where k = queue -> len) within a maximum * SQUARED distance of maxDist2 around key. Returns the corresponding * data records with their SQUARED distances in the search queue * (queue -> node [0] .. queue -> node [queue -> tail - 1]), with the * furthest neighbour at the queue head and queue->tail <= queue->len. * * This is a recursive function requiring params for the current node: * DataIdx is the data offset for records in the current node, which is * necessary for implied addressing into the leaf file. Org and size are * the origin and size of the current node's bounding box. * * This function accepts two optional callbacks for selective search: * /SHALL MODE ON * * -- If filter != NULL, filter->filtFunc() is called for each candidate * nearest neighbour along with the encapsulated state data provided * by the caller; records SHALL be accepted if filter->filtFunc() * returns nonzero. * * -- If attrib != NULL, attrib->findFunc() is called for each candidate * nearest neighbour that has been accepted by filter->filtFunc() and * will therefore be added to the search queue. * This function SHALL map a record attribute to its search queue node * via a lookup table ("attribute map") accessible via attrib->state, * and SHALL return a modifiable pointer to the corresponding search * queue node, or NULL if the attribute is absent. * This callback SHALL support replacing a previously enqueued record * with the candidate if they share the same attribute and the latter * is closer to the key; this avoids accumulating records with * duplicate attributes in the queue, which may be desirable in some * applications. Or not. Maybe not. * This callback SHALL also support updating the attribute map entry * for each resorted record when restoring the priority queue property * during insertions; the search queue nodes referenced in the * attribute map are directly modified by dereferencing the returned * pointers. * In addition, the callbacks attrib->addFunc() and attrib->delFunc() * SHALL insert new entries into the attribute map, or delete evicted * ones once the search queue is full. * * /SHALL MODE OFF * We use the word "SHALL" a lot here 'cos Lars likes it and it sounds * very imposing and professional and stuff. * * This function returns -1 on failure, else a positive value on success. * The maximum (squared) distance to the found neighours can be obtained * from the head of the search queue, queue -> node [0].dist2. * (Note this will only coincide with the return value if the queue is * full, so this value should not be relied on!). */ float OOC_Find1Nearest ( OOC_Octree *oct, OOC_Node *node, OOC_DataIdx dataIdx, const FVECT org, float size, const FVECT key, const OOC_SearchFilterCallback *filter, void *nnRec, OOC_DataIdx *nnIdx, float maxDist2 ); /* Find single nearest neighbour within max SQUARED distance maxDist2 * around key and copy corresponding record in nnRec, and return its * index in nnIdx (if nnIdx != NULL). This is an optimised version of * OOC_FindNearest() without a search queue */ int OOC_InitNearest ( OOC_SearchQueue *squeue, unsigned len, unsigned recSize ); /* Initialise NN search queue of length len and local buffa for records * of size recSize; precondition to calling OOC_FindNearest() */ void *OOC_GetNearest (const OOC_SearchQueue *queue, unsigned idx); /* Return pointer to record at pos idx in search queue after calling * OOC_FindNearest() */ #endif diff --git a/oococt.c b/oococt.c index 2348458..7d5f6f4 100644 --- a/oococt.c +++ b/oococt.c @@ -1,484 +1,484 @@ #ifndef lint static const char RCSid[] = "$Id: oococt.c,v 2.4 2017/08/14 21:12:10 rschregle Exp $"; #endif /* ====================================================================== Out-of-core octree data structure 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) ====================================================================== $Id: oococt.c,v 2.4 2017/08/14 21:12:10 rschregle Exp $ */ #if !defined(_WIN32) && !defined(_WIN64) || defined(PMAP_OOC) /* No Windoze support for now */ #include "oococt.h" #include "rtio.h" #include #include void OOC_Null (OOC_Octree *oct) /* Init empty octree */ { oct -> maxNodes = oct -> numNodes = oct -> leafMax = oct -> maxDepth = oct -> numData = oct -> size = oct -> recSize = oct -> mortonScale = 0; oct -> org [0] = oct -> org [1] = oct -> org [2] = oct -> bound [0] = oct -> bound [1] = oct -> bound [2] = 0; oct -> key = NULL; oct -> nodes = NULL; oct -> leafFile = NULL; oct -> leafFname = NULL; oct -> cache = NULL; } void OOC_Init (OOC_Octree *oct, unsigned recSize, FVECT org, RREAL size, RREAL *(*key)(const void*), FILE *leafFile, char *leafFname ) { oct -> maxNodes = oct -> numNodes = oct -> leafMax = oct -> maxDepth = oct -> numData = 0; VCOPY(oct -> org, org); oct -> size = oct -> bound[0] = oct -> bound[1] = oct -> bound[2] = size; VADD(oct -> bound, oct -> bound, org); oct -> mortonScale = size > 0 ? MORTON3D_MAX / size : 0; oct -> recSize = recSize; oct -> key = key; oct -> nodes = NULL; oct -> leafFile = leafFile; oct -> leafFname = leafFname; oct -> cache = NULL; /* Cache currently initialised externally */ } OOC_Node *NewNode (OOC_Octree *oct) /* Allocate new octree node, enlarge array if necessary. Return pointer to new node or NULL if failed. */ { OOC_Node *n = NULL; if (oct -> numNodes >= OOC_NODEIDX_MAX) { /* Node index overflow */ fprintf(stderr, "OOC_NewNode: node index overflow (numNodes = %d)\n", oct -> numNodes ); return NULL; } if (++oct -> numNodes > oct -> maxNodes) { /* Need to allocate root or enlarge array */ oct -> maxNodes += OOC_BLKSIZE / sizeof(OOC_Node); oct -> nodes = realloc(oct -> nodes, oct -> maxNodes * sizeof(OOC_Node) ); if (!oct -> nodes) { perror("OOC_NewNode: couldn't realloc() nodes"); return NULL; } } n = oct -> nodes + oct -> numNodes - 1; n -> node.num = n -> node.kid = n -> node.branch = n -> node.type = 0; return n; } int OOC_GetData (OOC_Octree *oct, OOC_DataIdx idx, void *rec) /* Reads indexed data record from leaf file and copies it to rec, else * returns -1 on failure */ { if (!oct || !rec) return -1; if (idx >= oct -> numData) { fprintf(stderr, "OOC_GetData: invalid data record index %u\n", idx); return -1; } if (oct -> cache) { /* Retrieve record from leaf file via I/O cache */ void *cachedRec; if (!(cachedRec = OOC_CacheData(oct -> cache, oct -> leafFile, idx))) return -1; /* It's safer to copy the record to the caller's local buffer as a * returned pointer may be invalidated by subsequent calls if the * page it points to is swapped out */ memcpy(rec, cachedRec, oct -> recSize); } else { /* No I/O caching; do (SLOW!) per-record I/O from leaf file */ const unsigned long pos = (unsigned long)oct -> recSize * idx; if (pread(fileno(oct -> leafFile), rec, oct -> recSize, pos) != oct -> recSize) { perror("OOC_GetData: failed seek/read in leaf file"); return -1; } } return 0; } int OOC_InBBox (const OOC_Octree *oct, const FVECT org, RREAL size, const FVECT key ) /* Check whether key lies within bounding box defined by org and size. * Compares Morton codes rather than the coordinates directly to account * for dicretisation error introduced by the former. */ { const MortonIdx keyZ = OOC_KEY2MORTON(key, oct); FVECT bbox = {size, size, size}; VADD(bbox, org, bbox); return keyZ > OOC_KEY2MORTON(org, oct) && keyZ < OOC_KEY2MORTON(bbox, oct); } int OOC_Branch (const OOC_Octree *oct, const FVECT org, RREAL size, const FVECT key, FVECT nuOrg ) /* Return index of octant containing key and corresponding origin if nuOrg * != NULL, or -1 if key lies outside all octants. Size is already assumed * to be halved, i.e. corresponding to the length of an octant axis. * Compares Morton codes rather than the coordinates directly to account * for dicretisation error introduced by the former. */ { int o; FVECT octOrg; for (o = 0; o < 8; o++) { VCOPY(octOrg, org); OOC_OCTORIGIN(octOrg, o, size); if (OOC_InBBox(oct, octOrg, size, key)) { if (nuOrg) VCOPY(nuOrg, octOrg); return o; } } return -1; } int OOC_BranchZ (const OOC_Octree *oct, const FVECT org, RREAL size, MortonIdx keyZ, FVECT nuOrg ) /* Optimised version of OOC_Branch() with precomputed key Morton code */ { int o; const FVECT cent = {size, size, size}; FVECT octOrg, bbox; for (o = 0; o < 8; o++) { VCOPY(octOrg, org); OOC_OCTORIGIN(octOrg, o, size); VADD(bbox, octOrg, cent); if (keyZ > OOC_KEY2MORTON(octOrg, oct) && keyZ < OOC_KEY2MORTON(bbox, oct) ) { if (nuOrg) VCOPY(nuOrg, octOrg); return o; } } return -1; } OOC_DataIdx OOC_GetKid (const OOC_Octree *oct, OOC_Node **node, unsigned kid) /* For a given octree node, return the sum of the data counters for kids * [0..k-1]. On return, the node either points to the k'th kid, or * NULL if the kid is nonexistent or the node is a leaf */ { unsigned k; OOC_Node *kidNode = OOC_KID1(oct, *node); /* NULL if leaf */ OOC_DataIdx dataIdx = 0; for (k = 0; k < kid; k++) { if (OOC_ISLEAF(*node)) /* Sum data counters of leaf octants */ dataIdx += (*node) -> leaf.num [k]; else /* Sum data counter of inner node's nonempty kids and advance * pointer to sibling */ if (OOC_ISBRANCH(*node, k)) { dataIdx += OOC_DATACNT(kidNode); kidNode++; } } /* Update node pointer to selected kid (or NULL if nonexistent or node is * a leaf) */ *node = OOC_ISBRANCH(*node, kid) ? kidNode : NULL; return dataIdx; } #ifdef DEBUG_OOC int OOC_Check (OOC_Octree *oct, const OOC_Node *node, FVECT org, RREAL size, OOC_DataIdx dataIdx ) /* Traverse tree & check for valid nodes; oct -> leafFile must be open; * return 0 if ok, otherwise -1 */ { unsigned k; FVECT kidOrg; const RREAL kidSize = size * 0.5; if (!oct || !node) return -1; if (OOC_ISLEAF(node)) { /* Node is a leaf; check records in each octant */ char rec [oct -> recSize]; /* Violates C89! */ OOC_OctCnt d; RREAL *key; for (k = 0; k < 8; k++) { VCOPY(kidOrg, org); OOC_OCTORIGIN(kidOrg, k, kidSize); for (d = node -> leaf.num [k]; d; d--, dataIdx++) { #ifdef DEBUG_OOC_CHECK fprintf(stderr, "OOC_Check: checking record %lu\n", (unsigned long)dataIdx ); #endif if (OOC_GetData(oct, dataIdx, rec)) { fprintf(stderr, "OOC_Check: failed getting record at %lu\n", (unsigned long)dataIdx ); return -1; } key = oct -> key(rec); if (!OOC_InBBox(oct, kidOrg, kidSize, key)) { fprintf(stderr, "OOC_Check: key [%f, %f, %f] (zIdx %020lu) " - "outside bbox [[%f-%f], [%f-%f], [%f-%f]] " + "outside bbox [ [%f, %f], [%f, %f], [%f, %f] ] " "in octant %d (should be %d)\n", key [0], key [1], key [2], OOC_KEY2MORTON(key, oct), kidOrg [0], kidOrg [0] + kidSize, kidOrg [1], kidOrg [1] + kidSize, kidOrg [2], kidOrg [2] + kidSize, k, OOC_Branch(oct, org, kidSize, key, NULL) ); /* return -1; */ } } } } else { /* Internal node; recurse on all kids */ const OOC_Node *kid = OOC_KID1(oct, node); OOC_DataIdx numData = 0; if (node -> node.kid >= oct -> numNodes) { fprintf(stderr, "OOC_Check: invalid node index %u\n", node -> node.kid ); return -1; } if (!node -> node.num) { fputs("OOC_Check: empty octree node\n", stderr); return -1; } for (k = 0; k < 8; k++) if (OOC_ISBRANCH(node, k)) { VCOPY(kidOrg, org); OOC_OCTORIGIN(kidOrg, k, kidSize); if (OOC_Check(oct, kid, kidOrg, kidSize, dataIdx)) return -1; dataIdx += OOC_DATACNT(kid); numData += OOC_DATACNT(kid); kid++; } if (node -> node.num != numData) { fprintf(stderr, "Parent/kid data count mismatch; expected %u, found %u\n", node -> node.num, numData ); return -1; } } return 0; } #endif int OOC_SaveOctree (const OOC_Octree *oct, FILE *out) /* Appends octree nodes to specified file along with metadata. Uses * RADIANCE's portable I/O routines. Returns 0 on success, else -1. Note * the internal representation of the nodes is platform specific as they * contain unions, therefore we use the portable I/O routines from the * RADIANCE library */ { int i; OOC_NodeIdx nc; OOC_Node *np = NULL; if (!oct || !oct->nodes || !oct->numData || !oct->numNodes || !out) { fputs("OOC_SaveOctree: empty octree\n", stderr); return -1; } /* Write octree origin, size, data record size, number of records, and * number of nodes */ for (i = 0; i < 3; i++) putflt(oct -> org [i], out); putflt(oct -> size, out); putint(oct -> recSize, sizeof(oct -> recSize), out); putint(oct -> numData, sizeof(oct -> numData), out); putint(oct -> numNodes, sizeof(oct -> numNodes), out); /* Write nodes to file */ for (nc = oct -> numNodes, np = oct -> nodes; nc; nc--, np++) { if (OOC_ISLEAF(np)) { /* Write leaf marker */ putint(-1, 1, out); /* Write leaf octant counters */ for (i = 0; i < 8; i++) putint(np -> leaf.num [i], sizeof(np -> leaf.num [i]), out); } else { /* Write node marker */ putint(0, 1, out); /* Write node, rounding field size up to nearest number of bytes */ putint(np -> node.kid, (OOC_NODEIDX_BITS + 7) >> 3, out); putint(np -> node.num, (OOC_DATAIDX_BITS + 7) >> 3, out); putint(np -> node.branch, 1, out); } if (ferror(out)) { fputs("OOC_SaveOctree: failed writing to node file", stderr); return -1; } } fflush(out); return 0; } int OOC_LoadOctree (OOC_Octree *oct, FILE *nodeFile, RREAL *(*key)(const void*), FILE *leafFile ) /* Load octree nodes and metadata from nodeFile and associate with records * in leafFile. Uses RADIANCE's portable I/O routines. Returns 0 on * success, else -1. */ { int i; OOC_NodeIdx nc; OOC_Node *np = NULL; if (!oct || !nodeFile) return -1; OOC_Null(oct); /* Read octree origin, size, record size, #records and #nodes */ for (i = 0; i < 3; i++) oct -> org [i] = getflt(nodeFile); oct -> size = getflt(nodeFile); oct -> bound [0] = oct -> bound [1] = oct -> bound [2] = oct -> size; VADD(oct -> bound, oct -> bound, oct -> org); oct -> mortonScale = MORTON3D_MAX / oct -> size; oct -> recSize = getint(sizeof(oct -> recSize), nodeFile); oct -> numData = getint(sizeof(oct -> numData), nodeFile); oct -> numNodes = getint(sizeof(oct -> numNodes), nodeFile); oct -> key = key; oct -> leafFile = leafFile; if (feof(nodeFile) || !oct -> numData || !oct -> numNodes) { fputs("OOC_LoadOctree: empty octree\n", stderr); return -1; } /* Allocate octree nodes */ if (!(oct -> nodes = calloc(oct -> numNodes, sizeof(OOC_Node)))) { fputs("OOC_LoadOctree: failed octree allocation\n", stderr); return -1; } /* Read nodes from file */ for (nc = oct -> numNodes, np = oct -> nodes; nc; nc--, np++) { OOC_CLEARNODE(np); /* Read node type marker */ if (getint(1, nodeFile)) { /* Read leaf octant counters and set node type */ for (i = 0; i < 8; i++) np -> leaf.num [i] = getint(sizeof(np -> leaf.num [i]), nodeFile ); OOC_SETLEAF(np); } else { /* Read node, rounding field size up to nearest number of bytes */ np -> node.kid = getint((OOC_NODEIDX_BITS + 7) >> 3, nodeFile); np -> node.num = getint((OOC_DATAIDX_BITS + 7) >> 3, nodeFile); np -> node.branch = getint(1, nodeFile); } if (ferror(nodeFile)) { fputs("OOC_LoadOctree: failed reading from node file\n", stderr); return -1; } } return 0; } void OOC_Delete (OOC_Octree *oct) /* Self-destruct */ { free(oct -> nodes); if (oct -> leafFile) fclose(oct -> leafFile); if (oct -> cache) OOC_DeleteCache(oct -> cache); } #endif /* NIX / PMAP_OOC */ diff --git a/pmapcontrib.c b/pmapcontrib.c index 8b031c8..b6fb66a 100644 --- a/pmapcontrib.c +++ b/pmapcontrib.c @@ -1,1388 +1,1498 @@ #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 $ */ #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; #ifdef PMAP_CONTRIB_DBG /* Sum/num of mRGBE encoding errors for wavelet coeffs */ static WAVELET_COEFF mrgbeDiffs = 0; static unsigned long mrgbeDiffCnt = 0; #endif extern void SDdisk2square(double sq[2], double diskx, double disky); 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 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 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; /* 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; /* 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 */ { 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_PADD4_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; 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++) for (j = 0; j < coeffDim; j++) if ((i >= WAVELET_PADD4_APPROXDIM || j >= WAVELET_PADD4_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 = PMAP_CONTRIB_XY2LIN(coeffDim, i, j); threshCoeffPtr++ -> coeff = (WAVELET_COEFF*)&( waveletMatrix [i] [j] ); } /* Num of nonzero detail coeffs in buffer, number actually expected */ numThresh = threshCoeffPtr - threshCoeffs; nNzDetailCoeffs = WAVELET_PADD4_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. */ for (i = j = coeffDim - 1; numThresh < nNzDetailCoeffs; numThresh++) { threshCoeffPtr -> idx = PMAP_CONTRIB_XY2LIN(coeffDim, i, j); threshCoeffPtr++ -> coeff = (WAVELET_COEFF*)(waveletMatrix [i][j]); } /* 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_PADD4_APPROXDIM-1] [0..WAVELET_PADD4_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++) { #ifdef PMAP_CONTRIB_DBG #if 0 /* Replace wavelet coeff with its linear index for debugging */ threshCoeffs [i].coeff [k] = threshCoeffs [i].idx = i + 1; #endif #endif 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); #if 0 if (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); } #else mrgbeDiffs += sqrt(dist2(decCoeff, threshCoeffs [i].coeff)); mrgbeDiffCnt++; #endif for (k = 0; k < 3; k++) { absCoeff = fabs(threshCoeffs [i].coeff [k]) > 1e-6 ? decCoeff [k] / threshCoeffs [i].coeff [k] : 0; #if 0 if (absCoeff > 1e-6 && fabs(absCoeff - 1) > 1) { sprintf(errmsg, "mRGBE encoding out of tolerance in encodeContribs()\n" "Encoded: [%.3g %.3g %.3g]\nDecoded: [%.3g %.3g %.3g]", threshCoeffs [i].coeff [0], threshCoeffs [i].coeff [1], threshCoeffs [i].coeff [2], decCoeff [0], decCoeff [1], decCoeff [2] ); error(CONSISTENCY, errmsg); } #endif if (absCoeff < 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); */ } } +#if 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; - + unsigned i, bin; + float maxBinDist2 [contrib -> nbins]; + /* 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 */ + + /* 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 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; + #if 0 + r2 = max(sqn -> dist2, (sqn + 1) -> dist2); + r2 = 0.25 * (pmap -> maxDist2 + r2 + 2 * sqrt(pmap -> maxDist2 * r2)); + #endif + + /* 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); + norm = 2 / (PI * PI); + #else + /* Normalise accumulated flux by search area PI * r^2, including + RADIANCE-specific lambertian factor PI */ +// norm = 1 / (PI * PI * r2); + norm = 1 / (PI * PI); + #endif + + memset(maxBinDist2, 0, contrib -> nbins * sizeof(float)); + for (i = 0, sqn = pmap -> squeue.node; + i < pmap -> squeue.tail; + i++, sqn++ + ) { + photon = getNearestPhoton(&pmap -> squeue, sqn -> idx); + bin = photonSrcBin(pmap, photon); + + if (sqn -> dist2 > maxBinDist2 [bin]) + maxBinDist2 [bin] = sqn -> dist2; + } + + /* Skip the extra photon */ + for (i = 0, sqn = pmap -> squeue.node; + i < pmap -> squeue.tail; + i++, sqn++ + ) { + /* Get photon's contribution to density estimate */ + photon = getNearestPhoton(&pmap -> squeue, sqn -> idx); + getPhotonFlux(photon, flux); + bin = photonSrcBin(pmap, photon); + r2 = maxBinDist2 [bin]; + if (r2 > FTINY) { + scalecolor(flux, norm / r2); + #ifdef PMAP_EPANECHNIKOV + /* Apply Epanechnikov kernel to photon flux based on distance */ + scalecolor(flux, 1 - sqn -> dist2 / r2); + #endif + addcolor(irrad, flux); + addcolor(contrib -> cbin [bin], flux); + if (!isnormal(contrib -> cbin [bin][0])) + error(CONSISTENCY, "Inf/NaN contribs"); + } + } + + return contrib; + } + +#else + 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, bin; + float maxBinDist2 [contrib -> nbins]; + + /* 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; 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++) { + for (i = 0, sqn = pmap -> squeue.node; + i < pmap -> squeue.tail; + i++, sqn++ + ) { /* Get photon's contribution to density estimate */ photon = getNearestPhoton(&pmap -> squeue, sqn -> idx); getPhotonFlux(photon, flux); + bin = photonSrcBin(pmap, photon); 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); - addcolor(contrib -> cbin [photonSrcBin(pmap, photon)], flux); + addcolor(contrib -> cbin [bin], flux); } return contrib; } - +#endif + 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); } } /* 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; } 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 p, numPreComp; unsigned i, j, k, coeffIdx, scDim, nBins, coeffDim, nCoeffs, nCompressedCoeffs, nNzDetailCoeffs; PhotonIdx pIdx; Photon photon; RAY ray; MODCONT *binnedContribs; COLR rgbe32 [WAVELET_PADD4_NUMAPPROX + 2]; LUENT *preCompContribNode; PhotonMap *preCompContribPmap, nuPmap; PreComputedContrib *preCompContrib; LUTAB lutInit = LU_SINIT(NULL, freePreCompContribNode); WaveletMatrix2 waveletMatrix; #ifdef PMAP_CONTRIB_DBG RAY dbgRay; #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; if (verbose) { sprintf(errmsg, "\nPrecomputing contributions for %ld photons\n", numPreComp ); eputs(errmsg); #if NIX fflush(stderr); #endif } photonRay(NULL, &ray, PRIMARY, NULL); ray.ro = NULL; for (p = 0; p < numPreComp; p++) { /* Get random photon from stratified distribution in source heap * to avoid duplicates and clustering */ pIdx = firstPhoton(pmap) + (unsigned long)( (p + 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 (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 */ binnedContribs = getPhotonContrib(pmap, &ray, ray.rcol); #ifdef PMAP_CONTRIB_DBG #if 0 /* Get all contribs from same photon for debugging */ /* Position must differ otherwise too many identical photon keys * will cause ooC octree to overflow! */ VCOPY(dbgRay.rop, photon.pos); getPhoton(pmap, 0, &photon); memcpy(&dbgRay, &ray, sizeof(RAY)); dbgRay.rsrc = photonSrcIdx(pmap, &photon); for (i = 0; i < 3; i++) dbgRay.ron [i] = photon.norm [i] / 127.0; binnedContribs = getPhotonContrib(pmap, &dbgRay, dbgRay.rcol); #endif #endif if (binnedContribs) { #if 0 if (!(p & PMAP_CNTUPDATE)) printf("Precomputing contribution photon %lu / %lu\n", p, numPreComp ); #endif /* 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; /* Set Shirley-Chiu square & wavelet matrix dimensions (number of bins resp. coefficients). */ preCompContrib -> scDim = scDim = sqrt( preCompContrib -> nBins = nBins = binnedContribs -> nbins ); if (scDim * scDim != nBins) /* Mkpmap shoulda took care of dis */ error(INTERNAL, "contribution bin count not " "integer square in preComputeContrib()" ); if (nBins > 1) { /* BINNING ENABLED; set up wavelet xform & compression. Get dimensions of wavelet coefficient matrix after boundary padding, and number of nonzero coefficients. The number of compressed (detail) coefficients is fixed and based on the number of NONZERO coefficients (minus approximations, see NOTE below) since zero coeffs are considered thresholded by default. !!! NOTE: THE APPROXIMATION COEFFICIENTS IN THE UPPER !!! LEFT OF THE MATRIX ARE _NEVER_ THRESHOLDED TO !!! MINIMISE COMPRESSION ARTEFACTS! THEY ARE ESSENTIAL !!! FOR RECONSTRUCTING THE ORIGINAL CONTRIBUTIONS. */ preCompContrib -> coeffDim = coeffDim = padWaveletXform2( NULL, NULL, scDim, &preCompContrib -> nNonZeroCoeffs ); preCompContrib -> nCoeffs = nCoeffs = coeffDim * coeffDim; nNzDetailCoeffs = WAVELET_PADD4_NUMDETAIL( preCompContrib -> nNonZeroCoeffs ); nCompressedCoeffs = (1 - compressRatio) * nNzDetailCoeffs; if (nCompressedCoeffs < 1) { error(WARNING, "maximum contribution compression, clamping number " "of wavelet coefficients in preComputeContrib()" ); nCompressedCoeffs = 1; } if (nCompressedCoeffs >= preCompContrib -> nCoeffs) { error(WARNING, "minimum contribution compression, clamping number " "of wavelet coefficients in preComputeContrib()" ); nCompressedCoeffs = nNzDetailCoeffs; } preCompContrib -> nCompressedCoeffs = nCompressedCoeffs; /* Lazily allocate primary and transposed wavelet * coefficient matrix */ preCompContrib -> waveletMatrix = waveletMatrix = allocWaveletMatrix2(coeffDim); preCompContrib -> tWaveletMatrix = allocWaveletMatrix2(coeffDim); if (!waveletMatrix || !preCompContrib -> tWaveletMatrix) error(SYSTEM, "out of memory allocating wavelet " "coefficient matrix in preComputeContrib()" ); /* 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 preComputeContrib()" ); /* Set size of compressed contributions, in bytes */ preCompContrib -> contribSize = PMAP_CONTRIB_ENCSIZE( nCompressedCoeffs ); } } else { /* LUT node already exists for this modifier; use its data */ 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; copy binned contribs row-by-row to * wavelet xform input matrix, optionally applying logarithm, * then compress & encode */ for (i = 0; i < scDim; i++) { memcpy(waveletMatrix [i], binnedContribs->cbin + PMAP_CONTRIB_XY2LIN(scDim, i, 0), scDim * WAVELET_COEFFSIZE ); #ifdef PMAP_CONTRIB_LOG for (j = 0; j < scDim; j++) for (k = 0; k < 3; k++) waveletMatrix [i] [j] [k] = PMAP_CONTRIB_VAL( waveletMatrix [i] [j] [k] ); #endif } if (encodeContribs(preCompContrib, compressRatio)) error(INTERNAL, "failed contribution " "compression/encoding in preComputeContrib()" ); /* Encode wavelet approx coeffs in the upper left of the * wavelet matrix to 32-bit RGBE. These are not thresholded * to minimise compression artefacts. */ for (i = 0; i < WAVELET_PADD4_APPROXDIM; i++) for (j = 0; j < WAVELET_PADD4_APPROXDIM; j++) { coeffIdx = PMAP_CONTRIB_XY2LIN( WAVELET_PADD4_APPROXDIM, i, j ); setcolr(rgbe32 [coeffIdx], fabs(waveletMatrix [i] [j] [0]), fabs(waveletMatrix [i] [j] [1]), fabs(waveletMatrix [i] [j] [2]) ); /* HACK: depending on the boundary extension mode, some * approx coeff 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 [coeffIdx] [k] = PMAP_CONTRIB_SET_RGBE32_SGN( rgbe32 [coeffIdx] [k], waveletMatrix [i] [j] [k] ); } /* Encode wavelet detail coeff range to 32-bit RGBE */ setcolr(rgbe32 [WAVELET_PADD4_NUMAPPROX], preCompContrib -> mrgbeRange.max [0], preCompContrib -> mrgbeRange.max [1], preCompContrib -> mrgbeRange.max [2] ); setcolr(rgbe32 [WAVELET_PADD4_NUMAPPROX + 1], preCompContrib -> mrgbeRange.min [0], preCompContrib -> mrgbeRange.min [1], preCompContrib -> mrgbeRange.min [2] ); /* Dump 32-bit RGBE approx coeffs followed by mRGBE range to * contribution heap file. NOTE: mRGBE range minimum and * maximum are reversed in the file to facilitate handling * the special (but pointless) case of a single compressed * coeff; if so, only the mRGBE maximum is dumped, since the * minimum is implicitly zero and can be omitted to save * space */ if (putbinary(rgbe32, sizeof(COLR), WAVELET_PADD4_NUMAPPROX + 1 + (nCompressedCoeffs > 1), preCompContribPmap -> contribHeap ) != WAVELET_PADD4_NUMAPPROX + 1 + (nCompressedCoeffs > 1) ) error(SYSTEM, "can't write wavelet approximation coefficients to " "contribution heap in preComputeContrib()" ); /* Now dump mRGBE encoded, compressed detail coefficients */ if (putbinary(preCompContrib -> mrgbeCoeffs, sizeof(mRGBE), nCompressedCoeffs, preCompContribPmap -> contribHeap ) != nCompressedCoeffs ) error(SYSTEM, "can't write wavelet detail coefficients to " "contribution heap in preComputeContrib()" ); } #if 0 /* HACK: signal newPhoton() to set precomputed photon's contribution source from ray -> rsrc */ /* XXX: DO WE STILL NEED THIS CRAP AFTER PRECOMP, SINCE CHILD PMAPS ARE PER-MODIFIER ANYWAY? */ preCompContribPmap -> lastContribSrc.srcIdx = -2; #endif /* 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 /* Bail out if no photons could be precomputed */ if (!pmap -> numPhotons) error(USER, "no contribution photons precomputed; try increasing -am" ); #ifdef PMAP_CONTRIB_DBG printf("*DBG* preComputeContrib: Avg mRGBE encoding diff = %g\n", mrgbeDiffs / mrgbeDiffCnt ); #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, scdim, 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 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(); 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, scdim); 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/pmapdata.c b/pmapdata.c index b58ae6a..4d24674 100644 --- a/pmapdata.c +++ b/pmapdata.c @@ -1,869 +1,873 @@ #ifndef lint static const char RCSid[] = "$Id: pmapdata.c,v 2.23 2021/04/14 11:26:25 rschregle Exp $"; #endif /* ========================================================================= Photon map types and high level interface to nearest neighbour lookups in underlying point cloud data structure. The default data structure is an in-core kd-tree (see pmapkdt.{h,c}). This can be overriden with the PMAP_OOC compiletime switch, which enables an out-of-core octree (see oococt.{h,c}). Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) (c) Fraunhofer Institute for Solar Energy Systems, supported by the German Research Foundation (DFG LU-204/10-2, "Fassadenintegrierte Regelsysteme" (FARESYS)) (c) Lucerne University of Applied Sciences and Arts, supported by the Swiss National Science Foundation (SNSF #147053, "Daylight Redirecting Components", SNSF #179067, "Light Fields for Spatio-Temporal Glare Assessment") (c) Tokyo University of Science, supported by the JSPS Grants-in-Aid for Scientific Research (KAKENHI JP19KK0115, "Three-Dimensional Light Flow") ========================================================================= $Id: pmapdata.c,v 2.23 2021/04/14 11:26:25 rschregle Exp $ */ #include "pmapdata.h" #include "pmapdens.h" #include "pmaprand.h" #include "pmapmat.h" #include "pmaproi.h" #include "otypes.h" #include "random.h" PhotonMap *photonMaps [NUM_PMAP_TYPES] = { NULL, NULL, NULL, NULL, NULL, NULL, NULL, #ifdef PMAP_PHOTONFLOW NULL, NULL, #endif 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 -> preCompContrib = 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 Windows? */ fdFlags = fcntl(fileno(pmap -> heap), F_GETFL); fcntl(fileno(pmap -> heap), F_SETFL, fdFlags | O_APPEND); #endif/* ftruncate(fileno(pmap -> heap), 0); */ } } void flushPhotonHeap (PhotonMap *pmap) { int fd; const unsigned long len = pmap -> heapBufLen * sizeof(Photon); if (!pmap) error(INTERNAL, "undefined photon map in flushPhotonHeap()"); if (!pmap -> heap || !pmap -> heapBuf) { /* Silently ignore undefined heap error(INTERNAL, "undefined heap in flushPhotonHeap()"); */ return; } /* Atomically seek and write block to heap */ /* !!! Unbuffered I/O via pwrite() avoids potential race conditions * !!! and buffer corruption which can occur with lseek()/fseek() * !!! followed by write()/fwrite(). */ fd = fileno(pmap -> heap); #ifdef DEBUG_PMAP sprintf(errmsg, "Proc %d: flushing %ld photons from pos %ld\n", getpid(), pmap -> heapBufLen, lseek(fd, 0, SEEK_END) / sizeof(Photon) ); eputs(errmsg); #endif /*if (pwrite(fd, pmap -> heapBuf, len, lseek(fd, 0, SEEK_END)) != len) */ if (write(fd, pmap -> heapBuf, len) != len) { error(SYSTEM, "failed append to heap file in flushPhotonHeap()"); /* Clean up temp file */ fclose(pmap -> heap); unlink(pmap -> heapFname); } #if NIX if (fsync(fd)) error(SYSTEM, "failed fsync in flushPhotonHeap()"); #endif pmap -> heapBufLen = 0; } #ifdef DEBUG_PMAP static int checkPhotonHeap (FILE *file) /* Check heap for nonsensical or duplicate photons */ { Photon p, lastp; int i, dup; rewind(file); memset(&lastp, 0, sizeof(lastp)); while (fread(&p, sizeof(p), 1, file)) { dup = 1; for (i = 0; i <= 2; i++) { if (p.pos [i] < thescene.cuorg [i] || p.pos [i] > thescene.cuorg [i] + thescene.cusize ) { sprintf(errmsg, "corrupt photon in heap at [%f, %f, %f]\n", p.pos [0], p.pos [1], p.pos [2] ); error(WARNING, errmsg); } dup &= p.pos [i] == lastp.pos [i]; } if (dup) { sprintf(errmsg, "consecutive duplicate photon in heap " "at [%f, %f, %f]\n", p.pos [0], p.pos [1], p.pos [2] ); error(WARNING, errmsg); } } return 0; } #endif int newPhoton (PhotonMap* pmap, const RAY* ray) { unsigned i; Photon photon; COLOR photonFlux; /* Account for distribution ratio */ if (!pmap || pmapRandom(pmap -> randState) > pmap -> distribRatio) return -1; /* Don't store on sources */ if (ray -> robj > -1 && islight(objptr(ray -> ro -> omod) -> otype)) return -1; /* Ignore photon if modifier in/outside exclude/include set */ if (ambincl != -1 && ray -> ro && ambincl != inset(ambset, ray->ro->omod) ) return -1; /* Store photon if within a region of interest (for ze Ecksperts!) */ if (!photonInROI(ray)) return -1; /* Adjust flux according to distribution ratio and ray weight */ copycolor(photonFlux, ray -> rcol); #if 0 /* Factored out ray -> rweight as deprecated (?) for pmap, and infact erroneously attenuates volume photon flux based on extinction, which is already factored in by photonParticipate() */ scalecolor( photonFlux, ray -> rweight / (pmap -> distribRatio ? pmap -> distribRatio : 1) ); #else scalecolor(photonFlux, 1.0 / (pmap -> distribRatio ? pmap -> distribRatio : 1) ); #endif setPhotonFlux(&photon, photonFlux); /* Set photon position and flags */ VCOPY(photon.pos, ray -> rop); photon.flags = 0; photon.caustic = PMAP_CAUSTICRAY(ray); /* Set photon's subprocess index (to consolidate contribution sources after photon distribution completes) */ photon.proc = PMAP_GETRAYPROC(ray); switch (pmap -> type) { 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; case PMAP_TYPE_TRANSIENT: #ifdef PMAP_PHOTONFLOW case PMAP_TYPE_TRANSLIGHTFLOW: #endif /* Set cumulative path length for transient photon, corresponding to its time of flight. The cumulative path length is obtained relative to the maximum path length photonMaxPathDist. The last intersection distance is subtracted as this isn't factored in until photonRay() generates a new photon ray. */ photon.aux.pathLen = photonMaxPathDist - ray -> rmax + ray -> rot; break; default: /* Set photon's path ID from ray serial num (supplemented by photon.proc below) */ photon.aux.pathID = ray -> rno; } /* Set normal */ for (i = 0; i <= 2; i++) photon.norm [i] = 127.0 * (isVolumePmap(pmap) #ifdef PMAP_PHOTONFLOW || isLightFlowPmap(pmap) #endif ? ray -> rdir [i] : ray -> ron [i] ); if (!pmap -> heapBuf) { /* Lazily allocate heap buffa */ #if NIX /* Randomise buffa size to temporally decorellate flushes in * multiprocessing mode */ srandom(randSeed + getpid()); pmap -> heapBufSize = PMAP_HEAPBUFSIZE * (0.5 + frandom()); #else /* Randomisation disabled for single processes on WIN; also useful * for reproducability during debugging */ pmap -> heapBufSize = PMAP_HEAPBUFSIZE; #endif if (!(pmap -> heapBuf = calloc(pmap -> heapBufSize, sizeof(Photon)))) error(SYSTEM, "failed heap buffer allocation in newPhoton"); pmap -> heapBufLen = 0; } /* Photon initialised; now append to heap buffa */ memcpy(pmap -> heapBuf + pmap -> heapBufLen, &photon, sizeof(Photon)); if (++pmap -> heapBufLen >= pmap -> heapBufSize) /* Heap buffa full, flush to heap file */ flushPhotonHeap(pmap); pmap -> numPhotons++; /* Print photon attributes */ if (printPhoton) /* Non-const kludge */ printPhoton((RAY*)ray, &photon, pmap); return 0; } void buildPhotonMap (PhotonMap *pmap, double *photonFlux, PhotonContribSourceIdx *contribSrcOfs, unsigned nproc ) { unsigned long n, nCheck = 0; unsigned i; Photon *p; COLOR flux; char nuHeapFname [sizeof(PMAP_TMPFNAME)]; FILE *nuHeap; /* Need double here to reduce summation errors */ double avgFlux [3] = {0, 0, 0}, CoG [3] = {0, 0, 0}, CoGdist = 0, avgPathLen = 0; FVECT d; if (!pmap) error(INTERNAL, "undefined photon map in buildPhotonMap()"); /* Get number of photons from heapfile size */ if (fseek(pmap -> heap, 0, SEEK_END) < 0) error(SYSTEM, "failed seek to end of photon heap in buildPhotonMap()"); pmap -> numPhotons = ftell(pmap -> heap) / sizeof(Photon); if (!pmap -> numPhotons) error(INTERNAL, "empty photon map in buildPhotonMap()"); if (!pmap -> heap) error(INTERNAL, "no heap in buildPhotonMap()"); #ifdef DEBUG_PMAP eputs("Checking photon heap consistency...\n"); checkPhotonHeap(pmap -> heap); sprintf(errmsg, "Heap contains %ld photons\n", pmap -> numPhotons); eputs(errmsg); #endif /* Allocate heap buffa */ if (!pmap -> heapBuf) { pmap -> heapBufSize = PMAP_HEAPBUFSIZE; pmap -> heapBuf = calloc(pmap -> heapBufSize, sizeof(Photon)); if (!pmap -> heapBuf) error(SYSTEM, "failed to allocate postprocessed photon heap in buildPhotonMap()" ); } /* We REALLY don't need yet another @%&*! heap just to hold the scaled * photons, but can't think of a quicker fix... */ mktemp(strcpy(nuHeapFname, PMAP_TMPFNAME)); if (!(nuHeap = fopen(nuHeapFname, "w+b"))) error(SYSTEM, "failed to open postprocessed photon heap in buildPhotonMap()" ); rewind(pmap -> heap); #ifdef DEBUG_PMAP eputs("Postprocessing photons...\n"); #endif while (!feof(pmap -> heap)) { #ifdef DEBUG_PMAP printf("Reading %lu at %lu... ", pmap -> heapBufSize, ftell(pmap->heap) ); #endif pmap -> heapBufLen = fread(pmap -> heapBuf, sizeof(Photon), pmap -> heapBufSize, pmap -> heap ); #ifdef DEBUG_PMAP printf("Got %lu\n", pmap -> heapBufLen); #endif if (ferror(pmap -> heap)) error(SYSTEM, "failed to read photon heap in buildPhotonMap()"); for (n = pmap -> heapBufLen, p = pmap -> heapBuf; n; n--, p++) { /* Update min and max pos and set photon flux */ for (i = 0; i <= 2; i++) { if (p -> pos [i] < pmap -> minPos [i]) pmap -> minPos [i] = p -> pos [i]; else if (p -> pos [i] > pmap -> maxPos [i]) pmap -> maxPos [i] = p -> pos [i]; /* Update centre of gravity with photon position */ CoG [i] += p -> pos [i]; } if (isContribPmap(pmap) && contribSrcOfs) /* Linearise contrib photon origin index from subprocess index * using the per-subprocess offsets in contribSrcOfs */ p -> aux.contribSrc += contribSrcOfs [p -> proc]; else if (isTransientPmap(pmap) #ifdef PMAP_PHOTONFLOW || isTransLightFlowPmap(pmap) #endif ) { /* Update min and max path length */ if (p -> aux.pathLen < pmap -> minPathLen) pmap -> minPathLen = p -> aux.pathLen; else if (p -> aux.pathLen > pmap -> maxPathLen) pmap -> maxPathLen = p -> aux.pathLen; /* Update avg photon path length */ avgPathLen += p -> aux.pathLen; } /* Scale photon's flux (hitherto normalised to 1 over RGB); in * case of a contrib photon map, this is done per light source, * and photonFlux is assumed to be an array */ getPhotonFlux(p, flux); if (photonFlux) { scalecolor(flux, + #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; } /* !!! Flush pending writes to new heap !!! */ fflush(nuHeap); #ifdef DEBUG_PMAP if (nCheck < pmap -> numPhotons) error(INTERNAL, "truncated photon heap in buildPhotonMap"); #endif /* Finalise average photon flux */ scalecolor(avgFlux, 1.0 / pmap -> numPhotons); copycolor(pmap -> photonFlux, avgFlux); /* Average photon positions to get centre of gravity */ for (i = 0; i < 3; i++) pmap -> CoG [i] = CoG [i] /= pmap -> numPhotons; /* Average photon path lengths */ pmap -> avgPathLen = avgPathLen /= pmap -> numPhotons; rewind(pmap -> heap); /* Compute average photon distance to centre of gravity */ while (!feof(pmap -> heap)) { pmap -> heapBufLen = fread(pmap -> heapBuf, sizeof(Photon), pmap -> heapBufSize, pmap -> heap ); for (n = pmap -> heapBufLen, p = pmap -> heapBuf; n; n--, p++) { VSUB(d, p -> pos, CoG); CoGdist += DOT(d, d); if (isTransientPmap(pmap) #ifdef PMAP_PHOTONFLOW || isTransLightFlowPmap(pmap) #endif ) { /* Add temporal distance (in units of photon path length) for transient photons */ d [0] = p -> aux.pathLen - avgPathLen; CoGdist += d [0] * d [0]; } } } pmap -> CoGdist = CoGdist /= pmap -> numPhotons; /* Swap heaps, discarding unscaled photons */ fclose(pmap -> heap); unlink(pmap -> heapFname); pmap -> heap = nuHeap; strcpy(pmap -> heapFname, nuHeapFname); #ifdef PMAP_OOC OOC_BuildPhotonMap(pmap, nproc); #else if (isTransientPmap(pmap) #ifdef PMAP_PHOTONFLOW || isTransLightFlowPmap(pmap) #endif ) kdT_TransBuildPhotonMap(pmap); else kdT_BuildPhotonMap(pmap); #endif /* Trash heap and its buffa */ free(pmap -> heapBuf); fclose(pmap -> heap); unlink(pmap -> heapFname); pmap -> heap = NULL; pmap -> heapBuf = NULL; } /* PHOTON MAP SEARCH ROUTINES ------------------------------------------ */ /* Dynamic max photon search radius increase and reduction factors */ #define PMAP_MAXDIST_INC 4 #define PMAP_MAXDIST_DEC 0.9 /* Num successful lookups before reducing in max search radius */ #define PMAP_MAXDIST_CNT 1000 /* Threshold below which we assume increasing max radius won't help */ #define PMAP_SHORT_LOOKUP_THRESH 1 /* Coefficient for adaptive maximum search radius */ #define PMAP_MAXDIST_COEFF 100 void findPhotons (PhotonMap* pmap, const RAY* ray) { int redo = 0; const RREAL *norm = ray -> ron, *pos = ray -> rop; if (!pmap -> squeue.len) { /* Lazy init priority queue */ #ifdef PMAP_OOC OOC_InitFindPhotons(pmap); #else kdT_InitFindPhotons(pmap); #endif pmap -> minGathered = pmap -> maxGather; pmap -> maxGathered = pmap -> minGather; pmap -> totalGathered = 0; pmap -> numLookups = pmap -> numShortLookups = 0; pmap -> shortLookupPct = 0; pmap -> minError = FHUGE; pmap -> maxError = -FHUGE; pmap -> rmsError = 0; /* SQUARED max search radius limit is based on avg photon distance to * centre of gravity, unless fixed by user (maxDistFix > 0) */ pmap -> maxDist0 = pmap -> maxDist2Limit = maxDistFix > 0 ? maxDistFix * maxDistFix : PMAP_MAXDIST_COEFF * pmap -> squeue.len * pmap -> CoGdist / pmap -> numPhotons; } do { pmap -> squeue.tail = 0; pmap -> maxDist2 = pmap -> maxDist0; if (isVolumePmap(pmap) #ifdef PMAP_PHOTONFLOW || isLightFlowPmap(pmap) && sphericalIrrad #endif ) { /* Volume photons are not associated with a normal or intersection point; null the normal and set origin as lookup point. */ /* If a lightflow photon map is used and sphericalIrrad = 1, the mean spherical irradiance is evaluated, so normals are irrelevant and not filtered. */ norm = NULL; pos = ray -> rorg; } /* XXX/NOTE: status returned by XXX_FindPhotons() is currently ignored; if no photons are found, an empty queue is returned under the assumption all photons are too distant to contribute significant flux. */ #ifdef PMAP_OOC OOC_FindPhotons(pmap, pos, norm); #else if (isTransientPmap(pmap) #ifdef PMAP_PHOTONFLOW || isTransLightFlowPmap(pmap) #endif ) kdT_TransFindPhotons(pmap, pos, norm); else kdT_FindPhotons(pmap, pos, norm); #endif #ifdef PMAP_LOOKUP_INFO fprintf(stderr, "%d/%d %s photons found within radius %.3f " "at (%.2f,%.2f,%.2f) on %s\n", pmap -> squeue.tail, pmap -> squeue.len, pmapName [pmap -> type], sqrt(pmap -> maxDist2), ray -> rop [0], ray -> rop [1], ray -> rop [2], ray -> ro ? ray -> ro -> oname : "" ); #endif if (pmap -> squeue.tail < pmap -> squeue.len * pmap->gatherTolerance) { /* Short lookup; too few photons found */ if (pmap -> squeue.tail > PMAP_SHORT_LOOKUP_THRESH) { /* Ignore short lookups which return fewer than * PMAP_SHORT_LOOKUP_THRESH photons under the assumption there * really are no photons in the vicinity, and increasing the max * search radius therefore won't help */ #ifdef PMAP_LOOKUP_WARN sprintf(errmsg, "%d/%d %s photons found at (%.2f,%.2f,%.2f) on %s", pmap -> squeue.tail, pmap->squeue.len, pmapName [pmap->type], ray -> rop [0], ray -> rop [1], ray -> rop [2], ray -> ro ? ray -> ro -> oname : "" ); error(WARNING, errmsg); #endif /* Bail out after warning if maxDist is fixed */ if (maxDistFix > 0) return; if (pmap -> maxDist0 < pmap -> maxDist2Limit) { /* Increase max search radius if below limit & redo search */ pmap -> maxDist0 *= PMAP_MAXDIST_INC; #ifdef PMAP_LOOKUP_REDO redo = 1; #endif #ifdef PMAP_LOOKUP_WARN sprintf(errmsg, redo ? "restarting photon lookup with max radius %.1e" : "max photon lookup radius adjusted to %.1e", pmap -> maxDist0 ); error(WARNING, errmsg); #endif } #ifdef PMAP_LOOKUP_REDO else { sprintf(errmsg, "max photon lookup radius clamped to %.1e", pmap -> maxDist0 ); error(WARNING, errmsg); } #endif } /* Reset successful lookup counter */ pmap -> numLookups = 0; } else { /* Skip search radius reduction if maxDist is fixed */ if (maxDistFix > 0) return; /* Increment successful lookup counter and reduce max search radius if * wraparound */ pmap -> numLookups = (pmap -> numLookups + 1) % PMAP_MAXDIST_CNT; if (!pmap -> numLookups) pmap -> maxDist0 *= PMAP_MAXDIST_DEC; redo = 0; } } while (redo); } Photon *find1Photon (PhotonMap *pmap, const RAY* ray, Photon *photon, void *photonIdx ) { /* Init (squared) search radius to avg photon dist to centre of gravity */ float maxDist2_0 = pmap -> CoGdist; int found = 0; #ifdef PMAP_LOOKUP_REDO #define REDO 1 #else #define REDO 0 #endif do { pmap -> maxDist2 = maxDist2_0; #ifdef PMAP_OOC found = OOC_Find1Photon(pmap, ray -> rop, ray -> ron, photon, (OOC_DataIdx*)photonIdx ); #else found = kdT_Find1Photon(pmap, ray -> rop, ray -> ron, photon); #endif if (found < 0) { /* Expand search radius to retry */ maxDist2_0 *= 2; #ifdef PMAP_LOOKUP_WARN sprintf(errmsg, "failed 1-NN photon lookup" #ifdef PMAP_LOOKUP_REDO ", retrying with search radius %.2f", maxDist2_0 #endif ); error(WARNING, errmsg); #endif } } while (REDO && found < 0); /* Return photon buffer containing valid photon, else NULL */ return found < 0 ? NULL : photon; } void getPhoton (PhotonMap *pmap, PhotonIdx idx, Photon *photon) { #ifdef PMAP_OOC if (OOC_GetPhoton(pmap, idx, photon)) #else if (kdT_GetPhoton(pmap, idx, photon)) #endif error(INTERNAL, "failed photon lookup"); } Photon *getNearestPhoton (const PhotonSearchQueue *squeue, PhotonIdx idx) { #ifdef PMAP_OOC return OOC_GetNearestPhoton(squeue, idx); #else return kdT_GetNearestPhoton(squeue, idx); #endif } PhotonIdx firstPhoton (const PhotonMap *pmap) { #ifdef PMAP_OOC return OOC_FirstPhoton(pmap); #else return kdT_FirstPhoton(pmap); #endif } /* PHOTON MAP CLEANUP ROUTINES ------------------------------------------ */ void deletePhotons (PhotonMap* pmap) { #ifdef PMAP_OOC OOC_Delete(&pmap -> store); #else kdT_Delete(&pmap -> store); #endif free(pmap -> squeue.node); free(pmap -> biasCompHist); pmap -> numPhotons = pmap -> minGather = pmap -> maxGather = pmap -> squeue.len = pmap -> squeue.tail = 0; #ifdef PMAP_PATHFILT if (pmap -> pathLUT) { lu_done(pmap -> pathLUT); free(pmap -> pathLUT); } if (pmap -> pathLUTKeys) { unsigned k; for (k = 0; k < pmap->squeue.len + 1; free(pmap->pathLUTKeys [k++])); free(pmap -> pathLUTKeys); } #endif /* Contrib stuff */ if (isContribPmap(pmap) && pmap -> preCompContribTab) { /* Parent contrib photon map; clean up child photon maps containing precomputed contribs via lu_done() */ lu_done(pmap -> preCompContribTab); free(pmap -> preCompContribTab); pmap -> preCompContribTab = NULL; } } diff --git a/pmapooc.h b/pmapooc.h index f1774da..d0e936d 100644 --- a/pmapooc.h +++ b/pmapooc.h @@ -1,89 +1,89 @@ /* ================================================================== 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") (c) Tokyo University of Science, supported by the JSPS Grants-in-Aid for Scientific Research (KAKENHI JP19KK0115, "Three-Dimensional Light Flow") ================================================================== $Id: pmapooc.h,v 1.1 2020/10/29 16:20:11 u-no-hoo Exp u-no-hoo $ */ #ifndef PMAPOOC_H #define PMAPOOC_H #include "oocnn.h" /* Suffixes for octree filenames */ /* #define PMAP_OOC_NODESUFFIX ".node" #define PMAP_OOC_HEAPSUFFIX ".heap" */ #define PMAP_OOC_LEAFSUFFIX ".leaf" /* Out-of-core octree constants */ #define PMAP_OOC_NUMBLK 32 /* Num blocks for external sort */ #define PMAP_OOC_BLKSIZE 1e8 /* Block size for external sort */ #define PMAP_OOC_LEAFMAX (OOC_OCTCNT_MAX) /* Max photons per leaf */ #define PMAP_OOC_MAXDEPTH (MORTON3D_BITS) /* Max octree depth */ typedef OOC_SearchQueueNode PhotonSearchQueueNode; typedef OOC_SearchQueue PhotonSearchQueue; typedef OOC_Octree PhotonStorage; typedef unsigned PhotonIdx; /* Forward declarations to break dependency loop with pmapdata.h */ struct PhotonMap; - + void OOC_BuildPhotonMap (struct PhotonMap *pmap, unsigned numProc); /* Build out-of-core octree pmap -> store from photons in unsorted * heapfile pmap -> heap and generate nodes and leaf file with prefix * pmap -> fileName. Photon map construction may be parallelised if * numProc > 1, if supported. The heap is destroyed on return. */ int OOC_SavePhotons (const struct PhotonMap *pmap, FILE *out); /* Save photons in out-of-core octree to file. Return -1 on error, * else 0 */ int OOC_LoadPhotons (struct PhotonMap *pmap, FILE *in); /* Load photons in out-of-core octree from file. Return -1 on error, * else 0 */ void OOC_InitFindPhotons (struct PhotonMap *pmap); /* Initialise NN search queue prior to calling kdT_FindPhotons() */ int OOC_FindPhotons ( struct PhotonMap* pmap, const FVECT pos, const FVECT norm ); /* Locate pmap -> squeue.len nearest photons to pos with similar normal * (NULL for volume photons) and return in search queue pmap -> squeue, * starting with the further photon at pmap -> squeue.node. Return -1 * if none found, else 0. */ int OOC_Find1Photon (struct PhotonMap* pmap, const FVECT pos, const FVECT norm, Photon *photon, OOC_DataIdx *photonIdx ); /* Locate single nearest photon to pos with similar normal. Return -1 * if none found, else 0. */ int OOC_GetPhoton (struct PhotonMap *pmap, PhotonIdx idx, Photon *photon); /* Retrieve photon referenced by idx from leaf file and return -1 on * error, else 0. */ Photon *OOC_GetNearestPhoton (const PhotonSearchQueue *squeue, PhotonIdx idx ); /* Retrieve photon from NN search queue after OOC_FindPhotons() */ PhotonIdx OOC_FirstPhoton (const struct PhotonMap* pmap); /* Return index to first photon in octree */ #endif diff --git a/pmcontrib2a.c b/pmcontrib2a.c new file mode 100644 index 0000000..2fe3b21 --- /dev/null +++ b/pmcontrib2a.c @@ -0,0 +1,677 @@ +#ifndef lint +static const char RCSid[] = "$Id$"; +#endif + +/* + ========================================================================= + Photon map routines for precomputed light source contributions. + This is the main photon distribution routine used by mkpmap. + + !!! THIS EXPERIMENTAL VERSION USES IMPORTANCE SAMPLING OF LIGHT SOURCE + !!! FLUX, SIMILAR TO distribPhotons(). CONSEQUENTLY, IT DOES NOT + !!! MAINTAIN A LIST OF PER-SOURCE FLUX, AND MAY UNDERESTIMATE SOURCES + !!! WITH LOW CONTRIBUTIONS. + + 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 "pmapdiag.h" +#include "pmaprand.h" +#include "pmapmat.h" +#include "pmaproi.h" +#include "pmapsrc.h" +#include "pmapio.h" +#include "otspecial.h" +#if NIX + #include + #include + #include +#endif + + + +/* Defs for photon emission counter array passed by sub-processes to parent + * via shared memory */ +typedef unsigned long PhotonContribCnt; + +/* Indices for photon emission counter array: num photons stored and num + * emitted per source */ +#define PHOTONCNT_NUMPHOT 0 +#define PHOTONCNT_NUMEMIT 1 +#define PHOTONCNT_SIZE 2 + + + +static void shmLock (int shmFile, int type) +/* Set/release record lock on shared memory file, + where type is either F_RDLCK, F_WRLCK or F_UNLCK*/ +{ + struct flock lock; + + lock.l_type = type; + lock.l_whence = SEEK_SET; + lock.l_start = lock.l_len = 0; + + /* F_SETLKW waits until file is released by other process. A loop handles + interrupt by a signal (error EINTR). */ + do + if (fcntl(shmFile, F_SETLKW, &lock) >= 0) + /* Okie dokie */ + return; + while (errno == EINTR); + + error(SYSTEM, + "failed to (un)lock shared memory in distribPhotonContrib()" + ); +} + + + +void distribPhotonContrib (PhotonMap* pmap, LUTAB *contribTab, + unsigned numContribs, int *contribMode, unsigned numProc +) +{ + EmissionMap emap; + char errmsg2 [128], shmFname [PMAP_TMPFNLEN]; + unsigned srcIdx, proc; + int shmFile, stat, pid; + double totalFlux = 0; /* Total emitted source flux */ + unsigned long *photonCnt; /* Photon counter array */ + FILE **contribSrcHeap = NULL; + char **contribSrcHeapFname = NULL; + PhotonContribSourceIdx *contribSrcOfs = NULL; + pid_t procPids [PMAP_MAXPROC]; + + if (!pmap) + error(USER, "no contribution photon map specified"); + + if (!nsources) + error(USER, "no light sources"); + + if (nsources > PMAP_MAXSRCIDX) + error(USER, "too many light sources"); + + if (!contribTab || !numContribs) + error(USER, "no modifiers specified for contribution photon map"); + + /* =================================================================== + * INITIALISATION - Set up emission and scattering funcs + * =================================================================== */ + emap.samples = NULL; + emap.src = NULL; + emap.maxPartitions = MAXSPART; + emap.partitions = (unsigned char*)malloc(emap.maxPartitions >> 1); + if (!emap.partitions) + error(USER, "can't allocate source partitions in distribPhotonContrib"); + + /* Initialise contrib photon map */ + initPhotonMap(pmap, PMAP_TYPE_CONTRIB); + initPmapContribTab(contribTab, contribMode); + pmap -> contribMode = *contribMode; + initPmapContrib(pmap); + initPhotonHeap(pmap); + initPhotonEmissionFuncs(); + initPhotonScatterFuncs(); + + /* Per-subprocess target counts */ + if (!(pmap -> distribTarget /= numProc)) + error(INTERNAL, "no photons to distribute in distribPhotonContrib"); + + /* Get photon ports from modifier list */ + getPhotonPorts(photonPortList); + + /* Get photon sensor modifiers */ + getPhotonSensors(photonSensorList); + + /* Get polyhedral regions of interest */ + getPolyROIs(pmapROImodList); + +#if NIX + /* Set up shared mem for photon counters (zeroed by ftruncate) */ + strcpy(shmFname, PMAP_TMPFNAME); + shmFile = mkstemp(shmFname); + + if (shmFile < 0 || ftruncate(shmFile, PHOTONCNT_SIZE) < 0) + error(SYSTEM, "failed shared mem init in distribPhotonContrib"); + + photonCnt = mmap(NULL, PHOTONCNT_SIZE, + PROT_READ | PROT_WRITE, MAP_SHARED, shmFile, 0 + ); + + if (photonCnt == MAP_FAILED) + error(SYSTEM, "failed mapping shared memory in distribPhotonContrib"); +#else + photonCnt [PHOTONCNT_NUMPHOT] = photonCnt [PHOTONCNT_NUMEMIT] = 0; +#endif /* NIX */ + + if (verbose) { + sprintf(errmsg, "\nIntegrating flux from %d sources", nsources); + + if (photonPorts) { + sprintf(errmsg2, " via %d ports", numPhotonPorts); + strcat(errmsg, errmsg2); + } + + strcat(errmsg, "\n"); + eputs(errmsg); + } + + /* ============================================================= + * FLUX INTEGRATION - Get total flux emitted from sources/ports + * ============================================================= */ + for (srcIdx = 0; srcIdx < nsources; srcIdx++) { + unsigned portCnt = 0; + const OBJREC *srcMod = findmaterial(source [srcIdx].so); + + /* Skip this source if its contributions are not sought */ + if (!lu_find(pmapContribTab, srcMod -> oname) -> data) { + sprintf(errmsg, "ignoring contributions from source %s", + source [srcIdx].so -> oname + ); + error(WARNING, errmsg); + continue; + } + + emap.src = source + srcIdx; + + do { /* Need at least one iteration if no ports! */ + emap.port = emap.src -> sflags & SDISTANT + ? photonPorts + portCnt + : NULL; + photonPartition [emap.src -> so -> otype] (&emap); + + if (verbose) { + sprintf(errmsg, "\tIntegrating flux from source %s ", + source [srcIdx].so -> oname + ); + + if (emap.port) { + sprintf(errmsg2, "via port %s ", + photonPorts [portCnt].so -> oname + ); + strcat(errmsg, errmsg2); + } + + sprintf(errmsg2, "(%lu partitions)\n", emap.numPartitions); + strcat(errmsg, errmsg2); + eputs(errmsg); +#if NIX + fflush(stderr); +#endif + } + + for (emap.partitionCnt = 0; + emap.partitionCnt < emap.numPartitions; + emap.partitionCnt++ + ) { + initPhotonEmission(&emap, pdfSamples); + totalFlux += colorAvg(emap.partFlux); + } + + portCnt++; + } while (portCnt < numPhotonPorts); + } + + if (totalFlux < FTINY) + error(USER, "zero flux from light sources"); + + /* Allocate & init per-subprocess contribution source heap files */ + contribSrcHeap = calloc(numProc, sizeof(FILE*)); + contribSrcHeapFname = calloc(numProc, sizeof(char*)); + contribSrcOfs = calloc(numProc, sizeof(PhotonContribSourceIdx)); + if (!contribSrcHeap || !contribSrcHeapFname || !contribSrcOfs) + error(SYSTEM, "failed contribution source heap allocation " + "in distribPhotonContrib()" + ); + + for (proc = 0; proc < numProc; proc++) { + contribSrcHeapFname [proc] = malloc(PMAP_TMPFNLEN); + if (!contribSrcHeapFname [proc]) + error(SYSTEM, "failed contribution source heap file allocation " + "in distribPhotonContrib()" + ); + + mktemp(strcpy(contribSrcHeapFname [proc], PMAP_TMPFNAME)); + if (!(contribSrcHeap [proc] = fopen(contribSrcHeapFname [proc], "w+b"))) + error(SYSTEM, "failed opening contribution source heap file " + "in distribPhotonContrib()" + ); + } + + /* Record start time for progress reports */ + repStartTime = time(NULL); + + if (verbose) { + sprintf(errmsg, "\nPhoton distribution @ %d procs\n", numProc); + eputs(errmsg); + } + + /* MAIN LOOP */ + for (proc = 0; proc < numProc; proc++) { +#if NIX + if (!(pid = fork())) { + /* SUBPROCESS ENTERS HERE; opened and mmapped files inherited */ +#else + if (1) { + /* No subprocess under Windoze */ +#endif + /* Local photon counters for this subprocess */ + unsigned passCnt = 0, prePassCnt = 0; + unsigned long preIncPhotonCnt, lastNumPhotons = 0, + numEmitted = 0, lastNumEmitted = 0; + /* Decorrelate shared photon counter updates to reduce + * contention by setting unique update intervals per subproc */ + const unsigned photonCntUpdate = PMAP_CNTUPDATE + proc; + + /* Seed RNGs from PID for decorellated photon distribution */ + pmapSeed(randSeed + proc, partState); + pmapSeed(randSeed + (proc + 1) % numProc, emitState); + pmapSeed(randSeed + (proc + 2) % numProc, cntState); + pmapSeed(randSeed + (proc + 3) % numProc, mediumState); + pmapSeed(randSeed + (proc + 4) % numProc, scatterState); + pmapSeed(randSeed + (proc + 5) % numProc, rouletteState); +#ifdef DEBUG_PMAP + /* Output child process PID after random delay to prevent corrupted + * console output due to race condition */ + usleep(1e6 * pmapRandom(rouletteState)); + fprintf(stderr, + "Proc %d: PID = %d (waiting 10 sec to attach debugger...)\n", + proc, getpid() + ); + /* Allow time for debugger to attach to child process */ + sleep(10); +#endif + /* ============================================================= + * 2-PASS PHOTON DISTRIBUTION + * Pass 1 (pre): emit fraction of target photon count + * Pass 2 (main): based on outcome of pass 1, estimate remaining + * number of photons to emit to approximate target + * count + * ============================================================= */ + while (passCnt < 2) { + double numEmit; + + if (!passCnt) { + /* INIT PASS 1 */ + if (++prePassCnt > maxPreDistrib) { + /* Skip if no photons contributed after sufficient + * iterations */ + sprintf(errmsg, + "proc %d: too many prepasses, no photons stored; " + "increase -apD?", proc + ); + error(USER, errmsg); + break; + } + + /* Num to emit is fraction of target count */ + numEmit = preDistrib * pmap -> distribTarget; + } + else { + /* INIT PASS 2 */ + /* Based on the outcome of the predistribution we can now + * figure out how many more photons we have to emit from + * the current source to meet the target count, + * distribTarget. This value is clamped to 0 in case + * the target has already been exceeded in pass 1. + * srcNumEmit and srcNumDistrib is the number of photons + * emitted and distributed (stored) from the current + * source in pass 1, respectively. */ + numEmit = numEmitted * (pmap -> numPhotons + ? max((float)pmap -> distribTarget / pmap -> numPhotons, 1.) - 1. + : 0. + ); + + if (!numEmit) + /* No photons left to distribute in main pass */ + break; + } + + /* PHOTON DISTRIBUTION LOOP */ + for (srcIdx = 0; srcIdx < nsources; srcIdx++) { + unsigned portCnt = 0; + emap.src = source + srcIdx; + + do { /* Need at least one iteration if no ports! */ + emap.port = emap.src -> sflags & SDISTANT + ? photonPorts + portCnt + : NULL; + photonPartition [emap.src -> so -> otype] (&emap); + + if (verbose && !proc) { + /* Output from subproc 0 only to avoid race condition + * on console I/O */ + if (!passCnt) + sprintf(errmsg, "\tPREPASS %d on source %s ", + prePassCnt, source [srcIdx].so -> oname + ); + else + sprintf(errmsg, "\tMAIN PASS on source %s ", + source [srcIdx].so -> oname + ); + + if (emap.port) { + sprintf(errmsg2, "via port %s ", + photonPorts [portCnt].so -> oname + ); + strcat(errmsg, errmsg2); + } + + sprintf(errmsg2, "(%lu partitions)\n", + emap.numPartitions + ); + strcat(errmsg, errmsg2); + eputs(errmsg); +#if NIX + fflush(stderr); +#endif + } + + for (emap.partitionCnt = 0; + emap.partitionCnt < emap.numPartitions; + emap.partitionCnt++ + ) { + double partNumEmit; + unsigned long partEmitCnt; + + /* Get photon origin within current source partishunn + * and build emission map */ + photonOrigin [emap.src -> so -> otype] (&emap); + initPhotonEmission(&emap, pdfSamples); + + /* Number of photons to emit from ziss partishunn; + * scale according to its normalised contribushunn to + * the emitted source flux */ + partNumEmit = numEmit * colorAvg(emap.partFlux) / + totalFlux; + partEmitCnt = (unsigned long)partNumEmit; + + /* Probabilistically account for fractional photons */ + if (pmapRandom(cntState) < partNumEmit - partEmitCnt) + partEmitCnt++; + + /* Integer counter avoids FP rounding errors during + * iteration */ + while (partEmitCnt--) { + RAY photonRay; + + /* Emit photon according to PDF (if any). If + * accepted, allocate associated contribution + * origin, and trace through scene until + * absorbed/leaked; emitPhoton() sets the emitting + * light source index in photonRay */ + /* NOTE: rejection sampling skips incrementing the + * emission counter (see below), thus compensating + * for the rejected photons by increasing the photon + * flux in proportion to the lower effective + * emission count. + * BUG: THIS INTERFERES WITH THE PROGRESS COUNTER + * REPORTED TO THE PARENT, AND WITH THE + * PREDISTRIBUTION PASS --> PHOTON DISTRIBUTION WILL + * FINISH EARLY, WITH FEWER PHOTONS THAN TARGETTED! + */ + if (!emitPhoton(&emap, &photonRay)) + continue; + + newPhotonContribSource(pmap, &photonRay, + contribSrcHeap [proc] + ); + + /* Update local emission counter */ + numEmitted++; + + /* Skip photon if it has an invalid bin. It will + implicitly contribute zero flux, so don't bother + tracing and storing it. However, it counts as + emitted to avoid bias (see enclosing while() loop + and local emission counter update above). */ + if (pmap -> lastContribSrc.srcBin < 0) + continue; + + /* Set subprocess index in photonRay for post- + * distrib contribution source index linearisation; + * this is propagated with the contrib source index + * in photonRay and set for photon hits by + * newPhoton() */ + PMAP_SETRAYPROC(&photonRay, proc); + tracePhoton(&photonRay); + + if (!(partEmitCnt && (numEmitted & photonCntUpdate))) { + /* Update global counters shared with siblings + in unique intervals to reduce overhead and + contention; + ALSO WHEN EMITTING FINAL PHOTON FROM THIS + PARTITION, OTHERWISE COUNTERS WILL SKIP LAST + UPDATE, BIASING THE PHOTON FLUX! */ + /* Shared photon counter file must be locked! */ + shmLock(shmFile, F_WRLCK); + /* Differentially increment number emitted using + * local counter numEmitted */ + preIncPhotonCnt = photonCnt [PHOTONCNT_NUMEMIT]; + photonCnt [PHOTONCNT_NUMEMIT] += + numEmitted - lastNumEmitted; + lastNumEmitted = numEmitted; + + /* Check overflow using pre-increment values */ + if (photonCnt [PHOTONCNT_NUMEMIT] < preIncPhotonCnt) { + sprintf(errmsg, "photon emission counter " + "overflow (was: %ld, is: %ld)", + preIncPhotonCnt, photonCnt [PHOTONCNT_NUMEMIT] + ); + error(INTERNAL, errmsg); + } + + /* Differentially increment photon counter */ + preIncPhotonCnt = photonCnt [PHOTONCNT_NUMPHOT]; + photonCnt [PHOTONCNT_NUMPHOT] += + pmap -> numPhotons - lastNumPhotons; + lastNumPhotons = pmap -> numPhotons; + + /* Check for photon counter overflow (this could + only happen before an emission counter overflow + if the scene has an absurdly high albedo and/or + very dense geometry) */ + if (photonCnt [PHOTONCNT_NUMPHOT] < preIncPhotonCnt) { + sprintf(errmsg, "photon counter overflow " + "(was: %ld, is: %ld)", preIncPhotonCnt, + photonCnt [PHOTONCNT_NUMPHOT] + ); + error(INTERNAL, errmsg); + } + + /* Release lock on shared photon counter file! */ + shmLock(shmFile, F_UNLCK); + } + } +#if !NIX + /* Synchronous progress report on Windoze */ + if (!proc && photonRepTime > 0 && + time(NULL) >= repLastTime + photonRepTime + ) { + unsigned s; + repComplete = pmap -> distribTarget * numProc; + repProgress = photonCnt [PHOTONCNT_NUMPHOT]; + + for (repEmitted = 0, s = 0; s < nsources; s++) + repEmitted += photonCnt [numEmitIdx]; + + pmapDistribReport(); + } +#endif + } + + portCnt++; + } while (portCnt < numPhotonPorts); + + if (!pmap -> numPhotons) { + /* Double predistrib factor in case no photons were stored + * for this source and redo pass 1 */ + preDistrib *= 2; + } + else { + /* Now do pass 2 */ + passCnt++; + } + } + } + + /* Flush heap buffa one final time to prevent data corruption */ + flushPhotonHeap(pmap); + /* Flush last contribution origin to origin heap file */ + newPhotonContribSource(pmap, NULL, contribSrcHeap [proc]); + /* Heap files closed automatically on exit + fclose(pmap -> heap); + fclose(orgHeap [proc]); */ +#ifdef DEBUG_PMAP + sprintf( + errmsg, "Proc %d total %ld photons\n", proc, pmap -> numPhotons + ); + eputs(errmsg); + fflush(stderr); +#endif +#ifdef PMAP_SIGUSR + signal(SIGUSR1, SIG_DFL); +#endif +#if NIX + /* Terminate subprocess */ + exit(0); +#endif + } + else { + /* PARENT PROCESS CONTINUES LOOP ITERATION HERE */ + if (pid < 0) + error(SYSTEM, + "failed to fork subprocess in distribPhotonContrib()" + ); + else + /* Saves child process IDs */ + procPids [proc] = pid; + } + } + +#if NIX + /* PARENT PROCESS CONTINUES HERE */ +#ifdef SIGCONT + /* Enable progress report signal handler */ + signal(SIGCONT, pmapDistribReport); +#endif + /* Wait for subprocesses to complete while reporting progress */ + proc = numProc; + while (proc) { + while (waitpid(-1, &stat, WNOHANG) > 0) { + /* Subprocess exited; check status */ + if (!WIFEXITED(stat) || WEXITSTATUS(stat)) { + /* Exited with error; terminate siblings, clean up & bail out */ + for (proc = 0; proc < numProc; proc++) + kill(procPids [proc], SIGKILL); + + /* Unmap shared memory, squish mapped file */ + munmap(photonCnt, PHOTONCNT_SIZE); + close(shmFile); + unlink(shmFname); + error(USER, "failed photon distribution"); + } + + --proc; + } + + /* Nod off for a bit and update progress */ + sleep(1); + + /* Asynchronous progress report from shared subprocess counters */ + shmLock(shmFile, F_RDLCK); + repProgress = photonCnt [PHOTONCNT_NUMPHOT]; + pmap -> numPhotons = photonCnt [PHOTONCNT_NUMPHOT]; + + for (repEmitted = 0, srcIdx = 0; srcIdx < nsources; srcIdx++) { + repEmitted += photonCnt [PHOTONCNT_NUMEMIT]; + } + + repComplete = pmap -> distribTarget * numProc; + shmLock(shmFile, F_UNLCK); + + if (photonRepTime > 0 && time(NULL) >= repLastTime + photonRepTime) + pmapDistribReport(); +#ifdef SIGCONT + else signal(SIGCONT, pmapDistribReport); +#endif + } +#endif /* NIX */ + + /* ================================================================ + * POST-DISTRIBUTION - Set photon flux and build kd-tree, etc. + * ================================================================ */ +#ifdef SIGCONT + /* Reset signal handler */ + signal(SIGCONT, SIG_DFL); +#endif + free(emap.samples); + + if (!pmap -> numPhotons) + error(USER, "empty contribution photon map"); + + /* Load per-subprocess contribution sources into pmap -> contribSrc */ + /* Dumb compilers apparently need the char** cast */ + pmap -> numContribSrc = buildContribSources(pmap, contribSrcHeap, + (char**)contribSrcHeapFname, contribSrcOfs, numProc + ); + if (!pmap -> numContribSrc) + error(INTERNAL, "no contribution sources in contribution photon map"); + + /* Set photon flux */ + totalFlux = photonCnt [PHOTONCNT_NUMEMIT] + ? totalFlux / photonCnt [PHOTONCNT_NUMEMIT] + : 0; +#if NIX + /* Photon counters no longer needed, unmap shared memory */ + munmap(photonCnt, PHOTONCNT_SIZE); + close(shmFile); + unlink(shmFname); +#else + free(photonCnt); +#endif + + if (verbose) { + eputs("\nBuilding contribution photon map...\n"); +#if NIX + fflush(stderr); +#endif + } + + /* Build underlying data structure; heap is destroyed. + XXX: Photon flux remains unchanged (i.e. normalised) in coefficient + mode, so totalFlux is omitted */ + buildPhotonMap(pmap, pmap -> contribMode ? &totalFlux : NULL, + contribSrcOfs, numProc + ); + + /* Precompute binned photon contributions, init lookup table pmap -> + preCompContribTab containing constituent pre-modifier photon maps */ + if (verbose) + eputs("\n"); + preComputeContrib(contribPmap); + + /* Free per-subprocess origin heap files */ + for (proc = 0; proc < numProc; proc++) + free(contribSrcHeapFname [proc]); + + free(contribSrcHeapFname); + free(contribSrcHeap); + free(contribSrcOfs); + + if (verbose) + eputs("\n"); +} + +#endif /* PMAP_CONTRIB */ diff --git a/pmcontrib4.c b/pmcontrib4.c index 37000a4..06f29b4 100644 --- a/pmcontrib4.c +++ b/pmcontrib4.c @@ -1,652 +1,800 @@ #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 "pmcontribcache.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" #include "random.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/coefficient mode if it doesn't match precomputed */ if (*pmapContribMode != preCompContribPmap -> contribMode) { sprintf(errmsg, "photon map contains precomputed %s, overriding rcontrib setting", preCompContribPmap -> contribMode ? "contributions" : "coefficients" ); error(WARNING, errmsg); *pmapContribMode = preCompContribPmap -> contribMode; } /* NOTE: preCompContribPmap -> preCompContrib is initialised by * loadPhotonMap() */ preCompContrib = (PreComputedContrib*)( preCompContribPmap -> preCompContrib ); /* Set number of bins and wavelet coefficients */ preCompContrib -> nBins = preCompContrib -> scDim * preCompContrib -> scDim; preCompContrib -> nCoeffs = preCompContrib -> coeffDim * preCompContrib -> coeffDim; /* Set wavelet coefficient filename */ preCompContrib -> waveletFname = malloc(strlen(dirName) + strlen(modName) + strlen(PMAP_CONTRIB_WAVELETFILE) ); if (!preCompContribPmap -> fileName) error(SYSTEM, "out of memory allocating filename in " "loadPreCompContrib()" ); snprintf(preCompContrib -> waveletFname, sizeof(dirName) + 5, PMAP_CONTRIB_WAVELETFILE, dirName, modName ); return 0; } void loadContribPhotonMap (PhotonMap *pmap, const char *fname) /* Load contribution pmap and its per-modifier precomputed children */ { LUTAB lutInit = LU_SINIT(NULL, freePreCompContribNode); /* Allocate & init LUT containing per-modifier child contrib pmaps */ if (!(pmap -> preCompContribTab = malloc(sizeof(lutInit)))) error(SYSTEM, "out of memory allocating precomputed contribution " "LUT in loadContribPhotonMap()" ); memcpy(pmap -> preCompContribTab, &lutInit, sizeof(lutInit)); /* NOTE: filename already set in loadPhotonMap() pmap -> fileName = (char*)fname; */ /* Init contrib photon map for light source contributions */ initPmapContrib(pmap); /* Load child contribution photon maps for each modifier in contribTab */ lu_doall(pmapContribTab, loadPreCompContrib, pmap); } static int decodeContribs (PreComputedContrib *preCompContrib) /* Decode and decompress mRGBE-encoded wavelet coefficients in preCompContrib -> mrgbeCoeffs, apply inverse wavelet transform and return decoded contributions in the wavelet coefficient matrix preCompContrib -> waveletMatrix. NOTE: THE WAVELET COEFFICIENT MATRIX IS ASSUMED TO BE ZEROED, with the wavelet approximation coefficients in the upper left of the matrix (i.e. preCompContrib -> waveletMatrix [0..approxDim-1][0..approxDim-1], where approxDim = WAVELET_PADD4_APPROXDIM). Returns 0 on success, else -1. */ { unsigned c, i, j, coeffIdx, scDim, coeffDim, nCoeffs, nCompressedCoeffs; WaveletMatrix2 waveletMatrix; mRGBERange *mrgbeRange; mRGBE *mrgbeCoeffs; DCOLOR fCoeff; if (!preCompContrib || !preCompContrib -> mrgbeCoeffs || !(waveletMatrix = preCompContrib -> waveletMatrix) || !preCompContrib -> tWaveletMatrix ) /* Should be initialised by caller! */ return -1; scDim = preCompContrib -> scDim; coeffDim = preCompContrib -> coeffDim; nCoeffs = preCompContrib -> nCoeffs; nCompressedCoeffs = preCompContrib -> nCompressedCoeffs; mrgbeCoeffs = preCompContrib -> mrgbeCoeffs; mrgbeRange = &preCompContrib -> mrgbeRange; /* Init mRGBE coefficient normalisation from range */ if (!mRGBEinit(mrgbeRange, mrgbeRange -> min, mrgbeRange -> max)) return -1; /* Decode mRGBE detail coeffs and populate wavelet coefficient matrix, based on the encoded incremental coefficient index. This omits the approximation coefficients in the upper left of the coefficient matrix (these should have already been set by the caller). NOTE: The detail coefficients in the matrix are assumed to be initialised to zero to account for those that were dropped by thresholding. */ for (c = coeffIdx = 0; c < nCompressedCoeffs; c++) { coeffIdx += mRGBEdecode(mrgbeCoeffs [c], mrgbeRange, fCoeff); i = PMAP_CONTRIB_LIN2X(coeffDim, coeffIdx); j = PMAP_CONTRIB_LIN2Y(coeffDim, coeffIdx); #ifdef PMAP_CONTRIB_DBG /* Check for invalid decoded coefficients */ if (i < WAVELET_PADD4_APPROXDIM && j < WAVELET_PADD4_APPROXDIM || coeffIdx >= nCoeffs ) error(CONSISTENCY, "wavelet coefficient index out of range in decodeContribs()" ); #endif copycolor(waveletMatrix [i] [j], fCoeff); } /* Do inverse 2D wavelet transform on preCompContrib -> waveletMatrix */ if (!padWaveletInvXform2(waveletMatrix, preCompContrib -> tWaveletMatrix, coeffDim, scDim ) ) error(INTERNAL, "failed inverse wavelet transform in decodeContribs()" ); +#if 0 #ifdef PMAP_CONTRIB_DBG for (i = 0; i < scDim; i++) { for (j = 0; j < scDim; j++) { #if 0 /* Dump decoded contribs for debugging */ printf("% 7.3g\t", colorAvg(waveletMatrix [i][j])); #endif /* Warn if coefficient is negative; this _can_ happen due to loss of precision by the mRGBE encoding */ if (min(min(waveletMatrix [i][j][0], waveletMatrix [i][j][1]), waveletMatrix [i][j][2] ) < 0 ) error(WARNING, "negative contributions in getPreCompContrib()"); } #if 0 putchar('\n'); } putchar('\n'); #else } #endif +#endif #endif return 0; } static void getPreCompContribByPhoton (const Photon *photon, OOC_DataIdx photonIdx, PreComputedContrib *preCompContrib, DCOLOR *binnedContribs ) /* Fetch and decode mRGBE-encoded wavelet coefficients for given photon from wavelet file at offset photonIdx, perform inverse wavelet xform, and return the reconstructed binned contributions in binnedContribs (which is assumed to be variable). Returns 0 on success, else -1. */ /* NOTE: binnedContribs IS ASSUMED TO POINT TO CACHE PAGE DATA! */ { unsigned i, j, k, coeffIdx, scDim, coeffDim, nCompressedCoeffs; COLOR fCoeff; COLR rgbe32 [WAVELET_PADD4_NUMAPPROX + 2]; WaveletMatrix2 waveletMatrix = NULL; if (!binnedContribs || !preCompContrib) /* Shouldn't happen */ error(INTERNAL, "contributions not initialised in getPreCompContrib()" ); if (preCompContrib -> nBins <= 1) /* No binning, so nothing to decode */ return; scDim = preCompContrib -> scDim; coeffDim = preCompContrib -> coeffDim; nCompressedCoeffs = preCompContrib -> nCompressedCoeffs; /* Lazily initialise preCompContrib */ if (!preCompContrib -> waveletFile) { /* Lazily open wavelet coefficient file */ preCompContrib -> waveletFile = fopen( preCompContrib -> waveletFname, "rb" ); if (!preCompContrib -> waveletFile) { sprintf(errmsg, "can't open wavelet coefficient file %s " "in getPreCompContrib()", preCompContrib -> waveletFname ); error(SYSTEM, errmsg); } /* Set record size of encoded coeffs in wavelet file */ preCompContrib -> contribSize = PMAP_CONTRIB_ENCSIZE( nCompressedCoeffs ); /* Lazily allocate buffer for mRGBE wavelet coeffs */ preCompContrib -> mrgbeCoeffs = calloc(nCompressedCoeffs, sizeof(mRGBE) ); if (!preCompContrib -> mrgbeCoeffs) error(SYSTEM, "out of memory allocating decoded wavelet " "coefficients in getPreCompContrib()" ); /* Lazily allocate primary and transposed wavelet matrices */ preCompContrib -> waveletMatrix = allocWaveletMatrix2(coeffDim); preCompContrib -> tWaveletMatrix = allocWaveletMatrix2(coeffDim); if (!preCompContrib -> waveletMatrix || !preCompContrib -> tWaveletMatrix ) error(SYSTEM, "out of memory allocating wavelet coefficient " "matrix in getPreCompContrib()" ); } /* Seek to photon index in wavelet file and read its associated coefficients */ if (fseek(preCompContrib -> waveletFile, photonIdx * preCompContrib -> contribSize, SEEK_SET ) < 0 ) error(SYSTEM, "can't seek in wavelet coefficient file " "in getPreCompContrib()" ); /* Read 32-bit encoded wavelet approximation coefficients and mRGBE * range; omit mRGBE minimum if only 1 compressed bin (=maximum * compression), since it's implicitly zero in this case. */ if (getbinary(rgbe32, sizeof(COLR), WAVELET_PADD4_NUMAPPROX + 1 + (nCompressedCoeffs > 1), preCompContrib -> waveletFile ) != WAVELET_PADD4_NUMAPPROX + 1 + (nCompressedCoeffs > 1) ) error(SYSTEM, "can't read approximation coefficients from " "wavelet coefficient file in getPreCompContrib()" ); if (getbinary(preCompContrib -> mrgbeCoeffs, sizeof(mRGBE), nCompressedCoeffs, preCompContrib -> waveletFile ) != nCompressedCoeffs ) error(SYSTEM, "can't read detail coefficients from " "wavelet coefficient file in getPreCompContrib()" ); /* Get mRGBE range (NOTE min/max are reversed in the coeff file) Note that direct assignment isn't possible as colr_color() returns float, not double. */ colr_color(fCoeff, rgbe32 [WAVELET_PADD4_NUMAPPROX]); copycolor(preCompContrib -> mrgbeRange.max, fCoeff); if (nCompressedCoeffs > 1) { colr_color(fCoeff, rgbe32 [WAVELET_PADD4_NUMAPPROX + 1]); copycolor(preCompContrib -> mrgbeRange.min, fCoeff); } else { /* Single compressed bin; mRGBE minimum is implicitly zero */ setcolor(preCompContrib -> mrgbeRange.min, 0, 0, 0); } /* Zero wavelet coefficient matrix and set approximation coefficients in * the upper left of the matrix, as expected by the decodeContribs(). */ waveletMatrix = preCompContrib -> waveletMatrix; zeroCoeffs2(waveletMatrix, coeffDim); for (i = 0; i < WAVELET_PADD4_APPROXDIM; i++) for (j = 0; j < WAVELET_PADD4_APPROXDIM; j++) { coeffIdx = PMAP_CONTRIB_XY2LIN(WAVELET_PADD4_APPROXDIM, i, j); /* Direct assignment to wavelet matrix via colr_color() isn't * possible as the latter returns float, not double */ colr_color(fCoeff, rgbe32 [coeffIdx]); /* HACK: depending on the boundary extension mode, some approx * coeffs may be NEGATIVE (!), so set sign from extra bit in * 32-bit RGBE mantissa. */ for (k = 0; k < 3; k++) waveletMatrix [i] [j] [k] = PMAP_CONTRIB_GET_RGBE32_SGN( rgbe32 [coeffIdx] [k], fCoeff [k] ); } /* All set, now decode mRGBE coeffs and invert wavelet transform */ if (decodeContribs(preCompContrib)) error(INTERNAL, "failed contribution decoding/decompression " "in getPreCompContrib()" ); /* Copy decoded contributions from wavelet coefficient matrix rows */ for (i = 0; i < scDim; i++) { #ifdef PMAP_CONTRIB_LOG for (j = 0; j < scDim; j++) for (k = 0; k < 3; k++) waveletMatrix [i] [j] [k] = PMAP_CONTRIB_IVAL( waveletMatrix [i] [j] [k] ); #endif memcpy(binnedContribs + PMAP_CONTRIB_XY2LIN(scDim, i, 0), waveletMatrix [i], scDim * WAVELET_COEFFSIZE ); } } typedef struct { const RAY *ray; RREAL rayCoeff [3]; COLORV *totalContrib; } PreCompContribRCData; +#if 0 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, k; if (!preCompContribPmap || !rcData || !(preCompContrib = (PreComputedContrib*)( preCompContribPmap -> preCompContrib )) ) { sprintf(errmsg, "missing precomputed contributions for modifier %s", modName ); error(INTERNAL, errmsg); } /* Find rcontrib's LUT entry for this modifier */ contribTabNode = lu_find(pmapContribTab, modName); if (!contribTabNode || !contribTabNode -> key || !(modCont = (MODCONT*)contribTabNode -> data) ) { sprintf(errmsg, "missing contribution LUT entry for modifier %s", modName ); error(INTERNAL, errmsg); } /* 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", modName ); error(SYSTEM, errmsg); } else { sprintf(errmsg, "bin count mismatch for modifier %s, resizing", 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 ) ) { if (preCompContrib -> nBins > 1) { /* BINNING ENABLED; lazily init contrib cache if necessary */ if (initContribCache(preCompContrib)) { /* CACHE ENABLED; fetch cache page for found photon */ DCOLOR *cacheBins; if (!getContribCache(preCompContrib, photonIdx, &cacheBins)) { /* PAGE NOT CACHED; decode contribs into new cache page */ getPreCompContribByPhoton(&photon, photonIdx, preCompContrib, cacheBins ); } else; /* PAGE CACHED; (re)use decoded contribs from cache! */ for (b = 0; b < preCompContrib -> nBins; b++) { /* Accumulate binned contributions into rcontrib's LUT entry for this modifier, weighted by ray coefficient. Also sum total contributions. */ /* NOTE: Using multcolor() on the decoded contribs would modify them in the cache, so use a temp variable here. */ for (k = 0; k < 3; k++) photonContrib [k] = cacheBins [b] [k] * rcData -> rayCoeff [k]; addcolor(modCont -> cbin [b], photonContrib); addcolor(rcData -> totalContrib, photonContrib); } } else { /* CACHE DISABLED; decode contribs into temp array on stack */ DCOLOR tempBins [preCompContrib -> nBins]; getPreCompContribByPhoton(&photon, photonIdx, preCompContrib, tempBins ); for (b = 0; b < preCompContrib -> nBins; b++) { /* Accumulate binned contributions into rcontrib's LUT entry for this modifier, weighted by ray coefficient. Also sum total contributions. */ for (k = 0; k < 3; k++) photonContrib [k] = tempBins [b] [k] * rcData -> rayCoeff [k]; addcolor(modCont -> cbin [b], photonContrib); addcolor(rcData -> totalContrib, photonContrib); } } } else { /* NO BINNING; get single contribution directly from photon, scale by ray coefficient and sum to total contribs */ getPhotonFlux(&photon, photonContrib); multcolor(photonContrib, rcData -> rayCoeff); addcolor(modCont -> cbin [0], photonContrib); addcolor(rcData -> totalContrib, photonContrib); } } return 0; } +#else + +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; + PhotonSearchQueue *squeue; + COLOR photonContrib; + unsigned i, b, k; + + if (!preCompContribPmap || !rcData || + !(preCompContrib = (PreComputedContrib*)( + preCompContribPmap -> preCompContrib + )) + ) { + sprintf(errmsg, "missing precomputed contributions for modifier %s", + modName + ); + error(INTERNAL, errmsg); + } + + /* Find rcontrib's LUT entry for this modifier */ + contribTabNode = lu_find(pmapContribTab, modName); + if (!contribTabNode || !contribTabNode -> key || + !(modCont = (MODCONT*)contribTabNode -> data) + ) { + sprintf(errmsg, "missing contribution LUT entry for modifier %s", + modName + ); + error(INTERNAL, errmsg); + } + + /* 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", modName + ); + error(SYSTEM, errmsg); + } + else { + sprintf(errmsg, "bin count mismatch for modifier %s, resizing", + modName + ); + error(WARNING, errmsg); + } + + modCont = (MODCONT*)contribTabNode -> data; + modCont -> nbins = preCompContrib -> nBins; + memset(modCont -> cbin, 0, sizeof(DCOLOR) * modCont -> nbins); + } + + if (!preCompContribPmap -> maxGather) + /* Zero bandwidth */ + return -1; + + /* Fetch closest photon and its contributions */ + findPhotons(preCompContribPmap, rcData -> ray); + squeue = &preCompContribPmap -> squeue; + for (i = 0; i < squeue -> tail; i++) { + photon = getNearestPhoton(squeue, squeue -> node [i].idx); + + if (preCompContrib -> nBins > 1) { + /* BINNING ENABLED; lazily init contrib cache if necessary */ + #if 0 + if (initContribCache(preCompContrib)) { + /* CACHE ENABLED; fetch cache page for found photon */ + DCOLOR *cacheBins; + + if (!getContribCache(preCompContrib, photonIdx, &cacheBins)) { + /* PAGE NOT CACHED; decode contribs into new cache page */ + getPreCompContribByPhoton(&photon, photonIdx, preCompContrib, + cacheBins + ); + } + else; /* PAGE CACHED; (re)use decoded contribs from cache! */ + + for (b = 0; b < preCompContrib -> nBins; b++) { + /* Accumulate binned contributions into rcontrib's LUT entry + for this modifier, weighted by ray coefficient. Also sum + total contributions. */ + /* NOTE: Using multcolor() on the decoded contribs would + modify them in the cache, so use a temp variable here. */ + for (k = 0; k < 3; k++) + photonContrib [k] = cacheBins [b] [k] * + rcData -> rayCoeff [k]; + + addcolor(modCont -> cbin [b], photonContrib); + addcolor(rcData -> totalContrib, photonContrib); + } + } + else + #endif + { + /* CACHE DISABLED; decode contribs into temp array on stack */ + DCOLOR tempBins [preCompContrib -> nBins]; + + getPreCompContribByPhoton(photon, squeue -> node [i].recIdx, + preCompContrib, tempBins + ); + + for (b = 0; b < preCompContrib -> nBins; b++) { + /* Accumulate binned contributions into rcontrib's LUT entry + for this modifier, weighted by ray coefficient. Also sum + total contributions. */ + for (k = 0; k < 3; k++) + photonContrib [k] = tempBins [b] [k] * + rcData -> rayCoeff [k] / squeue -> tail; + + addcolor(modCont -> cbin [b], photonContrib); + addcolor(rcData -> totalContrib, photonContrib); + } + } + } + else { + /* NO BINNING; get single contribution directly from photon, + scale by ray coefficient and sum to total contribs */ + getPhotonFlux(photon, photonContrib);\ + scalecolor(photonContrib, 1.0 / squeue -> tail); + multcolor(photonContrib, rcData -> rayCoeff); + addcolor(modCont -> cbin [0], photonContrib); + addcolor(rcData -> totalContrib, photonContrib); + } + } + + return 0; +} +#endif + 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 */