qbvhaccel.h

Go to the documentation of this file.

/***************************************************************************
 *   Copyright (C) 1998-2010 by authors (see AUTHORS.txt )                 *
 *                                                                         *
 *   This file is part of LuxRays.                                         *
 *                                                                         *
 *   LuxRays is free software; you can redistribute it and/or modify       *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 3 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   LuxRays is distributed in the hope that it will be useful,            *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>. *
 *                                                                         *
 *   LuxRays website: http://www.luxrender.net                             *
 ***************************************************************************/

#ifndef _LUXRAYS_QBVHACCEL_H
#define _LUXRAYS_QBVHACCEL_H

#include <string.h>
#include <xmmintrin.h>
#include <boost/cstdint.hpp>

#include "luxrays/luxrays.h"
#include "luxrays/core/acceleretor.h"

using boost::int32_t;

namespace luxrays {

#if defined(WIN32) && !defined(__CYGWIN__)
class __declspec(align(16)) QuadRay {
#else
class QuadRay {
#endif
public:
        QuadRay(const Ray &ray)
        {
                ox = _mm_set1_ps(ray.o.x);
                oy = _mm_set1_ps(ray.o.y);
                oz = _mm_set1_ps(ray.o.z);

                dx = _mm_set1_ps(ray.d.x);
                dy = _mm_set1_ps(ray.d.y);
                dz = _mm_set1_ps(ray.d.z);

                mint = _mm_set1_ps(ray.mint);
                maxt = _mm_set1_ps(ray.maxt);
        }

        __m128 ox, oy, oz;
        __m128 dx, dy, dz;
        mutable __m128 mint, maxt;
#if defined(WIN32) && !defined(__CYGWIN__)
};
#else
} __attribute__ ((aligned(16)));
#endif

class QuadTriangle : public Aligned16 {
public:

        QuadTriangle(const Triangle *tris, const Point *verts,
                        const unsigned int p1,
                        const unsigned int p2,
                        const unsigned int p3,
                        const unsigned int p4) {

                primitives[0] = p1;
                primitives[1] = p2;
                primitives[2] = p3;
                primitives[3] = p4;

                for (u_int i = 0; i < 4; ++i) {
                        const Triangle *t = &tris[primitives[i]];

                        reinterpret_cast<float *> (&origx)[i] = verts[t->v[0]].x;
                        reinterpret_cast<float *> (&origy)[i] = verts[t->v[0]].y;
                        reinterpret_cast<float *> (&origz)[i] = verts[t->v[0]].z;

                        reinterpret_cast<float *> (&edge1x)[i] = verts[t->v[1]].x - verts[t->v[0]].x;
                        reinterpret_cast<float *> (&edge1y)[i] = verts[t->v[1]].y - verts[t->v[0]].y;
                        reinterpret_cast<float *> (&edge1z)[i] = verts[t->v[1]].z - verts[t->v[0]].z;

                        reinterpret_cast<float *> (&edge2x)[i] = verts[t->v[2]].x - verts[t->v[0]].x;
                        reinterpret_cast<float *> (&edge2y)[i] = verts[t->v[2]].y - verts[t->v[0]].y;
                        reinterpret_cast<float *> (&edge2z)[i] = verts[t->v[2]].z - verts[t->v[0]].z;
                }
        }

        ~QuadTriangle() {
        }

        BBox WorldBound(const Triangle *tris, const Point *verts) const {
                return Union(
                                Union(tris[primitives[0]].WorldBound(verts), tris[primitives[1]].WorldBound(verts)),
                                Union(tris[primitives[2]].WorldBound(verts), tris[primitives[3]].WorldBound(verts)));
        }

        bool Intersect(const QuadRay &ray4, const Ray &ray, RayHit *rayHit) const {
                const __m128 zero = _mm_setzero_ps();
                const __m128 s1x = _mm_sub_ps(_mm_mul_ps(ray4.dy, edge2z),
                                _mm_mul_ps(ray4.dz, edge2y));
                const __m128 s1y = _mm_sub_ps(_mm_mul_ps(ray4.dz, edge2x),
                                _mm_mul_ps(ray4.dx, edge2z));
                const __m128 s1z = _mm_sub_ps(_mm_mul_ps(ray4.dx, edge2y),
                                _mm_mul_ps(ray4.dy, edge2x));
                const __m128 divisor = _mm_add_ps(_mm_mul_ps(s1x, edge1x),
                                _mm_add_ps(_mm_mul_ps(s1y, edge1y),
                                _mm_mul_ps(s1z, edge1z)));
                __m128 test = _mm_cmpneq_ps(divisor, zero);
                const __m128 inverse = _mm_div_ps(_mm_set_ps1(1.f), divisor);
                const __m128 dx = _mm_sub_ps(ray4.ox, origx);
                const __m128 dy = _mm_sub_ps(ray4.oy, origy);
                const __m128 dz = _mm_sub_ps(ray4.oz, origz);
                const __m128 b1 = _mm_mul_ps(_mm_add_ps(_mm_mul_ps(dx, s1x),
                                _mm_add_ps(_mm_mul_ps(dy, s1y), _mm_mul_ps(dz, s1z))),
                                inverse);
                test = _mm_and_ps(test, _mm_cmpge_ps(b1, zero));
                const __m128 s2x = _mm_sub_ps(_mm_mul_ps(dy, edge1z),
                                _mm_mul_ps(dz, edge1y));
                const __m128 s2y = _mm_sub_ps(_mm_mul_ps(dz, edge1x),
                                _mm_mul_ps(dx, edge1z));
                const __m128 s2z = _mm_sub_ps(_mm_mul_ps(dx, edge1y),
                                _mm_mul_ps(dy, edge1x));
                const __m128 b2 = _mm_mul_ps(_mm_add_ps(_mm_mul_ps(ray4.dx, s2x),
                                _mm_add_ps(_mm_mul_ps(ray4.dy, s2y), _mm_mul_ps(ray4.dz, s2z))),
                                inverse);
                const __m128 b0 = _mm_sub_ps(_mm_set1_ps(1.f),
                                _mm_add_ps(b1, b2));
                test = _mm_and_ps(test, _mm_and_ps(_mm_cmpge_ps(b2, zero),
                                _mm_cmpge_ps(b0, zero)));
                const __m128 t = _mm_mul_ps(_mm_add_ps(_mm_mul_ps(edge2x, s2x),
                                _mm_add_ps(_mm_mul_ps(edge2y, s2y),
                                _mm_mul_ps(edge2z, s2z))), inverse);
                test = _mm_and_ps(test,
                                _mm_and_ps(_mm_cmpgt_ps(t, ray4.mint),
                                _mm_cmplt_ps(t, ray4.maxt)));

                const int testmask = _mm_movemask_ps(test);             
                if (testmask == 0) return false;

                u_int hit;
                if ((testmask & (testmask - 1)) == 0) {
                        hit = UIntLog2(testmask);
                        ray.maxt = reinterpret_cast<const float *> (&t)[hit];
                } else {
                        for (u_int i = 0; i < 4; ++i) {
                                if (reinterpret_cast<const int *> (&test)[i] && reinterpret_cast<const float *> (&t)[i] < ray.maxt) {
                                        hit = i;
                                        ray.maxt = reinterpret_cast<const float *> (&t)[i];
                                }
                        }
                }

                ray4.maxt = _mm_set1_ps(ray.maxt);

                rayHit->t = ray.maxt;
                rayHit->b1 = reinterpret_cast<const float *> (&b1)[hit];
                rayHit->b2 = reinterpret_cast<const float *> (&b2)[hit];
                rayHit->index = primitives[hit];

                return true;
        }

        __m128 origx, origy, origz;
        __m128 edge1x, edge1y, edge1z;
        __m128 edge2x, edge2y, edge2z;
        unsigned int primitives[4];
};

// This code is based on Flexray by Anthony Pajot (anthony.pajot@alumni.enseeiht.fr)

#define NB_BINS 8

class QBVHNode {
public: 
        // The constant used to represent empty leaves. there would have been
        // a conflict with a normal leaf if there were 16 quads,
        // starting at 2^27 in the quads array... very improbable.
        // using MININT (0x80000000) can produce conflict when initializing a
        // QBVH with less than 4 vertices at the beginning :
        // the number of quads - 1 would give 0, and it would start at 0
        // in the quads array
        static const int32_t emptyLeafNode = 0xffffffff;
        
        __m128 bboxes[2][3];

        int32_t children[4];
        
        inline QBVHNode() {
                for (int i = 0; i < 3; ++i) {
                        bboxes[0][i] = _mm_set1_ps(INFINITY);
                        bboxes[1][i] = _mm_set1_ps(-INFINITY);
                }
                
                // All children are empty leaves by default
                for (int i = 0; i < 4; ++i)
                        children[i] = emptyLeafNode;
        }

        inline bool ChildIsLeaf(int i) const {
                return (children[i] < 0);
        }

        inline static bool IsLeaf(int32_t index) {
                return (index < 0);
        }

        inline bool LeafIsEmpty(int i) const {
                return (children[i] == emptyLeafNode);
        }

        inline static bool IsEmpty(int32_t index) {
                return (index == emptyLeafNode);
        }
        
        inline u_int NbQuadsInLeaf(int i) const {
                return static_cast<u_int>((children[i] >> 27) & 0xf) + 1;
        }

        inline static u_int NbQuadPrimitives(int32_t index) {
                return static_cast<u_int>((index >> 27) & 0xf) + 1;
        }
        
        inline u_int NbPrimitivesInLeaf(int i) const {
                return NbQuadsInLeaf(i) * 4;
        }

        inline u_int FirstQuadIndexForLeaf(int i) const {
                return children[i] & 0x07ffffff;
        }
        
        inline static u_int FirstQuadIndex(int32_t index) {
                return index & 0x07ffffff;
        }

        inline void InitializeLeaf(int i, u_int nbQuads, u_int firstQuadIndex) {
                // Take care to make a valid initialisation of the leaf.
                if (nbQuads == 0) {
                        children[i] = emptyLeafNode;
                } else {
                        // Put the negative sign in a plateform independent way
                        children[i] = 0x80000000;//-1L & ~(-1L >> 1L);
                        
                        children[i] |=  ((static_cast<int32_t>(nbQuads) - 1) & 0xf) << 27;

                        children[i] |= static_cast<int32_t>(firstQuadIndex) & 0x07ffffff;
                }
        }

        inline void SetBBox(int i, const BBox &bbox) {
                for (int axis = 0; axis < 3; ++axis) {
                        reinterpret_cast<float *>(&(bboxes[0][axis]))[i] = bbox.pMin[axis];
                        reinterpret_cast<float *>(&(bboxes[1][axis]))[i] = bbox.pMax[axis];
                }
        }

        
        inline int32_t BBoxIntersect(const QuadRay &ray4, const __m128 invDir[3],
                const int sign[3]) const {
                __m128 tMin = ray4.mint;
                __m128 tMax = ray4.maxt;

                // X coordinate
                tMin = _mm_max_ps(tMin, _mm_mul_ps(_mm_sub_ps(bboxes[sign[0]][0],
                                ray4.ox), invDir[0]));
                tMax = _mm_min_ps(tMax, _mm_mul_ps(_mm_sub_ps(bboxes[1 - sign[0]][0],
                                ray4.ox), invDir[0]));

                // Y coordinate
                tMin = _mm_max_ps(tMin, _mm_mul_ps(_mm_sub_ps(bboxes[sign[1]][1],
                                ray4.oy), invDir[1]));
                tMax = _mm_min_ps(tMax, _mm_mul_ps(_mm_sub_ps(bboxes[1 - sign[1]][1],
                                ray4.oy), invDir[1]));

                // Z coordinate
                tMin = _mm_max_ps(tMin, _mm_mul_ps(_mm_sub_ps(bboxes[sign[2]][2],
                                ray4.oz), invDir[2]));
                tMax = _mm_min_ps(tMax, _mm_mul_ps(_mm_sub_ps(bboxes[1 - sign[2]][2],
                                ray4.oz), invDir[2]));

                //return the visit flags
                return _mm_movemask_ps(_mm_cmpge_ps(tMax, tMin));
        }
};

/***************************************************/
class QBVHAccel  : public Accelerator {
public:
        QBVHAccel(const Context *context, u_int mp, u_int fst, u_int sf);

        ~QBVHAccel();

        BBox WorldBound() const { return worldBound; }

        AcceleratorType GetType() const { return ACCEL_QBVH; }
        void Init(const std::deque<Mesh *> &meshes, const unsigned int totalVertexCount,
                const unsigned int totalTriangleCount);

        const TriangleMeshID GetMeshID(const unsigned int index) const { return meshIDs[index]; }
        const TriangleMeshID *GetMeshIDTable() const { return meshIDs; }
        const TriangleID GetMeshTriangleID(const unsigned int index) const { return meshTriangleIDs[index]; }
        const TriangleID *GetMeshTriangleIDTable() const { return meshTriangleIDs; }

        bool Intersect(const Ray *ray, RayHit *hit) const;

        const TriangleMesh *GetPreprocessedMesh() const { return preprocessedMesh; }

        friend class MQBVHAccel;
        friend class OpenCLIntersectionDevice;

private:
        // A special initialization method used only by MQBVHAccel
        void Init(const Mesh *m);

        void BuildTree(u_int start, u_int end, u_int *primsIndexes,
                BBox *primsBboxes, Point *primsCentroids, const BBox &nodeBbox,
                const BBox &centroidsBbox, int32_t parentIndex, int32_t childIndex,
                int depth);
        
        void CreateTempLeaf(int32_t parentIndex, int32_t childIndex, u_int start, u_int end,
                const BBox &nodeBbox);

        inline int32_t CreateIntermediateNode(int32_t parentIndex, int32_t childIndex,
                const BBox &nodeBbox) {
                int32_t index = nNodes++; // increment after assignment
                if (nNodes >= maxNodes) {
                        QBVHNode *newNodes = AllocAligned<QBVHNode>(2 * maxNodes);
                        memcpy(newNodes, nodes, sizeof(QBVHNode) * maxNodes);
                        for (u_int i = 0; i < maxNodes; ++i)
                                newNodes[maxNodes + i] = QBVHNode();
                        FreeAligned(nodes);
                        nodes = newNodes;
                        maxNodes *= 2;
                }

                if (parentIndex >= 0) {
                        nodes[parentIndex].children[childIndex] = index;
                        nodes[parentIndex].SetBBox(childIndex, nodeBbox);
                }
                return index;
        }

        void PreSwizzle(int32_t nodeIndex, u_int *primsIndexes);

        void CreateSwizzledLeaf(int32_t parentIndex, int32_t childIndex, 
                u_int *primsIndexes);

        u_int nQuads;

        QuadTriangle *prims;

        QBVHNode *nodes;

        u_int nNodes, maxNodes;

        BBox worldBound;

        u_int fullSweepThreshold;

        u_int skipFactor;

        u_int maxPrimsPerLeaf;

        const Context *ctx;
        TriangleMesh *preprocessedMesh;
        const Mesh *mesh;
        TriangleMeshID *meshIDs;
        TriangleID *meshTriangleIDs;

        bool initialized;
};

}

#endif  /* _LUXRAYS_QBVHACCEL_H */

---