diff --git a/oocnn.c b/oocnn.c index c2310ff..e5de84a 100644 --- a/oocnn.c +++ b/oocnn.c @@ -1,476 +1,480 @@ #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, 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; /* Start with undefined max distance */ float d2max = -1; #ifdef PMAP_PHOTONFLOW /* 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; #ifdef PMAP_PHOTONFLOW 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 +#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; memcpy(OOC_GetNearest(squeue, sqn [kid].idx), rec, squeue -> recSize); #ifdef PMAP_PHOTONFLOW 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; #ifdef PMAP_PHOTONFLOW 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; #ifdef PMAP_PHOTONFLOW 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; memcpy(OOC_GetNearest(squeue, sqn[root].idx), rec, squeue->recSize); #ifdef PMAP_PHOTONFLOW 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 +#endif } } if (squeue -> tail >= squeue -> len) /* Search queue full; update max distance^2 from head node */ d2max = sqn [0].dist2; #if defined(PMAP_PHOTONFLOW) && 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; /* 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) /* 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, 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; /* 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, 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)) /* Update closest record and max SQUARED dist to key if it * passes filter */ if ((d2 = dist2(key, oct -> key(rec))) < maxDist2) { memcpy(nnRec, rec, oct -> recSize); maxDist2 = d2; } } } return maxDist2; } #endif /* NIX / PMAP_OOC */ diff --git a/oocnn.h b/oocnn.h index 04ce3e8..7d41c77 100644 --- a/oocnn.h +++ b/oocnn.h @@ -1,120 +1,123 @@ /* ========================================================================= 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_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 the SQUARED distance to furthest neighbour on - * success, else -1 on failure. */ + * 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, float maxDist2 ); /* Find single nearest neighbour within max SQUARED distance maxDist2 * around key and copy corresponding record in nnRec. 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/pmapkdt.c b/pmapkdt.c index 5cd6b96..e535f3d 100644 --- a/pmapkdt.c +++ b/pmapkdt.c @@ -1,776 +1,779 @@ /* ====================================================================== In-core kd-tree for photon map Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) (c) Fraunhofer Institute for Solar Energy Systems, supported by the German Research Foundation (DFG LU-204/10-2, "Fassadenintegrierte Regelsysteme" (FARESYS)) (c) Lucerne University of Applied Sciences and Arts, supported by the Swiss National Science Foundation (SNSF #147053, "Daylight Redirecting Components") (c) Tokyo University of Science, supported by the JSPS Grants-in-Aid for Scientific Research (KAKENHI JP19KK0115, "Three-Dimensional Light Flow") ====================================================================== $Id: pmapkdt.c,v 1.1 2020/11/10 01:10:57 u-no-hoo Exp u-no-hoo $ */ #include "pmapdata.h" /* Includes pmapkdt.h */ #include "pmapfilt.h" #include "source.h" #include "otspecial.h" #include "random.h" /* PHOTON MAP BUILD ROUTINES ------------------------------------------ */ void kdT_Null (PhotonKdTree *kdt) { kdt -> nodes = NULL; } static unsigned long kdT_MedianPartition ( const Photon *heap, unsigned long *heapIdx, unsigned long *heapXdi, unsigned long left, unsigned long right, unsigned dim ) /* Returns index to median in heap from indices left to right (inclusive) in dimension dim. The heap is partitioned relative to median using a quicksort algorithm. The heap indices in heapIdx are sorted rather than the heap itself. */ { const float *p; unsigned long l, r, lg2, n2, m, n = right - left + 1; unsigned d; /* Round down n to nearest power of 2 */ for (lg2 = 0, n2 = n; n2 > 1; n2 >>= 1, ++lg2); n2 = 1 << lg2; /* Determine median position; this takes into account the fact that only the last level in the heap can be partially empty, and that it fills from left to right */ m = left + ((n - n2) > (n2 >> 1) - 1 ? n2 - 1 : n - (n2 >> 1)); while (right > left) { /* Pivot node */ p = heap [heapIdx [right]].pos; l = left; r = right - 1; /* l & r converge, swapping elements out of order with respect to pivot node. Identical keys are resolved by cycling through dim. The convergence point is then the pivot's position. */ do { while (l <= r) { d = dim; while (heap [heapIdx [l]].pos [d] == p [d]) { d = (d + 1) % 3; if (d == dim) { /* Ignore dupes? */ error(WARNING, "duplicate keys in photon heap"); l++; break; } } if (heap [heapIdx [l]].pos [d] < p [d]) l++; else break; } while (r > l) { d = dim; while (heap [heapIdx [r]].pos [d] == p [d]) { d = (d + 1) % 3; if (d == dim) { /* Ignore dupes? */ error(WARNING, "duplicate keys in photon heap"); r--; break; } } if (heap [heapIdx [r]].pos [d] > p [d]) r--; else break; } /* Swap indices (not the nodes they point to) */ n2 = heapIdx [l]; heapIdx [l] = heapIdx [r]; heapIdx [r] = n2; /* Update reverse indices */ heapXdi [heapIdx [l]] = l; heapXdi [n2] = r; } while (l < r); /* Swap indices of convergence and pivot nodes */ heapIdx [r] = heapIdx [l]; heapIdx [l] = heapIdx [right]; heapIdx [right] = n2; /* Update reverse indices */ heapXdi [heapIdx [r]] = r; heapXdi [heapIdx [l]] = l; heapXdi [n2] = right; if (l >= m) right = l - 1; if (l <= m) left = l + 1; } /* Once left & right have converged at m, we have found the median */ return m; } static void kdT_Build ( Photon *heap, unsigned long *heapIdx, unsigned long *heapXdi, const float min [3], const float max [3], unsigned long left, unsigned long right, unsigned long root ) /* Recursive part of balancePhotons(..). Builds heap from subarray defined by indices left and right. min and max are the minimum resp. maximum photon positions in the array. root is the index of the current subtree's root, which corresponds to the median's 1-based index in the heap. heapIdx are the balanced heap indices. The heap is accessed indirectly through these. heapXdi are the reverse indices from the heap to heapIdx so that heapXdi [heapIdx [i]] = i. */ { float maxLeft [3], minRight [3]; Photon rootNode; unsigned d; /* Choose median for dimension with largest spread and partition accordingly */ const float d0 = max [0] - min [0], d1 = max [1] - min [1], d2 = max [2] - min [2]; const unsigned char dim = d0 > d1 ? d0 > d2 ? 0 : 2 : d1 > d2 ? 1 : 2; const unsigned long median = left == right ? left : kdT_MedianPartition(heap, heapIdx, heapXdi, left, right, dim); /* Place median at root of current subtree. This consists of swapping the median and the root nodes and updating the heap indices */ memcpy(&rootNode, heap + heapIdx [median], sizeof(Photon)); memcpy(heap + heapIdx [median], heap + root - 1, sizeof(Photon)); rootNode.discr = dim; memcpy(heap + root - 1, &rootNode, sizeof(Photon)); heapIdx [heapXdi [root - 1]] = heapIdx [median]; heapXdi [heapIdx [median]] = heapXdi [root - 1]; heapIdx [median] = root - 1; heapXdi [root - 1] = median; /* Update bounds for left and right subtrees and recurse on them */ for (d = 0; d <= 2; d++) if (d == dim) maxLeft [d] = minRight [d] = rootNode.pos [d]; else { maxLeft [d] = max [d]; minRight [d] = min [d]; } if (left < median) kdT_Build( heap, heapIdx, heapXdi, min, maxLeft, left, median - 1, root << 1 ); if (right > median) kdT_Build( heap, heapIdx, heapXdi, minRight, max, median + 1, right, (root << 1) + 1 ); } void kdT_BuildPhotonMap (struct PhotonMap *pmap) { Photon *nodes; unsigned long i; unsigned long *heapIdx, /* Photon index array */ *heapXdi; /* Reverse index to heapIdx */ /* Allocate kd-tree nodes and load photons from heap file */ if (!(nodes = calloc(pmap -> numPhotons, sizeof(Photon)))) error(SYSTEM, "failed in-core heap allocation in kdT_BuildPhotonMap"); rewind(pmap -> heap); i = fread(nodes, sizeof(Photon), pmap -> numPhotons, pmap -> heap); if (i != pmap -> numPhotons) error(SYSTEM, "failed loading photon heap in kdT_BuildPhotonMap"); pmap -> store.nodes = nodes; heapIdx = calloc(pmap -> numPhotons, sizeof(unsigned long)); heapXdi = calloc(pmap -> numPhotons, sizeof(unsigned long)); if (!heapIdx || !heapXdi) error(SYSTEM, "failed heap index allocation in kdT_BuildPhotonMap"); /* Initialize index arrays */ for (i = 0; i < pmap -> numPhotons; i++) heapXdi [i] = heapIdx [i] = i; /* Build kd-tree */ kdT_Build( nodes, heapIdx, heapXdi, pmap -> minPos, pmap -> maxPos, 0, pmap -> numPhotons - 1, 1 ); /* Cleanup */ free(heapIdx); free(heapXdi); } /* PHOTON MAP I/O ROUTINES ------------------------------------------ */ int kdT_SavePhotons (const struct PhotonMap *pmap, FILE *out) { unsigned long i, j; Photon *p = (Photon*)pmap -> store.nodes; for (i = 0; i < pmap -> numPhotons; i++, p++) { /* Write photon attributes */ for (j = 0; j < 3; j++) putflt(p -> pos [j], out); /* Bytewise dump otherwise we have portability probs */ for (j = 0; j < 3; j++) putint(p -> norm [j], 1, out); #ifdef PMAP_FLOAT_FLUX for (j = 0; j < 3; j++) putflt(p -> flux [j], out); #else for (j = 0; j < 4; j++) putint(p -> flux [j], 1, out); #endif putint(p -> primary, sizeof(p -> primary), out); putint(p -> flags, 1, out); if (ferror(out)) return -1; } return 0; } int kdT_LoadPhotons (struct PhotonMap *pmap, FILE *in) { unsigned long i, j; Photon *p; /* Allocate kd-tree based on initialised pmap -> numPhotons */ pmap -> store.nodes = calloc(sizeof(Photon), pmap -> numPhotons); if (!pmap -> store.nodes) error(SYSTEM, "failed kd-tree allocation in kdT_LoadPhotons"); /* Get photon attributes */ for (i = 0, p = pmap -> store.nodes; i < pmap -> numPhotons; i++, p++) { for (j = 0; j < 3; j++) p -> pos [j] = getflt(in); /* Bytewise grab otherwise we have portability probs */ for (j = 0; j < 3; j++) p -> norm [j] = getint(1, in); #ifdef PMAP_FLOAT_FLUX for (j = 0; j < 3; j++) p -> flux [j] = getflt(in); #else for (j = 0; j < 4; j++) p -> flux [j] = getint(1, in); #endif p -> primary = getint(sizeof(p -> primary), in); p -> flags = getint(1, in); if (feof(in)) return -1; } return 0; } /* PHOTON MAP SEARCH ROUTINES ------------------------------------------ */ void kdT_InitFindPhotons (struct PhotonMap *pmap) { pmap -> squeue.len = pmap -> maxGather + 1; pmap -> squeue.node = calloc( pmap -> squeue.len, sizeof(PhotonSearchQueueNode) ); if (!pmap -> squeue.node) error(SYSTEM, "can't allocate photon search queue"); #ifdef PMAP_PHOTONFLOW initPhotonPaths(pmap); #endif } #ifdef DEBUG_KDT_NN static int kdT_CheckSearchQueue ( kdT_SearchQueue *squeue, const float pos [3], unsigned root ) /* Priority queue sanity check */ { unsigned kid; const kdT_SearchQueueNode *sqn = squeue -> node; const Photon *photon; float d2, dv [3]; if (root < squeue -> tail) { photon = kdT_GetNearestPhoton(squeue, sqn [root].idx); VSUB(dv, pos, photon -> pos); d2 = DOT(dv, dv); 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 kdT_CheckSearchQueue(squeue, pos, kid); } if (++kid < squeue -> tail) { if (sqn [kid].dist2 > sqn [root].dist2) return -1; else return kdT_CheckSearchQueue(squeue, pos, kid); } } return 0; } #endif static float kdT_PutNearest ( kdT_SearchQueue *squeue, float d2, const Photon *photon, const kdT_SearchAttribCallback *attrib ) /* Insert new photon with SQUARED distance to query point into the search * priority queue, maintaining the most distant photon at the queue head. * If the queue is full, the new photon 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 photon attributes to their queue nodes * by calling attrib->findFunc() if attrib != NULL. It enables finding a * previously enqueued photon and replacing it if the new photon's SQUARED * distance d2 is lower. * * 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. */ { const Photon *tPhoton = NULL; kdT_SearchQueueNode *sqn = squeue -> node, **attribSqn = NULL, **tAttribSqn = NULL; unsigned root, kid, kid2; /* Start with undefined max distance */ - float d2max = -1; + float d2max = -1; #ifdef PMAP_PHOTONFLOW /* Find previously enqueued photon with same attribute if callback defined */ #ifdef PMAP_DEBUGPATHS printf("-------- Enqueue --------\n"); #endif attribSqn = (kdT_SearchQueueNode**)( attrib ? attrib -> findFunc(photon, 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 photon 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; #ifdef PMAP_PHOTONFLOW if (attrib) { /* Update root's attribute map entry to refer to kid */ tPhoton = kdT_GetNearestPhoton(squeue, sqn [root].idx); tAttribSqn = (kdT_SearchQueueNode**)( attrib -> findFunc(tPhoton, attrib -> state) ); #ifdef PMAP_DEBUGPATHS printf("update\n"); #endif if (!tAttribSqn) { /* Consistency error: a previously enqueued photon must * have an entry in the attribute map! */ error( CONSISTENCY, "kdT_PutNearest: corrupt attribute map" ); exit(-1); } *tAttribSqn = &sqn [kid]; } #endif kid = root; } else break; } /* Priority queue restored; assign new photon's queue entry */ - sqn [kid].dist2 = d2; + sqn [kid].dist2 = d2; sqn [kid].idx = (PhotonIdx)photon; #ifdef PMAP_PHOTONFLOW if (attrib) /* Add new photon to attribute map */ attrib -> addFunc(photon, &sqn [kid], attrib -> state); #endif } else { /* Search queue full or duplicate attribute in queue and new photon may be closer; set root to replaced photon, otherwise to queue head */ - root = attribSqn ? *attribSqn - sqn : 0; + root = attribSqn ? *attribSqn - sqn : 0; - if (d2 < sqn [root].dist2) { + if (d2 < sqn [root].dist2) { /* New photon closer than maximum at root; replace and resort towards * leaves (=search queue tail) */ -#ifdef PMAP_PHOTONFLOW +#ifdef PMAP_PHOTONFLOW if (attrib && !attribSqn) { /* New photon will replace root and has unique attribute; delete * root from the attribute map */ tPhoton = kdT_GetNearestPhoton(squeue, sqn [root].idx); attrib -> delFunc(tPhoton, 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; -#ifdef PMAP_PHOTONFLOW +#ifdef PMAP_PHOTONFLOW if (attrib) { /* Update kid's attribute map entry to refer to root */ tPhoton = kdT_GetNearestPhoton(squeue, sqn [kid].idx); tAttribSqn = (kdT_SearchQueueNode**)( attrib -> findFunc(tPhoton, attrib -> state) ); #ifdef PMAP_DEBUGPATHS printf("update\n"); - #endif + #endif if (!tAttribSqn) { /* Consistency error: a previously enqueued photon must - * have an entry in the attribute map! */ + * have an entry in the attribute map! */ error( CONSISTENCY, "kdT_PutNearest: corrupt attribute map" ); exit(-1); } *tAttribSqn = &sqn [root]; } #endif } else break; root = kid; } /* Priority queue restored; reassign head or replaced node */ sqn [root].dist2 = d2; sqn [root].idx = (PhotonIdx)photon; #ifdef PMAP_PHOTONFLOW if (attrib) { if (!attribSqn) /* Add new photon to attribute map */ attrib -> addFunc(photon, &sqn [root], attrib -> state); else { #ifdef PMAP_DEBUGPATHS attrib -> findFunc(photon, attrib -> state); printf("update\n"); #endif /* Update new photon'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_PHOTONFLOW) && defined (PMAP_DEBUGPATHS) if (attrib) /* Run sanity check on attribute map */ attrib -> checkFunc(attrib -> state); #endif - return d2max; + return d2max; } static void kdT_FindNearest ( PhotonMap *pmap, const float pos [3], const kdT_SearchFilterCallback *filter, const kdT_SearchAttribCallback *attrib, unsigned long node ) /* Recursive part of kdT_FindPhotons(). Locate pmap -> squeue.len nearest * neighbours which pass the specified filt (if not NULL) and return in * search queue starting at pmap -> squeue.node. */ { Photon *p = (Photon*)pmap -> store.nodes + node - 1; /* Signed distance to current photon's splitting plane */ float d = pos [p -> discr] - p -> pos [p -> discr], d2 = d * d, dv [3]; /* Search subtree closer to pos first; exclude other subtree if the distance to the splitting plane is greater than maxDist */ if (d < 0) { if (node << 1 <= pmap -> numPhotons) kdT_FindNearest(pmap, pos, filter, attrib, node << 1); if (d2 < pmap -> maxDist2 && node << 1 < pmap -> numPhotons) kdT_FindNearest(pmap, pos, filter, attrib, (node << 1) + 1); } else { if (node << 1 < pmap -> numPhotons) kdT_FindNearest(pmap, pos, filter, attrib, (node << 1) + 1); if (d2 < pmap -> maxDist2 && node << 1 <= pmap -> numPhotons) kdT_FindNearest(pmap, pos, filter, attrib, node << 1); } /* Bail out if photon is rejected by filter */ if (filter && !filter -> filtFunc(p, filter -> state)) return; /* Squared distance to current photon (note dist2() requires doubles) */ VSUB(dv, pos, p -> pos); d2 = DOT(dv, dv); /* Accept photon if closer than current max dist & add to search queue */ if (d2 < pmap -> maxDist2) { /* Insert record in search queue if it passes filter and its SQUARED * dist to key is below maxDist2 */ if ((d2 = kdT_PutNearest(&pmap -> squeue, d2, p, attrib)) >= 0) { /* Update maxDist2 if search queue is full */ pmap -> maxDist2 = d2; #ifdef DEBUG_KDT_NN if (kdT_CheckSearchQueue(&pmap -> squeue, pos, 0)) error(CONSISTENCY, "kdT_PutNearest: corrupt search queue"); #endif } } } int kdT_FindPhotons ( struct PhotonMap *pmap, const FVECT pos, const FVECT norm ) { kdT_SearchFilterCallback filt; kdT_SearchFilterState filtState; #ifdef PMAP_PHOTONFLOW kdT_SearchAttribCallback paths, *pathsPtr = NULL; #endif float p [3], n [3]; /* Photon pos & normal stored at lower precision */ VCOPY(p, pos); if (norm) VCOPY(n, norm); /* Set up filter callback */ filtState.pmap = pmap; filtState.norm = norm ? n : NULL; filt.state = &filtState; filt.filtFunc = filterPhoton; #ifdef PMAP_PHOTONFLOW if (isVolumePmap(pmap)) { /* Set up path ID callback */ paths.state = pmap; paths.findFunc = findPhotonPath; paths.addFunc = addPhotonPath; paths.delFunc = deletePhotonPath; paths.checkFunc = checkPhotonPaths; pathsPtr = &paths; resetPhotonPaths(pmap); } kdT_FindNearest(pmap, p, &filt, pathsPtr, 1); #else kdT_FindNearest(pmap, p, &filt, NULL, 1); #endif - + + /* Get (maximum distance)^2 from search queue head */ + pmap -> maxDist2 = pmap -> squeue.node [0].dist2; + /* Return success or failure (empty queue => none found) */ return pmap -> squeue.tail ? 0 : -1; } static void kdT_Find1Nearest ( PhotonMap *pmap, const float pos [3], const kdT_SearchFilterCallback *filter, Photon **photon, unsigned long node ) /* Recursive part of kdT_Find1Photon(). Locate single nearest neighbour to * pos with similar normal. Note that all heap and queue indices are * 1-based, but accesses to the arrays are 0-based! */ { Photon *p = (Photon*)pmap -> store.nodes + node - 1; /* Signed distance to current photon's splitting plane */ float d = pos [p -> discr] - p -> pos [p -> discr], d2 = d * d, dv [3]; /* Search subtree closer to pos first; exclude other subtree if the distance to the splitting plane is greater than maxDist */ if (d < 0) { if (node << 1 <= pmap -> numPhotons) kdT_Find1Nearest(pmap, pos, filter, photon, node << 1); if (d2 < pmap -> maxDist2 && node << 1 < pmap -> numPhotons) kdT_Find1Nearest(pmap, pos, filter, photon, (node << 1) + 1); } else { if (node << 1 < pmap -> numPhotons) kdT_Find1Nearest(pmap, pos, filter, photon, (node << 1) + 1); if (d2 < pmap -> maxDist2 && node << 1 <= pmap -> numPhotons) kdT_Find1Nearest(pmap, pos, filter, photon, node << 1); } /* Squared distance to current photon */ VSUB(dv, pos, p -> pos); d2 = DOT(dv, dv); if ( d2 < pmap -> maxDist2 && (!filter || filter -> filtFunc(p, filter -> state)) ) { /* Closest photon so far that passes filter */ pmap -> maxDist2 = d2; *photon = p; } } int kdT_Find1Photon ( struct PhotonMap *pmap, const FVECT pos, const FVECT norm, Photon *photon ) { kdT_SearchFilterCallback filt; kdT_SearchFilterState filtState; float p [3], n [3]; /* Photon pos & normal stored at lower precision */ VCOPY(p, pos); if (norm) VCOPY(n, norm); /* Set up filter callback */ filtState.pmap = pmap; filtState.norm = norm ? n : NULL; filt.state = &filtState; filt.filtFunc = filterPhoton; Photon *pnn = NULL; kdT_Find1Nearest(pmap, p, &filt, &pnn, 1); if (!pnn) /* No photon found => failed */ return -1; else { /* Copy found photon => successs */ memcpy(photon, pnn, sizeof(Photon)); return 0; } } /* PHOTON ACCESS ROUTINES ------------------------------------------ */ int kdT_GetPhoton ( const struct PhotonMap *pmap, PhotonIdx idx, Photon *photon ) { memcpy(photon, idx, sizeof(Photon)); return 0; } Photon *kdT_GetNearestPhoton (const PhotonSearchQueue *squeue, PhotonIdx idx) { return idx; } PhotonIdx kdT_FirstPhoton (const struct PhotonMap* pmap) { return pmap -> store.nodes; } void kdT_Delete (PhotonKdTree *kdt) { free(kdt -> nodes); kdt -> nodes = NULL; } diff --git a/pmapooc.c b/pmapooc.c index 155d69c..23079a7 100644 --- a/pmapooc.c +++ b/pmapooc.c @@ -1,340 +1,343 @@ /* ====================================================================== 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.c,v 1.2 2020/11/10 01:08:56 u-no-hoo Exp u-no-hoo $ + $Id: pmapooc.c,v 1.7 2021/03/22 23:00:00 rschregle Exp $ */ #include "pmapdata.h" /* Includes pmapooc.h */ #include "source.h" #include "otspecial.h" #include "oocsort.h" #include "oocbuild.h" /* PHOTON MAP BUILD ROUTINES ------------------------------------------ */ /* Returns photon position as sorting key for OOC_Octree & friends (notably * for Morton code generation). * !!! XXX: Uses type conversion from float via TEMPORARY storage; * !!! THIS IS NOT THREAD SAFE! * !!! RETURNED KEY PERSISTS ONLY IF COPIED BEFORE NEXT CALL! */ RREAL *OOC_PhotonKey (const void *p) { static FVECT photonPos; /* Temp storage for type conversion */ VCOPY(photonPos, ((Photon*)p) -> pos); return photonPos; } #ifdef DEBUG_OOC static void OOC_CheckKeys (FILE *file, const OOC_Octree *oct) { Photon p, lastp; RREAL *key; OOC_MortonIdx zIdx, lastzIdx = 0; rewind(file); memset(&lastp, 0, sizeof(lastp)); while (fread(&p, sizeof(p), 1, file) > 0) { key = OOC_PhotonKey(&p); zIdx = OOC_KEY2MORTON(key, oct); if (zIdx < lastzIdx) error(INTERNAL, "photons not sorted"); if (zIdx == lastzIdx) { sprintf(errmsg, "identical key %021ld " "for [%.9f, %.9f, %.9f]\tand [%.9f, %.9f, %.9f]", zIdx, lastp.pos [0], lastp.pos [1], lastp.pos [2], p.pos [0], p.pos [1], p.pos [2]); error(WARNING, errmsg); } lastzIdx = zIdx; memcpy(&lastp, &p, sizeof(p)); } } #endif void OOC_BuildPhotonMap (struct PhotonMap *pmap, unsigned numProc) { FILE *leafFile; char leafFname [1024]; FVECT d, octOrg; int i; RREAL octSize = 0; /* Determine octree size and origin from pmap extent and init octree */ VCOPY(octOrg, pmap -> minPos); VSUB(d, pmap -> maxPos, pmap -> minPos); for (i = 0; i < 3; i++) if (octSize < d [i]) octSize = d [i]; if (octSize < FTINY) error(INTERNAL, "zero octree size in OOC_BuildPhotonMap"); /* Derive leaf filename from photon map and open file */ strncpy(leafFname, pmap -> fileName, sizeof(leafFname)); strncat(leafFname, PMAP_OOC_LEAFSUFFIX, sizeof(leafFname) - strlen(leafFname) - 1); if (!(leafFile = fopen(leafFname, "w+b"))) error(SYSTEM, "failed to open leaf file in OOC_BuildPhotonMap"); #ifdef DEBUG_OOC eputs("Sorting photons by Morton code...\n"); #endif /* Sort photons in heapfile by Morton code and write to leaf file */ if (OOC_Sort(pmap -> heap, leafFile, PMAP_OOC_NUMBLK, PMAP_OOC_BLKSIZE, numProc, sizeof(Photon), octOrg, octSize, OOC_PhotonKey)) error(INTERNAL, "failed out-of-core photon sort in OOC_BuildPhotonMap"); /* Init and build octree */ OOC_Init(&pmap -> store, sizeof(Photon), octOrg, octSize, OOC_PhotonKey, leafFile); #ifdef DEBUG_OOC eputs("Checking leaf file consistency...\n"); OOC_CheckKeys(leafFile, &pmap -> store); eputs("Building out-of-core octree...\n"); #endif if (!OOC_Build(&pmap -> store, PMAP_OOC_LEAFMAX, PMAP_OOC_MAXDEPTH)) error(INTERNAL, "failed out-of-core octree build in OOC_BuildPhotonMap"); #ifdef DEBUG_OOC eputs("Checking out-of-core octree consistency...\n"); if (OOC_Check(&pmap -> store, OOC_ROOT(&pmap -> store), octOrg, octSize, 0)) error(INTERNAL, "inconsistent out-of-core octree; Time4Harakiri"); #endif } /* PHOTON MAP I/O ROUTINES ------------------------------------------ */ int OOC_SavePhotons (const struct PhotonMap *pmap, FILE *out) { return OOC_SaveOctree(&pmap -> store, out); } int OOC_LoadPhotons (struct PhotonMap *pmap, FILE *nodeFile) { FILE *leafFile; char leafFname [1024]; /* Derive leaf filename from photon map and open file */ strncpy(leafFname, pmap -> fileName, sizeof(leafFname)); strncat(leafFname, PMAP_OOC_LEAFSUFFIX, sizeof(leafFname) - strlen(leafFname) - 1); if (!(leafFile = fopen(leafFname, "r"))) error(SYSTEM, "failed to open leaf file in OOC_LoadPhotons"); if (OOC_LoadOctree(&pmap -> store, nodeFile, OOC_PhotonKey, leafFile)) return -1; #ifdef DEBUG_OOC /* Check octree for consistency */ if (OOC_Check( &pmap -> store, OOC_ROOT(&pmap -> store), pmap -> store.org, pmap -> store.size, 0 )) return -1; #endif - + return 0; } /* PHOTON MAP SEARCH ROUTINES ------------------------------------------ */ void OOC_InitFindPhotons (struct PhotonMap *pmap) { if (OOC_InitNearest( &pmap -> squeue, pmap -> maxGather + 1, sizeof(Photon) )) error(SYSTEM, "can't allocate photon search queue"); #ifdef PMAP_PHOTONFLOW if (isVolumePmap(pmap) && photonFlow) initPhotonPaths(pmap); #endif } static void OOC_InitPhotonCache (struct PhotonMap *pmap) /* Initialise OOC photon cache */ { static char warn = 1; if (!pmap -> store.cache && !pmap -> numDensity) { if (pmapCacheSize > 0) { const unsigned pageSize = pmapCachePageSize * pmap -> maxGather, numPages = pmapCacheSize / pageSize; /* Allocate & init I/O cache in octree */ pmap -> store.cache = malloc(sizeof(OOC_Cache)); if (!pmap -> store.cache || OOC_CacheInit( pmap -> store.cache, numPages, pageSize, sizeof(Photon) )) { error(SYSTEM, "failed OOC photon map cache init"); } } else if (warn) { error(WARNING, "OOC photon map cache DISABLED"); warn = 0; } } } int OOC_FindPhotons ( struct PhotonMap *pmap, const FVECT pos, const FVECT norm ) { OOC_SearchFilterCallback filt; OOC_SearchFilterState filtState; #ifdef PMAP_PHOTONFLOW OOC_SearchAttribCallback paths, *pathsPtr = NULL; #endif - float n [3]; + float res, n [3]; /* Lazily init OOC cache */ if (!pmap -> store.cache) OOC_InitPhotonCache(pmap); /* Set up filter callback */ if (norm) VCOPY(n, norm); filtState.pmap = pmap; filtState.norm = norm ? n : NULL; filt.state = &filtState; filt.filtFunc = filterPhoton; -#ifdef PMAP_PHOTONFLOW +#ifdef _PMAP_PHOTONFLOW if (isVolumePmap(pmap) && photonFlow) { /* Set up path ID callback */ paths.state = pmap; paths.findFunc = findPhotonPath; paths.addFunc = addPhotonPath; paths.delFunc = deletePhotonPath; paths.checkFunc = checkPhotonPaths; - pathsPtr = &paths; + pathsPtr = &paths; resetPhotonPaths(pmap); } - pmap -> maxDist2 = OOC_FindNearest( + res = OOC_FindNearest( &pmap -> store, OOC_ROOT(&pmap -> store), 0, pmap -> store.org, pmap -> store.size, pos, &filt, pathsPtr, &pmap -> squeue, pmap -> maxDist2 ); #else - pmap -> maxDist2 = OOC_FindNearest( + res = OOC_FindNearest( &pmap -> store, OOC_ROOT(&pmap -> store), 0, pmap -> store.org, pmap -> store.size, pos, &filt, NULL, &pmap -> squeue, pmap -> maxDist2 ); #endif /* PMAP_PHOTONFLOW */ - if (pmap -> maxDist2 < 0) + if (res < 0) error(INTERNAL, "failed k-NN photon lookup in OOC_FindPhotons"); - + + /* Get (maximum distance)^2 from search queue head */ + pmap -> maxDist2 = pmap -> squeue.node [0].dist2; + /* Return success or failure (empty queue => none found) */ return pmap -> squeue.tail ? 0 : -1; } int OOC_Find1Photon ( struct PhotonMap* pmap, const FVECT pos, const FVECT norm, Photon *photon ) { OOC_SearchFilterCallback filt; OOC_SearchFilterState filtState; float n [3], maxDist2; /* Lazily init OOC cache */ if (!pmap -> store.cache) OOC_InitPhotonCache(pmap); /* Set up filter callback */ if (norm) VCOPY(n, norm); filtState.pmap = pmap; filtState.norm = norm ? n : NULL; filt.state = &filtState; filt.filtFunc = filterPhoton; maxDist2 = OOC_Find1Nearest( &pmap -> store, OOC_ROOT(&pmap -> store), 0, pmap -> store.org, pmap -> store.size,pos, &filt, photon, pmap -> maxDist2 ); if (maxDist2 < 0) error(INTERNAL, "failed 1-NN photon lookup in OOC_Find1Photon"); if (maxDist2 >= pmap -> maxDist2) /* No photon found => failed */ return -1; else { /* Set photon distance => success */ pmap -> maxDist2 = maxDist2; return 0; } } /* PHOTON ACCESS ROUTINES ------------------------------------------ */ int OOC_GetPhoton (struct PhotonMap *pmap, PhotonIdx idx, Photon *photon) { return OOC_GetData(&pmap -> store, idx, photon); } Photon *OOC_GetNearestPhoton (const PhotonSearchQueue *squeue, PhotonIdx idx) { return OOC_GetNearest(squeue, idx); } PhotonIdx OOC_FirstPhoton (const struct PhotonMap* pmap) { return 0; }