Homec4science
Diffusion RADIANCE Photon Map (pmapmerge) Browse

RadiancePmap/pmapkdt.c
9bc999be9d6dpmapmerge

Last Change
Raw File

pmapkdt.c
View Options

/*
======================================================================
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
putbinary(&p -> aux, 1, sizeof(p -> aux), 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
getbinary(&p -> aux, 1, sizeof(p -> aux), 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_PATHFILT
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
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_PATHFILT
/* 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_PATHFILT
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_PATHFILT
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_PATHFILT
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_PATHFILT
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_PATHFILT
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_PATHFILT) && 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_PATHFILT
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.pos = pos;
filtState.norm = norm ? n : NULL;
filt.state = &filtState;
filt.filtFunc = filterPhoton;
/* No contribs with kd-tree */
filtState.srcMod = NULL;
#ifdef PMAP_PATHFILT
/* Set up path ID callback to prune photons from duplicate paths */
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 /* PMAP_PATHFILT */
/* 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;
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;
}
c4science · Help