Homec4science
Diffusion RADIANCE Photon Map (pmapmerge) Browse

RadiancePmap/pmaptkdt.c
7449a2b8e70epmapmerge

Last Change
Raw File

pmaptkdt.c
View Options

/*
======================================================================
In-core kd-tree for 4D TRANSIENT photon map. Based on pmapkdt.h with
extensions to handle photon path length (representing time of flight).
Roland Schregle (roland.schregle@{hslu.ch, gmail.com})
(c) Tokyo University of Science,
supported by the JSPS Grants-in-Aid for Scientific Research
(KAKENHI JP19KK0115, "Three-Dimensional Light Flow")
======================================================================
$Id$
*/
#include "pmapdata.h" /* Includes pmaptkdt.h */
#include "source.h"
#include "otspecial.h"
#include "random.h"
/* PHOTON MAP BUILD ROUTINES ------------------------------------------ */
static unsigned long kdT_TransMedianPartition (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 [0..3], where dim=3 is path length
(corresponding to time of flight). 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;
float pLen;
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;
pLen = heap [heapIdx [right]].aux.pathLen;
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 (d < 3
? heap [heapIdx [l]].pos [d] == p [d]
: heap [heapIdx [l]].aux.pathLen == pLen
) {
d = (d + 1) % 4;
if (d == dim) {
/* Ignore dupes? */
error(WARNING, "duplicate keys in photon heap");
l++;
break;
}
}
if (d < 3
? heap [heapIdx [l]].pos [d] < p [d]
: heap [heapIdx [l]].aux.pathLen < pLen
)
l++;
else break;
}
while (r > l) {
d = dim;
while (d < 3
? heap [heapIdx [r]].pos [d] == p [d]
: heap [heapIdx [r]].aux.pathLen == pLen
) {
d = (d + 1) % 4;
if (d == dim) {
/* Ignore dupes? */
error(WARNING, "duplicate keys in photon heap");
r--;
break;
}
}
if (d < 3
? heap [heapIdx [r]].pos [d] > p [d]
: heap [heapIdx [r]].aux.pathLen > pLen
)
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_TransBuild (Photon *heap,
unsigned long *heapIdx, unsigned long *heapXdi,
const float min [4], const float max [4],
unsigned long left, unsigned long right, unsigned long root
)
/* Recursive part of transBalancePhotons(..). Builds heap from subarray
defined by indices left and right. min and max are the minimum resp.
maximum photon positions and path lengths (corresponding to the photons'
time of flight) 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 [4], minRight [4], maxDist, dDist;
Photon rootNode;
unsigned d, dim;
unsigned long median;
/* Choose median for dimension with largest spread and partition
accordingly */
for (maxDist = 0, dim = 0, d = 0; d < 4; d++)
if ((dDist = max [d] - min [d]) > maxDist) {
maxDist = dDist;
dim = d;
}
median = left == right
? left
: kdT_TransMedianPartition(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 < 4; d++)
if (d == dim)
maxLeft [d] = minRight [d] = d < 3
? rootNode.pos [d] : rootNode.aux.pathLen;
else {
maxLeft [d] = max [d];
minRight [d] = min [d];
}
if (left < median)
kdT_TransBuild(heap, heapIdx, heapXdi, min, maxLeft,
left, median - 1, root << 1
);
if (right > median)
kdT_TransBuild(heap, heapIdx, heapXdi, minRight, max,
median + 1, right, (root << 1) + 1
);
}
void kdT_TransBuildPhotonMap (struct PhotonMap *pmap)
{
Photon *nodes;
unsigned long i;
unsigned long *heapIdx, /* Photon index array */
*heapXdi; /* Reverse index to heapIdx */
float min [4], max [4];
/* 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_TransBuildPhotonMap"
);
rewind(pmap -> heap);
i = fread(nodes, sizeof(Photon), pmap -> numPhotons, pmap -> heap);
if (i != pmap -> numPhotons)
error(SYSTEM,
"failed loading photon heap in kdT_TransBuildPhotonMap"
);
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_TransBuildPhotonMap"
);
/* Initialize index arrays */
for (i = 0; i < pmap -> numPhotons; i++)
heapXdi [i] = heapIdx [i] = i;
/* Build kd-tree */
VCOPY(min, pmap -> minPos);
min [3] = pmap -> minPathLen;
VCOPY(max, pmap -> maxPos);
max [3] = pmap -> maxPathLen;
kdT_TransBuild(nodes, heapIdx, heapXdi, min, max,
0, pmap -> numPhotons - 1, 1
);
/* Cleanup */
free(heapIdx);
free(heapXdi);
}
/* PHOTON MAP SEARCH ROUTINES ------------------------------------------ */
static void kdT_TransFindNearest (PhotonMap *pmap,
const float pos [4], const kdT_SearchFilterCallback *filter,
const kdT_SearchAttribCallback *attrib, unsigned long node
)
/* Recursive part of kdT_TransFindPhotons(). Locate pmap -> squeue.len
* nearest neighbours to 4D pos 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 -> discr < 3
? p -> pos [p -> discr]
: p -> aux.pathLen
),
d2 = d * d, dv [4];
/* 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_TransFindNearest(pmap, pos, filter, attrib, node << 1);
if (d2 < pmap -> maxDist2 && node << 1 < pmap -> numPhotons)
kdT_TransFindNearest(pmap, pos, filter, attrib, (node << 1) + 1);
}
else {
if (node << 1 < pmap -> numPhotons)
kdT_TransFindNearest(pmap, pos, filter, attrib, (node << 1) + 1);
if (d2 < pmap -> maxDist2 && node << 1 <= pmap -> numPhotons)
kdT_TransFindNearest(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 (dist2() would require doubles) */
VSUB(dv, pos, p -> pos);
d2 = DOT(dv, dv);
dv [3] = pos [3] - p -> aux.pathLen;
d2 += dv [3] * dv [3];
/* 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;
}
}
int kdT_TransFindPhotons (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 [4], n [3];
/* 4th key dimension: photon path length at given time in units of
* scene geometry for consistency with pos */
const float pathLen = pmap -> time * pmap -> velocity;
/* Photon pos & normal stored at lower precision; note path length is
* appended to position as 4th key dimension */
VCOPY(p, pos);
p [3] = pathLen;
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;
#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_TransFindNearest(pmap, p, &filt, pathsPtr, 1);
#else
kdT_TransFindNearest(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_TransFind1Nearest (PhotonMap *pmap,
const float pos [4], const kdT_SearchFilterCallback *filter,
Photon **photon, unsigned long node
)
/* Recursive part of kdT_Find1Photon(). Locate single nearest neighbour to
* 4D 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 -> discr < 3
? p -> pos [p -> discr]
: p -> aux.pathLen
),
d2 = d * d, dv [4];
/* 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_TransFind1Nearest(pmap, pos, filter, photon, node << 1);
if (d2 < pmap -> maxDist2 && node << 1 < pmap -> numPhotons)
kdT_TransFind1Nearest(pmap, pos, filter, photon, (node << 1) + 1);
}
else {
if (node << 1 < pmap -> numPhotons)
kdT_TransFind1Nearest(pmap, pos, filter, photon, (node << 1) + 1);
if (d2 < pmap -> maxDist2 && node << 1 <= pmap -> numPhotons)
kdT_TransFind1Nearest(pmap, pos, filter, photon, node << 1);
}
/* Squared distance to current photon (dist2() would require doubles) */
VSUB(dv, pos, p -> pos);
d2 = DOT(dv, dv);
dv [3] = pos [3] - p -> aux.pathLen;
d2 += dv [3] * dv [3];
if (d2 < pmap -> maxDist2 &&
(!filter || filter -> filtFunc(p, filter -> state))
) {
/* Closest photon so far that passes filter */
pmap -> maxDist2 = d2;
*photon = p;
}
}
int kdT_Find1TransPhoton (struct PhotonMap *pmap,
const FVECT pos, const FVECT norm, Photon *photon
)
{
kdT_SearchFilterCallback filt;
kdT_SearchFilterState filtState;
float p [4], n [3];
Photon *pnn = NULL;
/* 4th key dimension: photon path length at given time in units of
* scene geometry for consistency with pos */
const float pathLen = pmap -> time * pmap -> velocity;
/* Photon pos & normal stored at lower precision; note path length is
* appended to position as 4th key dimension */
VCOPY(p, pos);
p [3] = pathLen;
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;
kdT_TransFind1Nearest(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;
}
}
c4science · Help