diff --git a/pmapkdt.c b/pmapkdt.c index 5595468..5595b95 100644 --- a/pmapkdt.c +++ b/pmapkdt.c @@ -1,782 +1,782 @@ /* ====================================================================== 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", 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: 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 -> org, sizeof(p -> org), 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 -> org = getint(sizeof(p -> org), 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; #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].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; if (d2 < sqn [root].dist2) { /* New photon closer than maximum at root; replace and resort towards * leaves (=search queue tail) */ #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 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 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 [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; } 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; /* Get sought light source modifier from pmap -> lastContribOrg if precomputing contribution photon */ filtState.srcMod = isContribPmap(pmap) ? - findmaterial(source [ray -> rsrc].so) : NULL; + findmaterial(source [pmap -> lastContribOrg.srcIdx].so) : NULL; #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 /* 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; filtState.srcMod = 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; }