source: GTP/trunk/App/Demos/Illum/pathmap/KDTree.cpp @ 2197

Revision 2197, 29.3 KB checked in by szirmay, 18 years ago (diff)
RevLine 
[2197]1// Graphics Programming Methods
2// An Effective kd-tree Implementation
3// László Szécsi
4// 2003.
5
6// implementation of kd-tree methods
7
8#include "dxstdafx.h"
9#include "KDTree.hpp"
10#include <iostream>
11
12bool KDTree::traverse (const Ray& ray, HitRec& hitRec, float minT, float maxT)
13{
14        // handle dirs parallel to axis here in advance
15        // will not need to check later for every node;
16        Vector dirInverts;
17        for(int i=0; i<3; i++)
18                dirInverts[i] = 1.0f / ray.dir[i];
19
20//      hitRec.object   = NULL;
21        hitRec.isIntersect = false;
22        hitRec.depth = FLT_MAX;
23
24        float rayMin = minT;
25        float rayMax = maxT;
26
27        unsigned int tNode = 0;
28        signed int stackPointer = -1;
29
30        unsigned int leftChild;
31        unsigned int rightChild;
32
33        for(;;)
34        {
35                bool isLeaf = !followChildren(tNode, leftChild, rightChild);
36
37                if(isLeaf) // node is a leaf
38                {
39                        Intersectable** leafList = ((Intersectable***)nodeTable)[tNode];
40                        if(leafList) // leaf contains some objects
41                        {
42                                // first 4 bytes of leafList contains the length
43                                unsigned int leafListLength = (unsigned int)*leafList;
44                                Intersectable** leafListEnd = ++leafList + leafListLength;
45                                for (Intersectable** objectIterator = leafList; objectIterator < leafListEnd; objectIterator++) {
46                                        if((*objectIterator) == forbidden) continue;
47                                        float depth;
48                                        // check if this test was already done
49                                        if (ray.id != (*objectIterator)->lastTestedRayId)
50                                                (*objectIterator)->intersect (ray, depth, minT, maxT);
51                                        else
52                                                depth = (*objectIterator)->lastTestedRayResult.depth;
53                                        // intersection is acceptable only if it is within the node volume
54                                        if ((*objectIterator)->lastTestedRayResult.isIntersect &&
55                                                        depth < hitRec.depth &&
56                                                        depth > rayMin - epsilon &&
57                                                        depth < rayMax + epsilon)
58                                                hitRec = (*objectIterator)->lastTestedRayResult;
59                                }
60
61                                if(hitRec.isIntersect)
62                                        return true;
63                                else
64                                        if(stackPointer >= 0)
65                                        {
66                                                //have not found anything on this branch, pop the other one from stack
67                                                rayMin = traverseStack[stackPointer].rayMin;
68                                                rayMax = traverseStack[stackPointer].rayMax;
69                                                tNode = traverseStack[stackPointer].tNode;
70                                                stackPointer--;
71                                        }
72                                        else
73                                                return false;
74                        }
75                        else
76                        {
77                                //leaf is empty
78                                if(stackPointer >= 0)
79                                {
80                                        rayMin = traverseStack[stackPointer].rayMin;
81                                        rayMax = traverseStack[stackPointer].rayMax;
82                                        tNode = traverseStack[stackPointer].tNode;
83                                        stackPointer--;
84                                }
85                                else
86                                        return false;
87                        }
88                }
89                else //not a leaf, go deeper
90                {
91                        // gather cutting plane information from the nodeTable
92                        float nodeValue = nodeTable[tNode];
93                        unsigned char axis = *(unsigned int*)(&nodeValue) & 0x00000003;
94                        unsigned int bitfield = *(unsigned int*)(&nodeValue) & 0xfffffffc;
95                        float cutPlane = *(float*)(&bitfield);
96
97                        float origCoord = ray.origin[axis];
98                        float dirCoord = dirInverts[axis];
99
100            // calculate cutting plane - ray intersection distance 't'
101                        register float t = (cutPlane - origCoord) * dirCoord;
102
103                        // find the child volume nearer to the origin
104                        unsigned int nearNode = rightChild;
105                        unsigned int farNode = leftChild;
106                        if(origCoord + epsilon < cutPlane)
107                        {
108                                nearNode = leftChild;
109                                farNode = rightChild;
110                        }
111                        else
112                                if( origCoord < cutPlane + epsilon )            // origin is on the plane
113                                        if(     dirCoord < 0 )                                          // direction decides
114                                        {
115                                                tNode = leftChild; continue;
116                                        }
117                                        else
118                                        {
119                                                tNode = rightChild; continue;
120                                        }
121
122                        if(t <= 0 || rayMax + epsilon < t )
123                                // whole interval on near cell
124                                tNode = nearNode;
125                        else
126                        {
127                                if(t < rayMin - epsilon)
128                                        // whole interval on far cell
129                                        tNode = farNode;
130                                else
131                                {
132                                        // both cells
133                                        // push far branch to stack
134                                        stackPointer++;
135                                        traverseStack[stackPointer].rayMin = t;
136                                        traverseStack[stackPointer].rayMax = rayMax;
137                                        traverseStack[stackPointer].tNode = farNode;
138                                        // near branch is next to traverse
139                                        tNode = nearNode;
140                                        rayMax = t;
141                                }
142                        }
143                }
144        }
145        return false;
146}
147
148bool KDTree::traverseBackSide (const Ray& ray, HitRec& hitRec, float minT, float maxT)
149{
150        // handle dirs parallel to axis here in advance
151        // will not need to check later for every node;
152        Vector dirInverts;
153        for(int i=0; i<3; i++)
154                dirInverts[i] = 1.0f / ray.dir[i];
155
156//      hitRec.object   = NULL;
157        hitRec.isIntersect = false;
158        hitRec.depth = FLT_MAX;
159
160        float rayMin = minT;
161        float rayMax = maxT;
162
163        unsigned int tNode = 0;
164        signed int stackPointer = -1;
165
166        unsigned int leftChild;
167        unsigned int rightChild;
168
169        for(;;)
170        {
171                bool isLeaf = !followChildren(tNode, leftChild, rightChild);
172
173                if(isLeaf) // node is a leaf
174                {
175                        Intersectable** leafList = ((Intersectable***)nodeTable)[tNode];
176                        if(leafList) // leaf contains some objects
177                        {
178                                // first 4 bytes of leafList contains the length
179                                unsigned int leafListLength = (unsigned int)*leafList;
180                                Intersectable** leafListEnd = ++leafList + leafListLength;
181                                for (Intersectable** objectIterator = leafList; objectIterator < leafListEnd; objectIterator++) {
182//                                      if((*objectIterator) == forbidden) continue;
183                                        float depth;
184                                        // check if this test was already done
185                                        if (ray.id != (*objectIterator)->lastTestedRayId)
186                                                (*objectIterator)->intersectBackSide (ray, depth, minT, maxT);
187                                        else
188                                                depth = (*objectIterator)->lastTestedRayResult.depth;
189                                        // intersection is acceptable only if it is within the node volume
190                                        if ((*objectIterator)->lastTestedRayResult.isIntersect &&
191                                                        depth < hitRec.depth &&
192                                                        depth > rayMin - epsilon &&
193                                                        depth < rayMax + epsilon)
194                                                hitRec = (*objectIterator)->lastTestedRayResult;
195                                }
196
197                                if(hitRec.isIntersect)
198                                        return true;
199                                else
200                                        if(stackPointer >= 0)
201                                        {
202                                                //have not found anything on this branch, pop the other one from stack
203                                                rayMin = traverseStack[stackPointer].rayMin;
204                                                rayMax = traverseStack[stackPointer].rayMax;
205                                                tNode = traverseStack[stackPointer].tNode;
206                                                stackPointer--;
207                                        }
208                                        else
209                                                return false;
210                        }
211                        else
212                        {
213                                //leaf is empty
214                                if(stackPointer >= 0)
215                                {
216                                        rayMin = traverseStack[stackPointer].rayMin;
217                                        rayMax = traverseStack[stackPointer].rayMax;
218                                        tNode = traverseStack[stackPointer].tNode;
219                                        stackPointer--;
220                                }
221                                else
222                                        return false;
223                        }
224                }
225                else //not a leaf, go deeper
226                {
227                        // gather cutting plane information from the nodeTable
228                        float nodeValue = nodeTable[tNode];
229                        unsigned char axis = *(unsigned int*)(&nodeValue) & 0x00000003;
230                        unsigned int bitfield = *(unsigned int*)(&nodeValue) & 0xfffffffc;
231                        float cutPlane = *(float*)(&bitfield);
232
233                        float origCoord = ray.origin[axis];
234                        float dirCoord = dirInverts[axis];
235
236            // calculate cutting plane - ray intersection distance 't'
237                        register float t = (cutPlane - origCoord) * dirCoord;
238
239                        // find the child volume nearer to the origin
240                        unsigned int nearNode = rightChild;
241                        unsigned int farNode = leftChild;
242                        if(origCoord + epsilon < cutPlane)
243                        {
244                                nearNode = leftChild;
245                                farNode = rightChild;
246                        }
247                        else
248                                if( origCoord < cutPlane + epsilon )            // origin is on the plane
249                                        if(     dirCoord < 0 )                                          // direction decides
250                                        {
251                                                tNode = leftChild; continue;
252                                        }
253                                        else
254                                        {
255                                                tNode = rightChild; continue;
256                                        }
257
258                        if(t <= 0 || rayMax + epsilon < t )
259                                // whole interval on near cell
260                                tNode = nearNode;
261                        else
262                        {
263                                if(t < rayMin - epsilon)
264                                        // whole interval on far cell
265                                        tNode = farNode;
266                                else
267                                {
268                                        // both cells
269                                        // push far branch to stack
270                                        stackPointer++;
271                                        traverseStack[stackPointer].rayMin = t;
272                                        traverseStack[stackPointer].rayMax = rayMax;
273                                        traverseStack[stackPointer].tNode = farNode;
274                                        // near branch is next to traverse
275                                        tNode = nearNode;
276                                        rayMax = t;
277                                }
278                        }
279                }
280        }
281        return false;
282}
283
284void KDTree::deleteLeaves(unsigned int nodeID)
285{
286        unsigned int rc;
287        unsigned int lc;
288        if(!followChildren(nodeID, lc, rc))
289        {
290                Intersectable** leafList = *(Intersectable***)&nodeTable[nodeID];
291                if(leafList)
292                {
293                        delete leafList;
294                        return;
295                }
296                return;
297        }
298        deleteLeaves(lc);
299        deleteLeaves(rc);
300}
301
302bool KDTree::followChildren(unsigned int& parent, unsigned int &leftChild, unsigned int &rightChild)
303{
304        unsigned int subPos = parent & nCacheLineNodes;                 // position within cache line
305        unsigned int supPos = parent & ~nCacheLineNodes;                // cache line index
306        unsigned int leafBit = ((*(unsigned int *)&nodeTable[parent | nCacheLineNodes]) >> subPos) & 0x1;
307        if(leafBit)                                                     // node is a leaf, has no children
308                return false;                                   // leftChild & rightChild unchanged
309        subPos = (subPos << 1) + 1;                     // index of left child
310        if(subPos >= nCacheLineNodes)           // out of cache line, node is pointer node
311        {
312                parent = ((unsigned int*)nodeTable)[parent];    // follow pointer
313                subPos = parent & nCacheLineNodes;                              // same procedure with new parent node
314                supPos = parent & ~nCacheLineNodes;                             // referenced node cannot be leaf
315                subPos = (subPos << 1) + 1;                                             // or pointer node
316        }
317        leftChild = (parent & ~nCacheLineNodes) | subPos;       // recombine parent and node address
318        rightChild = (parent & ~nCacheLineNodes) | subPos + 1;
319        return true;
320}
321
322// retrieve would-have-been child nodes of a leaf node, if they are not in a last row
323bool KDTree::getFreeNodes(unsigned int leaf, unsigned int& leftFree, unsigned int& rightFree)
324{
325        unsigned int subPos = leaf & nCacheLineNodes;                   // position within cache line
326        unsigned int supPos = leaf & ~nCacheLineNodes;                  // cache line index
327        subPos = (subPos << 1) + 1;                     // index of left child position
328        if(subPos < nCacheLineNodes)            // within cache line, may have free children positions
329        {
330                unsigned int grandSubPos = ((subPos + 1) << 1) + 2;     // rightmost grandchild
331                if(grandSubPos >= nCacheLineNodes)                                      // children are on last row
332                        return false;                                                                   // not suitable for free nodes
333                leftFree = (leaf & ~nCacheLineNodes) | subPos;
334                rightFree = (leaf & ~nCacheLineNodes) | subPos + 1;
335                return true;
336        }
337        else
338                return false;                                   // leaf is on last row, has no free children
339}
340
341bool KDTree::isLeaf(unsigned int xnode)
342{
343        return ((*(unsigned int *)&nodeTable[xnode | nCacheLineNodes]) >> (xnode & nCacheLineNodes)) & 0x1;
344}
345
346float KDTree::getBoundValue(unsigned int * extremeIndex, unsigned char axis)
347{
348        //if Most Significant Bit of 'extremeIndex' is 1, take the maximum
349        //if Most Significant Bit of 'extremeIndex' is 0, take the minimum
350        unsigned int index = *extremeIndex;
351        return (index & 0x80000000)?
352                (objects[ index & 0x7fffffff]->bbox.maxPoint[axis]):
353                (objects[ index & 0x7fffffff]->bbox.minPoint[axis]);
354}
355
356void KDTree::build(unsigned int nodeID, unsigned int* boundsArray, unsigned int nObjects, BoundingBox& limits, unsigned char axisNmask)
357{
358        //separate the axisNmask parameter to axis and mask
359        unsigned char axis = axisNmask & 07;
360        unsigned char mask = axisNmask & 070;
361
362        //calculate some values for the cost function
363        float costBase;         // area of face perpendicular to axis
364        float costSteep;        // half of circumference of face perpendicular to axis
365                                                // surface area / 2 = costBase + costSteep * (length of edge parallel to axis)
366                                                // cost = #(objects within volume) * surface area / 2
367        switch(axis)
368        {
369        case 0:
370                costBase = (limits.maxPoint.z - limits.minPoint.z)
371                                        *       (limits.maxPoint.y - limits.minPoint.y);
372                costSteep = (limits.maxPoint.z - limits.minPoint.z)
373                                        +       (limits.maxPoint.y - limits.minPoint.y);
374                break;
375        case 1:
376                costBase = (limits.maxPoint.x - limits.minPoint.x)
377                                        *       (limits.maxPoint.z - limits.minPoint.z);
378                costSteep = (limits.maxPoint.x - limits.minPoint.x)
379                                        +       (limits.maxPoint.z - limits.minPoint.z);
380                break;
381        case 2:
382                costBase = (limits.maxPoint.x - limits.minPoint.x)
383                                        *       (limits.maxPoint.y - limits.minPoint.y);
384                costSteep = (limits.maxPoint.x - limits.minPoint.x)
385                                        +       (limits.maxPoint.y - limits.minPoint.y);
386                break;
387        }
388
389        // references to extrema (corresponding to 'axis') of node volume
390        float& limiMin = axis?((axis==1)?limits.minPoint.y:limits.minPoint.z):limits.minPoint.x;
391        float& limiMax = axis?((axis==1)?limits.maxPoint.y:limits.maxPoint.z):limits.maxPoint.x;
392
393        //sort the the extrema of the objects
394        //we will need this to find the best splitting plane
395        quickSort(boundsArray, boundsArray + (nObjects << 1) - 1, axis);
396
397        // the following code section could be used to find the spatial and object medians
398        // either to implement a simpler and less effective subdivision rule
399        // or to restrict optimum search to the median interval (between two median)
400        //***************************************************************************
401        // find the value for the spatial median.... easy
402        //      float spatialMedian = (limiMax + limiMin) / 2.0f;
403        //find the position in the array corresponding to the spatial median value
404        //      unsigned int* spatialMedianPosition = findBound(boundsArray, nObjects << 1, spatialMedian, axis);
405        //find object median position (middle of array, surprise-surprise), and value
406        //      unsigned int* objectMedianPosition = boundsArray + nObjects;
407        //      float objectMedian = getBoundValue(objectMedianPosition, axis);
408        //****************************************************************************
409
410        // find the value for the left and right "cutting off empty space" (COES) planes
411        // remember the object extrema are already sorted
412        float leftCoesCut = objects[boundsArray[0] & 0x7fffffff]->bbox.minPoint[axis];
413        float rightCoesCut = objects[boundsArray[(nObjects << 1)-1] & 0x7fffffff]->bbox.maxPoint[axis];
414
415        if(rightCoesCut > limiMax)
416                rightCoesCut = limiMax;
417        if(leftCoesCut < limiMin)
418                leftCoesCut = limiMin;
419
420        // find the estimated cost of ray casting after COES
421        float leftCoesCutCost = ((limiMax - leftCoesCut) * costSteep + costBase) * (nObjects);
422        float rightCoesCutCost = ((rightCoesCut - limiMin) * costSteep + costBase) * (nObjects);
423        // take care of numerical accuracy problems that arise when
424        // limiMin == limiMax           (possible and allowed)
425        if(leftCoesCutCost < 0) leftCoesCutCost = FLT_MAX;
426        if(rightCoesCutCost < 0) rightCoesCutCost = FLT_MAX;
427
428        // creep through array, evaluate cost for every possible cut,
429        // maintaining the counters for intersected and completely passed objects
430        int leftCount = 0;                              // objects on the left of current cut
431        int intersectedCount = 0;               // object intersected by current cut
432
433        // we may limit the search to a partition of the array, eg. the median interval
434        unsigned int* searchIntervalStart;
435        unsigned int* searchIntervalEnd;
436
437        // full range optimum search
438        searchIntervalStart = boundsArray + 1;
439        searchIntervalEnd = boundsArray + (nObjects << 1) - 1;
440
441        // these three variable identify the best cut found so far:
442        float minimumCostFound;
443        float minimumCostCut;
444        unsigned int* minimumCostPosition;
445
446        // is left or right COES better? the better is our best candidate to start with
447        int minimumIntersected = 0;
448        int minimumLeft = 0;
449        if(leftCoesCutCost < rightCoesCutCost)
450        {
451                minimumCostFound = leftCoesCutCost;
452                minimumCostPosition = boundsArray;
453                minimumCostCut = leftCoesCut;
454                minimumIntersected = 0;
455                minimumLeft = 0;
456        }
457        else
458        {
459                minimumCostFound = rightCoesCutCost;
460                minimumCostPosition = boundsArray + (nObjects << 1);
461                minimumCostCut = rightCoesCut;
462                minimumIntersected = 0;
463                minimumLeft = nObjects;
464        }
465
466        // creep up to the search interval... count objects
467        unsigned int* trace = boundsArray;
468        for(; trace< searchIntervalStart; trace++ )
469        {
470                if(*trace & 0x80000000)
471                {
472                        intersectedCount--;
473                        leftCount++;
474                }
475                else
476                        intersectedCount++;
477        }
478        // evaluate the cost function, if better is found keep it
479        for(; trace< searchIntervalEnd; trace++ )
480        {
481                float cut = getBoundValue(trace, axis);
482                float costIfCutHere;
483                if(*trace & 0x80000000)
484                {
485                        intersectedCount--;
486                        leftCount++;
487                        costIfCutHere =
488                                (leftCount+intersectedCount) * ((cut - limiMin) * costSteep + costBase) +
489                                (nObjects - leftCount) * ((limiMax - cut) * costSteep + costBase);
490                        if(leftCount == 0 || leftCount + intersectedCount == nObjects )
491                                costIfCutHere = FLT_MAX;
492                        if(costIfCutHere < minimumCostFound ){
493                                minimumCostFound = costIfCutHere;
494                                minimumCostPosition = trace + 1;
495                                minimumIntersected = intersectedCount;
496                                minimumLeft = leftCount;
497                                minimumCostCut = cut;
498                        }
499                }
500                else
501                {
502                        costIfCutHere =
503                                (leftCount+intersectedCount) * ((cut - limiMin) * costSteep + costBase) +
504                                (nObjects - leftCount) * ((limiMax - cut) * costSteep + costBase);
505                        if(leftCount == 0 || leftCount + intersectedCount == nObjects)
506                                costIfCutHere = FLT_MAX;
507                        if(costIfCutHere < minimumCostFound)
508                        {
509                                minimumCostFound = costIfCutHere;
510                                minimumCostPosition = trace;
511                                minimumIntersected = intersectedCount;
512                                minimumLeft = leftCount;
513                                minimumCostCut = cut;
514                        }
515                        intersectedCount++;
516                }
517
518        }
519
520        // automatic termination, based on the cost estimation
521        if(     minimumCostFound * 1.0001 >= nObjects * ((limiMax - limiMin) * costSteep + costBase))
522        {
523                // no appropriate splitting found, mark this axis as failed
524                mask |= (010 << axis);
525                if(mask != 070 && (nObjects != 1 || minimumCostFound < 10.0))
526                {
527                        //if there is still a hopeful axis, do not split here, go on with another direction
528                        do{             axis = (axis + 1)%3;
529                        }while(mask & (010 << axis));
530                        build(nodeID, boundsArray, nObjects, limits, mask | axis);
531                        return;
532                }
533                // if all 3 have axes have already failed, stop here, make a leaf
534                // save leaf length and objects list
535                Intersectable** leafList = new Intersectable*[nObjects+1];
536                Intersectable** leafLoader = leafList;
537                *(unsigned int*)leafLoader = nObjects;
538                leafLoader++;
539                for(unsigned int* bubo = boundsArray;bubo < boundsArray + (nObjects << 1);bubo++)
540                        if(!(*bubo & 0x80000000))
541                        {
542                                *leafLoader = *(objects + *bubo);
543                                leafLoader++;
544                        }
545                // make the node contain a pointer to the list
546                ((Intersectable ***)nodeTable)[nodeID] = leafList;
547                delete boundsArray;
548                //mark the node as a leaf
549                ((unsigned int*)nodeTable)[nodeID | nCacheLineNodes] |= (0x1 << (nodeID & nCacheLineNodes));
550                //add orphaned nodes to free node list
551                unsigned int lFree;
552                unsigned int rFree;
553                if(getFreeNodes(nodeID, lFree, rFree))
554                {
555                        freeNodes->insert(lFree);
556                        freeNodes->insert(rFree);
557                }
558                return;
559        }
560        // grant another chance to previously failed directions:
561        mask = 00;
562        // insert a back-pointer if necessary
563        nodeID = makePointer(nodeID);
564        // store the splitting plane value in the node
565        unsigned int bitfield = *(unsigned int*)(&minimumCostCut) & 0xfffffffc | axis;
566        nodeTable[nodeID] =  *(float*)(&bitfield);
567        // mark node as a non-leaf node
568        *(unsigned int*)&nodeTable[nodeID | nCacheLineNodes] &= ~(0x1 << (nodeID & nCacheLineNodes));
569
570        currentDepth++; if(currentDepth > depth) depth = currentDepth;
571
572        // allocate child nodes, pass the objects on to them
573        unsigned int leftChildMade;
574        unsigned int rightChildMade;
575
576        // handle the two branches, the one with less objects first
577        bool firstBranch = false;
578        bool rightBranchHasMoreObjects = (minimumLeft << 1) + minimumIntersected < nObjects;
579        followChildren(nodeID, leftChildMade, rightChildMade);
580        do
581        {
582                firstBranch = !firstBranch;
583                if(firstBranch && rightBranchHasMoreObjects
584                        || !firstBranch && !rightBranchHasMoreObjects)
585                {
586                        // left branch
587                        if(minimumLeft + minimumIntersected)
588                        {
589                                // branch has objects
590                                // create new object extrema array
591                                BoundingBox diver = limits;
592                                limiMax = minimumCostCut;
593                                unsigned int * leftChildBounds = new unsigned int[(minimumLeft + minimumIntersected) << 1];
594                                unsigned int* copycatLeft = leftChildBounds;
595                                unsigned int* snoop = boundsArray;
596                                for(; snoop < minimumCostPosition; snoop++)
597                                {
598                                        // for every mimimum left of cutting plane
599                                        if(!(*snoop & 0x80000000))
600                                        {
601                                                // add one for the mimimum
602                                                *copycatLeft = *snoop;
603                                                copycatLeft++;
604                                                // add one for the maximum
605                                                *copycatLeft = *snoop | 0x80000000;
606                                                copycatLeft++;
607                                        }
608                                }
609                                unsigned char turnedAxis = axis;
610                                do{             turnedAxis = (turnedAxis + 1)%3;
611                                }while(mask & (010 << turnedAxis));
612                                build(leftChildMade, leftChildBounds, minimumLeft + minimumIntersected, limits, mask | turnedAxis);
613                                limits = diver;
614                        }
615                        else
616                        {
617                                // make empty leaf
618                                unsigned int empty = leftChildMade;
619                                (unsigned int *&)nodeTable[empty] = 0x0;
620                                *(unsigned int*)&nodeTable[empty | nCacheLineNodes] |= 0x1 << (empty & nCacheLineNodes);
621                                // add orphaned nodes to free list
622                                unsigned int eleLeft;
623                                unsigned int eleRight;
624                                if(getFreeNodes(empty, eleLeft, eleRight))
625                                {
626                                        freeNodes->insert(eleLeft);
627                                        freeNodes->insert(eleRight);
628                                }
629                        }
630                }
631                else
632                {
633                        if(nObjects - minimumLeft)
634                        {
635                                BoundingBox diver = limits;
636                                limiMin = minimumCostCut;
637
638                                // create new bounds array
639                                unsigned int * rightChildBounds;
640                                rightChildBounds = new unsigned int[(nObjects - minimumLeft) << 1 ];
641
642                                unsigned int* copycatRight = rightChildBounds;
643                                unsigned int* snoop = minimumCostPosition;
644                                for(; snoop < boundsArray + (nObjects << 1); snoop++)
645                                {
646                                        // for every maximum right of cutting plane
647                                        if(*snoop & 0x80000000)
648                                        {
649                                                *copycatRight = *snoop & 0x7fffffff;
650                                                copycatRight++;
651                                                *copycatRight = *snoop;
652                                                copycatRight++;
653                                        }
654                                }
655
656                                unsigned int turnedAxis = axis;
657                                do{     turnedAxis = (turnedAxis + 1)%3;
658                                }while(mask & (010 << turnedAxis));
659                                build(rightChildMade, rightChildBounds, nObjects - minimumLeft, limits, mask | turnedAxis);
660                                limits = diver;
661                        }
662                        else
663                        {
664                                // empty leaf
665                                unsigned int empty = rightChildMade;
666                                (unsigned int *&)nodeTable[empty] = 0x0;
667                                *(unsigned int*)&nodeTable[empty | nCacheLineNodes] |= 0x1 << (empty & nCacheLineNodes);
668                                // orphaned nodes
669                                unsigned int eleLeft;
670                                unsigned int eleRight;
671                                if(getFreeNodes(empty, eleLeft, eleRight))
672                                {
673                                        freeNodes->insert(eleLeft);
674                                        freeNodes->insert(eleRight);
675                                }
676                        }
677                }
678        }while(firstBranch);
679        currentDepth--;
680        delete boundsArray;
681}
682
683void KDTree::quickSort(unsigned int *lo0, unsigned int* hi0, unsigned char axis)
684{
685        if ( hi0 > lo0)
686        {
687                unsigned int* lo = lo0;
688                unsigned int* hi = hi0;
689
690         // establish partition element as the midpoint of the array.
691                unsigned int* midPos = 
692                        lo0 + ( ((unsigned int)hi0 - (unsigned int)lo0) / sizeof(unsigned int) / 2);
693                unsigned int midIndex = *midPos;
694
695         // loop through the array until indices cross
696                while( lo <= hi )
697                {
698            // find the first element that is greater than or equal to
699            // the partition element starting from the left Index.
700                        // getBoundValue(lo, axis) < mid ) && *lo != midIndex)
701                        while( ( lo < hi0 ) && ( compare(*lo, midIndex, axis)))
702                                ++lo;
703
704            // find an element that is smaller than or equal to
705            // the partition element starting from the right Index.
706                        //( getBoundValue(hi, axis) > mid ) && *hi != midIndex)
707            while( ( hi > lo0 ) && ( compare(midIndex, *hi, axis)))
708               --hi;           
709                        // if the indices have not crossed, swap
710            if( lo <= hi )
711                        {
712                                unsigned int temp = *lo;
713                                *lo = *hi;
714                                *hi = temp;
715                                ++lo;
716                                --hi;
717                        }         
718                }
719        // If the right index has not reached the left side of array
720        // must now sort the left partition.
721        if( lo0 < hi )
722            quickSort( lo0, hi, axis);
723                // If the left index has not reached the right side of array
724        // must now sort the right partition.
725                if( lo < hi0 )
726            quickSort( lo, hi0, axis);
727        }
728}
729
730KDTree::KDTree (Intersectable** iObjects, int nObjects){
731        // take parameters
732        this->nObjects          = nObjects;
733        objects = iObjects;
734
735        // determine bounding box
736        createBoundingBox (iObjects, nObjects);
737
738        // determine error because of 21 bits mantissa representaion
739        float largest = 0.0f;
740        for(int i=0; i<3; i++)
741                if( boundingBox.maxPoint[i] > largest)
742                        largest = boundingBox.maxPoint[i];
743        for(int o=0; o<3; o++)
744                if( - boundingBox.minPoint[o] > largest)
745                        largest = - boundingBox.minPoint[o];
746        epsilon = largest / 524288.0f; // e = X * 2^19
747
748        // set up memory structure constants
749        nCacheLineBytes = KDTREE_CACHE_LINE_BYTES;
750        // a cache line has 2^n - 1 nodes, plus 4 bytes (1 node) for leaf bits
751        nCacheLineNodes = nCacheLineBytes / sizeof(float) - 1;
752
753        // size of cache line memory pool
754        maxNCacheLines = nObjects * 16 / nCacheLineNodes + 1;
755        // cache lines already allocated
756        nCacheLines = 1;
757
758        unsigned int maxPossibleFree = maxNCacheLines * nCacheLineNodes / 2;
759        // allocate free nodes heap
760        freeNodes = new Heap<unsigned int>(maxPossibleFree);
761
762        // allocate cache line memory pool (large size, will be shrinked after the build)
763        // nodes in cache line pool + offset correction
764        poolNodeTable = (float*)malloc(sizeof(float) *
765                ( maxNCacheLines * (nCacheLineNodes + 1) + nCacheLineNodes + 1 ));
766        nodeTable = (float*)(((unsigned int)poolNodeTable / nCacheLineBytes + 1) * nCacheLineBytes);
767
768        // collect the initial object extremes to pass on to the build algorithm
769        unsigned int* objectBoundaries = new unsigned int[nObjects << 1];
770        for(unsigned int fill=0;fill<nObjects;fill++){
771                objectBoundaries[fill << 1] = fill & 0x7fffffff;
772                objectBoundaries[(fill << 1) + 1] = fill | 0x80000000;
773        }
774
775        // GO! GO! GO!
776        depth = currentDepth = 1;
777        build(0, objectBoundaries, nObjects, boundingBox, 0);
778        traverseStack = new TraverseStack[depth];
779
780        // get rid of garbage
781        // decrease cache line memory pool to fit actual memory need
782        realloc(poolNodeTable, (nCacheLines + 1) * nCacheLineBytes);
783
784        delete freeNodes;
785}
786
787unsigned int* KDTree::findBound(unsigned int *extremeArrayStart, unsigned int nBounds, float loc, unsigned char axis)
788{
789        while(nBounds > 1)
790        {
791                unsigned int wBounds = nBounds >> 1;
792                if(getBoundValue(extremeArrayStart+wBounds, axis) < loc)
793                {
794                        extremeArrayStart += wBounds;
795                        nBounds -= wBounds;
796                }
797                else
798                        nBounds = wBounds;
799        }
800        return extremeArrayStart;
801
802}
803
804/*!this is more complicated then one would predict. rules are:
805- if values are significantly different, no problem
806- the minimum of a patch is smaller than its maximum
807- if a maximum and a minimum are near, the other extrema of the patches decide
808- if all four extrema are near, the 2 ends of a patch have to be simultaneusly
809smaller or bigger than the other two ends
810*/
811bool KDTree::compare(unsigned int aIndex, unsigned int bIndex, unsigned char axis)
812{
813        unsigned int caIndex = aIndex & 0x7fffffff;
814        float aBegin = objects[caIndex]->bbox.minPoint[axis];
815        float aEnd = objects[caIndex]->bbox.maxPoint[axis];
816        unsigned int cbIndex = bIndex & 0x7fffffff;
817        float bBegin = objects[cbIndex]->bbox.minPoint[axis];
818        float bEnd = objects[cbIndex]->bbox.maxPoint[axis];
819
820        if(aIndex & 0x80000000)
821                if(bIndex & 0x80000000)
822                {
823                        if(caIndex == cbIndex)
824                                return false;
825                        if(aEnd + epsilon < bEnd)
826                                return true;
827                        if(bEnd + epsilon < aEnd)
828                                return false;
829                        if(aBegin + epsilon < bBegin)
830                                return true;
831                        if(bBegin + epsilon < aBegin)
832                                return false;
833                        return caIndex < cbIndex;
834                }
835                else
836                {
837                        if(caIndex == cbIndex)
838                                return false;
839                        if(aEnd + epsilon < bBegin)
840                                return true;
841                        if(bBegin + epsilon < aEnd)
842                                return false;
843                        if(aBegin + epsilon < bEnd)
844                                return true;
845                        return caIndex < cbIndex;
846
847                }
848        else
849                if(bIndex & 0x80000000)
850                {
851                        if(caIndex == cbIndex)
852                                return true;
853                        if(bEnd + epsilon < aBegin)
854                                return false;
855                        if(aBegin + epsilon < bEnd)
856                                return true;
857                        if(bBegin + epsilon < aEnd)
858                                return false;
859                        return caIndex + epsilon < cbIndex;
860                }
861                else
862                {
863                        if(caIndex == cbIndex)
864                                return false;
865                        if(aBegin + epsilon < bBegin)
866                                return true;
867                        if(bBegin + epsilon < aBegin)
868                                return false;
869                        if(aEnd + epsilon < bEnd)
870                                return true;
871                        if(bEnd + epsilon < aEnd)
872                                return false;
873                        return caIndex < cbIndex;
874
875                }
876}
877
878inline unsigned int KDTree::makePointer(unsigned int originalNode)
879{
880        unsigned int subPos = originalNode & nCacheLineNodes;           // position within cache line
881        unsigned int supPos = originalNode & ~nCacheLineNodes;          // index of cache line
882        subPos = (subPos << 1) + 1;                                                                     // child position
883        if(subPos < nCacheLineNodes)                                                            // within cache line
884        {
885                return originalNode;
886        }
887        else                                                                                                            // find free node
888        {
889                unsigned int emptyNode;
890                if(!freeNodes->isEmpty())
891                {
892                        emptyNode = freeNodes->removeMin();
893                }
894                else
895                {
896                        //allocate a new cache line
897                        emptyNode = nCacheLines * (nCacheLineNodes + 1);        // root of new cache line
898                        nCacheLines++;
899                        if(nCacheLines > maxNCacheLines)
900                        {
901                                maxNCacheLines *= 1.1;
902                                float * largerPoolNodeTable = (float*)malloc(
903                                        sizeof(float) * ( maxNCacheLines * (nCacheLineNodes + 1) + nCacheLineNodes + 1 ));
904                                float* largerNodeTable = (float*)(((unsigned int)largerPoolNodeTable / nCacheLineBytes + 1) * nCacheLineBytes);
905                                memcpy(largerNodeTable, nodeTable, nCacheLines * nCacheLineBytes);
906                                nodeTable = largerNodeTable;
907                                free(poolNodeTable);
908                                poolNodeTable = largerPoolNodeTable;
909                        }
910                }
911                // mark as non-leaf (being last-row and non-leaf means pointer node)
912                *(unsigned int*)&nodeTable[originalNode | nCacheLineNodes] &= ~(0x1 << (originalNode & nCacheLineNodes));
913                // insert index of empty node
914                ((unsigned int*)nodeTable)[originalNode] = emptyNode;
915                return emptyNode;
916        }
917}
Note: See TracBrowser for help on using the repository browser.