kdtree.java

Sunflow是一个照片级的渲染系统

            nextOffset &= ~(3 << 30);
            switch (axis) {
                case 0:
                    max = bounds.getMaximum().x;
                    bounds.getMaximum().x = split;
                    v0 = dumpObj(nextOffset, vertOffset, maxN, bounds, file, mtlFile);
                    // restore and go to other side
                    bounds.getMaximum().x = max;
                    min = bounds.getMinimum().x;
                    bounds.getMinimum().x = split;
                    v0 = dumpObj(nextOffset + 2, v0, maxN, bounds, file, mtlFile);
                    bounds.getMinimum().x = min;
                    break;
                case 1 << 30:
                    max = bounds.getMaximum().y;
                    bounds.getMaximum().y = split;
                    v0 = dumpObj(nextOffset, vertOffset, maxN, bounds, file, mtlFile);
                    // restore and go to other side
                    bounds.getMaximum().y = max;
                    min = bounds.getMinimum().y;
                    bounds.getMinimum().y = split;
                    v0 = dumpObj(nextOffset + 2, v0, maxN, bounds, file, mtlFile);
                    bounds.getMinimum().y = min;
                    break;
                case 2 << 30:
                    max = bounds.getMaximum().z;
                    bounds.getMaximum().z = split;
                    v0 = dumpObj(nextOffset, vertOffset, maxN, bounds, file, mtlFile);
                    // restore and go to other side
                    bounds.getMaximum().z = max;
                    min = bounds.getMinimum().z;
                    bounds.getMinimum().z = split;
                    v0 = dumpObj(nextOffset + 2, v0, maxN, bounds, file, mtlFile);
                    // restore and go to other side
                    bounds.getMinimum().z = min;
                    break;
                default:
                    v0 = vertOffset;
                    break;
            }
            return v0;
        }
    }

    // type is encoded as 2 shifted bits
    private static final long CLOSED = 0L << 30;
    private static final long PLANAR = 1L << 30;
    private static final long OPENED = 2L << 30;
    private static final long TYPE_MASK = 3L << 30;

    // pack split values into a 64bit integer
    private static long pack(float split, long type, int axis, int object) {
        // pack float in sortable form
        int f = Float.floatToRawIntBits(split);
        int top = f ^ ((f >> 31) | 0x80000000);
        long p = ((long) top & 0xFFFFFFFFL) << 32;
        p |= type; // encode type as 2 bits
        p |= ((long) axis) << 28; // encode axis as 2 bits
        p |= (object & 0xFFFFFFFL); // pack object number
        return p;
    }

    private static int unpackObject(long p) {
        return (int) (p & 0xFFFFFFFL);
    }

    private static int unpackAxis(long p) {
        return (int) (p >>> 28) & 3;
    }

    private static long unpackSplitType(long p) {
        return p & TYPE_MASK;
    }

    private static float unpackSplit(long p) {
        int f = (int) ((p >>> 32) & 0xFFFFFFFFL);
        int m = ((f >>> 31) - 1) | 0x80000000;
        return Float.intBitsToFloat(f ^ m);
    }

    // radix sort on top 36 bits - returns sorted result
    private static void radix12(long[] splits, int n) {
        // allocate working memory
        final int[] hist = new int[2048];
        final long[] sorted = new long[n];
        // parallel histogramming pass
        for (int i = 0; i < n; i++) {
            long pi = splits[i];
            hist[0x000 + ((int) (pi >>> 28) & 0x1FF)]++;
            hist[0x200 + ((int) (pi >>> 37) & 0x1FF)]++;
            hist[0x400 + ((int) (pi >>> 46) & 0x1FF)]++;
            hist[0x600 + ((int) (pi >>> 55))]++;
        }

        // sum the histograms - each histogram entry records the number of
        // values preceding itself.
        {
            int sum0 = 0, sum1 = 0, sum2 = 0, sum3 = 0;
            int tsum;
            for (int i = 0; i < 512; i++) {
                tsum = hist[0x000 + i] + sum0;
                hist[0x000 + i] = sum0 - 1;
                sum0 = tsum;
                tsum = hist[0x200 + i] + sum1;
                hist[0x200 + i] = sum1 - 1;
                sum1 = tsum;
                tsum = hist[0x400 + i] + sum2;
                hist[0x400 + i] = sum2 - 1;
                sum2 = tsum;
                tsum = hist[0x600 + i] + sum3;
                hist[0x600 + i] = sum3 - 1;
                sum3 = tsum;
            }
        }

        // read/write histogram passes
        for (int i = 0; i < n; i++) {
            long pi = splits[i];
            int pos = (int) (pi >>> 28) & 0x1FF;
            sorted[++hist[0x000 + pos]] = pi;
        }
        for (int i = 0; i < n; i++) {
            long pi = sorted[i];
            int pos = (int) (pi >>> 37) & 0x1FF;
            splits[++hist[0x200 + pos]] = pi;
        }
        for (int i = 0; i < n; i++) {
            long pi = splits[i];
            int pos = (int) (pi >>> 46) & 0x1FF;
            sorted[++hist[0x400 + pos]] = pi;
        }
        for (int i = 0; i < n; i++) {
            long pi = sorted[i];
            int pos = (int) (pi >>> 55);
            splits[++hist[0x600 + pos]] = pi;
        }
    }

    private static class BuildTask {
        long[] splits;
        int numObjects;
        int n;
        byte[] leftRightTable;

        BuildTask(int numObjects) {
            splits = new long[6 * numObjects];
            this.numObjects = numObjects;
            n = 0;
            // 2 bits per object
            leftRightTable = new byte[(numObjects + 3) / 4];
        }

        BuildTask(int numObjects, BuildTask parent) {
            splits = new long[6 * numObjects];
            this.numObjects = numObjects;
            n = 0;
            leftRightTable = parent.leftRightTable;
        }
    }

    private void buildTree(float minx, float maxx, float miny, float maxy, float minz, float maxz, BuildTask task, int depth, IntArray tempTree, int offset, IntArray tempList, BuildStats stats) {
        // get node bounding box extents
        if (task.numObjects > maxPrims && depth < MAX_DEPTH) {
            float dx = maxx - minx;
            float dy = maxy - miny;
            float dz = maxz - minz;
            // search for best possible split
            float bestCost = INTERSECT_COST * task.numObjects;
            int bestAxis = -1;
            int bestOffsetStart = -1;
            int bestOffsetEnd = -1;
            float bestSplit = 0;
            boolean bestPlanarLeft = false;
            int bnl = 0, bnr = 0;
            // inverse area of the bounding box (factor of 2 ommitted)
            float area = (dx * dy + dy * dz + dz * dx);
            float ISECT_COST = INTERSECT_COST / area;
            // setup counts for each axis
            int[] nl = { 0, 0, 0 };
            int[] nr = { task.numObjects, task.numObjects, task.numObjects };
            // setup bounds for each axis
            float[] dp = { dy * dz, dz * dx, dx * dy };
            float[] ds = { dy + dz, dz + dx, dx + dy };
            float[] nodeMin = { minx, miny, minz };
            float[] nodeMax = { maxx, maxy, maxz };
            // search for best cost
            int nSplits = task.n;
            long[] splits = task.splits;
            byte[] lrtable = task.leftRightTable;
            for (int i = 0; i < nSplits;) {
                // extract current split
                long ptr = splits[i];
                float split = unpackSplit(ptr);
                int axis = unpackAxis(ptr);
                // mark current position
                int currentOffset = i;
                // count number of primitives start/stopping/lying on the
                // current plane
                int pClosed = 0, pPlanar = 0, pOpened = 0;
                long ptrMasked = ptr & (~TYPE_MASK & 0xFFFFFFFFF0000000L);
                long ptrClosed = ptrMasked | CLOSED;
                long ptrPlanar = ptrMasked | PLANAR;
                long ptrOpened = ptrMasked | OPENED;
                while (i < nSplits && (splits[i] & 0xFFFFFFFFF0000000L) == ptrClosed) {
                    int obj = unpackObject(splits[i]);
                    lrtable[obj >>> 2] = 0;
                    pClosed++;
                    i++;
                }
                while (i < nSplits && (splits[i] & 0xFFFFFFFFF0000000L) == ptrPlanar) {
                    int obj = unpackObject(splits[i]);
                    lrtable[obj >>> 2] = 0;
                    pPlanar++;
                    i++;
                }
                while (i < nSplits && (splits[i] & 0xFFFFFFFFF0000000L) == ptrOpened) {
                    int obj = unpackObject(splits[i]);
                    lrtable[obj >>> 2] = 0;
                    pOpened++;
                    i++;
                }
                // now we have summed all contributions from this plane
                nr[axis] -= pPlanar + pClosed;
                // compute cost
                if (split >= nodeMin[axis] && split <= nodeMax[axis]) {
                    // left and right surface area (factor of 2 ommitted)
                    float dl = split - nodeMin[axis];
                    float dr = nodeMax[axis] - split;
                    float lp = dp[axis] + dl * ds[axis];
                    float rp = dp[axis] + dr * ds[axis];
                    // planar prims go to smallest cell always
                    boolean planarLeft = dl < dr;
                    int numLeft = nl[axis] + (planarLeft ? pPlanar : 0);
                    int numRight = nr[axis] + (planarLeft ? 0 : pPlanar);
                    float eb = ((numLeft == 0 && dl > 0) || (numRight == 0 && dr > 0)) ? EMPTY_BONUS : 0;
                    float cost = TRAVERSAL_COST + ISECT_COST * (1 - eb) * (lp * numLeft + rp * numRight);
                    if (cost < bestCost) {
                        bestCost = cost;
                        bestAxis = axis;
                        bestSplit = split;
                        bestOffsetStart = currentOffset;
                        bestOffsetEnd = i;
                        bnl = numLeft;
                        bnr = numRight;
                        bestPlanarLeft = planarLeft;
                    }
                }
                // move objects left
                nl[axis] += pOpened + pPlanar;
            }
            // debug check for correctness of the scan
            for (int axis = 0; axis < 3; axis++) {
                int numLeft = nl[axis];
                int numRight = nr[axis];
                if (numLeft != task.numObjects || numRight != 0)
                    UI.printError(Module.ACCEL, "Didn't scan full range of objects @depth=%d. Left overs for axis %d: [L: %d] [R: %d]", depth, axis, numLeft, numRight);
            }
            // found best split?
            if (bestAxis != -1) {
                // allocate space for child nodes
                BuildTask taskL = new BuildTask(bnl, task);
                BuildTask taskR = new BuildTask(bnr, task);
                int lk = 0, rk = 0;
                for (int i = 0; i < bestOffsetStart; i++) {
                    long ptr = splits[i];
                    if (unpackAxis(ptr) == bestAxis) {
                        if (unpackSplitType(ptr) != CLOSED) {