Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F62647632
oocnn.c
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Tue, May 14, 14:21
Size
15 KB
Mime Type
text/x-c
Expires
Thu, May 16, 14:21 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
17672538
Attached To
R10977 RADIANCE Photon Map
oocnn.c
View Options
#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 <stdlib.h>
#include <string.h>
#include <math.h>
/* ========================= 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_PATHFILT
/* Find previously enqueued record with same attribute if callback
defined */
#ifdef PMAP_DEBUGPATHS
printf("-------- Enqueue --------\n");
#endif
attribSqn = (OOC_SearchQueueNode**)(
attrib ? attrib -> findFunc(rec, attrib -> state) : NULL
);
#ifdef PMAP_DEBUGPATHS
printf(attribSqn ? "found\n" : "not found\n");
#endif
#endif
if (!attribSqn && squeue -> tail < squeue -> len) {
/* No duplicate attribute in queue, and queue not full;
insert record at tail and resort towards head */
kid = squeue -> tail++;
while (kid) {
root = (kid - 1) >> 1;
/* Compare with parent and swap if necessary (bubble up),
else terminate */
if (d2 > sqn [root].dist2) {
sqn [kid].dist2 = sqn [root].dist2;
sqn [kid].idx = sqn [root].idx;
#ifdef PMAP_PATHFILT
if (attrib) {
/* Update root's attribute map entry to refer to kid */
tRec = OOC_GetNearest(squeue, sqn [root].idx);
tAttribSqn = (OOC_SearchQueueNode**)(
attrib -> findFunc(tRec, attrib -> state)
);
#ifdef PMAP_DEBUGPATHS
printf("update\n");
#endif
if (!tAttribSqn) {
/* Consistency error: a previously enqueued record must
* have an entry in the attribute map! */
fputs("OOC_PutNearest: corrupt attribute map\n", stderr);
exit(-1);
}
*tAttribSqn = &sqn [kid];
}
#endif
kid = root;
}
else break;
}
/* Priority queue restored; assign tail position as linear index into
* record buffa squeue -> nnRec and append record */
sqn [kid].dist2 = d2;
sqn [kid].idx = squeue -> tail - 1;
memcpy(OOC_GetNearest(squeue, sqn [kid].idx), rec, squeue -> recSize);
#ifdef PMAP_PATHFILT
if (attrib)
/* Add new record to attribute map */
attrib -> addFunc(rec, &sqn [kid], attrib -> state);
#endif
}
else {
/* Search queue full or duplicate attribute in queue and new node may
be closer; set root to replaced node, otherwise to queue head */
root = attribSqn ? *attribSqn - sqn : 0;
if (d2 < sqn [root].dist2) {
/* New record closer than maximum at root; replace and resort
towards leaves (=search queue tail) */
rootIdx = sqn [root].idx;
#ifdef PMAP_PATHFILT
if (attrib && !attribSqn) {
/* New record will replace root and has unique attribute; delete
* root from the attribute map */
tRec = OOC_GetNearest(squeue, sqn [root].idx);
attrib -> delFunc(tRec, attrib -> state);
}
#endif
while ((kid = (root << 1) + 1) < squeue -> tail) {
/* Compare with larger kid & swap if necessary (bubble down),
else terminate */
if (
(kid2 = (kid + 1)) < squeue -> tail &&
sqn [kid2].dist2 > sqn [kid].dist2
)
kid = kid2;
if (d2 < sqn [kid].dist2) {
sqn [root].dist2 = sqn [kid].dist2;
sqn [root].idx = sqn [kid].idx;
#ifdef PMAP_PATHFILT
if (attrib) {
/* Update kid's attribute map entry to refer to root */
tRec = OOC_GetNearest(squeue, sqn [kid].idx);
tAttribSqn = (OOC_SearchQueueNode**)(
attrib -> findFunc(tRec, attrib -> state)
);
#ifdef PMAP_DEBUGPATHS
printf("update\n");
#endif
if (!tAttribSqn) {
/* Consistency error: a previously enqueued record must
* have an entry in the attribute map! */
fputs("OOC_PutNearest: corrupt attribute map\n", stderr);
exit(-1);
}
*tAttribSqn = &sqn [root];
}
#endif
}
else break;
root = kid;
}
/* Reassign head's / replaced node's previous buffa index and
overwrite corresponding record */
sqn [root].dist2 = d2;
sqn [root].idx = rootIdx;
memcpy(OOC_GetNearest(squeue, sqn[root].idx), rec, squeue->recSize);
#ifdef PMAP_PATHFILT
if (attrib) {
if (!attribSqn)
/* Add new record to attribute map */
attrib -> addFunc(rec, &sqn [root], attrib -> state);
else {
#ifdef PMAP_DEBUGPATHS
attrib -> findFunc(rec, attrib -> state);
printf("update\n");
#endif
/* Update new record's existing attribute map entry */
*attribSqn = &sqn [root];
}
}
#endif
}
}
if (squeue -> tail >= squeue -> len)
/* Search queue full; update max distance^2 from head node */
d2max = sqn [0].dist2;
#if defined(PMAP_PATHFILT) && defined(PMAP_DEBUGPATHS)
if (attrib)
/* Run sanity check on attribute map */
attrib -> checkFunc(attrib -> state);
#endif
return d2max;
}
int OOC_InitNearest (OOC_SearchQueue *squeue, unsigned len, unsigned recSize)
{
squeue -> len = len;
squeue -> recSize = recSize;
squeue -> node = calloc(len, sizeof(OOC_SearchQueueNode));
squeue -> nnRec = calloc(len, recSize);
if (!squeue -> node || !squeue -> nnRec) {
perror("OOC_InitNearest: failed search queue allocation");
return -1;
}
return 0;
}
void *OOC_GetNearest (const OOC_SearchQueue *squeue, unsigned idx)
{
return squeue -> nnRec + idx * squeue -> recSize;
}
/* ============================ OCTREE STUFF ============================ */
static float OOC_BBoxDist2 (const FVECT bbOrg, RREAL bbSiz, const FVECT key)
/* Return minimum *SQUARED* distance between key and bounding box defined by
* bbOrg and bbSiz; a distance of 0 implies the key lies INSIDE the bbox */
{
/* Explicit comparison with bbox corners */
int i;
FVECT d;
for (i = 0; i < 3; i++) {
d [i] = key [i] - bbOrg [i];
d [i] = d [i] < 0 ? -d [i] : d [i] - bbSiz;
/* Set to 0 if inside bbox */
if (d [i] < 0)
d [i] = 0;
}
return DOT(d, d);
}
float OOC_FindNearest (
OOC_Octree *oct, OOC_Node *node, OOC_DataIdx dataIdx,
const FVECT org, float size, const FVECT key,
const OOC_SearchFilterCallback *filter,
const OOC_SearchAttribCallback *attrib,
OOC_SearchQueue *squeue, float maxDist2
)
{
const float kidSize = size * 0.5;
unsigned i, kid, kid0;
float d2;
char rec [oct -> recSize];
FVECT kidOrg;
OOC_DataIdx kidDataIdx, recIdx;
OOC_Node *kidNode;
if (!node)
/* No current node -- shouldn't happen */
return -1;
/* Start with suboctant closest to key */
for (kid0 = 0, i = 0; i < 3; i++)
kid0 |= (key [i] > org [i] + kidSize) << i;
for (i = 0; i < 8; i++) {
kid = kid0 ^ i;
kidNode = node;
kidDataIdx = dataIdx + OOC_GetKid(oct, &kidNode, kid);
/* Prune empty suboctant */
if ((!kidNode && !OOC_ISLEAF(node)) ||
(OOC_ISLEAF(node) && !node -> leaf.num [kid])
)
continue;
/* Set up suboctant */
VCOPY(kidOrg, org);
OOC_OCTORIGIN(kidOrg, kid, kidSize);
/* Prune suboctant if not overlapped by maxDist2 */
if (OOC_BBoxDist2(kidOrg, kidSize, key) > maxDist2)
continue;
if (kidNode) {
/* Internal node; recurse into non-empty suboctant */
maxDist2 = OOC_FindNearest(oct, kidNode, kidDataIdx, kidOrg,
kidSize, key, filter, attrib, squeue, maxDist2
);
if (maxDist2 < 0)
/* Bail out on error */
break;
}
else if (OOC_ISLEAF(node)) {
/* Leaf node; do linear check of all records in suboctant */
for (recIdx = kidDataIdx;
recIdx < kidDataIdx + node -> leaf.num [kid];
recIdx++
) {
if (OOC_GetData(oct, recIdx, rec))
return -1;
if (!filter || filter -> filtFunc(rec, filter -> state)) {
/* Insert record in search queue if it passes filter
* and its SQUARED dist to key is below maxDist2 */
if ((d2 = dist2(key, oct -> key(rec))) < maxDist2) {
if ((d2 = OOC_PutNearest(squeue, d2, rec, attrib)) >= 0)
/* Update maxDist2 if search queue is full */
maxDist2 = d2;
#ifdef DEBUG_OOC_NN
if (OOC_CheckSearchQueue(squeue, key, oct->key, 0)) {
fprintf(stderr,
"OOC_PutNearest: corrupt search queue\n"
);
return -1;
}
#endif
}
}
}
}
}
/* XXX: Note this will not correspond to the maximum distance in a
partially filled search queue, since it hasn't been finalised;
this distance should therefore be obtained from the search queue
head by the caller when this func returns! */
return maxDist2;
}
float OOC_Find1Nearest (
OOC_Octree *oct, OOC_Node *node, OOC_DataIdx dataIdx,
const FVECT org, float size, const FVECT key,
const OOC_SearchFilterCallback *filter, void *nnRec,
OOC_DataIdx *nnIdx, float maxDist2
)
{
const float kidSize = size * 0.5;
unsigned i, kid, kid0;
float d2;
char rec [oct -> recSize];
FVECT kidOrg;
OOC_DataIdx kidDataIdx, recIdx;
OOC_Node *kidNode;
if (!node)
/* No current node -- shouldn't happen */
return -1;
/* Start with suboctant closest to key */
for (kid0 = 0, i = 0; i < 3; i++)
kid0 |= (key [i] > org [i] + kidSize) << i;
for (i = 0; i < 8; i++) {
kid = kid0 ^ i;
kidNode = node;
kidDataIdx = dataIdx + OOC_GetKid(oct, &kidNode, kid);
/* Prune empty suboctant */
if ((!kidNode && !OOC_ISLEAF(node)) ||
(OOC_ISLEAF(node) && !node -> leaf.num [kid])
)
continue;
/* Set up suboctant */
VCOPY(kidOrg, org);
OOC_OCTORIGIN(kidOrg, kid, kidSize);
/* Prune suboctant if not overlapped by maxDist2 */
if (OOC_BBoxDist2(kidOrg, kidSize, key) > maxDist2)
continue;
if (kidNode) {
/* Internal node; recurse into non-empty suboctant */
maxDist2 = OOC_Find1Nearest(oct, kidNode, kidDataIdx, kidOrg,
kidSize, key, filter, nnRec, nnIdx, maxDist2
);
if (maxDist2 < 0)
/* Bail out on error */
break;
}
else if (OOC_ISLEAF(node))
/* Leaf node; do linear check of all records in suboctant */
for (recIdx = kidDataIdx;
recIdx < kidDataIdx + node -> leaf.num [kid];
recIdx++
) {
if (OOC_GetData(oct, recIdx, rec))
return -1;
if (!filter || filter -> filtFunc(rec, filter -> state))
/* If record passes filter, update closest record,
its index, and max SQUARED distance to its key */
if ((d2 = dist2(key, oct -> key(rec))) < maxDist2) {
memcpy(nnRec, rec, oct -> recSize);
maxDist2 = d2;
if (nnIdx)
*nnIdx = recIdx;
}
}
}
return maxDist2;
}
#endif /* NIX / PMAP_OOC */
Event Timeline
Log In to Comment