experiments

kernels/bvh/bvh.h

@@ -1132,10 +1132,16 @@ namespace embree
     };
 
     /*! BVHN Quantized Node */
-    struct __aligned(64) QuantizedNode : public BaseNode
+    struct __aligned(32) QuantizedNode : public BaseNode
     {
       using BaseNode::children;
-#if 0
+
+/* slower than explicit 8 or 16 bit quantization */
+//#define COMPACT10BIT
+
+#ifndef COMPACT10BIT
+
+#if 1
       typedef unsigned char T;
       static const T MIN_QUAN = 0;
       static const T MAX_QUAN = 255;
@@ -1145,6 +1151,12 @@ namespace embree
       static const T MAX_QUAN = 65535;
 #endif
 
+#else
+      typedef int T;
+      static const T MIN_QUAN = 0;
+      static const T MAX_QUAN = 1023;
+#endif
+
       struct Create2
       {
         template<typename BuildRecord>
@@ -1207,7 +1219,8 @@ namespace embree
                                          T lower_quant[N],
                                          T upper_quant[N],
                                          float &start,
-                                         float &scale)
+                                         float &scale,
+                                         const size_t dim)
       {
         const vbool<N> m_valid = lower != vfloat<N>(pos_inf);
         const float minF = reduce_min(lower);
@@ -1239,16 +1252,20 @@ namespace embree
         i_ceil_upper   = select(m_upper_correction,i_ceil_upper +1,i_ceil_upper);
 
         /* disable invalid lanes */
-        i_floor_lower = select(m_valid,i_floor_lower,MAX_QUAN);
+        i_floor_lower = select(m_valid,i_floor_lower,(int)MAX_QUAN);
         i_ceil_upper  = select(m_valid,i_ceil_upper ,0);
 
         /* store as uchar to memory */
+#ifdef COMPACT10BIT
+        i_floor_lower = vint<N>::loadu(lower_quant) | (i_floor_lower << (dim*10));
+        i_ceil_upper  = vint<N>::loadu(upper_quant) | (i_ceil_upper  << (dim*10));
+#endif
         vint<N>::store(lower_quant,i_floor_lower);
         vint<N>::store(upper_quant,i_ceil_upper);
         start = minF;
         scale = scale_diff;
 
-#if defined(DEBUG)
+#if defined(DEBUG) && !defined(COMPACT10BIT)
         vfloat<N> extract_lower( vint<N>::load(lower_quant) );
         vfloat<N> extract_upper( vint<N>::load(upper_quant) );
         vfloat<N> final_extract_lower = madd(extract_lower,scale_diff,minF);
@@ -1261,27 +1278,52 @@ namespace embree
       __forceinline void init(AlignedNode& node)
       {
         for (size_t i=0;i<N;i++) children[i] = emptyNode;
-        init_dim(node.lower_x,node.upper_x,lower_x,upper_x,start.x,scale.x);
-        init_dim(node.lower_y,node.upper_y,lower_y,upper_y,start.y,scale.y);
-        init_dim(node.lower_z,node.upper_z,lower_z,upper_z,start.z,scale.z);
+#ifdef COMPACT10BIT
+        lower_xyz = vint<N>(zero);
+        upper_xyz = vint<N>(zero);
+        init_dim(node.lower_x,node.upper_x,(T*)&lower_xyz,(T*)&upper_xyz,start.x,scale.x,0);
+        init_dim(node.lower_y,node.upper_y,(T*)&lower_xyz,(T*)&upper_xyz,start.y,scale.y,1);
+        init_dim(node.lower_z,node.upper_z,(T*)&lower_xyz,(T*)&upper_xyz,start.z,scale.z,2);
+#else
+        init_dim(node.lower_x,node.upper_x,lower_x,upper_x,start.x,scale.x,0);
+        init_dim(node.lower_y,node.upper_y,lower_y,upper_y,start.y,scale.y,1);
+        init_dim(node.lower_z,node.upper_z,lower_z,upper_z,start.z,scale.z,2);
+#endif
+
       }
 
+#ifdef COMPACT10BIT
+      __forceinline vfloat<N> dequantizeLowerX() const { return madd(vfloat<N>((lower_xyz >>  0) & 1023),scale.x,vfloat<N>(start.x)); }
+      __forceinline vfloat<N> dequantizeLowerY() const { return madd(vfloat<N>((lower_xyz >> 10) & 1023),scale.y,vfloat<N>(start.y)); }
+      __forceinline vfloat<N> dequantizeLowerZ() const { return madd(vfloat<N>((lower_xyz >> 20) & 1023),scale.z,vfloat<N>(start.z)); }
+      __forceinline vfloat<N> dequantizeUpperX() const { return madd(vfloat<N>((upper_xyz >>  0) & 1023),scale.x,vfloat<N>(start.x)); }
+      __forceinline vfloat<N> dequantizeUpperY() const { return madd(vfloat<N>((upper_xyz >> 10) & 1023),scale.y,vfloat<N>(start.y)); }
+      __forceinline vfloat<N> dequantizeUpperZ() const { return madd(vfloat<N>((upper_xyz >> 20) & 1023),scale.z,vfloat<N>(start.z)); }
+#else
       __forceinline vfloat<N> dequantizeLowerX() const { return madd(vfloat<N>(vint<N>::load(lower_x)),scale.x,vfloat<N>(start.x)); }
-
-      __forceinline vfloat<N> dequantizeUpperX() const { return madd(vfloat<N>(vint<N>::load(upper_x)),scale.x,vfloat<N>(start.x)); }
-
       __forceinline vfloat<N> dequantizeLowerY() const { return madd(vfloat<N>(vint<N>::load(lower_y)),scale.y,vfloat<N>(start.y)); }
-
-      __forceinline vfloat<N> dequantizeUpperY() const { return madd(vfloat<N>(vint<N>::load(upper_y)),scale.y,vfloat<N>(start.y)); }
-
       __forceinline vfloat<N> dequantizeLowerZ() const { return madd(vfloat<N>(vint<N>::load(lower_z)),scale.z,vfloat<N>(start.z)); }
-
+      __forceinline vfloat<N> dequantizeUpperX() const { return madd(vfloat<N>(vint<N>::load(upper_x)),scale.x,vfloat<N>(start.x)); }
+      __forceinline vfloat<N> dequantizeUpperY() const { return madd(vfloat<N>(vint<N>::load(upper_y)),scale.y,vfloat<N>(start.y)); }
       __forceinline vfloat<N> dequantizeUpperZ() const { return madd(vfloat<N>(vint<N>::load(upper_z)),scale.z,vfloat<N>(start.z)); }
+#endif
 
       template <int M>
       __forceinline vfloat<M> dequantize(const size_t offset) const { return vfloat<M>(vint<M>::loadu(all_planes+offset)); }
 
 #if defined(__AVX512F__)
+
+#ifdef COMPACT10BIT
+      template <int M>
+        __forceinline vfloat<M> dequantizeLowerUpperX(const vint<M> &p) const { return madd(vfloat<M>(permute(vint<M>::loadu(&lower_xyz) & 1023        ,p)),scale.x,vfloat<M>(start.x)); }
+
+      template <int M>
+        __forceinline vfloat<M> dequantizeLowerUpperY(const vint<M> &p) const { return madd(vfloat<M>(permute((vint<M>::loadu(&lower_xyz) >> 10) & 1023,p)),scale.y,vfloat<M>(start.y)); }
+
+      template <int M>
+        __forceinline vfloat<M> dequantizeLowerUpperZ(const vint<M> &p) const { return madd(vfloat<M>(permute((vint<M>::loadu(&lower_xyz) >> 20)       ,p)),scale.z,vfloat<M>(start.z)); }      
+
+#else
       template <int M>
       __forceinline vfloat<M> dequantizeLowerUpperX(const vint<M> &p) const { return madd(vfloat<M>(permute(vint<M>::load(lower_x),p)),scale.x,vfloat<M>(start.x)); }
 
@@ -1290,6 +1332,8 @@ namespace embree
 
       template <int M>
       __forceinline vfloat<M> dequantizeLowerUpperZ(const vint<M> &p) const { return madd(vfloat<M>(permute(vint<M>::load(lower_z),p)),scale.z,vfloat<M>(start.z)); }      
+#endif
+
 #endif
 
       union {
@@ -1302,6 +1346,12 @@ namespace embree
           T upper_z[N]; //!< 8bit discretized Z dimension of upper bounds of all N children
         };
         T all_planes[6*N];
+#ifdef COMPACT10BIT
+        struct {
+          vint<N> lower_xyz;
+          vint<N> upper_xyz;
+        };
+#endif
       };
 
       Vec3f start;

kernels/bvh/bvh8_factory.cpp

@@ -138,6 +138,8 @@ namespace embree
   DECLARE_SYMBOL2(Accel::Intersector16,BVH8Triangle4vIntersector16HybridPluecker);
   DECLARE_SYMBOL2(Accel::Intersector16,BVH8Triangle4iIntersector16HybridPluecker);
 
+  DECLARE_SYMBOL2(Accel::Intersector16,QBVH8Triangle4Intersector16HybridMoellerNoFilter);
+
   DECLARE_SYMBOL2(Accel::Intersector16,BVH8Triangle4vMBIntersector16HybridMoeller);
   DECLARE_SYMBOL2(Accel::Intersector16,BVH8Triangle4iMBIntersector16HybridMoeller);
   DECLARE_SYMBOL2(Accel::Intersector16,BVH8Triangle4vMBIntersector16HybridPluecker);
@@ -394,6 +396,8 @@ namespace embree
     IF_ENABLED_TRIS(SELECT_SYMBOL_INIT_AVX512KNL_AVX512SKX(features,BVH8Triangle4vIntersector16HybridPluecker));
     IF_ENABLED_TRIS(SELECT_SYMBOL_INIT_AVX512KNL_AVX512SKX(features,BVH8Triangle4iIntersector16HybridPluecker));
 
+    IF_ENABLED_TRIS(SELECT_SYMBOL_INIT_AVX512KNL_AVX512SKX(features,QBVH8Triangle4Intersector16HybridMoellerNoFilter));
+
     IF_ENABLED_TRIS(SELECT_SYMBOL_INIT_AVX512KNL_AVX512SKX(features,BVH8Triangle4vMBIntersector16HybridMoeller));
     IF_ENABLED_TRIS(SELECT_SYMBOL_INIT_AVX512KNL_AVX512SKX(features,BVH8Triangle4iMBIntersector16HybridMoeller));
     IF_ENABLED_TRIS(SELECT_SYMBOL_INIT_AVX512KNL_AVX512SKX(features,BVH8Triangle4vMBIntersector16HybridPluecker));
@@ -836,6 +840,9 @@ namespace embree
     Accel::Intersectors intersectors;
     intersectors.ptr = bvh;
     intersectors.intersector1 = QBVH8Triangle4Intersector1Moeller();
+#if defined (EMBREE_RAY_PACKETS)
+    intersectors.intersector16 = QBVH8Triangle4Intersector16HybridMoellerNoFilter();
+#endif
     return intersectors;
   }
 

kernels/bvh/bvh8_factory.h

@@ -201,6 +201,8 @@ namespace embree
     DEFINE_SYMBOL2(Accel::Intersector16,BVH8Triangle4vMBIntersector16HybridPluecker);
     DEFINE_SYMBOL2(Accel::Intersector16,BVH8Triangle4iMBIntersector16HybridPluecker);
 
+    DEFINE_SYMBOL2(Accel::Intersector16,QBVH8Triangle4Intersector16HybridMoellerNoFilter);
+
     DEFINE_SYMBOL2(Accel::Intersector16,BVH8Quad4vIntersector16HybridMoeller);
     DEFINE_SYMBOL2(Accel::Intersector16,BVH8Quad4vIntersector16HybridMoellerNoFilter);
     DEFINE_SYMBOL2(Accel::Intersector16,BVH8Quad4iIntersector16HybridMoeller);

kernels/bvh/bvh_intersector_hybrid.cpp

@@ -33,7 +33,7 @@
 #include "../geometry/object_intersector.h"
 
 #define SWITCH_DURING_DOWN_TRAVERSAL 1
-#define FORCE_SINGLE_MODE 0
+#define FORCE_SINGLE_MODE 1
 
 #define ENABLE_FAST_COHERENT_CODEPATHS 1
 
@@ -103,6 +103,10 @@ namespace embree
 
           /* select next child and push other children */
           BVHNNodeTraverser1<N,Nx,types>::traverseClosestHit(cur,mask,tNear,stackPtr,stackEnd);
+#if defined(__AVX512ER__)
+          //BVHNNodeTraverser1Permute<N,Nx,types>::traverseClosestHit(cur,mask,tNear,stackPtr,stackEnd,vray.octantIndex);
+#endif
+
         }
 
         /*! this is a leaf node */

kernels/bvh/bvh_intersector_hybrid.h

@@ -29,7 +29,7 @@ namespace embree
     class BVHNIntersectorKHybrid
     {
       /* right now AVX512KNL SIMD extension only for standard node types */
-      static const size_t Nx = types == BVH_AN1 ? vextend<N>::size : N;
+      static const size_t Nx = (types == BVH_AN1 | types == BVH_QN1) ? vextend<N>::size : N;
 
       /* shortcuts for frequently used types */
       typedef typename PrimitiveIntersectorK::Precalculations Precalculations;

kernels/bvh/bvh_intersector_hybrid16_bvh8.cpp

@@ -30,6 +30,8 @@ namespace embree
     IF_ENABLED_TRIS(DEFINE_INTERSECTOR16(BVH8Triangle4vIntersector16HybridPluecker,      BVHNIntersectorKHybrid<8 COMMA 16 COMMA BVH_AN1 COMMA true  COMMA ArrayIntersectorK_1<16 COMMA TriangleMvIntersectorKPluecker<SIMD_MODE(4) COMMA 16 COMMA true> > >));
     IF_ENABLED_TRIS(DEFINE_INTERSECTOR16(BVH8Triangle4iIntersector16HybridPluecker,      BVHNIntersectorKHybrid<8 COMMA 16 COMMA BVH_AN1 COMMA true  COMMA ArrayIntersectorK_1<16 COMMA TriangleMiIntersectorKPluecker<SIMD_MODE(4) COMMA 16 COMMA true > > >));
 
+    IF_ENABLED_TRIS(DEFINE_INTERSECTOR16(QBVH8Triangle4Intersector16HybridMoellerNoFilter,BVHNIntersectorKHybrid<8 COMMA 16 COMMA BVH_QN1 COMMA false COMMA ArrayIntersectorK_1<16 COMMA TriangleMIntersectorKMoeller  <SIMD_MODE(4) COMMA 16 COMMA false> > >));
+
     IF_ENABLED_TRIS(DEFINE_INTERSECTOR16(BVH8Triangle4vMBIntersector16HybridMoeller, BVHNIntersectorKHybrid<8 COMMA 16 COMMA BVH_AN2_AN4D COMMA false COMMA ArrayIntersectorK_1<16 COMMA TriangleMvMBIntersectorKMoeller <SIMD_MODE(4) COMMA 16 COMMA true> > >));
     IF_ENABLED_TRIS(DEFINE_INTERSECTOR16(BVH8Triangle4iMBIntersector16HybridMoeller, BVHNIntersectorKHybrid<8 COMMA 16 COMMA BVH_AN2_AN4D COMMA false COMMA ArrayIntersectorK_1<16 COMMA TriangleMiMBIntersectorKMoeller <SIMD_MODE(4) COMMA 16 COMMA true> > >));
     IF_ENABLED_TRIS(DEFINE_INTERSECTOR16(BVH8Triangle4vMBIntersector16HybridPluecker,BVHNIntersectorKHybrid<8 COMMA 16 COMMA BVH_AN2_AN4D COMMA true  COMMA ArrayIntersectorK_1<16 COMMA TriangleMvMBIntersectorKPluecker<SIMD_MODE(4) COMMA 16 COMMA true> > >));

kernels/bvh/bvh_intersector_node.h

@@ -79,6 +79,11 @@ namespace embree
         permX = select(vfloat<16>(dir.x) >= 0.0f, id, id2);
         permY = select(vfloat<16>(dir.y) >= 0.0f, id, id2);
         permZ = select(vfloat<16>(dir.z) >= 0.0f, id, id2);
+
+        octantIndex = \
+          (any(vfloat<16>(dir.x) < 0.0f) ? (size_t)1 : 0) |
+          (any(vfloat<16>(dir.y) < 0.0f) ? (size_t)2 : 0) |
+          (any(vfloat<16>(dir.z) < 0.0f) ? (size_t)4 : 0);
 #endif
       }
 
@@ -116,6 +121,7 @@ namespace embree
 #endif
 #if defined(__AVX512ER__) // KNL+
       vint16 permX, permY, permZ;
+      size_t octantIndex;
 #endif
 
       size_t nearX, nearY, nearZ;
@@ -927,6 +933,50 @@ namespace embree
     
 #endif
 
+
+    template<int N, int K>
+      __forceinline vbool<K> intersectNode(const typename BVHN<N>::QuantizedNode* node, size_t i, 
+                                           const Vec3vf<K>& org, const Vec3vf<K>& dir, const Vec3vf<K>& rdir, const Vec3vf<K>& org_rdir,
+                                           const vfloat<K>& tnear, const vfloat<K>& tfar, vfloat<K>& dist)
+ 
+    {
+      const vfloat8 start_x(node->start.x);
+      const vfloat8 scale_x(node->scale.x);
+      const vfloat8 lower_x = madd(node->dequantizeLowerX(),scale_x,start_x);
+      const vfloat8 upper_x = madd(node->dequantizeUpperX(),scale_x,start_x);
+      const vfloat8 start_y(node->start.y);
+      const vfloat8 scale_y(node->scale.y);
+      const vfloat8 lower_y = madd(node->dequantizeLowerX(),scale_y,start_y);
+      const vfloat8 upper_y = madd(node->dequantizeUpperX(),scale_y,start_y);
+      const vfloat8 start_z(node->start.z);
+      const vfloat8 scale_z(node->scale.z);
+      const vfloat8 lower_z = madd(node->dequantizeLowerX(),scale_z,start_z);
+      const vfloat8 upper_z = madd(node->dequantizeUpperX(),scale_z,start_z);
+
+#if defined(__AVX2__)
+      const vfloat<K> lclipMinX = msub(lower_x[i],rdir.x,org_rdir.x);
+      const vfloat<K> lclipMinY = msub(lower_y[i],rdir.y,org_rdir.y);
+      const vfloat<K> lclipMinZ = msub(lower_z[i],rdir.z,org_rdir.z);
+      const vfloat<K> lclipMaxX = msub(upper_x[i],rdir.x,org_rdir.x);
+      const vfloat<K> lclipMaxY = msub(upper_y[i],rdir.y,org_rdir.y);
+      const vfloat<K> lclipMaxZ = msub(upper_z[i],rdir.z,org_rdir.z);
+#else
+      const vfloat<K> lclipMinX = (lower_x[i] - org.x) * rdir.x;
+      const vfloat<K> lclipMinY = (lower_y[i] - org.y) * rdir.y;
+      const vfloat<K> lclipMinZ = (lower_z[i] - org.z) * rdir.z;
+      const vfloat<K> lclipMaxX = (upper_x[i] - org.x) * rdir.x;
+      const vfloat<K> lclipMaxY = (upper_y[i] - org.y) * rdir.y;
+      const vfloat<K> lclipMaxZ = (upper_z[i] - org.z) * rdir.z;
+#endif  
+
+      /* use mixed float/int min/max */
+      const vfloat<K> lnearP = maxi(min(lclipMinX, lclipMaxX), min(lclipMinY, lclipMaxY), min(lclipMinZ, lclipMaxZ));
+      const vfloat<K> lfarP  = mini(max(lclipMinX, lclipMaxX), max(lclipMinY, lclipMaxY), max(lclipMinZ, lclipMaxZ));
+      const vbool<K> lhit    = asInt(maxi(lnearP, tnear)) <= asInt(mini(lfarP, tfar));
+      dist = lnearP;
+      return lhit;
+    }
+
     //////////////////////////////////////////////////////////////////////////////////////
     // fast ray/BVHN::UnalignedNode intersection
     //////////////////////////////////////////////////////////////////////////////////////
@@ -1321,5 +1371,18 @@ namespace embree
         return true;
       }
     };
+
+    template<int N, int K>
+    struct BVHNNodeIntersectorK<N,K,BVH_QN1,false>
+    {
+      static __forceinline bool intersect(const typename BVHN<N>::NodeRef& node, const size_t i, const Vec3vf<K>& org, const Vec3vf<K>& dir, const Vec3vf<K>& rdir, const Vec3vf<K>& org_rdir,
+                                          const vfloat<K>& tnear, const vfloat<K>& tfar, const vfloat<K>& time, vfloat<K>& dist, vbool<K>& vmask)
+      {
+        vmask = intersectNode<N,K>(node.quantizedNode(),i,org,dir,rdir,org_rdir,tnear,tfar,dist);
+        return true;
+      }
+    };
+
+
   }
 }

kernels/bvh/bvh_traverser1.h

@@ -242,6 +242,7 @@ namespace embree
         }
         cur = toScalar(ptr);
 #endif        
+
       }
 
 #endif
@@ -679,5 +680,91 @@ namespace embree
     public:
       __forceinline explicit BVHNNodeTraverser1(const TravRay<N,Nx>& vray) : BVHNNodeTraverser1Transform<N, Nx, types, (bool)(types & BVH_FLAG_TRANSFORM_NODE)>(vray) {}
     };
+
+
+
+#if 0
+    // experimental stuff will be removed soon
+    /*! BVH regular node traversal for single rays. */
+    template<int N, int Nx, int types>
+    class BVHNNodeTraverser1Permute;
+
+    template<int Nx, int types>
+      class BVHNNodeTraverser1Permute<4,Nx,types>
+    {
+      typedef BVH4 BVH;
+      typedef BVH4::NodeRef NodeRef;
+      typedef BVH4::BaseNode BaseNode;
+
+    public:
+      static __forceinline void traverseClosestHit(NodeRef& cur,
+                                                   size_t mask,
+                                                   const vfloat<Nx>& tNear,
+                                                   StackItemT<NodeRef>*& stackPtr,
+                                                   StackItemT<NodeRef>* stackEnd,
+                                                   const unsigned int octantIndex)
+      {
+      }
+    };
+
+    template<int Nx, int types>
+      class BVHNNodeTraverser1Permute<8,Nx,types>
+    {
+      typedef BVH8 BVH;
+      typedef BVH8::NodeRef NodeRef;
+      typedef BVH8::BaseNode BaseNode;
+
+    public:
+      static __forceinline void traverseClosestHit(NodeRef& cur,
+                                                   size_t mask,
+                                                   const vfloat<Nx>& tNear,
+                                                   StackItemT<NodeRef>*& stackPtr,
+                                                   StackItemT<NodeRef>* stackEnd,
+                                                   const unsigned int octantIndex)
+      {
+#if defined(__AVX512ER__)   
+        assert(mask);
+        STAT3(normal.trav_hit_boxes[__popcnt(mask)],1,1,1);
+        const BVH::QuantizedNode* node = cur.quantizedNode();
+        const unsigned int gatherIDs = node->lookupTable[octantIndex];
+        PRINT(gatherIDs);
+        vuint16 gatherIDs_dist(gatherIDs);
+        gatherIDs_dist = gatherIDs_dist >> (vuint16( step ) + vuint16( step ) + vuint16( step ));
+        gatherIDs_dist &= 0x7;
+        vllong8 gatherIDs_nodes(zeroExtend32Bit(gatherIDs_dist));
+        vfloat16 dist = select(vbool16((int)mask),tNear,vfloat16(pos_inf));
+        dist = permute(dist,gatherIDs_dist);
+        vbool16 m_compact = dist != vfloat16(pos_inf);
+        vllong8 nodes = permute(vllong8::loadu(node->children),gatherIDs_nodes);
+        
+        dist  = vfloat16::compact(m_compact,dist);
+        nodes = vllong8::compact(vboold8(m_compact),nodes);
+
+        while (1)
+        {
+          assert(stackPtr < stackEnd);
+          //size_t r = __bscf(mask);
+          cur = toScalar(nodes);
+          //cur = node->child(r);
+          //assert(cur == node->child(r));
+          nodes = align_shift_right<1>(nodes,nodes);
+          cur.prefetch(types); 
+          //unsigned int d = *(unsigned int*)&tNear[r]; 
+          float dd = toScalar(dist);
+          //assert(dd == tNear[r]);
+          dist = align_shift_right<1>(dist,dist);
+
+          stackPtr->ptr = cur; 
+          *(float*)&stackPtr->dist = dd; 
+          stackPtr++;
+          assert(cur != BVH::emptyNode);
+          if (unlikely(mask == 0)) break;
+        }
+        stackPtr--;
+#endif
+      }
+    };
+#endif
+
   }
 }

tutorials/pathtracer/pathtracer_device.cpp

@@ -40,9 +40,28 @@ namespace embree {
 #define LEVEL_FACTOR    64.0f
 #define MAX_PATH_LENGTH  8
 
+#define FILTER_TEST 1
+
 bool g_subdiv_mode = false;
 unsigned int keyframeID = 0;
 
+struct GBufferEntry {
+  GBufferEntry() {}
+
+  Vec3fa diffuse_brdf;
+  Vec3fa direct_illum;
+  Vec3fa indirect_illum;
+  float distance;
+
+  __forceinline void clear() {
+    diffuse_brdf = Vec3fa(0.0f);
+    direct_illum = Vec3fa(0.0f);
+    indirect_illum = Vec3fa(0.0f);
+    distance = inf;
+  };
+
+};
+
 struct BRDF
 {
   float Ns;               /*< specular exponent */
@@ -885,6 +904,10 @@ Vec3fa g_accu_p;
 extern "C" bool g_changed;
 extern "C" int g_instancing_mode;
 
+/* g-buffer */
+GBufferEntry *gbuffer = nullptr;
+unsigned int gbuffer_width = 0;
+unsigned int gbuffer_height = 0;
 
 bool g_animation = true;
 bool g_use_smooth_normals = false;
@@ -1485,6 +1508,151 @@ Vec3fa renderPixelFunction(float x, float y, RandomSampler& sampler, const ISPCC
   return L;
 }
 
+
+
+  GBufferEntry renderGBufferFunction(float x, float y, RandomSampler& sampler, const ISPCCamera& camera, RayStats& stats)
+{
+  GBufferEntry gb;
+  gb.clear();
+
+  /* radiance accumulator and weight */
+  Vec3fa L = Vec3fa(0.0f);
+  Vec3fa Lw = Vec3fa(1.0f);
+  Medium medium = make_Medium_Vacuum();
+  float time = RandomSampler_get1D(sampler);
+
+  /* initialize ray */
+  RTCRay ray = RTCRay(Vec3fa(camera.xfm.p),
+                      Vec3fa(normalize(x*camera.xfm.l.vx + y*camera.xfm.l.vy + camera.xfm.l.vz)),0.0f,inf,time);
+
+  DifferentialGeometry dg;
+
+  /* iterative path tracer loop */
+  for (int i=0; i<MAX_PATH_LENGTH; i++)
+  {
+    /* terminate if contribution too low */
+    if (max(Lw.x,max(Lw.y,Lw.z)) < 0.01f)
+      break;
+
+    /* intersect ray with scene */
+    RTCIntersectContext context;
+    context.flags = (i == 0) ? g_iflags_coherent : g_iflags_incoherent;
+    rtcIntersect1Ex(g_scene,&context,ray);
+    RayStats_addRay(stats);
+    const Vec3fa wo = neg(ray.dir);
+
+    /* invoke environment lights if nothing hit */
+    if (ray.geomID == RTC_INVALID_GEOMETRY_ID)
+    {
+      //L = L + Lw*Vec3fa(1.0f);
+
+      /* iterate over all lights */
+      for (size_t i=0; i<g_ispc_scene->numLights; i++)
+      {
+        const Light* l = g_ispc_scene->lights[i];
+        Light_EvalRes le = l->eval(l,dg,ray.dir);
+        L = L + Lw*le.value;
+      }
+
+      break;
+    }
+    Vec3fa Ns = normalize(ray.Ng);
+
+    if (g_use_smooth_normals)
+      if (ray.geomID != RTC_INVALID_GEOMETRY_ID) // FIXME: workaround for ISPC bug, location reached with empty execution mask
+    {
+      Vec3fa dPdu,dPdv;
+      unsigned int geomID = ray.geomID; {
+        rtcInterpolate(g_scene,geomID,ray.primID,ray.u,ray.v,RTC_VERTEX_BUFFER0,nullptr,&dPdu.x,&dPdv.x,3);
+      }
+      Ns = normalize(cross(dPdv,dPdu));
+    }
+
+    /* compute differential geometry */
+    dg.geomID = ray.geomID;
+    dg.primID = ray.primID;
+    dg.u = ray.u;
+    dg.v = ray.v;
+    dg.P  = ray.org+ray.tfar*ray.dir;
+    dg.Ng = ray.Ng;
+    dg.Ns = Ns;
+    int materialID = postIntersect(ray,dg);
+    dg.Ng = face_forward(ray.dir,normalize(dg.Ng));
+    dg.Ns = face_forward(ray.dir,normalize(dg.Ns));
+
+    /*! Compute  simple volumetric effect. */
+    Vec3fa c = Vec3fa(1.0f);
+    const Vec3fa transmission = medium.transmission;
+    if (ne(transmission,Vec3fa(1.0f)))
+      c = c * pow(transmission,ray.tfar);
+
+    /* calculate BRDF */
+    BRDF brdf;
+    int numMaterials = g_ispc_scene->numMaterials;
+    ISPCMaterial** material_array = &g_ispc_scene->materials[0];
+    Material__preprocess(material_array,materialID,numMaterials,brdf,wo,dg,medium);
+
+    /* sample BRDF at hit point */
+    Sample3f wi1;
+    c = c * Material__sample(material_array,materialID,numMaterials,brdf,Lw, wo, dg, wi1, medium, RandomSampler_get2D(sampler));
+
+    /* fill gbuffer */
+    if (i == 0)
+    {
+      gb.distance = ray.tfar;      
+    }
+
+    /* iterate over lights */
+    for (size_t i=0; i<g_ispc_scene->numLights; i++)
+    {
+      const Light* l = g_ispc_scene->lights[i];
+      Light_SampleRes ls = l->sample(l,dg,RandomSampler_get2D(sampler));
+      if (ls.pdf <= 0.0f) continue;
+      RTCRay shadow = RTCRay(dg.P,ls.dir,dg.tnear_eps,ls.dist,time); shadow.transparency = Vec3fa(1.0f);
+      rtcOccluded1Ex(g_scene,&context,shadow);
+      RayStats_addShadowRay(stats);
+      //if (shadow.geomID != RTC_INVALID_GEOMETRY_ID) continue;
+      const Vec3fa mat_eval = Material__eval(material_array,materialID,numMaterials,brdf,wo,dg,ls.dir);
+
+      /* fill gbuffer */
+      if (i == 0)
+      {
+        gb.diffuse_brdf += mat_eval;
+      }
+
+      if (max(max(shadow.transparency.x,shadow.transparency.y),shadow.transparency.z) > 0.0f)
+      {
+        const Vec3fa illum = Lw*ls.weight*shadow.transparency;
+        L = L + illum*mat_eval;
+
+        /* fill gbuffer */
+        if (i == 0)
+        {
+          gb.direct_illum += illum;
+        }          
+        else
+        {
+          gb.indirect_illum += illum*mat_eval;          
+        }
+      }
+    }
+    if (i == 0)
+    {
+      gb.diffuse_brdf *= 1.0f / (float)g_ispc_scene->numLights;
+    }
+
+    if (wi1.pdf <= 1E-4f /* 0.0f */) break;
+    Lw = Lw*c/wi1.pdf;
+
+    /* setup secondary ray */
+    float sign = dot(wi1.v,dg.Ng) < 0.0f ? -1.0f : 1.0f;
+    dg.P = dg.P + sign*dg.tnear_eps*dg.Ng;
+    ray = RTCRay(dg.P,normalize(wi1.v),dg.tnear_eps,inf,time);
+  }
+  return gb;
+}
+
+
 /* task that renders a single screen tile */
 Vec3fa renderPixelStandard(float x, float y, const ISPCCamera& camera, RayStats& stats)
 {
@@ -1526,8 +1694,15 @@ void renderTileStandard(int taskIndex,
   for (unsigned int y=y0; y<y1; y++) for (unsigned int x=x0; x<x1; x++)
   {
     /* calculate pixel color */
+#if FILTER_TEST == 1
+    RandomSampler sampler;
+    RandomSampler_init(sampler, (int)x, (int)y, 0);
+    /* calculate pixel color */
+    float fx = x + RandomSampler_get1D(sampler);
+    float fy = y + RandomSampler_get1D(sampler);
+    gbuffer[y*width+x] = renderGBufferFunction(fx,fy,sampler,camera,g_stats[threadIndex]);
+#else
     Vec3fa color = renderPixelStandard((float)x,(float)y,camera,g_stats[threadIndex]);
-
     /* write color to framebuffer */
     Vec3fa accu_color = g_accu[y*width+x] + Vec3fa(color.x,color.y,color.z,1.0f); g_accu[y*width+x] = accu_color;
     float f = rcp(max(0.001f,accu_color.w));
@@ -1535,6 +1710,7 @@ void renderTileStandard(int taskIndex,
     unsigned int g = (unsigned int) (255.0f * clamp(accu_color.y*f,0.0f,1.0f));
     unsigned int b = (unsigned int) (255.0f * clamp(accu_color.z*f,0.0f,1.0f));
     pixels[y*width+x] = (b << 16) + (g << 8) + r;
+#endif
   }
 }
 
@@ -1638,12 +1814,61 @@ extern "C" void device_init (char* cfg)
 
 } // device_init
 
+  __forceinline float distance(int x, int y, int i, int j) {
+    const float dx = (float)(x-i);
+    const float dy = (float)(y-j);
+    return sqrtf(dx*dx + dy*dy);
+  }
+
+  template<class T>
+  __forceinline T gaussian(T x, float sigma) {
+    return exp(-(x*x)/(2.0f * sigma * sigma)) / (2.0f * M_PI * sigma * sigma);    
+  }
+
+  __forceinline Vec3fa getColor(int x, int y) {
+    return gbuffer[y*g_accu_width+x].direct_illum;
+  }
+Vec3fa applyBilateralFilter(int x, int y, int diameter, float sigmaI, float sigmaS) {
+  Vec3fa color(zero);
+  const int half = diameter / 2;
+
+#if 1
+  float wP = 0;
+  for(int i = 0; i < diameter; i++) {
+    for(int j = 0; j < diameter; j++) {
+      const int neighbor_x = x - (half - i);
+      const int neighbor_y = y - (half - j);
+      const float gi = gaussian(length(getColor(neighbor_x, neighbor_y) - getColor(x, y)), sigmaI);
+      const float gs = gaussian(distance(x, y, neighbor_x, neighbor_y), sigmaS);
+      const float w = gi * gs;
+      color += getColor(neighbor_x, neighbor_y) * w;
+      wP += w;
+    }
+  }
+  color *= 1.0f / wP;
+#else
+  int num = 0;
+  for(int i = 0; i < diameter; i++) {
+    for(int j = 0; j < diameter; j++) {
+      const int neighbor_x = x - (half - i);
+      const int neighbor_y = y - (half - j);
+      color += getColor(neighbor_x, neighbor_y);
+      num++;
+    }
+  }
+  color *= 1.0f / (float)num;
+
+#endif
+  return color;
+}
+
+
 /* called by the C++ code to render */
 extern "C" void device_render (int* pixels,
-                           const unsigned int width,
-                           const unsigned int height,
-                           const float time,
-                           const ISPCCamera& camera)
+                               const unsigned int width,
+                               const unsigned int height,
+                               const float time,
+                               const ISPCCamera& camera)
 {
   /* create scene */
   if (g_scene == nullptr) {
@@ -1662,6 +1887,14 @@ extern "C" void device_render (int* pixels,
       g_accu[i] = Vec3fa(0.0f);
   }
 
+  /* create gbuffer */
+  if (gbuffer_width != width || gbuffer_height != height) {
+    if (gbuffer) alignedFree(gbuffer);
+    gbuffer = (GBufferEntry*) alignedMalloc(width*height*sizeof(GBufferEntry));
+    gbuffer_width = width;
+    gbuffer_height = height;
+  }
+
   /* reset accumulator */
   bool camera_changed = g_changed || !g_accumulate; g_changed = false;
   camera_changed |= ne(g_accu_vx,camera.xfm.l.vx); g_accu_vx = camera.xfm.l.vx;
@@ -1680,8 +1913,10 @@ extern "C" void device_render (int* pixels,
       rtcCommit (g_scene);
     }
   }
+#if FILTER_TEST == 0
   else
     g_accu_count++;
+#endif
 
   /* render image */
   const int numTilesX = (width +TILE_SIZE_X-1)/TILE_SIZE_X;
@@ -1691,6 +1926,35 @@ extern "C" void device_render (int* pixels,
     for (size_t i=range.begin(); i<range.end(); i++)
       renderTileTask((int)i,threadIndex,pixels,width,height,time,camera,numTilesX,numTilesY);
   }); 
+
+#if FILTER_TEST == 1
+  // filter accumulation buffer
+  const int diameter = 4;
+  const float sigmaI = 12.0f;
+  const float sigmaS = 12.0f;
+
+  for(int y = diameter; y < gbuffer_height - diameter; y++) 
+    for(int x = diameter; x < gbuffer_width - diameter; x++) 
+    {
+#if 1
+      GBufferEntry &gb = gbuffer[y*gbuffer_width+x];
+      Vec3fa direct_illum = applyBilateralFilter(x, y, diameter, sigmaI, sigmaS);
+      Vec3fa color = gb.diffuse_brdf * direct_illum; // + gb.indirect_illum;
+      unsigned int r = (unsigned int) (255.0f * clamp(color.x,0.0f,1.0f));
+      unsigned int g = (unsigned int) (255.0f * clamp(color.y,0.0f,1.0f));
+      unsigned int b = (unsigned int) (255.0f * clamp(color.z,0.0f,1.0f));
+      pixels[y*gbuffer_width+x] = (b << 16) + (g << 8) + r;
+
+#else
+      Vec3fa color = applyBilateralFilter(x, y, diameter, sigmaI, sigmaS);
+      unsigned int r = (unsigned int) (255.0f * clamp(color.x,0.0f,1.0f));
+      unsigned int g = (unsigned int) (255.0f * clamp(color.y,0.0f,1.0f));
+      unsigned int b = (unsigned int) (255.0f * clamp(color.z,0.0f,1.0f));
+      pixels[y*g_accu_width+x] = (b << 16) + (g << 8) + r;
+#endif
+    }
+#endif
+  
   //rtcDebug();
 } // device_render