| 7 |
|
EGSR 2011. |
| 8 |
|
|
| 9 |
|
Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) |
| 10 |
< |
(c) Fraunhofer Institute for Solar Energy Systems, |
| 11 |
< |
Lucerne University of Applied Sciences & Arts |
| 10 |
> |
(c) Lucerne University of Applied Sciences and Arts, |
| 11 |
> |
supported by the Swiss National Science Foundation (SNSF, #147053) |
| 12 |
|
======================================================================= |
| 13 |
|
|
| 14 |
|
$Id$ |
| 17 |
|
|
| 18 |
|
#include "oococt.h" |
| 19 |
|
#include "oocsort.h" |
| 20 |
+ |
#include <stdlib.h> |
| 21 |
+ |
#include <string.h> |
| 22 |
|
|
| 23 |
|
|
| 24 |
|
/* Test for empty/full input queue, return pointer to head/tail */ |
| 26 |
|
#define QueueEmpty(q) (!(q) -> len) |
| 27 |
|
#define QueueHead(q) ((q) -> data + (q) -> head * (q) -> recSize) |
| 28 |
|
#define QueueTail(q) ((q) -> data + \ |
| 29 |
< |
(q) -> head + (q) -> len - 1) * (q) -> recSize) |
| 29 |
> |
((q)->head + (q)->len-1) % (q)->cap * (q)->recSize) |
| 30 |
|
|
| 31 |
|
|
| 32 |
< |
/* Input queue from heap */ |
| 33 |
< |
type struct { |
| 34 |
< |
void *data = NULL; |
| 35 |
< |
/* Queue head, length (from head), capacity and record size */ |
| 36 |
< |
unsigned head, len, cap, recSize; |
| 32 |
> |
|
| 33 |
> |
/* Input queue for bottom-up octree construction */ |
| 34 |
> |
typedef struct { |
| 35 |
> |
void *data; |
| 36 |
> |
unsigned head, len, cap, recSize; /* Queue head, length (from head), |
| 37 |
> |
capacity and record size */ |
| 38 |
> |
FILE *in; /* Input file for data records */ |
| 39 |
|
} OOC_BuildQueue; |
| 40 |
|
|
| 41 |
|
|
| 56 |
|
|
| 57 |
|
|
| 58 |
|
static int QueuePush (OOC_BuildQueue *q, const void *rec) |
| 59 |
< |
/* Append record to queue tail; returns new queue length or -1 on failure */ |
| 59 |
> |
/* Append record to queue tail; return new queue length or -1 on failure */ |
| 60 |
|
{ |
| 61 |
< |
int tail; |
| 58 |
< |
|
| 59 |
< |
if (!q || !rec) |
| 61 |
> |
if (!q || !rec || QueueFull(q)) |
| 62 |
|
return -1; |
| 61 |
– |
|
| 62 |
– |
if (q -> len >= q -> cap) |
| 63 |
– |
/* Queue full */ |
| 64 |
– |
return -1; |
| 63 |
|
|
| 64 |
< |
tail = (q -> head + q-> len++) % q -> cap; |
| 65 |
< |
memcpy(q -> data + tail * q -> recSize, rec, q -> recSize); |
| 66 |
< |
|
| 64 |
> |
++q->len; |
| 65 |
> |
memcpy(QueueTail(q), rec, q -> recSize); |
| 66 |
> |
|
| 67 |
|
return q -> len; |
| 68 |
|
} |
| 69 |
|
|
| 70 |
|
|
| 71 |
|
static int QueuePop (OOC_BuildQueue *q, void *rec) |
| 72 |
< |
/* Remove record from queue head; returns new queue length or -1 on failure */ |
| 72 |
> |
/* Remove record from queue head and return in rec if not NULL; return new |
| 73 |
> |
* queue length or -1 on failure */ |
| 74 |
|
{ |
| 75 |
< |
if (!q || !rec) |
| 75 |
> |
if (!q || QueueEmpty(q)) |
| 76 |
|
return -1; |
| 77 |
|
|
| 78 |
< |
if (!q -> len) |
| 79 |
< |
/* Queue empty */ |
| 80 |
< |
return -1; |
| 78 |
> |
/* Return head if rec != NULL */ |
| 79 |
> |
if (rec) |
| 80 |
> |
memcpy(rec, QueueHead(q), q -> recSize); |
| 81 |
|
|
| 83 |
– |
memcpy(rev, q -> data + q -> head * q -> recSize, q -> recSize); |
| 82 |
|
q -> head = (q -> head + 1) % q -> cap; |
| 83 |
|
|
| 84 |
|
return --q -> len; |
| 85 |
|
} |
| 86 |
|
|
| 87 |
|
|
| 88 |
< |
static OOC_Idx *OOC_BuildRecurse (OOC_Octree *oct, OOC_Node *node, |
| 89 |
< |
FVECT org, RREAL size, unsigned depth, |
| 90 |
< |
OOC_Queue *queue, FILE *heapFile, |
| 91 |
< |
FILE *nodeFile, FILE *leafFile) |
| 88 |
> |
static int QueueFill (OOC_BuildQueue *q) |
| 89 |
> |
/* Read records from q -> in until the queue is full; return queue |
| 90 |
> |
* length or -1 on failure */ |
| 91 |
> |
{ |
| 92 |
> |
static void *rec = NULL; |
| 93 |
> |
|
| 94 |
> |
if (!rec && !(rec = malloc(q -> recSize))) |
| 95 |
> |
return -1; |
| 96 |
> |
|
| 97 |
> |
while (!QueueFull(q) && !feof(q -> in) && |
| 98 |
> |
fread(rec, q -> recSize, 1, q -> in)) |
| 99 |
> |
QueuePush(q, rec); |
| 100 |
> |
|
| 101 |
> |
return q -> len; |
| 102 |
> |
} |
| 103 |
> |
|
| 104 |
> |
|
| 105 |
> |
|
| 106 |
> |
static OOC_DataIdx OOC_BuildRecurse (OOC_Octree *oct, OOC_Node* node, |
| 107 |
> |
FVECT org, RREAL size, unsigned depth, |
| 108 |
> |
OOC_BuildQueue *queue) |
| 109 |
|
/* Recursive part of OOC_Build(); insert records from input queue into |
| 110 |
< |
* octree node and subdivide if necessary. Returns number of records in |
| 111 |
< |
* subtree. */ |
| 110 |
> |
* octree node and subdivide into kids if necessary. Returns number of |
| 111 |
> |
* records in subtree or OOC_DATAIDX_ERR if failed. */ |
| 112 |
|
{ |
| 113 |
< |
if (QueueEmpty(queue) || !InNode(QueueHead(queue))) |
| 113 |
> |
int k; |
| 114 |
> |
const RREAL kidSize = size * 0.5; |
| 115 |
> |
|
| 116 |
> |
if (!oct || !node) |
| 117 |
> |
return OOC_DATAIDX_ERR; |
| 118 |
> |
|
| 119 |
> |
if (QueueEmpty(queue) || |
| 120 |
> |
!OOC_InBBox(oct, org, size, oct -> key(QueueHead(queue)))) |
| 121 |
|
/* Input exhausted or queue head outside node */ |
| 122 |
|
return 0; |
| 123 |
+ |
|
| 124 |
+ |
if (QueueFull(queue) && depth < oct -> maxDepth && |
| 125 |
+ |
OOC_InBBox(oct, org, size, oct -> key(QueueTail(queue)))) { |
| 126 |
+ |
/*************************** SUBDIVIDE NODE ************************* |
| 127 |
+ |
* At least leafMax + 1 records (since the queue is full) lie inside |
| 128 |
+ |
* the current node's bounding box, and maxDepth has not been reached |
| 129 |
+ |
* ==> subdivide this node. |
| 130 |
+ |
* (Note it suffices to only test the queue tail against the bounding |
| 131 |
+ |
* box, as the records are in Z-order) |
| 132 |
+ |
********************************************************************/ |
| 133 |
+ |
OOC_Node kid [8]; |
| 134 |
+ |
OOC_DataIdx dataCnt; |
| 135 |
+ |
FVECT kidOrg; |
| 136 |
+ |
|
| 137 |
+ |
#ifdef DEBUG_OOC_BUILD |
| 138 |
+ |
FVECT key; |
| 139 |
+ |
unsigned k2 = 0; |
| 140 |
+ |
#endif |
| 141 |
|
|
| 142 |
< |
if (InNode(QueueTail(queue)) && QueueFull(queue) && |
| 143 |
< |
depth < oct -> maxDepth) { |
| 144 |
< |
|
| 142 |
> |
/* We recurse on the nonempty kids first, then finalise their nodes so |
| 143 |
> |
* they are ordered consecutively, since the parent only indexes the 1st kid */ |
| 144 |
> |
for (k = 0; k < 8; k++) { |
| 145 |
> |
/* Clear kid node and get its octant origin */ |
| 146 |
> |
OOC_CLEARNODE(kid + k); |
| 147 |
> |
|
| 148 |
> |
VCOPY(kidOrg, org); |
| 149 |
> |
OOC_OCTORIGIN(kidOrg, k, kidSize); |
| 150 |
> |
|
| 151 |
> |
/* Recurse on kid and check for early bailout */ |
| 152 |
> |
if (OOC_BuildRecurse(oct, kid + k, kidOrg, kidSize, |
| 153 |
> |
depth + 1, queue) == OOC_DATAIDX_ERR) |
| 154 |
> |
return OOC_DATAIDX_ERR; |
| 155 |
> |
|
| 156 |
> |
#ifdef DEBUG_OOC_BUILD |
| 157 |
> |
if (!QueueEmpty(queue)) { |
| 158 |
> |
VCOPY(key, oct -> key(QueueHead(queue))); |
| 159 |
> |
k2 = OOC_Branch(oct, org, kidSize, key, NULL); |
| 160 |
> |
if (OOC_InBBox(oct, org, size, key) && k2 < k) { |
| 161 |
> |
fprintf(stderr, |
| 162 |
> |
"OOC_BuildRecurse, node subdiv: unsorted key [%f, %f, %f] with " |
| 163 |
> |
"octant %d (last %d with bbox [%f-%f, %f-%f, %f-%f])\n", |
| 164 |
> |
key [0], key [1], key [2], k2, k, |
| 165 |
> |
kidOrg [0], kidOrg [0] + kidSize, kidOrg [1], kidOrg [1] + kidSize, |
| 166 |
> |
kidOrg [2], kidOrg [2] + kidSize); |
| 167 |
> |
}} |
| 168 |
> |
#endif |
| 169 |
> |
} |
| 170 |
> |
|
| 171 |
> |
/* Now finalise consecutive kid nodes, skipping empty ones */ |
| 172 |
> |
for (k = 0; k < 8; k++) |
| 173 |
> |
if ((dataCnt = OOC_DATACNT(kid + k))) { |
| 174 |
> |
/* Nonzero kid ==> allocate and set node */ |
| 175 |
> |
if (!NewNode(oct)) { |
| 176 |
> |
fputs("OOC_BuildRecurse: failed to allocate new node\n", |
| 177 |
> |
stderr); |
| 178 |
> |
return OOC_DATAIDX_ERR; |
| 179 |
> |
} |
| 180 |
> |
OOC_SETROOT(oct, kid + k); |
| 181 |
> |
|
| 182 |
> |
/* Sum kid's data count to parent's and check for overflow */ |
| 183 |
> |
if ((dataCnt += node -> node.num) <= OOC_DATAIDX_MAX) |
| 184 |
> |
node -> node.num = dataCnt; |
| 185 |
> |
else { |
| 186 |
> |
fputs("OOC_BuildRecurse: data count overflow in node\n", |
| 187 |
> |
stderr); |
| 188 |
> |
return OOC_DATAIDX_ERR; |
| 189 |
> |
} |
| 190 |
> |
|
| 191 |
> |
/* Set kid index in parent (if first kid) and corresponding |
| 192 |
> |
* branch bit. The kid is the most recent node and thus at the |
| 193 |
> |
* end of the node array, which coincides with the current |
| 194 |
> |
* subtree root */ |
| 195 |
> |
if (!node -> node.branch) |
| 196 |
> |
node -> node.kid = OOC_ROOTIDX(oct); |
| 197 |
> |
|
| 198 |
> |
OOC_SETBRANCH(node, k); |
| 199 |
> |
} |
| 200 |
|
} |
| 201 |
|
else { |
| 202 |
+ |
/****************************** MAKE LEAF **************************** |
| 203 |
+ |
* Queue contains no more than leafMax records, queue tail lies |
| 204 |
+ |
* outside node's bounding box, or maxDepth reached |
| 205 |
+ |
* ==> make this node a leaf. |
| 206 |
+ |
*********************************************************************/ |
| 207 |
+ |
RREAL *key; |
| 208 |
+ |
|
| 209 |
+ |
#ifdef DEBUG_OOC_BUILD |
| 210 |
+ |
OOC_MortonIdx zIdx, lastzIdx = 0; |
| 211 |
+ |
FVECT /* key, */ |
| 212 |
+ |
lastKey, kidOrg; |
| 213 |
+ |
unsigned lastk = 0; |
| 214 |
+ |
#endif |
| 215 |
+ |
|
| 216 |
+ |
/* Mark as leaf (note it's been cleared by the parent call) */ |
| 217 |
+ |
OOC_SETLEAF(node); |
| 218 |
+ |
|
| 219 |
+ |
while (!QueueEmpty(queue) && |
| 220 |
+ |
OOC_InBBox(oct, org, size, (key = oct->key(QueueHead(queue))))) { |
| 221 |
+ |
/* Record lies inside leaf; increment data counter for octant |
| 222 |
+ |
* containing record. */ |
| 223 |
+ |
|
| 224 |
+ |
if ((k = OOC_Branch(oct, org, kidSize, key, NULL)) < 0) { |
| 225 |
+ |
/* Shouldn't happen, as key tested within bbox above */ |
| 226 |
+ |
fprintf(stderr, "OOC_BuildRecurse: buggered Morton code, " |
| 227 |
+ |
"disruption in space-time continuum?\n"); |
| 228 |
+ |
return OOC_DATAIDX_ERR; |
| 229 |
+ |
} |
| 230 |
+ |
|
| 231 |
+ |
if (node -> leaf.num [k] == OOC_OCTCNT_MAX) { |
| 232 |
+ |
/* Currently we're buggered here; merge records instead? */ |
| 233 |
+ |
fprintf(stderr, "OOC_BuildRecurse: data count overflow in " |
| 234 |
+ |
"leaf: depth = %d, count = %d\n", |
| 235 |
+ |
depth, node -> leaf.num [k]); |
| 236 |
+ |
return OOC_DATAIDX_ERR; |
| 237 |
+ |
} |
| 238 |
+ |
|
| 239 |
+ |
++node -> leaf.num [k]; |
| 240 |
+ |
|
| 241 |
+ |
#ifdef DEBUG_OOC_BUILD |
| 242 |
+ |
/* VCOPY(key, oct -> key(QueueHead(queue))); */ |
| 243 |
+ |
if ((zIdx = OOC_KEY2MORTON(key, oct)) < lastzIdx) { |
| 244 |
+ |
fprintf(stderr, "OOC_BuildRecurse, make leaf: unsorted zIdx %020ld for " |
| 245 |
+ |
"key [%f, %f, %f] (previous zIdx %020ld for " |
| 246 |
+ |
"key [%f, %f, %f]\n", zIdx, key [0], key [1], key [2], |
| 247 |
+ |
lastzIdx, lastKey [0], lastKey [1], lastKey [2]); |
| 248 |
+ |
} |
| 249 |
+ |
|
| 250 |
+ |
VCOPY(kidOrg, org); |
| 251 |
+ |
OOC_OCTORIGIN(kidOrg, k, kidSize); |
| 252 |
+ |
if (k < lastk || zIdx < lastzIdx) { |
| 253 |
+ |
fprintf(stderr, |
| 254 |
+ |
"OOC_BuildRecurse, make leaf: unsorted octant %d (last %d) with " |
| 255 |
+ |
"bbox [%f-%f, %f-%f, %f-%f] for key [%f, %f, %f] with zIdx %020ld " |
| 256 |
+ |
"(last [%f, %f, %f], zIdx %020ld)\n", |
| 257 |
+ |
k, lastk, kidOrg [0], kidOrg [0] + kidSize, |
| 258 |
+ |
kidOrg [1], kidOrg [1] + kidSize, kidOrg [2], kidOrg [2] + kidSize, |
| 259 |
+ |
key [0], key [1], key [2], zIdx, |
| 260 |
+ |
lastKey [0], lastKey [1], lastKey [2], lastzIdx); |
| 261 |
+ |
} |
| 262 |
+ |
lastk = k; |
| 263 |
+ |
lastzIdx = zIdx; |
| 264 |
+ |
VCOPY(lastKey, key); |
| 265 |
+ |
#endif |
| 266 |
+ |
|
| 267 |
+ |
/* Remove record from queue */ |
| 268 |
+ |
QueuePop(queue, NULL); |
| 269 |
+ |
} |
| 270 |
+ |
|
| 271 |
+ |
/* Refill queue for next node(s) */ |
| 272 |
+ |
if (QueueFill(queue) < 0) { |
| 273 |
+ |
fputs("OOC_Build: failed input queue fill\n", stderr); |
| 274 |
+ |
return OOC_DATAIDX_ERR; |
| 275 |
+ |
} |
| 276 |
+ |
} |
| 277 |
|
|
| 278 |
< |
} |
| 278 |
> |
return OOC_DATACNT(node); |
| 279 |
|
} |
| 280 |
+ |
|
| 281 |
|
|
| 282 |
|
|
| 283 |
< |
OOC_Octree *OOC_Build (OOC_Octree *oct, unsigned recSize, unsigned leafMax, |
| 284 |
< |
unsigned maxDepth, const char *heapFileName, |
| 285 |
< |
const char *nodeFileName, const char *leafFileName) |
| 286 |
< |
/* Build out-of-core octree from unsorted records of size recSize in |
| 287 |
< |
* heapFileName, such that leaves contain <= leafMax records, except those |
| 288 |
< |
* at maxDepth. Nodes and leaves are output separately to nodeFileName and leafFileName, |
| 289 |
< |
* respectively. Returns octree pointer on success, else NULL. */ |
| 283 |
> |
OOC_Octree *OOC_Build (OOC_Octree *oct, unsigned leafMax, unsigned maxDepth) |
| 284 |
> |
/* Bottom-up construction of out-of-core octree in postorder traversal. The |
| 285 |
> |
* octree oct is assumed to be initialised with its origin (oct -> org), |
| 286 |
> |
* size (oct -> size), key callback (oct -> key), and its associated leaf |
| 287 |
> |
* file (oct -> leafFile). |
| 288 |
> |
|
| 289 |
> |
* Records are read from the leafFile and assumed to be sorted in Z-order, |
| 290 |
> |
* which defines an octree leaf ordering. Leaves (terminal nodes) are |
| 291 |
> |
* constructed such that they contain <= leafMax records and have a maximum |
| 292 |
> |
* depth of maxDepth. |
| 293 |
> |
|
| 294 |
> |
* Note that the following limits apply: |
| 295 |
> |
* leafMax <= OOC_OCTCNT_MAX (see oococt.h) |
| 296 |
> |
* maxDepth <= OOC_MORTON_BITS (see oocsort.h) |
| 297 |
> |
|
| 298 |
> |
* The maxDepth limit arises from the fact that the Z-ordering has a limited |
| 299 |
> |
* resolution and will map node coordinates beyond a depth of |
| 300 |
> |
* OOC_MORTON_BITS to the same Z-index, causing nodes to be potentially read |
| 301 |
> |
* out of sequence and winding up in the wrong octree nodes. |
| 302 |
> |
|
| 303 |
> |
* On success, the octree pointer oct is returned, with the constructed |
| 304 |
> |
* nodes in oct -> nodes, and the node count in oct -> numNodes. On |
| 305 |
> |
* failure, NULL is returned. */ |
| 306 |
|
{ |
| 307 |
|
OOC_BuildQueue queue; |
| 308 |
< |
void *rec = malloc(recSize); |
| 122 |
< |
FILE *heapFile = fopen(heapFileName, "rb"), |
| 123 |
< |
*nodeFile = fopen(nodeFileName, "wb"), |
| 124 |
< |
*leafFile = fopen(leafFileName, "wb"); |
| 308 |
> |
OOC_Node root; |
| 309 |
|
|
| 310 |
< |
if (OOC_Init(oct, recSize, leafMax, maxDepth) < 0) { |
| 311 |
< |
perror("OOC_Build: failed octree init"); |
| 310 |
> |
if (!oct || !oct -> size) { |
| 311 |
> |
fputs("OOC_Build: octree not initialised", stderr); |
| 312 |
|
return NULL; |
| 313 |
|
} |
| 130 |
– |
|
| 131 |
– |
if (!heapFile) { |
| 132 |
– |
perror("OOC_Build: failed opening heap file"); |
| 133 |
– |
return NULL; |
| 134 |
– |
} |
| 314 |
|
|
| 315 |
< |
if (!nodeFile || !leafFile) { |
| 316 |
< |
perror("OOC_Build: failed opening node/leaf output file"); |
| 315 |
> |
if (!oct -> leafFile) { |
| 316 |
> |
fputs("OOC_Build: empty leaf file", stderr); |
| 317 |
|
return NULL; |
| 318 |
|
} |
| 319 |
+ |
|
| 320 |
+ |
oct -> leafMax = leafMax; |
| 321 |
+ |
oct -> maxDepth = maxDepth; |
| 322 |
+ |
queue.in = oct -> leafFile; |
| 323 |
+ |
rewind(queue.in); |
| 324 |
|
|
| 325 |
< |
/* Init and fill queue from heap file */ |
| 326 |
< |
if (!rec || !QueueInit(&queue, recSize, leafMax + 1)) { |
| 327 |
< |
perror("OOC_Build: failed input queue init"); |
| 325 |
> |
/* Init queue and fill from leaf file */ |
| 326 |
> |
if (!QueueInit(&queue, oct -> recSize, leafMax + 1) || |
| 327 |
> |
QueueFill(&queue) < 0) { |
| 328 |
> |
fputs("OOC_Build: failed input queue init\n", stderr); |
| 329 |
|
return NULL; |
| 330 |
|
} |
| 331 |
< |
|
| 332 |
< |
do { |
| 333 |
< |
if (feof(heapFile) || !fread(rec, recSize, 1, heapFile)) { |
| 334 |
< |
perror("OOC_Build: failed reading from heap file"); |
| 335 |
< |
|
| 151 |
< |
QueuePush(&queue, rec); |
| 152 |
< |
} while (!QueueFull(&queue)); |
| 153 |
< |
|
| 154 |
< |
if (!OOC_BuildRecurse(oct, oct -> root, queue, |
| 155 |
< |
heapFile, nodeFile, leafFile)) |
| 331 |
> |
|
| 332 |
> |
/* Clear octree root and recurse */ |
| 333 |
> |
OOC_CLEARNODE(&root); |
| 334 |
> |
if (OOC_BuildRecurse(oct, &root, oct -> org, oct -> size, 0, &queue) == |
| 335 |
> |
OOC_DATAIDX_ERR) |
| 336 |
|
return NULL; |
| 337 |
< |
|
| 338 |
< |
fclose(heapFile); |
| 339 |
< |
fclose(nodeFile); |
| 340 |
< |
fclose(leafFile); |
| 337 |
> |
|
| 338 |
> |
/* Finalise octree root */ |
| 339 |
> |
if (!NewNode(oct)) { |
| 340 |
> |
fputs("OOC_Build: failed to allocate octree root\n", stderr); |
| 341 |
> |
return NULL; |
| 342 |
> |
} |
| 343 |
> |
OOC_SETROOT(oct, &root); |
| 344 |
> |
/* Passing OOC_ROOT(oct) avoids annoying compiler warnings about (&root) |
| 345 |
> |
* always evaluating to true when calling OOC_DATAIDX() */ |
| 346 |
> |
oct -> numData = OOC_DATACNT(OOC_ROOT(oct)); |
| 347 |
|
|
| 348 |
|
return oct; |
| 349 |
|
} |
